YITP-22-98
Real tensor eigenvalue/vector distributions
of the Gaussian tensor model via a four-fermi theory

Eigenvalue distributions are important dynamical quantities in matrix models. The distributions give quantitative/qualitative insights into complex dynamical systems with strongly coupled degrees of freedom, such as atomic systems [1]. Topological transitions of the distributions [2, 3] give important insights into the properties of gauge theories and two-dimensional quantum gravity, such as phase structures in particular. The distributions are also important as one of the main tools to solve the matrix models [4].

Considering recent attentions to tensor models [5, 6, 7, 8] it would be interesting to study corresponding quantities and their roles in tensor models. While, in various applications [9]²²2See for example [10] for the tensor model perspectives on such applications., it is important to develop techniques to obtain eigenvalues/vectors [11, 12, 13] for given tensors, it is rather more important to understand statistical properties of tensor eigenvalues/vectors for ensembles of tensors in tensor models, where tensors are dynamical. In fact, not many results are known about their statistical properties: The expectation values of the numbers of real eigenvalues are computed for random real tensors [14, 15]; the sizes of the largest eigenvalues are estimated for random tensors [16]; Wigner semi-circle law in matrix models has been extended to tensor models [17].

A basic question about eigenvalues/vectors in tensor models is their distributions. In a previous study [18], the present author derived an exact formula³³3The formula is expressed with the confluent hypergeometric function of the second kind. for signed distributions of real eigenvectors for real symmetric order-three random tensors with the Gaussian distribution: Each eigenvector contributes to the distribution by $\pm 1$ depending on the sign of an associated Hessian matrix. An interesting intermediate step was that the problem was rewritten as the computation of a partition function of a four-fermi theory. This was achieved by Gaussian integrations over bosonic variables after rewriting the problem as a partition function of a system of bosonic and fermionic variables. The final formula was obtained by exactly computing the partition function of the four-fermi theory.

The above procedure of the previous paper can be extended to other kinds of eigenvalue/vector distributions. The purpose of the present paper is to extend it to the distribution of real eigenvalues/vectors for real symmetric order-three random tensors with the Gaussian distribution. In the previous paper, we only dealt with one couple of fermions, but, in this paper, we deal with $R$ replicas of fermions, and the eigenvector/value distribution is supposed to be obtained by an analytic continuation to $R=1/2$. We derive a four-fermi theory with the replicated fermions, which is more complex than that of the previous paper. We exactly compute the partition function of the four-fermi theory for some small-$N,R$ cases and show the precise agreement with Monte Carlo results. On the other hand, for large-$N$, we apply an approximation using the Schwinger-Dyson equation (See Appendix A), and obtain the eigenvalue/vector distribution by formally taking $R=1/2$. We find that the obtained approximate analytic expression has good agreement with Monte Carlo results, if we multiply an extra overall factor depending on $N$ to the expression, where $N$ denotes the range of the tensor indices.

This paper is organized as follows. In Section 2, we define the distribution of the real tensor eigenvalues/vectors of the Gaussian tensor model, and rewrite the distribution as a partition function of a system of bosonic and fermionic variables. In Section 3, by Gaussian integrations over the bosonic variables of the system, we obtain a four-fermi theory. In Section 4, we exactly compute the partition function of the four-fermi theory for some small-$N,R$ cases. In Section 5, we compute the partition function of the four-fermi theory for large-$N$ by applying an approximation using the Schwinger-Dyson equation for two-point functions, and obtain an approximate analytic expression of the eigenvalue/vector distribution by formally taking $R=1/2$. In Section 6, we perform some Monte Carlo simulations. The exact computations for small-$N,R$ are shown to precisely agree with the Monte Carlo results. On the other hand, we find that the approximate analytic expression of the eigenvalue/vector distribution agrees well with the Monte Carlo result, if it is multiplied by an extra factor depending on $N$. The final section is devoted to a summary and future problems. In Appendix C, we explain an overlap between the main result of Section 5 and a known result from random matrix theory.

In this paper, we restrict ourselves to the symmetric real tensors with three indices, $C_{abc}=C_{bac}=C_{bca}\in\mathbb{R}\ (a,b,c=1,2,\ldots,N)$, as the simplest case. There are several different definitions of tensor eigenvalues/vectors [11, 12, 13]. In this paper we employ the definition that real eigenvectors of a given $C$ are defined by $v$’s satisfying

$\displaystyle C_{abc}v_{b}v_{c}=v_{a},\ v\in\mathbb{R}^{N}.$

(1)

