# Symmetry-protected gates of Majorana qubits in a high-$T_{c}$ superconductor
platform

Matthew F. Lapa<EMAIL_ADDRESS>Kadanoff Center for Theoretical Physics,
University of Chicago, Chicago, Illinois 60637, USA  Meng Cheng
<EMAIL_ADDRESS>Department of Physics, Yale University, New Haven, CT
06520-8120, U.S.A  Yuxuan Wang<EMAIL_ADDRESS>Department of Physics,
University of Florida, Gainesville, FL 32611, USA

###### Abstract

We propose a platform for braiding Majorana non-Abelian anyons based on a
heterostructure between a $d$-wave high-$T_{c}$ superconductor and a quantum
spin-Hall insulator. It has been recently shown that such a setup for a
quantum spin-Hall insulator leads to a pair of Majorana zero modes at each
corner of the sample, and thus can be regarded as a higher-order topological
superconductor. We show that upon applying a Zeeman field in the region, these
Majorana modes split in space and can be manipulated for braiding processes by
tuning the field and pairing phase. We show that such a setup can achieve full
braiding, exchanging, and phase gates of the Majorana zero modes, all of which
are robust and protected by symmetries. Our analysis naturally includes
interaction effects and can be generalized to cases with fractional bulk
excitations. As many of the ingredients of our proposed platform have been
realized in recent experiments, our results provides a new route toward
topological quantum computation.

###### Contents

  1. I Introduction
  2. II Corner MZMs in a high-$T_{c}$ superconductor platform
    1. II.1 Emergent symmetries
  3. III Corner MZMs from bosonization
    1. III.1 Review of bosonization
    2. III.2 Derivation of the corner modes
  4. IV Symmetry-protected quantum gates of Majorana qubits
    1. IV.1 Manipulating Majorana corner modes via Zeeman field
      1. IV.1.1 Full braid via $2\pi$ rotation of the Zeeman field
      2. IV.1.2 Single exchange via $\pi$ rotation and flip of Zeeman field
      3. IV.1.3 Phase gate via generic rotation of Zeeman field
    2. IV.2 Clifford gates via manipulating Majorana modes within and across corners
  5. V Braiding parafermion modes
  6. VI Conclusion
  7. A Lattice model for second-order topological superconductor
  8. B Mode expansion of the bosonic fields
  9. C Bosonization in the $m=1$ case
    1. C.1 Important identities
    2. C.2 Definition of the fermion operators
    3. C.3 Derivation of anticommutation relations
    4. C.4 Correlation functions in the free theory
  10. D The variational approximation
  11. E Results from the variational approximation (Mathieu’s equation)
    1. E.1 Bound on $|E_{k_{1}}-E_{k_{2}}|$ for $H_{\text{eff}}$ and estimate of the splitting $\Delta E$
    2. E.2 Approximate form of $\chi_{k}(q)$ at large $\lambda$
    3. E.3 Calculating the Berry phase $\gamma_{k}$
  12. F Generalization to $m\neq 1$
  13. G Exact two-fold degeneracy of the domain wall model at $\delta=1/2$
  14. H Ground state degeneracy from the perspective of the ’t Hooft anomaly

##  I Introduction

In the past two decades, topological quantum computation has attracted great
interest in the condensed matter community. They key ingredient of this idea
is to encode and manipulate quantum information using non-Abelian anyons,
which are inherently non-local degrees of freedom and are thus immune to local
error at the hardware level. One of the most promising platform for the
physical realization of non-Abelian anyons is topological superconductors that
host Majorana Zero Modes (MZM) at boundaries and defects Read and Green
(2000); Kitaev (2001); Fu and Kane (2008, 2009); Lutchyn  _et al._ (2010);
Oreg  _et al._ (2010); Beenakker (2013); Leijnse and Flensberg (2012).
Adiabatic braiding and exchange of the MZMs generate Clifford gates in a
topologically protected manner Nayak and Wilczek (1996); Ivanov (2001a), and
implementations of such operations have been studied in various platforms
Alicea  _et al._ (2011); Sau  _et al._ (2011); van Heck  _et al._ (2012);
Hyart  _et al._ (2013); Hegde  _et al._ (2020). However, one drawback of the
Majorana platform is that the Clifford gates are not powerful enough to
achieve universal quantum computation Bravyi and Kitaev (2005). It is well
known that additional gates (i.e., the magic gate of a $\pi/8$ phase rotation)
must be supplemented to achieve universality, which however require non-
topological operations. A number of proposals to implement the magic gate in
the Majorana platform have been put forward Sau  _et al._ (2010); Jiang  _et
al._ (2011); Clarke  _et al._ (2016); Karzig  _et al._ (2016, 2019), most of
which rely on precise control over non-universal couplings to essentially
realize an arbitrary phase rotation (a notable example of a robust magic gate
using geometric decoupling was proposed in Ref. Karzig  _et al._ , 2016,
2019).

Recently, the concept of topological insulators and superconductors has been
generalized to higher-order topological insulators and superconductors.
Schindler  _et al._ (2018, 2018); Benalcazar  _et al._ (2017a, b); Zhang  _et
al._ (2013); You  _et al._ (2018); Kunst  _et al._ (2018); Song  _et al._
(2017); Geier  _et al._ (2018); Wang  _et al._ (2018a); Ezawa (2018); Khalaf
(2018); Matsugatani and Watanabe (2018); Lin and Hughes (2018); Dwivedi  _et
al._ (2018); Langbehn  _et al._ (2017); Parameswaran and Wan (2017); Wang  _et
al._ (2018b); Matsugatani and Watanabe (2018); Trifunovic and Brouwer (2019);
Tiwari  _et al._ (2019); Li  _et al._ (2019); You (2019); Călugăru  _et al._
(2019); Ahn and Yang (2020); Roy (2020); Wu  _et al._ (2020); Zhang  _et al._
(2019, 2020a); Vu  _et al._ (2020); Yan  _et al._ (2018a); Yan (2019); Roberts
_et al._ (2020a, b); Ono  _et al._ (2020); Ghorashi  _et al._ (2019, 2020).
Protected by crystalline symmetries Trifunovic and Brouwer (2019), higher-
order topological superconductors host MZMs at the corners in two spatial
dimensions and Majorana modes at the hinges or vertices in three dimensions.
With the flourishing ideas on the realization of higher-order topological
superconductors, it is natural to search for new possibilities of manipulating
Majorana modes using a higher-order topological superconductor. For example,
in a recent work Pahomi  _et al._ (2020), the authors proposed a protocol
through the manipulation of the Zeeman field and the pairing order parameter,
a full braid (corresponding to $\pi/2$ rotations) between a pair of MZMs can
be achieved (see also Ref. Luo  _et al._ , 2018). In another proposal Zhang
_et al._ (2020b, c), the authors showed that the exchange of MZMs can be
achieved through a multi-step process by tuning three independent Zeeman
fields, a protocol similar to that in a T-junction of superconducting
nanowires. Alicea  _et al._ (2011)

In this work, we propose a different setup in a higher-order topological
superconductor that allows for a much richer set of non-Abelian rotations in
the Hilbert space of Majoranas zero modes, including Clifford and symmetry-
protected phase gates for MZMs. Our proposed setup is based on several recent
works Yan  _et al._ (2018b); Wang  _et al._ (2018c) showing higher-order
topological superconductors can be achieved in a heterostructure involving a
(first-order) topological insulator and unconventional high-$T_{c}$
superconductors coupled via superconducting proximity effect. In particular,
we focus on a heterostructure between a $d$-wave high-$T_{c}$ superconductor,
for example the Bi based cuprate Bi2Sr2CaCu2O8+δ (BSCCO) that has recently
been realized in monolayers Yu  _et al._ (2019), and a quantum spin Hall
insulator, such as WTe2 Wu  _et al._ (2018). For a $d$-wave cuprate
superconductor, the pairing symmetry enforces the proximity-induced gap to
vanish along the certain directions. When such a pairing gap is induced on the
helical edge states of the underlying quantum spin Hall insulator, it creates
a Majorana mass domain wall at each corner, thus hosting _two_ MZMs. (To
achieve this one needs to suppress the single-particle tunneling between the
$d$-wave superconductor and the quantum Hall insulator. We discuss this in
Sec. II.) In the context of higher-order topology, the corner Majorana modes
are protected by mirror reflection symmetries together with time-reversal
symmetry and particle-hole symmetry. The mirror symmetries pin the MZMs at
high-symmetry directions, which form a Kramers pair.

For our purposes, however, the model-specific mirror symmetries are
unnecessary, and in fact intentionally broken by external control fields, so
that the corner MZMs can move along the edge. Instead, we identify two
emergent symmetries, an effective time-reversal and a chiral (an anti-unitary
charge conjugation) symmetry of the low-energy edge theory, which protect the
MZMs even when they are away from the mirror symmetric locations. Using a
bosonized edge theory, we determine the localization length and the excitation
gap of the MZMs in the presence of interaction effects, which are consistent
with the celebrated Kosterlitz-Thouless scaling for infinite systems.
Interestingly, even when the spatial profiles of the MZM become large and
overlap, their degeneracy remain protected by these symmetries. These
additional emergent symmetries also circumvent a no-go theorem Wölms  _et al._
(2016) that would have allowed local time-reversal-invariant perturbations to
spoil the universal non-Abelian Berry phases from braiding a Kramers pair of
MZMs.

The key additional ingredient in our platform is an in-plane Zeeman field. To
this end, we note that recently, heterostructures involving two dimensional
ferromagnets fabricated via molecular-beam epitaxy has already been shown
Kezilebieke  _et al._ (2020) to realize topological superconductivity Takei
_et al._ (2013). In the presence of a Zeeman field, the physical time-reversal
symmetry of the quantum spin Hall insulator is broken. However, we show that
the emergent _effective_ time-reversal and chiral symmetries are still intact,
protecting the MZMs. Since they are no longer Kramers partners, the MZMs can
split spatially. By tuning the Zeeman field, the position of the Majorana
modes can be manipulated. We show that this can be utilized to achieve various
non-Abelian rotations within the degenerate ground state subspace. First, we
show that rotating the in-plane Zeeman field by $2\pi$ is equivalent to a full
braid between the two MZMs, which is analogous to previous proposals. Second,
as the main result of this work, we demonstrate that by taking the in-plane
Zeeman field $\mathbf{B}$ through a “half-moon” contour in the $B_{x}$-$B_{y}$
plane that crosses $B=0$ (see Fig. 3), one can achieve an exchange process of
the two MZMs localized in the same corner, resulting in the hallmark non-
Abelian exchange statistics of the Ising anyons. Crucially, we show that the
non-Abelian Berry phase of this exchange process is protected by the
_physical_ time-reversal symmetry broken only by the Zeeman fields, robust
against local perturbations. Additionally, we show that dual to this process,
one can tune the phase of the complex superconducting order parameter along
one edge of the sample to go through the same “half-moon” contour in the
complex plane, and achieve the exchange of two MZMs from adjacent corners. The
Berry phase during this process is protected by the emergent chiral symmetry.
The combination of these two exchange processes realize the Clifford gates in
a qubit formed by four MZMs in two adjacent corners. Notably, a finite sample
of our setup realizes three qubits, with a set of Clifford gates available on
each edge. Third, we show that by going through a “slice of pie” contour (see
Fig. 4), the Zeeman field (and analogously the superconducting field) can
perform an arbitrary phase gate of the Majorana qubit, including the long-
sought-after “magic gate” for MZMs, crucial for universal topological quantum
computing. Remarkably, the Berry phases in this process are protected by U(1)
symmetries (which can be exact or emergent).

Our proposal has several advantages. First and foremost, the high-$T_{c}$
superconductor platform ensures a higher operating temperature, a larger
critical Zeeman field, and better localization of the Majorana modes. As we
mentioned, BSCCO and WTe2 are readily available 2d materials for $d$-wave
superconductivity and quantum spin Hall effect. Second, our protocols of
exchanging MZMs consist of simple manipulations of Zeeman or pairing fields,
which do not require physically moving around superconducting vortices or
tuning multiple parameters in each exchange process. Third, our setup can
achieve a universal phase gate protected by symmetries, including the $\pi/8$
magic gate, and thus holds promises for universal topological quantum
computation.

While the higher-order topological superconductor platform provides a feasible
realization of our proposal, our results are established within the framework
of the universal effective field theory description of the topological edge
states, which can then be straightforwardly adapted to other systems with the
same low-energy description. For instance, we note that the corner MZMs have
been shown to exist in similar platforms with iron-based high-$T_{c}$
superconductors Wang  _et al._ (2018c); Wu  _et al._ (2020). In general, all
of our results can be easily applied to MZMs realized at domain walls between
magnetic and superconducting regions on the edge of a quantum spin Hall
insulator.

Our analysis can also be directly extended to interacting topological phases
with fractional statistics in the bulk. We show that $\mathbb{Z}_{2m}$
parafermion modes Lindner  _et al._ (2012); Cheng (2012); Clarke  _et al._
(2013); Laubscher  _et al._ (2019, 2020) can be realized using a similar setup
with a fractional quantum spin Hall insulator. A key difference from the
Majorana case is that, here there are $m-1$ independent dynamical phases that
accompanies the non-Abelian Berry phases. Even thoughthe non-Abelian phase is
not topologically protected, for small $m$ it may be possible to precisely
control the time of operation to tune these dynamical phases to zero.
Interestingly, evidence for parafermions have been observed in a similar setup
with fractional quantum Hall states in the presence of superconductivity Gül
_et al._ (2020).

The remainder of this paper is organized as follows. In Sec. II we describe
the setup of our proposed platform that hosts pairs of MZMs at its corners. In
Sec. III we reformulate the derivation of the corner Majorana modes using a
bosonized language, which enables the inclusion of interaction effects and a
transparent interpretation of the non-Abelian Berry phases. As the main result
of this work, in Sec. IV we show that such a setup allows symmetry protected
Clifford gates and phase gates utilizing the MZMs by tuning an in-plane
magnetic field and the phase of the superconducting order parameter. In Sec. V
we generalize our setup to that with a fractional quantum spin Hall insulator
with $m\neq 1$, and show that the Berry phase accumulated using the same
protocol corresponds exactly to the exchange statistics of $\mathbb{Z}_{2m}$
parafermions.

We include various details in the Appendices. In Appendix A we present an
example of a lattice model for the setup that is a higher-order topological
superconductor protected by mirror symmetries and time-reversal symmetry. In
Appendices B, C, D, E, and F we present a detailed analysis of the
bosonization procedure for our setup and the non-Abelian Berry phases we
obtained in a more heuristic manner in the main text. In Appendices G and H we
prove the twofold ground state degeneracy at each corner corresponds to the
Majorana doublet, from both an operator algebra approach and ’t Hooft anomaly
perspective.

##  II Corner MZMs in a high-$T_{c}$ superconductor platform

Our platform is based on several recent proposals Yan  _et al._ (2018b); Wang
_et al._ (2018c) of higher-order topological superconductivity realized in a
heterostructure formed by a quantum spin Hall insulator (QSH) and a
high-$T_{c}$ $d$-wave superconductor, coupled via superconducting proximity
effect. As is well-known, a single-band $d$-wave superconductor hosts gapless
Bogoliubov quasiparticles with Dirac dispersion along the nodal (diagonal)
directions. For our purposes, the single-particle tunneling between the
$d$-wave superconductor and the QSH needs to be suppressed. This can be
achieved by taking advantage of the fact that the single-particle tunneling
and superconducting proximity effect have distinct spatial profiles: the
former effect is peaked at the nodal direction and vanishes along the $x$ and
$y$ directions, while for the latter it is the opposite. Thus, single-particle
tunneling can be effectively suppressed by geometrically separating the
diagonal portion QSH edge with the $d$-wave superconductor. We depict such a
setup in Fig. 1.

For specific lattice models such a phase can be classified as a topological
crystalline superconductor with time-reversal and mirror reflection
symmetries, which we analyze in Appendix A by applying recent results on
higher-order topological phases Trifunovic and Brouwer (2019). However, as we
will see below, our analysis actually does not rely on these symmetries, and
it is more general to start with a low-energy theory describing the edge modes
of a QSH, which we do below. For a full lattice model and its higher-order
topology, we refer the reader to Appendix A.

Figure 1: The schematics of the proposed setup between a QSH insulator and a
$d$-wave superconductor coupled via superconducting proximity effect. By
spatially separating the $d$-wave superconductor and the QSH in the nodal
direction, single-particle tunneling into the QSH edge can be suppressed. The
corner of the sample is subject to a Zeeman field and hosts two MZMs.

The existence of the corner Majorana pairs can be demonstrated by analyzing
the boundary states of the QSH and treating a superconducting gap $\Delta$ and
a Zeeman field $\mathbf{B}$ as perturbations. Consider a portion of the edge
near a corner of a QSH insulator shown in Fig. 1. A low-energy field theory
model of the QSH edge consists of a right-moving fermion $R(x)$ with spin up
(in the $z$-direction) and a left-moving fermion $L(x)$ with spin down (also
in the $z$-direction). These operators have standard anticommutation
relations, for example
$\\{R(x),R^{\dagger}(x^{\prime})\\}=\delta(x-x^{\prime})$, and they also obey
periodic boundary conditions. The kinetic energy for this system takes the
low-energy form

$H_{0}=-i\int dx\
\left(R^{\dagger}(x)\partial_{x}R(x)-L^{\dagger}(x)\partial_{x}L(x)\right)\ .$
(1)

Here we have assumed that both fermions have the same velocity which is set to
$1$.

The superconducting gap term, being a spin-singlet one, takes the low-energy
form on the edge

$H_{\text{SC}}=\int dx\ [\Delta(x)R(x)L(x)+\textrm{h.c.}]\ .$ (2)

Importantly, due to the $d$-wave pairing symmetry, the gap function is odd
under mirror reflection, and when projected onto the edge, $\Delta(x)$ is an
odd function (we choose $x$ to be along the diagonal direction). Finally, the
Zeeman term is projected to the edge as

$H_{\text{Z}}=\int dx\ BR^{\dagger}(x)L(x)+\textrm{h.c.}\ ,$ (3)

where $B\equiv B_{x}+iB_{y}=|B|e^{2i\tau}$. The full Hamiltonian is the sum of
all three of these terms, $H=H_{0}+H_{\text{SC}}+H_{\text{Z}}$.

