Symmetry-protected gates of Majorana qubits in a high-Tc higher-order topological superconductor platform

We propose a platform for braiding Majorana non-Abelian anyons based on a heterostructure between a d-wave high-Tc 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 arbitrary phase gates (including the π/8 magic gates) of the Majorana zero modes, all of which are robust and protected by symmetries. As many of the ingredients of our proposed platform have been realized in recent experiments, our results provide a new route toward universal topological quantum computation. Copyright M. F. Lapa et al. This work is licensed under the Creative Commons Attribution 4.0 International License. Published by the SciPost Foundation. Received 03-05-2021 Accepted 08-09-2021 Published 02-11-2021 Check for updates doi:10.21468/SciPostPhys.11.5.086


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 [1][2][3][4][5][6][7][8]. Adiabatic braiding and exchange of the MZMs generate Clifford gates in a topologically protected manner [9,10], and implementations of such operations have been studied in various platforms [11][12][13][14][15]. However, one drawback of the Majorana platform is that the Clifford gates are not powerful enough to achieve universal quantum computation [16]. It is well known that additional gates (i.e., the magic gate of a π/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 [17][18][19][20][21], 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. [20,21]).
Recently, the concept of topological insulators and superconductors has been generalized to higher-order topological insulators and superconductors . Protected by crystalline symmetries [39], 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 [56], the authors proposed a protocol through the manipulation of the Zeeman field and the pairing order parameter, a full braid (corresponding to π/2 rotations) between a pair of MZMs can be achieved (see also Ref. [57]). In another proposal [58,59], 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 [11].
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 [60,61] 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 Bi 2 Sr 2 CaCu 2 O 8+δ (BSCCO) that has recently been realized in monolayers [62], and a quantum spin Hall insulator, such as WTe 2 [63][64][65]. 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. In the context of higher-order topology, the corner Majorana modes are protected by mirror reflection symmetries together with time-reversal symmetry and particlehole 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 [66] 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 [67] to realize topological superconductivity [68]. 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 timereversal 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 vari-ous non-Abelian rotations within the degenerate ground state subspace. First, we show that rotating the in-plane Zeeman field by 2π 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 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. This includes 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), and hence are robust against random errors as long as the input for the phase angle is sufficiently precise.
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 WTe 2 are readily available 2d materials for d-wave superconductivity and quantum spin Hall effect. In particular, WTe 2 has been demonstrated [64,65] to have a U(1) spin axis needed for our purposes. 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 π/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 [46,61]. 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 2m parafermion modes [69][70][71][72][73] 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 though the non-Abelian phase is not topologically protected against unitary errors, 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 [74].
The remainder of this paper is organized as follows. In Sec. 2 we describe the setup of our proposed platform that hosts pairs of MZMs at its corners. In Sec. 3 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. 4 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. 5 we generalize our setup to that with a fractional quantum spin Hall insulator with m = 1, and show that the Berry phase accumulated using the same protocol corresponds exactly to the exchange statistics of 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.

Corner MZMs in a high-T c superconductor platform
Our platform is based on several recent proposals [60,61] 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 wellknown, 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, in which the corner region of the d-wave SC is rounded and spatially separated from the QSH layer. Alternatively, we note that nodeless d-wave superconductivity have been proposed for the high-T c monolayer superconductor FeSe/SrTiO 3 [75,76]. In addition, there are several proposals for corner pairs of MZMs with s-wave pairing [77,78]. All of these are free from the issue of single-particle tunneling.
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 [39]. 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.
The existence of the corner Majorana pairs can be demonstrated by analyzing the boundary states of the QSH and treating a superconducting gap ∆ and a Zeeman field 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(s) with spin up (in the z-direction) and a left-moving fermion L(s) with spin down (also in the z-direction). These operators have standard anticommutation relations, for example {R(s), R † (s )} = δ(s − s ), and they also obey periodic boundary conditions. The kinetic energy for this system takes the Figure 1: The schematics of the proposed setup between a QSH insulator and a dwave superconductor coupled via superconducting proximity effect. The corner of the sample is subject to a Zeeman field and hosts two MZMs. The thick green curve denotes the boundary of a d-wave superconductor, and the interior of the red curve denotes the Zeeman field region.

