Stability of Majorana Edge Zero Modes against Interactions

Tohru Koma

Department of Physics, Gakushuin University (retired), Mejiro, Toshima-ku,
Tokyo 171-8588, JAPAN

Abstract: We study an interacting Majorana chain with an open boundary
condition. In the case without interactions, the system shows a prototypical
Majorana edge zero mode in the sector of the ground state with a spectral gap
above the sector. We prove that both of the Majorana edge zero mode and the
non-vanishing spectral gap are stable against generic weak interactions whose
fermion parity is even. We also deal with the corresponding ladder systems,
and discuss the $\mathbb{Z}_{2}$ index for the Majorana edge zero modes.

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
  3. 3 Majorana-Hubbard chains
  4. 4 Majorana-Hubbard ladders and edge $\mathbb{Z}_{2}$ index
  5. 5 Explicit construction of Majorana edge zero modes
    1. 5.1 Exactly solvable cases
    2. 5.2 A series expansion method
  6. A Stability of the spectral gap
  7. B Proof of the existence of $\gamma$ satisfying $[H,\gamma]=0$
  8. C Locality of Majorana edge zero mode

##  1 Introduction

Majorana fermions have been in the focus of many researchers, so far. In
particular, Majorana edge zero modes [1] are the most typical ones in
condensed matter physics. Although they exhibit the peculiarity, i.e., one
half of a usual complex fermion, they are indeed realized at open boundaries
of tight binding models. However, the effects of interactions are still
unclear. Fidkowski and Kitaev [2] discussed the effects of interactions for
Majorana ladders from the point of view of topological classification. They
concluded that the topological invariant $\mathbb{Z}$ for the free fermion
classification is broken to $\mathbb{Z}_{8}$, due to the introduction of
generic interactions, under a certain assumption on the symmetry of the
Hamiltonian. Katsura, Schuricht and Takahashi [3] obtained the explicit forms
of the left and right Majorana edge zero modes for a one-parameter family of
interacting Majorana chains which are derived from a certain XYZ spin chain in
a magnetic field, via Jordan-Wigner transformation. Fendley [4] obtained the
zero mode operator for the exactly solvable XYZ chain. Moreover, the
parafermionic edge zero modes which are a generalization of the Majorana edge
zero mode were treated in [5, 6, 7]. For those zero modes, the dynamics [8, 9]
and the thermal effect [10] have been studied as well. See also [11, 12, 13,
17, 14, 15, 16, 18] for other several aspects of parafermions. However,
mathematical treatments of Majorana edge zero modes for interacting systems
are still rare.

In the present paper, we study an interacting Majorana chain and ladders with
an open boundary. In the case without interactions, these systems show
Majorana edge zero modes in the sector of the ground state with a spectral gap
above the sector. For the single chain, we prove that both of the Majorana
edge zero mode and the non-vanishing spectral gap are stable against generic
weak interactions that have even fermion parity. The results are summarized as
Theorem 3.1 in Sec. 3 below.

In Sec. 4, we deal with the ladders with $L$ legs for $L\geq 2$. Clearly, for
the $L$ copies of the noninteracting Majorana chains, there appear $L$
Majorana edge zero modes at the left edge. When $L$ is even, certain
interchain interactions lift all the degeneracy for these Majorana edge zero
modes. On the other hand, when $L$ is odd, there remains an unpaired Majorana
edge zero mode for any weak interaction that has even fermion parity. Namely,
the evenoddness of the Majorana edge zero modes is invariant against weak
perturbations with even parity. This implies that for adding Majorana edge
zero modes, the corresponding index is given by the additive group
$\mathbb{Z}_{2}$ in the sense of the stability against generic weak
perturbations.

Acknowledgements: I would like to thank Tomonari Mizoguchi for many useful
discussions. I also thank Hosho Katsura for many helpful comments.

##  2 Preliminaries

Consider first a prototypical Majorana chain whose Hamiltonian is given by [1,
2]

$H_{0}:=i\kappa\sum_{\ell=1}^{N-1}c_{2\ell-1}c_{2\ell}+i\sum_{\ell=1}^{N-1}c_{2\ell}c_{2\ell+1},$
(2.1)

where $\kappa$ is a real parameter, and $c_{i}$ is the Majorana fermion
operator at the site $i\in\\{1,2,\ldots,2N\\}$ with the length $2N$ of the
chain. The Majorana fermion operators satisfy $c_{i}^{\dagger}=c_{i}$ and obey
the anticommutation relations,

$\\{c_{i},c_{j}\\}=2\delta_{i,j}.$

We assume $0\leq|\kappa|<1$ so that there appears a spectral gap above the
ground state of $H_{0}$. Since the Hamiltonian $H_{0}$ does not contain
$c_{2N}$, $c_{2N}$ commutes with $H_{0}$. Namely, the right edge mode is
trivially given by $c_{2N}$.

In order to obtain the Majorana edge zero mode $\gamma_{0}$ at the left edge
for the Hamiltonian $H_{0}$, we set

$\gamma_{0}=\sum_{m=1}^{N}u_{m}c_{2m-1}.$ (2.2)

The coefficients $u_{m}$ are determined so as to satisfy the commutation
relation $[H_{0},\gamma_{0}]=0$. Clearly, when this relation is fulfilled, the
state $\gamma_{0}|0\rangle$ has the same energy for the ground state
$|0\rangle$, i.e., it is the Majorana zero mode. Note that

$[H_{0},\gamma_{0}]=i\kappa\sum_{\ell=1}^{N-1}\sum_{m=1}^{N}[c_{2\ell-1}c_{2\ell},c_{2m-1}]u_{m}+i\sum_{\ell=1}^{N-1}\sum_{m=1}^{N}[c_{2\ell}c_{2\ell+1},c_{2m-1}]u_{m}.$

The commutators in the summands are computed as follows:

$c_{2\ell-1}c_{2\ell}c_{2m-1}-c_{2m-1}c_{2\ell-1}c_{2\ell}=-2\delta_{\ell,m}c_{2\ell},$

and

$c_{2\ell}c_{2\ell+1}c_{2m-1}-c_{2m-1}c_{2\ell}c_{2\ell+1}=2\delta_{\ell+1,m}c_{2\ell}.$

Substituting these into the above the right-hand side, one has

$[H_{0},\gamma_{0}]=-2i\kappa\sum_{\ell=1}^{N-1}c_{2\ell}u_{\ell}+2i\sum_{\ell=1}^{N-1}c_{2\ell}u_{\ell+1}=0.$

This implies $u_{\ell+1}=\kappa u_{\ell}$ for $\ell=1,2,\ldots,N-1$.
Therefore, the coefficients $u_{\ell}$ are given by
$u_{\ell}=\kappa^{\ell-1}u_{1}$ with a constant $u_{1}$. As a result, one
obtains

$\gamma_{0}=\sum_{\ell=1}^{N}\kappa^{\ell-1}c_{2\ell-1},$ (2.3)

where we have set $u_{1}=1$ for simplicity. Since $|\kappa|<1$ from the
assumption, the mode is localized at the left edge of the chain. The
normalized operator is given by

$\hat{\gamma}_{0}:=\sqrt{\frac{1-\kappa^{2}}{1-\kappa^{2N}}}\gamma_{0},$ (2.4)

which satisfies $\hat{\gamma}_{0}^{2}=1$.

Let us consider the situation when the Hamiltonian $H_{0}$ is diagonalized in
terms of the usual complex fermions. The two Majorana edge modes are also
treated in terms of the corresponding complex single fermion. For this
purpose, we introduce a fermion operator,

$a_{\rm E}:=\frac{1}{2}(\hat{\gamma}_{0}+ic_{2N}).$ (2.5)

One can check $\\{a_{\rm E},a_{\rm E}^{\dagger}\\}=1$, and $a_{\rm E}^{2}=0$.
Let $|0\rangle$ be the ground state of $H_{0}$ which satisfies $a_{\rm
E}|0\rangle=0$. Then, the other ground state is given by $a_{\rm
E}^{\dagger}|0\rangle$. Clearly, one has $\langle 0|a_{\rm
E}^{\dagger}|0\rangle=0$. This implies that the sector of the ground state of
$H_{0}$ is spanned by $|0\rangle$ and $a_{\rm E}^{\dagger}|0\rangle$.

##  3 Majorana-Hubbard chains

Next we introduce interactions $V$ for the free Hamiltonian $H_{0}$ of (2.1).
More specifically, we consider the Majorana-Hubbard chains [19] whose generic
form of the Hamiltonians is given by

$H_{g}=H_{0}+gV$ (3.1)

with the interaction Hamiltonian,

$V=\sum_{i<j<k<\ell}K_{ijk\ell}c_{i}c_{j}c_{k}c_{\ell},$ (3.2)

where $g$ and $K_{ijk\ell}$ are real parameters. For simplicity, we assume
that the interaction $V$ is of finite range. In our approach, we can treat
more general interactions with even parity such as many-body interactions
consisting of six or eight Majorana fermion operators. We also assume that $V$
dose not contain $c_{2N}$. Therefore, $c_{2N}$ is still the right edge mode.

To begin with, we present a general argument about the Majorana edge zero
modes. When $g=0$ for the Hamiltonian $H_{g}=H_{0}+gV$, the sector of the
ground state is spanned by the ground state $|0\rangle$ and $a_{\rm
E}^{\dagger}|0\rangle$ as mentioned at the end of the preceding section. We
write $P(0)$ for the projection onto the sector of the ground state of
$H_{0}$.

In Appendix A, we prove that the spectral gap above the sector of the ground
state is stable against weak perturbations [20, 21, 22] under the assumption
that the sector of the ground state is two-fold degenerate. This degeneracy
will be proved below in this section. Namely, there exists a positive constant
$g_{\rm max}$ such that for $|g|\leq g_{\rm max}$, a non-vanishing spectral
gap above the sector of the ground state exists. Therefore, for $|g|\leq
g_{\rm max}$, there exists a unitary operator $U(g)$ such that the spectral
projection $P(g)$ onto the low energy sector for the Hamiltonian $H_{g}$ is
given by [23, 24]

$P(g)=U(g)P(0)U(g)^{\dagger}.$ (3.3)

These low energy states are given by

$|0,g\rangle:=U(g)|0\rangle\quad\mbox{and}\quad a_{\rm
E}(g)^{\dagger}|0,g\rangle,$

where have written $a_{\rm E}(g):=U(g)a_{\rm E}U(g)^{\dagger}$. In fact, from
the parity conservation, these two states are an eigenstate of $H_{g}$ and
orthogonal to each other. However, it is not clear whether or not the energies
are degenerate.

Let us consider the eigenvalue equation,

$H_{g}|0,g\rangle=E_{0}(g)|0,g\rangle,$ (3.4)

where $E_{0}(g)$ is the eigenenergy of $|0,g\rangle$. By multiplying the both
sides by $c_{2N}$ and using $[H_{g},c_{2N}]=0$, one has

$H_{g}c_{2N}|0,g\rangle=E_{0}(g)c_{2N}|0,g\rangle.$

Therefore, the state $c_{2N}|0,g\rangle$ is an eigenstate of $H_{g}$ with the
same eigenvalue $E_{0}(g)$. However, it has the opposite fermion parity to
that of $|0,g\rangle$. Since such a state is unique and must be equal to
$a_{\rm E}(g)^{\dagger}|0,g\rangle$. As a result, we have

$H_{g}a_{\rm E}(g)^{\dagger}|0,g\rangle=E_{0}(g)a_{\rm
E}(g)^{\dagger}|0,g\rangle.$

