Voltage-tunable Majorana bound states in time-reversal symmetric bilayer quantum spin Hall hybrid systems

TL;DR

This paper explores how a bilayer quantum spin Hall system coupled with an s-wave superconductor can host Majorana bound states, with tunable topological interfaces created by spatially varying gate voltages, revealing unique transport signatures.

Contribution

It demonstrates that Majorana bound states can be engineered in a homogeneous system via gate voltage control, offering a new method to realize topological superconductivity.

Findings

Majorana bound states appear at domain walls of pairing types.

Spatially dependent gate voltage can induce topological interfaces.

Distinct zero-bias anomalies signal Majorana states.

Abstract

We investigate hybrid structures based on a bilayer quantum spin Hall system in proximity to an s-wave superconductor as a platform to mimic time-reversal symmetric topological superconductors. In this bilayer setup, the induced pairing can be of intra- or inter-layer type, and domain walls of those different types of pairing potentials host Kramers partners (time-reversal conjugate pairs) of Majorana bound states. Interestingly, we discover that such topological interfaces providing Majorana bound states can also be achieved in an otherwise homogeneous system by a spatially dependent inter-layer gate voltage. This gate voltage causes the relative electron densities of the two layers to vary accordingly which suppresses the inter-layer pairing in regions with strong gate voltage. We identify particular transport signatures (zero-bias anomalies) in a five-terminal setup that are uniquely related to the presence of Kramers pairs of Majorana bound states.

Tables

Table 1. Table 1: List of allowed (a) and forbidden (f) scattering processes for a helical bilayer hybrid system assuming an incoming electron from the lhs in TI1. The first column defines the process, the second column the corresponding scattering coefficient ( SC ), and the last column the helicity constraint ( HC ) of this process.

scattering process	SC	HC
electron reflection in TI1	$r_{e1}$	f
electron reflection in TI2	$r_{e2}$	f
electron transmission in TI1 (EC)	$t_{e1}$	a
electron transmission in TI2 (CEC)	$t_{e2}$	a
Andreev reflection in TI1 (LAR)	$r_{h1}$	a
Andreev reflection in TI2 (CAR)	$r_{h2}$	a
hole transmission in TI1	$t_{h1}$	f
hole transmission in TI2	$t_{h2}$	f

Equations

$$H_k=\left(\begin{array}{cc}{H}_k^{\text{BL}}& \hat{\Delta }\\ {\hat{\Delta }}^{†}& -\mathcal{T}{H}_k^{\text{BL}}{\mathcal{T}}^{-1}\end{array}\right)=\left(\begin{array}{cc}{H}_k^{\text{BL}}& \hat{\Delta }\\ {\hat{\Delta }}^{†}& -H_k^{\text{BL}}\end{array}\right),$$

$$H_k^{\text{BL}}=\left(\begin{array}{cc}{H}_k^{\text{BHZ}}-{\mu }_1& h_k^t\\ h_k^t& H_k^{\text{BHZ}}-{\mu }_2\end{array}\right),$$

$$H_k^{\text{BHZ}}=\left(\begin{array}{cc}{h}_k& 0\\ 0& h_{-k}^{\ast }\end{array}\right),$$

$$h_k^{t\uparrow }=\frac{1}{2}\left(\begin{array}{cc}{t}_E& \alpha (\sin (k_x)+i\sin (k_y))\\ \alpha (\sin (k_x)-i\sin (k_y))& t_H\end{array}\right),$$

$$\hat{\Delta }=\left(\begin{array}{cc}{\Delta }_1& {\Delta }_c\\ {\Delta }_c& {\Delta }_2\end{array}\right),$$

$$\hat{H}=\frac{1}{2}\int dx\psi (x{)}^{†}\mathcal{H}\psi (x)\text{,}$$

$$\mathcal{H}=\hslash v_{\text{F}}\hat{k}({\zeta }_0{\tau }_0{s}_3)-\frac{{\mu }_1}{2}{\zeta }_3\left(1+{\tau }_3\right)s_0-\frac{{\mu }_2}{2}{\zeta }_3\left(1-{\tau }_3\right)s_0$$

$$-\frac{{\Delta }_1}{2}{\zeta }_2\left(1+{\tau }_3\right)s_2-\frac{{\Delta }_2}{2}{\zeta }_2\left(1-{\tau }_3\right)s_2-{\Delta }_c{\zeta }_2{\tau }_1{s}_2\text{,}$$

$$E$$

$$\Delta E(k=0)={\Delta }_1+{\Delta }_2-\sqrt{({\Delta }_1-{\Delta }_2{)}^2+4{\Delta }_c^2}\text{.}$$

$$\mathbb{N}=\text{sgn}({\Delta }_{+})\text{sgn}({\Delta }_{-})=\text{sgn}(\Delta +{\Delta }_c)\text{sgn}(\Delta -{\Delta }_c)\text{,}$$

$$\displaystyle E=\frac{\pm 1}{2}\Bigl{[}4(\Delta^{2}+\Delta_{c}^{2})+4k^{2}+V^{2}$$

$$\displaystyle\pm 4\sqrt{4\Delta^{2}\Delta_{c}^{2}+V^{2}(\Delta_{c}^{2}+k^{2})}\textit{ }\Bigr{]}^{1/2}\textit{,}$$

$$k_{E=0}$$

$$V_c=\pm 2\sqrt{{\Delta }_c^2-{\Delta }^2}.$$

$${\tilde{\Delta }}_{+}=\frac{1}{2}\sqrt{V^2+4{\Delta }^2}+{\Delta }_c\text{,}$$

$${\tilde{\Delta }}_{-}=\frac{1}{2}\sqrt{V^2+4{\Delta }^2}-{\Delta }_c\text{,}$$

