arXiv:yymm.nnnn

A superintegrable quantum field theory Marine De Clerck DAMTP, University of Cambridge, Cambridge, United Kingdom Oleg Evnin High Energy Physics Research Unit, Faculty of Science,
Chulalongkorn University, Bangkok, Thailand

Theoretische Natuurkunde, Vrije Universiteit Brussel (VUB)
&   International Solvay Institutes, Brussels, Belgium [email protected], [email protected] ABSTRACT

Gérard and Grellier proposed, under the name of the cubic Szegő equation, a remarkable classical field theory on a circle with a quartic Hamiltonian. The Lax integrability structure that emerges from their definition is so constraining that it allows for writing down an explicit general solution for prescribed initial data, and at the same time, the dynamics is highly nontrivial and involves turbulent energy transfer to arbitrarily short wavelengths. The quantum version of the same Hamiltonian is even more striking: not only the Hamiltonian itself, but also its associated conserved hierarchies display purely integer spectra, indicating a structure beyond ordinary quantum integrability. Here, we initiate a systematic study of this quantum system by presenting a mixture of analytic results and empirical observations on the structure of its eigenvalues and eigenvectors, conservation laws, ladder operators, etc.

1 Introduction

Our main concern here will be with the quantum Hamiltonian

$$H\equiv\frac{1}{2}\sum_{n+m=k+l}^{n,m,k,l\geq 0}a^{\dagger}_{n}a^{\dagger}_{m}a_{k}a_{l},$$ (1.1)

where $a^{\dagger}_{n}$ and $a_{n}$ with $n\geq 0$ are a set of usual bosonic creation-annihilation operators satisfying $[a_{n},a_{m}^{\dagger}]=\delta_{mn}$. The structure is that of a typical particle-number-conserving (nonrelativistic) quantum field theory with quartic interactions, as one may commonly see in theoretical condensed matter physics, yet the couplings between the modes are chosen in such a way as to impart the system a rich and elaborate integrability structure. Our goal is to investigate this structure.

The classical version of (1.1) is the Hamiltonian

$$H_{cl}\equiv\frac{1}{2}\sum_{n+m=k+l}^{n,m,k,l\geq 0}\bar{\alpha}_{n}\bar{\alpha}_{m}\alpha_{k}\alpha_{l}$$ (1.2)

for complex canonical coordinates $\alpha_{n}(t)$ whose conjugate momenta are $i\bar{\alpha}_{n}(t)$, with bars denoting the ordinary complex conjugation. The Hamiltonian (1.2) was introduced, studied and solved in a series of works by Patrick Gérard and Sandrine Grellier [1, 2, 3, 4], using position space notation that we shall explain below (as opposed to our mode space notation). We shall be referring to (1.1) and (1.2) as Gérard-Grellier (GG) Hamiltonians.

The main motivation for introducing and studying (1.2) was the topic of turbulent energy transfer from modes with long wavelength (small $n$) to those with short wavelength (large $n$). Such turbulent phenomena are an important topic in contemporary PDE mathematics, and are typically challenging to analyze – see [5] for a well-known example where a “weak weak” form of turbulence is established rigorously with considerable effort. The GG Hamiltonian (1.2) plays a very special role in this regard since it displays an array of turbulent phenomena that can be explicitly analyzed using integrability methods. We shall review the classical integrability of (1.2) in section 3 after spelling out in section 2 its place in the broader context of resonant Hamiltonian systems with their connections to weakly nonlinear equations appearing in physics.

While in-depth understanding is available for the integrable structure of the classical Hamiltonian (1.2) with its Lax pairs and hierarchies of conservation laws, this knowledge does not directly transfer to the quantum system, as naive quantization of the classical expressions without special attention to the ordering ambiguities does not recover the conserved charges of the quantum system. Perhaps the most striking feature of the quantum Hamiltonian (1.1), first pointed out in [6], is that its eigenvalue spectrum is made of integers. This feature is beyond any regular expectations, even for a quantum integrable system. Indeed, while a fundamental definition of quantum integrability is lacking [7, 8], the general expectation is that the energy levels of a typical quantum integrable system look like independent random numbers thrown on the real line [9], without level correlations or level repulsion. The perfect integers with rigid level spacings are clearly a different story, and such spectra are typically seen in superintegrable quantum-mechanical systems [10]. In view of this picture, we would like to refer to (1.1) as a superintegrable quantum field theory.

There is more to the story. In what follows, we shall present some approaches to constructing quantum conservation laws. One example is

$$H_{\min}=\hskip-5.69054pt\sum_{\begin{subarray}{c}n,m,k,l=1\\ n+m=k+l\end{subarray}}^{\infty}\hskip-5.69054pt\min(n,m,k,l)\ a^{\dagger}_{n}a^{\dagger}_{m}a_{k}a_{l}+\sum_{k=1}^{\infty}k^{2}a^{\dagger}_{k}a_{k},$$ (1.3)

first pointed out in [11]. This operator commutes with (1.1), and there will be further operators of this sort as well. All the conservation laws we find have purely integer spectra.

In what follows, we shall summarize our current understanding of the quantum Gérard-Grellier Hamiltonian, its hierarchy of conservation laws, and their spectra. This investigation will combine analytic and numerical tools.

2 The context of resonant Hamiltonian systems

Before focusing on the GG Hamiltonian (1.1), it could be good to step back for a moment and examine the more general system

$$H_{res}\equiv\frac{1}{2}\sum_{n+m=k+l}^{n,m,k,l\geq 0}C_{nmkl}\bar{\alpha}_{n}\bar{\alpha}_{m}\alpha_{k}\alpha_{l}$$ (2.1)

with arbitrary coupling coefficients $C$ satisfying $C_{nmkl}=C_{nmlk}=\bar{C}_{klnm}$. Evidently, setting $C$ to 1 leads to (1.2), while subsequent canonical quantization leads to (1.1). We refer to such Hamiltonians as resonant Hamiltonians because of the presence of the resonance condition $n+m=k+l$.

Hamiltonian systems of the form (2.1) emerge naturally from weakly nonlinear approximations to Hamiltonian PDEs whose linearized normal modes possess a perfectly resonant, equispaced spectrum of frequencies [12]. Examples include studies of nonlinear Schrödinger equation with harmonic potentials [13, 14, 15, 16, 17, 18, 19], which is a model for trapped Bose-Einstein condensates [20], as well as studies of waves in asymptotically anti-de Sitter spacetimes [21, 22, 23, 25, 26, 27], with and without gravitational interactions, in relation to the anti-de Sitter instability conjecture of Bizoń and Rostworowski [28]. These systems have also been studied in their own right, mostly as a dynamical arena for turbulent transfer of excitations from $\alpha_{n}$ with small $n$ to high $n$ (long wavelength to short wavelength in physical settings). This line of research includes the original articles [1, 2, 3, 4] that introduced and studied (1.2), as well as [29, 30, 31, 32, 33, 34]. Because of the presence of two conservation laws

$$N=\sum_{n=0}^{\infty}|\alpha_{n}|^{2},\qquad M=\sum_{n=1}^{\infty}n|\alpha_{n}|^{2},$$ (2.2)

transfer of excitations to higher $n$ is necessarily accompanied by transfer of excitations to low $n$, which is known as a ‘dual cascade’ [24], see [33, 34] for analytically tractable examples.

For a quick demonstration of how a Hamiltonian system of the form (2.1) can emerge from a realistic, physical PDE, consider the one-dimensional nonlinear Schrödinger equation with a harmonic potential:

$$i\,\frac{\partial\Psi}{\partial t}=\frac{1}{2}\left(-\frac{\partial^{2}}{\partial x^{2}}+x^{2}\right)\Psi+g|\Psi|^{2}\Psi.$$ (2.3)

With the nonlinearity turned off ($g=0$), the general solution is

$$\Psi=\sum_{n=0}^{\infty}\alpha_{n}\psi_{n}(x)e^{-iE_{n}t},\qquad E_{n}=n+\frac{1}{2},\qquad\frac{1}{2}\left(-\frac{\partial^{2}}{\partial x^{2}}+x^{2}\right)\psi_{n}=E_{n}\psi_{n},$$ (2.4)

with constants $\alpha_{n}$. At small nonzero coupling $g$, $\alpha_{n}$ are no longer constant in time and acquire slow variations. Assuming thus that $\alpha_{n}$ are functions of time, we can substitute (2.4) into (2.3) to obtain

$$i\,\frac{d\alpha_{n}}{dt}=g\sum_{k,l,m=0}^{\infty}C_{nmkl}\,\bar{\alpha}_{m}\alpha_{k}\alpha_{l}\,e^{i(E_{n}+E_{m}-E_{k}-E_{l})t},$$ (2.5)

where $C_{nmkl}=\int dx\,\psi_{n}\psi_{m}\psi_{k}\psi_{l}$. At small $g$, $\alpha_{n}$ vary slowly, while the last exponential factor oscillates on time scales of order 1. It is legitimate to expect that these ‘fast’ oscillations average out and only resonant terms with $E_{n}+E_{m}-E_{k}-E_{l}=n+m-k-l=0$ contribute substantially to the evolution. Technically, this provides a good approximation up to $t\sim 1/g$ (but not on longer timescales). The result of eliminating nonresonant terms is precisely the Hamiltonian equation of motion of (2.1), after $g$ has been absorbed into a redefinition of time. Related mathematically rigorous proofs protecting the accuracy of this resonant approximation can be found in [13, 17]. If one wants to obtain the GG Hamiltonian (1.2) via this resonant approximation process, one has to start with the so-called half-wave equation on a circle [35].

Only special values of the couplings $C$ arise from common physical PDEs, so an interesting question is how to engineer systems that produce desired values of couplings (including the ones we shall focus on in this paper). Ultracold atomic gases in harmonic traps provide a good setting, since harmonic potentials give rise to resonant structures necessary for (2.1). An outstanding problem is, however, to tune the couplings between the normal modes (for instance, by incorporating position-dependent nonlinearities). Another option is waves on a circle, but those waves will have to move in one direction, since the wavenumber index $n$ in $\alpha_{n}$ is nonnegative. Such ‘chiral waves’ are known to exist as edge states in some condensed matter system [36]. Once again, the interactions of individual normal modes will have to be differentially controlled.

Canonical quantization of (2.1) leads to the corresponding quantum Hamiltonian

$$H_{quant}=\frac{1}{2}\sum_{n+m=k+l}^{n,m,k,l\geq 0}C_{nmkl}a^{\dagger}_{n}a^{\dagger}_{m}a_{k}a_{l},$$ (2.6)

with $[a_{n},a_{m}^{\dagger}]=\delta_{mn}$. This Hamiltonian is strikingly simple to work with, at least at the level of numerical diagonalization, independent of the actual values of the couplings $C$. Indeed, quantization of (2.2) immediately leads to two quantum conservation laws

$$N=\sum_{n=0}^{\infty}a^{\dagger}_{n}a_{n},\qquad M=\sum_{n=1}^{\infty}n\,a^{\dagger}_{n}a_{n},$$ (2.7)

while the space of states is spanned by the Fock vectors: there is one such vector for each set of occupation numbers $\eta_{n}$, and it is obtained by acting on the vacuum $|0\rangle$ (annihilated by all $a_{n}$) with the creation operators $a^{\dagger}_{n}$ as $(a_{0}^{\dagger})^{\eta_{0}}(a_{1}^{\dagger})^{\eta_{1}}\cdots|0\rangle$. But there is only a finite number of choices of $\eta_{n}$ for each $N$ and $M$ (given by the number of integer partitions of $M$ into at most $N$ parts), while all matrix elements of $H_{quant}$ between Fock states with different $N$ or $M$ must vanish. For that reason, diagonalization of (2.6) is reduced to diagonalizing an infinite family of finite-sized matrices, which makes the eigenvalues directly accessible [6]. This is in fact how the integer eigenvalues of (1.1) were originally observed. Note that all of this has nothing to do with solvability of (2.1) and (2.6) and the structure is present for any values of the couplings, despite the fact that the corresponding classical dynamics may be arbitrarily complicated (chaotic, turbulent, etc). For special values of $C$ corresponding to quantum integrable systems, as in (1.1), there is much more structure still.

Just like its classical counterpart (2.1), the quantum Hamiltonian (2.6) arises as a weak coupling approximation to certain quantum field theories, and this is even more intuitive than the classical story [37, 38, 39, 40]. Indeed, when the classical normal modes have commensurate frequencies, the corresponding quanta have commensurate energies. As a result, quantization of free field theories with highly resonant spectra (nonrelativistic Schrödinger fields in harmonic traps [38, 39], relativistic fields in anti-de Sitter spacetimes [37, 40], etc) leads to highly degenerate multiparticle energy spectra. When weak interactions are turned on, these degenerate levels split, and $H_{quant}$ is the Hamiltonian one needs to diagonalize, in line with the standard Rayleigh–Schrödinger perturbation theory for degenerate spectra, to obtain the energy shifts that split the degenerate noninteracting energy levels. This picture has led to studies of the so-called lowest Landau level Hamiltonian, which belongs to the class (2.6) in relation to the physics of trapped ultracold atomic gases [41, 42, 43, 44, 45].

3 Classical integrability of the Gérard-Grellier Hamiltonian

Before proceeding to the quantum case (1.1), it would be wise to review the classical integrability theory of (1.2) developed in detail in the original works [1, 2, 3, 4]. In doing so, we will also attempt to adapt the presentation for the eyes of physicists.

The equations of motion corresponding to (1.2) read

$$i\dot{\alpha}_{n}=\sum_{m=0}^{\infty}\sum_{k=0}^{n+m}\bar{\alpha}_{m}\alpha_{k}\alpha_{n+m-k},$$ (3.1)

where dots will denote time derivatives. For some purposes, it is convenient to rewrite these equations in terms of the position space field $u(\theta)$ on the unit circle $\theta\in[0,2\pi)$:

$$u(\theta,t)\equiv\sum_{n=0}^{\infty}e^{in\theta}\alpha_{n}(t),$$ (3.2)

which yields

$$i\dot{u}=\Pi(|u|^{2}u),$$ (3.3)

with

$$\Pi\left(\sum_{n=-\infty}^{\infty}h_{n}e^{in\theta}\right)\equiv\sum_{n=0}^{\infty}h_{n}e^{in\theta}$$ (3.4)

being the Szegő projector. In (3.3), one can recognize the original cubic Szegő equation as introduced in [1]. Note that $u(\theta)$ naturally extends to the complex plane by writing $u(z)=\sum_{n=0}^{\infty}z^{n}\alpha_{n}(t)$, in which case $u(z=e^{i\theta})$ recovers the original definition $u(\theta)$. Many statements about the classical system are proved very effectively in this position space picture using properties of Szegő projectors, but in preparation for analyzing the quantum system, we will often stress the mode space representation, since, once quantized, the complex amplitudes $\alpha_{n}$ turn into a conventional creation-annihilation operator algebra, while the corresponding field operators $u(\theta)$ acquire rather unusual nonlocal commutation relations.

It was recognized early on [1] that (3.3) is a very special equation and in particular, it admits an infinite hierarchy of dynamically invariant manifolds formed by meromorphic $u(z)$ with a fixed number of poles in the complex plane. If the initial conditions are of this form, the evolution described by (3.3) will amount to moving the positions of the poles and changing their residues, but the number of singularities will not change. The poles thus act as solitons. In the mode space language, the invariant manifolds are defined by configurations of the form

$$\alpha_{n}(t)=\sum_{r=1}^{R}c_{r}(t)\,[p_{r}(t)]^{n}$$ (3.5)

at some prescribed value of $R$. Further inquiry into the structure of (3.3) has led to the discovery of its Lax pairs that we will briefly summarize immediately below. The original presentation of [1] revolves around the antilinear operator $H_{u}h=\Pi(u\bar{h})$, resulting in a Lax pair in which the first Lax operator is antilinear. This is used very wisely in the constructions of [1], but as it is rather unconventional, we will rely in our presentation on the linear Lax operator obtained by squaring $H_{u}$.

We proceed with presenting the Lax pair in the mode space language of (1.2). To this end, we introduce auxiliary ‘test’ vectors in the mode space $\vec{h}=(h_{0},h_{1},\cdots)$, and also, for future use, the special vector

$$\vec{1}=(1,0,0,\cdots).$$ (3.6)

The classical Lax pair $(\mathcal{L},\mathcal{M})$ is defined by its action on $\vec{h}$ as

$$(\mathcal{L}\vec{h})_{n}=\sum_{k,l=0}^{\infty}\alpha_{n+k}\bar{\alpha}_{k+l}h_{l},\qquad(\mathcal{M}\vec{h})_{n}=\sum_{m+p\geq n}^{m,p\geq 0}\alpha_{p}\bar{\alpha}_{m+p-n}h_{m}.$$ (3.7)

One can show that the relation

$$i\dot{\mathcal{L}}=[\mathcal{M},\mathcal{L}]$$ (3.8)

holds whenever the equations of motion (3.1) are satisfied. (There is an elegant proof of this statement in [1] using Szegő projectors, but for completeness, we provide, in Appendix A, a brute force proof in the mode space language used here.) This has all the usual consequences of the Lax theory, in particular, the conservation of an infinite tower of charges given by $\mathrm{Tr}[\mathcal{L}^{p}]$ for any $p$, whose time derivatives vanish by (3.8) and the cyclic property of the trace. But there is much more structure in this particular case than what is dictated by the general Lax theory.

One special property of the Lax operators (3.7) is

$$\mathcal{L}\vec{1}=\mathcal{M}\vec{1}.$$ (3.9)

As a result, one gets an extra Lax pair that provides an extension of the Lax structure (3.8). Namely, define the projector on $\vec{1}$:

