#  YITP-22-123
Exact analytic expressions of real tensor eigenvalue
distributions of Gaussian tensor model for small $N$

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

We obtain exact analytic expressions of real tensor eigenvalue/vector
distributions of real symmetric order-three tensors with Gaussian
distributions for $N\leq 8$. This is achieved by explicitly computing the
partition function of a zero-dimensional boson-fermion system with four-
interactions. The distributions are expressed by combinations of polynomial,
exponential and error functions as results of feasible complicated bosonic
integrals which appear after fermionic integrations. By extrapolating the
expressions and also using a previous result, we guess a large-$N$ expression.
The expressions are compared with Monte Carlo simulations, and precise and
good agreement are obtained with the exact and the large-$N$ expressions,
respectively. Understanding the feasibility of the integration is left for
future study, which would provide a general-$N$ analytic formula.

##  1 Introduction

Eigenvalue distributions are important dynamical quantities in studies of
matrix models. They model energy eigenvalue distributions of complex dynamical
systems [1]. They provide a major technique in solving matrix models [2].
Their topological properties differentiate phases of matrix models222For more
details, see for instance [3]., providing insights into the dynamics of gauge
theories [4, 5].

Recently, tensor models [6, 7, 8, 9] attract much attention in various
contexts [10]. While it is important to develop efficient techniques of
computing eigenvalues/vectors for certain tensors in various practical
applications of tensors [11], it is also interesting to study their
distributions for ensembles of tensors, because tensors are dynamical in
tensor models.

While there had already been a few interesting closely related results [12,
13, 14, 15], eigenvalue distributions in tensor models had not been explicitly
studied until very recently, when the present author derived some explicit
formulas [16, 17]. The general strategy employed in these studies was that the
problem was rewritten as computations of partition functions of zero-
dimensional fermion systems. In [16], an exact formula of signed distributions
of real eigenvalue/vector distributions for real symmetric order-three tensors
with Gaussian distributions was obtained, where each eigenvalue/vector
contributed to the distribution by $\pm 1$, depending on the sign of a Hessian
matrix associated to each eigenvalue/vector. In [17], real eigenvalue/vector
distributions for the same ensemble of tensors was studied. In this study,
however, only an approximate large-$N$ expression333$N$ denotes the dimension
of the index vector space of the tensor. was obtained by truncating a
Schwinger-Dyson equation, though some closely related exact results were also
obtained. While the functional form of the large-$N$ expression, which was a
Gaussian, agreed well with Monte Carlo simulations, the overall factor did not
and an improvement was expected.

In this paper, we rewrite the real eigenvalue/vector problem as a computation
of a partition function of a boson-fermion system with four-interactions,
instead of the fermion systems in the previous studies [16, 17]. We exactly
compute the real eigenvalue/vector distributions for real symmetric order-
three tensors with Gaussian distributions for $N\leq 8$, and extrapolate the
expressions to guess a large-$N$ expression, including the overall factor,
which was not correctly obtained by the approximation of the previous paper
[17]. We compare the exact expressions for small $N$ with Monte Carlo
simulations, and obtain precise agreement. We also compare the large-$N$
expression, and obtain good agreement.

This paper is organized as follows. In Section 2, we define what we compute,
the real eigenvalue/vector distributions for real symmetric order-three
tensors with Gaussian distributions. In Section 3, we rewrite the problem as
computation of a partition function of a zero-dimensional boson-fermion system
with four-interactions. In Section 4, we determine the general form which is
obtained by integrating out the fermions in the partition function. In Section
5, we perform the remaining bosonic integrals but one. In Section 6, we
explicitly perform the last bosonic integral for each case of $N\leq 8$, and
obtain the explicit exact expressions of the distributions. Interestingly, the
apparently complicated bosonic integrals are feasible and the final results
have simple expressions. In Section 7, we perform some Monte Carlo
simulations, and precise agreement with the exact expressions is obtained. In
Section 8, we perform an extrapolation of the exact results for small $N$ to
guess a large-$N$ expression. We compare it with Monte Carlo simulations, and
good agreement is obtained. The last section is devoted to a summary and
future prospects.

##  2 Real eigenvector/value distributions of tensors

In this paper, we restrict ourselves to the real symmetric order-three
tensors, i.e., $C_{abc}\in\mathbb{R}\ (a,b,c=1,2,\ldots,N)$, satisfying
$C_{abc}=C_{bac}=C_{bca}$, as the simplest case. As for eigenvalues/vectors of
tensors, there are a few similar but slightly different definitions in the
literature [18, 19, 20]. In this paper we employ the definition that the real
eigenvectors $v$ of a given $C$ are the non-zero solutions to

$\displaystyle C_{abc}v_{b}v_{c}=v_{a},\ (v\neq 0,\ v\in\mathbb{R}^{N}),$ (1)

where the repeated indices are assumed to be summed over, as will be assumed
in the rest of this paper. Then the distribution of the eigenvectors for a
given $C$ is given by

$\displaystyle\begin{split}\rho(v,C)&=\sum_{i=1}^{n_{C}}\delta^{N}(v-v^{i})\\\
&=|\det M|\prod_{a=1}^{N}\delta(v_{a}-C_{abc}v_{b}v_{c})\end{split}$ (2)

for $v\neq 0$ under the volume measure $d^{N}v=\prod_{a=1}^{N}dv_{a}$, where
$v^{i}\ (i=1,2,\ldots,n_{C})$ are all the solutions to (1), $|\cdot|$ denotes
the absolute value, and

$\displaystyle M_{ab}=\frac{\partial}{\partial
v_{a}}\left(v_{b}-C_{bcd}v_{c}v_{d}\right)=\delta_{ab}-2C_{abc}v_{c}.$ (3)

Here the absolute value of the determinant, $|{\rm det}M|$, is the Jacobian
associated to the change of the arguments of the $\delta$-functions performed
in (2).

When the tensor $C$ has a Gaussian distribution, the real eigenvector
distribution is given by