Here repeated indices are assumed to be summed over, as will also be assumed throughout this paper. Then the real eigenvector distribution for a given $C$ is given by

$\displaystyle\rho(v,C)=\sum_{i=1}^{N_{\rm sol}(C)}\prod_{a=1}^{N}\delta(v_{a}-v_{a}^{i}),$

(2)

where $v^{i}\ (i=1,2,\ldots,N_{\rm sol}(C))$ are all the solutions of (1).

The expression (2) can be rewritten as

$\displaystyle\rho(v,C)=\left|\det M\right|\,\prod_{a=1}^{N}\delta\left(v_{a}-C_{abc}v_{b}v_{c}\right),$

(3)

where $\det$ represents taking the determinant, $|\cdot|$ is to take the absolute value, and $M$ is a matrix defined by

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

(4)

Suppose the tensor $C$ is randomly distributed with the Gaussian distribution. Then the distribution of $v$ is given by

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

(5)

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

In the previous paper [18], we computed a similar quantity which has $\det M$ instead of $|\det M|$ in (3). The previous quantity was easier to compute, since it can simply be expressed by a couple of fermions by using the formula, $\det M=\int d{\bar{\psi}}d\psi\,e^{{\bar{\psi}}_{a}M_{ab}\psi_{b}}$ [19]. On the other hand, what makes the expression (3) difficult to deal with is that $|\det M|$ is not analytic. One way to turn it to an analytic expression is to consider a quantity,

$\displaystyle\rho(v,C,R,\epsilon)=\left(\det\left(M^{2}+\epsilon I\right)\right)^{R}\prod_{a=1}^{N}\delta\left(v_{a}-C_{abc}v_{b}v_{c}\right),$

(6)

where $\epsilon>0$ is a regularization parameter, and $I$ is the identity matrix, $I_{ab}=\delta_{ab}$. Note that the regularization parameter keeps the $O(N)$ invariance of the system⁴⁴4The system is invariant under the orthogonal transformation in the vector space associated to the lower index.. Then, for random $C$ with the Gaussian distribution, we have

$\displaystyle\rho(v,R,\epsilon)=\langle\rho(v,C,R,\epsilon)\rangle_{C}=A^{-1}\int_{\mathbb{R}^{D_{C}}}dC\,e^{-\alpha C^{2}}\left(\det\left(M^{2}+\epsilon I\right)\right)^{R}\prod_{a=1}^{N}\delta\left(v_{a}-C_{abc}v_{b}v_{c}\right).$

(7)

The expression corresponding to (5) can be obtained by putting $R=1/2$ and taking the $\epsilon\rightarrow+0$ limit:

$\displaystyle\rho(v)=\rho(v,1/2,+0).$

(8)

As we will find in Section 5, the regularization parameter $\epsilon$ is necessary for the large-$N$ approximate computation, to unambiguously obtain an expression which is valid in all the region of $v$: Directly starting from the expression with $\epsilon=0$ has a singularity and it is not clear how the expression should be extended to the whole region.

The real eigenvalues accompanied with real eigenvectors (Z-eigenvalues in the terminology of [11]) for a real symmetric order-three tensor $C$ are the $z$’s satisfying

$\displaystyle C_{abc}w_{b}w_{c}=z\,w_{a},$

(9)

where $z\in\mathbb{R}$, $w\in\mathbb{R}^{N}$with $|w|=1$ ($|w|=\sqrt{w_{a}w_{a}}$). The relation to (1) is $w=v/|v|$ and $z=1/|v|$. Therefore, once we obtain $\rho(v)$, the Z-eigenvalue distribution is obtained by

$\displaystyle\begin{split}\rho_{\rm eigenvalue}(z)&=\rho(v)S^{N-1}|v|^{N-1}\frac{d|v|}{dz}\\ &=\rho(1/z)S_{N-1}z^{-N-1},\end{split}$

(10)

where $S_{N-1}=2\pi^{N/2}/\Gamma(N/2)$ is the surface volume of a unit sphere in an $N$-dimensional space. To derive this expression, we have used that $\rho(v)$ actually depends only on $|v|$, because of the $O(N)$-invariance of (5), and have used $d^{N}v=S_{N-1}|v|^{N-1}d|v|$. Note that we have abusively written $1/z$ as the argument on the righthand side of (10) to represent an arbitrary vector of size $1/z$.

By using the well-known formulas, $\int_{\mathbb{R}}dx\,e^{ixy}=2\pi\delta(y)$ and $\int d{\bar{\psi}}d\psi\,e^{{\bar{\psi}}_{a}T_{ab}\psi_{b}}=\det T$ [19], the expression (7) can be rewritten as