$$P_{1}\vec{h}\equiv(\vec{1},\vec{h})\,\vec{1},$$ (3.10)

where the scalar products are defined by the evident formula

$$(\vec{h},\vec{g})=\sum_{n=0}^{\infty}\bar{h}_{n}g_{n}.$$ (3.11)

As $P_{1}$ is a time-independent operator, and in view of (3.9), we can write a trivial Lax equation

$$i\dot{P}_{1}=[\mathcal{L}-\mathcal{M},P_{1}],$$ (3.12)

where both sides are identically zero. At the same time, we can equivalently rewrite (3.8) as

$$i\dot{\mathcal{L}}=[\mathcal{L}-\mathcal{M},\mathcal{L}]$$ (3.13)

This shows that $\mathcal{L}$ and $P_{1}$ are compatible Lax operators (with the same Lax partner) and therefore, the trace of any product of powers of $\mathcal{L}$ and $P_{1}$ will be conserved by the usual Lax construction. In practice, this leads to the conservation of the following quantities:

$$I_{n}\equiv\mathrm{Tr}[\mathcal{L}^{n}P_{1}]=(\vec{1},\mathcal{L}^{n}\vec{1})=\hskip-8.53581pt\sum_{i_{1},\ldots,i_{2n-1}=0}^{\infty}\hskip-8.53581pt\bar{\alpha}_{i_{1}}\alpha_{i_{1}+i_{2}}\bar{\alpha}_{i_{2}+i_{3}}\alpha_{i_{3}+i_{4}}\cdots\bar{\alpha}_{i_{2n-2}+i_{2n-1}}\alpha_{i_{2n-1}}.$$ (3.14)

These coexist with the aforementioned standard Lax conservation laws

$$G_{n}\equiv\mathrm{Tr}[\mathcal{L}^{n}]=\hskip-8.53581pt\sum_{i,i_{1},\ldots,i_{2n-1}=0}^{\infty}\hskip-8.53581pt\bar{\alpha}_{i+i_{1}}\alpha_{i_{1}+i_{2}}\bar{\alpha}_{i_{2}+i_{3}}\alpha_{i_{3}+i_{4}}\cdots\bar{\alpha}_{i_{2n-2}+i_{2n-1}}\alpha_{i_{2n-1}+i},$$ (3.15)

which however find very little use in the analytic considerations of [1, 2, 3, 4]. The $I$-tower starts with $I_{1}=N$, while $I_{2}$ can be expressed through (1.2) and $N^{2}$. Similarly, the $G$-tower starts with $G_{1}=N+M$, while $G_{2}$ can be related to

$$\sum_{n+m=k+l}^{n,m,k,l\geq 0}[1+\min(n,m,k,l)]\,\bar{\alpha}_{n}\bar{\alpha}_{m}\alpha_{k}\alpha_{l}=2H_{cl}+H_{\min,cl}\,,$$

with

$$H_{\min,cl}\equiv\sum_{n+m=k+l}^{n,m,k,l\geq 0}\min(n,m,k,l)\,\bar{\alpha}_{n}\bar{\alpha}_{m}\alpha_{k}\alpha_{l}\,,$$ (3.16)

the classical analog of the quantum conserved operator $H_{\min}$ defined in (1.3).

There is another Lax-pair construction of a similar sort. Consider $P_{\alpha}$ defined by

$$P_{\alpha}\vec{h}=(\vec{\alpha},\vec{h})\,\vec{h},$$ (3.17)

where $\vec{\alpha}=(\alpha_{0},\alpha_{1},\cdots)$. By (3.1), we have

$$i\dot{P}_{\alpha}=[\mathcal{M},P_{\alpha}].$$ (3.18)

In other words, $P_{\alpha}$ is Lax-compatible with $\mathcal{L}$ (but not with $P_{1}$). This means we can freely combine $\mathcal{L}$ and $P_{\alpha}$, for example, by defining

$$\tilde{\mathcal{L}}=\mathcal{L}-P_{\alpha},$$ (3.19)

which would be called $K_{u}^{2}$ in the notation of [3]. Conservation laws may be built from $\tilde{\mathcal{L}}$, though they are not independent of the conservation laws built from $\mathcal{L}$. One will have the usual traces $\mathrm{Tr}[\tilde{\mathcal{L}}^{n}]$, as well as

$$\tilde{I}_{n}\equiv\mathrm{Tr}[\tilde{\mathcal{L}}^{n-1}P_{\alpha}]=(\vec{\alpha},\tilde{\mathcal{L}}^{n-1}\vec{\alpha})=\hskip-8.53581pt\sum_{i_{1},\ldots,i_{2n-1}=0}^{\infty}\hskip-8.53581pt\bar{\alpha}_{i_{1}}\alpha_{i_{1}+i_{2}+1}\bar{\alpha}_{i_{2}+i_{3}+1}\cdots\bar{\alpha}_{i_{2n-2}+i_{2n-1}+1}\alpha_{i_{2n-1}}.$$ (3.20)

To conclude this brief review of the classical integrability properties of the cubic Szegő equation, we mention that a general formula can be written for the evolution of arbitrary prescribed initial configurations. This is much more than what one expects in general from a Lax-integrable system. Translating the original results of [3] into the mode space language we are currently using, we introduce the initial data vector $\vec{\alpha}_{start}\equiv(\alpha_{0}(0),\alpha_{1}(0),\cdots)$ and the corresponding initial Lax operators $\mathcal{L}_{0}\equiv\mathcal{L}[\vec{\alpha}_{start}]$ and $\tilde{\mathcal{L}}_{0}\equiv\tilde{\mathcal{L}}[\vec{\alpha}_{start}]$, as well as the shift operator $S$ together with its conjugate $S^{\dagger}$:

$$S\vec{h}=(0,h_{0},h_{1},h_{2},\cdots),\qquad S^{\dagger}\vec{h}=(h_{1},h_{2},h_{3},\cdots).$$ (3.21)

Then,

$$\alpha_{n}(t)=([e^{-it\mathcal{L}_{0}}e^{it\tilde{\mathcal{L}}_{0}}S^{\dagger}]^{n}e^{-it\mathcal{L}_{0}}\vec{\alpha}_{start}|\vec{1}).$$ (3.22)

Further details can be found in [3, 4].

Lax-integrable deformations of the cubic Szegő equation have been constructed [29, 32], though only part of the analytic structure presented here survives when the equation has been deformed. Finally, we comment on the Lax pair corresponding to the classical evolution defined by (3.16), which we found in the course of our investigations. The operator $\mathcal{L}$ of (3.7) remains unchanged, while its Lax partner $\mathcal{M}$ must be replaced with

$$(\mathcal{M}\vec{h})_{n}=\sum_{m+p\geq n}^{m,p\geq 0}n\,\alpha_{p}\bar{\alpha}_{m+p-n}h_{m}-\hskip-8.5359pt\sum^{\min(n,m)\geq p}_{m,p\geq 0}\hskip-5.69046ptp\,\alpha_{n-p}\bar{\alpha}_{m-p}h_{m}-\hskip-2.84544pt\sum_{n+p\geq m}^{m\geq n}(m-n)\,\alpha_{p}\bar{\alpha}_{n+p-m}h_{m}$$ (3.23)

in order to accommodate the evolution generated by the Hamiltonian (3.16).

4 The quantum Gérard-Grellier Hamiltonian

We finally return to the quantum Hamiltonian (1.1), which will be the main subject of our study for the rest of this paper. As explained at the end of section 2, the space of states of this Hamiltonian is spanned by the Fock vectors

$$|\eta_{0},\eta_{1},\cdots\rangle\equiv\left(\prod_{k=0}^{\infty}\frac{(a^{\dagger}_{k})^{\eta_{k}}}{\sqrt{\eta_{k}!}}\right)|0,0,\cdots\rangle,\quad a_{k}|0,0,\cdots\rangle=0,\quad a_{k}^{\dagger}a_{k}|\eta_{0},\eta_{1},\cdots\rangle=\eta_{k}|\eta_{0},\eta_{1},\cdots\rangle.$$

Because the operators $N$ and $M$ defined in (2.7) commute with the Hamiltonian, the latter can only have nonvanishing matrix elements between the Fock states with the same values of $N$ (given by $\sum_{k}\eta_{k}$) and $M$ (given by $\sum_{k}k\eta_{k}$). But since there is only a finite number of such states within each $(N,M)$-block, one is left with independent finite-sized matrices. Diagonalizing these finite-sized matrices reveals integer eigenvalues, providing thereby an entry point [6] to the manifold analytic puzzles presented by (1.1). What kind of structure can we expect to underlie these patterns, also in view of the classical integrability properties reviewed in section 3?

The integer eigenvalues suggest that ladder operators must be present that convert eigenvectors to eigenvectors while changing the eigenvalues by integer shifts. We will indeed report some such operators at the end of this section, and more in section 9.

Quantization of the hierarchies of conservation laws given by (3.14) and (3.15) in the classical theory is a very natural outstanding question. Naive quantization based on replacing $\bar{\alpha}_{n}$ with $a_{n}^{\dagger}$ and $\alpha_{n}$ with $a_{n}$ fails for any obvious ordering prescription. For example, we could try to keep the order as in (3.14) and (3.15), or we could try to impose normal ordering, with all $a^{\dagger}$’s moved to the left and all $a$’s moved to the right. None of these prescriptions produce valid quantum conservation laws that commute with (1.1). In fact, we know that the classical conservation law $G_{2}$ can be rewritten in terms of a structure identical to the first term of the valid quantum conservation law (1.3). (We will provide a brute force proof that (1.3) commutes with (1.1) in section 6.) This makes us expect that other conservation laws will behave similarly: they will consist of the highest-order polynomial term that can be visualized as the corresponding classical conservation law with $(\bar{\alpha}_{n},\alpha_{n})$ replaced with $(a_{n}^{\dagger},a_{n})$, normal-ordered, plus lower-order polynomial quantum corrections. This is precisely the structure of (1.3), and it will be true of all other conservation laws we can explicitly construct. What is even more striking is that, for all quantum conservation laws we construct, we empirically find purely integer eigenvalue spectra. This will be reported in section 12, together with our knowledge about the rather peculiar quantum Lax pair, also obtained by introducing a quantum correction to the classical Lax pair (3.7).

Before proceeding with more in-depth studies of the GG Hamiltonian, it will be handy to assemble here a few useful operators and state the algebra they form. In much of our exposition, we will keep in the back of our mind the picture of diagonalizing simultaneously $H$ given by (1.1) and $H_{\min}$ given by (1.3). First, this will remove some of the possible degeneracies present when $H$ is diagonalized by itself. Second, there are nice algebraic relations that give this simultaneous diagonalization process a useful purpose.

We first remark that $a_{0}^{\dagger}$ commutes with $H_{\min}$, since the latter only involves modes number 1 and higher. For that reason, for any eigenvector $|\Psi\rangle$ of $H_{\min}$ in block $(N,M)$, we can take $a_{0}^{\dagger}|\Psi\rangle$ in block $(N+1,M)$ and it will be an eigenvector of $H_{\min}$ as well. Furthermore, we can introduce quantum shift operators $S$ and $S^{\dagger}$ analogous to (3.21):

$$S|\eta_{0},\eta_{1},\cdots\rangle=|0,\eta_{0},\eta_{1},\cdots\rangle,\qquad S^{\dagger}|\eta_{0},\eta_{1},\eta_{2}\cdots\rangle=|\eta_{1},\eta_{2}\cdots\rangle.$$ (4.1)

We define $S^{\dagger}$ to annihilate $|\eta_{0},\eta_{1},\cdots\rangle$ unless $\eta_{0}=0$. With this, we have

$$S^{\dagger}S=1,\qquad SS^{\dagger}=P_{0},\qquad P_{0}S=S,\qquad S^{\dagger}P_{0}=S^{\dagger},$$ (4.2)

where $P_{0}$ is the projector onto the null space of $a_{0}$, as well as the following useful relations:

$$SMS^{\dagger} =M-N,$$ (4.3)

$$S^{\dagger}HS =H,$$ (4.4)

$$H_{\min}S =S(2H+H_{\min}+2M+N).$$ (4.5)

The last relation means that, if we have simultaneously diagonalized $H$ and $H_{\min}$ in block $(N,M)$, and $|\Psi\rangle$ is one of their joint eigenvectors, the vector $S|\Psi\rangle$ residing in the block $(N,M+N)$ is an eigenvector of $H_{\min}$ with eigenvalue determined by (4.5).

The above picture suggests a recursive process for diagonalizing $H$ and $H_{\min}$: imagine we have already diagonalized these two operators in the blocks whose values of $N$ or $M$ are below the current block $(N,M)$. Then, take the eigenvectors from the block $(k,M-k)$, denoted as $V_{k,M-k}$, and transport them to the block ($k,M$) with $S$. In order to reach the block ($N,M$), one can apply $a_{0}^{\dagger}$ acting with it $N-k$ times. We obtain groups of (mutually) linearly independent vectors, all of which are eigenvectors of $H_{\min}$. Furthermore, the disjoint union of these vectors

$$\bigcup_{k=1}^{\min(N,M)}\frac{a_{0}^{\dagger N-k}}{\sqrt{(N-k)!}}SV_{k,M-k}$$ (4.6)

forms a complete basis in the block $(N,M)$. This follows directly from the linear independence of these states, combined with a counting argument relating the dimensions of the spaces $V_{k,M-k}$ and the size of the block $(N,M)$. The number of vectors in an ($N,M$) block is the number of integer partitions of $M$ into at most $N$ parts [6] denoted $p_{N}(M)$. This number-theoretic function satisfies the recursion relation

$$p_{N}(M)=p_{N-1}(M)+p_{N}(M-N).$$ (4.7)

(A partition of $M$ into at most $N$ parts is either a partition into exactly $N$ parts, so that no parts are zero and we can subtract 1 from each part, obtaining a partition of $M-N$ into at most $N$ parts, or it is a partition of $M$ into at most $N-1$ part.) Applying this relation recursively to the first term on the right-hand side of (4.7) yields

$$p_{N}(M)=\sum_{k=1}^{\min(N,M)}p_{k}(M-k)\,.$$ (4.8)

This shows that (4.6) contains precisely enough linearly independent vectors to be a basis for the block $(N,M)$. We have thus obtained, by this very simple process, an eigenbasis of $H_{\min}$ in the current block. The only issue is that, to restart the iterative procedure we need a joint eigenbasis of $H$ and $H_{\min}$, which may require re-diagonalization within the degenerate subspaces of $H_{\min}$. So we must inspect the spectrum and rediagonalize $H$ within any such degenerate subspaces.

We will see in sections 8 and 10 that, in application to the two highest subspaces of $H_{\min}$, this procedure gives an explicit construction of families of eigenvectors. The entry point of this construction is the joint top eigenvectors of $H$ and $H_{\min}$ that we will build in section 7. There is a relation between these top eigenspaces of $H_{\min}$ and classical invariant manifolds (3.5), as we shall explain in section 11. Before proceeding with that, we will construct, in section 5, bounds on $H$ and $H_{\min}$ that will help us identify these top subspaces. For future use we define, in addition to the operators introduced above, the following family of operators

$$J_{1}\equiv a_{0},\qquad J_{2n+1}=[J_{2n-1},H].$$ (4.9)

These are quantum analogs of the classical quantities with the same name in [1].

5 Energy bounds

The bounds on $H$ and $H_{\min}$ will play an essential role in our constructions of explicit families of eigenvectors below. They are also important more broadly in controlling the spectral properties of $H$ and $H_{\min}$, and will reveal useful algebraic structures. Consider first the following decomposition:

$$(N-1) (N+2M)-2H=\sum_{kl}(1+k+l)a^{\dagger}_{k}a^{\dagger}_{l}a_{k}a_{l}-2H$$

$$=\sum_{j=0}^{\infty}\left[(1+j)\sum_{k=0}^{j}a^{\dagger}_{k}a^{\dagger}_{j-k}a_{k}a_{j-k}-\sum_{k,l=0}^{j}a^{\dagger}_{k}a^{\dagger}_{j-k}a_{l}a_{j-l}\right]=\sum_{j=0}^{\infty}\sum_{\alpha=1}^{\lfloor{j/2}\rfloor}B_{j\alpha}^{\dagger}B_{j\alpha},$$ (5.1)

where

$$B_{j\alpha}\equiv\sum_{k=0}^{j}e^{(j,\alpha)}_{k}a_{k}a_{j-k},$$ (5.2)

with $e^{(j,\alpha)}_{k}$ being any set of vectors satisfying

$$e^{(j,\alpha)}_{k}=e^{(j,\alpha)}_{j-k},\qquad\sum_{k}e^{(j,\alpha)}_{k}=0,\qquad\sum_{\alpha=1}^{\lfloor{j/2}\rfloor}e^{(j,\alpha)}_{k}e^{(j,\alpha)}_{l}=(1+j)\frac{(\delta_{kl}+\delta_{k,j-l})}{2}-1.$$ (5.3)

(One can check that the bisymmetric matrix appearing on the right-hand side is semi-positive-definite and admits the stated decomposition.) Since the right-hand side of (5) is evidently positive-definite, we get the bound

$$2H\leq(N-1)(N+2M).$$ (5.4)

This bound is saturated on a vector $|\Psi\rangle$ if and only if

$$\forall j,\alpha:\quad B_{j\alpha}|\Psi\rangle=0.$$ (5.5)

We will use this and other similar formulas below to construct specific families of eigenvectors. Note that $(N-1)(N+2M)$ is nothing but the normal-ordered product of $N$ and $N+2M$.