$\displaystyle\rho(v)=A^{-1}\int_{\mathbb{R}^{\\#C}}dC\,e^{-\alpha
C^{2}}\rho(v,C),$ (4)

where $dC=\prod_{a\leq b\leq c=1}^{N}dC_{abc}$, $A=\int_{\\#C}dC\,e^{-\alpha
C^{2}}$, $C^{2}=C_{abc}C_{abc}$, $\alpha>0$, and $\\#C=N(N+1)(N+2)/6$, which
is the total number of the independent components of $C$.

It is worth commenting on the eigenvalue distribution corresponding to the
real eigenvector distribution above. An eigenvalue $\zeta$ accompanied with a
real eigenvector $v$ (a Z-eigenvalue in the terminology of [18]) is defined by

$\displaystyle C_{abc}w_{b}w_{c}=\zeta\,w_{a}\ (|w|=1,\ w\in\mathbb{R}^{N}),$
(5)

where $|w|=\sqrt{w_{a}w_{a}}$. Comparing with (1), we obtain the relation,

$\displaystyle\zeta=\frac{1}{|v|}.$ (6)

Therefore, by using (6) and the fact that $\rho(v)$ is actually a function of
$|v|$ because of the rotational symmetry of the distribution, we obtain the
real eigenvalue distribution,

$\displaystyle\rho_{\rm eig}(\zeta)=\rho(1/\zeta)S_{N-1}\zeta^{-N-1},$ (7)

where $1/\zeta$ in the argument abusively represents an arbitrary vector of
size $1/\zeta$, and $S_{N-1}=2\pi^{N/2}/\Gamma[N/2]$, the surface volume of a
unit sphere in an $N$-dimensional space.

##  3 A zero-dimensional boson-fermion system

In this section, we will rewrite (4) with (2) as a partition function of a
zero-dimensional boson-fermion system with four-interactions. An immediate
obstacle in doing this is the presence of an absolute value in (2), which is
not an analytic function. To rewrite it in an analytic form, we take

$\displaystyle|\det M|=\lim_{\epsilon\rightarrow+0}\frac{\det(M^{2}+\epsilon
I)}{\sqrt{\det(M^{2}+\epsilon I)}},$ (8)

where $I_{ab}=\delta_{ab}$ is an identity matrix, and the parameter $\epsilon$
is a positive small regularization parameter, which assures the convergence of
the integrals below.

As is well known, $\det(M^{2}+\epsilon I)$ can be rewritten as
$\det(M^{2}+\epsilon I)=\int
d{\bar{\psi}}d\psi\,e^{{\bar{\psi}}_{a}(M^{2})_{ab}\psi_{b}+\epsilon{\bar{\psi}}_{a}\psi_{a}}$
by introducing a fermion pair, ${\bar{\psi}}$ and $\psi$ [21]. However, the
exponent of this expression contains $C$ in a quadratic manner (see (3)) and
is difficult to handle when we perform an integration over $C$ as will be done
below. Therefore, as was done previously in [17], we introduce another fermion
pair, $\bar{\varphi},{\varphi}$, to rewrite the exponent in a form linear in
$C$:

$\displaystyle\det(M^{2}+\epsilon I)=(-1)^{N}\int d{\bar{\psi}}d\psi
d\bar{\varphi}d{\varphi}\,e^{-\bar{\varphi}{\varphi}-{\bar{\psi}}M{\varphi}-\bar{\varphi}M\psi-\epsilon{\bar{\psi}}\psi},$
(9)

where the contracted indices are suppressed for brevity:
$\bar{\varphi}{\varphi}=\bar{\varphi}_{a}{\varphi}_{a}$,
${\bar{\psi}}M{\varphi}={\bar{\psi}}_{a}M_{ab}{\varphi}_{b}$, and so on. The
equality can be shown by noting that
$\bar{\varphi}{\varphi}+{\bar{\psi}}M{\varphi}+\bar{\varphi}M\psi=(\bar{\varphi}+{\bar{\psi}}M)({\varphi}+M\psi)-{\bar{\psi}}M^{2}\psi$.
In a similar manner we obtain

$\displaystyle\frac{1}{\sqrt{\det(M^{2}+\epsilon
I)}}=\pi^{-N}\int_{\mathbb{R}^{2N}}d\phi d\sigma\,e^{-\sigma^{2}-2i\sigma
M\phi-\epsilon\phi^{2}},$ (10)

where $\sigma$ is a new bosonic variable, and
$\sigma^{2}=\sigma_{a}\sigma_{a}$, and so on.

By using (8), (9), (10) and a well-known formula,
$\int_{\mathbb{R}}dx\,e^{ipx}=2\pi\delta(p)$, (4) with (2) can be rewritten as

$\displaystyle\rho(v)=\lim_{\epsilon\rightarrow+0}A^{-1}(2\pi)^{-N}\pi^{-N}(-1)^{N}\int
dCd\lambda d\phi d\sigma d{\bar{\psi}}d\psi
d\bar{\varphi}d{\varphi}\,e^{S_{1}},$ (11)

where

$\displaystyle S_{1}=-\alpha
C^{2}+i\lambda_{a}(v_{a}-C_{abc}v_{b}v_{c})-\sigma^{2}-2i\sigma
M\phi-\epsilon\phi^{2}-\bar{\varphi}{\varphi}-{\bar{\psi}}M{\varphi}-\bar{\varphi}M\psi-\epsilon{\bar{\psi}}\psi$
(12)

with a new bosonic variable $\lambda_{a}\ (a=1,2,\ldots,N)$ to rewrite the
$\delta$-functions.

Next, let us perform the integrations over $C$ and $\lambda$ in the expression
(11). Let us first consider $C$. The terms containing $C$ in (12) are

$\displaystyle S_{C}=-\alpha
C^{2}-i\,C_{abc}\lambda_{a}v_{b}v_{c}+4iC_{abc}v_{a}\sigma_{b}\phi_{c}+2C_{abc}v_{a}{\bar{\psi}}_{b}{\varphi}_{c}+2C_{abc}v_{a}\bar{\varphi}_{b}\psi_{c}.$
(13)

Then we obtain

$\displaystyle\int_{\mathbb{R}^{\\#C}}dC\,e^{S_{C}}=A\,e^{\delta S_{C}},$ (14)

where

$\displaystyle\delta
S_{C}=\frac{1}{\alpha}\left(\sum_{s}\left(-\frac{i}{2}\lambda_{s_{a}}v_{s_{b}}v_{s_{c}}+2iv_{s_{a}}\sigma_{s_{b}}\phi_{s_{c}}+v_{s_{a}}{\bar{\psi}}_{s_{b}}{\varphi}_{s_{c}}+v_{s_{a}}\bar{\varphi}_{s_{b}}\psi_{s_{c}}\right)\right)^{2},$
(15)

where the summation is over all the permutations of $a,b,c$, which is
necessary because $C$ is a symmetric tensor. By explicitly expanding (15), we
obtain

$\displaystyle\begin{split}\delta
S_{C}&=-\frac{|v|^{4}}{12\alpha}\lambda_{a}B_{ab}\lambda_{b}-i\lambda_{a}(D_{a}+\tilde{D}_{a})+E_{1}+E_{2}+E_{3},\end{split}$
(16)

where

$\displaystyle\begin{split}B_{ab}&=\delta_{ab}+2\hat{v}_{a}\hat{v}_{b}=I_{\perp\,ab}+3I_{\parallel\,ab},\\\
D_{a}&=\frac{|v|^{3}}{3\alpha}\left({\bar{\psi}}_{\parallel}{\varphi}_{\parallel}{}+\bar{\varphi}_{\parallel}\psi_{\parallel}{}\right)\hat{v}_{a}+\frac{|v|^{3}}{3\alpha}\left({\bar{\psi}}_{a}{\varphi}_{\parallel}{}+{\bar{\psi}}_{\parallel}{\varphi}_{a}+\bar{\varphi}_{a}\psi_{\parallel}{}+\bar{\varphi}_{\parallel}\psi_{a}\right),\\\
\tilde{D}_{a}&=\frac{2i}{\alpha}v_{b}v_{c}\cdot\frac{1}{6}\sum_{s}v_{s_{a}}\sigma_{s_{b}}\phi_{s_{c}},\\\
E_{1}&=\frac{1}{\alpha}\left(\frac{1}{6}\sum_{s}\left({\bar{\psi}}_{s_{a}}{\varphi}_{s_{b}}v_{s_{c}}+\bar{\varphi}_{s_{a}}\psi_{s_{b}}v_{s_{c}}\right)\right)^{2},\\\
E_{2}&=-\frac{4}{\alpha}\left(\frac{1}{6}\sum_{s}v_{s_{a}}\sigma_{s_{b}}\phi_{s_{c}}\right)^{2},\\\
E_{3}&=\frac{4i}{\alpha}(v_{a}{\bar{\psi}}_{b}{\varphi}_{c}+v_{a}\bar{\varphi}_{b}\psi_{c})\cdot\frac{1}{6}\sum_{s}v_{s_{a}}\sigma_{s_{b}}\phi_{s_{c}}.\end{split}$
(17)

Here $\hat{v}=v/|v|$, and
${\bar{\psi}}_{\parallel}=\bar{\psi}_{a}\hat{v}_{a}$, etc., and
$I_{\parallel}$ and $I_{\perp}$ are respectively the projection matrices to
the parallel and transverse subspaces to $\hat{v}$.

Next let us perform the integration over $\lambda$. Picking up the terms
containing $\lambda$ (with no $C$) in (12) and (16), we obtain

$\displaystyle
S_{\lambda}=-\frac{|v|^{4}}{12\alpha}\lambda_{a}B_{ab}\lambda_{b}+i\lambda_{a}(v_{a}-D_{a}-\tilde{D}_{a}).$
(18)

Considering that $B$ is a sum of the projection matrices as in (17), we obtain

$\displaystyle\int_{\mathbb{R}^{N}}d\lambda\,e^{S_{\lambda}}=|v|^{-2N}(12\pi\alpha)^{\frac{N}{2}}(\det
B)^{-\frac{1}{2}}e^{\delta S_{\lambda}},$ (19)

where $\det B=3$ from (17), and

$\displaystyle\begin{split}\delta
S_{\lambda}&=-3\alpha|v|^{-4}(v-D-\tilde{D})_{a}B^{-1}_{ab}(v-D-\tilde{D})_{b}\\\
&=-\alpha|v|^{-2}+2\alpha|v|^{-3}D_{\parallel}+2\alpha|v|^{-3}\tilde{D}_{\parallel}\\\
&\ \ \ \ -3\alpha|v|^{-4}\left(D_{\perp}\cdot
D_{\perp}+\frac{1}{3}D_{\parallel}^{2}+2D_{\perp}\cdot\tilde{D}_{\perp}+\tilde{D}_{\perp}^{2}+\frac{2}{3}D_{\parallel}\tilde{D}_{\parallel}+\frac{1}{3}\tilde{D}_{\parallel}^{2}\right).\end{split}$
(20)

Here we have used $B^{-1}=I_{\perp}+\frac{1}{3}I_{\parallel}$, and
$D_{\perp}\cdot\tilde{D}_{\perp}=D_{\perp\,a}\tilde{D}_{\perp\,a}$, and so on,
where $X_{\perp}\ (X=D,\tilde{D})$ denotes vector $X$ projected to the
transverse subspace to $\hat{v}$. From (17) the projections of $D$ and
$\tilde{D}$ are more explicitly given by

$\displaystyle\begin{split}D_{\parallel}&=\frac{|v|^{3}}{\alpha}\left({\bar{\psi}}_{\parallel}{\varphi}_{\parallel}+\bar{\varphi}_{\parallel}\psi_{\parallel}\right),\\\
D_{\perp}&=\frac{|v|^{3}}{3\alpha}\left(\bar{\psi}_{\perp}{\varphi}_{\parallel}+{\bar{\psi}}_{\parallel}\varphi_{\perp}{}+\bar{\varphi}_{\perp}\psi_{\parallel}{}+\bar{\varphi}_{\parallel}\psi_{\perp}\right),\\\
\tilde{D}_{\parallel}&=\frac{2i|v|^{3}}{\alpha}\sigma_{\parallel}\phi_{\parallel},\\\
\tilde{D}_{\perp}&=\frac{2i|v|^{3}}{3\alpha}\left(\sigma_{\parallel}\phi_{\perp}+\sigma_{\perp}\phi_{\parallel}\right).\end{split}$
(21)

The first term in the second line of (20) is a constant contribution to the
potential. The second and the third terms in the same line are corrections to
the kinetic terms of the parallel components of the fermions and the bosons,
respectively. The other terms describe the four-interaction terms among the
bosons and the fermions. Collecting the results above, we obtain the following
intermediate expression,

$\displaystyle\rho(v)=\lim_{\epsilon\rightarrow+0}3^{\frac{N-1}{2}}\pi^{-\frac{3N}{2}}\alpha^{\frac{N}{2}}|v|^{-2N}e^{-\frac{\alpha}{|v|^{2}}}(-1)^{N}\int
d\phi d\sigma d{\bar{\psi}}d\psi
d\bar{\varphi}d{\varphi}\,e^{\tilde{K}_{F}+\tilde{K}_{B}+\tilde{V}_{F}+\tilde{V}_{B}+\tilde{V}_{BF}},$
(22)

where

$\displaystyle\begin{split}\tilde{K}_{F}&=-\bar{\varphi}_{\perp}\cdot\varphi_{\perp}-\bar{\psi}_{\perp}\cdot\varphi_{\perp}-\bar{\varphi}_{\perp}\cdot\psi_{\perp}-\epsilon\bar{\psi}_{\perp}\cdot\psi_{\perp}-\bar{\varphi}_{\parallel}{\varphi}_{\parallel}+{\bar{\psi}}_{\parallel}{\varphi}_{\parallel}+\bar{\varphi}_{\parallel}\psi_{\parallel}-\epsilon{\bar{\psi}}_{\parallel}\psi_{\parallel},\\\
\tilde{K}_{B}&=-{\sigma_{\perp}}^{2}-2i{\sigma_{\perp}}\cdot{\phi_{\perp}}-\epsilon{\phi_{\perp}}^{2}-{\sigma_{\parallel}}^{2}+2i{\sigma_{\parallel}}{\phi_{\parallel}}-\epsilon{\phi_{\parallel}}^{2},\\\
\tilde{V}_{F}&=E_{1}-3\alpha|v|^{-4}\left(D_{\perp}\cdot
D_{\perp}+\frac{1}{3}D_{\parallel}^{2}\right),\\\
\tilde{V}_{B}&=E_{2}-3\alpha|v|^{-4}\left(\tilde{D}_{\perp}^{2}+\frac{1}{3}\tilde{D}_{\parallel}^{2}\right),\\\
\tilde{V}_{BF}&=E_{3}-3\alpha|v|^{-4}\left(2D_{\perp}\cdot\tilde{D}_{\perp}+\frac{2}{3}D_{\parallel}\tilde{D}_{\parallel}\right).\end{split}$
(23)

Here $\tilde{K}_{F},\tilde{V}_{F}$ contain only the fermions,
$\tilde{K}_{B},\tilde{V}_{B}$ only the bosons, and $\tilde{V}_{BF}$ is a
mixture. Note that the quadratic terms of the parallel components in
$\tilde{K}_{F}$ and $\tilde{K}_{B}$ have been corrected by the aforementioned
terms from the second line of (20).

Now we want to compute $\tilde{V}_{B},\tilde{V}_{F},\tilde{V}_{BF}$ more
explicitly. The computation of $\tilde{V}_{F}$ in (23) is essentially the same
as the derivation of the four-fermi interaction terms in [17] for $R=1$. By
noting that two of the terms in [17] do not appear, because
$\bar{\psi}_{\perp}\cdot\bar{\psi}_{\perp}=\bar{\varphi}_{\perp}\cdot\bar{\varphi}_{\perp}=0$
(for $R=1$), we obtain

$\displaystyle\begin{split}\tilde{V}_{F}&=-\frac{|v|^{2}}{6\alpha}\left((\bar{\psi}_{\perp}\cdot\varphi_{\perp})^{2}+(\bar{\varphi}_{\perp}\cdot\psi_{\perp})^{2}+2\bar{\psi}_{\perp}\cdot\bar{\varphi}_{\perp}\varphi_{\perp}\cdot\psi_{\perp}+2\bar{\psi}_{\perp}\cdot\psi_{\perp}{}\bar{\varphi}_{\perp}\cdot\varphi_{\perp}{}\right).\end{split}$
(24)

What seems surprising in this expression is that the parallel components of
the fermions cancel out from the interactions after all.

As for $\tilde{V}_{B}$ and $\tilde{V}_{BF}$ in (23), by using the results of
the computations of $E_{2}$ and $E_{3}$ in Appendix A, we obtain

$\displaystyle\begin{split}\tilde{V}_{B}&=-\frac{2|v|^{2}}{3\alpha}\left({\sigma_{\perp}}^{2}{\phi_{\perp}}^{2}+({\sigma_{\perp}}\cdot{\phi_{\perp}})^{2}\right),\\\
\tilde{V}_{BF}&=\frac{2i|v|^{2}}{3\alpha}\left(\bar{\psi}_{\perp}\cdot{\sigma_{\perp}}\,\varphi_{\perp}\cdot{\phi_{\perp}}+\bar{\varphi}_{\perp}\cdot{\sigma_{\perp}}\,\psi_{\perp}\cdot{\phi_{\perp}}+\bar{\psi}_{\perp}\cdot{\phi_{\perp}}\,\varphi_{\perp}\cdot{\sigma_{\perp}}+\bar{\varphi}_{\perp}\cdot{\phi_{\perp}}\,\psi_{\perp}\cdot{\sigma_{\perp}}\right).\end{split}$
(25)

We again find the surprising fact that no parallel components appear in
$\tilde{V}_{B}$ and $\tilde{V}_{BF}$.

Because the parallel components only exist in $\tilde{K}_{B},\tilde{K}_{F}$
and do not interact, these can trivially be integrated out, that generates the
overall factors of $\pi$ and $-1$ for the bosons and fermions, respectively.
Therefore we finally obtain

$\displaystyle\rho(v)=\lim_{\epsilon\rightarrow+0}3^{\frac{N-1}{2}}\pi^{-\frac{3N}{2}+1}\alpha^{\frac{N}{2}}|v|^{-2N}e^{-\frac{\alpha}{|v|^{2}}}(-1)^{N-1}\int
d{\phi_{\perp}}d{\sigma_{\perp}}d\bar{\psi}_{\perp}d\psi_{\perp}d\bar{\varphi}_{\perp}d\varphi_{\perp}\,e^{\tilde{K}^{\perp}_{F}+\tilde{K}^{\perp}_{B}+\tilde{V}_{F}+\tilde{V}_{B}+\tilde{V}_{BF}},$
(26)

where

$\displaystyle\begin{split}\tilde{K}^{\perp}_{F}&=-\bar{\varphi}_{\perp}\cdot\varphi_{\perp}-\bar{\psi}_{\perp}\cdot\varphi_{\perp}-\bar{\varphi}_{\perp}\cdot\psi_{\perp}-\epsilon\bar{\psi}_{\perp}\cdot\psi_{\perp},\\\
\tilde{K}^{\perp}_{B}&=-{\sigma_{\perp}}^{2}-2i{\sigma_{\perp}}\cdot{\phi_{\perp}}-\epsilon{\phi_{\perp}}^{2},\end{split}$
(27)

and $\tilde{V}_{F},\tilde{V}_{B},\tilde{V}_{BF}$ are given in (24) and (25).

##  4 Integrations over fermions

In the expression (26) the variables projected to the transverse directions,
i.e., ${\bar{\psi}}_{\perp}$, etc., just represent $N-1$-dimensional degrees
of freedom with no other restrictions. Therefore we can simply regard them as
$N-1$-dimensional variables from the beginning. The external variable, the
vector $v$, only appears in the coupling constants of the four-interactions
through $|v|$. Therefore we can simply write

$\displaystyle\rho(v)=\lim_{\epsilon\rightarrow+0}3^{\frac{N-1}{2}}\pi^{-\frac{3N}{2}+1}\alpha^{\frac{N}{2}}|v|^{-2N}e^{-\frac{\alpha}{|v|^{2}}}(-1)^{N-1}\int_{N^{*}}d\phi
d\sigma d{\bar{\psi}}d\psi
d\bar{\varphi}d{\varphi}\,e^{K_{F}+K_{B}+V_{F}+V_{B}+V_{BF}},$ (28)

where all the variables are $N-1$-dimensional, and

$\displaystyle\begin{split}K_{F}&=-\bar{\varphi}\cdot{\varphi}-{\bar{\psi}}\cdot{\varphi}-\bar{\varphi}\cdot\psi-\epsilon{\bar{\psi}}\cdot\psi,\\\
K_{B}&=-\sigma^{2}-2i\sigma\cdot\phi-\epsilon\phi^{2},\\\
V_{F}&=-\frac{|v|^{2}}{6\alpha}\left(({\bar{\psi}}\cdot{\varphi})^{2}+(\bar{\varphi}\cdot\psi)^{2}+2{\bar{\psi}}\cdot\bar{\varphi}\,{\varphi}\cdot\psi+2{\bar{\psi}}\cdot\psi\,\bar{\varphi}\cdot{\varphi}{}\right),\\\
V_{B}&=-\frac{2|v|^{2}}{3\alpha}\left(\sigma^{2}\phi^{2}+(\sigma\cdot\phi)^{2}\right),\\\
V_{BF}&=\frac{2i|v|^{2}}{3\alpha}\left({\bar{\psi}}\cdot\sigma\,{\varphi}\cdot\phi+\bar{\varphi}\cdot\sigma\,\psi\cdot\phi+{\bar{\psi}}\cdot\phi\,{\varphi}\cdot\sigma+\bar{\varphi}\cdot\phi\,\psi\cdot\sigma\right).\end{split}$
(29)

Note that, for simplicity, we are abusively using the same notations of the
variables as those in Section 3 with a different dimension, and, to indicate
this difference, the symbol $N^{*}$ is attached to the integral symbol in
(28).

To compute $\rho(v)$ explicitly, we first perform the fermionic integrations.
The integrand in (28) can be expanded in the fermions. A useful property in
this expansion is that the expansion in $V_{BF}$ stops at the fourth order:

$\displaystyle e^{V_{BF}}=\sum_{n=0}^{4}\frac{1}{n!}\left(V_{BF}\right)^{n}.$
(30)

This can be proven as follows. The fermions in $V_{BF}$ are projected to
$\phi$ or $\sigma$. Therefore $V_{BF}$ contains only eight independent
fermions in total, i.e.,
${\bar{\psi}}\cdot\phi,{\bar{\psi}}\cdot\sigma,\psi\cdot\phi,\psi\cdot\sigma,\bar{\varphi}\cdot\phi,\bar{\varphi}\cdot\sigma,{\varphi}\cdot\phi,{\varphi}\cdot\sigma$.
Since each term of $V_{BF}$ contains two of these fermions and products of
more than eight of these fermions vanish, we obtain

$\displaystyle(V_{BF})^{n}=0\hbox{ for }n>4.$ (31)

We also notice that, when $\sigma=\pm\phi$, only four of these fermions are
independent. Therefore,

$\displaystyle(V_{BF})^{n}=0\hbox{ for }n>2,\hbox{ when }\sigma=\pm\phi.$ (32)

Now we want to determine the functional forms of the fermionic integrals of
each summand in (30):

$\displaystyle\beta_{n}=\frac{(-1)^{N-1}}{n!}\int_{N^{*}}d{\bar{\psi}}d\psi
d\bar{\varphi}d{\varphi}\,(V_{BF})^{n}e^{K_{F}+V_{F}}\ \ \ (n\leq 4).$ (33)

From the form of $V_{BF}$ and the $O(N-1)$ symmetry, $\beta_{n}$ should be a
polynomial function of $\sigma^{2},\phi^{2},\sigma\cdot\phi$, and its order in
$\sigma,\phi$ should be $2n$. In addition, each term of the polynomial
function should contain equal numbers of $\sigma$ and $\phi$, and the
polynomial function should also be invariant under the interchange
$\sigma\leftrightarrow\phi$ as a whole.

From the above considerations, we uniquely obtain for $n=1$

$\displaystyle\beta_{1}=a_{1}\sigma\cdot\phi.$ (34)

Here $a_{1}$ is a function of $\epsilon$ and $|v|^{2}/\alpha$, and all $a_{i}$
below are also so.

For $n=2$, we have two possibilities,

$\displaystyle\beta_{2}=a_{2}(\sigma\cdot\phi)^{2}+a_{3}\sigma^{2}\phi^{2}.$
(35)

As for $n=3$, $\beta_{3}$ must vanish for $\sigma=\pm\phi$ due to (32).
Considering also the other conditions mentioned above, we uniquely obtain

$\displaystyle\beta_{3}$
$\displaystyle=a_{4}\left(\sigma^{2}\phi^{2}\sigma\cdot\phi-(\sigma\cdot\phi)^{3}\right).$
(36)

As for $n=4$, we need not only (32) but also the following property: Taking
the derivatives of $(V_{BF})^{4}$ with respect to $\sigma,\phi$ three times or
less must vanish, when $\sigma=\pm\phi$. This can be proven as follows.
Obviously, $(V_{BF})^{4}$ is proportional to the product of the eight
projected fermions, ${\bar{\psi}}\cdot\phi,{\bar{\psi}}\cdot\sigma,\ldots$.
After taking the derivatives of this three times or less with respect to
$\sigma,\phi$, each term contains at least five of the eight projected
fermions. Therefore all these terms vanish for $\sigma=\pm\phi$ because of the
same reason for (32). Now using this property, we can uniquely determine

$\displaystyle\beta_{4}=a_{5}\left(\sigma^{2}\phi^{2}-(\sigma\cdot\phi)^{2}\right)^{2}.$
(37)

Collecting all the results above, we conclude that the fermionic integrations
of the partition function has the following general form,

$\displaystyle\begin{split}&\int_{N^{*}}d{\bar{\psi}}d\psi
d\bar{\varphi}d{\varphi}\,e^{K_{F}+V_{F}+V_{BF}}\\\
&=a_{0}+a_{1}\sigma\cdot\phi+a_{2}(\sigma\cdot\phi)^{2}+a_{3}\sigma^{2}\phi^{2}+a_{4}\left(\sigma^{2}\phi^{2}\sigma\cdot\phi-(\sigma\cdot\phi)^{3}\right)+a_{5}\left(\sigma^{2}\phi^{2}-(\sigma\cdot\phi)^{2}\right)^{2}.\end{split}$
(38)

Here $a_{i}$ are generally functions of $\epsilon$ and $|v|^{2}/\alpha$, but
$\epsilon$ in them turn out to simply disappear in the $\epsilon\rightarrow+0$
limit in (28) without any essential roles. Therefore in the following sections
we ignore it by just putting $\epsilon=0$ in $a_{i}$.

##  5 Integrations over bosons

Assuming the form (38) of the fermionic integrations, we perform the bosonic
integrations but one in this section.

Let us first generalize the bosonic kinetic term by introducing new parameters
$\epsilon_{i}\ (i=1,2,3)$:

$\displaystyle
K^{\epsilon}_{B}=-\epsilon_{1}\sigma^{2}-2i\epsilon_{2}\,\sigma\cdot\phi-\epsilon_{3}\phi^{2},$
(39)

where $\epsilon_{1},\epsilon_{3}>0$ are assumed for the convergence of the
bosonic integration. The original kinetic term in (29) corresponds to
$\epsilon_{1}=\epsilon_{2}=1,\ \epsilon_{3}=\epsilon$.

Using this new kinetic term, (28), and (38), $\rho(v)$ can be expressed as

$\displaystyle\begin{split}&\rho(v)=3^{\frac{N-1}{2}}\pi^{-\frac{3N}{2}+1}\alpha^{\frac{N}{2}}|v|^{-2N}e^{-\frac{\alpha}{|v|^{2}}}\,G_{N},\end{split}$
(40)

where

$\displaystyle\begin{split}G_{N}=\left(a_{0}+a_{1}D_{2}+a_{2}D_{2}^{2}+a_{3}D_{1}+a_{4}(D_{1}D_{2}-D_{2}^{3})+a_{5}(D_{1}-D_{2}^{2})^{2}\right)\left.\int_{N^{*}}d\sigma
d\phi\,e^{K^{\epsilon}_{B}+V_{B}}\right|_{\epsilon_{1}=\epsilon_{2}=1\atop\epsilon_{3}=+0},\end{split}$
(41)

with $D_{1},D_{2}$ being the following partial derivative operators,

$\displaystyle\begin{split}D_{1}&=\frac{\partial^{2}}{\partial\epsilon_{1}\partial\epsilon_{3}},\\\
D_{2}&=-\frac{1}{2i}\frac{\partial}{\partial\epsilon_{2}}.\end{split}$ (42)

The bosonic integration in (41) does not seem to be fully integrable, but it
has a simpler expression. Since $\sigma$ appears at most quadratically in
$K_{B}^{\epsilon}+V_{B}$, the $\sigma$ integration can be performed:

$\displaystyle\int_{N^{*}}d\sigma\,e^{K^{\epsilon}_{B}+V_{B}}=\pi^{\frac{N-1}{2}}(\epsilon_{1}+8z\phi^{2})^{-\frac{1}{2}}(\epsilon_{1}+4z\phi^{2})^{\frac{N-2}{2}}\exp\left(-\frac{(\epsilon_{2}^{2}+\epsilon_{3}(\epsilon_{1}+8z\phi^{2}))\phi^{2}}{\epsilon_{1}+8z\phi^{2}}\right),$
(43)

where for brevity we have introduced

$\displaystyle\begin{split}z=\frac{|v|^{2}}{6\alpha}.\end{split}$ (44)

(43) does not depend on the angular directions of $\phi$, and therefore the
integration over these produces the spherical volume,
$2\pi^{(N-1)/2}|\phi|^{N-2}/\Gamma[(N-1)/2]$. Then, after applying the
derivative operators in (41) and performing a replacement of variable,
$|\phi|=\sqrt{x}/\sqrt{1-8zx}$, we obtain

$\displaystyle\begin{split}G_{N}=\frac{\pi^{N-1}}{4\,\Gamma\left[\frac{N-1}{2}\right]}\int_{0}^{\frac{1}{8z}}&dx\,e^{-x}x^{\frac{N-3}{2}}(1-4zx)^{-\frac{N+2}{2}}\\\
&\cdot(4a_{0}+2(-2ia_{1}+a_{2}+(N-1)a_{3})x\\\ &\ \ \
+(8iza_{1}-(3+4z)(a_{2}+a_{3})-i(N-2)a_{4})x^{2}+8z(a_{2}+a_{3})x^{3}),\end{split}$
(45)

where we have used

$\displaystyle
64z^{2}a_{0}-8iza_{1}+(4z-1)a_{2}+(-1-4z+8(N-1)z)a_{3}-i(N-2)a_{4}+N(N-2)a_{5}=0$
(46)

to delete $a_{5}$. As we will see in Section 6, (46) holds for all the cases
we consider (namely, $N\leq 8$). In fact this relation is essentially
important, because otherwise the integrand has an extra factor $1/(1-8zx)$,
and the possibility of the integration over $x$ which will be performed in
Section 6 would become unclear.

##  6 Analytic expressions

In this section, we explicitly compute the integration (45) for each case of
$N\leq 8$ to obtain exact analytic expressions. What is surprising is that the
seemingly difficult integrations can be done explicitly. The apparent reason
is that the integrand of (45) turns out to be a sum of a total derivative and
a simple integrable function.

###  6.1 $N=1$

This case is trivial, since the integration in (28) can just be ignored. We
obtain

$\displaystyle\rho_{N=1}(v)=\pi^{-\frac{1}{2}}\alpha^{\frac{1}{2}}|v|^{-2}e^{-\frac{\alpha}{|v|^{2}}}.$
(47)

Since $v=1/C$ for $N=1$ from (1), $\rho_{N=1}(v)\,dv$ is indeed equivalent to
$\alpha^{1/2}\pi^{-1/2}e^{-\alpha C^{2}}dC$, which is the Gaussian
distribution of $C$.

###  6.2 $N=2$

In this case, since the variables in (28) are all one-dimensional, we can
ignore the indices. Then we immediately obtain

$\displaystyle\begin{split}V^{N=2}_{F}&=-4z\,{\bar{\psi}}\psi\bar{\varphi}{\varphi},\\\
V^{N=2}_{B}&=-8z\,\sigma^{2}\phi^{2},\\\
V^{N=2}_{BF}&=8iz\,\sigma\phi\left({\bar{\psi}}{\varphi}+\bar{\varphi}\psi\right).\end{split}$
(48)

Hence, by explicitly expanding the integrand, we obtain,

$\displaystyle\int_{N^{*}}d{\bar{\psi}}d\psi
d\bar{\varphi}d{\varphi}\,e^{K_{F}^{N=2}+V_{F}^{N=2}+V_{BF}^{N=2}}=1+4z-16iz\sigma\phi-64z^{2}\sigma^{2}\phi^{2}.$
(49)

This determines

$\displaystyle\begin{split}&a_{0}=1+4z,\\\ &a_{1}=-16iz,\\\
&a_{2}=-64z^{2},\\\ &\hbox{Others}=0,\end{split}$ (50)

by comparing with (38). This indeed satisfies (46). Putting (50) into (45), we
obtain

$\displaystyle\begin{split}G_{N=2}&=\sqrt{\pi}\int_{0}^{\frac{1}{8z}}dx\,e^{-x}x^{-\frac{1}{2}}(1+4z-8zx)\\\
&=\sqrt{\pi}\left((1+4z)\gamma\left[\frac{1}{2},\frac{1}{8z}\right]-8z\gamma\left[\frac{3}{2},\frac{1}{8z}\right]\right)\\\
&=\sqrt{\pi}\left(\gamma\left[\frac{1}{2},\frac{1}{8z}\right]+\sqrt{8z}\,e^{-\frac{1}{8z}}\right),\end{split}$
(51)

where the lower incomplete gamma function $\gamma[\cdot,\cdot]$ is defined by

$\displaystyle\gamma\left[a,y\right]=\int_{0}^{y}dt\,t^{a-1}\,e^{-t},$ (52)

and we have used its property,

$\displaystyle\gamma\left[a+1,y\right]=a\gamma\left[a,y\right]-y^{a}e^{-y}.$
(53)

The lower incomplete gamma function with index $1/2$ in (51) is related to the
error function by

$\displaystyle\gamma\left[\frac{1}{2},y\right]=\sqrt{\pi}\,\hbox{erf}\left[\sqrt{y}\right].$
(54)

Therefore, from (44) and (51), $G_{N=2}$ is represented by a combination of
polynomial, exponential, and error functions of $|v|$. This is common to the
other cases of $N$ shown below.

###  6.3 $N=3$

The fermionic integration of (38) seems too complicated to perform by hand.
Rather, we use a Mathematica package for Grassmann variables [22]. The result
is

$\displaystyle\begin{split}a_{0}&=1+4z+28z^{2},\\\ a_{1}&=-16i(z+2z^{2}),\\\
a_{2}&=-32(3z^{2}+2z^{3}),\\\ a_{3}&=-32(-z^{2}+6z^{3}),\\\
a_{4}&=-256iz^{3},\\\ a_{5}&=256z^{4}.\end{split}$ (55)

This indeed satisfies (46). Then, by putting this into (45), we obtain

$\displaystyle\begin{split}G_{N=3}=\pi^{2}\int_{0}^{\frac{1}{8z}}dx\,\frac{e^{-x}}{(1-4zx)^{\frac{5}{2}}}\big{(}&1+4z+28z^{2}-16z(1+3z+14z^{2})x\\\
&+16z^{2}(5+16z+16z^{2})x^{2}-128z^{3}(1+4z)x^{3}\big{)}.\end{split}$ (56)

A surprising fact is that the integrand in (56) is actually a total
derivative, and we therefore obtain

$\displaystyle\begin{split}G_{N=3}&=-\pi^{2}\int_{0}^{\frac{1}{8z}}dx\,\frac{d}{dx}\left(\frac{e^{-x}}{(1-4zx)^{\frac{3}{2}}}\left(1-2z-4z(3+4z)x+32z^{2}(1+4z)x^{2}\right)\right)\\\
&=\pi^{2}\left(1-2z+4\sqrt{2}ze^{-\frac{1}{8z}}\right).\end{split}$ (57)

###  6.4 Larger odd $N$

The strategy taken in Section 6.3 can be generalized for larger odd $N$ in the
following manner. We first compute the $a_{i}\,(i=0,1,\ldots,5)$, which are
listed for $N=5,7$ in Appendix B, by using the aforementioned Mathematica
package. They indeed satisfy (46). Then, by putting them into (45), we obtain
$G_{N}$. Similarly to the case of $N=3$ in 6.3, what we find is that the
integrand of $G_{N}$ is a total derivative of the following form for $N=5,7$:

$\displaystyle G_{N:{\rm
odd}}=\frac{\pi^{N-1}}{4\,\Gamma\left[\frac{N-1}{2}\right]}\int_{0}^{\frac{1}{8z}}$
$\displaystyle
dx\,\frac{d}{dx}\frac{e^{-x}\sum_{n=0}^{\frac{N+1}{2}}b_{n}x^{n}}{(1-4zx)^{\frac{N}{2}}},$
(58)

where $b_{n}\,(n=0,1,\ldots,(N+1)/2)$ are some polynomial functions of $z$.
The explicit forms of $b_{n}$ are given in Appendix B. Therefore, we obtain

$\displaystyle G_{N:{\rm
odd}}=\frac{\pi^{N-1}}{4\,\Gamma\left[\frac{N-1}{2}\right]}\left(2^{\frac{N}{2}}e^{-\frac{1}{8z}}\sum_{n=0}^{\frac{N+1}{2}}\frac{b_{n}}{(8z)^{n}}-b_{0}\right).$
(59)

By using the $b_{i}$ in Appendix B, the explicit expressions of $G_{N}$ for
$N=5,7$ are given by

$\displaystyle\begin{split}G_{N=5}&=\pi^{4}\left(1-12z+12z^{2}+\sqrt{2}e^{-\frac{1}{8z}}(1+12z+12z^{2})\right),\\\
G_{N=7}&=\pi^{6}\left(1-30z+180z^{2}-120z^{3}+\frac{\sqrt{2}e^{-\frac{1}{8z}}}{8z}(1+8z+120z^{2}-480z^{3}+2640z^{4})\right).\end{split}$
(60)

Let us lastly comment on our equipment. The computations were done on a
machine which had a Xeon W2295 (3.0GHz, 18 cores), 128GB DDR4 memory, and
Ubuntu 20 as OS. The computation of $a_{i}$ quickly takes longer time as $N$
becomes larger. $G_{N=8}$, which appears in the next subsection, was the
largest feasible case, while we failed to obtain $G_{N=9}$ seemingly because
of a memory shortage.

###  6.5 Larger even $N$

The difference from the odd case of Section 6.4 is that the integrand of
$G_{N}$ is a sum of a total derivative and a simple integrable term:

$\displaystyle
G_{N}=\frac{\pi^{N-1}}{4\,\Gamma\left[\frac{N-1}{2}\right]}\int_{0}^{\frac{1}{8z}}$
$\displaystyle
dx\,\left(c_{0}\,x^{-\frac{1}{2}}e^{-x}+\frac{d}{dx}\frac{x^{\frac{1}{2}}e^{-x}\sum_{n=0}^{\frac{N}{2}}b_{n}x^{n}}{(1-4zx)^{\frac{N}{2}}}\right).$
(61)

Then, by doing the integration, we obtain

$\displaystyle
G_{N}=\frac{\pi^{N-1}}{4\,\Gamma\left[\frac{N-1}{2}\right]}\left(c_{0}\gamma\left[\frac{1}{2},\frac{1}{8z}\right]+2^{\frac{N}{2}}e^{-\frac{1}{8z}}\sum_{n=0}^{\frac{N}{2}}\frac{b_{n}}{(8z)^{n+\frac{1}{2}}}\right).$
(62)

The lists of $a_{i},b_{i},c_{0}$ for $N=4,6,8$ are given in Appendix C.
$a_{i}$ indeed satisfy (46). By putting the values of $b_{i},c_{0}$ we obtain
the explicit forms of $G_{N}$ as

$\displaystyle\begin{split}G_{N=4}&=\pi^{\frac{5}{2}}\left(6\sqrt{2}e^{-\frac{1}{8z}}\sqrt{z}(1+2z)+(1-6z)\,\gamma\left[\frac{1}{2},\frac{1}{8z}\right]\right),\\\
G_{N=6}&=\pi^{\frac{9}{2}}\left(\frac{2\sqrt{2}e^{-\frac{1}{8z}}(1+15z+180z^{3})}{3\sqrt{z}}+(1-20z+60z^{2})\,\gamma\left[\frac{1}{2},\frac{1}{8z}\right]\right),\\\
G_{N=8}&=\pi^{\frac{13}{2}}\Bigg{(}\frac{\sqrt{2}e^{-\frac{1}{8z}}(1+210z^{2}-2100z^{3}+12600z^{4}+25200z^{5})}{15z^{\frac{3}{2}}}\\\
&\hskip
113.81102pt+(1-42z+420z^{2}-840z^{3})\,\gamma\left[\frac{1}{2},\frac{1}{8z}\right]\Bigg{)}.\end{split}$
(63)

##  7 Comparison with Monte Carlo simulations

In this section, we compare the results in Section 6 with Monte Carlo
simulations. The procedure is basically the same as that used in [16, 17]. To
make this paper self-contained, however, we review the method below.

The eigenvector equation (1) is a system of polynomial equations and it can be
solved by an appropriate polynomial equation solver, unless $N$ is too large.
We use Mathematica 13 for this purpose. It gives generally complex solutions
to the equation (1), and we pick up only real ones, since we are only counting
real eigenvectors (or Z-eigenvalues). Whether this method covers all the real
solutions or not can be checked by whether the number of generally complex
solutions obtained for each $C$ by a polynomial equation solver agrees with
the known number $2^{N}-1$ [20]444In fact, for large $N$, Mathematica 13 seems
to miss a few solutions for some $C$. We have not pursued the reason for that,
but the missing portion is $\sim 10^{-4}$ even for our largest case of $N=16$,
and is statistically irrelevant in the present study..

With the above method of solving the eigenvector equation (1), our procedure
of Monte Carlo simulation is given as follows.

  * •

Randomly generate a real symmetric tensor $C$ with components,
$C_{ijk}=\sigma/\sqrt{d(i,j,k)}\ (1\leq i\leq j\leq k\leq N)$, where $\sigma$
has a normal distribution with mean value zero and standard deviation one, and
$d(i,j,k)$ is a degeneracy factor,

$\displaystyle\begin{split}d(i,j,k)=\left\\{\begin{array}[]{ll}1,&i=j=k\\\
3,&i=j\neq k,\,i\neq j=k,\,k=i\neq j\\\ 6,&i\neq k\neq j\neq
i\end{array}\right..\end{split}$ (64)

This random generation of $C$ corresponds to $\alpha=1/2$ in (4), since

$\displaystyle\begin{split}C^{2}=C_{abc}C_{abc}=\sum_{i\leq j\leq
k=1}^{N}d(i,j,k)\,C_{ijk}^{2}\end{split}$ (65)

due to $C$ being a symmetric tensor.

  * •

Compute the real eigenvectors of a generated $C$ by the aforementioned method.

  * •

Store each size $|v|$ of all the real eigenvectors.

  * •

Repeat the above processes.

By this procedure, we obtain a sequence of sizes of real eigenvectors,
$|v|_{i}\ (i=1,2,\ldots,L)$. Then the size distribution of real eigenvectors
can be obtained by

$\displaystyle\rho^{\rm MC}_{\rm size}((k+1/2)\delta v)=\frac{1}{\delta
vN_{C}}\sum_{i=1}^{L}\theta(k\delta v<|v|_{i}\leq(k+1)\delta v),$ (66)

where $\delta v$ is a bin size, $N_{C}$ is the total number of randomly
generated $C$, $k=0,1,2,\ldots$, and $\theta(\cdot)$ is a support function
which takes 1 if the inequality of the argument is satisfied, but zero
otherwise. For the comparison with the analytical results obtained in Section
6, (66) should be compared with the size distribution

$\displaystyle\rho_{\rm size}(|v|)=\rho(v)S_{N-1}|v|^{N-1},$ (67)

where $S_{N-1}=2\pi^{N/2}/\Gamma[N/2]$ denotes the surface volume of a unit
sphere in an $N$-dimensional space, and the vector $v$ in the argument of
$\rho(\cdot)$ is an arbitrary vector of size $|v|$ due to the rotational
symmetry.

Since a Z-eigenvalue is related with the size of an eigenvector by the
relation (6), the real (or Z-) eigenvalue distribution is given by

$\displaystyle\rho^{\rm MC}_{\rm eig}((k+1/2)\delta\zeta)=\frac{1}{\delta\zeta
N_{C}}\sum_{i=1}^{L}\theta(k\delta\zeta<1/|v|_{i}\leq(k+1)\delta\zeta),$ (68)

where $\delta\zeta$ is a bin size. This quantity should be compared with (7).

Figure 1: The results of the Monte Carlo simulation for $N=8,N_{C}=10000$ are
compared with the analytic expressions. Left: The eigenvector size
distribution. The Monte Carlo result (66) with $\delta v=0.03$ (dots) is
compared with (67) using (40) and $G_{N=8}$ in (63) (solid line). Right: The
eigenvalue distribution. The Monte Carlo result (68) with $\delta\zeta=0.03$
(dots) is compared with the analytic expression through (7).

In Figure 1, the Monte Carlo results are compared with the analytic
expressions for $N=8$. They agree precisely. Similar precise agreement has
been obtained for all $N\leq 8$.

##  8 Extrapolation to general $N$

In this section, we point out a few patterns which exist in the expressions of
the distributions derived for small $N$ in Section 6, and guess an
extrapolation to general $N$. A motivation for doing this is to improve the
large-$N$ expression previously obtained by an approximation using a
Schwinger-Dyson equation [17]. The issue of the previous result was that,
while the functional form agreed well with the numerical simulation for
large-$N$, the overall factor did not. In this section, using the
extrapolation, we will guess the overall factor for general $N$ and find good
agreement with Monte Carlo simulations.

After a thought one notices that $G_{N}$ for even $N$ (namely, $N=2,4,6,8$) in
Section 6 can be expressed by the following general form:

$\displaystyle\begin{split}G_{N:{\rm
even}}=&\pi^{N-\frac{3}{2}}z^{\frac{N-1}{2}}H_{N-1}\left[\frac{1}{2\sqrt{z}}\right]\gamma\left[\frac{1}{2},\frac{1}{8z}\right]\\\
&+\sqrt{2}\pi^{N-\frac{3}{2}}\frac{N!}{\frac{N}{2}!}z^{\frac{N-1}{2}}e^{-\frac{1}{8z}}\left(1+\frac{d_{1}}{z}+\frac{d_{2}}{z^{2}}+\cdots+\frac{d_{N-3}}{z^{N-3}}\right),\end{split}$
(69)

where $H_{n}[\cdot]$ are Hermite polynomials, $d_{i}$ are some coefficients
generally depending on $N$, and specifically

$\displaystyle d_{1}=\frac{1+(-1)^{\frac{N}{2}}}{4}.$ (70)

As for $d_{i}\ (i\geq 2)$, we could not find reasonably simple functions of
$N$.

Let us discuss the real eigenvalue distribution for large-$N$, assuming (69)
with (70). Because of the relations (6) and (44), and the fact that the major
part of the distribution is around $\zeta\sim 0$ as in Figure 1, we are
interested in a $1/z$ expansion of (69) with (70). By explicitly doing this
one obtains

$\displaystyle
G_{N}\sim\sqrt{2}\pi^{N-\frac{3}{2}}\frac{\Gamma[N+1]}{\Gamma\left[\frac{N}{2}+1\right]}z^{\frac{N-1}{2}}\left(1+\frac{1}{8z}+\cdots\right),$
(71)

where the factorials are replaced by Gamma functions to also be applicable for
odd $N$ below. By combining with the previous result in [17] that the
large-$N$ eigenvalue distribution is given by a Gaussian function of $|v|$,
one could assume that the expansion in $1/z$ of (71) comes from the
expression,

$\displaystyle
G_{N}\sim\sqrt{2}\pi^{N-\frac{3}{2}}\frac{\Gamma[N+1]}{\Gamma\left[\frac{N}{2}+1\right]}z^{\frac{N-1}{2}}e^{\frac{1}{8z}}.$
(72)

Putting this into (7) using (6), (40) and (44), we obtain

$\displaystyle\rho_{\rm eig}(\zeta)\sim
2^{-\frac{N}{2}+2}\alpha^{\frac{1}{2}}\pi^{-\frac{1}{2}}\frac{\Gamma[N+1]}{\Gamma\left[\frac{N}{2}+1\right]\Gamma\left[\frac{N}{2}\right]}e^{-\frac{\alpha}{4}\zeta^{2}},$
(73)

which indeed is a Gaussian distribution. While the coefficient $\alpha/4$ in
the exponent indeed agrees with the previous result in [17], the overall
factor is different.

Figure 2:  The guessed large-$N$ formula (73) (solid line) is compared with
the Monte Carlo simulation (68) (dots). Left:
$N=14,\delta\zeta=0.03,N_{C}=400$. Middle: $N=15,\delta\zeta=0.06,N_{C}=100$.
Right: $N=16,\delta\zeta=0.1,N_{C}=72$.

In Figure 2 we compare the large-$N$ expression (73) with the Monte Carlo
simulations for $N=14,15,16$.555With our machine power (See the last part of
Section 6.4), it took one full day for the computation of the $N=16$ case.
They agree well, and the agreement seems to become slightly better, as $N$
becomes larger.

In the above discussion, we have not used the odd $N$ cases. The reason is
that we could not find simple expressions valid across all the odd $N$ cases
obtained in Section 6. We merely notice

$\displaystyle G_{N:{\rm
odd}}=\pi^{N-1}z^{\frac{N-1}{2}}H_{N-1}\left(\frac{1}{2\sqrt{z}}\right)+d_{0}^{\prime}\,e^{-\frac{1}{8z}}\left(1+\frac{d^{\prime}_{1}}{z}+\cdots+\frac{d^{\prime}_{N-3}}{z^{N-3}}\right),$
(74)

where we do not have any simple expressions for $d^{\prime}_{i}$.

Figure 3:  The coefficients $e_{0},e_{1}$ of the expansion (75) are plotted
against $N$. $\bar{e}_{0}=\sqrt{2}\pi^{N-3/2}\Gamma[N]/\Gamma[N/2+1]$ (see
(71)). The values of $e_{0},e_{1}$ for odd $N$ seem to approach those for even
$N$, as $N$ becomes larger.

Rather, we find that the large-$z$ behavior of $G_{N:{\rm odd}}$ approaches
that of $G_{N:{\rm even}}$, as $N$ becomes larger. To see this, let us
introduce the parameters of the $1/z$ expansion as

$\displaystyle G_{N}\sim
e_{0}\,z^{\frac{N-1}{2}}\left(1+\frac{e_{1}}{z}+\cdots\right).$ (75)

As in (71), for even $N$, we have

$\displaystyle
e_{0}=\frac{\sqrt{2}\pi^{N-\frac{3}{2}}\Gamma[N+1]}{\Gamma\left[\frac{N}{2}+1\right]},\
e_{1}=\frac{1}{8}.$ (76)

In Figure 3, we plot $e_{0},e_{1}$ for each $N$, which are extracted by
performing the $1/z$ expansions of the explicit expressions in Section 6.
Indeed $e_{0},e_{1}$ for odd $N$ seem to approach those for even $N$, as $N$
becomes larger. This implies that the large-$N$ formula (73) can commonly be
used for both even and odd $N$.

##  9 Summary and future prospects

In this paper we have explicitly computed the real eigenvalue/vector
distributions for real symmetric order-three tensors with Gaussian
distributions for $N\leq 8$. This has been achieved by rewriting the problem
as the computation of a partition function of a zero-dimensional boson-fermion
system with four-interactions. We have extrapolated the exact expressions for
$N\leq 8$ to guess a large-$N$ expression. Monte Carlo simulations have been
performed, and precise agreement and good agreement have been obtained for the
exact small-$N$ expressions and the large-$N$ expression, respectively.

It seems surprising that the complicated integrations of the partition
function can be performed exactly with the final simple results expressed by
polynomial, exponential, and error functions. It is suspected that this is
true for any $N$, but deeper understanding of this aspect needed for such
generalization is left for future study. Considering the final simple forms,
it would be possible that there exist much easier ways to compute the
distributions than what has been done in this paper.

There would exist various directions to pursue by a similar strategy of this
paper. One obvious direction is the generalization of tensors to those with
different orders and non-Gaussian distributions. It would also be interesting
to pursue the possibility of topological changes of the distributions in
tensor models as in matrix models. In fact, some tensor systems seemingly
suggest such possibilities [23, 24]. It would also be possible to compute
correlations among eigenvalues/vectors.

## Acknowledgements

The author is supported in part by JSPS KAKENHI Grant No.19K03825.

##  Appendix Appendix A Computations of $E_{2}$ and $E_{3}$

From (17),

$\displaystyle E_{2}=-\frac{4}{\alpha}A_{1},$ (77)

with

$\displaystyle
A_{1}=\left(\frac{1}{6}\sum_{s}v_{s_{a}}\sigma_{s_{b}}\phi_{s_{c}}\right)^{2}.$
(78)

By explicitly writing down the sum and noting that one of the symmetrization
is not necessary, we obtain

$\displaystyle\begin{split}A_{1}&=\frac{1}{6}v_{a}\sigma_{b}\phi_{c}(v_{a}\sigma_{b}\phi_{c}+v_{b}\sigma_{c}\phi_{a}+v_{c}\sigma_{a}\phi_{b}+v_{b}\sigma_{a}\phi_{c}+v_{a}\sigma_{c}\phi_{b}+v_{c}\sigma_{b}\phi_{a})\\\
&=\frac{1}{6}\left(v^{2}\sigma^{2}\phi^{2}+v\cdot\sigma\,\sigma\cdot\phi\,\phi\cdot
v+v\cdot\phi\,\sigma\cdot v\,\phi\cdot\sigma+v\cdot\sigma\,\sigma\cdot
v\,\phi^{2}+v^{2}\,\sigma\cdot\phi\,\phi\cdot\sigma+v\cdot\phi\,\sigma^{2}\,\phi\cdot
v\right)\\\
&=\frac{|v|^{2}}{6}\left(\sigma^{2}\phi^{2}+2{\phi_{\parallel}}{\sigma_{\parallel}}\,\sigma\cdot\phi+{\sigma_{\parallel}}^{2}\phi^{2}+{\phi_{\parallel}}^{2}\sigma^{2}+(\sigma\cdot\phi)^{2}\right)\\\
&=\frac{|v|^{2}}{6}\left({\sigma_{\perp}}^{2}{\phi_{\perp}}^{2}+({\phi_{\perp}}\cdot{\sigma_{\perp}})^{2}+2{\sigma_{\parallel}}^{2}{\phi_{\perp}}^{2}+2{\sigma_{\perp}}^{2}{\phi_{\parallel}}^{2}+4{\sigma_{\parallel}}{\phi_{\parallel}}\,{\phi_{\perp}}\cdot{\sigma_{\perp}}+6{\sigma_{\parallel}}^{2}{\phi_{\parallel}}^{2}\right),\end{split}$
(79)

where we have used $v\cdot\phi=|v|{\phi_{\parallel}}$, etc., and
$\phi\cdot\sigma={\phi_{\parallel}}{\sigma_{\parallel}}+{\phi_{\perp}}\cdot{\sigma_{\perp}}$,
etc..

Let us next compute

$\displaystyle E_{3}=\frac{4i}{\alpha}A_{2},$ (80)

where

$\displaystyle A_{2}$
$\displaystyle=(v_{a}{\bar{\psi}}_{b}{\varphi}_{c}+v_{a}\bar{\varphi}_{b}\psi_{c})\,\frac{1}{6}\sum_{s}v_{s_{a}}\sigma_{s_{b}}\phi_{s_{c}},$
(81)

from (17). By explicitly expanding the sum, we obtain

$\displaystyle\begin{split}A_{2}&=v_{a}{\bar{\psi}}_{b}{\varphi}_{c}(v_{a}\sigma_{b}\phi_{c}+v_{b}\sigma_{c}\phi_{a}+v_{c}\sigma_{a}\phi_{b}+v_{b}\sigma_{a}\phi_{c}+v_{a}\sigma_{c}\phi_{b}+v_{c}\sigma_{b}\phi_{a})+(\psi\leftrightarrow{\varphi})\\\
&=\frac{|v|^{2}}{6}\big{(}\bar{\psi}_{\perp}\cdot{\sigma_{\perp}}\,\varphi_{\perp}\cdot{\phi_{\perp}}+\bar{\psi}_{\perp}\cdot{\phi_{\perp}}\,\varphi_{\perp}\cdot{\sigma_{\perp}}+\bar{\varphi}_{\perp}\cdot{\sigma_{\perp}}\,\psi_{\perp}\cdot{\phi_{\perp}}+\bar{\varphi}_{\perp}\cdot{\phi_{\perp}}\,\psi_{\perp}\cdot{\sigma_{\perp}}\\\
&\ \ \ \ \ \ \ \ \
+2{\bar{\psi}}_{\parallel}{\sigma_{\parallel}}\varphi_{\perp}\cdot{\phi_{\perp}}+2{\phi_{\parallel}}{\bar{\psi}}_{\parallel}\varphi_{\perp}\cdot{\sigma_{\perp}}+2{\varphi}_{\parallel}{\phi_{\parallel}}\bar{\psi}_{\perp}\cdot{\sigma_{\perp}}+2{\varphi}_{\parallel}{\sigma_{\parallel}}\bar{\psi}_{\perp}\cdot{\phi_{\perp}}\\\
&\ \ \ \ \ \ \ \ \
+2\bar{\varphi}_{\parallel}{\sigma_{\parallel}}\psi_{\perp}\cdot{\phi_{\perp}}+2{\phi_{\parallel}}\bar{\varphi}_{\parallel}\psi_{\perp}\cdot{\sigma_{\perp}}+2\psi_{\parallel}{\phi_{\parallel}}\bar{\varphi}_{\perp}\cdot{\sigma_{\perp}}+2\psi_{\parallel}{\sigma_{\parallel}}\bar{\varphi}_{\perp}\cdot{\phi_{\perp}}\\\
&\ \ \ \ \ \ \ \ \
+6{\bar{\psi}}_{\parallel}{\varphi}_{\parallel}{\phi_{\parallel}}{\sigma_{\parallel}}+6\bar{\varphi}_{\parallel}\psi_{\parallel}{\phi_{\parallel}}{\sigma_{\parallel}}\big{)},\end{split}$
(82)

where we have taken similar computational steps as those for $A_{1}$.

##  Appendix Appendix B Lists of $a_{i}$ and $b_{i}$ for $N=5,7$

For $N=5$,

$\displaystyle\begin{split}a_{0}&=1-8z+120z^{2}+480z^{3}+2640z^{4},\\\
a_{1}&=-16i(z-6z^{2}+60z^{3}+120z^{4}),\\\
a_{2}&=-32(3z^{2}-10z^{3}+100z^{4}+200z^{5}),\\\
a_{3}&=-32(-z^{2}+10z^{3}+20z^{4}+280z^{5}),\\\
a_{4}&=-256i(z^{3}+20z^{5}),\\\ a_{5}&=256(z^{4}+4z^{5}+28z^{6}),\end{split}$
(83)

and

$\displaystyle\begin{split}b_{0}&=-4(1-12z+12z^{2}),\\\
b_{1}&=4(-1+22z-132z^{2}+120z^{3}),\\\ b_{2}&=48(z-8z^{2}+60z^{3}+80z^{4}),\\\
b_{3}&=-128(z^{2}+60z^{4}+240z^{5}).\end{split}$ (84)

For $N=7$,

$\displaystyle\begin{split}a_{0}&=1-36z+660z^{2}-3360z^{3}+25200z^{4}+100800z^{5}+504000z^{6},\\\
a_{1}&=-16i(z-30z^{2}+440z^{3}-1680z^{4}+8400z^{5}+16800z^{6}),\\\
a_{2}&=-32(3z^{2}-70z^{3}+840z^{4}-1680z^{5}+11760z^{6}+30240z^{7}),\\\
a_{3}&=-32(-z^{2}+30z^{3}-280z^{4}+1680z^{5}+5040z^{6}+36960z^{7}),\\\
a_{4}&=-256i(z^{3}-16z^{4}+168z^{5}+1680z^{7}),\\\
a_{5}&=256(z^{4}-8z^{5}+120z^{6}+480z^{7}+2640z^{8}),\end{split}$ (85)

and

$\displaystyle\begin{split}b_{0}&=8(-1+30z-180z^{2}+120z^{3}),\\\
b_{1}&=-8(1-44z+600z^{2}-2640z^{3}+1680z^{4}),\\\
b_{2}&=4(-1+58z-1160z^{2}+9360z^{3}-28560z^{4}+16800z^{5}),\\\
b_{3}&=16(3z-100z^{2}+1480z^{3}-6720z^{4}+25200z^{5}+33600z^{6}),\\\
b_{4}&=-128(z^{2}-20z^{3}+280z^{4}+8400z^{6}+33600z^{7}).\end{split}$ (86)

##  Appendix Appendix C Lists of $a_{i},b_{i},c_{0}$ for $N=4,6,8$

For $N=4$,

$\displaystyle\begin{split}a_{0}&=1+60z^{2}+240z^{3},\\\
a_{1}&=-16i(z+20z^{3}),\\\ a_{2}&=-32(3z^{2}+2z^{3}+24z^{4}),\\\
a_{3}&=-32(-z^{2}+6z^{3}+32z^{4}),\\\ a_{4}&=-256i(z^{3}+2z^{4}),\\\
a_{5}&=256(z^{4}+4z^{5}),\end{split}$ (87)

and

$\displaystyle\begin{split}c_{0}&=-2(-1+6z),\\\ b_{0}&=4(-1+6z),\\\
b_{1}&=16(3z-2z^{2}+40z^{3}),\\\
b_{2}&=-128(z^{2}+4z^{3}+28z^{4}).\end{split}$ (88)

For $N=6$,

$\displaystyle\begin{split}a_{0}&=1-20z+280z^{2}+8400z^{4}+33600z^{5},\\\
a_{1}&=-16i(z-16z^{2}+168z^{3}+1680z^{5}),\\\
a_{2}&=-32(3z^{2}-34z^{3}+300z^{4}+360z^{5}+2400z^{6}),\\\
a_{3}&=-32(-z^{2}+18z^{3}-60z^{4}+600z^{5}+2880z^{6}),\\\
a_{4}&=-256i(z^{3}-6z^{4}+60z^{5}+120z^{6}),\\\
a_{5}&=256(z^{4}+60z^{6}+240z^{7}),\end{split}$ (89)

and

$\displaystyle\begin{split}c_{0}&=3(1-20z+60z^{2}),\\\
b_{0}&=-6(1-20z+60z^{2}),\\\ b_{1}&=4(-1+18z)(1-20z+60z^{2}),\\\
b_{2}&=16(3z-56z^{2}+540z^{3}-240z^{4}+3360z^{5}),\\\
b_{3}&=-128(z^{2}-8z^{3}+120z^{4}+480z^{5}+2640z^{6}).\end{split}$ (90)

For $N=8$,

$\displaystyle\begin{split}a_{0}&=1-56z+1428z^{2}-15120z^{3}+105840z^{4}+2116800z^{6}+8467200z^{7},\\\
a_{1}&=-16i(z-48z^{2}+1020z^{3}-8640z^{4}+45360z^{5}+302400z^{7}),\\\
a_{2}&=-32(3z^{2}-118z^{3}+2080z^{4}-12880z^{5}+58800z^{6}+84000z^{7}+470400z^{8}),\\\
a_{3}&=-32(-z^{2}+46z^{3}-760z^{4}+6160z^{5}-8400z^{6}+117600z^{7}+537600z^{8}),\\\
a_{4}&=-256i(z^{3}-30z^{4}+440z^{5}-1680z^{6}+8400z^{7}+16800z^{8}),\\\
a_{5}&=256(z^{4}-20z^{5}+280z^{6}+8400z^{8}+33600z^{9}),\end{split}$ (91)

and

$\displaystyle\begin{split}c_{0}&=-(15/2)(-1+42z-420z^{2}+840z^{3}),\\\
b_{0}&=15(-1+42z-420z^{2}+840z^{3}),\\\
b_{1}&=-10(-1+24z)(-1+42z-420z^{2}+840z^{3}),\\\
b_{2}&=4(1-40z+360z^{2})(-1+42z-420z^{2}+840z^{3}),\\\
b_{3}&=48(z-52z^{2}+1140z^{3}-10360z^{4}+50400z^{5}-16800z^{6}+201600z^{7}),\\\
b_{4}&=-128(z^{2}-36z^{3}+660z^{4}-3360z^{5}+25200z^{6}+100800z^{7}+504000z^{8}).\end{split}$
(92)

## References

  * [1] E. Wigner, “On the Distribution of the Roots of Certain Symmetric Matrices,” Annals of Mathematics 67 (2), 325-327 (1958) https://doi.org/10.2307/1970008.
  * [2] E. Brezin, C. Itzykson, G. Parisi and J. B. Zuber, “Planar Diagrams,” Commun. Math. Phys. 59, 35 (1978) doi:10.1007/BF01614153
  * [3] B. Eynard, Counting surfaces : CRM Aisenstadt chair lectures, Prog. Math. Phys. 70, Birkhäuser, 2016.
  * [4] 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
  * [5] 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
  * [6] 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
  * [7] N. Sasakura, “Tensor model for gravity and orientability of manifold,” Mod. Phys. Lett. A 6, 2613-2624 (1991) doi:10.1142/S0217732391003055
  * [8] 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
  * [9] 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]].
  * [10] M. Ouerfelli, V. Rivasseau and M. Tamaazousti, “The Tensor Track VII: From Quantum Gravity to Artificial Intelligence,” [arXiv:2205.10326 [hep-th]].
  * [11] L. Qi, H. Chen, Y, Chen, Tensor Eigenvalues and Their Applications, Springer, Singapore, 2018.
  * [12] P. Breiding, “The expected number of eigenvalues of a real gaussian tensor,” SIAM Journal on Applied Algebra and Geometry 1, 254-271 (2017).
  * [13] P. Breiding, “How many eigenvalues of a random symmetric tensor are real?,” Transactions of the American Mathematical Society 372, 7857-7887 (2019).
  * [14] O. Evnin, “Melonic dominance and the largest eigenvalue of a large random tensor,” Lett. Math. Phys. 111, 66 (2021) doi:10.1007/s11005-021-01407-z [arXiv:2003.11220 [math-ph]].
  * [15] R. Gurau, “On the generalization of the Wigner semicircle law to real symmetric tensors,” [arXiv:2004.02660 [math-ph]].
  * [16] N. Sasakura, “Signed distributions of real tensor eigenvectors of Gaussian tensor model via a four-fermi theory,” [arXiv:2208.08837 [hep-th]].
  * [17] N. Sasakura, “Real tensor eigenvalue/vector distributions of the Gaussian tensor model via a four-fermi theory,” [arXiv:2209.07032 [hep-th]].
  * [18] L. Qi, “Eigenvalues of a real supersymmetric tensor,” Journal of Symbolic Computation 40, 1302-1324 (2005).
  * [19] L.H. Lim, “Singular Values and Eigenvalues of Tensors: A Variational Approach,” in Proceedings of the IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP ’05), Vol. 1 (2005), pp. 129–132.
  * [20] D. Cartwright and B. Sturmfels, “The number of eigenvalues of a tensor,” Linear algebra and its applications 438, 942-952 (2013).
  * [21] J. Zinn-Justin, Quantum Field Theory and Critical Phenomena, (Clarendon Press, Oxford, 1989).
  * [22] M. Headrick, J Michelson, grassmann.m, https://people.brandeis.edu/~headrick/Mathematica/.
  * [23] 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]].
  * [24] N. Sasakura, “Splitting-merging transitions in a tensor-vectors system in exact large-$N$ limits,” [arXiv:2206.12017 [hep-th]].