This Hamiltonian can be diagonalized in the standard way by constructing
lowering operators of the form $\mathcal{O}_{\eta}=\int
dx\,\left\\{\eta_{1}(x)R(x)+\eta_{2}(x)R^{\dagger}(x)+\eta_{3}(x)L(x)+\eta_{4}(x)L^{\dagger}(x)\right\\}$
that satisfy $[H,\mathcal{O}_{\eta}]=-E\mathcal{O}_{\eta}$ with a non-negative
energy $E\geq 0$. Indeed, imposing this relation leads to the usual
Bogoliubov-de Gennes equations for the “spinor”
$\eta(x)=(\eta_{1}(x),\eta_{2}(x),\eta_{3}(x),\eta_{4}(x))^{T}$. Without loss
of generality, taking $\Delta$ and $B$ as real (their constant phases can be
absorbed into the definition of $R$ and $L$), we get

$\left(i\Gamma_{1}\partial_{x}+\Delta(x)\Gamma_{2}+B\Gamma_{13}\right)\eta(x)=E\eta(x)\
,$ (4)

where $\Gamma_{1},\Gamma_{2}$, and $\Gamma_{13}$ are $4\times 4$ matrices
defined as

$\displaystyle\Gamma_{1}$ $\displaystyle=$ $\displaystyle
s_{z}\otimes\mathbb{I}$ (5a) $\displaystyle\Gamma_{2}$ $\displaystyle=$
$\displaystyle s_{y}\otimes\tau_{y}$ (5b) $\displaystyle\Gamma_{13}$
$\displaystyle=$ $\displaystyle s_{x}\otimes\tau_{z}\ ,$ (5c)

and where $s_{x,y,z}$ and $\tau_{x,y,z}$ are the Pauli matrices. We note that
$\\{\Gamma_{1},\Gamma_{2}\\}=\\{\Gamma_{1},\Gamma_{13}\\}=0$, but
$[\Gamma_{2},\Gamma_{13}]=0$. This means that the superconducting and
ferromagnetic mass terms compete with each other.

To analyze this system we can use the fact that $[\Gamma_{2},\Gamma_{13}]=0$
to rotate to a basis in which $\Gamma_{2}$ and $\Gamma_{13}$ are both
diagonal. The required unitary matrix $U$ is given by

$U=\frac{1}{2}\begin{pmatrix}-1&-1&-1&1\\\ -1&1&1&1\\\ -1&-1&1&-1\\\
1&-1&1&1\end{pmatrix}\ .$ (6)

If we define a new spinor $\tilde{\eta}(x)$ via
$\eta(x)=U^{\dagger}\tilde{\eta}(x)$, then we find that $\tilde{\eta}(x)$
satisfies

$\left(i\tilde{\Gamma}_{1}\partial_{x}+\Delta(x)\tilde{\Gamma}_{2}+B\tilde{\Gamma}_{13}\right)\tilde{\eta}(x)=E\tilde{\eta}(x)\
,$ (7)

with

$\displaystyle\tilde{\Gamma}_{1}$ $\displaystyle=$ $\displaystyle
s_{x}\otimes\tau_{z}$ (8a) $\displaystyle\tilde{\Gamma}_{2}$ $\displaystyle=$
$\displaystyle s_{z}\otimes\mathbb{I}$ (8b) $\displaystyle\tilde{\Gamma}_{13}$
$\displaystyle=$ $\displaystyle s_{z}\otimes\tau_{z}\ .$ (8c)

The key property of this new equation for $\tilde{\eta}(x)$ is that it breaks
up into two decoupled $2\times 2$ Dirac equations with masses equal to
$M_{\pm}(x)=\Delta(x)\pm B$. Just as in the Jackiw-Rebbi model, a fermion zero
mode is associated with each domain wall in $M_{+}(x)$ and $M_{-}(x)$, i.e.,
where $\Delta(x)=B$ or $\Delta(x)=-B$ (see Fig. 2). For the profile given by
solid lines in Fig. 2, there is a mass domain wall in $M_{-}(x)$ (marked by
the blue dot to the left), and the zero-energy solution is

$\eta(x)=\frac{1}{2}\begin{pmatrix}i\\\ -i\\\ 1\\\
1\end{pmatrix}e^{-\int_{x_{0}}^{x}ds\,|M_{-}(s)|}$ (9)

where $x_{-}$ is the location where $\Delta=B$. Another MZM, located at
$M_{+}(x)=0$, marked by the blue dot to the right in Fig. 2, can be similarly
obtained. It is straightforward to verify that both solutions for
$\mathcal{O}_{\eta}$ are Hermitian, and correspond to Majorana fermions.
Therefore we find that, for odd $\Delta(x)$, there exist a pair of MZMs
separated by a length $\ell$, the length of the region where $|\Delta|=B$. If
$B=0$, the two MZMs overlap in space, and form a Dirac zero mode. Indeed, they
form a Kramers doublet required by the time-reversal symmetry.

Interestingly, we note that as one tunes the Zeeman field in a given direction
through zero, the two MZMs swaps positions. We can also consider an
alternative configuration in which $\Delta$ is a constant and $B(x)$ changes
sign. This is relevant for a given edge of the QSH with opposite Zeeman field
applied to the two corners it connects, which we will discuss in Sec IV.2. By
the same token, MZMs are nucleated at the nodes of $\Delta\pm B(x)=0$. These
two MZMs switch positions when $\Delta$ is tuned through zero. Later we will
build upon this observation and propose a protocol for non-Abelian braiding of
the Majorana modes.

Figure 2: A spatially varying superconducting mass $\Delta(x)$ (blue curve)
and a constant magnetic field $B$ (solid red line). The dashed red line is the
curve $-B$. For a QSH edge with this configuration of $\Delta(x)$ and $B$,
single Majorana fermions (represented by the blue dots) are localized at the
points where $\Delta(x)=B$ and where $\Delta(x)=-B$. The central region where
$B\geq|\Delta(x)|$ has a length $\ell$.

###  II.1 Emergent symmetries

As we discuss in Appendix A in a specific lattice model, the corner Majorana
modes are protected by crystalline (mirror) symmetries of the bulk theory.
However, these symmetries are rather restrictive for manipulating of the
Majoranas, and in experimental realizations of QSH these symmetries may not be
present anyway. Here we show that fortunately there are several _emergent_
symmetries in the edge theory that protect the MZMs and allow for additional
perturbations to be included.

The edge theory including both the pairing field $\Delta$ and a Zeeman field
$\mathbf{B}$ can be written in first quantized BdG form as

$\displaystyle H=$ $\displaystyle
ks_{z}\tau_{0}+\operatorname{Re}\Delta(x)s_{y}\tau_{y}+\operatorname{Im}\Delta(x)s_{y}\tau_{x}$
$\displaystyle+B_{x}s_{x}\tau_{z}+B_{y}s_{y}\tau_{0},$ (10)

where, e.g. $s_{z}\tau_{0}:=s_{z}\otimes\tau_{0}$. Here $s_{z}=\pm 1$
distinguishes two counter-propagating modes, which transform and couple to
external field like physical spin, which we will refer to as such. We define
$\mathbf{B}=(B_{x},B_{y})$ coupling to edge modes in such a way as the “in-
plane” fields, while it is understood that they may not lie in the plane of
the two-dimensional system.

First we give the full emergent symmetries when $\Delta$ and $B$ are absent.
There are two U(1) subgroups: U(1)s generated by spin:

$U_{s,\phi}=\exp(is_{z}\tau_{z}\phi/2),$ (11)

and U(1)c generated by charge

$\displaystyle U_{c,\theta}=\exp(i\tau_{z}\theta/2).$ (12)

Note that $U_{c,2\pi}=U_{s,2\pi}=-s_{0}\tau_{0}$ is the fermion parity
symmetry. The theory also enjoys a number of discrete symmetries. It is
invariant under the time-reversal symmetry

$\mathcal{T}=is_{y}K,~{}~{}{\mathcal{T}H(k)\mathcal{T}^{-1}=H(-k),}$ (13)

where $K$ is the complex conjugation. In addition, there is a chiral symmetry

$\mathcal{C}=s_{y},~{}~{}{\mathcal{C}H(k)\mathcal{C}^{-1}=-H(k).}$ (14)

As we will see in the next Section, $\mathcal{C}$ is an anti-unitary charge-
conjugation symmetry for many-body states. Ludwig (2015) The BdG Hamiltonian
(10) also has a particle-hole symmetry

$\displaystyle\mathcal{P}=\tau_{x}K,~{}~{}\mathcal{P}H(k)\mathcal{P}^{-1}=-H(-k)^{*},$
(15)

which is not a physical symmetry, but rather a redundancy of the BdG
formalism.

Now we consider the full Hamiltonian (10). The in-plane Zeeman field $B_{x,y}$
breaks the spin rotation and the time-reversal symmetries, but is invariant
under the composite symmetry

$\tilde{\mathcal{T}}\equiv U_{s,\pi}\mathcal{T}=s_{x}\tau_{z}K,$ (16)

which is a symmetry of (10) for a uniform phase of $\Delta$ (taken to be real
without loss of generality). Such an anti-unitary symmetry which squares to
one, along with $\bar{\mathcal{C}}$ places the edge theory in class BDI, which
admits a $\mathbb{Z}$ classification, corresponding to a winding number that
equals the number of symmetry protected MZMs. Therefore, the two MZMs at a
given corner can be viewed as being protected by $\tilde{\mathcal{T}}$. In
addition, we note that $U_{s,2\pi}=-s_{0}\tau_{0}$, the fermion parity, is
obviously still a symmetry.

In the opposite situation in which the magnitude of the unidirectional Zeeman
field (say ${\bf B}=B_{x}\hat{x}$) is spatial dependent and has a domain wall
and $\max(B)>|\Delta|$, two of the Majorana modes from different corners move
to the edge connecting the two corners. In this case, the composite symmetry
$\tilde{\mathcal{T}}$ does not protect MZMs from different corners from
hybridizing (their winding numbers under BDI are opposite), unless additional
symmetries exist. Similar to the spin symmetry, the pairing field breaks both
U(1)c and the chiral symmetry, but preserves their combination

$\displaystyle\tilde{\mathcal{C}}\equiv U_{c,\pi}\mathcal{C}={s_{y}\tau_{z}}.$
(17)

With the composite chiral symmetry $\mathcal{\bar{C}}$, the edge theory
additionally belongs to class AIII, which admits another $\mathbb{Z}$
classification. Such a classification protects the two corner Majorana modes
overlapping on the edge – as can be verified from Eq. (9) and its counterpart
for the other corner, they carry opposite quantum numbers of the unitary
operator $\mathcal{\tilde{C}P\tilde{T}}=s_{z}\tau_{x}$.

So far we have identified the symmetries at the level of the effective BdG
Hamiltonian for the edge states. Typically some of the symmetries are not
exact in the microscopic theory. For instance, in the bulk theory of the QSH,
in general due to the spin-orbit coupling, the spin of edge states $L(x)$ and
$R(x)$ depends on momentum and on location of the edge. However, U(1)s is
realized as an emergent symmetry at low energies as long as the pairing gap is
much greater than the specific spin-orbit coupling that causes momentum and
position dependent spin texture.

Similarly, while ${\mathcal{C}}$ is an exact symmetry of our lattice model for
QSH in Appendix A, a generic QSH insulator is not particle-hole symmetric.
However, for the edge theory, ${\mathcal{C}}$ emerges as an approximate
symmetry as long as the chemical potential is tuned to the crossing point of
the helical edge states. For a 2d system, this can be experimentally achieved
via gating.

Finally, we note that while our analysis of the emergent symmetries for a
given corner and for a given edge appear quite different within the BdG
formalism, as we shall see, within the field-theoretical approach, the
treatments for a given corner and for a given edge are completely symmetric.
In fact, the T-duality of the compact free boson theory, the (1+1)d version
particle-vortex duality Senthil  _et al._ (2019), relates the two symmetries
$\tilde{\mathcal{C}}$ and $\tilde{\mathcal{T}}$.

##  III Corner MZMs from bosonization

We now switch to a bosonization description of the edge of a QSH insulator and
the resulting ground state degeneracy representing the MZMs in the presence of
a proximity SC field. The bosonized treatment has two advantages. First,
interactions can be easily incorporated by turning on a Luttinger parameter
$K\neq 1$ and by generalizing to a fractional QSH (FQSH) state. In particular,
we obtain the scaling behavior of localization length of the Majorana zero
modes upon varying the Luttinger parameter. Second, the calculation of the
non-Abelian Berry phases are rather transparent in the bosonized formalism,
which has an analog of the Berry phases in a 1d lattice.

For the sake of generality, in the bosonized theory we replace the QSH
insulator with a $\nu=1/m$ fractional quantum spin Hall (FQSH) state Levin and
Stern (2009, 2012) and include a Luttinger parameter $K$ to capture
interaction effects. The non-interacting QSH state we have focused on thus far
corresponds to the special case with $m=1$ and $K=1$. We note that in a recent
work Chua  _et al._ (2020), the authors developed a similar bosonization
apporach to Majorana zero modes for a non-interacting open system.

###  III.1 Review of bosonization

The edge of the FQSH state with an emergent $U(1)_{s}$ symmetry can be
described by two bosonic fields $\phi_{\uparrow}(x)$ and
$\phi_{\downarrow}(x)$ obeying the commutation relations

$\displaystyle\left[\phi_{\uparrow}(x),\partial_{x^{\prime}}\phi_{\uparrow}(x^{\prime})\right]$
$\displaystyle=$ $\displaystyle\frac{2\pi i}{m}\delta(x-x^{\prime})$ (18a)
$\displaystyle\left[\phi_{\downarrow}(x),\partial_{x^{\prime}}\phi_{\downarrow}(x)\right]$
$\displaystyle=$ $\displaystyle-\frac{2\pi i}{m}\delta(x-x^{\prime})$ (18b)
$\displaystyle\left[\phi_{\uparrow}(x),\phi_{\downarrow}(x^{\prime})\right]$
$\displaystyle=$ $\displaystyle 0\ .$ (18c)

In the K-matrix formalism for edges of Chern-Simons theories, this system
corresponds to the matrix $K=m\sigma^{z}$. Both fields
$\phi_{\uparrow/\downarrow}(x)$ are defined to have compactification radius
$2\pi$. This means that all physical operators must be invariant under the
shift $\phi_{\uparrow}(x)\to\phi_{\uparrow}(x)+2\pi$, and likewise for
$\phi_{\downarrow}(x)$. Then the allowed operators containing zero derivatives
of these fields must be built from exponentials of the form
$e^{in\phi_{\uparrow/\downarrow}(x)}$ for some integer $n\in\mathbb{Z}$.

The charge density current, and spin operator for the edge are defined to be

$\displaystyle\rho(x)=$
$\displaystyle\frac{1}{2\pi}\left(\partial_{x}\phi_{\uparrow}(x)+\partial_{x}\phi_{\downarrow}(x)\right)$
$\displaystyle j(x)=$
$\displaystyle\frac{1}{2\pi}\left(\partial_{x}\phi_{\uparrow}(x)-\partial_{x}\phi_{\downarrow}(x)\right)$
$\displaystyle s(x)=$
$\displaystyle\frac{1}{4\pi}\left(\partial_{x}\phi_{\uparrow}(x)-\partial_{x}\phi_{\downarrow}(x)\right)\
.$ (19)

We then find that the right- and left-moving electron operators for the free
fermion case are given by

$\displaystyle R(x)$ $\displaystyle\sim$
$\displaystyle\frac{1}{\sqrt{\ell}}:e^{-im\phi_{\uparrow}(x)}:$ (20a)
$\displaystyle L(x)$ $\displaystyle\sim$
$\displaystyle\frac{1}{\sqrt{\ell}}:e^{im\phi_{\downarrow}(x)}:\ ,$ (20b)

where $\ell$ is the length of the system. We present the details of the
bosonization dictionary in Appendices B and C for $m=1$ and in Appendix F for
$m\neq 1$.

With this definition we find that acting with $R(x)$ or $L(x)$ lowers the
total charge by one unit, as expected for an operator that annihilates a
single electron. In addition, the anticommutation relation
$\\{R(x),R(x^{\prime})\\}=0$ (which should be obeyed by any fermionic
operator) follows from the fact that $m$ is an odd integer.

The basic kinetic energy term for the bosonic fields
$\phi_{\uparrow/\downarrow}(x)$ takes the form

$H_{0}=\frac{1}{2\pi}\int dx\
\left[\frac{1}{2}\sum_{\sigma=\uparrow,\downarrow}(\partial_{x}\phi_{\sigma}(x))^{2}+g\partial_{x}\phi_{\uparrow}(x)\partial_{x}\phi_{\downarrow}(x)\right]\
.$ (21)

Here we have also incorporated a density-density interaction (with coupling
constant $g$) between the spin up and spin down fermions111Note that
$\frac{1}{2\pi}\partial_{x}\phi_{\sigma}(x)$ is the density operator for
excitations with spin $\sigma\in\\{\uparrow,\downarrow\\}$.. We see that $g>0$
corresponds to a repulsive interaction, while $g<0$ corresponds to an
attractive interaction.

For the domain wall configurations that we study in this paper, it is
convenient to introduce new non-chiral fields $\varphi(x)$ and $\vartheta(x)$
defined as

$\displaystyle\varphi(x)$ $\displaystyle=$
$\displaystyle\frac{m}{2}\left(\phi_{\uparrow}(x)+\phi_{\downarrow}(x)\right)$
(22a) $\displaystyle\vartheta(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\phi_{\uparrow}(x)-\phi_{\downarrow}(x)\right)\
,$ (22b)

which satisfy the commutation relation

$\displaystyle[\varphi(x),\partial_{x^{\prime}}\vartheta(x^{\prime})]=\pi
i\delta(x-x^{\prime}).$ (23)

In terms of these fields we find that

$\displaystyle\rho(x)=$ $\displaystyle\frac{1}{\pi m}\partial_{x}\varphi(x)$
$\displaystyle j(x)=$ $\displaystyle\frac{1}{\pi}\partial_{x}\vartheta(x)$
$\displaystyle s(x)=$ $\displaystyle\frac{1}{2\pi}\partial_{x}\vartheta(x),$
(24)

