Fractional Topological Superconductivity and Parafermion Corner States

We consider a system of weakly coupled Rashba nanowires in the strong spin-orbit interaction (SOI) regime. The nanowires are arranged into two tunnel-coupled layers proximitized by a top and bottom superconductor such that the superconducting phase difference between them is $\pi$. We show that in such a system strong electron-electron interactions can stabilize a helical topological superconducting phase hosting Kramers partners of $\mathbb{Z}_{2m}$ parafermion edge modes, where $m$ is an odd integer determined by the position of the chemical potential. Furthermore, upon turning on a weak in-plane magnetic field, the system is driven into a second-order topological superconducting phase hosting zero-energy $\mathbb{Z}_{2m}$ parafermion bound states localized at two opposite corners of a rectangular sample. As a special case, zero-energy Majorana corner states emerge in the non-interacting limit $m=1$, where the chemical potential is tuned to the SOI energy of the single nanowires.

Introduction.-The search for topological phases of matter has generated an enormous amount of research. Motivated by the discovery and classification of topological insulators (TIs) and topological superconductors (TSCs), the field has been driven by the desire to access phases with increasingly exotic properties. In particular, it has been found that the effects of strong electronelectron interactions can lead to exotic fractionalized phases, which are considered particularly interesting due to their potential use for topological quantum computation. However, only one-dimensional (1D) systems allow for an analytically tractable description of such strong interactions via the bosonization formalism. In order to study strongly interacting systems in more than one dimension, one therefore resorts to the so-called coupledwire approach, where two-or three-dimensional systems are built up from weakly coupled 1D channels, such as nanowires. This approach has proven to be exceptionally fruitful in accessing the fractional counterparts of several well-known topological phases such as fractional quantum Hall states [1][2][3], fractional TIs and TSCs [4][5][6][7][8][9][10][11][12][13][14][15], as well as fractional spin liquids [16][17][18].
Recently, a lot of interest has been raised by the generalization of conventional TIs/TSCs to so-called higher-order TIs/TSCs . While a conventional d-dimensional TI/TSC exhibits (d − 1)-dimensional gapless boundary modes, a d-dimensional nth-order TI/TSC hosts gapless modes at its (d−n)-dimensional boundaries. While electron-electron interactions have been taken into account in a few cases [40,41], the main focus was on noninteracting systems, in particular neglecting the possible existence of exotic fractional phases supporting emergent parafermions. This raises the question whether a coupled-wire approach can be used to extend the class of higher-order topological phases to the fractional regime. In this work, we show that this is indeed possible and explicitly construct a two-dimensional (2D) fractional second-order TSC exhibiting exotic parafermion corner states.
Our model consists of two layers of coupled Rashba nanowires with proximity-induced superconductivity of a FIG. 1. The setup consists of two layers of coupled Rashba nanowires where the index τ = 1 (τ =1) denotes the upper (lower) layer. The strength of the Rashba SOI associated with propagation along the x direction is given by α (−α) for the upper (lower) layer. Both layers are brought into proximity to an s-wave bulk superconductor such that there is a phase difference of π between them. In addition, the two layers are strongly coupled by interlayer tunneling of strength Γ. Neighboring nanowires of the same layer are weakly coupled via a spin-conserving hopping term of strength tz, via a spin-flip hopping term of strength β (−β) associated with Rashba SOI along the z direction, and via a crossed-Andreev superconducting term of strength ∆c (−∆c), where the last two terms are again of opposite sign for the two layers. Finally, a weak in-plane magnetic field of strength B is applied. We show that in specific regions of parameter space, this system hosts two zero-energy corner states (here denoted by γ1 and γ2) at two opposite corners of a rectangular sample. In the non-interacting case, these states are Majorana zero modes, whereas strong electron-electron interactions lead to Z2m parafermion corner states, where m is an odd integer depending on the position of the chemical potential.
phase difference of π between the upper and lower layers, see Fig. 1. In a first step, we show that in the presence of strong electron-electron interactions, such a setup exhibits a helical topological superconducting phase with gapless helical Z 2m parafermion edge modes propagating along the edges. Here, m is an odd integer determined by the position of the chemical potential µ. In the special case m = 1, where µ is tuned to the SOI energy of the single nanowires, Majorana edge modes emerge even in the non-interacting regime. At lower densities, the fractional regime m > 1 emerges in the presence of strong electronelectron interactions as the SOI and Fermi wavevectors get commensurable.
In a second step, we include a small time-reversal breaking perturbation in the form of a weak in-plane magnetic field to gap out the helical edge modes. For a finite rectangular sample, we find Z 2m parafermions localized at two opposite corners of the system depending on the direction of the magnetic field, which places our model in the class of 2D fractional second-order TSCs.
Model.-We consider two layers of coupled Rashba nanowires proximitized by bulk s-wave superconductors, see Fig. 1. Each nanowire of length L is modeled by a free-particle Hamiltonian Here, ψ † nτ σ (x) [ψ nτ σ (x)] creates (destroys) an electron at position x in the n-th wire in the layer τ ∈ {1,1} of spin σ ∈ {1,1}, where we define the spin quantization axis along the SOI direction. The Rashba coefficient α is taken to be of equal magnitude for all nanowires, but of opposite sign for the two layers. The SOI energy associated with propagation along the nanowire is E so = 2 k 2 so /(2m) for k so = mα/ 2 , and the chemical potential µ is defined relative to E so . The proximityinduced superconductivity is described by where we have set the phase difference between the two superconductors to π. This can, for example, be realized by the Josephson-junction setup shown in the inset of Fig. 1, where the phase difference between the two superconductors is adjusted by controlling the magnetic flux through the superconducting loop [58,59]. Alternatively, a thin insulating layer of randomly oriented magnetic impurities [60] could be placed between one of the layers and the corresponding superconductor such that the phase difference of π arises due to spin-flip tunneling via the impurities [61][62][63][64]. Furthermore, the two layers are coupled by interlayer tunneling of the form such that the total Hamiltonian describing an effective double nanowire (DNW) composed of two strongly cou-pled nanowires from different layers is given by H n = H 0,n + H ∆,n + H Γ,n . Finally, the DNWs are weakly coupled via a spin-conserving hopping term of strength t z , via a spin-flip hopping term of strength β (−β) associated with Rashba SOI along the z direction as well as via a crossed-Andreev superconducting term of strength ∆ c (−∆ c ), where the last two terms are again of opposite sign for the two layers. Here, |t z |, |β|, |∆ c | ≪ |∆|, |Γ|.
The interwire Hamiltonian can then be written as The total Hamiltonian is now given by H = n H n + H ⊥ , which in momentum space takes the form H = Here, τ i , η i , and σ i for i ∈ {x, y, z} are Pauli matrices acting in layer, particle-hole, and spin space, respectively, and a z is the spacing between neighboring nanowires. The system belongs to the symmetry class DIII [65] with time-reversal (particle-hole) symmetry given by T = iσ y K (P = η x K).
Helical topological superconducting phase.-We now demonstrate that the system can be brought into a helical topological superconducting phase hosting two counterpropagating Z 2m parafermion edge modes in the presence of strong electron-electron interactions. For this, we follow the method developed before for fractional TIs [4,8]: First, we solve the DNW Hamiltonian H n and demonstrate that, due to the interplay between ∆ and Γ, the elementary excitations are given by gapless Z 2m parafermion modes. We note that, in contrast to Refs. [4,8], there are two competing gap-opening mechanisms, such that when the system is brought close to the critical point Γ ≈ ∆, again half of the modes are left gapless. Second, we include weak hoppings between DNWs to gap out the parafermion modes in the bulk but leave Kramers pairs of gapless parafermion modes at the edges of the system. Again, if β and ∆ c counterbalance each other, the edge modes propagating along the x axis are perfectly localized at the outermost DNWs. Importantly, the topological phase is robust against deviations from these fine-tuned points, which will, however, lead to increased localization lengths of the edge states.
For illustrative purposes, we first consider the noninteracting regime with m = 1 and set µ = 0. To treat the DNW Hamiltonian H n , we linearize the spectra of the single nanowires around the Fermi points [66] as , L nτ σ (x) vary slowly on the scale of k −1 so and the Fermi momenta are given by k rτ σ F = (στ + r)k so [67]. We note that upon a change of basis defined byL nκν = (L nκν − iκνL nκν )/ √ 2,R nκν = (R nκν − iκνR nκν )/ √ 2, H n takes a block-diagonal form, while the structure of the Fermi momenta remains unchanged. For ∆ = 0, the exterior branchesR nκκ andL nκκ are fully gapped by superconductivity, whereas the interior branchesL nκκ andR nκκ have two competing gap-opening mechanisms given by interlayer tunneling and superconductivity. In the following, we thus focus on the interior branches only and tune the system to the critical point ∆ = Γ. In the new basis, the superconducting and tunneling term take the form H Γ,n = iΓ κ dxR † nκκLnκκ + H.c., , where the two decoupled sectors labeled by κ are related by time-reversal symmetry. Focusing on the first sector (corresponding to κ = 1), we find two gapless counterpropagating Majorana modes per DNW that can be written as Next, we add small interwire hopping terms [see Eq. (4)], where we set t z = 0 for simplicity. Focusing on the lowenergy sector spanned by the states given in Eq. (6), H ⊥ takes a form similar to a Kitaev chain [68] of coupled 1D modes, where N is the number of DNWs. At the special point ∆ c = β, the modes χ L11 and χ RN 1 do not enter H ⊥ and, thus, stay gapless in contrast to all other bulk modes. Obviously, the same is true for their time-reversal partners χ R11 and χ LN1 . Thus, the system is in a helical topological superconducting phase with Kramers partners of gapless Majorana modes propagating along the edges. Even though this result was derived using a considerable amount of fine-tuning, the topological properties of the system remain qualitatively identical for a broad range of parameters as long as the bulk gap does not close. In particular, our results do not change if a small t z is included, see Fig. 2(a). If, on the other hand, the system is infinite along the z axis and finite along the x axis, we apply the standard procedure of matching decaying eigenfunctions [69] to find a Kramers pair of gapless Majorana edge modes propagating perpendicular to the DNWs, see the Supplemental Material (SM) [70] for details and again Fig. 2(a) for a numerical confirmation.
Fractional helical topological superconducting phase.-Now we focus on the fractional counterpart of the helical superconducting phase discussed above. We In the presence of a small in-plane magnetic field, ∆Z > 0, we find Majorana bound states localized at two opposite corners of the system. The inset shows the spectrum confirming that these two states (red dots) are indeed at zero energy. The numerical parameters are N = 100, µ = 0, ksoL = 85, Γ/Eso ≈ 0.6, ∆/Eso ≈ 0.55, tz/Eso ≈ 0.01, β/Eso ≈ 0.28, ∆c/Eso ≈ 0.11, and, in (b), ∆Z /Eso ≈ 0.07 and φ = −π/16. We note that the found topological phases are stable against disorder and do not rely on spatial symmetries. tune the chemical potential to a fractional value µ/E so = −1+1/m 2 , where m is an odd integer. The new Fermi momenta are now given by k rτ σ F = (τ σ + r/m)k so . For m > 1, the interlayer tunneling term given in Eq. (3) no longer conserves momentum and is therefore suppressed. However, for the special values of chemical potential introduced above, momentumconserving terms can be constructed by including backscattering terms arising from electron-electron interactions [2,71]. These terms are given byH Γ, Similarly, we can write down a dressed superconducting termH ∆,n = ∆ κ,ν κν dx , where g B is the strength of a single backscattering process caused by electron-electron interactions, are assumed to be large [15,50,51,72]. In order to treat the interacting Hamiltonian analytically, we adapt a bosonized language [71]:R nκν (x) = e iφ1nκν (x) , L nκν (x) = e iφ1 nκν (x) for bosonic fields φ rnκν (x) satisfying standard non-local commutation relations. The dressed superconducting and tunneling terms can be simplified by introducing new bosonic operators where v is the Fermi velocity and we focus on the special values of Luttinger liquid (LL) parameters K nκν = 1/m. Again, half the modes are fully gapped by superconductivity, while for the other modes superconductivity and interlayer tunneling compete. Introducing conjugate fields ϕ nκ = (η 1nκκ − η1 nκκ − π/2)/(2 √ m), θ nκ = (η 1nκκ + η1 nκκ )/(2 √ m), the competing part of the above Hamiltonian can be rewritten as ForΓ =∆, this Hamiltonian corresponds to two timereversed copies of a well-known self-dual sine-Gordon model [73,74]. For m = 1, we thus expect to find a single gapless Majorana mode per time-reversal sector, which is consistent with our analysis of the non-interacting regime in the previous section. To study the more general case, we start by noting that for our choice of LL parameters, the competing terms have the same scaling dimension, which allows us to explicitly study the properties of the system along the self-dual line. For m > 1, however, the superconducting and tunneling terms are irrelevant to first order in the renormalization group (RG) analysis, suggesting a flow to a trivial LL fixed point. To resolve this issue, Ref. [72] argued that upon including a third-order term in the RG equations, a multicritical fixed point is encountered, which in our case separates a gapless phase, a phase dominated by superconductivity, and a phase dominated by interlayer tunneling. Such a fixed point has been shown to be described by a Z 2m parafermion theory [74], which means in our case that there are two bulk Z 2m parafermion modes related by time-reversal symmetry residing within each DNW. We now refermionize the above model in order to obtain an explicit expression for specific primary fields [74] of these parafermion theories. In particular, we define new composite chiral fermion operatorsψ commute with the superconducting and tunneling term, and the same is true for their Kramers partners χ rnκ , which prompts us to identify them as the ψ m primary fields of the Z 2m parafermion theories describing each DNW. Note that these fields are local in terms of electrons, which makes them particularly convenient to handle.
Similar to the non-interacting case, we introduce dressed interwire couplings for m > 1, which now couple theR nκν fields. Assuming that the interwire terms are relevant [70] and repeating the analysis of the integer case for the modes given in Eq. (10), we find that the bulk of the system is fully gapped, while there is a Kramers pair of gapless modes propagating along the edges of a finite sample. These modes correspond to ψ m primary fields of a Z 2m parafermion theory. However, it is expected [42] that there are indeed two full Z 2m parafermion theories residing at the edges of the system.
Majorana and parafermion corner states.-We now show that in the presence of a weak inplane magnetic field, the system enters a secondorder topological superconducting phase.
Let us start from the (non-interacting) Zeeman Hamiltonian H Z = ∆ Z n,τ,σ,σ ′ dx ψ † nτ σ [cos(φ)(σ x ) σσ ′ + sin(φ)(σ z ) σσ ′ ]ψ nτ σ ′ . For m > 1, momentum-conserving terms are once again constructed by including suitable backscattering processes, such that the dressed term then couples theR nκν fields. In the following, we focus on the regime where the magnetic field strength ∆ Z is small enough not to modify the bulk structure. However, as time-reversal symmetry is broken, the helical edge modes are gapped out. Assuming that the system size is large such that far away from the corners, all four edges can be treated independently, we calculate the projection of H Z onto the edge states for all four edges, see the SM [70]. If we label the edges of a rectangular sample by an index p = 0, ..., 3 in counterclockwise order starting from the bottom edge, the projection of the Zeeman Hamiltonian onto the edge p is given by where we have defined ϕ p = pπ/2 and γ y is a Pauli matrix acting on the low-energy subspace spanned, in this order, by the low-energy edge mode belonging to the time-reversal sector κ = 1 and its Kramers partner belonging to the sector κ =1. This shows that the mass term changes sign at two corners of the system. Explicitly, the sign change occurs at two diagonally opposite corners of the sample depending on the direction of the magnetic field. For φ ∈ (0, π/2) ∪ (π, 3π/2) [φ ∈ (π/2, π) ∪ (3π/2, 2π)] the sign change occurs at the top-left and bottom-right (top-right and bottom-left) corners. In the spirit of a Jackiw-Rebbi model [75], there are bound states at the corners where the mass term changes sign, which in our case inherit the exotic properties of the propagating modes and thus can be identified as zero-energy Z 2m parafermion corner states. Again, while this result was derived for the local ψ m fields, we expect that our arguments generalize to the full set of Z 2m primary fields. In the non-interacting limit m = 1, we find zero-energy Majorana corner states, which is verified numerically in Fig. 2(b).
Conclusions.-We have studied a system consisting of two layers of coupled Rashba nanowires in the presence of interlayer tunneling and proximity-induced superconductivity of a phase difference of π between the layers. We have shown that strong electron-electron interactions can stabilize a helical topological superconducting phase exhibiting Kramers partners of gapless Z 2m parafermion edge modes. Upon turning on a small in-plane magnetic field, the system enters a second-order topological superconducting phase hosting exotic zero-energy parafermion bound states at two corners of a rectangular sample depending on the direction of the magnetic field.
Our analytical approach is limited to the perturbative regime. However, if the system parameters are increased resulting in non-perturbative gaps, the parafermions will still be present as long as the bulk gap is not closed. Such non-perturbative regimes could be accessed numerically. Thus, our model provides a proof of principle for the existence of helical fractional TSCs and second-order fractional TSCs and we envision it to be a representative of a more general class of systems exhibiting the same parafermionic features. In this appendix, we explicitly write down the dressed interwire terms coupling the gapless parafermion modes found to reside within each DNW (see the main text). Let us start from the non-interacting case m = 1. Focusing on the interior branchesL nκκ ,R nκκ defined in the main text, the interwire term for t z = 0 reads From the m = 1 case, momentum-conserving terms can be constructed for m > 1 by including backscattering processes in a similar way as was discussed for interlayer hopping and superconductivity in the main text. Explicitly, we define the dressed terms as where again k = (m − 1)/2. As can easily be checked by changing to bosonic fields, the two competing terms in H ⊥ have the same structure as the two competing DNW terms introduced in the main text. Therefore, the arguments for the relevance of H ⊥ can be directly adapted from the corresponding arguments for H n . In terms of the composite fermions defined in the main text, the interwire Hamiltonian reads from where we can repeat the analysis of the non-interacting case.

