How can a Bogoliubov transformation be implemented numerically?

Question

I'm working with an antiferromagnetic spin model on a bipartite lattice, which I am trying to analyze via Holstein-Primakoff spin-wave theory. Because the theory has two sublattices, there are naturally two boson flavors $a,b$. In Fourier space, I've found an exact form for the Hamiltonian in terms of the bosons:

$$\sum_{\vec{k}}\begin{pmatrix} a_{k}^{\dagger} && b_{k}^{\dagger} && a_{-k} && b_{-k} \end{pmatrix} \begin{pmatrix} c && w^{*} && 0 && x^{*}\\ w && c && x && 0\\ 0 && x^{*} && c && w^{*}\\ x && 0 && w && c \end{pmatrix} \begin{pmatrix} a_{k}\\ b_{k}\\ a_{-k}^{\dagger}\\ b_{-k}^{\dagger} \end{pmatrix} $$

I have exact expressions for the parameters $c, x$, and $w$, and they depend on $k$. The trouble is that I'm not aware of a simple way to implement a Bogoliubov transformation for Hamiltonians of this type - there are some papers on spin-wave theory which address similar problems, but the examples I have seen do not have the same structure of couplings, so those closed forms do not apply.

For my purposes, I don't need a closed form - just a numerical procedure for extracting the band structure. Is there such an approach, in general? Of course, I'm happy to get a closed form for the band structure, but all I want is a well-defined numerical procedure. The difficulty is that Bogoliubov transformations are non-unitary, so brute-forcing the relevant equations seems computationally expensive. Besides, it would be useful to know how to do this in general.

Answer 1

The diagonalization of general quadratic bosonic Hamiltonians is discussed in J.H.P. Colpa, Physica A: Statistical Mechanics and its Applications 93, 327 (1978). It's a problem commonly encountered in spin-wave theories for noncollinear magnets with a higher number of sublattices, for example.

The central results are as follows. If you can write your Hamiltonians in the form you have, i.e. $$ H = \vec{\Psi}_k^\dagger h(k) \vec{\Psi}_k $$ where $\vec{\Psi}^\dagger_k=\left( a_1^\dagger(k), \dots, a_n^\dagger(k), a_1(-k),\dots, a_n(-k)\right)$, and $h(k)$ is a $2n\times 2n$ positive semi-definite matrix, then you can generally perform a so-called paraunitary diagonalization by (numerically) diagonalizing $\sigma_3 h(k)$, where $$ \sigma_3 = \begin{pmatrix} I_{n\times n} && 0\\ 0 && -I_{n\times n}\\ \end{pmatrix}$$ and $I_{n\times n}$ is the $n\times n$ identity-matrix. The structure of $\sigma_3$ ensures that bosonic commutation relations are preserved. (In the fermionic case, $\sigma_3$ can be replaced by the $2n\times 2n$ unit matrix.) Because of particle-hole symmetry in the enlarged space (in your case $4\times 4$ rather than $2\times 2$), eivenvalues show up in pairs $\pm \epsilon_n$. The physical eigenvalues are likely the positive ones. At least that tends to be the case in spin-wave treatments where we're interested in excitations.

So in a numerical approach, you can simply diagonalize $\sigma_3 h(k)$ for each $k$ and select eigenstates with positive eigenvalues. Before plotting the band structure it may be useful to sort the eigenvalues.