$$\mathbb{N}=\text{sgn}({\tilde{\Delta }}_{-})=\text{sgn}(\frac{1}{2}\sqrt{V^2+4{\Delta }^2}-{\Delta }_c)\text{.}$$

$${\psi }_{\text{I}}(l_1)$$

$${\psi }_{\text{II}}(0)$$

$${\psi }_{\text{III}}(l_2)$$

$${\psi }_{\text{I}}(x)$$

$${\psi }_{\text{II}}(x)$$

$${\psi }_{\text{III}}(x)$$

$${\psi }_{\text{IV}}(x)$$

$$\left(\begin{array}{c}{I}_1\\ I_2\\ I_3\\ I_4\end{array}\right)=\left[\begin{array}{cccc}{G}_1^{\prime }& G_{\text{CAR}}^{12}& -G_{\text{EC}}^{13}& -G_{\text{CEC}}^{14}\\ G_{\text{CAR}}^{21}& G_2^{\prime }& -G_{\text{CEC}}^{23}& -G_{\text{EC}}^{24}\\ -G_{\text{EC}}^{31}& -G_{\text{CEC}}^{32}& G_3^{\prime }& G_{\text{CAR}}^{34}\\ -G_{\text{CEC}}^{41}& -G_{\text{EC}}^{42}& G_{\text{CAR}}^{43}& G_4^{\prime }\end{array}\right]\left(\begin{array}{c}{V}_1\\ V_2\\ V_3\\ V_4\end{array}\right)\text{,}$$

$$G_1^{\prime }=$$

$$G_2^{\prime }=$$

$$G_3^{\prime }=$$

$$G_4^{\prime }=$$

$$\frac{\mathbf{dI}}{\mathbf{dV}}=\frac{e^2}{h}$$

$$\left[\begin{array}{cccc}{G}_1(e{V}_1)& R_{\text{CAR}}^{12}(e{V}_2)& -T_{\text{EC}}^{13}(e{V}_3)& -T_{\text{CEC}}^{14}(e{V}_4)\\ R_{\text{CAR}}^{21}(e{V}_1)& G_2(e{V}_2)& -T_{\text{CEC}}^{23}(e{V}_3)& -T_{\text{EC}}^{24}(e{V}_4)\\ -T_{\text{EC}}^{31}(e{V}_1)& -T_{\text{CEC}}^{32}(e{V}_2)& G_3(e{V}_3)& R_{\text{CAR}}^{34}(e{V}_4)\\ -T_{\text{CEC}}^{41}(e{V}_1)& -T_{\text{EC}}^{42}(e{V}_2)& R_{\text{CAR}}^{43}(e{V}_3)& G_4(e{V}_4)\end{array}\right]$$

$$G_1=$$

$$G_2=$$

$$G_3=$$

$$G_4=$$

$$R_{\text{LAR}}^{\text{E}=0}=$$

$$R_{\text{CAR}}^{\text{E}=0}=$$

Full text

I Introduction

One of the fundamental prerequisites of topological quantum computing dassarma (TQC) is based to the hardware exhibiting quasi-particles with non-Abelian exchange statistics. The primary example of such non-Abelian quasi-particles is provided by Majorana bound states (MBS). They are zero energy excitations with anyonic braiding properties. There is a plethora of theoretical proposals Read2000 ; kitaev2001 ; fu2008 ; fu2009 ; dassarma2010 ; oppen2010 ; alicea2010 ; alicea2012 ; beenakker2013 , predicting favorable conditions for hosting MBS. Moreover, experimental evidence of MBS has been reported in several systems kouwen2012 ; Das2012 ; marcus2016 ; marcus22016 ; Zhang2018 .

In particular, superconducting hybrid structures have been proposed to host MBS. Those hybrids can be realized in various ways, for instance, by coupling topological insulators (TIs) Kane2009b ; beenakker2009 ; trauz2014 ; trauz22014 ; recher2015 ; keidel2018 , semiconducting nanowires dassarma2010 ; oppen2010 ; lossklino2012 ; lee2012 ; aguado2012 ; lossklino2013 ; simon2013 , or chains of adatoms bernevig2013 ; loss2013 ; mcdonald2014 ; yazdani2014 ; franke2015 ; meyer2016 to an ordinary $s$-wave superconductor (SC) via the proximity effect. For our proposal, it is important to be able to attain a decent proximity effect between a quantum spin Hall (QSH) insulator – a two-dimensional (2D) TI – and an $s$-wave SC. Fortunately, this task has already been achieved in the laboratory in QSH insulators based on HgTe in proximity to Al or Nb molenkamp2014 ; kouwen2015 ; yacoby2017 ; molenkamp2017 ; Deacon2017 .

Inspired by the work of Klinovaja, Yacoby, and Loss Klino2014 , we investigate bilayer QSH systems in proximity to an ordinary $s$-wave SC (see Fig. 1 for a schematic). In this setup, we assume competing SC pairings, namely the direct pairing within the individual layers and the crossed pairing between the layers. The latter pairing can be interpreted as a Cooper pair splitting between the layers recher2002 ; fisher2002 ; Sato2010 ; recher2013 . Importantly, a significant inter-layer pairing as compared to the intra-layer pairing can only be realized in the presence of interactions. This requirement can be viewed as a consequence of a deeper theoretical concept, that interactions are necessary for the development of time-reversal symmetric topological superconductivity in low spatial dimensions flensberg2016a ; flensberg2016b ; haim2018 .

It has been shown before Klino2014 that Kramers pairs of MBS are located at the boundaries between two regions in space where the two different pairings (direct and crossed) dominate, respectively. Remarkably, we demonstrate below how an additional inter-layer gate voltage can provide an experimental knob for locally tuning the domain boundary of the superconducting pairing potentials. This approach works because the inter-layer gate voltage effectively suppresses the crossed pairing.

