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.
First of all, at zero temperature not all of the electrons in the BCS ground sate (GS) form Cooper pairs. One way to think about this is that deep inside the Fermi surface, there is no available states for scattering events, including the phonon-mediated ones, and thus we wouldn't have the effective attraction to form the Cooper pairs.
This can also be seen from the coefficients $u_k$ and $v_k$ in the postulated BCS GS given in the original post. When $k$ is way below the Fermi surface, then $|\xi_k|>>\Delta$, with $\xi_k<0$. In this case, $|u_k|\to 0$ and $v_k\to 1$, i.e. the states far below the Fermi surface are almost fully occupied by electrons created by $c^{\dagger}_{k\uparrow}c^{\dagger}_{-k\downarrow}$. These electrons are created in pairs, but they are NOT Cooper pairs!!! The pair creation operator does not tell us anything regarding whether the pair is a Cooper or just a normal pair. Only when it is close to the Fermi surface, and thus possible to have phonon-mediated scattering, is the Cooper pair creation possible.
In summary, the BCS GS looks like the following: deep down inside the Fermi surface, it's just normal electrons. Close to the Fermi surface, it's mostly a superposition $\left(u_k+v_kc^{\dagger}_{k\uparrow}c^{\dagger}_{-k\downarrow}\right)|0\rangle$ of empty state, with amplitude $u_k$, and occupied Cooper pairs, with amplitude $v_k$.
EDIT:
After I submitted the above answer, I realized that I missed out the other part of your question, which is about the excitations of the BCS GS and here it goes:
You mentioned that "the excitations of the BCS state must be created or destroyed in pairs". This is not correct. The excitations of the BCS GS are fundamentally different from Cooper pairs or the normal electrons. Instead, it is the quasi-particle created by
$$\gamma^{\dagger}_{p\uparrow}=u^*_pc^{\dagger}_{p\uparrow}-v^*_pc_{-p\downarrow}$$
Not surprisingly, this is just the Bogolyubov transformation used to diagonalize the BCS mean field hamiltonian and this why it is the legitimate excitations of the GS. The corresponding eigenvalue is exactly
$$E_k=\sqrt{\xi_k^2+\Delta^2}$$
This is where the gap $\Delta$ to the BCS excitation comes from. This kind of excitation is slightly less intuitive compared to objects like single electrons or Cooper pairs. It has spin-1/2, but it doesn't carry well-defined charge, whereas a Cooper pair has spin 0 and charge $2e$.
Also, You can easily check that
$$\gamma_{p\uparrow}|BCS\rangle=\gamma_{-p\downarrow}|BCS\rangle=0$$
which makes perfect sense since the annihilation operator for the excitation annihilates the GS.
Best Answer
The localization of Majorana zero modes has a well-defined meaning: consider a Kitaev chain with two ends. Because of the zero modes, there are two nearly degenerate ground states, let us call them $|0\rangle$ and $|1\rangle$, which have opposite fermion parities. They are localized as "single-particle wavefunctions" in the following sense: if one computes the matrix element $\langle 1|c^\dagger(x)|0\rangle$ where $c^\dagger(x)$ is the creation operator for fermions, the result is an exponentially decaying function of $x$ away from the edge. This definition works even when the system is interacting. Intuitively it means that the weight for creating a single fermion excitation is localized near the edge, and in the bulk there is a finite gap to the single particle excitations.