Energy level splitting for weakly interacting bosons in a harmonic trap

TL;DR

This paper investigates how weak contact interactions cause energy level splitting in two-dimensional trapped bosons, revealing patterns related to quantum chaos and symmetries, with implications for understanding turbulence and quantum integrability.

Contribution

It provides a detailed analysis of energy level splitting patterns for weakly interacting bosons in a harmonic trap, connecting quantum resonant systems to classical turbulence studies.

Findings

Energy levels split and form patterns related to symmetries.

Wigner-Dyson statistics observed in level spacings.

Explicit description of the highest energy states from each unperturbed level.

Abstract

We consider identical quantum bosons with weak contact interactions in a two-dimensional isotropic harmonic trap. When the interactions are turned off, the energy levels are equidistant and highly degenerate. At linear order in the coupling parameter, these degenerate levels split, and we study the patterns of this splitting. It turns out that the problem is mathematically identical to diagonalizing the quantum resonant system of the two-dimensional Gross-Pitaevskii equation, whose classical counterpart has been previously studied in the mathematical literature on turbulence. Our purpose is to explore the implications of the symmetries and energy bounds of this resonant system, previously studied for the classical case, for the quantum level splitting. Simplifications in computing the splitting spectrum numerically result from exploiting the symmetries. The highest energy state emanating from each unperturbed level is explicitly described by our analytics. We furthermore discuss the energy level spacing distributions in the spirit of quantum chaos theory. After separating the eigenvalues into blocks with respect to the known conservation laws, we observe the Wigner-Dyson statistics within specific large blocks, which leaves little room for further integrable structures in the problem beyond the symmetries that are already explicitly known.

Equations

$$\mathcal{H}={\mathcal{H}}_0+g{\mathcal{H}}_{int},$$

$${\mathcal{H}}_0=\frac{1}{2}\int (\nabla {\Psi }^{†}\cdot \nabla \Psi +(x^2+y^2){\Psi }^{†}\Psi )dxdy,$$

$${\mathcal{H}}_{int}=\pi \int {\Psi }^{†2}{\Psi }^{2}dxdy.$$

$$[{\Psi }^{†}(x,y),\Psi (x^{\prime },y^{\prime })]=-\delta (x-x^{\prime })\delta (y-y^{\prime }).$$

$$\Psi (x,y)=\sum _{n,m}{\alpha }_{nm}{\psi }_{nm}(x,y),$$

