define the averaged separation from the system center $$\sqrt{\langle x^2 \rangle} \equiv
\sqrt{
\left(\sum_{\sigma,l} (l-l_c)^2 \vert{\mbox{\boldmath $v$}}_{\sigma,l}\vert^2\right)
/
\left(\sum_{\sigma,l} \vert{\mbox{\boldmath $v$}}_{\sigma,l}\vert^2\right)
}.$$ If the modes localize to the system ends, $\sqrt{\langle x^2 \rangle}$ are close to $l_c$, the maximum value it can take. Note that, when $\epsilon^-$ and $\epsilon^+$ are numerically degenerate ($\epsilon^+ \lesssim 10^{-12}$ in our work), any linear combination of the two eigenvectors that correspond to these eigenvalues would be obtained as the eigenvector corresponding to $\epsilon^+$, so within our BdG calculation we do not directly observe eigenmodes that are localized at only one of the ends of the system.

In the density-matrix renormalization group (DMRG) simulations of the same model, [@Stoudenmire2011; @Tezuka2012] however, a pair of Majorana modes localized at each end of the system have been obtained. We believe that the DMRG calculation, with limitation on the entanglement entropy between the two ends of the system, automatically chooses less entangled degenerate ground states, which are connected to each other by operating either of the two Majorana operators, for the subspaces with even and odd numbers of fermions. We find that the localization of ${\mbox{\boldmath $u$}}_{\sigma,l}$ and ${\mbox{\boldmath $v$}}_{\sigma,l}$, in terms of $\sqrt{\langle x^2\rangle}$, precisely corresponds to 