Is BdG transformation unitarily equivalent? Consider the fermionic Hamiltonian of the form
$$H = \sum_{i,j=1}^n\left( \alpha_{ij} c^\dagger_i c_j + \frac12 \gamma_{ij} c^\dagger_i c^\dagger_j + \frac12 \gamma^*_{ji}c_ic_j \right).$$
In BdG formalism, we use the transformation of the form
$$a_i = A_{ij} c_j + B_{ij} c_j^\dagger, \quad i=1,\ldots, n,$$
with the matrix $A_{ij}$ and $B_{ij}$ is defined suitably to preserve the anticommutation relation and to make the Hamiltonian
$$H = \sum_{i=1}^n E_i \left( a^\dagger_i a_i - \frac12\right) + \mathrm{const}$$
in the simple form. This transformation is always possible (Theorem 38 of arXiv:0908.0787).
However, finding the transformed Hamiltonian is not the end of the problem. What we are really interested is the eigenvalues and eigenstates of $H$. To find them, a natural idea is the following:

*

*Try to find the ground state, characterized by the condition $a_i |\psi\rangle = 0$ for each $i=1,\ldots, n$.


*Obtain the excited states by acting $a^\dagger_i$ on $|\psi\rangle$.
However, this approach gives me two questions:

*

*Is the characterization in 1. well-defined? In other words, does there exist $|\psi\rangle$ that satisfies $a_i |\psi\rangle = 0$? Also, it would be much better if there is a simpler way to obtain the ground state, by not solving $a_i |\psi\rangle = 0$.


*Can all the eigenvectors be obtained by procedure 2.?
It would be best if there is a unitary operator $U$ such that $U a_i U^{-1} = c_i$. In this case, the answers of 1. and 2. is automatically positive. Does such $U$ exist? If so, how can one find $U$?
Any help will be appreciated, although it partially answers my question.
 A: Indeed, it is true that one can find a unitary $U$ to relate $a$ to $c$. To find this unitary, it is most convenient to write everything in the Majorana basis: $c_i=\frac{1}{2}(\chi_i+i\eta_i)$, where $\chi_i^2=\eta_i^2=1$ are Majorana operators. Call them $\gamma_a, a=1,2,\dots,2n$ （not to be confused with the $\gamma_{ij}$ parameters in your Hamiltonian), then the Hamiltonian can be rewritten as $H=\frac{i}{2}\Gamma^T A\Gamma$, where $A$ is real skew-symmetric matrix, and $\Gamma$ is the column vector $(\gamma_1,\gamma_2,\dots,\gamma_{2n})^T$.  Then the Bogoliubov transform means that one finds a SO$(2n)$ transformation $V$ on $\gamma_a$'s: $\Gamma'=V\Gamma$ ($\Gamma'$ is basically equivalent to your $a$, written in Majorana basis), such that $A$ is "diagonalized" into the following form:
$$
VAV^T=\begin{pmatrix}
0 & \theta_1 & & &\\
-\theta_1 & 0 & & &\\
 &          &\ddots & &\\
 & & & 0 & \theta_n\\
& & & -\theta_n & 0
\end{pmatrix}
$$
Here $\theta_i>0$ (assuming no zero eigenvalues for simplicity). These $\theta_i$'s are basically the energy $E_i$'s, but let me use a different notation for clarity. The diagonalization can always be done. Then we define a unitary transformation
$$
U=\prod_{j=1}^n \left(\cos\frac{\theta_j}{2}+\sin\frac{\theta_j}{2}\Gamma'_{2j-1}\Gamma'_{2j}\right)
$$
This $U$ satisfies $U\Gamma U^\dagger = \Gamma'$. There is a more compact expression for $U$:
$$
U=\exp \left( \frac{1}{4} \Gamma^T D \Gamma\right)
$$
where the matrix $D$ satisfies $e^D=V$.
