#  YITP-22-66
Splitting-merging transitions in tensor-vectors systems
in exact large-$N$ limits

Naoki<EMAIL_ADDRESS>
Yukawa Institute for Theoretical Physics, Kyoto University,
and
CGPQI, Yukawa Institute for Theoretical Physics, Kyoto University,
Kitashirakawa, Sakyo-ku, Kyoto 606-8502, Japan

Matrix models have phase transitions in which distributions of variables
change topologically like the Gross-Witten-Wadia transition. In a recent
study, similar splitting-merging behavior of distributions of dynamical
variables was observed in a tensor-vectors system by numerical simulations. In
this paper, we study the system exactly in some large-$N$ limits, in which the
distributions are discrete sets of configurations rather than continuous. We
find cascades of first-order phase transitions for fixed tensors, and first-
and second-order phase transitions for random tensors, being characterized by
breaking patterns of replica symmetries. The system is of interest across
three different subjects at least: The splitting dynamics plays essential
roles in emergence of classical spacetimes in a tensor model of quantum
gravity; The splitting dynamics automatically detects the rank of a tensor in
the tensor rank decomposition in data analysis; The system provides a variant
of the spherical $p$-spin model for spin glasses with a new non-trivial
parameter. We discuss some implications of the results from these
perspectives. The results are compared with some numerical simulations to
check the large-$N$ convergence and the assumptions made in the analysis.

##  1 Introduction

The purpose of this paper is to better understand the dynamics of the
dynamical system defined by the partition function,

$\displaystyle
Z_{N,R}(\beta,C)=\int_{I}d\phi\,e^{-\beta(C-\phi\phi\phi)^{2}},$ (1)

where $\beta$ is an inverse temperature, $C$ denotes a real symmetric tensor
of order three, $C_{abc}=C_{bac}=C_{bca}\ (a,b,c=1,2,\ldots,N)$, the
integration variables are $\phi_{a}^{i}\ (a=1,2,\ldots,N,\ i=1,2,\ldots,R)$,
the integration region is $I=\mathbb{R}^{NR}$, the integration measure is
$d\phi=\prod_{a=1}^{N}\prod_{i=1}^{R}d\phi_{a}^{i}$, and we use a short-hand
notation,

$\displaystyle(C-\phi\phi\phi)^{2}=\left(C_{abc}-\sum_{i=1}^{R}\phi_{a}^{i}\phi_{b}^{i}\phi_{c}^{i}\right)\left(C_{abc}-\sum_{i=1}^{R}\phi_{a}^{i}\phi_{b}^{i}\phi_{c}^{i}\right),$
(2)

where pairwise repeated lower indices are assumed to be summed over, as is
assumed throughout this paper. On the other hand summations over the upper
indices have to be always explicitly indicated. Note that, in addition to an
$O(N)$ symmetry with respect to the lower indices, the system (1) is invariant
under relabeling the upper index of $\phi_{a}^{i}$ (namely, invariant under
$\phi^{i}_{a}\leftrightarrow\phi^{j}_{a}$). In this paper, we call it real
replica symmetry, distinguishing it from the replica symmetry which appears
later in the replica trick.

One of the motivations to study the system (1) comes from a tensor model in
the Hamiltonian formalism, which we call canonical tensor model (CTM) [1, 2].
Tensor models were originally introduced as a generalization of the matrix
models, which are successful in describing two-dimensional quantum gravity,
with a hope to extend the success to higher dimensions [3, 4, 5, 6]. However,
these tensor models do not generate macroscopic spacetimes, suffering from
dominance of singular objects like branched polymers [7, 8]. A motivation of
considering CTM is to overcome the issue by introducing a temporal direction
into tensor models, trying to follow the success of the causal dynamical
triangulation over the dynamical triangulation in generating macroscopic
spacetimes [9], where the former has a temporal direction, while the latter
not. Indeed, in [10], it was explicitly shown by numerical simulations that a
wave function of CTM seems to have a twofold phase structure, and classical
spacetimes emerge in the one which we call the classical phase.222The other
was called the quantum phase. Here the transition between the two phases can
be characterized by splitting-merging transitions of distributions of the
dynamical variables, like those in the matrix counterparts, such as the Gross-
Witten-Wadia transition [11, 12] and the transitions among multi-cut large-$N$
solutions [13]. However, the simulations were not convincing enough to
conclude whether the phases are really different (in some large-$N$ limits),
not just separated by crossovers. The system (1) is a part of this wave
function333The wave function has an integral expression, which is a multi-
variable extension of the integral representation of Airy functions. Its
integrand is complex, but, by taking only the modus of the integrand, one
obtains the system (1). In [10] the wave function, of complex values, was
analyzed by the re-weighting method of the Monte Carlo simulations, in which
the system (1) played the role of a statistical system with a positive weight
in the method., and plays essential roles in the dynamics of the above
transition. A result of this paper is that, at least in the large-$N$ limits
we consider, the system (1) has sharp phase transitions characterized by
breaking patterns of the real/genuine replica symmetries. The phases which
appear at large values of $\beta$ correspond to the classical phase found in
the previous paper.444 In fact, we will find cascades of first-order phase
transitions in this paper, meaning that what was called the classical phase in
[10] could be a collection of phases in general. This, however, must be taken
with caution, since the large-$N$ limits we take in this paper is different
from that in the tensor model, which should be $R\propto N^{2}$ (See the last
section).

Another motivation of studying (1) comes from that the exponent in (1) can be
used as a cost function of the tensor rank decomposition [14, 15, 16, 17],
which is an important technique in data analysis.555See [18] for some
developments to data analysis techniques from random tensor studies. In the
present case of a real symmetric tensor $C$ of order three, a real symmetric
tensor rank decomposition is defined to find $\phi_{a}^{i}\in\mathbb{R}$
satsfying

$\displaystyle C_{abc}=\sum_{i=1}^{R}\phi_{a}^{i}\phi_{b}^{i}\phi_{c}^{i},$
(3)

which is equivalent to require the cost function to vanish. The minimum value
of $R$ which realizes this decomposition for a $C$ is called the rank of $C$.
The decomposition (3) is a sort of an extension of the singular value
decomposition of the matrix to the tensor, but the hardness is largely
different [19]. While a matrix can be decomposed by straightforward
procedures, a practical method for the tensor case is to optimize
$\phi_{a}^{i}$ so that the cost function $(C-\phi\phi\phi)^{2}$ be minimized
or vanish [20]. Here one of the difficult issues is that we do not have prior
knowledge of an appropriate value of $R$ for a given $C$: If we take a too
large $R$ for the optimization, the decomposition will be overfitting, and, if
a too small $R$ is taken, we will miss some properties of $C$. What was found
in [10] and is interesting in the system (1) is that, when $\beta$ is taken
large enough, the dominant configurations of $\phi_{a}^{i}$ are such that they
are separated into two parts, the dominant and minor parts666Here the upper
index of $\phi_{a}^{i}$ has been relabeled without loss of generality.,

$\displaystyle\begin{split}C_{abc}&=\sum_{i=1}^{R_{c}}\phi_{a}^{i}\phi_{b}^{i}\phi_{c}^{i}+\sum_{i=R_{c}+1}^{R}\phi_{a}^{i}\phi_{b}^{i}\phi_{c}^{i},\\\
\phi_{a}^{i}&\not\sim 0\ (i=1,2,\ldots,R_{c}),\ \ \ \phi_{a}^{i}\sim 0\
(i=R_{c}+1,R_{c}+2,\ldots,R).\end{split}$ (4)

In fact, in the examples considered in [10], the values of $R_{c}$ agreed with
the ranks of $C$ (or very near values in a few large system cases). Therefore
the system (1) seems to have an intrinsic dynamics which automatically detects
an appropriate rank for a given $C$. A purpose of this paper is to study this
interesting property, which was found in the previous numerical simulation, by
an exact method. In the long run, understanding the system (1) would provide
some solutions to the long-standing issues in the tensor rank decomposition.

Another interesting link of the system (1) to physics is spin glasses. It can
be regarded as a variant of spherical $p$-spin model for spin glasses [21,
22], which is defined by the following Hamiltonian and a constraint,

$\displaystyle\begin{split}&H=C_{a_{1}a_{2}\cdots
a_{p}}\phi_{a_{1}}\phi_{a_{2}}\cdots\phi_{a_{p}},\\\
&\phi_{a}\phi_{a}=N,\end{split}$ (5)

where the tensor $C_{a_{1}a_{2}\cdots a_{p}}$ is assumed to take random
numbers777The simplest realization would be a normal distribution.. Thus the
system (1) can be regarded as a multi-real-replica extension of the spherical
$p$-spin model. Note that we are interested in finite $R$, while a replica
number is taken to vanish in the replica trick. Considering the connection to
the tensor rank decomposition, it should be a non-trivial question how (1)
behaves in $R$.

The non-triviality of the $R$-dependence can also be seen in another way. Note
that our system (1) does not have the constraint in (5), which prohibits
$\phi_{a}^{i}$ to run away to infinity, assuring the stability of the model.
Though the exponent in (1) is semidefinite in our case, it is a non-trivial
question whether (1) is finite or not, because the exponent contains flat
directions, such as $\phi_{a}^{i}=-\phi_{a}^{j}$, which extends to infinity.
This question about the finiteness was systematically studied mainly by
numerical methods in [23]888There are also some closely related studies [24,
25, 26]., and it was checked/conjectured that the system is finite only for
$R\lesssim(N+1)(N+2)/2$. In this paper, we are free from this instability,
because we consider large-$N$ limits with finite $R$.

This paper is organized as follows. In Section 2, we introduce a large-$N$
limit with finite $R$ of the system (1) for fixed $C$, and derive an exact
expression of the free energy in the limit. In Section 3, we consider a few
representative examples of $C$, and study the free energy by explicitly
computing minima of the expression obtained in Section 2. As $\beta$ is
increased, the system undergoes a cascade of first-order phase transitions, in
which the number of non-zero $\psi^{i}$, which appear later, increases one by
one, eventually reaches a phase in which the number agrees with the tensor
rank of $C$, and stays there. In Section 4, we consider the case of random
$C$. To incorporate the random $C$, we employ the replica trick up to one-step
replica symmetry breaking (1RSB) as in the case of the spherical $p$-spin
model. When $\beta$ is small, the system is in the replica symmetric phase. As
$\beta$ is increased, the system undergoes a first- or second-order phase
transition, depending on a parameter we introduce, and enters the 1RSB spin-
glass phase. In Section 5, we perform some numerical simulations and compare
with the exact results.