$\displaystyle\rho(v,R,\epsilon)=A^{-1}(2\pi)^{-N}\int dCd\lambda d\bar{\psi}d\psi\,e^{S_{1}},$

(11)

where

$\displaystyle S_{1}=-\alpha C^{2}+i\lambda_{a}(v_{a}-C_{abc}v_{b}v_{c})+{\bar{\psi}}_{a}^{i}M^{2}_{ab}\psi_{bi}+\epsilon{\bar{\psi}}_{a}^{i}\psi_{ai}.$

(12)

Here we have assumed $R$ to be integer, $\bar{\psi}_{a}^{i},\psi_{ai}\ (a=1,2,\ldots,N;\ i=1,2,\ldots,R)$ are fermionic, and $\lambda_{a}\ (a=1,2,\ldots,N)$ are bosonic with the integration region $\mathbb{R}^{N}$. Note that the paired new indices (namely, $i$) are also assumed to be summed over in the above expression, as will be assumed throughout this paper. Note also that the system (11) is invariant under the following $GL(R)$ transformation concerning the new index:

$\displaystyle\begin{split}\psi^{\prime}{}_{ai}&=G{}_{i}{}^{i^{\prime}}\psi_{ai^{\prime}},\ \bar{\psi}^{\prime}{}_{a}^{i}=\bar{\psi}^{i^{\prime}}_{a}G^{-1}{}_{i^{\prime}}{}^{i},\end{split}$

(13)

for $G\in GL(R)$. Another comment is that the appearance of the imaginary number in (12) does not violate the reality of the integral, since the integral is symmetric under $\lambda\rightarrow-\lambda$.

In the similar manner as the procedure taken in the previous paper [18], we want to first integrate over the bosonic variables, $C$ and $\lambda$, to obtain a fermionic theory. However, $M^{2}$ contains $C$ in a quadratic manner and is not straightforward to deal with. To circumvent this matter, let us introduce another pair of fermions $\bar{\varphi},\varphi$ to rewrite (12) into an expression linear in $M$:

$\displaystyle\rho(v,R,\epsilon)=A^{-1}(2\pi)^{-N}(-1)^{NR}\int dCd\lambda d\bar{\psi}d\psi d\bar{\varphi}d\varphi\,e^{S_{2}},$

(14)

where

$\displaystyle\begin{split}S_{2}=-\alpha C^{2}+&i\lambda_{a}(v_{a}-C_{abc}v_{b}v_{c})-\bar{\varphi}_{a}^{i}\varphi_{ai}+\epsilon\bar{\psi}_{a}^{i}\psi_{ai}-\bar{\psi}_{a}^{i}M_{ab}\varphi_{bi}-\bar{\varphi}_{a}^{i}M_{ab}\psi_{bi}.\end{split}$

(15)

The equivalence between (11) and (14) can be shown by

$\displaystyle\begin{split}\int d\bar{\varphi}d\varphi\,e^{-\bar{\varphi}_{a}\varphi_{a}-\bar{\psi}_{a}M_{ab}\varphi_{b}-\bar{\varphi}_{a}M_{ab}\psi_{b}}&=\int d\bar{\varphi}d\varphi\,e^{-(\bar{\varphi}_{a}+\bar{\psi}_{b}M_{ba})(\varphi_{a}+M_{ac}\psi_{c})+\bar{\psi}_{a}M^{2}_{ab}\psi_{b}}\\ &=(-1)^{N}e^{\bar{\psi}_{a}M^{2}_{ab}\psi_{b}}.\end{split}$

(16)

Now let us first integrate over $C$, assuming that the final result does not depend on this change of the order of the integrations. The terms containing $C$ in (15) are

$\displaystyle S_{C}=-\alpha C^{2}-iC_{abc}\lambda_{a}v_{b}v_{c}+2C_{abc}\bar{\psi}_{a}^{i}\varphi_{bi}v_{c}+2C_{abc}\bar{\varphi}_{a}^{i}\psi_{bi}v_{c}.$

(17)

Therefore the integration over $C$ is just a Gaussian integration and we obtain

$\displaystyle\int_{\mathbb{R}^{D_{C}}}dCe^{S_{C}}=A\,e^{\delta S_{C}},$

(18)

where

$\displaystyle\begin{split}\delta S_{C}&=\frac{1}{\alpha}\left(\frac{1}{6}\sum_{\sigma}\left(-\frac{i}{2}\lambda_{\sigma_{a}}v_{\sigma_{b}}v_{\sigma_{c}}+\bar{\psi}^{i}_{\sigma_{a}}\varphi_{\sigma_{b}i}v_{\sigma_{c}}+\bar{\varphi}^{i}_{\sigma_{a}}\psi_{\sigma_{b}i}v_{\sigma_{c}}\right)\right)^{2}.\end{split}$