and that $H_{0}$ can be rewritten in the form

$H_{0}=\frac{v}{2\pi}\int dx\
\left[\frac{1}{mK}(\partial_{x}\varphi(x))^{2}+mK(\partial_{x}\vartheta(x))^{2}\right]\
,$ (25)

where the renormalized velocity $v$ and Luttinger parameter $K$ are related to
the coupling constant $g$ as

$\displaystyle v$ $\displaystyle=$ $\displaystyle\frac{1}{m}\sqrt{1-g^{2}}$
(26a) $\displaystyle K$ $\displaystyle=$ $\displaystyle\sqrt{\frac{1-g}{1+g}}\
.$ (26b)

Note the singularity in $K$ at $g=-1$ and the zeros in $v$ and $K$ at $g=1$.
In addition, we have $K<1$ for repulsive interactions and $K>1$ for attractive
interactions, while $K=1$ in the absence of interactions ($g=0$). For later
use it is convenient to combine $m$ and $K$ into a modified Luttinger
parameter

$K^{\prime}=mK\ ,$ (27)

as it is this modified Luttinger parameter that actually appears in $H_{0}$.

In this bosonized formalism, a superconducting mass term takes the form

$\displaystyle\Delta R^{\dagger}(x)L^{\dagger}(x)+\text{h.c.}{\propto}$
$\displaystyle\cos\left[2m\vartheta(x)+2\rho\right]\ ,$ (28)

where $\rho$ is the superconducting phase. Similarly, a ferromagnetic mass
term takes the form

$\displaystyle(B_{x}+iB_{y})R^{\dagger}(x)L(x)+\text{h.c.}{\propto}$
$\displaystyle\cos\left[2\varphi(x)+2\tau\right],$ (29)

where $B_{x}+iB_{y}=|B|e^{2i\tau}$. A more rigorous derivation of the mass
terms in terms of boson fields is done using the mode expansion, as we
describe below and in Appendix D.

It is instructive to see how the bosonic variables $\varphi$ and $\vartheta$
transform under the emergent symmetries identified in the previous section:

$\displaystyle U_{s,\phi}:~{}$
$\displaystyle\vartheta\to\vartheta,~{}\varphi\to\varphi-\phi/2$
$\displaystyle U_{c,\theta}:~{}$
$\displaystyle\varphi\to\varphi,~{}\vartheta\to\vartheta-\theta/2$
$\displaystyle\mathcal{T}:~{}$
$\displaystyle\vartheta\to-\vartheta,~{}\varphi\to\varphi+\pi/2$
$\displaystyle\mathcal{C}:~{}$
$\displaystyle\varphi\to-\varphi,~{}{\vartheta\to\vartheta-\pi/2}.$ (30)

Despite their different forms for the BdG Hamiltonian, at the field theory
level both $\mathcal{C}$ and $\mathcal{T}$ are antiunitary symmetries, since
each flips the sign on one of the dual fields $\vartheta$ and $\varphi$. In
particular $\mathcal{C}$ flips the sign of the charge but not the current,
which is thus an anti-unitary charge conjugation symmetry.

The composite symmetries $\tilde{\mathcal{T}}=U_{s,\pi}\mathcal{T}$ and
$\tilde{\mathcal{C}}=U_{c,\pi}\mathcal{C}$ now become

$\displaystyle\tilde{\mathcal{T}}:~{}$
$\displaystyle\vartheta\to-\vartheta,~{}\varphi\to\varphi$
$\displaystyle\tilde{\mathcal{C}}:~{}$
$\displaystyle\varphi\to-\varphi,~{}\vartheta\to\vartheta.$ (31)

Thus under $\mathcal{\tilde{T}}$ ($\mathcal{\tilde{C}}$), the Zeeman (SC) term
remains invariant. Under the T-duality of the free boson theory
$\varphi\leftrightarrow\vartheta$ (which is an emergent symmetry when $K=1$),
$\tilde{\mathcal{T}}$ and $\tilde{\mathcal{C}}$ are exchanged, as well as the
Zeeman and the SC terms.

In our analysis below we will mainly invoke $U_{s,\phi}$,
$\mathcal{\tilde{T}}$ for the MZMs localized at a given corner, and due to the
T-duality, the results directly carry over to the case of MZMs overlapping on
an edge.

###  III.2 Derivation of the corner modes

We analyze the states hosted by a corner region with Zeeman fields sandwiched
between two superconducting regions with opposite pairing gap related by the
$d$-wave symmetry. To simplify the calculations, let us fix the length of the
magnetic region to be $\ell$, within which the Zeeman field is a constant, and
take the limit in which the superconducting gap $|\Delta|\to\infty$ in the
superconducting region and thus the $\vartheta$ at the two ends of the
magnetic region are completely pinned. Without loss of generality, we take

$\vartheta(0)=0\mod\frac{\pi}{m},~{}~{}~{}~{}\vartheta(\ell)=\frac{\pi}{2m}\mod\frac{\pi}{m}.$
(32)

In the magnetic region, the Hamiltonian is given by

$\displaystyle H=\int_{0}^{\mathcal{\ell}}dx\ \Big{\\{}\frac{v}{2\pi}$
$\displaystyle\left[\frac{1}{mK}(\partial_{x}\varphi(x))^{2}+mK(\partial_{x}\vartheta(x))^{2}\right]$
$\displaystyle+b\cos[2\varphi(x)+2\tau]\Big{\\}}\ ,$ (33)

where $b$ is a coupling constant induced by $B$. To analyze the low-energy
spectrum of this Hamiltonian it is helpful to perform a mode expansions for
$\varphi$ and $\vartheta$:

$\displaystyle\varphi(x)$
$\displaystyle=mq-\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}\cos(\kappa_{n}x)(b_{n}+b^{\dagger}_{n})$ (34a)
$\displaystyle\vartheta(x)$ $\displaystyle=\left(p+\frac{1}{2}\right)\frac{\pi
x}{m\ell}$ $\displaystyle+i\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}\sin(\kappa_{n}x)(b_{n}-b^{\dagger}_{n})\ ,$ (34b)

where $\kappa_{n}=\frac{\pi n}{\ell}$, and where we included the dimensionless
ultraviolet cutoff $\epsilon$ to control the oscillator sums. Removing the
cutoff corresponds to taking $\epsilon\to 0$, and one can check that in this
limit the fields obey the correct commutation relations in Eq. (23) (see
Appendix C.3). One can also see that the field $\vartheta(x)$ obeys the
boundary conditions from Eq. (32) with quantized winding number
$p\in\mathbb{Z}$. Importantly, from Eq. (23) we have the commutation relation

$\left[q,p\right]=i,~{}[b_{n},b_{n}^{\dagger}]=1,$ (35)

indicating $q$ and $p$ are conjugate variables. As a result of the
quantization of $p$, $q$ is a compact variable with $q\sim q+2\pi$.

In the absence of the Zeeman field, the eigenstates of $H_{0}$ are labeled by
the winding number $p$ and the occupation number for a set of new
quasiparticle modes:

$H_{0}=\frac{\pi
vK}{2m\mathcal{\ell}}\left(p+\frac{1}{2}\right)^{2}+v\sum_{n}\kappa_{n}a_{n}^{\dagger}a_{n}+\text{const.}$
(36)

The operators $\\{a_{n}\\}$ are related to $\\{b_{n}\\}$ via a Bogoliubov
transformation

$a_{n}=\cosh(\eta)b_{n}+\sinh(\eta)b^{\dagger}_{n}\ ,~{}~{}~{}e^{-2\eta}=K.$
(37)

which we discuss in details in Appendix D. The Fock space structure is
guaranteed by $\mathcal{\tilde{T}}$ symmetry; for example a symmetry breaking
Zeeman term $\sim B_{z}\partial_{x}\vartheta$ term would condense the
quasiparticles in the ground state.

The Zeeman field term makes the quasiparticles massive, and further increases
the quasiparticle gap. Therefore, for low-energy states we only need to focus
on the Fock vacuum sector and consider the $q$ and $p$ modes. The effective
Hamitonian is given by

$H_{\rm eff}=\alpha\tilde{p}^{2}-\beta\cos(2mq+2\tau)$ (38)

in which $\tilde{p}\equiv p+1/2$, and $\alpha,\beta$ are coupling constants
renormalized by quasiparticle fluctuations. In Appendices D and F, using a
variational approximation, we derive the coefficients for $K^{\prime}<2$, and
for the Zeeman field within the range

$\frac{1}{a}\left(\frac{a}{\ell}\right)^{2-K^{\prime}}\ll B\ll\frac{1}{a},$
(39)

where $a$ is a short-distance cutoff, e.g., given by the underlying lattice.
The results are given by

$\displaystyle\alpha$ $\displaystyle=\frac{vK^{\prime}}{2\pi m^{2}\ell}$
$\displaystyle\beta$ $\displaystyle\sim
B^{\frac{2}{2-K^{\prime}}}a^{\frac{2K^{\prime}-2}{2-K^{\prime}}}\ell.$ (40)

Using Eq. (39), it is straightforward to see that here
$\beta\gg\alpha$.222Interestingly, we note that for the free fermion case with
$K=1$, the prefactor $\beta$ of the cosine term is actually proportional to
$B^{2}$ rather than $B$.

The Schrödinger equation with the Hamiltonian in Eq. (38) is known as the
Mathieu’s equation Weinstein and Keller (1985); Dunne and Unsal (2016); Simon
(1984), and can be viewed as the equation of motion of a single particle in a
1d ring modeled by a periodic lattice potential. According to Bloch’s theorem,
the eigenstates $\ket{\psi_{k}}$ of this Hamiltonian is labeled by lattice
momenta $k$ 333This is not to be confused with the _actual_ lattice momenta of
our QSH/$d$-SC system., i.e.
$e^{i\frac{\pi}{m}\tilde{p}}\ket{\psi_{k}}=e^{i\frac{\pi}{m}k}\ket{\psi_{k}}$,
which take quantized values inside the Brillouin zone. The lattice constant is
$\pi/m$, and thus $k\in[-m,m)$. Since $\tilde{p}$ takes half-integer quantized
values, so does $k$. The offset $1/2$ in the quantization of $\tilde{p}$ is
analogous to the effect of a magnetic flux through the lattice ring, causing a
twisted boundary condition and the same amount of offset in the lattice
momenta $k$. Therefore, we have

$k\in(\mathbb{Z}+1/2)\cap(-m,m).$ (41)

Under $\mathcal{\tilde{T}}$, $\tilde{p},k\to-\tilde{p},-k$. From Eqs. (24) and
(34), we see that physically the tunneling current and spin quantum numbers
are directly related to the $k$ via

$\displaystyle
J=2S=\int_{0}^{\ell}\frac{dx}{m\pi}\partial_{x}\vartheta(x)=\frac{k}{m}\mod\frac{1}{m}.$
(42)

The eigenstates of this Hamiltonian form energy bands labeled by the band
index and lattice momenta $\\{k\\}$. Each band consists of $2m$ states. We
will focus on the lowest band.

In the limit of large $\ell$ or large $B$, we have $\beta\gg\alpha$ and the 1d
lattice is in a flat-band limit. As we show in Appendix E.1, in the limit
$\beta\gg\alpha$ (which is the same as Eq. (39)) the bandwidth is
exponentially suppressed as

$\Delta E\sim\mathrm{const}\times\exp\left(-\frac{\ell}{\xi}\right).$ (43)

Here the correlation length $\xi$ for $K^{\prime}<2$ is expressed as

$\xi=\sqrt{\frac{\alpha}{\beta}}\ell\sim\left(\frac{1}{Ba}\right)^{\frac{1}{2-K^{\prime}}}a,$
(44)

following the familiar Kosterlitz-Thouless scaling behavior. The $2m$ states
in the lowest band are approximately degenerate. This is the same $2m$
degeneracy given by a pair of $\mathbb{Z}_{2m}$ parafermions Lindner  _et al._
(2012); Cheng (2012); Vaezi (2013); Clarke  _et al._ (2013); Barkeshli  _et
al._ (2013a, b).

In particular for $m=1$, the ground state degeneracy corresponds to a pair of
MZMs, consistent with what we found using the BdG formalism. The exponential
supression of the hybridization energy of the parafermions indicates that
these modes are exponentially localized in space, consistent with the results
we obtained for the free fermion case. Indeed, for $m=1,K=1$ we restore the
familiar result $\xi\sim 1/B$; see Eq. (9).

In the flat band limit $\beta\gg\alpha$, the gap separating the MZMs from
excited states can be obtained by approximating Eq. (38) by expanding the
cosine potential. We find for the excitation gap

$\Delta_{\rm ex}\sim\sqrt{\alpha\beta}=\sqrt{K}B(Ba)^{\frac{K-1}{2-K}}.$ (45)

For the free fermion system, this expression reduces to the Zeeman energy
$\sim B$.

In the opposite limit $\alpha\gg\beta$, the Mathieu’s equation is in the weak
potential limit, and the band dispersion is similar to that of a free
particle. In particular, as $B$ goes to zero, the band gap closes at the BZ
boundary ($\pm 1$ for $m=1$) and the spectrum restores the parabolic
dispersion. If either of the states in the lowest bands reside at the
Brillouin zone boundary, the MZMs will be “poisoned” by excited states.

Fortunately, in the presence of the twist boundary condition causing
$k\in\mathbb{Z}+1/2$, the quantized lattice momentum $k$ do not take values at
the BZ boundary, and the Majorana states in the lowest band remain degenerate
and separated from higher bands. This is consistent with our findings in the
previous Section. The size of their spatial profile is $\mathcal{O}(\ell)$.
The analysis here further shows that the excitation gap results from the
quantization of lattice momentum, and is given by

$\Delta_{\rm ex}\sim v/\ell.$ (46)

We prove the two-fold degeneracy for a general $\beta/\alpha$ more rigorously
in Appendix G.

Interestingly, such a robust ground state degeneracy has been elucidated
Gaiotto  _et al._ (2017) from the perspective of a mixed ’t Hooft anomaly
between time-reversal symmetry (corresponding to our generalized time-reversal
$\tilde{\mathcal{T}}$) and fermion parity symmetry $\varphi\to\varphi+\pi$
(generated by our $U_{s,2\pi}$) of field theories with a $\Theta$-term at
$\Theta=\pi$. The anomaly ensures that independent of basis choice, one of the
two classical symmetries is represented as a double cover at the quantum
level, leading to the two-fold degeneracy. The two states are related by time-
reversal and differ by fermion parity quantum numbers. We present the proof of
the degeneracy in Appendices H in a way that reveals a clear analogy with
Appendix D of Ref. Gaiotto  _et al._ , 2017.

We end this section by noting that in the alternative configuration when the
Zeeman fields at two different corners are antiparallel, the two MZMs can be
obtained in a dual bosonized theory, related to our discussion above by
$\vartheta\leftrightarrow\varphi$ and
$\mathcal{\tilde{T}}\leftrightarrow\mathcal{\tilde{C}}$.

##  IV Symmetry-protected quantum gates of Majorana qubits

In this Section we focus on the $m=1$ case where the degenerate corner states
correspond to a pair of MZMs. As we showed in the fermionic language, these
two MZMs are located at $\Delta(x)=\pm B$ and switch position when the Zeeman
field is flipped. We now show that when the in-plane Zeeman field is _rotated_
back, the full process induces an non-Abelian Berry phase that is the same as
exchanging two-dimensional MZMs (or Ising anyons). Furthermore, we show that
by tuning the Zeeman field as well as the superconducting order parameter, one
can realize all Clifford gates and universal phase gates on the ground state
qubit, protected by the emergent symmetries of the theory.

As is well-known, with a fixed fermion parity, four MZMs form a two-level
qubit. This is realized by two adjacent corners of the higher-order
topological superconductor platform. The Clifford gates consist of exchanging
the Majorana modes both within the same corner and across different corners,
the protocol of which we discuss below.

###  IV.1 Manipulating Majorana corner modes via Zeeman field

####  IV.1.1 Full braid via $2\pi$ rotation of the Zeeman field

Before we discuss the exchange process of the Majorana modes, let us first
consider an adiabatic process involving a full $2\pi$ rotation of the in-plane
Zeeman field and compute the Berry phase. In our conventions this corresponds
to keeping the magnitude $B$ of the magnetic field fixed while tuning the
angular parameter $\tau$ from $0$ to $\pi$ in Eq. (38). We will show that this
correspond to a full braid of two Majoranas at a given corner.

Let $|\psi_{k}(B,\tau)\rangle$ be the ground state of $H_{\rm eff}(B,\tau)$
with lattice momentum $k$. According to Eq. (41), $k=\pm 1/2$. In other words,

$e^{i\pi\tilde{p}}|\psi_{k}(B,\tau)\rangle=e^{i\pi
k}|\psi_{k}(B,\tau)\rangle=\pm i|\psi_{k}(B,\tau)\rangle.$ (47)

In the 1d lattice interpretation, the parameter rotation angle of the Zeeman
field $2\tau$ corresponds to a displacement of the periodic potential by an
amount of $\tau$, and thus the eigenstate can be expressed via a translation
operator:

$|\psi_{k}(B,\tau)\rangle\sim e^{-i\tau\tilde{p}}|\psi_{k}(B,0)\rangle.$ (48)

However, from Eq. (47), the right hand side is not single-valued upon a full
$2\pi$ rotation $(\tau=\pi)$ of the Zeeman field. For an unambiguous
calculation of the Berry phase, one should choose the phases of the states
$|\psi_{k}(B,\tau)\rangle$ so that these states are single-valued functions of
$B$ and $\tau$ (defined modulo $\pi$) in the region of the parameter space
that is of interest for the Berry phase calculation.

This issue can be addressed by adding to each ground state a c-number phase
factor $e^{i\tau k}$ corresponding to their lattice momenta

$|\psi_{k}(B,\tau)\rangle=e^{i{\tau
k}}e^{-i{\tau\tilde{p}}}|\psi_{k}(B,0)\rangle\ .$ (49)

From Eq. (47), this choice ensures that $|\psi_{k}(B,\tau)\rangle$ returns to
itself when we wind the Zeeman field by $2\pi$ (translating $\tau$ by $\pi$),
i.e., we have