We present analytical results based on an effective model of helical edge states with proximity induced pairing, and corroborate our findings by numerical data on a microscopic lattice model of the 2D bilayer system. We start by investigating the spectral properties of the hybrid system. Based on the two complementary models, we examine the emergence of MBS at a spatial interface between two inequivalent Bogoliubov-de-Gennes (BdG) band structures.

Subsequently, we consider a five-terminal setup which allows us to distinguish all transport processes (electron reflection, Andreev reflection, electron cotunneling, and crossed Andreev reflection) by multi-terminal transport characteristics. Calculating the conductance matrix following the approaches by Blonder, Tinkham, and Klapwijk (BTK) BTK as well as Büttiker buettiker , we find – in the presence of MBS due to a domain wall in the superconducting pairing – dips (in the local conductance) and peaks (in the nonlocal conductance). By tuning the inter-layer gate voltage, this picture can even be inverted, i.e. we can find peaks for the local and dips for the nonlocal conductance. We carefully explain below to which extent the zero-bias features in the conductance are indeed related to the emergence of Kramers pairs of MBS.

The article is organized as follows: In Sec. II, we present the full model in a tight-binding and effective low energy Hamiltonian description. In Sec. III, we examine the spectral properties of this model, where we show two different gap closing mechanisms. In Sec. IV, we explain the appearance of MBS on a spatial interface between two different band gaps in the BdG picture. In Sec. V, we present zero bias anomalies in local and non-local differential conductances. Finally, we conclude in Sec. VI. Some technical details are moved to App. A and App. B.

II Model

We study a hybrid quantum system consisting of two identical HgTe quantum wellsrecher2012 in the QSH regime coupled to an $s$-wave superconductor (see Fig. 1). Additionally, we assume that an inter-layer gate voltage allows us to tune the relative electron densities between the layers. This is a difficult task in experiments because the SC screens the influence of the gates. Hence, the geometry of SC and charging gates have to be chosen in a clever way. We compare our effective low energy results to exact numerical data on a microscopic bulk lattice model of the setup which, in reciprocal space, is described by the Bogoliubov de Gennes (BdG) Hamiltonian

[TABLE]

where the matrix structure is in Nambu space, $\mathcal{T}$ denotes the time reversal operator, and the second equality follows from time-reversal symmetry (TRS). The bilayer Hamiltonian $H_{k}^{\text{BL}}$ describing two identical, tunnel-coupled HgTe quantum wells has the following form in the layer pseudo-spin space:

[TABLE]

where the inter-layer potentials $\mu_{1,2}$ are crucial for some of the physics discussed in this work. In Eq. (2), $H_{k}^{\text{BHZ}}$ is assumed to be the Bloch Hamiltonian of the standard Bernevig-Hughes-Zhang modelbhz2006 on a square lattice with lattice constant $a_{\text{x,y}}=a=1$, i.e.

[TABLE]

where one of the time-reversal conjugate (Kramers) blocks is given by $h_{k}=A(\sin(k_{x})\sigma_{x}-\sin(k_{y})\sigma_{y})-(B\sigma_{z}+D\sigma_{0})(4-2\cos(k_{x})-2\cos(k_{y}))+m\sigma_{z}$, and $\sigma_{i}$ are the standard Pauli matrices in orbital space. In the basis $(|E,\uparrow>,|H,\uparrow>,|E,\downarrow>|H,\downarrow>)$ the tunnelling Hamiltonian $h_{k}^{t}$ is diagonal in spin space. Its matrix structure in orbital space is for spin up given by

[TABLE]

and from TRS we have $\left(h_{k}^{t\uparrow}\right)^{*}=h_{-k}^{t\downarrow}$. Finally, the pairing matrix $\hat{\Delta}$ has the following structure in the layer pseudo-spin space

[TABLE]

where $\Delta_{c}$ denotes the crossed Andreev (inter-layer) pairingrecher2002 ; fisher2002 ; Sato2010 ; feinberg2000 ; recher2013 ; braggio2019 from Cooper pairs delocalized over the two layers while $\Delta_{1}$ and $\Delta_{2}$ denote (intra-layer) pairing within the layers $1$ and $2$, respectively.

If the superconducting pairings $\Delta_{1,2,\text{c}}$ as well as the inter-layer voltage $V$ are much smaller than the bulk insulating gap of the quantum wells, we expect an effective edge description along the lines of Ref. Klino2014, to be a good approximation. There, the insulating bulk of the setup is neglected, thus only considering superconducting pairings in the helical edge states of the QSH bilayer system. For simplicity, we neglect direct tunneling of electrons between the layers, which is reasonable for $\alpha/a<\Delta_{1,2,\text{c}}$, and because the tunnelling decays exponentially with the distance of the two samples.

For the effective one-dimensional edge description along the $x$-axis, the edge modes are linearized around the Fermi-momentum $k_{F}$, and the fermion operators are expressed in terms of slowly varying left and right moving fields $\psi(x)=(R_{1}(x),L_{1}(x),R_{2}(x),L_{2}(x)$, $R_{1}^{\dagger}(x),L_{1}^{\dagger}(x),R_{2}^{\dagger}(x),L_{2}^{\dagger}(x))^{T}$. The resulting non interacting Hamiltonian can be written as

[TABLE]

with the Hamiltonian density

[TABLE]

