Time-Reversal Invariant Topological Superconductivity in Planar Josephson Bijunction

We consider a Josephson bijunction consisting of a thin $SIS$ $\pi$-Josephson junction sandwiched between two-dimensional semiconducting layers with strong Rashba spin-orbit interaction. Each of these layers forms an $SNS$ junction due to proximity-induced superconductivity. The $SIS$ junction is assumed to be thin enough such that the two Rashba layers are tunnel-coupled. We show that, by tuning external gates, this system can be controllably brought into a time-reversal invariant topological superconducting phase with a Kramers pair of Majorana bound states being localized at the end of the normal region for a large parameter phase space. In particular, in the strong spin-orbit interaction limit, the topological phase can be accessed already in the regime of small tunneling amplitudes.

In contrast, time-reversal invariant (TRI) topological superconductors hosting Kramers pairs of MBSs  are still lacking experimental evidence. Such schemes are attractive since they avoid magnetic fields and their detrimental effect on superconductivity. Many such proposals rely on unconventional superconductors or can only be achieved with strong electron-electron interactions [56,57]. A major drawback of them is that they do not allow for an easy in-situ control over the topological phase.
It is the goal of the present work to close this gap. For this we propose a Josephson bijunction which can be brought into the TRI topological superconducting phase via the control of experimentally accessible parameters, thereby hosting multiple Kramers pairs of MBSs, see Fig. 1. The setup consists of a superconductor-insulatorsuperconductor (SIS) Josephson junction sandwiched between semiconducting layers with Rashba spin-orbit interaction (SOI) [see Fig. 1(a)]. The superconductors are assumed to be s-wave and to have a phase difference φ controlled by applying a magnetic flux or generated by introducing an additional layer of randomly oriented magnetic impurities [58][59][60][61][62]. Due to the proximity effect this leads to the formation of two tunnel- coupled superconductor-normal-superconductor (SN S) Josephson junctions of width W and phase difference φ [see Fig. 1(a)]. This setup leads to an intricate interplay between the formation of Andreev bound state bands (ABSBs) in both SN S junctions on one hand and the hybridization of the two Rashba layers on the other hand, with striking consequences. In particular, we find a periodic closing and reopening of the topological gap in the spectrum of the ABSBs as a function of k so W , where k so is the SOI momentum which can be tuned electrically via gates [63][64][65]. Interestingly, in the strong SOI limit, the TRI topological phase can be reached for a relatively small tunneling amplitude. We believe that with the recent advances in the fabrication and study of bilayer systems [66][67][68][69] and van der Waals heterostructures [70][71][72][73][74], the setup proposed here lies within experimental reach and could be a promising route in realizing Kramers pairs of MBSs without the need of magnetic fields.
Model. The Josephson bijunction consists of two semiconducting layers labeled by τ = ±1. Due to the structural asymmetry the Rashba SOI vectors α τ are naturally antiparallel and aligned along the z axis, which is normal to the layers. The layers are modeled by where the in-plane momentum operator, and σ i are the Pauli matrices acting in spin space. The operator ψ † τ σ (r) creates an electron with spin projection σ along the z axis in the τ -layer at the position r = (x, y). The SOI energy E τ,so = ̵ h 2 k 2 so,τ 2m, with the SOI momentum k so,τ = mα τ ̵ h 2 , is the energy difference between the bottom of the Rashba bands and the degeneracy point at k = 0, from which µ τ is measured.
The regions that are in contact with the bulk superconductors become superconducting as well [20,25,26,28] and the proximitized Rashba layers themselves form an SN S junction with a normal region of width W [75]. The superconducting regions are described by [76] where ∆ τ (y) = ∆e isgn(y)φ 2 θ( y − W 2) 2, ∆ > 0 is the strength of the induced superconducting pairing and φ is the phase difference between the superconducting regions inherited from the parent superconductors. We also allow for electron tunneling between the two layers, where we assume the tunneling amplitude to be uniform in the x direction such thatΓ(y) = Γθ(W 2− y )+Γ ′ θ( y − W 2), and without loss of generality Γ > 0 (Γ ′ ≥ 0) for the coupling between the two layers [83]. The total Hamiltonian modeling the setup is then given by First, we study the system numerically for two geometries [see the Supplemental Material (SM) [84] for details]: (i) semi-infinite geometry -the system is translationally invariant along the direction of the junction, i.e., along the x-direction, which allows us to parametrize the eigenstates by the good quantum number k x ; (ii) finite geometry -the system is finite in both the x-and ydirection. Second, to gain a better physical understanding, we then treat the problem also analytically.
Spectrum of Andreev Bound State Bands. In order to better understand the evolution of the ABSBs, it is instructive to consider the effect of a potential barrier of height V at the interfaces between the superconducting and the normal regions (y = ±W 2 ) [85]. If V is very large, the interfaces are intransparent such that there are two independent superconducting regions and an isolated normal region, which consists of two tunnel-coupled quantum wires. Since the physics of the normal region in the confined direction (y-direction) is equivalent to a particle in a box, the spectrum consists of quantized transverse subbands labeled by an index n, where the spacing between them depends on the width W and roughly behaves as 1 W 2 . Out of these subbands only the ones crossing the chemical potential or lying within the superconducting gap are important for us. By increasing W the subband spacing is reduced and as a result more of these will be shifted into the superconducting gap [see Fig. 2a]. The probability density of the n-th subband, ψ (n) (y) 2 , has n peaks, which allows us to determine n numerically from the shape of the wavefunction. Upon lowering V , the electron and hole bands get cou- pled with corresponding anti-crossings and, as a result, ABSBs form [see Fig. 2b]. This picture qualitatively explains why the number of ABSBs increases with W , and that generically the ABSBs have contributions coming from different transverse subbands of the normal region [see Fig. 2]. For Γ ≠ 0, the ABSB spreads in both layers. More concretely, when a new subband is entering the superconducting gap, the ABSB highest in energy will obtain mainly contributions from this additional subband around k x = 0, while at larger momenta k x it has mainly contributions from a subband with smaller n. As W is further increased the same behavior can be observed for all ABSBs [see Fig. 2]. For the remainder we will put V = 0 and work with ABSBs.
Kramers pairs of MBSs. If the phase difference between the superconductors is φ = π, the system is timereversal invariant and is in the DIII symmetry class and has a Z 2 number classification [86]. The time-reversal (particle-hole) symmetry operator is given by Θ = iσ 2 K (P = η 1 K), where η i are the Pauli matrices acting in particle-hole space, and K is the complex conjugation operator.
To investigate topological phase transitions, we look for gap closings of the ABSBs at k x = 0, which only happens at particular combinations of k so W and Γ ∆. An analytical expression for the gap closing condition can be derived in the semi-infinite geometry in the strong SOI limit, i.e. Γ, ∆ ≪ E so , for α 1 = −α1 and Γ ′ = 0 (see SM [84] for details). In this limit, we find a remarkably simple expression for the gap closing condition We stress that it is essential for obtaining this result to go beyond standard linearization and retain band curvature effects. Thus, for a fixed ratio Γ ∆ there exist multiple solutions, whenever the cotangent is positive [see Fig.  3a]. In the interval where the cotangent is negative, e.g. for k so W ∈ (π 2, π), no gap closing occurs and the ABSBs are always gapped [see Fig. 3a]. Interestingly, in this region of the phase diagram the topological phase can in principal be achieved with an infinitesimal tunneling amplitude Γ. Relaxing the condition of exactly opposite SOI in the two layers, one finds that the functional form of the gap closing lines remain essentially unchanged [see Fig. 3]. However, for weak SOI, ∆ ∼ E so , the functional form deviates from the cotangent behavior, and, more importantly, there are no longer regions where no gap closing takes place [see Fig. 3b]. Note that the position of the gap closing lines depend on the product k so W , and it is thus possible to tune into a particular phase by changing the SOI α, e.g. by applying electric fields via external gates.
The closing and reopening of the gap of the ABSBs is accompanied by the emergence of a Kramers pairs of MBSs at each end of the normal region [see Fig. 4], thus these gap closing points mark topological phase transitions. Due to their topological nature they are robust against potential disorder, as we checked numerically. Importantly, the emergence of the MBSs does neither rely on the fine-tuned point α 1 = α1 nor on Γ ′ = 0. It turns out that MBSs exist for both regimes Γ > Γ ′ and Γ < Γ ′ [see Fig. 4], which implies that the relative strength Γ Γ ′ is not crucial to reach the topological phase. However, it is essential that Γ is non-zero.
As has become clear from the discussion above, for fixed coupling strength Γ and increasing width W , the lowest lying ABSB will change its 'orbital content' around k x = 0 and will have contributions from transverse subbands with larger n while the contributions from subbands with smaller n are 'pushed' to higher values of k x [see pair of MBSs [see Fig. 5]. The wavefunctions of the MBSs show that they are localized on one end of the normal region as is expected. In the y-direction, however, the states that emerge for larger W have n peaks, revealing that they stem from the topological phase transition of the corresponding subband.
Having discussed the physics around k x = 0, it is also important to look at finite momenta, as the localization length of the MBSs is determined by the smallest gap. Gaps at higher momenta, where the band has contributions from lower subbands, become smaller with increasing W and thus, the localization length of MBSs coming from these subbands increases [87,88]. This delocalization leads to a larger overlap of their wavefunction and thus they split in energy [see Fig. 5]. Typically, when a new Kramers pair of MBSs emerges, the MBSs coming from subbands with lower n are already split due their overlap, and there is only one Kramers pair of MBSs at each end that is at zero-energy. Thus, disorder induced coupling between an even number of Kramers pairs is suppressed since they are no longer energetically degenerate. As we have verified numerically, the phases with an even number of Kramers pairs of MBSs are indeed quite robust against potential disorder. Only in the limit where the junction is much longer in the x-direction than any of the localization lengths, the MBSs on opposite ends do not overlap, and one can always find a term that gaps out all MBSs without violating any symmetries. Strictly speaking only the phases with an odd number of Kramers pair of MBSs are topologically protected and thus stable against any symmetry-respecting perturbations.
Breaking time-reversal symmetry. If TRS is broken by a Zeeman field along the x-direction, the system behaves similarly to what is known in topological double nanowires setups [47]. Increasing ∆ Z continuously from zero brings the system into a trivial phase until a critical field value is reached and the gap of the ABSBs closes. Increasing ∆ Z further leads to a reopening of the gap and the system enters a topological phase with a total of two MBSs, one MBS at each end of the normal region, see SM [84]. Importantly, there exist regions where the topological phase can be reached with relatively small Zeeman fields [see In the main text we present numerical results for the semi-infinite and the finite geometries in which the Rashba bilayer system is assumed to have a rectangular shape. In this section we explicitly discretize the total Hamiltonian H defined by Eqs. (1)-(3) of the main text.

Semi-infinite geometry
In the semi-infinite geometry, we assume, without loss of generality, that the system is translationally invariant along the x and finite along the y direction with the length L y = (N y − 1)a, where N y is the number of lattice sites in y-direction and a the lattice constant. The width of the left superconducting region is given by L 1 = (N 1 − 1)a and the width of the normal region is W = (N W + 1)a. The total Hamiltonian for the semi-infinite geometry is given bȳ The operator c † kxτ σn creates an electron with momentum k x and spin projection σ (along z-axis) in the layer τ at the lattice site m. Here, t is the amplitude for a hopping process between two neighboring lattice sites used to set the effective mass as t = ̵ h 2 (2ma 2 ). The spin-flip hopping amplitudeα is related to the SOI parameter byα = α 2a. The spin-orbit energy E so is given by E so =α 2 t [1,2]. We use this Hamiltonian to obtain the spectra and wavefunction in Fig. (2) in the main text, and to determine the position of the gap closing of the ABSBs.

Finite rectangular geometry
In the finite geometry we assume the system to be finite in both x and y directions and of size L x ×L y = (N x −1)(N y − 1)a 2 . The indices (n, m) label sites in the two-dimensional square lattice with n ∈ {1, . . . , N x } and m ∈ {1, . . . , N y }. The total HamiltonianH =H 1 +H1 +H Γ +H Γ ′ +H D in the finite geometry is then given bȳ The operator c † τ σnm creates an electron with spin projection σ in the layer τ at the lattice site (n, m). Note thatH ′ in the semi-infinite geometry can be obtained fromH by applying the Fourier transformation in the x direction. We use this Hamiltonian to obtain the spectra and MBS wavefunction in Figs. (4) and (5) in the main text, as well as to identify the different topological phases.

APPENDIX B: HYBRIDIZATION OF TRANSVERSE SUBBANDS
As discussed in the main text the 'orbital content' of the ABSBs depends on the momentum k x and changes as a function of W . When W is increased more transverse subbands of the normal region move into the superconducting gap and as a consequence more ABSBs are formed. Generally, the ABSBs in the tunnel coupled Josephson bijunction are gapped. In the semi-infinite geometry, the gap only closes at k x = 0 (the momentum along the junction) for certain combinations of k so W and Γ ∆. Examining the wavefunction of the ABS at a fixed momentum reveals which transverse subband is dominating at that point. In order to better understand the topological phase transitions of several transverse subbands, we consider the lowest ABSB at k x = 0, Γ ∆ = 1.5, and when W is varied from W = 26 to W = 51 [see Fig. 1]. For W = 26 the dominant contribution is coming from the second transverse subband. As W is increased, the gap closes at W = 31 and then reopens as W is increased further. This corresponds to the topological phase transition of the second transverse subband. Since the width of the normal region is increased, a new ABSB forms and the second lowest ABSB then has a dominant contribution from the third transverse subband. At some point the second and third transverse subband hybridize, such that the dominant contribution of the lowest ABSB is switched (at W = 40) and the situation is now opposite. Increasing W further, the gap of the lowest ABSB then closes at W = 50, which corresponds to the topological phase transition of the third transverse subband. While for the lowest ABSB the main contribution at k x = 0 is coming from higher transverse subbands as W is increased, the contribution of the lower transverse subbands appears at finite momenta [see Fig. 2]. The local minima at these finite momenta, where the band has dominant contributions from lower subbands, determines the localization lengths of potential bound states [3]. For W = 41, the n = 1 and the n = 2 subbands are in the topological phase [see Fig. 1], and as the gaps around k x = 0.17 (n = 1) and around k x = 0.12 (n = 2) are quite small [see Fig. 2], the localization length of the MBS become quite large. We note that these localization lengths can be decreased if zigzag-shaped boundaries of the SNS interface are considered [4]. This explains the behavior discussed in Fig. 5 in the main text.

APPENDIX C: ABS WAVEFUNCTION AND TOPOLOGICAL PHASE TRANSITION CRITERION
In this section we give a detailed derivation of the ABS wavefunction present in the coupled Josephson junction setup. As explained in the main text, we assume the system to be translationally invariant along the x direction (along the junction), and finite in the y direction. We then focus on the subgap states at k x = 0 with energy E < ∆ and the time-reversal invariant case φ = π. The BdG equations for this problem then read whereΓ(y) = Γθ(W 2 − y ), and ∆(y) = sgn(y)∆θ( y − W 2) and without loss of generality we assume Γ > 0 and ∆ > 0, as discussed in the main text. The wavefunctionΦ E (y) is given byΦ , v E,1↑ (y)]. Since these two sectors are decoupled the spectrum is two-fold degenerate, and it is enough to focus on one sector in order to derive the zero-energy condition. Without loss of generality, we solve the BdG equations for Φ (1) E (y) and we drop from now on the superscript. First, we perform a unitary transformation U = e iπσ 2 4 , which maps σ 1 → σ 3 and σ 3 → −σ 1 . With these preparations the reduced BdG equations read which are then solved in the normal region y < W 2, in the left superconducting region y < −W 2, and in the right superconducting region y > W 2. The general solution in the normal region is given by where λ ± = 1 ± 1 + Γ 2 Analogously one finds the wavefunction in the right superconducting region, which decays for y → ∞, The total ABS wavefunction then reads Φ(y) = Φ (L) (y)θ(−y − W 2) + Φ (N ) (y)θ(W 2 − y ) + Φ (R) (y)θ(y − W 2). One has to impose the boundary conditions at the two interfaces, which can be recast in matrix form as where a = (a 1 , . . . , a 16 ) is the vector of linear coefficients. Non-trivial solutions exist only if the matrix M is singular,i.e. det(M ) = 0. The solution of this equation then yields the condition on the system parameters (k so W , Γ, ∆) under which a zero-energy ABS exists. This is the topological phase transition criterion we are looking for. As can be seen above, the presence of the square root factors makes it impossible to find an analytical closed form expression. However, in the strong SOI limit, i.e. Γ, ∆ ≪ E so , these square roots can be expanded and thus the matrix M and the equation Eq. (12) are expanded in the strong SOI limit, and can be brought into form In doing so, we find Thus for Γ > 0, Eq. (13) can only be fulfilled if One can see that the system only has zero-energy bound states when the cotangent is positive, since throughout the derivation it was assumed that Γ > 0.

APPENDIX C: EFFECT OF BREAKING TIME-REVERSAL SYMMETRY
In the main text we presented the results for the TRI topological superconducting phase and briefly mentioned the effect of breaking time-reversal symmetry (TRS). In this section we give a more detailed discussion of the TRS broken phase.
In our system, TRS can be broken by either deviating from a π-junction, i.e φ = π + δφ, or by applying a Zeeman field. The TRS breaking term coming from the deviation from a π-junction, acts in the same way on all states independently of their spin-structure. This results in trivially gapping out the MBS, and the system is always in a trivial phase.
We take into account the effect of applying a Zeeman field along the junction (the x-direction), by adding the term where∆ Z (y) = ∆ Z θ(W 2 − y ), and we assume without loss of generality that ∆ Z ≥ 0. We find that the Zeeman field leads to a closing of the gap in the ABSBs at k x = 0 for some critical fields. Starting with ∆ Z = 0 for some fixed width W , one starts in the TRI topological phases discussed in the main text, which then become trivial under the application of the Zeeman field, since the Kramers pairs of MBS are gapped out. Increasing ∆ Z furhter until the ABSB gap closes and reopens the first time, the system enters a topological phase with one MBS at each end of the junction [see Fig. 3]. After the second gap closing and reopening the system has an even number of MBS, which are not protected against disorder, and indeed we observe that they split as a function of disorder strength, and hence, this region is also trivial.