The Hamiltonian $H_{\min}$ may be visualized in terms of a ‘layered cake’ construction: we first take $H$ and shift all the mode indices of the $a$’s by +1, then by +2, +3 and so on, and add up the results, obtaining

$$H_{\min}=\sum_{n+m=k+l}^{n,m,k,l\geq 0}\sum_{s=1}^{\infty}a^{\dagger}_{n+s}a^{\dagger}_{m+s}a_{k+s}a_{l+s}+\sum_{k=1}^{\infty}k^{2}a^{\dagger}_{k}a_{k}.$$ (5.6)

Indeed, a given term $a^{\dagger}_{\bar{n}}a^{\dagger}_{\bar{m}}a_{\bar{k}}a_{\bar{l}}$ will be present in the above sum for $s=1..\min(\bar{n},\bar{m},\bar{k},\bar{l})$, yielding the coefficient in (1.3) upon summation over $s$. We can then apply the decomposition (5) within this ‘layered cake’ construction:

$$\sum_{n+m=k+l}^{n,m,k,l\geq 0}a^{\dagger}_{n+s}a^{\dagger}_{m+s}a_{k+s}a_{l+s}=\sum_{j=0}^{\infty}(1+j)\sum_{k=0}^{j}a^{\dagger}_{k+s}a^{\dagger}_{j-k+s}a_{k+s}a_{j-k+s}-\sum_{j=0}^{\infty}\sum_{\alpha=1}^{\lfloor{j/2}\rfloor}B_{js\alpha}^{\dagger}B_{js\alpha},$$ (5.7)

where

$$B_{js\alpha}\equiv\sum_{k=0}^{j}e^{(j,\alpha)}_{k}a_{k+s}a_{j-k+s}.$$ (5.8)

Now,

$$\sum_{s=1}^{\infty}\sum_{j=0}^{\infty}(1+j)\sum_{k=0}^{j}a^{\dagger}_{k+s}a^{\dagger}_{j-k+s}a_{k+s}a_{j-k+s}=\sum_{s=1}^{\infty}\sum_{k,l=0}^{\infty}(1+k+l)a^{\dagger}_{k+s}a^{\dagger}_{l+s}a_{k+s}a_{l+s}$$ (5.9)

$$=\sum_{s=1}^{\infty}\sum_{k,l=s}^{\infty}(1+k+l-2s)a^{\dagger}_{k}a^{\dagger}_{l}a_{k}a_{l}=\sum_{k,l=1}^{\infty}a^{\dagger}_{k}a^{\dagger}_{l}a_{k}a_{l}\sum_{s=1}^{\min(k,l)}(1+k+l-2s).$$ (5.10)

For the last sum, we get

$$\sum_{s=1}^{\min(k,l)}(1+k+l-2s) =(1+k+l)\min(k,l)-(\min(k,l)+1)\min(k,l)$$ (5.11)

$$=[k+l-\min(k,l)]\min(k,l)=\max(k,l)\min(k,l)=kl.$$

Furthermore,

$$\sum_{k,l=1}^{\infty}kl\,a^{\dagger}_{k}a^{\dagger}_{l}a_{k}a_{l}+\sum_{k=1}^{\infty}k^{2}a^{\dagger}_{k}a_{k}=M^{2}$$ (5.12)

Putting everything together,

$$H_{\min}=M^{2}-\sum_{s=1}^{\infty}\sum_{j=0}^{\infty}\sum_{\alpha=1}^{\lfloor{j/2}\rfloor}B_{js\alpha}^{\dagger}B_{js\alpha}.$$ (5.13)

Not all of the conditions $B_{js\alpha}|\Psi\rangle=0$ are independent. For example, enforcing this condition for $s=1$ and all $j$ and $\alpha$ automatically enforces it for all higher values of $s$. For such vectors $H_{\min}=M^{2}$. If, on the other hand, the conditions hold for $s=2$ but not for $s=1$, they hold automatically for all $s\geq 2$, and only terms with $s=1$ are left in (5.13). These terms can furthermore be simplified using $B_{js\alpha}|\Psi\rangle=0$ with $s\geq 2$ since that should eliminate subsets of terms in $B_{j1\alpha}$. We will use such subspaces defined by $B_{js\alpha}|\Psi\rangle=0$ for diagonalizing $H$ and $H_{\min}$.

In (5), we used a complete-basis decompositions of quadratic forms, ending up with somewhat awkward expressions in terms of the eigenvectors $e_{k}$. In many ways, it is more convenient to take advantage of the formula

$$(1+j)\sum_{k=0}^{j}\bar{x}_{k}x_{k}-\sum_{k,l=0}^{j}\bar{x}_{k}x_{l}=\frac{1}{2}\sum_{k,l=0}^{j}(\bar{x}_{k}x_{k}+\bar{x}_{l}x_{l}-\bar{x}_{k}x_{l}-\bar{x}_{l}x_{k})=\frac{1}{2}\sum_{k,l=0}^{j}|x_{k}-x_{l}|^{2},$$ (5.14)

and write instead of (5)

$$2H=(N-1)(N+2M)-\frac{1}{2}\sum_{j=0}^{\infty}\sum_{k,l=0}^{j}(a_{k}a_{j-k}-a_{l}a_{j-l})^{\dagger}(a_{k}a_{j-k}-a_{l}a_{j-l}).$$ (5.15)

An analogous representation of the bound (5.13) is

$$H_{\min}=M^{2}-\frac{1}{2}\sum_{s=1}^{\infty}\sum_{j=0}^{\infty}\sum_{k,l=0}^{j}(a_{k+s}a_{j-k+s}-a_{l+s}a_{j-l+s})^{\dagger}(a_{k+s}a_{j-k+s}-a_{l+s}a_{j-l+s}).$$ (5.16)

Introducing $\tilde{k}=k+s$, $\tilde{l}=l+s$, $\tilde{j}=j+2s$ and dropping tildes, we get

$$H_{\min}=M^{2}-\frac{1}{2}\sum_{s=1}^{\infty}\sum_{j=2s}^{\infty}\sum_{k,l=s}^{j-s}(a_{k}a_{j-k}-a_{l}a_{j-l})^{\dagger}(a_{k}a_{j-k}-a_{l}a_{j-l}).$$ (5.17)

The summand is now $s$-independent, so the sums can be simplified. The index range conditions $s\leq k\leq j-s$ and $s\leq l\leq j-s$ can be recast as $s\leq k$, $s\leq j-k$, $s\leq l$, $s\leq j-l$ or, equivalently, $s\leq\min(k,j-k,l,j-l)$, and then

$$H_{\min}=M^{2}-\frac{1}{2}\sum_{j=0}^{\infty}\sum_{k,l=0}^{j}\min(k,j-k,l,j-l)(a_{k}a_{j-k}-a_{l}a_{j-l})^{\dagger}(a_{k}a_{j-k}-a_{l}a_{j-l}).$$ (5.18)

It will be economical for some purposes to switch to thinking in terms of the following two Hamiltonians:

$$H_{0}=\frac{1}{2}\sum_{j=0}^{\infty}\sum_{k,l=0}^{j}D^{\dagger}_{jkl}D_{jkl},$$ (5.19)

and

$$H_{1}=\frac{1}{2}\sum_{j=0}^{\infty}\sum_{k,l=0}^{j}\min(k,j-k,l,j-l)D^{\dagger}_{jkl}D_{jkl},$$ (5.20)

with

$$D_{jkl}\equiv a_{k}a_{j-k}-a_{l}a_{j-l}.$$ (5.21)

Evidently,

$$2H=(N-1)(N+2M)-H_{0},\qquad H_{\min}=M^{2}-H_{1}.$$ (5.22)

An advantage is that $H_{1}$ is purely quartic, in contrast with (1.3) that has a quadratic ‘quantum’ piece. Viewed from the vantage point of the above formulas, this quadratic piece has a trivial origin: it comes from normal-ordering $M^{2}$. We will rely on $H_{0}$ and $H_{1}$ in the next section, and prove that they commute, which will automatically imply that $[H,H_{\min}]=0$.

We remark that sum-of-squares decompositions like (5.19) and (5.20) are very typical of the ‘factorization method’ and supersymmetric quantum mechanics [46]. However, since $D$ and $D^{\dagger}$ do not commute for different index values, it is not immediately obvious how to use it for constructing a full solution.

6 A brute force proof that $[H,H_{\min}]=0$

Since we use simultaneous diagonalization of $H$ and $H_{\min}$ in our constructions of the eigenvectors below, the property $[H,H_{\min}]=0$ is essential for our considerations. The Hamiltonian $H_{\min}$ was originally guessed in [11] as an outcome of a numerical procedure that algorithmically constructs conservation laws polynomial in the creation-annihilation operators as explicit matrices with the individual $(N,M)$-blocks of the Hilbert space. It is easy to verify within any given $(N,M)$-block that the two operators indeed commute, but what we need is a general analytic proof.

While we expect that, eventually, a systematic understanding of the conservation laws of $H$ will emerge in the spirit of Lax theory, we are far from that point, and a practical solution is to provide a brute force proof that $[H,H_{\min}]=0$. It is more convenient to prove instead that $[H_{0},H_{1}]=0$, with the definitions (5.19-5.20), which would automatically imply $[H,H_{\min}]=0$ in view of (5.22). We provide a proof below. It is essential for our subsequent considerations, but rather bulky and tedious. The rest of our treatment will only use the final result, but not the details of this proof, so the details can be freely skipped.

We remark that we could in principle throw away the top polynomial piece of $[H_{0},H_{1}]$ after normal ordering, because that piece must agree with the classical computation, and in the classical theory $H$ and $H_{\min}$ Poisson-commute from the general Lax theory described in section 3. We opt, however, for a completely mechanical proof that follows through detailed bookkeeping of all the terms emerging from the commutation and proceeds to demonstrate that the result is zero.

We introduce the shorthand $\sum_{\{jkl\}}\equiv\sum_{j=0}^{\infty}\sum_{k,l=0}^{j}$. Then,

$$[H_{0},H_{1}]=\frac{1}{4}\sum_{\{jkl\}}\sum_{\{j^{\prime}k^{\prime}l^{\prime}\}}\mathrm{min}_{j^{\prime}k^{\prime}l^{\prime}}[D^{\dagger}_{jkl}D_{jkl},D^{\dagger}_{j^{\prime}k^{\prime}l^{\prime}}D_{j^{\prime}k^{\prime}l^{\prime}}],$$ (6.1)

with the definition

$$\mathrm{min}_{jkl}\equiv\min(k,j-k,l,j-l).$$ (6.2)

We have

$$[D^{\dagger}_{jkl}D_{jkl},D^{\dagger}_{j^{\prime}k^{\prime}l^{\prime}}D_{j^{\prime}k^{\prime}l^{\prime}}] =D^{\dagger}_{jkl}D_{jkl}D^{\dagger}_{j^{\prime}k^{\prime}l^{\prime}}D_{j^{\prime}k^{\prime}l^{\prime}}-D^{\dagger}_{j^{\prime}k^{\prime}l^{\prime}}D_{j^{\prime}k^{\prime}l^{\prime}}D^{\dagger}_{jkl}D_{jkl}$$

$$=D^{\dagger}_{jkl}[D_{jkl},D^{\dagger}_{j^{\prime}k^{\prime}l^{\prime}}]D_{j^{\prime}k^{\prime}l^{\prime}}-D^{\dagger}_{j^{\prime}k^{\prime}l^{\prime}}[D_{j^{\prime}k^{\prime}l^{\prime}},D^{\dagger}_{jkl}]D_{jkl}.$$

So,

$$[H_{0},H_{1}]=\sum_{\{jkl\}}\sum_{\{j^{\prime}k^{\prime}l^{\prime}\}}\mathrm{min}_{j^{\prime}k^{\prime}l^{\prime}}\left(D^{\dagger}_{jkl}[a_{k}a_{j-k},a^{\dagger}_{k^{\prime}}a^{\dagger}_{j^{\prime}-k^{\prime}}]D_{j^{\prime}k^{\prime}l^{\prime}}-D^{\dagger}_{j^{\prime}k^{\prime}l^{\prime}}[a_{k^{\prime}}a_{j^{\prime}-k^{\prime}},a^{\dagger}_{k}a^{\dagger}_{j-k}]D_{jkl}\right),$$

where we have used the antisymmetry of $D_{jkl}$ with respect to permuting $k$ and $l$. Now,

$$[a_{k}a_{j-k},a^{\dagger}_{k^{\prime}}a^{\dagger}_{j^{\prime}-k^{\prime}}] =\delta_{kk^{\prime}}a_{j-k}a^{\dagger}_{j^{\prime}-k^{\prime}}+\delta_{k,j^{\prime}-k^{\prime}}a_{j-k}a^{\dagger}_{k^{\prime}}+\delta_{j-k,k^{\prime}}a_{k}a^{\dagger}_{j^{\prime}-k^{\prime}}+\delta_{j-k,j^{\prime}-k^{\prime}}a_{k}a^{\dagger}_{k^{\prime}}$$

$$=\delta_{kk^{\prime}}a^{\dagger}_{j^{\prime}-k^{\prime}}a_{j-k}+\delta_{k,j^{\prime}-k^{\prime}}a^{\dagger}_{k^{\prime}}a_{j-k}+\delta_{j-k,k^{\prime}}a^{\dagger}_{j^{\prime}-k^{\prime}}a_{k}+\delta_{j-k,j^{\prime}-k^{\prime}}a^{\dagger}_{k^{\prime}}a_{k}$$

$$+2\delta_{kk^{\prime}}\delta_{jj^{\prime}}+2\delta_{k,j-k^{\prime}}\delta_{jj^{\prime}}.$$

Plugging this back in and keeping in mind that $j-k$ can be interchanged with $k$ using the symmetry of $D_{jkl}$, and similarly for the primed indices, we get

$$\begin{split}\frac{[H_{0},H_{1}]}{4}&=\sum_{\{jkl\}}\sum_{\{j^{\prime}k^{\prime}l^{\prime}\}}\mathrm{min}_{j^{\prime}k^{\prime}l^{\prime}}\big[D^{\dagger}_{jkl}(\delta_{kk^{\prime}}a^{\dagger}_{j^{\prime}-k^{\prime}}a_{j-k}+\delta_{kk^{\prime}}\delta_{jj^{\prime}})D_{j^{\prime}k^{\prime}l^{\prime}}\\ &\hskip 85.35826pt-D^{\dagger}_{j^{\prime}k^{\prime}l^{\prime}}(\delta_{kk^{\prime}}a^{\dagger}_{j-k}a_{j^{\prime}-k^{\prime}}+\delta_{kk^{\prime}}\delta_{jj^{\prime}})D_{jkl}\big]\\ &=\sum_{j,j^{\prime}=0}^{\infty}\sum_{l=0}^{j}\sum_{l^{\prime}=0}^{j^{\prime}}\sum_{k=0}^{\min(j,j^{\prime})}\mathrm{min}_{j^{\prime}kl^{\prime}}\big[D^{\dagger}_{jkl}a^{\dagger}_{j^{\prime}-k}a_{j-k}D_{j^{\prime}kl^{\prime}}-(jl\leftrightarrow j^{\prime}l^{\prime})\big]\\ &\hskip 56.9055pt+\sum_{j=0}^{\infty}\sum_{k,l,l^{\prime}=0}^{j}\mathrm{min}_{jkl^{\prime}}\big[D^{\dagger}_{jkl}D_{jkl^{\prime}}-D^{\dagger}_{jkl^{\prime}}D_{jkl}\big].\end{split}$$ (6.3)

Take the last line:

$$\sum_{j=0}^{\infty}\sum_{k,l,l^{\prime}=0}^{j}\mathrm{min}_{jkl^{\prime}}\big[D^{\dagger}_{jkl}D_{jkl^{\prime}}-D^{\dagger}_{jkl^{\prime}}D_{jkl}\big]=\sum_{j=0}^{\infty}\sum_{k,l,l^{\prime}=0}^{j}\big[\mathrm{min}_{jkl^{\prime}}-\mathrm{min}_{jkl}\big]D^{\dagger}_{jkl}D_{jkl^{\prime}}$$

$$=\sum_{j=0}^{\infty}\sum_{k,l,l^{\prime}=0}^{j}\big[\mathrm{min}_{jkl^{\prime}}-\mathrm{min}_{jkl}\big]\big(a^{\dagger}_{k}a^{\dagger}_{j-k}a_{k}a_{j-k}-a^{\dagger}_{l}a^{\dagger}_{j-l}a_{k}a_{j-k}-a^{\dagger}_{k}a^{\dagger}_{j-k}a_{l^{\prime}}a_{j-l^{\prime}}+a^{\dagger}_{l}a^{\dagger}_{j-l}a_{l^{\prime}}a_{j-l^{\prime}}\big).$$

The first term in the round brackets does not depend on $l$ and $l^{\prime}$, while $\sum_{ll^{\prime}}(\mathrm{min}_{jkl^{\prime}}-\mathrm{min}_{jkl})$ is evidently 0, thereby eliminating the first term in the round brackets. For the second term in the round brackets, we do the substitution $k\to l^{\prime}\to l\to k$, which gives

$$\sum_{j=0}^{\infty}\sum_{k,l,l^{\prime}=0}^{j}\big[\mathrm{min}_{jll^{\prime}}-\mathrm{min}_{jkl^{\prime}}\big]a^{\dagger}_{k}a^{\dagger}_{j-k}a_{l^{\prime}}a_{j-l^{\prime}}.$$

This combines nicely with the third term in the round brackets, leaving

$$\sum_{j=0}^{\infty}\sum_{k,l,l^{\prime}=0}^{j}\mathrm{min}_{jkl^{\prime}}\big[D^{\dagger}_{jkl}D_{jkl^{\prime}}-D^{\dagger}_{jkl^{\prime}}D_{jkl}\big]$$

