the single-particle components of the Hamiltonian, $$\begin{aligned}
\hat{H_{\sigma\sigma}} = -\frac{t}{2}\sum_{l=0}^{L-2}\sum_{\sigma=\uparrow,\downarrow}
(\hat c_{\sigma,l}^\dag \hat c_{\sigma,l+1} + \mathrm{h.c.})\nonumber\\
+\sum_{l=0}^{L-1} \sum_{\sigma=\uparrow,\downarrow}(t-\mu+\epsilon_{l}
+(-)^\sigma (\hat{{\mbox{\boldmath $B$}}}\cdot\hat{{\mbox{\boldmath $z$}}})\Gamma
) \hat n_{\sigma,l},\end{aligned}$$ $$\begin{aligned}
\hat{H_{\downarrow\uparrow}} &=& \hat{H_{\uparrow\downarrow}}^\dag =
\frac{\alpha}{2} \sum_{l=0}^{L-2} 
(\hat c_{\downarrow,l}^\dag \hat c_{\uparrow,l+1}
-\hat c_{\downarrow,l+1}^\dag \hat c_{\uparrow,l})\\
&-&((\hat{{\mbox{\boldmath $B$}}}\cdot\hat{{\mbox{\boldmath $x$}}}) + i(\hat{{\mbox{\boldmath $B$}}}\cdot\hat{{\mbox{\boldmath $y$}}}))\Gamma
\sum_{l=0}^{L-1}\hat c_{\downarrow,l}^\dag \hat c_{\uparrow,l},\end{aligned}$$ in which $(-)^\sigma = \delta_{\sigma\downarrow}-\delta_{\sigma\uparrow}$. In the following we call the $4L$-dimensional matrix in the left hand side of eq. (\[eqn:BdG\]) $M_\mathrm{BdG}$.

We work in the limit of low temperature $T\rightarrow 0$. We obtain the particle distribution by $$\begin{aligned}
n_{\sigma,l} &=& \sum_q\left[f(\epsilon_q) |u_{\sigma,l}|^2 + f(-\epsilon_q) |v_{\sigma,l}|^2\right],
\label{eqn:density}\end{aligned}$$ in which $f(\epsilon)\equiv (1+\mathrm{e}^{\epsilon/(k_\mathrm{B}T)})
\rightarrow \Theta(-\epsilon)$ ($T\rightarrow 0$) is the Fermi distribution function, with $\Theta(x)$ being the step function. The total number of fermions with spin $\sigma$ is obtained as $N_\sigma = \sum_l n_{\sigma,l}$.

The sum in eq. (\[eqn:density\]) is taken over all $q$. The matrix $M_\mathrm{BdG}$ is Hermitian, and has pairs of positive and negative eigenvalues with equal absolute values, because if $M_\mathrm{BdG}({\mbox{\boldmath $u$}}_\uparrow, {\mbox{\boldmath $u$}}_\downarrow, {\mbox{\boldmath $v$}}_\uparrow, {\mbox{\boldmath $v$}}_\downarrow)^\mathrm{T}
= \epsilon({\mbox{\boldmath $u$}}_\uparrow, {\mbox{\boldmath $u$}}_\downarrow, {\mbox{\boldmath $v$}}_\uparrow, {\mbox{\boldmath $v$}}_\downarrow)^\mathrm{T}$, we have $M_\mathrm{BdG}({\mbox{\boldmath $v$}}_\uparrow, {\mbox{\boldmath $v$}}_\downarrow, {\mbox{\boldmath $u$}}_\uparrow, {\mbox{\boldmath $u$}}_\downarrow)^\mathrm{T}
= -\epsilon({\mbox{\boldmath $v$}}_\uparrow, {\mbox{\boldmath $v$}}_\downarrow, {\mbox{\boldmath $u$}}_\uparrow, {\mbox{\boldmath $u$}}_\downarrow)^\mathrm{T}.$ For a positive 