(19)

Here the summation over $\sigma$ is over all the permutations of $a,b,c$, the necessity of which comes from that $C$ is a symmetric tensor. Explicitly expanding the expression (19), we obtain

$\displaystyle\delta S_{C}=-\frac{|v|^{4}}{12\alpha}\lambda_{a}B_{ab}\lambda_{b}-i\lambda_{a}D_{a}+E$

(20)

with

	$\displaystyle B_{ab}=\delta_{ab}+2\hat{v}_{a}\hat{v}_{b}=I_{\perp\,ab}+3I_{\parallel\,ab},$		(21)
	$\displaystyle D_{a}=\frac{|v|^{3}}{3\alpha}\left({\bar{\psi}}_{\parallel}^{i}{\varphi}_{\parallel}{}_{i}+\bar{\varphi}_{\parallel}^{i}\psi_{\parallel}{}_{i}\right)\hat{v}_{a}+\frac{|v|^{3}}{3\alpha}\left({\bar{\psi}}_{a}^{i}{\varphi}_{\parallel}{}_{i}+{\bar{\psi}}_{\parallel}^{i}{\varphi}_{ai}+\bar{\varphi}_{a}^{i}\psi_{\parallel}{}_{i}+\bar{\varphi}_{\parallel}^{i}\psi_{ai}\right),$		(22)
	$\displaystyle E=\frac{1}{\alpha}\left(\frac{1}{6}\sum_{\sigma}\left({\bar{\psi}}_{\sigma_{a}}^{i}{\varphi}_{\sigma_{b}i}v_{\sigma_{c}}+\bar{\varphi}_{\sigma_{a}}^{i}\psi_{\sigma_{b}i}v_{\sigma_{c}}\right)\right)^{2},$		(23)

where $\hat{v}_{a}=v_{a}/|v|$, and $\parallel$ denotes the projection to $\hat{v}$, namely, $\psi_{\parallel}=\psi_{a}\hat{v}_{a}$, and so on. As above, the matrix $B$ can be rewritten as a linear sum of the projection operators, where $I_{\parallel}$ is the projection operator to the one-dimensional linear subspace spanned by $\hat{v}$, and $I_{\perp}$ is to the subspace transverse to $\hat{v}$.

From (15) and (20), the remaining terms containing $\lambda$ are

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

(24)

Integrating over $\lambda$ gives

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

(25)

with

$\displaystyle\begin{split}\delta S_{\lambda}&=-3\alpha|v|^{-4}(v_{a}-D_{a})B^{-1}_{ab}(v_{b}-D_{b})\\ &=-3\alpha|v|^{-4}\left(D_{\perp}\cdot D_{\perp}+\frac{1}{3}(|v|-D_{\parallel})^{2}\right)\\ &=-\frac{\alpha}{|v|^{2}}+2\left({\bar{\psi}}_{\parallel}^{i}{\varphi}_{\parallel}^{i}+\bar{\varphi}_{\parallel}^{i}\psi_{\parallel}^{i}\right)-3\alpha|v|^{-4}\left(D_{\perp}\cdot D_{\perp}+\frac{1}{3}D_{\parallel}^{2}\right),\end{split}$

(26)

where we have used $B^{-1}=I_{\perp}+\frac{1}{3}I_{\parallel}$, and from (22) the projections of $D$ to the parallel and transverse directions to $\hat{v}$ are respectively given by

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

(27)

While the last term of (26) gives some four-fermi interaction terms, the second term gives some corrections to the quadratic terms in the parallel direction.

Collecting all the results above and using $\det B=3$, we obtain a four-fermi theory,

$\displaystyle\rho(v,R,\epsilon)=3^{\frac{N-1}{2}}\pi^{-\frac{N}{2}}\alpha^{\frac{N}{2}}|v|^{-2N}e^{-\frac{\alpha}{|v|^{2}}}(-1)^{NR}\int d{\bar{\psi}}d\psi d\bar{\varphi}d{\varphi}\,e^{K+K_{\parallel}+V},$

(28)

where the quadratic term $K,\,K_{\parallel}$ and the four-fermi interaction terms $V$ are respectively given by