# Open-orbit induced low field extremely large magnetoresistance in
graphene/h-BN superlattices

Zihao Wang<EMAIL_ADDRESS>Institute for Functional Intelligent
Materials, NUS, Block S9, Level 9, 4 Science Drive, Singapore 117544,
Singapore  Department of Materials Science and Engineering, National
University of Singapore, 117544, Singapore  Pablo M. Perez-Piskunow
<EMAIL_ADDRESS>Catalan Institute of Nanoscience and
Nanotechnology - ICN2, (CSIC and BIST), Campus UAB, Bellaterra, Barcelona
08193, Spain  Calvin Pei Yu Wong  Institute of Materials Research and
Engineering (IMRE), Agency for Science Technology and Research (A*STAR), 2
Fusionopolis Way, Innovis #08-03, Singapore 128634, Republic of Singapore
Matthew Holwill  Department of Physics and Astronomy, University of
Manchester, Oxford Road, Manchester M13 9PL, UK  Jiawei Liu  Institute for
Functional Intelligent Materials, NUS, Block S9, Level 9, 4 Science Drive,
Singapore 117544, Singapore  Wei Fu  Institute of Materials Research and
Engineering (IMRE), Agency for Science Technology and Research (A*STAR), 2
Fusionopolis Way, Innovis #08-03, Singapore 128634, Republic of Singapore
Junxiong Hu  Department of Physics, National University of Singapore, 2
Science Drive 3, Singapore 117551, Republic of Singapore  T Taniguchi
National Institute for Materials Science, Ibaraki 305-0044, Japan  K Watanabe
National Institute for Materials Science, Ibaraki 305-0044, Japan  Ariando
Ariando  Department of Physics, National University of Singapore, 2 Science
Drive 3, Singapore 117551, Republic of Singapore  Lin Li  Institute for
Functional Intelligent Materials, NUS, Block S9, Level 9, 4 Science Drive,
Singapore 117544, Singapore  Kuan Eng Johnson Goh  Institute of Materials
Research and Engineering (IMRE), Agency for Science Technology and Research
(A*STAR), 2 Fusionopolis Way, Innovis #08-03, Singapore 128634, Republic of
Singapore  Department of Physics, National University of Singapore, 2 Science
Drive 3, Singapore 117551, Republic of Singapore  Stephan Roche
<EMAIL_ADDRESS>Catalan Institute of Nanoscience and Nanotechnology -
ICN2, (CSIC and BIST), Campus UAB, Bellaterra, Barcelona 08193, Spain  ICREA
Institucio Catalana de Recerca i Estudis Avancats, Barcelona 08010, Spain
Jeil Jung<EMAIL_ADDRESS>Department of Physics, University of Seoul,
Seoul 02504, Korea  Department of Smart Cities, University of Seoul, Seoul
02504, Korea  Konstantin Novoselov<EMAIL_ADDRESS>Institute for Functional
Intelligent Materials, NUS, Block S9, Level 9, 4 Science Drive, Singapore
117544, Singapore  Department of Materials Science and Engineering, National
University of Singapore, 117544, Singapore  Nicolas Leconte
<EMAIL_ADDRESS>Department of Physics, University of Seoul, Seoul
02504, Korea