$|\psi_{k}(B,\tau+\pi)\rangle=|\psi_{k}(B,\tau)\rangle\ .$ (50)

This is not the only possible choice, and it is well-known that the Berry
phases that we obtain are invariant under any redefinition
$|\psi_{k}(B,\tau)\rangle\to|\tilde{\psi}_{k}(B,\tau)\rangle=e^{i\theta_{k}(B,\tau)}|\psi_{k}(B,\tau)\rangle$,
provided that the new states $|\tilde{\psi}_{k}(B,\tau)\rangle$ are also
single-valued functions of $B$ and $\tau$ in the relevant region of the
parameter space.

We now compute the Berry phases $\gamma_{k}$ picked up by the states
$|\psi_{k}(B,\tau)\rangle$ during the $2\pi$ rotation of the Zeeman field. The
Berry phases $\gamma_{k}$ are given by the standard formula

$\gamma_{k}=i\int_{0}^{\pi}d\tau\
\langle\psi_{k}(B,\tau)|\partial_{\tau}|\psi_{k}(B,\tau)\rangle\ ,$ (51)

and so using Eqs. (49) we find that

$\gamma_{k}={\pi}\langle\psi_{k}(B,0)|\tilde{p}|\psi_{k}(B,0)\rangle-{\pi k}\
.$ (52)

Intuitively, the first term evaluates the average momentum of the Bloch
states, and the second term lattice momentum. As we mentioned, depending on
the depth of the periodic potential in Eq. (38) there are two important limits
– the tight-binding limit ($\beta\gg\alpha$) and weak periodic potential limit
($\alpha\gg\beta$). In the first limit, Bloch states are approximately a
linear superposition of bound states each at a potential minimum, while in the
second limit, Bloch states are approximately plane-wave states. Therefore,
heuristically we find that in former limit the average momentum
$\langle\tilde{p}\rangle$ approaches zero, while in the latter
$\langle\tilde{p}\rangle$ approaches the lattice momentum $k$. Thus we have in
the tight-binding limit,

$\displaystyle\gamma_{k}=-\pi k,~{}~{}~{}k=\pm\frac{1}{2},$ (53)

which can be understood as coming from “dragging” the Bloch state by a lattice
constant. In the opposite limit, the Berry phase vanishes, i.e., the lattice
potential is so weak that translating it does not induce a significant change
in the wave function.

Here focus on the tight-binding limit ($\beta\gg\alpha$). This is the same
condition as (39), i.e., $Ba\gg(a/\ell)^{2-K}$. In Appendix E.3 we directly
compute the Berry phase for the state $|\psi_{k}(B,0)\rangle$ using our
analysis of that state based on Mathieu’s equation. There we show that the
deviations of the Berry phase from the approximate result in Eq. (53) is
indeed exponentially small in $\ell$. This result is topological, in the sense
that the Berry phase does not depend on parameters such as $B$ and $K$ up to
exponentially small corrections.

Noting that $k=\pm 1/2$, the result in (53) matches exactly the non-Abelian
Berry phases accrued during a full $2\pi$ braiding of two MZMs. Ivanov (2001b)
We note that the full braiding of the Majorana modes have also been proposed
in a similar higher-order topological superconductor platform in Ref. Pahomi
_et al._ (2020).

####  IV.1.2 Single exchange via $\pi$ rotation and flip of Zeeman field

We now show that owing to symmetries of the system, one can also perform a
_single exchange_ of two Majoranas using a different adiabatic process that
also involves only the external Zeeman field. To motivate this process, recall
from the previous subsection that the Berry phase for a $2\pi$ rotation of the
Zeeman field within the $x$-$y$ plane is equal (at large $\ell$) to the Berry
phase for a full braid (double exchange) of the fractional quasiparticles
localized near the ends of the FM region. In this subsection we show that this
Berry phase, and the adiabatic process itself, can be split into two equal
contributions in a symmetry-protected manner, such that each contribution on
its own yields the Berry phase for a single exchange of fractional
quasiparticles. The Berry phase $\tilde{\gamma}_{k}$ for this process is then
given by _half_ of the value $\gamma_{k}$ for the full braid,

$\tilde{\gamma}_{k}=-\frac{\pi k}{2},~{}~{}~{}k=\pm\frac{1}{2}\ .$ (54)

We find that $\tilde{\gamma}_{\pm 1/2}$ for the two ground states differ by
$\pi/2$, and this is exactly the relative Berry phase expected for a single
exchange of two MZMs Ivanov (2001a).

To achieve this, we consider “half-moon” paths of the Zeeman field
$\mathbf{B}=(B_{x},B_{y})$ denoted in Fig. 3. This path consists of a half
circle from $\tau=0$ to $\tau=\pi/2$, and a straight line with $B_{y}=0$
sweeping the $\mathbf{B}$ field back to the initial configuration passing
through the origin. Crucially, along the arc $\mathbf{B}$ must be in the
tight-binding regime, i.e., $Ba\gg(a/\ell)^{2-K}$. As we shall see below, the
precise shape of the arc does not matter, as long as it is in the flat-band
limit. We denote such a contour transversed counterclockwise by
$\mathcal{C}_{B}^{\prime}$, and its image under inversion in the $B_{x}-B_{y}$
plane, transversed clockwise, by $\mathcal{C}_{B}^{\prime\prime}$.

Let $|\psi_{k}(\mathbf{B})\rangle=|\psi_{k}(B_{x},B_{y})\rangle$ be the ground
state of the Hamiltonian in the sector with “lattice momentum” $k$. In our
previous notation we had $B_{x}=B\cos(2\tau)$ and $B_{y}=B\sin(2\tau)$, and so
$|\psi_{k}(\mathbf{B})\rangle$ can be identified with the state
$|\psi_{k}(B,\tau)\rangle$ that we defined in Eq. (49). The Berry phase
$\gamma_{k}$ for the full $2\pi$ rotation of $\mathbf{B}$ can be written as
the line integral

$\gamma_{k}=i\oint_{\mathcal{C}_{B}}d\mathbf{B}\cdot\langle\psi_{k}(\mathbf{B})|\nabla_{\mathbf{B}}|\Psi_{k}(\mathbf{B})\rangle\
,$ (55)

where $\mathcal{C}_{B}$ is the circular contour of radius $B$ centered at the
origin of the $B_{x}$-$B_{y}$ plane.

The integral expression for $\gamma_{k}$ can be split into two contributions
as

$\displaystyle\gamma_{k}$
$\displaystyle=i\oint_{\mathcal{C}^{\prime}_{B}}d\mathbf{B}\cdot\langle\psi_{k}(\mathbf{B})|\nabla_{\mathbf{B}}|\Psi_{k}(\mathbf{B})\rangle$
$\displaystyle\quad+i\oint_{\mathcal{C}^{\prime\prime}_{B}}d\mathbf{B}\cdot\langle\psi_{k}(\mathbf{B})|\nabla_{\mathbf{B}}|\Psi_{k}(\mathbf{B})\rangle\
$ $\displaystyle\equiv\gamma_{k}^{\prime}+\gamma^{\prime\prime}_{k}.$ (56)

Here $\gamma_{k}^{\prime}$ and $\gamma_{k}^{\prime\prime}$ are the
contributions to the total Berry phase from the half-moon paths
$\mathcal{C}^{\prime}_{B}$ and $\mathcal{C}^{\prime\prime}_{B}$, respectively.

Figure 3: The “half-moon”-shaped contours $\mathcal{C}^{\prime}_{B}$ (blue)
and $\mathcal{C}^{\prime\prime}_{B}$ (red) in the $B_{x}$-$B_{y}$ plane. The
magnetic field in the semicircle portion is in the tight-binding regime
$Ba\gg(a/\ell)^{2-K}$. The Berry phase for either one of these paths is equal
(up to exponentially small corrections) to the known Berry phase for a single
braid of fractional quasiparticles.

An important prerequisite for a well-defined Berry phase is that the system
remains gapped during the process. This is indeed true for the half-moon
contour, since the ground state qubit is always energetically separated from
the excited states: on the outer arc, the corner region is gapped by the
Zeeman field (see Eq. (45)). Near the origin, the Zeeman field vanishes but a
finite size gap still exists (see Eq. (46)). In addition, in order for the
dynamical phases to cancel, the two ground states should remain degenerate,
which is guaranteed by the $\mathcal{\tilde{T}}$ symmetry for any value of
$B$, including at $B=0$.

Provided that the Zeeman field is the only odd component under the
$\mathcal{T}$ symmetry (the physical time-reversal symmetry preserved by the
quantum spin Hall and $d$-wave superconductor; not to be confused with
$\tilde{\mathcal{T}}$), the contour
$\mathcal{C}^{\prime}_{B}\to-\mathcal{C}^{\prime\prime}_{B}$ under
$\mathcal{T}$ (which reverses both the Zeeman field and the orientation of the
contour). The Berry phase, is obviously odd under time-reversal, and
therefore,

$\displaystyle\gamma_{k}^{\prime}-\gamma_{k}^{\prime\prime}=$ $\displaystyle
i\oint_{\mathcal{C}^{\prime}_{B}}d\mathbf{B}\cdot\langle\psi_{k}(\mathbf{B})|\nabla_{\mathbf{B}}|\Psi_{k}(\mathbf{B})\rangle$
$\displaystyle+i\oint_{-\mathcal{C}^{\prime\prime}_{B}}d\mathbf{B}\cdot\langle\psi_{k}(\mathbf{B})|\nabla_{\mathbf{B}}|\Psi_{k}(\mathbf{B})\rangle\
$ $\displaystyle=$ $\displaystyle\ 0.$ (57)

Combining Eqs. (56, 57), we see that

$\gamma_{k}^{\prime}=\gamma_{k}^{\prime\prime}=\tilde{\gamma}_{k}=-\frac{\pi
k}{2},~{}~{}~{}k=\pm\frac{1}{2},$ (58)

precisely the Berry phase during an exchange process of the Majoranas.

Let us summarize the conditions required to have the quantized value of the
Berry phase:

  1. 1.

The Berry phase is robust to small deformations of the arc as long as the
flat-band condition is maintained.

  2. 2.

The magnetic field must go through zero in the straight-line segment.
($B_{y}=0$ in our case).

  3. 3.

The system must have the $\mathcal{T}$ symmetry in the absence of the Zeeman
field, and the $\tilde{\mathcal{T}}=U_{s,\pi}\mathcal{T}$ symmetry in its
presence. This means that the phase difference between the two superconducting
regions right outside the corner region must be $\pi$, which is naturally
realized in our setup with a $d$-wave superconductor.

####  IV.1.3 Phase gate via generic rotation of Zeeman field

Figure 4: The “slice-of-pie” countour in the $B_{x}$-$B_{y}$ plane, which
achieves a phase gate for the two Majorana modes at a given corner. The
magnetic field in the semicircle portion is in the tight-binding regime
$Ba\gg(a/\ell)^{2-K}$.

In this last subsection we build on the idea of the previous subsection and
show that it is possible to obtain a continuous family of Berry phase values
by taking the system along a “slice of pie” path in the parameter space of the
external magnetic field $\mathbf{B}=(B_{x},B_{y})$, shown in Fig. 4. The
specific path that we consider is as follows. We start with $B_{x}=B>0$ and
$B_{y}=0$. In the first part of the path we rotate the Zeeman field
counterclockwise by an angle of $2\tau_{0}$, ending up at
$B_{x}=B\cos(2\tau_{0})$ and $B_{y}=B\sin(2\tau_{0})$. In the second part of
the path we traverse the straight line segment from
$\mathbf{B}=(B\cos(2\tau_{0}),B\sin(2\tau_{0}))$ to the origin
$\mathbf{B}=(0,0)$. Finally, in the third part of the path we traverse the
straight segment from the origin $\mathbf{B}=(0,0)$ back to our starting point
$\mathbf{B}=(B,0)$. Crucially, we assume that the first part of the path,
namely the curved segment that is traversed at constant magnitude
$|\mathbf{B}|=B$, is taken in the tight-binding regime. We will also assume
that in the absence of $\mathbf{B}$, the theory has U(1)s symmetry.

To calculate the Berry phase $\gamma_{k}(\tau_{0})$ in this process, we first
consider the contribution from the two straight line paths,
$\mathcal{C}_{B}(0)$ and $\mathcal{C}_{B}(\tau_{0})$. Since the only term in
the system that violates the spin rotation symmetry $U_{s,\theta}$ in Eq. (11)
is the Zeeman field, and that the two paths $\mathcal{C}_{B}(0)$ and
$-\mathcal{C}_{B}(\tau)$ are related by the U(1)s symmetry, the total
contribution to the Berry phase along the straight paths
$\mathcal{C}_{B}(0)+\mathcal{C}_{B}(\tau)$ is zero.

Then $\gamma_{k}(\tau_{0})$ is exactly equal to the contribution from the
curved part of the path, and so

$\gamma_{k}(\tau_{0})=i\int_{0}^{\tau_{0}}d\tau\
\langle\psi_{k}(B,\tau)|\partial_{\tau}|\psi_{k}(B,\tau)\rangle\ .$ (59)

If we evaluate this using our variational approximation in the large $B$
regime (following the ideas from earlier in this section), then we find the
total Berry phase as

$\gamma_{k}(\tau_{0})=-k\tau_{0},~{}~{}~{}k=\pm\frac{1}{2}.$ (60)

We note that this argument can also be applied to the symmetry protection for
the exchange process. The result (60) is protected by the spin-rotation
symmetry U(1)s when there is no Zeeman field.

Again let us summarize the conditions required to have the quantized value of
the Berry phase:

  1. 1.

The Berry phase is robust to small deformations of the arc as long as the
flat-band condition is maintained.

  2. 2.

The two straight segments must be precisely related by a $2\tau_{0}$ rotation
in the $B_{x}-B_{y}$ plane.

  3. 3.

The system must have U(1)s symmetry when there is no Zeeman field, and the
$\tilde{\mathcal{T}}=U_{s,\pi}\mathcal{T}$ symmetry in its presence to protect
the ground state degeneracy.

Due to the ground state degeneracy it is also possible to choose as the
initial state a superposition of $k=\pm 1/2$. In the next subsection we
discuss such a situation where we choose a different basis for initial states.

###  IV.2 Clifford gates via manipulating Majorana modes within and across
corners

Figure 5: A qubit formed by four MZMs from adjacent corners. By tuning the
Zeeman field $\mathbf{B}$ and the superconducting order $\Delta$ one can
achieve symmetry protected Clifford gates.

In this Subsection we consider a configuration of two adjacent corners subject
to antiparallel Zeeman fields and the edge between them are gapped by SC
order, which we depict in Fig. 5. From Eq. (24), the tunneling current through
each corner is given by Eq. (42):

$J=\int_{0}^{\ell}\frac{dx}{\pi}\partial_{x}\vartheta(x)=\frac{k}{m}\mod\frac{1}{m}$
(61)

corresponding to the fermion parity at the corner

${(-)^{F}\equiv e^{i\pi k}.}$ (62)

Since quasiparticles can tunnel between corners through the edge, only the
combined fermion parity of the two corners is conserved. For a given parity
(say even), such a configuration with two corners and one edge form a single
qubit, with the two energy levels distinguished by the parity at a given
corner, which we label as

$\displaystyle|\uparrow\rangle\equiv$
$\displaystyle|k=\frac{1}{2},k^{\prime}=\frac{1}{2}\rangle,$
$\displaystyle|\downarrow\rangle\equiv$
$\displaystyle|k=-\frac{1}{2},k^{\prime}=-\frac{1}{2}\rangle,$ (63)

where $k$ and $k^{\prime}$ are the respective lattice momenta in the two
corners. With this notation, the exchange operation in either corner leads to
a Berry phase represented by

$\tilde{\gamma}=e^{i\frac{\pi\sigma_{z}}{4}},$ (64)

where $\sigma_{z}$ is the Pauli matrix in the Hilbert space of (63).

To realize Clifford gates, one additionally needs to achieve the non-Abelian
unitary operator

$\bar{\gamma}=e^{i\frac{\pi\sigma_{x}}{4}}.$ (65)

In Ref. Ivanov, 2001a this can be achieved by swapping different sets of
Majorana pairs. Similarly, here we show that $\bar{\gamma}$ is achieved by
manipulating two Majorana modes _across_ different corners.

As we discussed in Sec. II.1, with antiparallel in-plane Zeeman fields in the
two corners the edge region is described by a theory dual to the one for the
corner regions, with Majorana modes protected instead by $\mathcal{\tilde{C}}$
symmetry. According to (60), such a duality is simply the usual
$\vartheta\leftrightarrow\varphi$ duality in bosonization. In this basis, the
two states forming the qubit are then eigenstates of

$Q=\int_{\rm edge}\frac{dx}{\pi}\partial_{x}\varphi(x),$ (66)

which is the fermion parity in the superconducting edge. Within the subspace
of the ground state qubit, the charge eigenstate is a superposition between
the two different eigenstates for tunneling current $J$, thus we can rewrite
$Q$ as (which is time-reversal invariant)

$Q=|\uparrow\rangle\langle\downarrow|+|\downarrow\rangle\langle\uparrow|=\sigma_{x},$
(67)

and its eigenstates are labeled by $\langle\sigma_{x}\rangle=\pm 1$.

In order to induce Berry phases, we can similarly design contours in the
complex plane of $\Delta$ similar to that of $B=B_{x}+iB_{y}$. Assuming that
the system size is much larger than the superconducting coherence length, one
can tune the pairing fields for different edges independently. Due to the
$\vartheta\leftrightarrow\varphi$ duality, one can straightforwardly obtain
that a “half-moon” contour leads to the non-Abelian phase given in (65).
Notice here this value is topological and protected by the dual
$\mathcal{\tilde{C}}$ symmetry, in which, as we showed $\mathcal{C}$ is an
emergent symmetry guaranteed by properly gating the sample to charge
neutrality.

In addition, it is straightforward to see that the phase gate operation can be
realized for Majorana’s across different corners by taking a “slice of pie”
contour (analog of that in Fig. 4) of an angle $2\rho_{0}$ in $\Delta$, which
is protected by the $U(1)_{c}$ symmetry. Thus we have two types of phase gates
available, namely,