low-energy form
Here we have assumed that both fermions have the same velocity which is set to 1. (Note that, throughout this work, we will use s and s for coordinates along the edge of the 2D sample.) The superconducting gap term, being a spin-singlet one, takes the low-energy form on the edge Importantly, due to the d-wave pairing symmetry, the gap function is odd under mirror reflection, and when projected onto the edge, ∆(s) is an odd function (we choose the origin at the corner). Finally, the Zeeman term is projected to the edge as where B ≡ B x + iB y = |B|e 2iτ . The full Hamiltonian is the sum of all three of these terms, This Hamiltonian can be diagonalized in the standard way by constructing lowering op- with a non-negative energy E ≥ 0. Indeed, imposing this relation leads to the usual Bogoliubov-de Gennes equations for the "spinor" η(s) = (η 1 (s), η 2 (s), η 3 (s), η 4 (s)) T . Without loss of generality, taking ∆ and B as real (their constant phases can be absorbed into the definition of R and L), we get where Γ 1 , Γ 2 , and Γ 13 are 4 × 4 matrices defined as and where s x, y,z and τ x, y,z are the Pauli matrices. We note that {Γ 1 , Γ 2 } = {Γ 1 , Γ 13 } = 0, but [Γ 2 , Γ 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 [Γ 2 , Γ 13 ] = 0 to rotate to a basis in which Γ 2 and Γ 13 are both diagonal. The required unitary matrix U is given by If we define a new spinorη(s) via η(s) = U †η (s), then we find thatη(s) satisfies The key property of this new equation forη(s) is that it breaks up into two decoupled 2 × 2 Dirac equations with masses equal to M ± (s) = ∆(s) ± B. Just as in the Jackiw-Rebbi model, a fermion zero mode is associated with each domain wall in M + (s) and M − (s), i.e., where ∆(s) = B or ∆(s) = −B (see Fig. 2). For the profile given by solid lines in Fig. 2, there is a mass domain wall in M − (s) (marked by the blue dot to the left), and the zero-energy solution is where s 0 is the location where ∆ = B. Another MZM, located at M + (s) = 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 O η are Hermitian, and correspond to Majorana fermions. Therefore we find that, for odd ∆(s), there exist a pair of MZMs separated by a length , the length of the region where |∆| < 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 swap positions. We can also consider an alternative configuration in which ∆ is a constant and B(s) 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 4.2. By the same token, MZMs are nucleated at the nodes of ∆ ± B(s) = 0. These two MZMs switch positions when ∆ is tuned through zero. Later we will build upon this observation and propose a protocol for non-Abelian braiding of the Majorana modes.

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 ∆ and a Zeeman field B can be written in first quantized BdG form as where, e.g. s z τ 0 := s z ⊗τ 0 . Here s z = ±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 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 ∆ and B are absent. There are two U(1) subgroups: U(1) s generated by spin: and U(1) c generated by charge Note that U c,2π = U s,2π = −s 0 τ 0 is the fermion parity symmetry. The theory also enjoys a number of discrete symmetries. It is invariant under the time-reversal symmetry where K is the complex conjugation. In addition, there is a chiral symmetry As we will see in the next Section, C is an anti-unitary charge-conjugation symmetry for manybody states. [79] The BdG Hamiltonian (10) also has a particle-hole symmetry 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 symmetrỹ which is a symmetry of (10) for a uniform phase of ∆ (taken to be real without loss of generality). Such an anti-unitary symmetry which squares to one, along with C places the edge theory in class BDI, which admits a 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 byT . In addition, we note that U s,2π = −s 0 τ 0 , the fermion parity, is obviously still a symmetry.
In the opposite situation in which the magnitude of the unidirectional Zeeman field (say B = B xx ) is spatial dependent and has a domain wall and max(B) > |∆|, two of the Majorana modes from different corners move to the edge connecting the two corners. In this case, the composite symmetryT 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 combinationC With the composite chiral symmetryC, the edge theory additionally belongs to class AIII, which admits another 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 operatorCPT = s z τ 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 Rashba-type spin-orbit coupling, the spin of edge states L(s) and R(s) may depend 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. For the quantum spin Hall material WTe 2 , however, we note that recent theoretical [64] and experimental works [65] have shown that indeed there exists a spin axis for the edge states and a U(1) s symmetry at the microscopic level.
Similarly, while 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, 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 [80], relates the two symmetriesC andT .

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 = 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 ν = 1/m fractional quantum spin Hall (FQSH) state [81,82] 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 [83], the authors developed a similar bosonization apporach to Majorana zero modes for a noninteracting open system.

Review of bosonization
The edge of the FQSH state with an emergent U(1) s symmetry can be described by two bosonic fields φ ↑ (s) and φ ↓ (s) obeying the commutation relations In the K-matrix formalism for edges of Chern-Simons theories, this system corresponds to the matrix K = mσ z . Both fields φ ↑/↓ (s) are defined to have compactification radius 2π. This means that all physical operators must be invariant under the shift φ ↑ (s) → φ ↑ (s) + 2π, and likewise for φ ↓ (s). Then the allowed operators containing zero derivatives of these fields must be built from exponentials of the form e inφ ↑/↓ (s) for some integer n ∈ . The charge density current, and spin operator for the edge are defined to be We then find that the right-and left-moving electron operators for the free fermion case are given by where 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 = 1.
With this definition we find that acting with R(s) or L(s) lowers the total charge by one unit, as expected for an operator that annihilates a single electron. In addition, the anticommutation relation {R(s), R(s )} = 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 φ ↑/↓ (s) takes the form Here we have also incorporated a density-density interaction (with coupling constant g) between the spin up and spin down fermions. 1 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 ϕ(s) and ϑ(s) defined as which satisfy the commutation relation In terms of these fields we find that and that H 0 can be rewritten in the form where the renormalized velocity v and Luttinger parameter K are related to the coupling con- 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 as it is this modified Luttinger parameter that actually appears in H 0 . In this bosonized formalism, a superconducting mass term takes the form where ρ is the superconducting phase. Similarly, a ferromagnetic mass term takes the form where B x + iB y = |B|e 2iτ . 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 ϕ and ϑ transform under the emergent symmetries identified in the previous section: Despite their different forms for the BdG Hamiltonian, at the field theory level both C and T are antiunitary symmetries, since each flips the sign on one of the dual fields ϑ and ϕ. In particular C flips the sign of the charge but not the current, which is thus an anti-unitary charge conjugation symmetry. The composite symmetriesT = U s,π T andC = U c,π C now becomẽ Thus underT (C), the Zeeman (SC) term remains invariant. Under the T-duality of the free boson theory ϕ ↔ ϑ (which is an emergent symmetry when K = 1),T andC are exchanged, as well as the Zeeman and the SC terms.
In our analysis below we will mainly invoke U s,φ ,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.

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 , within which the Zeeman field is a constant, and take the limit in which the superconducting gap |∆| → ∞ in the superconducting region and thus the ϑ at the two ends of the magnetic region are completely pinned. Without loss of generality, we take In the magnetic region, the Hamiltonian is given by 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 ϕ and ϑ: where κ n = πn , and where we included the dimensionless ultraviolet cutoff ε to control the oscillator sums. Removing the cutoff corresponds to taking ε → 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 ϑ(s) obeys the boundary conditions from Eq. (32) with quantized winding number p ∈ . Importantly, from Eq. (23) we have the commutation relation indicating q and p are conjugate variables. As a result of the quantization of p, q is a compact variable with q ∼ q + 2π.
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: The operators {a n } are related to {b n } via a Bogoliubov transformation which we discuss in details in Appendix D. The Fock space structure is guaranteed byT symmetry; for example a symmetry breaking Zeeman term ∼ B z ∂ x ϑ 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 in whichp ≡ p + 1/2, and α, β are coupling constants renormalized by quasiparticle fluctuations. In Appendices D and F, using a variational approximation, we derive the coefficients for K < 2, and for a strong Zeeman field within the range where a is a short-distance cutoff, e.g., given by the underlying lattice. The results are given by Using Eq. (38), it is straightforward to see that here β α. 2 The Schrödinger equation with the Hamiltonian in Eq. (37) is known as the Mathieu's equation [84][85][86], 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 |ψ k 〉 of this Hamiltonian is labeled by lattice momenta k, 3 i.e. e i π mp |ψ k 〉 = e i π m k |ψ k 〉, which take quantized values inside the Brillouin zone. The lattice constant is π/m, and thus k ∈ [−m, m). Sincep takes half-integer quantized values, so does k. The offset 1/2 in the quantization ofp 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 UnderT ,p, k → −p, −k. From Eqs. (24) and (33), we see that physically the tunneling current and spin quantum numbers are directly related to the k via 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 or large B, we have β α and the 1d lattice is in a flat-band limit. As we show in Appendix E.1, in the limit β α (which is the same as Eq. (38)) the bandwidth is exponentially suppressed as Here the correlation length ξ for K < 2 is expressed as 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 2m parafermions [69][70][71][87][88][89].
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 ξ ∼ 1/B; see Eq. (9).
In the flat band limit β α, the gap separating the MZMs from excited states can be obtained by approximating Eq. (37) by expanding the cosine potential. We find for the excitation gap For the free fermion system, this expression reduces to the Zeeman energy ∼ B.
In the opposite limit with B 1 a a 2−K , 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 (±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 ∈ + 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 O( ). The analysis here further shows that the excitation gap results from the quantization of lattice momentum, and is given by We prove the two-fold degeneracy for a general β/α more rigorously in Appendix G. Interestingly, such a robust ground state degeneracy has been elucidated [90] from the perspective of a mixed 't Hooft anomaly between time-reversal symmetry (corresponding to our generalized time-reversalT ) and fermion parity symmetry ϕ → ϕ + π (generated by our U s,2π ) of field theories with a Θ-term at Θ = π. 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. [90].
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 ϑ ↔ ϕ andT ↔C.

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 ∆(s) = ±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 system. 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.
While we use terms such as "braiding" and "exchanging" the MZMs in what follows, it should be emphasized that they are distinct from their counterparts realized by directly moving the anyonic excitations. For example, in our "braiding" process, the physical locations of the MZMs are not changed, and in our "exchange" process, the MZM's exchange locations but they do not remain well separated during the process. More precisely, we use these terms to mean that the resulting non-Abelian Berry phase, when protected by symmetries, is identical to those from the braiding and exchanging anyons. In this sense our proposal is closely connected to holonomic quantum computing [91].

Full braid via 2π 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π 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 τ from 0 to π in Eq. (37). We will show that this correspond to a full braid of two Majoranas at a given corner.
Let |ψ k (B, τ)〉 be the ground state of H eff (B, τ) with lattice momentum k. According to Eq. (40), k = ±1/2. In other words, In the 1d lattice interpretation, the parameter rotation angle of the Zeeman field 2τ corresponds to a displacement of the periodic potential by an amount of τ, and thus the eigenstate can be expressed via a translation operator: However, from Eq. (46), the right hand side is not single-valued upon a full 2π rotation (τ = π) of the Zeeman field. For an unambiguous calculation of the Berry phase, one should choose the phases of the states |ψ k (B, τ)〉 so that these states are single-valued functions of B and τ (defined modulo π) 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τk corresponding to their lattice momenta From Eq. (46), this choice ensures that |ψ k (B, τ)〉 returns to itself when we wind the Zeeman field by 2π (translating τ by π), i.e., we have This is not the only possible choice, and it is well-known that the Berry phases that we obtain are invariant under any redefinition provided that the new states |ψ k (B, τ)〉 are also single-valued functions of B and τ in the relevant region of the parameter space. We now compute the Berry phases γ k picked up by the states |ψ k (B, τ)〉 during the 2π rotation of the Zeeman field. The Berry phases γ k are given by the standard formula and so using Eqs. (48) we find that 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. (37) there are two important limits -the tight-binding limit (β α) and weak periodic potential limit (α β). 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 〈p〉 approaches zero, while in the latter 〈p〉 approaches the lattice momentum k. Thus we have in the tight-binding limit, 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 (β α). This is the same condition as (38), i.e., Ba (a/ ) 2−K . In Appendix E.3 we directly compute the Berry phase for the state |ψ k (B, 0)〉 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. (52) is indeed exponentially small in . 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 = ±1/2, the result in (52) matches exactly the non-Abelian Berry phases accrued during a full 2π braiding of two MZMs. [92] We note that the full braiding of the Majorana modes have also been proposed in a similar higher-order topological superconductor platform in Ref. [56].

Single exchange via π 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π rotation of the Zeeman field within the x-y plane is equal (at large ) 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γ k for this process is then given by half of the value γ k for the full braid, We find thatγ ±1/2 for the two ground states differ by π/2, and this is exactly the relative Berry phase expected for a single exchange of two MZMs [10].
To achieve this, we consider "half-moon" paths of the Zeeman field B = (B x , B y ) denoted in Fig. 3. This path consists of a half circle from τ = 0 to τ = π/2, and a straight line with B y = 0 sweeping the B field back to the initial configuration passing through the origin. Crucially, along the arc B must be in the tight-binding regime, i.e., Ba (a/ ) 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 C B , and its image under inversion in the B x − B y plane, transversed clockwise, by C B ".
Let |ψ k (B)〉 = |ψ k (B x , B y )〉 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τ) and B y = B sin(2τ), and so |ψ k (B)〉 can be identified with the state |ψ k (B, τ)〉 that we defined in Eq. (48). The Berry phase γ k for the full 2π rotation of B can be written as the line integral where C B is the circular contour of radius B centered at the origin of the B x -B y plane. The integral expression for γ k can be split into two contributions as Here γ k and γ k " are the contributions to the total Berry phase from the half-moon paths C B and C" B , respectively. 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. (44)). Near the origin, the Zeeman field vanishes but a finite size gap still exists (see Eq. (45)). In addition, in order for the dynamical phases to cancel, the two ground states should remain degenerate, which is guaranteed by theT symmetry for any value of B, including at B = 0.
Provided that the Zeeman field is the only odd component under the T symmetry (the physical time-reversal symmetry preserved by the quantum spin Hall and d-wave superconductor; not to be confused withT ), the contour C B → −C" B under T (which reverses both the Zeeman field and the orientation of the contour). The Berry phase is obviously odd under time-reversal, and therefore, Combining Eqs. (55,56), we see that  3. The system must have the T symmetry in the absence of the Zeeman field, and thẽ T = U s,π T symmetry in its presence. This means that the phase difference between the two superconducting regions right outside the corner region must be π, which is naturally realized in our setup with a d-wave superconductor.