$$\frac{1}{2}(-$$

$$-i{\partial }_{\phi }{\psi }_{nm}=m{\psi }_{nm},$$

$$[{\alpha }_{nm}^{†},{\alpha }_{n^{\prime }{m}^{\prime }}]=-{\delta }_{n{n}^{\prime }}{\delta }_{m{m}^{\prime }}.$$

$${\mathcal{H}}_0=\sum _{nm}n{\alpha }_{nm}^{†}{\alpha }_{nm},$$

$$|\{{\eta }_{nm}\}⟩=\prod _{nm}\frac{({\alpha }_{nm}^{†}{)}^{{\eta }_{nm}}}{\sqrt{{\eta }_{nm}!}}|0⟩.$$

$${\alpha }_{nm}^{†}{\alpha }_{nm}|\{\eta \}⟩={\eta }_{nm}|\{\eta \}⟩,$$

$$E_{\{\eta \}}=\sum _{nm}n{\eta }_{nm}.$$

$${\psi }_{nm}^{\text{mbox}norm.}=\sqrt{\frac{(\frac{1}{2}(n-|m|))!}{(\frac{1}{2}(n+|m|))!}}\frac{r^{|m|}}{\sqrt{\pi }}{L}_{\frac{n-|m|}{2}}^{|m|}(r^2)e^{-r^2/2}{e}^{im\varphi },$$

$${\psi }_{nm}=(-1{)}^{\frac{1}{2}(m-|m|)}{\psi }_{nm}^{\text{mbox}norm.}.$$

$$L_n^{\mu }(\rho )=\sum _{k=0}^n(-1{)}^k\frac{(n+\mu )!}{k!(n-k)!(k+\mu )!}{\rho }^k,$$

$$L_n^{\mu }(\rho )=(-1{)}^{\mu }\frac{(n+\mu )!}{n!}{\rho }^{-\mu }{L}_{\mu +n}^{-\mu }(\rho ).$$

$$\sqrt{\frac{((n-|m|)/2)!}{((n+|m|)/2)!}}{r}^{|m|}{L}_{\frac{n-|m|}{2}}^{|m|}(r^2)e^{-r^2/2}$$

$$=(-1{)}^{\frac{1}{2}(m-|m|)}\sqrt{\frac{((n-m)/2)!}{((n+m)/2)!}}{r}^m{L}_{\frac{n-m}{2}}^m(r^2)e^{-r^2/2},$$

$${\psi }_{nm}=\sqrt{\frac{(\frac{1}{2}(n-m))!}{(\frac{1}{2}(n+m))!}}\frac{r^m}{\sqrt{\pi }}{L}_{\frac{n-m}{2}}^m(r^2)e^{-r^2/2}{e}^{im\varphi }.$$

$${\mathcal{H}}_{int}=$$

$$C_{n_1{n}_2{n}_3{n}_4}^{m_1{m}_2{m}_3{m}_4}=2\pi \int {\psi }_{n_1{m}_1}^{\ast }{\psi }_{n_2{m}_2}^{\ast }{\psi }_{n_3{m}_3}{\psi }_{n_4{m}_4}rdrd\varphi$$

$$C_{0000}^{0000}=1.$$

$$⟨\{\eta \}|{\mathcal{H}}_{int}|\{{\eta }^{\prime }\}⟩,$$

$$N_{\{\eta \}}=\sum _{nm}{\eta }_{nm}.$$

$${\tilde{E}}_I=E+g{\epsilon }_I.$$

$${\mathcal{H}}_{res}=\frac{1}{2}\sum _{\begin{array}{c}{n}_1+n_2=n_3+n_4\\ m_1+m_2=m_3+m_4\end{array}}{C}_{n_1{n}_2{n}_3{n}_4}^{m_1{m}_2{m}_3{m}_4}{\alpha }_{n_1{m}_1}^{†}{\alpha }_{n_2{m}_2}^{†}{\alpha }_{n_3{m}_3}{\alpha }_{n_4{m}_4}.$$

$${\mathcal{H}}_{res}=\frac{1}{2}\sum _{\begin{array}{c}{n}_1+n_2=n_3+n_4\\ m_1+m_2=m_3+m_4\end{array}}{C}_{n_1{n}_2{n}_3{n}_4}^{m_1{m}_2{m}_3{m}_4}{\alpha }_{n_1{m}_1}^{\ast }{\alpha }_{n_2{m}_2}^{\ast }{\alpha }_{n_3{m}_3}{\alpha }_{n_4{m}_4}$$

$$i\frac{d{\alpha }_{nm}}{dt}=\sum _{\begin{array}{c}n+n_1=n_2+n_3\\ m+m_1=m_2+m_3\end{array}}{C}_{n{n}_1{n}_2{n}_3}^{m{m}_1{m}_2{m}_3}{\alpha }_{n_1{m}_1}^{\ast }{\alpha }_{n_2{m}_2}{\alpha }_{n_3{m}_3}.$$

$$N$$

$$E$$

$$M$$

$$Z_{+}$$

$$Z_{-}$$

$$W$$

$${\mathcal{H}}_{LLL}=\frac{1}{2}\sum _{\begin{array}{c}{n}_1,n_2,n_3,n_4=0\\ n_1+n_2=n_3+n_4\end{array}}^{\infty }{C}_{n_1{n}_2{n}_3{n}_4}{\alpha }_{n_1}^{†}{\alpha }_{n_2}^{†}{\alpha }_{n_3}{\alpha }_{n_4}.$$

Full text

I Introduction

Energy levels of quantum identical interacting bosons in harmonic traps have often been studied in the literature split1 ; split2 ; split3 ; split4 ; split5 ; split6 ; split7 ; split8 ; split9 ; breathing , motivated in particular by the physics of cold atomic gases. Contact interactions between the bosons commonly appear in such studies as the simplest possible choice for the two-particle potential that is expected to retain realistic features. Our purpose in this article is to report on rich mathematical structures emerging from studies of bosons with weak contact interactions in a two-dimensional isotropic harmonic trap, and the consequences these structures have for the makeup of the energy spectrum.

Most of the past publications dealing with related systems have focused on particular small numbers of bosons. For instance, split1 reports an exact solution of the problem involving two bosons in a harmonic trap with contact interactions of arbitrary strength. A notable exception is breathing , which studies effects of symmetries on the spectrum of an arbitrary number of atoms with finite strength contact interactions in a harmonic trap. Our perspective will be quite similar, but the strength of the interactions will be assumed small (more precisely, our results are valid at linear order in the coupling parameter). This allows for more extensive analytic exploration of the structure of the spectrum.

If the interactions are turned off altogether, the eigenstates of the multi-boson system are simply constructed by having the bosons independently occupy the energy levels of the harmonic potential, and the total energy is a sum of contributions of the individual noninteracting bosons. Since the energies of harmonic potential eigenstates are integer in appropriate units, one ends up with integer energy eigenvalues for the multi-boson system, and the degeneracies of the energy levels are very high, since there are many ways to partition a given integer amount of energy into the integer energies of individual bosons. (These level degeneracies furthermore grow without bound as the energy is increased.) When weak contact interactions between the bosons are turned on, the highly degenerate free-boson energy levels split, and the pattern of the splitting at leading order in the interaction strength is computed by diagonalizing the interaction Hamiltonian within each degenerate unperturbed energy level, which is a standard quantum-mechanical procedure. These level splitting patterns will be the main subject of our analysis.

If one is interested in the splitting pattern of a given unperturbed level, diagonalizing a finite-sized matrix of the interaction Hamiltonian matrix elements is all one has to do. This is straightforward for numerical evaluation, yet little can be said about the general structure of the spectrum, since one matrix has to be diagonalized for each unperturbed level, and the sizes of the matrices (determined by the unperturbed level degeneracies) grow without bound as one moves to higher energies. However, specifically for bosons with contact interactions in an isotropic two-dimensional harmonic trap, the diagonalization problems involved display rich symmetry structures that relate the energy shifts of different unperturbed levels, and impose significant constraints on the patterns in the spectrum.

A key observation that makes our analysis possible is that the diagonalization problem involved in finding the energy shifts is identical to diagonalizing the quantized version of the resonant system of the two-dimensional Gross-Pitaevskii equation in an isotropic harmonic trap. The Gross-Pitaevskii equation (which physically describes the condensed regime of our interacting bosons) is a partial differential equation (PDE) that can be approximated in the weakly nonlinear regime by the corresponding resonant Hamiltonian system, an approach commonly taken in the mathematical literature on turbulence FGH ; continuous . The enhanced symmetries of this classical resonant system have been presented explicitly in the past, see continuous . Quantum resonant systems, on the other hand, have been considered in quantres from a perspective geared toward quantum chaos studies. Our strategy here is to translate the classical symmetries and energy bounds of continuous to the corresponding quantum resonant system, and to relate the results to the energy level splitting of bosons with weak contact interactions.

Another issue that naturally comes to mind is whether there are symmetries in the problem beyond the ones we use. To shed light on this question, we turn to quantum chaos theory GMW ; haake ; DKPR , which purports that integrable and chaotic systems are characterized by qualitatively different distributions of distances between neighboring energy levels. For integrable systems, the distribution is Poissonian btint , so that, as far as the level spacing statistics is concerned, the energy levels appear random and completely uncorrelated. Chaotic systems, on the other hand, are associated with a phenomenon known as level repulsion, captured by the Wigner-Dyson distribution of the level spacings BGS . In our case, we already know a few explicit operators that commute with the Hamiltonian, and one must factor this knowledge into the study of the level spacing statistics. If the statistics is built indiscriminately for all energy levels, one observes the Poisson distribution typical of integrable systems. This is merely a reflection of the fact that eigenvalues split into blocks with respect to the conserved quantities, and without correlations between eigenvalues in the different blocks, the joint level spacing statistics mimics that of random energy levels DKPR . What is more relevant is to plot the level spacing distributions separately within such blocks. Once we do that taking into account all the conserved quantitites of continuous , most of the large blocks display Wigner-Dyson-like, rather than Poissonian, distributions. This suggests that no integrability is to be expected, and further analytic structures, if any, must be very subtle. (We do, however, find Poisson-like distributions in a very restricted part of the spectrum characterized by rapid rotation.)

The paper is organized as follows: In section II, we review bosons with contact interactions in harmonic traps, and their energy level splitting at linear order in the coupling parameter, and then explain the connection between these considerations and quantum resonant systems. In section III, we demonstrate the use of our techniques by focusing on specific multiplets in the so-called Lowest Landau Level (LLL) sector, for which the analysis simplifies. In section IV, we generalize the analysis of section III to arbitrary energy levels. We conclude with a summary and discussion.

II Trapped bosons, energy

shifts and resonant systems

We consider the second-quantized representation of identical bosons with contact interactions of strength $g$ in an isotropic harmonic potential:

[TABLE]

Here, $\Psi(x,y)$ is a quantum nonrelativistic bosonic field satisfying the commutation relations

[TABLE]

We shall be focusing on small values of the coupling parameter, $g\ll 1$. It will be convenient to further assume $g\geq 0$ to simplify the wording (for example, to refer to the highest rather than lowest energy state within each fine structure multiplet), though our derivations are equally valid for small negative couplings. The factor of $\pi$ is inserted on the right-hand side of (3) for future convenience and may be absorbed, if desired, into a redefinition of $g$.

We note that known subtleties exist with defining the operator product at the same point in (3), but these subtleties will not be relevant for our treatment. The issue is commonly stated in the language of first quantization, where the naive product given in (3) corresponds to the interparticle potential $V({\bf r}_{i}-{\bf r}_{j})\sim\delta({\bf r}_{i}-{\bf r}_{j})$ between particles number $i$ and $j$ of an $N$-body system. It is known that wavefunctions of multiparticle systems with contact interactions possess singularities at coincident particle positions ${\bf r}_{i}={\bf r}_{j}$, originating from the singularities of the Laplacian Green’s function in more than one spacial dimension. Since one has to define products of the wave functions and the potential to formulate the Schrödinger equation, naive $\delta$-functions in the potential are not acceptable and need to be replaced by a modification that can be multiplied by wavefunctions with specific singularities at coincident particle positions. This is clearly explained in split1 for two particles in three spatial dimensions, where the wavefunctions are allowed to have $1/r$ singularities and the $\delta({\bf r})$ has to be replaced by the operator $\delta({\bf r})\partial_{r}r$. Note that, if applied to any nonsingular function, this operator acts exactly like a naive $\delta$-function, but its action is also defined for functions with $1/r$ singularities. This also explains why this subtlety is irrelevant for our considerations: as we shall deal with the first order of perturbation theory in $g$, we only need to compute the matrix elements of the interaction Hamiltonian between noninteracting eigenstates. But the noninteracting eigenstates are smooth functions, without any singularities at coincident particle positions, and hence using naive $\delta$-functions, and naive products of second-quantized fields in (3), is just as good as using the correct regulated expressions, and evidently we shall not encounter any divergences in evaluating matrix elements of the naive expression (3).

Setting for a moment $g=0$ in (1), one obtains noninteracting bosons in a harmonic trap, and the energy spectrum descends directly from the one-particle spectrum by a simple addition of energies. One can decompose $\Psi$ in terms of the harmonic oscillator eigenfunctions $\psi_{nm}$ as

[TABLE]

with $\psi_{nm}$ satisfying

[TABLE]

where $\varphi$ is the polar angle in the $(x,y)$-plane. These states thus carry $n+1$ units of energy and $m$ units of angular momentum, with $n$ being a nonnegative integer and $m\in\{-n,-n+2,\dots,n-2,n\}$. The creation-annihilation operators $\alpha^{\dagger}_{nm}$ and $\alpha_{nm}$ satisfy the standard commutation relations

[TABLE]

The noninteracting Hamiltonian $\mathcal{H}_{0}$ is then expressed as

[TABLE]

where we have subtracted the irrelevant vacuum energy contribution by replacing $n+1$ with $n$. The eigenstates of $\mathcal{H}_{0}$ are simply given by the Fock basis states, generated by acting on the vacuum state $|0\rangle$ (one has $\alpha_{nm}|0\rangle=0$) with the creation operators $\alpha^{\dagger}_{nm}$ to produce a state with the set of occupation numbers $\{\eta_{nm}\}$, one for each one-particle mode labelled by $n$ and $m$:

[TABLE]

These satisfy, for any $n$ and $m$,

[TABLE]

and hence they are eigenstates of (8) with eigenvalues

[TABLE]

There are of course many ways to generate the same value of $E$ by combining the integer numbers $n$ and $\eta_{nm}$. Correspondingly, the energy levels of (8) are highly degenerate (and the degeneracies grow without bound as one moves to higher energies).

A note is in order on the set of eigenfunctions to be used in the decomposition (5). A common choice for the orthonormal set of two-dimensional harmonic oscillator energy eigenfunctions is given (in polar coordinates) by

[TABLE]

see, e.g., dahl . Here, $L_{n}^{\alpha}$ are the generalized Laguerre polynomials. These functions are, of course, defined up to phase factors, which are completely irrelevant in the noninteracting theory. When considering interactions, however, different choices of the phase factors in (12) will lead to equivalent theories, but may be more or less efficient from the standpoint of simplifying the algebraic expressions and making their structure more apparent. For the purposes of our algebra it turns out beneficial to introduce an extra sign factor in the definition of the eigenfunctions as follows:

[TABLE]

With this sign factor, the wavefunctions can be conveniently transformed in a manner identical to the considerations of BBCE2 using identities for Laguerre polynomials. Using the explicit expressions for the Laguerre polynomials,

[TABLE]

and remembering that factorials of negative numbers are infinite, we obtain for every integer $\mu$

[TABLE]

Therefore, one can see that introducing the extra sign factors amounts to using eigenfunctions that can be written as the original $\psi^{\mbox{\tiny norm.}}_{nm}$ but without absolute values,

[TABLE]

or

[TABLE]

In the following, we shall use these expressions in the decompositions (5), which in fact brings us in accord with the conventions of continuous and lets us conveniently reuse the mathematical structures developed there.

We are now in a position to study weak contact interactions. To this end, we substitute the decomposition (5) in the interaction Hamiltonian (3) to obtain

[TABLE]

The sum only contains terms satisfying $m_{1}+m_{2}=m_{3}+m_{4}$ as an immediate consequence of the angular momentum conservation by contact interactions. We have furthermore introduced the interaction coefficients defined by

[TABLE]

whose properties will play a key role in our analysis. Note that, with the conventions we have adopted,

[TABLE]

A standard approach to small perturbations of quantum dynamics in confining potentials, known as the Rayleigh-Schrödinger perturbation theory, is to analyze the corrections to energy levels and eigenstates as power series in the interaction strength $g$. If the unperturbed levels are degenerate, as they are in the case at hand, level splitting will be in general induced by perturbations. At linear order in $g$, this level splitting is analyzed by computing the eigenvalues $\varepsilon_{I}$ of the matrix

[TABLE]

where $|\{\eta\}\rangle$ and $|\{\eta^{\prime}\}\rangle$ are two unperturbed Fock states of the form (9) with the same value of the unperturbed energy given by (11). Since $\mathcal{H}_{int}$ conserves the number of particles, $|\{\eta\}\rangle$ and $|\{\eta^{\prime}\}\rangle$ must also contain the same total number of particles

[TABLE]

For any given $N$ and $E$, (21) is a finite-sized numerical matrix with the entries expressed through the interaction coefficients (19). Once its eigenvalues $\varepsilon_{I}$ have been found, the corresponding perturbed energy levels are given at order $g$ by

[TABLE]

We note in passing that the structure of the unperturbed levels, and the diagonalization problem arising here at linear order in $g$, parallel closely what one would have encountered if treating quantum relativistic interacting fields in Anti-de Sitter spacetime (a brief summary can be found in madagascar ). This is not a coincidence, since nonrelativistic bosonic fields in harmonic potentials arise systematically through taking nonrelativistic limits of field systems in Anti-de Sitter spacetime BEL ; BEF . (Energy levels of quantum interacting fields in Anti-de Sitter spacetime have recently been considered from a different perspective in BSS .)

Because the energy carried by $|\{\eta\}\rangle$ and $|\{\eta^{\prime}\}\rangle$ in (21) is the same, the two annihilation operators in $\mathcal{H}_{int}$ must remove the same total amount of energy as what the two creation operators add, i.e., only terms with $n_{1}+n_{2}=n_{3}+n_{4}$ in (18) may contribute in the matrix elements (21). This results in a simplification that is straightforward, but has welcome analytic consequences. Namely, for the purposes of computing the shifted energy levels (23), one may replace $\mathcal{H}_{int}$ in (21) by $\mathcal{H}_{res}$ defined by

[TABLE]

We note that the classical system corresponding to (24) is described by the Hamiltonian

[TABLE]

for complex-valued dynamical variables $\alpha_{nm}(t)$ and $\alpha^{*}_{nm}(t)$ with the symplectic form $i\sum_{nm}d\alpha^{*}_{nm}\wedge d\alpha_{nm}$. The corresponding equations of motion are

[TABLE]

This classical Hamiltonian has been studied in the literature continuous ; BBCE2 as the resonant approximation to the Gross-Pitaevskii equation, which is the classical limit of (1) – see also FGH where the same resonant Hamiltonian emerges from approximating a different related PDE. Similar resonant Hamiltonian systems emerge as approximations to other physically motivated PDEs, see BEL ; BEF ; CF ; BHP ; FPU ; CEV ; BMR . It is not surprising that a close relation exists between applying the resonant approximation in a classical theory and the first order of the Rayleigh-Schrödinger perturbation theory for the quantum version of the same problem, since both approaches succeed in approximating the time evolution of the perturbed system on time scales of order $1/g$ at small values of the coupling parameter $g$.

Studies of the classical resonant system (25-26) have produced a large set of conserved quantities (a list can be found in continuous ). Since the conserved quantities are bilinear in $\alpha_{nm}$ and $\alpha^{*}_{nm}$, their quantization is straightforward and does not incur ordering ambiguities. This results in a set of operators commuting with the quantum Hamiltonian (24) that can be presented in terms of the following combinations

[TABLE]

The last three operators are not Hermitian, and their Hermitian conjugates $Z_{+}^{\dagger}$, $Z_{-}^{\dagger}$ and $W^{\dagger}$ should also be included. We do not explicitly give the summation ranges of $n$ and $m$, but it is understood, here and elsewhere, that summations run over all possible values of $n$ and $m$ for which the creation-annihilation operators in the summand correspond to existing modes. As a reminder, $\alpha_{nm}$ corresponds to an actual oscillator mode if $n\geq 0$, $|m|\leq n$, and $n-m$ is even.

While the conservation of $N$, $E$ and $M$ is straightforwardly seen from the general structure of the resonant Hamiltonian (24), the conservation of $Z_{+}$, $Z_{-}$ and $W$ relies on the specific form of the interaction coefficients $C$ given by (19) and is not immediately obvious. These conservation laws were established in FGH where the resonant system (25-26) was derived for bosons with contact interactions without a harmonic trap. The reason the same resonant system is of relevance both with and without a harmonic trap is that, specifically in two dimensions and for contact interactions, the two situations are mapped into each other using the so-called ‘lens transform,’ also known as the pseudo-conformal compactification Carles ; Tao . Thus, as explicitly pointed out in continuous , the same conservation laws are respected in the presence of a harmonic trap. As the original derivations of the conservation laws in the mathematical literature would have taken us rather far outside our present context, we feel it beneficial to present elementary proofs, which we give in Appendix A. The essence of these proofs is that the harmonic oscillator mode functions (17) are expressed through the Laguerre polynomials, and identities for the Laguerre polynomials imply that the interaction coefficients defined by (19) satisfy certain finite difference equations with respect to their mode number indices. These finite difference equations, in turn, imply the conservation of $Z_{+}$, $Z_{-}$ and $W$.

What do the conserved quantities (27-32) tell us about the diagonalization of (24), and consequently about the energy shifts (23)? First of all, since $N$, $E$ and $M$ commute not only with $\mathcal{H}_{res}$ but also with each other, the four operators can be diagonalized simultaneously, grouping the eigenvalues of $\mathcal{H}_{res}$ into $(N,E,M)$-blocks, labelled by the integer eigenvalues of $N$, $E$ and $M$ (which equal $\sum\eta_{nm}$, $\sum n\,\eta_{nm}$ and $\sum m\,\eta_{nm}$, respectively). In fact, since $N$, $E$ and $M$ are diagonal in the Fock basis, the block-diagonal structure of $\mathcal{H}_{res}$ is seen directly by considering its Fock basis matrix elements. Evidently, there are finitely many states for a given triplet $(N,E,M)$. Each such block describes energy shifts (23) of the unperturbed level of energy $E$ in the sector with $N$ particles carrying a total of $M$ units of angular momentum. Now, the action of $Z_{+}$ moves any state in an $(N,E,M)$-block to a state in the $(N,E+1,M+1)$-block, the action of $Z_{-}$ moves any state in an $(N,E,M)$-block to a state in the $(N,E+1,M-1)$-block, and the action of $W$ moves any state in an $(N,E,M)$-block to a state in the $(N,E+2,M)$-block, as depicted in Fig. 1.

Since $Z_{+}$, $Z_{-}$ and $W$ commute with $\mathcal{H}_{res}$, their action will copy the eigenvalues from lower to higher $(N,E,M)$-blocks. A distinctive signature of this structure is that many energy shifts $\varepsilon_{I}$ in (23) will be identical for different unperturbed energy levels, and thus infinitely many integer energy differences of the noninteracting boson spectrum will survive at linear order in $g$. This is in fact a reflection of the general pattern valid for any $g$, translated into the resonant system (24). Indeed, the center-of-mass motion is known to separate for bosons with arbitrary translation-invariant interactions confined in a harmonic trap (see, e.g., bbbb ). The energy of the center-of-mass, given by an energy eigenvalue of a two-dimensional harmonic oscillator, is simply added to the energy of the relative motion. The two-dimensional harmonic oscillator admits two independent raising operators increasing the energy value by 1 unit, which corresponds in the language of the resonant system to the action of $Z_{+}$ and $Z_{-}$. Additional symmetry enhancement exists specifically for bosons with contact interactions in two spatial dimensions, where the symmetry is extended to the full Schrödinger group Niederer ; OFN . This results in the existence of the Pitaevskii-Rosch breathing mode and the corresponding independent raising operator increasing the energy by 2 units breathing . This corresponds to the action of $W$. We note that the raising operators of the original system are not symmetries in the standard Hamiltonian sense (they have nontrivial commutation relations with the Hamiltonian), but they are translated into ordinary Hamiltonian symmetries (operators commuting with the Hamiltonian) in the language of the resonant system (24).

While the symmetries of the resonant system (24) originate from the known symmetry structures of the original problem (1), their operation has more powerful consequences. Indeed, under normal circumstances, symmetries allow generating new energy eigenstates from a known energy eigenstate, but do not simplify the process of finding additional energy states outside a given symmetry multiplet. In our case, the spectrum of (24) splits into finite-sized blocks (originating from the degenerate energy levels of noninteracting bosons), which allows for computing parts of the spectrum exactly, independently of other parts of the spectrum. Then, for example, instead of computing all eigenvalues in blocks up to a certain energy, one can simply diagonalize a single sufficiently large block, and then recover the lower blocks by repeatedly acting with $Z_{\pm}^{\dagger}$ and $W^{\dagger}$ on the explicitly found eigenvectors. (Alternatively, one can diagonalize operators like $Z_{+}Z_{+}^{\dagger}$ simultaneously with $\mathcal{H}_{res}$, and this will allow for recovering energy eigenvalues from lower blocks, as we shall briefly describe in section IV.) Separation of energy eigenstates with respect to the action of the symmetries (30-32) will also play a crucial role in our analysis of the level spacing statistics. Furthermore, a number of mathematical results are known for the classical counterpart (25) of our system (24), such as bounds on the classical Hamiltonian, and they have implications for the quantum problem we are considering. Our goal for the rest of this treatment will be to systematically explore these issues, and to comment on the eigenvalue statistics in hope of addressing the question whether the set of symmetries we have used is complete. As considerations for general $(N,E,M)$-blocks become rather involved, we shall start by focusing on blocks with $E=M$, known as the Lowest Landau Level (LLL), where the algebra simplifies considerably, and all essential concepts and techniques may be more transparently demonstrated. We shall thereafter proceed with our analysis of general $(N,E,M)$-blocks.

III The LLL sector

III.1 Classical and quantum LLL truncation

It is well-known continuous ; GT ; BBCE ; GGT that the classical resonant system (25-26) can be consistently truncated to the set of modes satisfying $n=m$, the so-called Lowest Landau Level (LLL) sector. Setting all modes with $n\neq m$ to zero in the initial state guarantees that they never get excited by the evolution equation (26). (The same set of LLL modes has often appeared in more phenomenological studies of rapidly rotating trapped Bose-Einstein condensates fetter .)

Classical truncations of this sort in general have no direct implications for the quantum theory, since the uncertainty principle makes it impossible to set canonical variables to zero. The LLL truncation happens to have a direct translation to the quantum system (24), nonetheless, for reasons that are essentially kinematical. Consider an $(N,E,M)$-block of Fock states with $N$ particles, $E$ units of energy and $M$ units of angular momentum, and impose $E=M$. Then, the occupation numbers $\eta_{nm}$ in the corresponding Fock states must vanish for $n\neq m$. Indeed, only modes with $m\leq n$ exist, and having any modes with $m<n$ excited guarantees that the total angular momentum $M$ is less than the total energy $E$, in contradiction with our assumption. Thus, any block with $E=M$ is entirely composed of LLL states. Creation-annihilation operators corresponding to non-LLL modes do not contribute to matrix elements (21) between any two states in an $E=M$ block. We shall refer to the $E=M$ blocks as ‘LLL blocks’ for obvious reasons.

In view of the above picture, for the analysis of level splitting in the LLL blocks, the full resonant Hamiltonian (24) can be replaced by the following simpler LLL Hamiltonian

[TABLE]

Here, we have relabelled $\alpha_{n_{i}m_{i}}$ with $n_{i}=m_{i}$ as simply $\alpha_{n_{i}}$. The LLL interaction coefficients can be directly evaluated continuous as a special case of (19) resulting in the following simple expression:

[TABLE]

The classical system corresponding to the LLL Hamiltonian (33) has been studied as an approximation to the Gross-Pitaevskii equation for Bose-Einstein condensates continuous ; GT ; BBCE ; GGT , and possesses many special properties, for example it admits an invariant manifold where the nonlinear equations can be solved exactly, and the solutions show interesting long-term return behaviors BBCE . Spatial positions of zeros of the wavefunctions (known as ‘vortices’) corresponding to some solutions may also be analyzed to a great extent GGT . This system is a representative of a very large class of partially solvable resonant systems developed in AO , which are of the form (33) but with different choices of the interaction coefficients $C$, and which share many of the special dynamical features we have just mentioned.

The quantum Hamiltonian (33-34), which is what is of interest for us here, inherits the following conserved quantities from (27-32):

[TABLE]

The commutators of these operators with $\mathcal{H}_{LLL}$ vanish, and the commutators among themselves and with $Z^{\dagger}$ vanish except for

[TABLE]

Energy spectra of systems of the form (33) with general interaction coefficients $C$ have been studied numerically in quantres . We shall now examine how the specific symmetry structures emerging for $C$ given by (34) influence the eigenvalue patterns.

III.2 Structure of the Hamiltonian blocks

The structure of the eigenvalue problem for the LLL Hamiltonian (33) is inherited as a simplified version from what has been described in the previous section for the full resonant Hamiltonian (24). On the other hand, the generalities of diagonalizing Hamiltonians of the form (33) with arbitrary interaction coefficients $C$ have been spelled out in quantres , and we shall closely follow that treatment.

One starts with defining the LLL Fock basis as

[TABLE]

such that

[TABLE]

for any $k$, and $\eta_{k}$ are nonnegative integers. The eigenvalues of $N$ and $E$ in this basis are evidently

[TABLE]

The Hamiltonian (33) has nonvanishing matrix elements between $|\eta_{0},\eta_{1},\dots\rangle$ and $|\eta^{\prime}_{0},\eta^{\prime}_{1},\dots\rangle$ only for two sets of occupation numbers $\{\eta_{k}\}$ and $\{\eta^{\prime}_{k}\}$ having the same values of $N$ and $E$. Thus, the Hamiltonian is block-diagonal in the Fock basis, where the blocks are labeled by the nonnegative numbers $(N,E)$. The number of states in an $(N,E)$-block is given by the number of integer partitions of $E$ into at most $N$ parts quantres , a well known number-theoretic function usually denoted as $p_{N}(E)$. One thus has to diagonalize finite-sized $p_{N}(E)\times p_{N}(E)$ matrices to get the eigenvalues of (33) within each $(N,E)$-block.

Up to this point, our discussion of the diagonalization has been generic and did not make any reference to the specific form of the interaction coefficients given by (34). For this specific form of the interaction coefficients, an extra conserved operator $Z$ given by (35) and its Hermitian conjugate enter the game. These operators act as raising and lowering operators for $E$: namely, acting on a state in an $(N,E)$-block, $Z$ produces a state in the $(N,E+1)$-block, and $Z^{\dagger}$, a state in the $(N,E-1)$-block. Since $Z$ and $Z^{\dagger}$ commute with $\mathcal{H}_{LLL}$, eigenvectors of the latter are mapped into eigenvectors by this action, while their eigenvalues remain intact.

We are thus brought to the first qualitative conclusion of our analysis. By the action of $Z$, eigenvalues are recursively copied from lower $(N,E)$-blocks to higher $(N,E)$-blocks. Thus, each eigenvalue is present in infinitely many copies. For any two blocks $(N,E_{1})$ and $(N,E_{2})$ with $E_{1}<E_{2}$, the $\mathcal{H}_{LLL}$ eigenvalues of the first block are a subset of the eigenvalues of the second block. The energy shifts $\varepsilon_{I}$ in (23) for the two corresponding unperturbed levels are evidently in the same relation.

How do new eigenvalues emerge in this picture as we move to higher values of $E$? The action of $Z$ copies all the $p_{N}(E)$ eigenvalues of the $(N,E)$-block into the $(N,E+1)$-block. The latter block has $p_{N}(E+1)>p_{N}(E)$ eigenvalues, and the excess eigenvalues must correspond to eigenvectors annihilated by $Z^{\dagger}$. Indeed, if they were not annihilated by $Z^{\dagger}$, the action of $Z^{\dagger}$ would have produced an eigenvector the $(N,E)$-block with the same eigenvalue, in contradiction with our assumption that the eigenvalue is new. Thus, any new eigenvalues in the $(N,E+1)$-block compared to the $(N,E)$-block must come from vectors $|\Psi\rangle$ in the kernel of $Z^{\dagger}$,

[TABLE]

In order for the whole picture to be consistent, the dimension of this kernel within the $(N,E+1)$-block must precisely account for the difference in the dimensions of the $(N,E)$- and $(N,E+1)$-blocks:

[TABLE]

To see this structure more explicitly, note first that $Z$ cannot annihilate a nonvacuum state. Indeed, by (36),

[TABLE]

Thus, by acting with $Z$ on a complete basis $|\Psi_{I}\rangle$ in the $(N,E)$-block containing $p_{N}(E)$ elements, we obtain $p_{N}(E)$ linearly independent vectors $Z|\Psi_{I}\rangle$ in the $(N,E+1)$-block. Consider the orthogonal complement of the subspace spanned by these vectors within the $(N,E+1)$-block, which by construction has $p_{N}(E+1)-p_{N}(E)$ dimensions. Any of the vectors $|\Phi\rangle$ in this orthogonal complement must be annihilated by $Z^{\dagger}$ since, on the one hand, $Z^{\dagger}|\Phi\rangle$ belongs to the $(N,E)$-block, and on the other hand it is orthogonal to all states in the $(N,E)$-block, since $\langle\Psi_{I}|Z^{\dagger}|\Phi\rangle=0$ as a consequence of $|\Phi\rangle$ being orthogonal to $Z|\Psi_{I}\rangle$. Therefore any such $|\Phi\rangle$ is in the kernel of $Z^{\dagger}$, while none of the vectors $Z|\Psi_{I}\rangle$ can be in the kernel of $Z^{\dagger}$ as $|Z^{\dagger}Z|\Psi_{I}\rangle|^{2}>0$ by (36). We have thereby established (41).

The picture we have presented leads to simplifications in generating the spectrum of $\mathcal{H}_{LLL}$ numerically. Instead of individually diagonalizing all the $(N,E)$-blocks, as one does for generic resonant systems quantres , we may simply choose one block with a sufficiently large $E$, and diagonalize explicitly only this block. All the eigenvalues of the blocks with the same $N$ and lower $E$ are contained in this block due to the inheritance property we have described. These eigenvalues of the lower blocks can be extracted by either repeatedly acting on the eigenvectors of the given $(N,E)$-block with $Z^{\dagger}$, or by diagonalizing $ZZ^{\dagger}$ simultaneously with $\mathcal{H}_{LLL}$ in the given $(N,E)$-block. Indeed, the eigenvectors in the kernel (40) will be annihilated by $ZZ^{\dagger}$, while nonzero eigenvalues of $ZZ^{\dagger}$ will mark the energy eigenstates that have descended from the lower values of $E$ via the action of $Z$. The eigenvalues of $ZZ^{\dagger}$ can be obtained from the commutation relations (36), and they depend on the value of $E$ in the block where a given energy eigenvalue first appears before being transported to the current block by the action of $Z$. This is why inspecting the eigenvalues of $ZZ^{\dagger}$ allows for a recovery of the lower $(N,E)$-blocks from a given block.

One particular eigenvalue stands out prominently in the inheritance process we have described. One can start with the $(N,0)$-block that contains only one state, $|N,0,0,\ldots\rangle$, which is of course an eigenstate of $\mathcal{H}_{LLL}$,

[TABLE]

The corresponding eigenvalue $N(N-1)/2$ will be propagated by the action of $Z$ to all the higher $(N,E)$-blocks, where it will correspond to the vector $Z^{E}|N,0,0,\ldots\rangle$. We shall now derive bounds on $\mathcal{H}_{LLL}$, which make it clear that this eigenvalue is always the highest one in the spectrum, while the other eigenvalues of an $(N,E)$-blocks lie between [math] and $N(N-1)/2$.

III.3 Bounds on the LLL Hamiltonian

In establishing bounds on $\mathcal{H}_{LLL}$, the following explicit representation of the resonant summation in (33) will prove convenient:

[TABLE]

(Resonant summations in this form have been extensively employed in meloturb .) From this representation and (34), one may immediately write

[TABLE]

where

[TABLE]

Thus, $\mathcal{H}_{LLL}$ is manifestly a nonnegative operator and all of its eigenvalues are nonnegative as claimed under (43).

We then proceed to establish the upper bound on $\mathcal{H}_{LLL}$, which is slightly less straightforward. Consider

[TABLE]

Similarly to the previous proof, we wish to write $\mathcal{H}_{inv}$ as a sum of manifestly nonnegative terms. Using $[\alpha^{\dagger}_{l},\alpha_{k}]=-\delta_{kl}$, we first note that

[TABLE]

We then insert into this representation of $N(N-1)/2$ a tautological decomposition of unity in terms of the binomial identity

[TABLE]

With these preliminaries, one can straightforwardly write

[TABLE]

where

[TABLE]

This proves that $\mathcal{H}_{LLL}$ is bounded from above by $N(N-1)/2$.

To summarize, all the eigenvalues $\mathcal{H}_{LLL}$ within an $(N,E)$ block lie between 0 and the largest eigenvalue $N(N-1)/2$, which corresponds to the state $Z^{E}|N,0,0,\ldots\rangle$ as we have previously claimed.

III.4 Level spacing statistics

We have seen in the above treatment that the known symmetries of the LLL Hamiltonian have powerful implications for the structure of its spectrum. A question naturally arises whether the symmetries generated by (35) are all there is, or there are extra conserved operators in addition to $N$, $E$ and $Z$. While there is no algorithmic way to construct extra conserved operators or prove their absence, quantum chaos theory GMW ; haake provides an attractive set of tools to shed light on this question. The bulk of quantum chaos theory revolves around a set of conjectures derived from extensive numerical experimentation, and we will thus not have water-tight theorems at our disposal. Nonetheless, we feel that the indicators supplied by this approach are of great practical use in our context.

A key tenet of quantum chaos theory is that the quantum spectra of classically chaotic systems have distances between neighboring energy levels that obey a qualitatively different statistics from the quantum spectra of classically integrable systems. More specifically, distances between energy eigenvalues for a generic integrable system follow btint the Poisson distribution

[TABLE]

Here, $s$ are distances between neighboring levels of a properly normalized (more technically, ‘unfolded’) energy level sequence that we shall discuss below. The distribution (52) is the same as the distribution of distances between points randomly thrown on a line with a unit mean density. For a chaotic system, on the other hand, one expects BGS the Wigner-Dyson distribution, which is well approximated for practical purposes by the ‘Wigner surmise’

[TABLE]

This is the statistics of distances between eigenvalues of real symmetric random matrices with independent identically distributed entries. Note that (53) vanishes at zero energy level separation, a phenomenon known as ‘level repulsion.’

We have to clarify what is meant by the ‘unfolding’ procedure used to generate the normalized level distances $s$ from a given set of energy levels. The need for unfolding is most easily understood in the original context of the quantum chaos theory, that is, the energy spectra of chaotic quantum-mechanical systems with a few degrees of freedom, such as quantum billiards BGS . In this situation, there are infinitely many energy levels and the distances between the neighboring levels may become progressively longer or shorter as one moves to higher energies. Thus, one is not even guaranteed to have the distribution of plain energy distances converge to a limit as more and more energy levels are included. What one must do is differentially rescale the spectrum so that on intervals containing many energy levels the mean density is 1 in any part of the spectrum. Then, the common distribution function for level distances will exist, and it is the level distance distributions in thus unfolded level sequences for which the quantum chaos conjectures are formulated. A contemporary discussion of unfolding can be found in GMW ; abdulmagd ; luukko .

Unfolding may seem more subtle for inherently finite samples of energy levels, such as the levels within the $(N,E)$-blocks of the LLL Hamiltonian. Nonetheless, the following simple-minded definition quantres will prove perfectly effective for our purposes. Consider a sample with $K$ eigenvalues denoted $E_{I}$ with $I=1,\ldots,K$. We shall assume $K$ to be reasonably large (comparable to 1000 in our numerics). One can choose an integer $\Delta$ close to $\sqrt{K}$, and differentially stretch the spectrum so that the level density is 1 on intervals containing $2\Delta$ energy levels in any part of the spectrum. To this end, define the raw unfolded sequence of level spacings

[TABLE]

with $I=\Delta+1,\Delta+2,\ldots,K-\Delta$. The division by $(E_{I+\Delta}-E_{I-\Delta})$ factors out the mean level density over intervals containing $2\Delta$ adjacent levels and centered on the point of observation. We then compute the average of this raw sequence $\bar{s}=(\sum_{I=\Delta+1}^{K-\Delta}s_{I}^{(raw)})/(K-2\Delta)$, and define the normalized unfolded sequence

[TABLE]

whose mean by construction equals 1. This is our definition of $s$ for comparing with (52) and (53). While the smoothing parameter $\Delta$ and the bin size for plotting distribution histograms have to be chosen by hand (though $\sqrt{K}$ is the natural scale for both), there are no free parameters in the distributions (52) and (53) and thus no fitting involved in the comparison after the histograms have been plotted.

We note that the chaotic Wigner-Dyson distribution, which we have approximately represented by the Wigner surmise (53), is only expected to emerge if the spectrum is separated into blocks according to the values of all pairwise commuting conserved operators (the analogs of globally defined classical integrals of motion in involution). The underlying philosophy is that for each specific value of such commuting conserved operators, one essentially has an independent Hamiltonian system, and the spectra of these individual Hamiltonian systems are not expected to be correlated. Then, lumping together the different blocks for computing the level spacing statistics would wash out the level repulsion phenomenon, which is responsible for the shape of (53) and requires correlations between the energy levels (see for example section 3.1 of DKPR ).

This is precisely what one could see in the numerical analysis of the LLL Hamiltonian in quantres . In that analysis, guided by the general structure of the eigenvalue problems for resonant Hamiltonians, the eigenvalues were split into $(N,E)$-blocks. Within such blocks, a Poissonian distribution of unfolded level spacings was seen, and the question was raised about its implications. We now answer this question by noticing that $N$ and $E$ are not the only quantities commuting with each other and $\mathcal{H}_{LLL}$, as $ZZ^{\dagger}$ commutes with all the three. One must then separate the LLL Hamiltonian eigenvalues into blocks not only according to their values of $N$ and $E$, but also according to the values of $ZZ^{\dagger}$.

In practice, because of the structure of the eigenvalue problem of $\mathcal{H}_{LLL}$ described in section III.2, it is sufficient to consider the eigenvectors of $\mathcal{H}_{LLL}$ annihilated by $ZZ^{\dagger}$, i.e., the eigenvectors in the kernel of $Z^{\dagger}$. Indeed, all other $ZZ^{\dagger}$-blocks are simply inherited through repeated action of $Z$ on the $Z^{\dagger}$-kernel of some other $(N,E)$-block with a lower value of $E$. So the variety of statistics present in the problem is fully exhausted by looking at $Z^{\dagger}$-kernels only.

We have implemented this procedure for the $N=E=27$ block that has already been considered in quantres . If the joint level spacing statistics in the entire block is plotted, one sees a Poissonian distribution, as in quantres . If the levels outside the $Z^{\dagger}$-kernel are excised, one obtains the distribution in Fig. 2. This distribution is very far from the Poissonian shape, and close to the Wigner-Dyson shape, which is a strong argument against integrability. (For an integrable system, one expects the distribution to be Poissonian in all large blocks, along the lines of the Berry-Tabor results btint .) One may notice that the actual distribution in Fig. 2 is slightly shifted to the left compared to the Wigner-Dyson curve. This shift, in fact, increases when one moves to other blocks away from $N=M$. There, the level repulsion decreases, and one gets what is known as crossover behaviors. For $N\ll M$, the distribution acquires Poissonian features. This does not change our conclusions on integrability, but creates room for some subtle structures. The actual distributions are given in Appendix B.

III.5 Generalizations

We briefly comment on possible generalizations of the material of this section that take us outside the immediate scope of the present article. Resonant systems of the form (33), with various assignments of the interaction coefficients $C$, have appeared in studies of weakly nonlinear problems emerging from a range of physical applications. Apart from cold atomic gases and Bose-Einstein condensates, there are applications in high-energy and gravitational physics, in particular due to connections to the topic of instability of the Anti-de Sitter spacetime BR ; rev2 , and dynamics of related nonlinear wave equations in highly symmetric spacetimes CF ; BEL . Such resonant systems are also of relevance to the Hartree equation in harmonic potentials BEF .

A number of resonant systems emerging from such weakly nonlinear considerations CF ; BEL ; BEF ; BBCE2 possess features closely reminiscent of what we have used above to analyze the resonant LLL Hamiltonian. In AO , these features have been distilled into a formulation of a very large class of resonant systems that, in particular, admit an analog of the $Z$-conservation that has played a crucial role in the derivations of this section. In quantres , numerical analysis of the quantum version of these systems was performed, demonstrating a number of distinctive traits, such as the eigenvalue inheritance from lower to higher values of $E$ in the $(N,E)$-blocks, and the maximal eigenvalue within each block being equal to $N(N-1)/2$. We have seen in this section how these observations can be proved for the LLL system. Analogous proofs will hold for other systems in the large class formulated in AO . In particular, the $Z$-conservation will directly translate into eigenvalue inheritance, while an upper bound on the eigenvalues will follow (at least for systems with positive interaction coefficients) from a generalization of the summation identity (49) that plays a defining role in the constructions of AO . This opens up space for applications of the material of this section, in a slightly modified form, to a number of very different physical problems.

IV Generic energy levels

As we turn to generic $(N,E,M)$-blocks of the resonant Hamiltonian (24), we essentially have to retrace the steps taken in the previous section for the LLL blocks with $E=M$, but the amount of technical effort required increases. A large part of the reason is that the interaction coefficients (19) no longer have a simple closed form like (34) and are expressed through the Laguerre polynomials as

[TABLE]

There is also a much bigger set of conserved operators (27-32) that must be taken into account. (Note that integrals of quadruple products of Laguerre polynomials of the form (56) have many special properties. For example, for the radially symmetric case $m_{i}=0$, there is appreciable mathematical literature on positivity lagpos and combinatorial interpretations lagcomb1 ; lagcomb2 of such integrals.)

As for the LLL case, there is inheritance of eigenvalues of $\mathcal{H}_{res}$ from blocks with lower $E$ and $M$ to higher blocks via the action of $Z_{\pm}$ and $W$ as depicted on Fig. 1. New eigenvalues, not found in the lower blocks, always enter through the joint kernel of $Z_{\pm}^{\dagger}$ and $W^{\dagger}$, i.e., the subspace of states $|\Psi\rangle$ satisfying

[TABLE]

Similarly to the previous section, acting on the state where all the $N$ particles occupy the lowest energy mode with different arrangements of $Z_{\pm}$ and $W$ results in states of energy $N(N-1)/2$ in each of the higher $(N,E,M)$-blocks. Unlike the LLL sector, these states may have high multiplicities within their $(N,E,M)$-blocks as there are many inequivalent ways to reach the same $(N,E,M)$-block acting with different arrangements of $Z_{\pm}$ and $W$. As for the LLL sector, $N(N-1)/2$ is the highest energy state within each $(N,E,M)$-block as a consequence of bounds on the resonant Hamiltonian (24).

Establishing bounds on $\mathcal{H}_{res}$ is considerably more involved than for the case of $\mathcal{H}_{LLL}$, and we give the proofs in Appendix C. These proofs retrace the mathematical literature FGH ; continuous , but in a language more accessible to physicists. The result is that all the eigenvalues in general $(N,E,M)$-blocks lie between [math] and the maximal eigenvalue $N(N-1)/2$, as in the previous section.

The conserved operators $N$, $E$, $M$, $Z_{\pm}$, $Z_{\pm}^{\dagger}$, $W$ and $W^{\dagger}$ obey the following commutator algebra, where we only give the nonzero commutators explicitly:

[TABLE]

A convenient combination is

[TABLE]

satisfying $[Z_{\pm},\tilde{W}^{\dagger}]=0$. The operators $Z_{+}Z_{+}^{\dagger}$, $Z_{-}Z_{-}^{\dagger}$ and $\tilde{W}\tilde{W}^{\dagger}$ commute among themselves and with $N$, $E$, $M$ and $\mathcal{H}_{res}$. In fact, if one chooses a large $(N,E,M)$-block and diagonalizes $\mathcal{H}_{res}$, $Z_{+}Z_{+}^{\dagger}$, $Z_{-}Z_{-}^{\dagger}$ and $\tilde{W}\tilde{W}^{\dagger}$ simultaneously within this block, this allows for recovering the energy eigenvalues of all the lower blocks from which one can reach the given block by the moves in Fig. 1. Indeed, the kernel (57), containing all the newly added energy eigenvalues not found in the lower blocks, evidently corresponds to zero eigenvalues of $Z_{+}Z_{+}^{\dagger}$, $Z_{-}Z_{-}^{\dagger}$ and $\tilde{W}\tilde{W}^{\dagger}$, while all the other ‘descendant’ energy eigenvalues that have arrived in the current block from the lower blocks are marked by positive eigenvalues of $Z_{+}Z_{+}^{\dagger}$, $Z_{-}Z_{-}^{\dagger}$ and $\tilde{W}\tilde{W}^{\dagger}$. These latter eigenvalues are easily obtained from the commutation relations (58), and encode the number of moves in Fig. 1 necessary to reach the current block from the block where a particular energy eigenvalue first appeared.

In implementing our analysis of the level spacing statistics, as in section III.4, we must likewise separate the eigenvalues of $\mathcal{H}_{res}$ according to all the six mutually commuting conserved quantities $N$, $E$, $M$, $Z_{+}Z_{+}^{\dagger}$, $Z_{-}Z_{-}^{\dagger}$ and $\tilde{W}\tilde{W}^{\dagger}$. In practice, it is sufficient to consider zero values of the last three quantities, i.e., the eigenstates in the joint kernel (57). Eigenvalues outside this kernel in a given $(N,E,M)$-block always descend through the action of $Z_{\pm}$ and $W$ from the kernel eigenvalues of a lower $(N,E,M)$-block. With this picture in mind, we have plotted the unfolded level spacing distribution in the $(Z_{\pm}^{\dagger},W^{\dagger})$-kernel of the $(N,E,M)=(7,22,16)$ block, given in Fig. 3. The distribution is Wignerian and suggests that no symmetry structures beyond (27-32) are to be expected. We give two extra level spacing distributions for different blocks in Appendix B. In this case, unlike the LLL sector, the distribution is very close to the Wigner-Dyson curve for all the blocks we have studied.

V Discussion

We have analyzed the problem of energy level splitting for quantum bosons in a two-dimensional isotropic harmonic trap under the influence of weak contact interactions. The problem is shown to be identical (at linear order in the interaction strength) to diagonalizing the quantized resonant system of the Gross-Pitaevskii equation. Known symmetries of this resonant system impart a highly constrained structure to the energy spectrum. The action of the symmetries copies the energy shifts from lower to higher unperturbed energy levels, so that infinitely many exactly integer energy differences characteristic of noninteracting bosons in a harmonic trap survive (this is a reflection of the exact pattern in the spectrum valid at finite coupling). New values of energy shifts (in addition to the ones having descended from lower unperturbed levels via the action of the symmetries) always enter through explicitly described kernels of the symmetry generators. These symmetry structures can be used to optimize the numerical construction of eigenstate energies. The highest energy level splitting off the unperturbed level with $E$ units of energy and $N$ particles is explicitly given as

[TABLE]

while all the others lie between $E$ and $\tilde{E}_{max}$.

We have furthermore analyzed the level spacing statistics in the fine structure emanating from individual unperturbed levels under the influence of weak contact interactions. The purpose is to use the standard lore of quantum chaos theory to shed light on the question of possible integrability (or, more generally, extra symmetry structures) in the resonant system of the Gross-Pitaevskii equation (whose quantized version describes the energy level splitting). Speculations about possible integrability of this resonant system (and its close relatives) have been voiced in the literature FGH ; CF ; quantres , in particular due to some intriguing parallels seen between these systems and the cubic Szegő equation GG , which is known to be integrable. Our analysis of the level spacing statistics suggests that no integrability is to be expected. Outside the LLL sector (blocks with the maximal amount of rotation), we observe Wignerian level spacing distributions, characteristic of strongly chaotic systems. Within the LLL sector, there is a range of behaviors ranging from nearly Wignerian, to crossover, to Poissonian, which may be alluding to some very subtle analytic structure, but not to integrability.

A few classical systems possess symmetries closely analogous to those of the Gross-Pitaevskii equation, which is the classical version of (1), both at the level of the original PDE, and at the level of extra symmetry enhancements in the resonant approximation. We specifically mention the conformally coupled cubic wave equation on the 3-sphere CF and the Hartree equation in a harmonic trap in four spatial dimensions BEF . Similarly, if one replaces the $|\Psi|^{4}$ interaction that has characterized our considerations by $|\Psi|^{6}$, there are systems in lower numbers of spatial dimensions fennell ; quintic that respect symmetry structures of the sort underlying our current study. There is additionally a very large class of resonant systems explicitly constructed in AO in a way that guarantees the presence of the same type of symmetries. For all the systems mentioned, one expects that quantization will generate energy level patterns closely analogous to what we have observed in the present treatment.

We finally comment on the ‘quantum Gross-Pitaevskii equations’ QGP , a hierarchy of equations inspired by variational tensor network techniques and interpolating between the classical Gross-Pitaevskii equation and quantum bosons with contact interactions described by (1). Since the operation of enhanced symmetries in the weak coupling regime has by now been seen both in the Gross-Pitaevskii equation continuous ; GT ; BBCE ; BBCE2 ; GGT and in the full quantum system in our present work, one expects it to have manifestations within the entire quantum Gross-Pitaevskii hierarchy. It would be interesting to investigate this subject further.

Acknowledgements.

We thank Piotr Bizoń for a discussion on bounds for resonant Hamiltonians. This research has been supported by FWO-Vlaanderen (projects G044016N and G006918N), by Vrije Universiteit Brussel through the Strategic Research Program “High-Energy Physics,” and by CUniverse research promotion project (CUAASC) at Chulalongkorn University. M. D. C. is supported by a PhD fellowship from the Research Foundation Flanders (FWO).

Appendix A Conservation of $Z_{+}$, $Z_{-}$ and $W$

The three proofs in this appendix closely retrace section 3 and the appendix of AO . One has to prove that $[Z_{\pm},\mathcal{H}_{res}]=[W,\mathcal{H}_{res}]=0$. Evaluating the commutators directly, which is identical to the time-differentiation in section 3 of AO , produces expressions that vanish provided that certain four-term identities are satisfied by the interaction coefficients (19). These identities are then proved using the following properties of the Laguerre polynomials:

[TABLE]

A.1 Conservation of $Z_{+}$

Direct computation of the commutator $[Z_{+},\mathcal{H}_{res}]$ produces an expression that vanishes if

[TABLE]

whenever $n_{3}=n+n_{1}-n_{2}-1$ and $m_{3}=m+m_{1}-m_{2}-1$. Multiplying by

[TABLE]

and substituting the interaction coefficients (56), we write

[TABLE]

where we have used $n_{3}=n+n_{1}-n_{2}-1$ and $m_{3}=m+m_{1}-m_{2}-1$. We apply the identities (62-63) to the first two terms, while for the last two terms, we use (62), obtaining

[TABLE]

Using the remaining identity (61), we note that this expression is a total derivative:

[TABLE]

The Laguerre polynomials of the form $L_{\frac{n-m}{2}}^{m}$ do not have any powers of $\rho$ below $\rho^{-m}$, which ensures this last expression is zero, completing the proof.

A.2 Conservation of $Z_{-}$

The proof is analogous to the one above. Direct computation of the commutator $[Z_{-},\mathcal{H}_{res}]$ produces an expression that vanishes if

[TABLE]

whenever $n_{3}=n+n_{1}-n_{2}-1$ and $m_{3}=m+m_{1}-m_{2}+1$. Multiplying by (65) and substituting (56), we get

[TABLE]

We apply identity (61) on the first two terms. The last two terms can be rewritten by applying (62) first, followed by (63) and using again (62), which is equivalent to

[TABLE]

Afterwards, we group all the terms to isolate a total derivative, as in the previous proof

[TABLE]

The expression is zero, which proves the conservation of $Z_{-}$.

A.3 Conservation of $W$

Direct computation of the commutator $[W,\mathcal{H}_{res}]$ produces an expression that vanishes if

[TABLE]

whenever $n+n_{1}=n_{2}+n_{3}-2$ and $m+m_{1}=m_{2}+m_{3}$. Multiplying by (65) and substituting (56), we get

[TABLE]

We now apply the identity (63) to the four terms, yielding

[TABLE]

where we have used $n+n_{1}=n_{2}+n_{3}-2$ and $m+m_{1}=m_{2}+m_{3}$. Finally, we apply (62) to the second and third term, and identify a total derivative using (61), giving

[TABLE]

The Laguerre polynomials of the form $L_{\frac{n-m}{2}}^{m}$ do not have any powers of $\rho$ below $\rho^{-m}$, which ensures that the last expression is zero and $W$ is conserved.

Appendix B More information on the level spacing statistics

In this appendix, we present a few extra plots of level spacing statistics that are not crucial for our conclusions in the main text (that integrability is unlikely), but give a broader perspective, and also bring forth some interesting details.

As we remarked in the main text in section III, the LLL sector is characterized by distributions that are not exactly Wignerian, but somewhat shifted towards the Poissonian curve. Such distributions are often described as ‘crossover’ behaviors, and characterize chaotic systems with relatively weak chaotic dynamics (for example, small chaotic perturbations of integrable systems). Observing these distributions makes integrability unlikely, but leaves some room for more subtle features. The distributions appear close to Wignerian when $N\approx M$ and morph in the direction of the Poissonian curve as one moves away from that region. If $N\ll M$, Poissonian behaviors are observed. This is depicted in Fig. 4 for two particular blocks.

Away from the LLL sector, we have consistently observed Wignerian distributions, exemplified by two specific blocks in Fig. 5. In particular, the level repulsion is clearly visible as the leftmost point of each plot fits the Wignerian curve very accurately.

Appendix C Bounds on the two-dimensional resonant Hamiltonian

Proofs of bounds on the resonant Hamiltonian (24) are given in FGH , but in a language that may be not immediately accessible for readers with physics backgrounds. The purpose of this appendix is to recast these proofs in a way that is more straightforward, even if less rigorous. We also need to translate the proofs of FGH , given for the classical theory, into statements about the quantum resonant Hamiltonian (24).

C.1 Positivity

Consider

[TABLE]

with the normalized harmonic oscillator energy eigenstates $\psi_{nm}$ given by (17) and ${\bf r}\equiv(x,y)$. This expression is, of course, simply the Heisenberg operator of the field $\Psi$. We use the angle-like notation $\theta$ for the time variable in this representation to emphasize that the harmonic oscillator evolution (and hence the above Heisenberg operator) is exactly periodic. With this notation, the resonant Hamiltonian (24) with the interaction coefficients (19) is written as

[TABLE]

where we have used $\int_{-\pi}^{\pi}e^{i(n_{1}+n_{2}-n_{3}-n_{4})\theta}d\theta=2\pi\delta_{n_{1}+n_{2},n_{3}+n_{4}}$. Since this expression is an integral of the manifestly positive operator $\Psi^{\dagger 2}\Psi^{2}$, $\mathcal{H}_{res}$ is positive.

C.2 Upper bound

Proving the upper bound $\mathcal{H}_{res}\leq N(N-1)/2$ will require more effort than the positivity proof. The classical version of this statement is proved in FGH and turns out equivalent to constructing bounds on the so-called Strichartz norms. The latter problem has been extensively studied in the mathematical literature starting with strnorm1 ; strnorm2 , however, sharp bounds with explicit numerical coefficients of the sort we need here appeared only rather recently strnorm3 ; strnorm4 . We find the proofs of strnorm4 particularly graceful and easy to adapt for our purposes. Arguments based on the heat kernel, as in FGH ; strnorm3 , are not straightforwardly translated to the quantum Hamiltonian (24). Our presentation below essentially retraces the proof of strnorm4 in a physicist’s language, and extends it to the quantum case.

The evolution of (76) is expressed through the harmonic oscillator propagator as

[TABLE]

We note that $\Psi({\bf r},\theta+\pi)=-\Psi(-{\bf r},\theta)$, which lets us reduce the range of integration in (77) as

[TABLE]

Then,

[TABLE]

where we have reverted to the notation $\Psi({\bf r})\equiv\Psi({\bf r},0)$. Now introduce ${\bf k}={\bf r}/\sin\theta$ and $\tau=\cot\theta$ so that

[TABLE]

We note that the transformations we have just performed are closely linked to the connection between harmonic oscillator and free motion known as the lens transform Carles ; Tao . One may notice the appearance of the characteristic energy-momentum conservation of a free particle in the $\delta$-functions of (81). It may seem surprising that the ‘energies’ and ‘momenta’ in the $\delta$-functions are expressed through the coordinates ${\bf r}_{i}$, but in fact the expression is invariant under Fourier transform, as emphasized in FGH , so the momentum-space analog of (81) is given by the same expression, but with ${\bf r}_{i}$ replaced by the corresponding momenta ${\bf k}_{i}$ and the conservation of energy given by the standard expression. Note that the said invariance of $\mathcal{H}_{res}$ under the Fourier transform of $\Psi$ is a direct consequence of (79) being an integral of a periodic function, whose integration range can be freely shifted, and (78) being invariant under a shift of $\theta$ by $\pi/2$ combined with a Fourier transform (which is the well-known emergence of Fourier transforms at a quarter-period in the evolution of the quantum harmonic oscillator).

Back to our problem, introduce ${\bf R}=({\bf r}_{1}+{\bf r}_{2})/2$, ${\bf s}=({\bf r}_{1}-{\bf r}_{2})/2$, ${\bf R}^{\prime}=({\bf r}_{3}+{\bf r}_{4})/2$, ${\bf s}^{\prime}=({\bf r}_{3}-{\bf r}_{4})/2$ so that $d{\bf r}_{1}d{\bf r}_{2}d{\bf r}_{3}d{\bf r}_{4}=4d{\bf R}\,d{\bf s}\,d{\bf R}^{\prime}\,d{\bf s}^{\prime}$ and

[TABLE]

We now examine the action of the following operator $\hat{S}$ on a general (complex number valued) test function $h({\bf R},{\bf s})$:

[TABLE]

In polar coordinates $(s,\sigma)$ on the ${\bf s}$-plane, any function can be expanded in a Fourier series as

[TABLE]

Then,

[TABLE]

Thus, $\hat{S}$ is simply the orthogonal projector on the zero angular momentum subspace ($m=0$). But any projector is bounded from above by the identity operator, and hence, for any $h({\bf R},{\bf s})$,

[TABLE]

Consider the expectation value $\langle\Phi|\mathcal{H}_{res}|\Phi\rangle$ with an arbitrary state $|\Phi\rangle$. Insert the decomposition of unity in terms of a complete basis of states $|\Phi^{\prime}\rangle$ into $\mathcal{H}_{res}$ given by (82) between $\Psi^{\dagger}({\bf R}+{\bf s})\,\Psi^{\dagger}({\bf R}-{\bf s})$ and $\Psi({\bf R}^{\prime}+{\bf s}^{\prime})\,\Psi({\bf R}^{\prime}-{\bf s}^{\prime})$ to obtain

[TABLE]

with $h_{\Phi\Phi^{\prime}}({\bf R},{\bf s})=\langle\Phi^{\prime}|\Psi({\bf R}+{\bf s})\Psi({\bf R}-{\bf s})|\Phi\rangle$. By (86), for any $|\Phi\rangle$,

[TABLE]

Converting back from ${\bf R}$ and ${\bf s}$ to ${\bf r}_{1}$ and ${\bf r}_{2}$, we get

[TABLE]

since $N=\int\Psi^{\dagger}({\bf r})\Psi({\bf r})d{\bf r}$ and $[\Psi^{\dagger}({\bf r}^{\prime}),\Psi({\bf r})]=-\delta({\bf r}^{\prime}-{\bf r})$. We have thus proved that $\mathcal{H}_{res}$ is bounded from above by $N(N-1)/2$, as intended.

We conclude with some remarks on the proof of the upper bound we have just presented. Our proof is admittedly rather naive from a mathematical perspective. First, we have not been careful in specifying the functional spaces to which the various functions belong in deriving (86), though a mathematically rigorous derivation of the same bound can be found in strnorm4 . Second, in the quantum context, everything becomes even more subtle, as rigorous analysis would require a much more accurate definition of the space of states and dealing with issues of convergence in the decomposition of unity in (87), etc. We believe that our derivations capture the essence of the problem, however. Importantly, the end result is just a constraint on the diagonalization of finite-sized matrices for the problem of level splitting found in the main text, which is perfectly well-defined in terms of elementary mathematics. This constraint (the eigenvalues lying between [math] and $N(N-1)/2$) is furthermore empirically validated by our numerical diagonalization for concrete energy levels.

The derivations presented above are phrased through the representation of the resonant Hamiltonian given by (81), which is illuminating, but rather unintuitive in the context of a harmonic trap problem. It would be interesting to construct a more direct proof in the Laguerre polynomial basis, but we are not aware of an easy way to do this. If one only keeps the modes with $m=0$ in the resonant system, the positivity results for integrals of products of Laguerre polynomials lagpos and the summation identities of BBCE2 allow for a direct generalization of the elementary upper bound proof in section III.3. This proof cannot be immediately adapted, however, to the full resonant system, which has both positive and negative interaction coefficients.