$$=-\sum_{j=0}^{\infty}\sum_{k,l,l^{\prime}=0}^{j}\big[\mathrm{min}_{jll^{\prime}}-\mathrm{min}_{jkl}\big]a^{\dagger}_{k}a^{\dagger}_{j-k}a_{l^{\prime}}a_{j-l^{\prime}}+\sum_{j=0}^{\infty}\sum_{k,l,l^{\prime}=0}^{j}\big[\mathrm{min}_{jkl^{\prime}}-\mathrm{min}_{jkl}]a^{\dagger}_{l}a^{\dagger}_{j-l}a_{l^{\prime}}a_{j-l^{\prime}}.$$

Finally, we interchange $k$ and $l$ in the first sum, making the two sums cancel each other. We have thus proved that the last line in (6.3) vanishes. Note that, up to this point, we have not used the explicit form of $\mathrm{min}_{jkl}$ but only its index permutation symmetries.

We now turn to the first line of the last representation in (6.3):

$$\sum_{j,j^{\prime}=0}^{\infty}\sum_{l=0}^{j}\sum_{l^{\prime}=0}^{j^{\prime}}\sum_{k=0}^{\min(j,j^{\prime})}(\mathrm{min}_{j^{\prime}kl^{\prime}}-\mathrm{min}_{jkl})D^{\dagger}_{jkl}a^{\dagger}_{j^{\prime}-k}a_{j-k}D_{j^{\prime}kl^{\prime}}$$

$$=\sum_{j,j^{\prime}=0}^{\infty}\sum_{l=0}^{j}\sum_{l^{\prime}=0}^{j^{\prime}}\sum_{k=0}^{\min(j,j^{\prime})}(\mathrm{min}_{j^{\prime}kl^{\prime}}-\mathrm{min}_{jkl})(a^{\dagger}_{k}a^{\dagger}_{j-k}a^{\dagger}_{j^{\prime}-k}a_{j-k}a_{k}a_{j^{\prime}-k}-a^{\dagger}_{k}a^{\dagger}_{j-k}a^{\dagger}_{j^{\prime}-k}a_{j-k}a_{l^{\prime}}a_{j^{\prime}-l^{\prime}}$$

$$\hskip 142.26378pt-a^{\dagger}_{l}a^{\dagger}_{j-l}a^{\dagger}_{j^{\prime}-k}a_{j-k}a_{k}a_{j^{\prime}-k}+a^{\dagger}_{l}a^{\dagger}_{j-l}a^{\dagger}_{j^{\prime}-k}a_{j-k}a_{l^{\prime}}a_{j^{\prime}-l^{\prime}})$$

$$=\sum_{k=0}^{\infty}\sum_{j,j^{\prime}=k}^{\infty}\sum_{l=0}^{j}\sum_{l^{\prime}=0}^{j^{\prime}}(\mathrm{min}_{j^{\prime}kl^{\prime}}-\mathrm{min}_{jkl})(a^{\dagger}_{k}a^{\dagger}_{j-k}a^{\dagger}_{j^{\prime}-k}a_{j-k}a_{k}a_{j^{\prime}-k}-a^{\dagger}_{k}a^{\dagger}_{j-k}a^{\dagger}_{j^{\prime}-k}a_{j-k}a_{l^{\prime}}a_{j^{\prime}-l^{\prime}}$$

$$\hskip 142.26378pt-a^{\dagger}_{l}a^{\dagger}_{j-l}a^{\dagger}_{j^{\prime}-k}a_{j-k}a_{k}a_{j^{\prime}-k}+a^{\dagger}_{l}a^{\dagger}_{j-l}a^{\dagger}_{j^{\prime}-k}a_{j-k}a_{l^{\prime}}a_{j^{\prime}-l^{\prime}})$$

$$=\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j^{\prime}+k}(\mathrm{min}_{j^{\prime}+k,kl^{\prime}}-\mathrm{min}_{j+k,kl})[a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{k}a_{j}a_{j^{\prime}}-a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{l^{\prime}}a_{j}a_{j^{\prime}+k-l^{\prime}}$$

$$\hskip 142.26378pt-a^{\dagger}_{l}a^{\dagger}_{j^{\prime}}a^{\dagger}_{j+k-l}a_{k}a_{j}a_{j^{\prime}}+a^{\dagger}_{l}a^{\dagger}_{j^{\prime}}a^{\dagger}_{j+k-l}a_{l^{\prime}}a_{j}a_{j^{\prime}+k-l^{\prime}}].$$ (6.4)

The first term in the square brackets vanishes using the permutation $jl\leftrightarrow j^{\prime}l^{\prime}$. Then, take the third term:

$$\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j^{\prime}+k}(\mathrm{min}_{j^{\prime}+k,kl^{\prime}}-\mathrm{min}_{j+k,kl})a^{\dagger}_{l}a^{\dagger}_{j^{\prime}}a^{\dagger}_{j+k-l}a_{k}a_{j}a_{j^{\prime}}$$

$$=\sum_{j,j^{\prime},l=0}^{\infty}\sum_{k=\max(0,l-j)}^{\infty}\sum_{l^{\prime}=0}^{j^{\prime}+k}(\mathrm{min}_{j^{\prime}+k,kl^{\prime}}-\mathrm{min}_{j+k,kl})a^{\dagger}_{l}a^{\dagger}_{j^{\prime}}a^{\dagger}_{j+k-l}a_{k}a_{j}a_{j^{\prime}}$$

$$=\sum_{j,j^{\prime},l=0}^{\infty}\sum_{k=\max(j-l,0)}^{\infty}\sum_{l^{\prime}=0}^{j^{\prime}+k+l-j}(\mathrm{min}_{j^{\prime}+k+l-j,k+l-j,l^{\prime}}-\mathrm{min}_{l+k,k+l-j,l})a^{\dagger}_{l}a^{\dagger}_{j^{\prime}}a^{\dagger}_{k}a_{k+l-j}a_{j}a_{j^{\prime}}$$

$$=\sum_{k,l,j^{\prime}=0}^{\infty}\sum_{j=0}^{k+l}\sum_{l^{\prime}=0}^{j^{\prime}+k+l-j}(\mathrm{min}_{j^{\prime}+k+l-j,k+l-j,l^{\prime}}-\mathrm{min}_{l+k,k+l-j,l})a^{\dagger}_{l}a^{\dagger}_{j^{\prime}}a^{\dagger}_{k}a_{k+l-j}a_{j}a_{j^{\prime}}$$

$$=\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j+j^{\prime}+k-l}(\mathrm{min}_{j+j^{\prime}+k-l,j+k-l,l^{\prime}}-\mathrm{min}_{j+k,j+k-l,j})a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{j+k-l}a_{l}a_{j^{\prime}}$$

$$=\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j+j^{\prime}+k-l}(\mathrm{min}_{j+j^{\prime}+k-l,l^{\prime}j^{\prime}}-\mathrm{min}_{j+k,kl})a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{j+k-l}a_{l}a_{j^{\prime}}$$

$$=\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j^{\prime}+l}(\mathrm{min}_{j^{\prime}+l,ll^{\prime}}-\mathrm{min}_{j+k,kl})a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{j+k-l}a_{l}a_{j^{\prime}}$$

$$=\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l^{\prime}=0}^{j^{\prime}+k}\sum_{l=0}^{j+l^{\prime}}(\mathrm{min}_{j+l^{\prime},ll^{\prime}}-\mathrm{min}_{j^{\prime}+k,kl^{\prime}})a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{l^{\prime}}a_{j}a_{j^{\prime}+k-l^{\prime}}$$

$$=\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l^{\prime}=0}^{j^{\prime}+k}(jl^{\prime}-(j+l^{\prime}+1)\mathrm{min}_{j^{\prime}+k,kl^{\prime}})a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{l^{\prime}}a_{j}a_{j^{\prime}+k-l^{\prime}}$$

$$=\frac{1}{2}\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l^{\prime}=0}^{j^{\prime}+k}[j(j^{\prime}+k)-(2j+j^{\prime}+k+2)\mathrm{min}_{j^{\prime}+k,kl^{\prime}}]a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{l^{\prime}}a_{j}a_{j^{\prime}+k-l^{\prime}},$$

where we have used

$$\sum_{k=0}^{j}\mathrm{min}_{jkl}=\sum_{k=0}^{j}\min(k,j-k,l,j-l)=l(j-l),$$ (6.5)

and then symmetrized with respect to $l^{\prime}\leftrightarrow j^{\prime}+k-l^{\prime}$. The second term in (6.4) gives by direct evaluation

$$\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j^{\prime}+k}(\mathrm{min}_{j^{\prime}+k,kl^{\prime}}-\mathrm{min}_{j+k,kl})a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{l^{\prime}}a_{j}a_{j^{\prime}+k-l^{\prime}}$$

$$=\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l^{\prime}=0}^{j^{\prime}+k}[(j+k+1)\mathrm{min}_{j^{\prime}+k,kl^{\prime}}-jk]a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{l^{\prime}}a_{j}a_{j^{\prime}+k-l^{\prime}}.$$

Adding this up with the third term we get,

$$\frac{1}{2}\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l^{\prime}=0}^{j^{\prime}+k}[j(j^{\prime}-k)-(j^{\prime}-k)\mathrm{min}_{j^{\prime}+k,kl^{\prime}}]a^{\dagger}_{k}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a_{l^{\prime}}a_{j}a_{j^{\prime}+k-l^{\prime}},$$ (6.6)

but this expression equal minus itself under the interchange of $k$ and $j^{\prime}$, and hence the sum of the second and third terms in (6.4) vanishes.

Finally, take the last term in (6.4) and process it in a manner similar to the third term, so that the indices of the three $a^{\dagger}$’s are $k$, $j$ and $j^{\prime}$:

$$\sum_{k,j,j^{\prime}=0}^{\infty}\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j^{\prime}+k}(\mathrm{min}_{j^{\prime}+k,kl^{\prime}}-\mathrm{min}_{j+k,kl})a^{\dagger}_{l}a^{\dagger}_{j^{\prime}}a^{\dagger}_{j+k-l}a_{l^{\prime}}a_{j}a_{j^{\prime}+k-l^{\prime}}$$

$$=\sum_{j,j^{\prime},l=0}^{\infty}\sum_{k=\max(0,l-j)}^{\infty}\sum_{l^{\prime}=0}^{j^{\prime}+k}(\mathrm{min}_{j^{\prime}+k,kl^{\prime}}-\mathrm{min}_{j+k,kl})a^{\dagger}_{l}a^{\dagger}_{j^{\prime}}a^{\dagger}_{j+k-l}a_{l^{\prime}}a_{j}a_{j^{\prime}+k-l^{\prime}}$$

$$=\sum_{j,j^{\prime},l=0}^{\infty}\sum_{k=\max(j-l,0)}^{\infty}\sum_{l^{\prime}=0}^{j^{\prime}+k+l-j}(\mathrm{min}_{j^{\prime}+k+l-j,k+l-j,l^{\prime}}-\mathrm{min}_{k+l,k+l-j,l})a^{\dagger}_{l}a^{\dagger}_{j^{\prime}}a^{\dagger}_{k}a_{l^{\prime}}a_{j}a_{j^{\prime}+l+k-j-l^{\prime}}$$

$$=\sum_{k,l,j^{\prime}=0}^{\infty}\sum_{j=0}^{l+k}\sum_{l^{\prime}=0}^{j^{\prime}+k+l-j}(\mathrm{min}_{j^{\prime}+k+l-j,k+l-j,l^{\prime}}-\mathrm{min}_{k+l,k+l-j,l})a^{\dagger}_{l}a^{\dagger}_{j^{\prime}}a^{\dagger}_{k}a_{l^{\prime}}a_{j}a_{j^{\prime}+l+k-j-l^{\prime}}$$

$$=\sum_{k,j,j^{\prime}=0}^{\infty}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a^{\dagger}_{k}\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j+j^{\prime}+k-l}(\mathrm{min}_{j+j^{\prime}+k-l,j+k-l,l^{\prime}}-\mathrm{min}_{j+k,j+k-l,j})a_{l}a_{l^{\prime}}a_{j+j^{\prime}+k-l-l^{\prime}}$$

$$=\sum_{k,j,j^{\prime}=0}^{\infty}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a^{\dagger}_{k}\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j+j^{\prime}+k-l}(\mathrm{min}_{j+j^{\prime}+k-l,j^{\prime}l^{\prime}}-\mathrm{min}_{j+k,lk})a_{l}a_{l^{\prime}}a_{j+j^{\prime}+k-l-l^{\prime}}$$

$$\equiv\sum_{k,j,j^{\prime}=0}^{\infty}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a^{\dagger}_{k}A_{kjj^{\prime}}=\frac{1}{3}\sum_{k,j,j^{\prime}=0}^{\infty}a^{\dagger}_{j}a^{\dagger}_{j^{\prime}}a^{\dagger}_{k}(A_{kjj^{\prime}}+A_{kj^{\prime}j}+A_{j^{\prime}jk}),$$

with

$$A_{kjj^{\prime}}\equiv\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j^{\prime}+l}(\mathrm{min}_{j^{\prime}+l,ll^{\prime}}-\mathrm{min}_{j+k,lk})a_{j+k-l}a_{l^{\prime}}a_{j^{\prime}+l-l^{\prime}}.$$ (6.7)

One can easily verify using computer algebra for any specific $k$, $j$ and $j^{\prime}$ that $A_{kjj^{\prime}}+A_{kj^{\prime}j}+A_{j^{\prime}jk}=0$, but proving this directly appears challenging because of the complexity of the summation regions. This difficulty can be bypassed by converting sums into integrals using generating functions, and we present a complete proof in Appendix B. With this, we have established that the last term in the square brackets of (6.4) gives a vanishing contribution, while all the other terms have been proved earlier to give vanishing contributions as well, hence $[H_{0},H_{1}]=0$, and therefore $[H,H_{\min}]=0$.

7 The highest energy states

Each $(N,M)$-block contains one eigenvector with the highest possible eigenvalues of $H$ and $H_{\min}$ defined by the bounds (5.15) and (5.18). These eigenvalues are $(N-1)(N+2M)/2$ and $M^{2}$, respectively. The corresponding vector $|\psi^{\max}_{N,M}\rangle$ must satisfy, by (5.15) and (5.18),

$$a_{k}a_{j-k}|\psi^{\max}_{N,M}\rangle=a_{l}a_{j-l}|\psi^{\max}_{N,M}\rangle$$ (7.1)

for all $j$, $k$ and $l$. We note that the constraints (7.1) are preserved by the action of the following operator

$$Z=\sum_{n=0}^{\infty}(n+1)a^{\dagger}_{n+1}a_{n}\,,$$ (7.2)

which moves states from the block $(N,M)$ to $(N,M+1)$. If we assume that $\ket{\psi}$ satisfies (7.1), then $Z\ket{\psi}$ also satisfies (7.1):

$$a_{j-k}a_{k}Z|\psi\rangle =Za_{j-k}a_{k}|\psi\rangle+\left[(1-\delta_{k0})ka_{j-k}a_{k-1}+(1-\delta_{j-k,0})(j-k)a_{j-k-1}a_{k}\right]|\psi\rangle$$

$$=Za_{j-l}a_{l}|\psi\rangle+ja_{j-1}a_{0}|\psi\rangle=a_{j-l}a_{l}Z|\psi\rangle+(ja_{j-1}a_{0}-ja_{j-1}a_{0})|\psi\rangle=a_{j-l}a_{l}Z|\psi\rangle.$$

A possible strategy is therefore to solve the constraints (7.1) in low $M$ blocks and use the $Z$-operator to transport the vectors to higher $M$ blocks. For this purpose, we consider the blocks with $M=0$, which contain a single state $|N,0,...\rangle$. Such a state evidently satisfies (7.1). Hence, at least one state that saturates the bounds (5.15) and (5.18) can be found in each $(N,M)$ block by transporting the single state in block ($N,0$) to the block $(N,M)$ by the repeated action of the $Z$-operator:

$$\ket{\psi^{\max}_{N,M}}=\frac{Z^{M}|N,0,\cdots\rangle}{\sqrt{\langle N,0,\cdots|Z^{\dagger M}Z^{M}|N,0,...\rangle}}$$ (7.3)

Note that the kernel of $Z$ is trivial. Therefore, the states (7.3) define genuine eigenstates of $H$ and $H_{\min}$. The corresponding eigenvalues are $(N-1)(N+2M)/2$ and $M^{2}$, since they saturate (5.15) and (5.18). The denominator in (7.3) is relatively straightforward to compute using the identity [40]

$$Z^{\dagger m}Z^{m}=\sum_{k=0}^{m}\frac{m!^{2}}{k!(m-k)!^{2}}\frac{(k-1+2\hat{M}+\hat{N})!}{(-1+2\hat{M}+\hat{N})!}Z^{m-k}Z^{\dagger m-k}.$$ (7.4)

Acting with this on the state $|N,0,...\rangle$, the only term that contributes is $k=m$, such that

$$\ket{{\psi}^{\max}_{N,M}}=\sqrt{\frac{(N-1)!}{M!(N+M-1)!}}Z^{M}|N,0,...\rangle.$$ (7.5)