Phase gate via generic rotation of Zeeman field
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 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θ 0 , ending up at B x = B cos(2θ 0 ) and B y = B sin(2θ 0 ). In the second part of the path we traverse the straight line segment from B = (B cos(2θ 0 ), B sin(2θ 0 )) to the origin B = (0, 0). Finally, in the third part of the path we traverse the straight segment from the origin B = (0, 0) back to our starting point B = (B, 0). Crucially, we assume that the first part of the path, namely the curved segment that is traversed at constant magnitude |B| = B, is taken in the tight-binding regime. We will also assume that in the absence of B, the theory has U(1) s symmetry.
To calculate the Berry phase γ k (θ 0 ) in this process, we first consider the contribution from the two straight line paths, C B (0) and C B (θ 0 ). Since the only term in the system that violates the spin rotation symmetry U s,θ in Eq. (11) is the Zeeman field, and that the two paths C B (0) and −C B (θ 0 ) are related by the U(1) s symmetry, the total contribution to the Berry phase along the straight paths C B (0) + C B (θ 0 ) is zero.
Then γ k (θ 0 ) is exactly equal to the contribution from the curved part of the path, and so 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 We note that this argument can also be applied to the symmetry protection for the exchange process. The result (59) 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. The Berry phase is robust to small deformations of the arc as long as the flat-band condition is maintained.
2. In regions outside Ba (a/ ) 2−K or Ba (a/ ) 2−K , The tracks of magnetic field must be precisely related by a 2θ 0 rotation in the B x − B y plane. 3. The system must have U(1) s symmetry when there is no Zeeman field, and theT = U s,π T symmetry in its presence to protect the ground state degeneracy. Here it is achieved by applying the Zeeman field in the direction perpendicular to the spin axis.
Due to the ground state degeneracy it is also possible to choose as the initial state a superposition of k = ±1/2. In the next subsection we discuss such a situation where we choose a different basis for initial states.