Thus, the two states are degenerate. Further, this can be rewritten as

$[H_{g},a_{\rm E}(g)^{\dagger}]|0,g\rangle=0.$

The operator $a_{\rm E}(g)^{\dagger}$ in the commutator can be written

$a_{\rm E}(g)^{\dagger}=U(g)a_{\rm
E}^{\dagger}U(g)^{\dagger}=\frac{1}{2}U(g)(\hat{\gamma}_{0}-ic_{2N})U(g)^{\dagger}=\frac{1}{2}U(g)\hat{\gamma}_{0}U(g)^{\dagger}-\frac{i}{2}c_{2N}$

from the definitions and the assumption $[H_{g},c_{2N}]=0$ on the Hamiltonian
$H_{g}$. These yield

$[H_{g},U(g)\hat{\gamma}_{0}U(g)^{\dagger}]|0,g\rangle=0.$ (3.5)

Thus, the operator $U(g)\hat{\gamma}_{0}U(g)^{\dagger}$ creates the Majorana
edge zero mode above the ground state $|0,g\rangle$ of the Hamiltonian
$H_{g}$. Namely, the Majorana edge zero mode is stable against the
perturbations.

The above result (3.5) suggests that the operator relation
$[H_{g},U(g)\hat{\gamma}_{0}U(g)^{\dagger}]=0$ holds without acting on the
state $|0,g\rangle$. In fact, Goldstein and Chamon [25] proved that there
exists a non-trivial operator $\gamma$ with odd fermion parity such that it
satisfies $[H,\gamma]=0$ for a generic Hamiltonian $H$ which consists of an
odd number of Majorana fermion operators. If the operator $\gamma$ satisfies
$[H,\gamma]=0$, then $[H,\gamma^{\dagger}]=0$. This yields
$[H,\gamma+\gamma^{\dagger}]=0$. Therefore, one can take $\gamma$ to satisfy
$\gamma^{\dagger}=\gamma$. The proof of the existence of $\gamma$ is given in
Appendix B.

Since the present Hamiltonian $H_{g}$ consists of the odd number of the
Majorana fermions, there exists an operator $\hat{\gamma}_{g}$ such that
$[H_{g},\hat{\gamma}_{g}]=0$. Clearly, this operator $\hat{\gamma}_{g}$ also
creates the Majorana edge zero mode at the left edge. Therefore, we have

$U(g)\hat{\gamma}_{0}U(g)^{\dagger}|0,g\rangle=\hat{\gamma}_{g}|0,g\rangle$
(3.6)

for a small coupling constant $g$. The locality of the Majorana edge mode is
proved in Appendix C. Thus, our main results are summarized as follows:

###### Theorem 3.1

There exists a positive constant $g_{\rm max}$ such that for $|g|\leq g_{\rm
max}$, the ground state of the Hamiltonian $H_{g}$ of (3.1) is two-fold
degenerate with a non-vanishing spectral gap above the sector of the ground
state. One of the two ground-state vectors can be taken to be a Majorana edge
mode which is localized near the left edge.

##  4 Majorana-Hubbard ladders and edge $\mathbb{Z}_{2}$ index

In this section, we deal with Majorana ladders [2] which are a system of $L$
Majorana chains with interchain interactions. Here, $L\geq 2$ is the number of
the legs of the ladder. As we will show below, we can expect that the
evenoddness of the Majorana left edge modes is stable against generic weak
perturbations with even fermion parity.

Consider first the case with $L=2$. The Hamiltonian is given by

$H_{0}^{(2)}:=i\sum_{j=1}^{2}\sum_{\ell=1}^{N-1}(\kappa
c_{2\ell-1,j}c_{2\ell,j}+c_{2\ell,j}c_{2\ell+1,j}),$ (4.1)

which is the sum of the two copies of the Majorana chain of (2.1). The two
copies of the set of the Majorana operators satisfy the anti-commutation
relations,

$\\{c_{m,i},c_{n,j}\\}=2\delta_{m,n}\delta_{i,j},\quad\mbox{for
}m,n=1,2,\ldots,2N-1;\ i,j=1,2.$ (4.2)

The two operators, $c_{2N,j}$ $j=1,2$, are the trivial right edge modes.
Clearly, the two left Majorana edge modes are given by

$\gamma_{0,j}=\sum_{\ell=1}^{N}\kappa^{\ell-1}c_{2\ell-1,j}\quad\mbox{for
}j=1,2.$ (4.3)

The normalized operators are given by

$\hat{\gamma}_{0,j}:=\sqrt{\frac{1-\kappa^{2}}{1-\kappa^{2N}}}\gamma_{0,j},$
(4.4)

for $j=1,2$, which satisfy $(\hat{\gamma}_{0,j})^{2}=1$ for $j=1,2$.

Let us consider a very simple interchain interaction,

$V=i\hat{\gamma}_{0,1}\hat{\gamma}_{0,2},$ (4.5)

which is instructive to understand the difference between the ladder with two
legs and the single chain. Clearly, this interaction has even fermion parity
and is of short range.

In order to diagonalize the total Hamiltonian $H_{g}=H_{0}^{(2)}+gV$ in terms
of usual complex fermions, we introduce a complex fermion operator,

$\chi:=\frac{1}{2}(\hat{\gamma}_{0,1}+i\hat{\gamma}_{0,2}).$ (4.6)

Then, one has

$\\{\chi,\chi^{\dagger}\\}=1$ (4.7)

and

$i\hat{\gamma}_{0,1}\hat{\gamma}_{0,2}=2\chi^{\dagger}\chi-1.$ (4.8)

This implies that the two eigenvalues of $V$ are given by $\pm 1$. Similarly,
for the right edge modes, we can introduce the corresponding complex fermion
$(c_{2N,1}+ic_{2N,2})/2$. Therefore, the degeneracy of the left zero modes is
lifted by the perturbation $V$ even for any small coupling $g$.

Clearly, this result can be extended to the generic ladders with even number
legs. On the other hand, for odd number legs, there must remain a unpaired
Majorana fermion. Thus, the evenoddness of Majorana edge modes can be expected
to be the $\mathbb{Z}_{2}$ invariant under weak perturbations with even
parity.

##  5 Explicit construction of Majorana edge zero modes

###  5.1 Exactly solvable cases

Clearly, from the result $[H_{0},\gamma_{0}]=0$ of Sec. 2, one has
$[H_{g},\gamma_{0}]=[H_{0},\gamma_{0}]+g[V,\gamma_{0}]=g[V,\gamma_{0}]$.
Therefore, if $[V,\gamma_{0}]=0$, then $\gamma_{0}$ is the desired solution.

Since the operator $\gamma_{0}$ consists of $c_{i}$ at only the odd site
$i=2\ell-1$, $\gamma_{0}$ commutes with an interaction $V$ which consists of
$c_{j}$ at only the even site $j=2k$. Namely, the interaction $V$ has the
form,

$V=\sum_{i<j<k<\ell}K_{2i,2j,2k,2\ell}c_{2i}c_{2j}c_{2k}c_{2\ell}.$

In a sense, it may be said that this is a trivial case.

In order to present other examples, we introduce an operator,

$b_{\ell}:=\kappa c_{2\ell-1}-c_{2\ell+1}.$

Note that

$\displaystyle\\{b_{\ell},c_{2\ell-1}+\kappa c_{2\ell+1}\\}$ $\displaystyle=$
$\displaystyle(\kappa c_{2\ell-1}-c_{2\ell+1})(c_{2\ell-1}+\kappa
c_{2\ell+1})+(c_{2\ell-1}+\kappa c_{2\ell+1})(\kappa c_{2\ell-1}-c_{2\ell+1})$
$\displaystyle=$
$\displaystyle\kappa+\kappa^{2}c_{2\ell-1}c_{2\ell+1}-c_{2\ell+1}c_{2\ell-1}-\kappa+\kappa-
c_{2\ell-1}c_{2\ell+1}+\kappa^{2}c_{2\ell+1}c_{2\ell-1}-\kappa$
$\displaystyle=$ $\displaystyle 0.$

Using this identity, we have

$\displaystyle\\{b_{\ell},\gamma_{0}\\}$ $\displaystyle=$
$\displaystyle\sum_{m=1}\\{b_{\ell},\kappa^{m-1}c_{2m-1}\\}$ $\displaystyle=$
$\displaystyle\\{b_{\ell},\kappa^{\ell-1}c_{2\ell-1}+\kappa^{\ell}c_{2\ell+1}\\}$
$\displaystyle=$ $\displaystyle\kappa^{\ell-1}\\{b_{\ell},c_{2\ell-1}+\kappa
c_{2\ell+1}\\}=0.$

As an interaction $V$, let us consider

$V=\sum_{\ell=1}b_{\ell}c_{2\ell}c_{2\ell+2}c_{2\ell+4}=\sum_{\ell=1}(\kappa
c_{2\ell-1}-c_{2\ell+1})c_{2\ell}c_{2\ell+2}c_{2\ell+4}.$

From the above observations, we have $[\gamma_{0},V]=0$.

Similarly, we can consider

$V=\sum_{\ell=1}[b_{\ell}b_{\ell+1}c_{2\ell}c_{2\ell+2}+\mbox{h.c.}].$

This also satisfies $[\gamma_{0},V]=0$. In passing, we note that

$\displaystyle b_{\ell}b_{\ell+1}c_{2\ell}c_{2\ell+2}+\mbox{h.c.}$
$\displaystyle=$ $\displaystyle
b_{\ell}b_{\ell+1}c_{2\ell}c_{2\ell+2}+c_{2\ell+2}c_{2\ell}b_{\ell+1}b_{\ell}$
$\displaystyle=$
$\displaystyle(b_{\ell}b_{\ell+1}-b_{\ell+1}b_{\ell})c_{2\ell}c_{2\ell+2}$
$\displaystyle=$ $\displaystyle 2(\kappa^{2}c_{2\ell-1}c_{2\ell+1}-\kappa
c_{2\ell-1}c_{2\ell+3}+c_{2\ell+1}c_{2\ell+3})c_{2\ell}c_{2\ell+2},$

where we have used

$\displaystyle b_{\ell}b_{\ell+1}-b_{\ell+1}b_{\ell}$ $\displaystyle=$
$\displaystyle(\kappa c_{2\ell-1}-c_{2\ell+1})(\kappa
c_{2\ell+1}-c_{2\ell+3})-(\kappa c_{2\ell+1}-c_{2\ell+3})(\kappa
c_{2\ell-1}-c_{2\ell+1})$ $\displaystyle=$ $\displaystyle
2\kappa^{2}c_{2\ell-1}c_{2\ell+1}-2\kappa
c_{2\ell-1}c_{2\ell+3}+2c_{2\ell+1}c_{2\ell+3}.$

###  5.2 A series expansion method

When the coupling constant $g$ of the interaction is weak compared to the
spectral gap above the ground-state sector of the unperturbed Hamiltonian
$H_{0}$, a series expansion method can be expected to work well [21, 26].

In order to find $\gamma$ satisfying $[H_{0}+gV,\gamma]=0$ for a given
interaction $V$ with a weak coupling $g$, we expand $\gamma$ into power series
of $g$ as follows:

$\gamma=\sum_{n=0}g^{n}\gamma_{n}.$

Substituting this into $[H_{0},\gamma]=-g[V,\gamma]$, one has

$\sum_{n=0}g^{n}[H_{0},\gamma_{n}]=-\sum_{n=0}g^{n+1}[V,\gamma_{n}].$