where the Pauli matrices $s_{i}$ ($\zeta_{i}$) [$\tau_{i}$] act on spin (electron-hole) [upper-lower QSH system] space and $\hat{k}=-i\hbar\partial_{x}$ is the momentum operator in real space representation. By the symmetric choice of $\mu_{1,2}=\pm V/2$, we parametrize the appearance of the chemical potentials in Eq. (II) with a inter-layer gate voltage $V$ corresponding to the chemical potential difference $\mu_{1}-\mu_{2}$. We allow the superconducting pairing potentials $\Delta_{1,2}(x)$, $\Delta_{c}(x)$, and the chemical potential offsets $\mu_{1,2}(x)$ to depend on the spatial coordinate $x$. In the following, we put $\hbar=v_{F}=1$.

Interestingly, our model refers to a superconducting hybrid system that preserves TRS. Hence, it can be classified by a topological quantum number $\mathbb{N}$ that is completely determined by the Fermi-surface properties, more precisely by the sign of the order parameter at the Fermi-surfaces in 1D zhang2010 . We employ this type of classification in the next section.

The characterization and detection of MBS are the main tasks of this article. It turns out that transport properties of hybrid bilayer structures are particularly suitable for the identification of MBS. Hence, we investigate below multi-terminal setups in which we assume that (i) the two layers of the bilayer QSH system are coupled to separate electron reservoirs and (ii) all superconducting regions are grounded. In Fig. 2, we schematically illustrate two setups of that kind that will be further analyzed below. In particular, both figures illustrate N-S-S’-N junctions in which the N regions are characterized by gapless helical edge states (with the same helicity in the two layers, for concreteness) and the S/S’ regions are superconducting regions in which either different pairing amplitudes dominate or different inter-layer gate voltages are assumed.

III Spectroscopic properties

In this section, we discuss the spectroscopic features of our model. We identify BdG band inversions that can be generated by a variation of $\Delta_{c}$ with respect to $\Delta_{1,2}$ (at fixed $\mu_{1,2}$) or by a variation of $V=\mu_{1}-\mu_{2}$ with respect to $\Delta_{c}$ (at fixed $\Delta_{c}>\Delta_{1,2}$). Interestingly, Majorana bound states emerge in hybrid structures of QSH bilayers at interfaces between regions in space with different band orderings. Before we analyze such interfaces and the resulting bound states, we discuss different processes of tuning the band ordering of the BdG excitation spectrum in our system. The general idea is to identify possibilities of gap closings and reopenings as a function of adjustable parameters of the model. Some of the possibilities have been realized before (presented in Sec. III.1)Klino2014 , others are not yet mentioned in the literature (presented in Sec. III.2).

III.1 Gap closing by tuning pairing potentials

First, we consider the effective edge Hamiltonian of Eq. (II) under the choice $\mu_{1,2}=0$, i.e. in the absence of an inter-layer voltage. We look at the translational-invariant case in which the operator $\hat{k}$ can be replaced by the wave vector $k$ in Eq. (II). Furthermore, we assume (for simplicity) that all spatially homogeneous pairing terms are purely real and positive. Under this choice of parameters, the energy eigenvalues of the BdG Hamiltonian are given by

[TABLE]

We note that each energy level has a twofold (Kramers) degeneracy because of TRS. To identify the points in parameter space, where the band ordering is changing, we first consider the gap closing conditions of the spectrum. In this model, the gap closing happens at $k=0$. Therefore, it is sufficient to look at the energy difference $\Delta E(k=0)$ of the lowest ”conduction” band and the highest energy state of the ”valence” band, i.e.

[TABLE]

It is easy to see that the gap closes at $\Delta_{c}^{2}=\Delta_{1}\Delta_{2}$. The dispersion relation for different choices of system parameters is shown by the thick blue and green lines in Fig. 3. Evidently, a band inversion can be generated by a relative change of $\Delta_{c}$ with respect to $\Delta_{1,2}$.

The TRS ensures the supercondcuting pairings to be real and without loss of generality, we assume them to be positiv. By taking $\Delta_{1}=\Delta_{2}=\Delta$, we can identify two relevant pairing potentials at the Fermi-surface, namely $\Delta_{+}=\Delta+\Delta_{c}$ and $\Delta_{-}=\Delta-\Delta_{c}$. Accordingly, we are able to calculate a topological invariant for this superconductor that preserves TRS zhang2010 , i.e.

[TABLE]

where sgn($\Delta_{+})$=sgn($\Delta+\Delta_{c})$ is always greater than zero. Hence, for $\Delta_{c}>\Delta$, the invariant multiplies to $\mathbb{N}=-1$. In such hybrid structures, a Kramers pair of Majorana bound states appears at the boundary of two regions in space with a different topological invariant $\mathbb{N}$ as we demonstrate below.

We have corroborated our results based on the effective edge Hamiltonian of Eq. (II) by numerical simulations on the full 2D lattice model specified in Eq. (1) accounting for possible effects of the finite bulk energy gap on the edge physics. Our results fully confirm the gap closing and gap reopening behaviour from the competition of $\Delta_{c}$ and $\Delta_{1},\Delta_{2}$. The numerical simulations are shown by the transparent lines in Fig. 3.

As one of the main results of our present work, we show next that an inter-layer voltage $V$ allows us to switch between the different SC pairings thus enabling the flexible tunability of domain walls and MBS bound to them.

III.2 Gap closing by tuning inter-layer voltage

In this section, we demonstrate that a gap closing can also be driven by a voltage, where we use the variation of the relative chemical potential parameters ($V=\mu_{1}-\mu_{2}$, where $\mu_{1,2}=\pm V/2$) of the two layers. On the basis of a spatially homogeneous version of the Hamiltonian stated in Eq. (II) with the simplification $\Delta_{1}=\Delta_{2}\equiv\Delta$, we obtain the dispersion relation

[TABLE]