$\tilde{\gamma}(\tau_{0})=e^{i{\tau_{0}\sigma_{z}}},~{}~{}~{}\bar{\gamma}(\rho_{0})=e^{i{\rho_{0}\sigma_{x}}}.$
(68)

We note that in Ref. Wölms  _et al._ , 2016 the authors pointed out that a
qubit made out of a Kramers doublet of MZMs can be subject to a non-Abelian
Berry phase in the presence of an adiabatic local perturbation without lifting
the Kramers degeneracy, unless they carry distinct quantum numbers. In our
case at $B=0$, the Majorna zero modes overlap in space and form a Kramers
pair. However, the two states associated with the corner MZMs are
distinguished by their “lattice momenta” $k\in\\{-1/2,1/2\\}$ and their
fractional (topological) spins (see Eq. (42)) $S=\pm k/2$, and hence are
protected by symmetry from local perturbations.

Finally, note here that so far our platform has only involved two edge-sharing
corners of a semi-infinite sample. By simple math, a $d$-wave superconductor
setup produces four such corners, corresponding to three qubits (with a fixed
fermion parity for the sample). With the protocol above, one can realize a set
of Clifford gates on each edge, leading to a richer set of quantum gates in
the enlarged Hilbert space.

##  V Braiding parafermion modes

In Sec. IV we have completely focused on the $m=1$ case, in which the ground
states are MZMs. It is straightforward to generalize our full braid, exchange,
and phase gates to a generic $m$. For example, via a half-moon contour in
${B}$, we obtain

$\tilde{\gamma}_{k}=-\frac{\pi k}{2m},~{}~{}~{}k=[-m,m)\cap(\mathbb{Z}+1/2)\
.$ (69)

This result is exactly the exchange statistics of $\mathbb{Z}_{2m}$
parafermions. Lindner  _et al._ (2012); Cheng (2012); Vaezi (2013); Clarke
_et al._ (2013); Barkeshli  _et al._ (2013a, b)

However, unless $\beta\gg\alpha$, the $2m$ eigenstates in the lowest band are
not degenerate. In the half-moon contour, this indicates that on the straight
line portion through $B=0$ of the contour, even though the $2m$ states remain
separated from the other excited states, the topological Berry phase cannot be
separated from a $k$-dependent dynamical phase that is non-universal. The $2m$
states come in $m$ pairs, leading to $m-1$ independent relative dynamical
phases.

For small $m$, it may be possible to precisely control the system parameters
and the duration of the exchange process such that it is a common period for
all $m-1$ modes, which eliminates the dynamical phases. However, the result is
not topologically protected.

##  VI Conclusion

In this work we have proposed a platform for topological quantum computing
based on a heterostructure between a high-$T_{c}$ $d$-wave superconductor and
a quantum spin Hall insulator, which can be regarded as a higher-order
topological superconductor. We demonstrated that, via tuning the a Zeeman
field applied to the corner region and the superconducting order parameter,
such a setup can realize non-Abelian Clifford gates of Majorana qubits that
are protected by time-reversal and charge conjugation symmetries, as well as
phase gates protected by $U(1)$ symmetries. Within our analysis, interaction
effects and generalization to a fractional quantum spin Hall states can
naturally be incorporated.

In our proposed setup, the $d$-wave superconductor ensures a large critical
temperature and a large critical field, making the manipulation of Majorana’s
via an external Zeeman field easier to realize in experiments. Recent
advancements in low-dimensional materials have made the key components of the
heterostructure, including $d$-wave superconductors (and its monolayer version
Yu  _et al._ (2019)), quantum spin Hall insulators Wu  _et al._ (2018) and
two-dimensional ferromagnets Kezilebieke  _et al._ (2020), readily available.
Other than the specific combination of ingredients in our proposal, our
theoretical analysis is based on low-energy effective field theories, which
can be easily adapted to other topological materials with magnetism and
superconductivity. We note that recently signatures of parafermions have been
observed in a similar setup with a fractional quantum Hall insulator Gül  _et
al._ (2020). It will be extremely interesting to see if one can demonstrate
and manipulate these non-Abelian anyons using the protocols we propose in this
work.

###### Acknowledgements.

We thank A. Jahin, A. Tiwari and X.-X. Zhang for useful discussions. M.F.L.
acknowledges the support of the Kadanoff Center for Theoretical Physics at the
University of Chicago. M.F.L. is also supported by the Simons Collaboration on
Ultra-Quantum Matter, which is a grant from the Simons Foundation (651440).
M.C. acknowledges support from NSF under award number DMR-1846109 and the
Alfred P. Sloan foundation. Y.W. is supported by startup funds at the
University of Florida.

##  Appendix A Lattice model for second-order topological superconductor

In this Appendix we present a lattice model for the second-order topological
superconductor given by a quantum spin Hall (QSH) insulator with a proximity
effect induced $d$-wave superconudcting ($d$-SC) gap, given by $H_{0}+H_{\rm
SC}=\int
d\mathbf{k}\Psi^{\dagger}(\mathbf{k})\left(\mathcal{H}_{0}+\mathcal{H}_{\rm
SC}\right)(\mathbf{k})\Psi(\mathbf{k})$, where
$\Psi^{\dagger}=(\psi^{\dagger}(\mathbf{k}),\psi(-\mathbf{k}))$ and

$\displaystyle\mathcal{H}_{0}(\mathbf{k})+\mathcal{H}_{\rm SC}(\mathbf{k})=$
$\displaystyle\sin k_{a}s_{z}\sigma_{z}\tau_{0}+\sin
k_{b}s_{0}\sigma_{y}\tau_{z}$ $\displaystyle+(\cos k_{a}+\cos
k_{b}+m-1)s_{0}\sigma_{x}\tau_{z}$ $\displaystyle+\Delta\sin k_{x}\sin
k_{y}s_{y}\sigma_{0}\tau_{y}.$ (70)

Here $s_{x,y,z}$ denotes the spin degree of freedom, $\sigma_{x,y,z}$ is a
band index, and $\tau_{x,y,z}$ are Pauli matrices in the Nambu space. The
first three terms describes the normal state, which is a quantum spin Hall
insulator, and the last term is a pairing term of $d$-wave symmetry, coming
from the proximity effect with a high-$T_{c}$ superconductor. The Hamiltonian
has a time-reversal symmetry given by $\mathcal{T}=is_{y}K$, diagonal mirror
symmetries $\mathcal{M}_{a}=s_{x}\sigma_{z}\tau_{y}$,
$\mathcal{M}_{b}=s_{y}\sigma_{y}\tau_{x}$, and a particle-hole symmetry
${\mathcal{P}}=\tau_{x}K$.

Such a Hamiltonian has been studied in Refs. Yan  _et al._ , 2018b; Wang  _et
al._ , 2018c as a second-order topological superconductor protected by time-
reversal symmetry and a $C_{4}$ rotation symmetry. For our purposes, we will
instead rely on the mirror reflection symmetries $\mathcal{M}_{a,b}$. Higher-
order topological crystalline insulators and superconductors with mirror
symmetries have been classified in Ref. Trifunovic and Brouwer, 2019 based on
a $K$-theory analysis by Shiozaki and Sato Shiozaki and Sato (2014). In our
case, the reflection symmetry anticommutes with both time-reversal and
particle-hole conjugation. According to the terminology in Ref. Trifunovic and
Brouwer, 2019, it belongs to symmetry class DIII${}^{\mathcal{M}_{++}}$, which
in terms of second-order topology admits a $\mathbb{Z}_{2}$ classification in
2d. In the nontrivial phase symmetric corners of the sample host a pair of
MZMs that form Kramer doublet and have the _same_ mirror eigenvalue. Here the
$\mathbb{Z}_{2}$ classification is an _intrinsic_ bulk property. As such, one
cannot remove the Majorana doublet by modifying the boundary termination
without breaking the symmetry. For example, one can glue a 1d time-reversal
invariant topological superconductor on one of the edges, upon coupling to the
bulk, this gaps out the corner Majorana doublet, but this procedure
necessarily violates mirror symmetry.

We also consider an in-plane Zeeman field, either applied throughout the bulk
or only near the corners, given by $H_{\rm Z}=\int
d\mathbf{k}\Psi^{\dagger}(\mathbf{k})\mathcal{H}_{\rm
Z}(\mathbf{k})\Psi(\mathbf{k})$ where

$\displaystyle\mathcal{H}_{\rm Z}=$ $\displaystyle
B_{x}s_{x}\sigma_{0}\tau_{z}+B_{y}s_{y}\sigma_{0}\tau_{0}.$ (71)

The Zeeman field breaks both $\mathcal{T}$ and $\mathcal{M}_{a,b}$, but
preserves the composite symmetry $\mathcal{T}\mathcal{M}_{a,b}$. Importantly,
the other Zeeman term $\sim B_{z}s_{z}\sigma_{0}\tau_{z}$ is odd under this
action and is forbidden. Together with the particle-hole symmetry
${\mathcal{P}}=\tau_{x}K$, the Hamiltonian $H=H_{0}+H_{\rm Z}$ preserves
composite chiral (anti)symmetries
$\mathcal{PM}_{a,b}=s_{z}\sigma_{z}\tau_{z}$, which anticommutes with
${\mathcal{P}}$. Such a phase belongs to class D${}^{\mathcal{PM}_{-}}$, which
also admits a $\mathbb{Z}_{2}$ classification. This indicates the MZMs in the
absence of $H_{B}$ remains robust, despite time-reversal symmetry being
broken. However, as pointed out in Ref. Trifunovic and Brouwer, 2019, this
$\mathbb{Z}_{2}$ invariant is _extrinsic_. In fact, as we show in the main
text, the Majorana modes can be moved (or evem removed) by symmetric boundary
perturbations.

Experimentally, a $d$-SC/QSH heterostructure can be achieved by stacking
$d$-wave high-$T_{c}$ superconductor BSCCO and quantum spin Hall insulator
WTe2. Of course, depending material details and the geometry of stacking, such
a heterostucture may not realize the mirror symmetries we specified above.
However, as we show in the main text, the Majorana modes can be protected by
other emergent on-site symmetries.

##  Appendix B Mode expansion of the bosonic fields

In this appendix we explain in more detail the mode expansions for
$\varphi(x)$ and $\vartheta(x)$ that we use in our analysis in this paper. To
obtain the mode expansion for $\vartheta(x)$, we first identify a complete set
of functions of $x\in(0,\ell)$ that also obey the boundary conditions Eq.
(32), and then we expand $\vartheta(x)$ as a series in these functions with
operator-valued coefficients. We then expand $\varphi(x)$ in terms of a
complementary set of functions (also with operator-valued coefficients) in
such a way that $\vartheta(x)$ and $\varphi(x)$ obey the correct commutation
relations.

In our case the operator-valued coefficients that appear in the mode
expansions consist of zero mode operators $q$ and $p$ and a set of oscillator
raising and lowering operators $b_{n}$ and $b^{\dagger}_{n}$, with
$n\in\mathbb{N}\setminus\\{0\\}$ (i.e., we have an oscillator variable for
each integer $n\geq 1$). These operators obey the standard commutation
relations $[q,p]=i$ and
$[b_{n},b^{\dagger}_{n^{\prime}}]=\delta_{nn^{\prime}}$ (with all other
commutators vanishing). In addition, the zero mode operator $q$ is a compact
variable and is defined modulo $2\pi$, while its conjugate momentum $p$ is
defined to have integer eigenvalues. This means that the Hilbert space
$\mathcal{H}_{\text{zm}}$ associated with the zero mode operators $q$ and $p$
is spanned by the states $|s\rangle$, $s\in\mathbb{Z}$, which are eigenstates
of $p$, $p|s\rangle=s|s\rangle$, and with $e^{\pm iq}|s\rangle=|s\pm
1\rangle$. We can also define a basis $|q\rangle$ of eigenstates of $q$, with
$\langle q|s\rangle=\frac{e^{iqs}}{\sqrt{2\pi}}$, and we have $\langle
q|q^{\prime}\rangle=\delta_{2\pi}(q-q^{\prime})$, where
$\delta_{2\pi}(q-q^{\prime})=\sum_{s\in\mathbb{Z}}\frac{1}{2\pi}e^{i(q-q^{\prime})s}$
is the $2\pi$-periodic delta function.

In terms of these operators, the mode expansions for $\varphi(x)$ and
$\vartheta(x)$ take the form

$\displaystyle\varphi(x)$
$\displaystyle=mq-\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}\cos(\kappa_{n}x)(b_{n}+b^{\dagger}_{n})$ (72a)
$\displaystyle\vartheta(x)$ $\displaystyle=(p+\delta)\frac{\pi x}{m\ell}$
$\displaystyle+i\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}\sin(\kappa_{n}x)(b_{n}-b^{\dagger}_{n})\ ,$ (72b)

where $\kappa_{n}=\frac{\pi n}{\ell}$. Here we have set the twist for
$\vartheta$ at the two ends as a generic $\delta$; in the setup discussed in
the main text, we have $\delta=1/2$. The exponential factor $\epsilon$ is a
dimensionless ultraviolet cutoff, which we will discuss in details for $m=1$
(Appendix C) and $m\neq 1$ (Appendix F). We can see that this cutoff serves to
control the oscillator sums at high momenta $\kappa_{n}$. Removing the cutoff
corresponds to taking $a\to 0$, and one can check that in this limit the
fields obey the correct commutation relations
$[\vartheta(x),\partial_{x^{\prime}}\varphi(x^{\prime})]=[\varphi(x),\partial_{x^{\prime}}\vartheta(x^{\prime})]=\pi
i\delta(x-x^{\prime})$ (these follow from Eqs. (18) and the definition of
$\varphi(x)$ and $\vartheta(x)$ in terms of $\phi_{\uparrow/\downarrow}(x)$).

Finally, as in the main text, it is convenient to define the shifted zero mode
momentum operator $\tilde{p}$ via

$\tilde{p}=p+\delta\ .$ (73)

This will be useful because almost all of our expressions will involve the
shifted momentum $\tilde{p}$ instead of the original momentum $p$. Note that,
since $p$ is defined to have integer eigenvalues, the eigenvalues of
$\tilde{p}$ lie in the set $\mathbb{Z}+\delta$ (the integers shifted by
$\delta$).

##  Appendix C Bosonization in the $m=1$ case

In this appendix we explain how to carefully define the fermionic operators
$R(x)$ and $L(x)$ in terms of the bosonic fields $\phi_{\uparrow}(x)$ and
$\phi_{\downarrow}(x)$ in the non-fractional case with $m=1$. Specifically, we
define $R(x)$ and $L(x)$ as _normal-ordered exponentials_ of
$\phi_{\uparrow}(x)$ and $\phi_{\downarrow}(x)$, and with a dimensionful
prefactor that depends on the length $\ell$. We then show that the operators
$R(x)$ and $L(x)$ constructed in this way actually do obey the standard
anticommutation relations of fermionic fields.

###  C.1 Important identities

There are two basic identities that we will use repeatedly in the derivations
in this appendix, and so we record them here for reference. Let $X$ and $Y$ be
any two operators such that their commutator $[X,Y]$ is a c-number. Then we
have

$e^{X}e^{Y}=e^{Y}e^{X}e^{[X,Y]}$ (74)

and

$e^{X}e^{Y}=e^{X+Y}e^{\frac{1}{2}[X,Y]}\ .$ (75)

###  C.2 Definition of the fermion operators

We now present the definition of the operators $R(x)$ and $L(x)$. We start
with the mode expansions for the fields $\phi_{\uparrow}(x)$ and
$\phi_{\downarrow}(x)$, which take the form (recall that we take $m=1$)

$\displaystyle\phi_{\uparrow}(x)$ $\displaystyle=q+\tilde{p}\frac{\pi
x}{\ell}-\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}\left(e^{-i\kappa_{n}x}b_{n}+\text{h.c.}\right)$ (76a)
$\displaystyle\phi_{\downarrow}(x)$ $\displaystyle=q-\tilde{p}\frac{\pi
x}{\ell}-\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}\left(e^{i\kappa_{n}x}b_{n}+\text{h.c.}\right)\ .$ (76b)

In addition, recall that $\kappa_{n}=\frac{\pi n}{\ell}$ and that
$\epsilon=\frac{\pi a}{\ell}$. For later use, we note here that
$\phi_{\downarrow}(x)=\phi_{\uparrow}(-x)$ (this relation actually holds for
any $m$ and not just $m=1$).444It should be clear that there is no problem
with plugging a negative value of the position coordinate into our mode
expansions for $\phi_{\uparrow}(x)$ and $\phi_{\downarrow}(x)$.

Given these mode expansions, our definition of the fermion operators $R(x)$
and $L(x)$ is as follows. First, for any operator $\mathcal{O}$ of the form

$\mathcal{O}=Aq+\sum_{n=1}^{\infty}B_{n}b^{\dagger}_{n}+\sum_{n=1}^{\infty}C_{n}b_{n}+Dp\
,$ (77)

we define the normal-ordered exponential $:e^{\mathcal{O}}:$ by

$:e^{\mathcal{O}}:\ =\
e^{Aq}e^{\sum_{n=1}^{\infty}B_{n}b^{\dagger}_{n}}e^{\sum_{n=1}^{\infty}C_{n}b_{n}}e^{Dp}\
.$ (78)

Then our definition of $R(x)$ and $L(x)$ is

$\displaystyle R(x)$ $\displaystyle=$ $\displaystyle\frac{e^{i\delta\frac{\pi
x}{\ell}}}{\sqrt{2\ell}}:e^{-i\phi_{\uparrow}(x)}:$ (79a) $\displaystyle L(x)$
$\displaystyle=$ $\displaystyle\frac{e^{i\delta\frac{\pi
x}{\ell}}}{\sqrt{2\ell}}:e^{i\phi_{\downarrow}(x)}:\ .$ (79b)

In other words, the fermionic operators $R(x)$ and $L(x)$ are defined in terms
of normal-ordered exponentials of $\phi_{\uparrow}(x)$ and
$\phi_{\downarrow}(x)$, with an additional prefactor proportional to
$\ell^{-\frac{1}{2}}$. This prefactor ensures that the fermionic operators
have the correct units and anticommutation relations, as we show below. We can
also use the definition of the normal-ordered exponential to write out these
operators in more detail. For example, we find that