It turns out that the bounds (5.15) and (5.18) that lead to (7.1) are saturated by a single state in each $(N,M)$ block, given by (7.5). The proof for uniqueness relies on showing that all the Fock basis coefficients of a state satisfying (7.1) are related, so that $\ket{{\psi}^{\max}_{N,M}}$ is unique up to normalization. Take any two Fock basis vectors $\ket{e_{1}}$ and $\ket{e_{2}}$. There exist two products of annihilation operators $L_{1}$ and $L_{2}$ that turn them into the same vector (up to normalization) $L_{1}\ket{e_{1}}\sim L_{2}\ket{e_{2}}$ in some lower $(N^{\prime},M^{\prime})$-block. Since both vectors started in the block $(N,M)$ and have been mapped to the block $(N^{\prime},M^{\prime})$, the number of annihilation operators and the sum of their indices must be equal in both $L_{1}$ and $L_{2}$. The claim is that for a vector $\ket{\psi}$ satisfying (7.1), one has $L_{1}\ket{\psi}=L_{2}\ket{\psi}$. A constructive approach to verify this claim, is to consider extreme values $a_{j-k}a_{k}$ (with $j-k>k$) that are present in one of the $L_{i}$ but not in both, and exchange them for $a_{j-k-1}a_{k+1}$ using (7.1). One can repeat this on either side, depending on which side has the largest (or smallest) mode number in annihilation operators acting on $\ket{\psi}$. The procedure stops when the products of annihilation operators on both sides coincide, which is what we wanted to prove. The equation $L_{1}\ket{\psi}=L_{2}\ket{\psi}$ tells us in particular that the coefficients of the vectors $\ket{e_{1}}$ and $\ket{e_{2}}$ in the expansion of $\ket{\psi}$ are related. Since this argument applies to any pair of Fock vectors $\ket{e_{1}}$ and $\ket{e_{2}}$, all coefficients in $\ket{\psi}$ are related and (7.1) has a unique solution in each ($N,M$)-block.

8 Eigenstates within the top eigenspace of $H_{\min}$

In the previous section, we constructed the joint top eigenvector of $H$ and $H_{\min}$ in every ($N,M$)-block. We shall now explore the structure of the remaining eigenvectors of $H$ that reside in the top eigenspace of $H_{\min}$.

Consider the subspace of states in an $(N,M)$-block annihilated by all $B_{js\alpha}$ for $s\geq 1$ and call it $V^{(1)}_{NM}$. From (5.20) and (5.22), it follows that their eigenvalue under $H_{\min}$ equals $M^{2}$, which saturates the upper bound for $H_{\min}$ (but not necessarily for $H$). Similarly to (7.1), these states satisfy

$$a_{k+1}a_{j-k+1}|\phi^{(1)}_{N,M}\rangle=a_{l+1}a_{j-l+1}|\phi^{(1)}_{N,M}\rangle,$$ (8.1)

for any $j,k$ and $l$. Following the same line of thought as in the uniqueness proof at the end of the previous section, one can show that the dimension of this subspace is at most $\min(N,M)$. Indeed, the constraints (8.1) can be used to relate coefficients of Fock vectors that have the same number of modes $a^{\dagger}_{0}$, while vectors with a different eigenvalue for $a^{\dagger}_{0}a_{0}$ cannot be related using these constraints. This leaves $\min(N,M)$ coefficients in the Fock expansion unconstrained.

The idea is now to construct a basis for $V^{(1)}_{NM}$ and diagonalize $H$ in that $H_{\min}$ eigenspace. Since $H$ and $H_{\min}$ commute and any eigenstate outside of $V^{(1)}_{NM}$ has a lower eigenvalue for $H_{\min}$, this subspace is invariant under $H$ as well. (Note that the highest energy state derived in the previous section belongs to $V^{(1)}_{NM}$.) The key observation that allows us to construct a basis for $V^{(1)}_{NM}$ is the fact that the constraints (7.1) and (8.1) are related by the action of the shift operator $S$. Therefore, any shifted highest energy state (7.5) would satisfy (8.1). A basis for $V^{(1)}_{NM}$ can thus be constructed by considering the states (7.5) in the blocks $(k,M-k)$, applying the shift operator, and then populating mode 0 by acting with $a_{0}^{\dagger}$ repeatedly to reach the target $(N,M)$-block:

$$|\phi^{(1,k)}_{N,M}\rangle\equiv\frac{a_{0}^{\dagger(N-k)}}{\sqrt{(N-k)!}}S|\psi^{\max}_{k,M-k}\rangle.$$ (8.2)

Note that $k$ ranges from 1 to $\min(N,M)$, as expected by the above counting argument. Within each block, the states (8.2) are linearly independent as they contain a different number of modes $a^{\dagger}_{0}$. In conclusion, the states (8.2) provide us with a complete basis for the top eigenspace of $H_{\min}$.

The next step is to compute the matrix elements of $H$ in the basis (8.2) and diagonalize the resulting matrix. Fortunately, since a quartic resonant Hamiltonian can at most create or destroy one zero mode, $H$ is tridiagonal in the basis (8.2). Therefore, it is enough to focus on cases $k^{\prime}=k$ and $k^{\prime}=k+1$, in combination with the hermiticity of $H$. It will be useful to decompose $H$ according to the number of appearances of the zero mode. Let us first focus on the case $k^{\prime}=k$. The only terms in $H$ that contribute in that case have an identical number of $a^{\dagger}_{0}$ and $a_{0}$. These terms can be conveniently written as

$$H\to\frac{1}{2}\sum_{n+m=k+l}^{n,m,k,l\geq 0}a^{\dagger}_{n+1}a^{\dagger}_{m+1}a_{k+1}a_{l+1}+2a^{\dagger}_{0}a_{0}\hat{N}-\frac{3}{2}\left(a^{\dagger}_{0}a_{0}\right)^{2}-\frac{1}{2}a^{\dagger}_{0}a_{0}\,,$$ (8.3)

such that

$$\bra{\phi^{(1,k)}_{N,M}}H\ket{\phi^{(1,k)}_{N,M}} =\bra{\psi^{\max}_{k,M-k}}H\ket{\psi^{\max}_{k,M-k}}+2(N-k)N-\frac{3}{2}\left(N-k\right)^{2}-\frac{1}{2}(N-k)$$

$$=\left(\frac{N(N-1)}{2}-M\right)+\left(N+M+1\right)k-2k^{2}.$$ (8.4)

We now turn to the case $k^{\prime}=k+1$, where the terms in $H$ need to have a single $a_{0}$ present, and no other creation-annihilation operator with index 0, to give a nonzero contribution:

$$\bra{\phi^{(1,k+1)}_{N,M}}H\ket{\phi^{(1,k)}_{N,M}} =\frac{\bra{\psi^{\max}_{k+1,M-k-1}}S^{\dagger}a_{0}^{(N-k-1)}Ha_{0}^{\dagger(N-k)}S\ket{\psi^{\max}_{k,M-k}}}{\sqrt{(N-k-1)!(N-k)!}}\,,$$

$$=\frac{\bra{\psi^{\max}_{k+1,M-k-1}}S^{\dagger}a_{0}^{(N-k)}\left(\sum_{n,m=1}^{\infty}a^{\dagger}_{n}a^{\dagger}_{m}a_{n+m}\right)a_{0}^{\dagger(N-k)}S\ket{\psi^{\max}_{k,M-k}}}{\sqrt{(N-k-1)!(N-k)!}}\,,$$

$$=\sqrt{N-k}\bra{\psi^{\max}_{k+1,M-k-1}}K_{+}\ket{\psi^{\max}_{k,M-k}}\,,$$ (8.5)

with

$$K_{+}=\sum_{n,m=0}^{\infty}a^{\dagger}_{n}a^{\dagger}_{m}a_{n+m+1}\,.$$ (8.6)

Using the constraints (7.1) on the bra-state, the action of this operator $K_{+}$ on the highest-energy states can be simplified to

$$\bra{\phi^{(1,k+1)}_{N,M}}H\ket{\phi^{(1,k)}_{N,M}} =\sqrt{N-k}\bra{\psi^{\max}_{k+1,M-k-1}}a^{\dagger}_{0}\sum_{j=0}^{\infty}(j+1)a^{\dagger}_{j}a_{j+1}\ket{\psi^{\max}_{k,M-k}}\,,$$

$$=\sqrt{N-k}\bra{\psi^{\max}_{k+1,M-k-1}}a^{\dagger}_{0}Z^{\dagger}\ket{\psi^{\max}_{k,M-k}}\,.$$ (8.7)

It is simplest to consider the action of this product of operators on the bra-state. Since $a_{0}$ commutes with the constraints (7.1), this operator transports the highest energy state from the block $(N,M)$ to the highest energy state in block $(N-1,M)$. Moreover, $a_{0}$ and $Z$ commute so that the action of $a_{0}$ in (8.7) can be straightforwardly computed using the expression (7.5)

$$a_{0}\ket{\psi^{\max}_{N,M}}=\sqrt{\frac{N(N-1)}{N+M-1}}\ket{\psi^{\max}_{N-1,M}}\,.$$ (8.8)

In addition, we established above that $Z$ maps the highest energy state (7.5) in block $(N,M)$ to the highest energy state in block $(N,M+1)$. Using (7.5), one can compute the proportionality factor as

$$Z\ket{\psi^{\max}_{N,M}}=\sqrt{(M+1)(N+M)}\ket{\psi^{\max}_{N,M+1}}.$$ (8.9)

Substituting (8.9) and (8.8) in (8.7), the off-diagonal matrix elements become

$$\bra{\phi^{(1,k+1)}_{N,M}}H\ket{\phi^{(1,k)}_{N,M}}=\sqrt{(N-k)(M-k)(k+1)k}=\bra{\phi^{(1,k)}_{N,M}}H\ket{\phi^{(1,k+1)}_{N,M}}.$$ (8.10)

This completes the computation of the matrix elements of $H$ in the top eigenbasis of $H_{\min}$.

Tridiagonal matrices analogous to the ones we have arrived at, with nonzero matrix elements given by (8.4) and (8.10), have been diagonalized in [47]. Remarkably, their spectrum follows a square integer sequence. This property therefore also defines the eigenvalues of $H$ in the top $H_{\min}$ eigenspace. In fact, all the eigenvalues of $H$ that we will be able to bring under some degree of analytic control are parts of square integer sequences of this sort. We now repeat the argument of [47], adapted to our needs, for completeness.

The claim is that the eigenvalues of the tridiagonal matrix given by (8.4) and (8.10) are

$$E^{(1,k)}=\left(\frac{N+M+2-2k}{2}\right)^{2}+c^{(1)}\quad\text{for}\quad k=1,\dots,\min(N,M),$$ (8.11)

with $c^{(1)}\equiv N(N-1)/2-M-(N-M)^{2}/4$ and the superscript refers to the top eigenspace of $H_{\min}$. In the following, we shall focus on deriving the square integer part of (8.11) and therefore subtract the constant part $c^{(1)}$ from the diagonal elements (8.4) in our treatment below.

Before delving into the proof, it is convenient to rewrite the matrix elements in terms of the dimensionality of the basis $n\equiv\min(N,M)$ and the parameter $a\equiv|N-M|$, instead of $N$ and $M$. The following relations hold:

$$\max(N,M)=a+n\quad\text{and}\quad N+M=\min(N,M)+\max(N,M)=2n+a.$$ (8.12)

Note that either $N=n$ and $M=a+n$ or $M=n$ and $N=a+n$. The diagonal elements (8.4), with the constant $c^{(1)}$ subtracted, can then be written as

$$a_{k}=-2k^{2}+(2n+a+1)k+a^{2}/4\quad\text{for}\quad k=1,\cdots,n\,,$$ (8.13)

while the off-diagonal terms become

$$b_{k}=\sqrt{(n-k)(n+a-k)(k+1)k}\quad\text{for}\quad k=1,\cdots,n-1.$$ (8.14)

Together, they make up a matrix that we will refer to as $H_{n}$. We now want to show that the eigenvalues of this tridiagonal matrix $H_{n}$ are

$$\tilde{E}_{k}=\left(\frac{a}{2}+k\right)^{2}\quad\text{for}\quad k=1,\dots,n.$$ (8.15)

The proof works by induction [47]. The statement is trivially true at $n=1$, since in that case $a_{1}=\tilde{E}_{1}$. We aim to show that the statement is true at $n$ assuming it holds at $n-1$, keeping $a$ fixed. At level $n$ we consider the decomposition $H_{n}\equiv\left(n+{a}/{2}\right)^{2}I-C_{n}$ and notice that the matrix $C_{n}$ can be written as the product $C_{n}\equiv A_{n}A_{n}^{T}$ with $A_{n}$ a lower bidiagonal matrix. Its diagonal elements $h_{k}$ and lower subdiagonal elements $r_{k}$ are respectively given by

$$h_{k}=\sqrt{(n-k)(n+a-k)} \quad\text{for}\quad k=1,\cdots,n,$$ (8.16)

$$r_{k}=-\sqrt{(k+1)k}\quad\text{for}\quad k=1,\cdots,n-1.$$ (8.17)

The key observation is that $D_{n}\equiv\left(n+{a}/{2}\right)^{2}I-A_{n}^{T}A_{n}$, where the matrix $A_{n}$ and $A_{n}^{T}$ have swapped their positions, decomposes as

$$D_{n}=\begin{pmatrix}&&&0\\ &H_{n-1}&&\vdots\\ &&&0\\ 0&\cdots&0&\left(n+{a}/{2}\right)^{2}\end{pmatrix}.$$ (8.18)

By use of the singular value decomposition of $A_{n}$, it is straightforward to show that $A_{n}^{T}A_{n}$ and $A_{n}A_{n}^{T}$ share the same nonzero eigenvalues. We therefore find that $H_{n-1}$ and $H_{n}$ display the same $n-1$ lowest eigenvalues (at fixed value of $a$) and the additional eigenvalue for $H_{n}$ is $\tilde{E}_{n}$. This completes the proof. Note that the above construction with the swapped order of $A$ and $A^{T}$ in the partner matrices $C$ and $D$ is strongly reminiscent of the factorization method in quantum mechanics [46], with its connections to supersymmetric quantum mechanics.

The eigenvectors of $H_{n}$ can be derived in an analogous manner. We first focus on the highest energy vector, with energy $\tilde{E}_{n}$. This corresponds to the null direction of $C_{n}$, which satisfies $A_{n}^{T}v=0$. This condition is solved by

$$v_{k+1}=-\frac{h_{k}}{r_{k}}v_{k}\quad\text{for}\quad k=1,\cdots,n-1\,,$$ (8.19)

where we used that $h_{n}=0$. This uniquely specifies the highest energy eigenvector of $H_{n}$. The normalized components of this vector in the basis (8.2) are

$$v^{(1,1)}_{k}=\sqrt{\frac{n!(a+n)!}{(2n+a-1)!}}\sqrt{\frac{1}{k}\binom{n-1}{k-1}\binom{n+a-1}{k-1}}=\sqrt{\frac{M!N!}{(N+M-1)!}}\sqrt{\frac{1}{k}\binom{N-1}{k-1}\binom{M-1}{k-1}}\,.$$

The resulting energy eigenvector is precisely the highest energy state $\ket{{\psi}^{\max}_{N,M}}$ derived in the previous section.

In order to find the other eigenvectors, we use the product decomposition (8.18). We first focus on the $(n-1)$-dimensional space in which $H_{n-1}$ is defined. Following identical steps as above, we find the highest energy state $v^{(n-1)}$ of this matrix using (8.19) and replacing $n\rightarrow n-1$. This state, embedded into an $n$-dimensional space by appending a 0 to its column vector representation, is also an eigenstate of $D_{n}$. It is straightforward to show that, if $v^{(n-1)}$ is an eigenvector of $A_{n}^{T}A_{n}$, then $A_{n}v^{(n-1)}$ is an eigenvector of $A_{n}A_{n}^{T}$ with the same eigenvalue. The resulting vector $A_{n}v^{(n-1)}$ is an eigenvector of $H_{n}$ with eigenvalue $\tilde{E}_{n-1}$, as required. This process can be repeated for any of the eigenvalues $\tilde{E}_{k}$.

The result of this section is complete analytic understanding of the eigenvalues and eigenvectors of $H$ in the top eigenspace of $H_{\min}$. (For later reference, we shall denote the normalized eigenvectors associated with the energies (8.11) as $\ket{\psi^{(1,k)}_{N,M}}$.) Under the shift operator, these joint eigenstates are mapped to eigenstates of $H_{\min}$ with distinct eigenvalues. These eigenvalues can be found using (4.5). Applying (4.5) to the eigenstates with energies (8.11) considered in block $(N,M-N$) yields the following sequence of $H_{\min}$ eigenvalues in block $(N,M)$

$$E_{\min}=\left(M-k+1\right)^{2}+\left(k-1\right)^{2}\,,$$ (8.20)

with $k=1,\dots,\min(N,M-N)$. These $H_{\min}$ eigenvalues are preserved when transporting the shifted states to higher $N$ blocks using $a^{\dagger}_{0}$. For each $k$, one could try to construct a subspace in the block $(N,M)$ in which to diagonalize $H$ following the construction (4.6). Of course, we are not guaranteed to already know all the states in lower blocks that would map to the eigenvalues (8.20) under the shift $S$, and $H$ might therefore not leave the constructed subspaces invariant. Suppose we could complete the states $\ket{\psi^{(1,k)}_{N,M}}$ to a joint eigenbasis of $H$ and $H_{\min}$ in lower blocks and construct (4.6). Following similar steps as in the previous section, we would find that the Hamiltonian $H$ has nonzero matrix elements in two cases only. First, the matrix elements between two same states can be computed using (8.3) and is generally nonzero. For a pair of distinct states with the same number of zero modes turned on in (4.6), the respective contribution vanishes by orthogonality of states in the spaces $V_{NM}$. Second, pairs of states (4.6) that differ by one mode $a^{\dagger}_{0}$ are related by $H$ if and only if the states from which they were derived using the shift operator have a nonzero overlap through the operator $K_{+}$, as in (8.5). All other matrix elements vanish.