where $V=\mu_{1}-\mu_{2}$. The gap closing condition can be derived by solving Eq. (III.2) for $k$ at zero energy, resulting in the solutions

[TABLE]

To find the true gap closing condition, we have to identify (within this set of possible solutions) the ones that are purely real. Proper solutions are found for two choices of inter-layer voltage, i.e.

[TABLE]

We coin the values of $V$ that fulfil Eq. (13) the critical $V_{\text{c}}$. They correspond to gap closings at wave vectors equal to zero. The choice of a non symmetric inter-layer voltage leads to a closing at finite $k$.

The dispersion relation for different values of the transition parameters $V$ is shown in Fig. 4, where we assume $\Delta_{c}^{2}>\Delta^{2}$. Evidently, a band inversion also happens under this choice of the variation of parameters. Again, we verify our results by numerical calculations on the full 2D lattice model given by Eq. (1). The transparent lines in Fig. 4 complement with numerical data from microscopic simulations of this voltage driven transition the low energy model for various values of the inter-layer voltage $V$ and confirm the voltage induced gap closing at the transition to the trivial superconducting phase with increasing inter-layer voltage $V$.

In the presence of $V$, the BdG dispersion relation allows us to define two effective pairings at the Fermi surface, i.e.

[TABLE]

where $\tilde{\Delta}_{+}$ is greater than zero, because $V$ is chosen to be real. Hence, the topological invariant simplifies to

[TABLE]

Evidently, the sign of the invariant changes at $V=V_{c}$. By interfacing domains with $V<V_{c}$ and $V>V_{c}$, we again find Kramers pairs of MBS that are exponentially localized at the domain walls with spatially oscillating tails. We elaborate on these bound states in the following section.

IV Majorana bound states

In this section, we discuss the emergence of Kramers pairs of MBS at the interface of the two S-S’ junctions in Fig. 2 (in the two layers). These bound states appear at junctions between sectors with different band orderings. As described in the previous section, there are two possibilities to achieve this task: (i) by a variation of $\Delta_{c}$ with respect to $\Delta_{1,2}$ and (ii) by a variation of $V=\mu_{1}-\mu_{2}$ at fixed $\Delta_{c}>\Delta_{1,2}$. We present results for both cases below. Technically, we look at a scattering problem of a N-S-S’-N junction in both cases. Since our system under consideration is time-reversal symmetric, we always obtain Kramers pairs of MBS. However, in our formalism, we employ scattering theory to identify a fingerprint of a MBS in a given scattering state. This approach typically yields one of the two Kramers partners for a particular choice of the scattering state. To identify the other Kramers partner, we need to look at the time-reversed scattering state in which direction of motion and spin are reversed. The corresponding scattering states are specified in the next section and App. A. The resulting bound states are then found by solving the secular equation that derives from the continuity conditions of the spinors and identifying solutions at zero excitation energy.

IV.1 MBS by variation of crossed pairing

The N-S-S’-N junction under consideration is naturally divided into four regions in space which we name I-IV (from left to right). At each interface a continuity condition of the wave function has to be fulfilled because our Hamiltonian describes a set of first-order differential equations. On the basis of the coordinate system defined in Fig. 2 (a), we label the spinor $\psi(x)$ accordingly $\psi_{\text{I}}(x)-\psi_{\text{IV}}(x)$ and take into account three different continuity conditions, i.e.

[TABLE]

In region I and IV, we consider free propagating helical edge states in the bilayer setup in the absence of superconducting pairing. In region II, $\Delta_{1}=\Delta_{2}=\Delta$ is finite but $\Delta_{c}=0$. In region III, we tune the band inversion by $\Delta_{c}$ with respect to $\Delta_{1}=\Delta_{2}=\Delta$. All calculations are done for a scattering state with an incoming electron approaching the junction from region I in TI1. Similar bound states can be determined with other choices of scattering states. For better illustration, we only show the bound states in region II and III. In Fig. 5, we plot the absolute square of the wave function of a solution to the secular equation at zero energy for different choices of $\Delta_{c}$ in region III with respect to finite $\Delta_{1}=\Delta_{2}=\Delta$. As expected, the MBS show an exponential decay with a localization length inversely proportional to the induced superconducting gaps ($\lambda_{\text{MBS}}\propto 1/\Delta E(k=0)$) in both regions, respectively. Realistic values for a superconducting localization length are discussed in the Conclusion.

IV.2 MBS by variation of inter-layer voltage

A similar analysis can be done for the setup shown in Fig. 2 (b). Again, we divide the N-S-S’-N junction into four regions in space labeled I-IV. Similar to the previous section, we investigate a scattering event with an incoming electron approaching the junction from region I in TI1. Now, the superconducting pairing potentials are chosen homogeneously but a step-like variation of the inter-layer voltage discriminates region II from region III. For simplicity, we put $V=0$ in all regions in space and just vary $V$ in region III. The resulting bound states are plotted in Fig. 6. They emerge for $V>V_{\text{c}}$ in accordance to the band inversion properties discussed in Sec. III.2. Interestingly, the bound states shown in Fig. 6 slightly differ from their counterparts in Fig. 5. They are not centered precisely at $x=0$ but rather shifted to the right hand-side. Additionally, they develop oscillations (as a function of space) for larger steps in the inter-layer voltage. The oscillations arise only for $V>2\Delta_{c}$ and their strength is independent of the normal pairing term $\Delta$. Since the two edges are connected only by the crossed Andreev pairing, a symmetric shift of $\mu_{1,2}$ closes the gap at $k=0$. For $V>2\Delta_{c}$, the Fermi-momentum is away from $k=0$, which leads to a real (oscillating) part in the momentum (see Eq. 12). This can be seen in analogy to a spin orbit coupling induced gap opening for a one dimensional nanowire.lossklino2012 For more details, we refer to App. B.

