How can a Bogoliubov transformation be implemented numerically? 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.
 A: 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.
