# The Edge of Random Tensor Eigenvalues with Deviation

Nicolas Delporte<EMAIL_ADDRESS>Okinawa Institute of Science and
Technology Graduate University, 1919-1, Tancha, Onna, Kunigami District,
Okinawa 904-0495, Japan. Naoki Sasakura<EMAIL_ADDRESS>Yukawa
Institute for Theoretical Physics, Kyoto University,Kitashirakawa, Sakyo-ku,
Kyoto 606-8502, Japan.  CGPQI, Yukawa Institute for Theoretical Physics, Kyoto
University,Kitashirakawa, Sakyo-ku, Kyoto 606-8502, Japan.

###### Abstract

The largest eigenvalue of random tensors is an important feature of systems
involving disorder, equivalent to the ground state energy of glassy systems or
to the injective norm of quantum states. For symmetric Gaussian random tensors
of order $3$ and of size $N$, in the presence of a Gaussian noise, continuing
the work [1], we compare the edge of the genuine eigenvalue distribution with
that of the signed distribution, using field theoretic methods at large $N$
combined with earlier rigorous results of [2]. We support our claims with
Monte Carlo simulations. We find two critical points as the variance of the
noise increases, the first of which corresponding to the emergence of an
outlier from the main part of the spectrum and the second where this outlier
merges with the corresponding trivial eigenvalue and they both become complex.

##  1 Introduction

While random matrices have become today a major tool in sciences, through a
tremendous range of applicability 111[3] is an excellent introductory
reference to grasp the breadth of topics intimately related to random
matrices, from chaos, RNA folding and QCD to random surfaces, knot theory and
free probability., random tensors, omnipresent in classical and quantum
networks, are still lacking the coherence that binds the matrices. Besides
their use in data sciences, random tensors also serve as a discrete approach
to quantum gravity [4, 5, 6] (a recent review is [7]) generating simplicial
complexes from Feynman diagram expansions, attempting to generalize in higher
dimensions the success of random matrices to produce and solve models of
random surfaces with a large $N$ expansion that discriminates surfaces of
different topology (a modern treatment is done in [8]).

With tensors of size $N$ and of any order $p\geq 3$, an important milestone
was achieved in [9, 10, 11] that set up an action invariant under a tensor
product $G^{\otimes p}$, for some compact Lie group $G$, leading to an
analogous large $N$ expansion 222The dependence in $N$ is coming from the
choice of representation of the group $G$ or from the dimension of $G$.,
separating this time classes of graphs of different _degree_. The leading
order was shown to consist of _melonic_ graphs, enjoying a simplifying
recursive structure that can be solved explicitly in that limit. This melonic
limit turned out in fact to form a strong universality class competing with a
two-dimensional Brownian sphere phase obtained with matrix-like interactions
(and a more singular intermediate branched polymer phase) [12]. Escaping the
melonic universality class is still an open problem.

Matrices are often analysed through the lens of their spectral properties
[13], revealing aspects of universality and integrability. For tensors,
corresponding spectra are not as straightforward, since firstly there exist
several definitions of eigenvalues [14] and secondly, the number of the
eigenvalues is typically exponential in the size of the tensor [15]. Other
approaches involve studying the invariants under the tensor symmetry group
[16], although they grow factorially with the number of tensors contracted
together [17, 18]. In fact, most tensor problems that generalize the linear
systems of equations and inequalities of matrix problems (real eigenvalues,
spectral norm, the best rank-1 approximation, etc.) are NP-hard [19]. This
offers another reason motivating the search of efficient approaches in the
study of random tensors.

For a symmetric real tensor $C_{i_{1}\dots i_{p}}$ $(1\leq
i_{1},\dots,i_{p}\leq N)$ of order $p$, a Z-eigenpair [20] is a real solution
to the equations

$\displaystyle\sum_{1\leq i_{2}\dots i_{p}\leq N}C_{i_{1}\dots
i_{p}}w_{i_{2}}\dots w_{i_{p}}=zw_{i_{1}}\,,\quad(1\leq i_{1}\leq
N,z\in\mathbb{R},w\in\mathbb{R}^{N},\absolutevalue{w}=1),$ (1)

that is they form the zero set of polynomials on $S_{N-1}$. When the tensor
has random Gaussian entries, the polynomials are called Kostlan polynomials
(see [21] for the expectation value of the cardinality of the set, its
variance at large $N$, as well as an exhaustive chronology of the results,
originating with Kac [22] and Rice [23]; [24] offers a different approach
providing an exact counting of the average number of roots, from which
asymptotic expressions at large $N$ are however difficult to extract). The
problem is identical to the determination of the critical points of pure
spherical $p$-spin models first solved with replicas [25] and TAP equations
[26]. In [2] random matrix techniques allowed to obtain rigorously asymptotics
and finite $N$ expressions of the expected number of critical points (its
normalized logarithm is called in the spin glass community the “complexity”)
and the asymptotic expected distribution of critical points of given finite
index. Fluctuations around the ground state were described in [27] as
following a Gumbel distribution. We refer to [28] for a comprehensive overview
of the history in understanding the critical points of high-dimensional random
landscapes. The scaling of the largest real eigenvalue $z_{max}\sim\sqrt{N}$
was also obtained from combinatorial methods [29].

Offering a playground that connects complexity theory, glasses, neural
networks, etc. [30], many variations of the spherical $p$-spins have been
studied: mixed, Ising, with external field…333See [31] for a concise guide to
the recent literature. The last case corresponds to a Hamiltonian of the form

$H=-\sum_{1\leq i_{1}\dots i_{p}\leq N}J_{i_{1}\dots i_{p}}x_{i_{1}}\dots
x_{i_{p}}-\sum_{i=1}^{N}h_{i}x_{i},\quad\absolutevalue{x}^{2}=N,\quad\expectationvalue{J_{i_{1}\dots
i_{p}}^{2}}=\frac{J^{2}}{pN^{p-1}}\,,$ (2)

where the external field $(h_{i})_{1\leq i\leq N}$ can be deterministic or
random, of norm or standard deviation $h$. It was shown [32, 33] that for $h$
large enough, the number of critical points is trivial (there is only one
maximum and one minimum) whereas after an explicit threshold, it becomes
exponential. In the presence of a quadratic potential [34], the density and
distribution of stationary points grow more and more complex, as the deviation
from a unit mass scales with $N^{-1/6}$ or $N^{-1/2}$.

Rescaling $z=1/\absolutevalue{v}^{p-2},w=v/\absolutevalue{v}$ in (1) and
adding a Gaussian random deviation $\nu=(\nu_{i})_{1\leq i\leq N}$ of variance
$2\beta$ leads to the following eigenvector equation

$\displaystyle\sum_{1\leq i_{2}\dots i_{p}\leq N}C_{i_{1}\dots
i_{p}}v_{i_{2}}\dots v_{i_{p}}=v_{i_{1}}+\nu_{i_{1}}\,,\quad(1\leq i_{1}\leq
N,v\neq 0,v\in\mathbb{R}^{N})\,.$ (3)

The deviation can arise as a noise or an error as one tries to solve the
eigenvector equation. It bears resemblance to the case of an external
potential, although the outcomes and the appropriate large $N$ behaviour turn
out different. We observe that as the variance of the noise increases, the
second largest eigenvalue emerges from the non-trivial main body of the
spectrum and later merges with the corresponding trivial first largest
eigenvalue, in a sense becoming an outlier.

In the case of random matrices, it is well-known [35] that outliers escape the
limiting spectral distribution of a random matrix, the noise, in the presence
of an added deterministic matrix, the signal, when the eigenvalues of the
signal matrix are above a given threshold, corresponding to the BBP transition
[36]. In other words, information about the signal can be recovered, when it
is sufficiently strong compared to the noise. It was also shown [36] that the
fluctuations of the largest eigenvalue differ from the Tracy-Widom law of the
Wigner ensemble when the signal is detectable.

Another example of a hard tensor problem is the spiked tensor model [37],
asking for the detectability and recovery threshold of a given signal vector,
hidden in a tensor to which a random tensor is added. In contrast with the
matrix case that possess an easy (with power iteration) and an impossible to
recover the signal phase, separated by the BBP transition, it seems that the
tensor case contains an intermediate “hard” regime, above the information-
theoretically possible threshold, where there is no polynomial time algorithm
that recovers the signal [38] (see however [39] for some promising methods
based on tensor invariants).

Focusing on tensors of order 3, one of the authors recovered the distribution
of [2] with a field theoretic approach at large [40] and finite $N$ [41],
reviving the earlier techniques that were approximating the absolute value of
the determinant with the determinant (what below we call the signed
distribution) or the use of replicas (see e.g. [42]). In [1], similar
techniques were extended to the situation where the tensor has a non-zero
background value 444In the case where the background was a rank one tensor,
the detectability threshold of [38] was recovered. and a Gaussian noise was
introduced in the eigenvector equation, work that we will pursue. More
recently, by matching the edge (signalling the largest eigenvalue) and the
critical point of the signed distribution to those of the genuine spectral
density, the case was made for some useful features that are contained in the
signed spectrum, easier to compute exactly [43].

The aim of this program is to promote the use of techniques inspired from
field theory in the study of random tensors that are suitable to understand
large size phenomena, combined with rigorous methods from probability theory,
and compare further the limits of what can be learned about the original
problem from some approximations that nonetheless lead to simple and explicit
formulas. The underlying method is to write an effective action for quantum
fields, that symmetries and a correct decomposition of the degrees of freedom
can simplify. Schwinger-Dyson equations for two-point and higher order
correlators often reduce the apparent difficulty of the problem.

