D ec 2 01 9 Majorana and parafermion corner states from two coupled sheets of bilayer graphene

We consider a setup consisting of two coupled sheets of bilayer graphene in the regime of strong spin-orbit interaction, where electrostatic confinement is used to create an array of effective quantum wires. We show that for suitable interwire couplings the system supports a topological insulator phase exhibiting Kramers partners of gapless helical edge states, while the additional presence of a small in-plane magnetic field and weak proximity-induced superconductivity leads to the emergence of zero-energy Majorana corner states at all four corners of a rectangular sample, indicating the transition to a second-order topological superconducting phase. The presence of strong electronelectron interactions is shown to promote the above phases to their exotic fractional counterparts. In particular, we find that the system supports a fractional topological insulator phase exhibiting fractionally charged gapless edge states and a fractional second-order topological superconducting phase exhibiting zero-energy Z2m parafermion corner states, where m is an odd integer determined by the position of the chemical potential.


I. INTRODUCTION
Over the last few decades, topological phases of quantum matter have been the subject of extensive studies, both in theory and in experiments. In particular, a lot of work has been dedicated to the description and classification of topological insulators (TIs) and topological superconductors (TSCs) in various spatial dimensions [1][2][3]. Recently, the generalization of conventional TIs and TSCs to higher-order TIs and TSCs has attracted strong interest . While conventional ddimensional TIs and TSCs exhibit gapless edge states at their (d − 1)-dimensional boundaries, nth-order ddimensional TIs or TSCs exhibit gapless edge states at their (d − n)-dimensional boundaries.
In the search for suitable platforms to realize topologically non-trivial physics, graphene and graphenebased systems [32,33] such as carbon nanotubes and bilayer graphene (BLG) have attracted particular attention. While the unusual low-energy properties of these systems make them interesting in their own right, they also have been proposed to support topologically nontrivial phases of matter, hosting, e.g., gapless edge states or localized Majorana zero modes [34][35][36][37][38][39][40][41][42][43]. Unfortunately, most of these proposals require strong spin-orbit interaction (SOI) as a crucial ingredient, whereas SOI is weak in standard graphene [44]. In the last few years, however, considerable experimental progress in creating van der Waals heterostructures has made it possible to induce strong SOI in graphene by proximity to transition metal dichalcogenides (TMDs) [45][46][47][48][49][50][51][52][53][54][55][56], which has led to renewed interest in graphene-based systems as promising candidates to realize topologically non-trivial systems in the lab. These considerations, together with the recent interest in higher-order topological phases of matter, have prompted us to devise a graphene-based system realizing second-order topological superconducting phases. In particular, we consider an array of coupled quantum wires arising in bilayer graphene due to electrostatic confinement [57][58][59][60][61], see Fig. 1. The combined FIG. 1. The model consists of two coupled sheets of ABstacked bilayer graphene subject to electrostatic confinement, such that effective 1D wires arise at domain walls between gates set to opposite voltages ±V0/2. The upper left panel shows the spectrum of an effective wire localized at one of the domain walls with the in-gap states highlighted in red. Note that due to SOI each of the bands is split into two shifted copies. The light green lines indicate the values of chemical potential that will be of interest in the remainder of this paper. We now consider an array of such effective wires, where a unit cell is defined as consisting of four wires, see the dashed box. The wires are weakly coupled via a layer-conserving hopping term between neighboring wires within the same unit cell (between neighboring wires belonging to different unit cells) of strength ty,τ (t ′ y,τ ), as well as via an inter-bilayer hopping term of strength tz. Note that in order to introduce a hierarchy of interwire terms, the wires are arranged in an armchair-like order. In particular, the setup shown here leads to ty,1 ≈ t ′ y,1 < tz < t ′ y,1 ≈ t y,1 , as the strength of the hopping terms naturally decreases with the separation of the wires. effects of competing interwire hopping terms, an in-plane magnetic field, and proximity-induced superconductivity lead to the formation of zero-energy corner states at all four corners of a rectangular sample. In the noninteracting case, these corner states are Majorana bound states. The major benefit of studying an array of coupled wires, however, is the additional possibility of including the effects of strong electron-electron interactions in an analytically tractable way [62][63][64][65][66][67][68][69][70][71][72][73]. Using bosonization techniques, we show that suitable interactions can drive the system into a fractional phase exhibiting zeroenergy Z 2m parafermion corner states for an odd integer m, placing our model in the class of fractional secondorder TSCs.
This paper is organized as follows. In Sec. II, we describe a model for our setup, which consists of an array of coupled wires arising in bilayer graphene due to electrostatic confinement. In Sec. III A, we show that this system is a topological insulator for a certain range of parameters. Section III B extends this result to the interacting case, showing that the system supports fractionally charged edge states for suitable values of chemical potential and sufficiently strong electron-electron interactions. In Sec. IV A, we then include suitable superconducting and magnetic perturbations to gap out the helical edge states found previously and show that, in the non-interacting case, the system is driven into a secondorder topological superconducting phase with Majorana corner states at all four corners of a rectangular sample. In Sec. IV B, we again extend our analysis to the interacting case and show that the system can be driven into a phase hosting exotic Z 2m parafermion corner states, where m is an odd integer depending on the chemical potential. We summarize our results in Sec. V.