Clifford and phase gates via manipulating Majorana modes within and across corners
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. (41): corresponding to the fermion parity at the corner 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 where k and k are the respective lattice momenta in the two corners. With this notation, the exchange operation in either corner leads to a Berry phase represented bỹ where σ z is the Pauli matrix in the Hilbert space of (62).
To realize Clifford gates, one additionally needs to achieve the non-Abelian unitary operatorγ = e i πσ x 4 .
In Ref. [10] this can be achieved by swapping different sets of Majorana pairs. Similarly, here we show thatγ is achieved by manipulating two Majorana modes across different corners. As we discussed in Sec. 2.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 byC symmetry. According to (59), such a duality is simply the usual ϑ ↔ ϕ duality in bosonization. In this basis, the two states forming the qubit are then eigenstates of 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) and its eigenstates are labeled by 〈σ x 〉 = ±1.
In order to induce Berry phases, we can similarly design contours in the complex plane of ∆ 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 ϑ ↔ ϕ duality, one can straightforwardly obtain that a "half-moon" contour leads to the non-Abelian phase given in (64). Notice here this value is topological and protected by the dualC symmetry, in which, as we showed 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ρ 0 in ∆, which is protected by the U(1) c symmetry. Thus we have two types of phase gates available, namely,γ which for θ 0 = ρ 0 = π/8 correspond to the magic gates [16]. We note that in Ref. [66] 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 ∈ {−1/2, 1/2} and their fractional spins (see Eq. (41)) S = ±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.

Braiding parafermion modes
In Sec. 4 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 obtaiñ This result is exactly the exchange statistics of 2m parafermions. [69][70][71][87][88][89] However, unless β α, 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 eliminate dynamical phases by precisely controlling the system parameters and the duration of the exchange process such that it is a common period for all m − 1 modes. However, the result is not topological protected against unitary errors induced by imperfect cancelation of dynamical phases.

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 (including the π/8 magic gates needed for universal topological quantum computing) 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 [62]), quantum spin Hall insulators [63] and two-dimensional ferromagnets [67], 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 [74]. 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.

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 SC = dkΨ † (k) (H 0 + H SC ) (k)Ψ(k), where Ψ † = (ψ † (k), ψ(−k)) and H 0 (k) + H SC (k) = sin k x s z σ z τ 0 + sin k y s 0 σ y τ z + (cos k x + cos k y + m − 1)s 0 σ x τ z + ∆ sin k x sin k y s y σ 0 τ y .
Here s x, y,z denotes the spin degree of freedom, σ x, y,z is a band index, and τ 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 T = is y K, diagonal mirror symmetries M a = s x σ z τ y , M b = s y σ y τ x , and a particle-hole symmetry P = τ x K. Such a Hamiltonian has been studied in Refs. [60,61] 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 M a,b . Higher-order topological crystalline insulators and superconductors with mirror symmetries have been classified in Ref. [39] based on a K-theory analysis by Shiozaki and Sato [93]. In our case, the reflection symmetry anticommutes with both time-reversal and particle-hole conjugation. According to the terminology in Ref. [39], it belongs to symmetry class DIII M ++ , which in terms of second-order topology admits a 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 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 The Zeeman field breaks both T and M a,b , but preserves the composite symmetry T M a,b . Importantly, the other Zeeman term ∼ B z s z σ 0 τ z is odd under this action and is forbidden. Together with the particle-hole symmetry P = τ x K, the Hamiltonian H = H 0 + H Z preserves composite chiral (anti)symmetries PM a,b = s z σ z τ z , which anticommutes with P. Such a phase belongs to class D PM − , which also admits a 2 classification. This indicates the MZMs in the absence of H Z remains robust, despite time-reversal symmetry being broken. However, as pointed out in Ref. [39], this 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 WTe 2 . 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.

B Mode expansion of the bosonic fields
In this appendix we explain in more detail the mode expansions for ϕ(s) and ϑ(s) that we use in our analysis in this paper. To obtain the mode expansion for ϑ(s), we first identify a complete set of functions of x ∈ (0, ) that also obey the boundary conditions Eq. (32), and then we expand ϑ(s) as a series in these functions with operator-valued coefficients. We then expand ϕ(s) in terms of a complementary set of functions (also with operator-valued coefficients) in such a way that ϑ(s) and ϕ(s) 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 † n , with n ∈ \ {0} (i.e., we have an oscillator variable for each integer n ≥ 1). These operators obey the standard commutation relations [q, p] = i and [b n , b † n ] = δ nn (with all other commutators vanishing). In addition, the zero mode operator q is a compact variable and is defined modulo 2π, while its conjugate momentum p is defined to have integer eigenvalues. This means that the Hilbert space H zm associated with the zero mode operators q and p is spanned by the states |s〉, s ∈ , which are eigenstates of p, p|s〉 = s|s〉, and with e ±iq |s〉 = |s±1〉. We can also define a basis |q〉 of eigenstates of q, with 〈q|s〉 = e iqs 2π , and we have 〈q|q 〉 = δ 2π (q − q ), where δ 2π (q − q ) = s∈ 1 2π e i(q−q )s is the 2π-periodic delta function. In terms of these operators, the mode expansions for ϕ(s) and ϑ(s) take the form where κ n = πn . Here we have set the twist for ϑ at the two ends as a generic δ; in the setup discussed in the main text, we have δ = 1/2. The exponential factor ε is a dimensionless ultraviolet cutoff, which we will discuss in details for m = 1 (Appendix C) and m = 1 (Appendix F). We can see that this cutoff serves to control the oscillator sums at high momenta κ n . Removing the cutoff corresponds to taking a → 0, and one can check that in this limit the fields obey the correct commutation relations [ϑ(s), ∂ s ϕ(s )] = [ϕ(s), ∂ s ϑ(s )] = πiδ(s − s ) (these follow from Eqs. (18) and the definition of ϕ(s) and ϑ(s) in terms of φ ↑/↓ (s)).
Finally, as in the main text, it is convenient to define the shifted zero mode momentum operatorp viap = p + δ .
This will be useful because almost all of our expressions will involve the shifted momentump instead of the original momentum p. Note that, since p is defined to have integer eigenvalues, the eigenvalues ofp lie in the set + δ (the integers shifted by δ).

