Bogoliubov Transformation with Complex Hamiltonian Consider the following Hamiltonian:
$$H=\sum_k \begin{pmatrix}a_k^\dagger & b_k \end{pmatrix}
\begin{pmatrix}\omega_0 & \Omega f_k \\ \Omega f_k^* & \pm \omega_0\end{pmatrix} \begin{pmatrix}a_k \\\ b_k^\dagger\end{pmatrix}\tag{1}$$
for bosonic operators ($+$) or fermionic operators ($-$). The standard way to do  Bogoliubov transformations is to use the transformations:
$$M_{\text{boson}}=\begin{pmatrix} \cosh(\theta) & \sinh(\theta)\\ \sinh(\theta)&\cosh(\theta)\end{pmatrix},\quad M_{\text{fermion}}=\begin{pmatrix} \cos(\theta) & \sin(\theta)\\ -\sin(\theta)&\cos(\theta)\end{pmatrix}$$
However, in this case these won't work as they will give complex values of $\theta$, and to ensure that our (anti-)commutators remain intact we need $\theta$ to be real. 
Thus my question is: How do we generalize the Bogoliubov to solve problems of the form of (1)?
This question is based of this one: Bogoliubov transformation with a slight twist 
 A: There is always a bottom-line answer to this question: write the complex boson/fermion in terms of real boson/fermion ($a=a_R+i a_I$, etc), plug it in, and then diagonalize it by orthogonal matrices. This is probably the more natural way to do it for particle non-conserving systems.
If one insists on doing it in terms of complex boson/fermion, it's still possible, but many of the time annoying. This is because one (generically) also need to transform within the real and the imaginary part of the variables, which forces one to double the size of the matrix to include $( a,a^{\dagger},b,b^{\dagger})^T$ all together, like the Nambu spinor when one solves for superconductors' mean-field Hamiltonian. The annoying part is that one needs to take care of the redundancy in the matrix components.
A: There are two methods to tackling this problem:

*

*As pointed out in Yen-Ta Huang's answer and also in this Everett You's answer (EY16) to this related question we can split the creation and annihilation operators into a real and an imaginary part.

*As hinted at in (Capri, 2002; pg448) we can generalize the Bogoliubov transform to work with complex Hamiltonians.

Here I will do a simple example with the following fermionic Hamiltonian:
$$H=\varepsilon c_1^\dagger c_1+\varepsilon c_2^\dagger c_2+\lambda i(c_1^\dagger c_2^\dagger-c_2c_1)\tag{1}$$
Method 1
We let:
$$c_j=a_j+i b_j\quad \text{for}\quad j=1,2 \tag{2}$$
where $a_j^\dagger=a_j$ and $b_j^\dagger=b_j$. As shown in EY16 for $a_j$ and $b_j$ we have the following commutation relations
$$\{a_j,a_j\}=\{b_j,b_j\}=1$$
$$\{a_1,a_2\}=\{b_1,b_2\}=\{a_i,b_j\}=0$$
Thus subbing (2) into (1) we get that (after some algebra):
$$H=2i (\varepsilon a_1b_1+\varepsilon a_2 b_2+\lambda a_1 a_2-\lambda b_1 b_2)$$
$$=2i\begin{pmatrix} a_1 &b_2 \end{pmatrix}\begin{pmatrix} \varepsilon & \lambda \\ \lambda &\varepsilon \end{pmatrix}\begin{pmatrix} b_1 \\ a_2\end{pmatrix}$$
As explained in EY16 a Bogoliubov transformation of $a_j$ and $b_j$ is an orthogonal transformation in the case of fermions. Thus if we let:
$$\begin{pmatrix} b_1 \\ a_2\end{pmatrix}=\begin{pmatrix} \cos(\theta) & \sin(\theta)\\ -\sin(\theta) & \cos(\theta)\end{pmatrix} \begin{pmatrix} e_1 \\ d_2\end{pmatrix}$$
$$\begin{pmatrix} a_1 \\ b_2\end{pmatrix}=\begin{pmatrix} \cos(\theta) & \sin(\theta)\\ -\sin(\theta) & \cos(\theta)\end{pmatrix} \begin{pmatrix} d_1 \\ e_2\end{pmatrix}$$
with the new fermionic creation and annihilation operators being given by $f_j=d_j+ie_j$ with an appropriate choice of $\theta$ this will diagonalize the Hamiltonian
Method 2
In method 2 we simply generalize the Bogoliubov transformation. Consider the transformation:
$$f_j=u_jc_j+v_j c_j^\dagger$$
we are needing to enforce the conditions that:
$$\{f_i,f_j\}=0, \quad \{ f_i,f_j^\dagger\}=\delta_{ij}$$
If we do this we get that we need:
$$u_1v_2+u_2v_1=0\tag{3}$$
and
$$|u_j|^2+|v_j|^2=1\tag{4}$$
(4) implies that we have:
$$u_j=\cos(\theta_j) e^{i\phi_j^u}\quad v_j=\sin(\theta_j) e^{i\phi_j^v}$$
whilst with these (3) implies that:
$$\cos(\theta_1)\sin(\theta_2)=-\cos(\theta_2) \sin(\theta_1),\quad \phi_1^u+\phi_2^v=\phi_2^u+\phi_1^v$$
Putting these together the general Bololiubov transformation of fermionic operators is:
$$e^{i\tilde \phi_1} \begin{pmatrix}e^{i\tilde \phi_2} \cos(\theta_p) & e^{i\tilde \phi_3}\sin(\theta_p)\\ -e^{-i\tilde \phi_3}\sin(\theta_p) & e^{-i\tilde \phi_2}\cos(\theta_p) \end{pmatrix}$$
The standard method of the Bololiubov transformation can then be followed with this.
For reference the general Bololiubov transformation for bosons is (according to my calcuations:
$$e^{i\tilde \phi_1} \begin{pmatrix}e^{i\tilde \phi_2} \cosh(\theta_p) & e^{i\tilde \phi_3}\sinh(\theta_p)\\ e^{-i\tilde \phi_3}\sinh(\theta_p) & e^{-i\tilde \phi_2}\cosh(\theta_p) \end{pmatrix}$$