II. MODEL
We consider a system consisting of two coupled sheets of AB-stacked bilayer graphene as shown in Fig. 1. Each layer of graphene is a honeycomb lattice consisting of two non-equivalent atoms A and B coupled by a hopping element t. Each sheet of bilayer graphene is then composed of two layers of graphene coupled by a hopping amplitude t ⊥ between the A atoms of the first layer and the B atoms of the second layer. The effective Hamiltonian for a single sheet of BLG in the momentum space is then given by where λ i , γ i , and η i for i ∈ {x, y, z} are Pauli matrices acting in valley, sublattice, and layer space, respectively, and v F is the Fermi velocity for electrons in graphene. Within each sheet, electrostatic confinement can be used to form effective 1D wires. In particular, we consider creating one-dimensional domain walls between gates set to opposite voltages ±V 0 /2. This leads to propagating onedimensional states localized around the region where the voltage changes sign [57]. We then consider an array of such effective wires, where the wires are arranged in an armchair-like order as shown in Fig. 1. For later convenience, we define a unit cell as consisting of four wires (two belonging to the upper layer, and two belonging to the lower layer). As such, each wire can be labeled by a unit cell index n as well as two indices (ν, τ ), where ν ∈ {1,1} denotes the position within the unit cell and τ ∈ {1,1} denotes the layer.
The case of a single effective wire without SOI [57] as well as in the presence of curvature-induced SOI [61] has been thoroughly analyzed in previous works. However, even though the curvature-induced SOI is considerably larger than the intrinsic SOI of standard graphene, it is still relatively small [74][75][76][77][78][79][80][81]. In order to access a regime with stronger SOI and avoid the need for curvature, we consider a van der Waals heterostructure combining layers of graphene and a TMD [45][46][47][48][49][50][51][52][53][54][55][56]. In this case, the proximity-induced SOI is of the form H so = αλ z σ z + α R (λ z γ x σ y − γ y σ x ), where σ i for i ∈ {x, y, z} is a Pauli matrix acting in spin space [47,50]. While the predicted values for α and α R vary across the literature and depend on the specific TMD that is used, we find that in our case, the Rashba-like term proportional to α R is suppressed by strong interlayer tunneling t ⊥ , which is why we focus on the term proportional to α [82].
In the following, we consider a step function potential of strength V 0 /2, where, without loss of generality, we assume V 0 > 0 and focus on the case where the direction of confinement is along the armchair direction of the graphene lattice. Adapting the results of Refs. [57,61] to our setup to include the SOI, we find that the bulk gap of the spectrum is given by V 0 , while there are eight in-gap modes per effective wire, see the spectrum shown in Fig. 1. Explicitly, the energies of the in-gap modes for wire ν in the layer τ are given by Here, λ ∈ {1,1} is the valley index, σ ∈ {1,1} is the spin projection onto the z axis that is determined by the SOI of strength α, which we take to be of equal magnitude but of opposite sign for the two sheets of BLG, and κ ∈ {1,1} is an additional subband index. For now, we tune the chemical potential in layer τ to µ τ = τ [V 0 /(2 √ 2) + α], while an alternative choice is described in Appendix A. Note that for these values of chemical potential, the sector κ = 1 (κ =1) corresponds to modes with small Fermi momenta close to k x = 0 (large Fermi momenta far away from k x = 0). As an additional simplification, we note that for α ≪ V 0 , the above energy spectrum is approximately linear for small k x (i.e., in the sector κ = 1) such that To proceed, we work in the regime of strong spinorbit interaction, which allows us to linearize the spectrum of each channel around the respective Fermi point.
The electron operator for the n-th unit cell can then be represented as Ψ n = κ,ν,τ,σ (R nκντ σ e ik (ντ )κντ σ are slowly right-moving [left-moving] fields and k λκντ σ F are the respective Fermi momenta. Note that here and in the following we implicitly assume the presence of weak but finite intervalley scattering, which is why we do no longer consider the valley degree of freedom to be a good quantum number but treat modes differing only in their valley index as right and left moving modes of the same species. The effective Hamiltonian describing the uncoupled wires can be written as where we have made use of the fact that the Fermi velocities of the different branches for a fixed κ are approximately the same given that α ≪ V 0 . Explicitly, we note from Eq. (3) that, for the values of chemical potential of interest to us, we approximately have We will now couple neighboring wires in various ways. Neglecting all fast oscillating terms, we consider the interwire Hamiltonian with 0 ≤ t y,τ , t ′ y,τ , t z ≪ α. Here, t y,τ (t ′ y,τ ) is a spinconserving intralayer hopping element between neighboring wires within the same unit cell (between neighboring wires belonging to different unit cells) and t z is a spinconserving hopping element between neighboring wires belonging to different layers. The strength of these hopping amplitudes can be controlled by varying the interwire distance as well as the strength and the shape of the confinement potential.
Furthermore, if a superconducting TMD such as NbSe 2 is used, superconductivity will be induced in the graphene bilayers [53]. The corresponding effective Hamiltonian then reads Additionally, we consider the effect of an in-plane Zeeman field along the x direction. Combined with intervalley scattering, which we assume to be present in the system with broken translational invariance, this term takes the form Finally, the total Hamiltonian is defined as In the remainder of this paper, we will focus on the regime ∆ sc , ∆ Z ≪ t y,τ , t ′ y,τ , t z such that the superconducting and Zeeman term can be treated as weak perturbations to the interwire terms. Numerically, however, our analysis can be extended to the nonperturbative regime, confirming that the found topological properties persist as long as the bulk gap is not closed.

A. Non-interacting case
In this section, we demonstrate that our model supports a TI phase with Kramers partners of gapless edge states propagating along the edges of a large but finite sample. For this, we set ∆ Z = ∆ sc = 0. As can immediately be verified from Eqs. (5)-(7), the branches with κ =1 are trivially gapped by interlayer hopping. As such, we focus on the sector κ = 1 in the following. For this sector, we find that the branches R n1ντ (ντ ) and L n1ντ (ντ ) are trivially gapped by interlayer hopping, whereas the different hopping processes compete for the branches R n1ντ (ντ ) and L n1ντ (ντ ) , see Fig. 2. In the following, we are interested in the regime t y,1 ≈ t ′ y,1 < t z < t ′ y,1 ≈ t y,1 . Such a hierarchy is natural for an armchair-like arrangement of the effective wires as shown in Fig. 1, as the strength of the hopping terms can be expected to decrease with the For simplicity, only the sector κ = 1 is shown. The branches R n1ντ (ντ ) and L n1ντ (ντ ) are trivially gapped by interlayer hopping of strength tz. For the branches R n1ντ (ντ ) and L n1ντ (ντ ) , on the other hand, the interlayer hopping term competes with intracell hopping of strength ty,τ and intercell hopping of strength t ′ y,τ .
separation of the wires. For simplicity, let us assume that t y,1 = t ′ y,1 = 0 and t ′ y,1 = t y,1 . By direct inspection of Eqs. (5)- (7), we find that the bulk of the system is fully gapped.
In order to find edge states in a system that is finite along the y direction and consists of N unit cells, we note that the parameter regime of our interest is in the same part of the topological phase diagram as t z ≪ t ′ y,1 , see Appendix B. In this limit, we find that the two modes R 11111 and L 11111 (R N 1111 and L N 1111 ) at the left (right) edge of the system stay gapless. The presence of these helical edge states is not affected by deviations from the above fine-tuned point as long as the bulk gap does not close.
Let us now assume that the system is finite along the x direction and infinite along the y direction. We apply the standard procedure of matching decaying eigenfunctions to find edge states propagating along the y direction [83]. The projection of the Hamiltonian H = H 0 + H ⊥ onto the sector κ = 1 can then be written in momentum space as H = ky dx Ψ † ky H(k y )Ψ ky , with the Hamiltonian density H(k y ) given by Here, ν i , τ i , σ i , and ρ i for i ∈ {x, y, z} are Pauli matrices acting in wire, layer, spin, and right/left mover space, respectively, and a y is the size of a unit cell in the y direction. Next, we focus on k y = 0 and a single edge of the system at x = 0. In order to satisfy vanishing boundary conditions, we require R ky 1ντ σ (0) = −L ky1ντ σ (0). From this condition, we find that, given t y,1 > t z , there are two exponentially decaying solutions localized to the edge of the system. These are given by in the basis of Ψ ky=0 , where we defined a = e −x/ξ1 and b = e −x/ξ2 e ikF x with ξ 1 = v 1 /(t y,1 −t z ) and ξ 2 = v 1 /t z . It is straightforward to verify that these edge states are Kramers partners and related by time-reversal via Φ − = −iσ y ρ x KΦ + , where K denotes the complex conjugation.
Putting together all of the above results, we conclude that our system is in a topological insulator phase with a Kramers pair of gapless edge states running along the edges of a large but finite sample.

B. Interacting case
Let us now address the construction of the fractional counterpart of the above phase. For this, we tune the chemical potential in layer τ to µ τ = τ [V 0 /(2 √ 2)+ α/m], where m is an odd integer and m = 1 reproduces the noninteracting case discussed above. Again, the interlayer hopping term given in Eq. (7) trivially gaps out the sector κ =1 corresponding to large Fermi momenta. Therefore, we again focus on the sector κ = 1.
As a first step, we note that for m > 1 the hopping processes between neighboring wires belonging to the same layer [see Eqs. (5) and (6)] no longer conserve momentum. However, momentum-conserving terms can be constructed by including single-electron backscattering processes arising from strong electron-electron interactions, see Fig. 3 for a graphical illustration in the case m = 3. Explicitly, the dressed interwire terms are given by Here, t  In the following, let us assume that the above terms flow to strong coupling in a renormalization group (RG) sense. This can always be achieved if their bare coupling constants are sufficiently large or if their scaling dimensions are the lowest ones among all possible competing terms. The original interlayer hopping term given in Eq. (7) does not commute with the above terms and therefore cannot order simultaneously. Instead, the interlayer term that commutes with the above terms is, to lowest order, given by the dressed term where again t . Following the same arguments as above, we assume that this term flows to strong coupling. The total interwire Hamiltonian in the interacting case is now defined as H In order to facilitate the analytical description of the interacting system, we introduce bosonic fields φ 1n1ντ σ (x) and φ1 n1ντ σ (x) defined as R n1ντ σ (x) = e iφ1n1ντσ (x) and L n1ντ σ (x) = e iφ1 n1ντ σ (x) . The fields φ rn1ντ σ (x) satisfy the non-local com- With this choice, R n1ντ σ and L n1ντ σ satisfy the proper fermionic anticommutation relations among themselves, while the commutation relations between different species can be satisfied by an appropriate choice of Klein factors [84], which we will not explicitly include here. The dressed interwire terms given in Eqs. (12)- (14) can be simplified by introducing new bosonic operators η rn1ντ σ (x) = m+1 2 φ rn1ντ σ (x) − m−1 2 φr n1ντ σ (x). The new fields obey the commutation relations [η rn1ντ σ (x), η r ′ n ′ 1ν ′ τ ′ σ ′ (x ′ )] = irmπδ rr ′ δ nn ′ δ νν ′ δ τ τ ′ δ σσ ′ sgn(x − x ′ ) and, for m > 1, they carry fractional charge e/m [65]. In terms of these new fields, the dressed interwire terms take the form (17) We note that the bulk of the system is now fully gapped, while for a system consisting of N unit cells the modes η 111111 and η1 11111 (η 1N 1111 and η1 N 1111 ) at the left edge (right edge) of the system stay gapless. These edge states carry fractional charges e/m, as expected for a fractional TI.
In order to study the emerging fractional edge states further, let us define new composite chiral fermion operators R from which we recover the non-interacting case for m = 1.
Indeed, H (m) ⊥ has the exact same form as in the noninteracting case, except that R n1ντ σ (L n1ντ σ ) is replaced by R    As in the non-interacting case, we also consider a semiinfinite geometry where the system is finite along the x direction and infinite along the y direction. By introducing the Fourier transforms R is the velocity of the composite fields.
By continuity, we therefore find that our system hosts a Kramers pair of fractionally charged gapless edge states running along the edges of a large but finite sample, which allows us to identify our system as a fractional topological insulator. This means that we can write an effective edge theory in terms of two conjugate bosonic fields η 1 and η1 with [η r (l), η r ′ (l ′ )] = irmπδ rr ′ sgn(l − l ′ ), where l is an edge coordinate which is defined mod 2[L + (N − 1)a y ] and runs along the edge of the sample in the counterclockwise direction [17].

A. Non-interacting case
In this section, we show that the terms H sc and H Z [see Eqs. (8) and (9)] can drive the system into a secondorder topological superconducting phase. Again, we start by treating the non-interacting case. Importantly, we consider ∆ sc and ∆ Z to be small enough not to modify the bulk gap structure. However, they may modify the low-energy behavior of the system by gapping out the helical edge states found above. This statement is confirmed explicitly by considering the effective lowenergy edge theory. We assume that the system size is sufficiently large such that far away from the corners all four edges can be treated independently. Crucially, we find that the Zeeman term does not open a gap at k y = 0 in the spectrum of the edges states propagating along the y direction. This can be verified explicitly by using the form of the edge state wave functions given by Eq. (11), for which we find Φ + |H Z |Φ − = Φ + |H Z |Φ + = Φ − |H Z |Φ − = 0. Alternatively, one can arrive at the same conclusion in a more general way by exploiting the symmetries of the system. Indeed, at k y = 0 the system has an additional symmetry represented by the operator O = ν z τ y σ z ρ x , which anticommutes with the Hamiltonian H(k y = 0) [see Eq. (10)] for t y,1 + t ′ y,1 = t y,1 + t ′ y,1 . In addition, Φ ± defined in Eq. (11) are eigenstates of O: OΦ ± = Φ ± . Furthermore, we find {H Z , O} = 0, which then implies Φ + |H Z |Φ − = 0. All other matrix elements are trivially zero as the edge states Φ ± are eigenstates of σ z , showing that the magnetic term indeed does not open a gap in the spectrum of the edge states propagating along the y direction. Therefore, these edge states can only be gapped by superconductivity, whereas in the spectrum of the edge states propagating along the x direction both mechanisms can in principle open a gap. If we choose |∆ Z | > |∆ sc |, the magnetic term dominates over the superconducting one such that it is responsible for gapping the edge states propagating along the x direction. In analogy to previous works studying domain walls between competing gapping mechanisms in systems with helical edge states [85][86][87], we find localized Majorana zero modes at the domain walls between the regions where the superconducting/magnetic term dominates, which in this case means at all four corners of the system. However, in contrast to previous works, we apply both the superconducting as well as the magnetic term throughout the entire system. Figure 4 verifies our results numerically. Importantly, our numerical analysis confirms that the corner states are robust against small deviations from the fine-tuned point t y,1 + t ′ y,1 = t y,1 + t ′ y,1 . In addition, we confirmed numerically that the zero-energy corner states are robust against disorder that breaks all spatial symmetries but neither closes the bulk nor the edge state gaps. This confirms that the Majorana corner states are protected purely by the particle-hole symmetry enforced by superconductivity, while the spatial symmetry O is not playing a crucial role.

B. Interacting case
The above results can be extended rather straightforwardly to the interacting case. In order to gap out the edge states given in Eq. (21), the Zeeman term as well as the superconducting term need to be dressed by interactions in the standard way. To lowest order, the terms that can open a gap in the edge state spectrum are given by with ∆ . Again, these terms have exactly the same form as in the noninteracting case, except that R n1ντ σ (L n1ντ σ ) is replaced by R n1ντ σ ). We assume that the above terms are relevant (either due to their bare coupling constants or their scaling dimension) but considerably smaller than the interwire terms entering H (m) ⊥ . In this case, these additional terms gap out the edge states without changing the bulk gap structure. We can now repeat the above symmetry argument to find that the magnetic term does not open a gap in the spectrum of the edge states propagating along the y direction; this gap is opened by superconductivity only. On the other hand, the edge states propagating along the x direction are gapped by the Zeeman term for |∆ (m) sc |. We are thus effectively dealing with domain walls occurring naturally at the corners of a fractional 2D TI, despite the fact that the superconducting and magnetic terms are uniform and act both simultaneously on the entire system. Given this analogy to domain walls, we can follow Refs. [85][86][87][88] and show that every domain wall between a region gapped by superconductivity and a region gapped by a magnetic field hosts a zero-energy parafermion bound state that is spatially localized to the domain wall, i.e., to the corner of the sample. To make this statement explicit in terms of the fields considered here, we rewrite the left and right moving fields η 1 and η1 describing the low-energy edge theory in terms of conjugate fields ϕ = (η 1 − η1)/(2m) and θ = (η 1 + η1)/(2m) with [ϕ(l), θ(l ′ )] = iπ 2m sgn(l − l ′ ). The dressed superconducting and magnetic terms given in Eqs. (22) and (23) projected onto the low-energy part of the spectrum now take the form H dl cos(2mϕ). Let us now label our edges of the system by s ∈ {0, ..., 3} starting from the right edge of the sample and proceeding in counterclockwise order. In the strong-coupling regime, we find that along the x (y) edges we have ϕ i = π m (p i +1/2) [θ i = π m (q j +1/2)] for p i , q j ∈ Z, where i ∈ {0, 2}, j ∈ {1, 3} label the respective edge. These operators satisfy [p i (l), q j (l ′ )] = im 2π sgn(l − l ′ ). If we label the corners of a rectangular sample by v ∈ {0, ..., 3} in counterclockwise order starting from the corner between edges 0 and 1, we can define operators acting locally on the corners as γ 2k = e iπ(p 2k −q 2k+1 )/m , γ 2k+1 = e iπ(p 2k+2 −q 2k+1 )/m for k ∈ {0, 1}. These operators commute with the Hamiltonian as they act on domain walls between segments gapped by competing mechanisms and satisfy Z 2m parafermionic commutation As such, we find a single zero-energy Z 2m parafermion corner state per corner.