We plan our discussion as follows. We first start in Sec. 2 by deriving the
exact signed eigenvalue distribution in the presence of disorder at finite $N$
using a quartic fermionic theory. The section 3 uses a quartic action with a
combination of fermions and bosons to rewrite the absolute value of the
determinant, that we analyze in detail. In the large $N$ limit, we can
factorize the transverse (noting that the transverse part is independent of
the disorder) and parallel parts of the action. Writing the large $N$
Schwinger-Dyson equations together with a supersymmetry, we obtain the
limiting expression of the genuine spectral distribution. In Sec. 4, we
compare the asymptotics of the edge of the spectrum of both the signed and
genuine distributions, for the former using a Lefschetz thimble decomposition
and for the latter, asymptotic expressions of the large $N$ action, giving the
same as the signed case, extending thus the results of [1]. We conclude
briefly in Sec. 5 presenting some open questions. Three appendices close our
paper. App. A containing details on our simulations and App. B giving more
numerical support for the behaviour of the edge of the genuine spectrum, using
the finite $N$ formulas for the spectral density from [2]. Finally, relying on
the BRST invariance of the signed distribution, App. C explains that the total
number of eigenvectors is still exponential in $N$, despite their apparent
decrease as the variance $\beta$ of the noise grows.

##  2 Signed distribution

Let us quickly remind the setting of [1]. For the simplicity of our
discussions, we consider a tensor $C$ of order three ($p=3$). The eigenvector
distribution $\rho$ is given by 555Our convention is to sum over repeated
indices.

$\displaystyle\rho(v,C,\nu)=\sum_{i=1}^{n_{sol}}\delta^{(N)}(v-v^{i})=\absolutevalue{\det
M}\prod_{a=1}^{N}\delta(v_{a}-C_{abc}v_{b}v_{c}+\nu_{a})\,,$ (4)
$\displaystyle M_{ab}=\delta_{ab}-2C_{abc}v_{c}\,,$ (5)

where $\nu_{a}$ is a Gaussian random vector of variance $2\beta$ and
$\\{v^{i}\\}_{1\leq i\leq n_{\text{sol}}}$ solve the eigenvector equation (3).
Additionally, we take the tensor $C$ to be a fully-symmetric Gaussian random
variable with $\\#C=N(N+1)(N+2)/6$ independent components such that the
averaged distribution is