9 Ladder operators

We have observed numerically that the square integer structure displayed by the top eigenspace of $H_{\min}$ is replicated in other eigenspaces. A key component in the diagonalization process of $H$ is the behavior of $K_{+}$ defined by (8.6) on joint eigenstates of $H$ and $H_{\min}$, as we have just discussed. In this section, we show that $K_{+}$ has in fact an even more important role in the structure of eigenstates and eigenvalues: it is a ladder operator for the combination $2H+H_{\min}$.

Introducing the Hermitian conjugate of $K_{-}\equiv K_{+}^{\dagger}$, the following commutation relation holds:

$$[2H+H_{\min},K_{\pm}]=\pm K_{\pm}.$$ (9.1)

Hence, $K_{+}$ ($K_{-}$) raises (lowers) the eigenvalues of $2H+H_{\min}$ by one unit. Note that the part of this commutation relation quintic in creation-annihilation operators cancels out between the contributions from $H$ and $H_{\min}$. The remaining cubic part is precisely the right-hand side of (9.1). This automatically implies that the classical quantities corresponding to $K_{\pm}$ are conservation laws for the dynamics induced by $2H+H_{\min}$, while their function as ladder operators in the quantum theory arises purely from quantum corrections to the classical Poisson brackets and the quantum correction in $H_{\min}$.

To verify (9.1), we will prove an equivalent statement:

$$[H_{\min},\tilde{J}_{3}]=0,\qquad\tilde{J}_{3}\equiv\sum_{n,m=1}^{\infty}a^{\dagger}_{m+n}a_{m}a_{n}.$$ (9.2)

This commutation relation follows straightforwardly by first observing that

$$J_{3}\equiv\sum_{n,m=0}^{\infty}a^{\dagger}_{m+n}a_{m}a_{n}$$

commutes with $H_{\min}$. From (4.9), $J_{3}$ is defined by the commutator of $a_{0}$ and $H$, both of which commute with $H_{\min}$. The difference between $\tilde{J}_{3}$ and $J_{3}$ involves $Na_{0}$ and a cubic operator made of $a_{0}$ and $a_{0}^{\dagger}$, all of which commute with $H_{\min}$. This proves (9.2). To relate (9.2) to (9.1), we notice that shifting the mode numbers up by one unit $a_{n}\rightarrow a_{n+1}$ (and the same for $a^{\dagger}$) sends $K_{-}$ to $\tilde{J}_{3}$ and $2H+H_{\min}$ to $H_{\min}$, up to extra terms involving $N$ and $M$ that are described by (4.5). One can straightforwardly commute these extra terms with $K_{-}$, producing the nonzero right-hand-side in (9.1).

This ladder operator has a complementary action to the operator $a^{\dagger}_{0}$ which acts as a ladder operator for $H_{\min}$ by copying its eigenvalues from block $(N,M)$ to $(N+1,M)$, as well as the shift operator $S$ whose action on the eigenvectors is determined by (4.5). Of course, all the operators obtained by commuting either $K_{\pm}$ or $a^{\dagger}_{0}$ with any of the conservation laws of the model trivially define additional (higher polynomial) ladder operators.

10 Eigenstates within lower eigenspaces of $H_{\min}$

Having understood the structure of the top eigenspace of $H_{\min}$, we would like to follow similar steps for the $H_{\min}$ eigenspace associated to the next eigenvalue in the sequence (8.20),

$$E^{(N,M)}_{min}=(M-1)^{2}+1.$$ (10.1)

To this end, we consider the states

$$\ket{\phi^{(2,k)}_{N,M}}\equiv\frac{a_{0}^{\dagger(N-k-1)}}{\sqrt{(N\!-\!k\!-\!1)!}}S\ket{\psi^{(1,2)}_{k+1,M-k-1}}$$ (10.2)

where $\ket{\psi^{(1,2)}_{N,M}}$ is the next-to-highest energy eigenvectors of $H$ (constructed as $A_{n}v^{(n-1)}$ in the previous section) in the top eigenspace of $H_{\min}$ in block $(N,M)$. The states (10.2) exist for $1\leq k\leq\min(N-1,M-3)$. Equation (8.20) shows that the states (10.2) all share the same $H_{\min}$ eigenvalue equal to (10.1) and we aim to diagonalize the Hamiltonian in this subspace. However, in order to perform this diagonalization, a complete basis of the eigenspace at energy (10.1) is needed. In contrast to section 8, we do not have a counting argument showing that the states (10.2) form a complete basis for the eigenspace of $H_{\min}$ at value (10.1). However, as we shall see below, the Hamiltonian $H$ indeed leaves the subspace spanned by the basis (10.2) invariant due to the simple action of $K_{+}$ on states $\ket{\psi^{(1,2)}_{N,M}}$, see (10.7).

Evaluating $H$ in the basis (10.2) once more yields a tridiagonal matrix, with eigenvalues following the square integer sequence:

$$E^{(2,k)}=\left(\frac{N+M-2-2k}{2}\right)^{2}+\frac{N(N+1)}{2}-M-\frac{(M-N-2)^{2}}{4},$$ (10.3)

with $k=1,\cdots,\min(N-1,M-3)$. To show this, we construct the matrix elements of $H$ in (10.2). The computation of the diagonal elements is identical to section 8. It now takes as input the energy of the states $|\psi^{(1,2)}_{k+1,M-k-1}\rangle$,

$$E^{(1,2)}_{k}=(k-1)\left(M-1-\frac{k}{2}\right),$$ (10.4)

so that

$$\bra{\phi^{(2,k)}_{N,M}}H\ket{\phi^{(2,k)}_{N,M}}=-2k^{2}+(N+M-3)k+\frac{N(N+1)}{2}-M.$$ (10.5)

For the off-diagonal term with $k^{\prime}=k+1$, we write

$$\bra{\phi^{(2,k+1)}_{N,M}}H\ket{\phi^{(2,k)}_{N,M}}=\sqrt{N-k-1}\bra{\psi^{(1,2)}_{k+2,M-k-2}}K_{+}\ket{\psi^{(1,2)}_{k+1,M-k-1}}=\bra{\phi^{(2,k)}_{N,M}}H\ket{\phi^{(2,k+1)}_{N,M}}\,,$$ (10.6)

The operator $K_{+}$ can be shown to have a simple action on the states $\ket{\psi^{(1,2)}_{n,m}}$:

$$K_{+}\ket{\psi^{(1,2)}_{N,M}}=\sqrt{N(N-1)(M-2)}\ket{\psi^{(1,2)}_{N+1,M-1}}\,.$$ (10.7)

The steps to show (10.7) are straightforward but tedious, and we summarize this derivation in Appendix C. In light of the previous section, the map (10.7) is only possible by the subtle relation between the eigenvalue of the pair of states with respect to $2H+H_{\min}$, which have to differ by exactly one unit to be consistent with (9.1).

We can then write (10.6) as

$$\bra{\phi^{(2,k+1)}_{N,M}}H\ket{\phi^{(2,k)}_{N,M}}=\sqrt{(N-k-1)(M-k-3)(k+1)k}.$$ (10.8)

with $k=1,\cdots,\min(M-3,N-1)-1$. The resulting matrix elements (10.5) and (10.8) have a form identical to (8.4) and (8.10). This matrix can diagonalized in an analogous manner to section 8, now with $n=\min(N-1,M-3)$ and $a=|M-N-2|$. All subsequent steps are identical and one finds the square integer spectrum (10.3), with eigenvectors, for the operator $H$ in the $H_{\min}$ eigenspace at level (10.1).

We can now shift the joint (subspace) eigenbasis $\ket{\psi^{(2,k)}_{N,M-N}}$ of $H$ and $H_{\min}$ using $S$ and find that the resulting $H_{\min}$ eigenvalues in block $(N,M)$ belong to the sequence

$$E_{\min}=(M-k-1)^{2}+(k+1)^{2},$$ (10.9)

with $k=1,\dots,\min(N-1,M-N-3)$. Note that there are generally at least two different types of states that are mapped to the $E^{(N,M)}_{\min}=(M-i)^{2}+i^{2}$ eigenspace (with $i>1$) under the shift: the states $\ket{\psi^{(1,1+i)}_{N,M-N}}$ and the states $\ket{\psi^{(2,i-1)}_{N,M-N}}$. This means, using (4.5), that these states share the same eigenvalue under $2H+H_{\min}$. This level is therefore degenerate, and $K_{+}$ will generically not be diagonal with respect to these two types of joint eigenvectors of $H$ and $H_{\min}$. This fact represents an obstacle for straightforwardly applying the above method for diagonalizing $H$ in further lower eigenspaces of $H_{\min}$. Indeed, the Hamiltonian $H$ in a basis composed of shifted joint eigenvectors will no longer be tridiagonal. Nonetheless, we have empirically found that the square integer sequences continue to emerge, suggesting that there exists a preferred basis for $K_{+}$ in which the $H$ eigenblocks acquire the tridiagonal structure of the sort encountered above.

We conclude this section with closed-form expressions for some of the numerically observed sequences of eigenvalues of $H$. At level $E_{\min}=\left(M-i\right)^{2}+i^{2}$, with $i>0$ there are at least $i$ distinct square integer sequences parameterized by $n=0,\dots,i-1$ that follow the pattern

$$E^{i,n}_{k}= \left(\frac{N+M-3i+n+1-2k}{2}\right)^{2}+\frac{(N+n)(N+n+1)}{2}$$

$$\hskip 71.13188pt-M+\frac{(1+n)(n-2i+2)}{2}-\left(\frac{M-N-n-i-1}{2}\right)^{2}\,,$$ (10.10)

for $k=1,\dots,\min(N-i+n,M-2i-1)$. This formula reproduces in particular the energies (10.3) when $i=1$. In addition, at level $E_{\min}=(M-3-i)^{2}+(i+3)^{2}-2(i+2)$ for $i>0$ there are at least $i$ distinct square integer sequences parameterized by $n=0,\dots,i-1$ that are distributed according to the sequences:

$$E^{i,n}_{k}= \left(\frac{N+M-3i+n-5-2k}{2}\right)^{2}+\frac{(N+n+1)(N+n+2)}{2}$$

$$\hskip 71.13188pt-M+\frac{(n+1)(n-2i-2)}{2}-\left(\frac{M-N-n-i-5}{2}\right)^{2}\,,$$ (10.11)

for $k=1,\dots,\min(N-i+n-1,M-2i-6)$.

11 Relation to classical invariant manifolds

The families of eigenstates derived above are reminiscent of the integrable sectors in the resonant systems studied in [39, 40]. These system were previously shown to have a single (three-dimensional) invariant manifold [15, 27, 30] and the exactly solvable set of Hamiltonian eigenstates can be combined into coherent-like superpositions that reproduce the classical dynamics in the semiclassical limit [39, 40]. This connection between special subsets of Hamiltonian eigenstates and classical dynamically invariant manifolds can furthermore be effectively visualized using phase space techniques [48]. There is a considerable resemblance between this story outlined in [39, 40] and the special families of eigenvectors of the GG Hamiltonian constructed above, hence one may expect the eigenstates of $H$ to include sectors that encode the classical dynamics within the invariant manifolds (3.5).

The simplest correspondence of this sort is the relation between the top eigenstates of $H$ within each $(N,M)$-block and the ansatz (3.5) with $R=1$. To see this, we examine the expectation values of the mode occupation number $a^{\dagger}_{n}a_{n}$ within these states, which is an analog of the classical observable $|\alpha_{n}|^{2}$. For the highest energy state (7.5), one can show that

$$a_{n}|\psi^{\max}_{N,M}\rangle=\sqrt{N(N-1)}\,\,\sqrt{\frac{M!(N+M-n-2)!}{(M-n)!(N+M-1)!}}\,\,|\psi^{\max}_{N-1,M-n}\rangle\,.$$ (11.1)

In the semiclassical limit, where $N,M\gg n$ are large while keeping their ratio $d^{2}\equiv M/N$ fixed, one then gets

$$\langle\psi^{\max}_{N,M}|a^{\dagger}_{n}a_{n}|\psi^{\max}_{N,M}\rangle=\frac{N}{1+d^{2}}\left(\frac{d^{2}}{1+d^{2}}\right)^{n}.$$ (11.2)

This has exactly the same form as the observable $|\alpha_{n}|^{2}$ for $\alpha_{n}$ within the ansatz (3.5) at $R=1$. The energy of the solutions in this lowest manifold, denoted by $\mathcal{M}(1)$ in [1], likewise saturates the classical energy bound [1].

We expect this connection with classical invariant manifolds to extend to more complicated families of eigenvectors, though we will not explore this in detail here. To clarify our intuition, we recall that, according to [1], in addition to the invariant manifolds (3.5), denoted there $\mathcal{M}(R)$, one has the invariant manifolds

$$\alpha_{n}(t)=c(t)\,\delta_{n,0}+\sum_{r=1}^{R}c_{r}(t)\,[p_{r}(t)]^{n},$$ (11.3)

denoted $\tilde{\mathcal{M}}(R)$. There is an interlacing embedding of these manifolds into each other. The manifold $\mathcal{M}(R)$ is evidently a submanifold of $\tilde{\mathcal{M}}(R)$ obtained by setting $c=0$. On the other hand, the manifold $\tilde{\mathcal{M}}(R)$ is a submanifold of $\mathcal{M}(R+1)$ obtained by setting $p_{R+1}=0$, $c_{R+1}=c$. The action of the classical shift operator (3.21) moves configurations up this invariant manifold ladder. Namely, if we apply the shift to a configuration (3.5), the result still satisfies (3.5), with redefined $c_{r}$ and $p_{r}$, except for $n=0$, since now $\alpha_{0}=0$. In other words, we obtain a configuration that fits the ansatz (3.5), except for the shifted $\alpha_{0}$. But this is precisely the definition (11.3), so the configuration is in $\tilde{\mathcal{M}}(R)$. Furthermore, since $\tilde{\mathcal{M}}(R)\in\mathcal{M}(R+1)$, subsequent action of the shift operator will push the configuration to $\tilde{\mathcal{M}}(R+1)$ and so on.

As explained above, the top eigenstates correspond to configurations in $\mathcal{M}(1)$ — more precisely, one will have to take superpositions of the top eigenstates in different nearby blocks to obtain a state localized around a specific classical configuration $\alpha_{n}=cp^{n}$. We have furthermore been repeatedly obtaining other families of eigenvectors by acting on these top vectors (and their descendants) with the quantum shift operator (4.1). This strongly suggests that these lower families correspond to higher invariant manifolds described above. One would need, however, a more thorough and systematic understanding of the descendant eigenstate families to explore this correspondence comprehensively.

12 Quantum Lax pair and conservation laws

Quantum Lax pairs [49] are not a very common subject for discussion, since in general, they are both difficult to construct because of operator ordering problems and do not provide the same empowerment of the formalism and automatic construction of conservation laws as in the corresponding classical theory. They have found, however, some applications to quantum Calogero systems [50], see [51, 52].

For the quantum GG Hamiltonian, we have been able to identify a curious Lax pair structure that we shall describe here. As in the classical case, introduce an auxiliary space $\vec{h}\equiv(h_{0},h_{1},\ldots)$ of number-valued vectors in the mode space. The Lax operators will be infinite matrices acting on these vectors, while the entries of these matrices will be Hilbert space operators made of $a$ and $a^{\dagger}$. The specific expressions are

$$(\mathcal{L}\vec{h})_{n}=\sum_{k,l=0}^{\infty}a_{l+k}^{\dagger}a_{k+n}h_{l}+nh_{n},\qquad(\mathcal{M}\vec{h})_{n}=\sum_{l+m\geq n}^{l,m\geq 0}a_{l+m-n}^{\dagger}a_{l}h_{m}.$$ (12.1)

Except for the second ‘quantum’ term in $\mathcal{L}$, this replicates the structure of the classical Lax pair (3.7). Introducing the identity matrix $I$ in the space of $\vec{h}$, so that $I\vec{h}=\vec{h}$, the evolution of $\mathcal{L}$ is

$$i\dot{\mathcal{L}}=[\mathcal{L},HI]\,,$$ (12.2)

and the Lax pair condition can be written simply as

$$[HI+\mathcal{M},\mathcal{L}]=0.$$ (12.3)

To prove by brute force that (12.3) is valid, we can evaluate the commutators directly and then rewrite the result in terms of normal-ordered operators. The top-order piece, which is quartic in the creation-annihilation operators has to vanish by the classical computation in Appendix A, since all the terms are normal-ordered, no commutations of $a$ and $a^{\dagger}$ are necessary, and the classical and quantum computations will agree line-by-line. The remainder consists of two quadratic contributions from

$$([\mathcal{M},\mathcal{L}]h)_{n}= \sum_{l+m\geq n}^{k,k^{\prime},l,m\geq 0}a_{l+m-n}^{\dagger}a_{l}a_{k^{\prime}+k}^{\dagger}a_{k+m}h_{k^{\prime}}+\sum_{l+m\geq n}^{l,m\geq 0}m\,a_{l+m-n}^{\dagger}a_{l}h_{m}$$

$$-\sum_{k^{\prime}+m\geq l}^{k,k^{\prime},m,l\geq 0}a_{l+k}^{\dagger}a_{k+n}a_{k^{\prime}+m-l}^{\dagger}a_{k^{\prime}}h_{m}-n\sum_{l+m\geq n}^{l,m\geq 0}a_{l+m-n}^{\dagger}a_{l}h_{m}.$$ (12.4)