V. CONCLUSIONS
We have considered a model based on two coupled sheets of bilayer graphene in the strong SOI regime. Electrostatic confinement is used to create effective 1D quantum wires, which are then tunnel-coupled in various ways. For a certain range of parameters, this system can be brought into a topological insulator phase characterized by the presence of a Kramers pair of gapless helical edge states. Furthermore, a small in-plane magnetic field and weak proximity-induced superconductivity drive the systems into a second-order topological superconducting phase with zero-energy Majorana corner states at all four corners of a rectangular sample. Even more interestingly, the fact that we are dealing with effective 1D systems allows us to take into account the effects of strong electron-electron interactions in an analytically tractable way. Using a bosonization approach, we have shown that for sufficiently strong electron-electron interactions and suitable values of chemical potential, the system can be brought into a fractional topological insulator phase with fractionally charged gapless helical edge states as well as into a fractional second-order topological superconducting phase hosting exotic Z 2m parafermion corner states, where m is an odd integer determined by the position of the chemical potential.
In particular, we envision the strong SOI in the graphene bilayers to be induced by proximity to a few layers of NbSe 2 , which at the same time also induces superconductivity into the system. The recent progress in fabricating van der Waals heterostructures puts such a setup well into experimental reach. From a more general perspective, we therefore believe that our system demonstrates the potential use of electrostatically generated arrays of effective quantum wires in bilayer graphene as designer platforms to realize topologically non-trivial physics.
On the other hand, we note that while gated bilayer graphene turns out to be a particularly convenient platform to realize the model proposed here, our results can The bulk gap closing line at tz = (ty,1 + t ′ y,1 )(t y,1 + t ′ y,1 ) marks the transition from the topologically trivial phase to the TI phase with a Kramers pair of gapless edge states. (b) Phase diagram of the second-order phase as a function of ∆Z and ∆sc. Importantly, we consider both ∆sc and ∆Z to be smaller than the bulk gap Egap, such that the bulk gap is never closed by these terms. However, the edge gap closes and reopens for ∆sc = ∆Z, corresponding to the phase transition between the topologically trivial phase and the second-order TSC (SOTSC) phase with one Majorana corner state per corner. The green dot at ∆Z = ∆sc = 0 corresponds to the first-order phase, while the light green line for ∆sc = 0 and ∆Z > 0 indicates the phase with partially gapped edge states (gapped along the x direction but gapless along the y direction). The other parameters are the same as in Fig. 4. readily be adapted to different realizations of coupled 1D wires with similar low-energy properties. In this Appendix, we comment on an alternative realization of the second-order TSC phase. For this, we tune the chemical potential to µ τ = V 0 /(2 √ 2) + τ α instead of the values chosen in the main text. By replacing κτ →κ, we note that the sectorκ now once again corresponds to the modes with small Fermi momenta close to k x = 0. The spectrum for the sectorκ = 1 is identical to the one shown in Fig. 2 and can be analyzed in the same way as in the main text. However, the sectorκ =1 has to be treated differently in this case. Indeed, the interlayer hopping term [see Eq. (7)] is not able to open a gap for the sectorκ =1 anymore, as it couples right (left) with right (left) moving modes. Therefore, the first-order topological insulator phase is not present in this case, as the bulk of the system is not fully gapped out. However, once superconductivity is taken into account, the gapless modes in the sectorκ =1 are trivially gapped out by superconductivity, while the sectorκ = 1 can be treated in the exact same way as before. Therefore, we find a second-order TSC phase with the same properties as in the main text.

