Consider a set of \(N\) second-quantized creation operators \[\vec{a} = [\hat{a}_1, \hat{a}_2, \cdots, \hat{a}_N]\] which satisfy the canonical (anti-)commutation relations (CCR), \[[\hat{a}_i, \hat{a}_j^{\dagger}]_{\pm} = \delta_{ij},\] depending on whether they describe bosons or fermions. Suppose we are interested in a system described by a Hamiltonian that is quadratic in these operators, \(H=\sum_{ij} \hat{a}_i^{\dagger}h_{ij}\hat{a}_j + \frac{1}{2} \hat{a}_i^{\dagger}\Delta_{ij}\hat{a}^{\dagger}_j +\frac{1}{2}\hat{a}_i \Delta_{ij}^* \hat{a}_j\). Generically, this may be cast into a matrix equation like so (upto some constant terms),
\[ H = \frac{1}{2} \begin{bmatrix} \vec{a}^\dagger & \vec{a} \end{bmatrix} \begin{bmatrix} h &\Delta \\ \pm \Delta^{\dagger} & \pm h^{T} \end{bmatrix} \begin{bmatrix} \vec{a} \\ \vec{a}^\dagger \end{bmatrix} = \frac{1}{2} \vec{A}^{\dagger} H_{BdG} \vec{A}, \] where \(h\) is generally an \(N \times N\) Hermitian matrix and \(\Delta\) is an \(N \times N\) anti-symmetric matrix describing pairing terms. The positive (negative) sign corresponds to bosons (fermions) due to the CCR. Note that we have introduced a stacked ‘Nambu’ vector \(\vec{A} = [\hat{a}_1, \hat{a}_2, \cdots, \hat{a}_N, \hat{a}^{\dagger}_1, \hat{a}^{\dagger}_2, \cdots, \hat{a}^{\dagger}_N]\) such that the CCR is now written as follows
\[ [A_i, A_j^{\dagger}]_{\pm} = \Sigma_{ij} \]
where \(\Sigma_{ij} = \begin{bmatrix} I_N & 0_N \\ 0_N & \pm I_N \end{bmatrix}\) since \(A\) consists of both creation and annihilation operators.
In such a case, it is convenient to look for the energy spectrum of the system by finding a new set of operators, \[B = [\hat{b}_1, \hat{b}_2, \cdots, \hat{b}_N, \hat{b}^{\dagger}_1, \hat{b}^{\dagger}_2, \cdots, \hat{b}^{\dagger}_N],\] related to \(A\) through a linear transformation \(T\) (note that we interpret \(\vec{A}\) and \(\vec{A}^{\dagger}\) as column and row vectors, respectively): \[\vec{A} = T \vec{B} \quad \implies \quad \vec{A}^\dagger = \vec{B}^\dagger T^\dagger.\]
such that \(H = \frac{1}{2} \vec{B}^{\dagger}T^{\dagger} H_{BdG} T\vec{B} = \frac{1}{2} \vec{B}^{\dagger}\Lambda\vec{B} = \sum_{i} \Lambda_{i} \hat{B}^{\dagger}_i\hat{B}_i\) where \(\Lambda\) is diagonal. The transformation \(T\) must chosen such that \(\vec{B}\) also satisfies the CCR in order to ensure a physical interpretation. This is the general framework of Bogoliubov transformations.
Such a transformation necessarily mixes the creation and annihilation operators and invites the interpretation of quasi-particles that have characteristics of holes and particles created by the original operators. In some places, this procedure can be referred to as ‘diagonalizing’ the system, i.e, solving the free part, but such a notion typically adds to the confusion because \(T^{\dagger}H_{BdG}T\) is not a similarity transform in general.
We now introduce \(Q = T^{-1}\) such that \(\vec{B} = Q \vec{A}\). The constraint on \(Q\) follows from a simple calculation enforcing the CCR.
\[ \begin{align} [Q_{im}A_m, A_n^{\dagger} Q^{\dagger}_{nj}]_{\pm} &= \Sigma_{ij} \\ Q_{im} Q^{\dagger}_{nj} [A_m, A_n^{\dagger}]_{\pm} = \Sigma_{ij} \\ Q_{im}\Sigma_{mn}Q_{nj}^{\dagger} = \Sigma_{ij} \end{align} \]
which means that \(Q \, \Sigma \, Q^\dagger = \Sigma.\) This tells us that for bosons, \(Q \in \mathrm{Sp}(2N, \mathbb{C})\) the complex symplectic group, while for fermions, \(Q \in \mathrm{U}(2N)\), the unitary group. Note that this holds for the inverse transformation \(T = Q^{-1}\) as well (i.e, \(T \Sigma T^{\dagger} = \Sigma\)). Perhaps rather notably, for fermions this remains a similarity transform, so the eigenvalues may simply be obtained by unitarily diagonalizing \(H_{BdG}\), however, for bosons, it is necessarily different.
We now need to identify an appropriate \(T\) such that \(T^{\dagger}H_{BdG}T = \Lambda\) is diagonal. This can be done with some simple manipulations as follows. Since we have \(\Sigma^2 = I_{2N}\), it follows that \(T^{\dagger} = \Sigma T^{-1} \Sigma\). Then, \[ \begin{align} T^{\dagger} H_{BdG} T &= \Lambda\\ (\Sigma T^{-1} \Sigma) H_{BdG} T &= \Lambda \\ T^{-1} \Sigma H_{BdG} T &= \Sigma \Lambda\\ \Sigma H_{BdG} T &= T \Sigma \Lambda \end{align} \]
Since \(\Lambda = \text{diag}(E_1, \dots, E_N, E_{N+1} \dots, E_{2N})\) is diagonal, we know that \(\Sigma \Lambda = \tilde\Lambda = \text{diag}(E_1, \dots, E_N, \pm E_{N+1}, \dots, \pm E_{2N})\) remains diagonal. Writing \(T = [\phi_1 \ \dots \phi_{2N}]\) in terms of its column vectors, we now have:
\[ \begin{equation} \Sigma H_{BdG} \phi_i = \tilde \Lambda_i \phi_i \end{equation} \]
which is a normal eigenvalue equation where we diagonalize the matrix, \[ \Sigma H_{BdG} = \begin{bmatrix} h &\Delta \\ -\Delta^{\dagger} & -h^{T} \end{bmatrix}\] for both bosons and fermions. We can actually infer something about the spectrum in general using the structure of \(H_{BdG}\). In both cases, we have \(\Omega (\Sigma H_{BdG}) \Omega = -(\Sigma H_{BdG})^*\) where,
\[ \Omega = \begin{bmatrix} 0_N & I_N \\ I_N & 0_N \end{bmatrix}.\]
This immediately means that for every eigenvector \(\vec{v} = [\vec{v}_1, \vec{v}_2]\) with eigenvalue \(\epsilon\), there is a corresponding pair \(\vec{v}' = [\vec{v}_2^*, \vec{v}_1^*]\) with eigenvalue \(-\epsilon\) in \(\tilde \Lambda\). The physical spectrum is then given by \(\Lambda = \Sigma \tilde \Lambda\).
In the case of fermions, \(\Lambda = \tilde \Lambda\), so the ground state already has all \(N\) negative energy states occupied, with the remaining \(N\) positive energies indicating physical excitations. The structure in the eigenvalues is interpreted as a manifestation of particle-hole symmetry in the original Hamiltonian \(H_{BdG}\). For bosons on the other hand, the \(\Sigma\) flips back the negative half of the spectrum, seemingly resulting in a doubly degenerate spectrum of positive energies in \(\Lambda\). This degeneracy, however, simply results from doubling the size of the Hamiltonian while constructing \(H_{BdG}\) and there are only \(N\) physically distinct excitations.
Before wrapping up, let us consider the case when a transformation may be valid for both bosons and fermions, i.e, it is both symplectic and unitary over \(2N \times 2N\) matrices. Suppose we consider a general Bogoliubov transformation as follows (since \(\vec{a} \mapsto M \vec{a} \implies \vec{a}^{\dagger} \mapsto M^* \vec{a}^{\dagger}\), considering both cases as column vectors as they are stacked in \(\vec{A}\)),
\[ Q = \begin{bmatrix} A & B \\ B^* & A^* \end{bmatrix} \]
From the bosonic and fermionic constraints, we have \(Q\Sigma_{-} Q^{\dagger} = \Sigma_{-}\) and \(Q Q^{\dagger} = I_{2N}\), respectively. Then \(Q\Sigma_{-} = \Sigma_{-} Q\), which immediately enforces \(B = 0_N\). The fermionic constraint then implies unitarity on the remaining blocks, \(AA^{\dagger} = I_N\). We then have,
\[ Q = \begin{bmatrix} A & 0_N \\ 0_N & A^* \end{bmatrix}, \hspace{1cm} A \in U(N) \]
which is a unitary transformation acting purely at the single particle level without mixing the creation and annihilation operators. This is perfectly consistent, since one would expect that the only transformations that are not affected by the exchange statistics must be a basis change of the single particle orbitals.