BdG & para-Hermitian diagonalization

How to diagonalize a quadratic bosonic Hamiltonian that contains anomalous (pairing) terms ${\hat{b}}^{†}{\hat{b}}^{†}$ and $\hat{b}\hat{b}$ — the situation one always lands in after a Holstein–Primakoff expansion of a magnet, or for any Bogoliubov problem with a condensate. The naive instinct is to stack the operators into a vector and diagonalize the coefficient matrix like an ordinary Hermitian eigenproblem. That is wrong, and the reason is structural rather than numerical: the transformation that relates old and new bosons must preserve the bosonic commutation relations, which is not unitarity but para-unitarity — unitarity with respect to an indefinite metric $g=\mathrm{diag}(1_n,-1_n)$. Get the metric right and two good things follow at once: the physical magnon frequencies come out correctly, and the positive-definiteness of the coefficient matrix is exactly the dynamical stability test of the reference state. This note is the linear-algebra core that LSWT invokes (it is its §3–4). The fermionic superconducting BdG problem is a different object — same doubling idea, but a positive- definite (Euclidean) metric and a particle–hole charge-conjugation constraint, so eigenvalues come in $\pm E$ pairs — and I do not develop it here.

The quadratic bosonic Hamiltonian and the BdG form

After expanding to quadratic order one has, at each wavevector $\mathbf{k}$, bosons $b_{\mathbf{k},1},\dots ,b_{\mathbf{k},n}$ (sublattice/band index) with a Hamiltonian containing both normal ($b^{†}b$) and anomalous ($b^{†}{b}^{†}$, $bb$) terms. Collect them in a Nambu-like spinor

$${\mathbf{\beta }}_{\mathbf{k}}=(b_{\mathbf{k},1},\dots ,b_{\mathbf{k},n},b_{-\mathbf{k},1}^{†},\dots ,b_{-\mathbf{k},n}^{†}{)}^{\top }(\text{length }2n),$$

so that (dropping a constant zero-point shift)

$$\hat{H}=\frac{1}{2}\sum _{\mathbf{k}}{\mathbf{\beta }}_{\mathbf{k}}^{†}\mathcal{H}(\mathbf{k}){\mathbf{\beta }}_{\mathbf{k}},\mathcal{H}(\mathbf{k})=\left(\begin{array}{cc}A(\mathbf{k})& B(\mathbf{k})\\ B^{†}(\mathbf{k})& A^{\top }(-\mathbf{k})\end{array}\right).$$

$\mathcal{H}(\mathbf{k})$ is the BdG (Bogoliubov–de Gennes) matrix, $2n\times 2n$ and Hermitian, $\mathcal{H}={\mathcal{H}}^{†}$. Hermiticity of $\hat{H}$ forces $A=A^{†}$; bosonic exchange symmetry of the pairing block forces $B=B^{\top }$; and the lower-right block is fixed to $A^{\top }(-\mathbf{k})$ by the doubling. Crucially, the doubled basis is redundant — the lower $n$ components are the adjoints of (the $-\mathbf{k}$ partners of) the upper $n$ — so $\mathcal{H}$ obeys a particle–hole constraint

$${\tau }_x{\mathcal{H}}^{\ast }(-\mathbf{k}){\tau }_x=\mathcal{H}(\mathbf{k}),{\tau }_x=\left(\begin{array}{cc}0& 1_n\\ 1_n& 0\end{array}\right),$$

which is what guarantees the spectrum below comes in $\pm \omega$ pairs and that the positive and negative branches carry the same physics.

Why Hermitian diagonalization is the wrong eigenproblem

The commutation relations in the doubled basis are not the identity. With $[b_a,b_b^{†}]={\delta }_{ab}$ one finds

$$[{\mathbf{\beta }}_{\mathbf{k}},{\mathbf{\beta }}_{\mathbf{k}}^{†}]=g,g=\left(\begin{array}{cc}{1}_n& 0\\ 0& -1_n\end{array}\right)=\mathrm{diag}(1_n,-1_n),$$

the indefinite (para-unitary) metric: the second block of operators are creation operators, and their commutator carries the opposite sign. A change of basis $\mathbf{\beta }=T\mathbf{\gamma }$ preserves the bosonic algebra of the new operators $\mathbf{\gamma }$ iff

$$Tg{T}^{†}=g(\text{para-unitarity, }T\in U(n,n)),$$

