# Wave Function Structure in Two-Body Random Matrix Ensembles

## Abstract

We study the structure of eigenstates in two-body interaction random matrix ensembles and find significant deviations from random matrix theory expectations. The deviations are most prominent in the tails of the spectral density and indicate localization of the eigenstates in Fock space. Using ideas related to scar theory we derive an analytical formula that relates fluctuations in wave function intensities to fluctuations of the two-body interaction matrix elements. Numerical results for many-body fermion systems agree well with the theoretical predictions.

PACS numbers: 24.10.Cn, 05.45.+b, 24.60.Ky, 05.30.-d

Random matrix theory (RMT) has become a powerful tool for describing statistical properties of wave functions and energy levels in complex quantum systems . The use of a two-body random matrix ensemble (TBRE)  is of particular interest for many-body systems since classical RMT implies the presence of $`k`$-body forces ($`k>2`$) and gives the unphysical semicircle as the spectral density. The TBRE displays the same spectral fluctuations as classical RMT while its Gaussian spectral density agrees well with nuclear shell model calculations. The situation is not so clear for the structure of wave functions in TBRE. Recent results show that ground states of shell model Hamiltonians with two-body random interactions favor certain quantum numbers  and thereby differ considerably from RMT expectations.

In this letter we examine the structure of wave functions in TBRE, and compare numerical results with theoretical predictions. Such a study is not only interesting on its own but is also motivated by the ongoing importance of TBRE for nuclear  and mesoscopic physics . We recall that deviations from RMT indicate some degree of wave function non-ergodicity and are related to phenomena like Fock space localization in many-body systems , and scars of periodic orbits  or invariant manifolds  in classically chaotic systems.

To quantify the degree of localization of a given wave function, it is useful to introduce the notion of an inverse participation ratio (IPR). Thus, let $`D`$ be the total dimension of the relevant Fock subspace, let $`|b`$ label the single-particle basis states (Fock states), and let $`|\alpha `$ represent the eigenstates of the Hamiltonian. Then the overlap intensities

$$P_{\alpha b}=|\alpha |b|^2$$
(1)   
are the squares of the expansion coefficients, and have mean value $`1/D`$. The IPR of eigenstate $`|\alpha `$ is defined as the first nontrivial moment of the intensity distribution, namely the ratio of the mean squared $`P_{\alpha b}`$ to the square of the mean:

$$\mathrm{IPR}_\alpha =\frac{\frac{1}{D}\underset{b=1}{\overset{D}{}}P_{\alpha b}^2}{\left(\frac{1}{D}_{b=1}^DP_{\alpha b}\right)^2}=D\underset{b=1}{\overset{D}{}}P_{\alpha b}^2.$$
(2)   
The IPR measures the inverse fraction of Fock states that participate in building up the full wave function $`|\alpha `$, i.e. $`\mathrm{IPR}_\alpha =1`$ for a wave function that has equal overlaps $`P_{\alpha b}`$ with all basis states, and $`\mathrm{IPR}_\alpha =D`$ for the other extreme of a wave function that is composed entirely of one basis state. While complete information about wave function ergodicity is contained in the full distribution of intensities $`P_{\alpha b}`$, the IPR serves as a very useful one-number measure of the degree of Fock-space localization. In RMT, the $`P_{\alpha b}`$ are given (in the large-$`D`$ limit) by squares of Gaussian random variables, in accordance with the Porter-Thomas distribution, leading to $`\mathrm{IPR}_{\mathrm{RMT}}=3`$ for real wave functions. For finite $`D`$, the IPR is slightly below its asymptotic value (with the deviation falling off as $`1/D`$), but is still uniform over the entire spectrum.

A very simple two-body interaction model, however, already displays behavior which is qualitatively different from this naive RMT expectation. Let the Hamiltonian be given by

$$H=\frac{1}{2}\underset{i,j,k,l}{}V_{ijkl}a_i^{}a_j^{}a_ka_l,$$
(3)   
where the single-particle indices $`i`$, $`j`$, $`k`$, $`l`$ run from $`1`$ to $`M`$ (the number of available single-particle states) and the $`V_{ijkl}`$ are Gaussian random variables with unit variance. The operators $`a_j^{}`$ and $`a_j`$ create and annihilate a fermion in the single-particle state labeled by $`j`$, respectively, and obey the usual anti-commutation rules. The dimension of the $`N`$-particle Fock subspace is given by $`D=\left(\genfrac{}{}{0pt}{}{M}{N}\right)`$. We notice that this model contains no explicit one-body terms and therefore does not display Anderson-type localization effects. On physical grounds, we assume that the two-body interaction is generated by a real symmetric potential, $`V(r_1,r_2)=V(r_2,r_1)`$ and choose real single-particle wave functions. These conditions lead to the constraints $`V_{ikjl}=V_{ljki}=V_{jilk}=V_{ijkl}`$ and reduce the number of independent variables describing a given realization (Eq. 3) of the ensemble. While such correlations between matrix elements will determine factors of $`2`$ in the calculation below, they do not qualitatively affect our results. The choice of bosons rather than fermions also does not qualitatively change the localization behavior.