$\displaystyle R(x)=$
$\displaystyle\frac{1}{\sqrt{2\ell}}e^{-iq}\text{exp}\left\\{i\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}e^{i\kappa_{n}x}b^{\dagger}_{n}\right\\}$
$\displaystyle\times\text{exp}\left\\{i\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}e^{-i\kappa_{n}x}b_{n}\right\\}e^{-ip\frac{\pi x}{\ell}}\ .$
(80)

We now mention a few important properties of the operators $R(x)$ and $L(x)$.
First, these operators, as we have defined them above, are functions of the
ultraviolet cutoff $a$, although we have not indicated this dependence in our
notation. Later we will show that these operators behave exactly like
fermionic fields in the $a\to 0$ limit. We also note that, in our open
geometry (and with our choice of boundary conditions), the fields $R(x)$ and
$L(x)$ are not independent but are actually related by the identity

$L(x)=e^{i\frac{\pi x}{\ell}}R^{\dagger}(-x)\ .$ (81)

This identity can be derived by taking the Hermitian conjugate of our
expression for $R(-x)$, and by using the rearrangement identity

$e^{iq}e^{-ip\frac{\pi x}{\ell}}=e^{-ip\frac{\pi x}{\ell}}e^{iq}e^{i\frac{\pi
x}{\ell}}\ ,$ (82)

which can be derived using Eq. (74). For our setup, however, we will focus on
the interval $(0,\ell)$, in which $L$ and $R^{\dagger}$ can be treated as
independent fields.

###  C.3 Derivation of anticommutation relations

We now show that, in the limit $\epsilon\to 0$, the operators $R(x)$ and
$L(x)$ that we defined actually do obey the standard anticommutation relations
for fermionic fields. We start by deriving the anticommutator
$\\{R(x),R(y)\\}$ between the right-moving field at two different points $x$
and $y$. For this calculation we first define four quantities $(1)$, $(2)$,
$(3)$, and $(4)$ via

$\displaystyle(1)$ $\displaystyle=$ $\displaystyle e^{-ip\frac{\pi x}{\ell}}$
(83a) $\displaystyle(2)$ $\displaystyle=$ $\displaystyle e^{-iq}$ (83b)
$\displaystyle(3)$ $\displaystyle=$
$\displaystyle\text{exp}\left\\{i\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}e^{-i\kappa_{n}x}b_{n}\right\\}$ (83c) $\displaystyle(4)$
$\displaystyle=$
$\displaystyle\text{exp}\left\\{i\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}e^{i\kappa_{n}y}b^{\dagger}_{n}\right\\}\ .$ (83d)

Then using Eq. (74) we find that

$(1)(2)=(2)(1)e^{i\frac{\pi x}{\ell}}$ (84)

and

$\displaystyle(3)(4)$ $\displaystyle=$
$\displaystyle(4)(3)e^{-\sum_{n=1}^{\infty}\frac{e^{-\epsilon
n}}{n}e^{-i\kappa_{n}(x-y)}}$ (85) $\displaystyle=$
$\displaystyle(4)(3)e^{\ln\left[1-e^{-\epsilon}e^{-i\frac{\pi}{\ell}(x-y)}\right]}$
$\displaystyle=$
$\displaystyle(4)(3)\left[1-e^{-\epsilon}e^{-i\frac{\pi}{\ell}(x-y)}\right]\
,$

where we used the infinite series $\ln(1-z)=-\sum_{n=1}^{\infty}z^{n}/n$
(valid for $|z|<1$) to get from the first to the second line. Putting these
results together yields the formula

$R(x)R(y)=\frac{1}{2\ell}:e^{-i\phi_{\uparrow}(x)-i\phi_{\uparrow}(y)}:\left[e^{i\frac{\pi
x}{\ell}}-e^{-\epsilon}e^{i\frac{\pi y}{\ell}}\right]\ .$ (86)

By examining the term in square brackets, which tends to $e^{i\frac{\pi
x}{\ell}}-e^{i\frac{\pi y}{\ell}}$ as $a\to 0$, we can see that

$\lim_{a\to 0}\\{R(x),R(y)\\}=0\ \ \forall\ \ x,\ y\ ,$ (87)

which is the expected anticommutator for a fermionic field with itself.

Next, we consider the anticommutator of $R(x)$ with $R^{\dagger}(y)$. For this
calculation we define the operator $A(x)$ by

$A(x)=\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}e^{i\kappa_{n}x}b^{\dagger}_{n}\ ,$ (88)

and so

$A^{\dagger}(x)=\sum_{n=1}^{\infty}\frac{e^{-\frac{\epsilon
n}{2}}}{\sqrt{n}}e^{-i\kappa_{n}x}b_{n}\ .$ (89)

In terms of this operator we can rewrite $R(x)$ and $R^{\dagger}(y)$ as

$\displaystyle R(x)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\ell}}e^{-iq}e^{iA(x)}e^{iA^{\dagger}(x)}e^{-ip\frac{\pi
x}{\ell}}$ (90) $\displaystyle R^{\dagger}(y)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\ell}}e^{ip\frac{\pi
y}{\ell}}e^{-iA(y)}e^{-iA^{\dagger}(y)}e^{-iq}\ .$ (91)

Then, using similar rearrangement identities as in our previous calculation
(using Eq. (74) again), we obtain the formulas

$\displaystyle R(x)R^{\dagger}(y)$ $\displaystyle=$
$\displaystyle\frac{1}{2\ell}e^{-ip\frac{\pi}{\ell}(x-y)}e^{iA(x)-iA(y)}e^{iA^{\dagger}(x)-iA^{\dagger}(y)}$
(92) $\displaystyle\times$
$\displaystyle\frac{e^{-i\frac{\pi}{\ell}(x-y)}}{1-e^{-\epsilon}e^{-i\frac{\pi}{\ell}(x-y)}}$

and

$\displaystyle R^{\dagger}(y)R(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2\ell}e^{-ip\frac{\pi}{\ell}(x-y)}e^{iA(x)-iA(y)}e^{iA^{\dagger}(x)-iA^{\dagger}(y)}$
(93) $\displaystyle\times$
$\displaystyle\frac{1}{1-e^{-\epsilon}e^{i\frac{\pi}{\ell}(x-y)}}\ .$

For the anticommutator we then find the formula

$\displaystyle\\{R(x),R^{\dagger}(y)\\}$ $\displaystyle=$ $\displaystyle
e^{-ip\frac{\pi}{\ell}(x-y)}e^{iA(x)-iA(y)}e^{iA^{\dagger}(x)-iA^{\dagger}(y)}$
(94) $\displaystyle\times$ $\displaystyle d(x-y;\epsilon)\ ,$

where we defined the function $d(x-y;a)$ by

$d(x-y;\epsilon)=\frac{1}{2\ell}\left[\frac{e^{-i\frac{\pi}{\ell}(x-y)}}{1-e^{-\epsilon}e^{-i\frac{\pi}{\ell}(x-y)}}+\frac{1}{1-e^{-\epsilon}e^{i\frac{\pi}{\ell}(x-y)}}\right]\
.$ (95)

We now examine the properties of the function $d(x-y;a)$. For $a>0$ we can
expand the denominators in $d(x-y;a)$ as geometric series to obtain

$\displaystyle
d(x-y;\epsilon)=\frac{1}{2\ell}\Bigg{[}e^{-i\frac{\pi}{\ell}(x-y)}$
$\displaystyle\sum_{n=0}^{\infty}e^{-\epsilon n}e^{-i\frac{\pi n}{\ell}(x-y)}$
$\displaystyle+$ $\displaystyle\ \sum_{n=0}^{\infty}e^{-\epsilon
n}e^{i\frac{\pi n}{\ell}(x-y)}\Bigg{]}\ .$ (96)

From this expression we can see that

$\lim_{\epsilon\to 0}d(x-y;\epsilon)=\delta_{2\ell}(x-y)\ ,$ (97)

where

$\delta_{2\ell}(x-y)=\frac{1}{2\ell}\sum_{n\in\mathbb{Z}}e^{i\frac{2\pi
n}{2\ell}(x-y)}$ (98)

is the $2\ell$-periodic delta function. Since this function is zero for $x\neq
y$ modulo $2\ell$, and since all of the prefactors from Eq. (94) are equal to
1 when $x=y$ modulo $2\ell$, we then find that

$\lim_{a\to 0}\\{R(x),R^{\dagger}(y)\\}=\delta_{2\ell}(x-y)\ .$ (99)

Therefore we have proven that, in the limit $\epsilon\to 0$, the operator
$R(x)$ obeys the standard anticommutation relations for a fermionic field
operator. Since the anti-commutation relation (99) holds down to the smallest
length scales, one can identify the ultraviolet cutoff $\epsilon$ in the
bosonic theory as that in the fermionic theory, which is given by the lattice
constant $a$ as

$\epsilon=\frac{\pi a}{\ell}\ .$ (100)

For the left-moving field $L(x)$, we can use our results for $R(x)$ and the
relation (81) between $L(x)$ and $R(x)$ to immediately conclude that $L(x)$
also obeys the standard anticommutation relations for a fermionic field
operator.

Finally, there is one more interesting anticommutation relation that we can
obtain for our system with open boundary conditions. Since in this case the
left- and right-moving fermionic fields are not independent, we find that, for
$x,y\in(0,\ell)$,

$\displaystyle\\{R(x),L(y)\\}$ $\displaystyle=$ $\displaystyle e^{i\frac{\pi
y}{\ell}}\\{R(x),R^{\dagger}(-y)\\}$ (101) $\displaystyle=$ $\displaystyle
e^{i\frac{\pi y}{\ell}}\delta_{2\ell}(x+y)$ $\displaystyle=$ $\displaystyle 0\
,$

where the last line holds since $x+y\neq 0$ for $x,y\in(0,\ell)$. For
$x,y\in(0,\ell)$ we also have

$\displaystyle\\{R(x),L^{\dagger}(y)\\}$ $\displaystyle=$ $\displaystyle
e^{-i\frac{\pi y}{\ell}}\\{R(x),R(-y)\\}$ (102) $\displaystyle=$
$\displaystyle 0\ .$

Therefore, for our system with open boundary conditions, we find that the
left- and right-moving fermionic fields already have the correct
anticommutation relations, and we do not need to include any extra _Klein
factors_ to ensure that $R(x)$ and $L(y)$ (and $R(x)$ and $L^{\dagger}(y)$)
anticommute.

###  C.4 Correlation functions in the free theory

To complete this section we present the formulas for the two-point correlation
functions of $R(x)$ and $L(x)$ in the free theory before adding any
perturbation terms. The form of these correlation functions will complete the
demonstration that we have correctly constructed the fermionic operators from
the bosonic fields.

We consider the free bosonic vacuum $|0\rangle$ that satisfies $p|0\rangle=0$
and $b_{n}|0\rangle=0$ for all $n\in\\{1,2,3,\dots\\}$. Using our previous
rearrangement of the product $R^{\dagger}(y)R(x)$, we find that in this state
we have

$\langle
0|R^{\dagger}(y)R(x)|0\rangle=\frac{1}{2\ell}\frac{1}{\left[1-e^{-\epsilon}e^{i\frac{\pi}{\ell}(x-y)}\right]}\
,$ (103)

and by taking the complex conjugate we find that

$\langle
0|R^{\dagger}(x)R(y)|0\rangle=\frac{1}{2\ell}\frac{1}{\left[1-e^{-\epsilon}e^{-i\frac{\pi}{\ell}(x-y)}\right]}\
.$ (104)

To check that this formula makes sense, we can investigate its behavior in the
bulk of the system, which corresponds to taking the limit $\ell\to\infty$
while keeping $x-y$, $\delta$. We also hold the ultraviolet cutoff $a$ fixed
in this limit (although it is safe to take it to zero at this point). Then in
this limit we find that

$\langle
0|R^{\dagger}(x)R(y)|0\rangle\to\frac{1}{2\pi}\frac{1}{\left[i(x-y)+a\right]}\
,$ (105)

which is the correct bulk correlation function of a free right-moving fermion
in one spatial dimension (with the correct normalization).

We can now do a similar calculation for the left-moving fermion. Using the
relation between $R(x)$ and $L(x)$, we first find that

$\langle 0|L^{\dagger}(x)L(y)|0\rangle=e^{-i\frac{\pi}{\ell}(x-y)}\langle
0|R(-x)R^{\dagger}(-y)|0\rangle\ .$ (106)

Then, using our previous rearrangement of $\langle
0|R(x)R^{\dagger}(y)|0\rangle$, we find that

$\langle
0|L^{\dagger}(x)L(y)|0\rangle=\frac{1}{2\ell}\frac{1}{\left[1-e^{-\epsilon}e^{i\frac{\pi}{\ell}(x-y)}\right]}\
.$ (107)

If we now take the bulk limit then in this case we find that

$\langle
0|L^{\dagger}(x)L(y)|0\rangle\to\frac{1}{2\pi}\frac{1}{\left[-i(x-y)+a\right]}\
,$ (108)

which is the correct bulk correlation function of a free left-moving fermion
in one spatial dimension.

##  Appendix D The variational approximation

In this appendix we explain the variational approximation that we use to study
the ground states of the domain wall Hamiltonian from Eq. (33) of Sec. III.2
of the main text. This variational approximation leads us to an effective
Hamiltonian that only involves the zero mode operators $q$ and $p$ from the
mode expansions of $\varphi(x)$ and $\vartheta(x)$, and in the later
appendices we analyze this effective Hamiltonian in detail and use it to make
predictions for the physical properties of the original domain wall model.

The variational method is familiar from quantum mechanics. It allows one to
obtain information about the ground state of a system by making sufficiently
clever guesses for the form of the ground state wave function. Here we apply
this method to study the ground state of a quantum field theory, and in this
setting there are additional complications associated with divergences present
in a quantum field theory without a proper cutoff. Therefore, we perform our
variational calculation for the domain wall model with a finite ultraviolet
cutoff $a$. Then, at the end of the calculation, we consider the system with a
small but finite value of $a$, as in condensed matter systems it is sensible
to keep a finite ultraviolet cutoff $a$ (which can be intuitively thought of
as being related to the scale of the crystal lattice).

Our starting point is the full Hamiltonian for the domain wall model. We
denote this Hamiltonian by $H(a)$ to indicate that we are working with a
finite ultraviolet cutoff $a>0$. As in Sec. II, this Hamiltonian takes the
form

$H(a)=H_{0}(a)-B\int_{0}^{\ell}dx\
\left(R^{\dagger}(x)L(x)+\text{h.c.}\right)\ .$ (109)

Here, $R(x)$ and $L(x)$ are the fermionic operators from Eq. (79) that we
constructed from the bosonic fields $\phi_{\uparrow}(x)$ and
$\phi_{\downarrow}(x)$. For simplicity we also assume that the magnetic field
points along the positive $x$-axis, so that $\tau=0$ in our previous notation.
Finally, $H_{0}(a)$ is the part of the Hamiltonian that contains the kinetic
energy term and the density-density interactions.

To prepare for our variational approximation, we first write out the term
$H_{0}(a)$ using our mode expansions for $\varphi(x)$ and $\vartheta(x)$. We
find that

$\displaystyle H_{0}(a)$ $\displaystyle=\frac{vK\pi}{2\ell}\tilde{p}^{2}$
$\displaystyle+\frac{v}{4}\left(\frac{1}{K}+K\right)\sum_{n=1}^{\infty}e^{-\epsilon
n}\kappa_{n}(b^{\dagger}_{n}b_{n}+b_{n}b^{\dagger}_{n})$
$\displaystyle+\frac{v}{4}\left(\frac{1}{K}-K\right)\sum_{n=1}^{\infty}e^{-\epsilon
n}\kappa_{n}(b_{n}b_{n}+b^{\dagger}_{n}b^{\dagger}_{n})\ .$ (110)

The oscillator part of $H_{0}(a)$ can be diagonalized by making a Bogoliubov
transformation to new oscillator variables $a_{n}$ defined by

$a_{n}=\cosh(\eta)b_{n}+\sinh(\eta)b^{\dagger}_{n}\ ,$ (111)

where the real parameter $\eta$ is related to $K$ as

$e^{-2\eta}=K.$ (112)

In terms of these new variables, we find that

$\displaystyle H_{0}(a)$ $\displaystyle=\frac{vK\pi}{2\ell}\tilde{p}^{2}$
$\displaystyle+\frac{v}{2}\sum_{n=1}^{\infty}e^{-\epsilon
n}\kappa_{n}(a^{\dagger}_{n}a_{n}+a_{n}a^{\dagger}_{n})\ .$ (113)

For later use we also note the reverse Bogoliubov transformation,

$b_{n}=\cosh(\eta)a_{n}-\sinh(\eta)a^{\dagger}_{n}\ ,$ (114)

which allows us to express $b_{n}$ in terms of $a_{n}$ and $a^{\dagger}_{n}$.

As we mentioned in the main text, the Zeeman term $\sim B\int R^{\dagger}L$
gaps out the region and makes the oscillator modes massive. To find a suitable
trial state for the oscillator part of the Hilbert space, an additional
Bogoliubov transformation is needed. To this end we introduce yet another set
of oscillator variables, which we denote by $\tilde{a}_{n}$ (with Hermitian
conjugates $\tilde{a}_{n}^{\dagger}$). These will be related to the $a_{n}$
oscillators via the Bogoliubov transformation

$a_{n}=\cosh(\zeta_{n})\tilde{a}_{n}-\sinh(\zeta_{n})\tilde{a}^{\dagger}_{n}\
,$ (115)

where we have allowed the parameter $\zeta_{n}\in\mathbb{R}$ that determines
the transformation to depend on the index $n$. Using these new variables, we
can rewrite $H_{0}(a)$ in the form

$\displaystyle H_{0}(a)$ $\displaystyle=\frac{vK\pi}{2\ell}\tilde{p}^{2}$
$\displaystyle+\frac{v}{2}\sum_{n=1}^{\infty}e^{-\epsilon
n}\kappa_{n}\cosh(2\zeta_{n})(\tilde{a}^{\dagger}_{n}\tilde{a}_{n}+\tilde{a}_{n}\tilde{a}^{\dagger}_{n})$
$\displaystyle-\frac{v}{2}\sum_{n=1}^{\infty}e^{-\epsilon
n}\kappa_{n}\sinh(2\zeta_{n})(\tilde{a}_{n}\tilde{a}_{n}+\tilde{a}^{\dagger}_{n}\tilde{a}^{\dagger}_{n})\
.$ (116)