$\displaystyle\begin{split}K&=\left(\begin{array}[]{cc}\bar{\psi}_{\perp}^{i}&\bar{\varphi}_{\perp}^{i}\end{array}\right)\left(\begin{array}[]{cc}\epsilon&-1\cr-1&-1\end{array}\right)\left(\begin{array}[]{c}\psi_{\perp}{}_{i}\cr\varphi_{\perp}{}_{i}\end{array}\right),\\ K_{\parallel}&=\left(\begin{array}[]{cc}{\bar{\psi}}_{\parallel}^{i}&\bar{\varphi}_{\parallel}^{i}\end{array}\right)\left(\begin{array}[]{cc}\epsilon&1\cr 1&-1\end{array}\right)\left(\begin{array}[]{c}\psi_{\parallel}{}_{i}\cr{\varphi}_{\parallel}{}_{i}\end{array}\right),\\ V&=-3\alpha|v|^{-4}\left(D_{\perp}\cdot D_{\perp}+\frac{1}{3}D_{\parallel}^{2}\right)+E\\ &=-\frac{|v|^{2}}{6\alpha}\Big{(}\bar{\psi}_{\perp}^{i}\cdot\bar{\psi}_{\perp}^{j}\varphi_{\perp}{}_{i}\cdot\varphi_{\perp}{}_{j}+\bar{\psi}_{\perp}^{i}\cdot\varphi_{\perp}{}_{j}\bar{\psi}_{\perp}^{j}\cdot\varphi_{\perp}{}_{i}+\bar{\varphi}_{\perp}^{i}\cdot\bar{\varphi}_{\perp}^{j}\psi_{\perp}{}_{i}\cdot\psi_{\perp}{}_{j}\\ &\ \ \ \ \ \ \ \ \ \ \ \ +\bar{\varphi}_{\perp}^{i}\cdot\psi_{\perp}{}_{j}\bar{\varphi}_{\perp}^{j}\cdot\psi_{\perp}{}_{i}+2\bar{\psi}_{\perp}^{i}\cdot\bar{\varphi}_{\perp}^{j}\varphi_{\perp}{}_{i}\cdot\psi_{\perp}{}_{j}+2\bar{\psi}_{\perp}^{i}\cdot\psi_{\perp}{}_{j}\bar{\varphi}_{\perp}^{j}\cdot\varphi_{\perp}{}_{i}\Big{)}.\end{split}$

(29)

What is surprising in $V$ in (29) is that the parallel components are not contained. Therefore, the parallel components are free theory and can trivially be integrated over:

$\displaystyle\int d{\bar{\psi}}_{\parallel}d\psi_{\parallel}d\bar{\varphi}_{\parallel}d{\varphi}_{\parallel}\,e^{K_{\parallel}}=(-1)^{R},$

(30)

where we have already taken the $\epsilon\rightarrow+0$ limit, since this is smooth. Then we finally obtain a four-fermi theory containing only the transverse components,

$\displaystyle\rho(v,R,\epsilon)=3^{\frac{N-1}{2}}\pi^{-\frac{N}{2}}\alpha^{\frac{N}{2}}|v|^{-2N}e^{-\frac{\alpha}{|v|^{2}}}(-1)^{(N-1)R}\int d\bar{\psi}_{\perp}d\psi_{\perp}d\bar{\varphi}_{\perp}d\varphi_{\perp}\,e^{K+V}.$

(31)

Simple crosschecks of the formula, (31) with (29), can be performed by exactly computing it. In the $N=1$ case, we have no transverse directions, and hence (31) gives

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

(32)

Indeed, from that the eigenvector is given by $v=1/C$ for $N=1$, the Gaussian distribution of $C$ is the same as (32), because $\pi^{-\frac{1}{2}}\alpha^{\frac{1}{2}}e^{-\alpha C^{2}}dC=\rho_{N=1}(v,R,0)dv$.

Another obvious consequence is that, for integer $R$, the four-fermi theory (31) can in principle be computed by expanding the exponential in $K,V$:

$\displaystyle\int d\bar{\psi}_{\perp}d\psi_{\perp}d\bar{\varphi}_{\perp}d\varphi_{\perp}\,e^{K+V}=\sum_{n=0}^{(N-1)R}\frac{1}{n!(2(N-1)R-2n)!}\int d\bar{\psi}_{\perp}d\psi_{\perp}d\bar{\varphi}_{\perp}d\varphi_{\perp}\,V^{n}\,K^{2(N-1)R-2n}.$

(33)

This formula comes from the fact that each term of the integrand must be a product of exactly $4(N-1)R$ fermions for the fermionic integral to be non-zero [19]. The $\epsilon\rightarrow+0$ limit of (33) is straightforward, and, for integer $R$, (31) generally has the form,