Fig. 1(top) shows a smoothed ensemble-averaged spectral density

$$\rho (E)=\mathrm{Tr}\delta (EH)=\underset{\alpha =1}{\overset{D}{}}\delta (EE_\alpha ),$$
(4)   
for $`M=13`$ and $`N=6`$, while in Fig. 1(bottom) we plot $`\mathrm{IPR}/\mathrm{IPR}_{\mathrm{RMT}}1`$ as a function of energy. Obviously, we observe strong deviations from RMT behavior at the edges of the spectrum. This implies that wave function intensities do not have a Porter-Thomas distribution.

To understand this surprising result, we may adapt a formalism previously used successfully to understand wave function scars and other types of anomalous quantum localization behavior in single-particle systems. Let $`\rho _b(E)`$ be the local density of states (strength function) of the basis state $`b`$:

$$\rho _b(E)=b|\delta (EH)|b=\underset{\alpha =1}{\overset{D}{}}P_{\alpha b}\delta (EE_\alpha );$$
(5)   
it is given by the Fourier transform of the autocorrelation function

$$A_b(t)=b|e^{iHt}|b.$$
(6)   
Let us assume that $`A(t)`$ displays only two different time scales: the initial decay time $`T_{\mathrm{decay}}`$ of the Fock state $`|b`$ due to interactions, and the Heisenberg time $`T_HDT_{\mathrm{decay}}`$ (i.e. $`\mathrm{}`$ over the mean level spacing) at which individual eigenlevels are resolved. Following the initial decay, random long-time recurrences in $`A_b(t)`$ can be shown to be convolved with the short-time behavior . In the energy domain this produces random oscillations $`f_b(E)`$ multiplying a smooth envelope $`\rho _b^{\mathrm{sm}}(E)`$ given by the short-time dynamics:

$$\rho _b(E)=f_b(E)\rho _b^{\mathrm{sm}}(E).$$
(7)   
Here

$$\underset{b=1}{\overset{D}{}}\rho _b^{\mathrm{sm}}(E)=\rho ^{\mathrm{sm}}(E),$$
(8)   
while $`f_b(E)`$ is a fluctuating function with mean value of unity:

$$f_b(E)=\frac{\underset{\alpha =1}{\overset{D}{}}r_{\alpha b}\delta (EE_\alpha )}{\rho ^{\mathrm{sm}}(E)}.$$
(9)   
The $`r_{\alpha b}`$ are random $`\chi ^2`$ variables with mean value one. Then, substituting into Eq. 5, we have the individual wave function intensities given by

$$P_{\alpha b}=r_{\alpha b}\frac{\rho _b^{\mathrm{sm}}(E_\alpha )}{\rho ^{\mathrm{sm}}(E_\alpha )}.$$
(10)   
Both $`r_{\alpha b}`$ and $`\rho _b^{\mathrm{sm}}`$ have $`b`$-dependent fluctuations, which are uncorrelated under our assumption of separated time scales. Then using Eq. 10 we can express the IPR (Eq. 2) as ($`\delta \rho _b^{\mathrm{sm}}(E)\rho _b^{\mathrm{sm}}(E)\rho _b^{\mathrm{sm}}(E)`$)

$`\mathrm{IPR}_\alpha `$ $`=`$ $`{\displaystyle \frac{r_{\alpha b}^2}{r_{\alpha b}^2}}\times {\displaystyle \frac{\rho _b^{\mathrm{sm}}(E_\alpha )^2}{\rho _b^{\mathrm{sm}}(E_\alpha )^2}}`$ (11)   
$`=`$ $`\mathrm{IPR}_{\mathrm{RMT}}\left(1+{\displaystyle \frac{\delta \rho _b^{\mathrm{sm}}(E_\alpha )^2}{\rho _b^{\mathrm{sm}}(E_\alpha )^2}}\right),`$ (12)   
where all averages are over the Fock basis index $`b`$, i.e. $`\mathrm{}D^1_{b=1}^D`$.

In the limit of many particles (and many holes) $`N`$, $`MN1`$ the spectrum approaches a Gaussian shape 