###### Abstract

We report intriguing and hitherto overlooked low-field room temperature
extremely large magnetoresistance (XMR) patterns in graphene/hexagonal boron
nitride (h-BN) superlattices that emerge due to the existence of open orbits
within each miniband. This finding is set against the backdrop of the
experimental discovery of the Hofstadter butterfly in moiré superlattices,
which has sparked considerable interest in the fractal quantum Hall regime. To
cope with the challenge of deciphering the low magnetic field dynamics of
moiré minibands, we utilize a novel semi-classical calculation method,
grounded in zero-field Fermi contours, to predict the nontrivial behavior of
the Landau-level spectrum. This is compared with fully quantum simulations,
enabling an in-depth and contrasted analysis of transport measurements in
high-quality graphene-hBN superlattices. Our results not only highlight the
primary observation of the open-orbit induced XMR in this system but also shed
new light on other intricate phenomena. These include the nuances of single
miniband dynamics, evident through Lifshitz transitions, and the complex
interplay of semiclassical and quantum effects between these minibands.
Specifically, we document transport anomalies linked to trigonal warping, a
semiclassical deviation from the expected linear characteristics of Landau
levels, and magnetic breakdown phenomena indicative of quantum tunneling, all
effects jointly contributing to the intricacies of a rich electronic landscape
uncovered at low magnetic fields.

