case the BdG equation has two eigenstates with $\epsilon\sim 0$.

Let us consider the $(2L)$-th and $(2L+1)$-th smallest eigenvalues, $\epsilon^-$ and $\epsilon^+$, which satisfy $-\epsilon^- = \epsilon^+$ because eigenvalues appear in pairs with same absolute value and opposite signs as mentioned above. We can only have an even number of vanishing eigenvalues, and if we have them $\epsilon^-$ and $\epsilon^+$ should be included in them; otherwise $\epsilon^- < 0 < \epsilon^+$.

For $\mu\ll\{-t,\min(\epsilon_{l})\}$, the number of fermions in the system is negligible. As $\mu$ is increased, $N_\sigma$ increases, with $N_\downarrow > N_\uparrow$ for $\Gamma > 0$. $\epsilon^+$ initially decreases linearly in $\mu$, reflecting the linear decrease of the required energy to add a single particle in the system. When it is closer to zero, however, the value of $\epsilon_\pm$ approaches more slowly to zero, especially for a smaller system. When we fit the decrease by a function of the shape $\exp{(-c\vert\mu-\mu_0\vert)}$, the exponent $c$ is roughly in proportion to the system size $L$. This also suggests that the modes corresponding to $\epsilon_\pm$ are spatially localized.

### Detection of end Majorana fermions

From the eigenvector of $M_\mathrm{BdG}$ corresponding to the eigenvalue $\epsilon^+$, $({\mbox{\boldmath $u$}}_\uparrow, {\mbox{\boldmath $u$}}_\downarrow, {\mbox{\boldmath $v$}}_\uparrow, {\mbox{\boldmath $v$}}_\downarrow)^\mathrm{T}$, we 