C Bosonization in the m = 1 case
In this appendix we explain how to carefully define the fermionic operators R(s) and L(s) in terms of the bosonic fields φ ↑ (s) and φ ↓ (s) in the non-fractional case with m = 1. Specifically, we define R(s) and L(s) as normal-ordered exponentials of φ ↑ (s) and φ ↓ (s), and with a dimensionful prefactor that depends on the length . We then show that the operators R(s) and L(s) 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 and e X e Y = e X +Y e

C.2 Definition of the fermion operators
We now present the definition of the operators R(s) and L(s). We start with the mode expansions for the fields φ ↑ (s) and φ ↓ (s), which take the form (recall that we take m = 1) In addition, recall that κ n = πn and that ε = πa . For later use, we note here that φ ↓ (s)=φ ↑ (−x) (this relation actually holds for any m and not just m = 1). 4 Given these mode expansions, our definition of the fermion operators R(s) and L(s) is as follows. First, for any operator O of the form we define the normal-ordered exponential : e O : by Then our definition of R(s) and L(s) is In other words, the fermionic operators R(s) and L(s) are defined in terms of normal-ordered exponentials of φ ↑ (s) and φ ↓ (s), with an additional prefactor proportional to − 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 We now mention a few important properties of the operators R(s) and L(s). 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 → 0 limit. We also note that, in our open geometry (and with our choice of boundary conditions), the fields R(s) and L(s) are not independent but are actually related by the identity This identity can be derived by taking the Hermitian conjugate of our expression for R(−x), and by using the rearrangement identity which can be derived using Eq. (73). For our setup, however, we will focus on the interval (0, ), in which L and R † can be treated as independent fields.