We also find that $b_{n}$ is related to $\tilde{a}_{n}$ via the relation

$b_{n}=\cosh(\eta+\zeta_{n})\tilde{a}_{n}-\sinh(\eta+\zeta_{n})\tilde{a}^{\dagger}_{n}\
,$ (117)

which follows from identities for the hyperbolic trigonometric functions.

We now discuss our choice of variational trial state. Let $|0,\zeta\rangle$ be
the Fock vacuum state annihilated by all the $\tilde{a}_{n}$,

$\tilde{a}_{n}|0,\zeta\rangle=0\ \forall\ n\ .$ (118)

Note also that, as we are now working in terms of the $\tilde{a}_{n}$
oscillator variables, the full Hilbert space $\mathcal{H}_{\text{tot}}$ of our
domain wall model is equal to the tensor product

$\mathcal{H}_{\text{tot}}=\mathcal{H}_{\text{zm}}\otimes\mathcal{H}_{\text{F}}\
,$ (119)

where $\mathcal{H}_{\text{F}}$ is the Fock space generated by the action of
the raising operators $\tilde{a}_{n}^{\dagger}$ on the Fock vacuum
$|0,\zeta\rangle$, and $\mathcal{H}_{\text{zm}}$ is the Hilbert space for the
zero modes ($q$ and $p$ act within $\mathcal{H}_{\text{zm}}$). The trial
ground state $|\Psi\rangle$ that we consider respects the tensor product
structure of the Hilbert space and it takes the tensor product form

$|\Psi\rangle=|\psi\rangle\otimes|0,\zeta\rangle\ ,$ (120)

where $|\psi\rangle$ is a state in the zero mode Hilbert space
$\mathcal{H}_{\text{zm}}$ and $|0,\zeta\rangle$ is the Fock vacuum for the
$\tilde{a}_{n}$ variables.

The nontrivial part of our variational calculation is the problem of finding
the parameters $\zeta_{n}$ and the zero mode state
$|\psi\rangle\in\mathcal{H}_{\text{zm}}$ that minimize the energy expectation
value $\langle\Psi|H(a)|\Psi\rangle$. In fact, we will not carry out this
optimization procedure completely on the $\zeta_{n}$ parameters. Instead, we
will use a heuristic argument to obtain the behavior of the energy expectation
value with the correct choice of $\zeta_{n}$.

To proceed with the variational calculation we need to compute the expectation
value $\langle\Psi|H(a)|\Psi\rangle$ and then consider this expectation value
in the small $a$ limit. We now present this calculation, omitting many of the
details since the required manipulations are similar to the ones we used in
Appendix C to prove the bosonization formulas. For the kinetic term we find
that

$\langle\Psi|H_{0}(a)|\Psi\rangle=\frac{vK\pi}{2\ell}\langle\psi|\tilde{p}^{2}|\psi\rangle+\frac{v}{2}\sum_{n=1}^{\infty}e^{-\epsilon
n}\kappa_{n}\cosh(2\zeta_{n})\ ,$ (121)

where the second term here is the vacuum energy for the $\tilde{a}_{n}$
oscillators. For the Zeeman term we find that

$\langle\Psi|R^{\dagger}(x)L(x)|\Psi\rangle=\frac{1}{2}f(x;a;\zeta)\langle\psi|e^{i2q}|\psi\rangle\
,$ (122)

where the function $f(x;a;\zeta)$ is given by

$f(x;a;\zeta)=\frac{1}{\ell}\frac{1}{1-e^{-\epsilon}}\left[\frac{1-e^{-\epsilon-i\frac{2\pi
x}{\ell}}}{1-e^{-\epsilon+i\frac{2\pi
x}{\ell}}}\right]^{\frac{1}{2}}e^{i\frac{2\pi
x}{\ell}}e^{-\sum_{n=1}^{\infty}\frac{e^{-\epsilon
n}}{n}e^{-2(\eta+\zeta_{n})}(1+\cos(2\kappa_{n}x))}\ .$ (123)

We note that the summation in the exponent of the last factor

$\sum_{n=1}^{\infty}\frac{e^{-\epsilon n}}{n}e^{-2(\eta+\zeta_{n})}$ (124)

is logarithmical, with the series effectively truncated by the factors
$e^{-\epsilon n}$ and $e^{-2\zeta_{n}}$. The first factor $e^{-\epsilon n}$,
with $\epsilon=\pi a/\ell$, is a ultraviolet cutoff for the mode number,
$n\lesssim\ell/a$. The second factor $e^{-2\zeta_{n}}$, on the other hand, is
due to the quasiparticle mass gap $\Delta_{Z}$ induced by the Zeeman field
$B$. Obviously we have

$\Delta_{Z}=B\textrm{ for free fermions,}$ (125)

but in general $\Delta_{Z}$ is renormalized by interactions and is a parameter
determined by the choice of $\\{\zeta_{n}\\}$. The $e^{-2\zeta_{n}}$ factor
effectively serves as an infrared cutoff for the mode number $n$, as for high
enough modes $\kappa_{n}\gg\Delta_{Z}$ the effect of the Zeeman field can be
neglected. This means that the correct choice of $\zeta_{n}$ sets a lower
limit of summation, $n\gtrsim\Delta_{Z}\ell$. Therefore in the regime

$a\ll 1/\Delta_{Z}\ll\ell,$ (126)

the following summation can be approximated by

$\displaystyle\sum_{n=1}^{\infty}\frac{e^{-\epsilon
n}}{n}e^{-2(\eta+\zeta_{n})}\approx-K\ln(\Delta_{Z}a),$ (127)

where we have used Eq. (112). Plugging (127) into (123) we obtain

$f(x;a;\zeta)=\Delta_{Z}^{K}a^{K-1}f_{0}(x;a;\zeta),$ (128)

where $f_{0}$ is a well-behaved $\mathcal{O}(1)$ function. With this choice,
we then define an energy scale $\beta$ via

$\beta=B\int_{0}^{\ell}dx\ f(x;a;\zeta)\ ,$ (129)

and we find that at small $a$ we have an extensive behavior

$\beta\sim B\Delta_{Z}^{K}a^{K-1}\ell.$ (130)

For a free fermion system with $K$=1, this energy is $\sim B^{2}\ell$. Using
these results, we can now complete our calculation of
$\langle\Psi|H(a)|\Psi\rangle$. We find that

$\langle\Psi|H(a)|\Psi\rangle=\alpha\langle\psi|\tilde{p}^{2}|\psi\rangle-\beta\langle\psi|\cos(2q)|\psi\rangle\
,$ (131)

where the coefficients $\alpha$ and $\beta$ are given by

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle\frac{vK\pi}{2\ell}$
(132a) $\displaystyle\beta$ $\displaystyle\sim$ $\displaystyle
B\Delta_{Z}^{K}a^{K-1}\ell\ ,$ (132b)

and we have omitted the $c$-number vacuum energy term for the $\tilde{a}_{n}$
oscillators. This result tells us that for our variational approximation the
zero mode state $|\psi\rangle$ should be chosen to be the lowest energy state
of the effective zero mode Hamiltonian

$H_{\text{eff}}=\alpha\tilde{p}^{2}-\beta\cos(2q).$ (133)

We notice that in the limit $\beta\gg\alpha$, this Hamiltonian describes an
approximate harmonic oscillator, with energy level spacing given by
$\sim\sqrt{\alpha\beta}$. These energy levels form a Fock space of zero modes,
now with mass $\sim\sqrt{\alpha\beta}$. We can identify the mass scale of
former zero modes with that of the oscillator modes:

$\Delta_{Z}\sim\sqrt{\alpha\beta}.$ (134)

Self-consistency between Eqs. (132) and (134) fixes the scale of $\Delta_{Z}$
as

$\Delta_{Z}\sim B\left(Ba\right)^{\frac{K-1}{2-K}}.$ (135)

Notably, we indeed recover $\Delta_{Z}=B$ for the free fermion case $K=1$.
With Eq. (135), the condition (126) translates to

$\frac{1}{a}\left(\frac{a}{\ell}\right)^{2-K}\ll B\ll\frac{1}{a},$ (136)

which requires $K<2$. For $K\geq 2$ a separate variational ansatz is needed,
which we postpone to future studies. As a sanity check, we see that the
condition $\beta\gg\alpha$ we needed is precisely one of the conditions in Eq.
(136). We note that the condition $K<2$ is consistent with Kosterlitz
renormalization group results on the sine-Gordon model in infininte spacetime,
under which the cosine term is a relevant perturbation.

We can rewrite the parameters $\alpha$ and $\beta$ as

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle\frac{vK\pi}{2\ell}$
(137a) $\displaystyle\beta$ $\displaystyle\sim$ $\displaystyle
B^{\frac{2}{2-K}}a^{\frac{2K-2}{2-K}}\ell.$ (137b)

As we discuss in the next appendix, $H_{\text{eff}}$ is closely related to
Mathieu’s equation, and so we can use known results on that equation to study
$H_{\text{eff}}$ and solve our variational problem. In that appendix we
present error estimates for various quantities, and those error estimates are
_exponentially_ small in $\ell$. The key to obtaining that scaling for the
error estimates is the fact that the parameters $\alpha$ and $\beta$ in
$H_{\text{eff}}$ satisfy the relation

$\lambda^{2}\equiv\frac{\beta}{\alpha}\sim
B^{\frac{2}{2-K}}a^{\frac{2K-2}{2-K}}\ell^{2}\ .$ (138)

It is convenient to define a correlation length such that
$\lambda\propto\ell/\xi$, and we have

$\xi/a\sim\left(\frac{1}{Ba}\right)^{\frac{1}{2-K}},$ (139)

which diverges at $K\to 2$ obeying the familiar Kosterlitz-Thouless scaling
behavior.

The full Hamiltonian $H(a)$ and the effective Hamiltonian $H_{\text{eff}}$
both have a $\mathbb{Z}_{2}$ symmetry generated by the operator $e^{i\pi p}$
(the symmetry is $\mathbb{Z}_{2}$ because $p$ has integer eigenvalues). This
means that the Hilbert space of the domain wall model is broken up into two
different sectors, where the states in each sector have opposite eigenvalues
($\pm 1$) of $e^{i\pi p}$. It also means that we should carry out our
variational calculation separately in each sector to study the ground state of
the Hamiltonian within each sector.

We can gain a more physical understanding of this $\mathbb{Z}_{2}$ symmetry by
noting that $e^{i\pi p}$ is proportional to the parity $e^{i\pi S}$ of the
total spin $S$ in the FM region, since $S$ is given explicitly by

$\displaystyle S$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{\ell}dx\
\left[\partial_{x}\phi_{\uparrow}(x)-\partial_{x}\phi_{\downarrow}(x)\right]$
(140) $\displaystyle=$ $\displaystyle p+\delta$ $\displaystyle=$
$\displaystyle\tilde{p}\ .$

From a physical point of view this makes sense since the spin parity $e^{i\pi
S}$ commutes with the Hamiltonian in the FM region. Therefore in what follows
it is convenient for us to label the two sectors of the Hilbert space by a
“lattice momentum” $k\in[-1,1)\cap(\mathbb{Z}+\delta)$, such that any given
state is an eigenstate of $e^{i\pi\tilde{p}}$ with eigenvalue $e^{i\pi
k}$.555Note that $k$ should lie in $[-1,1)$ because this is the first
Brillouin zone for a lattice with lattice spacing equal to $\pi$. Note that,
since $p$ has integer eigenvalues, there are only two such values of $k$ in
the set $[-1,1)\cap(\mathbb{Z}+\delta)$. For example, in the case of
$\delta=\frac{1}{2}$, which is our main interest, we have
$k\in\\{-\tfrac{1}{2},\tfrac{1}{2}\\}$.

Our variational method can be used to study the lowest energy states of $H(a)$
with all possible eigenvalues of $e^{i\pi\tilde{p}}$ (i.e., all possible
values of the lattice momentum $k$). In particular, we are interested in
estimating the energy splitting between the lowest energy states of $H(a)$
with different $e^{i\pi\tilde{p}}$ eigenvalues. We therefore define separate
variational trial states
$|\Psi_{k}\rangle=|\psi_{k}\rangle\otimes|0,\zeta\rangle$ for each allowed
value of $k$, where $|\psi_{k}\rangle$ should be chosen to be the ground state
of $H_{\text{eff}}$ in the sector of $\mathcal{H}_{\text{zm}}$ with
$e^{i\pi\tilde{p}}=e^{i\pi k}$.

Let $E_{k}$ be the energy of the ground state $|\psi_{k}\rangle$ of
$H_{\text{eff}}$ in the sector with $e^{i\pi\tilde{p}}=e^{i\pi k}$,

$H_{\text{eff}}|\psi_{k}\rangle=E_{k}|\psi_{k}\rangle\ \ \ ,\ \ \
e^{i\pi\tilde{p}}|\psi_{k}\rangle=e^{i\pi k}|\psi_{k}\rangle\ ,$ (141)

and let $k_{1}$ and $k_{2}$ be the two allowed values of $k$ in the set
$[-1,1)\cap(\mathbb{Z}+\delta)$. Our variational estimate for the energy
splitting $\Delta E$ between the ground states of the domain wall model in the
sectors with $e^{i\pi\tilde{p}}=e^{i\pi k_{1}}$ and $e^{i\pi\tilde{p}}=e^{i\pi
k_{2}}$ is then given by

$\displaystyle\Delta E$
$\displaystyle\approx\Big{|}\Big{(}\langle\Psi_{k_{1}}|H(a)|\Psi_{k_{1}}\rangle-\langle\Psi_{k_{2}}|H(a)|\Psi_{k_{2}}\rangle\Big{)}\Big{|}$
$\displaystyle=|E_{k_{1}}-E_{k_{2}}|\ .$ (142)

In the next appendix we use known results on Mathieu’s differential equation
to show that $|E_{k_{1}}-E_{k_{2}}|$ is exponentially small in $\ell$. Then
our variational approximation predicts that the ground states of the domain
wall model with different $e^{i\pi\tilde{p}}$ eigenvalues are very nearly
degenerate for large $\ell$. In addition, for the special case where
$\delta=\frac{1}{2}$, which is our main interest in this paper, we show in
Appendices G and H that the full domain wall Hamiltonian $H(a)$ has an
additional discrete symmetry that guarantees that every eigenstate of $H(a)$
has a partner with the exact same energy, and so the ground state of of $H(a)$
is exactly degenerate for any $\ell$ (not just approximately degenerate for
large $\ell$).

##  Appendix E Results from the variational approximation (Mathieu’s equation)

In this appendix we present our main results on the domain wall model that we
described in Sec. III.2. We first apply known mathematical results on
Mathieu’s equation to understand the ground states of the effective zero mode
Hamiltonian $H_{\text{eff}}$. We then use these results in our variational
approximation to obtain nontrivial predictions for certain properties of the
domain wall model, including the finite-size splitting of the nearly
degenerate ground states (when $\delta\neq 1/2$ – there is an exact degeneracy
at $\delta=\frac{1}{2}$ that we discuss in Appendices G and H), the
correlation functions of the fermionic operators, and the Berry phase for
certain adiabatic processes involving the external magnetic field.

###  E.1 Bound on $|E_{k_{1}}-E_{k_{2}}|$ for $H_{\text{eff}}$ and estimate of
the splitting $\Delta E$

We first explain how known results on Mathieu’s differential equation can be
used to bound the difference $|E_{k_{1}}-E_{k_{2}}|$ between the energies of
the lowest energy states of $H_{\text{eff}}$ in the two sectors with different
$e^{i\pi\tilde{p}}$ eigenvalue. We start by explaining the relation between
$H_{\text{eff}}$ and Mathieu’s differential equation. We first note that, by
construction, $p$ has integer eigenvalues. It follows from this that all
states in the zero mode Hilbert space $\mathcal{H}_{\text{zm}}$ are invariant
under the action of $e^{i2\pi p}$, which is the operator that translates $q$
by $2\pi$, $e^{i2\pi p}qe^{-i2\pi p}=q+2\pi$. Therefore, if $|\psi\rangle$ is
any state in $\mathcal{H}_{\text{zm}}$, then its wave function
$\psi(q)=\langle q|\psi\rangle$ is $2\pi$-periodic in $q$,
$\psi(q+2\pi)=\psi(q)$. Next, as we discussed in the previous appendix,
$H_{\text{eff}}$ also commutes with the operator $e^{i\pi p}$ that translates
$q$ by $\pi$, $e^{i\pi p}qe^{-i\pi p}=q+\pi$. Accordingly, all eigenstates of
$H_{\text{eff}}$ can be chosen to be eigenstates of $e^{i\pi p}$, as we have
discussed (and we actually labeled states by their eigenvalue of the closely
related operator $e^{i\pi\tilde{p}}$).

Let $|\psi\rangle$ be an eigenstate of $H_{\text{eff}}$ with energy $E$. Then
the wave function $\psi(q)$ satisfies the Schrodinger equation

$\alpha\left(-i\frac{d}{dq}+\delta\right)^{2}\psi(q)-\beta\cos(2q)\psi(q)=E\psi(q)\
,$ (143)

where $p$ has become the differential operator $-i\frac{d}{dq}$. If we define
a new wave function $\chi(q)$ by

$\psi(q)=e^{-i\delta q}\chi(q)\ ,$ (144)

then we find that $\chi(q)$ satisfies

$-\alpha\chi^{\prime\prime}(q)-\beta\cos(2q)\chi(q)=E\chi(q)\ ,$ (145)

where $\chi^{\prime}(q)=\frac{d\chi(q)}{dq}$. This equation can be brought
into the standard form of Mathieu’s equation by dividing through by
$\alpha\neq 0$ to obtain

$-\chi^{\prime\prime}(q)-\lambda^{2}\cos(2q)\chi(q)=\mathcal{E}\chi(q)\ ,$
(146)

where we remind $\lambda^{2}=\beta/\alpha$ and $\mathcal{E}=E/\alpha$. In
addition, the $2\pi$-periodicity of $\psi(q)$ implies that $\chi(q)$ obeys the
periodicity condition