The zeroth order is given by

$[H_{0},\gamma_{0}]=0.$

This is already fulfilled. The $n$-th order is given by

$[H_{0},\gamma_{n}]=-[V,\gamma_{n-1}]$ (5.1)

for $n=1,2,\ldots$.

As a first demonstration, we consider a very simple case
$V=c_{1}c_{2}c_{3}c_{4}$. It consists of only a single term, but it can be
expected to represent one of essential, effective interactions because the
Majorana edge mode is localized at the edge.

Let us consider the right-hand side of (5.1) in the case of $n=1$. Note that

$\displaystyle[V,\gamma_{0}]$ $\displaystyle=$
$\displaystyle[c_{1}c_{2}c_{3}c_{4},c_{1}]+\kappa[c_{1}c_{2}c_{3}c_{4},c_{3}]$
$\displaystyle=$ $\displaystyle-2c_{2}c_{3}c_{4}-2\kappa c_{1}c_{2}c_{4}.$

Therefore, in order to find $\gamma_{1}$, it is sufficient to solve

$[H_{0},\gamma_{1}]=2c_{2}c_{3}c_{4}+2\kappa c_{1}c_{2}c_{4}$ (5.2)

from (5.1). Before searching for a solution $\gamma_{1}$, we remark that there
are many solutions even restricting $\gamma_{1}$ to a sum of terms of three
Majorana fermions. In fact, the quantity $(\gamma_{0}H_{0}+H_{0}\gamma_{0})/2$
commutes with $H_{0}$, and is a sum of terms of three Majorana fermions.

Since $\gamma_{0}$ has the form of the power series of $\kappa$, we similarly
set

$\gamma_{1}=\sum_{j=0}\kappa^{j}\gamma_{1}^{(j)}.$ (5.3)

Note that

$\displaystyle[H_{0},\gamma_{1}]$ $\displaystyle=$
$\displaystyle\kappa[H_{0,1},\gamma_{1}]+[H_{0,0},\gamma_{1}]$ (5.4)
$\displaystyle=$
$\displaystyle\sum_{j=0}\kappa^{j+1}[H_{0,1},\gamma_{1}^{(j)}]+\sum_{j=0}\kappa^{j}[H_{0,0},\gamma_{1}^{(j)}],$

where we have written

$H_{0}=\kappa H_{0,1}+H_{0,0}$

with

$H_{0,1}:=i\sum_{\ell=1}c_{2\ell-1}c_{2\ell}$

and

$H_{0,0}:=i\sum_{\ell=1}c_{2\ell}c_{2\ell+1}.$

Therefore, from (5.2) and (5.4), we have

$[H_{0,0},\gamma_{1}^{(0)}]=2c_{2}c_{3}c_{4},$ (5.5)
$[H_{0,1},\gamma_{1}^{(0)}]+[H_{0,0},\gamma_{1}^{(1)}]=2c_{1}c_{2}c_{4},$
(5.6)

and

$[H_{0,1},\gamma_{1}^{(j)}]+[H_{0,0},\gamma_{1}^{(j+1)}]=0\mbox{\ \ for \
}j=1,2,\ldots.$ (5.7)

In order to find the solution $\gamma_{1}^{(0)}$ of (5.5), we try
$\gamma_{1}^{(0)}=-ic_{2}c_{3}c_{5}$ as a candidate. Actually, we have

$\displaystyle[H_{0,0},\gamma_{1}^{(0)}]$ $\displaystyle=$
$\displaystyle\sum_{\ell=1}[c_{2\ell}c_{2\ell+1},c_{2}c_{3}c_{5}]$
$\displaystyle=$
$\displaystyle[c_{2}c_{3},c_{2}c_{3}c_{5}]+[c_{4}c_{5},c_{2}c_{3}c_{5}]$
$\displaystyle=$ $\displaystyle 2c_{2}c_{3}c_{4}.$

Thus, $\gamma_{1}^{(0)}=-ic_{2}c_{3}c_{5}$ is one of the solutions of (5.5).

Next, in order to obtain $\gamma_{1}^{(1)}$, let us compute the first
commutator $[H_{0,1},\gamma_{1}^{(0)}]$ in the left-hand side of (5.6). Note
that

$\displaystyle[H_{0,1},\gamma_{1}^{(0)}]$ $\displaystyle=$
$\displaystyle\sum_{\ell=1}[c_{2\ell-1}c_{2\ell},c_{2}c_{3}c_{5}]$
$\displaystyle=$
$\displaystyle[c_{1}c_{2},c_{2}c_{3}c_{5}]+[c_{3}c_{4},c_{2}c_{3}c_{5}]+[c_{5}c_{6},c_{2}c_{3}c_{5}]$
$\displaystyle=$ $\displaystyle
2c_{1}c_{3}c_{5}-2c_{2}c_{4}c_{5}-2c_{2}c_{3}c_{6}.$

Substituting this into the left-hand side of (5.6), one has

$[H_{0,0},\gamma_{1}^{(1)}]=2c_{1}c_{2}c_{4}-2c_{1}c_{3}c_{5}+2c_{2}c_{4}c_{5}+2c_{2}c_{3}c_{6}.$
(5.8)

We prepare the following:

$[H_{0,0},(-i)c_{3}c_{4}c_{5}]=\sum_{\ell=1}[c_{2\ell}c_{2\ell+1},c_{3}c_{4}c_{5}]=[c_{2}c_{3},c_{3}c_{4}c_{5}]+[c_{4}c_{5},c_{3}c_{4}c_{5}]=2c_{2}c_{4}c_{5},$
$[H_{0,0},(-i)c_{2}c_{3}c_{7}]=\sum_{\ell=1}[c_{2\ell}c_{2\ell+1},c_{2}c_{3}c_{7}]=[c_{2}c_{3},c_{2}c_{3}c_{7}]+[c_{6}c_{7},c_{2}c_{3}c_{7}]=2c_{2}c_{3}c_{6},$
$[H_{0,0},(-i)c_{1}c_{3}c_{4}]=\sum_{\ell=1}[c_{2\ell}c_{2\ell+1},c_{1}c_{3}c_{4}]=[c_{2}c_{3},c_{1}c_{3}c_{4}]+[c_{4}c_{5},c_{1}c_{3}c_{4}]=2c_{1}c_{2}c_{4}-2c_{1}c_{3}c_{5},$

and

$[H_{0,0},(-i)c_{1}c_{2}c_{5}]=\sum_{\ell=1}[c_{2\ell}c_{2\ell+1},c_{1}c_{2}c_{5}]=[c_{2}c_{3},c_{1}c_{2}c_{5}]+[c_{4}c_{5},c_{1}c_{2}c_{5}]=-2c_{1}c_{3}c_{5}+2c_{1}c_{2}c_{4}.$

Therefore, we have a solution,

$\gamma_{1}^{(1)}=-ic_{3}c_{4}c_{5}-ic_{2}c_{3}c_{7}-i\lambda_{1}^{(1)}c_{1}c_{3}c_{4}-i(1-\lambda_{1}^{(1)})c_{1}c_{2}c_{5},$

where $\lambda_{1}^{(1)}$ is a real parameter.

Further, $\gamma_{1}^{(2)}$ is determined by

$[H_{0,0},\gamma_{1}^{(2)}]=-[H_{0,1},\gamma_{1}^{(1)}]$ (5.9)

from the equation (5.7) for $j=1$. The right-hand side is calculated as

$\displaystyle-[H_{0,1},\gamma_{1}^{(1)}]$ $\displaystyle=$
$\displaystyle-\sum_{\ell=1}[c_{2\ell-1}c_{2\ell},c_{3}c_{4}c_{5}+c_{2}c_{3}c_{7}+\lambda_{1}^{(1)}c_{1}c_{3}c_{4}+(1-\lambda_{1}^{(1)})c_{1}c_{2}c_{5}]$
$\displaystyle=$ $\displaystyle
2c_{3}c_{4}c_{6}-2c_{1}c_{3}c_{7}+2c_{2}c_{4}c_{7}+2c_{2}c_{3}c_{8}+2\lambda_{1}^{(1)}c_{2}c_{3}c_{4}+2(1-\lambda_{1}^{(1)})c_{1}c_{2}c_{6}.$

Note that

$[H_{0,0},(-i)c_{2}c_{3}c_{9}]=\sum_{\ell=1}[c_{2\ell}c_{2\ell+1},c_{2}c_{3}c_{9}]=2c_{2}c_{3}c_{8},$
(5.11)

and

$[H_{0,0},(-i)c_{1}c_{2}c_{7}]=[c_{2}c_{3},c_{1}c_{2}c_{7}]+[c_{6}c_{7},c_{1}c_{2}c_{7}]=-2c_{1}c_{3}c_{7}+2c_{1}c_{2}c_{6}.$
(5.12)

By choosing $\lambda_{1}^{(1)}=0$, we can obtain the three terms in the right-
hand side of (5.2), and the rest are $2c_{3}c_{4}c_{6}$ and
$2c_{2}c_{4}c_{7}$.

In order to obtain the corresponding contributions to the rest of the two
terms, we prepare the following:

$I_{246}:=[H_{0,0},(-i)c_{2}c_{4}c_{6}]=-2c_{3}c_{4}c_{6}-2c_{2}c_{5}c_{6}-2c_{2}c_{4}c_{7},$
$I_{257}:=[H_{0,0},(-i)c_{2}c_{5}c_{7}]=-2c_{3}c_{5}c_{7}+2c_{2}c_{5}c_{6}+2c_{2}c_{4}c_{7},$
$I_{347}:=[H_{0,0},(-i)c_{3}c_{4}c_{7}]=2c_{3}c_{4}c_{6}-2c_{3}c_{5}c_{7}+2c_{2}c_{4}c_{7},$

and

$I_{356}:=[H_{0,0},(-i)c_{3}c_{5}c_{6}]=2c_{3}c_{4}c_{6}-2c_{3}c_{5}c_{7}+2c_{2}c_{5}c_{6}.$

From these results, one has

$2I_{246}+I_{257}-2I_{347}+I_{356}=-6c_{3}c_{4}c_{6}-6c_{2}c_{4}c_{7}.$

This implies

$[H_{0,0},(\frac{2}{3}ic_{2}c_{4}c_{6}+\frac{1}{3}ic_{2}c_{5}c_{7}-\frac{2}{3}ic_{3}c_{4}c_{7}+\frac{1}{3}ic_{3}c_{5}c_{6})]=2c_{3}c_{4}c_{6}+2c_{2}c_{4}c_{7}.$

Combining this, (5.9), (5.2), (5.11) and (5.12), we obtain

$\gamma_{1}^{(2)}=\frac{2}{3}ic_{2}c_{4}c_{6}+\frac{1}{3}ic_{2}c_{5}c_{7}-\frac{2}{3}ic_{3}c_{4}c_{7}+\frac{1}{3}ic_{3}c_{5}c_{6}-ic_{1}c_{2}c_{7}-ic_{2}c_{3}c_{9}.$

Thus, we can obtain the approximate solution of the edge mode by the series
expansion method.

##  Appendix A Stability of the spectral gap

In this appendix, we prove that the gap above the two-fold degenerate ground
state of the unperturbed Hamiltonian $H_{0}$ does not close under weak
interactions. As proved in Sec. 3, if there exists a non-vanishing spectral
gap above the sector of the ground state which are spanned by the ground-state
vector and the Majorana edge mode, then these two energies are always
degenerate. We assume that for a sufficiently small $|g|$, there exists a
spectral gap $\Delta\tilde{E}>0$ above the two-fold degenerate ground state
which continuously connects with the spectral gap in the case of $g=0$. Under
this assumption, we will estimate $|g|$ and $\Delta\tilde{E}$.