##  2 Large-$N$ limits with finite $R$ and fixed $C$

The strategy of this section is that we dip $C$ of a finite dimension into a
large-$N$ system. More precisely, let us introduce a parameter $n$, and assume
that $C$ takes non-zero values only within this sub-dimension:

$\displaystyle\begin{split}&C_{a_{1}a_{2}a_{3}}\neq 0,\ \hbox{only if}\
\forall a_{i}\in\\{1,2,\ldots,n\\},\\\
&C_{a_{1}a_{2}a_{3}}=0,\hbox{otherwise}.\end{split}$ (6)

We assume that $n$ is kept finite in the large-$N$ limit (However, we will
later consider $n\propto N$ in the large-$N$ limit for random $C$). Let us
introduce new variables to separate $\phi_{a}^{i}$ into two parts:

$\displaystyle\begin{split}&\psi_{a}^{i}=\phi_{a}^{i}\ (a=1,2,\ldots,n),\\\
&\tilde{\psi}_{a}^{i}=\phi_{a}^{i}\ (a=n+1,n+2,\ldots,N),\end{split}$ (7)

for $\forall i=1,2,\ldots,R$. With these variables, the exponent
$\beta(C-\phi\phi\phi)^{2}$ of (1) can be rewritten as

$\displaystyle\begin{split}S_{\psi\tilde{\psi}}(\beta,C)&=\beta(C-\psi\psi\psi)^{2}+3\beta\sum_{i,j=1}^{R}(\psi^{i}\cdot\psi^{j})^{2}\,\tilde{\psi}^{i}\cdot\tilde{\psi}^{j}\\\
&\ \ \
+3\beta\sum_{i,j=1}^{R}\psi^{i}\cdot\psi^{j}\,(\tilde{\psi}^{i}\cdot\tilde{\psi}^{j})^{2}+\beta\sum_{i,j=1}^{R}(\tilde{\psi}^{i}\cdot\tilde{\psi}^{j})^{3},\end{split}$
(8)

where $\cdot$ denotes the inner product,
$\psi^{i}\cdot\psi^{j}=\psi_{a}^{i}\psi_{a}^{j}$, and $C$ abusively denotes
the $n$-sub-dimensional part. $S_{\psi\tilde{\psi}}(\beta,C)$ is invariant
under the $SO(N-n)$ transformation with respect to the lower index of
$\tilde{\psi}_{a}^{i}$, and we can factor out the degrees of freedom by
introducing the following new variable, which is called overlap in spin glass
theory [22],

$\displaystyle\tilde{Q}_{ij}=\tilde{\psi}^{i}\cdot\tilde{\psi}^{j},\
(i,j=1,2,\ldots,R).$ (9)

This constraint can be embedded into the system by rewriting the partition
function (1) as

$\displaystyle Z_{N,R}(\beta,C)=\int d\psi d\tilde{\psi}d\lambda
d\tilde{Q}\,e^{-S_{\psi\tilde{Q}}(\beta,C)+i\sum_{i,j=1}^{R}\lambda_{ij}(\tilde{Q}_{ij}-\tilde{\psi}^{i}\cdot\tilde{\psi}^{j})},$
(10)

where we have ignored an irrelevant overall factor, and

$\displaystyle\begin{split}S_{\psi\tilde{Q}}(\beta,C)=\beta(C-\psi\psi\psi)^{2}+3\beta\sum_{i,j=1}^{R}(\psi^{i}\cdot\psi^{j})^{2}\tilde{Q}_{ij}+3\beta\sum_{i,j=1}^{R}\psi^{i}\cdot\psi^{j}\,(\tilde{Q}_{ij})^{2}+\beta\sum_{i,j=1}^{R}(\tilde{Q}_{ij})^{3}.\end{split}$
(11)

Integrating over $\tilde{\psi}$ in(10) generates a new term $-(N-n)/2\
\log\mathrm{det}\lambda$ in the exponent, and then by assuming large-$N$ and
carrying out the $\lambda$ integration by taking the saddle point, we obtain

$\displaystyle S^{\rm
eff}_{\psi\tilde{Q}}(\bar{\beta},C)=(N-n)\left(S_{\psi\tilde{Q}}(\bar{\beta},C)-\frac{1}{2}\ln\mathrm{det}\tilde{Q}\right).$
(12)

where we have introduced $\bar{\beta}$ by $\beta=(N-n)\,\bar{\beta}$.
Therefore, in the large-$N$ limit, the free energy of the system is given by

$\displaystyle\begin{split}\bar{\beta}F(\bar{\beta},C)&=-\lim_{N\rightarrow\infty}\frac{1}{N}\log
Z(\beta,C)\\\ &=\min_{\psi,\tilde{Q}\geq
0}\left(S_{\psi\tilde{Q}}(\bar{\beta},C)-\frac{1}{2}\ln\mathrm{det}\tilde{Q}\right),\end{split}$
(13)

where $\tilde{Q}\geq 0$ represents that $\tilde{Q}$ is constrained to be a
positive semidefinite matrix999Namely, all the eigenvalues are zero or
positive. due to (9).

In general, when the parameters $\bar{\beta}$ and $C$ are gradually changed,
the free energy (13) will undergo various first-order phase transitions with
finite jumps of the minimum in $\psi,\tilde{Q}$. As will explicitly be shown
in some representative cases in Section 3, the splitting between the dominant
and minor parts (4) occurs as the results of the first-order phase
transitions.

##  3 Some representative examples for fixed $C$

In this section, we explicitly compute (13) for some $C$’s as a function of
$\bar{\beta}$. The minimization must be performed under the condition that
$\tilde{Q}$ is a positive semidefinite matrix. A convenient way to implement
it is to parameterize $\tilde{Q}$ as

$\displaystyle\tilde{Q}_{ij}=v^{i}\cdot v^{j},$ (14)

where $v^{i}\ (i=1,2,\ldots,R)$ are $R$-dimensional real vectors which have
non-zero elements only in a triangular part:

$\displaystyle v^{i}_{a}\left\\{\begin{array}[]{l}\neq 0,\ \hbox{only if
}1\leq a\leq i,\\\ =0,\ \hbox{for }i<a\leq R.\end{array}\right.$ (17)

Unless $n$ and $R$ are very large, the minimum (13) can be obtained by
repeating the minimization process many times enough, starting from random
values of $\psi$ and $v$.

###  3.1 $n=1$

Though $C=C_{111}$ cannot even be called a tensor for $n=1$, this simplest
case still shows the splitting (4), as we will see below. For $n=1$, the
“tensor rank decomposition” (3) is just a scalar equation,

$\displaystyle C=\sum_{i=1}^{R}(\psi^{i})^{3},$ (18)

and, for $R>1$, it has a continuously infinite number of solutions extending
to infinity. Considering the abundance of the solutions, the dominance shown
below of the particular splitted configurations (3) in the system (1) is
rather surprising. This dominance should be a non-trivial effect of the
integration volume.

Figure 1: Left two panels: $n=1,\ R=1,\ C=1$ are taken. The left of them plots
the free energy (13) against $\bar{\beta}$. There is a first-order phase
transition at $\bar{\beta}\sim 0.8$. The right of them plots the value of
$\psi$, which has a jump from zero to a finite value. The right two panels:
Similar plots for $n=1,\ R=2,\ C=1$. Only one of $\psi^{1},\psi^{2}$ gets
finite at $\bar{\beta}>\bar{\beta}_{c}$.

Two examples with $R=1,2$ each are shown in Figure 1. In both cases, there are
first-order phase transitions at $\bar{\beta}=\bar{\beta}_{c}\sim 0.8$. In the
left example with $R=1$, the only $\psi^{1}$ makes a jump from zero to a
finite value, while, in the right example with $R=2$, only one of the
$\psi^{i}$, say $\psi^{1}$, makes a jump, but the other, $\psi^{2}$, stays
zero: The splitting (3) is realized in the phase at
$\bar{\beta}>\bar{\beta}_{c}$. Note that $\psi^{2}$, corresponding to the
minor part in (4), is exactly zero in this large-$N$ limit. One can also check
that the above phenomenon is universal for any $R\geq 2$: Only one of
$\psi^{i}$ gets non-zero values and the others stay zero at
$\bar{\beta}>\bar{\beta}_{c}$.

The above transition can be restated as a real replica symmetry breaking. In
the phase at $\bar{\beta}<\bar{\beta}_{c}$, the real replica symmetry, namely,
the symmetric group $S_{R}$ interchanging $\phi^{i}\ (i=1,2,\ldots,R)$, is
unbroken because of $\forall\psi^{i}=0$. In the other phase at
$\bar{\beta}>\bar{\beta}_{c}$, it is spontaneously broken,

$\displaystyle S_{R}\rightarrow S_{R-1},\ \hbox{at
}\bar{\beta}=\bar{\beta}_{c},$ (19)

because one of $\psi^{i}$ gets finite.

The same symmetric structure can also be checked for $\tilde{Q}$.

Figure 2: An example of $n=2,\ R=4$. $C_{111}=C_{222}=1,\ C_{122}=-2$ is taken
for $C$. From the left to the right panels, the free energy (13),
$\bar{\beta}d\bar{\beta}F/d\bar{\beta}$, and
$|\psi^{i}|=\sqrt{\psi^{i}_{a}\psi_{a}^{i}}$ are plotted against
$\bar{\beta}$, respectively.

###  3.2 $n=2$

The typical ranks101010The space of $C$ (with a normalization) can be
classified by the ranks of $C$. Typical ranks are the ranks which appear with
a finite measure in the space. In other words, when a $C$ is randomly chosen,
the rank of $C$ will be one of the typical ranks, each appearing with a non-
zero probability. of the $n=2$ real symmetric tensors are 2 and 3 [27].
Therefore, depending on the choice of $C$, we typically find 2 or 3 first-
order phase transition points for $R\geq 3$, as $\beta$ is changed. Figure 2
shows an example with $\hbox{rank}(C)=3$ and $R=4$. We indeed find three
first-order phase transitions. The number of non-zero vectors of $\psi^{i}$
increases one by one from zero to three, as the system undergoes each first-
order phase transition when $\bar{\beta}$ is increased. One vector keeps
vanishing how large $\bar{\beta}$ is taken, because $R-\hbox{rank}(C)=1$ in
this case. One can check that the splitting, as in (3), is universal for any
$R\geq 4$ when $\bar{\beta}$ is taken large enough. In this phase, the number
$\hbox{rank}(C)$ of $\psi^{i}$ takes non-zero values, while the others are
exactly zero. Thus there are typically the following two possible patterns of
real replica symmetry breaking for $n=2$:

$\displaystyle\begin{split}&S_{R}\rightarrow S_{R-1}\rightarrow S_{R-2},\hbox{
for rank}(C)=2,\\\ &S_{R}\rightarrow S_{R-1}\rightarrow S_{R-2}\rightarrow
S_{R-3},\hbox{ for rank}(C)=3.\end{split}$ (20)

Figure 3:  The error of the tensor rank decomposition, $|C-\psi\psi\psi|/|C|$,
is plotted against $\bar{\beta}$ for the same example as in Figure 2.

It would be instructive to see how the tensor rank decomposition of $C$ is
performed as $\bar{\beta}$ is increased. As shown in Figure 3, the tensor rank
decomposition is improved at each time the first-order phase transitions
occurs.

###  3.3 Lie-group symmetric $C$

Figure 4: An example with $n=3,\ R=5$ and an $SO(2)$ invariant $C$ with
$\alpha=0.5$ in (21). In the right panel, all the non-zero $\psi^{i}$ have the
same size, and are degenerate in the plot.

For $n\geq 3$, $C$ can be taken Lie-group symmetric. As an illustrative
example, let us consider the simplest case of $n=3$ and $SO(2)$-invariant $C$.
More precisely, $C$ is given by

$\displaystyle
C_{abc}=\hbox{const.}\,e^{-\alpha(m_{a}^{2}+m_{b}^{2}+m_{c}^{2})}\int_{0}^{2\pi}d\theta\,f_{a}f_{b}f_{c},$
(21)

where

$\displaystyle f_{1}=\frac{1}{\sqrt{2}},\ f_{2}=\cos(\theta),\
f_{3}=\sin(\theta),$ (22)

$m_{a}$ are the angular frequencies of the functions $f_{a}$, namely,
$m_{1}=0,m_{2}=m_{3}=1$. The parameter $\alpha$ has been introduced to make
smooth the sharp frequency cutoff, at $m=1$ in this case, by choosing
$\alpha\sim O(1)$. The overall factor const. is a normalization factor for
$C_{abc}C_{abc}=1$. It is easy to check that $C$ is invariant under an
arbitrary $SO(2)$ rotation between $f_{2}$ and $f_{3}$, which corresponds to a
shift of $\theta$.

Figure 4 shows the free energy and so on for this case. There exist four
first-order phase transitions, as $\bar{\beta}$ is increased. At each time a
transition undergoes, the number of non-zero $\psi^{i}$ increases by one, and
the maximum number is four, which is the rank of $C$.

The symmetry breaking pattern of this case is more interesting than the
previous cases. Because of the $SO(2)$ symmetry of $C$, the system initially
has $SO(2)\times S_{R}$ symmetry. Since the non-zero $\psi^{i}$ in each phase
are all different from each other, one would suspect that the fate of the real
replica symmetry breaking would be the same as the previous cases. However, it
is easy to check in each phase that the set of non-zero $\psi^{i}$ are
invariant under a discrete subgroup of $SO(2)$. For example, in the phase with
two non-zero $\psi^{i}$s, say, $\psi^{1},\psi^{2}$, they are related by
$\psi^{1}=\hbox{Rot}(\pi)\psi^{2}$, where $\hbox{Rot}(\pi)$ is the $SO(2)$
rotation by angle $\pi$ (See Figure 5). This is similar in the other phases,
with three non-zero $\psi^{i}$’s being related by $\hbox{Rot}(\pi/3)$ and four
non-zero $\psi^{i}$’s related by $\hbox{Rot}(\pi/4)$. In all, the symmetry
breaking pattern is given by

$\displaystyle SO(2)\times S_{R}\rightarrow SO(2)\times S_{R-1}\rightarrow
Z_{2}\times S_{R-2}\rightarrow Z_{3}\times S_{R-3}\rightarrow Z_{4}\times
S_{R-4},$ (23)

where $Z_{n}$ denotes the cyclic group. As illustrated in Figure 5, one can
check that non-zero $\psi^{i}$’s form discretized $S^{1}$, and the lattice
spacing becomes finer, as $\bar{\beta}$ is increased. Note that, since there
is a finite jump of configurations at each transition, any direct connections
(subgroup structre, etc.) do not generally exist between the consecutive
breaking patterns.

Figure 5: Illustration of the development of the symmetry patterns in (23).
Non-zero $\psi^{i}$’s form discretized $S^{1}$, the lattice spacing of which
becomes finer as $\bar{\beta}$ is increased. The dots illustrate the locations
of non-zero $\psi^{i}$’s, and the dashed line an imaginary $S^{1}$. With
increasing $\bar{\beta}$, an $S^{1}$ gradually emerges by the discrete steps
of the first-order phase transitions.

##  4 Random $C$

In this section we will consider the cases with random values of $C$. The
motivation is to understand the dynamics of the system (1) for general values
of $C$, rather than for some particular values as studied in Section 3.
Following the successes in the study of spin glasses, the main strategy we
take is to apply the replica trick,

$\displaystyle\bar{\beta}F=-\left.\frac{\partial}{\partial
T}\lim_{N\rightarrow\infty}\frac{1}{NR}\log\left\langle
Z_{N,R}(\beta,C)^{T}\right\rangle_{C}\right|_{T=0},$ (24)

where $\bar{\beta}=\beta/N$, $\langle\cdot\rangle_{C}$ denotes the average
over the random distribution of $C$, and $T$ is the replica number. As in (6),
we restrict the range of non-zero values of $C$ as $C_{abc}\ (a,b,c\leq n)$,
and assume them to be distributed by the normal distribution:

$\displaystyle\langle{\cal O}\rangle_{C}=A\int\prod_{a\leq b\leq
c=1}^{n}dC_{abc}\,{\cal O}\,e^{-\alpha C_{abc}C_{abc}},$ (25)

where $\alpha$ is a positive number, and $A$ is a normalization factor for
$\langle 1\rangle_{C}=1$.

Let us first rewrite the power $T$ in (24) by introducing $T$ replicas of
$\phi_{a}^{i}$:

$\displaystyle\left\langle Z_{N,R}(\beta,C)^{T}\right\rangle_{C}=A\int
dC\,e^{-\alpha C^{2}}\,Z_{N,R}(\beta,C)^{T}=A\int
dC\int_{I^{T}}\prod_{t=0}^{T-1}\prod_{a,i=1}^{N,R}d\phi^{it}_{a}\,e^{-S_{T}},$
(26)

where

$\displaystyle\begin{split}S_{T}=\alpha\,C_{abc}C_{abc}+\beta\sum_{t=0}^{T-1}\left(C_{abc}-\sum_{i=1}^{R}\phi_{a}^{it}\phi_{b}^{it}\phi_{c}^{it}\right)^{2},\end{split}$
(27)

and we have introduced an additional upper index $t\ (t=0,1,\ldots,T-1)$ for
$\phi$, which starts from zero for later convenience. With the same spirit as
(7), we divide $\phi$ into two parts,

$\displaystyle\begin{split}&\psi_{a}^{it}=\phi_{a}^{it}\ (a=1,2,\ldots,n),\\\
&\tilde{\psi}_{a}^{it}=\phi_{a}^{it}\ (a=n+1,n+2,\ldots,N).\end{split}$ (28)

After a straightforward computation, we obtain

$\displaystyle\begin{split}S_{T}&=(\alpha+\beta
T)\left(C_{abc}-\frac{\beta}{\alpha+\beta
T}\sum_{i,t}\psi_{a}^{it}\psi_{b}^{it}\psi_{b}^{it}\right)^{2}-\frac{\beta^{2}}{\alpha+\beta
T}\sum_{i,i^{\prime},t,t^{\prime}}\left(\psi^{it}\cdot\psi^{i^{\prime}t^{\prime}}\right)^{3}\\\
&+\beta\sum_{i,i^{\prime},t}\left(\psi^{it}\cdot\psi^{i^{\prime}t}\right)^{3}+3\beta\sum_{i,i^{\prime},t}\left(\psi^{it}\cdot\psi^{i^{\prime}t}\right)^{2}\left(\tilde{\psi}^{it}\cdot\tilde{\psi}^{i^{\prime}t}\right)\\\
&+3\beta\sum_{i,i^{\prime},t}\left(\psi^{it}\cdot\psi^{i^{\prime}t}\right)\left(\tilde{\psi}^{it}\cdot\tilde{\psi}^{i^{\prime}t}\right)^{2}+\beta\sum_{i,i^{\prime},t}\left(\tilde{\psi}^{it}\cdot\tilde{\psi}^{i^{\prime}t}\right)^{3},\end{split}$
(29)

where the ranges of sums have been omitted for brevity, since they are
obvious.

The first term of (29) can be integrated over $C$, and this generates

$\displaystyle\frac{\\#C}{2}\log\left(\frac{\alpha+\beta T}{\alpha}\right)$
(30)

as an additional term to the exponent, where $\\#C=n(n+1)(n+2)/6$, namely, the
number of independent elements of $C$. To compute the other terms, let us
introduce the overlaps,

$\displaystyle\begin{split}&Q_{iti^{\prime}t^{\prime}}=\psi^{it}\cdot\psi^{i^{\prime}t^{\prime}},\\\
&\tilde{Q}_{iti^{\prime}t^{\prime}}=\tilde{\psi}^{it}\cdot\tilde{\psi}^{i^{\prime}t^{\prime}},\end{split}$
(31)

as in (9). The same procedure as before generates similar logarithmic terms as
in (12),

$\displaystyle-\frac{n}{2}\log\det Q-\frac{N-n}{2}\log\det\tilde{Q},$ (32)

where $Q$ and $\tilde{Q}$ are regarded as $RT\times RT$ matrices in the
determinants. By assembling the above expressions, we obtain

$\displaystyle\begin{split}S_{T}^{\rm
eff}(\beta)=&\frac{\\#C}{2}\log\left(1+\frac{\beta}{\alpha}T\right)-\frac{\beta^{2}}{\alpha+\beta
T}\sum_{i,i^{\prime},t,t^{\prime}}\left(Q^{iti^{\prime}t^{\prime}}\right)^{3}+\beta\sum_{i,i^{\prime},t}\left(Q^{iti^{\prime}t}\right)^{3}\\\
&+3\beta\sum_{i,i^{\prime},t}\left(Q^{iti^{\prime}t}\right)^{2}\left(\tilde{Q}^{iti^{\prime}t}\right)+3\beta\sum_{i,i^{\prime},t}\left(Q^{iti^{\prime}t}\right)\left(\tilde{Q}^{iti^{\prime}t}\right)^{2}+\beta\sum_{i,i^{\prime},t}\left(\tilde{Q}^{iti^{\prime}t}\right)^{3}\\\
&-\frac{n}{2}\log\det Q-\frac{N-n}{2}\log\det\tilde{Q}.\end{split}$ (33)

For later use, it is more convenient to perform a rescaling of variables in
(33): $Q=(N/\beta)^{1/3}\bar{Q},\ \tilde{Q}=(N/\beta)^{1/3}\bar{\tilde{Q}}$.
Then we obtain

$\displaystyle\begin{split}S_{T}^{\rm
eff}(\beta)=&\frac{\\#C}{2}\log\left(1+\frac{\beta}{\alpha}T\right)+\frac{NRT}{6}\log\beta-\frac{N\beta}{\alpha+\beta
T}\sum_{i,i^{\prime},t,t^{\prime}}\left(\bar{Q}^{iti^{\prime}t^{\prime}}\right)^{3}+N\sum_{i,i^{\prime},t}\left(\bar{Q}^{iti^{\prime}t}\right)^{3}\\\
&+3N\sum_{i,i^{\prime},t}\left(\bar{Q}^{iti^{\prime}t}\right)^{2}\left(\bar{\tilde{Q}}^{iti^{\prime}t}\right)+3N\sum_{i,i^{\prime},t}\left(\bar{Q}^{iti^{\prime}t}\right)\left(\bar{\tilde{Q}}^{iti^{\prime}t}\right)^{2}+N\sum_{i,i^{\prime},t}\left(\bar{\tilde{Q}}^{iti^{\prime}t}\right)^{3}\\\
&-\frac{n}{2}\log\det\bar{Q}-\frac{N-n}{2}\log\det\bar{\tilde{Q}},\end{split}$
(34)

where we have ignored an unimportant constant shift.

Further computations depend on assumptions made on replica symmetry breaking.
Below we consider only two possibilities, replica symmetric (RS) and one-step
replica symmetry breaking (1RSB). A reason to consider only these
possibilities is that 1RSB is known to be exact [21] in the spherical $p$-spin
model. In fact numerical simulations, which we will show in Section 5, seem to
support this assumption at least in the large-$N$ limit of this paper.

###  4.1 $R=1$ with random $C$

Since $R=1$, we can simply write $\bar{Q}^{tt^{\prime}}$ in place of
$\bar{Q}^{iti^{\prime}t^{\prime}}$, neglecting the real replica index. Under
the assumption of 1RSB [21, 22], the overlap $\bar{Q}$ is assumed to have the
form,

$\displaystyle\bar{Q}^{tt^{\prime}}=\delta_{\lfloor t/M\rfloor,\lfloor
t^{\prime}/M\rfloor}I_{t\bmod M,t^{\prime}\bmod
M}(q_{0},q_{1})+q_{2}(1-\delta_{\lfloor t/M\rfloor,\lfloor
t^{\prime}/M\rfloor}),$ (35)

where $q_{i}$ are new variables, $\lfloor\cdot\rfloor$ is the floor function,
$a\bmod b$ denotes $a$ modulo $b$, and $I$ is an $M\times M$ matrix with
components,

$\displaystyle
I_{tt^{\prime}}(q_{0},q_{1})=q_{0}\,\delta_{tt^{\prime}}+q_{1}(1-\delta_{tt^{\prime}}).$
(36)

Then it is straightforward to obtain

$\displaystyle\begin{split}&\frac{1}{T}\sum_{t,t^{\prime}}(\bar{Q}^{tt^{\prime}})^{3}=(T-M)q_{2}^{3}+(M-1)q_{1}^{3}+q_{0}^{3},\\\
&\frac{1}{T}\sum_{t}(\bar{Q}^{tt})^{3}=q_{0}^{3},\\\
&\log\det\bar{Q}=\log\left(q_{0}+(M-1)q_{1}+(T-M)q_{2}\right)+\left(\frac{T}{M}-1\right)\log\left(q_{0}-q_{1}+M(q_{1}-q_{2})\right)\\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+\left(T-\frac{T}{M}\right)\log(q_{0}-q_{1}).\end{split}$ (37)

These expressions are also assumed for $\bar{\tilde{Q}}$, with replacements
$q_{i}\rightarrow\tilde{q}_{i},\ M\rightarrow\tilde{M}$.

By putting (37) and the corresponding expressions of $\bar{\tilde{Q}}$ to
(34), and applying the replica trick formula in (24), one obtains

$\displaystyle\begin{split}\bar{\beta}&F^{1RSB}=\frac{\\#C}{2N}\frac{\beta}{\alpha}+\frac{1}{6}\log\beta-\frac{\beta}{\alpha}\left(-Mq_{2}^{3}+(M-1)q_{1}^{3}+q_{0}^{3}\right)+q_{0}^{3}+3q_{0}^{2}\tilde{q}_{0}+3q_{0}\tilde{q}_{0}^{2}+\tilde{q}_{0}^{3}\\\
&-\frac{\bar{n}}{2}\left(\left(1-\frac{1}{M}\right)\log(q_{0}-q_{1})+\frac{1}{M}\log(q_{0}-q_{1}+M(q_{1}-q_{2}))+\frac{q_{2}}{q_{0}+(M-1)q_{1}-Mq_{2}}\right)\\\
&-\frac{1-\bar{n}}{2}\left(\left(1-\frac{1}{\tilde{M}}\right)\log(\tilde{q}_{0}-\tilde{q}_{1})+\frac{1}{\tilde{M}}\log(\tilde{q}_{0}-\tilde{q}_{1}+\tilde{M}(\tilde{q}_{1}-\tilde{q}_{2}))+\frac{\tilde{q}_{2}}{\tilde{q}_{0}+(\tilde{M}-1)\tilde{q}_{1}-\tilde{M}\tilde{q}_{2}}\right),\end{split}$
(38)

where we have introduced $\bar{n}=n/N$.

Here let us comment on the ranges of the parameters which appear in the above
expressions (See [22] for more details). From the physical point of view, the
overlaps of configurations in the same state will not be less than those in
different states. Therefore

$\displaystyle q_{0}\geq q_{1}\geq q_{2},\
\tilde{q}_{0}\geq\tilde{q}_{1}\geq\tilde{q}_{2}.$ (39)

In addition, the computations of
$\langle\langle\psi_{a_{1}}\cdots\psi_{a_{k}}\rangle\langle\psi_{a_{1}}\cdots\psi_{a_{k}}\rangle\rangle_{C}$
under the 1RSB assumption (and the same for $\tilde{\psi}$) leads to a
probability distribution of the overlap,

$\displaystyle\langle
P(q)\rangle_{C}=\frac{M-1}{T-1}\delta(q-q_{1})+\frac{T-M}{T-1}\delta(q-q_{2}),$
(40)

and a similar one for $\tilde{q}_{i}$. Then the positivity of this probability
distribution in the $T\rightarrow 0$ limit requires

$\displaystyle 0\leq M\leq 1,\ 0\leq\tilde{M}\leq 1.$ (41)

One can check that (39) and (41) assures the arguments of the logarithms and
the denominators in (38) do not become less than zero.

For simplicity of discussions, in the rest of this paper, we assume
$q_{2},\tilde{q}_{2}=0$, since this is known to be true for the $p$-spin
spherical model [21]. The consistency of our analysis of this section and the
numerical results of Section 5 will support this simplification.

When $\beta$ is small enough, namely, in high temperatures, the system will be
in RS phase. The RS expression of the free energy can be obtained by putting
$q_{1}=\tilde{q}_{1}=0$ (and $q_{2}=\tilde{q}_{2}=0$) to the 1RSB expression
(38):

$\displaystyle\begin{split}\bar{\beta}F^{RS}=&\frac{\\#C}{2N}\frac{\beta}{\alpha}+\frac{1}{6}\log\beta+\left(1-\frac{\beta}{\alpha}\right)q_{0}^{3}+3q_{0}^{2}\tilde{q}_{0}+3q_{0}\tilde{q}_{0}^{2}+\tilde{q}_{0}^{3}-\frac{\bar{n}}{2}\log
q_{0}-\frac{1-\bar{n}}{2}\log\tilde{q}_{0}.\end{split}$ (42)

Note that $M,\tilde{M}$ have disappeared from the expression, since they have
no roles for $q_{1}=\tilde{q}_{1}=q_{2}=\tilde{q}_{2}=0$ in (35).

The free energy of the system is obtained by searching for the stationary
points of these expressions (38) and (42) in terms of the variables,
$q_{i},\tilde{q}_{i},M,\tilde{M}$.

####  4.1.1 $R=1$, $\bar{n}=1$ with random $C$

In this subsection, let us perform concrete analysis on the case with $R=1,\
\bar{n}=1$ and random $C$ by using the expressions above. Because of
$\bar{n}=1$, we can simply ignore the terms with $\tilde{q}_{i}$. Therefore,
from (38) and (42), we have

$\displaystyle\bar{\beta}F^{RS}_{R=1\,\bar{n}=1}=$
$\displaystyle\frac{\\#C}{2N}\frac{\beta}{\alpha}+\frac{1}{6}\log\beta+\left(1-\frac{\beta}{\alpha}\right)q_{0}^{3}-\frac{1}{2}\log
q_{0},$ (43)

and

$\displaystyle\begin{split}\bar{\beta}F^{1RSB}_{R=1\,\bar{n}=1}=&\frac{\\#C}{2N}\frac{\beta}{\alpha}+\frac{1}{6}\log\beta-\frac{\beta}{\alpha}\left((M-1)q_{1}^{3}+q_{0}^{3}\right)+q_{0}^{3}\\\
&-\frac{1}{2M}\left((M-1)\log(q_{0}-q_{1})+\log(q_{0}-q_{1}+Mq_{1})\right),\end{split}$
(44)

where we have put $q_{2}=0$ following the spherical $p$-spin model [21].

Let us first consider the RS case (43). By solving
$\frac{\partial\bar{\beta}F_{R=1\,\bar{n}=1}^{RS}}{\partial q_{0}}=0$, one
obtains

$\displaystyle\bar{\beta}F^{RS}_{R=1\,\bar{n}=1}=\frac{\\#C}{2N}\frac{\beta}{\alpha}+\frac{1}{6}\log\beta+\frac{1}{6}+\frac{1}{6}\log
6+\frac{1}{6}\log\left(1-\frac{\beta}{\alpha}\right).$ (45)

Because of the divergence at $\beta=\alpha$, the RS phase cannot exist for all
$\beta$ and there must be a critical point $\beta_{c}<\alpha$, by which the RS
phase is bounded.

A comment is in order. The appearance of this diverging point can be traced to
the appearance of a curious indefinite coefficient $1-\beta/\alpha$ in (43).
This looks confusing, because the starting expression (27) is positive
semidefinite. However, after integrating out $C$, we obtain (33), which
contains a negative coefficient. One can check that this expression is still
positive semidefinite for $T\geq 1$, but after taking $T=0$ of the replica
trick, it becomes indefinite. Therefore the presence of a phase transition is
an interplay between integrating out and the analytic continuation.

We assume that the RS solution is connected to the 1RSB solution at
$\beta=\beta_{c}<\alpha$ with the requirement that the free energy is
continuous at this point. The stationary condition of (44) is given by

$\displaystyle\frac{\partial\bar{\beta}F^{1RSB}_{R=1\,\bar{n}=1}}{\partial
q_{0}}=\frac{\partial\bar{\beta}F^{1RSB}_{R=1\,\bar{n}=1}}{\partial
q_{1}}=\frac{\partial\bar{\beta}F^{1RSB}_{R=1\,\bar{n}=1}}{\partial M}=0.$
(46)

These equations, however, cannot be fully solved algebraically, because of the
logarithmic terms in the last equation. However, using Mathematica, the first
two equations can analytically be solved, in which $q_{0},q_{1}$ are
analytically expressed by $M$. Then the last equation can be numerically
solved to obtain $M$. With this procedure, one can compute the free energy for
each value of $\alpha,\beta$.

Figure 6: The free energy and its first and second derivatives for $R=1,\
\bar{n}=1$ and random $C$. $\alpha=1/2$ is taken without loss of generality.
$F^{\prime}$ denotes the free energy with the subtraction of the first term of
(43) and (44). There is a second order phase transition with a finite jump of
the second derivative at $\beta_{c}\sim 0.407\ (\log_{10}\beta_{c}\sim-0.39)$.

Figure 6 shows the free energy up to the second derivatives. Here, we have
subtracted the first term in (43) and (44) and have denoted the subtracted
free energy by $F^{\prime}$, because this first term behaves in $\sim N^{2}$
and is much larger than the other terms. Note that this term is linear in
$\beta$ and therefore does not appear in the second derivative. To be
concrete, we take $\alpha=1/2$ without loss of generality. The transition
point can be determined by requiring the continuity of the free energy,
$F^{RS}_{R=1\,\bar{n}=1}=F^{1RSB}_{R=1\,\bar{n}=1}$, at $\beta=\beta_{c}$, and
the result is $\beta_{c}\sim 0.407\ (\log_{10}0.407\sim-0.39)$. The system is
in the RS phase for $\beta<\beta_{c}$, and in the 1RSB phase for
$\beta>\beta_{c}$. The transition is second-order with a finite jump of the
second derivative of the free energy with respect to $\beta$. The reason for
the continuity of the first derivative of the free energy at the transition
point can be attributed to the continuity of $q_{0}$ and $M=1$, while $q_{1}$
makes a jump, as shown in Figure 7.

Figure 7: The behavior of the configurations for the same case as in Figure 6.
$q_{0}$ is continuous, while $q_{1}$ makes a jump at the transition point.

####  4.1.2 $R=1$, $\bar{n}=1/2$ with random $C$

In this subsection, we consider a case with $R=1$, $\bar{n}<1$, and random
$C$. To be concrete, we set $\bar{n}=1/2$ as a representative case of
$\bar{n}<1$. The main difference from the $\bar{n}=1$ case is that we find a
first-order phase transition in place of the second-order one for $\bar{n}=1$.

Let us first reduce the general form of the free energy given in (38), since
it seems enough to assume $q_{2},\tilde{q}_{1},\tilde{q}_{2}=0$ for the
present analysis. Here the assumption $\tilde{q}_{1},\tilde{q}_{2}=0$ of RS
for the variable $\bar{\tilde{Q}}$ would be natural, because only its diagonal
components appear in the interaction terms of (34). By putting these
assumptions and $\bar{n}=1/2$ to (38), we obtain

$\displaystyle\begin{split}\bar{\beta}&F^{1RSB}_{R=1\,\bar{n}=1/2}=\frac{\\#C}{2N}\frac{\beta}{\alpha}+\frac{1}{6}\log\beta-\frac{\beta}{\alpha}\left((M-1)q_{1}^{3}+q_{0}^{3}\right)+q_{0}^{3}+3q_{0}^{2}\tilde{q}_{0}+3q_{0}\tilde{q}_{0}^{2}+\tilde{q}_{0}^{3}\\\
&-\frac{1}{4M}\left(\left(M-1\right)\log(q_{0}-q_{1})+\log(q_{0}-q_{1}+Mq_{1})\right)-\frac{1}{4}\log\tilde{q}_{0}.\end{split}$
(47)

The RS case can simply be obtained by putting $q_{1}=0$ to (47):

$\displaystyle\begin{split}\bar{\beta}&F^{RS}_{R=1\,\bar{n}=1/2}=\frac{\\#C}{2N}\frac{\beta}{\alpha}+\frac{1}{6}\log\beta+\left(1-\frac{\beta}{\alpha}\right)q_{0}^{3}+3q_{0}^{2}\tilde{q}_{0}+3q_{0}\tilde{q}_{0}^{2}+\tilde{q}_{0}^{3}-\frac{1}{4}\log
q_{0}-\frac{1}{4}\log\tilde{q}_{0}.\end{split}$ (48)

In the RS case, the equations,
$\partial\bar{\beta}F^{RS}_{R=1\,\bar{n}=1/2}/\partial
q_{0}=\partial\bar{\beta}F^{RS}_{R=1\,\bar{n}=1/2}/\partial\tilde{q}_{0}=0$,
can be solved algebraically for $q_{0},\tilde{q}_{0}$, and the free energy can
be obtained as a function of $\beta/\alpha$, though the expression is too
complicated to explicitly write down here.

As for the 1RSB case, we have to fully resort to numerical analysis. To be
concrete, let us take $\alpha=1/2$ without loss of generality. Then we can
find that there is a phase transition between RS and 1RSB phases at
$\beta=\beta_{c}\sim 0.58\ (\log_{10}\beta_{c}\sim-0.24)$, as is shown in
Figure 8. In the figure, the first terms of (47) and (48) are subtracted in
showing the free energy and its first derivative. The transition is first-
order, as can be clearly seen in the plots (see also Figure 9).

Figure 8: The free energy and its first and second derivatives for $R=1,\
\bar{n}=1/2$. Without loss of generality, $\alpha=1/2$ is taken to be
concrete. A first-order phase transition exists at $\beta=\beta_{c}\sim 0.58\
(\log_{10}\beta_{c}\sim-0.24)$. Here $F^{\prime}$ denotes the free energy with
the subtraction of the first term.

Figure 9: The configurations are plotted for the same case as in Figure 8.
$q_{0}$ and $q_{1}$ are shown by solid and dashed lines, respectively, in the
left panel. The middle and right panels show $\tilde{q}_{0}$ and $M$,
respectively. The configurations jump at the phase transition point.

###  4.2 $R>1,\ \bar{n}=1$ with random $C$

To analyze the $R>1,\ \bar{n}=1$ case we assume the following form for
$\bar{Q}$:

$\displaystyle\begin{split}\bar{Q}^{iti^{\prime}t^{\prime}}&=\delta_{\lfloor
t/M\rfloor,\lfloor t^{\prime}/M\rfloor}\left(\delta_{ii^{\prime}}I_{t\bmod
M,t^{\prime}\bmod M}(q_{0},q_{1})+(1-\delta_{ii^{\prime}})I_{t\bmod
M,t^{\prime}\bmod M}(q_{0}^{\prime},q_{1}^{\prime})\right),\end{split}$ (49)

were the matrix $I$ is defined in (36), and we have already assumed that the
$q_{2},q_{2}^{\prime}$ terms are not needed as in the previous cases. Here
note that we are assuming that the real replica symmetry concerning the
$i,i^{\prime}$ indices of $\bar{Q}^{ii^{\prime}tt^{\prime}}$ is not broken,
while the replica symmetry can be broken in general concerning $t,t^{\prime}$.

Putting some formula shown in Appendix B to (34) and using (24), we obtain the
free energy as

$\displaystyle\begin{split}\bar{\beta}F_{R>1,\bar{n}=1}^{1RSB}&=\frac{\\#C}{2NR}\frac{\beta}{\alpha}+\frac{1}{6}\log\beta-\frac{\beta}{\alpha}\left((M-1)q_{1}^{3}+q_{0}^{3}\right)+q_{0}^{3}\\\
&+(R-1)\left(-\frac{\beta}{\alpha}\left((M-1)q^{\prime}_{1}{}^{3}+q^{\prime}_{0}{}^{3}\right)+q^{\prime}_{0}{}^{3}\right)\\\
&-\frac{1}{2MR}\Big{[}\log\left\\{q_{0}+(R-1)q_{0}^{\prime}+(M-1)(q_{1}+(R-1)q_{1}^{\prime})\right\\}\\\
&\ \ \ \ \ \ \ \ \ \
+(M-1)\log\left\\{q_{0}+(R-1)q_{0}^{\prime}-(q_{1}+(R-1)q_{1}^{\prime})\right\\}\\\
&\ \ \ \ \ \ \ \ \ \
+(R-1)\log\left\\{q_{0}-q_{0}^{\prime}+(M-1)(q_{1}-q_{1}^{\prime})\right\\}\\\
&\ \ \ \ \ \ \ \ \ \
+(M-1)(R-1)\log\left\\{q_{0}-q_{0}^{\prime}-(q_{1}-q_{1}^{\prime})\right\\}\Big{]}.\end{split}$
(50)

The RS case can be obtained by putting $q_{1}=q_{1}^{\prime}=0$ as

$\displaystyle\begin{split}\bar{\beta}F_{R>1,\bar{n}=1}^{RS}&=\frac{\\#C}{2NR}\frac{\beta}{\alpha}+\frac{1}{6}\log\beta+\left(1-\frac{\beta}{\alpha}\right)q_{0}^{3}+(R-1)\left(1-\frac{\beta}{\alpha}\right)q^{\prime}_{0}{}^{3}\\\
&-\frac{1}{2R}\Big{[}\log\left(q_{0}+(R-1)q_{0}^{\prime}\right)+(R-1)\log\left(q_{0}-q_{0}^{\prime}\right)\Big{]}.\end{split}$
(51)

In a similar manner as in the previous sections, one can numerically analyze
the stationary solutions of (50) and (51). Then what we find is that there are
only solutions with $q_{0}^{\prime}=q_{1}^{\prime}=0$. Therefore the present
case reduces to the case with $R=1,\ \bar{n}=1$ case in Section 4.1.1, as (50)
and (51) reduce to (44) and (43), respectively.

###  4.3 $R>1,\ \bar{n}<1$

In this case, we have $\bar{\tilde{Q}}$ in addition. Since only the diagonal
components of $\bar{\tilde{Q}}$ with respect to $tt^{\prime}$ indices couple
with $\bar{Q}$ in (34), it would be reasonable to assume the following
diagonal form for $\bar{\tilde{Q}}$:

$\displaystyle\bar{\tilde{Q}}^{iti^{\prime}t^{\prime}}=\delta_{tt^{\prime}}\left(q_{0}\delta_{ii^{\prime}}+(1-\delta_{ii^{\prime}})q_{0}^{\prime}\right).$
(52)

Then the interaction terms are given by

$\displaystyle\begin{split}&3\sum_{i,i^{\prime},t}\left(\bar{Q}^{iti^{\prime}t}\right)^{2}\left(\bar{\tilde{Q}}^{iti^{\prime}t}\right)+3\sum_{i,i^{\prime},t}\left(\bar{Q}^{iti^{\prime}t}\right)\left(\bar{\tilde{Q}}^{iti^{\prime}t}\right)^{2}+\sum_{i,i^{\prime},t}\left(\bar{\tilde{Q}}^{iti^{\prime}t}\right)^{3}\\\
&=RT\left(3\left(q_{0}^{2}\tilde{q}_{0}+(R-1){q^{\prime}_{0}}^{2}\tilde{q}^{\prime}_{0}\right)+3\left(q_{0}\tilde{q}_{0}^{2}+(R-1)q^{\prime}_{0}\tilde{q}^{\prime}_{0}{}^{2}\right)+\left(\tilde{q}_{0}^{3}+(R-1)\tilde{q}^{\prime}_{0}{}^{3}\right)\right),\end{split}$
(53)

and the determinant term of $\bar{\tilde{Q}}$ is given by

$\displaystyle\log\det\bar{\tilde{Q}}=T\log(\tilde{q}_{0}+(R-1)\tilde{q}_{0}^{\prime})+(R-1)T\log(\tilde{q}_{0}-\tilde{q}_{0}^{\prime}).$
(54)

By appropriately including the parameter $\bar{n}$, the expression of the free
energy can straightforwardly be obtained from (50), (53) and (54) in the same
manner as before. Then we numerically searched for the stationary solutions.
However, all the solutions have
$q_{0}^{\prime}=q_{1}^{\prime}=\tilde{q}_{0}^{\prime}=0$. This concludes that
the present case also reduces to the $R=1$ case.

##  5 Comparisons with numerical simulations

In this section, we perform some numerical simulations, and compare with the
exact results. By this we can check the exact results obtained based on some
assumptions in the previous sections, and can also check the large-$N$
convergence explicitly. The latter will be useful for future study, since this
provides us some information on how much we can rely on numerical simulations
in the analysis of our system.

We apply Hamiltonian Monte Carlo method [28] to the system (1), accompanied
with parallel tempering [29] across different values of $\beta$. We compute
$\beta\langle s\rangle$ and $\beta^{2}\langle
s^{2}\rangle_{c}=\beta^{2}(\langle s^{2}\rangle-\langle s\rangle^{2})$ with
$s=\beta(Q-\phi\phi\phi)^{2}$ to compare with $\beta\frac{\partial\beta
F}{\partial\beta},\ -\beta^{2}\frac{\partial^{2}\beta F}{\partial\beta^{2}}$
of the exact results, respectively, where $\langle\cdot\rangle$ denotes the
statistical average in the Monte Carlo simulations. We also compute the
$\beta$ dependencies of the configurations. The error estimates are performed
by the Jackknife method.

As for random $C$, every independent component of $C$ is generated by the
normal distribution with mean zero and standard deviation $1/\sqrt{m}$, where
$m$ is the multiplicity of the components, namely,
$m(C_{iii})=1,m(C_{iij})=3,m(C_{ijk})=6$ for different $i,j,k$111111This
rescaling is necessary, because the $SO(N)$ invariant expression has the form,
$C_{abc}C_{abc}=C_{iii}{}^{2}+3\,C_{iij}{}^{2}+6\,C_{ijk}{}^{2}$ due to $C$
being a symmetric tensor.. This distribution corresponds to taking
$\alpha=1/2$, which was indeed taken in the plots shown in the previous
sections. In general the details of the results depend on each particular
randomly generated $C$, but we find that, when $N$ is large enough, such
sample dependence seems irrelevant at least for the data we will show, except
for one quantity in Section 5.4. We will actually see more or less good
agreement with the exact results, even if we take one sample, except for the
quantity in Section 5.4.

The machine had a Xeon W2295 (3.0GHz, 18 cores), 128GB DDR4 memory, and Ubuntu
20 as OS. The program was written in C++ with the use of pthread for
parallelization. The leapfrog numbers were taken around a hundred. The numbers
of the samples of each run are $10^{4}$-$10^{6}$. Every run typically took for
3-20 hours121212The longest run was about 50 hours, which was for $N=70,\
R=2$. with active use of parallelization, without serious tunings of the speed
of the program.

We consider some representative cases in the following subsections.

###  5.1 $n=2,\ R=4$, fixed $C$

This case corresponds to the exact results in Section 3.2.

Figure 10 shows the comparison with the numerical simulations. The convergence
of the free energy (first derivative) to the exact result is rather slow (the
left panel), and it would be difficult to conclude the presence of multiple
first-order transitions solely from the data. On the other hand, from the plot
of the configurations (the right panel), one would better be able to conclude
there are jumps, implying first-order phase transitions.

Figure 10: The case with $n=2,\ R=4$ with fixed $C$. Comparisons between the
exact results of Section 3.2 and the numerical simulations shown by dots with
error bars. $N=20$ results are shown blue, and $N=60$ shown red.

###  5.2 $\bar{n}=1,\ R=1$, random $C$

This corresponds to the exact results in Section 4.1.1.

As for the free energy, to kill the first term in (43) and (44), which is
$O(N^{2})$ larger than the other terms, we compare the second derivative of
the free energy with the numerical results131313In fact, we tried to compare
the first derivative of the free energy with the simulation data by
subtracting the first term by hand. However, the data still seemed to be
affected by some unknown quantities smaller than $O(N^{2})$ but significant,
and we could not perform any reliable comparisons., which are shown in the
left two panels of Figure 11. We see that the numerical results seem to
converge to the exact result, though the speed of convergence is rather slow.
Especially, the errors are very large in the right region (1RSB) of the
transition point, which is supposed to be glassy.

In the right panel of Figure 11, we compare $q_{0}$ and
$(\beta/N)^{1/3}\langle\psi\cdot\psi\rangle$ (Note the change of the
normalization for (34)), and find a very good agreement between them. In this
sense, the situation seems to be the same as in Section 5.1, namely, the
numerical simulations are more reliable on configurations than on the free
energy.

Figure 11: The case with $\bar{n}=1,\ R=1$ with random $C$. The exact results
of Section 4.1.1 are shown by solid lines, and the numerical simulations by
dots with error bars. $N=30$ results are shown blue, and $N=50$ shown red.
Left two panels: The second derivative of the free energy. Right panel:
$q_{0}$ (solid line) and $(\beta/N)^{1/3}\langle\psi\cdot\psi\rangle$ (dots
with error bars).

###  5.3 $\bar{n}=1/2,\ R=1$, random $C$

This corresponds to the exact results in Section 4.1.2. The analysis is almost
the same as in Section 5.2 but with an additional data for $\tilde{q}_{0}$.
The results are plotted in Figure 12. The convergence is much slower than the
previous case. Especially, the convergence of the second derivative of the
free energy is not so good, even raising some doubts about the exact results,
but the convergence of the configurations seems to support the exact results.

Figure 12: The case with $\bar{n}=1/2,\ R=1$ and random $C$. The simulation
data are colored: $N=30$ with blue, $N=50$ with back, and $N=70$ with red.
Left: The second derivative of the free energy. Middle:
$(\beta/N)^{1/3}\langle\psi\cdot\psi\rangle$ compared with $q_{0}$ (solid)
($q_{1}$ is shown dashed). Right:
$(\beta/N)^{1/3}\langle\tilde{\psi}\cdot\tilde{\psi}\rangle$ compared with
$\tilde{q}_{0}$.

###  5.4 $R=2$, random $C$

As discussed in Sections 4.2 and 4.3, the $R>1$ cases reduce to the $R=1$
cases discussed in Sections 4.1.1 and 4.1.2. We took $R=2$ and $\bar{n}=1,1/2$
for numerical simulations, and compared with the exact results. As before,
there are good matches for $q_{0},q_{1}$, while the agreement on the free
energies is good in the RS region, but is difficult to compare because of
large errors in the 1RSB region. The reduction to $R=1$ is based on that
$q_{0}^{\prime}=q_{1}^{\prime}=\tilde{q}_{0}^{\prime}=0$ obtained in Sections
4.1.1 and 4.1.2. As shown in Figure 13, $\langle\psi^{1}\cdot\psi^{2}\rangle$
takes non-zero values in the 1RSB region and are dependent on samples. We
expect that the average of these values over samples should vanish in the
large $N$ limit to match $q_{0}^{\prime}=q_{1}^{\prime}=0$, but the number of
the data is too small to derive such a conclusion. On the other hand, as shown
in the right panel of Figure 13,
$\langle\tilde{\psi}^{1}\cdot\tilde{\psi}^{2}\rangle$ is small and would be
consistent with zero, in agreement $\tilde{q}_{0}^{\prime}=0$.

Figure 13: Left panel: $\bar{n}=1,\ R=2$ with random $C$.
$(\beta/N)^{1/3}\langle\psi^{1}\cdot\psi^{2}\rangle$ is plotted for $N=30$
(blue), $N=50$ (black) and $N=70$ (red, orange). Two data are plotted to see
the sample dependence for $N=70$. Right panel: $\bar{n}=1/2,\ R=2$ with random
$C$. $(\beta/N)^{1/3}\langle\tilde{\psi}^{1}\cdot\tilde{\psi}^{2}\rangle$ is
plotted for $N=30$ (blue), $N=50$ (black) and $N=70$ (red).

##  6 Summary and discussions

In matrix models, there are phase transitions in which distributions of
dynamical variables change topologically, like in Gross-Witten-Wadia
transition [11, 12] and in the transitions among the large-$N$ limit multi-cut
solutions [13]. In a recent study [10], similar splitting-merging behavior of
dynamical variables was observed in a tensor-vectors system by numerical
simulations, but the results were not convincing enough to characterize it and
conclude whether this is a phase transition or just a crossover. In this
paper, we have studied the system exactly in some large-$N$ limits, and have
found cascades of first-order phase transitions for fixed tensors, and a
second- or first- order phase transition for random tensors, applying the
replica trick for the random cases. These phases can be distinguished by
breaking patterns of real replica symmetries for fixed tensors, and those of
replica symmetries for random tensors, respectively. We have also performed
some numerical simulations to compare with the exact results: We have found
consistent results, which support the assumptions made in the derivation of
the exact results; We have also found rather slow convergence of numerical
data toward sharp transitions of the exact results, which implies the
necessity of our cautious attitude toward numerical simulations of our system
in future study.

As will be mentioned in the final paragraph, the large-$N$ limit of this paper
is different from that in the tensor model. However, the results of this paper
suggest the following interesting possibilities. The presence of the sharp
phase transitions implies that what we called the quantum and classical phases
in the previous study [10] of the tensor model wave function could actually be
different phases. This means that classical spacetimes could emerge through
phase transitions in the tensor model. A new finding of this paper is that the
transition rather consists of multiple first-order phase transitions.
Therefore the classical phase, where there emerge classical spacetimes, could
not be just a single phase but rather a collection of phases. This aspect can
explicitly be seen in the $SO(2)$-invariant example given in Section 3.3: We
see that the number of discrete points forming $S^{1}$ increases one by one as
the system undergoes the phase transitions, as in Figure 5. It would be
straightforward to check similar matters in more general cases $S^{n}$ with
$SO(n+1)$ symmetry [10], and it would also be more interesting to consider
general spacetimes like black holes in future study, using the general
procedure of constructing tensors corresponding to geometries developed in
[30, 10].

Here we would like to stress that the phase transitions in the previous
paragraph are all first-order and should be linked solely to the emergence of
classical spacetimes, which should not be mixed with emergence of continuum
physics in them. The latter is a different thing which could be checked by
analyzing fluctuation modes around the emergent classical spacetimes.

Another potential application of our system is to the tensor rank
decomposition [14, 15, 16, 17] in data analysis. This is an important
technique to extract information from tensors of real-life data in data
analysis, but an efficient procedure is yet to be established because of its
hardness [19]. An interesting insight from our model is that the problems
should be common across the tensor rank decomposition and the spin glasses,
since our model can be seen bilaterally as a cost function of the former and a
generalization of the spherical $p$-spin model for the latter. In this sense,
this paper can be seen as a single instance of using a technique developed in
spin glass theory for a problem in data analysis (See also [18]). Though
restricted to the large-$N$ limit with finite $R$ in this paper, the
confirmation of the splitting phenomenon (4) as a result of phase transitions
seems meaningful for the tensor rank decomposition, since the splitting gives
a tensor rank decomposition with an automatically determined rank for an
arbitrarily given tensor. Though this phenomenon cannot be used immediately
for the actual process of the tensor rank decomposition due to its
inefficiency, studying our model and the phenomenon more deeply could lead to
new useful procedures in the future.

As argued in Section 1, our system (1) with random $C$ can be considered to be
a variant of the spherical $p$-spin model with a new non-trivial parameter
$R$. In this paper’s limit of large-$N$ with finite $R$, however, the RS/1RSB
analyses have reduced to the same as $R=1$, and this seemed to be consistent
with the numerical results, as shown in Sections 4 and 5. While this suggests
the possibility that our system is just included in the same universality
class as the standard one, we have found large fluctuations of the overlaps
between real replicas in Section 5.4. This fact would suggest that, for our
system to be properly treated in mean field theory, we would need to consider
many numbers of real replicas by taking large-$R$ as well. Such large-$N,R$
limits could open new possibilities.

In fact the large-$N$ limits with finite $R$ we considered in this paper has
the disadvantage that the distribution of dynamical variables is a finite set
of configurations. If it were continuous, the dynamics would be more
interesting, like topological transitions as in the corresponding phenomena in
matrix models. This disadvantage could be dealt with by simultaneously taking
the large $R$ limit, which increases the number of overlap variables in
correlation with $N$. In fact, there are some suggestions that urge us to look
into this direction:

  * •

In [23], it was argued that the system (1) is convergent, only if
$R\lesssim(N+1)(N+2)/2$. Therefore it would be interesting to take a large-$N$
limit with $R=N^{2}\,r$ with fixed $r$, and study the dependence on $r$. Some
intereresting phenomena are expected to occur near $r=1/2$.

  * •

In the tensor model, the hermiticity of the Hamiltonian constraint requires
$R=(N+2)(N+3)/4$ [31]141414 As for the derivation of this value of $R$, it
would be easier to see an appendix of [10] for a short summary with the same
normalization as is used in this paper. When the wave function is absolute-
squared for probability, the value which matters is given by $2R=(N+2)(N+3)/2$
[25], which curiously reproduces $r=1/2$ as an interesting case.. This urges
us to look into the same kind of limit as above.

  * •

In [24], a different but closely related model was treated, and it was found
that there is a transit parameter region, $N\lesssim R\lesssim N^{2}/2$, in
which dominant graphs gradually change. This would suggest that not only the
above one but also $R=N\,r^{\prime}$ with fixed $r^{\prime}$ could provide an
interesting large-$N$ limit.

  * •

From the view point of the tensor rank decomposition, Alexander-Hirschowitz
theorem[32] tells that, with some exceptions, the complex symmetric general
rank of order-three symmetric tensors is given by $\lceil(N+1)(N+2)/6\rceil$,
and combining with the argument in [33], the symmetric real rank of a
generally given real tensor of order three should be smaller than or equal to
$3\lceil(N+1)(N+2)/6\rceil$. This number curiously agrees well with the number
given in the first item above, and gives a motivation to study the same
large-$N$ limit.

  * •

The results in Figure 13 imply that $\langle\psi^{1}\cdot\psi^{2}\rangle$
could vanish in the 1RSB region only after the large-$N$ limit with fixed $R$
has been taken. Since this vanishing is the main cause for the reduction of
$R>1$ to $R=1$, we would expect more interesting dynamics to appear, if this
vanishing is changed in the simultaneous large-$N,R$ limits proposed above.

It is stimulating that the different perspectives above actually point to the
similar large-$N,R$ limits. We hope we can report some progress in this
direction in future studies.

## Acknowledgements

The work of N.S. is supported in part by JSPS KAKENHI Grant No.19K03825.

##  Appendix Appendix A Rank of $SO(2)$ invariant $C$

A fact used in Section 3.3 is that the rank of $SO(2)$ invariant $C$ of $n=3$
is given by four. This is a special case of the following general formula:

$\displaystyle\hbox{Rank}(C^{SO(2)})=\frac{3n-1}{2}.$ (55)

This statement is equivalent to that $r=(3n-1)/2$ is the minimum number for
the following summation expression of the integral to hold:

$\displaystyle\begin{split}\int_{0}^{2\pi}d\theta
f_{a}f_{b}f_{c}=\sum_{j=1}^{r}A_{a}^{j}A_{b}^{j}A_{c}^{j},\end{split}$ (56)

where

$\displaystyle\\{f_{a}|a=1,2,\ldots,2\Lambda+1\\}=\\{\frac{1}{\sqrt{2}},\cos\theta,\sin\theta,\cos
2\theta,\cdots,\sin(\Lambda\theta)\\},$ (57)

with $n=2\Lambda+1$.

To prove (56) let us recall what is called quadrature [34]. Quadrature is an
approximation to an integral in terms of a summation, which however gives
exact values if the integrand is contained in a particular set of functions.
In the present case, we want a quadrature,

$\displaystyle\int_{0}^{2\pi}d\theta\
g(\theta)=\sum_{j=1}^{r}w_{j}\,g(\theta_{j}),$ (58)

where $g(\theta)$ is contained in a certain set of periodic functions of
period $2\pi$, $\theta_{j}$ are quadrature nodes, and $w_{j}$ are quadrature
weights. Let us consider the set of $g(\theta)$ to be given by

$\displaystyle\\{1,\cos\theta,\sin\theta,\cos 2\theta,\ldots,\sin M\theta\\}.$
(59)

Then it is straightforward to show that the nodes can be taken as

$\displaystyle\theta_{j}=\frac{2\pi j}{M+1},\ (j=0,1,\ldots,M),$ (60)

with weights $w_{j}=2\pi/(M+1)$, because

$\displaystyle\frac{2\pi}{M+1}\sum_{j=0}^{M}e^{i\frac{2\pi
kj}{M+1}}=\delta_{k,0}=\int_{0}^{2\pi}d\theta e^{ik\theta},$ (61)

for integer $0\leq k\leq M$. Here the number of nodes, $r=M+1$, is the least,
because $M+1$ independent conditions must be satisfied.

In our present case, $g$ is given by a product of three functions in (57).
Therefore $M=3\Lambda$, and hence $r=M+1=(3n-1)/2$, which indeed agrees with
(55). By applying (58) to (56), we obtain

$\displaystyle
A_{a}^{j}=\left(\frac{2\pi}{3\Lambda+1}\right)^{\frac{1}{3}}f_{a}(\theta_{j}).$
(62)

##  Appendix Appendix B Details of the derivation of (50)

In this appendix, we show some details of the derivation of the free energy
(50).

From the expression (49), we find

$\displaystyle\begin{split}\sum_{ii^{\prime}tt^{\prime}}(\bar{Q}^{iti^{\prime}t^{\prime}})^{3}&=R\sum_{tt^{\prime}}(\bar{Q}^{1t1t^{\prime}})^{3}+R(R-1)\sum_{tt^{\prime}}(\bar{Q}^{1t2t^{\prime}})^{3}\\\
&=RT(q_{0}^{3}+(M-1)q_{1}^{3})+R(R-1)T((q_{0}^{\prime})^{3}+(M-1)(q_{1}^{\prime})^{3}).\end{split}$
(63)

We also obtain

$\displaystyle\begin{split}\sum_{ii^{\prime}t}(\bar{Q}^{iti^{\prime}t})^{3}&=R\sum_{t}(\bar{Q}^{1t1t})^{3}+R(R-1)\sum_{t}(\bar{Q}^{1t2t})^{3}\\\
&=RTq_{0}^{3}+R(R-1)T(q_{0}^{\prime})^{3}.\end{split}$ (64)

Computation of $\det(\bar{Q})$ is also straightforward but a little involved.
First we note that $\bar{Q}$ has the form,

$\displaystyle\bar{Q}^{iti^{\prime}t^{\prime}}=A_{tt^{\prime}}\delta_{ii^{\prime}}+A^{\prime}_{tt^{\prime}}(1-\delta_{ii^{\prime}}),$
(65)

where

$\displaystyle A_{tt^{\prime}}=\delta_{\lfloor t/M\rfloor,\lfloor
t^{\prime}/M\rfloor}I_{t\bmod M,t^{\prime}\bmod M}(q_{0},q_{1}),$ (66)

and similarly for $A^{\prime}$ with $q_{0}^{\prime}$ and $q_{1}^{\prime}$.
Concerning the indices $i,i^{\prime}$, $\bar{Q}$ in (65) has the
“eigenvalues”, $A+(R-1)A^{\prime}$ and $A-A^{\prime}$ with degeneracies $1$
and $R-1$, respectively. On the other hand, $A$ has the eigenvalues,
$q_{0}+(M-1)q_{1}$ and $q_{0}-q_{1}$ with degeneracies $T/M$ and $(M-1)T/M$,
respectively. Combining these, we obtain the following list of eigenvalues of
$\bar{Q}$:

$\displaystyle\begin{split}&q_{0}+(R-1)q_{0}^{\prime}+(M-1)(q_{1}+(R-1)q_{1}^{\prime}),\
\ \hbox{deg}=T/M,\\\ &q_{0}+(R-1)q_{0}^{\prime}-(q_{1}+(R-1)q_{1}^{\prime}),\
\ \hbox{deg}=(M-1)T/M,\\\ &q_{0}-q_{0}^{\prime}+(M-1)(q_{1}-q_{1}^{\prime}),\
\ \hbox{deg}=(R-1)T/M,\\\ &q_{0}-q_{0}^{\prime}-(q_{1}-q_{1}^{\prime}),\ \
\hbox{deg}=(M-1)(R-1)T/M.\end{split}$ (67)

This determines the terms coming from $\det\bar{Q}$.

## References

  * [1] N. Sasakura, “Canonical tensor models with local time,” Int. J. Mod. Phys. A 27, 1250020 (2012) doi:10.1142/S0217751X12500200 [arXiv:1111.2790 [hep-th]].
  * [2] N. Sasakura, “Uniqueness of canonical tensor model with local time,” Int. J. Mod. Phys. A 27, 1250096 (2012) doi:10.1142/S0217751X12500960 [arXiv:1203.0421 [hep-th]].
  * [3] J. Ambjorn, B. Durhuus and T. Jonsson, “Three-dimensional simplicial quantum gravity and generalized matrix models,” Mod. Phys. Lett. A 6, 1133-1146 (1991) doi:10.1142/S0217732391001184
  * [4] N. Sasakura, “Tensor model for gravity and orientability of manifold,” Mod. Phys. Lett. A 6, 2613-2624 (1991) doi:10.1142/S0217732391003055
  * [5] N. Godfrey and M. Gross, “Simplicial quantum gravity in more than two-dimensions,” Phys. Rev. D 43, 1749-1753 (1991) doi:10.1103/PhysRevD.43.R1749
  * [6] R. Gurau, “Colored Group Field Theory,” Commun. Math. Phys. 304, 69-93 (2011) doi:10.1007/s00220-011-1226-9 [arXiv:0907.2582 [hep-th]].
  * [7] V. Bonzom, R. Gurau, A. Riello and V. Rivasseau, “Critical behavior of colored tensor models in the large N limit,” Nucl. Phys. B 853, 174-195 (2011) doi:10.1016/j.nuclphysb.2011.07.022 [arXiv:1105.3122 [hep-th]].
  * [8] R. Gurau and J. P. Ryan, “Colored Tensor Models - a review,” SIGMA 8, 020 (2012) doi:10.3842/SIGMA.2012.020 [arXiv:1109.4812 [hep-th]].
  * [9] J. Ambjorn, J. Jurkiewicz and R. Loll, “Emergence of a 4-D world from causal quantum gravity,” Phys. Rev. Lett. 93, 131301 (2004) doi:10.1103/PhysRevLett.93.131301 [arXiv:hep-th/0404156 [hep-th]].
  * [10] T. Kawano and N. Sasakura, “Emergence of Lie group symmetric classical spacetimes in the canonical tensor model,” PTEP 2022, no.4, 043A01 (2022) doi:10.1093/ptep/ptac045 [arXiv:2109.09896 [hep-th]].
  * [11] D. J. Gross and E. Witten, “Possible Third Order Phase Transition in the Large N Lattice Gauge Theory,” Phys. Rev. D 21, 446-453 (1980) doi:10.1103/PhysRevD.21.446
  * [12] S. R. Wadia, “$N$ = Infinity Phase Transition in a Class of Exactly Soluble Model Lattice Gauge Theories,” Phys. Lett. B 93, 403-410 (1980) doi:10.1016/0370-2693(80)90353-6
  * [13] B. Eynard, Counting Surfaces, Prog.Math.Phys. 70, Springer, 2016, doi:10.1007/978-3-7643-8797-6
  * [14] F. L. Hitchcock, “The expression of a tensor or a polyadic as a sum of products,” Journal of Mathematics and Physics 6 no. 1-4, (1927) 164–189. http://dx.doi.org/10.1002/sapm192761164.
  * [15] J. D. Carroll and J.-J. Chang, “Analysis of individual differences in multidimensional scaling via an n-way generalization of “eckart-young” decomposition,” Psychometrika 35 no. 3, (Sep, 1970) 283–319. https://doi.org/10.1007/BF02310791.
  * [16] Landsberg, J. M., Tensors: Geometry and Applications.  American Mathematical Society, Providence, 2012.
  * [17] P. Comon, “Tensors: a Brief Introduction,” IEEE Signal Processing Magazine 31 no. 3, (May, 2014) 44–53. https://hal.archives-ouvertes.fr/hal-00923279.
  * [18] M. Ouerfelli, V. Rivasseau and M. Tamaazousti, “The Tensor Track VII: From Quantum Gravity to Artificial Intelligence,” [arXiv:2205.10326 [hep-th]].
  * [19] C. J. Hillar and L. Lim, “Most tensor problems are NP-Hard”, Journal of the ACM. 60 (6) (2013): 1-39. arXiv:0911.1393. doi:10.1145/2512329
  * [20] W. Hackbusch, Tensor Spaces and Numerical Tensor Calculus, Springer-Verlag Berlin Heidelberg 2012, https://doi.org/10.1007/978-3-642-28027-6.
  * [21] A. Crisanti and H.-J. Sommers, “The spherical p-spin interaction spin glass model: the statics”, Z. Phys. B 87, 341 (1992).
  * [22] T. Castellani and A. Cavagna, “Spin-glass theory for pedestrians”, J. Stat. Mech.: Theo. Exp. 2005, 05012 [arXiv: cond-mat/0505032].
  * [23] D. Obster and N. Sasakura, “Counting Tensor Rank Decompositions,” Universe 7, no.8, 302 (2021) doi:10.3390/universe7080302 [arXiv:2107.10237 [gr-qc]].
  * [24] L. Lionni and N. Sasakura, “A random matrix model with non-pairwise contracted indices,” PTEP 2019, no.7, 073A01 (2019) doi:10.1093/ptep/ptz057 [arXiv:1903.05944 [hep-th]].
  * [25] N. Sasakura and S. Takeuchi, “Numerical and analytical analyses of a matrix model with non-pairwise contracted indices,” Eur. Phys. J. C 80, no.2, 118 (2020) doi:10.1140/epjc/s10052-019-7591-9 [arXiv:1907.06137 [hep-th]].
  * [26] D. Obster and N. Sasakura, “Phases of a matrix model with non-pairwise index contractions,” PTEP 2020, no.7, 073B06 (2020) doi:10.1093/ptep/ptaa085 [arXiv:2004.03152 [hep-th]].
  * [27] P. Comon and G. Ottaviani, “On the typical rank of real binary forms,” Linear and Multilinear Algebra 60, 657 (2012) https://doi.org/10.1080/03081087.2011.624097, [arXiv:0909.4865 [math.AG]].
  * [28] R. Neal, Handbook of Markov Chain Monte Carlo, 113 (2011). doi:10.1201/b10905 [arXiv:1206.1901 [stat.CO]].
  * [29] D. J. Earl and M. W. Deem, “Parallel tempering: Theory, applications, and new perspectives”. Phys. Chem. Chem. Phys. 7: 3910 (2005) doi:10.1039/b509983h.
  * [30] T. Kawano, D. Obster and N. Sasakura, “Canonical tensor model through data analysis: Dimensions, topologies, and geometries,” Phys. Rev. D 97, no.12, 124061 (2018) doi:10.1103/PhysRevD.97.124061 [arXiv:1805.04800 [hep-th]].
  * [31] N. Sasakura, “Quantum canonical tensor model and an exact wave function,” Int. J. Mod. Phys. A 28, 1350111 (2013) doi:10.1142/S0217751X1350111X [arXiv:1305.6389 [hep-th]].
  * [32] J. Alexander and A. Hirschowitz, “Polynomial interpolation in several variables,” J. Algebraic Geom., 4, 201-222 (1995).
  * [33] E. Ballico, “An upper bound for the real tensor rank and the real symmetric tensor rank in terms of the complex ranks,” Linear and Multilinear Algebra 62, 1546-1552 (2014) doi:10.1080/03081087.2013.839671
  * [34] L. N. Trefethen, Approximation Theory and Approximation Practice, Extended Edition, Other Titles in Applied Mathematics, SIAM Textbooks, 2019, https://doi.org/10.1137/1.9781611975949.