V Transport properties

In this section, we aim to identify signatures of MBS that are observable in transport measurements. Similar to the previous section, we look at different types of N-S-S’-N junctions and calculate the corresponding scattering problem. Importantly, helicity puts tight constraints on the allowed scattering processes. For an incoming electron from region I in TI1 we have to take into account eight different scattering processes. In Tab. 1, we list them all and state which ones are allowed/forbidden by helicity. The allowed scattering processes are given by specific names: Electron cotunneling (EC), crossed electron cotunneling (CEC), local Andreev reflection (LAR), and crossed Andreev reflection (CAR). We stress that the reflection process of CAR is due to the one-dimensional system shown in Fig. 2. For the transport picture related to Fig. 1 one can see this process also as transmission of a hole to TI2.

The complete scattering states are composed of four terms (for the four different regions). They take the explicit forms

[TABLE]

where the momenta in the different regions in space are stated in App. A. In Eq. (V), the vectors $\vec{\mathbf{e_{i}}}$ in the two N regions are 8-dimensional Euclidean basis vectors. The vectors $\vec{\mathbf{u}_{i}}$ and $\vec{\mathbf{v}_{i}}$ that appear in the scattering states for the regions S and S’ are too cumbersome to spell them out explicitly. We specify them in App. A. The continuity conditions stated in Eq. (IV.1) have to hold again. This set of equations allows us to calculate all coefficients of the scattering states (Eq. V), in particular, the transmission and reflection coefficients, $t_{i}$ and $r_{i}$, respectively.

The spin helical leads in our setup provide a stringent selectiveness in scattering processes listed in Tab. 1. In a five-terminal setup, in which the left N region and the right N region are each separately contacted by two electron reservoirs in each TI and the regions S and S’ are grounded, we investigate the (local and non-local) differential conductance $dI_{i}/dV_{j}$. The indices $i$ and $j$ run from 1 to 4 (over the numbering of the four electron reservoirs connected to TI1 and TI2). The current $I_{i}$ corresponds to the current that leaves or enters reservoir $i$ and the voltage $V_{j}$ corresponds to a bias applied locally to reservoir $j$.

Let us, for concreteness, specify the allowed scattering processes for an incoming right-moving electron in TI1: (i) Electron cotunnelling (EC) from reservoir 1 (connected to the lhs of TI1) to reservoir 3 (connected to the rhs of TI1). This process is reflected by a finite net current $G_{\text{EC}}^{31}(V_{1}-V_{3})$ between reservoirs 1 and 3; (ii) crossed electron cotunnelling (CEC) from reservoir 1 to reservoir 4 (connected to the rhs of TI2). This process is reflected by a finite net current $G_{\text{CEC}}^{41}(V_{1}-V_{4})$ between reservoirs 1 and 4; (iii) local Andreev reflection (LAR) at reservoir 1 with corresponding current $2G_{\text{LAR}}V_{1}$; and (iii) crossed Andreev reflection (CAR) between reservoirs 1 and reservoir 2 (connected to the lhs of TI2) with corresponding current $G_{\text{CAR}}^{21}(V_{1}+V_{2})$ Falci .

If we continue the thought for any possible incoming particle at any desired reservoir, the full multi-terminal current-voltage relation can be written in compact form as

[TABLE]

where the diagonal terms of the conductance matrix are given by

[TABLE]

Along the lines of BTK BTK , we can relate the non-linear conductance $dI_{i}/dV_{j}$ of our hybrid system (for simplicity, at zero temperature $T=0$) to corresponding reflection and transmission probabilities through

[TABLE]

with $R=|r|^{2}$, $T=|t|^{2}$, and the diagonal terms

[TABLE]

Equation (19) nicely illustrates the beauty of the helical bilayer system. The five-terminal setup allows us to measure individually the transmission/reflection probabilities of each scattering process. This remarkable feature solely stems from the particular constraints due to helicity summarized in Tab. 1. In the following sections, we argue that this particular advantage of our system (compared to non-helical analogues) gives rise to unique transport signatures that can be related to the presence/absence of (Kramers partners of) MBS at the S-S’ interface.

V.1 Transport in the absence of inter-layer voltage

Since MBS are bound states at zero excitation energy $E=0$, we first focus on the corresponding scattering coefficients. In particular, we try to identify zero-bias anomalies in the structure illustrated in Fig. 2 (a). Fortunately, it is possible to derive and display analytical results for all relevant transmission/reflection probabilities. Under the parameter choice $\mu_{1}=\mu_{2}=0$ and $\Delta_{1}=\Delta_{2}=\Delta$ and $l_{1}<0<l_{2}$, we obtain

[TABLE]

with $N=\left(e^{2\Delta l_{1}}+e^{2l_{2}(\Delta-\Delta_{c})}\right)\left(e^{2\Delta l_{1}}+e^{2l_{2}(\Delta+\Delta_{c})}\right)$. Let us first interpret these results before we discuss concrete signatures thereof in the differential conductance. Evidently, for $\Delta_{\text{c}}=0$, the probabilities $R_{\text{CAR}}$ and $T_{\text{CEC}}$ vanish because the two QSH layers are totally disconnected in this case. Then, regions II and III are identical which is the reason why the finite probabilities $R_{\text{LAR}}$ and $T_{\text{EC}}$ follow \tanh^{2}\big{(}(|l_{1}|+l_{2})\Delta\big{)} and 1-\tanh^{2}\big{(}(|l_{1}|+l_{2})\Delta\big{)}, respectively. For increasing superconducting region, LAR is getting more dominant and EC is suppressed.