$\displaystyle\rho(v)=A^{-1}\int_{\mathbb{R}^{N}}\differential\nu\,e^{-\absolutevalue{\nu}^{2}/4\beta}\int_{\mathbb{R}^{\\#C}}\differential
C\,e^{-\alpha C^{2}}\rho(v,C,\nu)\,,$ (6) $\displaystyle C^{2}=\sum_{1\leq
a,b,c\leq N}C_{abc}C_{abc}\,,\quad\differential C=\prod_{1\leq a\leq b\leq
c\leq N}\differential C_{abc}\,,\quad\differential\nu=\prod_{1\leq i\leq
N}\differential\nu_{i}\quad\alpha,\beta>0,$ (7)

and the normalisation

$A=\int_{\mathbb{R}^{N}}\differential\nu\,e^{-\absolutevalue{\nu}^{2}/4\beta}\int_{\mathbb{R}^{\\#C}}\differential
C\,e^{-\alpha C^{2}}\,.$ (8)

Given the rescaling freedom in equations (1) and (3), the eigenvector
distribution $\rho(v)$ is related to the eigenvalue distribution
$\tilde{\rho}(z)$ by:

$\displaystyle\tilde{\rho}(z)\differential
z=\rho(\absolutevalue{v}=1/z)\absolutevalue{v}^{N-1}S_{N-1}\differential\absolutevalue{v}\,,$

with the spherical volume $S_{N-1}=2\pi^{N/2}/\Gamma(N/2)$, using also the
invariance of the measure of the tensor under $O(N)$.

We will first compute the signed spectral distribution following the steps of
[40] :

$\displaystyle\begin{split}\rho_{s}(v)&=\expectationvalue{\sum_{i=1}^{n_{sol}}\text{sign}(M(v^{i}))\delta^{(N)}(v-v^{i})}\\\
&=A^{-1}\int\differential\nu\differential
C\exp(-\absolutevalue{\nu}^{2}/4\beta-\alpha C^{2})\det
M(v)\prod_{a=1}^{N}\delta(v_{a}-C_{abc}v_{b}v_{c}+\nu_{a})\\\
&=A^{-1}(2\pi)^{-N}\int\differential
C\differential\psi\differential\bar{\psi}\differential\lambda\,e^{S}\end{split}$
(9)

where, after integrating over the deviation $\nu$, the action is

$\displaystyle S=-\alpha
C^{2}-\beta\lambda^{2}+i\lambda_{a}(v_{a}-C_{abc}v_{b}v_{c})+\bar{\psi}_{a}(\delta_{ab}-2C_{abc}v_{c})\psi_{b}\,,$
(10)

and $\lambda^{2}=\sum_{a=1}^{N}\lambda_{a}^{2}$. For any $N$-dimensional non-
singular matrix $M$, any vector $\phi\in\mathbb{R}^{N}$ and $N$-dimensional
Grassmann variables $\psi,\bar{\psi}$, we recall the useful formulas

$\int_{-\infty}^{+\infty}\differential\lambda\,e^{i\lambda
x}=2\pi\delta(x)\,,\quad\int\differential\phi\,e^{-\phi
M\phi}=\frac{1}{\sqrt{\det(M/\pi)}}\,,\quad\int\differential\bar{\psi}\differential\psi\,e^{\bar{\psi}M\psi}=\det
M\,.$ (11)

After the integration over $C$ that cancels the normalizing constant $A$, the
action becomes

$-\beta\lambda^{2}+i\lambda\cdot
v+\bar{\psi}\cdot\psi+\frac{1}{\alpha}\left(\frac{1}{6}\sum_{\sigma}\left(\frac{i}{2}\lambda_{\sigma_{a}}v_{\sigma_{b}}v_{\sigma_{c}}+\bar{\psi}_{\sigma_{a}}\psi_{\sigma_{b}}v_{\sigma_{c}}\right)\right)^{2}\,,$
(12)

where the last sum is over all permutations
$(\sigma_{a},\sigma_{b},\sigma_{c})$ of the indices $(a,b,c)$, since the
tensor $C$ is fully symmetric. Gathering the terms containing the Lagrange
multiplier $\lambda$, we find

$\displaystyle
S_{\lambda}=-\frac{1}{12\alpha}B_{ab}\lambda_{a}\lambda_{b}+i\lambda_{a}D_{a},$
(13) $\displaystyle
B_{ab}=(12\alpha\beta+\absolutevalue{v}^{4})\delta_{ab}+2\absolutevalue{v}^{2}v_{a}v_{b},$
(14) $\displaystyle D_{a}=v_{a}+\frac{1}{3\alpha}(\bar{\psi}_{a}\psi\cdot
v\absolutevalue{v}^{2}+\bar{\psi}\cdot
v\psi_{a}\absolutevalue{v}^{2}+\bar{\psi}\cdot v\psi\cdot vv_{a}),$ (15)

that can again be integrated exactly. It is convenient to project the
Grassmann variables on the unit vector $\hat{v}_{a}=v_{a}/\absolutevalue{v}$:
$\psi_{\parallel}=\psi\cdot\hat{v}$ and
$\psi_{\perp}=\psi-\hat{v}\psi_{\parallel}$. It is then easy to write the
inverse of the covariance matrix $B$

$\displaystyle\begin{split}B^{-1}_{ab}&=\frac{1}{12\alpha\beta+\absolutevalue{v}^{4}}\mathds{1}_{\perp
ab}+\frac{1}{12\alpha\beta+3\absolutevalue{v}^{4}}\ \mathds{1}_{\parallel
ab}\\\
&=\frac{1}{12\alpha\beta+\absolutevalue{v}^{4}}\delta_{ab}-\frac{2\absolutevalue{v}^{2}}{(12\alpha\beta+3\absolutevalue{v}^{4})(12\alpha\beta+\absolutevalue{v}^{4})}v_{a}v_{b}\,.\end{split}$
(16)

After integrating over $\lambda$, we obtain the effective action

$\displaystyle\begin{split}S_{\psi}=&-\frac{1}{2}\ln(3\absolutevalue{v}^{4}+12\alpha\beta)-\frac{N-1}{2}\ln(\absolutevalue{v}^{4}+12\alpha\beta)-\frac{\alpha\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+4\alpha\beta}\\\
&+\bar{\psi}_{\perp}\cdot\psi_{\perp}-\frac{\absolutevalue{v}^{4}-4\alpha\beta}{\absolutevalue{v}^{4}+4\alpha\beta}\bar{\psi}_{\parallel}\psi_{\parallel}-\frac{8\beta\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+12\alpha\beta}\bar{\psi}_{\perp}\cdot\psi_{\perp}\bar{\psi}_{\parallel}\psi_{\parallel}-\frac{\absolutevalue{v}^{2}}{6\alpha}(\bar{\psi}_{\perp}\cdot\psi_{\perp})^{2},\end{split}$
(17)

where the first two terms arise from the determinant of the covariance $B$.
Note that the quartic term mixing parallel and transverse components is new
compared to [40]. Doing next the integral over the parallel components gives

$\displaystyle\int\differential\bar{\psi}_{\parallel}\differential\psi_{\parallel}\,e^{S_{\psi}}=-\left(\frac{\absolutevalue{v}^{4}-4\alpha\beta}{\absolutevalue{v}^{4}+4\alpha\beta}+\frac{8\beta\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+12\alpha\beta}\bar{\psi}_{\perp}\cdot\psi_{\perp}\right)e^{S_{\psi_{\perp}}(k=1)+S_{0}},$
(18) $\displaystyle
S_{\psi_{\perp}}(k)=k\bar{\psi}_{\perp}\cdot\psi_{\perp}-\frac{\absolutevalue{v}^{2}}{6\alpha}(\bar{\psi}_{\perp}\cdot\psi_{\perp})^{2},$
(19) $\displaystyle
S_{0}=-\frac{1}{2}\ln(3\absolutevalue{v}^{4}+12\alpha\beta)-\frac{N-1}{2}\ln(\absolutevalue{v}^{4}+12\alpha\beta)-\frac{\alpha\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+4\alpha\beta},$
(20)

and the integral over the transverse components can also be done as follows
[40]

$\displaystyle\begin{split}Z_{\perp}(k)=&\int\differential\bar{\psi}_{\perp}\differential\psi_{\perp}\,e^{S_{\psi_{\perp}}(k)}\\\
=&\sum_{n\geq
0}\frac{1}{n!}\left(-\frac{\absolutevalue{v}^{2}}{6\alpha}\right)^{n}\left(\frac{\partial}{\partial
k}\right)^{2n}\underbrace{\int\differential\bar{\psi}_{\perp}\differential\psi_{\perp}\,e^{k\bar{\psi}_{\perp}\cdot\psi_{\perp}}}_{=k^{N-1}}\\\
=&\left(\frac{3\alpha}{2\absolutevalue{v}^{2}}\right)^{1-N/2}kU\left(1-\frac{N}{2},\frac{3}{2},\frac{3\alpha
k^{2}}{2\absolutevalue{v}^{2}}\right)\,.\end{split}$ (21)

Collecting all the integration constants, the signed spectral distribution
writes as

$\displaystyle\rho_{s}(v)=-\left(\frac{\sqrt{12\pi\alpha}}{2\pi}\right)^{N}\left(\frac{\absolutevalue{v}^{4}-4\alpha\beta}{\absolutevalue{v}^{4}+4\alpha\beta}+\frac{8\beta\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+12\alpha\beta}\partial_{k}\right)Z_{\perp}(k)\rvert_{k=1}e^{S_{0}},$
(22) $\displaystyle
Z_{\perp}(k)=k\left(\frac{3\alpha}{2\absolutevalue{v}^{2}}\right)^{1-N/2}U\left(1-\frac{N}{2},\frac{3}{2},\frac{3\alpha
k^{2}}{2\absolutevalue{v}^{2}}\right)\,,$ (23)

with $U$ being a confluent hypergeometric function of the second kind.

For the distribution of the absolute value $\absolutevalue{v}$, one looks at
$\rho^{\text{size}}_{s}(v)=\rho_{s}(\absolutevalue{v})\absolutevalue{v}^{N-1}S_{N-1}$,

$\rho^{\text{size}}_{s}(v)=-\frac{3\alpha
2^{N/2}v^{2N-3}}{\Gamma(N/2)}e^{S_{0}}\left(\frac{\absolutevalue{v}^{4}-4\alpha\beta}{\absolutevalue{v}^{4}+4\alpha\beta}+\frac{8\beta\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+12\alpha\beta}\partial_{k}\right)kU\left(1-\frac{N}{2},\frac{3}{2},\frac{3\alpha
k^{2}}{2\absolutevalue{v}^{2}}\right)\bigg{|}_{k=1}\,.$ (24)

We compare in Fig. 1 the result of our analytic computation to Monte Carlo
simulations (cf. Appendix A for the details on the simulations), for which we
see excellent agreement. We also noted numerically that the integral of
$\rho^{\text{size}}_{s}$ is independent of $\beta$ (see App. C for the
explanation of this independence.).

Figure 1: Comparison of the analytic expression (in continuous blue) of the
signed eigenvalue distributions for $\beta=0.1,0.069/N^{2},0.001/N^{2}$ with
the numerical simulation (in dotted red, with associated numerical errors),
for $N=12$ and 10000 iterations).

##  3 Genuine distribution

In this section, we will look at the genuine distribution eq. (4), following
[41]. We will first obtain the field theory corresponding to the genuine
distribution in Section 3.1. Then in Section 3.2 we will compute the field
theory partition function with approximations that we detail. We will
crosscheck our result with Monte Carlo simulations finding excellent
agreement.

###  3.1 Corresponding field theory and its symmetries

The absolute value of the determinant in (4) can be written as

$\absolutevalue{\det M}=\lim_{\epsilon\to 0^{+}}\frac{\det(M^{2}+\epsilon
I)}{\sqrt{\det(M^{2}+\epsilon I)}},$ (25)

and introducing the scalar $(\phi_{a})_{1\leq a\leq N}$ and Grassmann
variables $(\psi_{a},\bar{\psi}_{a},\varphi_{a},\bar{\varphi}_{a})_{1\leq
a\leq N}$, respectively for the denominator and the numerator, and after a
Hubbard-Stratonovich procedure that decomposes the quadratic terms in $C$ into
linear ones with the introduction of an extra Gaussian vector
$(\sigma_{a})_{1\leq a\leq N}$, we obtain the following integral

$\rho(v,\alpha,\beta)=\lim_{\epsilon\to
0^{+}}\frac{(-1)^{N}}{A(2\pi^{2})^{N}}\int\differential
C\differential\lambda\differential\phi\differential\sigma\differential\bar{\psi}\differential\psi\differential\bar{\varphi}\differential\varphi\,e^{S},$
(26) $S=-\alpha
C^{2}-\beta\lambda^{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\,,$
(27)

where the quadratic term for the Lagrange multiplier $\lambda$ comes from the
integration over the Gaussian noise.

The integration over $C$ gives, like in [41]:

$\displaystyle\begin{split}\delta
S_{C}&=-\frac{\absolutevalue{v}^{4}}{12\alpha}\lambda_{a}B_{ab}\lambda_{b}-i\lambda_{a}(D_{a}+\tilde{D}_{a})+E_{1}+E_{2}+E_{3}\,,\\\
B_{ab}&=\delta_{ab}\left(1+\frac{12\alpha\beta}{\absolutevalue{v}^{4}}\right)+2\hat{v}_{a}\hat{v}_{b}=\mathds{1}_{\perp
ab}\left(1+\frac{12\alpha\beta}{\absolutevalue{v}^{4}}\right)+\mathds{1}_{\parallel
ab}\left(3+\frac{12\alpha\beta}{\absolutevalue{v}^{4}}\right)\,,\\\
D_{a}&=D_{\parallel}\hat{v}_{a}+D_{\perp a}\,,\\\
\tilde{D}_{a}&=\tilde{D}_{\parallel}\hat{v}_{a}+\tilde{D}_{\perp a}\,,\\\
E_{1}&=\frac{1}{\alpha}\left(\frac{1}{6}\sum_{s}(\bar{\psi}_{s_{a}}\varphi_{s_{a}}v_{s_{a}}+\bar{\varphi}_{s_{a}}\psi_{s_{a}}v_{s_{a}}\right)^{2}\,,\\\
E_{2}&=-\frac{4}{\alpha}\left(\frac{1}{6}\sum_{s}v_{s_{a}}\sigma_{s_{a}}\phi_{s_{a}}\right)^{2}\,,\\\
E_{3}&=\frac{4i}{\alpha}(v_{a}\bar{\psi}_{b}\varphi_{c}+v_{a}\bar{\varphi}_{b}\psi_{c})\left(\frac{1}{6}\sum_{s}v_{s_{a}}\sigma_{s_{a}}\phi_{s_{a}}\right)\,,\end{split}$
(28)

where $\hat{v}=v/\absolutevalue{v}$. The summation over $s$ encodes a sum over
the symmetric permutations of all indices $a,b,c$. The parallel and transverse
components to $v$ of $D,\tilde{D}$ are given by

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

After integrating over $\lambda$, we obtain:

$\displaystyle\begin{split}\delta
S_{\lambda}=&-\frac{\alpha\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+4\alpha\beta}+\frac{2\alpha\absolutevalue{v}}{\absolutevalue{v}^{4}+4\alpha\beta}(D_{\parallel}+\tilde{D}_{\parallel})\\\
&-\left[\frac{3\alpha}{\absolutevalue{v}^{4}+12\alpha\beta}(D_{\perp}\cdot
D_{\perp}+\tilde{D}_{\perp}\cdot\tilde{D}_{\perp}+2D_{\perp}\cdot\tilde{D}_{\perp})+\frac{\alpha}{\absolutevalue{v}^{4}+4\alpha\beta}(D_{\parallel}^{2}+\tilde{D}_{\parallel}^{2}+2D_{\parallel}\tilde{D}_{\parallel})\right]\,.\end{split}$
(30)

Putting it all together, the density rewrites as a zero-dimensional quartic
field theory of bosons and fermions:

$\displaystyle\rho(v,\alpha,\beta)$ $\displaystyle=\lim_{\epsilon\to
0^{+}}\frac{3^{(N-1)/2}\alpha^{N/2}}{\pi^{N/2}}\frac{1}{(\absolutevalue{v}^{4}+4\alpha\beta)^{1/2}(\absolutevalue{v}^{4}+12\alpha\beta)^{(N-1)/2}}\exp(-\frac{\alpha\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+4\alpha\beta})Z_{FT}\,,$
(31) $\displaystyle Z_{FT}=\frac{(-1)^{N}}{\pi^{N}}\int\differential
B\differential F\exp(K_{\parallel}+\delta V+S_{\perp})\,,$ (32)

where $\int dBdF$ stands for the integration over all the bosonic and
fermionic variables. We have combined the terms quadratic in parallel
components into $K_{\parallel}$:

$K_{\parallel}=\begin{pmatrix}\bar{\psi}_{\parallel}&\bar{\varphi}_{\parallel}\end{pmatrix}\begin{pmatrix}\epsilon&\gamma_{1}\\\
\gamma_{1}&-1\end{pmatrix}\begin{pmatrix}\psi_{\parallel}\\\
\varphi_{\parallel}\end{pmatrix}+\begin{pmatrix}\phi_{\parallel}&\sigma_{\parallel}\end{pmatrix}\begin{pmatrix}-\epsilon&\gamma_{1}i\\\
\gamma_{1}i&-1\end{pmatrix}\begin{pmatrix}\phi_{\parallel}\\\
\sigma_{\parallel}\end{pmatrix}$ (33)

and the mixing or quartic terms in the parallel components into $\delta V$:

$\displaystyle\begin{split}\delta
V=&-\gamma_{2}\bar{\psi}_{\parallel}\psi_{\parallel}\bar{\varphi}_{\parallel}\varphi_{\parallel}+2\gamma_{2}\left(i(\bar{\psi}_{\parallel}\varphi_{\parallel}+\bar{\varphi}_{\parallel}\psi_{\parallel})\sigma_{\parallel}\phi_{\parallel}-(\sigma_{\parallel}\phi_{\parallel})^{2}\right)\\\
&+2\gamma_{3}\left(i(\bar{\psi}_{\perp}\varphi_{\parallel}+\bar{\psi}_{\parallel}\varphi_{\perp}+\bar{\varphi}_{\perp}\psi_{\parallel}+\bar{\varphi}_{\parallel}\psi_{\perp})\cdot(\sigma_{\parallel}\phi_{\perp}+\sigma_{\perp}\phi_{\parallel})-(\sigma_{\parallel}\phi_{\perp}+\sigma_{\perp}\phi_{\parallel})^{2}\right)\\\
&+\gamma_{3}\left(\bar{\psi}_{\perp}\cdot\bar{\varphi}_{\perp}\psi_{\parallel}\varphi_{\parallel}+\psi_{\perp}\cdot\varphi_{\perp}\bar{\psi}_{\parallel}\bar{\varphi}_{\parallel}-\bar{\psi}_{\perp}\cdot\psi_{\perp}\bar{\varphi}_{\parallel}\varphi_{\parallel}-\bar{\psi}_{\perp}\cdot\varphi_{\perp}\bar{\psi}_{\parallel}\varphi_{\parallel}-\bar{\varphi}_{\perp}\cdot\varphi_{\perp}\bar{\psi}_{\parallel}\psi_{\parallel}-\bar{\varphi}_{\perp}\cdot\psi_{\perp}\bar{\varphi}_{\parallel}\psi_{\parallel}\right),\end{split}$
(34)

introducing the constants:

$\gamma_{1}=\frac{\absolutevalue{v}^{4}-4\alpha\beta}{\absolutevalue{v}^{4}+4\alpha\beta}\,,\quad\gamma_{2}=\frac{8\beta\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+4\alpha\beta}\,,\quad\gamma_{3}=\frac{8\beta\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+12\alpha\beta}\,.$
(35)

The terms involving only transverse components are collected in $S_{\perp}$
(notice that it is independent of $\beta$), with the quadratic terms into
$K_{\perp}$, the quartic bosonic, fermionic and sharing bosonic and fermionic
transverse components respectively into $V_{B}$, $V_{F}$ and $V_{BF}$:

$\displaystyle\begin{split}&S_{\perp}=K_{\perp}+V_{F}+V_{B}+V_{BF}\,,\\\ &\ \
\
K_{\perp}=\begin{pmatrix}\bar{\psi}_{\perp}&\bar{\varphi}_{\perp}\end{pmatrix}\begin{pmatrix}\epsilon&-1\\\
-1&-1\end{pmatrix}\begin{pmatrix}\psi_{\perp}\\\
\varphi_{\perp}\end{pmatrix}+\begin{pmatrix}\phi_{\perp}&\sigma_{\perp}\end{pmatrix}\begin{pmatrix}-\epsilon&-i\\\
-i&-1\end{pmatrix}\begin{pmatrix}\phi_{\perp}\\\
\sigma_{\perp}\end{pmatrix}\,,\\\ &\ \ \
V_{B}=-\frac{2\absolutevalue{v}^{2}}{3\alpha}\left(\sigma_{\perp}^{2}\phi_{\perp}^{2}+(\sigma_{\perp}\cdot\phi_{\perp})^{2}\right)\,,\\\
&\ \ \
V_{F}=-\frac{\absolutevalue{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)\,,\\\
&\ \ \
V_{BF}=\frac{2i\absolutevalue{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}$
(36)

#### Symmetries.

Two symmetries are important to understand the large $N$ limit of the system:
the global $O(1)\times O(N-1)$ symmetry, leading to the factorization between
the parallel and transverse components action in that limit, and the graded
symmetry exchanging the bosonic and the fermionic variables,

$\displaystyle\delta_{1}\psi_{a}=-\phi_{a}\,,\quad\delta_{1}\phi_{a}=\frac{1}{2}\bar{\psi}_{a}\,,\quad\delta_{1}\varphi_{a}=i\sigma_{a}\,,\quad\delta_{1}\sigma_{a}=\frac{i}{2}\bar{\varphi}_{a}\,,$
(37)
$\displaystyle\delta_{2}\bar{\psi}_{a}=\phi_{a}\,,\quad\delta_{2}\phi_{a}=\frac{1}{2}\psi_{a}\,,\quad\delta_{2}\bar{\varphi}_{a}=-i\sigma_{a}\,,\quad\delta_{2}\sigma_{a}=\frac{i}{2}\varphi_{a}\,.$
(38)

Assumed unbroken 666We have verified that those equalities hold for $N=1$. It
is also not a proper supersymmetry since the number of bosons differ from the
number of fermions and there is no space dependence [44]., this symmetry
implies the following relations between the expectation values of the
fermionic and the bosonic fields, which will be assumed in Section 3.2:

$\displaystyle\begin{split}&\expectationvalue{\bar{\psi}_{\perp}\cdot\psi_{\perp}}=2\expectationvalue{\phi_{\perp}\cdot\phi_{\perp}},\quad\expectationvalue{\bar{\varphi}_{\perp}\cdot\psi_{\perp}}=-2i\expectationvalue{\phi_{\perp}\cdot\sigma_{\perp}},\\\
&\expectationvalue{\bar{\varphi}_{\perp}\cdot\varphi_{\perp}}=-2\expectationvalue{\sigma_{\perp}\cdot\sigma_{\perp}},\quad\expectationvalue{\bar{\psi}_{\perp}\cdot\varphi_{\perp}}=-2i\expectationvalue{\sigma_{\perp}\cdot\phi_{\perp}}.\end{split}$
(39)

###  3.2 Computation

In this section, we will obtain the genuine distribution by explicitly
computing the partition function (32) by an approximation, which is expected
to be valid for large $N$. As shown in Section 3.1, the field theory can be
regarded as an interacting system between the parallel and transverse sectors.
For large $N$ the transverse sector dominates, and therefore we can neglect
the parallel components to determine the dynamics of the transverse ones by
their self-interaction. In a second step, the dynamics of the parallel sector
can be computed in the background of the transverse sector, where the
interaction between the parallel and transverse sectors is considered. More
explicitly, in the large $N$ limit, the partition function factorizes into the
form,

$Z_{FT}\sim Z_{FT}^{\perp}\cdot Z_{FT}^{\parallel}(Q^{*}),$ (40)

where $Z_{FT}^{\perp}$ denotes the partition function determined solely by the
transverse sector, and $Z_{FT}^{\parallel}(Q^{*})$ denotes that of the
parallel sector with $Q^{*}$ being the background of the transverse sector.

An important fact in our approximation is that $Z_{FT}^{\perp}$ does not
depend on $\beta$, while the $\beta$-dependence appears in
$Z_{FT}^{\parallel}(Q^{*})$. Therefore we can obtain $Z_{FT}^{\perp}$ using
some former results without the noise. We will first review the Schwinger-
Dyson (SD) method, which was employed in [45]. This determines
$Z_{FT}^{\perp}$ as well as the background $Q^{*}$ in the large $N$ limit. We
will also use the exact formula given in [2] for better approximations at
finite $N$.

#### Transverse sector.

The transverse sector is defined by the partition function,

$Z_{FT}^{\perp}=\frac{(-1)^{N-1}}{\pi^{N-1}}\int\differential
B_{\perp}\differential F_{\perp}e^{S_{\perp}},$ (41)

where $S_{\perp}$ is given in (36). Note that it does not contain $\beta$.
Neglecting the subleading terms in the potentials $V_{F}$ and $V_{BF}$ and
assuming that the $O(N-1)$ symmetry remains unbroken in the large $N$ limit,
we can write an effective transverse action in terms of the two-point
correlators:

$\displaystyle S_{\perp}(Q)$ $\displaystyle=\frac{N-1}{2}\left(\epsilon
Q_{11}-Q_{22}-2Q_{12}-x(Q_{12}^{2}+Q_{11}Q_{22})-\log(-\det
Q)\right)+\mathcal{O}(1)\,,$ (42)
$Q=\frac{1}{N-1}\begin{pmatrix}\expectationvalue{\bar{\psi}_{\perp}\cdot\psi_{\perp}}&\expectationvalue{\bar{\psi}_{\perp}\cdot\varphi_{\perp}}\\\
\expectationvalue{\bar{\varphi}_{\perp}\cdot\psi_{\perp}}&\expectationvalue{\bar{\varphi}_{\perp}\cdot\varphi_{\perp}}\end{pmatrix}\,,\quad
x=\frac{(N-1)\absolutevalue{v}^{2}}{3\alpha}\,,$ (43)

where we used that $Q_{12}=Q_{21}$,777This can be derived from (39). and
defined the variable $x$. The logarithmic term arises as a Jacobian of the
change of variables from the fermionic fields to the bilinear $Q$ (see [46],
section B4) 888We introduced a sign since $\det Q^{*}<0$..

Notice that bosons and fermions have the same leading order action, so the
saddle point is also supersymmetric, as was already observed in [47]. Similar
graded symmetries were used in earlier treatments of systems with rough
landscapes, when approximating the absolute value of the determinant by the
determinant itself. It is now understood that this supersymmetry holds if
there is at most one step replica symmetry breaking, but it is broken if
higher, when the number of metastable states is very sensitive to external
perturbations, like in the Sherrington-Kirkpatrick model [48].

#### Large-$N$ Solution.

Since the effective transverse action is independent of $\beta$ and
contributions from the parallel sector are subleading in $N$, we can use the
saddle point analysis of [45]. At large $N$, the leading order saddle point of
$S_{\perp}$ that reproduces the free theory at $v=0$ is singular at the
critical point $x=1/4$, and expanded at small $\epsilon$, is given by:

  * •

$0<x<1/4$:

$\displaystyle\begin{split}&Q_{11}^{*}=\frac{-\sqrt{1-4x}+1}{2x\sqrt{1-4x}}+\mathcal{O}(\epsilon)\,,\\\
&Q_{12}^{*}=\frac{1-\sqrt{1-4x}-4x}{2x\sqrt{1-4x}}+\mathcal{O}(\epsilon)\,,\\\
&Q_{22}^{*}=\mathcal{O}(\epsilon^{2})\,,\end{split}$ (44)

  * •

$x>1/4$:

$\displaystyle\begin{split}&Q_{11}^{*}=\frac{\sqrt{4x-1}}{2\sqrt{\epsilon}x}-\frac{1}{2x}+\mathcal{O}(\sqrt{\epsilon})\,,\\\
&Q_{12}^{*}=-\frac{1}{2x}+\frac{\sqrt{\epsilon}}{2x\sqrt{4x-1}}+\mathcal{O}(\epsilon^{3/2})\,,\\\
&Q_{22}^{*}=-\frac{\sqrt{\epsilon}\sqrt{4x-1}}{2x}+\frac{\epsilon}{2x}+\mathcal{O}(\epsilon^{3/2})\,.\end{split}$
(45)

When $\epsilon\to 0^{+}$, the effective action evaluates to:

$S_{\perp}(Q^{*})=\begin{cases}\frac{(N-1)}{4}\left(2+\log
16+\frac{1-\sqrt{1-4x}}{x}-4\log\frac{1-\sqrt{1-4x}}{x}\right)\,,&0<x<1/4\,,\\\
\frac{(N-1)}{4x}\left(1+2x+2x\log x\right)\,,&x>1/4\,.\end{cases}$ (46)

Thus the partition function in the large-$N$ is given by

$Z_{FT}^{\perp}\sim e^{S_{\perp}(Q^{*})-(N-1)},$ (47)

where the subtraction serves to ensure the normalization $Z_{FT}^{\perp}=1$
for $|v|=0\ (x=0)$.

The random matrix point of view for this action is as arising from eigenvalues
of the Hessian of the spherical $p$-spin at a given stationary point (in fact
the minima, that are dominating the critical points of any fixed degree at
that level) [2], where the two $x$ regions of (46) are lying respectively
outside or inside the bulk of the large $N$ spectrum for the first and second
line of eq. (46) (see [32] for a pedagogical account).

#### Parallel sector.

Now we want to compute the remaining part of the partition function (32):

$Z_{FT}^{\parallel}(Q^{*})=-\frac{1}{\pi}\int\differential
B_{\parallel}\differential F_{\parallel}\,e^{K_{\parallel}+\delta V},$ (48)

where all the transverse fields are replaced by the two-field correlators
$Q^{*}$ obtained above. To deal with the parallel components, we can do a
Hubbard-Stratonovich transformation with a Gaussian variable $g$ to split the
quartic terms of the first line of eq. (34)

$\frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty}\differential
ge^{-g^{2}+\sqrt{2\gamma_{2}}(\bar{\psi}_{\parallel}\varphi_{\parallel}+\bar{\varphi}_{\parallel}\psi_{\parallel}+2i\sigma_{\parallel}\phi_{\parallel})g}$
(49)

and integrate out the now Gaussian parallel components of the fermionic and
bosonic fields (the last ones bringing each a factor $\sqrt{\pi}$).

In the large $N$ limit, the first term of the second line in $\delta V$ is
subleading, mixing parallel and transverse components of the bosonic and
fermionic variables, as well as the first two of the third line. In that
limit, the partition function factorizes since the resulting determinant
factor after integrating the Gaussian parallel components is subleading in $N$
with respect to the purely transverse action (42), leading to

$Z^{\parallel}_{FT}(Q^{*})=\frac{-1}{\sqrt{\pi}}\int_{-\infty}^{+\infty}\differential
ge^{-g^{2}}\lim_{\epsilon\rightarrow
0^{+}}\left[\big{(}\epsilon-\gamma_{3}(N-1)Q^{*}_{22}\big{)}\big{(}1+\gamma_{3}(N-1)Q^{*}_{11}\big{)}+\left(\gamma_{1}-(N-1)\gamma_{3}Q^{*}_{12}+\sqrt{2\gamma_{2}}g\right)^{2}\right]^{1/2}\,.$
(50)

We evaluated this integral numerically with Mathematica, in the limits
$N\to\infty$ and $\epsilon\to 0^{+}$ (in that order).

#### Finite $N$ case.

For finite $N$, we can take advantage of the exact formula given in [2].999See
[45] for a field theoretical derivation. From (31), (40), the independence of
$Z_{FT}^{\perp}$ on $\beta$, and $Z_{FT}^{\parallel}(Q^{*})=1$ for $\beta=0$,
we find

$\rho^{\rm
size}(v,\alpha,\beta)\sim\left(\frac{|v|^{4}}{|v|^{4}+4\alpha\beta}\right)^{\frac{1}{2}}\left(\frac{|v|^{4}}{|v|^{4}+12\alpha\beta}\right)^{\frac{N-1}{2}}\exp\left(-\frac{\alpha|v|^{2}}{|v|^{4}+4\alpha\beta}+\frac{\alpha}{|v|^{2}}\right)Z_{FT}^{\parallel}(Q^{*})\,\rho^{\rm
size}(v,\alpha,\beta=0),$ (51)

where $\rho^{\rm size}(v,\alpha,\beta=0)$ is the genuine distribution without
the noise. By properly changing the variables in [2] and considering that
their normalization of the Gaussian tensor corresponds to $\alpha=1/2$, we
obtain

$\rho^{\rm
size}(|v|,\alpha=1/2,\beta=0)=2\sqrt{N}|v|^{-2}(p-1)^{(N-1)/2}\exp(-(p-2)/(4(p-1)|v|^{2}))\rho_{N}(x),$
(52)

where $\rho_{N}(\cdot)$ is given in Section 7.2 of [2],
$x=-\sqrt{p/(2(p-1))}/(\sqrt{N}|v|)$, and $p=3$ (see Appendix B for more
explanation). Note that the formula (51) still depends on the SD approximation
for $Q^{*}$, which is valid in the large-$N$ limit.

#### Summary.

We have crosschecked the expression (51) with Monte Carlo simulations (see
Appendix A for the details). The agreement is shown in Fig. 2 for different
values of $\beta$. On the middle plot, we see the emergence of the outlier and
on the rightmost plot, its merging with the trivial eigenvector. One also
notices that as $\beta$ increases, the spectrum seems to concentrates at
smaller values of eigenvectors. In Appendix C, we show that although the
number of eigenvectors decreases with $\beta$, it still grows exponentially
with $N$, being bounded from below by the integral of the signed density,
itself independent of $\beta$.

Figure 2: Genuine eigenvector distributions for
$\beta=0.1,0.069/N^{2},0.001/N^{2}$ (from left to right), the analytic
expression (in continuous blue) and Monte Carlo simulation (in dotted red)
with the associated numerical error, for $N=10$ and 10000 iterations.

##  4 The largest eigenvalue

###  4.1 Signed distribution

In order to determine the location of the smallest non-trivial eigenvector
$v_{\text{edge}}$, we look at the dominant behaviour in the exponential of the
distribution $\rho^{\text{size}}_{s}(v)$. We first analyse $Z_{\perp}(k)$,
following [43]. One can express the confluent hypergeometric function $U$ with
Hermite polynomials, introducing $Ny^{2}=3\alpha/(2\absolutevalue{v}^{2})$,

$\displaystyle U\left(1-\frac{N}{2},\frac{3}{2},\frac{3\alpha
k^{2}}{2\absolutevalue{v}^{2}}\right)=\frac{\Gamma(\frac{N}{2})\Gamma(\frac{N+1}{2})}{\Gamma(N)}\frac{1}{\sqrt{\pi}}\frac{1}{\sqrt{N}yk}H_{N-1}(\sqrt{N}yk)\,.$
(53)

Those polynomials obey the relations

$\displaystyle H^{\prime}_{N}(y)=2NH_{N-1}(y)\,,\quad H_{N}(y)=\frac{N!}{2\pi
i}\oint_{\mathcal{C}}\differential t\frac{e^{2ty-t^{2}}}{t^{N+1}}\,,$ (54)

where the contour $\mathcal{C}$ winds counterclockwise around the origin. The
signed density can then be expressed as

$\displaystyle\begin{split}\rho^{\text{size}}_{s}(v)=&\frac{3\alpha
2^{N/2}v^{2N-3}}{\sqrt{\pi}}\frac{\Gamma(\frac{N+1}{2})}{\Gamma(N)}e^{S_{0}}\left(\frac{\absolutevalue{v}^{4}-4\alpha\beta}{\absolutevalue{v}^{4}+4\alpha\beta}+\frac{8\beta\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+12\alpha\beta}\partial_{k}\right)\frac{1}{\sqrt{N}y}H_{N-1}(\sqrt{N}yk)\bigg{|}_{k=1}\\\
\sim&\frac{3\alpha
2^{N/2}v^{2N-3}}{\sqrt{\pi}}\frac{\Gamma(\frac{N+1}{2})}{\Gamma(N)}e^{S_{0}}\left(\frac{\absolutevalue{v}^{4}-4\alpha\beta}{\absolutevalue{v}^{4}+4\alpha\beta}\frac{N^{-N/2}}{y}+\frac{16\beta\absolutevalue{v}^{2}}{\absolutevalue{v}^{4}+12\alpha\beta}N^{(2-N)/2}\right)\frac{(N-1)!}{2\pi
i}\oint\differential\tau e^{Nf(y,\tau)}\end{split}$ (55)

in the large $N$ limit, using the variable $\tau=t/\sqrt{N}$ and

$f(y,\tau)=2\tau y-\tau^{2}-\log\tau\,.$ (56)

The saddle points of $f(y,\tau)$ are located at

$\tau_{\pm}=\frac{1}{2}(y\pm\sqrt{y^{2}-2})$ (57)

with a critical point at $y_{c}=\sqrt{2}$. The contour integral can be
deformed into an integral over the Lefschetz thimbles101010See for instance
Section 3 of [49] for the steepest descent method in terms of Lefschetz
thimbles. originating from the different saddle points $\tau_{*}$, obeying the
flow equation,

$\derivative{\tau(t)}{t}=-\overline{\partialderivative{f(y,\tau)}{\tau}}\,,\quad\lim_{t\to-\infty}\tau(t)=\tau_{*},$
(58)

and the saddles that actually contribute are those whose upward flow,
determined by

$\derivative{\tau(t)}{t}=\overline{\partialderivative{f(y,\tau)}{\tau}}\,,\quad\lim_{t\to-\infty}\tau(t)=\tau_{*}\,,$
(59)

intersect the original contour $\mathcal{C}$. As in [43], only the $\tau_{-}$
saddle point contributes in the $y>y_{c}$ region, that corresponds to the
small $v$ region, while both $\tau_{\pm}$ are relevant in the region
$y<y_{c}$.

In order to track the displacement of the edge of the density, we gather the
leading factors in $N$, we have for $y>y_{c}$

$\displaystyle\begin{split}\rho^{\text{size}}_{s}(v)\sim&\exp(N\left[f(y,\tau_{-})-\frac{1}{2}-\frac{1}{2}\log(1+\frac{16\bar{\beta}y^{4}}{3\alpha})-\frac{6\alpha^{2}y^{2}}{9\alpha^{2}+16\alpha\bar{\beta}y^{4}}\right]),\end{split}$
(60)

or also taking into account the contribution of the second saddle point for
$y<y_{c}$,

$\displaystyle\begin{split}h(y)&=\lim_{N\to\infty}\frac{1}{N}\log\rho^{\text{size}}_{s}(y)\\\
&=\real
f(y,\tau_{-})-\frac{1}{2}-\frac{1}{2}\log(1+\frac{16\bar{\beta}y^{4}}{3\alpha})-\frac{6\alpha^{2}y^{2}}{9\alpha^{2}+16\alpha\bar{\beta}y^{4}}.\end{split}$
(61)

We observe three different regimes in the roots of the exponent of
distribution, represented on Fig 3, and the corresponding plots of $h(y)$ on
Fig 4, for $\bar{\beta}=0.003,0.045,0.095$ with $\alpha=1/2$.

Figure 3: The different real roots of the asymptotics of the signed density as
a function of the rescaled noise variance $\bar{\beta}$ ($\alpha=1/2$).

Figure 4: The asymptotics of the logarithm of the signed density as
$\bar{\beta}$ increases, from left to right $\bar{\beta}=0.003,0.045,0.095$
($\alpha=1/2$).

There is a first branching point at the location of the critical point
$y_{c}=\sqrt{2}$, where $\bar{\beta}_{1}$ solves

$0=1-\frac{12\alpha}{9\alpha+64\bar{\beta}_{1}}+\frac{1}{2}\log
2-\frac{1}{2}\log\left(1+\frac{64\bar{\beta}_{1}}{3\alpha}\right)\,,\quad\bar{\beta}_{1}=3\alpha/64\,,$
(62)

giving for $\alpha=1/2$, $\bar{\beta}_{1}=3/128\approx 0.023$. An expansion
near the critical point indicates the following behaviour for
$\bar{\beta}<\bar{\beta}_{1}$

$\absolutevalue{y-y_{c}}\sim\frac{128\sqrt{2}}{9}\absolutevalue{\bar{\beta}-\bar{\beta}_{1}}^{2}\,,\quad
y>y_{c}\,,$ (63)

while on the other side $\bar{\beta}>\bar{\beta}_{1}$, we have two different
exponents

$\absolutevalue{y-y_{c}}\sim\begin{cases}\frac{512\sqrt{2}}{9}\absolutevalue{\bar{\beta}-\bar{\beta}_{1}}^{2}\,,&y>y_{c}\,,\\\
\frac{2^{7/4}}{\sqrt{3}}\absolutevalue{\bar{\beta}-\bar{\beta}_{1}}\,,&y<y_{c}\,.\end{cases}$
(64)

The second branching point is located where the two largest real eigenvalues,
corresponding to the two maxima on the right in the middle picture of Fig. 4,
merge and become complex. Here the locations of the two maxima can be obtained
by solving the extremal condition,

$\frac{\partial h(y)}{\partial y}=0.$ (65)

The solution of this extremal condition satisfying $y>\sqrt{2}$ and
$\bar{\beta}>0$ is given by

$\bar{\beta}=\frac{3\alpha(-3+2y^{2}+2y\sqrt{-2+y^{2}})}{16y^{4}}.$ (66)

Then it is straightforward to see that the $\bar{\beta}$ in (66) takes the
maximum value $\bar{\beta}_{2}$ at $y=\sqrt{3}$. In other words, for
$\bar{\beta}>\bar{\beta}_{2}$, there are no solutions to (65) satisfying
$y>\sqrt{2}$. For $\alpha=1/2$, $\bar{\beta}_{2}=(3+2\sqrt{3})/96\approx
0.067334$. One can also check that $h(y)=0$ for $\bar{\beta}=\beta_{2},\
y=\sqrt{3}$, implying that the two maxima indeed satisfy $h(y)=0$ as well.

###  4.2 Genuine distribution

The leading order genuine density is given by

$\displaystyle\begin{split}\rho^{\text{size}}(v)=&\frac{2\cdot
3^{(N-1)/2}\alpha^{N/2}}{\Gamma(N/2)}\frac{v^{N-1}}{(v^{4}+4\alpha\beta)^{1/2}(v^{4}+12\alpha\beta)^{(N-1)/2}}\exp(-\frac{\alpha
v^{2}}{v^{4}+4\alpha\beta})\\\ &\exp(\frac{N-1}{4}\left(2+\log
16+\frac{1-\sqrt{1-4x}}{x}-4\log(1-\sqrt{1-4x})+4\log
x\right)-(N-1))\,,\end{split}$ (67)

with $x=(N-1)v^{2}/3\alpha$ and in the second line, the approximation of the
large $N$ field theory partition function was taken, subtracting
$S_{eff}(0)=N-1$ to ensure that it is appropriately normalized. 111111We also
multiplied by a factor $\pi^{N}$ to account for the two scalar fields $\phi$
and $\sigma$. The large $N$ behaviour of the genuine distribution is then

$\displaystyle\begin{split}\log\rho^{\text{size}}(x)=&\frac{N}{4}\bigg{(}-\frac{12\alpha
x}{4\bar{\beta}+9\alpha x^{2}}+\log
16+\frac{1-\sqrt{1-4x}}{x}-4\log(1-\sqrt{1-4x})+4\log x\\\ &+2\log 9\alpha
x-2\log(3\alpha(4\bar{\beta}+3\alpha
x^{2}))\bigg{)}+\mathcal{O}(1)\,,\end{split}$ (68)

of which the real part is the same as eq. (61).

##  5 Conclusion

To summarize, extending the series of works [45, 41, 1, 40, 43], we have
obtained the exact signed distribution and the large $N$ genuine distribution
of the Z-eigenvalues of Gaussian real random fully-symmetric tensors of order
3, in the presence of a Gaussian noise, using respectively, a quartic
fermionic theory and an interacting bosonic and fermionic theory, focusing on
the location of the largest eigenvalues. Generally, the noise concentrated the
spectrum at smaller eigenvectors as it grew. We observed that the edge of the
two distributions (signed and genuine) obeyed the same equation, and we found
two critical points as the variance of the noise increased: the first,
corresponding to the edge of the semi-circle of the random matrix integral
[2], occurred at the emergence of an outlier (the second largest eigenvalue)
and the second happened when the (trivial) largest eigenvalue merged with the
outlier. At the first transition point, the behavior of the location of the
edge of the distribution as a function of the size of the Gaussian noise
changed in a discrete manner that we determined. We have supported our
analytical expressions with numerical simulations.

Some interesting directions to pursue are, the corresponding problem for
higher rank tensors or in the line of [40, 50] the associated largest
eigenvalue of complex tensors, that is important for measures of entanglement
entropy of random quantum states, tensors with more general symmetries or with
correlated entries, in order to investigate the existence of universality
classes distinct from those of random matrices.

Another intriguing problem stems from viewing the random tensor as the
adjacency “matrix” of a hypergraph [51]. Understanding the leading eigenvalues
of the tensor and their correlations can teach us about the asymptotic
stability of the system on such random networks, analogous to the May-Wigner
instabilities (e.g. [52] for a different way to incorporate non-linearities)
or to turbulence [53]. Hopefully, one could return from the knowledge of
relevant eigenvectors to the study of simplicial complexes built from tensors.

More broadly, large $N$ corrections, fluctuations of the edge of the tensor
eigenspectrum, second and higher order moments of the distribution in the
spirit of [54, 55] are an important program that we leave for the future.

## Aknowledgements

We are grateful for the help and support provided by the Scientific Computing
and Data Analysis section of Core Facilities at OIST that allowed part of the
simulations. ND acknowledges the organizers of the YITP-ExU long-term workshop
Quantum Information, Quantum Matter and Quantum Gravity (QIMG2023), Yukawa
Institute for Theoretical Physics, Kyoto University, during which part of this
work was accomplished. We would also like to thank R. Toriumi for the
discussions at the initial stage of this work. NS is supported in part by JSPS
KAKENHI Grant No.19K03825 and ND is supported in part by JSPS KAKENHI Grant
No.21K20340.

##  Appendix A Numerical simulations

Let us explain the setting of our numerical simulations. Similarly to [56, 45,
41], we generate a random tensor $\sigma_{ijk}$ with components normally
distributed as a centered Gaussian of unit variance $\mathcal{N}(0,1)$, and
multiply with the appropriate combinatorial factor that accounts for the
symmetries of the indices:

$\displaystyle c_{ijk}=\frac{\sigma_{ijk}}{\sqrt{d(i,j,k)}}$ (69)
$\displaystyle d(i,j,k)=\begin{cases}1,\quad i=j=k,\\\ 3,\quad i=j\neq
k,\text{or }i\neq j=k,\text{or }i=k\neq j,\\\ 6,\quad i\neq j\neq k\neq
i,\end{cases}$ (70)

ensuring that the measure arises from the norm

$c^{2}=\sum_{1\leq a,b,c\leq N}c_{abc}^{2}=\sum_{1\leq i\leq j\leq k\leq
N}d(i,j,k)c_{ijk}^{2}$ (71)

Additionally, we generate a random Gaussian vector $\\{\delta v_{i}\\}_{1\leq
i\leq N}$ of variance $2\beta$ to account for the width of the Lagrange
multiplier $\lambda$. We solve the eigenvector equation

$v_{i}=c_{ijk}v_{j}v_{k}+\delta v_{i}$ (72)

and keep the real solutions, storing their absolute value $\absolutevalue{v}$
and the sign of the determinant of the Jacobian matrix $M_{ab}$ for the signed
case. The distribution of the absolute value $\absolutevalue{v}$ without and
with the weight of the sign of the determinant of $M_{ab}$ are represented
above in Fig. 1 and Fig 2.

##  Appendix B Finite $N$ distribution and interpretation of the solutions of
$h(y)=0$

We first elaborate on how a finite $N$ expression for the eigenvector
distribution can be obtained from [2]. We start from a rather trivial
expression,

$\rho^{\rm size}(v,\alpha,\beta)=\frac{\rho^{\rm
size}(v,\alpha,\beta)}{\rho^{\rm size}(v,\alpha,\beta=0)}\rho^{\rm
size}(v,\alpha,\beta=0),$ (73)

and take the advantage that the exact expression of $\rho^{\rm
size}(v,\alpha=1/2,\beta=0)$ for arbitrary $N$ is given in [2]. The relation
between their variable $u$ and our $\absolutevalue{v}$ is given by

$\frac{1}{|v|}=-\sqrt{N}u,$ (74)

which can be derived by comparing the stationary condition of the Hamiltonian
(energy) in [2] and our eigenvector equation (see an appendix of [41] for an
explicit derivation). Our $\rho^{\rm size}(|v|,\alpha=1/2,\beta=0)$ is equal
to the derivative of (2.9) of[2] with respect to $|v|$, considering that they
are the number density and the number, respectively. Then, by using
(7.4)-(7.8) of [2], we obtain the exact expression of $\rho^{\rm
size}(|v|,\alpha=1/2,\beta=0)$.121212More explicitly, $\rho^{\rm
size}(|v|,\alpha=1/2,\beta=0)=2\sqrt{N}|v|^{-2}(p-1)^{(N-1)/2}\exp(-(p-2)/(4(p-1)|v|^{2}))\rho_{N}(x)$,
where $\rho_{N}(\cdot)$ is that of [2], $x=-\sqrt{p/(2(p-1))}/(\sqrt{N}|v|)$,
and $p=3$.

Next let us estimate the first factor of (73). From (31) and using the
approximation (40), we find

$\frac{\rho^{\rm size}(v,\alpha,\beta)}{\rho^{\rm
size}(v,\alpha,\beta=0)}\sim\left(\frac{|v|^{4}}{|v|^{4}+4\alpha\beta}\right)^{\frac{1}{2}}\left(\frac{|v|^{4}}{|v|^{4}+12\alpha\beta}\right)^{\frac{N-1}{2}}\exp\left(-\frac{\alpha|v|^{2}}{|v|^{4}+4\alpha\beta}+\frac{\alpha}{|v|^{2}}\right)\frac{Z_{FT}^{\parallel}}{Z_{FT}^{\parallel}(\beta=0)},$
(75)

where we have in particular used the fact that $Z_{FT}^{\perp}$ is free from
$\beta$. Note that $Z_{FT}^{\parallel}$ can approximately be expressed by (50)
with $Q_{ij}$ as done in the main text.

Let us argue now that the roots of $h(y)=0$ correspond to single eigenvectors
which are isolated from the main part of the distributions. Having to deal
with single eigenvectors, we could not rely on the large $N$ expressions. As
shown in Fig. 3, the number of solutions of $h(y)=0$ changes depending on the
value of $\bar{\beta}$: two solutions for $\bar{\beta}<\bar{\beta}_{1}$, three
for $\bar{\beta}_{1}<\bar{\beta}<\bar{\beta}_{2}$, and one for
$\beta_{2}<\bar{\beta}$. Since $h(y)>0$ in the region $y<y_{\rm smallest}$,
where $y_{\rm smallest}$ denotes the smallest root of $h(y)=0$, it is clear
that $y_{\rm smallest}$ is the edge of the distribution in the large-$N$
limit. On the other hand, the other roots of $y$ are not edges of any regions
of $h(y)>0$, and the interpretation of these solutions was not so clear.

By generalizing the arguments of [43], the locations of the eigenvectors can
be estimated by solving

$\int_{0}^{|v_{n}|}\differential z\,\rho^{\rm
size}(z,\alpha,\beta)=n-\frac{1}{2},$ (76)

where $\absolutevalue{v_{n}}$ denotes the estimation of the absolute value of
the $n$-th smallest eigenvector, and we can use the above result for
$\rho^{\rm size}(z,\alpha,\beta)$ for $\alpha=1/2$.

In Fig. 5, the locations of the smallest eigenvectors up to $n=6$ are plotted
in $y$ against $\bar{\beta}$ for $N=50$ by solving (76). The value $N=50$
still seems too small for a definite conclusion, but, comparing with Fig. 3,
the behavior of the smallest two eigenvectors seems consistent with what the
roots of $h(y)=0$ imply: the first two smallest eigenvectors can be identified
with the isolated roots found for
$\bar{\beta}_{1}<\bar{\beta}<\bar{\beta}_{2}$, and they “disappear” after
$\bar{\beta}=\bar{\beta}_{2}$. (More precisely, in the region
$\bar{\beta}_{2}<\bar{\beta}$, the smallest eigenvector is succeeded by the
third one.)

Figure 5: The solutions to (76) up to $n=6$ are plotted in $y$ against
$\bar{\beta}$ for $N=50$. The lines of Fig. 3 are also drawn on the figure for
comparison.

##  Appendix C A comment on the large $\beta$ limit

As is shown in an example on Fig. 6, the effect of taking larger $\beta$
generally reduces the total number of eigenvectors. The question we want
discuss here is whether the number vanishes in the large $\beta$ limit. We
will show that there remain an exponentially large number of them.

Figure 6: The behavior of the mean number of eigenvectors against $\beta$ for
$N=10$. The number decreases, as the noise size $\beta$ becomes larger.

Since the difference between the genuine and signed distributions is to take
the absolute value of the determinant or not, the mean number of eigenvectors
is bounded from below by the signed distribution as

$\langle\\#{\rm Eigenvectors}\rangle=\int_{\mathbb{R}^{N}}\differential
v\,\rho(v)\geq\left|\int_{\mathbb{R}^{N}}\differential v\,\rho_{s}(v)\right|.$
(77)

In fact, $\int_{\mathbb{R}^{N}}\differential v\,\rho_{s}(v)$ is constant
against $\beta$, as will be shown below. Therefore, if it is non-zero,
$\langle\\#{\rm Eigenvectors}\rangle$ is bounded from below by a positive
number even for $\beta\rightarrow\infty$.

Let us prove the constancy of $\int_{\mathbb{R}^{N}}\differential
v\,\rho_{s}(v)$ against $\beta$. We first observe that since we integrate over
$v$, $v$ can be regarded as a dynamical variable of the system. We introduce
the following BRST symmetry $Q_{\rm BRST}$ for the variables
$v,\lambda,\bar{\psi},\psi$:

${Q_{BRST}\,}v_{a}=\bar{\psi}_{a},\ {Q_{BRST}\,}\lambda_{a}=0,\
{Q_{BRST}\,}\bar{\psi}_{a}=0,\ {Q_{BRST}\,}\psi_{a}=i\lambda_{a},$ (78)

which indeed satisfies $Q_{BRST}^{2}=0$. Then we observe

${Q_{BRST}\,}\left(\left(v_{a}-C_{abc}v_{b}v_{c}\right)\psi_{a}\right)=i\lambda_{a}(v_{a}-C_{abc}v_{b}v_{c})+\bar{\psi}_{a}(\delta_{ab}-2C_{abc}v_{c})\psi_{b},$
(79)

which agrees with the latter two terms of (10). This implies that those terms
form a BRST invariant action. In addition, we can further observe

$i{Q_{BRST}\,}\left(\lambda_{a}\psi_{a}\right)=-\lambda^{2},$

which means that the second term with $\beta$ in (10) is in fact BRST exact.
By a standard argument, the partition function

$\int\differential
v\differential\lambda\differential\bar{\psi}\differential\psi\,e^{-\beta\lambda^{2}+i\lambda_{a}(v_{a}-C_{abc}v_{b}v_{c})+\bar{\psi}_{a}(\delta_{ab}-2C_{abc}v_{c})\psi_{b}}$
(80)

is invariant under the change of $\beta$, which proves the invariance of the
integral of the signed distribution $\int_{\mathbb{R}^{N}}\differential
v\rho_{s}(v)$ under the change of $\beta$.

It is not difficult to compute the explicit value of the bound (77). The
integration of $\rho_{s}(v)$ over $v$ can most easily be computed for
$\beta=0$. Putting $\beta=0$ into (24), we obtain

$\displaystyle\begin{split}\int_{0}^{\infty}\differential|v|\,\rho_{s}^{\text{size}}(|v|,\beta=0)&=\frac{\sqrt{3}\alpha
2^{N/2}}{\Gamma(N/2)}\int_{0}^{\infty}\differential|v|\
|v|^{-3}e^{-\alpha/|v|^{2}}U\left(1-N/2,3/2,3\alpha/(2|v|^{2})\right)\cr&=\left\\{\begin{array}[]{ll}-1+\frac{(-1)^{N/2}2^{-N/2+1}\Gamma(N)}{\sqrt{3}\Gamma(N/2)\Gamma(N/2+1)}{}_{2}F_{1}(1/2,1;N/2+1;2/3)&\hbox{for
even }N,\\\ -1&\hbox{for odd }N,\end{array}\right.\end{split}$ (81)

where ${}_{2}F_{1}(\cdot,\cdot;\cdot;\cdot)$ is the hypergeometric function.
In fact, $\rho_{s}^{\text{size}}(|v|,\beta=0)$ is missing the trivial solution
$v=0$, as can be seen in its expression. Since the $v=0$ solution always
exists, and ${\rm det}M=1$ for the solution, one should add 1 to the result of
(81), which cancels the first terms of (81). 131313The main reason why we have
to count the trivial solution $v=0$ for $\beta=0$ is that the trivial solution
is moved to a $v\neq 0$ solution for $\beta>0$, and therefore is counted in
the integral $\int_{0}^{\infty}d|v|\,\rho_{s}(|v|)$. Hence for the integral to
be constant against $\beta$ for $\beta\geq 0$, the trivial solution must be
counted. Thus we obtain a bound,

$\langle\\#{\rm
Eigenvectors}\rangle\geq\frac{2^{-N/2+1}\Gamma(N)}{\sqrt{3}\Gamma(N/2)\Gamma(N/2+1)}{}_{2}F_{1}(1/2,1;N/2+1;2/3)\sim\frac{2^{(N+1)/2}}{\sqrt{3\pi
N}}$ (82)

for even $N$ and the asymptotic behavior for large $N$ has also been
estimated. For odd $N$, we cannot derive a meaningful bound from (77), but
(82) could be extended to the odd case too, since $\langle\\#{\rm
Eigenvectors}\rangle$ is numerically a smooth function of $N$, no matter if
$N$ is odd or even. We see that, even for $\beta\rightarrow\infty$, there
remain an exponentially large number of eigenvectors.

## References

  * [1] N. Sasakura, _Real eigenvector distributions of random tensors with backgrounds and random deviations_ , _Progress of Theoretical and Experimental Physics_ 2023 (2023) 123A01.
  * [2] A. Auffinger, G.B. Arous and J. Černỳ, _Random matrices and complexity of spin glasses_ , _Communications on Pure and Applied Mathematics_ 66 (2013) 165.
  * [3] G. Akemann, J. Baik and P. Di Francesco, _The Oxford handbook of random matrix theory_ , Oxford University Press (2011).
  * [4] J. Ambjørn, B. Durhuus and T. Jonsson, _Three-dimensional simplicial quantum gravity and generalized matrix models_ , _Modern Physics Letters A_ 6 (1991) 1133.
  * [5] M. Gross, _Tensor models and simplicial quantum gravity in $>2-D$_, _Nuclear Physics B-Proceedings Supplements_ 25 (1992) 144.
  * [6] N. Sasakura, _Tensor model for gravity and orientability of manifold_ , _Modern Physics Letters A_ 6 (1991) 2613.
  * [7] R. Gurau and V. Rivasseau, _Quantum gravity and random tensors_ , _Contribution to the Poincaré Seminar, December 16th 2023_ (2024)  [arXiv:2401.13510].
  * [8] D. Anninos and B. Mühlmann, _Notes on matrix models (matrix musings)_ , _Journal of Statistical Mechanics: Theory and Experiment_ 2020 (2020) 083109.
  * [9] R. Gurau, _The $1/N$ expansion of colored tensor models_, _Annales Henri Poincare_ 12 (2011) 829.
  * [10] R. Gurau and V. Rivasseau, _The $1/N$ expansion of colored tensor models in arbitrary dimension_, _EPL (Europhysics Letters)_ 95 (2011) .
  * [11] R. Gurau, _The complete $1/N$ expansion of colored tensor models in arbitrary dimension_, _Annales Henri Poincaré_ 13 (2011) .
  * [12] L. Lionni, _Colored Discrete Spaces: Higher Dimensional Combinatorial Maps and Quantum Gravity_ , Springer International Publishing (2018), 10.1007/978-3-319-96023-4.
  * [13] M.L. Mehta, _Random matrices_ , Elsevier (2004).
  * [14] L. Qi, H. Chen and Y. Chen, _Tensor eigenvalues and their applications_ , vol. 39, Springer (2018).
  * [15] D. Cartwright and B. Sturmfels, _The number of eigenvalues of a tensor_ , _Linear Algebra and its Applications_ 438 (2013) 942–952.
  * [16] R. Gurau, _On the generalization of the Wigner semicircle law to real symmetric tensors_ , 2004.02660.
  * [17] J.B. Geloun, _On the counting tensor model observables as $U(N)$ and $O(N)$ classical invariants_, _Contribution to the Corfu Summer Institute 2019_ (2020)  [2005.01773].
  * [18] I.R. Klebanov, F. Popov and G. Tarnopolsky, _TASI lectures on large $N$ tensor models_, _arXiv:1808.09434_ (2018) .
  * [19] C.J. Hillar and L.-H. Lim, _Most tensor problems are NP-hard_ , _Journal of the ACM (JACM)_ 60 (2013) 1.
  * [20] L. Qi, _Eigenvalues of a real supersymmetric tensor_ , _Journal of Symbolic Computation_ 40 (2005) 1302.
  * [21] E. Subag, _Concentration for the zero set of random polynomial systems_ , _arXiv:2303.11924_ .
  * [22] M. Kac, _On the average number of real roots of a random algebraic equation (ii)_ , _Proceedings of the London Mathematical Society_ 2 (1948) 390.
  * [23] S.O. Rice, _Mathematical analysis of random noise_ , _The Bell System Technical Journal_ 23 (1944) 282.
  * [24] P. Breiding, _How many eigenvalues of a random symmetric tensor are real?_ , _Transactions of the American Mathematical Society_ 372 (2019) 7857.
  * [25] A. Crisanti and H.-J. Sommers, _The spherical p-spin interaction spin glass model: the statics_ , _Zeitschrift für Physik B Condensed Matter_ 87 (1992) 341.
  * [26] A. Crisanti and H.-J. Sommers, _Thouless-Anderson-Palmer approach to the spherical p-spin spin glass model_ , _Journal de Physique I_ 5 (1995) 805.
  * [27] E. Subag and O. Zeitouni, _The extremal process of critical points of the pure p-spin spherical spin glass model_ , _Probability theory and related fields_ 168 (2017) 773.
  * [28] V. Ros and Y.V. Fyodorov, _The high-d landscapes paradigm: spin-glasses, and beyond_ , _Spin Glass Theory and Far Beyond: Replica Symmetry Breaking after 40 Years_ (2022) 24.
  * [29] O. Evnin, _Melonic dominance and the largest eigenvalue of a large random tensor_ , _Lett. Math. Phys._ 111 (2021) 66 [2003.11220].
  * [30] G. Parisi, _Nobel lecture: Multiple equilibria_ , _Reviews of Modern Physics_ 95 (2023) 030501.
  * [31] V. Piccolo, _Topological complexity of spiked random polynomials and finite-rank spherical integrals_ , _arXiv:2312.12323_ .
  * [32] Y. Fyodorov, _High-dimensional random fields and random matrix theory_ , _Markov Processes and Related Fields_ 21 (2015) 483.
  * [33] D. Belius, J. Černỳ, S. Nakajima and M.A. Schmidt, _Triviality of the geometry of mixed p-spin spherical hamiltonians with external field_ , _Journal of Statistical Physics_ 186 (2022) 12.
  * [34] J. Grela and B.A. Khoruzhenko, _Glass-like transition described by toppling of stability hierarchy_ , _Journal of Physics A: Mathematical and Theoretical_ 55 (2022) 154001.
  * [35] S. Péché, _Deformed ensembles of random matrices_ , in _Proceedings of the International Congress of Mathematicians, Seoul_ , vol. 3, pp. 1059–1174, 2014.
  * [36] J. Baik, G. Ben Arous and S. Péché, _Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices_ , _The Annals of Probability_ 33 (2005) 1643.
  * [37] E. Richard and A. Montanari, _A statistical model for tensor PCA_ , _Advances in neural information processing systems_ 27 (2014) .
  * [38] G.B. Arous, S. Mei, A. Montanari and M. Nica, _The landscape of the spiked tensor model_ , _Communications on Pure and Applied Mathematics_ 72 (2019) 2282.
  * [39] M. Ouerfelli, M. Tamaazousti and V. Rivasseau, _Selective multiple power iteration: from tensor PCA to gradient-based exploration of landscapes_ , _The European Physical Journal Special Topics_ (2023) 1.
  * [40] N. Sasakura, _Signed eigenvalue/vector distribution of complex order-three random tensor_ , _Progress of Theoretical and Experimental Physics_ (2024) .
  * [41] N. Sasakura, _Exact analytic expressions of real tensor eigenvalue distributions of gaussian tensor model for small $N$_, _Journal of Mathematical Physics_ 64 (2023) .
  * [42] J. Kurchan, _Replica trick to calculate means of absolute values: applications to stochastic equations_ , _Journal of Physics A: Mathematical and General_ 24 (1991) 4969.
  * [43] M.R. Kloos and N. Sasakura, _Usefulness of signed eigenvalue/vector distributions of random tensors_ , 2403.12427.
  * [44] G. Parisi, _Computing the number of metastable states in infinite-range models_ , _Les Houches_ 85 (2006) .
  * [45] N. Sasakura, _Real tensor eigenvalue/vector distributions of the Gaussian tensor model via a four-fermi theory_ , _Progress of Theoretical and Experimental Physics_ 2023 (2022) 013A02.
  * [46] F. Zamponi, _Mean field theory of spin glasses_ , 1008.4844.
  * [47] A. Cavagna, J.P. Garrahan and I. Giardina, _Quenched complexity of the mean-field p-spin spherical model with external magnetic field_ , _Journal of Physics A: Mathematical and General_ 32 (1999) 711.
  * [48] I. Giardina, _Course 10: Metastable states in glassy systems_ , _Les Houches_ 85 (2007) 373.
  * [49] E. Witten, _Analytic Continuation Of Chern-Simons Theory_ , _AMS/IP Stud. Adv. Math._ 50 (2011) 347 [1001.2933].
  * [50] S. Dartois and B. McKenna, _Injective norm of real and complex random tensors I: From spin glasses to geometric entanglement_ , 2404.03627.
  * [51] J. Li, X. Wu, J. Lü and L. Lei, _Enhancing predictive accuracy in social contagion dynamics via directed hypergraph structures_ , _Communications Physics_ 7 (2024) 129.
  * [52] Y.V. Fyodorov and B.A. Khoruzhenko, _Nonlinear analogue of the May-Wigner instability transition_ , _Proceedings of the National Academy of Sciences_ 113 (2016) 6827.
  * [53] S. Dartois, O. Evnin, L. Lionni, V. Rivasseau and G. Valette, _Melonic turbulence_ , _Communications in Mathematical Physics_ 374 (2020) 1179.
  * [54] E. Subag, _The complexity of spherical p-spin models - a second moment approach_ , _Annals of Probability_ 45 (2017) 3385.
  * [55] V. Ros, G. Ben Arous, G. Biroli and C. Cammarota, _Complex energy landscapes in spiked-tensor and simple glassy models: Ruggedness, arrangements of local minima, and phase transitions_ , _Physical Review X_ 9 (2019) .
  * [56] N. Sasakura, _Signed distributions of real tensor eigenvectors of Gaussian tensor model via a four-fermi theory_ , _Phys. Lett. B_ 836 (2023) 137618 [2208.08837].