$\displaystyle\rho(v,R,0)=3^{\frac{N-1}{2}}\pi^{-\frac{N}{2}}\alpha^{\frac{N}{2}}|v|^{-2N}e^{-\frac{\alpha}{|v|^{2}}}{\cal L}_{N,R},$

(34)

where ${\cal L}_{N,R}$ is a polynomial function of $x=|v|^{2}/(6\alpha)$, ${\cal L}_{N,R}=1+c_{1}x+c_{2}x^{2}+\cdots+c_{(N-1)R}x^{(N-1)R}$, where $c_{i}$ are some coefficients depending on $N,R$. Then, as in (8), $\rho(v)$ could be obtained by taking the analytic continuation of ${\cal L}_{N,R}$ to $R=1/2$.

However, the exact computation of ${\cal L}_{N,R}$ for general $N,R$ seems to be a difficult task. It becomes more and more cumbersome very quickly, as $N,R$ increase. Still it is doable for some small $N,R$ by using computers. By a Mathematica package for Grassmann variables [20], we obtain

$\displaystyle\begin{split}{\cal L}_{N=2,R=1}&=1+4x,\\ {\cal L}_{N=3,R=1}&=1+4x+28x^{2},\\ {\cal L}_{N=2,R=2}&=1+24x+48x^{2},\\ {\cal L}_{N=3,R=2}&=1+40x+552x^{2}+1248x^{3}+5136x^{4}.\end{split}$

(35)

The formula (34) with (35) will be checked by Monte Carlo simulations in Section 6.

As was demonstrated in the last part of Section 4, the expression (31) can in principle be computed for general integer $N,R$, and $R=1/2$ could be taken by analytic continuation. However, this does not seem straightforward in practice. Therefore, in this section, we apply the Schwinger-Dyson equation to the four-fermi theory, assuming the factorization of the four-point correlation functions into the products of two-point functions. A self-contained brief explanation of the approximation is given in Appendix A.

The four-fermi theory (31) has the following symmetry, $O(N-1)\times GL(R)$, for the transverse variables $X,\,\bar{X}\ (X=\psi_{\perp},\varphi_{\perp})$:

$\displaystyle\begin{split}X^{\prime}{}_{ai}&=g_{1}{}_{a}{}^{a^{\prime}}g_{2}{}_{i}{}^{i^{\prime}}X_{a^{\prime}i^{\prime}},\ \bar{X}^{\prime}{}_{a}^{i}=g_{1}{}_{a}{}^{a^{\prime}}\bar{X}^{i^{\prime}}_{a^{\prime}}g_{2}^{-1}{}_{i^{\prime}}{}^{i},\end{split}$

(36)

where $g_{1}$ belongs to the defining representation of $O(N-1)$ in the subspace transverse to $\hat{v}$, and $g_{2}$ belongs to the defining representation of $GL(R)$. Assuming that these symmetries do not break down in the large-$N$ limit, we can assume the following forms of the two point functions:

$\displaystyle\begin{split}\langle\bar{\psi}_{\perp}{}_{a}^{i}\psi_{\perp}{}_{bj}\rangle&=Q_{11}I_{\perp\,ab}\,\delta^{i}_{j},\\ \langle\bar{\psi}_{\perp}{}_{a}^{i}\varphi_{\perp}{}_{bj}\rangle&=Q_{12}I_{\perp\,ab}\,\delta^{i}_{j},\\ \langle\bar{\varphi}_{\perp}{}_{a}^{i}\psi_{\perp}{}_{bj}\rangle&=Q_{21}I_{\perp\,ab}\,\delta^{i}_{j},\\ \langle\bar{\varphi}_{\perp}{}_{a}^{i}\varphi_{\perp}{}_{bj}\rangle&=Q_{22}I_{\perp\,ab}\,\delta^{i}_{j},\\ \hbox{Others}&=0,\end{split}$

(37)

where $Q_{ij}$ are some numbers. From (69) in Appendix A, the effective action written in terms of $Q_{ij}$ is given by

$\displaystyle\begin{split}S_{\rm eff}&=\langle K+V\rangle-(N-1)R\log(-\det Q)\\ &=(N-1)R\left(\epsilon Q_{11}-Q_{12}-Q_{21}-Q_{22}-\frac{|v|^{2}(N-1)}{6\alpha}\left(Q_{12}^{2}+Q_{21}^{2}+2Q_{11}Q_{22}\right)-\log(-\det Q)\right)\end{split}$

(38)