not $T{T}^{†}=1$. So the allowed transformations form the non-compact group $U(n,n)$, not $U(2n)$. The Heisenberg equation of motion makes the consequence concrete:

$$i{\partial }_t\mathbf{\beta }=[\mathbf{\beta },\hat{H}]=g\mathcal{H}\mathbf{\beta },$$

so the dynamical matrix whose eigenvalues are the physical frequencies is $g\mathcal{H}$, not $\mathcal{H}$. $g\mathcal{H}$ is non-Hermitian (it is the product of two Hermitian matrices), yet because $\mathcal{H}$ is Hermitian its spectrum is constrained — real when $\mathcal{H}$ is positive-definite (below), and otherwise complex-conjugate pairs that signal instability. Diagonalizing $\mathcal{H}$ itself, as an ordinary Hermitian matrix, would silently replace $g$ by $1$ and return eigenvalues of $\mathcal{H}$, which are neither the frequencies nor anything physical: it is the right matrix paired with the wrong metric. This is the single most common error in hand-rolled spin-wave codes.

The Bogoliubov transformation as a para-unitary diagonalization

The goal is a para-unitary $T$ ($Tg{T}^{†}=g$) bringing $\hat{H}$ to a sum of independent oscillators,

$$T^{†}\mathcal{H}T=\mathcal{E}=\mathrm{diag}({\omega }_1,\dots ,{\omega }_n,{\omega }_1,\dots ,{\omega }_n),\hat{H}=\sum _{\mathbf{k}}\sum _{m=1}^n{\omega }_m(\mathbf{k})({\gamma }_{\mathbf{k}m}^{†}{\gamma }_{\mathbf{k}m}+\frac{1}{2}),$$

with $\mathbf{\gamma }=T^{-1}\mathbf{\beta }$ the Bogoliubov quasiparticles (magnons) and ${\omega }_m>0$ the magnon dispersion. Equivalently, because $T^{-1}=g{T}^{†}g$ for a para-unitary $T$, the diagonalization is an eigenproblem of the non-Hermitian dynamical matrix,

$$g\mathcal{H}T=T(g\mathcal{E}),g\mathcal{E}=\mathrm{diag}({\omega }_1,\dots ,{\omega }_n,-{\omega }_1,\dots ,-{\omega }_n).$$

The columns of $T$ are the right eigenvectors of $g\mathcal{H}$ and the paraunitary normalization $Tg{T}^{†}=g$ fixes their (indefinite) lengths — the positive-norm eigenvectors are the magnon annihilation modes, the negative-norm ones their conjugate partners, and the $\pm \omega$ pairing is the particle–hole constraint above. Solving $g\mathcal{H}$ directly works but is numerically delicate: $g\mathcal{H}$ is non-Hermitian, its eigenvectors are not orthogonal in the ordinary sense, and one must sort eigenvalues into $\pm$ branches and re-normalize against $g$ by hand. Colpa’s method below replaces this with two Hermitian operations and makes the stability test fall out for free.

Positive-definiteness is the stability test

Before any diagonalization, the meaningful question is whether the chosen reference state (the classical magnetic ground state the bosons fluctuate around) is actually a local minimum. In the quadratic theory the energy of a fluctuation is $\frac{1}{2}{\mathbf{\beta }}^{†}\mathcal{H}\mathbf{\beta }$, so the reference state is dynamically stable iff $\mathcal{H}(\mathbf{k})\succ 0$ (positive-definite) for every $\mathbf{k}$. If $\mathcal{H}$ has a zero eigenvalue there is a flat (Goldstone) direction — fine at $\Gamma$ for a continuous-symmetry-broken magnet, expected. If $\mathcal{H}$ acquires a negative eigenvalue at some $\mathbf{k}$, the energy is unbounded below along that fluctuation: the assumed order is unstable (a competing modulated/incommensurate state lowers the energy), and correspondingly $g\mathcal{H}$ develops complex-conjugate frequency pairs — imaginary magnon energies. The positive-definiteness of $\mathcal{H}$ is therefore not a precondition one hopes for but the diagnostic itself: a failed Cholesky factorization at a given $\mathbf{k}$ is the code telling you the magnetic structure is wrong there. (This is the bosonic analogue of the imaginary-mode diagnostic in lattice dynamics — same logic, different metric.)