Appendix B: Phase diagram in the non-interacting case
In the main text, we assumed t y,1 = t ′ y,1 = 0 and t ′ y,1 = t y,1 for analytical simplicity. However, we argued that the topological properties of the system stay qualitatively the same in an extended region of parameter space. In this Appendix, we confirm this statement by calculating a condition for the closing of the bulk gap for a more general choice of hopping amplitudes. We start by considering the first-order phase with ∆ sc = ∆ Z = 0. To simplify matters, we focus on the case of t y,τ , t ′ y,τ , t z ≥ 0, but our analysis can easily be extended to account for negative values of the hopping amplitudes as well. The bulk Hamiltonian is obtained from Eq. (10) upon replacing −i∂ x → k x , and has time-reversal symmetry expressed by T = iσ y ρ x K, where K denotes the complex conjugation. Therefore, our system belongs to the symmetry class AII [89]. As long as t ′ y,1 , t y,1 > t y,1 , t ′ y,1 , the bulk gap can only close at k x = k y = 0, where the eigenenergies are explicitly given by E 2,±,± = ± t y,1 + t y,1 + t ′ y,1 + t ′ y,1 ± (t y,1 − t y,1 + t ′ y,1 − t ′ y,1 ) 2 + 4t 2 z /2. (B2) Thus, we find that the bulk of the system is fully gapped except for the values t z = 0 or t 2 z = (t y,1 +t ′ y,1 )(t y,1 +t ′ y,1 ) These conditions define the potential boundaries between topologically non-equivalent phases. In our case, we can identify the region corresponding to the topologically non-trivial phase by checking for the existence of gapless edge states: In the main text, we argued that the system hosts a Kramers pair of gapless edge states for t y,1 = t ′ y,1 = 0 and t y,1 = t ′ y,1 ≫ t z . Therefore, we identify the region of the phase diagram for which t 2 z < (t y,1 + t ′ y,1 )(t y,1 + t ′ y,1 ) as the topologically nontrivial one. This is visualized in Fig. 5(a).
Let us now turn to the second-order phase. In the presence of a magnetic field, ∆ Z = 0, time-reversal symmetry is broken, while superconductivity, ∆ sc = 0, enforces particle-hole symmetry. This places our system in the symmetry class D [89]. Assuming that ∆ sc and ∆ Z are much smaller than the bulk gap, these terms cannot result in a closing of the bulk gap. However, we find that the edge gap closes at the points |∆ sc | = |∆ Z |, which separates the topologically non-trivial second-order phase with one zero-energy Majorana corner state per corner from the trivial one with no corner states. The corresponding phase diagram is shown in Fig. 5(b).