$$\rho ^{\mathrm{sm}}(E)=\frac{D}{\sqrt{2\pi E_0^2}}\mathrm{exp}(E^2/2E_0^2),$$
(13)   
where $`E_0^2=D^1\mathrm{Tr}H^2`$ is given by the mean sum of squares of matrix elements in any given row of the Hamiltonian. The same arguments lead to the vanishing of higher-order cumulants for the individual strength functions $`\rho _b^{\mathrm{sm}}(E)`$, so each of these should also have a Gaussian shape, but with centroid $`c_b=H_{bb}`$ and variance $`v_b=_b^{}H_{bb^{}}^2`$. Due to these fluctuations in the centroids and widths as one goes through the different basis states $`b`$, we have

$`\delta \rho _b^{\mathrm{sm}}(E)`$ $`=`$ $`{\displaystyle \frac{\rho _b^{\mathrm{sm}}(E)}{c_b}}\delta c_b+{\displaystyle \frac{\rho _b^{\mathrm{sm}}(E)}{v_b}}\delta v_b`$ (14)   
$`=`$ $`\left[{\displaystyle \frac{E}{E_0}}{\displaystyle \frac{\delta c_b}{E_0}}+\left({\displaystyle \frac{E^2}{E_0^2}}1\right){\displaystyle \frac{\delta v_b}{2E_0^2}}\right]\rho ^{\mathrm{sm}}(E),`$ (15)   
where in the second line we have used the Gaussian form of Eq. 13. Substituting into Eq. 11, we obtain the general form:

$$\frac{\mathrm{IPR}(E)}{\mathrm{IPR}_{\mathrm{RMT}}}1=\frac{(\delta c_b)^2}{E_0^2}\left(\frac{E^2}{E_0^2}\right)+\frac{(\delta v_b)^2}{4E_0^4}\left(\frac{E^2}{E_0^2}1\right)^2,$$
(16)   
valid of course not only for two-body interactions but for any ergodic Hamiltonian with a Gaussian density of states. The quantities $`(\delta c_b)^2`$ and $`(\delta v_b)^2`$ do depend on the parameters of the model, but remarkably the IPR behavior is always given to leading order by a quartic polynomial in energy. From this point of view, it is easy to understand the enhancement in fluctuations near the edge of the spectrum (i.e. for large $`|E|`$): two Gaussians differing only slightly in their centroid or width may look almost identical in the bulk, but the relative difference increases dramatically as one moves into the tail of the spectrum.

For our simple model (Eq. 3) one may easily compute the coefficients in Eq. 16. The squared width $`v_b=E_0^2`$ of the full spectrum is given by the sum of squares of entries in one row of the Hamiltonian and equals the number of independent terms in Eq. 3 that couple any given Fock state to other Fock states, times the mean squared value ($`=1`$) of each such term. A simple counting argument then shows $`E_0^2=2\left(\genfrac{}{}{0pt}{}{N}{2}\right)\left[1+2(MN)+\left(\genfrac{}{}{0pt}{}{MN}{2}\right)\right]`$. Because all the contributions are independent $`\chi ^2`$ variables of mean $`1`$ and variance $`2`$, the variance in the sum is given by twice the number of contributions: $`(\delta v_b)^2=2E_0^2`$. Finally, the variance in the centroid is given by the number of terms in the Hamiltonian (Eq. 3) which contribute to each diagonal element in the Fock basis, namely $`(\delta c_b)^2=\left(\genfrac{}{}{0pt}{}{N}{2}\right)`$. Then

$`{\displaystyle \frac{\mathrm{IPR}(E)}{\mathrm{IPR}_{\mathrm{RMT}}}}1`$ $`=`$ $`\left[1+2(MN)+\left({\displaystyle \genfrac{}{}{0pt}{}{MN}{2}}\right)\right]^1`$ (17)   
$`\times `$ $`\left[\left({\displaystyle \frac{E^2}{E_0^2}}\right)+\left[2\left({\displaystyle \genfrac{}{}{0pt}{}{N}{2}}\right)\right]^1\left({\displaystyle \frac{E^2}{E_0^2}}1\right)^2\right]`$ (18)   
$``$ $`{\displaystyle \frac{2}{M^2}}\left[{\displaystyle \frac{E^2}{E_0^2}}+{\displaystyle \frac{1}{N^2}}\left({\displaystyle \frac{E^2}{E_0^2}}1\right)^2\right],`$ (19)   
where in the last line we have taken the dilute many-particle limit $`MN1`$. The fluctuations are always strongest near the edge of the spectrum; in particular close to the ground state ($`E_{\mathrm{gs}}^22E_0^2\mathrm{ln}D2E_0^2N\mathrm{ln}(M/N)`$), we find