In order to separate the edge modes $\gamma_{0}$ and $c_{2N}$ from the bulk
variables so that the two edge variables are independent of the bulk
variables, we introduce new variables,
$\tilde{c}_{3},\tilde{c_{5}},\ldots,\tilde{c}_{2N-1}$:

$\mathcal{N}_{2\ell+1}\tilde{c}_{2\ell+1}=\sum_{m=1}^{\ell}\kappa^{m-1}c_{2m-1}-\frac{\sum_{m=1}^{\ell}\kappa^{2(m-1)}}{\kappa^{\ell}}c_{2\ell+1}$
(A.1)

for $\ell=1,2,\ldots,N-1$, where $\mathcal{N}_{2\ell+1}$ is the normalization
constant so that $\tilde{c}_{2\ell+1}$ satisfies
$(\tilde{c}_{2\ell+1})^{2}=1$. We recall the expression of the left edge mode,

$\gamma_{0}=\sum_{\ell=1}^{N}\kappa^{\ell-1}c_{2\ell-1}.$ (A.2)

From these expressions, one has

$\\{\gamma_{0},\tilde{c}_{2\ell+1}\\}=0$ (A.3)

for $\ell=1,2,\ldots,N-1$, and

$\\{\tilde{c}_{2\ell+1},\tilde{c}_{2\ell^{\prime}+1}\\}=0\quad\mbox{for \ \
}\ell\neq\ell^{\prime}.$ (A.4)

Since $\mathcal{N}_{2\ell+1}\sim\kappa^{-\ell}$ for a small $\kappa$, the
operator $\tilde{c}_{2\ell+1}$ is still localized around the site $2\ell+1$.

To rewrite the Hamiltonian $H_{0}$ in terms of the variables
$\tilde{c}_{2\ell+1}$, we note that

$H_{0}=i\kappa\sum_{\ell=1}^{N-1}c_{2\ell-1}c_{2\ell}+i\sum_{\ell=1}^{N-1}c_{2\ell}c_{2\ell+1}=i\sum_{\ell=1}^{N-1}(\kappa
c_{2\ell-1}-c_{2\ell+1})c_{2\ell}.$ (A.5)

From (A.1), one has

$\mathcal{N}_{2\ell+1}\tilde{c}_{2\ell+1}=\mathcal{N}_{2\ell-1}\tilde{c}_{2\ell-1}+\frac{\sum_{m=1}^{\ell}\kappa^{2(m-1)}}{\kappa^{\ell}}(\kappa
c_{2\ell-1}-c_{2\ell+1})$ (A.6)

for $\ell=2,3,\ldots,N-1$. This yields

$\kappa
c_{2\ell-1}-c_{2\ell+1}=\frac{\kappa^{\ell}}{\sum_{m=1}^{\ell}\kappa^{2(m-1)}}(\mathcal{N}_{2\ell+1}\tilde{c}_{2\ell+1}-\mathcal{N}_{2\ell-1}\tilde{c}_{2\ell-1})$
(A.7)

for $\ell=2,3,\ldots,N-1$. Substituting this and (A.1) for $\ell=1$ into the
above expression (A.5) of the Hamiltonian $H_{0}$, we have

$H_{0}=i\kappa\mathcal{N}_{3}\tilde{c}_{3}\tilde{c}_{2}+i\sum_{\ell=2}^{N-1}\frac{\kappa^{\ell}}{\sum_{m=1}^{\ell}\kappa^{2(m-1)}}(\mathcal{N}_{2\ell+1}\tilde{c}_{2\ell+1}-\mathcal{N}_{2\ell-1}\tilde{c}_{2\ell-1})\tilde{c}_{2\ell},$
(A.8)

where we have written $\tilde{c}_{2\ell}=c_{2\ell}$ for $\ell=1,2,\ldots,N-1$.
Thus, the Hamiltonian $H_{0}$ can be written in terms of the variables,
$\tilde{c}_{2},\tilde{c}_{3},\ldots,\tilde{c}_{2N-1}$. Although the new set of
the variables is given by the $2N$ variables,
$\\{\hat{\gamma}_{0},\tilde{c}_{2},\tilde{c}_{3},\ldots,\tilde{c}_{2N-1},c_{2N}\\}$,
the expression of the Hamiltonian $H_{0}$ is independent of the Majorana edge
zero modes, $\hat{\gamma}_{0}$ and $c_{2N}$. Namely, the two sets of the
variables are separated in the expression of $H_{0}$.

By using a pure imaginary, antisymmetric matrix $\tilde{A}$, the Hamiltonian
$H_{0}$ can be written in the form,

$H_{0}=\sum_{i,j\in\\{2,3,\ldots,2N-1\\}}\tilde{c}_{i}\tilde{A}_{i,j}\tilde{c}_{j}.$
(A.9)

By the construction, the Hamiltonian $H_{0}$ has the unique ground state with
the non-vanishing spectral gap above it, except for the Majorana edge zero
modes. In the following, we will treat these two ground states separately.

Following Hastings [20], we introduce two copies of the Majorana fermions, and
two copies of the Hamiltonian ${H}_{0}$ with opposite signs so that the total
Hamiltonian is given by

$\tilde{H}_{0}:=\sum_{i,j\in\\{2,3,\ldots,2N-1\\}}\tilde{c}_{i}\tilde{A}_{i,j}\tilde{c}_{j}-\sum_{i,j\in\\{2,3,\ldots,2N-1\\}}\tilde{d}_{i}\tilde{A}_{i,j}\tilde{d}_{j},$
(A.10)

where the Majorana operator $\tilde{d}_{j}$ is a copy of $\tilde{c}_{j}$. This
Hamiltonian has also the unique ground state with the non-vanishing spectral
gap above it. We write

$\tilde{H}_{0}=\begin{pmatrix}\tilde{c},&\tilde{d}\end{pmatrix}\begin{pmatrix}\tilde{A}&0\\\
0&-\tilde{A}\end{pmatrix}\begin{pmatrix}\tilde{c}\\\ \tilde{d}\end{pmatrix},$

and $P_{\pm}(\tilde{A})$ for the projection onto the positive and negative
energies for $\tilde{A}$, respectively. We also introduce [20]

$\hat{U}:=\frac{1}{\sqrt{2}}\begin{pmatrix}1&is(\tilde{A})\\\
is(\tilde{A})&1\end{pmatrix},$

where $s(\tilde{A}):=P_{+}(\tilde{A})-P_{-}(\tilde{A})$. Then, one has

$\hat{U}^{\dagger}\begin{pmatrix}\tilde{A}&0\\\
0&-\tilde{A}\end{pmatrix}\hat{U}=\begin{pmatrix}0&i|\tilde{A}|\\\
-i|\tilde{A}|&0\end{pmatrix}.$

We define a pair of Majorana fermions,

$\begin{pmatrix}\eta^{\rm R}\\\ \eta^{\rm
I}\end{pmatrix}:=\hat{U}^{\dagger}\begin{pmatrix}\tilde{c}\\\
\tilde{d}\end{pmatrix}.$ (A.11)

From the above observation, one has

$\tilde{H}_{0}=\begin{pmatrix}\eta^{\rm R},&\eta^{\rm
I}\end{pmatrix}\begin{pmatrix}0&i|\tilde{A}|\\\
-i|\tilde{A}|&0\end{pmatrix}\begin{pmatrix}\eta^{\rm R}\\\ \eta^{\rm
I}\end{pmatrix}=2i\eta^{\rm R}|\tilde{A}|\eta^{\rm I}.$

Further, we define a complex fermion,

$\eta_{i}:=\frac{1}{2}(\eta_{i}^{\rm R}+i\eta_{i}^{\rm I}),$ (A.12)

for $i=2,3,\ldots,2N-1$. Then,

$\tilde{H}_{0}=2\eta^{\dagger}|\tilde{A}|\eta,$ (A.13)

where we have dropped the constant in the right-hand side. Since there is a
spectral gap at the zero energy in the spectrum of the matrix $\tilde{A}$, we
have

$\tilde{H}_{0}=2\eta^{\dagger}|\tilde{A}|\eta\geq\Delta
E_{0}\eta^{\dagger}\eta$ (A.14)

with a constant $\Delta E_{0}>0$.

Next, consider the interaction term $V$ of the Hamiltonian $H_{g}$ of (3.1).
The interaction $V$ is written in terms of the operators,
$c_{1},c_{2},\ldots,c_{2N-1}$. Our aim is to express the interaction $V$ in
terms of the operators, $\eta$ and $a_{\rm E}$.

For this purpose, we first want to express the operators,
$c_{1},c_{3},\ldots,c_{2N-1}$, in terms of
$\gamma_{0},\tilde{c}_{3},\tilde{c}_{5},\ldots,\tilde{c}_{2N-1}$. From the
expression (A.2) of $\gamma_{0}$, one has

$c_{1}=\gamma_{0}-(\kappa c_{3}+\kappa^{2}c_{5}+\kappa^{3}c_{7}+\cdots).$
(A.15)

Combining this with (A.1) for $\ell=1$, one obtains

$\displaystyle c_{3}$ $\displaystyle=$ $\displaystyle\kappa
c_{1}-\kappa\mathcal{N}_{3}\tilde{c}_{3}$ (A.16) $\displaystyle=$
$\displaystyle\kappa\gamma_{0}-\kappa(\kappa
c_{3}+\kappa^{2}c_{5}+\kappa^{3}c_{7}+\cdots)-\kappa\mathcal{N}_{3}\tilde{c}_{3}$
$\displaystyle=$
$\displaystyle\kappa\gamma_{0}-\kappa^{2}c_{3}-\kappa(\kappa^{2}c_{5}+\kappa^{3}c_{7}+\cdots)-\kappa\mathcal{N}_{3}\tilde{c}_{3}.$

Therefore,

$c_{3}=\frac{\kappa}{1+\kappa^{2}}[\gamma_{0}-(\kappa^{2}c_{5}+\kappa^{3}c_{7}+\cdots)]-\frac{\kappa}{1+\kappa^{2}}\mathcal{N}_{3}\tilde{c}_{3}.$
(A.17)

By using this expression, we can eliminate $c_{3}$ from the right-hand side of
(A.15). Similarly, from (A.1) for $\ell=2$, one has

$\displaystyle\mathcal{N}_{5}\tilde{c}_{5}$ $\displaystyle=$ $\displaystyle
c_{1}+\kappa c_{3}-\frac{1+\kappa^{2}}{\kappa^{2}}c_{5}$ $\displaystyle=$
$\displaystyle\gamma_{0}-\kappa^{2}c_{5}-(\kappa^{3}c_{7}+\kappa^{4}c_{9}+\cdots)-\frac{1+\kappa^{2}}{\kappa^{2}}c_{5}.$

Therefore,

$c_{5}=\frac{\kappa^{2}}{1+\kappa^{2}+\kappa^{4}}[\gamma_{0}-(\kappa^{3}c_{7}+\kappa^{4}c_{9}+\cdots)]-\frac{\kappa^{2}}{1+\kappa^{2}+\kappa^{4}}\mathcal{N}_{5}\tilde{c}_{5}.$
(A.18)