Suggested keywords

###### pacs:

33.15.Ta

## Introduction

The revelation of a fractal electronic spectrum in graphene superimposed on
hexagonal boron nitride (hBN) Dean  _et al._ (2013); Ponomarenko  _et al._
(2013a); Hunt  _et al._ (2013a) marked the inception of a decade-long research
odyssey into the properties of the Hofstadter butterfly Hofstadter (1976).
This research endeavor has primarily revolved around the strategic
introduction of superlattices in pristine graphene, leading to the intricate
division of its intrinsic $\pi-$electron bands into a series of moirÃ©
minibands Dean  _et al._ (2013); Ponomarenko  _et al._ (2013a); Hunt  _et al._
(2013a). The fractal’s self-similar nature, as posited in the original
Hofstadter model Hofstadter (1976), suggests an endless emergence of spectrum
replica features. However, this pursuit is constrained physically by some
fundamental mechanisms that are disregarded in the idealized mathematical
fractal model. These mechanisms, either resulting from the dynamics within
individual minibands or from interactions between different minibands, present
distinct experimental phenomena acting on high-energy minibands at low
magnetic fields. Over the past decade, research on graphene’s Landau levels
has predominantly concentrated on high-field and low-energy regions. Indeed,
because the onset of the quantum hall effect is influenced by disorder, high-
energy minibands have been challenging to map clearly due to high requirements
in measurement precision and due to the too large samples that are needed for
theoretical simulation at low fields. Woods  _et al._ ; Amet  _et al._ (2013);
Hunt  _et al._ (2013a); Wang  _et al._ (2016); Yankowitz  _et al._ (2014);
Hunt  _et al._ (2013b); Ponomarenko  _et al._ (2013b); Wallbank  _et al._
(2013); Chen  _et al._ (2014); Ni  _et al._ (2019); Song  _et al._ (2013);
Ribeiro-Palau  _et al._ (2018); San-Jose  _et al._ (2014); Kim  _et al._
(2018); Lin and Ni (2019); Long  _et al._ (2022).