APPENDIX B: EDGE MODES PROPAGATING ALONG THE z DIRECTION
To confirm the existence of helical edge modes propagating along the z direction, we assume that the system is finite along the x direction and infinite along the z direction and apply the standard procedure of matching decaying eigenfunctions. Once all terms in the Hamiltonian are dressed by suitable backscattering processes (see the main text as well as Appendix A), it is convenient to work directly in terms of the composite fermionsψ (m) nκν for general m. Changing to momentum space along the z axis, we write the problem in terms of the Fourier-transformed fieldsψ (m) kz κν and linearize the spectrum around the Fermi points [1], kzκν (x)] are again slowly varying right-moving (left-moving) fields. The total Hamiltonian separates into a part corresponding to the exterior branches and a part corresponding to the interior branches [2] given by H int = i vκ z ν z ∂ x +βsin(k z a z )κ z ν x +Γκ z ν y + [∆ + cos(k z a z )∆ c ]κ z η y ν y (S5) in the basisΨ Here, κ i and ν i for i ∈ {x, y, z} are Pauli matrices, and the two sectors labeled by κ are related by time-reversal symmetry. As in the main text, we focus on the regimeΓ ≈∆ andβ,∆ c ≪Γ,∆, and assume that all terms in the Hamiltonian are RG-relevant for all m. For∆,∆ c > 0, we then find Kramers pairs of zero-energy solutions at k z a z = π which are exponentially localized to the system edges at x = 0, L. To demonstrate this, we consider the Hamiltonians Kramers partners of zero-energy wave functions at k x = 0, which we label as Φ 0,± for the bottom edge (p = 0) and Φ 2,± for the top edge (p = 2) of the system in correspondence with the labeling of the edges used in the main text. We note that the DNW Hamiltonian exhibits an additional symmetry corresponding to the operator O 1 = η y ν z , which anticommutes with both interlayer tunneling as well as superconductivity. Furthermore, O 1 commutes with the particle-hole symmetry operatorP. For the edge states Φ Hence, we find where γ y is a Pauli matrix acting on the low-energy subspace spanned by Φ (m) 0,± . In order to calculate the effective Hamiltonian for the top edge, we note that our system is invariant under rotation around the y axis by an angle π, which leads to (S16) We now treat the edges along the z direction, where we use the properly normalized zero-energy wave functions found in Appendix B and label them as Φ (S17) Therefore, we arrive at Again, the effective Hamiltonian for the right edge can be obtained by exploiting the two-fold rotation symmetry of the system, which gives us Combining the above results, we arrive at the effective Hamiltonian given in Eq. (11) of the main text. Following Ref. [3], we conclude that there exist zero-energy bound states at the corners where the mass term changes sign. Importantly, we note that this argument is independent of any gauge choice. If one would naively multiply an arbitrary phase factor to a solution on a particular edge, the time-reversal relation between the two Kramers partners at this edge changes, while the corresponding relations stay unmodified for all other edges, which would then contradict the idea of the solutions being connected to form a single set of counterpropagating edge modes.