By using this expression, we can eliminate $c_{5}$ from both of the right-hand
sides of (A.15) and (A.17). Therefore, by repeating this procedure, the
operators, $c_{1},c_{3},\ldots,c_{2N-1}$, can be expressed in terms of
$\gamma_{0},\tilde{c}_{3},\tilde{c}_{5},\ldots,\tilde{c}_{2N-1}$. Since
$c_{2\ell}=\tilde{c}_{2\ell}$, the interaction potential $V$ can be expressed
in terms of the operators,
$\gamma_{0},\tilde{c}_{2},\tilde{c}_{3},\tilde{c}_{4},\ldots,\tilde{c}_{2N-2},\tilde{c}_{2N-1}$.
Then, the resulting terms in the interaction decay by exponential law with the
distance about these operators because the interaction $V$ which is written in
terms of $c_{1},c_{2},\ldots,c_{2N-1}$ is assumed to be finite range. Further,
from (A.11) and (A.12), one has

$\tilde{c}_{i}=\frac{1}{\sqrt{2}}\Bigl{[}\eta_{i}+\eta_{i}^{\dagger}-\sum_{j=2}^{2N-1}s(\tilde{A})_{i,j}(\eta_{j}-\eta_{j}^{\dagger})\Bigr{]}$
(A.19)

for $i=2,3,\ldots,2N-1$. We recall $\hat{\gamma}_{0}=a_{\rm E}+a_{\rm
E}^{\dagger}$. From these observations, one notices that the interaction $V$
can be written in terms of $a_{\rm E},\eta_{2},\eta_{3},\ldots,\eta_{2N-1}$.
In the following, we denote the interaction by

$\tilde{V}:={V}(a_{\rm E},\eta_{2},\eta_{3},\ldots,\eta_{2N-1}).$ (A.20)

The total Hamiltonian which we consider is given by

$\tilde{H}_{g}:=\tilde{H}_{0}+g\tilde{V},$ (A.21)

where $\tilde{H}_{0}$ is given by (A.13). From the definitions (A.9) and
(A.10), this right-hand side is written

$\tilde{H}_{0}+g\tilde{V}=H_{0}+gV-\tilde{d}\tilde{A}\tilde{d}.$

Clearly, the first and second terms in the right-hand side are those of the
original Hamiltonian $H_{g}$, and the third term is the additional Hamiltonian
whose spectrum has the gap above the unique ground state. Thus, it is enough
to show that the spectrum of the Hamiltonian $\tilde{H}_{g}$ has a non-
vanishing spectral gap above the two low-energy states.

As mentioned at the beginning of this section, we assume that for a small
$|g|$, there exists a spectral gap $\Delta\tilde{E}>0$ above the two-fold
degenerate ground state which continuously connects with the spectral gap in
the case of $g=0$. Both of the estimates of $|g|$ and $\Delta\tilde{E}$ will
be obtained below.

We write $|\tilde{0}\rangle$ for the ground state of $\tilde{H}_{0}$ of (A.13)
which satisfies $\eta_{i}|\tilde{0}\rangle=0$ for all $i=2,3,\ldots,2N-1$ and
$a_{\rm E}|\tilde{0}\rangle=0$ for the operator $a_{\rm E}$ of the Majorana
edge zero mode. Then, the other ground state is given by $a_{\rm
E}^{\dagger}|\tilde{0}\rangle$. We also write $\tilde{P}_{0}(0)$ for the
projection onto the sector of these two ground states. Under the above
assumption on the spectral gap $\Delta\tilde{E}>0$ above the two low-energy
states of $\tilde{H}_{g}$, there exists a unitary operator $\tilde{U}(g)$ such
that [23]

$\tilde{U}^{\dagger}(g)\tilde{P}_{0}(g)\tilde{U}(g)=\tilde{P}_{0}(0),$ (A.22)

where $\tilde{P}_{0}(g)$ is the projection onto the sector of the two low-
energy states of $\tilde{H}_{g}$.

Note that

$\tilde{H}_{g}=[1-\tilde{P}_{0}(g)]\tilde{H}_{g}[1-\tilde{P}_{0}(g)]+\tilde{P}_{0}(g)\tilde{H}_{g}\tilde{P}_{0}(g).$

Therefore, one has

$\displaystyle\tilde{H}(g)$ $\displaystyle:=$
$\displaystyle\tilde{U}^{\dagger}(g)\tilde{H}_{g}\tilde{U}(g)$
$\displaystyle=$
$\displaystyle[1-\tilde{P}_{0}(0)]\tilde{U}^{\dagger}(g)\tilde{H}_{g}\tilde{U}(g)[1-\tilde{P}_{0}(0)]+\tilde{P}_{0}(0)\tilde{U}^{\dagger}(g)\tilde{H}_{g}\tilde{U}(g)\tilde{P}_{0}(0).$

This implies that $\tilde{P}_{0}(0)$ is the projection onto the two low-energy
states for the transformed Hamiltonian $\tilde{H}(g)$. In addition, the vector
$|\tilde{0}\rangle$ is the eigenvector of $\tilde{H}(g)$. In fact, one has

$\langle\tilde{0}|a_{\rm
E}\tilde{U}^{\dagger}(g)\tilde{H}_{g}\tilde{U}(g)|\tilde{0}\rangle=0$

because $\tilde{U}^{\dagger}(g)\tilde{H}_{g}\tilde{U}(g)$ has even fermion
parity by the construction [24] of the unitary operator $\tilde{U}(g)$. Thus,
it is sufficient to show that there exists a non-vanishing spectral gap above
the two low-energy states for this Hamiltonian $\tilde{H}(g)$.

We write

$\tilde{H}_{0}=\sum_{Z\subset\Lambda}\tilde{H}_{0,Z},$

where $\Lambda=\\{1,2,\ldots,2N-1\\}$, and $\tilde{H}_{0,Z}$ is the local
Hamiltonian whose explicit form is given by

$\tilde{H}_{0,Z}=2\eta_{i}^{\dagger}|\tilde{A}|_{i,j}\eta_{j}$

with $Z=\\{i,j\\}$ for $i,j=2,3,\ldots,2N-1$. Clearly, one has

$\tilde{H}_{0,Z}|\tilde{0}\rangle=0$ (A.23)

for all $Z$. Similarly, we can write the interaction $\tilde{V}$ in the sum of
the local interactions $\tilde{V}_{Z}$,

$\tilde{V}=\sum_{Z\subset\Lambda}\tilde{V}_{Z},$

where the support of $a_{\rm E}$ is given by the single site $\\{1\\}$. We
write

$\tilde{H}_{0}(g):=\tilde{U}^{\dagger}(g)\tilde{H}_{0}\tilde{U}(g)=\sum_{Z\subset\Lambda}\tilde{U}^{\dagger}(g)\tilde{H}_{0,Z}\tilde{U}(g)$
(A.24)

and

$\tilde{V}(g):=\tilde{U}^{\dagger}(g)\tilde{V}\tilde{U}(g)=\sum_{Z\subset\Lambda}\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g).$
(A.25)

Obviously,

$\tilde{H}(g)=\tilde{H}_{0}(g)+g\tilde{V}(g).$ (A.26)

Following Hastings [20], we introduce

$\mathcal{R}_{g}(\cdots):=\int_{-\infty}^{+\infty}\tilde{\tau}_{g,t}(\cdots)w(t)dt,$

where

$\tilde{\tau}_{g,t}(\cdots):=\exp[it\tilde{H}(g)](\cdots)\exp[-it\tilde{H}(g)],$

and $w(t)$ is a function such that its Fourier transform $\hat{w}$ satisfies
$\hat{w}(0)=1$ and $\hat{w}(\omega)=0$ for $|\omega|\geq\Delta\tilde{E}$. From
the condition $\hat{w}(0)=1$, one has $\mathcal{R}_{g}(1)=1$ and

$\mathcal{R}_{g}(\tilde{H}(g))=\tilde{H}(g).$ (A.27)

Since the Hamiltonian $\tilde{H}(g)$ has the spectral gap $\Delta\tilde{E}>0$
above the sector of the two low-energy states from the assumption on the
Hamiltonian $\tilde{H}_{g}$, the condition $\hat{w}(\omega)=0$ for
$|\omega|\geq\Delta\tilde{E}$ yields

$[1-\tilde{P}_{0}(0)]\mathcal{R}_{g}(\mathcal{A})\tilde{P}_{0}(0)=0$ (A.28)

for any operator $\mathcal{A}$. As a concrete function, we use the function
$w(t)$ of the equation (2.1) in Lemma 2.3 in [24]. From (A.24), (A.25), (A.26)
and (A.27), we have

$\tilde{H}(g)=\sum_{Z\subset\Lambda}\mathcal{R}_{g}(\tilde{U}^{\dagger}(g)\tilde{H}_{0,Z}\tilde{U}(g))+g\sum_{Z\subset\Lambda}\mathcal{R}_{g}(\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g)).$
(A.29)

Consider first the summand in the second sum in the right-hand side of (A.29).
Note that

$\displaystyle\mathcal{R}_{g}(\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g))|\tilde{0}\rangle$
$\displaystyle=$
$\displaystyle\tilde{P}_{0}(0)\mathcal{R}_{g}(\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g))|\tilde{0}\rangle$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}dt\;w(t)\tilde{P}_{0}(0)e^{it\tilde{H}(g)}\tilde{P}_{0}(0)\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g)|\tilde{0}\rangle
e^{-it\tilde{E}_{0}(g)}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}dt\;w(t)|\tilde{0}\rangle\langle\tilde{0}|\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g)|\tilde{0}\rangle$
$\displaystyle=$
$\displaystyle|\tilde{0}\rangle\langle\tilde{0}|\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g)|\tilde{0}\rangle,$

where $\tilde{E}_{0}(g)$ is the energy eigenvalue of $\tilde{H}(g)$ for the
eigenvector $|\tilde{0}\rangle$, and we have used (A.28) and $\hat{w}(0)=1$;
we have also used $\langle\tilde{0}|a_{\rm
E}\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g)|\tilde{0}\rangle$=0, which
holds because the operator $\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g)$
has even fermion parity. This can be written

$\mathcal{R}_{g}(:\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g):)|\tilde{0}\rangle=0,$
(A.30)

where $:(\cdots):$ is the subtraction of the expectation value with respect to
the state $|\tilde{0}\rangle$ which is defined by

$:\mathcal{A}:=\mathcal{A}-\langle\tilde{0}|\mathcal{A}|\tilde{0}\rangle$
(A.31)

for an operator $\mathcal{A}$. We write

$\tilde{\mathcal{W}}_{Z}^{(1)}(g):=g\mathcal{R}_{g}(:\tilde{U}^{\dagger}(g)\tilde{V}_{Z}\tilde{U}(g):).$
(A.32)

Next consider the summand in the first sum in the right-hand side of (A.29).
Note that

$\displaystyle\tilde{U}^{\dagger}(g)\tilde{H}_{0,Z}\tilde{U}(g)-\tilde{H}_{0,Z}$
$\displaystyle=$
$\displaystyle\int_{0}^{g}dg^{\prime}\frac{d}{dg^{\prime}}\tilde{U}^{\dagger}(g^{\prime})\tilde{H}_{0,Z}\tilde{U}(g^{\prime})$
(A.33) $\displaystyle=$
$\displaystyle\int_{0}^{g}dg^{\prime}\tilde{U}^{\dagger}(g^{\prime})i[\tilde{H}_{0,Z},\tilde{D}(g^{\prime})]\tilde{U}(g^{\prime}),$

where we have used [24]

$\frac{d}{dg}\tilde{U}(g)=i\tilde{D}(g)\tilde{U}(g)$