First, there are terms that emerge from bringing the quartic contributions to $[\mathcal{M},\mathcal{L}]$ into a normal-ordered form; second, there are terms arising from the quantum piece $nh_{n}$ in the $\mathcal{L}$-operator of (12.1). We must show that these two quadratic pieces cancel each other. We consider the quartic contributions from the first and third terms in (12.4) in turn. The first term gives the following quadratic contribution:

$$\sum_{l+m\geq n}^{k,k^{\prime},l,m\geq 0}a_{l+m-n}^{\dagger}a_{l}a_{k^{\prime}+k}^{\dagger}a_{k+m}h_{k^{\prime}} \to\sum_{l+m\geq n}^{k,k^{\prime},l,m\geq 0}a_{l+m-n}^{\dagger}\delta_{l,k^{\prime}+k}a_{k+m}h_{k^{\prime}}=\sum_{k+k^{\prime}+m\geq n}^{k,k^{\prime},m\geq 0}a_{k+k^{\prime}+m-n}^{\dagger}a_{k+m}h_{k^{\prime}}$$

$$=\sum_{\begin{subarray}{c}k^{\prime}+j\geq n\\ j\geq m\end{subarray}}^{j,k^{\prime},m\geq 0}a_{k^{\prime}+j-n}^{\dagger}a_{j}h_{k^{\prime}}=\sum_{k^{\prime}+j\geq n}^{j,k^{\prime}\geq 0}j\,a_{k^{\prime}+j-n}^{\dagger}a_{j}h_{k^{\prime}}.$$ (12.5)

First, we remove the sum over $l$ using the Kronecker $\delta$, then we introduce the index $j=k+m$ to swap the sum over $k$ for a sum over $j$. The label $m$ no longer appears in the summand, and therefore the respective sum is straightforwardly computed. We can perform similar steps for the third term in (12.4):

$$\sum_{k^{\prime}+m\geq l}^{k,k^{\prime},m,l\geq 0}a_{l+k}^{\dagger}a_{k+n}a_{k^{\prime}+m-l}^{\dagger}a_{k^{\prime}}h_{m} \to\sum_{k^{\prime}+m\geq l}^{k,k^{\prime},m,l\geq 0}a_{l+k}^{\dagger}\delta_{k+n,k^{\prime}+m-l}a_{k^{\prime}}h_{m}=\sum_{\begin{subarray}{c}k^{\prime}+m\geq n\\ k^{\prime}+m\geq n+k\end{subarray}}^{k,k^{\prime},m\geq 0}a_{k^{\prime}+m-n}^{\dagger}a_{k^{\prime}}h_{m}$$

$$=\sum_{k^{\prime}+m\geq n}^{k^{\prime},m\geq 0}(k^{\prime}+m-n)a_{k^{\prime}+m-n}^{\dagger}a_{k^{\prime}}h_{m},$$ (12.6)

where we performed the sum over $l$ and this eliminated all dependence on the label $k$ in the summand. Subtracting (12.6) from (12.5), the remaining terms precisely cancel the second and fourth contributions in (12.4). We have therefore explicitly shown that (12.3) holds at the quantum level, with the modified Lax pair (12.1).

It is usually impossible to directly construct quantum conservation laws from such quantum Lax pairs as $\mathrm{Tr}\mathcal{L}^{n}$, since the lack of commutativity of Hilbert space operators upsets the cyclic property of the trace [49]. However, in our case, there is a trick that closely replicates the construction of classical conservation laws in (3.14) originally discovered in [1], which is different from the usual trace-based construction. Namely, we first introduce the vector $\vec{1}\equiv(1,0,0,0,\ldots)$ and notice the identity

$$\mathcal{L}\vec{1}=\mathcal{M}\vec{1}.$$ (12.7)

Then we introduce the operators $(\vec{1},\mathcal{L}^{n}\vec{1})$ and observe that they commute with the Hamiltonian because

$$[H,(\vec{1},\mathcal{L}^{n}\vec{1})]=(\vec{1},[HI,\mathcal{L}^{n}]\vec{1})=(\vec{1},(\mathcal{L}^{n}\mathcal{M}-\mathcal{M}\mathcal{L}^{n})\vec{1})=0.$$ (12.8)

As the operators $(\vec{1},\mathcal{L}^{n}\vec{1})$ formally commute with the Hamiltonian, they provide an infinite family of conservation laws. There is unfortunately a complication with this construction, however, that precludes its immediate utilization for practical purposes. The operators obtained in this manner are not normal-ordered, and attempting to normal-order them, one gets contributions in the form of lower order polynomials in $a$ and $a^{\dagger}$ with infinite coefficients. It is natural to expect that these infinities are themselves lower-order conservation laws and may be systematically subtracted, yielding finite, well-defined conserved operators. Namely, one can truncate all sums at a finite large mode number and add a combination of lower-order conservation laws with divergent coefficients depending on this cutoff. With suitable finetuning, the total expression will have a well-defined limit as the cutoff is removed, providing a conservation law. We have verified using computer algebra that this can be done for the first few members of the tower. (We found the FORM software [53] useful for this purpose, due to its built-in ability to effectively handle large polynomial expressions, possibly made of noncommuting variables.) We do not have, however, a systematic theory of such subtractions, and hence no way to employ the Lax pair (12.1) to construct conservation laws effectively.

In the course of our numerical experimentation, we have come across another construction of conservation laws that appears to work at all orders. The formula is very simple and gives the conserved quantities as

$$J_{2n}=[J_{2n+1},a_{0}^{\dagger}],$$ (12.9)

where the sequence $J_{2n+1}$ is defined recursively by (4.9). Note that $J_{2}$ is related to $N$, while $J_{4}$ is related to $H$. The higher $J$’s are independent polynomial conservation laws. To illustrate the sort of quantum corrections that are generated by (12.9), we provide an explicit expression for $J_{6}$:

$$J_{6} =2\hskip-8.5359pt\sum_{k,l,m,n,p=0}^{\infty}\hskip-8.5359pta^{\dagger}_{k}a^{\dagger}_{l+m}a^{\dagger}_{n+p}a_{k+l}a_{m+n}a_{p}+2\hskip-2.84544pt\sum_{k,l,m=0}^{\infty}(k+l+m)\,a^{\dagger}_{k}a^{\dagger}_{l+m}a_{k+l}a_{m}$$ (12.10)

$$\hskip 113.81102pt+2H_{\min}+(2N+6)H+N^{3}+N^{2}+2MN+2M.$$

One can, of course, throw away the second line since it is expressed through the lower conservation laws. The sextic piece in the first line matches the structure of the classical conserved quantities (3.14), but there is a nontrivial quartic correction.

We have not been able to prove analytically that (12.9) are conserved, though we have run extensive checks that the corresponding matrices commute with the Hamiltonian within individual $(N,M)$-blocks. It is not difficult to understand, however, why the corresponding classical quantities are conserved as they are related to the classical tower (3.14). This follows immediately from the expression for the classical quantities $J^{cl}_{2n+1}\equiv i\left\{J^{cl}_{2n-1},H_{cl}\right\}$ corresponding to the operators (4.9), with $J^{cl}_{1}\equiv\alpha_{0}$. In terms of the classical modes, they are

$$J^{cl}_{2n+1}=\hskip-8.53581pt\sum_{i_{1},\ldots,i_{2n}=0}^{\infty}\hskip-8.53581pt\alpha_{i_{1}}\bar{\alpha}_{i_{1}+i_{2}}\alpha_{i_{2}+i_{3}}\cdots\bar{\alpha}_{i_{2n-1}+i_{2n}}\alpha_{i_{2n}}.$$ (12.11)

Then, one straightforwardly finds the Poisson bracket

$$i\left\{J^{cl}_{2n+1},\bar{\alpha}_{0}\right\}=\sum_{n=0}^{\infty}\frac{\partial J^{cl}_{2n+1}}{\partial\alpha_{n}}\frac{\partial\bar{\alpha}_{0}}{\partial\bar{\alpha}_{n}}=\sum_{i=0}^{n}I_{n-i}I_{i}$$ (12.12)

with $I_{0}\equiv 1$, which is indeed conserved. It is remarkable that applying the correspondence principle to the equivalent representation (12.12) of the tower (3.20) to obtain (12.9) leads to quantum conservation laws without running into operator ordering issues.

Additionally, one expects an extra quantum tower related to the classical traces of powers of the Lax operator given by (3.15). That tower starts with $M$ of (2.7), and we also know its quartic member (1.3) that we discovered empirically in [11] and also provided a brute force proof that it is conserved in section 6. Additionally, we have recovered by trial-and-error as well as guided by the corresponding classical conservation law $G_{3}$ of (3.15) the following sextic conserved operator:

$$\sum_{t,n,l,m,k,p,=0}^{\infty}a^{\dagger}_{t+n}a^{\dagger}_{l+m}a^{\dagger}_{k+p}a_{n+l}a_{m+k}a_{p+t}+6\sum_{t,n,l,m=0}^{\infty}ta^{\dagger}_{t+n}a^{\dagger}_{l+m}a_{n+l}a_{m+t}+\sum_{n=0}^{\infty}n^{3}a^{\dagger}_{n}a_{n}.$$ (12.13)

We do not know how to extend this tower to higher orders, nor is it clear what relation it bears to the quantum Lax pair (12.1).

Finally, perhaps the most remarkable property of all, which we once again only know from numerical observations, is that all eigenvalues of the conservation laws are integers, not only for $H$ and $H_{\min}$, but also for (12.9) and (12.13).

13 Parallels with Calogero and Benjamin-Ono systems

The presence of integer eigenvalues for the GG Hamiltonian and its partner conservation laws alludes to the structures typical of Calogero systems [50, 51, 52], the simplest of which describes $N$ particles on a line with $1/r^{2}$ two-body interactions:

$$H_{Calogero}=\sum_{i=1}^{N}\frac{p_{i}^{2}}{2}+\sum_{i<j}\frac{g}{(x_{i}-x_{j})^{2}}.$$ (13.1)

Even beyond the mere observation of integer eigenvalues, the peculiar construction of conservation laws in (3.14) that does not utilize traces finds its parallels in considerations of Calogero particles [52]. Despite these intriguing similarities, the GG Hamiltonian is in the form of a field theory rather than a quantum-mechanical many-body system, while straightforward second quantization of the Calogero system does not produce anything resembling the GG Hamiltonian [54].

More practical parallels can be seen between the structures revealed by our present exposition and the quantum Benjamin-Ono system [55, 56]. At the classical level, similarity of integrability structures of the Benjamin-Ono and cubic Szegő equations has been revealed in [57, 58] and has led to an explicit analytic solution of the Benjamin-Ono equation analogous to (3.22) – see [59, 60] for further generalizations. The quantum Benjamin-Ono Hamiltonian defined in [55] can be recast in a notation compatible with (1.1) as

$$H_{qBO}=\sum_{m,n=1}^{\infty}\sqrt{\alpha^{3}\,mn(m+n)}\left(a^{\dagger}_{m}a^{\dagger}_{n}a_{m+n}+a^{\dagger}_{m+n}a_{n}a_{m}\right)+\alpha(\alpha-1)\sum_{n=1}^{\infty}n^{2}a^{\dagger}_{n}a_{n},$$ (13.2)

where $\alpha$ is a parameter. (Compared to [55], we have changed the normalization of the creation-annihilation operators in order to give them standard commutation relations.) This Hamiltonian conserves $M$ defined by (2.7), but not the particle number $N$. It possesses a quantum Lax pair [55] leading to a construction of conservation laws reminiscent of (12.8), without encountering the divergences seen in the previous section for the GG Hamiltonian.

A complete solution for the joint eigenstates of the Hamiltonian (13.2) and its hierarchy of conservation laws was given in [55]. The wavefunctions are expressed through Jack’s polynomials in the variables $x_{i}$ defined by decomposing the classical variables $\alpha_{n}$ corresponding to the quantum operators $a_{n}$ as

$$\alpha_{n}=\frac{1}{\sqrt{n\alpha}}\left(x_{1}^{n}+x_{2}^{n}+x_{3}^{n}+\cdots\right),$$ (13.3)

which is reminiscent of the formula (3.5) for the invariant manifolds of the GG Hamiltonian. The eigenstates of the Hamiltonian and other conserved operators are indexed by partitions and the corresponding eigenvalues are expressed by elegant formulas in terms of part lengths [55].

Despite these tantalizing parallels between the structures observed in the GG Hamiltonian and the solution of the quantum Benjamin-Ono Hamiltonian in [55], we have not been able to make immediate use of these parallels. As already mentioned, the construction of conservation laws based on the quantum Lax pair for the GG Hamiltonian runs into divergences when attempting to normal-order the resulting operators. One can furthermore examine the numerically computed integer spectra within specific $(N,M)$-blocks and try to guess the corresponding formula in terms of part lengths. This does not appear possible, however, even for the very simple blocks with $N=2$ that contain only one nonzero eigenvalue [6]. Intuitively, one should keep in mind the mode $a_{0}$ absent in the Benjamin-Ono Hamiltonian but present in the GG Hamiltonian. In the Fock states, the occupation number of this mode $\eta_{0}$ encodes the number of parts $N-\eta_{0}$ in the corresponding partition. Apparently, the eigenvalue formulas must be sensitive to the number of parts, and not only to the part lengths as in [55].

14 Conclusions

We have undertaken studies of integrability properties of the Hamiltonian (1.1) and its partner Hamiltonian (1.3), relying on a mixture of analytic methods and numerical experimentation, and revealing rich and surprising algebraic structures. While we have not been able to solve the model in the sense of giving an explicit construction of all of its eigenfunctions and the corresponding eigenvalues, a number of prominent patterns have been identified. This includes the decompositions (5.15) and (5.17) providing energy bounds, explicit construction of many families of eigenvectors in sections 7, 8 and 10, and alluding to their connections to classical invariant manifolds in section 11, the ladder operator (8.6), the quantum Lax pair (12.1) capable of generating higher conservation law, but requiring subtraction of divergent terms arising in the course of normal ordering, and empirically discovered conservation laws (12.9) and (12.13).

We hope that further studies of this system will supply the missing pieces of the puzzle, and similarities with the quantum Benjamin-Ono equation fully solved in [55] are particularly promising in this regard. One ingredient that has been missing from our construction is the $r$-matrix theory, which has been systematically used in some cases for generating the Hamiltonian eigenstates, cf. the considerations for quantum bosons with contact interactions in [61, 62]. There is also a distinction between the Lax pairs (3.7) and (12.1) that our considerations revolve around, formulated as operators acting on an infinite-dimensional auxiliary space (this construction goes back to the original presentation due to Gérard and Grellier in [1]), and finite-sized Lax matrices depending on a spectral parameter commonly seen in treatments of inverse scattering theory. It could be useful to find reformulations of our algebraic structures that would bridge this gap with other approaches.

The Gérard-Grellier Hamiltonian forms a very special point in the space of theories (justifying the term ‘superintegrable’ we have been employing to characterize it) and admits deformations preserving its Lax integrability [29, 32]. Once its quantum integrability has been brought fully under control, those deformations will emerge as natural further targets for study.

Appendix A The classical Lax pair

To verify that the classical Lax pair

$$(\mathcal{L}\vec{h})_{n}=\sum_{k,l=0}^{\infty}\alpha_{n+k}\bar{\alpha}_{k+l}h_{l},\qquad(\mathcal{M}\vec{h})_{n}=\sum_{m+p\geq n}^{m,p\geq 0}\alpha_{p}\bar{\alpha}_{p+m-n}h_{m}.$$ (A.1)

works properly, we must prove that

$$i(\dot{\mathcal{L}}\vec{h})_{n}-(\mathcal{M}\mathcal{L}\vec{h})_{n}+(\mathcal{L}\mathcal{M}\vec{h})_{n}=0.$$ (A.2)

Using the classical equation of motion

$$i\dot{\alpha}_{n}=\sum_{m=0}^{\infty}\sum_{p=0}^{n+m}\bar{\alpha}_{m}\alpha_{p}\alpha_{n+m-p}\,,$$ (A.3)

we have

$$i(\dot{\mathcal{L}}\vec{h})_{n}=-\sum_{k,l,m=0}^{\infty}\sum_{p=0}^{k+l+m}\bar{\alpha}_{p}\bar{\alpha}_{k+l+m-p}\alpha_{m}\alpha_{k+n}h_{l}+\sum_{k,l,m=0}^{\infty}\sum_{p=0}^{k+n+m}\bar{\alpha}_{l+k}\bar{\alpha}_{m}\alpha_{p}\alpha_{k+n+m-p}h_{l}.$$ (A.4)

On the other hand,

$$(\mathcal{M}\mathcal{L}\vec{h})_{n}=\sum_{k,l=0}^{\infty}\sum_{m+p\geq n}^{m,p\geq 0}\bar{\alpha}_{m+p-n}\bar{\alpha}_{l+k}\alpha_{p}\alpha_{k+m}h_{l},$$ (A.5)

$$(\mathcal{L}\mathcal{M}\vec{h})_{n}=\sum_{k,l=0}^{\infty}\sum_{m+p\geq l}^{m,p\geq 0}\bar{\alpha}_{l+k}\bar{\alpha}_{m+p-l}\alpha_{k+n}\alpha_{p}h_{m},$$ (A.6)

We first relabel the indices in (A.6) to bring the summand to the same form as in the first term of (A.4). Namely, we rename $m$ into $l$, $p$ into $m$, and $l$ into $p$, obtaining

$$(\mathcal{L}\mathcal{M}\vec{h})_{n}=\sum_{k,p=0}^{\infty}\sum_{l+m\geq p}^{l,m\geq 0}\bar{\alpha}_{p+k}\bar{\alpha}_{l+m-p}\alpha_{k+n}\alpha_{m}h_{l}.$$ (A.7)