$\chi(q+2\pi)=e^{i2\pi\delta}\chi(q)\ .$ (147)

To apply known results from the study of the Mathieu’s equation, we need to
study the behavior of $\chi(q)$ under translations by $\pi$, which is the
period of the potential $\cos(2q)$ that appears in the equation. This behavior
will depend on the eigenvalue of a given state under the action of the
operator $e^{i\pi p}$. In particular, for a state $|\psi_{k}\rangle$ that
satisfies $e^{i\pi p}|\psi_{k}\rangle=e^{i\pi(k-\delta)}|\psi_{k}\rangle$ (so
that $e^{i\pi\tilde{p}}|\psi_{k}\rangle=e^{i\pi k}|\psi_{k}\rangle$), we find
that the corresponding function $\chi_{k}(q)=e^{i\delta
q}\psi_{k}(q)=e^{i\delta q}\langle q|\psi_{k}\rangle$ satisfies the
periodicity condition

$\displaystyle\chi_{k}(q+\pi)$ $\displaystyle=$ $\displaystyle
e^{i\delta(q+\pi)}\psi_{k}(q+\pi)$ (148) $\displaystyle=$ $\displaystyle
e^{i\pi k}\chi_{k}(q)\ ,$

and this simple relation explains why we chose to label our states by their
eigenvalue of $e^{i\pi\tilde{p}}$ instead of their eigenvalue of $e^{i\pi p}$.

It is known from Floquet theory (similar to Bloch’s theorem from condensed
matter physics), that the spectrum of the Mathieu operator
$-\frac{d^{2}}{dq^{2}}-\lambda^{2}\cos(2q)$ is divided into distinct energy
bands. In addition, the eigenfunctions within each energy band are labeled by
a wave number $k\in[-1,1)$, which corresponds to the Brillouin zone of a one-
dimensional lattice with period $\pi$. An eigenfunction $\chi_{k}(q)$
characterized by the wave number $k$ obeys exactly the periodicity condition
from Eq. (148). From this we see that the lowest energy state of
$H_{\text{eff}}$ in the sector with $e^{i\pi\tilde{p}}=e^{i\pi k}$ corresponds
exactly to the eigenfunction labeled by $k$ within the lowest band of the
spectrum of $-\frac{d^{2}}{dq^{2}}-\lambda^{2}\cos(2q)$. Therefore, the energy
splitting $|E_{k_{1}}-E_{k_{2}}|$ between the two lowest energy states of
$H_{\text{eff}}$ with different $e^{i\pi\tilde{p}}$ eigenvalues is certainly
less than $\alpha$ times the width $|W_{0}(\lambda)|$ of the lowest band of
the Mathieu operator $-\frac{d^{2}}{dq^{2}}-\lambda^{2}\cos(2q)$ (we multiply
by $\alpha$ because $E=\alpha\mathcal{E}$).

An asymptotic formula for the width $|W_{0}(\lambda)|$ at large $\lambda$ was
obtained in Ref. Weinstein and Keller, 1985 (see also Ref. Dunne and Unsal,
2016 for a convenient summary of the properties of the spectrum of the Mathieu
operator). It takes the form666In Ref. Weinstein and Keller, 1985 the
bandwidth $|W_{0}(\lambda)|$ was denoted by $|B_{0}(\lambda)|$, but we use
$|W_{0}(\lambda)|$ here to avoid confusion with the magnetic field in our
problem.

$|W_{0}(\lambda)|=\frac{2^{\frac{19}{4}}}{\pi^{\frac{1}{2}}}\lambda^{\frac{3}{2}}e^{-\lambda\sqrt{8}}\left[1+O(\lambda^{-\frac{1}{2}})\right]\
.$ (149)

The key feature of this formula is the factor of $e^{-\lambda\sqrt{8}}$. The
presence of this factor implies that, when $\lambda$ is large, the width
$|W_{0}(\lambda)|$ of the lowest band is exponentially small in $\lambda$. Now
for our model (which we obtained from our variational approximation), this
means that the splitting $|E_{k_{1}}-E_{k_{2}}|$ is exponentially small in
$\ell$,

$|E_{k_{1}}-E_{k_{2}}|\lesssim\text{constant}\times e^{-\frac{\ell}{\xi}}\ ,$
(150)

where $\xi\equiv\lambda\sqrt{8}$ is the correlation length given by Eq.
(139)). Thus, our variational approximation predicts that for large $\ell$ the
energy splitting $\Delta E$ of the two ground states in our domain wall model
is exponentially small in the length $\ell$ of the FM region. For the free
fermion case with $K=1$, the correlation length is given by $\xi\sim 1/B$,
consistent with the decaying behavior from solving the Dirac equation with a
mass domain wall.

Finally, we close this section by noting that the case we are most interested
in in this paper is the special case where $\delta=\frac{1}{2}$. In Appendices
G and H we will show that in this case the two ground states of our domain
wall model are _exactly_ degenerate, and not just approximately degenerate as
we have predicted here for a general $\delta$.

###  E.2 Approximate form of $\chi_{k}(q)$ at large $\lambda$

For the Berry phase calculation later in this appendix we will need to
understand the form of the eigenfunctions $\chi_{k}(q)$ of the Mathieu
operator in the limit of large $\lambda$ (we referred to this as the “tight-
binding” limit in the main text). Therefore, in this subsection we review some
known facts about $\chi_{k}(q)$ in this limit.

When $\lambda$ is large, the eigenfunctions of the Mathieu operator in its
lowest band are well-approximated by a weighted sum of Gaussians localized in
each valley of the $\cos(2q)$ potential (see the proof of Theorem 1 in Ref.
Simon, 1984). These approximate eigenfunctions can be constructed as follows.
We first expand $\cos(2q)$ to order $q^{2}$ about its minimum at $q=0$ and
study the resulting approximate Mathieu operator near $q=0$. Up to a constant,
we find the operator $-\frac{d^{2}}{dq^{2}}+2\lambda^{2}q^{2}$, and it is
well-known that the lowest energy eigenfunction of this operator is a Gaussian
of the form

$\chi_{0}(q)=\left(\frac{\lambda\sqrt{2}}{\pi}\right)^{\frac{1}{4}}e^{-\frac{\lambda}{\sqrt{2}}q^{2}}\
,$ (151)

where we have chosen the coefficient so that $\int_{-\infty}^{\infty}dq\
|\chi_{0}(q)|^{2}=1$.

Let $\chi_{k}(q)$ be the eigenfunction in the lowest band of the Mathieu
operator and obeying the periodicity condition
$\chi_{k}(q+\pi)=e^{ik\pi}\chi_{k}(q)$. At large $\lambda$, this eigenfunction
is given approximately by the periodic sum

$\chi_{k}(q)=\frac{1}{\sqrt{2}}\sum_{n\in\mathbb{Z}}e^{ikn\pi}\chi_{0}(q-n\pi)\
,$ (152)

which contains all translations of $\chi_{0}(q)$ by integer multiples of
$\pi$, with the translation by $n\pi$ accompanied by the $k$-dependent phase
factor $e^{ikn\pi}$. The factor of $\sqrt{2}$ is included here so that
$\chi_{k}(q)$ obeys the normalization condition777In our problem the original
wave functions are defined for $q\in[0,2\pi)$. This explains our extra factor
of $\sqrt{2}$ as compared with Ref. Simon, 1984, where the wave functions were
normalized for integration over one period of the periodic potential (which is
$\pi$ in our case).

$\int_{0}^{2\pi}dq\
\overline{\chi_{k}(q)}\chi_{k}(q)=1+O(e^{-\frac{\pi^{2}}{2\sqrt{2}}\lambda})\
,$ (153)

where the integral is restricted to $[0,2\pi)$ because this is the physical
range of $q$ in our problem. To understand the error estimate here, note that
the overlap of $\chi_{0}(q)$ and $\chi_{0}(q-d)$ is exponentially small in
$\lambda$,

$\int_{-\infty}^{\infty}dq\
\chi_{0}(q)\chi_{0}(q-d)=e^{-\frac{d^{2}}{2\sqrt{2}}\lambda}\ .$ (154)

This means that the dominant contribution to $\int_{0}^{2\pi}dq\
\overline{\chi_{k}(q)}\chi_{k}(q)$ comes from the overlap between Gaussians in
the same position, while the overlap between Gaussians that are offset by some
amount accounts for the error term. The smallest possible offset is equal to
the period $\pi$ of the $\cos(2q)$ potential, and so the error estimate in our
expression for $\int_{0}^{2\pi}dq\ \overline{\chi_{k}(q)}\chi_{k}(q)$ follows
from taking $d=\pi$ in Eq. (154). Finally, for later use we remind the reader
that for our model the parameter $\lambda$ is proportional to $\xi/\ell$,
where the correlation length $\xi$ was defined in Eq. (139).

###  E.3 Calculating the Berry phase $\gamma_{k}$

We now calculate the Berry phase $\gamma_{k}$ in Eq. (52) using properties of
the eigenstates of our effective zero mode Hamiltonian $H_{\text{eff}}$. To
simplify the notation we denote $|\psi_{k}(B,0)\rangle$ by $|\psi_{k}\rangle$
in what follows.

To start, we note that
$\langle\psi_{k}|\tilde{p}|\psi_{k}\rangle=\langle\psi_{k}|p|\psi_{k}\rangle+\delta$,
and so we focus on evaluating $\langle\psi_{k}|p|\psi_{k}\rangle$. For this
matrix element we have

$\displaystyle\langle\psi_{k}|p|\psi_{k}\rangle$
$\displaystyle=-i\int_{0}^{2\pi}dq\
\overline{\psi_{k}(q)}\frac{d}{dq}\psi_{k}(q)$
$\displaystyle=-\delta\int_{0}^{2\pi}dq\ |\chi_{k}(q)|^{2}-i\int_{0}^{2\pi}dq\
\overline{\chi_{k}(q)}\frac{d}{dq}\chi_{k}(q)\ ,$ (155)

where we remind the reader that $\psi_{k}(q)=e^{-i\delta q}\chi_{k}(q)$. We
then use the approximate form (152) of the wave functions $\chi_{k}(q)$ at
large $\ell$ (more precisely, at large $\lambda$) to find that

$\int_{0}^{2\pi}dq\
|\chi_{k}(q)|^{2}=1+O(e^{-\frac{\pi^{2}}{8}\frac{\ell}{\xi}})$ (156)

and

$-i\int_{0}^{2\pi}dq\
\overline{\chi_{k}(q)}\frac{d}{dq}\chi_{k}(q)=O(e^{-\frac{\pi^{2}}{8}\frac{\ell}{\xi}})\
,$ (157)

where the error terms on these estimates come from the calculation of the
overlap of two shifted Gaussians (recall Eq. (154)). Therefore our final
result is that

${\langle\psi_{k}|p|\psi_{k}\rangle}=-\delta+O(e^{-\frac{\pi^{2}}{8}\frac{\ell}{\xi}})\
,$ (158)

and so we find that the Berry phase $\gamma_{k}$ is given (within our
variational approximation) by

$\gamma_{k}=-k\pi+O(e^{-\frac{\pi^{2}}{8}\frac{\ell}{\xi}})\ .$ (159)

The most interesting aspect of our result for $\gamma_{k}$ is that, for
$\ell\gg\xi$, the Berry phase is equal to the _topological_ value

$\displaystyle\gamma_{k,\text{top}}$ $\displaystyle=$ $\displaystyle-k\pi\ ,$
(160)

up to corrections that are exponentially small in the length $\ell$ (which is
also the separation between the fractional quasiparticles at the ends of the
FM region). These exponentially small corrections to topological Berry phases
are always expected in finite size systems, but they are very rarely
calculated explicitly. Our ability to capture these corrections here is a
significant demonstration of the power of our variational method.

##  Appendix F Generalization to $m\neq 1$

In this appendix we briefly explain the generalization of our results to the
fractional case of $m>1$ (i.e., a domain wall configuration at the boundary of
a fractional quantum spin Hall system). Recall from Appendix C that in the
$m=1$ case we were able to precisely construct bosonized fermion operators
$R(x)$ and $L(x)$ that obey the correct anticommutation relations of fermion
field operators. In contrast to that result, in the $m>1$ case we are not
aware of a precise construction of bosonized fermion operators $R(x)$ and
$L(x)$ that exactly obey the correct anticommutation relations. One possible
guess in this case is to define $R(x)$ and $L(x)$ via

$\displaystyle R(x)$ $\displaystyle=$ $\displaystyle\frac{e^{i\delta\frac{\pi
x}{\ell}}}{\sqrt{2\ell}}:e^{-im\phi_{\uparrow}(x)}:$ (161a) $\displaystyle
L(x)$ $\displaystyle=$ $\displaystyle\frac{e^{i\delta\frac{\pi
x}{\ell}}}{\sqrt{2\ell}}:e^{im\phi_{\downarrow}(x)}:\ .$ (161b)

With these definitions one still finds that $\\{R(x),R(y)\\}=0$ and
$\\{L(x),L(y)\\}=0$. However, the other anticommutators no longer exactly
match the expected answer for fermionic operators. For example, in the limit
of $\epsilon\to 0$, $\\{R(x),R^{\dagger}(y)\\}\neq\delta_{2\ell}(x-y)$ but is
instead equal to some more complicated distribution.888This fact about the
bosonized fermion operators in the fractional case does not seem to be widely
known. At least, we are not aware of any discussion of it in the literature.
Heuristically, the deviation between $\\{R(x),R^{\dagger}(y)\\}$ and
$\delta_{2\ell}(x-y)$ is due to a short length scale of the strongly
interacting system above which interacting fermion systems develops
topological order, which we can identify as the ultraviolet cutoff $\epsilon$
in the mode expansion (72) of the boson fields.

Because of this issue, in this appendix only we adopt a less precise (but
commonly used) definition of the bosonized fermion operators. Specifically, we
define $R(x)$ and $L(x)$ via

$\displaystyle R(x)$ $\displaystyle\sim$
$\displaystyle\frac{1}{\sqrt{2a}}e^{-im\phi_{\uparrow}(x)}$ (162a)
$\displaystyle L(x)$ $\displaystyle\sim$
$\displaystyle\frac{1}{\sqrt{2a}}e^{im\phi_{\downarrow}(x)}\ ,$ (162b)

where we have not used any normal-ordering prescription, and where we used the
ultraviolet cutoff $a$ (instead of the infrared cutoff $\ell$) to obtain the
correct dimensions. Loosely speaking, using this definition we have
$\\{R(x),R^{\dagger}(0)\\}=0$ if $x\neq 0$, and $\\{R(x),R^{\dagger}(0)\\}\sim
1/a\to\infty$, similar to a $\delta$-function.

We again carry out a variational calculation using a trial state
$|\Psi\rangle=|\psi\rangle\otimes|0,\zeta\rangle$, where
$\zeta=(\zeta_{1},\zeta_{2},\dots)$ is again chosen so that the expectation
value of the magnetic field term in the state $|\Psi\rangle$ is extensive. In
particular, in this case we find that

$\langle\Psi|R^{\dagger}(x)L(x)|\Psi\rangle\sim\frac{1}{2}f_{m}(x;a;\zeta)\langle\psi|e^{i2mq}|\psi\rangle\
,$ (163)

where the function $f_{m}(x;a;\zeta)$ is given by

$\displaystyle f_{m}(x;a;\zeta)=$ $\displaystyle{\frac{1}{a}e^{i\frac{\pi
mx}{\ell}}e^{\frac{i}{m}\sum_{n=1}^{\infty}\frac{e^{-\epsilon
n}}{n}\sin(2\kappa_{n}x)}}$ $\displaystyle\times
e^{-\sum_{n=1}^{\infty}\frac{e^{-\epsilon
n}}{n}e^{-2(\eta^{\prime}+\zeta_{n})}(1+\cos(2\kappa_{n}x))}\ ,$ (164)

where now

$e^{-2\eta^{\prime}}=K^{\prime}\equiv mK\ .$ (165)

Since the first line in Eq. (164) is equal to $1/a$ times a pure phase factor
(i.e., a complex number of unit modulus), we find that by performing the
summation in the exponent of the last factor of Eq. (164) with the same
variational scheme in Appendix D,

$f_{m}(x;a;\zeta)=\Delta_{Z}^{K^{\prime}}a^{K^{\prime}-1}f_{m,0}(x;a;\zeta),$
(166)

where $f_{m,0}(x;a;\zeta)$ is an order one quantity. This result is very
similar to the $m=1$ case from Appendix D, with the important difference that
$K$ is replaced by $K^{\prime}=mK$.

In this way we find that $|\psi\rangle$ should again be chosen to be the
ground state of an effective zero mode Hamiltonian, and in this case this zero
mode Hamiltonian takes the form

$H_{\text{eff}}=\alpha\tilde{p}^{2}-\beta\cos(2mq)\ .$ (167)

Following the self-consistency relation in Appendix D, we have

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\frac{vK^{\prime}\pi}{2\ell}$ (168a) $\displaystyle\beta$
$\displaystyle\sim$ $\displaystyle
B^{\frac{2}{2-K^{\prime}}}a^{\frac{2K^{\prime}-2}{2-K^{\prime}}}\ell.$ (168b)

We see that we again have $\alpha\propto\frac{1}{\ell}$ and
$\beta\propto\ell$, and so we again have

$\lambda^{2}=\frac{\beta}{\alpha}\propto\ell^{2}\ .$ (169)

The main difference between the analysis in this case and the analysis in the
$m=1$ case is that $H_{\text{eff}}$ (and the full domain wall Hamiltonian
$H(a)$) have a $\mathbb{Z}_{2m}$ symmetry instead of a $\mathbb{Z}_{2}$
symmetry. This symmetry is generated by the operator $e^{i\frac{\pi p}{m}}$,
and it can again be related to the conservation of the parity of the spin $S$
in the FM region. Indeed, in this case we have

$S=\frac{1}{2\pi}\int_{0}^{\ell}dx\
\left[\partial_{x}\phi_{\uparrow}(x)-\partial_{x}\phi_{\downarrow}(x)\right]=\frac{\tilde{p}}{m}\
,$ (170)

and the Hamiltonian commutes with $e^{i\pi S}=e^{i\frac{\pi\tilde{p}}{m}}$.
The Hilbert space of the model breaks up into sectors labeled by the different