$$\frac{\mathrm{IPR}(E_{\mathrm{gs}})}{\mathrm{IPR}_{\mathrm{RMT}}}1=\frac{4}{M^2}\left[N\mathrm{ln}\frac{M}{N}+2\mathrm{ln}^2\frac{M}{N}\right].$$
(20)   
Let us compare the IPR prediction of Eq. 18 with the numerical data presented in Fig. 1(bottom). While the prediction qualitatively describes the correct trend of the IPR as a function of energy, it is not in quantitative agreement with the data and fails by as much as a factor of $`2`$ at the edge of the spectrum. This is not very surprising, given that the actual behavior of the spectral density $`\rho ^{\mathrm{sm}}(E)`$ is far from Gaussian for the parameters we have chosen. Although it is known that the spectrum approaches a Gaussian form in the many-particle limit , the number of particles that can be simulated numerically is not nearly large enough for the Gaussian form to be a good quantitative approximation, particularly in the tail (see Fig. 1(top)). A good way to measure deviations from the Gaussian form is to compute the fourth cumulant divided by the square of the second cumulant $`E_0^2`$: for a Gaussian this quantity vanishes while for our system it equals $`0.7`$, actually closer to the semicircle value of $`1`$.

To obtain quantitatively valid predictions, we must correct for deviations from the Gaussian shape. The ansatz

$$\rho ^{\mathrm{sm}}(E)=\frac{16ϵ\left(\frac{E}{E_0}\right)^2+ϵ\left(\frac{E}{E_0}\right)^4}{13ϵ}\frac{D}{\sqrt{2\pi E_0^2}}\mathrm{exp}\left(\frac{E^2}{2E_0^2}\right)$$
(21)   
allows for nonzero higher cumulants while keeping the number of states and the width $`E_0^2`$ fixed. The parameter $`ϵ`$ is obtained by a least squares fit of the numerical data to this form. A derivation starting from the first line of Eq. 14 yields an IPR that depends on the energy through a rational function instead of the quartic polynomial of Eq. 18. As we can see in Fig. 1(bottom), this leads to a surprisingly good quantitative prediction for the IPR behavior, given that we are applying perturbation theory around a Gaussian shape for a spectrum that in reality is very far from the Gaussian limit.

Gaussian spectral densities have been observed in nuclear shell model calculations with realistic  or random two-body interactions . This is mainly due to the presence of spin and isospin degrees of freedom, which cause important correlations between Hamiltonian matrix elements. Let us therefore consider adding spin to the Hamiltonian (Eq. 3), so that $`a_j^{}`$ and $`a_j`$ now create and annihilate a fermion in the single-particle state labeled by $`j(n_j,s_j)`$, with $`n_j`$ and $`s_j=\pm 1/2`$ denoting the orbital and spin quantum numbers. In what follows we consider $`N`$ fermions in a shell of fixed total spin $`S=_{j=1}^Ns_j=0`$. We assume that the random matrix elements $`V_{ijkl}`$ depend only on the orbital quantum numbers; this reduces the number of independent matrix elements considerably. The (sparse) two-body interaction matrix is constructed using a code similar to the one described in Ref..

The spectral density of this TBRE agrees very well with a Gaussian, and the IPR is predicted by Eq. 16. However, the determination of the quantities $`(\delta c_b)^2`$ and $`(\delta v_b)^2`$ is more difficult since the spin degree of freedom makes the required counting of matrix elements a non-trivial task. Alternatively, one may obtain these quantities directly from the numerically generated Hamiltonians. Fig. 2 shows that the numerical data and theoretical results are in good agreement. This confirms the validity of the simple theory derived in this letter.

We note that our results are in qualitative agreement with nuclear shell model calculations using realistic interactions  and with calculations in atomic physics . Our calculation shows that deviations from ergodicity near the ground state do not require the presence of one-body terms in the Hamiltonian, though such terms may of course enhance the degree of localization.

In summary we have studied numerically and analytically the wave function structure in many-body fermion systems with random two-body interactions. Near the edge of the spectrum, wave function intensities of this two-body random ensemble exhibit fluctuations that deviate strongly from random matrix theory predictions, while good agreement is obtained in the bulk of the spectrum. The numerical results agree well with the theoretical prediction that is derived from arguments used in scar theory. In particular, we have presented a simple formula that relates fluctuations of the wave function intensities to fluctuations of the two-body matrix elements of the Hamiltonian.

We thank George F. Bertsch for suggesting this study and for several very useful discussions. This research was supported by the DOE under Grant DE-FG-06-90ER40561.