I define $H\sim\psi_{\alpha}^{\dagger}\xi_{\alpha\beta}\psi_{\beta}$,
where I forget the sum/integrals and all these boring staff. I also
define $\xi_{\alpha\beta}\equiv\boldsymbol{\xi}\cdot\boldsymbol{\sigma}=\xi_{0}+\xi_{x}\sigma_{x}+\xi_{y}\sigma_{y}+\xi_{z}\sigma_{z}$
to have the most generic one-body Hamiltonian written in a compact
form. The one body Hamiltonian then reads, in matrix notation
$$
H\sim\left(\begin{array}{cc}
\psi_{\uparrow}^{\dagger} & \psi{}_{\downarrow}^{\dagger}\end{array}\right)\left(\begin{array}{cc}
\xi_{0}+\xi_{z} & \xi_{x}-\mathbf{i}\xi_{y}\\
\xi_{x}+\mathbf{i}\xi_{y} & \xi_{0}-\xi_{z}
\end{array}\right)\left(\begin{array}{c}
\psi_{\uparrow}\\
\psi_{\downarrow}
\end{array}\right)
$$
as can be easily check.
One now wants to add the particle-hole double space (Nambu space).
One uses that (the anti-commutation relation)
$$
\psi_{\alpha}^{\dagger}\xi_{\alpha\beta}\psi_{\beta}=-\xi_{\alpha\beta}\psi_{\beta}\psi{}_{\alpha}^{\dagger}+\delta_{\alpha\beta}\xi_{\alpha\beta}=-\psi_{\beta}\left(\xi_{\alpha\beta}\right)^{T}\psi{}_{\alpha}^{\dagger}+\delta_{\alpha\beta}\xi_{\alpha\beta}
$$
and you get the unavoidable trace over the one-body energy. This nevertheless
renormalizes your energy in a standard way, and one usually drops
off this extra term. We thus get
$$
H\sim\dfrac{1}{2}\left(\begin{array}{cccc}
\psi{}_{\uparrow}^{\dagger} & \psi{}_{\downarrow}^{\dagger} & \psi_{\uparrow} & \psi_{\downarrow}\end{array}\right)\left(\begin{array}{cc}
\boldsymbol{\xi}\cdot\boldsymbol{\sigma} & -\mathbf{i}\sigma_{y}\Delta\\
\mathbf{i}\sigma_{y}\Delta & -\left(\boldsymbol{\xi}\cdot\boldsymbol{\sigma}\right)^{T}
\end{array}\right)\left(\begin{array}{c}
\psi_{\uparrow}\\
\psi_{\downarrow}\\
\psi{}_{\uparrow}^{\dagger}\\
\psi{}_{\downarrow}^{\dagger}
\end{array}\right)-\dfrac{1}{2}\text{Tr}\left\{ \boldsymbol{\xi}\cdot\boldsymbol{\sigma}\right\}
$$
in a mixed notation (block matrix in the middle, full vectors on the
edge). Note the only important thing here that $\left(\boldsymbol{\xi}\cdot\boldsymbol{\sigma}\right)^{T}=\left(\boldsymbol{\xi}\cdot\boldsymbol{\sigma}\right)^{\ast}$
and only the $\sigma_{y}$ component changes sign (look at your $\alpha p\sigma_{y}\tau_{z}$
term in the BdG Hamiltonian)
Your ordering convention is found by an obvious change of basis from
mine. Then you choose a representation for the tensor product and
you're done. One more time, you can not avoid the final trace term,
but most of the people forget to discuss it. It has almmost no role,
except when you want to describe some effects related to the phase
transition of superconductivity (for instance, to correctly write
the free energy, you need it).
One more thing: the Hamiltonian you gave is a bit famous at the moment
for hosting Majorana fermions. If you diagonalise the spin-part, you
end up with a $p$-wave effective superconductivity at low energy.
It is a standard procedure that you can find in any good quantum mechanical book. I suggest Messiah, since it's the best I know ! Everyone should check the Messiah before asking on this web-site ;-)
I call $\Phi$ the 4-spinor, it is defined as
$$\Phi\equiv\left(\begin{array}{c}
u\\
v
\end{array}\right)\;;\;\Phi^{\dagger}\equiv\left(\begin{array}{cc}
u^{\dagger} & v^{\dagger}\end{array}\right)\;;\; u\equiv\left(\begin{array}{c}
u_{\uparrow}\\
u_{\downarrow}
\end{array}\right)$$
in my favorite representation.
The Schrödinger equation then reads
$$\mathbb{H}\Phi=\mathbf{i}\hbar\dfrac{\partial\Phi}{\partial t}$$
with $\mathbb{H}$ identified with your matrix. When the problem is time-independent, the usual substitution $\Phi=e^{-\mathbf{i}\varepsilon t/\hbar}\Psi$ gives back your BdG equation. Then you calculate
$$-\mathbf{i}\hbar\dfrac{\partial\Phi^{\dagger}}{\partial t}=\Phi^{\dagger}\mathbb{H}^{\dagger}$$
for the adjoint spinor. You can verify that $\mathbb{H}^{\dagger}$ is the Hermitian conjugate of the $\mathbb{H}$ indeed (transpose of the complex conjugate).
Then you define the probability of particle as $\rho_{p}=\left|\Phi\right|=\Phi^{\dagger}\Phi$ and you calculate
$$\dfrac{\partial\rho_{p}}{\partial t}=\dfrac{\partial\Phi^{\dagger}}{\partial t}\Phi+\Phi^{\dagger}\dfrac{\partial\Phi}{\partial t}$$
and you inject the previous relations. By virtue of the conservation of the number of particle, you must obtain something like $\partial_{t}\rho_{p}+\nabla\cdot\mathbf{j}_{p}=0$, so you next have to identify the current-probability $\mathbf{j}_{p}$ as some divergence in your expression for $\partial_{t}\rho_{p}$. For a s-wave superconductor and when $H_{0}=p^{2}/2m=-\hbar^{2}\nabla^{2}/2m$, I found
$$\mathbf{j}_{p}=\dfrac{\hbar}{m}\left[\Im\left\{ u^{\dagger}\nabla u-v^{\dagger}\nabla v\right\} \right]$$
which sounds correct. This result does not depend on the spin-interaction. Nevertheless, you may choose a $H_{0}=v_{0}\left(\boldsymbol{\sigma\cdot p}\right)$ and you will obtain something different of course (in that case the identification of the divergence of the current is straightforward, so I let you finding the details in the Messiah ... section about the Dirac equation ... )
The trick is that charge are not conserved. Indeed, if you reproduce the same calculation for the density of charge $\rho_{e}=\Phi^{\dagger}\tau_{3}\Phi=u^{\dagger}u-v^{\dagger}v$ you should find something like
$$\dfrac{\partial\rho_{e}}{\partial t}+\nabla\cdot\mathbf{j}_{e}=\dfrac{2e}{\hbar}\left[\Delta^{\ast}\left(v^{\dagger}\sigma_{2}u\right)+\Delta\left(u^{\dagger}\sigma_{2}v\right)\right]$$
with $\mathbf{j}_{e}=\dfrac{e\hbar}{m}\Im\left\{ u\nabla u^{\dagger}+v\nabla v^{\dagger}\right\} $, and the result does not depend on the spin interaction neither. The first paper to discuss this problem was the seminal paper by Tinkham, Blonder and Klapwick (I can send it to you if you have some difficulties to obtain it, follow the links in my profiles and you may find my email address out, but first check around you: library, colleagues, ...)
Only for the spin current might you find that a non-diagonal $H_{0}$ will induce fancy effects. The spin-current might be defined from the spin-density $\rho_{s}^{j}=\Phi^{\dagger}\sigma_{j}\Phi=u^{\dagger}\sigma_{j}u+v^{\dagger}\sigma_{j}v$ and you calculate the time-derivative as above. The current is not conserved one more time. If you really insist I may find the time to write out a generic expression for the spin-current for a generic Hamiltonian.
Please ask for more details if the above explanations are not sufficient.
Important Edit: Contrary to what is said above, the charge current is obviously conserved in a superconductor, provided the gap parameter is treated self-consistently. See
for instance.
Best Answer
Great question. One should view the diagonalization of quadratic fermion models as arising from the representation theory of the spin group and its Lie algebra. Because your model is quadratic, the Hamiltonian lies in $\mathfrak{spin}_4$, which can be identified with $\mathfrak{so}_4$ by pushing forward via the universal covering map $$\gamma_*: \underbrace{\mathfrak{spin}_{2n}}_\text{quadratic models}\xrightarrow{\simeq} \underbrace{\mathfrak{so}_{2n}}_\text{skew-symmetric matrices}$$
$\mathfrak{so}_{2n}$ is the algebra of real antisymmetric matrices, which is precisely the statement of particle-hole symmetry. Therefore, "particle-hole symmetry" is already built into your model, not because of anything physical, but because it is quadratic. I would study general characteristics of the eigenvalue problem for $\mathfrak{so}_{n}$, where all of the properties you mention will become clearer and less obfuscated by physical data (see https://en.wikipedia.org/wiki/Skew-symmetric_matrix#Spectral_theory).