Here, by fabricating ultraclean graphene superlattice samples that minimize
disorder and conducting ultra-high-precision transport measurements that
mitigate thermal noise, we have obtained a clear experimental mapping in the
low-field and high-energy regions. Meanwhile, theoretically using an
innovative semi-classical calculation technique - which exploits zero-field
calculations to predict behavior under finite magnetic fields, hence bypassing
the impracticality of too large samples required for full quantum simulations
- we embark on a systematic and comprehensive analysis of the electronic
behavior in the same regime. We unravel a knowledge frontier that goes beyond
the limits of fractality, entering a region where Lifshitz transitions play a
pivotal role. Specifically, the experimental manifestation of open orbits at
saddle point vHS is seminally identified as XMR persisting up to room
temperature, and the Landau level splitting between the hole-side Secondary
Dirac Point (h-SDP) and the neighboring vHS is attributed to trigonal warping
captured by the lattice-reconstructed tight-binding model. We finally also
identify a semi-classical deviation from the linear characteristics typically
expected for Landau levels, as evidenced in the Wannier representation of the
Hofstadter diagram, and observe signatures of magnetic breakdown, indicative
of quantum tunneling phenomena.

## G-hBN: full spectrum

Figure 1: Full spectrum magnetic-field dependence of electronic states in
G-hBN. a, $\sigma_{\text{xx}}$ (top) and $R_{\text{xy}}$ (bottom) transport
measurements. The black dashed line represents the position of SDP for hole
and electron side,respectively. The region enclosed by grey dashed lines
labels the region where Landau level splitting occurs. The region enclosed by
the purple dashed lines points to the place where the non-linear Wannier
diagram appears. The four vHS positions, manifested as conductivity dips and
Rxy sign reversal, are denoted by A, B, C, and D. b,. Fully quantum ($Q$)
numerical charge-carrier-dependent capacitance calculation using Lanczos
recursion. c, Charge-carrier dependent semi-classical ($SC$) LL predictions
using the zero-field Fermi contour behavior. Each color represents the LLs
belonging to a specific minibands.

In Fig. 1a, we present the full spectrum of the high-resolution transport
$\sigma_{\text{xx}}$ and $R_{\text{xy}}$ measurements on the graphene-hBN
superlattice. In addition to the usual Hofstadter butterfly signatures, we
observe a variety of new features in the $[0-5]$T regime at the onset of the
Quantum Hall Effect that have been overlooked in past experiments. In the
$R_{\text{xy}}$ map, the positions of the vHS denoted as A, B, C, and D are
determined by a sign reversal, corresponding to the conductivity dips in the
$\sigma_{\text{xx}}$ map. In the $R_{\text{xy}}$ map, the blue double arrows
mark the distance between the vHS and SDP on both sides, indicating the
asymmetry between the electron and hole side that leads to markedly different
observations. Both to the left of vHS A and around vHS D, a surprising
reshaping of the linear Landau levels is observed in the regions enclosed by
purple dashed lines. Indeed, a non-linearity occurs along with their
broadening as the magnetic field increases. Additionally, between the vHS A
and the h-SDP in the region enclosed by dashed grey lines, the linear Landau
levels split. All these phenomena revolve around vHS, a type of Lifshitz
transition, or within regions where bands overlap, and thus cannot be
explained by the single isolated band theory used for the Hofstadter
butterfly. In the right panels, both the quantum (Fig. 1b) and the semi-
classical (Fig. 1c) charge carrier density-dependent pictures are showcased
through the integration of the energy-dependent DOS from Fig. 2d. The semi-
classical simulation agrees well with the non-linear Landau levels observed in
Fig. 1a. Furthermore, the crossing of Landau levels belonging to different
minibands is also clearly illustrated. In the zero-field semi-classical
picture, all four vHS are vertical in the $V_{g}(B)$ map, yet this is markedly
different from the experiments. Indeed, the reason for the tilting at A and D
will be rationalized later. To this end, we need to first review how our semi-
classical approach works.

## Semi-classical approach

Figure 2:  Semi-classical simulations highlighting Lifshitz transitions and
miniband overlap. Panels a and b focus on a square-lattice single band toy
model, c and d correspond to the realistic G-hBN model containing multiple
minibands. a, Electronic band structure and its corresponding DOS. b, Semi-
classical predictions of the LL (white dots) with the full quantum
calculations in the background (dark colors correspond to a large number of
states while light colors correspond to the gaps separating the LL). The
dashed lines correspond to the positions of vHS where the black line based on
the semi-classical model is vertical and the red line drawn based on the
quantum simulations deviates to the right when the field is increasing. c, 3D
and 2D electronic band structure and DOS, indicating the electron-hole
asymmetry, and the existence of a variety of minibands potentially separated
by a gap or overlapping with each other. d, the semi-classical LL predictions
where we only plot the first 10 LLs for each energy. The different colors
represent the LLs belonging to different minibands from c. The green arrows
indicate the position of the saddle point vHS separating the electron and hole
character of each miniband.

To properly understand the capability and limitation of our semi-classical
method, pivotal for revealing the intra- and inter-miniband dynamics in the
moire superlattice, we first illustrate its applicability on a single band
square lattice model (see Methods for details). The model’s 3D bandstructure
and DOS (shown in Fig. 2a) exhibit a vHS where the nature of this single band
model changes from electrons to hole character, denoted by vHS1 followed by
vHS2 where distinct Fermi contours coalesce into one. Employing a Marching
Squares algorithm, we uniquely associate each eigenenergy with a specific
Fermi contour. We then compute the area enclosed by each of these contours and
semi-classically relate it to a magnetic length and magnetic field to obtain a
semi-classical prediction for the LL, grounded in the bands’ behavior in the
absence of a magnetic field (see Methods for details). In Fig. 2b, we finally
overlay these predictions with the quantum treatment results in a charge-
carrier density representation by integration of the DOS. It is important to
note that the semi-classical lines predict the existence of states rather than
mirroring the absence of states within gaps separating the Landau levels that
follow the well-known $n/n_{0}=t\phi/\phi_{0}+s$ Diophantine relation between
the normalized charge carrier density $n/n_{0}$ and the normalized magnetic
flux $\phi/\phi_{0}$. Consequently, these semi-classical predictions should be
contrasted with the quantum depiction’s dark-colored regions where the number
of states peaks.

Overall, we observe a noteworthy congruence between the semi-classical
approach premised on zero-field calculations and quantum calculations where
the magnetic field is explicitly incorporated via the Peierls phase
correction. Two regimes emerge where the concordance between semi-classical
and quantum methodologies appears to diminish. At vHS1, where the band
transitions from hole to electron character, the semi-classical LL predictions
cannot be resolved by our algorithm, unlike with the quantum approach that
shows intersecting and overlapping Landau level fans originating from each of
the band edges. At vHS2 where the white dots (dark regions) split towards
lower values of $n$ into doubled features due to the Lifshitz transition the
semi-classical perspective anticipates this occurring at a constant $n$ value
(dashed dark line), irrespective of the magnetic field, while the quantum view
indicates a shift towards higher $n$ values with increasing magnetic field $B$
(dashed red line). The latter discrepancy arises because our semi-classical
approach does not directly incorporate magnetic field corrections, and hence
cannot capture magnetic field-induced band renormalization.