in the leading order of $N-1$. Here the minus sign in the logarithm is appropriate for the current computation, since the solution of $Q_{ij}$ below gives $\det Q<0$. In the derivation of (38), we have used the factorizations into two-point functions and (37), such as

$\displaystyle\begin{split}\langle\bar{\psi}_{\perp}^{i}\cdot\varphi_{\perp}{}_{j}\bar{\psi}_{\perp}^{j}\cdot\varphi_{\perp}{}_{i}\rangle&=\langle\bar{\psi}_{\perp}^{i}{}_{a}\varphi_{\perp}{}_{aj}\rangle{}\langle\bar{\psi}_{\perp}^{j}{}_{b}\varphi_{\perp}{}_{bi}\rangle-\langle\bar{\psi}_{\perp}^{i}{}_{a}\varphi_{\perp}{}_{bi}\rangle{}\langle\bar{\psi}_{\perp}^{j}{}_{b}\varphi_{\perp}{}_{aj}\rangle\\ &=(N-1)^{2}RQ_{12}^{2}-(N-1)R^{2}Q_{12}^{2},\end{split}$

(39)

and have taken only the leading order terms in $N-1$.

The stationary equations,

$\displaystyle\frac{\partial S_{\rm eff}}{\partial Q_{11}}=\frac{\partial S_{\rm eff}}{\partial Q_{12}}=\frac{\partial S_{\rm eff}}{\partial Q_{21}}=\frac{\partial S_{\rm eff}}{\partial Q_{22}}=0,$

(40)

have four independent solutions, and, among them, we should choose a solution which reproduces the $v=0$ limit (this is a free theory limit, since $V=0$) given by

$\displaystyle Q_{11}=\frac{1}{1+\epsilon},\ Q_{12}=Q_{21}=\frac{-1}{1+\epsilon},\ Q_{22}=\frac{\epsilon}{1+\epsilon},\hbox{ for }v=0.$

(41)

The explicit expression of the solution is given in Appendix B with $x=v^{2}(N-1)/(3\alpha)$. The $\epsilon\rightarrow+0$ limit has a singularity at $x=1/4$, and for each region, we obtain for $\epsilon\rightarrow+0$

• •

  $0<x<1/4$

  $\displaystyle\begin{split}&Q_{11}=\frac{-\sqrt{1-4x}+1}{2x\sqrt{1-4x}},\\ &Q_{12}=Q_{21}=\frac{1-\sqrt{1-4x}-4x}{2x\sqrt{1-4x}},\\ &Q_{22}=0,\end{split}$

  (42)

• •

  $x>1/4$

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

  (43)

The latter solution is not well-defined in $\epsilon\rightarrow+0$, but the $S_{\rm eff}$ for those solutions has well-defined limits in both regions:

$\displaystyle S^{\epsilon=+0}_{\rm eff}(x)=\left\{\begin{array}[]{ll}\frac{(N-1)R}{2}\left(2+\log(16)+\frac{1-\sqrt{1-4x}}{x}-4\log\left(1-\sqrt{1-4x}\right)+4\log(x)\right)&\hbox{ for }0<x<\frac{1}{4},\\ \frac{(N-1)R}{2x}(1+2x+2x\log(x))&\hbox{ for }\frac{1}{4}<x.\end{array}\right.$

(46)

Then from (71) we obtain for large-$N$

$\displaystyle\rho(v,R,+0)=3^{\frac{N-1}{2}}\pi^{-\frac{N}{2}}\alpha^{\frac{N}{2}}|v|^{-2N}e^{-\frac{\alpha}{|v|^{2}}}e^{S^{\epsilon=+0}_{\rm eff}(x)-2(N-1)R},$

(47)

where the subtraction in the exponent is due to $S^{\epsilon=+0}_{\rm eff}(0)=2(N-1)R$. Note that, in the above computation, the limiting process of $\epsilon\rightarrow+0$ is independent from taking $R=1/2$, and therefore the result does not depend on the order of them.

A delicate matter concerning taking $R=1/2$ is that $R$ in the above computation is assumed to be positive integers, since $R$ is the number of replicas. Therefore we do not know whether this expression can truly be analytically extended to $R=1/2$. This is generally a difficult question to answer, which is often raised in applying the replica trick. However, as is shown in Appendix C, the above result (47) with (46) can be shown to agree with a large-$N$ expression computed by random matrix theory in another context.

Finally, we want to argue that only the expression at $x>1/4$ in (46) practically matters, when $N$ is large. The reason is as follows. Because of the exponential damping factor $e^{-\frac{\alpha}{|v|^{2}}}$, $\rho(v,R,+0)$ takes significant values only in the region $|v|\gtrsim\sqrt{\alpha}$. This means that the relevant region satisfies $x=|v|^{2}(N-1)/(3\alpha)\gtrsim(N-1)/3>1/4$ for large $N$. Therefore, we can practically use $S_{\rm eff}^{\epsilon=+0}$ at $x>1/4$ for all the region of $x$ for large-$N$. By doing this, we finally obtain, for large-$N$ and $R=1/2$,

$\displaystyle\rho(v)=\rho(v,1/2,+0)=(N-1)^{\frac{N-1}{2}}e^{-\frac{N-1}{2}}\pi^{-\frac{N}{2}}\alpha^{\frac{1}{2}}|v|^{-N-1}e^{-\frac{\alpha}{4|v|^{2}}}.$

(48)

The result (48) is for the distribution of vector $v$. For comparison with Monte Carlo simulations, it is more convenient to transform it to the distribution of the size $|v|$. Since $d^{N}v=S_{N-1}|v|^{N-1}d|v|$ with $S_{N-1}=2\pi^{N/2}/\Gamma(N/2)$ being the surface volume of a unit sphere in $N$-dimensions, we obtain for large-$N$

$\displaystyle\rho_{\rm size}(|v|)=\frac{2(N-1)^{\frac{N-1}{2}}e^{-\frac{N-1}{2}}\alpha^{\frac{1}{2}}}{\Gamma\left(\frac{N}{2}\right)}|v|^{-2}e^{-\frac{\alpha}{4|v|^{2}}}.$

(49)

One can also transform the result to the real eigenvalue distribution (Z-eigenvalues in the terminology of [11]), using (10): For large-$N$,

$\displaystyle\rho_{\rm eigevalue}(z)=\frac{2(N-1)^{\frac{N-1}{2}}e^{-\frac{N-1}{2}}\alpha^{\frac{1}{2}}}{\Gamma\left(\frac{N}{2}\right)}e^{-\frac{\alpha}{4}z^{2}}.$

(50)

Interestingly, the real eigenvalue distribution is just given by the Gaussian function in this approximation.

In this section, we perform some Monte Carlo simulations, and compare the results with the exact small-$N$ and the approximate large-$N$ results, which are respectively obtained in Sections 4 and 5.

The eigenvector equations (1) are a system of polynomial equations, which can be solved by appropriate polynomial equation solvers. We use Mathematica 12 for this purpose. It generally produces all the solutions including complex ones, and we pick up all the real ones out of them. Whether all the real eigenvectors are covered can be guaranteed by checking whether the number of all the generally complex solutions agrees with the known formula, $2^{N}-1$ [13].

With the above way of obtaining real eigenvectors, our procedure of numerical simulations can be summarized as follows.

• •

  Randomly generate $C$ with the Gaussian distribution: Each independent component is randomly generated by $C_{ijk}=\sigma_{ijk}/\sqrt{d(i,j,k)},\ (i\leq j\leq k)$. Here $\sigma_{ijk}$ is a random number following the normal distribution of mean value zero and standard deviation one. $d(i,j,k)$ is a degeneracy factor defined by

  $\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}$

  (51)

  The above distribution corresponds to choosing $\alpha=1/2$ in (5), since

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

  (52)

  due to $C$ being a symmetric tensor.

• •

  Pick up all the real eigenvectors by the way mentioned above.

• •

  Store the size $|v|$ and the value of $\det M$ for each of all the real eigenvectors.

• •

  Repeat the above processes.

By the above repeating procedure we obtain a data of $(|v_{i}|,\det M_{i})\ (i=1,2,\ldots,L)$, where $L$ is the total number of real eigenvectors generated. For the general case of $R$, we define

$\displaystyle\rho_{R}^{\rm sim}((k+1/2)\delta v)=\frac{1}{\delta vN_{C}}\sum_{i=1}^{L}|\det M|^{2R-1}\,\theta(k\delta v<|v_{i}|\leq(k+1)\delta v),$

(53)

where $N_{C}$ denotes the total number of randomly generated $C$, $\delta v$ is a bin size, $k=0,1,2,\ldots$, and $\theta$ is a support function which takes 1 if the inequality of the argument is satisfied, but zero otherwise. This quantity corresponds to $\rho(v,R,0)S_{N-1}|v|^{N-1}$ with $\alpha=1/2$ of the analytical result by a derivation similar to that of (49). The $-1$ in the exponent of $|\det M|$ in (53) comes from the consideration of the measure associated to the delta functions, namely, the difference between (2) and (3).