C.3 Derivation of anticommutation relations
We now show that, in the limit ε → 0, the operators R(s) and L(s) that we defined actually do obey the standard anticommutation relations for fermionic fields. We start by deriving the anticommutator {R(s), 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 where we used the infinite series ln(1 − z) = − ∞ n=1 z n /n (valid for |z| < 1) to get from the first to the second line. Putting these results together yields the formula By examining the term in square brackets, which tends to e i πx − e i π y as a → 0, we can see that which is the expected anticommutator for a fermionic field with itself. Next, we consider the anticommutator of R(s) with R † ( y). For this calculation we define the operator A(s) by and so In terms of this operator we can rewrite R(s) and R † ( y) as Then, using similar rearrangement identities as in our previous calculation (using Eq. (73) again), we obtain the formulas and For the anticommutator we then find the formula where we defined the function d(x − y; a) by 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 From this expression we can see that where is the 2 -periodic delta function. Since this function is zero for x = y modulo 2 , and since all of the prefactors from Eq. (93) are equal to 1 when x = y modulo 2 , we then find that Therefore we have proven that, in the limit ε → 0, the operator R(s) obeys the standard anticommutation relations for a fermionic field operator. Since the anti-commutation relation (97) holds down to the smallest length scales, one can identify the ultraviolet cutoff ε in the bosonic theory as that in the fermionic theory, which is given by the lattice constant a as For the left-moving field L(s), we can use our results for R(s) and the relation (80) between L(s) and R(s) to immediately conclude that L(s) 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 ∈ (0, ), where the last line holds since x + y = 0 for x, y ∈ (0, ). For x, y ∈ (0, ) we also have Therefore, for our system with open boundary conditions, we find that the left-and rightmoving fermionic fields already have the correct anticommutation relations, and we do not need to include any extra Klein factors to ensure that R(s) and L( y) (and R(s) and L † ( 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(s) and L(s) 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〉 that satisfies p|0〉 = 0 and b n |0〉 = 0 for all n ∈ {1, 2, 3, . . . }. Using our previous rearrangement of the product R † ( y)R(s), we find that in this state we have and by taking the complex conjugate we find that To check that this formula makes sense, we can investigate its behavior in the bulk of the system, which corresponds to taking the limit → ∞ while keeping x − y, δ. 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 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(s) and L(s), we first find that Then, using our previous rearrangement of 〈0|R(s)R † ( y)|0〉, we find that If we now take the bulk limit then in this case we find that which is the correct bulk correlation function of a free left-moving fermion in one spatial dimension.

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. (33a) of Sec. 3.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 ϕ(s) and ϑ(s), 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. 2, this Hamiltonian takes the form Here, R(s) and L(s) are the fermionic operators from Eq. (78) that we constructed from the bosonic fields φ ↑ (s) and φ ↓ (s). For simplicity we also assume that the magnetic field points along the positive x-axis, so that τ = 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 ϕ(s) and ϑ(s). We find that The oscillator part of H 0 (a) can be diagonalized by making a Bogoliubov transformation to new oscillator variables a n defined by where the real parameter η is related to K as In terms of these new variables, we find that e −εn κ n (a † n a n + a n a † n ) .
For later use we also note the reverse Bogoliubov transformation, which allows us to express b n in terms of a n and a † n . As we mentioned in the main text, the Zeeman term ∼ B R † 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ã n (with Hermitian conjugatesã † n ). These will be related to the a n oscillators via the Bogoliubov transformation a n = cosh(ζ n )ã n − sinh(ζ n )ã † n , where we have allowed the parameter ζ n ∈ that determines the transformation to depend on the index n. Using these new variables, we can rewrite H 0 (a) in the form We also find that b n is related toã n via the relation which follows from identities for the hyperbolic trigonometric functions. We now discuss our choice of variational trial state. Let |0, ζ〉 be the Fock vacuum state annihilated by all theã n ,ã Note also that, as we are now working in terms of theã n oscillator variables, the full Hilbert space H tot of our domain wall model is equal to the tensor product where H F is the Fock space generated by the action of the raising operatorsã † n on the Fock vacuum |0, ζ〉, and H zm is the Hilbert space for the zero modes (q and p act within H zm ). The trial ground state |Ψ〉 that we consider respects the tensor product structure of the Hilbert space and it takes the tensor product form where |ψ〉 is a state in the zero mode Hilbert space H zm and |0, ζ〉 is the Fock vacuum for thẽ a n variables. The nontrivial part of our variational calculation is the problem of finding the parameters ζ n and the zero mode state |ψ〉 ∈ H zm that minimize the energy expectation value 〈Ψ|H(a)|Ψ〉. In fact, we will not carry out this optimization procedure completely on the ζ n parameters. Instead, we will use a heuristic argument to obtain the behavior of the energy expectation value with the correct choice of ζ n .
To proceed with the variational calculation we need to compute the expectation value 〈Ψ|H(a)|Ψ〉 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 where the second term here is the vacuum energy for theã n oscillators. For the Zeeman term we find that where the function f (x; a; ζ) is given by e −εn n e −2(η+ζ n ) (1+cos(2κ n x)) .
We note that the summation in the exponent of the last factor ∞ n=1 e −εn n e −2(η+ζ n ) is logarithmical, with the series effectively truncated by the factors e −εn and e −2ζ n . The first factor e −εn , with ε = πa/ , is a ultraviolet cutoff for the mode number, n /a. The second factor e −2ζ n , on the other hand, is due to the quasiparticle mass gap ∆ Z induced by the Zeeman field B. Obviously we have ∆ Z = B for free fermions, but in general ∆ Z is renormalized by interactions and is a parameter determined by the choice of {ζ n }. The e −2ζ n factor effectively serves as an infrared cutoff for the mode number n, as for high enough modes κ n ∆ Z the effect of the Zeeman field can be neglected. This means that the correct choice of ζ n sets a lower limit of summation, n ∆ Z . Therefore in the regime the following summation can be approximated by where we have used Eq. (110). Plugging (124) into (120) we obtain where f 0 is a well-behaved O(1) function. With this choice, we then define an energy scale β via and we find that at small a we have an extensive behavior For a free fermion system with K=1, this energy is ∼ B 2 . Using these results, we can now complete our calculation of 〈Ψ|H(a)|Ψ〉. We find that where the coefficients α and β are given by and we have omitted the c-number vacuum energy term for theã n oscillators. This result tells us that for our variational approximation the zero mode state |ψ〉 should be chosen to be the lowest energy state of the effective zero mode Hamiltonian We notice that in the limit β α, this Hamiltonian describes an approximate harmonic oscillator, with energy level spacing given by ∼ αβ. These energy levels form a Fock space of zero modes, now with mass ∼ αβ. We can identify the mass scale of the former zero modes with that of the oscillator modes: Self-consistency between Eqs. (129) and (131) fixes the scale of ∆ Z as Notably, we indeed recover ∆ Z = B for the free fermion case K = 1. With Eq. (132), the condition (123) translates to 1 a which requires K < 2. For K ≥ 2 a separate variational ansatz is needed, which we postpone to future studies. As a sanity check, we see that the condition β α we needed is precisely one of the conditions in Eq. (133). 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 α and β as As we discuss in the next appendix, H eff is closely related to Mathieu's equation, and so we can use known results on that equation to study H eff and solve our variational problem. In that appendix we present error estimates for various quantities, and those error estimates are exponentially small in . The key to obtaining that scaling for the error estimates is the fact that the parameters α and β in H eff satisfy the relation It is convenient to define a correlation length such that λ ∝ /ξ, and we have ξ/a ∼ 1 Ba which diverges at K → 2 obeying the familiar Kosterlitz-Thouless scaling behavior. The full Hamiltonian H(a) and the effective Hamiltonian H eff both have a 2 symmetry generated by the operator e iπp (the symmetry is 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 (±1) of e iπ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 2 symmetry by noting that e iπp is proportional to the parity e iπS of the total spin S in the FM region, since S is given explicitly by From a physical point of view this makes sense since the spin parity e iπ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 ∈ [−1, 1) ∩ ( + δ), such that any given state is an eigenstate of e iπp with eigenvalue e iπk . 5 Note that, since p has integer eigenvalues, there are only two such values of k in the set [−1, 1) ∩ ( + δ). For example, in the case of δ = 1 2 , which is our main interest, we have k ∈ {− 1 2 , 1 2 }. Our variational method can be used to study the lowest energy states of H(a) with all possible eigenvalues of e iπ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πp eigenvalues. We therefore define separate variational trial states |Ψ k 〉 = |ψ k 〉 ⊗ |0, ζ〉 for each allowed value of k, where |ψ k 〉 should be chosen to be the ground state of H eff in the sector of H zm with e iπp = e iπk .
Let E k be the energy of the ground state |ψ k 〉 of H eff in the sector with e iπp = e iπk , and let k 1 and k 2 be the two allowed values of k in the set [−1, 1) ∩ ( + δ). Our variational estimate for the energy splitting ∆E between the ground states of the domain wall model in the sectors with e iπp = e iπk 1 and e iπp = e iπk 2 is then given by 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 . Then our variational approximation predicts that the ground states of the domain wall model with different e iπp eigenvalues are very nearly degenerate for large . In addition, for the special case where δ = 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 (not just approximately degenerate for large ).

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. 3.2. We first apply known mathematical results on Mathieu's equation to understand the ground states of the effective zero mode Hamiltonian H 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 δ = 1/2 -there is an exact degeneracy at δ = 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 eff and estimate of the splitting ∆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 eff in the two sectors with different e iπp eigenvalue. We start by explaining the relation between H 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 H zm are invariant under the action of e i2πp , which is the operator that translates q by 2π, e i2πp qe −i2πp = q + 2π. Therefore, if |ψ〉 is any state in H zm , then its wave function ψ(q) = 〈q|ψ〉 is 2π-periodic in q, ψ(q + 2π) = ψ(q). Next, as we discussed in the previous appendix, H eff also commutes with the operator e iπp that translates q by π, e iπp qe −iπp = q + π. Accordingly, all eigenstates of H eff can be chosen to be eigenstates of e iπp , as we have discussed (and we actually labeled states by their eigenvalue of the closely related operator e iπp ).
Let |ψ〉 be an eigenstate of H eff with energy E. Then the wave function ψ(q) satisfies the Schrodinger equation where p has become the differential operator −i d dq . If we define a new wave function χ(q) by then we find that χ(q) satisfies where dq . This equation can be brought into the standard form of Mathieu's equation by dividing through by α = 0 to obtain where we remind λ 2 = β/α and E = E/α. In addition, the 2π-periodicity of ψ(q) implies that χ(q) obeys the periodicity condition To apply known results from the study of the Mathieu's equation, we need to study the behavior of χ(q) under translations by π, 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πp . In particular, for a state |ψ k 〉 that satisfies e iπp |ψ k 〉 = e iπ(k−δ) |ψ k 〉 (so that e iπp |ψ k 〉 = e iπk |ψ k 〉), we find that the corresponding function χ k (q) = e iδq ψ k (q) = e iδq 〈q|ψ k 〉 satisfies the periodicity condition and this simple relation explains why we chose to label our states by their eigenvalue of e iπp instead of their eigenvalue of e iπp . It is known from Floquet theory (similar to Bloch's theorem from condensed matter physics), that the spectrum of the Mathieu operator − d 2 dq 2 −λ 2 cos(2q) is divided into distinct energy bands. In addition, the eigenfunctions within each energy band are labeled by a wave number k ∈ [−1, 1), which corresponds to the Brillouin zone of a one-dimensional lattice with period π. An eigenfunction χ k (q) characterized by the wave number k obeys exactly the periodicity condition from Eq. (145). From this we see that the lowest energy state of H eff in the sector with e iπp = e iπk corresponds exactly to the eigenfunction labeled by k within the lowest band of the spectrum of − d 2 dq 2 − λ 2 cos(2q). Therefore, the energy splitting |E k 1 − E k 2 | between the two lowest energy states of H eff with different e iπp eigenvalues is certainly less than α times the width |W 0 (λ)| of the lowest band of the Mathieu operator − d 2 dq 2 − λ 2 cos(2q) (we multiply by α because E = αE).
An asymptotic formula for the width |W 0 (λ)| at large λ was obtained in Ref. [84] (see also Ref. [85] for a convenient summary of the properties of the spectrum of the Mathieu operator). It takes the form 6 The key feature of this formula is the factor of e −λ 8 . The presence of this factor implies that, when λ is large, the width |W 0 (λ)| of the lowest band is exponentially small in λ. 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 , where ξ ≡ λ 8 is the correlation length given by Eq. (136)). Thus, our variational approximation predicts that for large the energy splitting ∆E of the two ground states in our domain wall model is exponentially small in the length of the FM region. For the free fermion case with K = 1, the correlation length is given by ξ ∼ 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 δ = 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 δ.

E.2 Approximate form of χ k (q) at large λ
For the Berry phase calculation later in this appendix we will need to understand the form of the eigenfunctions χ k (q) of the Mathieu operator in the limit of large λ (we referred to this as the "tight-binding" limit in the main text). Therefore, in this subsection we review some known facts about χ k (q) in this limit.
When λ is large, the eigenfunctions of the Mathieu operator in its lowest band are wellapproximated by a weighted sum of Gaussians localized in each valley of the cos(2q) potential (see the proof of Theorem 1 in Ref. [86]). 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 − d 2 dq 2 + 2λ 2 q 2 , and it is well-known that the lowest energy eigenfunction of this operator is a Gaussian of the form where we have chosen the coefficient so that ∞ −∞ dq |χ 0 (q)| 2 = 1. Let χ k (q) be the eigenfunction in the lowest band of the Mathieu operator and obeying the periodicity condition χ k (q + π) = e ikπ χ k (q). At large λ, this eigenfunction is given approximately by the periodic sum which contains all translations of χ 0 (q) by integer multiples of π, with the translation by nπ accompanied by the k-dependent phase factor e iknπ . The factor of 2 is included here so that χ k (q) obeys the normalization condition 7 where the integral is restricted to [0, 2π) because this is the physical range of q in our problem. To understand the error estimate here, note that the overlap of χ 0 (q) and χ 0 (q − d) is exponentially small in λ, This means that the dominant contribution to 2π 0 dq χ k (q)χ 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 π of the cos(2q) potential, and so the error estimate in our expression for 2π 0 dq χ k (q)χ k (q) follows from taking d = π in Eq. (151). Finally, for later use we remind the reader that for our model the parameter λ is proportional to ξ/ , where the correlation length ξ was defined in Eq. (136).

E.3 Calculating the Berry phase γ k
We now calculate the Berry phase γ k in Eq. (51) using properties of the eigenstates of our effective zero mode Hamiltonian H eff . To simplify the notation we denote |ψ k (B, 0)〉 by |ψ k 〉 in what follows.
To start, we note that 〈ψ k |p|ψ k 〉 = 〈ψ k |p|ψ k 〉+δ, and so we focus on evaluating 〈ψ k |p|ψ k 〉. For this matrix element we have where we remind the reader that ψ k (q) = e −iδq χ k (q). We then use the approximate form (149) of the wave functions χ k (q) at large (more precisely, at large λ) to find that and where the error terms on these estimates come from the calculation of the overlap of two shifted Gaussians (recall Eq. (151)). Therefore our final result is that and so we find that the Berry phase γ k is given (within our variational approximation) by The most interesting aspect of our result for γ k is that, for ξ, the Berry phase is equal to the topological value up to corrections that are exponentially small in the length (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.

F Generalization to m = 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(s) and L(s) 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(s) and L(s) that exactly obey the correct anticommutation relations. One possible guess in this case is to define R(s) and L(s) via With these definitions one still finds that {R(s), R( y)} = 0 and {L(s), L( y)} = 0. However, the other anticommutators no longer exactly match the expected answer for fermionic operators. For example, in the limit of ε → 0, {R(s), R † ( y)} = δ 2 (x − y) but is instead equal to some more complicated distribution. 8 Heuristically, the deviation between {R(s), R † ( y)} and δ 2 (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 ε in the mode expansion (71) 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(s) and L(s) via where we have not used any normal-ordering prescription, and where we used the ultraviolet cutoff a (instead of the infrared cutoff ) to obtain the correct dimensions. Loosely speaking, using this definition we have {R(s), R † (0)} = 0 if x = 0, and {R(s), R † (0)} ∼ 1/a → ∞, similar to a δ-function.
We again carry out a variational calculation using a trial state |Ψ〉 = |ψ〉 ⊗ |0, ζ〉, where ζ = (ζ 1 , ζ 2 , . . . ) is again chosen so that the expectation value of the magnetic field term in the state |Ψ〉 is extensive. In particular, in this case we find that where the function f m (x; a; ζ) is given by where now Since the first line in Eq. (161) 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. (161) with the same variational scheme in Appendix D, where f m,0 (x; a; ζ) 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 = mK.
In this way we find that |ψ〉 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 Following the self-consistency relation in Appendix D, we have We see that we again have α ∝ 1 and β ∝ , and so we again have The main difference between the analysis in this case and the analysis in the m = 1 case is that H eff (and the full domain wall Hamiltonian H(a)) have a 2m symmetry instead of a 2 symmetry. This symmetry is generated by the operator e i π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 |ψ k 〉 should be chosen to be the ground state of H eff in the sector with e i πp m = e i πk m . By again exploiting the connection to the Mathieu's equation, 9 we find that 〈q|ψ k 〉 = ψ k (q) = e −iδq χ k (q), where now the function χ k (q) should be chosen to be the eigenfunction in the lowest band of the operator − d 2 dq 2 −λ 2 cos(2mq) that also satisfies the periodicity condition All of our previous results can now be carried over to this case. The only difference is that there are now small changes in the asymptotic formula for the width |W 0 (λ)| of the lowest band of the Mathieu operator, and the approximate form of the eigenfunction χ k (q) in the tight-binding regime of large λ. These quantities are now given by and where now Note that the new factors of m in |W 0 (λ)| can be understood from the expression for |W 0 (λ)| at m = 1 by making the change of variables q = mq in the Mathieu's equation with m = 1. Also, the factor of 1/ 2m in the expression for χ k (q) is again present to ensure approximate normalization when integrated over the interval [0, 2π), which is 2m times larger than the period π/m of cos(2mq). Using these new formulas we again predict (in the tight-binding regime) an exponentially small splitting between the ground state energies E k of H(a) in sectors with different values of k, where the new correlation length ξ m is of the same order as the correlation length in the m = 1 case. Finally, we find that the Berry phase γ k associated with the full 2π rotation of the in-plane magnetic field is given approximately by The main difference compared to the integer case is the presence of the factor of 1/m, indicating a fractional value for the Berry phase. We again find exponential suppression of the corrections to this topological value. One important point for this Berry phase calculation is that we now choose the phase of the state |ψ k (B, τ)〉 according to the formula and this choice will ensure that the states are single-valued along the path that we take through the parameter space. In particular, with this choice we will again have |ψ k (B,τ+π)〉=|ψ k (B,τ)〉.

G Exact two-fold degeneracy of the domain wall model at δ = 1/2
Our main interest in this paper is domain wall configurations in which the central FM region is surrounded by two SC regions with opposite signs of the superconducting mass ∆(s). In this case, the central FM region is described by our domain wall model with the parameter value δ = 1 2 . In this appendix we show that in this situation the domain wall Hamiltonian H(a) has an exact two-fold degeneracy of all of its eigenstates (and this holds for any value of the integer m). We explain this symmetry structure in the particular case that the in-plane magnetic field B points along the positive x-axis, as the Hamiltonian with a rotated B is unitarily equivalent to this case (and so the structure of the energy spectrum will be the same).
We define the operator Γ 2 by its action on the operators q, p, b n , and b † n that appear in the mode expansions of the bosonic fields ϕ(s) and ϑ(s) in our model. As we mentioned above, we also choose Γ 2 to be anti-unitary. We define Γ 2 in such a way that it squares to the identity operator, and we define its actions on q, p, b n , and b † n as: Therefore, Γ 2 only acts nontrivially on p. However, it can also act on other expressions by complex conjugation since it is anti-unitary. We also note that Γ 2p Γ 2 = −p, wherep = p + 1 2 at δ = 1 2 . With these definitions one can easily see that Γ 2 ϕ(s)Γ 2 = ϕ(s) and Γ 2 ϑ(s)Γ 2 = −ϑ(s). From Eq. (59) of the main text, this Γ 2 operator is precisely the antiunitary time-reversal sym-metryT .
These relations in turn imply that Γ 2 φ ↑ (s)Γ 2 = φ ↓ (s) and Γ 2 φ ↓ (s)Γ 2 = φ ↑ (s). Finally, these relations imply that the bosonized fermion operators R(s) and L(s) satisfy and so we find that Γ 2 does indeed commute with the domain wall Hamiltonian H(a) at δ = 1 2 . Finally, we investigate the interplay between Γ 2 and Γ 1 . We have = e −iπ(−p−1) , and so Γ 2 anticommutes with Γ 1 . This completes our demonstration of the three properties of Γ 2 that we stated above. As we mentioned above, this then implies an exact two-fold degeneracy of all of the eigenstates of the domain wall Hamiltonian H(a).

H Ground state degeneracy from the perspective of the 't Hooft anomaly
In this Appendix we show that the ground state degeneracy due to the Majorana pair in each corner can be viewed as a consequence of a mixed 't Hooft anomaly between the generalized time-reversal symmetry U s,π T and fermion parity symmetry U s,2π . We begin with the partition function of the corner region in terms of the boson fields ϑ and ϕ (we set v = 1 and keep a generic m), given by Z = DϕDϑe iS , where subject to the spatial boundary condition The parameter α is rather unusual and absent from most literature on bosonization, which we will explain and determine shortly. Recall that the first term arises from the insertion of complete sets of conjugate coherent states |ϕ〉 and |π〉 ≡ |∂ x (ϑ + α)〉, which gives the matrix element x 〈ϕ(x, t + d t)|π(x, t)〉〈π(x, t)|ϕ(x, t)〉 = exp d t i∂ x (ϑ(s) + α(s))∂ t ϕ(s) π ds . (182) Indeed, it is straightforward to verify that [ϕ(s), π(s )] = iπδ(s − s ).
In the above we have used and from this the parameter α(s) can be determined by noticing the Hilbert space constraint of the compactification ϕ(s) ∼ ϕ(s) + 2πm. This requires that Given the spatial boundary condition Eq. (181), we can choose Note that this procedure is essentially the same as the one adopted in Eq. (48). After integrating the ϑ field, this leads to the action Notice that compared to the usual sine-Gordon model, we have an additional term with Θ = π. After integrating over x ∈ [0, ), this is precisely a Θ-term in a 1d quantum field theory. The partition function Z = e iS has two symmetries, a generalized time-reversalT under which ϕ → ϕ, t → −t, and a translation ϕ → ϕ + π (for m = 1 this is fermion parity U s,2π ). In particular, the former symmetry is only realized at Θ = 0, π in the presence of periodic temporal boundary conditions. As pointed out in Ref. [90], such the theory Θ = π admits a 't Hooft anomaly between the two symmetries. To this end, we couple the spin up and down fermions with a gauge field ±A s , via ∂ ϕ → ∂ ϕ − A s , and we show that gauge invariance and time-reversal symmetry are incompatible.
After integrating out spatially oscillatory modes, we have With the cosine term, the gauge group is lowered from U (1) under such a transformation. Here the incompatibility of 2m and time-reversal of the quantum theory is characteristic of a 't Hooft anomaly. In general, the ground state degeneracy due to the 't Hooft anomaly can be proven by contradiction. Suppose there is a unique ground state, and then due to time-reversal symmetry of the partition function Z[A s ], the ground state must carry zero charge under the gauge field, since the (temporal) gauge field A s is odd under time-reversal. However, if so, the ground state path integral could not admit a gauge anomaly, since being charge neutral it would not respond to any gauge transformation. Therefore, the ground state must be degenerate.
In this special case of m = 1, the symmetry properties of Z[A s ] from Eq. (191) can be captured by [90] Z[A s ] ∼ exp i d tA s /2 + exp −i d tA s /2 , (192) which indicates that the ground state is two-fold degenerate in the absence of the background gauge field. Each ground state carries a fractional charge ± 1 2 , and therefore, the gauge group is represented projectively, or equivalently as a double cover. While classically the 2 timereversal symmetry and the 2 gauge symmetry combines to 4 , at a quantum level the symmetry group is 8 .
Recalling that spin-1/2 fermions are charged ±1 objects under A s , we conclude that the two ground states have spin S = ± 1 4 . This indeed agrees with the results from the main text using (41)