with the self-adjoint operator $\tilde{D}(g)$. Therefore, we have

$\mathcal{R}_{g}(\tilde{U}^{\dagger}(g)\tilde{H}_{0,Z}\tilde{U}(g))=\mathcal{R}_{g}(\tilde{H}_{0,Z})+\int_{0}^{g}dg^{\prime}\mathcal{R}_{g}(\tilde{U}^{\dagger}(g^{\prime})i[\tilde{H}_{0,Z},\tilde{D}(g^{\prime})]\tilde{U}(g^{\prime})).$
(A.34)

The second term in the right-hand side satisfies

$\int_{0}^{g}dg^{\prime}\mathcal{R}_{g}(:\tilde{U}^{\dagger}(g^{\prime})i[\tilde{H}_{0,Z},\tilde{D}(g^{\prime})]\tilde{U}(g^{\prime}):)|\tilde{0}\rangle=0$

in the same way as in the above. Therefore, we write

$\tilde{\mathcal{W}}_{Z}^{(2)}(g):=\int_{0}^{g}dg^{\prime}\mathcal{R}_{g}(:\tilde{U}^{\dagger}(g^{\prime})i[\tilde{H}_{0,Z},\tilde{D}(g^{\prime})]\tilde{U}(g^{\prime}):).$
(A.35)

In order to deal with the first term $\mathcal{R}_{g}(\tilde{H}_{0,Z})$ in the
right-hand side of (A.34), we note that

$\displaystyle
e^{it\tilde{H}(g)}\tilde{H}_{0,Z}e^{-it\tilde{H}(g)}-\tilde{H}_{0,Z}$
$\displaystyle=$
$\displaystyle\int_{0}^{t}dt^{\prime}\frac{d}{dt^{\prime}}e^{it^{\prime}\tilde{H}(g)}\tilde{H}_{0,Z}e^{-it^{\prime}\tilde{H}(g)}$
(A.36) $\displaystyle=$
$\displaystyle\int_{0}^{t}dt^{\prime}\;e^{it^{\prime}\tilde{H}(g)}i[\tilde{H}(g),\tilde{H}_{0,Z}]e^{-it^{\prime}\tilde{H}(g)}.$

On the other hand, the above relation (A.33) yields

$\displaystyle\tilde{H}(g)$ $\displaystyle=$
$\displaystyle\tilde{U}^{\dagger}(g)\tilde{H}_{0}\tilde{U}(g)+g\tilde{U}^{\dagger}(g)\tilde{V}\tilde{U}(g)$
$\displaystyle=$
$\displaystyle\tilde{H}_{0}+\int_{0}^{g}dg^{\prime}\;\tilde{U}^{\dagger}(g^{\prime})i[\tilde{H}_{0},\tilde{D}(g^{\prime})]\tilde{U}(g^{\prime})+g\tilde{U}^{\dagger}(g)\tilde{V}\tilde{U}(g).$

Substituting this into the right-hand side of (A.36), we have

$\displaystyle e^{it\tilde{H}(g)}\tilde{H}_{0,Z}e^{-it\tilde{H}(g)}$
$\displaystyle=$
$\displaystyle\tilde{H}_{0,Z}+\int_{0}^{t}dt^{\prime}\;e^{it^{\prime}\tilde{H}(g)}i[\tilde{H}_{0},\tilde{H}_{0,Z}]e^{-it^{\prime}\tilde{H}(g)}+\tilde{\mathcal{M}}_{t,Z}(g),$
(A.37)

where

$\tilde{\mathcal{M}}_{t,Z}(g):=\int_{0}^{t}dt^{\prime}\;\tilde{\tau}_{g,t^{\prime}}\left(ig[\tilde{U}^{\dagger}(g)\tilde{V}\tilde{U}(g),\tilde{H}_{0,Z}]+\int_{0}^{g}dg^{\prime}\;i[\tilde{U}^{\dagger}(g^{\prime})i[\tilde{H}_{0},\tilde{D}(g^{\prime})]\tilde{U}(g^{\prime}),\tilde{H}_{0,Z}]\right).$
(A.38)

From the expression of (A.37), one notices that

$\tilde{\mathcal{M}}_{t,Z}(g)|\tilde{0}\rangle=0.$ (A.39)

Actually, $|\tilde{0}\rangle$ is the eigenvector of $\tilde{H}(g)$, and
$\tilde{H}_{0,Z}|\tilde{0}\rangle=0$. In addition, we obtain

$\sum_{Z\subset\Lambda}\mathcal{R}_{g}(\tilde{H}_{0,Z})=\tilde{H}_{0}+\sum_{Z\subset\Lambda}\int_{-\infty}^{+\infty}dt\;w(t)\tilde{\mathcal{M}}_{t,Z}(g).$
(A.40)

We write

$\tilde{\mathcal{W}}_{Z}^{(3)}(g):=\int_{-\infty}^{+\infty}dt\;w(t)\tilde{\mathcal{M}}_{t,Z}(g).$
(A.41)

Consequently, from (A.29), (A.32), (A.34), (A.35), (A.40) and (A.41), we
obtain the desired expression of the Hamiltonian $\tilde{H}(g)$,

$\tilde{H}(g)=\tilde{H}_{0}+\tilde{\mathcal{W}}(g),$ (A.42)

with

$\tilde{\mathcal{W}}(g):=\sum_{i=1}^{3}\sum_{Z\subset\Lambda}\tilde{\mathcal{W}}_{Z}^{(i)}(g),$
(A.43)

where we have dropped the constant term. The operators
$\tilde{\mathcal{W}}_{Z}^{(i)}(g)$ satisfy

$\tilde{\mathcal{W}}_{Z}^{(i)}(g)\rightarrow 0\quad\mbox{as \ \ }g\rightarrow
0,$ (A.44)

and

$\tilde{\mathcal{W}}_{Z}^{(i)}(g)|\tilde{0}\rangle=0$ (A.45)

for all $Z$ and $i=1,2,3$. Therefore, in the same way as in [22], we obtain

$[\tilde{\mathcal{W}}(g)]^{2}\leq(\tilde{g})^{2}(\tilde{\mathfrak{N}})^{2},$
(A.46)

where $\tilde{g}$ is a positive constant which satisfies $\tilde{g}\rightarrow
0$ as $g\rightarrow 0$, and

$\tilde{\mathfrak{N}}:=a_{\rm E}^{\dagger}a_{\rm
E}+\sum_{i=2}^{2N-1}\eta_{i}^{\dagger}\eta_{i}.$ (A.47)

By using the bound (A.14), we have

$\displaystyle[\tilde{\mathcal{W}}(g)]^{2}$ $\displaystyle\leq$
$\displaystyle(\tilde{g})^{2}(\tilde{\mathfrak{N}})^{2}$ (A.48)
$\displaystyle=$ $\displaystyle(\tilde{g})^{2}\Bigl{(}a_{\rm
E}^{\dagger}a_{\rm E}+\sum_{i}\eta_{i}^{\dagger}\eta_{i}\Bigr{)}^{2}$
$\displaystyle\leq$ $\displaystyle(\tilde{g})^{2}\Bigl{(}1+\frac{1}{\Delta
E_{0}}\tilde{H}_{0}\Bigr{)}^{2}$ $\displaystyle\leq$ $\displaystyle
2(\tilde{g})^{2}\Bigl{[}1+\frac{1}{(\Delta
E_{0})^{2}}(\tilde{H}_{0})^{2}\Bigr{]}.$

By the construction of the Hamiltonian
$\tilde{H}(g)=\tilde{H}_{0}+\tilde{\mathcal{W}}(g)$, one has

$\langle\tilde{0}|\tilde{H}(g)|\tilde{0}\rangle=\langle\tilde{0}|\tilde{H}_{0}|\tilde{0}\rangle+\langle\tilde{0}|\tilde{\mathcal{W}}(g)|\tilde{0}\rangle=0.$

This implies that the ground state energy $\tilde{E}_{0}(g)$ of $\tilde{H}(g)$
satisfies $\tilde{E}_{0}(g)\leq 0$.

Let $\tilde{\Psi}_{1}$ be an excited state of $\tilde{H}(g)$ such that the
energy eigenvalue $\tilde{E}_{1}(g)$ does not connect continuously with the
ground-state energy zero of the unperturbed Hamiltonian $\tilde{H}_{0}$ when
varying the coupling constant $g$ from non-zero to zero. From
$\tilde{H}(g)\tilde{\Psi}_{1}=\tilde{E}_{1}(g)\tilde{\Psi}_{1}$ with the
normalization $\|\tilde{\Psi}_{1}\|=1$, one has

$\displaystyle[\tilde{E}_{1}(g)]^{2}$ $\displaystyle=$
$\displaystyle\langle\tilde{\Psi}_{1},[\tilde{H}(g)]^{2}\tilde{\Psi}_{1}\rangle$
(A.49) $\displaystyle=$
$\displaystyle\langle\tilde{\Psi}_{1},(\tilde{H}_{0})^{2}\tilde{\Psi}_{1}\rangle+\langle\tilde{\Psi}_{1},[\tilde{\mathcal{W}}(g)]^{2}\tilde{\Psi}_{1}\rangle+\langle\tilde{\Psi}_{1},\tilde{H}_{0}\tilde{\mathcal{W}}(g)\tilde{\Psi}_{1}\rangle+\langle\tilde{\Psi}_{1},\tilde{\mathcal{W}}(g)\tilde{H}_{0}\tilde{\Psi}_{1}\rangle$
$\displaystyle\geq$
$\displaystyle\langle\tilde{\Psi}_{1},(\tilde{H}_{0})^{2}\tilde{\Psi}_{1}\rangle+\langle\tilde{\Psi}_{1},\tilde{H}_{0}\tilde{\mathcal{W}}(g)\tilde{\Psi}_{1}\rangle+\langle\tilde{\Psi}_{1},\tilde{\mathcal{W}}(g)\tilde{H}_{0}\tilde{\Psi}_{1}\rangle.$

The second and third terms in the right-hand side can be estimated as follows:

$\displaystyle\bigl{|}{\rm
Re}\langle\tilde{\Psi}_{1},\tilde{H}_{0}\tilde{\mathcal{W}}(g)\tilde{\Psi}_{1}\rangle\bigr{|}$
$\displaystyle\leq$
$\displaystyle\sqrt{\langle\tilde{\Psi}_{1},(\tilde{H}_{0})^{2}\tilde{\Psi}_{1}\rangle\langle\tilde{\Psi}_{1},[\tilde{\mathcal{W}}(g)]^{2}\tilde{\Psi}_{1}\rangle}$
(A.50) $\displaystyle\leq$ $\displaystyle\frac{\sqrt{2}\tilde{g}}{\Delta
E_{0}}\sqrt{\langle\tilde{\Psi}_{1},(\tilde{H}_{0})^{2}\tilde{\Psi}_{1}\rangle\langle\tilde{\Psi}_{1},[(\Delta
E_{0})^{2}+(\tilde{H}_{0})^{2}]\tilde{\Psi}_{1}\rangle}$ $\displaystyle\leq$
$\displaystyle\sqrt{2}\tilde{g}\Delta E_{0}+\frac{\sqrt{2}\tilde{g}}{\Delta
E_{0}}\langle\tilde{\Psi}_{1},(\tilde{H}_{0})^{2}\tilde{\Psi}_{1}\rangle,$

where we have used the bound (A.48). Substituting this into the above right-
hand side, we have

