# Getting the eigenvalues of a quadratic boson Hamiltonian numerically

I have the following quadratic Hamiltonian (of boson type):

$$\hat{H}=\epsilon b^\dagger b -v(b^\dagger b^\dagger + b b )$$, where both $$\epsilon$$ and $$v$$ are real parameters. The operators $$(b,b^\dagger)$$ Are the bosonic operators and they obey the usual commutation relations (specific for the bosonic case).

My task is to find the eigenvalues of this Hamiltonian, while working in the bases of the occupation number: $$|n\rangle =\frac{(b^\dagger)^n}{\sqrt{n}}|0\rangle$$, with $$n=0,1,2,3,\dots$$.

Finding numerically (e.g. in Mathematica) the eigenvalues of $$H$$ means to get the Hamiltonian in a matrix form, then just perform the calculation given by the so-called secular determinant $$\text{det}[H-\lambda I_n]=0$$.

The matrix form of H, given in terms of occupation number representation. Thus, the basis is given by: $$|0\rangle,|1\rangle,2\rangle \dots$$

Writing the matrix elements of $$H$$ as $$H_{ij}=\langle i|H|j\rangle$$, with the indices $$i,j$$ going through the entire basis, it is possible to obtain an $$n\times n$$ dimensional matrix which will be further introduced into the secular determinant in order to get its eigenvalues.

However, since the dimension of the space can be arbitrary ($$n$$ can be practically infinite), one cannot just diagonalize the matrix numerically for the entire space, so a truncation must be done.

So, if for example we stop at $$n=4$$, we only have the basis:

$$|0\rangle, |1\rangle, |2\rangle, |3\rangle$$

What is interesting is that, from the expression of the initial Hamiltonian, any matrix element: $$\langle i|H|j\rangle$$ will be non-zero only if the states $$i$$ and $$j$$ differ by 0 or 2 (since each of the three boson operator products will act on orthonormal states and produce results of the form $$c_ic_j\langle i| j \rangle)$$. I hope I’ve explained this clearly enough, if not, just try to compute $$\langle 1|H|2\rangle$$ and $$\langle 2|H|2\rangle$$.

In other words, the matrix associated with $$H$$ will be a tridiagonal matrix.

Because of this, there can be two cases:

• We take the sequence of bosonic states: $$|0\rangle,|2\rangle,|4\rangle,\dots$$
• We take the sequence of bosonic states: $$|1\rangle,|3\rangle,|5\rangle,\dots$$

Now, my approach (I’m using Mathematica by the way) was to take each case separately, and start with a low $$n$$, then keep increasing it by one unit. At each $$n$$, I try to compute by hand the matrix elements $$\langle i|H|j\rangle$$ (with $$(i,j)$$ being $$0,2,4,\dots$$ or $$1,3,5,\dots$$). Putting these elements in Mathematica and then just get the solutions of the secular determinant, I am able to obtain the $$\lambda$$’s.

Is my approach correct? Is this a proper way to obtain eigenvalues numerically?

Thank you in advance. (And I apologize in advance if my understanding of finding the eigenvalues of a bosonic Hamiltonian is wrong).

• I believe you already know that the spectrum of the above Hamiltonian is known exactly (can be computed using various approaches, simplest perhaps is through Bogoliubov transformation). Is the motivation for the question about numerics? Your said algorithm seems fine (infact you have exploited $\mathcal{Z}_2$ symmetry to permutationally transform the truncated matrix in occupation number basis which is pentadiagonal matrix in conventional ordering of basis into a block diagonal matrix with tridiagonal entries). – Sunyam Jul 1 at 15:38

What you want to do is use the $$\mathfrak{sp}(2,\mathbb{R})\sim \mathfrak{su}(1,1)$$ as your Hamiltonian can be expressed as combinations of the generators \begin{align} \hat K_{+}=\frac{1}{4}\hat a^\dagger \hat a^\dagger\, ,\qquad \hat K_{-}=\frac{1}{4}\hat a \hat a\, ,\qquad \hat K_0=\frac{1}{4}(\hat a^\dagger\hat a +\hat a\hat a^\dagger)\, , \end{align} and then make an $$SU(1,1)$$ transformation of your $$H$$ to eliminate cross terms. Specifically your Hamiltonian can be expressed in terms of the Hermitian operators $$\hat K_0$$ and $$K_x=\frac{1}{4} \left(\hat a^\dagger \hat a^\dagger+\hat a \hat a\right)$$ so an $$SU(1,1)$$ "rotation" $$U(\theta)=e^{i\theta \hat K_y}$$ might do the trick although depending on the ratio of $$v/\epsilon$$.

Indeed you can verify for yourself that a $$U$$ of the above form will transform $$\hat K_0$$ into a combo of $$\hat K_0$$ and $$\hat K_x$$, and will transform $$\hat K_x$$ into a combo of $$\hat K_0$$ and $$\hat K_x$$.

Note that, if $$v>> \epsilon$$, then you cannot expect this will work as clearly the eigenstates of $$H$$ would be "close to" the eigenstates of $$K_x$$, which are not normalizable: even a truncation scheme will fail. Truncation will only work in the "perturbative" approach where $$v<< \epsilon$$ else you might have to use an extremely large state space (and it might not even converge, as stated above).

The difference in the regimes is geometrically linked to two possible ways a plane can cut a hyperboloid. In one of these the result is a closed curve but in the other it's an open curve, and it's not possible to go between the two regimes by a continuous transformation. The boundary between the two regimes depends on the ratio $$v/\epsilon$$.

If this works then you will have to compute elements of the form $$\langle m\vert U\vert m'\rangle$$. There are $$SU(1,1)$$ group functions for which there are closed form expressions, given for instance in terms of Jacobi polynomials as per Eq.(3.20 b) of

Ui, H., 1970. Clebsch - Gordan formulas of the SU (1, 1) group, Progress of Theoretical Physics, 44 (3), pp.689 - 702

Similar expressions are also found in the appropriate chapter of

Perelomov, A., 2012. Generalized coherent states and their applications. Springer Science & Business Media

or

Vilenkin, N.J. and Klimyk, A.U., 2013. Representation of Lie groups and special functions: Volume 3: Classical and quantum groups and special functions (Vol. 75). Springer Science & Business Media

although I maintain the exposition of Ui is the clearest to follow.