These initial evaluations using a toy-model underscore the efficacy of the
semi-classical approach, prompting us to apply the same methodology to zero-
degree aligned G-hBN, as depicted on the right side of Fig. 2, where,
similarly, two kinds of vHS are expected to occur within each miniband.
Firstly, each of the minibands represented in different colors in the 3D band
structure in Fig. 2c contain an electron-to-hole character transition. These
occur at the saddle points indicated by the green arrows. This transition
becomes more apparent from the semi-classical LL predictions in Fig 2d (using
the same single color for each miniband as in Fig 2c) showing similar behavior
as in vHS1 from Fig 2b using the toy model. Secondly, the jumps in the number
of Fermi contours at each energy in the bottom panel of Fig 2d as counted by
our algorithm, suggest the occurrence of Lifshitz transitions where the
topology of a miniband band changes. In addition to these effects that match
with the single-band toy-model effects, we finally also observe in Fig 2c and
d several regions where minibands are overlapping (e.g. between the green-
shaded bands on the hole side or between the red-shaded bands on the electron
side), a requirement for the observation of inter-miniband phenomena.
Encouraged by this diverse array of features captured by the semi-classical
approach, we will delve in Figs. 3, 4 and 5 into specific regions of the
spectrum to explain the experimental features observed in Fig. 1. We
categorize them depending on the distinction if they stem from intra-miniband
Lifshitz transitions or from the interaction between different minibands.

## Intra-miniband Lifshitz physics

### Open orbits

Figure 3:  Open-orbit induced XMR persisting up to room temperature. a,
$R_{\text{xx}}$ Rxx at low magnetic fields. The positive XMR occurs at the
positions A, B, C, D and E where open orbits exist within each miniband.
Inset: Fermi contours for each miniband showing the existence of these open
orbits whenever a transition from hole to electron character of the bands
occurs at a vHS in the system. We repeat the Fermi contours in a $3\times 3$
tiling where we use a divergent colormap for better visualization. Red color
corresponds to hole-type carriers and blue color to electron-type carriers
putting the saddle point at the white colour. The dashed lines indicate as a
guide to the eye the path corresponding to an open orbit. We note that for our
choice in commensuration cell, the primary Dirac cones are folded back onto
the $\Gamma$-point. Only one band out of 2 energy-degenerate ones is shown
here. b, The linear fitting of the XMR of the peak resisstance values for E, A
and B. c, $R_{\text{XX}}$ at $0$ and $1$T at room temperature where the bottom
panel weakly shows the residual XMR signal at E. Much stronger signals are
observed up to $200$K, as illustrated in the Supplemental material.

In Fig. 3a, we measure the XMR at low magnetic fields from $0$T to $0.65$T
where Shubnikov de Haas oscillations are not yet significant. At the saddle
point vHS indicated by A, B, C, D and E, the resistivity increases
significantly with magnetic field. The corresponding positive XMR behaviors
for hole-side vHS are shown with linear fitting to magnetic fields in Figs.
3b. Since the distances between the SDP and vHS on this hole side are larger
than those the electron side (as already highlighted in Fig. 1a), the
amplitude of the peaks corresponding to the vHS are less influenced by the
more significant SDF peak. The values of XMR at $1.6$K achieve $10000\%$,
$3000\%$ and $800\%$ at $0.65$T for vHS E, A, and B, respectively. Due to C
and D being closer to the SDP, their positive XMR is more prone to smearing
when the temperature increases. Therefore we focus on vHS E when checking the
temperature dependence to conclude on an XMR that can persist up to room
temperature, as demonstrated in Figs. 3c. The room temperature signal at E is
very weak, but we illustrate in the Supplemental material how it strongly
survives up to $200$ K.

The linearity of the XMR is a significant hallmark symbolizing the emergence
of open orbits Takatsu  _et al._ (2013); Zhang  _et al._ (2019); Chambers
(1990); Pippard (1989); Lifshits  _et al._ (1973); Abrikosov (1988). Indeed,
the gate voltage at which each of these occurs perfectly matches with the
saddle point vHS that indicate a transition from electron to hole character
Markiewicz (1994). These transitions give rise to open orbits at zero field in
reciprocal space as illustrated by the insets Figs. 3a where the open orbits
are demarcated by dashed black lines within the white regions of the colormap.
The positive XMR in open orbits can be rationalized as follows. Because a
magnetic field induces a force on a moving electron, its momentum changes, and
thus its path in real space also gets affected. If electrons are on closed
orbits in reciprocal space, the change in moment following the k-path is
periodic, and thus overall, the net change in motion in real space is
negligible. Electrons move almost freely through the sample unaffected by the
magnetic field. However, if they are on a closed orbit in reciprocal space,
the change in momentum will be mostly in one primary direction. The magnetic
field thus continually modifies the momentum of an electron at this energy in
a specific direction. This leads to an increase in scattering and less
predictable electron motion. The electrons don’t move freely through the
sample and can experience more scattering events. This in turn explains the
positive XMR. By relating the observed XMR to the occurrences of open orbits,
we provide the first documentation of such features in layered moiré systems.
Notably, the observation of open orbits does not require quantum mechanical
treatment for validation. When comparing the energy at which each of these
open orbits occur in Figs. 1a and c we nonetheless note that the zero-field
semi-classical approach does not capture their weak field-dependence due to
the omittance of magnetic field-induced band renormalization. To strengthen
the attribution of the XMR to open orbits, nine other similar samples with
graphene doubly aligned to both sides with h-BN at various twist angles are
provided in the Supplemental material showing similar observations.

### Trigonal-warping

Figure 4: Trigonal warping transport signatures at negative gate voltages. a,
$R_{\text{xx}}$, $R_{\text{xy}}$, $\frac{d(1/\sigma_{\text{xx}})}{dB}$
measurements and the sketch based on $\frac{d(1/\sigma_{\text{xx}})}{dB}$
showing splitting of Landau levels. b Semi-classical confirmation of splitting
where we plot the $30$ first LLs for each energy. c, sketch based on b. d,
Fermi contours showing the split originates in the trigonal-warping induced
merging of Fermi pockets indicated by green arrows where the relaxed model
shows stronger features than when using the rigid model parameters.

In Fig. 4a, our attention is directed to the very low-B regime at negative
gate voltage corresponding to the region demarcated by grey dashed lines in
Fig. 1a where the current ultra-high resolution measurements exhibit
surprising LL level splittings at the SDP: initially, two lines split into
four, eventually evolving into eight lines. This pattern can be elucidated by
the semi-classical predictions in Fig. 4d that show that at lower energy in
the red-colored region of the map, three separate Fermi pockets exist around
the K-point (see green arrows). These pockets merge into a single pocket as
the energy slightly increases. The relaxed parameters that are used throughout
this paper show clearer signatures of this effect compared to the rigid model
parameters indicative of the importance of lattice reconstruction in this
system. This fusion of contours that happens at a Lifshitz transtion has been
associated in our toy-model with a reduction in the number of LLs where two
pockets merged into one. Given that here three pockets coalesce into one, it
logically follows that a single LL should split into three levels, a
hypothesis supported by the semi-classical LL predictions in Fig. 4b and the
corresponding sketch in Fig. 4c. From an experimental perspective, the
progression from 2 (approximately 3) observable bands to 4 (approximately 6)
and then to 8 (approximately 9) visible bands implies that, in reality, two
pockets merge initially before joining with the third pocket. This sequence of
events could potentially be attributed to minor non-uniform strain within the
sample. The uncertainty on the exact values stems from the difficulty of
extracting very fine features. As a side note, a similar observation can be
made again here in terms of the semi-classical approach, namely that it
ignores any magnetic field-induced renormalization of bands, missing the
experimentally observed gate voltage-dependence of the splitting behavior.

The origin of these three pockets can be attributed to trigonal warping
KostyaScience, defined as a deformation of the circular Fermi cones around the
K-points that tends to be trigonal due to the lattice symmetries. Such
trigonal warping usually has a substantial effect on the bandstructure of
graphene, hence explaining its observation through ARPES measurements Ortix
_et al._ (2012). Its electronic effects on the contrary are often subtle
making it challenging for transport measurements to detect them. Here by
matching the observation of these subtle features in the experiments with the
trigonal warping-induced deformation of this higher energy moire mini band
beyond the h-SDP, we provide the first experimental confirmation of trigonal
warping effects by means of transport measurements.

## Inter-miniband physics

### Non-linear Landau levels

Figure 5: Non-linear features and magnetic breakdown on the electron side. a,
$R_{\text{xx}}$ measurements with focus on positive gate voltage beyond the
e-SDP showing non-linear features emerging at $5.5n/n_{0}$ both towards the
left and right sides of the spectrum. b, Similar non-linear $R_{\text{xy}}$
features showing red-green color alternation indicating overlap between hole
and electron character charge carriers. c, Semi-classical (SC) Landau level
predictions confirming the non-linear behavior of the experimental features
including the overlapping hole and electron character part of the bands. We
only plot the ten first LLs for figure clarity. We mirror the positive
magnetic field calculations to obtain the corresponding negative field bands
for easier comparison with a and b. The inset shows a zoom on the relevant
bands from Fig. 2c highlighting the multitude of overlapping bands around D.

Turning our attention here to inter-miniband phenomena, we focus on the
measurements in Fig. 5 corresponding to the purple dashed regions from Fig. 1a
to investigate the departure from non-linear Landau levels in the Wannier
representation of the Hofstadter butterfly. Although older works have observed
such signatures, no clear origin has been provided Lu  _et al._ (2020).The
electron side features giving the clearest $R_{\text{xx}}$ and $R_{\text{xy}}$
signatures in Figs. 5a and b respectively, we draw from the alternating colors
corresponding to hole and electron charge carriers in $R_{\text{xy}}$ that
miniband overlap is happening here. Fig. 5c confirms this attribution through
the dense crossing between brown, orange, pink, and purple bands in the
$[4-7]n/n_{0}$ range to 7 n/n0 (brown, orange, pink, and possibly purple
colors in c). The colors match the colors from the minibands in the inset.

Because no magnetic field or quantum effect is explicitly included in this
semi-classical approach, the non-linear effect must be explained without
including quantum effects such as magnetic tunneling at the origin of magnetic
breakdown. Reminding that the position of the Landau levels in the Wannier
representation depends solely on the number of electrons contained within a
band, the specific shape of each band does not affect the linear dependence.
However, when two (mini-)bands coexist at the same energy, the filling of the
two bands happens concurrently, hence while filling the first band, some
electrons will start filling the second band hence disturbing the linear
behavior of the first band (and vice-versa).

### Magnetic breakdown

Beyond the non-linear effects in Fig. 5 that are already captured at zero-
field, our semi-classical approach cannot capture quantum tunneling between
different bands that can happen at finite magnetic field. However, the
experimental maps show clear broadening of the non-linear features with
magnetic field when alternating between the overlapping electron and hole
character bands. This is commonly referred to as magnetic breakdown Cohen and
Falicov (1961); Falicov and Sievert (1965); Pippard (1965); Blount (1962). The
quantum results for this low B-field regime are not shown, as computational
constraints hinder the resolution of similar features, likely due to the
requirement of exceedingly large samples for very low magnetic fields and
extensive cyclotron orbits. Yet, the magnetic treatment in Fig. 1b shows that
this approach captures the experimental magnetic broadening breakdown effects
in the $[1-6]$T regime hence overcoming this limitation of the semi-classical
approach and confirming the quantumness of these smoothened bands. We finally
note that experimentally, nonlinear features and magnetic breakdown are
particularly pronounced in double moire systems at low-magnetic field due to
hBN encapsulation of the graphene layer, as illustrated in the Supplemental
material, rationalized by the fact that the additional superlattice features
increase the occurrence of mini-band overlap at the root of this peculiar
behavior. As a side note, in the semi-classical picture from Fig. 5c, the
brown band originating at $4n/n_{0}$ also displays a signature of trigonal
warping-induced Lifshitz transitions causing the LLs to split around
$4.5n/n_{0}$ as was observed on the hole side in Fig. 4. However, the effect
is much weaker and cannot be resolved in the current measurement.

## Conclusions and Outlook

From a fundamental point of view, our analysis of XMR measurements in
graphene/hBN superlattices provides compelling evidence of fundamental
limitations to accessing endless intricacies of the fractal quantum spectrum
in moire superlattices, as highlighted through seminal observations of
phenomena that occur in the very low B-field regime before the quantum Hall
effect can fully develop. From an application point of view, we particularly
want to highlight the open orbit-induced XMR signatures that survive up to
room temperature, hence providing a promising playground to engineer similar
behavior at lower gate voltages for use in magnetic field sensors, MRAM and
logic gates Gopinadhan  _et al._ (2013); Gani  _et al._ (2020); Cho and Fuhrer
(2008); Friedman  _et al._ (2010); Liao  _et al._ (2012). From a
methodological point of view, the theoretical analysis of high-quality G-hBN
samples puts in perspective the strengths and limitations inherent to semi-
classical methodologies when juxtaposed with quantum mechanical calculations.
The utility of semi-classical approaches, particularly evident in the
elucidation of zero-field bands and Lifshitz transitions, offers a substantial
contribution to understanding complex physical phenomena. In addition to the
open orbit-induced XMR, these semi-classical calculations notably explain the
emergence of trigonal warping-driven transport anomalies and shed light on the
notable deviations from the linear B-field behavior expected for Landau levels
in the Wannier representation. The full quantum picture in turn allows to
capture magnetic breakdown due to tunneling between different bands. We expect
the semi-classical approach to be useful for a variety of layered materials
where Lifshitz transitions within single minibands and semi-classical and
quantum interactions between different minibands can trigger a variety of
exotic properties.

In the current understanding, the XMR associated with open orbits is primarily
interpreted through a reciprocal space framework of electronic band
structures. Our preliminary investigations, not detailed in this manuscript,
indicate that in real space these open orbits correspond to localized states
at the corners of the moiré pattern in this 2D material. This contrasts with
the more spatially extended states linked to closed orbits. Introduction of a
magnetic field results in a notable modulation of these spatial
characteristics: states associated with open orbits exhibit a tendency towards
delocalization, whereas those corresponding to closed orbits show enhanced
localization due to early onset of cyclotron orbits. The observed negative XMR
in the latter can still be coherently explained within the paradigm of weak
localization as the localization of states is very limited at these small
magnetic fields. Conversely, for open orbit states, the observed positive XMR
potentially aligns with a model involving edge states. These edge states are
hypothesized to interfere with each other, leading to an increase in
resistivity. Future studies are planned to substantiate this hypothesis by
directly visualizing, in real space, the current paths that correspond to open
orbits in reciprocal space.