More interesting physics arises, of course, for finite $\Delta_{\text{c}}$. Remarkably, the length scales $l_{1}$ and $l_{2}$ enter into the transport coefficients in a non-trivial way. By a comparison of the exponentials in the transport coefficients in Eqs. (21)-(24), we identify a new condition for a transport phase transition, given by $\Delta_{c}=(1+\frac{|l_{1}|}{l_{2}})\Delta$. At this point, all transport processes at $E=0$ have the same probability $1/4$. For larger values of $\Delta_{\text{c}}$, we find vanishing LAR, EC, CEC, and a perfect CAR at $E=0$. Therefore, we see that this perfect CAR is in one-to-one relation with the presence of MBS in our system. The resulting local ($dI_{1}/dV_{1}$) and non-local ($dI_{1}/dV_{2}$) differential conductances as a function of biases $V_{1}$ and $V_{2}$, respectively, are shown in Fig. 7.

Evidently, the pronounced zero bias dip (peak) in the local (non-local) differential conductance can serve as a clear signature of the presence of Kramers pairs of MBS in our system. However, it is interesting to realize that for transport properties the condition for zero bias features is slightly altered compared to the spectroscopic properties discussed in the previous section. This peculiarity has been realized in simpler hybrid structures based on single-layer quantum spin Hall systems before Fleckenstein . While for the spectroscopic appearance of zero energy bound states, the criterion $\Delta_{\text{c}}>\Delta$ is relevant (for the choice $\Delta_{1}=\Delta_{2}=\Delta$), the zero bias features in transport only emerge if $\Delta_{\text{c}}>(1+\frac{|l_{1}|}{l_{2}})\Delta$. In the latter case, the different weights of spatial inhomogeneities matter with respect to the length scales $l_{1}$ and $l_{2}$.

V.2 Transport in the presence of inter-layer voltage

For finite inter-layer voltage, our analytical results for the transport probabilities are too lengthy to be explicitly written down here. Therefore, we first graphically display the relevant reflection ($R_{\text{LAR}}$ and $R_{\text{CAR}}$) probabilities in Fig. 8 as a function of the excitation energy before we turn to the presentation of the corresponding differential conductance.

At the critical inter-layer voltage defined by $V_{\text{c}}$ via Eq. (13) we realize a qualitative change in the transmission/reflection properties. For $V<V_{\text{c}}$, at small excitation energies, transport is dominated by CAR because of our choice of pairing potentials $\Delta_{c}=2\Delta$ in Fig. 8. This behavior drastically changes for $V>V_{\text{c}}$. Then, LAR dominates for $E/\Delta\approx 0$. This comes from the fact that the energetic distance between the Fermi-energies of the two layers ($\mu_{1,2}$) is greater than the gap opened by the crossed pairing. The oscillations at larger excitation energies $E/\Delta\geq 1$ are Fabry-Pérot resonances due to quasiparticle interference. The width of the zero bias peaks or dips can be varied by the system size. The larger the superconducting regions, the smaller is the width of the feature at zero excitation energy.

The corresponding differential conductances are shown in Fig. 9.

Interestingly, we can identify a dip in the non-local conductance $dI_{1}/dV_{2}$ in the regime of large inter-layer voltage in which a MBS emerges at the step in $V$. This dip should be put in contrast to the peak that we observe in Fig. 7 where the MBS is tuned by a step-like variation of $\Delta_{c}$. Hence, we observe that zero-bias anomalies (peaks or dips) signal the presence of MBS in the five-terminal setup.

V.3 Length dependence of the MBS features

We have introduced in Sec. V.1 an effective transport phase transition weighted by the length scales $l_{1}$ and $l_{2}$. Now we show that the length scales of region II ($|l_{1}|$) and the length of region III ($l_{2}$) can affect the MBS and its conductance features. Therefore, we take a parameter choice where the MBS and its signatures already exist and vary the length of region III. We show in Fig.  10 the resulting local conductance with the related MBS for three sizes of region III.

Interestingly, we find that with the decrease of region III the dip of the conductance vanishes. Coincidently, the MBS gets shifted to the left interface, but still peaks at $x=0$. The reason for the length dependencies is that the scattering process of a particle at zero energy is not only influenced by the magnitude of the potential but also by the length of the acting potential. The new condition of the emergence of a zero bias anomaly on a N-S-S’-N junction is given by the strength of the order parameter times the corresponding length.

VI Summary and Conclusion

In this article, we have compared the low energy model of a bilayer hybrid system of helical edge states with induced superconducting pairing to numerical data on a microscopic tight-binding model, confirming the validity of the effective model, and the competition between a direct pairing and a crossed pairing superconducting gap. Going beyond previous work on a similar setup, we have introduced an inter-layer voltage that can suppress the crossed pairing driving the system from the crossed pairing phase to the direct pairing phase by merely tuning the voltage. This feature represents a simple experimental knob for tuning the topological phase of the system. With the help of a spatially varying inter-layer voltage, domains with different effective pairings hosting Kramers pairs of MBS are shown to be achievable. Additionally, we have discussed a five-terminal setup for a proposed transport experiment. We have shown that the emergence of the MBS is in accordance with zero bias anomalies in local and non-local differential conductance. The zero bias behaviour is characterized by peaks or dips in the differential conductance. Finally, we have also addressed the question of how the lengths of the different scattering regions can affect transport signatures and the emergence of Kramers pairs of MBS in distinct ways.

Let us conclude with a brief discussion of the feasibility of our proposal. We require that the distance $d$ of the two QSH layers should be smaller than the coherence length $\lambda$ of the superconductor (e.g. for Al $\lambda\approx 1.6\mu$m). Moreover, inter-edge tunnelling should be avoided. In a previous work on bilayer QSH insulators, some of us have shown that inter-layer tunnelling almost vanishes for a barrier thickness of $d>10$nm recher2012 . Thus, there is a window of interlayer distance $d$ that is favorable for our idea.