$[\tilde{E}_{1}(g)]^{2}+2\sqrt{2}\tilde{g}\Delta
E_{0}\geq\left[1-\frac{2\sqrt{2}\tilde{g}}{\Delta
E_{0}}\right]\langle\tilde{\Psi}_{1},(\tilde{H}_{0})^{2}\tilde{\Psi}_{1}\rangle.$

Further, using the bound (A.48) again, we obtain

$\displaystyle{2(\tilde{g})^{2}}\frac{[\tilde{E}_{1}(g)]^{2}+2\sqrt{2}\tilde{g}\Delta
E_{0}}{{(\Delta E_{0})^{2}}-{2\sqrt{2}\tilde{g}}{\Delta
E_{0}}}+2(\tilde{g})^{2}$ $\displaystyle\geq$
$\displaystyle\langle\tilde{\Psi}_{1},[\tilde{\mathcal{W}}(g)]^{2}\tilde{\Psi}_{1}\rangle$
(A.51) $\displaystyle=$
$\displaystyle\langle\tilde{\Psi}_{1},[\tilde{H}(g)-\tilde{H}_{0}]^{2}\tilde{\Psi}_{1}\rangle$
$\displaystyle=$
$\displaystyle\langle\tilde{\Psi}_{1},[\tilde{E}_{1}(g)-\tilde{H}_{0}]^{2}\tilde{\Psi}_{1}\rangle.$

Since we can assume $\tilde{E}_{1}(g)\leq\Delta E_{0}$, we have the bound,

$\langle\tilde{\Psi}_{1},[\tilde{E}_{1}(g)-\tilde{H}_{0}]^{2}\tilde{\Psi}_{1}\rangle\leq\frac{4\Delta
E_{0}}{\Delta E_{0}-2\sqrt{2}\tilde{g}}(\tilde{g})^{2}.$ (A.52)

When $\tilde{E}_{1}(g)$ satisfies $\Delta E_{0}\geq\tilde{E}_{1}(g)>\Delta
E_{0}/2$ for a small $|g|$, one has

$\frac{1}{2}\Delta
E_{0}>\sqrt{\langle\tilde{\Psi}_{1},[\tilde{E}_{1}(g)-\tilde{H}_{0}]^{2}\tilde{\Psi}_{1}\rangle}$

from the spectrum of $\tilde{H}_{0}$. Therefore, from the above bound (A.52),
we obtain

$\frac{1}{2}\Delta E_{0}<\Delta E_{0}-\left(\frac{4\Delta E_{0}}{\Delta
E_{0}-2\sqrt{2}\tilde{g}}\right)^{1/2}\tilde{g}\leq\tilde{E}_{1}(g)\leq\Delta
E_{0}$ (A.53)

for a small $|g|$. This implies that there exists a small positive constant
$g_{\rm max}$ such that the bound (A.53) can be realized for any coupling
constant $g$ which satisfies $|g|\leq g_{\rm max}$.

##  Appendix B Proof of the existence of $\gamma$ satisfying $[H,\gamma]=0$

In this appendix, we give a proof of the existence of a Majorana operator
$\gamma$ which commutes with the Hamiltonian $H$, i.e., $[H,\gamma]=0$,
following the paper [25] by Goldstein and Chamon. The Hamiltonian $H$ has even
fermion parity and consists of the operators, $c_{1},c_{2},\ldots,c_{2N-1}$,
whose number is odd $2N-1$.

Let us consider all the monomials which have odd fermion parity, and consist
of the operators, $c_{1},c_{2},\ldots,c_{2N-1}$. There are in total
$\sum_{m=1}^{N}{2N-1\choose 2m-1}=2^{2N-2}$ such operators. We denote them by
$C_{a}$ with indices $a$. We can write the operators $C_{a}$ in the form,

$C_{a}:=(i)^{n_{a}(n_{a}-1)/2}c_{i_{1}(a)}c_{i_{2}(a)}\cdots
c_{i_{n_{a}}(a)},$ (B.1)

where $n_{a}$ is the number of $c$-Majorana fermions in the product, and for
each monomial with index $a$,
$I(a):=\\{i_{1}(a),i_{2}(a),\ldots,i_{n_{a}}(a)\\}$ is the list of the sites
which satisfy conditions, $i_{j}(a)\in\\{1,2,\cdots,2N-1\\}$ for
$j\in\\{1,2,\ldots,n_{a}\\}$ and $i_{k}(a)<i_{\ell}(a)$ for $k<\ell$. Of
course, $n_{a}$ is odd and satisfies $n_{a}\leq 2N-1$. Then, one has

$C_{a}^{\dagger}=C_{a}\quad\mbox{and}\quad(C_{a})^{2}=1.$

Clearly, the operator $\gamma$ is written as

$\gamma=\sum_{a}\alpha_{a}C_{a}$ (B.2)

with coefficients $\alpha_{a}$.

The product of two $C_{a_{1}}$ and $C_{a_{2}}$ gives

$C_{a_{1}}C_{a_{2}}=(i)^{s(a_{1},a_{2})}C_{a_{3}}$

with some phase factor, where the set $I(a_{3})$ of the sites in $C_{a_{3}}$
satisfies $I(a_{3})=I(a_{1})\cup I(a_{2})\backslash I(a_{1})\cap I(a_{2})$. We
introduce the usual trace inner product for two operators, $A$ and $B$:

$(A,B):=\frac{1}{2^{2N}}{\rm Tr}A^{\dagger}B.$

Then, one has

$(C_{a},C_{b})=\delta_{a,b}$ (B.3)

for all $C_{a}$ and $C_{b}$. Actually,

${\rm Tr}\;c_{j_{1}}c_{j_{2}}\cdots c_{j_{n}}=0$

for any monomial $c_{j_{1}}c_{j_{2}}\cdots c_{j_{n}}$ with $n\geq 1$. In the
case of $n={\rm odd}$, this is true because of the odd fermion parity. When
$n={\rm even}$, one has

${\rm Tr}\;c_{j_{1}}c_{j_{2}}\cdots c_{j_{n}}={\rm Tr}\;c_{j_{2}}\cdots
c_{j_{n}}c_{j_{1}}=(-1)^{n-1}{\rm Tr}\;c_{j_{1}}c_{j_{2}}\cdots
c_{j_{n}}=-{\rm Tr}\;c_{j_{1}}c_{j_{2}}\cdots c_{j_{n}}$

by using the cyclic property of trace and the commutation relations. This
implies the vanishing of the trace.

Let us consider

$H_{a,b}:=(C_{a},[H,C_{b}])$ (B.4)

for the Hamiltonian $H$ with even fermion parity. This is a $2^{2N-2}\times
2^{2N-2}$ matrix. Further, this is pure imaginary and antisymmetric. Actually,
one has

$\overline{H_{a,b}}=\overline{(C_{a},[H,C_{b}])}=([H,C_{b}],C_{a})=2^{-2N}{\rm
Tr}(-[H,C_{b}])C_{a}=-2^{-2N}{\rm Tr}\;C_{a}[H,C_{b}]=-H_{a,b}.$

Thus, $H_{a,b}$ is pure imaginary. Further, the cyclic property of trace
yields

${\rm Tr}\;C_{a}[H,C_{b}]={\rm Tr}\;(C_{a}HC_{b}-C_{a}C_{b}H)={\rm
Tr}\;C_{b}C_{a}H-{\rm Tr}\;C_{b}HC_{a}.=-{\rm Tr}\;C_{b}[H,C_{a}]$

This implies $H_{a,b}=-H_{b,a}$, i.e., $H_{a,b}$ is antisymmetric.

One can easily show $[H,c_{1}c_{2}\cdots c_{2N-1}]=0$. Therefore, the matrix
$H_{a,b}$ has the corresponding zero eigenvalue. Since the matrix $H_{a,b}$ is
pure imaginary and antisymmetric with the even dimension $2^{2N-2}$, there
must exist at least one more zero eigenvalue. Namely, there exists an
eigenvector $\alpha_{b}$ such that it satisfies the eigenvalue equation,

$0=\sum_{b}H_{a,b}\alpha_{b}=\sum_{b}(C_{a},[H,C_{b}])\alpha_{b}=(C_{a},[H,\gamma]),$
(B.5)

where we have written

$\gamma=\sum_{b}\alpha_{b}C_{b}.$ (B.6)

The above equation must hold for any $C_{a}$. Therefore, one obtains
$[H,\gamma]=0$, and $\gamma\neq c_{1}c_{2}\cdots c_{2N-1}$.

##  Appendix C Locality of Majorana edge zero mode

In this appendix, we prove the locality of the operator
$U(g)\hat{\gamma}_{0}U(g)^{\dagger}$, which creates the Majorana edge zero
mode above the ground state. From (2.3) and (2.4), it is sufficient to prove
the locality of

$U(g)c_{2\ell-1}U(g)^{\dagger}$

for a not large $\ell\in\\{1,2,\ldots,N\\}$. For this purpose, we will use a
local approximation [28, 24, 27].

To begin with, we recall that the Hamiltonian $H_{g}$ does not contain the
operator $c_{2N}$ at the right edge. This is crucial for the following
argument. Note that

$(c_{2m-1}c_{2m})c_{2m-1}(c_{2m}c_{2m-1})=-c_{2m-1}\quad\mbox{for \
}m=1,2,\ldots,N,$
$(c_{2m-1}c_{2m})c_{2m}(c_{2m}c_{2m-1})=-c_{2m}\quad\mbox{for \
}m=1,2,\ldots,N-1,$
$(c_{2m-1}c_{2N})c_{2m-1}(c_{2N}c_{2m-1})=-c_{2m-1}\quad\mbox{for \
}m=1,2,\ldots,N-1,$

and

$(c_{2m-1}c_{2N})c_{2m-1}c_{2m}(c_{2N}c_{2m-1})=-c_{2m-1}c_{2m}\quad\mbox{for
\ }m=1,2,\ldots,N-1.$

Namely, these unitary transformations change the sign of the operators. Having
these relations in mind, we introduce a unitary transformation,

$\mathcal{U}_{\geq
m}(\sigma_{m},\sigma_{m+1},\ldots,\sigma_{N}):=\mathcal{U}_{m}(\sigma_{m})\mathcal{U}_{m+1}(\sigma_{m+1})\cdots\mathcal{U}_{N}(\sigma_{N}),$

for $m=1,2,\ldots,N$, where

$\mathcal{U}_{m}(\sigma_{m}):=\begin{cases}1,&\sigma_{m}=1;\\\
c_{2m}c_{2m-1},&\sigma_{m}=2;\\\ c_{2N}c_{2m-1},&\sigma_{m}=3;\\\
c_{2m}c_{2m-1}c_{2N}c_{2m-1},&\sigma_{m}=4,\end{cases}$

for $m=1,2,\ldots,N-1$, and

$\mathcal{U}_{N}(\sigma_{N}):=\begin{cases}1,&\sigma_{N}=1;\\\
c_{2N}c_{2N-1},&\sigma_{N}=2.\end{cases}$

By definition, the operator $\mathcal{U}_{\geq
m}(\sigma_{m},\sigma_{m+1},\ldots,\sigma_{N})$ has even fermion parity. For
any operator $A$, we define

$\Pi_{\geq
2m-1}(A):=\frac{1}{2}\left(\frac{1}{4}\right)^{N-m}\sum_{\sigma_{N}=1,2}\;\sum_{\sigma_{m},\sigma_{m+1},\ldots,\sigma_{N-1}}\mathcal{U}_{\geq
m}^{\dagger}(\sigma_{m},\ldots,\sigma_{N})A\mathcal{U}_{\geq
m}(\sigma_{m},\ldots,\sigma_{N}).$ (C.1)