## Methods

### Experiments

#### Device Fabrication

Monolayer graphene and hexagonal boron nitride with a thickness of $\sim 20$
nm were mechanically exfoliated on SiO2/Si substrates and identified under an
optical microscope. The flakes with a well-defined long and straight edge were
selectively utilized. The hBN-graphene-hBN heterostructures were assembled
using standard dry-transfer technique with PMMA-coated PDMS stamps. To obtain
the single-aligned samples, the straight edges of the top hBN and graphene
were aligned to zero. Simultaneously, the bottom hBN layer was intentionally
rotated to 15 degrees to prevent the formation of unintended superlattices. By
contrast, for doubly-aligned devices, all the selected straight edges of the
flakes were aligned in unison. All the heterostructures were pre-examined with
Raman spectrum to ensure the precise alignment. Subsequently, standard
electron beam lithography was carried out to define a Hall-bar geometry using
a PMMA mask, followed by CHF3/O2 ion etching. Cr/Au ($2/70$ nm) was thermally
deposited to form edge contacts.

#### Transport measurements

Transport measurements were performed in a dilution refrigerator with a base
temperature of $30$ mK. We used standard low-frequency lock-in techniques with
an excitation frequency of about $6-30$ Hz and an excitation current of about
$2-5$ nA. The current flowing through the sample was amplified by a current
pre-amplifier and measured by the lock-in amplifier. The four-probe voltage
was amplified by a voltage preamplifier at $\times 1,000$ gain and measured by
another lock-in amplifier. The Second Order Low Pass Filter was applied with
the R1 = $512\Omega$, C1 = $4.7$ nF, R2 = $816\Omega$, C2 = $2.2$ nF.

### Simulations

#### Semi-classical and quantum methods

In the semi-classical method, we associate each eigenenergy with a specific
Fermi contour within a specific band using a Marching Square algorithm as
illustrated in the Supplemental material. We then integrate the area of these
contours at zero magnetic field for each energy to obtain a semi-classical
prediction of each LL for this energy. The area in the case of graphene can be
linked to the magnetic length $l_{B}$ as

$A(E,N)=2\pi\frac{1}{l_{B}(E,N)^{2}}\left(N+\frac{1}{2}\right)$ (1)

where $N$ is the Landau level number. Conversely, we have

$l_{B}(E,N)=\sqrt{\frac{2\pi}{A(E)}\left(N+\frac{1}{2}\right)}.$ (2)

The magnetic length is then related to the magnetic field in Tesla when
$l_{B}$ is given in nm by

$B(E,N)=\left(\frac{2.6}{l_{B}(E,N)}\right)^{2}.$ (3)

By obtaining the magnetic field corresponding to the area for each energy and
each value of $N$, we obtain the dots that are plotted in Figs. 2d. For the
quantum calculations, the density of states (DOS) maps are obtained using
Lanczos recursion or KPM Fan  _et al._ (2021) (see Supplemental material for
details) on large systems containing on average 100 million atoms. This allows
to probe very fine meV resolution features. Quantum capacitance maps are
calculated using the equations outlined in the Supplemental material. The
transformation method to go from Hofstadter energy-dependent maps Hofstadter
(1976) to Wannier charge carrier density n-dependent maps Wannier (1978) for
both the semi-classical and quantum results is also given in the Supplemental
material. An open-source Python package that extends Varjas  _et al._ (2020)
the KPM capabilities of KWANT Groth  _et al._ (2014) is available for download
Pablo Piskunow (2023).

#### Tight-binding models

For the toy-model calculations from Fig. 2a and b, we engineer a square
lattice model that induces a change in band topology causing Lifshitz
transitions as well as a transition from the hole to electron character of the
charge carriers within this single band model. The first and second-neighbor
interactions (see Supplemental for sketch) are given by $t_{1}=-1$ [$t$],
$t_{2}=-0.5$ [$t$], $t_{3}=0.5$ [$t$] and $t_{4}=0.5$ [$t$]. For the G-hBN
calculations, because the accuracy of the TB model is critical for our
intended purpose of making qualitative comparisons with experiment, we use a
DFT-inspired TB model as introduced in Ref. Leconte and Jung (2020) based on a
relaxed-model continuum parametrization DaSilva  _et al._ (2015). For the
quantum treatment of the perpendicular magnetic field, the Peierls phase
correction to the hopping terms is included following the recipe introduced in
Ref. Cresti (2021).

###### Acknowledgements.

This work was supported by the Korean NRF through the Grant
NRF-2021R1A6A3A13045898 (N.L.) and by the Samsung Science and Technology
Foundation Grant No. SSTF-BAA1802-06 (J.J.). We acknowledge computational
support from KISTI Grant No. KSC-2022-CRE-0514 and by the resources of Urban
Big data and AI Institute (UBAI) at UOS. J.J. also acknowledges support by the
Korean Ministry of Land, Infrastructure and Transport (MOLIT) from the
Innovative Talent Education Program for Smart Cities. S.R acknowledges grant
PCI2021-122035-2A-2 funded by MCIN/AEI/10.13039/501100011033 and European
Union NextGenerationEU/PRTR, funding from the European Union’s Horizon 2020
research and innovation programme under grant No 881603 (Graphene Flagship),
and the support from Departament de Recerca i Universitats de la Generalitat
de Catalunya. ICN2 is funded by the CERCA Programme/Generalitat de Catalunya
and supported by the Severo Ochoa Centres of Excellence programme, Grant
CEX2021-001214-S, funded by MCIN/AEI/10.13039.501100011033.

## References

  * Dean  _et al._ (2013) C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone, and P. Kim, Nature 497, 598 (2013).
  * Ponomarenko  _et al._ (2013a) L. A. Ponomarenko, R. V. Gorbachev, G. L. Yu, D. C. Elias, R. Jalil, A. A. Patel, A. Mishchenko, A. S. Mayorov, C. R. Woods, J. R. Wallbank, M. Mucha-Kruczynski, B. A. Piot, M. Potemski, I. V. Grigorieva, K. S. Novoselov, F. Guinea, V. I. Fal’ko, and A. K. Geim, Nature 497, 594 (2013a).
  * Hunt  _et al._ (2013a) B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and R. C. Ashoori,  (2013a), 10.1126/science.1237240.
  * Hofstadter (1976) D. R. Hofstadter, Phys. Rev. B 14, 2239 (1976).
  * (5) C. R. Woods, L. Britnell, A. Eckmann, R. S. Ma, J. C. Lu, H. M. Guo, X. Lin, G. L. Yu, Y. Cao, R. V. Gorbachev, A. V. Kretinin, J. Park, L. A. Ponomarenko, M. I. Katsnelson, Y. N. Gornostyrev, K. Watanabe, T. Taniguchi, C. Casiraghi, H.-J. Gao, A. K. Geim, and K. S. Novoselov, Nat. Phys. 10, 451 456.
  * Amet  _et al._ (2013) F. Amet, J. R. Williams, K. Watanabe, T. Taniguchi, and D. Goldhaber-Gordon, Physical Review Letters 110, 216601 (2013).
  * Wang  _et al._ (2016) E. Wang, X. Lu, S. Ding, W. Yao, M. Yan, G. Wan, K. Deng, S. Wang, G. Chen, L. Ma, J. Jung, A. V. Fedorov, Y. Zhang, G. Zhang, and S. Zhou, Nature Physics 12 (2016), 10.1038/nphys3856.
  * Yankowitz  _et al._ (2014) M. Yankowitz, J. Xue, and B. J. LeRoy, Journal of Physics: Condensed Matter 26, 303201 (2014).
  * Hunt  _et al._ (2013b) B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and R. C. Ashoori, Science 340, 1427 (2013b), https://www.science.org/doi/pdf/10.1126/science.1237240 .
  * Ponomarenko  _et al._ (2013b) L. Ponomarenko, R. Gorbachev, G. Yu, D. Elias, R. Jalil, A. Patel, A. Mishchenko, A. Mayorov, C. Woods, J. Wallbank, M. Mucha-Kruczynski, B. Piot, M. Potemski, I. Grigorieva, K. Novoselov, F. Guinea, V. Fal’ko, and A. Geim, Nature 497, 594 (2013b).
  * Wallbank  _et al._ (2013) J. R. Wallbank, A. A. Patel, M. Mucha-Kruczynski, A. K. Geim, and V. I. Fal’ko, Physical Review B 87, 245408 (2013).
  * Chen  _et al._ (2014) Z. G. Chen, Z. Shi, W. Yang, X. Lu, Y. Lai, H. Yan, F. Wang, G. Zhang, and Z. Li, Nature Communications 5, 1 (2014).
  * Ni  _et al._ (2019) G. X. Ni, H. Wang, B.-Y. Jiang, L. X. Chen, Y. Du, Z. Y. Sun, M. D. Goldflam, A. J. Frenzel, X. M. Xie, M. M. Fogler, and D. N. Basov, Nature Communications 10, 4360 (2019).
  * Song  _et al._ (2013) J. C. W. Song, A. V. Shytov, and L. S. Levitov, Phys. Rev. Lett. 111, 266801 (2013).
  * Ribeiro-Palau  _et al._ (2018) R. Ribeiro-Palau, C. Zhang, K. Watanabe, T. Taniguchi, J. Hone, and C. R. Dean, Science 361, 690 (2018).
  * San-Jose  _et al._ (2014) P. San-Jose, A. Gutiérrez-Rubio, M. Sturla, and F. Guinea, Phys. Rev. B 90, 075428 (2014).
  * Kim  _et al._ (2018) H. Kim, N. Leconte, B. L. Chittari, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Jung, and S. Jung, Nano Letters 18, 7732 (2018).
  * Lin and Ni (2019) X. Lin and J. Ni, Phys. Rev. B 100, 195413 (2019).
  * Long  _et al._ (2022) M. Long, P. A. Pantaleón, Z. Zhan, F. Guinea, J. Ángel Silva-Guillén, and S. Yuan, npj Computational Materials 8, 73 (2022).
  * Takatsu  _et al._ (2013) H. Takatsu, J. J. Ishikawa, S. Yonezawa, H. Yoshino, T. Shishidou, T. Oguchi, K. Murata, and Y. Maeno, Phys. Rev. Lett. 111, 056601 (2013).
  * Zhang  _et al._ (2019) S. Zhang, Q. Wu, Y. Liu, and O. V. Yazyev, Phys. Rev. B 99, 035142 (2019).
  * Chambers (1990) R. G. Chambers,  _Electrons in Metals and Semiconductors_ (Chapman and Hall, London, 1990).
  * Pippard (1989) A. B. Pippard,  _Magnetoresistance in Metals_ (Cambridge University Press, New York, 1989).
  * Lifshits  _et al._ (1973) I. M. Lifshits, M. Y. Azbel, and M. I. Kaganov,  _Electron Theory of Metals_ (Consultants Bureau, New York, 1973).
  * Abrikosov (1988) A. A. Abrikosov,  _Fundamentals of the Theory of Metals_ (North Holland, Amsterdam, 1988) reprint.
  * Markiewicz (1994) R. S. Markiewicz, Journal of Physics: Condensed Matter 6, 3059 (1994).
  * Ortix  _et al._ (2012) C. Ortix, L. Yang, and J. van den Brink, Phys. Rev. B 86, 081405 (2012).
  * Lu  _et al._ (2020) X. Lu, J. Tang, J. R. Wallbank, S. Wang, C. Shen, S. Wu, P. Chen, W. Yang, J. Zhang, K. Watanabe, T. Taniguchi, R. Yang, D. Shi, D. K. Efetov, V. I. Fal’ko, and G. Zhang, Phys. Rev. B 102, 045409 (2020).
  * Cohen and Falicov (1961) M. H. Cohen and L. M. Falicov, Phys. Rev. Lett. 7, 231 (1961).
  * Falicov and Sievert (1965) L. M. Falicov and P. R. Sievert, Phys. Rev. 138, A88 (1965).
  * Pippard (1965) A. B. Pippard, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 287, 165 (1965), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1965.0174 .
  * Blount (1962) E. I. Blount, Phys. Rev. 126, 1636 (1962).
  * Gopinadhan  _et al._ (2013) K. Gopinadhan, Y. J. Shin, I. Yudhistira, J. Niu, and H. Yang, Phys. Rev. B 88, 195429 (2013).
  * Gani  _et al._ (2020) M. Gani, K. A. Shah, S. A. Parah, and P. Misra, Physics Letters A 384, 126171 (2020).
  * Cho and Fuhrer (2008) S. Cho and M. S. Fuhrer, Phys. Rev. B 77, 081402 (2008).
  * Friedman  _et al._ (2010) A. L. Friedman, J. L. Tedesco, P. M. Campbell, J. C. Culbertson, E. Aifer, F. K. Perkins, R. L. Myers-Ward, J. K. Hite, C. R. Eddy, Jr, G. G. Jernigan, and D. K. Gaskill, Nano Lett. 10, 3962 (2010).
  * Liao  _et al._ (2012) Z.-M. Liao, H.-C. Wu, S. Kumar, G. S. Duesberg, Y.-B. Zhou, G. L. W. Cross, I. V. Shvets, and D.-P. Yu, Adv. Mater. 24, 1862 (2012).
  * Fan  _et al._ (2021) Z. Fan, J. H. Garcia, A. W. Cummings, J. E. Barrios-Vargas, M. Panhans, A. Harju, F. Ortmann, and S. Roche, Physics Reports 903, 1 (2021), linear scaling quantum transport methodologies.
  * Wannier (1978) G. H. Wannier, Physica Status Solidi (b) 88, 757 (1978).
  * Varjas  _et al._ (2020) D. Varjas, M. Fruchart, A. R. Akhmerov, and P. M. Perez-Piskunow, Phys. Rev. Res. 2, 013229 (2020).
  * Groth  _et al._ (2014) C. W. Groth, M. Wimmer, A. R. Akhmerov, and X. Waintal, New Journal of Physics 16, 063065 (2014).
  * Pablo Piskunow (2023) Pablo Piskunow, “Kpm-tools documentation,”  (2023), accessed: 2023-12-03.
  * Leconte and Jung (2020) N. Leconte and J. Jung, 2D Materials 7, 031005 (2020).
  * DaSilva  _et al._ (2015) A. M. DaSilva, J. Jung, S. Adam, and A. H. MacDonald, Phys. Rev. B 91, 245422 (2015).
  * Cresti (2021) A. Cresti, Phys. Rev. B 103, 045402 (2021).