We thank C. Fleckenstein, F. Keidel, J. Klinovaja, D. Loss, H. W. Schumacher and P. Silvestrov for fruitful discussions. We acknowledge financial support by the DFG (SPP1666, SFB1170 “ToCoTronics”, SFB1143 “Correlated Magnetism: From Frustration to Topology” ID 247310070 , AS 327/5-1, RE 2978/7-1 and the DACH project on “Majorana and parafermions in topological insulators”), the ENB Graduate school on “Topological Insulators”, and the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter (EXC 2147, project-id 39085490).

Appendix A MBS wave functions in the absence of an inter-layer gate voltage

To get more insight into the shape of the MBS, we examine here the structure of the wave vectors in more detail. First, we take into account the setup shown in Fig. 2 (a). Region I implies the free TI edge states of both samples with right and left moving electrons or holes following the dispersion $E=\pm k_{e,h}^{1,2}$ ($\hbar=v_{\text{F}}=1$). The related wave vectors for electrons and holes in region II, from $x=l_{1}$ ($l_{1}<0$) to $x=0$, with energies in the superconducting gap are given by

[TABLE]

We identify an exponential decay into region II with a decay length inversely proportional to the superconducting gap $\lambda_{II}=1/\Delta$ at zero energy.

In region III from $x=0$ to $x=l_{2}$, where the crossed pairing enters the game, we obtain the following wave vectors:

[TABLE]

Again, we identify an exponential decay into region III, inversely proportional to the gap $|\Delta-\Delta_{c}|$, illustrated in Fig. 3. Since we consider the MBS, which has zero excitation energy we find a vanishing real part $\operatorname{Re}\left[k_{\pm}\right]=0$. This implies that there is no oscillating part in the wave function. The imaginary part of the wave vectors at zero energy vanishes for $\Delta_{\text{c}}=\Delta$, which is also an indication for the phase transition. The exponential decaying MBS in Fig. 5 has a perfect symmetric shape for $\Delta_{c}=2\Delta$, which comes from the fact, that in this case the energy gaps in both regions have the same magnitude and conclusively the same decay length.

Assuming a symmetric length of region II and III around $x=0$ ($l_{1}=-l_{2}=-L$) and employing continuity conditions, we find the zero energy wave function, up to a normalization factor:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $A=\frac{2e^{2L(2\Delta+\Delta_{c})}}{\left(e^{4\Delta L}+e^{2\Delta_{c}L}\right)\left(e^{2L(2\Delta+\Delta_{c})}+1\right)}$, $\tilde{L}_{1}=L+x$, $\tilde{L}_{2}=2L-x$ and $\tilde{L}_{3}=3L-x$. The beauty of the scattering wavefunction approach can also be seen in Eqs. (28) and (29), by comparing the process related entries on both sides.

Assuming the incident electron is a right mover coming from TI1, which restricts us to only four allowed scattering states, mentioned in the main text. Those processes are related to the four non vanishing entries of the wavefunction. The resulting MBS arises at the interface at $x=0$. For simplicity, we only display scattering states in region II and region III $(\Psi_{\text{II}}(x)$ and $\Psi_{\text{III}}(x))$. The resulting bound state has not yet the properties of a MBS. To create the MBS we first write the bound state in the compact form

[TABLE]

Next we superpose the bound state Eq. (31) with its charge conjugated partner. Note that the charge conjugation operator must satisfy the condition

[TABLE]

which is fulfilled by $\mathcal{C}=U_{c}\mathcal{K}=\zeta_{1}\mathcal{K}$, where $\mathcal{K}$ denotes the complex conjugation and $U_{c}$ is a unitary transformation. Conclusively, we can construct two MBS by the superpositions

[TABLE]

It turns out, that the charge conjugated wave function $\mathcal{C}\Psi(x)$ is exactly the wave function for the bound state when we assume an incoming hole in TI1. Interestingly, we obtain another MBS which is related to a symmetry, which we call sample symmetry $\mathcal{S}=\tau_{1}$. This symmetry changes the sample space if it acts on the Hamiltonian, so that there are two additional MBSs, namely

[TABLE]

These MBS correspond to bound states if we assume an incoming electron in TI2 with the superposition of its charge conjugated wave function. Those states are orthogonal to each other. Finally, we find with the TRS operator $\mathcal{T}=is_{2}\mathcal{K}$, eight different MBS located at the interfaces. They can be specified as

[TABLE]

Appendix B MBS wave functions in the presence of an inter-layer gate voltage

Here, we display the momenta for the region with the inter-layer gate voltage and examine the wave vectors for zero energy. We use the setup in Fig. 2 (b), where we assume $\Delta_{c}>\Delta$ in region II and an additional inter-layer gate voltage in region III. Again, region I and IV represent the free edge states in TI1 and TI2. The related wave vectors in region II are given by Eq. (26) and were already discussed in App. A. In region III, we vary the inter-layer gate voltage. The resulting momenta in Eq. (V) are given by

[TABLE]

where the inter-layer gate voltage is defined as $V=\mu_{1}-\mu_{2}$. The creation of the scattering state for an incoming rightmoving electron in TI1 is done in the same manner as in App. A. We get more information of the MBS in Fig. 6, by looking at the momenta in region III. Fig. 11 shows the imaginary and real parts of Eq. (36) at zero energy in dependence of $V/\Delta_{c}$ under the assumption that $\Delta_{\text{c}}=2\Delta$.

This result explains, why we find in region III an oscillating part, together with an exponential decay.