This is the average over the unitary transformations. When an operator $A$ is
written in a product form $A=A^{\prime}A^{\prime\prime}$ which satisfies ${\rm
supp}A^{\prime}\subset\\{1,2,\ldots,2m-2\\}$ and ${\rm
supp}A^{\prime\prime}\subset\\{2m-1,2m,\ldots,2N-1\\}$, one has $\Pi_{\geq
2m-1}(A)=A^{\prime}\Pi_{\geq 2m-1}(A^{\prime\prime})$. If the operator $A\neq
1$ satisfies the condition ${\rm supp}A\subset\\{2m-1,2m,\ldots,2N-1\\}$, then
$\Pi_{\geq 2m-1}(A)=0$. Of course, if $A=1$, then $\Pi_{\geq 2m-1}(A)=1$.

For the operator $U(g)c_{2\ell-1}U(g)^{\dagger}$, the local approximation is
given by

$\Pi_{\geq 2m-1}(U(g)c_{2\ell-1}U(g)^{\dagger}).$

From the definition(C.1) of $\Pi_{\geq 2m-1}(\cdots)$, the difference between
the two operators can be written

$\displaystyle U(g)c_{2\ell-1}U(g)^{\dagger}-\Pi_{\geq
2m-1}(U(g)c_{2\ell-1}U(g)^{\dagger})$ (C.2) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{1}{4}\right)^{N-m}\sum_{\sigma_{m},\ldots,\sigma_{N}}\mathcal{U}_{\geq
m}^{\dagger}(\sigma_{m},\ldots,\sigma_{N})\left[\mathcal{U}_{\geq
m}(\sigma_{m},\ldots,\sigma_{N}),U(g)c_{2\ell-1}U(g)^{\dagger}\right]$

by using the commutator. In order to evaluate the commutator, we recall the
following facts [24]: The unitary operator $U(g)$ is the solution of the
Schrödinger-type equation,

$-i\frac{d}{dg}U(g)=D(g)U(g),$ (C.3)

with the initial condition $U(0)=1$, where the self-adjoint operator $D(g)$ is
written

$D(g)=\int_{-\infty}^{+\infty}dt\;W(t)\;e^{itH_{g}}Ve^{-itH_{g}}$ (C.4)

with a certain real-valued function $W(t)$ which shows subexponential decay
[24]. The operator $D(g)$ is correspond to the Hamiltonian of the Schrödinger-
type equation. The range of the corresponding interaction of $D(g)$ is
determined by the interaction Hamiltonian $V$ of the present Hamiltonian
$H_{g}$ as seen in the right-hand side of $(\ref{Dg})$. When $V$ is finite-
range or exponentially decaying interaction, $D(g)$ shows subexponentially
decaying interaction [24]. Under the assumption on the range of the
interaction $V$, the following Lieb-Robinson bound holds [24, 22]:

$\|[U(g)AU(g)^{\dagger},B]\|\leq\|A\|\|B\|\sum_{i\in{\rm supp}A,j\in{\rm
supp}B}F(|i-j|)$ (C.5)

for any operator $A$ with odd fermion parity, and any operator $B$ with even
fermion parity, where $F$ is a subexponentially decaying function. Here, it is
crucial that at least one of the operators, $A$ and $B$, must have even
fermion parity [29, 27, 22]. Actually, if the two operators, $A$ and $B$,
anticommute with each other for large distance, then the commutator in the
Lieb-Robinson bound must be replaced by the anticommutator [29]. Combining the
Lieb-Robinson bound (C.5) with (C.2), we have

$\displaystyle\|U(g)c_{2\ell-1}U(g)^{\dagger}-\Pi_{\geq
2m-1}(U(g)c_{2\ell-1}U(g)^{\dagger})\|$ (C.6) $\displaystyle\leq$
$\displaystyle\frac{1}{2}\left(\frac{1}{4}\right)^{N-m}\sum_{\sigma_{m},\ldots,\sigma_{N}}\left\|\left[\mathcal{U}_{\geq
m}(\sigma_{m},\ldots,\sigma_{N}),U(g)c_{2\ell-1}U(g)^{\dagger}\right]\right\|$
$\displaystyle\leq$ $\displaystyle\sum_{n\geq
2m-1}F(|2\ell-1-n|)\leq\tilde{F}(|\ell-m|)$

for $\ell\leq m$, where we have used that the operator $\mathcal{U}_{\geq
m}(\sigma_{m},\ldots,\sigma_{N})$ has even fermion parity as mentioned above,
and $\tilde{F}$ is a subexponentially decaying function. This implies the
locality of $U(g)c_{2\ell-1}U(g)^{\dagger}$. In fact, we have the local
decomposition, [24]

$U(g)c_{2\ell-1}U(g)^{\dagger}=\sum_{m=1}^{\infty}\Delta_{m}(c_{2\ell-1}),$
(C.7)

where

$\Delta_{1}(c_{2\ell-1})=\Pi_{\geq 1}(U(g)c_{2\ell-1}U(g)^{\dagger}),$

and

$\Delta_{m}(c_{2\ell-1})=\Pi_{\geq
2m-1}(U(g)c_{2\ell-1}U(g)^{\dagger})-\Pi_{\geq
2m-3}(U(g)c_{2\ell-1}U(g)^{\dagger})$

for $m\geq 2$. From the above bound (C.6), one has

$\displaystyle\|\Delta_{m}(c_{2\ell-1})\|$ $\displaystyle=$
$\displaystyle\|\Pi_{\geq 2m-1}(U(g)c_{2\ell-1}U(g)^{\dagger})-\Pi_{\geq
2m-3}(U(g)c_{2\ell-1}U(g)^{\dagger})\|$ $\displaystyle\leq$
$\displaystyle\|\Pi_{\geq
2m-1}(U(g)c_{2\ell-1}U(g)^{\dagger})-U(g)c_{2\ell-1}U(g)^{\dagger}\|$
$\displaystyle+$ $\displaystyle\|U(g)c_{2\ell-1}U(g)^{\dagger}-\Pi_{\geq
2m-3}(U(g)c_{2\ell-1}U(g)^{\dagger})\|$ $\displaystyle\leq$
$\displaystyle\tilde{F}(|\ell-m|)+\tilde{F}(|\ell-m+1|)$

for a large $m$. This implies a subexponentially decay.

## References

  * [1] A. Yu Kitaev, Unpaired Majorana fermions in quantum wires, Phys.-Usp. 44, 131 (2001).
  * [2] L. Fidkowski and A. Kitaev, Effects of interactions on the topological classification of free fermion systems, Phys. Rev. B 81, 134509 (2010).
  * [3] H. Katsura, D. Schuricht, M. Takahashi, Exact ground states and topological order in interacting Kitaev/Majorana chains, Phys. Rev. B 92, 115137 (2015).
  * [4] P. Fendley, Strong zero modes and eigenstate phase transitions in the XYZ/interacting Majorana chain, J. Phys. A: Math. Theor. 49 (30), 30LT01 (2016).
  * [5] P. Fendley, Parafermionic edge zero modes in $\mathbb{Z}_{n}$-invariant spin chains, J. Stat. Mech. P11020 (2012).
  * [6] A. S. Jermyn, R. S. K. Mong, J. Alicea, and P. Fendley, Stability of zero modes in parafermion chains, Phys. Rev. B 90, 165106 (2014).
  * [7] A. Alexandradinata, N. Regnault, Chen Fang, Matthew J. Gilbert, B. Andrei Bernevig, Parafermionic phases with symmetry-breaking and topological order, Phys. Rev. B 94, 125103 (2016).
  * [8] D. J. Yates, A. G. Abanov, A. Mitra, Lifetime of almost strong edge-mode operators in one dimensional, interacting, symmetry protected topological phases, Phys. Rev. Lett. 124, 206803 (2020).
  * [9] D. J. Yates, A. G. Abanov, A. Mitra, Dynamics of almost strong edge modes in spin chains away from integrability, Phys. Rev. B 102, 195419 (2020).
  * [10] D. V. Else, P. Fendley, J. Kemp, C. Nayak, Prethermal Strong Zero Modes and Topological Qubits, Phys. Rev. X 7, 041062 (2017).
  * [11] W.-T. Xu, G.-M. Zhang, Matrix product states for topological phases with parafermions, Phys. Rev. B 95, 195122 (2017).
  * [12] W.-T. Xu, G.-M. Zhang, Classifying parafermionic gapped phases using matrix product states, Phys. Rev. B 97, 035160 (2018).
  * [13] W.-T. Xu, G.-M. Zhang, Enriched classification of parafermionic gapped phases with time reversal symmetry, Phys. Rev. B 97, 125135 (2018).
  * [14] A. Chew, D. F. Mross, J. Alicea, Fermionized parafermions and symmetry-enriched Majorana modes, Phys. Rev. B 98, 085143 (2018).
  * [15] A. Calzona, T. Meng, M. Sassetti, T. L. Schmidt, $\mathbb{Z}_{4}$ parafermions in one-dimensional fermionic lattices, Phys. Rev. B 98, 201110 (2018).
  * [16] R. W. Bomantara, $\mathbb{Z}_{4}$ parafermion $\pm\pi/2$ modes in interacting periodically driven superconducting chain, Phys. Rev. B 104, 121410 (2021).
  * [17] F. C. Alcaraz, R. A. Pimenta, Free-parafermionic Z(N) and free-fermionic XY quantum chains, Phys. Rev. E 104, 054121 (2021).
  * [18] J. Wouters, F. Hassler, H. Katsura, D. Schuricht Phase diagram of an extended parafermion chain, SciPost Phys. Core 5, 008 (2022).
  * [19] A. Rahmani, M. Franz, Interacting Majorana Fermions, Rep. Prog. Phys. 82, 084501 (2019).
  * [20] M. B. Hastings, The stability of free fermi Hamiltonians, J. Math. Phys. 60, 042201 (2019).
  * [21] W. De Roeck, and M. Salmhofer, Persistence of exponential decay and spectral gaps for interacting fermions, Commun. Math. Phys. 365, 773–796 (2019).
  * [22] T. Koma, Stability of the Spectral Gap for Lattice Fermions, arXiv:2005.04548.
  * [23] T. Kato, Perturbation Theory for Linear Operators, 2nd edn. (Springer, Berlin, Heidelberg, New York, 1980).
  * [24] S. Bachmann, S. Michalakis, B. Nachtergaele, and R. Sims, Automorphic Equivalence within Gapped Phases of Quantum Lattice Systems, Commun. Math. Phys. 309, 835–871 (2012).
  * [25] G. Goldstein, and C. Chamon, Exact Zero Modes in Closed Systems of Interacting Fermions, Phys. Rev. B 86, 115122 (2012).
  * [26] A. Giuliani, V. Mastropietro, and M. Porta, Universality of the Hall Conductivity in Interacting Electron Systems, Commun. Math. Phys. 349, 1107–1161 (2017).
  * [27] B. Nachtergaele, R. Sims, and A. Young, Lieb-Robinson bounds, the spectral flow, and stability of the spectral gap for lattice fermion systems, arXiv:1705.08553.
  * [28] S. Bravyi, M. Hastings, and F. Verstraete, Lieb-Robinson bounds and the generation of correlations and Topological quantum order, Phys. Rev. Lett. 97, 050401 (2006).
  * [29] M. B. Hastings, and T. Koma, Spectral Gap and Exponential Decay of Correlations, Commun. Math. Phys. 265, 781–804 (2006).