Colpa’s Cholesky route (the positive-definite case)

When $\mathcal{H}(\mathbf{k})\succ 0$, Colpa’s algorithm diagonalizes the para-Hermitian problem using only Hermitian, numerically stable steps:

1. Cholesky-factor the BdG matrix. Positive-definiteness guarantees a unique upper-triangular $\mathbf{K}$ with $$\mathcal{H}={\mathbf{K}}^{†}\mathbf{K}.$$ (The factorization failing — a non-positive pivot — is precisely the stability failure above; so the very first step is the test.)
2. Form and diagonalize a Hermitian matrix. The matrix $$\mathbf{W}\equiv \mathbf{K}g{\mathbf{K}}^{†}$$ is Hermitian (since $g^{†}=g$), so it has a unitary eigendecomposition ${\mathbf{U}}^{†}\mathbf{WU}=\mathbf{L}$, $\mathbf{L}$ real diagonal, $\mathbf{U}$ unitary. This is an ordinary Hermitian eigenproblem — LAPACK zheev, well-conditioned.
3. Read off the physical frequencies. $\mathbf{W}=\mathbf{K}g{\mathbf{K}}^{†}$ and $g\mathcal{H}=g{\mathbf{K}}^{†}\mathbf{K}$ share the same spectrum by the $\mathrm{eig}(AB)=\mathrm{eig}(BA)$ identity (take $A=\mathbf{K}$, $B=g{\mathbf{K}}^{†}$; $\mathbf{K}$ is invertible, so the full spectra coincide, not merely the nonzero part). The real eigenvalues $\mathbf{L}$ of the Hermitian $\mathbf{W}$ are therefore exactly the eigenvalues of the dynamical matrix $g\mathcal{H}$ — the set $\{{\omega }_1,\dots ,{\omega }_n,-{\omega }_1,\dots ,-{\omega }_n\}$. The physical magnon frequencies are $g$ times the eigenvalues of $\mathbf{W}$, $$\mathcal{E}=g\mathbf{L}=\mathrm{diag}({\omega }_1,\dots ,{\omega }_n,{\omega }_1,\dots ,{\omega }_n)\succ 0,$$ with the diagonal of $g\mathbf{L}$ arranged so the first $n$ (positive-norm) entries are the positive frequencies ${\omega }_m>0$ and the last $n$ their partners. Positive-definiteness of $\mathcal{H}$ guarantees every ${\omega }_m$ is real and positive.
4. Assemble the para-unitary transformation. The Bogoliubov matrix is $$T={\mathbf{K}}^{-1}\mathbf{U}{\mathcal{E}}^{1/2},$$ which one checks satisfies both $T^{†}\mathcal{H}T=\mathcal{E}$ (so $\hat{H}$ is diagonal) and $Tg{T}^{†}=g$ (so the new operators are bona fide bosons). The ${\mathcal{E}}^{1/2}$ factor is exactly the paraunitary normalization that the bare eigenvectors of the non-Hermitian $g\mathcal{H}$ would otherwise require by hand.

The whole procedure is two Hermitian factorizations (a Cholesky and an eigendecomposition) plus a triangular solve — fast, stable, and self-checking, which is why it is the standard back end in spin-wave packages (SpinW, and equivalents). The $u,v$ Bogoliubov coefficients of the classic two-sublattice antiferromagnet, $\tanh 2{\theta }_{\mathbf{k}}=|B_{\mathbf{k}}|/A_{\mathbf{k}}$, ${\omega }_{\mathbf{k}}=\sqrt{A_{\mathbf{k}}^2-|B_{\mathbf{k}}{|}^2}$, are the $n=1$ special case of all of the above — and the condition $A_{\mathbf{k}}>|B_{\mathbf{k}}|$ for real ${\omega }_{\mathbf{k}}$ is just $\mathcal{H}\succ 0$ written out.

Prerequisites

bosonic second quantization · Hermitian eigenproblems · Cholesky factorization

Builds toward: LSWT (this is its §3–4)

Key references

• The Colpa diagonalization — J. H. P. Colpa, Physica A 93, 327 (1978).
• Bogoliubov transformation for bosons with anomalous terms — R. M. White, M. Sparks & I. Ortenburger, Phys. Rev. 139, A450 (1965).
• The original transformation — N. N. Bogoliubov, J. Phys. (USSR) 11, 23 (1947).