We then introduce $\tilde{p}=p+k$ and drop tildes to obtain

$$(\mathcal{L}\mathcal{M}\vec{h})_{n}=\sum_{k=0}^{\infty}\sum_{p=k}^{\infty}\sum_{l+m\geq p-k}^{l,m\geq 0}\bar{\alpha}_{p}\bar{\alpha}_{k+l+m-p}\alpha_{k+n}\alpha_{m}h_{l}.$$ (A.8)

At this point, the summand is identical to the first term of (A.4), but the summation range can be expressed more effectively. The summation range is defined by the inequalities $p\geq k$, $l,m\geq 0$, $l+m\geq p-k$, which are more compactly written as $l,m\geq 0$, $k\leq p\leq k+l+m$, yielding

$$(\mathcal{L}\mathcal{M}\vec{h})_{n}=\sum_{k,l,m=0}^{\infty}\sum_{p=k}^{k+l+m}\bar{\alpha}_{p}\bar{\alpha}_{k+l+m-p}\alpha_{k+n}\alpha_{m}h_{l}.$$ (A.9)

Then, we turn to (A.5) and attempt to redefine the indices to make the summand identical to the second term of (A.4). This is accomplished by introducing $\tilde{m}=m+p-n$, in terms of which the index range conditions $m,p\geq 0$, $m+p\geq n$ become $\tilde{m}\geq p-n$, $p\geq 0$, $\tilde{m}\geq 0$, or even better, $\tilde{m}\geq 0$, $0\leq p\leq\tilde{m}+n$. Dropping the tildes, we obtain

$$(\mathcal{M}\mathcal{L}\vec{h})_{n}=\sum_{k,l,m=0}^{\infty}\sum_{p=0}^{n+m}\bar{\alpha}_{m}\bar{\alpha}_{l+k}\alpha_{p}\alpha_{k+m+n-p}h_{l}.$$ (A.10)

Finally, with (A.4), (A.9) and (A.10),

$$(\dot{\mathcal{L}}\vec{h})_{n}-(\mathcal{M}\mathcal{L}\vec{h})_{n}+(\mathcal{L}\mathcal{M}\vec{h})_{n}= -\sum_{k,l,m=0}^{\infty}\sum_{p=0}^{k-1}\bar{\alpha}_{p}\bar{\alpha}_{k+l+m-p}\alpha_{m}\alpha_{k+n}h_{l}$$ (A.11)

$$+\sum_{k,l,m=0}^{\infty}\sum_{p=n+m+1}^{k+n+m}\bar{\alpha}_{l+k}\bar{\alpha}_{m}\alpha_{p}\alpha_{k+n+m-p}h_{l}.$$

In the second line, redefine $p$ to $n+m+k-p$, after that interchange $m$ and $p$, after that interchange the order of summations over $k$ and $m$ and define $\tilde{k}=k+p-m$, thereafter dropping the tildes and interchanging the order of summations of $k$ and $p$:

$$\sum_{k,l,m=0}^{\infty}\sum_{p=n+m+1}^{k+n+m}\bar{\alpha}_{l+k}\bar{\alpha}_{m}\alpha_{p}\alpha_{k+n+m-p}h_{l}=\sum_{k,l,m=0}^{\infty}\sum_{p=0}^{k-1}\bar{\alpha}_{l+k}\bar{\alpha}_{m}\alpha_{p}\alpha_{k+n+m-p}h_{l}$$

$$\hskip 28.45274pt=\sum_{k,l,p=0}^{\infty}\sum_{m=0}^{k-1}\bar{\alpha}_{l+k}\bar{\alpha}_{p}\alpha_{m}\alpha_{k+n+p-m}h_{l}=\sum_{m,l,p=0}^{\infty}\sum_{k=m+1}^{\infty}\bar{\alpha}_{l+k}\bar{\alpha}_{p}\alpha_{m}\alpha_{k+n+p-m}h_{l}$$

$$\hskip 56.9055pt=\sum_{m,l,p=0}^{\infty}\sum_{\tilde{k}=p+1}^{\infty}\bar{\alpha}_{l+\tilde{k}+m-p}\bar{\alpha}_{p}\alpha_{m}\alpha_{\tilde{k}+n}h_{l}=\sum_{k,l,m=0}^{\infty}\sum_{p=0}^{k-1}\bar{\alpha}_{k+l+m-p}\bar{\alpha}_{p}\alpha_{m}\alpha_{k+n}h_{l}.$$

Thus, the first and second lines of (A.11) are identical with a relative minus sign, and hence (A.2) holds, ensuring the validity of the classical Lax pair.

Appendix B Some further details for the proof that $[H,H_{\min}]=0$

To complete the proof in section 6 that $[H,H_{\min}]=0$, we need to demonstrate that $A$ defined by (6.7) satisfies $A_{kjj^{\prime}}+A_{kj^{\prime}j}+A_{j^{\prime}jk}=0$. To do so, we introduce the position space field $u(z)$ according to

$$a_{n}=\frac{1}{2\pi i}\oint\frac{dz}{z^{n+1}}u(z),\qquad u(z)\equiv\sum_{n=0}^{\infty}a_{n}z^{n},$$ (B.1)

so that

$$\begin{split}&A_{kjj^{\prime}}=\frac{1}{(2\pi i)^{3}}\oint\frac{dz}{z}\oint\frac{dv}{v}\oint\frac{dw}{w}A_{kjj^{\prime}}(z,v,w),\\ &A(z,v,w)\equiv\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j^{\prime}+l}(\mathrm{min}_{j^{\prime}+l,ll^{\prime}}-\mathrm{min}_{j+k,lk})z^{l^{\prime}}v^{j+k-l}w^{j^{\prime}+l-l^{\prime}}.\end{split}$$ (B.2)

To evaluate $A_{kjj^{\prime}}(z,v,w)$, we need $\sum_{m=0}^{n}x^{m}=(1-x^{n+1})/(1-x)$, together with

$$\sum_{m=0}^{k+l}\mathrm{min}_{k+l,lm}x^{m}\equiv\sum_{m=0}^{k+l}\min(k,l,m,k+l-m)x^{m}=\sum_{m=0}^{k+l}\sum_{s=1}^{\min(k,l,m,k+l-m)}x^{m}$$

$$=\sum_{s=1}^{\min(k,l)}\sum_{m=s}^{k+l-s}x^{m}=\sum_{s=1}^{\min(k,l)}x^{s}\sum_{m=0}^{k+l-2s}x^{m}=\sum_{s=1}^{\min(k,l)}\frac{x^{s}-x^{k+l-s+1}}{1-x}$$

$$=\sum_{s=0}^{\min(k,l)-1}\frac{x^{s+1}-x^{k+l-s}}{1-x}=\frac{1}{1-x}\left[x\frac{1-x^{\min(k,l)}}{1-x}-x^{k+l}\frac{1-x^{-\min(k,l)}}{1-1/x}\right]$$

$$=\frac{x}{(1-x)^{2}}(1-x^{\min(k,l)}-x^{\max(k,l)}+x^{k+l})=\frac{x(1-x^{\min(k,l)})(1-x^{\max(k,l)})}{(1-x)^{2}}$$

$$\hskip 85.35826pt=\frac{x(1-x^{k})(1-x^{l})}{(1-x)^{2}}.$$

Then,

$$\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j^{\prime}+l}\mathrm{min}_{j^{\prime}+l,ll^{\prime}}z^{l^{\prime}}v^{j+k-l}w^{j^{\prime}+l-l^{\prime}}=v^{j+k}w^{j^{\prime}}\sum_{l=0}^{j+k}\frac{w^{l}}{v^{l}}\frac{z(1-z^{j^{\prime}}/w^{j^{\prime}})(1-z^{l}/w^{l})}{w(1-z/w)^{2}}$$

$$=\frac{zv^{j+k}w^{j^{\prime}+1}}{(w-z)^{2}}(1-z^{j^{\prime}}/w^{j^{\prime}})\sum_{l=0}^{j+k}\left(\frac{w^{l}}{v^{l}}-\frac{z^{l}}{v^{l}}\right)$$

$$=\frac{zwv^{j+k}}{(w-z)^{2}}(w^{j^{\prime}}-z^{j^{\prime}})\left(\frac{1-(w/v)^{j+k+1}}{1-w/v}-\frac{1-(z/v)^{j+k+1}}{1-z/v}\right)$$

$$=\frac{zw}{(w-z)^{2}}(w^{j^{\prime}}-z^{j^{\prime}})\left(\frac{v^{j+k+1}-w^{j+k+1}}{v-w}-\frac{v^{j+k+1}-z^{j+k+1}}{v-z}\right)$$

$$=\frac{v^{j+k+1}zw(w^{j^{\prime}}-z^{j^{\prime}})}{(w-z)(v-w)(v-z)}-\frac{zw^{j+k+2}(w^{j^{\prime}}-z^{j^{\prime}})}{(w-z)^{2}(v-w)}+\frac{z^{j+k+2}w(w^{j^{\prime}}-z^{j^{\prime}})}{(w-z)^{2}(v-z)}.$$ (B.3)

Similarly,

$$\sum_{l=0}^{j+k}\sum_{l^{\prime}=0}^{j^{\prime}+l}\mathrm{min}_{j+k,lk}z^{l^{\prime}}v^{j+k-l}w^{j^{\prime}+l-l^{\prime}}=v^{j+k}w^{j^{\prime}}\sum_{l=0}^{j+k}\mathrm{min}_{j+k,lk}\frac{w^{l}}{v^{l}}\frac{1-(z/w)^{j^{\prime}+l+1}}{1-z/w}$$

$$=\frac{v^{j+k}w^{j^{\prime}+1}}{w-z}\sum_{l=0}^{j+k}\mathrm{min}_{j+k,lk}\frac{w^{l}}{v^{l}}-\frac{v^{j+k}z^{j^{\prime}+1}}{w-z}\sum_{l=0}^{j+k}\mathrm{min}_{j+k,lk}\frac{z^{l}}{v^{l}}$$

$$=\frac{vw^{j^{\prime}+2}(v^{j}-w^{j})(v^{k}-w^{k})}{(w-z)(v-w)^{2}}-\frac{vz^{j^{\prime}+2}(v^{j}-z^{j})(v^{k}-z^{k})}{(w-z)(v-z)^{2}}$$ (B.4)

By subtracting (B.4) from (B.3), we obtain $A_{kjj^{\prime}}(z,v,w)$. Since the $(z,v,w)$-integration in (B.2) is fully symmetric under all permutations of $z$, $v$ and $w$, while $A_{kjj^{\prime}}(z,v,w)$ is only symmetric with respect to interchanging $z$ and $w$ it is wise to equivalently consider the symmetrized combination

$$A^{\mathrm{symm}}_{kjj^{\prime}}(z,v,w)\equiv\frac{1}{3}\left[A_{kjj^{\prime}}(z,v,w)+A_{kjj^{\prime}}(v,z,w)+A_{kjj^{\prime}}(z,w,v)\right]$$

$$=\frac{zvw}{3(w-v)(w-z)(v-z)}\left[P_{kjj^{\prime}}(v,z)+P_{kjj^{\prime}}(z,w)+P_{kjj^{\prime}}(w,v)\right]$$

$$-\frac{wv^{2}}{3(v-z)(v-w)^{2}}\left[P_{kjj^{\prime}}(v,w)+P_{kjj^{\prime}}(w,v)\right]-\frac{zw^{2}}{3(w-z)^{2}(w-v)}\left[P_{kjj^{\prime}}(z,w)+P_{kjj^{\prime}}(w,z)\right]$$

$$\hskip 85.35826pt+\frac{vz^{2}}{3(w-z)(v-z)^{2}}\left[P_{kjj^{\prime}}(z,v)+P_{kjj^{\prime}}(v,z)\right],$$

with $P_{abc}(x,y)\equiv 2x^{a+b}y^{c}-x^{b+c}y^{a}-x^{a+c}y^{b}$. But since $P_{abc}+P_{acb}+P_{cba}=0$, $A^{\mathrm{symm}}_{kjj^{\prime}}(z,v,w)+A^{\mathrm{symm}}_{kj^{\prime}j}(z,v,w)+A^{\mathrm{symm}}_{j^{\prime}jk}(z,v,w)=0$, and hence by (B.2), $A_{kjj^{\prime}}+A_{kj^{\prime}j}+A_{j^{\prime}jk}=0$. This proves that the last term in the square brackets of (6.4) gives a vanishing contribution, while all the other terms have been proven earlier to give vanishing contributions as well, hence $[H_{0},H_{1}]=0$, and therefore $[H,H_{\min}]=0$.

Appendix C Action of $K_{\pm}$

In this appendix we prove that $K_{+}$ (8.6) maps between next-to-highest eigenvectors in the top $H_{\min}$ eigenbasis, i.e.,

$$K_{+}\ket{\psi^{(1,2)}_{N,M}}=\sqrt{N(N-1)(M-2)}\ket{\psi^{(1,2)}_{N+1,M-1}}\,,$$ (C.1)

with

$$\ket{\psi^{(1,2)}_{N,M}}=\sum_{k=1}^{\min(N,M)}v^{(1,2)}_{k,{\scriptscriptstyle{NM}}}\frac{a_{0}^{\dagger(N-k)}}{\sqrt{(N-k)!}}S|\psi^{\max}_{k,M-k}\rangle\,.$$ (C.2)

We start by noting that $K_{+}$ commutes with $a^{\dagger}_{0}$ and transporting $K_{+}$ through the shift operator gives

$$K_{+}S=SK_{+}^{(2)}+2a^{\dagger}_{0}S\tilde{Z}^{\dagger}+a^{\dagger 2}_{0}Sa_{0}$$ (C.3)

with

$$K_{+}^{(2)}\equiv\sum_{n,m=0}^{\infty}a^{\dagger}_{n}a^{\dagger}_{m}a_{n+m+2}\,,\qquad\tilde{Z}^{\dagger}\equiv\sum_{n=0}^{\infty}a^{\dagger}_{n}a_{n+1}\,.$$ (C.4)

The operator $K_{+}$ is hence found to be tridiagonal in the basis $|\phi^{(1,k)}_{N,M}\rangle$ given by (8.2). Let us compute each contribution separately. First, it is straightforward to show by direct computation that

$$K^{(2)}_{+}\ket{\psi^{\max}_{N,M}}=\sqrt{\frac{M(M-1)N(N+1)}{N+M-1}}\ket{\psi^{\max}_{N+1,M-2}}$$ (C.5)

using the explicit expression for the highest energy state (7.5) and the following commutation relations

$$[K^{(2)}_{+},Z]=2K_{+}\,,\qquad[K_{+},Z]=J_{3}^{\dagger}\,\qquad\text{and}\qquad[J_{3}^{\dagger},Z]=0\,.$$ (C.6)

The second contribution in (C.3) can be found using

$$[\tilde{Z}^{\dagger},Z^{m}]=mNZ^{m-1}\,,$$ (C.7)

which leads to

$$\tilde{Z}^{\dagger}\ket{\psi^{\max}_{N,M}}=N\sqrt{\frac{M}{N+M-1}}\ket{\psi^{\max}_{N,M-1}}\,.$$ (C.8)

Finally, the third contribution to (C.3) is easily found using (8.8), which we copy here for ease

$$a_{0}\ket{\psi^{\max}_{N,M}}=\sqrt{\frac{N(N-1)}{N+M-1}}\ket{\psi^{\max}_{N-1,M}}\,.$$ (C.9)

Combining (C.3), (C.5), (C.8) and (C.9) one finds

$$K_{+}\ket{\psi^{(1,2)}_{N,M}}= \sum_{k=1}^{\min(N,M-2)}v^{(1,2)}_{k,{\scriptscriptstyle{NM}}}\sqrt{\frac{(M-k)(M-k-1)k(k+1)}{M-1}}\frac{a_{0}^{\dagger(N-k)}}{\sqrt{(N-k)!}}S\ket{\psi^{\max}_{k+1,M-k-2}}\,,$$

$$+ \sum_{k=1}^{\min(N,M-1)}2v^{(1,2)}_{k,{\scriptscriptstyle{NM}}}k\sqrt{\frac{M-k}{M-1}}\frac{a_{0}^{\dagger(N-k+1)}}{\sqrt{(N-k)!}}S\ket{\psi^{\max}_{k,M-k-1}}\,,$$

$$+ \sum_{k=2}^{\min(N,M)}v^{(1,2)}_{k,{\scriptscriptstyle{NM}}}\sqrt{\frac{k(k-1)}{M-1}}\frac{a_{0}^{\dagger(N-k+2)}}{\sqrt{(N-k)!}}S\ket{\psi^{\max}_{k-1,M-k}}\,.$$ (C.10)

The coefficients $v^{(1,2)}_{k,{\scriptscriptstyle{NM}}}$ can be computed following the recursion steps detailed at the end of section 8:

$$v^{(1,2)}_{k,{\scriptscriptstyle{NM}}} =\frac{(2n+a-1)k-n(n+a)}{\sqrt{k(2n+a-1)(2n+a-3)}}\sqrt{\frac{(n-2)!(n+a-2)!}{(2n+a-4)!}}\sqrt{\binom{n-1}{k-1}\binom{n+a-1}{k-1}}$$

$$=\frac{(N+M-1)k-NM}{\sqrt{k(N+M-1)(N+M-3)}}\sqrt{\frac{(N-2)!(M-2)!}{(N+M-4)!}}\sqrt{\binom{N-1}{k-1}\binom{M-1}{k-1}}\,.$$ (C.11)

Substituting (C.11) in (C.10) and rearranging terms leads to (C.1).