Solution (source code)

= Solution

For each complex fermion define two <Majorana fermion operators>
$$
c_{2j-1}=a_j+a_j^\dagger,
\qquad
c_{2j}=-i(a_j-a_j^\dagger).
$$
They obey $c_r^\dagger=c_r$, $\{c_r,c_s\}=2\delta_{rs}$, and $a_j=(c_{2j-1}+ic_{2j})/2$. Substituting these inverse relations expands every hopping and pairing monomial as a bilinear in the $c_r$. Diagonal terms $c_r^2=1$ contribute only a constant.

For $r\ne s$, $c_rc_s$ is anti-Hermitian. Hermiticity of the original <quadratic fermion Hamiltonian> therefore makes its coefficient purely imaginary, so it can be written $iA_{rs}$ with $A_{rs}$ real. Since
$$
\sum_{rs}A_{rs}c_rc_s
=\frac12\sum_{rs}(A_{rs}-A_{sr})c_rc_s
+\frac12\operatorname{tr}A,
$$
the symmetric part again changes only the constant. Absorbing conventional factors into $A$ gives
$$
H=i\sum_{rs}A_{rs}c_rc_s+\text{constant},
\qquad A^T=-A,
\qquad A\in M_{2n}(\mathbb R).
$$