Nonequilibrium Electrical, Thermal and Spin Transport in Open Quantum Systems of Topological Superconductors, Semiconductors and Metals

Abstract

We study nonequilibrium transport in various open quantum systems whose systems and leads/baths are made of topological superconductors (TSs), semiconductors, and metals. Using quantum Langevin equations and Green’s function method, we derive exact expressions for steady-state electrical, thermal, and spin current at the junctions between a system and leads. We validate these current expressions by comparing them with the results from direct time-evolution simulations. We then show how an electrical current injected in TS wires divides into two parts carried by single electronic excitations and Cooper pairs. We further show ballistic thermal transport in an open TS wire in the topological phase under temperature or voltage bias. The thermal current values grow significantly near the topological phase transition, where thermal conductance displays a sharp quantized peak as predicted earlier. We relate the quantized thermal conductance to the zero-frequency thermoelectric transmission coefficient of the open TS wire. We also observe a large thermoelectric current near the topological transition of the TS wires. We introduce a differential spin conductance which displays a quantized zero-bias peak at zero temperature for a spinful TS wire in the topological phase. The role of superconducting baths in transport is demonstrated by thoroughly examining the features of zero-temperature differential electrical conductance and thermal conductance in open systems with TS baths. Our new thermoelectric and spin transport findings in various two-terminal geometries are beneficial to the present challenges in probing the emergence of Majorana quasi-particles in experiments.

1 Introduction

The study of quantum transport in superconducting materials has attracted much attention in recent years due to its applicability in detecting intriguing topological phases of superconductors [1, 2]. Several recent experiments based on electrical transport measurements have strongly suggested possible evidence of Majorana bound states (MBSs), which are exotic quasiparticle excitations localized at the edges of one-dimensional (1D) topological superconductors (TSs) [3,4,5,6,7]. While the electrical transport measurements can detect the emergence of protected zero modes of the MBSs, the thermal transport measurements are further useful in directly probing these chargeless Majorana quasiparticles [2, 8, 9]. There has been a surge of interest in quantized thermal conductances for detecting fractionally charged and neutral modes [10,11,12]. The MBSs in experimentally realized TSs made of hybrid nanowires combining spin-orbit coupled semiconductor and superconductor materials have a non-zero spin polarization [13, 14]. Therefore, in such systems’ spin transport is expected to provide valuable information about the systems’ MBSs and topological properties [15,16,17]. The thermal and spin transport measurements in probing the emergence of Majorana fermions can become essential tools in the present scenario when the detection of MBSs through electrical current peak signals is not unambiguous [18,19,20,21,22,23] because such electrical signals can also be produced in these devices due to other than Majorana fermions, such as, other quantum states that are not Majoranas and imperfections in the nanowire [24, 25].

On the theoretical side, while there is a vast number of theoretical studies on electrical transport in systems with TSs [24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41], the thermal [2, 8, 9, 42,43,44,45] and mostly spin transport [46, 47] in these systems are much less explored. One of this paper’s primary goals is to give a detailed and unified description of electrical, thermal, and spin transport in different devices made of TSs. The quantum transport in such devices are generally studied utilizing the generalized scattering theory [8, 9, 48, 49], the Keldysh formalism [28, 30, 35, 44, 50], and the quantum Langevin equations & Green’s function (LEGF) method [25, 33, 41]. These theoretical techniques are mainly employed to calculate electrical currents and differential conductances, which are supposed to show non-trivial behaviour in superconductors’ topological phases. We apply the LEGF method to develop a unified electrical, thermal, and spin transport description in 1D superconducting devices. The LEGF method is useful in the physical understanding of nonequilibrium processes as it provides a nice picture of the baths (generating the bias) as sources of thermal and quantum noise and dissipation. The last feature is also helpful in deriving fluctuations in the current (e.g., current-current correlators) applying this method.

The LEGF method is an open-quantum system formulation within the Heisenberg representation of quantum mechanics to study nonequilibrium quantum transport in mesoscale and nanoscale devices [51,52,53,54,55]. This method which is based on a direct solution of the Heisenberg equations of system and bath variables, has been applied extensively in calculating steady-state [54, 55] and time-dependent [53, 56] electrical, thermal, and optical transport properties in devices consisting of noninteracting baths (e.g., metals, harmonic oscillators) [25, 33, 41, 57,58,59]. An extension of this method for superconducting baths lacks to date. Such an extension is necessary for dealing with the Josephson effect within an open-quantum system framework [40] and various Majorana braiding schemes proposed using TSs [1, 26]. In this paper, we consider hybrid devices of X-Y-Z configuration where X, Y, and Z can be made of TSs, normal metal (N), and spin-orbit coupled semiconductor (SM) in the presence of a magnetic field [40]. Here, the X and Z wires act as leads/baths in thermal equilibrium, and the Y wire is the system through which the transport is happening. Searching the signatures of TS wire leads in quantum transport when these TS wires are in the topological phase is another aspect of our study. To this end, we consider both the Kitaev chain [60] and the experimentally realizable 1D semiconductor-superconductor heterostructures [61, 62] as a TS lead.

It is worth noting that we have already performed a direct time-evolution study of electrical current at both the junctions of various hybrid devices of X-Y-Z configuration in [40]. We have detected a persistent and oscillating electrical current at both junctions of a TS-N-TS device, even in the absence of any phase or voltage or thermal bias when multiple MBSs and/or Andreev bound states (ABSs) within the bulk-gap are present near the junctions. Moreover, the amplitude and period of the oscillating current strongly depend on the middle N wire’s initial conditions indicating the absence of thermalization. Therefore, the presence of such bound states (both MBS and ABS) localized near the junctions prohibits the full hybrid devices from attaining a unique nonequilibrium steady state (NESS) at a long time. It can be noted that the initial-condition dependence of supercurrents in a phase-biased superconducting nanojunction of topologically trivial BCS superconductors has also been investigated in Refs. [63,64,65,66].

Interestingly, it has been demonstrated that the presence of MBSs amplifies the amplitude of zero-bias oscillating currents at the junctions compared to the amplitude of the same generated solely due to ABS. The generalized LEGF method can not be applied to study transport in those devices which do not have unique NESS. On the other hand, thermalization in tandem with a unique NESS can be achieved by tuning the system parameters in TS-N-Z and TS-TS-Z devices where TS is a Kitaev chain, and Z is either an N wire or a Kitaev chain at the topological phase transition point (TP). A unique NESS is reached since the Z lead’s energy spectrum becomes gapless in such situations. There is no longer any bound state from the middle N/TS wire, and the middle wire gets equilibrated with the boundary wire(s). Similarly, we achieve thermalization and unique NESS in TS-SM-Z and TS-TS-Z devices with TS wires made of semiconductor-superconductor heterostructures when the Z wire is either an SM or a TS at TP, and there is no mid-gap state in the spectrum of boundary wires. Hence one can apply generalized LEGF method to TS-N-Z, TS-SM-Z, and TS-TS-Z devices once these systems attain steady-state.

1.1 Overview and New Findings

In Sect. 2, we introduce the Hamiltonian of different models of TS, N, and SM wires and describe their statistical properties when they act a lead/bath. We particularly emphasize a general and detailed description of the TS wires as leads since this is the highlight of our present study. Next, we apply the generalized LEGF method in Sect. 3 for deriving analytical expressions for the steady-state electrical, thermal and spin currents in various devices with TS leads. We also evaluate differential electrical and thermal conductances, which are incredibly convenient in identifying TS leads’ role. This section is divided into two parts, (a) devices with Kitaev chains and (b) those with semiconductor-superconductor heterostructures. In Sect. 4, we employ these expressions of currents and differential conductances to calculate several impressive results. We first validate our steady-state current formulas by comparing them to the long-time currents obtained from the direct time-evolution simulation. We then show how electrical current injected in one junction separates into different parts of charge current carried by single electronic excitations and Cooper pairs inside TS wires of an N-TS-N. We further discuss several interesting features of thermal current and linear-response thermal conductance in the N-TS-N device. The properties of a sharp quantized peak of thermal conductance near the topological phase transition are especially highlighted. We relate the quantized thermal conductance to the zero-frequency thermoelectric transmission coefficient of the open TS wire. We observe a large thermoelectric current near the topological transition of the TS wires, which might be potentially applicable. We then discuss an interesting electrical current asymmetry in a TS-N-N device with spatial asymmetry in tunneling rates. We here introduce a differential spin conductance which displays a quantized zero-bias peak at zero temperature for a spinful TS wire in the topological phase. The role of superconducting leads/baths in transport is shown by thoroughly examining the zero-temperature differential electrical conductances (DECs) and thermal conductances in various devices. We conclude the paper’s central part in Sect. 5 by providing an outlook and problems of interest shortly. Eight appendices are further added to include analytical expressions and our method’s various details for the interested readers.

To summarize, we develop a generalized LEGF formalism which can be applied to study nonequilibrium electrical, thermal, and spin transport in an X-Y-Z device made of superconducting leads as well as metallic leads; this is a technical advancement as the previous use of LEGF was only applicable to the devices with metallic leads. We note that the use of superconducting leads is of present interest as this can give better insight for Majorana detection [36, 38]. On the application side of this generalized LEGF method, apart from verifying some of the already known results using this independent technique, we find several new features. (i) In a TS-TS-N device where TS is a Kitaev chain, we observe a quantized peak in thermal conductance near the TP of the middle TS wire only when the left superconducting lead is also at the TP. This feature of thermal conductance might have a potential application in detecting the topological phase of TS lead via thermal conductance measurement. (ii) We find a sizeable thermoelectric current in an N-TS-N device near the TP of the middle TS wire. Experimental measurements of such a large thermoelectric current or conductance in these devices would be much easier than the relatively small thermal conductance peak near the phase transition. Thus, the thermoelectric current or conductance might be a better probe to detect the TS wires’ topological phase transition experimentally. (iii) We show quantized zero-bias peaks of differential spin conductance in SM-TS-SM devices when the middle TS wire of semiconductor-superconductor heterostructures wire is in the topological phase. Further, we verify that the zero-bias peaks of differential spin conductance are robust against disorder. So, the detection of zero-temperature zero-bias peaks of differential spin conductance via spin tunneling spectroscopy may open up a new avenue to detect elusive Majorana fermions which emerge at the edges of TS wire in the topological phase. In the following sections, we will discuss this generalized LEGF method and its applications in detail.

2 Topological Superconductor Leads and Their Statistical Properties

We consider a hybrid device of X-Y-Z configuration consisting of a finite wire Y whose left and right terminals (ends) are connected to the semi-infinite wire X and Z respectively. We here treat these semi-infinite X and Z wires as leads (baths) and impose canonical or grand-canonical equilibrium for their statistical properties from the beginning before they are connected to the middle Y wire. Hereafter, we choose X, Y, and Z wires to be an N or an SM or a TS wire. We are particularly interested in treating the X and Z wires made of TSs. Below, we introduce the mean-field Hamiltonians of two different 1D TS models and discuss their statistical properties in thermal equilibrium. These models are (a) the Kitaev chain of a spinless p-wave superconductor [60] and (b) the Majorana wire of semiconductor-superconductor heterostructure. The latter one is a spinful TS engineered with a Rasbha spin-orbit coupled semiconductor nanowire proximity coupled to an s-wave superconductor in the presence of a magnetic field along the direction of the wire [61, 62]. We also study the Hamiltonians of N and SM wires as the limiting case of the Kitaev chain and Majorana wire, respectively, and their statistical properties in equilibrium.

We here write the Hamiltonian of the TS wires and those of N and SM wires in a matrix format using a ‘double-fermion’ basis, convenient for our nonequilibrium transport analysis with the quantum Langevin equations. Let us start by writing a most general 1D (mean-field) Hamiltonian [67] of fermions as follows

$$\begin{aligned} \frac{H}{\hbar }=\sum _{l'=1}^L \psi _{l'}^\dagger U_{l'} \psi _{l'} +\sum _{l'=1}^{L-1} (\psi _{l'+1}^\dagger V_{l'} \psi _{l'}+\psi _{l'}^\dagger V_{l'}^\dagger \psi _{l'+1}), \end{aligned}$$

(1)

where $l'=1,\dots ,L$ are the lattice sites along the wire, and $\psi _{l'}$ ($\psi _{l'}^\dagger$) is a column (row) vector of fermion annihilation/creation operators at the $l'$-th lattice cite. This Hamiltonian is defined over an LM dimensional Hilbert space ${\mathbb {H}}=h_1\otimes h_2 \otimes \dots \otimes h_L$ where $h_{l'}$ is an M dimensional local Hilbert space defined at $l'$-th site. Clearly $\psi _{l'}$ is an M dimensional column vector, whereas $U_{l'}$ and $V_{l'}$ are $M \times M$ dimensional matrices. In case of superconductors, $\psi _{l'}$ consists of both electron annihilation and creation operators. The generalized Hamiltonian (1) can accommodate both the Kitaev chain and the Majorana wire for some particular choices of $\psi$, U and V.

2.1 Kitaev Chain

To write the Kitaev chain in the above generalized form, we introduce the position space Nambu spinor $\psi _{l'}=(c_{l'},c_{l'}^\dagger )^T$, which implies that the dimension of local Hilbert space ($h_{l'}$) is two ($M=2$). With the above form of $\psi _{l'}$ and the following choice of $U_{l'}$ and $V_{l'}$,

$$\begin{aligned} U_{l'}= \frac{1}{2}\begin{bmatrix} -\epsilon &{}0\\ 0&{} \epsilon \\ \end{bmatrix}, \, \text {and} ~ ~ V_{l'}= \frac{1}{2}\begin{bmatrix} -\gamma &{} \Delta \\ -\Delta ^* &{} \gamma \\ \end{bmatrix}, \end{aligned}$$

the above Hamiltonian H (1) represents that of the Kitaev chain which is denoted by

$$\begin{aligned} \frac{H_{\mathrm{K}}}{\hbar }=-\gamma \sum _{l'=1}^{L-1}(c^{\dagger }_{l'}c_{l'+1}+c^{\dagger }_{l'+1}c_{l'})-\epsilon \sum _{l'=1}^{L}(c^{\dagger }_{l'}c_{l'}-\frac{1}{2}) - \sum _{l'=1}^{L-1}( \Delta c^{\dagger }_{l'}c^{\dagger }_{l'+1}+ \Delta ^* c_{l'+1}c_{l'}), \nonumber \\ \end{aligned}$$

(2)

where, $\gamma$ is hopping, $\epsilon$ is the on-site energy, and $\Delta$ denotes superconducting pairing potential. The parameters $\gamma ,\epsilon$ and $\Delta$ have dimension of frequency, and we assume them to be real hereafter. In order to study nonequilibrium transport using LEGF, we further introduce the following generalized basis:

$$\begin{aligned} \mathbf{a}\equiv & {} [ a_1, a_2, \dots , a_{2l'-1},a_{2l'},\dots , a_{2L-1}, a_{2L}]^T = [c_1, c^{\dagger }_1, \dots ,c_{l'},c_{l'}^\dagger ,\dots , c_{L}, c^{\dagger }_{L}]^T.\nonumber \\ \end{aligned}$$

(3)

Clearly, $a_{2l'}=a^{\dagger }_{2l'-1}$. In the above basis, $H_{\mathrm{K}}$ can be written in a quadratic form as follows, $H_{\mathrm{K}}=\frac{\hbar }{2}{} \mathbf{a}^{\dagger }{\mathcal {K}}\mathbf{a}=\frac{\hbar }{2}\sum _{l,m}{\mathcal {K}}_{lm} a^{\dagger }_la_m$, where ${\mathcal {K}}$ is an $2L \times 2L$ Hermitian matrix, and $l,m =1,\dots , 2L$ [40, 68]. It can be noted that the index l in $a_l$ (or $a_l^\dagger$) does not represent the actual physical site of the wire. For a given l, one can define a map to the physical site $l'$ of spinless fermions as: $l'=(l+1)/2$ for odd values of l, and $l'=(l/2)$ for even values of l. In the presence of pairing ($\Delta \ne 0$), the superconducting wire undergoes a topological phase transition as $\epsilon$ is tuned across $2\gamma$. The wire is in a topological phase for $\mid \epsilon \mid <2 \mid \gamma \mid$, and the TS wire hosts two spatially-localized MBSs at the opposite ends of the wire for a relatively long wire. The wire transits to a topologically trivial phase (non-topological phase) without the MBSs for $\mid \epsilon \mid >2\mid \gamma \mid$. The topological phase transition near $\mid \epsilon \mid =2 \mid \gamma \mid$ is also accompanied by a bulk-gap closing in its energy dispersion. The superconducting wire has a bulk-gap in its spectrum both in the topologically non-trivial and trivial phases, and the gap vanishes at the topological phase transition around $\mid \epsilon \mid =2 \mid \gamma \mid$. The two phases of the Kitaev chain can be characterized unambiguously by the quantized value of the geometric phase, namely the Pancharatnam–Zak phase, which acts as a topological invariant for such 1D systems [69, 70]. The Pancharatnam–Zak phase’s values are $\pi$ and 0, respectively in the topological and non-topological phases.

The Hamiltonian matrix ${\mathcal {K}}$ can be diagonalized by solving the Hermitian eigenvalue problem

$$\begin{aligned} {\mathcal {K}} U_r= \omega _r U_r, \end{aligned}$$

(4)

where $\omega _r$ and $U_r$ are the r-th eigenvalue and corresponding eigenfunction of the Kitaev Hamiltonian. Since ${\mathcal {K}}$ is a Hermitian matrix in Nambu basis, it satisfies following property [68]

$$\begin{aligned} \Sigma \, {\mathcal {K}} \, \Sigma =-{\mathcal {K}}^{ *}, \end{aligned}$$

(5)

where, $\Sigma = \mathbb {1}_{L} \otimes \sigma _x$ and $\mathbb {1}_{L}$ is an identity matrix of size L. When ${\mathcal {K}}$ is real, all the components of $U_r$ are also real. From Eq. 5, it follows

$$\begin{aligned} {\mathcal {K}} V_r = -\omega _r \, V_r, \end{aligned}$$

where, $V_r=\Sigma U^{*}_r$. Thus, the vector $V_r$ is also be an eigenvector with the eigenvalue $-\omega _r$. We group the eigenvalues of ${\mathcal {K}}$ into pairs ($\pm \,\omega _r$) with $\omega _r > 0$ . The eigenvectors of ${\mathcal {K}}$ obey the completeness relation,

$$\begin{aligned} \sum _{r>0}\left( U_r U_r^{\dagger } + V_r V_r^{ \dagger } \right) =\mathbb {1}_{2L}, \end{aligned}$$

(6)

where, the notation $r>0$ means that the sum is limited to the positive eigenvalues. From the properties of normalized eigenvectors, it follows

$$\begin{aligned} U_r^{\dagger }U_s=\delta _{rs}\,, ~~V_r^{\dagger }V_s=\delta _{rs}\, ,~~ U_r^{\dagger }V_r=0. \end{aligned}$$

(7)

By using the aforesaid properties of eigenvectors and eigenvalues, one can express ${\mathcal {K}}$ in the following form

$$\begin{aligned} {\mathcal {K}}= \sum _{r>0} \omega _r (U_r U^{\dagger }_r- V_r V_r^{ \dagger }). \end{aligned}$$

(8)

Sometimes it is more convenient to express (6) and (7) by the components of $U_r$ and $V_r$. Thus, denoting $U_r$ by

$$\begin{aligned} U_r= ( \phi _r(1), \psi _r(1), \dots , \phi _r(L ), \psi _r(L ) )^T, \end{aligned}$$

one can rewrite (6) as

$$\begin{aligned}&\sum _{r>0}( \phi _r (i') \phi _r^{ *}(j')+\psi _r^{ *}(i') \psi _r(j'))=\delta _{i'j'}\,, \nonumber \\&\sum _{r>0}(\phi _r(i') \psi _r^{ *}(j')+\psi _r^{ *}(i') \phi _r(j'))=0. \end{aligned}$$

(9)

Applying the expressions (8) and (9), $H_{\mathrm{K}}$ (2) can be expressed in a diagonal form:

$$\begin{aligned} H_{\mathrm{K}}= & {} \sum _{r>0} \hbar \omega _r (q_r^{ \dagger }\, q_r-q_r \,q_r^{\dagger }) = \sum _{r>0} \hbar \omega _r (2q_r^{ \dagger }\, q_r-1), \end{aligned}$$

where, the fermionic quasiparticle destruction operators $q_r$ are defined as

$$\begin{aligned} q_r = U_r^\dagger \mathbf{a}= \sum _{i'=1}^{L}( \phi _r^{ *}(i')c_{i'} +\psi _r^{ *}(i') c_{i'}^\dagger ). \end{aligned}$$

(10)

The ground state energy of $H_{\mathrm{K}}$ is $E_g=-\sum _{r>0} \hbar \omega _r$, which corresponds to the quasiparticle vacuum $\mid \varnothing \rangle$. It can be shown that $q_r$ and $q_r^{ \dagger }$ satisfy anticommutation relations (e.g., $\{q_r,q_s^{ \dagger }\}=\delta _{r s}$) that indicate fermionic nature of the quasiparticles. Evidently, this Hamiltonian can be diagonalized in terms of these Bogoliubov quasiparticles, which are linear superpositions of the excitations of negatively charged electrons and positively charged electron holes. These quasiparticle creation operators acting on $\mid \varnothing \rangle$ create many-particle states. Moreover, any second-quantized fermionic operator defined on this Hilbert space automatically takes care of the Pauli exclusion principle. Using (9), one can easily invert (10) to get back

$$\begin{aligned} c_{j'}=\sum _{r>0} ( \phi _r(j') q_r + \psi ^{ *}_r(j') q_r^{\dagger }). \end{aligned}$$

(11)

We here use a semi-infinite Kitaev chain to model the TS leads/baths for X or Z wire. We assume that the wire is in thermal equilibrium at temperature T and chemical potential $\mu$ before connecting it to the middle wire. It is now convenient for the superconducting leads to perform a gauge transformation such that the chemical potential does not explicitly appear in the leads’ thermal density matrix. Under such gauge transformation, the chemical-potential differences instead enter in our calculation through time-dependent phases in the tunnel couplings of the superconducting leads to the middle wire [37]. Therefore, the quasiparticle modes of the superconducting leads/baths in thermal equilibrium satisfy the following relations:

$$\begin{aligned} \langle q_r^{ \dagger } \, q_s \rangle = f(\omega _r, T) \, \delta _{r s}, ~\langle q_r \, q_s^{ \dagger } \rangle = f(-\omega _r, T) \, \delta _{r s}, \end{aligned}$$

(12)

where $\langle .. \rangle$ denotes equilibrium expectation with thermal density matrix. All other expectations like $\langle q_r^{ \dagger }\, q_s^{ \dagger } \rangle$, $\langle q_r \, q_s \rangle$ are always zero. Here, $f(\omega ,T)=1/(\exp [\hbar \omega /k_BT] +1)$ describes the equilibrium distribution of the fermionic quasiparticles of the baths. Here, we emphasize that the Fermi distribution does not capture the contribution of non-Abelian, zero-energy Majorana quasiparticles, which we assume being noninteracting, do not thermalize at temperature T. Nevertheless, the presence of such Majorana quasiparticles are expected to enter in our transport analysis through their contributions in the tunneling/scattering matrix of the transport coefficients. Using relations (11) and (12), one can readily evaluate the equilibrium correlations of particle creation and annihilation operators of the Kitaev chain. The detailed expressions of these correlation matrices are given in Appendix C.

The Hamiltonian $H_{\mathrm{K}}$ (2) reduces to that of a spinless N wire in the absence of pairing ($\Delta =0$). The derivation of a semi-infinite N lead’s normal modes and their statistical properties are straight-forward [54], and are not reproduced here. In the N bath case, we explicitly include the chemical potential in its thermal density matrix, and the tunneling matrix between the N bath and middle wire is now time-independent. If an N bath is kept at temperature T and chemical potential $\mu$, its thermal density matrix for normal modes can be derived from Eq. 12 after substituting $(\hbar \omega _r-\mu )$ in the place of $\hbar \omega _r$. One can also find the equilibrium correlations in terms of particle creation and annihilation operators $c_{l'}^{\dagger },c_{l'}$.

2.2 Majorana Wire

The Hamiltonian of the Majorana wire can also be dealt in the same manner as the Kitaev chain. Nevertheless, the local Hilbert space now is four dimensional ($M=4$), and the Nambu spinor reads as $\psi _{l'}=(c_{l' \uparrow },c_{l' \downarrow },c_{l' \uparrow }^\dagger ,c_{l' \downarrow }^\dagger )^T$, where $c_{l' \sigma }$ annihilates an electron with spin $\sigma =\uparrow ,\downarrow$ at the lattice site $l'$. With the following forms of $U_{l'}$ and $V_{l'}$,

$$\begin{aligned} U_{l'}= \begin{bmatrix} (\epsilon - \gamma ) &{} B &{} 0 &{} -\Delta \\ B&{} -(\epsilon - \gamma ) &{} \Delta &{} 0 \\ 0 &{} \Delta ^* &{} (\epsilon - \gamma ) &{} -B \\ -\Delta ^* &{} 0 &{} -B &{} -(\epsilon - \gamma ) \end{bmatrix},~ \end{aligned}$$

,

$$\begin{aligned} V_{l'}= \frac{1}{2}\begin{bmatrix} -\gamma &{} -\zeta &{} 0 &{}0 \\ \zeta &{} -\gamma &{} 0 &{} 0 \\ 0 &{} 0 &{} \gamma &{} \zeta \\ 0 &{} 0 &{} -\zeta &{} \gamma \end{bmatrix}, \end{aligned}$$

the generalized Hamiltonian, H (1) reduces to the Majorana wire Hamiltonian:

$$\begin{aligned} \frac{H_{\mathrm{M}}}{\hbar }= & {} \sum _{l'=1}^{L-1}\big [-\gamma \sum _{\sigma =\uparrow , \downarrow }(c^{\dag }_{l',\sigma } c_{l'+1,\sigma }+ c^{\dag }_{l'+1,\sigma } c_{l',\sigma })+\zeta (c^{\dag }_{l'+1,\uparrow }c_{l',\downarrow }-c^{\dag }_{l'+1,\downarrow }c_{l',\uparrow } \nonumber \\&+c^{\dag }_{l',\downarrow }c_{l'+1,\uparrow }-c^{\dag }_{l',\uparrow }c_{l'+1,\downarrow })\big ] +2\sum _{l'=1}^{L}\big [(\epsilon -\gamma )\sum _{\sigma =\uparrow , \downarrow }(c^{\dag }_{l',\sigma }c_{l',\sigma }-\frac{1}{2}) \nonumber \\&+B (c^{\dag }_{l',\uparrow }c_{l',\downarrow }+c^{\dag }_{l',\downarrow }c_{l',\uparrow })- (\Delta c^{\dag }_{l',\uparrow }c^{\dag }_{l',\downarrow }+\Delta ^* c_{l',\downarrow } c_{l',\uparrow }) \big ], \end{aligned}$$

(13)

where, $\epsilon$ represents the on-site energy, $\Delta$ is the proximity induced s-wave superconducting pairing potential, B is the magnetic field applied along the direction of wire (say, x-axis), $\gamma$ is the hopping, and $\zeta$ is the Rashba spin-orbit coupling strength. We again assume all the parameters, which have dimension of frequency, to be real. We define x, y, z components of total spin by $\sigma _x^T=\frac{\hbar }{2}\sum _{l'=1}^{L}(c^{\dag }_{l',\uparrow }c_{l',\downarrow }+c^{\dag }_{l',\downarrow }c_{l',\uparrow }), \sigma _y^T=\frac{\hbar }{2} \sum _{l'=1}^{L}i(c^{\dag }_{l',\downarrow }c_{l',\uparrow }-c^{\dag }_{l',\uparrow }c_{l',\downarrow })$, and $\sigma _z^T=\frac{\hbar }{2} \sum _{l'=1}^{L}(c^{\dag }_{l',\uparrow }c_{l',\uparrow }-c^{\dag }_{l',\downarrow }c_{l',\downarrow })$, respectively. We have $[\sigma _x^T,H_{\mathrm{M}}] \ne 0$ for $\zeta \ne 0$, $[\sigma _y^T,H_{\mathrm{M}}] \ne 0$ for $B \ne 0$, and $[\sigma _z^T,H_{\mathrm{M}}] \ne 0$ for $\zeta$ or $B\ne 0$. We later discuss consequences of the above commutation relations (conservation laws) on spin transport in Majorana wires.

Further, we introduce a generalized basis for the Majorana wire:

$$\begin{aligned}&\mathbf{b} \equiv [ b_1, b_2, b_3, b_4, \dots , b_{4L-3}, b_{4L-2}, b_{4L-1}, b_{4L}]^T \nonumber \\&\quad =[c_{1 \uparrow }, c_{1 \downarrow }, c^{\dagger }_{1 \uparrow }, c^{\dagger }_{1 \downarrow } , \dots , c_{L \uparrow }, c_{L \downarrow }, c^{\dagger }_{L \uparrow }, c^{\dagger }_{L \downarrow }]^T. \end{aligned}$$

(14)

In this basis, the Majorana wire Hamiltonian can also be written in a quadratic form as follows, $H_{\mathrm{M}}=\frac{\hbar }{2}\mathbf{b}^{\dagger }{\mathcal {M}}\mathbf{b}=\frac{\hbar }{2}\sum _{l,m}{\mathcal {M}}_{lm}b^{\dagger }_lb_m$, where ${\mathcal {M}}_{lm}$ is an $4L \times 4L$ square matrix and $l,m =1,\dots ,4L$. These $b_l$ operators are mutually related by following relations: $b_{4l'-3}^\dagger =b_{4l'-1}$ and $b_{4l'-2}^\dagger =b_{4l'}$. For this case, the index l of operator $b_l$ can be related to the actual physical lattice site $l'$ applying a set of rules: (i) $l'=l/4$ if l is even and divisible by 4, (ii) $l'=(l+2)/4$ if l is even but not divisible by 4, (iii) $l'=(l+1)/4$ if l is odd and $(l+1)$ is divisible by 4, (iv) $l'=(l+3)/4$ if l is odd but $(l+1)$ is not divisible by 4. The Majorana wire undergoes a topological phase transition at a certain critical magnetic field, $B_c = \sqrt{\Delta ^2+\epsilon ^2}$ [26, 61, 62]. For an applied magnetic field $B>B_c$, this heterostructure is driven into a chiral p-wave topological superconducting phase supporting two zero-energy MBSs at the two ends of the nanowire. However, in the opposite limit $B<B_c$, this system remains in a non-topological phase, and does not host MBS at the edges.

Like the Kitaev chain, the Majorana wire Hamiltonian matrix ${\mathcal {M}}$ can also be diagonalized using the Bogoliubov quasiparticles after solving the Hermitian eigenvalue problem:

$$\begin{aligned} {\mathcal {M}} U_r= \omega _r U_r, \end{aligned}$$

(15)

where, $\omega _r$ the r-th eigenvalue and $U_r$ is the corresponding eigenfunction of the Majorana wire Hamiltonian. In the Nambu spin basis, $U_r$ is represented by an 4L dimensional column vector. The eigenvalues of ${\mathcal {M}}$ can also be grouped into pairs ($\pm \,\omega _r$) with $\omega _r > 0$. If $V_r$ represents an eigenvector with an eigenvalue $-\omega _r$, then $U_r$ and $V_r$ satisfy completeness relation:

$$\begin{aligned} \sum _{r>0}\left( U_r U_r^{\dagger } + V_r V_r^{ \dagger } \right) =\mathbb {1}_{4L}\, , \end{aligned}$$

(16)

and orthogonality relations similar to (7). Hence, we write ${\mathcal {M}}$ as

$$\begin{aligned} {\mathcal {M}}= \sum _{r>0} \omega _r (U_r U^{\dagger }_r- V_r V_r^{ \dagger })\, . \end{aligned}$$

(17)

To define the quasiparticle excitations for the Majorana wire ($H_{\mathrm{M}}$), we first express the eigenvectors in terms of their components:

$$\begin{aligned} U_r= & {} \left[ \Phi _{r }(1), \Psi _{r}(1), \dots , \Phi _{r }(L), \Psi _{r }(L)\right] ^T \, ,\nonumber \\ V_r= & {} \left[ \Psi _{r }^*(1), \Phi _{r}^*(1), \dots , \Psi _{r }^*(L), \Phi _{r }^*(L)\right] ^T\, , \end{aligned}$$

(18)

where $\Phi _r(l')$, and $\Psi _r(l')$ are two component objects defined as $\Phi _r(l')=(\phi _{r \uparrow }(l'), \phi _{r \downarrow }(l'))$, and $\Psi _r(l')=(\psi _{r \uparrow }(l'), \psi _{r \downarrow }(l'))$. Now, we define the fermionic Bogoliubov quasiparticles as

$$\begin{aligned} q_r= & {} U_r^\dagger \mathbf{b}= \sum _{i'=1}^{L} \, \sum _{\sigma = \uparrow , \downarrow } ( \phi _{r \sigma }^{ *}(i')c_{i' \sigma }+\psi _{r \sigma } ^{ *}(i') c_{i' \sigma }^\dagger )\, , \nonumber \\ q_r^{\dagger }= & {} \mathbf{b}^\dagger U_r = \sum _{i'=1}^{L} \sum _{\sigma = \uparrow , \downarrow } ( \phi _{r \sigma }(i')c_{i' \sigma }^\dagger +\psi _{r \sigma }(i') c_{i' \sigma })\, . \end{aligned}$$

(19)

We can also express the electron operators using the Bogoliubov quasiparticles by inverting the above relations (19) while utilizing orthonormality relations of eigenvectors:

$$\begin{aligned} c_{j' \sigma }= & {} \sum _{r>0} ( \phi _{r \sigma } (j') q_r + \psi ^{ *}_{r \sigma } (j') q_r^{\dagger }), \nonumber \\ c_{j' \sigma }^{\dagger }= & {} \sum _{r>0} ( \phi ^{ *}_{r \sigma } (j') q_r^{\dagger } + \psi _{r \sigma }(j')q_r), \end{aligned}$$

(20)

where $\sigma =\uparrow ,\downarrow$. Applying the anti-commutation relations of quasiparticles, it is easy to express the Majorana wire Hamiltonian in a diagonal form:

$$\begin{aligned} H_{\mathrm{M}}&= \sum _{r>0} \hbar \omega _r (2q_r^{ \dagger }\, q_r-1). \end{aligned}$$

Like the Kitaev chain leads, we include the chemical-potential differences as a time-dependent phase in the Majorana wire lead’s tunnel couplings to the middle wire. The thermal density matrix of quasiparticle modes of a semi-infinite Majorana wire bath kept at a temperature T takes exactly similar forms as Eq. 12. In this case, it is assumed that the chemical potentials for up and down spins are the same. This density matrix can be used to calculate the equilibrium correlations in terms of electrons’ creation and annihilation operators $c_{j' \sigma }^{\dagger }, c_{j' \sigma }$. The detailed expressions of equilibrium correlations are given in Appendix E.

The Hamiltonian of an SM in the presence of a magnetic field used in our study can be obtained from $H_{\mathrm{M}}$ (13) by dropping the superconducting pairing $\Delta =0$. Like an N bath, we explicitly include the chemical potential ($\mu$) in the thermal density matrix for an SM bath. The thermal density matrix for its quasiparticle modes can also be derived from Eq. 12 by substituting $(\hbar \omega _r-\mu )$ in the place of $\hbar \omega _r$.

Fig. 1

Schematic of open quantum systems of X-Y-Z configuration where X, Y, Z are made of topological superconductor (TS) of Kitaev chain and normal metal (N) [top], and TS of Majorana wire and semiconductor (SM) [bottom]. The finite middle Y wire acts a system, and the semi-infinite boundary X, Z wires are leads/baths, which are kept at some chemical potential $(\mu _{\mathrm{X}},\mu _{\mathrm{Z}})$ and temperature $(T_{\mathrm{X}}, T_{\mathrm{Z}})$ from the beginning. The individual components of the devices are coupled via tunneling $\gamma _{\mathrm{XY}},\gamma _{\mathrm{YZ}}$ and spin-orbit coupling $\zeta _{\mathrm{XY}},\zeta _{\mathrm{YZ}}$

3 Quantum Langevin Equations and Steady-State Transport

In this section, we discuss the procedure of calculating steady-state nonequilibrium transport properties using a generalized LEGF method. As mentioned before, the hybrid device X-Y-Z consists of three separate wires X, Y, and Z. Total length of the device is $L=L_{\mathrm{X}}+L_{\mathrm{Y}}+L_{\mathrm{Z}}$, and $L_{\alpha '}$ is the length of $\alpha '$ wire for $\alpha '= X,Y,Z$. We assume that both $L_{\mathrm{X}}$ and $L_{\mathrm{Z}}$ are much greater than $L_{\mathrm{Y}}$, thus both X and Z wires can be treated as baths connected to Y. The first site of the Y wire is connected to the $L_{\mathrm{X}}$-th site of left bath, and the $L_{\mathrm{Y}}$-th site of Y wire is connected to the first site of right bath (see Fig. 1). The Hamiltonian of the full hybrid device, which is made of X, Y, Z wires and two contacts, is given by

$$\begin{aligned} H^{\mathrm{F}}=H^{\mathrm{X}}+H^{\mathrm{Y}}+H^{\mathrm{Z}}+H^{\mathrm{XY}}+H^{\mathrm{YZ}}, \end{aligned}$$

(21)

where $H^{\alpha '}$ is the Hamiltonian corresponding to the $\alpha '$ wire, and the contact Hamiltonians for X-Y and Y-Z junctions are $H^{\mathrm{XY}}$ and $H^{\mathrm{YZ}}$, respectively. The exact form of $H^{\alpha '}$ and $H^{\alpha ' \beta '}$ should be chosen according to the type of wire and junction. In our study, $H^{\alpha '}$ can be Hamiltonian of a Kitaev or a Majorana chain of TS or an N or an SM, which we have introduced in the previous section. We choose $H^{\alpha ' \beta '}$ to be either an N junction or an SM junction respectively for our study with the Kitaev chain and the Majorana wire. As discussed in the previous section, $\alpha '$ wire Hamiltonian ($H^{\alpha '}$) can be written in the matrix form, which is denoted by ${\mathcal {H}}^{\alpha '}$.

Here, we generalize the LEGF method so that it can also be applied to the devices of our interest, i.e., a hybrid junction with superconducting (topological/non-topological) leads. Let us briefly indicate the steps leading to generalized quantum Langevin equations of motion. We assume that the leads/baths (X and Z wires) are disconnected from the middle wire (Y) at time $t \le t_0$. Each bath is assumed to be in thermal equilibrium characterized by its temperature $T_\alpha$ and chemical potential $\mu _\alpha$ for $\alpha ={X,Z}$. It should be noted that we do not explicitly include chemical potentials in the thermal density matrix of superconducting baths. Nevertheless, we do have chemical potentials for an N or an SM bath. Here we make some critical approximations to make analytical progress. We assume the mean-field superconducting pair potential, $\Delta$ of the TS wires remains fixed during time-evolution. We have here avoided any self-consistent evolution of $\Delta$ for weak nonequilibrium boundary conditions. This approximation also implies no time-evolution of the superconducting phase of the middle TS wire¹. Further we assume that the temperature scale of the superconductor (T) is much smaller than the critical temperature of the superconductor ($T_c$), so we can approximate $\Delta (T) \sim \Delta (T=0)$ within this low temperature regime. Effect of higher temperature will be manifested as a suppression of conductance peak height [38]. In principle, any nonequilibrium boundary condition such as a voltage or a temperature bias would influence the pairing of the middle TS wire, which is then needed to be determined within a self-consistent mean-field approach in nonequilibrium. However, the qualitative features of our results are not affected as the existence of superconducting gap and the gap closing phenomenon are crucial for such features [35]. It will be an interesting problem to model such nonequilibrium scenario self-consistently and study the feedback mechanism that stabilizes all the mean-field parameters.

In the following, we are mostly interested in studying time-independent transport in hybrid junctions featuring unique nonequilibrium steady-state. Thus, we do not consider transport in a device with two superconducting leads kept at a non-zero chemical potential difference ($\mu _{\mathrm{X}} -\mu _{\mathrm{Z}} \ne 0$), which results in time-dependent steady-state. Nevertheless, for a system with one TS lead and another N/SM lead, we choose to bias the N/SM lead with a non-zero chemical potential while keeping the superconducting one at zero chemical potential to avoid unnecessary complication that may arise from the shifting of quasiparticle energy levels due to the chemical potential. In all the devices of type TS-TS-N/TS-TS-SM, we are mainly interested in calculating different types differential conductances at Y-Z junction. These are local quantities which are not affected by the supercurrent that may develop at X-Y junction due to the phase difference, $\delta \phi =\phi _{\mathrm{X}} -\phi _{\mathrm{Y}}$ when the complex pairing potentials of two superconductors are $\Delta _{\mathrm{X}}=\mid \Delta _{\mathrm{X}} \mid e^{i \phi _{\mathrm{X}}}$ and $\Delta _{\mathrm{Y}}=\mid \Delta _{\mathrm{Y}} \mid e^{i \phi _{\mathrm{Y}}}$ respectively. So, for simplicity we choose the pair potentials of both the superconductors to be real as the phase difference has no role in the differential conductances calculated for right TS-N/SM junction. However for TS-TS-TS at TP, relative phases of pairing potentials play some important role, but we avoided that complicacy by choosing them to be real as this particular case is used only to prove the applicability of NEGF for such junction.

In the subsequent subsections, we study charge, energy, and spin transport in the devices of X-Y-Z configuration in which at least one component is either the Kitaev chain or the Majorana wire. We first discuss the LEGF method and find expressions for steady-state electrical and thermal current for the Kitaev chains. Later, we sketch the approach and derive expressions for electrical and spin current in steady-state for the Majorana wire. We highlight the important steps and provide relatively concise formulas in the main text, and include lengthy expressions and details in the appendices. Due to the complexity of spin-orbit coupled Majorana wire, we do not repeat the tedious calculation of energy current for such a system. However, our method for energy current can readily be applied for the Majorana wire.

3.1 Electrical and Thermal Current: Kitaev Chain

For a hybrid device with spinless TS wires, we choose normal metallic contacts with tunneling rates $\gamma _{\alpha ' \beta '}$, whose Hamiltonians read as

$$\begin{aligned} H^{\alpha ' \beta '}=-\hbar \gamma _{\alpha ' \beta '} (c_{l'}^{\dagger }c_{l'+1}+c_{l'+1}^{\dagger }c_{l'}), \end{aligned}$$

(22)

where $\alpha ' \beta '=\mathrm{XY}$, $l'=L_{\mathrm{X}}$, and $\alpha ' \beta '=\mathrm{YZ}$, $l'=L_{\mathrm{X}}+L_{\mathrm{Y}}$ respectively for X-Y and Y-Z junction. At time $t=t_0$, we connect the baths to the opposite ends of the middle Y wire and look for the steady-state properties of the Y wire at later time t ($t>> t_0+\tau$; $\tau$ is some characteristic time scale for reaching the steady state). As we need to consider the time evolution of total X-Y-Z device, it is convenient to use the full generalized basis: $\mathbf{a} \equiv [a_1, a_2, \dots , a_{2L-1}, a_{2L}]^T$. For $t > t_0$, the Heisenberg equations of motion for the Y wire variables are given by

$$\begin{aligned} {\dot{a}}_l= & {} -i \sum _{m=2L_{\mathrm{X}}+1}^{2L_{\mathrm{XY}}} {\mathcal {H}}^{\mathrm{Y}}_{lm} a_m -i\gamma _{\mathrm{XY}}\sum _{k= 2L_{\mathrm{X}}-1}^{2L_{\mathrm{X}}} (-1)^k a_l \, \delta _{l,k+2}-i\gamma _{\mathrm{YZ}}\sum _{k= 2L_{\mathrm{X}}+1}^{2L_{\mathrm{X}}+2} (-1)^k a_l \, \delta _{l,k-2},\nonumber \\ \end{aligned}$$

(23)

where, $l = 2L_{\mathrm{X}}+1,\dots ,2L_{\mathrm{XY}}$ and $L_{\mathrm{XY}}=L_{\mathrm{X}}+L_{\mathrm{Y}}$. Similarly, the Heisenberg equations for the creation and annihilation operators of X and Z wires read as

$$\begin{aligned} {\dot{a}}_l = -i \sum _{m=1}^{2L_{\mathrm{X}}} {\mathcal {H}}^{\mathrm{X}}_{lm} a_m +i\gamma _{\mathrm{XY}}\sum _{k= 2L_{\mathrm{X}}+1}^{2L_{\mathrm{X}}+2} (-1)^{k+1} a_l \, \delta _{l,k-2} \, , \end{aligned}$$

(24)

for $l = 1, \dots , 2L_{\mathrm{X}}$, and

$$\begin{aligned} {\dot{a}}_l = -i \sum _{m=2L_{\mathrm{XY}}+1}^{2L} {\mathcal {H}}^{\mathrm{Z}}_{lm} a_m +i\gamma _{\mathrm{YZ}} \sum _{k= 2L_{\mathrm{XY}}-1}^{2L_{\mathrm{XY}}} (-1)^{k+1} a_l \, \delta _{l,k+2} \, , \end{aligned}$$

(25)

for $l = 2L_{\mathrm{XY}}+1,\dots ,2L$. The Eqs. 24, 25 are coupled, inhomogeneous, first-order differential equations that can be formally solved for the boundary bath operators by using the retarded Green’s function. Substituting these solutions for the bath operators into Eq. 23, one can rewrite Eq. 23 as Eq. A11, a generalized quantum Langevin equation (see Appendix A for a detailed discussion). The quantum Langevin equations of Y wire variables can be solved in the frequency domain using the Fourier transformation. However, the application of Fourier transformation is reliable only for the systems with unique NESS such that a memory of the initial state of the middle wire is irrelevant. It is worth mentioning here that the Fourier transform method is not applicable in the presence of bound states which prevent equilibration, and one needs to solve Eq. A11 numerically to examine the time evolution in such a case [40, 54]. Therefore, we solve the quantum Langevin equations of the middle wire by Fourier transformation only when a unique NESS is reached. To this end, we first consider the limit $t_0 \rightarrow -\infty$. The Fourier transform of the wire variables are defined as ${\tilde{a}}_l (\omega )=\frac{1}{2\pi }\int _{-\infty }^{\infty } dt \, a_l(t)\, e^{i \omega t}$. We get the following steady-state solutions for ${\tilde{a}}_l(\omega )$ after taking Fourier transform of the quantum Langevin equations of the Y wire variables, $a_l(t)$ (A11):

$$\begin{aligned} {\tilde{a}}_l(\omega )= & {} \sum _{m=2L_{\mathrm{X}}+1}^{2L_{\mathrm{XY}}} {\tilde{G}}^+_{l,m}(\omega )\left( \sum _{k=1,2} {\tilde{\eta }}_{k}^{\mathrm{X}} (\omega )\, \delta _{m,2L_{\mathrm{X}}+k} +\sum _{k=1,2} {\tilde{\eta }}_{k}^{\mathrm{Z}} (\omega )\, \delta _{m,2L_{\mathrm{XY}}+k -2} \right) ,\nonumber \\ \end{aligned}$$

(26)

where $l=2L_{\mathrm{X}}+1,\dots ,2L_{\mathrm{XY}}$. Here, ${\tilde{\eta }}^{\mathrm{X}}_{1,2}(\omega )$ and ${\tilde{\eta }}^{\mathrm{Z}}_{1,2}(\omega )$ are the noise terms arising in the process of integrating out the variables of X and Z bath, respectively (see Appendix A for definition). These noise terms keep track of the nonequilibrium boundary conditions across the middle wire, which we impose in the beginning through the boundary wires. The retarded Green’s function ${\tilde{G}}^+(\omega )$ of the full system in the Fourier domain is defined as

$$\begin{aligned} {\tilde{G}}^+(\omega )={(\omega \mathbb {1}_{2L_{\mathrm{Y}}} - {\mathcal {H}}^{\mathrm{Y}}- {\tilde{\Sigma }}_{\mathrm{X}}^+(\omega ) -{\tilde{\Sigma }}^+_{\mathrm{Z}}(\omega ))^{-1}} ={(\omega \mathbb {1}_{2L_{\mathrm{Y}}} - \tilde{{\mathcal {H}}}^{\mathrm{Y}})}^{-1}\,, \end{aligned}$$

(27)

where ${\tilde{\Sigma }}^+_{\mathrm{X,Z}}$ are the self-energy corrections to the Y wire Hamiltonian originated from its interactions to the respective baths. The effective Hamiltonian matrix of the Y wire which can be non-Hermitian, is given by $\tilde{{\mathcal {H}}}^{\mathrm{Y}}={\mathcal {H}}^{\mathrm{Y}}+ {\tilde{\Sigma }}^+_{\mathrm{X}}(\omega ) +{\tilde{\Sigma }}^+_{\mathrm{Z}}(\omega )$. These ${\tilde{\Sigma }}_{\mathrm{X,Z}}^+$ are square matrices of dimension $2L_{\mathrm{Y}} \times 2L_{\mathrm{Y}}$. The components of the self-energy terms ${\tilde{\Sigma }}^+_{\mathrm{X,Z}}$ are as following:

$$\begin{aligned}{}[{\tilde{\Sigma }}_{\mathrm{X}}^+(\omega )]_{lm}= & {} \gamma _{\mathrm{XY}} ^2 \sum _{k=2L_{\mathrm{X}}-1}^{2L_{\mathrm{X}}} {{\tilde{G}}^{\mathrm{X +}}_{k,k}(\omega )}\, \delta _{l,k+2}\, \delta _{l,m} \nonumber \\&-\gamma _{\mathrm{XY}} ^2 \sum _{\begin{array}{c} k,k'=2L_{\mathrm{X}}-1 \\ k \ne k' \end{array}}^{2L_{\mathrm{X}}}{{\tilde{G}}^{\mathrm{X +}}_{k,k'}(\omega )}\, \delta _{l,k+2}\, \delta _{m,k'+2} \, ,\nonumber \\ {[{\tilde{\Sigma }}^+_{\mathrm{Z}}(\omega )]}_{lm}= & {} \gamma _{\mathrm{YZ}}^2 \sum _{k=1}^{2}{{\tilde{G}}^{\mathrm{Z+}}_{k,k}(\omega )}\,\delta _{l,2L_{\mathrm{XY}}-2+k}\, \delta _{l,m} \nonumber \\&-\gamma _{\mathrm{YZ}} ^2 \sum _{\begin{array}{c} k,k'=1 \\ k \ne k' \end{array}}^{2}{{\tilde{G}}^{\mathrm{Z +}}_{k,k'}(\omega )}\, \delta _{l,2L_{\mathrm{XY}}-2+k}\, \delta _{m,2L_{\mathrm{XY}}-2+k'} \, , \end{aligned}$$

(28)

where $l,m= 2L_{\mathrm{X}}+1,\dots ,2L_{\mathrm{XY}}$. Here, ${\tilde{G}}^{\alpha +}_{l,m}(\omega )$ is the retarded Green’s function of isolated bath wires ($\alpha ={X,Z}$). For ${\tilde{G}}^{\mathrm{X +}}_{l,m}(\omega )$, $l,m =2L_{\mathrm{X}}-1,2L_{\mathrm{X}}$ correspond to the right most site (i.e., $L_{\mathrm{X}}$-th site) of the X reservoir, whereas in case of ${\tilde{G}}^{\mathrm{Z +}}_{l,m}(\omega )$, $l,m=1,2$ correspond to the left most site (i.e., $L_{\mathrm{XY}}+1$-th site) of the Z reservoir in the full X-Y-Z device. The detailed definitions and expressions of ${\tilde{G}}^{\alpha +}_{l,m}(\omega )$ are given in Appendix B. Since $\tilde{{\mathcal {H}}}^{\mathrm{Y}}$ is a block diagonal matrix, numerical values of ${\tilde{G}}^+(\omega )$ can be calculated by inverting $(\omega \mathbb {1}_{2L_{\mathrm{Y}}} -\tilde{{\mathcal {H}}}^{\mathrm{Y}})$.

We further write steady-state solutions for some of the bath variables ${\tilde{a}}_l(\omega )$ (A3) defined at the edges of the baths. For example, these ${\tilde{a}}_{2L_{\mathrm{X}}-1}(\omega )$ and ${\tilde{a}}_{2L_{\mathrm{X}}}(\omega )$ for the left X bath read as

$$\begin{aligned} {\tilde{a}}_{2L_{\mathrm{X}}-1}(\omega )\gamma _{\mathrm{XY}}= & {} -{\tilde{\eta }}^{\mathrm{X}}_1(\omega )- \sum _{m=2L_{\mathrm{X}}+1}^{2L_{\mathrm{X}}+2}{[{\tilde{\Sigma }}^+_{\mathrm{X}}(\omega )]}_{2L_{\mathrm{X}}+1,m} \, {\tilde{a}}_{m}(\omega ), \\ {\tilde{a}}_{2L_{\mathrm{X}}}(\omega )\gamma _{\mathrm{XY}}= & {} {\tilde{\eta }}^{\mathrm{X}}_2(\omega )+ \sum _{m=2L_{\mathrm{X}}+1}^{2L_{\mathrm{X}}+2}{[{\tilde{\Sigma }}^+_{\mathrm{X}}(\omega )]}_{2L_{\mathrm{X}}+2,m} \, {\tilde{a}}_{m}(\omega ) \, . \end{aligned}$$

These boundary variables of the baths would be useful in evaluating the transport coefficients through the middle wire, which we discuss below.

For these hybrid devices with spinless particle, the transport coefficients of interest are electrical (charge) and thermal (energy) conductance. The electrical conductance measurements in such devices are sensitive to the emergence of the Majorana zero modes at the edges of the Kitaev chains. However, the charge neutrality of the Majorana quasiparticles poses a challenge to unambiguous detection of such topologically protected modes through electrical conductance. It is rather interesting to probe these charge-neutral modes through the thermal transport, which we also evaluate here [8, 10,11,12, 43, 44]. Since the total number of particles (spinless electrons) is not conserved for the TS wires, the particle current of spinless electrons is, in general, not well-defined inside the TS wires. Nevertheless, the total particle number (of spinless electrons) is conserved for N wires or N junctions, and we mostly define charge currents carried by spinless electrons at those segments. Using the conservation of particles (of spinless electrons) at the junctions, we describe the charge/electrical current across the links after multiplying particle current by electron’s charge e:

$$\begin{aligned} J^e_{\alpha ' \beta '}=ie\gamma _{ \alpha ' \beta '}\langle (c_{l'}^{\dagger }c_{l'+1}- c_{l'+1}^{\dagger }c_{l'})\rangle , \end{aligned}$$

(29)

where again $\alpha ' \beta '=\mathrm{XY}$, $l'=L_{\mathrm{X}}$, and $\alpha ' \beta '=\mathrm{YZ}$, $l'=L_{\mathrm{XY}}$ respectively for the X-Y and Y-Z junction. The expectation $\langle .. \rangle$ denotes averaging over the initial density matrix of the baths. To find the electrical current using Eq. 29, we first evaluate the noise-noise correlations for the baths. In Appendix C, we outline the procedure of finding the noise-noise correlations, and we also list there all the noise-noise correlations relevant for the calculation of electrical current at the junctions. Next, we rewrite the currents in Eq. 29 in the following compact form:

$$\begin{aligned} J^e_{\mathrm{XY}}= & {} - 2e \gamma _{\mathrm{XY}} \, \text {Im}[\langle a_{2L_{\mathrm{X}}+1}^\dagger (t) a_{2L_{\mathrm{X}}-1} (t) \rangle ]\, , \end{aligned}$$

(30)

$$\begin{aligned} J^e_{\mathrm{YZ }}= & {} - 2e \gamma _{\mathrm{YZ}} \, \text {Im}[\langle a_{2L_{\mathrm{XY}}+1}^\dagger (t) a_{2L_{\mathrm{XY}}-1} (t) \rangle ] \, . \end{aligned}$$

(31)

First we take the Fourier transformation of the bath and wire variables in Eqs. 30 and 31, then those variables are substituted with Eq. 26 and the Fourier transformed version of Eq. A3. Using the noise-noise correlations (C3–C5), we finally obtain the analytical expressions for $J^e_{\mathrm{XY}}$ and $J^e_{\mathrm{YZ}}$ (see Appendix F). Due to the absence of conservation of total number of particles inside a superconductor, we observe that $J^e_{\mathrm{XY}} \ne J^e_{\mathrm{YZ}}$ for N-TS-N devices for arbitrary nonequilibrium boundary conditions [33]. In Sect. 4, we show an interesting conversion of a part of the injected electrical current to a Cooper pair current inside the middle TS wires. Nevertheless, $J^e_{\mathrm{XY}}=J^e_{\mathrm{YZ}}$ for a symmetric bias, $\mu _{\mathrm{X}}=-\mu _{\mathrm{Z}}$, in an N-TS-N device. The emergence of a zero-energy MBS is expected to manifest a quantized zero-bias peak of height $2e^2/h$ in the zero-temperature DEC. For N-TS-N devices placed under a symmetric bias, we define DEC as $2\left[ \frac{dJ^e_{\mathrm{XY}}}{dV}\right]$ or $2\left[ \frac{dJ^e_{\mathrm{YZ}}}{dV}\right]$, where $\mu _{\mathrm{X}}-\mu _{\mathrm{Z}}=eV$². Writing $\mu _{\mathrm{X}}=-\mu _{\mathrm{Z}}=\mu$, we can also express the DEC in such a device as $e\left[ \frac{dJ^e_{\mathrm{XY}}}{d\mu }\right]$ or $e\left[ \frac{dJ^e_{\mathrm{YZ}}}{d \mu }\right]$. For TS-N-N ans TS-TS-N devices, we apply an asymmetrical bias by setting $\mu _{\mathrm{X}}=0$ and $\mu _{\mathrm{Z}}= eV$, and change eV to find DEC. For such systems, we are only interested in the DEC at the Y-Z junction, which is defined as $\left[ \frac{dJ^e_{\mathrm{YZ}}}{d V}\right]$.

The expressions for $J^e_{\mathrm{XY}}$ and $J^e_{\mathrm{YZ}}$ can be written in a simple and neat Landauer current form in the presence of a temperature bias and zero chemical potentials, $\mu _{\mathrm{X}}=\mu _{\mathrm{Z}}=0$. There, $J^e_{\mathrm{XY}}$ (F1) and $J^e_{\mathrm{YZ}}$ (F2) are simplified as a product of a frequency-dependent transmission coefficient ${\mathcal {T}}(\omega )$ and a difference between the Fermi functions of the boundary leads. We note that generally $J^e_{\mathrm{XY}} \ne J^e_{\mathrm{YZ}}$ for a TS-TS-N and a TS-TS-TS at TP in the above scenario. However, the expression F3 implies that the $J^e_{\mathrm{XY}}=J^e_{\mathrm{YZ}}$ for an N-TS-N device in the above limit of bias. Thus, we get the following expression for the thermoelectric current generated by a temperature bias from the X and Z bath with temperature $T_{\mathrm{X}}$ and $T_{\mathrm{Z}}$, respectively:

$$\begin{aligned} J^e_{\mathrm{XY}}=J^e_{\mathrm{YZ}}=\int _{- \infty }^{ \infty }d \omega \,e\, {\mathcal {T}}(\omega ) \left( f(\omega ,T_{\mathrm{X}}) -f(\omega ,T_{\mathrm{Z}})\right) =\frac{\pi ^2 k_B^2 eT}{3 \hbar ^2}\left[ \frac{ \partial {\mathcal {T}}}{\partial \omega } \right] _{\omega =0} \Delta T,\nonumber \\ \end{aligned}$$

(32)

where ${\mathcal {T}}(\omega )~={\mathcal {T}}^1_{\mathrm{XY}}(\omega )+{\mathcal {T}}^2_{\mathrm{XY}}(\omega )$, whose explicit forms are given in Appendix F. We obtain the last expression in the above equation by making a linear response expansion for small temperature differences $\Delta T \ll T$, where $T_{\mathrm{X}}=T+\Delta T/2$ and $T_{\mathrm{Z}}=T-\Delta T/2$. As discussed above, all over the calculation we assume T itself is also small so that we can ignore the effect of T on the superconducting pair potential. We show later that $\left[ \frac{ \partial {\mathcal {T}}}{\partial \omega } \right] _{\omega =0}$ shows a sharp dip near the topological phase transition of the middle TS wires. We here propose to experimentally probe the topological phase transition in TS wires by measuring such a large dip in the thermoelectric current.

While electrical currents and differential conductances are extensively explored in an N-TS junction for the search of elusive Majorana fermions, the thermal/energy currents are relatively less studied in such junction [8, 43,44,45]. However, recent experiments [10,11,12] have suggested it might be possible to probe such very small thermal conductances arising from a few conducting channels rather accurately. Motivated by these developments, we derive the expressions of energy currents and linear-response thermal conductance in hybrid devices made of Kitaev chains. Further motivation stems from the fact that while electrical current is not the same across such mean-field models of TS wires, the energy current remains the same across the TS wires, which we explicitly demonstrate in Appendix H. Using the continuity equation for the conserved energy across the junction between two wires, we derive the following expressions for energy current at X-Y and Y-Z junctions:

$$\begin{aligned} J^{u}_{\mathrm{XY}}= & {} {\bar{J}}^q_{\mathrm{XY}}+{\bar{J}}^p_{\mathrm{XY}} \nonumber \\= & {} 2\hbar \gamma _{\mathrm{XY}} \left( \gamma _{\mathrm{Y}} \mathrm{Im}[\langle c^\dagger _{L_{\mathrm{X}}+2}(t) c_{L_{\mathrm{X}}} (t)\rangle ]+ \Delta _{\mathrm{Y}}\mathrm{Im}[\langle c^\dagger _{L_{\mathrm{X}}+2}(t) c^\dagger _{L_{\mathrm{X}}} (t)\rangle ] \right) \nonumber \\&+ 2 \hbar \gamma _{\mathrm{XY}} \epsilon _{\mathrm{Y}} \mathrm{Im}[\langle c^\dagger _{L_{\mathrm{X}}+1}(t) c_{L_{\mathrm{X}}} (t)\rangle ] \, , \end{aligned}$$

(33)

$$\begin{aligned} J^u_{\mathrm{YZ}}= & {} {\bar{J}}^q_{\mathrm{YZ}}+{\bar{J}}^p_{\mathrm{YZ}} \nonumber \\= & {} 2 \hbar \gamma _{\mathrm{YZ}} \left( \gamma _{\mathrm{Y}} \mathrm{Im}[\langle c^\dagger _{L_{\mathrm{XY}}+1}(t) c_{L_{\mathrm{XY}-1}} (t)\rangle ]-\Delta _{\mathrm{Y}} \mathrm{Im}[\langle c^\dagger _{L_{\mathrm{XY}}+1}(t) c^\dagger _{L_{\mathrm{XY}-1}} (t)\rangle ]\right) \nonumber \\&+ 2\hbar \gamma _{\mathrm{YZ}} \epsilon _{\mathrm{Y}}\mathrm{Im}[\langle c^\dagger _{L_{\mathrm{XY}}+1}(t) c_{L_{\mathrm{XY}}} (t)\rangle ] \, , \end{aligned}$$

(34)

assuming the pairing potential $\Delta _{\mathrm{Y}}$ is real. The first parts within the big parenthesis in both (33) and (34) denote ${\bar{J}}^q_{\alpha ' \beta '}$, and the other parts of the expressions represent ${\bar{J}}^p_{\alpha ' \beta '}$, which is zero when the on-site energy of the Y wire is zero. As shown in the above expressions, ${\bar{J}}^q_{\alpha ' \beta '}$ can be separated into two parts; one is electronic part which has explicit dependence on the hopping parameter, and another is Cooper pair part which is explicitly related to the pairing potential. Like the electrical currents, we can again find explicit expressions of the energy currents in the steady state by using the steady-state solutions of the variables appearing in Eqs. 33 and 34. Due to the conservation of total energy in the middle wire in our all studied models, the energy current remains the same across the middle wire including for a TS wire in an N-TS-N device. The detailed expression of the energy current is given in Appendix G. When the boundary lead wires are kept at a finite temperature bias and at zero chemical potential ($\mu _{\mathrm{X}}=\mu _{\mathrm{Z}}=0$), the expression of $J^{u}_{\mathrm{XY}}$ and $J^{u}_{\mathrm{YZ}}$ (G3) can be expressed in a simple Landauer current form as a multiplication of a frequency-dependent transmission coefficient with a difference between the Fermi functions of the boundary baths [54, 55, 57].

In the linear response regime, we can further simplify the Landauer form of energy current by assuming $T_{\mathrm{X}}=T+\Delta T/2$ and $T_{\mathrm{Z}}=T-\Delta T/2$ along with $\mu _{\mathrm{X}}=\mu _{\mathrm{Z}}=0$, where the temperature difference $\Delta T=T_{\mathrm{X}}-T_{\mathrm{Z}}$ is much small compared to the mean temperature T. By expanding the Fermi functions about T, we then write for the energy current:

$$\begin{aligned} J^{u}_{\mathrm{XY}}= J^{u}_{\mathrm{YZ}}= & {} \int _{- \infty }^{ \infty } d \omega \, \hbar {\mathcal {A}}(\omega ) \left( f(\omega ,T_{\mathrm{X}}) -f(\omega ,T_{\mathrm{Z}})\right) \nonumber \\= & {} \int _{- \infty }^{ \infty }d \omega \,\hbar \left( {\mathcal {A}}(0)+ \left[ \frac{ \partial {\mathcal {A}}}{\partial \omega } \right] _{\omega =0}\omega \right) \left[ \frac{\partial f}{\partial T}\right] _{T}\Delta T \nonumber \\= & {} \frac{\pi ^2 k_B^2 }{3 \hbar }\left[ \frac{ \partial {\mathcal {A}}}{\partial \omega } \right] _{\omega =0} T \Delta T, \end{aligned}$$

(35)

where ${\mathcal {A}}(\omega )={\mathcal {A}}^1(\omega )+{\mathcal {A}}^2(\omega )$ (Appendix G) is the frequency-dependent transmission coefficient of energy across the middle Y wire due to a temperature bias. We now define linear-response thermal conductance $G_{\mathrm{T}}$ as

$$\begin{aligned} G_{\mathrm{T}}=\frac{J^{u}_{\mathrm{XY}}}{\Delta T}=\frac{\pi ^2 k_B^2T }{3 \hbar }\left[ \frac{ \partial {\mathcal {A}}}{\partial \omega } \right] _{\omega =0}, \end{aligned}$$

(36)

where we can identify $[\partial {\mathcal {A}}/\partial \omega ]_{\omega =0} ={\mathcal {T}}(\omega =0)$ in an N-TS-N device and other devices with TS leads. We also notice that ${\mathcal {A}}(\omega )= \omega {\mathcal {T}}(\omega )$ (32) for all values of $\omega$ when $\epsilon _{\mathrm{Y}}=0$ in an N-TS-N device. The last result is probably due to the existence of two independent Majorana (conduction) channels without any scattering between them when $\epsilon _{\mathrm{Y}}=0$ [33]. $G_{\mathrm{T}}$ shows some interesting features across the TP in an N-TS-N device. We discuss the properties of steady-state energy current and $G_{\mathrm{T}}$ for different hybrid devices in Sect. 4.2.

3.2 Electrical and Spin Current: Majorana Wire

Next, we explore steady-state quantum transport in spinful models of superconductors and semiconductors. Apart from the electrical currents, we intend to find the features of spin current in such devices, which have at least one part made of Majorana wires. Motivated by the experimental set-ups in [3, 6], we now consider that the tunneling Hamiltonians also include spin-orbit couplings. Thus, the tunneling Hamiltonians for the X-Y and Y-Z junctions are the following:

$$\begin{aligned} H^{\alpha ' \beta '}= & {} -\hbar \gamma _{\alpha ' \beta '} \sum _{\sigma =\uparrow , \downarrow } (c^{\dag }_{l',\sigma } c_{l'+1,\sigma } + c^{\dag }_{l'+1,\sigma } c_{l',\sigma } ) \nonumber \\&+ \hbar \zeta _{\alpha ' \beta ' }(c^{\dag }_{l'+1,\uparrow }c_{l',\downarrow }-c^{\dag }_{l'+1,\downarrow }c_{l',\uparrow }+h.c.)\, , \end{aligned}$$

(37)

where $\alpha ' \beta ' =\mathrm{XY}$, $l'=L_{\mathrm{X}}$, and $\alpha ' \beta ' = \mathrm{YZ}$, $l'=L_{\mathrm{XY}}$ respectively for the X-Y and Y-Z junction. Here, $\gamma _{\alpha ' \beta '}$ represents tunneling rate, and $\zeta _{\alpha ' \beta ' }$ represents the strength of Rashba spin-orbit coupling at the tunnel junction.

The basic structure of the generalized LEGF is same for Kitaev and Majorana wire lead. However, the calculation of steady-state electrical and spin currents at the junctions of these devices made of Majorana/SM wire is a bit cumbersome due to the presence of spin degrees of freedom and spin-orbit coupling in the Majorana and SM wires. So without going into much details, we highlight some of the main steps to find the steady-state solutions of the Heisenberg’s equations of motion in Appendix D. We also explain the method of calculating noise-noise correlations for Majorana wire leads in Appendix E.

We define the total electrical current of both spin components of electrons at the junctions between the wires. Again, the total particle density of electrons at the junctions is conserved, and we use the continuity equations to write the charge currents of electrons across the junctions. We immediately get the following expressions for electrical current from the X bath to the Y wire and from the Y wire to the Z bath, respectively:

$$\begin{aligned} {\mathcal {J}}^e_{\mathrm{XY}}= & {} -2e\gamma _{\mathrm{XY}} \sum _{\sigma = \uparrow , \downarrow }\mathrm{Im}[\langle c^{\dagger }_{L_{\mathrm{X}}+1, \sigma }(t)\, c_{L_{\mathrm{X}}, \sigma }(t) \rangle ] \nonumber \\&-2e \zeta _{\mathrm{XY}} \sum _{\begin{array}{c} \sigma ,\sigma '=\uparrow ,\downarrow \\ \sigma \ne \sigma ' \end{array}}(-1)^{r(\sigma )}\mathrm{Im}[\langle c^{\dagger }_{L_{\mathrm{X}}+1, \sigma }(t)\, c_{L_{\mathrm{X}}, \sigma ' }(t) \rangle ], \end{aligned}$$

(38)

$$\begin{aligned} {\mathcal {J}}^e_{\mathrm{YZ}}= & {} -2e\gamma _{\mathrm{YZ}} \sum _{\sigma = \uparrow , \downarrow }\mathrm{Im}[\langle c^{\dagger }_{L_{\mathrm{XY}}+1, \sigma }(t)\, c_{L_{\mathrm{XY}}, \sigma }(t) \rangle ] \nonumber \\&-2e \zeta _{\mathrm{YZ}} \sum _{\begin{array}{c} \sigma ,\sigma '=\uparrow ,\downarrow \\ \sigma \ne \sigma ' \end{array}}(-1)^{r(\sigma )}\mathrm{Im}[\langle c^{\dagger }_{L_{\mathrm{XY}}+1, \sigma }(t)\, c_{L_{\mathrm{XY}}, \sigma ' }(t) \rangle ], \end{aligned}$$

(39)

where, $r(\uparrow ) (r(\downarrow ))=1(2)$. These expressions can be directly employed to study time evolution of electrical currents at the junctions after numerically solving the equations of motion for annihilation and creation operators of the full device. To calculate junction currents at NESS using the generalized LEGF, we first take the Fourier transformation of the expressions (38) and (39), then substitute the variables (${\tilde{b}}_l(\omega )$) of Y wire (D4) and the baths into it. Finally, plugging the noise-noise correlations, we can derive the electrical currents at the junctions in the steady state. Like the TS-N-N and TS-TS-N devices, we set $\mu _{\mathrm{X}}=0$ and $\mu _{\mathrm{Z}}= eV$, and tune eV to find DEC in TS-SM-SM and TS-TS-SM devices. For such systems, we are again interested in the DEC at Y-Z junction, which is defined as $\left[ \frac{d{\mathcal {J}}^e_{\mathrm{YZ}}}{d V}\right]$.

Recently, the spin transport has been experimentally explored in semiconductor -superconductor hybrid devices [17]. In order to calculate spin transport in devices with the Majorana wires, we need to define spin currents at the junctions. The expressions for spin currents at the junctions can be obtained by applying the continuity equation for the local spin density around the junctions. We define local spin density operator for the x-component of spin at site $l'$ as $\sigma _x^{l'}=\frac{\hbar }{2}(c^{\dag }_{l',\uparrow }c_{l',\downarrow }+c^{\dag }_{l',\downarrow }c_{l',\uparrow })$. Employing the continuity equations at the terminal sites of the Y wire ($l'=L_{\mathrm{X}}+1,L_{\mathrm{XY}}$), we derive below the expressions of x-component of spin current from the X bath to the Y wire and from the Y wire to the Z bath, respectively:

$$\begin{aligned} {\mathcal {J}}^{s_x}_{\mathrm{XY}}= & {} -\hbar \gamma _{\mathrm{XY}} \sum _{\begin{array}{c} \sigma ,\sigma '=\uparrow ,\downarrow \\ \sigma \ne \sigma ' \end{array}} \mathrm{Im}[\langle c^{\dagger }_{L_{\mathrm{X}}+1, \sigma }(t)\, c_{L_{\mathrm{X}}, \sigma '}(t) \rangle ] \nonumber \\&+ \hbar \zeta _{\mathrm{XY}} \sum _{\sigma =\uparrow ,\downarrow } (-1)^{r(\sigma )}\mathrm{Im}[\langle c^{\dagger }_{L_{\mathrm{X}}+1, \sigma }(t)\, c_{L_{\mathrm{X}}, \sigma }(t) \rangle ], \end{aligned}$$

(40)

$$\begin{aligned} {\mathcal {J}}^{s_x}_{\mathrm{YZ}}= & {} -\hbar \gamma _{\mathrm{YZ}} \sum _{\begin{array}{c} \sigma ,\sigma '=\uparrow ,\downarrow \\ \sigma \ne \sigma ' \end{array}} \mathrm{Im}[\langle c^{\dagger }_{L_{\mathrm{XY}}+1, \sigma }(t)\, c_{L_{\mathrm{XY}}, \sigma '}(t) \rangle ] \nonumber \\&+ \hbar \zeta _{\mathrm{YZ}} \sum _{\sigma =\uparrow ,\downarrow } (-1)^{r(\sigma )}\mathrm{Im}[\langle c^{\dagger }_{L_{\mathrm{XY}}+1, \sigma }(t)\, c_{L_{\mathrm{XY}}, \sigma }(t) \rangle ], \end{aligned}$$

(41)

where, $r(\uparrow ) (r(\downarrow ))=1(2)$. We later apply the above expressions (41,41) to find the expectation value of spin currents in different devices made of Majorana wires. Inspired by the extensive applications of DEC in probing the topological state of superconductors in experiments, we here introduce differential spin conductance (DSC), which we write as $\left[ \frac{d{\mathcal {J}}^{s_x}_{\mathrm{YZ}}}{d V}\right]$ for the Y-Z junction. In Sect. 4.4, we discuss some novel features of zero-temperature DSC in various spinful TS wire devices by tuning $\mu _{\mathrm{Z}}$ from $-eV$ to eV while keeping $\mu _{\mathrm{X}}=0$.

4 Results and Discussion

In this section, we discuss many new results that we obtained from our steady-state electrical, thermal, and spin current formulas which have been derived using the generalized LEGF method in the previous section. We have earlier argued that the generalized LEGF technique is applicable to those systems which reach a unique NESS in the long-time limit $t \rightarrow \infty$. For example, the LEGF can easily be applied to an N-TS-N device because such a junction satisfies the requirements for the unique NESS [33, 41]. However, an X-Y-Z configuration with a superconducting wire for X lead does not attain a unique NESS unless the Z wire is an N/SM or a topological superconductor at TP. The basic criterion for achieving a unique NESS in these devices is the absence of a bound state with energy within the full system’s bulk energy gap. There is no bound state in the hybrid device when the energy spectrum of Z wire is gapless for an N/SM or a TS at TP, thus the currents at the junctions in the long-time limit ($t \rightarrow \infty$) do not depend on the initial density matrix of the middle wire indicating a unique NESS [40]. Below, we validate our steady-state current expressions by comparing them with the currents’ long-time values calculated directly through numerically solving the Heisenberg equations of the full device using the time-dependent Green’s function techniques [40, 54].

It should be noted that we choose some arbitrary initial conditions (e.g., local density) for the Y wire in our direct time-evolution numerics [40]. For example, we can choose $n_{l'}$ as the initial number of spinless fermions at $l'$-th site of the Y wire for spinless models. Since every site of a spinful system like an SM/Majorana wire can be filled by two spins; $n_{l',\sigma }$ represents the number of fermions with spin $\sigma$ for spinful models. In our first principle calculation, we numerically evaluate some dynamical quantities like electric and thermal currents. We infer that the system has reached a unique NESS if these dynamical quantities become independent of time and the initial densities of Y wire in the long-time limit. In Appendix H, for the shake of completeness, we briefly discuss the aforesaid numerical method for studying the time evolution of the currents in an X-Y-Z device made of Kitaev chain and N wires.

The rest of this section is divided in the four subsections to discuss (1) charge current carried by Cooper pairs in an N-TS-N device where the TS is a Kitaev chain and in an SM-TS-SM device with a Majorana wire, (2) thermal and thermoelectric currents and conductances, (3) electrical current and DEC in hybrid devices with Kitaev chain leads, and (4) spin current, DEC and DSC in hybrid devices with Majorana wire leads.

Fig. 2

Comparison of electrical currents $J^e_{\mathrm{XY}}(t), J^e_{\mathrm{YZ}}(t)$ obtained from the first-principle/direct time-evolution numerics and the generalized LEGF method (dashed lines) for an N-TS-N device made of a Kitaev chain. The initial numbers of spinless electrons $(n_{l'})$ at the middle TS wire used for the time-evolution numerics are indicated on the headings. In both panels, $L_{\mathrm{X}}=L_{\mathrm{Z}}=2000, L_{\mathrm{Y}}=10$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=\gamma =1$, $\gamma _{\mathrm{Y}}=0.5$, $\Delta _{\mathrm{X}}=\Delta _{\mathrm{Z}}=0,\Delta _{\mathrm{Y}}=0.15$, $\epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Z}}=0$, $\epsilon _{\mathrm{Y}}=0.01$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.25$, $T_{\mathrm{X}}=0.02,T_{\mathrm{Z}}=0.02$ and $\mu _{\mathrm{X}}=0.2,\mu _{\mathrm{Z}}=-0.4$. All above parameters except lengths are in units of $\gamma$

4.1 Charge Current Carried by Cooper Pairs in N-TS-N and SM-TS-SM

First, we validate the generalized LEGF method for an N-TS-N and an SM-TS-SM device by comparing the steady-state currents with the results obtained by numerically evaluating the time-evolution of the full devices. The first-principle/direct numerics shows that the long-time transport in the N-TS-N (SM-TS-SM) device is independent of the initial conditions for the finite TS wire when the energy band of the N (SM) wire is wider than that of the TS wire. Therefore, the results obtained from the time-evolution numerics and the generalized LEGF should match with each other in such scenarios. In Fig. 2, we compare the LEGF formulas with the time-evolution numerics. We find very good agreement between the two values at a long time. Moreover, the values of electrical currents at a long time do not depend on the initial density of the middle TS wire as shown in panels of Fig. 2a, b. Here, our time-evolution numerics shows that the currents at X-Y and Y-Z junctions reach constant non-zero values respectively, 0.02049 and 0.02751 for both the initial conditions of the Y wire. On the other hand, the steady-state calculation yields $J^e_{\mathrm{XY}}=0.020492$ and $J^e_{\mathrm{YZ}}=0.02751$. Thus, both the calculations match excellently. We further observe such a good matching for other sets of parameters as well as for an SM-TS-SM device where the LEGF values of ${\mathcal {J}}^e_{\mathrm{XY}}=0.018851$ and ${\mathcal {J}}^e_{\mathrm{YZ}}=0.020826$ do agree with the time-evolution numerics as shown in Fig. 3. Here the numerical values of the junction currents are in units of $\gamma e$. These results confirm the validity of generalized LEGF approach for steady-state electrical transport in N-TS-N and SM-TS-SM devices.

Fig. 3

Comparison of electrical currents ${\mathcal {J}}^e_{\mathrm{XY}}(t), {\mathcal {J}}^e_{\mathrm{YZ}}(t)$ obtained from the first-principle/direct time-evolution numerics and the generalized LEGF method (dashed lines) for an SM-TS-SM device of a Majorana wire. The initial density $n_{l',\sigma }=c^{\dagger }_{l',\sigma }c_{l',\sigma }$ with $\sigma =\uparrow ,\downarrow$ of the middle TS wire for the direct time-evolution numerics is indicated on the heading of the top panel. In both panels, $L_{\mathrm{X}}=L_{\mathrm{Z}}=1500, L_{\mathrm{Y}}=10$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=\gamma =1$, $\gamma _{\mathrm{Y}}=0.5$, $\Delta _{\mathrm{Y}}=0.15,\Delta _{\mathrm{Y}}=\Delta _{\mathrm{Z}}=0, \zeta _{\mathrm{X}}=\zeta _{\mathrm{Y}}=\zeta _{\mathrm{Z}}=\zeta _{\mathrm{XY}}=\zeta _{\mathrm{YZ}}=0.2, B_{\mathrm{X}}=B_{\mathrm{Y}}=B_{\mathrm{Z}}=0.3$, $\epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Z}}=0, \epsilon _{\mathrm{Y}}=0.01$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.25$, $T_{\mathrm{X}}=T_{\mathrm{Z}}=0.02$ and $\mu _{\mathrm{X}}=0.2,\mu _{\mathrm{Z}}=-0.4$. All above parameters except lengths are in units of $\gamma$

Fig. 4

Dependence of different components of the electrical currents on position of the lattice site inside the middle superconducting wire of a N-TS-N device made of a Kitaev chain. The parameters are $L_{\mathrm{Y}}=20$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=\gamma =1$, $\gamma _{\mathrm{Y}}=0.5$, $\Delta _{\mathrm{Y}}=0.1$, $\epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Z}}=0,\epsilon _{\mathrm{Y}}=0.05$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.25$, $T_{\mathrm{X}}=T_{\mathrm{Z}}=0.02$ and $\mu _{\mathrm{X}}=0.2,\mu _{\mathrm{Z}}=-0.4$. All above parameters except lengths are in units of $\gamma$

We find from Fig. 2 (Fig. 3) that $J^e_{\mathrm{XY}} \ne J^e_{\mathrm{YZ}}$ for an N-TS-N device [33] (or ${\mathcal {J}}^e_{\mathrm{XY}} \ne {\mathcal {J}}^e_{\mathrm{YZ}}$ for an SM-TS-SM device [25]), since the total and local electron number are not conserved inside a superconductor modeled by a mean-field Hamiltonian as ours in this paper. The violation of local and global electron number conservation can be expressed, respectively, as $[c^{\dagger }_{l'}c_{l'},H^{\mathrm{Y}}] \ne 0$ and $[\sum _{l'=L_{\mathrm{X}}+1}^{L_{\mathrm{X}}+L_{\mathrm{Y}}}c^{\dagger }_{l'}c_{l'},H^{\mathrm{Y}}] \ne 0$ when $\Delta _{\mathrm{Y}} \ne 0$ for a spinless Y wire. Due to the local as well global non-conservation of electron number inside a superconductor, the electron current density is same neither locally nor globally, i.e., the incoming electron current density is not equal to the outgoing electron current density in a two-terminal geometry. We can rather say that such superconductor modelled by a mean-field Hamiltonian acts as a reservoir for its electrons. Thus, the single-electron charge current entering into the superconductor from the left edge may not be equal to the single-electron charge current coming out from the right edge of the superconductor. We note that the tunneling Hamiltonians (22,37) do not contain any pairing term; thus, the junction current, that is passing between the boundary wire (lead) and the middle wire, is solely carried by single electrons. Nevertheless, Cooper pairs play a significant role in the electrical, thermal, and spin transport inside the superconductors. We discuss below the contribution of the Cooper pairs to the electrical current inside a superconductor.

For devices like N-TS-N and SM-TS-SM, some interesting phenomena can be revealed when we investigate the charge currents inside the superconductors made of either a Kitaev chain or a Majorana wire. Inside a Kitaev chain, the total charge current can be written as a sum of currents carried by (i) single electrons ($J^e$) and (ii) pairs of electrons or Cooper pairs ($J^\Delta$). The latter is the charge current generated due to the motion of the Cooper pairs. It can be demonstrated using the continuity equation for charge density at $l'$-th site of the Y wire made of a Kitaev chain that the charge current going into $l'$-th site is $J^{\mathrm{in}}_{l'}=J^e_{l'-1,l'}-J^\Delta _{l'-1,l'}$, and the charge current coming out from $l'$-th site is $J^{\mathrm{out}}_{l'}=J^e_{l',l'+1}+J^\Delta _{l',l'+1}$, where $J^e_{l',l'+1}$ is the single-electron current from site $l'$ to $l'+1$, and $J^\Delta _{l',l'+1}$ is the Copper pair contribution to the charge current between site $l'$ to $l'+1$. Considering the continuity equation for charge density at $l'$-th site ($\rho _{l'}=c_{l'}^\dagger c_{l'}$), expressions for $J^e_{l',l'+1}$ and $J^\Delta _{l',l'+1}$ are obtained as:

$$\begin{aligned}&J^e_{l',l'+1}=-2e\gamma _{\mathrm{Y}} \mathrm{Im}[ \langle c^\dagger _{l'+1}(t) c_{l'} (t)\rangle ]\, , \nonumber \\&J^\Delta _{l',l'+1}=2e\Delta _{\mathrm{Y}} \mathrm{Im}[ \langle c^\dagger _{l'}(t) c^\dagger _{l'+1} (t)\rangle ]\,, \end{aligned}$$

(42)

where $\gamma _{\mathrm{Y}}$ and $\Delta _{\mathrm{Y}}$ are the hopping and the real pairing potential of the Y wire respectively. From the conservation of total electrical charges at any (e.g., $l'$-th) site inside the TS wire, we have $J^{\mathrm{in}}_{l'}=J^{\mathrm{out}}_{l'}$ for $l'=2,\dots ,L_{\mathrm{Y}}-1$, which implies

$$\begin{aligned} J^e_{l'-1,l'}-J^\Delta _{l'-1,l'}=J^e_{l',l'+1}+J^\Delta _{l',l'+1}. \end{aligned}$$

(43)

Now, $J^e_{\mathrm{XY}}$ is the incoming electrical current for the first site ($l'=1$), where as $J^e_{\mathrm{YZ}}$ is the outgoing electrical current for the last site ($l'=L_{\mathrm{Y}}$). Thus, we further have $J^e_{\mathrm{XY}}=J^{\mathrm{out}}_{1}$ and $J^{\mathrm{in}}_{L_{\mathrm{Y}}}=J^e_{\mathrm{YZ}}$. For notational convenience, we use $J^{\mathrm{in}}_{1}$ and $J^{\mathrm{out}}_{L_{\mathrm{Y}}}$ to represent $J^e_{\mathrm{XY}}$ and $J^e_{\mathrm{YZ}}$ respectively. Now, the conservation of electrical charge yields

$$\begin{aligned} J^e_{\mathrm{XY}}=J^e_{\mathrm{YZ}}+2\sum _{l'=1}^{L_{\mathrm{Y}}-1}J^\Delta _{l',l'+1}, \end{aligned}$$

(44)

which relates the electrical currents at the left and right junction. In Fig. 4, we plot different components of the electrical currents which are carried by single electrons ($J^e_{l',l'+1}$) and Cooper pairs ($J^\Delta _{l',l'+1}$) for each bond $(l',l'+1)$, along with the incoming current ($J^{\mathrm{in}}_{l'}$) and the outgoing current ($J^{\mathrm{out}}_{l'}$) at each lattice sites of the middle wire. We depict $J^e_{l',l'+1}$ and $J^\Delta _{l',l'+1}$ by triangle and diamond symbols respectively, and these are placed at the middle of $l'$-th and $(l'+1)$-th sites as these are associated with $(l',l'+1)$ bond. On the other hand, $J^{\mathrm{in}}_{l'}$ and $J^{\mathrm{out}}_{l'}$ which are respectively represented by circle and ‘+’ symbols, are placed just at the position of $l'$-th lattice site.

4.2 Thermal and Thermoelectric Currents and Conductances

We first discuss the thermal (thermoelectric) current and thermal conductance of a Kitaev chain in an N-TS-N device. In Fig. 14, we show that the values of $J^{u}_{\mathrm{XY}}(t)$ and $J^{u}_{\mathrm{YZ}}(t)$ due to a voltage bias become equal at long times for any initialization of the middle wire which confirms our prediction in Sect. 3.1 based on the conservation of total energy for a TS or an N middle wire. Both for a voltage bias and a thermal bias, the energy transport seems to be ballistic (i.e., $J^{u}_{\mathrm{XY}}$ and $J^{u}_{\mathrm{YZ}}$ are independent of $L_{\mathrm{Y}}$) in an N-TS-N device for the TS wire in a topological phase³. In Table. 1, we show the values of $J^{u}_{\mathrm{XY}}~(=J^{u}_{\mathrm{YZ}})$ for different $L_{\mathrm{Y}}$ in an N-TS-N device with a finite temperature bias (second column) and a finite voltage bias (third column) when the middle TS wire is in a topological phase. The values of $J^{u}_{\mathrm{XY}}~(=J^{u}_{\mathrm{YZ}})$ do not seem to vary with $L_{\mathrm{Y}}$ for longer $L_{\mathrm{Y}}$ within the numerical precision in our study. We note that we take a relatively large bias to achieve better numerical precision in showing the ballistic thermal transport in the middle TS wire’s topological regime. Such a large bias includes the contributions of the above-gap modes in the thermal transport. We also find a ballistic thermal transport near the TS wire’s topological phase transition, as shown in the last column of Table. 1. In comparison to the topological phase, the values of thermal current are relatively higher near the TP (even for a smaller bias), where the thermal conductance rapidly changes with $\epsilon _{\mathrm{Y}}$. We discuss it below.

Table 1 Length dependence of the energy current $J^{u}_{\mathrm{XY}}, J^{u}_{\mathrm{YZ}}$ in an N-TS-N device made of a Kitaev chain

Fig. 5

Thermal conductance $G_{\mathrm{T}}$, zero-frequency thermoelectric transmission ${\mathcal {T}}(\omega =0)$ and thermoelectric conductance $G_{\mathrm{TE}}$ of the middle TS wire as a function of onsite energy $\epsilon _{\mathrm{Y}}$ in an N-TS-N and a TS-TS-N device made of a Kitaev chain. In all panels, $L_{\mathrm{Y}}=100$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Y}}=\gamma _{\mathrm{Z}}=\gamma =1$, $\epsilon _{\mathrm{Z}}=0, \Delta _{\mathrm{Z}}=0$. Further, $\epsilon _{\mathrm{X}}=\Delta _{\mathrm{X}}=0, \Delta _{\mathrm{Y}}=0.3, \gamma _{\mathrm{YZ}}=\gamma _{\mathrm{XY}}=[0.05,0.15,0.25,0.35,0.5,1]$ from left to right curves in panels a, c, $\Delta _{\mathrm{X}}= \Delta _{\mathrm{Y}}=0.3, \gamma _{\mathrm{YZ}}=\gamma _{\mathrm{XY}}=0.15$ in panel b, and $\epsilon _{\mathrm{X}}=\Delta _{\mathrm{X}}=0, \gamma _{\mathrm{YZ}}=\gamma _{\mathrm{XY}}=0.15$ in panel d. All above parameters except lengths are in units of $\gamma$

We notice that the linear-response thermal conductance $G_{\mathrm{T}}$ shows a sharp peak near the TP as $\epsilon _{\mathrm{Y}}$ is tuned through a topological phase transition in the middle TS wire. The value of the sharp peak in $G_{\mathrm{T}}$ at the transition is quantized, and its value is $\pi ^2k_B^2T/6h$, which was predicted earlier in [2, 8, 9]. Such a behavior of $G_{\mathrm{T}}$ across the TP has been proposed for the detection of the topological phase transition in TS wires. Our open quantum-system description of thermal transport confirms the quantized height, which we depict in Fig. 5a. Nevertheless, we further observe that the height of sharp peak in $G_{\mathrm{T}}$ at the transition is sensitive to the tunneling rates $\gamma _{\mathrm{XY}}$ and $\gamma _{\mathrm{YZ}}$ for the open TS wires as shown in Fig. 5a. The value of $\epsilon _{\mathrm{Y}}$, where the sharp peak in $G_{\mathrm{T}}$ appears in the open TS wire, seems to move towards $2\gamma _{\mathrm{Y}}$ (the TP for an isolated TS wire) as $\gamma _{\mathrm{XY}},\gamma _{\mathrm{YZ}}$ approach $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Y}}=\gamma _{\mathrm{Z}}$.

We can also relate the quantized value of $G_{\mathrm{T}}$ to the zero-frequency thermoelectric transmission coefficient (${\mathcal {T}}(\omega =0)$) across the open TS wire under the temperature bias as described in Sect. 3.1. In Fig. 5c, we display $G_{\mathrm{T}}/(\pi ^2k_B^2T/6h)=4\pi {\mathcal {T}}(\omega =0)$ for different tunneling rates. We further observe a large thermoelectric current and a huge dip in the thermoelectric conductance $G_{\mathrm{TE}}=J^{e}_{\mathrm{XY}}/\Delta T=(\pi ^2k_B^2eT/3\hbar ^2)\left[ \frac{ \partial {\mathcal {T}}}{\partial \omega } \right] _{\omega =0}$ near the topological phase transition of the middle TS wire in an N-TS-N device. The height of the dip in $G_{\mathrm{TE}}$ seems to depend on the values of pairing amplitude of the TS wires, and the height increases with a decrease in $\Delta _{\mathrm{Y}}$ as shown in Fig. 5d. The large values of thermoelectricity may generate potential applications of these TS wires. Experimental measurements of such a large dip in the thermoelectric current or conductance in these devices would be much easier than that of the relatively small thermal conductance peak near the phase transition. Therefore, the thermoelectric current or conductance might be a better probe to detect the TS wires’ topological phase transition experimentally.

We inspect the features of $G_{\mathrm{T}}$ in various devices with TS leads to unveil the role of TS leads in transport. In Fig. 5b, we show $G_{\mathrm{T}}$ from the Y-Z junction with $\epsilon _{\mathrm{Y}}$ in a TS-TS-N device by varying $\epsilon _{\mathrm{X}}$ of the TS lead. Surprisingly, we observe a quantized peak in $G_{\mathrm{T}}$ near the TP of the middle TS wire only when the TS lead is also at the TP. When the TS lead is away from the TP, the peak in $G_{\mathrm{T}}$ near the TP of the middle TS wire disappears. Therefore in Fig. 5b, we have plotted just one quantized peak in $G_{\mathrm{T}}$ corresponding to $\epsilon _{\mathrm{X}}=2\gamma _{\mathrm{X}}$. To our opinion, a topological transition of the TS leads can also be probed by the thermal conductance measurements.

4.3 Electrical Current and Differential Conductance in Devices with Kitaev Chain Leads

Different hybrid systems with one or multiple TS leads have been investigated in the recent years to probe emergence of Majorana quasiparticles as well as for efficient quantum devices [26, 36,37,38,39,40, 71, 71]. In Appendix H, we show validity of steady-state transport in TS-N-Z and TS-TS-Z with Z=N and TS at TP by comparing the LEGF results with the direct time-evolution numerics at long time. We observe an interesting electrical current asymmetry in a TS-N-N device with spatial asymmetry (broken parity). The spatial asymmetry in such devices can be engineered by creating different tunneling rates at the X-Y and Y-Z junctions. It is generally expected to have different forward and backward currents under the reversal of bias when there is spatial asymmetry in nonlinear models. Such a difference in currents (rectification) between the forward and reversed bias is generated due to a variation in the distributions of inelastically scattered modes under bias reversal, which is intrinsically related to the nonlinearity. However, the quantum transport for our noninteracting model of N wires and mean-field model of TS wires is expected to be linear. Therefore, a rectification in electrical current generally is not expected in our TS-N-N devices even in the presence of different tunneling rates. Nevertheless, we find a large change in the steady-state electrical currents when we reverse the tunneling rates keeping all other parameters including the bias unaltered. For example, we find $J^e_{\mathrm{XY}}=J^e_{\mathrm{YZ}}=-0.014245$ (in units of $\gamma e$) for $\gamma _{\mathrm{XY}}=0.125, \gamma _{\mathrm{YZ}}=0.25$, and $J^e_{\mathrm{XY}}=J^e_{\mathrm{YZ}}=-0.006628$ (in units of $\gamma e$) for $\gamma _{\mathrm{XY}}=0.25, \gamma _{\mathrm{YZ}}=0.125$, where $L_{\mathrm{Y}}=3$, and the other relevant parameters in the units of $\gamma$ for both cases are $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=\gamma =1, \gamma _{\mathrm{Y}}=0.5, \Delta _{\mathrm{X}}=0.3,\Delta _{\mathrm{Y}}=\Delta _{\mathrm{Z}}=0, \epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Z}}=0,\epsilon _{\mathrm{Y}}=0.05, T_{\mathrm{X}}=T_{\mathrm{Z}}=0.02$ and $\mu _{\mathrm{X}}=0,\mu _{\mathrm{Z}}=0.5$. Such a change in electrical current is related to different hybridization of the Majorana quasiparticle at the TS wires near X-Y junction for different strength of the tunneling rate at that junction. The electrical currents at the both junctions are the same for each set of tunneling rates as expected for a middle N wire. Nevertheless, the electrical currents at the junctions do not change when the bias, e.g., the temperature of the leads, is reversed keeping the tunneling rates fixed at the junctions. This clarifies no true rectification in these hybrid devices with a mean-field model of TS wires.

Fig. 6

Zero-temperature differential electrical conductance $dJ^e_{\mathrm{YZ}}/dV$ at the Y-Z junction of an N-TS-N device (a) and a TS-TS-N device (b), where the TS wires are made of a Kitaev chain. In both panels, $L_{\mathrm{Y}}=50$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=1.0,\gamma _{\mathrm{Y}}=0.5$, $\Delta _{\mathrm{Y}}=0.15, \Delta _{\mathrm{Z}}=0$, $\epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Y}}=\epsilon _{\mathrm{Z}}=0$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.15$, and $\mu _{\mathrm{X}}=0$. Also, $\Delta _{\mathrm{X}}=0$ in a, and $\Delta _{\mathrm{X}}=0.3$ in b. All above parameters except lengths are in units of $\gamma$

Fig. 7

Zero-temperature differential electrical conductance $dJ^e_{\mathrm{YZ}}/dV$ at the Y-Z junction of an N-N-N device (a) and a TS-N-N device made of a Kitaev chain (b). In both panels, $L_{\mathrm{Y}}=50$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=1.0,\gamma _{\mathrm{Y}}=0.5$, $\Delta _{\mathrm{Y}}=\Delta _{\mathrm{Z}}=0$, $\epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Y}}=\epsilon _{\mathrm{Z}}=0$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.15$, and $\mu _{\mathrm{X}}=0$. Also, $\Delta _{\mathrm{X}}=0$ in a, and $\Delta _{\mathrm{X}}=0.3$ in b. All above parameters except lengths are in units of $\gamma$

To identify the unique role of a TS lead compared to an N lead, we also investigate zero-temperature DEC in different devices. Such a zero-temperature DEC has been earlier proposed and applied to detect the emergence of MBS in various junctions of TS wires such as an N-TS (or SM-TS) and an N-TS-N (or SM-TS-SM) devices [3,4,5]. A quantized zero-bias peak of height $2e^2/h$ within the superconducting pairing gap appears in the zero-temperature DEC when the middle TS wire is in a perfect topological phase. The zero-bias peak in the zero-temperature DEC disappears in the topologically trivial phase of the TS wires. Here, we further apply the zero-temperature DEC to quantify TS leads’ role in our different hybrid devices. In Fig. 6, we compare the zero-temperature DEC at the right TS-N junction of a TS-TS-N device to that of an N-TS-N device. While the height of the zero-bias DEC peak is the same in both cases, the height of the DEC peaks above the superconducting pairing gap is $2e^2/h$ for a TS-TS-N device in contrast to $e^2/h$ for an N-TS-N device. This is an intriguing feature as the DEC properties at the Y-Z junction are mostly expected to depend on local properties (e.g., the density of states) of the Y and Z wires for a TS middle wire, which does not conserve the particle number. Therefore, our results indicate that while the features of zero-bias DEC are mostly determined by the local properties of the Y and Z wires, the above-gap DEC peaks are controlled by the properties of both leads (X and Y wires) nonlocally [8]. The height $2e^2/h$ of the above-gap DEC peaks is mainly due to the superconducting pairing of the X wire, which we confirm by keeping the X wire in a topologically trivial phase. For superconducting leads, both electron-type and hole-type quasiparticle excitations contribute to the density of states for the energy spectrum above the superconducting gap. In contrast, for metallic leads, only electrons contribute to the density of states. Moreover, for the energy range (E), which is much greater than the pairing gap of the superconducting leads ($E>>\Delta$), the quasiparticle density of states for the superconductor is almost double of the density of states for the normal metal at that energy E [72]. A larger value of density of states for superconducting lead, which is almost two times that of the metallic bath, is the reason for observing above-gap DEC with an approximate height of around $2e^2/h$. Further, the width of the finite-voltage above-gap DEC peaks is mainly controlled by the couplings $\gamma _{\mathrm{XY}},\gamma _{\mathrm{YZ}}$.

In Fig. 7, we further compare the zero-temperature DEC at the right N-N junction from a TS-N-N device to that of an N-N-N device. We again find that the DEC peaks’ height is $2e^2/h$ for a TS-N-N device in contrast to $e^2/h$ for an N-N-N device. We also notice a weak zero-bias DEC peak in a TS-N-N device whose height is almost $2e^2/h$. The zero-bias peak emerges due to the MBS in the TS wires across the left junction. Therefore, the TS lead’s MBS has a signature in the right link of a coherent device. To clarify the role of the topological phase of the TS leads, we check the zero-temperature DEC in a TS-N-N device when the TS lead is in a topologically trivial phase. We observe that most of the zero-temperature DEC peaks disappear within the left TS lead’s bulk-gap in a trivial phase. The height of the zero-bias peak is also less than $e^2/h$ for the TS lead in a trivial phase.

4.4 Spin Current, Differential Electrical and Spin Conductances in Devices with Majorana Wires

Fig. 8

Comparison of x-component of spin currents ${\mathcal {J}}^{s_x}_{\mathrm{XY}}(t)$ and ${\mathcal {J}}^{s_x}_{\mathrm{YZ}}(t)$ at both junctions obtained from the first-principle/direct time-evolution numerics (full lines) and the generalized LEGF method (dashed lines) for a TS-SM-SM device made of a Majorana wire. The initial density $n_{l',\sigma }=c^{\dagger }_{l',\sigma }c_{l',\sigma }$ with $\sigma =\uparrow ,\downarrow$ of the middle SM wire for the direct time-evolution numerics is indicated on the heading of the top panel. In both panels, $L_{\mathrm{X}}=L_{\mathrm{Z}}=900, L_{\mathrm{Y}}=3$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=\gamma =1$, $\gamma _{\mathrm{Y}}=0.5$, $\Delta _{\mathrm{X}}=0.3,\Delta _{\mathrm{Y}}=\Delta _{\mathrm{Z}}=0, \zeta _{\mathrm{X}}=\zeta _{\mathrm{Y}}=\zeta _{\mathrm{Z}}=\zeta _{\mathrm{XY}}=\zeta _{\mathrm{YZ}}=0.2, B_{\mathrm{X}}=B_{\mathrm{Y}}=B_{\mathrm{Z}}=0.4$, $\epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Z}}=0, \epsilon _{\mathrm{Y}}=0.05$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.25$, $T_{\mathrm{X}}=T_{\mathrm{Z}}=0.02$ and $\mu _{\mathrm{X}}=0,\mu _{\mathrm{Z}}=0.5$. All above parameters except lengths are in units of $\gamma$

We now confirm the validity of steady-state spin current expression derived in the earlier section. For the Majorana wire, the total spin polarization along x-axis is not conserved in the presence of spin-orbit coupling (see Sect. 2.2). Therefore, the x-component of spin current does not remain the same at the left and right junctions of a semiconductor middle wire, i.e., ${\mathcal {J}}^{s_x}_{\mathrm{XY}} \ne {\mathcal {J}}^{s_x}_{\mathrm{YZ}}$. In Fig. 8, we plot the ${\mathcal {J}}^{s_x}_{\mathrm{XY}}(t)$ and ${\mathcal {J}}^{s_x}_{\mathrm{YZ}}(t)$ calculated using the first-principle/direct time-evolution numerics and the generalized LEGF method in a TS-SM-SM device. We find good agreement between the steady-state values of ${\mathcal {J}}^{s_x}_{\mathrm{XY}}=-0.021369$ and ${\mathcal {J}}^{s_x}_{\mathrm{YZ}}= -0.023747$ (in units of $\gamma \hbar /2$) from the generalized LEGF method, and the long-time values of ${\mathcal {J}}^{s_x}_{\mathrm{XY}}(\gamma t=800)=-0.021369$ and ${\mathcal {J}}^{s_x}_{\mathrm{YZ}}(\gamma t=800)=-0.023746$ from the direct time-evolution numerics. We further investigate the steady-state spin current in an SM-TS-SM device, which was experimentally explored recently in [17]. We get from the LEGF, ${\mathcal {J}}^{s_x}_{\mathrm{XY}}=0.004686$ and ${\mathcal {J}}^{s_x}_{\mathrm{YZ}}=-0.01072$ (in units of $\gamma \hbar /2$) for an SM-TS-SM device with $L_{\mathrm{Y}}=10$, $\gamma _{\mathrm{Y}}=1$, $\Delta _{\mathrm{X}}=0, \Delta _{\mathrm{Y}}=0.3$, and $\mu _{\mathrm{X}}=0.2,\mu _{\mathrm{Z}}=0.5$ (other parameters are the same as Fig. 8). The long-time values from the direct time-evolution numerics are ${\mathcal {J}}^{s_x}_{\mathrm{XY}}(\gamma t=1475)=0.004592$ and ${\mathcal {J}}^{s_x}_{\mathrm{YZ}}(\gamma t=1475)=-0.01047$, which show some deviations between the two methods for longer $L_{\mathrm{Y}}$.

Fig. 9

Zero-temperature differential electrical conductance $d{\mathcal {J}}^e_{\mathrm{YZ}}/dV$ at the Y-Z junction of an SM-TS-SM device (a), and a TS-TS-SM device (b), where the TS wires are the Majorana wire. In both panels, $L_{\mathrm{Y}}=40$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=\gamma _{\mathrm{Y}}=1.0$, $\Delta _{\mathrm{Y}}=0.3, \Delta _{\mathrm{Z}}=0, \zeta _{\mathrm{X}}=\zeta _{\mathrm{Y}}=\zeta _{\mathrm{Z}}=\zeta _{\mathrm{XY}}=\zeta _{\mathrm{YZ}}=0.2, B_{\mathrm{X}}=B_{\mathrm{Y}}=B_{\mathrm{Z}}=0.4$, $\epsilon _{\mathrm{Y}}=0,\epsilon _{\mathrm{Z}}=1$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.25$, and $\mu _{\mathrm{X}}=0$. Also, $\Delta _{\mathrm{X}}=0, \epsilon _{\mathrm{X}}=1$ in a, and $\Delta _{\mathrm{X}}=0.3, \epsilon _{\mathrm{X}}=0$ in b. All above parameters except lengths are in units of $\gamma$

Next, we clarify the role of TS leads made of Majorana wires. For this, we again compare the zero-temperature DEC in a TS-SM-SM device to an SM-SM-SM device and a TS-TS-SM device to an SM-TS-SM device. In Fig. 9, we show the zero-temperature DEC in TS-TS-SM and SM-TS-SM devices where the TS wires are kept in the topological phase. In contrast to the Kitaev chains, the difference in the finite bias DEC height above the bulk-gap between a TS-TS-SM device and an SM-TS-SM device is relatively less. It is probably due to the more structured density of states of the spin-orbit coupled wires. Nevertheless, there is a clear zero-bias DEC peak of height $2e^2/h$ in a TS-SM-SM device, which is absent in an SM-SM-SM device. We show it in Fig. 10. Therefore, the TS leads in the topological phase do have a signature in transport even for the spinful model of TSs.

Fig. 10

Zero-temperature differential electrical conductance $d{\mathcal {J}}^e_{\mathrm{YZ}}/dV$ at the Y-Z junction of an SM-SM-SM device (a), and a TS-SM-SM device (b), where the TS wires are the Majorana wire. In both panels, $L_{\mathrm{Y}}=40$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=\gamma _{\mathrm{Y}}=1.0$, $\Delta _{\mathrm{Y}}=\Delta _{\mathrm{Z}}=0, \zeta _{\mathrm{X}}=\zeta _{\mathrm{Y}}=\zeta _{\mathrm{Z}}=\zeta _{\mathrm{XY}}=\zeta _{\mathrm{YZ}}=0.2, B_{\mathrm{X}}=B_{\mathrm{Y}}=B_{\mathrm{Z}}=0.4$, $\epsilon _{\mathrm{Y}}=\epsilon _{\mathrm{Z}}=1$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.25$, and $\mu _{\mathrm{X}}=0$. Also, $\Delta _{\mathrm{X}}=0, \epsilon _{\mathrm{X}}=1$ in a, and $\Delta _{\mathrm{X}}=0.3, \epsilon _{\mathrm{X}}=0$ in b. All above parameters except lengths are in units of $\gamma$

In Sect. 3.2, we introduce a definition of DSC for spinful TS wire junctions. We here discuss the features of zero-temperature DSC in SM-TS-SM, TS-TS-SM, and TS-SM-SM devices. We show in Fig. 11b, c that the value of zero-bias DSC peak at zero temperature is also quantized in the above devices, while Y wire is in the topological phase. The quantized value of $d{\mathcal {J}}^{s_x}_{\mathrm{YZ}}/dV$ at zero bias is 2 in the unit of $(e/h)(\hbar /2) = (e/4 \pi )$. The quantized value of zero-bias DSC reminds us that of DEC in such devices. We note that the other components of DSC (e.g., y and z components) are not quantized in our devices. We observe the quantization of DSC is mainly along the direction of polarization of the electrons in the boundary wires (e.g., Z wire for the DSC at Y-Z junction) [73]. The quantization of zero-bias DSC at zero temperature is a consequence of Majorana zero modes that developed at the edges of the middle TS wire. It is a local phenomenon as it depends mainly on the localized Majorana zero mode and its tunnel coupling to the boundary wire. We have verified that the quantization of the zero-temperature zero-bias peak height of DSC is robust against disorder in onsite energies of the middle TS wire (check Fig. 11) and the strength of the tunneling rates. It should be noted that neither the electrical charge nor the total spin components is conserved inside a topological superconductor like Majorana wire (TS with spin), however zero-bias peaks in both DSC and DEC signify the presence of Majorana zero modes at the edges of the Majorana wire.

For comparison, we further depict the zero-temperature DSC in an SM-SM-SM device in Fig. 11a, which does not show a quantized zero-bias peak of height 2. Therefore, the quantized DSC along with the DEC can be employed to probe the emergence of MBSs in engineered semiconductor-superconductor heterostructures. In Fig. 11d, we also notice a sharp zero-bias peak in zero-temperature DSC from a TS-SM-SM device when the boundary X wire is in the topological phase, and the height of the zero-bias peak is almost quantized to 2. Nevertheless, there are finite DSC peaks within the superconducting gap of the boundary X wire in Fig. 11d, and these peaks are due to the middle semiconductor wire as also seen for zero-temperature DEC in Fig. 10b.

Fig. 11

Zero-temperature differential spin conductance $d{\mathcal {J}}^{s_x}_{\mathrm{YZ}}/dV$ at the Y-Z junction of an SM-SM-SM device (a), an SM-TS-SM device (b), a TS-TS-SM device (c), and a TS-SM-SM device (d), where the TS wires are the Majorana wire. The middle Y wire in all the devices has an onsite disorder with a uniform distribution of zero mean and variance 0.4, and we show the DSC signals after disorder averaging over 20 realizations. In all panels, $L_{\mathrm{Y}}=40$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=\gamma _{\mathrm{Y}}=1.0$, $\Delta _{\mathrm{Z}}=0, \zeta _{\mathrm{X}}=\zeta _{\mathrm{Y}}=\zeta _{\mathrm{Z}}=\zeta _{\mathrm{XY}}=\zeta _{\mathrm{YZ}}=0.2, B_{\mathrm{X}}=B_{\mathrm{Y}}=B_{\mathrm{Z}}=0.4$, $\epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Y}}=\epsilon _{\mathrm{Z}}=0$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.25$, and $\mu _{\mathrm{X}}=0$. Also, $\Delta _{\mathrm{X}}=\Delta _{\mathrm{Y}}=0$ in a, $\Delta _{\mathrm{X}}=0,\Delta _{\mathrm{Y}}=0.3$ in b, $\Delta _{\mathrm{X}}=\Delta _{\mathrm{Y}}=0.3,$ in c, $\Delta _{\mathrm{X}}=0.3,\Delta _{\mathrm{Y}}=0$ in d. All above parameters except lengths are in units of $\gamma$

5 Conclusion

This paper develops a unified open-quantum system description of nonequilibrium electrical, thermal, and spin transport in various devices whose systems and leads/baths are made of different TS wires. We demonstrate that the quantum LEGF method suits perfectly to derive neat expressions of different steady-state currents. We here notably extend the applications of the LEGF method to the baths of TS wires. We mainly try to reveal several significant thermal and spin transport features in different junctions of TS, SM, and N wires, which are either less explored or have not been investigated. While the spin transport in the junctions of TS and SM wires are rarely investigated in the past, a systematic study of spin transport can be fruitful in disclosing necessary signatures of Majorana quasiparticles (such as spin polarizations of Majoranas [14]) and the topological nature of the systems [17, 46, 47]. The obtained expression of spin current here, is expected to be useful for further investigation of spin transport. We hope the quantized zero-bias peak in zero-temperature DSC would be tested in engineered spinful TS wire junctions [17] for experimentally detecting topological phases of Majorana wires. Another interesting quantity that may be studied in this framework is the local compressibility which shows divergence at topological phase transition and can be an interesting probe to detect the same [74]. This compressibility while combined with thermal transport may reveal new Majorana quantization [75].

The open-quantum system description of transport is more appropriate in determining the quantum materials’ topological signatures as it incorporates the bath-induced dissipation in the topological materials (and generates an effective non-Hermitian notion) akin to the engineered TS devices in experiments [3,4,5,6,7, 17]. Such a description is also required to explain the discrepancy between the experimentally measured and theoretically proposed (using transport in isolated systems) height of the zero-bias DEC manifesting emergence of the MBSs and the characteristics of the quantized peak in thermal conductance indicating the topological phase transition. For example, the position of the quantized peak in thermal conductance as a function of Fermi energy of the middle TS wires shifts with the system-bath coupling.

The expressions of electrical, thermal, and spin current, can be further explored to investigate thermoelectric, magnetoelectric, and thermomagnetic transport properties, especially in the linear response regime. Such studies may show further quantization of different unexplored transport quantities.

The LEGF method suits well to extend the above calculation to find current fluctuation across its average value due to thermal and quantum noises. For this, we need to obtain two-point correlators of current operators at the different instants of time, which can be analyzed to investigate correlations of the nonequilibrium current (shot noise) at the same or different junctions [28, 44, 45] due to the quantization of the charge and energy of the carriers.

Appendices

Appendix A: Quantum Langevin Equations for Y Wire Connected to the Kitaev Chain Leads

The Eqs. 24 and 25 in the main text can be solved using single-particle retarded Green’s function of the isolated $\alpha$ lead/bath wire ($\alpha = X, Z$):

$$\begin{aligned} G^+_\alpha (\tau )=-i e^{-i{\mathcal {H}}^{\alpha } \tau } \theta (\tau ), \end{aligned}$$

(A1)

where both $G^+_\alpha$ and ${\mathcal {H}}^{\alpha }$ are matrices, and $\theta (\tau )$ is the Heaviside step function. This matrix Green’s function is related to the operator valued Green’s function as ${\mathcal {G}}^+_\alpha (\tau )=\sum _{l,m}[{G}^+_\alpha (\tau )]_{l,m}a^{\dagger }_la_m$, where $[{G}^+_\alpha (\tau )]_{l,m} \equiv [- i e^{-i\mathcal {H}^{\alpha } \tau } ]_{lm} \theta (\tau )$. Here, $l,m= 1,\dots , 2L_{\mathrm{X}}$ for $\alpha =X$, and $l,m=2L_{\mathrm{XY}}+1, \dots , 2L$ for $\alpha =Z$. These bath Green’s functions are the solution of the following equation:

$$\begin{aligned} \left( i \frac{\partial }{\partial t} -{\mathcal {H}}^{\alpha } \right) G_{\alpha }^+(t -t')= \delta (t-t')\mathbb {1}_{2L_{\alpha }}. \end{aligned}$$

(A2)

Using these Green’s functions (A1), we formally solve the Heisenberg equations of the boundary wires in Eqs. 24, 25, and the solutions are for $t>t_0$:

$$\begin{aligned} a_{k}(t)= & {} i \sum _{l=1}^{2L_{\mathrm{X}}} {[G_{\mathrm{X}}^+(t -t_0)]}_{k,l} \, a_{l}(t_0) -\gamma _{\mathrm{XY}} \int _{t_0}^{t} dt' \sum _{l=2L_{\mathrm{X}}-1}^{2L_{\mathrm{X}}} (-1)^{l+1} {[G_{\mathrm{X}}^+(t -t')]}_{k,l} a_{l+2}(t') \, , \nonumber \\ a_{k'}(t)= & {} i \sum _{l=2L_{\mathrm{XY}}+1}^{2L} {[G_{\mathrm{Z}}^+(t -t_0)]}_{k'-2L_{\mathrm{XY}},l-2L_{\mathrm{XY}}} \, a_{l}(t_0) \nonumber \\&-\gamma _{\mathrm{YZ}} \int _{t_0}^{t} dt' \sum _{l=1}^2 (-1)^{l+1} {[G_{\mathrm{Z}}^+(t -t')]}_{k'-2L_{\mathrm{XY}},l} \, a_{2 L_{\mathrm{XY}}-2+l}(t') \, , \end{aligned}$$

(A3)

where $k=2L_{\mathrm{X}}-1,2L_{\mathrm{X}}$ and $k'=2L_{\mathrm{XY}}+1 ,2L_{\mathrm{XY}}+2$. $G_{\mathrm{X}}^+(t -t_0)$ and $G_{\mathrm{Z}}^+(t -t_0)$ are the single-particle retarded Green’s functions of the X and Z baths, respectively. Now, plugging the formal solutions of the X wire variables (A3) into Eq. 23, we get the quantum Langevin equation of Y wire variables at the X-Y junction:

$$\begin{aligned} {\dot{a}}_{k}(t)= -i \sum _{l=2L_{\mathrm{X}}+1}^{2L_{\mathrm{XY}}} \mathcal {H}^\mathrm{Y}_{k,l} \, a_l(t) -i \eta _{k-2L_{\mathrm{X}}}^{\mathrm{X}}(t) - i \int _{t_0}^t dt' \sum _{l=2L_{\mathrm{X}}+1}^{2L_{\mathrm{X}}+2}{[\Sigma _{\mathrm{X}}^+(t -t')]}_{k,l}\, a_{l}(t'),\nonumber \\ \end{aligned}$$

(A4)

for $k=2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2$, where

$$\begin{aligned} \eta _1^{\mathrm{X}}(t)= & {} -i\gamma _{\mathrm{XY}} \sum _{l=1}^{2L_{\mathrm{X}}}{[G_{\mathrm{X}}^+(t -t_0)]}_{2L_{\mathrm{X}}-1,l} \, a_{l}(t_0), \nonumber \\ \eta _2^{\mathrm{X}}(t)= & {} i\gamma _{\mathrm{XY}} \sum _{l=1}^{2L_{\mathrm{X}}}{[G_{\mathrm{X}}^+(t -t_0)]}_{2L_{\mathrm{X}},l} \, a_{l}(t_0), \end{aligned}$$

(A5)

are the noise terms from the X bath, and the dissipative terms generated in the Y wire due to the coupling to the X bath are:

$$\begin{aligned}{}[\Sigma _{\mathrm{X}}^+(t)]_{l,m}= & {} \gamma _{\mathrm{XY}}^2\sum _{k=2L_{\mathrm{X}}-1}^{2L_{\mathrm{X}}}[G_{\mathrm{X}}^+(t)]_{k,k} \, \delta _{l,k+2}\, \delta _{m,l} \nonumber \\&-\gamma _{\mathrm{XY}}^2 \sum _{\begin{array}{c} k,k'=2L_{\mathrm{X}}-1 \\ k \ne k' \end{array}}^{2L_{\mathrm{X}}} [G_{\mathrm{X}}^+(t)]_{k,k'} \delta _{l,k+2}\delta _{m,k'+2} . \end{aligned}$$

(A6)

Similarly, the quantum Langevin equations for the Y wire variables at Y-Z junction are:

$$\begin{aligned} {\dot{a}}_{k}(t)= & {} -i \sum _{j=2L_{\mathrm{X}}+1}^{2L_{\mathrm{XY}}} {\mathcal {H}}^{\mathrm{Y}}_{k,j} \, a_j(t) -i \eta _{k-2L_{\mathrm{X}}+2}^{\mathrm{Z}}(t) - i \int _{t_0}^t dt' \sum _{l=2L_{\mathrm{XY}}-1}^{2L_{\mathrm{XY}}} {[\Sigma _{\mathrm{Z}}^+(t -t')]}_{k,l}\, a_{l}(t'),\nonumber \\ \end{aligned}$$

(A7)

for $k=2L_{\mathrm{XY}}-1,2L_{\mathrm{XY}}$, where

$$\begin{aligned} \eta _{1}^{\mathrm{Z}}(t)= & {} -i\gamma _{\mathrm{YZ}} \sum _{l=2L_{\mathrm{XY}}+1}^{2L} {[G_{\mathrm{Z}}^+(t -t_0)]}_{1,l-2L_{\mathrm{XY}}} \, a_{l}(t_0), \nonumber \\ \eta _{2}^{\mathrm{Z}}(t)= & {} i\gamma _{\mathrm{YZ}} \sum _{l=2L_{\mathrm{XY}}+1}^{2L} {[G_{\mathrm{Z}}^+(t -t_0)]}_{2,l-2L_{\mathrm{XY}}} \, a_{l}(t_0), \end{aligned}$$

(A8)

are the noises from the Z bath, and

$$\begin{aligned}{}[\Sigma _{\mathrm{Z}}^+(t)]_{l,m}= & {} \gamma _{\mathrm{YZ}}^2 \sum _{k=1}^{2}[G_{\mathrm{Z}}^+(t)]_{k,k} \, \delta _{l,2L_{\mathrm{XY}}-2+k}\delta _{m,l} \nonumber \\&-\gamma _{\mathrm{YZ}}^2 \sum _{\begin{array}{c} k,k'=1 \\ k \ne k' \end{array}}^{2}[G_{\mathrm{Z}}^+(t)]_{k,k'} \, \delta _{l,2L_{\mathrm{XY}}-2+k}\delta _{m,2L_{\mathrm{XY}}-2+k'}, \end{aligned}$$

(A9)

are the dissipative terms due to the Z bath. For the internal sites of the Y wire, the Heisenberg equations take the following form:

$$\begin{aligned} {\dot{a}}_{l}(t)=-i \sum _{j=2L_{\mathrm{X}}+1}^{2L_{\mathrm{XY}}} {\mathcal {H}}^{\mathrm{Y}}_{l,j} \, a_j(t), \end{aligned}$$

(A10)

for $l=2L_{\mathrm{X}}+3,\dots ,2L_{\mathrm{XY}}-2$. We finally rewrite Eqs. A4, A7, and A10 in a compact form:

$$\begin{aligned} {\dot{a}}_{l}(t)= & {} -i \sum _{m=2L_{\mathrm{X}}+1}^{2L_{\mathrm{XY}}} {\mathcal {H}}^{\mathrm{Y}}_{l,m} a_m(t) -i \sum _{k=1}^2 \left( \delta _{l,2L_X+k}\, \eta _{k}^{\mathrm{X}}(t)+\delta _{l,2L_{\mathrm{XY}}-2+k}\, \eta _{k}^{\mathrm{Z}}(t) \right) \nonumber \\&-i(\delta _{l,2L_X+1}+\delta _{l,2L_X+2})\int _{t_0}^t dt' \sum _{m=2L_X+1}^{2L_X+2} {[\Sigma _{\mathrm{X}}^+(t -t')]}_{l,m}\, a_{m}(t') \nonumber \\&-i(\delta _{l,2L_{\mathrm{XY}}-1}+\delta _{l,2L_{\mathrm{XY}}})\int _{t_0}^t dt' \sum _{m=2L_{\mathrm{XY}}-1}^{2L_{\mathrm{XY}}} {[\Sigma _{\mathrm{Z}}^+(t -t')]}_{l,m}\, a_{m}(t') , \end{aligned}$$

(A11)

where, $l \in [2L_{\mathrm{X}}+1,2L_{\mathrm{XY}}]$. The Eq. A11 is in a form of generalized quantum Langevin equations [54]. The Eq. A11 can be solved analytically by using the Fourier transformation for a particular class of systems, which have a unique NESS. Nevertheless, to apply the Fourier transformation on Eq. A11, we first need to define the Fourier transformation of the noise, and the dissipative/self-energy terms: ${\tilde{\eta }}_{1,2}^{\alpha }(\omega )=\frac{1}{2\pi }\int _{-\infty }^{\infty } dt \,\eta _{1,2}^{\alpha }(t)\, e^{i \omega t}$ and ${\tilde{\Sigma }}_{\alpha }^+(\omega )=\int _{-\infty }^{\infty } dt \, \Sigma _{\alpha }^+(t)\, e^{i \omega t}$ for $\alpha =X,Z$. We note that the steady-state solutions for the Y wire variables in frequency domain are given in Eq. 26 of the main text.

Appendix B: Green’s Function of an Isolated Kitaev Chain and Its Properties

Let us consider an isolated $\alpha$ wire lead $(\alpha =X,Z)$ of the Kitaev chain represented by the Hamiltonian $H^{\alpha }$. The r-th eigenvalue of $H^{\alpha }$ and the corresponding eigenfunction are respectively denoted by $\hbar \omega _r^\alpha$ and $U_r^\alpha = [ \phi _r^\alpha (1), \psi _r^\alpha (1), \dots , \phi _r^\alpha (L ), \psi _r^\alpha (L )]^T$. As we have already discussed, the Hamiltonian matrix ${\mathcal {H}}^\alpha$ for the isolated $\alpha$ wire can be expressed as,

$$\begin{aligned} {\mathcal {H}}^\alpha = \sum _{r>0} \omega _r^\alpha (U^\alpha _r U_r^{\alpha \,\dagger }- V^\alpha _r V_r^{\alpha \, \dagger }). \end{aligned}$$

(B1)

In Appendix A, we have introduced the single-particle retarded Green’s function of $\alpha$ wire using Eq. A2. Next, we rewrite the bath Green’s function by using the form of matrix ${\mathcal {H}}^\alpha$ as

$$\begin{aligned} G^+_{\alpha }(t-t')=-i \theta (t-t') \sum _{r>0}\left( U^\alpha _r U_r^{\alpha \,\dagger } e^{-i \omega _r^\alpha (t-t')} +V^\alpha _r V_r^{\alpha \, \dagger } e^{i \omega _r^\alpha (t-t')} \right) . \end{aligned}$$

(B2)

It is easy to check the following properties of the Green’s function:

$$\begin{aligned}&{[G_\alpha ^+(\tau )]}_{2i',2j'}=-{[G_\alpha ^+(\tau )]}_{2i'-1,2j'-1}^*\, , \nonumber \\&{[G_\alpha ^+(\tau )]}_{2i',2j'-1}=-{[G_\alpha ^+(\tau )]}_{2i'-1,2j'}^* \, , \end{aligned}$$

(B3)

for $i',j'= 1,\dots , L_\alpha$. The Fourier transformation of Green’s function is defined as: ${{\tilde{G}}^{\alpha +}_{l,m}(\omega )}=\frac{1}{2 \pi }\int _{-\infty }^{\infty } d\tau \, { [G_\alpha ^+(\tau )]}_{l,m}\, e^{i \omega \tau }$. Hence, the retarded Green’s function in the frequency domain satisfies the following relation:

$$\begin{aligned}&{{\tilde{G}}^{\alpha +}_{2i',2j'}(\omega )}=-{[{\tilde{G}}^{\alpha +}_{2i'-1,2j'-1}(-\omega )]}^*, \nonumber \\&{{\tilde{G}}^{\alpha +}_{2i',2j'-1}(\omega )}=-{[{\tilde{G}}^{\alpha +}_{2i'-1,2j'}(-\omega )]}^*\, . \end{aligned}$$

(B4)

In the frequency domain, components of the bath Green’s functions are expressed as

$$\begin{aligned} {{\tilde{G}}^{\alpha +}_{2i'-1,2j'-1}(\omega )}= & {} \sum _{r>0} \left[ \frac{\phi _r^\alpha (i')\phi ^{\alpha *}_r(j')}{\omega -\omega _r^\alpha } +\frac{\psi ^{\alpha *}_r(i')\psi _r^\alpha (j')}{\omega +\omega _r^\alpha } \right] \nonumber \\&- i \pi \sum _{r>0} \left[ \phi _r^\alpha (i')\phi ^{\alpha *}_r(j')\delta ( \omega -\omega _r^\alpha ) + \psi ^{\alpha *}_r(i')\psi _r^\alpha (j') \delta ( \omega +\omega _r^\alpha )\right] , \nonumber \\ {{\tilde{G}}^{\alpha +}_{2i',2j'}(\omega )}= & {} \sum _{r>0} \left[ \frac{\psi _r^\alpha (i')\psi ^{\alpha *}_r(j')}{\omega -\omega _r^\alpha } +\frac{\phi ^{\alpha *}_r(i')\phi _r^\alpha (j')}{\omega +\omega _r^\alpha } \right] \nonumber \\&- i \pi \sum _{r>0} \left[ \psi _r^\alpha (i')\psi ^{\alpha *}_r(j')\delta ( \omega -\omega _r^\alpha ) + \phi ^{\alpha *}_r(i')\phi _r^\alpha (j') \delta ( \omega +\omega _r^\alpha )\right] , \nonumber \\ {{\tilde{G}}^{\alpha +}_{2i',2j'-1}(\omega )}= & {} \sum _{r>0} \left[ \frac{\psi _r^\alpha (i')\phi ^{\alpha *}_r(j')}{\omega -\omega _r^\alpha } +\frac{\phi ^{\alpha *}_r(i')\psi _r^\alpha (j')}{\omega +\omega _r^\alpha } \right] \nonumber \\&- i \pi \sum _{r>0} \left[ \psi _r^\alpha (i')\phi ^{\alpha *}_r(j')\delta ( \omega -\omega _r^\alpha ) + \phi ^{\alpha *}_r(i')\psi _r^\alpha (j')\delta ( \omega +\omega _r^\alpha )\right] , \nonumber \\ {{\tilde{G}}^{\alpha +}_{2i'-1,2j'}(\omega )}= & {} \sum _{r>0} \left[ \frac{\phi _r^\alpha (i')\psi ^{\alpha *}_r(j')}{\omega -\omega _r^\alpha } +\frac{\psi ^{\alpha *}_r(i')\phi _r^\alpha (j')}{\omega +\omega _r^\alpha } \right] \nonumber \\&- i \pi \sum _{r>0} \left[ \phi _r^\alpha (i')\psi ^{\alpha *}_r(j')\delta ( \omega -\omega _r^\alpha ) + \psi ^{\alpha *}_r(i')\phi _r^\alpha (j')\delta ( \omega +\omega _r^\alpha )\right] .\nonumber \\ \end{aligned}$$

(B5)

We observe that $\{ \phi _r^\alpha (i') \}$ and $\{ \psi _r^\alpha (j') \}$ are real for a real superconducting gap $\Delta _\alpha$. Thus, we can write

$$\begin{aligned} \text {Im}[ {\tilde{G}}^{\alpha +}_{2i'-1,2j'-1}(\omega )]= & {} -\pi \sum _{r>0} \left[ \phi _r^\alpha (i')\phi ^{\alpha *}_r(j')\delta ( \omega -\omega _r^\alpha ) + \psi ^{\alpha *}_r(i')\psi _r^\alpha (j') \delta ( \omega +\omega _r^\alpha )\right] ,\nonumber \\ \text {Im}[ {\tilde{G}}^{\alpha +}_{2i',2j'}(\omega )]= & {} -\pi \sum _{r>0} \left[ \psi _r^\alpha (i')\psi ^{\alpha *}_r(j') \delta ( \omega -\omega _r^\alpha ) + \phi ^{\alpha *}_r(i')\phi _r^\alpha (j') \delta ( \omega +\omega _r^\alpha )\right] , \nonumber \\ \text {Im}[{\tilde{G}}^{\alpha +}_{2i',2j'-1}(\omega ) ]= & {} -\pi \sum _{r>0} \left[ \psi _r^\alpha (i')\phi ^{\alpha *}_r(j') \delta ( \omega -\omega _r^\alpha ) + \phi ^{\alpha *}_r(i')\psi _r^\alpha (j') \delta ( \omega +\omega _r^\alpha )\right] , \nonumber \\ \text {Im}[{\tilde{G}}^{\alpha +}_{2i'-1,2j'}(\omega )]= & {} - \pi \sum _{r>0} \left[ \phi _r^\alpha (i')\psi ^{\alpha *}_r(j') \delta ( \omega -\omega _r^\alpha ) + \psi ^{\alpha *}_r(i')\phi _r^\alpha (j') \delta ( \omega +\omega _r^\alpha )\right] .\nonumber \\ \end{aligned}$$

(B6)

The above relations (B6) are of significant importance, and we subsequently use them in deriving the noise-noise correlations of superconducting leads with a real superconducting gap.

Appendix C: Noise-Noise Correlations for Kitaev Chain Leads

Employing the Eqs. 11 and 12, we calculate the equilibrium correlation functions for a Kitaev chain lead. For such an isolated $\alpha$ wire bath ($\alpha = X, Z$), the equilibrium correlations are given by,

$$\begin{aligned} \langle a_{2i'-1}^\dagger \, a_{2j'-1} \rangle= & {} \langle c_{i'}^{\dagger }\, c_{j'} \rangle = \sum _{r>0}\left[ \phi ^{\alpha *}_r(i') \phi ^\alpha _r(j') f(\omega _r^\alpha , T_\alpha )+ \psi _r^\alpha (i') \psi _r^{\alpha *}(j')f(-\omega _r^\alpha , T_\alpha )\right] , \nonumber \\ \langle a_{2i'-1} \, a_{2j'-1}^\dagger \rangle= & {} \langle c_{i'}\, c_{j'}^\dagger \rangle = \sum _{r>0}\left[ \psi _r^{\alpha *}(i') \psi _r^\alpha (j') f(\omega _r^\alpha , T_\alpha )+\phi _r^\alpha (i') \phi _r^{\alpha *}(j') f(-\omega _r^\alpha ,T_\alpha )\right] , \nonumber \\ \langle a_{2i'-1} \, a_{2j'-1} \rangle= & {} \langle c_{i'} \, c_{j'} \rangle = \sum _{r>0}\left[ \psi _r^{\alpha *}(i') \phi _r^\alpha (j')f(\omega _r^\alpha ,T_\alpha )+\phi _r^\alpha (i') \psi _r^{\alpha *}(j')f(-\omega _r^\alpha ,T_\alpha )\right] , \nonumber \\ \langle a_{2i'-1}^\dagger \, a_{2j'-1}^\dagger \rangle= & {} \langle c_{i'}^{ \dagger } \, c_{j'}^{ \dagger } \rangle = \sum _{r>0}\left[ \phi _r^{\alpha *}(i') \psi _r^\alpha (j')f(\omega _r^\alpha , T_\alpha )+\psi _r^\alpha (i') \phi _r^{\alpha *}(j')f(-\omega _r^\alpha , T_\alpha )\right] ,\nonumber \\ \end{aligned}$$

(C1)

where $i',j'=1,\dots ,L_{\alpha }$. Here, $T_\alpha$ is the temperature of the $\alpha$ wire bath, and $\hbar \omega _r^\alpha$ is the r-th eigenvalue, and $\{ \phi _r^{\alpha },\psi _r^{\alpha }\}$ are the components of corresponding eigenfunction. Moreover, these relations are valid only when the baths are kept at zero chemical potential. The other possible two-point correlations between $a_l$ variables can easily be derived employing the relation: $a_{2i'}=a_{2i'-1}^\dagger$. Using the definitions of noises (A5), we now find the noise-noise correlation in time domain:

$$\begin{aligned} \langle \eta _1^{\mathrm{X} \dagger }(t) \eta _1^{\mathrm{X}}(t')\rangle= & {} \gamma _{\mathrm{XY}}^2 \theta (t-t_0)\theta (t'-t_0) \sum _{r>0} \left( \phi _r^{\mathrm{X}}(L_{\mathrm{X}}) \phi _r^{\mathrm{X *}}(L_{\mathrm{X}}) e^{i \omega _r^{\mathrm{X}} (t-t')} f(\omega _r^{\mathrm{X}},T_{\mathrm{X}}) \right. \nonumber \\&\left. +\psi _r^{\mathrm{X} *}(L_{\mathrm{X}}) \psi _r^{\mathrm{X}}(L_{\mathrm{X}}) e^{-i \omega _r^{\mathrm{X}}(t-t')} f(-\omega _r^{\mathrm{X}},T_{\mathrm{X}}) \right) . \end{aligned}$$

(C2)

The properties of the noises can be written in a convenient form in the frequency domain. Thus, we now convert the above noise-noise correlation to frequency domain, and we take the Fourier transformation (by letting $t_0 \rightarrow -\infty$)

$$\begin{aligned} \langle {\tilde{\eta }}_1^{\mathrm{X \dagger }}(\omega ) {\tilde{\eta }}_1^{\mathrm{X}}(\omega ')\rangle= & {} -{\left( \frac{1}{2 \pi }\right) }^2\int _{-\infty }^{\infty } dt e^{-i \omega t} \int _{-\infty }^{\infty } dt' e^{i \omega ' t'}\langle \eta _1^{\mathrm{X \dagger }}(t) \eta _1^{\mathrm{X}}(t')\rangle \nonumber \\= & {} \gamma _{\mathrm{XY}}^2 \, \delta (\omega -\omega ')\, \sum _{r>0} \left( \phi _r^{\mathrm{X}}(L_{\mathrm{X}}) \phi _r^{\mathrm{X*}}(L_{\mathrm{X}})\delta (\omega -\omega _r^{\mathrm{X}})f(\omega ,T_{\mathrm{X}}) \right. \nonumber \\&\left. +\psi _r^{\mathrm{X*}}(L_{\mathrm{X}}) \psi _r^{\mathrm{X}}(L_{\mathrm{X}}) \delta (\omega +\omega _r^{\mathrm{X}})f(\omega ,T_{\mathrm{X}}) \right) \nonumber \\= & {} - \frac{\gamma _{\mathrm{XY}}^2}{\pi } \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}}-1,2L_{\mathrm{X}}-1}(\omega )] f(\omega ,T_{\mathrm{X}}) \, \delta _{\omega ,\omega '}\, , \end{aligned}$$

(C3)

where, $\delta _{\omega ,\omega '} \equiv \delta (\omega -\omega ')$. We have applied (B6) in the last line of the above expression, and this substitution is valid when $\{ \phi _r^{\mathrm{X}}(i'),\psi _r^{\mathrm{X}}(j')\}$ are real, which is true for the N and the Kitaev chain Hamiltonian with a real superconducting gap ($\Delta _{\mathrm{X}}$). We similarly express all other noise-noise correlations from the X bath (an N or a Kitaev chain with real $\Delta _{\mathrm{X}}$) in the compact forms by applying (B6):

$$\begin{aligned} \langle {\tilde{\eta }}_2^{\mathrm{X \dagger }}(\omega ) {\tilde{\eta }}_2^{\mathrm{X}}(\omega ')\rangle= & {} - \frac{\gamma _{\mathrm{XY}}^2}{\pi } \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}},2L_{\mathrm{X}}}(\omega )]f(\omega , T_{\mathrm{X}})\delta _{\omega ,\omega '} \, , \nonumber \\ \langle {\tilde{\eta }}_1^{\mathrm{X \dagger }}(\omega ) {\tilde{\eta }}_2^{\mathrm{X}}(\omega ')\rangle= & {} \frac{\gamma _{\mathrm{XY}}^2}{\pi } \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}},2L_{\mathrm{X}}-1}(\omega )]f(\omega , T_{\mathrm{X}})\delta _{\omega ,\omega '}\, , \nonumber \\ \langle {\tilde{\eta }}_2^{\mathrm{X \dagger }}(\omega ) {\tilde{\eta }}_1^{\mathrm{X}}(\omega ')\rangle= & {} \frac{\gamma _{\mathrm{XY}}^2}{\pi } \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}}-1,2L_{\mathrm{X}}}(\omega )]f(\omega , T_{\mathrm{X}})\delta _{\omega ,\omega '}. \end{aligned}$$

(C4)

The noise-noise correlations for the Z bath can also be written in the following compact forms provided the bath is an N or a Kitaev chain with a real $\Delta _{\mathrm{Z}}$:

$$\begin{aligned} \langle {\tilde{\eta }}_1^{\mathrm{Z \dagger }}(\omega ) {\tilde{\eta }}_1^{\mathrm{Z}}(\omega ')\rangle= & {} -\frac{ \gamma _\mathrm{YZ}^2}{\pi } \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{1,1}(\omega )]f(\omega , T_{\mathrm{Z}})\delta _{\omega ,\omega '} \, , \nonumber \\ \langle {\tilde{\eta }}_2^{\mathrm{Z \dagger }}(\omega ) {\tilde{\eta }}_2^{\mathrm{Z}}(\omega ')\rangle= & {} - \frac{ \gamma _{\mathrm{YZ}}^2}{\pi } \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{2,2}(\omega )]f(\omega , T_{\mathrm{Z}})\delta _{\omega ,\omega '} \, , \nonumber \\ \langle {\tilde{\eta }}_1^{\mathrm{Z \dagger }}(\omega ) {\tilde{\eta }}_2^{\mathrm{Z}}(\omega ')\rangle= & {} \frac{ \gamma _{\mathrm{YZ}}^2}{\pi } \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{2,1}(\omega )]f(\omega , T_{\mathrm{Z}})\delta _{\omega ,\omega '} \, , \nonumber \\ \langle {\tilde{\eta }}_2^{\mathrm{Z \dagger }}(\omega ) {\tilde{\eta }}_1^{\mathrm{Z}}(\omega ')\rangle= & {} - \frac{ \gamma _{\mathrm{YZ}}^2}{\pi } \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{1,2}(\omega )]f(\omega , T_{\mathrm{Z}})\delta _{\omega ,\omega '}. \end{aligned}$$

(C5)

The above noise-noise correlations (C3–C5) are in the form of fluctuation-dissipation relations. It is clear from these relations that we need to find the boundary Green’s function of the leads to evaluate the noise-noise correlations. We know an exact analytical form of these boundary Green’s functions for an N bath [54, 57]. However, we use numerical methods to find the boundary Green’s function of semiconducting and superconducting leads. In this paper, we apply the highly-convergent iterative method of Lopez Sancho et al. [76] to calculate the boundary Green’s function.

We note that the above noise-noise correlations, which are proportional to the off-diagonal terms of boundary Green’s function, are identically zero for an N lead. These off-diagonal terms only exist for the TS leads in the presence of superconducting pairing. Since we here often take a non-zero chemical potential for the N baths, the Fermi functions used in the aforementioned noise-noise correlations need some modifications. When the metallic $\alpha$ wire bath is kept at a chemical potential $\mu _\alpha$, the Fermi functions in the $\langle {\tilde{\eta }}_1^{ \alpha \dagger }(\omega ) {\tilde{\eta }}_1^{\alpha }(\omega ')\rangle$ and $\langle {\tilde{\eta }}_2^{\alpha \dagger }(\omega ) {\tilde{\eta }}_2^{\mathrm{\alpha }}(\omega ')\rangle$ should be modified as $f(\omega -\frac{\mu _\alpha }{\hbar }, T_{\alpha })$ and $f(\omega +\frac{\mu _\alpha }{\hbar }, T_{\alpha })$, respectively. Thus, the expressions of the junction currents would also be modified accordingly. We further clarify that we use an exact analytical expression for the boundary Green’s function for an N lead [54, 57].

Appendix D: Application of LEGF to a Device Made of Majorana Wire

Here, for a X-Y-Z device with L lattice sites, we use the following generalized basis: $\mathbf{b} \equiv [ b_1, b_2, b_3, b_4, \dots , b_{4L-3}, b_{4L-2}, b_{4L-1}, b_{4L}]^T$. For time $t > t_0$, the Heisenberg equations of motion of the Y wire variables are

$$\begin{aligned} {\dot{b}}_l= & {} -i \sum _{m=4L_{\mathrm{X}}+1}^{4L_{\mathrm{XY}}} {\mathcal {H}}^{\mathrm{Y}}_{lm} b_m +i\gamma _{\mathrm{XY}}\sum _{m=4L_{\mathrm{X}}+1}^{4L_{\mathrm{X}}+2} b_{m-4} \, \delta _{l,m}-i\gamma _{\mathrm{XY}}\sum _{m=4L_{\mathrm{X}}+3}^{4L_{\mathrm{X}}+4} \, b_{m-4} \delta _{l,m} \nonumber \\&+i\gamma _{\mathrm{YZ}}\sum _{m=4L_{\mathrm{XY}}-3}^{4L_{\mathrm{XY}}-2} b_{m+4} \, \delta _{l,m} -i\gamma _{\mathrm{YZ}}\sum _{m=4L_{\mathrm{XY}}-1}^{4L_{\mathrm{XY}}} \, b_{m+4} \delta _{l,m} \nonumber \\&- i\zeta _{\mathrm{XY}} \big ( b_{4L_{\mathrm{X}}-2}\, \delta _{l,4L_{\mathrm{X}}+1} -b_{4L_{\mathrm{X}}-3}\, \delta _{l,4L_{\mathrm{X}}+2}-b_{4L_{\mathrm{X}}}\, \delta _{l,4L_{\mathrm{X}}+3}+b_{4L_{\mathrm{X}}-1}\, \delta _{l,4L_{\mathrm{X}}+4}\big ) \nonumber \\&+i\zeta _{\mathrm{YZ}}\big (b_{4L_{\mathrm{XY}}+2}\, \delta _{l,4L_{\mathrm{XY}}-3} -b_{4L_{\mathrm{XY}}+1}\, \delta _{l,4L_{\mathrm{XY}}-2} \nonumber \\&- b_{4L_{\mathrm{XY}}+4}\, \delta _{l,4L_{XY}-1}+ b_{4L_{\mathrm{XY}}+3}\, \delta _{l,4L_{XY}}\big ), \end{aligned}$$

(D1)

where, $l = 4L_{\mathrm{X}}+1,\dots ,4L_{\mathrm{XY}}$. The equations of motion for the bath (X and Z wires) variables read as

$$\begin{aligned} {\dot{b}}_l= & {} -i \sum _{m=1}^{4L_{\mathrm{X}}} {\mathcal {H}}^{\mathrm{X}}_{lm} b_m +i\gamma _{\mathrm{XY}}\sum _{m=4L_{\mathrm{X}}-3}^{4L_{\mathrm{X}}-2} b_{m+4} \, \delta _{l,m}-i\gamma _{\mathrm{XY}}\sum _{m=4L_{\mathrm{X}}-1}^{4L_{\mathrm{X}}} \, b_{m+4} \delta _{l,m} \nonumber \\&+ i\zeta _{\mathrm{XY}}\left( b_{4L_{\mathrm{X}}+2}\, \delta _{l,4L_{\mathrm{X}}-3}-b_{4L_{\mathrm{X}}+1}\, \delta _{l,4L_{\mathrm{X}}-2}-b_{4L_{\mathrm{X}}+4}\, \delta _{l,4L_{\mathrm{X}}-1}+b_{4L_{\mathrm{X}}+3}\, \delta _{l,4L_{\mathrm{X}}}\right) , \qquad \nonumber \\ \end{aligned}$$

(D2)

for $l = 1, \dots , 4L_{\mathrm{X}}$, and

$$\begin{aligned} {\dot{b}}_l= & {} -i \sum _{m=4L_{\mathrm{XY}}+1}^{4L} {\mathcal {H}}^{\mathrm{Z}}_{lm} b_m +i\gamma _{\mathrm{YZ}}\sum _{m=4L_{\mathrm{XY}}+1}^{4L_{\mathrm{XY}}+2} b_{m-4} \, \delta _{l,m}-i\gamma _{\mathrm{YZ}}\sum _{m=4L_{\mathrm{XY}}+3}^{4L_{\mathrm{XY}}+4} \, b_{m-4} \delta _{l,m} \nonumber \\&-i\zeta _{\mathrm{YZ}}\left( b_{4L_{\mathrm{XY}}-2}\, \delta _{l,4L_{XY}+1}- b_{4L_{\mathrm{XY}}-3}\, \delta _{l,4L_{XY}+2} \right. \nonumber \\&\left. -b_{4L_{\mathrm{XY}}}\, \delta _{l,4L_{\mathrm{XY}}+3}+b_{4L_{\mathrm{XY}}-1}\, \delta _{l,4L_{\mathrm{XY}}+4}\right) , \qquad \end{aligned}$$

(D3)

for $l = 4L_{\mathrm{XY}}+1,\dots ,4L$. Like the Kitaev wire in the previous subsection 3.1, Eqs. D2 and D3 can be solved using the single-particle retarded Green’s function of the isolated bath wire. In the first step, the above equations for the variables of X and Z near the junctions are solved using the Green’s function to find formal solutions of $b_l(t)$ in the integral form in time domain. In the next step, these formal solutions of bath variables are Fourier transformed assuming steady state at long time, and the bath modes, ${\tilde{b}}_l(\omega )$ are written in terms of their noises and Y wire variables.

It is noteworthy that the matrix Green’s function for $\alpha$ Majorana/semiconductor wire with $L_\alpha$ sites is represented by a square matrix of dimension $4L_\alpha \times 4L_\alpha$ in $\mathbf{{b}}$ basis as $[{G}^+_\alpha (\tau )]_{l,m} \equiv [-i e^{-i{\mathcal {H}}^{ \alpha } \tau } ]_{lm} \theta (\tau )$, where $l,m=1,\dots ,4L_\alpha$ and ${\mathcal {H}}^{ \alpha }$ is the corresponding wire Hamiltonian. Substituting the formal solutions of Eq. D2 and Eq. D3 into Eq. D1, we get a set of generalized quantum Langevin equations for Y wire variables in which the contributions of the baths (boundary wires) enters into the equations via the self-energy corrections and the noise terms. Due to the presence of spin-orbit coupling in the SM junctions, both the self-energy and noise terms also depend on the spin-orbit coupling in addition to its dependency on the tunnel coupling. The steady-state solutions of the aforesaid Y wire equations are possible through the Fourier transformation method when the middle wire attains a unique NESS. In such cases, the solutions of the variables of Y wire, ${\tilde{b}}_l(\omega )$ read (for $l=4L_{\mathrm{X}}+1,\dots ,4L_{\mathrm{XY}}$)

$$\begin{aligned} {\tilde{b}}_l(\omega )= & {} \sum _{m=4L_{\mathrm{X}}+1}^{4L_{\mathrm{XY}}} {\tilde{G}}^+_{l,m}(\omega )\left( {\tilde{\xi }}_{1 \uparrow }^{\mathrm{X}} (\omega )\, \delta _{m,4L_{\mathrm{X}}+1}+{\tilde{\xi }}_{1 \downarrow }^{\mathrm{X}} (\omega )\, \delta _{m,4L_{\mathrm{X}}+2} \right. \nonumber \\&+{\tilde{\xi }}_{2 \uparrow }^{\mathrm{X}} (\omega ) \, \delta _{m,4L_{\mathrm{X}}+3} +{\tilde{\xi }}_{2 \downarrow }^{\mathrm{X}} (\omega ) \, \delta _{m,4L_{\mathrm{X}}+4} +{\tilde{\xi }}_{1 \uparrow }^{\mathrm{Z}} (\omega )\, \delta _{m,4L_{\mathrm{XY}} -3} \nonumber \\&\left. +{\tilde{\xi }}_{1 \downarrow }^{\mathrm{Z}} (\omega )\, \delta _{m,4L_{\mathrm{XY}} -2} +{\tilde{\xi }}_{2 \uparrow }^{\mathrm{Z}} (\omega ) \, \delta _{m,4L_{\mathrm{XY}}-1} +{\tilde{\xi }}_{2 \downarrow }^{\mathrm{Z}} (\omega ) \, \delta _{m,4L_{\mathrm{XY}}}\right) , \end{aligned}$$

(D4)

where ${\tilde{\xi }}^{\alpha }_{1 \sigma }(\omega )$ and ${\tilde{\xi }}^{\alpha }_{2 \sigma }(\omega )$ ($\sigma = \uparrow ,\downarrow$) are the Fourier transform of the noise terms associated with $\alpha$= X, Z bath (see E2 for definition), and ${\tilde{G}}^+_{l,m}(\omega )$ is the Fourier transform of the retarded Green’s function of the full system which, in this case, is a square matrix of dimension $4L_{\mathrm{Y}} \times 4L_{\mathrm{Y}}$. This retarded Green’s function is defined in a similar fashion to the relation (27), and the effective non-Hermitian Hamiltonian of Y wire reads $\tilde{{\mathcal {H}}}^{\mathrm{Y}}={\mathcal {H}}^{\mathrm{Y}}+ {\tilde{\Sigma }}^+_{\mathrm{X}}(\omega ) +{\tilde{\Sigma }}^+_{\mathrm{Z}}(\omega )$. Here, ${\mathcal {H}}^{\mathrm{Y}}$ is the Hamiltonian of the finite Majorana/SM wire. The structure of ${\tilde{\Sigma }}^+_{\mathrm{X/Z}}(\omega )$ is determined by the tunnelling Hamiltonian (37) and the Hamiltonian of the lead made of either a Majorana or an SM wire (13). The self-energy correction term associated with the X lead can be written as:

$$\begin{aligned}{}[{\tilde{\Sigma }}_{\mathrm{X}}^+(\omega )]_{lm}= & {} \sum _{k=4L_{\mathrm{X}}-3}^{4L_{\mathrm{X}}}\left[ \delta _{l,4L_{\mathrm{X}}+1} \left( \gamma _{\mathrm{XY}} \, {{\tilde{G}}^{\mathrm{X +}}_{4L_{\mathrm{X}}-3,k}(\omega )}-\zeta _{\mathrm{XY}} {{\tilde{G}}^{\mathrm{X +}}_{4L_{\mathrm{X}}-2,k}(\omega )}\right) \right. \nonumber \\&+\delta _{l,4L_{\mathrm{X}}+2} \left( \gamma _{\mathrm{XY}} \, {{\tilde{G}}^{\mathrm{X +}}_{4L_{\mathrm{X}}-2,k}(\omega )}+\zeta _{\mathrm{XY}} {{\tilde{G}}^{\mathrm{X +}}_{4L_{\mathrm{X}}-3,k}(\omega )}\right) \nonumber \\&-\delta _{l,4L_{\mathrm{X}}+3} \left( \gamma _{\mathrm{XY}} \, {{\tilde{G}}^{\mathrm{X +}}_{4L_{\mathrm{X}}-1,k}(\omega )}-\zeta _{\mathrm{XY}} {{\tilde{G}}^{\mathrm{X +}}_{4L_{\mathrm{X}},k}(\omega )}\right) \nonumber \\&\left. -\delta _{l,4L_{\mathrm{X}}+4} \left( \gamma _{\mathrm{XY}} \, {{\tilde{G}}^{\mathrm{X +}}_{4L_{\mathrm{X}},k}(\omega )}+\zeta _{\mathrm{XY}} {{\tilde{G}}^{\mathrm{X +}}_{4L_{\mathrm{X}}-1,k}(\omega )}\right) \right] \Gamma ^{\mathrm{X}}_{km}\, ,\nonumber \\ \mathrm{where}, \, \Gamma ^{\mathrm{X}}_{km}= & {} \left( {\mathcal {X}}_{k m}^{\mathrm{X}} \gamma _{\mathrm{XY}}+{\mathcal {Y}}_{k m}^{\mathrm{X}}\zeta _{\mathrm{XY}}\right) ,\nonumber \\ {\mathcal {X}}_{k m}^{\mathrm{X}}= & {} \sum _{s=4L_{\mathrm{X}}-3}^{4L_{\mathrm{X}}-2} \delta _{k,s} \,\delta _{m,s+4}-\sum _{s=4L_{\mathrm{X}}-1}^{4L_{\mathrm{X}}} \delta _{k,s} \,\delta _{m,s+4}\, , \nonumber \\ \text {and}, \,{\mathcal {Y}}_{k m}^{\mathrm{X}}= & {} \delta _{k,4L_{\mathrm{X}}-3} \,\delta _{m,4L_{\mathrm{X}}+2}-\delta _{k,4L_{\mathrm{X}}-2} \, \delta _{m,4L_{\mathrm{X}}+1} \nonumber \\&-\delta _{k,4L_{\mathrm{X}}-1}\, \delta _{m,4L_{\mathrm{X}}+4}+\delta _{k,4L_{\mathrm{X}}} \, \delta _{m,4L_{\mathrm{X}}+3}\, . \end{aligned}$$

(D5)

Similarly, the self-energy correction term for the Z bath reads

$$\begin{aligned}{}[{\tilde{\Sigma }}_{\mathrm{Z}}^+(\omega )]_{lm}= & {} \sum _{k=1}^{4} \left[ \delta _{l,4L_{\mathrm{XY}}-3} \left( \gamma _{\mathrm{YZ}} \, {{\tilde{G}}^{\mathrm{Z +}}_{1,k}(\omega )}+\zeta _{\mathrm{YZ}} {{\tilde{G}}^{\mathrm{Z +}}_{2,k}(\omega )}\right) \right. \nonumber \\&+\delta _{l,4L_{\mathrm{XY}}-2} \left( \gamma _{\mathrm{YZ}} \, {{\tilde{G}}^{\mathrm{Z +}}_{2,k}(\omega )}-\zeta _{\mathrm{YZ}} {{\tilde{G}}^{\mathrm{Z +}}_{1,k}(\omega )}\right) \nonumber \\&-\delta _{l,4L_{\mathrm{XY}}-1} \left( \gamma _{\mathrm{YZ}} \, {{\tilde{G}}^{\mathrm{Z +}}_{3,k}(\omega )}+\zeta _{\mathrm{YZ}} {{\tilde{G}}^{\mathrm{Z +}}_{4,k}(\omega )}\right) \nonumber \\&\left. -\delta _{l,4L_{\mathrm{XY}}} \left( \gamma _{\mathrm{YZ}} \, {{\tilde{G}}^{\mathrm{Z +}}_{4,k}(\omega )}-\zeta _{\mathrm{YZ}} {{\tilde{G}}^{\mathrm{Z +}}_{3,k}(\omega )}\right) \right] \Gamma ^{\mathrm{Z}}_{km}\, , \nonumber \\ \mathrm{where},\, \Gamma ^{\mathrm{Z}}_{km}= & {} \left( {\mathcal {X}}_{k m}^{\mathrm{Z }}\gamma _{\mathrm{YZ}}+ {\mathcal {Y}}_{k m}^{\mathrm{Z}}\zeta _{\mathrm{YZ}}\right) , \nonumber \\ ~{\mathcal {X}}_{k m}^{\mathrm{Z}}= & {} \sum _{s=1}^{2} \delta _{k,s} \,\delta _{m,4L_{\mathrm{XY}}-4+s}-\sum _{s=3}^{4} \delta _{k,s} \,\delta _{m,4L_{\mathrm{XY}}-4+s}\, , \nonumber \\ \mathrm{and}, \, {\mathcal {Y}}_{k m}^{\mathrm{Z}}= & {} -\delta _{k,1} \,\delta _{m,4L_{\mathrm{XY}}-2}+\delta _{k,2} \, \delta _{m,4L_{\mathrm{XY}}-3}+\delta _{k,3}\, \delta _{m,4L_{\mathrm{XY}}}-\delta _{k,4} \, \delta _{m,4L_{\mathrm{XY}}-1}.\nonumber \\ \end{aligned}$$

(D6)

It is worth-noting that for ${\tilde{G}}^{\mathrm{X +}}_{l,m}(\omega )$, $l,m =4L_{\mathrm{X}}-3,\dots ,4L_{\mathrm{X}}$ correspond to the right most site (i.e., $L_{\mathrm{X}}$-th site) of the X reservoir, whereas in case of ${\tilde{G}}^{\mathrm{Z +}}_{l,m}(\omega )$, $l,m=1,\dots ,4$ correspond to the left most site (i.e., $L_{\mathrm{XY}}+1$-th site) of the Z reservoir in the full X-Y-Z device. As we have observed earlier, the noise-noise correlations are essential for finding the currents at the junctions. We outline the derivations of such correlations for the Majorana wire bath in Appendix E.

Appendix E: Noise-Noise Correlations for Majorana Wire Leads

The noise-noise correlations for the Majorana wire leads have several new types of terms in comparison to the Kitaev chain in the presence of two types of spin degrees of freedom in the former. To this end, we first calculate the equilibrium correlation functions by applying Eq. 20. For an isolated $\alpha$ bath wire ($\alpha =X,Z$), whose r-th eigenvalue is $\hbar \omega _r^\alpha$, and $\{ \phi _{r \uparrow }^{\alpha },\phi _{r \downarrow }^{\alpha },\psi _{r \uparrow }^{\alpha },\psi _{r \downarrow }^{\alpha }\}$ are the components of the corresponding eigenfunction, the equilibrium correlation functions at temperature $T_\alpha$ are:

$$\begin{aligned} \langle c_{i' \sigma }^{\dagger }\, c_{j' \rho } \rangle= & {} \sum _{r>0}\left[ \phi ^{\alpha * }_{r \sigma }(i')\phi ^\alpha _{r \rho }(j')f(\omega _r^\alpha , T_\alpha )+\psi _{r \sigma }^\alpha (i') \psi _{r \rho }^{\alpha *}(j') f(-\omega _r^\alpha , T_\alpha )\right] , \nonumber \\ \langle c_{i' \sigma }\, c_{j' \rho }^\dagger \rangle= & {} \sum _{r>0}\left[ \psi _{r \sigma }^{\alpha *}(i') \psi _{r \rho }^\alpha (j') f(\omega _r^\alpha , T_\alpha )+\phi _{r \sigma }^\alpha (i') \phi _{r \rho }^{\alpha *}(j') f(-\omega _r^\alpha ,T_\alpha )\right] , \nonumber \\ \langle c_{i' \sigma } \, c_{j' \rho } \rangle= & {} \sum _{r>0}\left[ \psi _{r \sigma }^{\alpha *}(i') \phi _{r \rho }^\alpha (j')f(\omega _r^\alpha ,T_\alpha )+\phi _{r \sigma }^\alpha (i') \psi _{r \rho }^{\alpha *}(j')f(-\omega _r^\alpha ,T_\alpha )\right] , \nonumber \\ \langle c_{i' \sigma }^{ \dagger } \, c_{j'\rho }^{ \dagger } \rangle= & {} \sum _{r>0}\left[ \phi _{r \sigma }^{\alpha *}(i') \psi _{r \rho }^\alpha (j')f(\omega _r^\alpha , T_\alpha )+\psi _{r \sigma }^\alpha (i') \phi _{r \rho }^{\alpha *}(j')f(-\omega _r^\alpha , T_\alpha )\right] , \qquad \end{aligned}$$

(E1)

where $i',j'=1,\dots ,L_{\alpha }$ and $\sigma ,\rho =\uparrow , \downarrow$. The noise terms related to the X and Z baths are defined as follows:

$$\begin{aligned} \xi ^{\mathrm{X}}_{k \uparrow }(t)= \eta ^{\mathrm{X}}_{k \uparrow }(t)-{\bar{\zeta }}_{XY}\eta ^{\mathrm{X}}_{k \downarrow }(t) ,&\, \xi ^{\mathrm{X}}_{k \downarrow }(t)= \eta ^{\mathrm{X}}_{k \downarrow }(t)+{\bar{\zeta }}_{XY}\eta ^{\mathrm{X}}_{k \uparrow }(t) \, , \nonumber \\ \xi ^{\mathrm{Z}}_{k \uparrow }(t)= \eta ^{\mathrm{Z}}_{k \uparrow }(t)+{\bar{\zeta }}_{YZ}\eta ^{\mathrm{Z}}_{k \downarrow }(t) ,&\, \xi ^{\mathrm{Z}}_{k \downarrow }(t)= \eta ^{\mathrm{Z}}_{k \downarrow }(t)-{\bar{\zeta }}_{YZ}\eta ^{\mathrm{Z}}_{k \uparrow }(t) \, , \end{aligned}$$

(E2)

where $k=1, 2$ and ${\bar{\zeta }}_{\alpha \beta }=\zeta _{\alpha \beta }/\gamma _{\alpha \beta }$. These $\eta ^\alpha _{k \sigma }$ terms are quite similar to the noise terms of the Kitaev chain leads, and they read as

$$\begin{aligned} \eta _{k \uparrow }^{\mathrm{X}}(t)= & {} i\gamma _{\mathrm{XY}} (-1)^k \sum _{l=1}^{4L_{\mathrm{X}}}{[G_{\mathrm{X}}^+(t -t_0)]}_{4L_{\mathrm{X}}+2k-5,l} \, b_{l}(t_0)\, ,\nonumber \\ \eta _{k \downarrow }^{\mathrm{X}}(t)= & {} i \gamma _{\mathrm{XY}} (-1)^k \sum _{l=1}^{4L_{\mathrm{X}}}{[G_{\mathrm{X}}^+(t -t_0)]}_{4L_{\mathrm{X}}+2k-4,l} \, b_{l}(t_0) \, ,\nonumber \\ \eta _{k \uparrow }^{\mathrm{Z}}(t)= & {} i \gamma _{\mathrm{YZ}} (-1)^k \sum _{l=4L_{\mathrm{XY}}+1}^{4L}{[G_{\mathrm{Z}}^+(t -t_0)]}_{2k-1,l-4L_{\mathrm{XY}}} \, b_{l}(t_0)\, ,\nonumber \\ \eta _{k \downarrow }^{\mathrm{Z}}(t)= & {} i\gamma _{\mathrm{YZ}} (-1)^k \sum _{l=4L_{\mathrm{XY}}+1}^{4L}{[G_{\mathrm{Z}}^+(t -t_0)]}_{2k,l-4L_{\mathrm{XY}}} \, b_{l}(t_0)\, , \end{aligned}$$

(E3)

where ${[G_{\alpha }^+(t -t_0)]}_{l,m}$ are the components of bath Green’s functions, and $b_l$’s are operators in the generalized basis $\mathbf{b}$ (14). These bath Green’s functions can be expanded in terms of eigenfunctions of the Hamiltonian of the bath wire. For example, some components of ${[G_{\alpha }^+(t -t_0)]}_{l,m}$ are given below, and others can be derived identically.

$$\begin{aligned} {[G_\alpha ^+(\tau )]}_{4i'-3,4j'-3}= & {} - i \theta (\tau ) \sum _{r>0} \left( \phi _{r\uparrow }^\alpha (i' ) \phi _{r \uparrow }^{\alpha *}(j') e^{-i\, \omega ^\alpha _r \tau }+\psi ^{\alpha *}_{r \uparrow }(i') \psi _{r \uparrow }^\alpha (j') e^{i\, \omega _r^\alpha \tau }\right) , \\ {[G_\alpha ^+(\tau )]}_{4i'-2,4j'-2}= & {} - i \theta (\tau ) \sum _{r>0} \left( \phi _{r\downarrow }^\alpha (i' ) \phi _{r \downarrow }^{\alpha *}(j') e^{-i\, \omega ^\alpha _r \tau }+\psi ^{\alpha *}_{r \downarrow }(i') \psi _{r \downarrow }^\alpha (j') e^{i\, \omega _r^\alpha \tau }\right) , \\ {[G_\alpha ^+(\tau )]}_{4i'-1,4j'-1}= & {} - i \theta (\tau ) \sum _{r>0} \left( \psi _{r \uparrow }^\alpha (i') \psi _{r \uparrow }^{\alpha *}(j') e^{-i \, \omega _r^\alpha \tau }+\phi ^{\alpha *}_{r \uparrow }(i') \phi _{r \uparrow }^\alpha (j') e^{i\, \omega _r^\alpha \tau } \right) , \\ {[G_\alpha ^+(\tau )]}_{4i',4j'}= & {} - i \theta (\tau ) \sum _{r>0} \left( \psi _{r \downarrow }^\alpha (i') \psi _{r \downarrow }^{\alpha *}(j') e^{-i \, \omega _r^\alpha \tau }+\phi ^{\alpha *}_{r \downarrow }(i') \phi _{r \downarrow }^\alpha (j') e^{i\, \omega _r^\alpha \tau } \right) , \\ {[G_\alpha ^+(\tau )]}_{4i'-3,4j'-2}= & {} - i \theta (\tau ) \sum _{r>0} \left( \phi _{r\uparrow }^\alpha (i' ) \phi _{r \downarrow }^{\alpha *}(j') e^{-i\, \omega ^\alpha _r \tau }+\psi ^{\alpha *}_{r \uparrow }(i') \psi _{r \downarrow }^\alpha (j') e^{i \, \omega _r^\alpha \tau }\right) , \\ {[G_\alpha ^+(\tau )]}_{4i'-3,4j'-1}= & {} - i \theta (\tau ) \sum _{r>0} \left( \phi _{r \uparrow }^\alpha (i') \psi _{r \uparrow }^{\alpha *}(j') e^{-i \, \omega _{r \uparrow }^\alpha \tau }+\psi ^{\alpha *}_{r \uparrow }(i') \phi _{r \uparrow }^\alpha (j') e^{i \, \omega _r^\alpha \tau }\right) . \end{aligned}$$

In our steady-state transport analysis, we do not apply ${[G_{\alpha }^+(t -t_0)]}_{l,m}$ directly, we rather take the Fourier transformation of these bath Green’s functions using

$$\begin{aligned} {{\tilde{G}}^{\alpha +}_{l,m}(\omega )}=\frac{1}{2 \pi }\int _{-\infty }^{\infty } d\tau \, { [G_\alpha ^+(\tau )]}_{l,m}\, e^{i \omega \tau }. \end{aligned}$$

Since the steady-state current calculation requires the noise-noise correlations in the Fourier domain, we further define the Fourier transform of the aforementioned functions as ${\tilde{\eta }}_{k \sigma }^{\alpha }(\omega )=\frac{1}{2\pi }\int _{-\infty }^{\infty } dt \,\eta _{k \sigma }^{\alpha }(t)\, e^{i \omega t}$ and ${\tilde{\xi }}_{k \sigma }^{\alpha }(\omega )=\frac{1}{2\pi }\int _{-\infty }^{\infty } dt \,\xi _{k \sigma }^{\alpha }(t)\, e^{i \omega t}$.

Like the Kitaev chain leads in Appendix C, the noise-noise correlations for the Majorana wires are initially calculated in time domain using E2, E3, E1, and the orthonormality relations between the eigenvectors of the bath’s Hamiltonian. In the next step, we perform the Fourier transformation of the noise-noise correlations. Finally, we substitute ${{\tilde{G}}^{\alpha +}_{l,m}(\omega )}$ into the Fourier transformed noise-noise correlation to rewrite them in a simplified form. This calculation is tedious but straightforward. Hence, we here only mention the important results, which are the correlations between ${\tilde{\eta }}_{k \sigma }^{\alpha }(\omega )$’s. These last correlations, in turn, can be applied to derive the required noise-noise correlations $\langle {\tilde{\xi }}_{k \sigma }^{\alpha \dagger }(\omega ) {\tilde{\xi }}_{m \sigma }^{\alpha }(\omega ') \rangle$. The correlations, $\langle {\tilde{\eta }}_{k \sigma }^{\alpha \dagger }(\omega ) {\tilde{\eta }}_{m \sigma }^{\alpha }(\omega ') \rangle$, are expressed in a compact notation as follows:

$$\begin{aligned} \langle {\tilde{\eta }}_{k \sigma }^{\mathrm{X \dagger }}(\omega ) {\tilde{\eta }}_{m \sigma '}^{\mathrm{X}}(\omega ')\rangle= & {} (-1)^{k+m+1} \frac{ \gamma _{\mathrm{XY}}^2}{\pi } \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{4L_{\mathrm{X}}-6+2m+r(\sigma '),4L_{\mathrm{X}}-6+2k+r(\sigma )}(\omega )]\nonumber \\&\times f(\omega , T_{\mathrm{X}})\delta _{\omega ,\omega '} \, , \nonumber \\ \langle {\tilde{\eta }}_{k \sigma }^{\mathrm{Z \dagger }}(\omega ) {\tilde{\eta }}_{m \sigma '}^{\mathrm{Z}}(\omega ')\rangle= & {} (-1)^{k+m+1} \frac{ \gamma _{\mathrm{YZ}}^2}{\pi } \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{2m-2+r(\sigma '),2k-2+r(\sigma )}(\omega )]f(\omega , T_{\mathrm{Z}})\delta _{\omega ,\omega '},\nonumber \\ \end{aligned}$$

(E4)

where $\sigma , \sigma ' =\uparrow , \downarrow$ and $r(\uparrow )(r(\downarrow ))=1(2)$. It is easy to find all the noise-noise correlations using E2 and E4. The noise-noise correlations depend linearly on the boundary Green’s functions of the baths. We apply the same iterative method of Ref. [76] to numerically calculate these boundary Green’s functions for the Majorana wire and the SM lead.

For an SM lead, all correlations of type $\langle {\tilde{\eta }}_{k \sigma }^{ \alpha \dagger }(\omega ) {\tilde{\eta }}_{m \tau }^{\alpha }(\omega ')\rangle$ are zero when $k \ne m$. Like an N bath, an SM bath can also be kept at a non-zero chemical potential. When an SM lead is kept at chemical potential $\mu _\alpha$, the Eq. E4 needs modification, and the Fermi functions in the $\langle {\tilde{\eta }}_{1 \sigma } ^{ \alpha \dagger }(\omega ) {\tilde{\eta }}_{1 \tau }^{\alpha }(\omega ')\rangle$ and $\langle {\tilde{\eta }}_{2 \sigma } ^{\alpha \dagger }(\omega ) {\tilde{\eta }}_{2 \tau }^{\mathrm{\alpha }}(\omega ')\rangle$ should be replaced by $f(\omega -\frac{\mu _\alpha }{\hbar }, T_{\alpha })$ and $f(\omega +\frac{\mu _\alpha }{\hbar }, T_{\alpha })$, respectively.

Appendix F: Expressions of Electrical Currents $J^{e}_{\mathrm{XY}}$ and $J^{e}_{\mathrm{YZ}}$

The full expression of the electrical current from the X bath to the Y wire can be separated in three different parts. When the chemical potentials of the boundary leads are zero, only the first part survives in the presence of a temperature bias from the boundary leads. We also need to set the chemical potential $(\mu _{\mathrm{X}},\mu _{\mathrm{Z}})$ to zero for the TS bath(s) of various devices in the following current expressions.

$$\begin{aligned} J^{e}_{\mathrm{XY}}= & {} \int _{- \infty }^{ \infty }d \omega \,e\left( {\mathcal {T}}^1_{\mathrm{XY}}(\omega )+{\mathcal {T}}^2_{\mathrm{XY}}(\omega ) \right) (f(\omega -\frac{\mu _{\mathrm{X}}}{\hbar },T_{\mathrm{X}}) -f(\omega -\frac{\mu _{\mathrm{Z}}}{\hbar },T_{\mathrm{Z}})) \nonumber \\&+\int _{- \infty }^{ \infty }d \omega \,e\, {\mathcal {T}}^2_{\mathrm{XY}}(\omega ) \,(f(\omega +\frac{\mu _{\mathrm{X}}}{\hbar },T_{\mathrm{X}}) -f(\omega -\frac{\mu _{\mathrm{X}}}{\hbar },T_{\mathrm{X}})) \nonumber \\&+\int _{- \infty }^{ \infty }d \omega \,e\, {\mathcal {T}}^3_{\mathrm{XY}}(\omega )\, (f(\omega +\frac{\mu _{\mathrm{Z}}}{\hbar },T_{\mathrm{Z}}) -f(\omega -\frac{\mu _{\mathrm{Z}}}{\hbar },T_{\mathrm{Z}})), \end{aligned}$$

(F1)

where,

$$\begin{aligned} {\mathcal {T}}^1_{\mathrm{XY}}(\omega )= & {} \left( -\frac{2 \gamma _{\mathrm{XY}}^2}{\pi }\right) \sum _{l=2L_{\mathrm{X}}-1}^{2L_{\mathrm{X}}} (-1)^{l+1} \text {Im} \left[ [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,l}^* \, \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}}-1,l-2}(\omega )] \right. \\&+ [{\tilde{\Sigma }}_{\mathrm{X}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,l}^* \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}}-1,l-2}(\omega )] \\&+ \left. [{\tilde{\Sigma }}_{\mathrm{X}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1}^* [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+2,l} \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{l-2,2L_{\mathrm{X}}-1}(\omega )] \right] , \\ {\mathcal {T}}^2_{\mathrm{XY}}(\omega )= & {} \left( \frac{2 \gamma _{\mathrm{XY}}^2}{\pi }\right) \sum _{l=2L_{\mathrm{X}}-1}^{2L_{\mathrm{X}}} (-1)^{l+1} \text {Im} \left[ [{\tilde{\Sigma }}_{\mathrm{X}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2}\right. \\&\times [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,l}^* \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}},l-2}(\omega )] \\&+\left. [{\tilde{\Sigma }}_{\mathrm{X}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2}^* [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+2,l} \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{l-2,2L_{\mathrm{X}}}(\omega )] \right] , \end{aligned}$$

and,

$$\begin{aligned} {\mathcal {T}}^3_{\mathrm{XY}}(\omega )= & {} \left( \frac{2 \gamma _{\mathrm{YZ}}^2}{\pi }\right) \sum _{l=2L_{\mathrm{XY}}-1}^{2L_{\mathrm{XY}}} (-1)^{l+1} \text {Im} \left[ [{\tilde{\Sigma }}_{\mathrm{X}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1}[{\tilde{G}}^+(\omega )]^*_{2L_{\mathrm{X}}+1,2L_{\mathrm{XY}}} \right. \\&\times [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,l} \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{l-2L_{\mathrm{XY}}+2,2}(\omega )] \\&+ \left. [{\tilde{\Sigma }}_{\mathrm{X}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{XY}}}^* [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+2,l} \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{l-2L_{\mathrm{XY}}+2,2}(\omega )] \right] . \end{aligned}$$

We can similarly simplify the electrical current from the Y wire to the Z bath:

$$\begin{aligned} J^{e}_{\mathrm{YZ }}= & {} \int _{- \infty }^{ \infty }d \omega \,e\left( {\mathcal {T}}^1_{\mathrm{YZ}}(\omega )+{\mathcal {T}}^2_{\mathrm{YZ}}(\omega ) \right) (f(\omega -\frac{\mu _{\mathrm{X}}}{\hbar },T_{\mathrm{X}}) -f(\omega -\frac{\mu _{\mathrm{Z}}}{\hbar },T_{\mathrm{Z}})) \nonumber \\&+\int _{- \infty }^{ \infty }d \omega \,e\, {\mathcal {T}}^2_{\mathrm{YZ}}(\omega ) \,(f(\omega +\frac{\mu _{\mathrm{X}}}{\hbar },T_{\mathrm{X}}) -f(\omega -\frac{\mu _{\mathrm{X}}}{\hbar },T_{\mathrm{X}})) \nonumber \\&+\int _{- \infty }^{ \infty }d \omega \,e\, {\mathcal {T}}^3_{\mathrm{YZ}}(\omega )\, (f(\omega +\frac{\mu _{\mathrm{Z}}}{\hbar },T_{\mathrm{Z}}) -f(\omega -\frac{\mu _{\mathrm{Z}}}{\hbar },T_{\mathrm{Z}})), \end{aligned}$$

(F2)

where,

$$\begin{aligned} {\mathcal {T}}^1_{\mathrm{YZ}}(\omega )= & {} \left( \frac{2 \gamma _{\mathrm{XY}}^2}{\pi }\right) \sum _{l=2L_{\mathrm{XY}}-1}^{2L_{\mathrm{XY}}} (-1)^{l+1} \text {Im} \left[ [{\tilde{\Sigma }}_{\mathrm{Z}}^+(\omega )]_{1,1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}}-1,2L_{\mathrm{X}}+1}\right. \\&\times [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}}-1,l}^* \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}}-1,l-2}(\omega )] \\&\left. + [{\tilde{\Sigma }}_{\mathrm{Z}}^+(\omega )]_{1,2}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}}-1,2L_{\mathrm{X}}+1}^* [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}},l} \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{l-2,2L_{\mathrm{X}}-1}(\omega )] \right] , \\ {\mathcal {T}}^2_{\mathrm{YZ}}(\omega )= & {} \left( -\frac{2 \gamma _{\mathrm{XY}}^2}{\pi }\right) \sum _{l=2L_{\mathrm{XY}}-1}^{2L_{\mathrm{XY}}} (-1)^{l+1} \text {Im} \left[ [{\tilde{\Sigma }}_{\mathrm{Z}}^+(\omega )]_{1,1}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}}-1,2L_{\mathrm{X}}+2} \right. \\&\times [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}}-1,l}^* \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}},l-2}(\omega )] \\&+[{\tilde{\Sigma }}_{\mathrm{Z}}^+(\omega )]_{1,2}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}}-1,2L_{\mathrm{X}}+2}^* [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}},l} \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{l-2,2L_{\mathrm{X}}}(\omega )] , \end{aligned}$$

and,

$$\begin{aligned} {\mathcal {T}}^3_{\mathrm{YZ}}(\omega )= & {} \left( -\frac{2 \gamma _{\mathrm{YZ}}^2}{\pi }\right) \sum _{l=2L_{\mathrm{XY}}-1}^{2L_{\mathrm{XY}}} (-1)^{l+1} \text {Im} \left[ [{\tilde{\Sigma }}_{\mathrm{Z}}^+(\omega )]_{1,1}[{\tilde{G}}^+(\omega )]^*_{2L_{\mathrm{XY}}-1,2L_{\mathrm{XY}}} \right. \\&\times [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}}-1,l} \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{l-2L_{\mathrm{XY}}+2,2}(\omega )] \\&\left. + [{\tilde{\Sigma }}_{\mathrm{Z}}^+(\omega )]_{1,2}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}}-1,2L_{\mathrm{XY}}}^*[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{XY}},l} \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{l-2L_{\mathrm{XY}}+2,2}(\omega )] \right] . \end{aligned}$$

We note that for TS-TS-N/TS-N-N devices with bias $\mu _{\mathrm{X}}=0, \mu _{\mathrm{Z}}=\mu$, the above expression becomes independent of ${\mathcal {T}}^2_{\mathrm{YZ}}(\omega )$ whose coefficient returns zero value for such a biasing.

We can further simplify the aforesaid expression (F2) for an N-TS-N device by considering $J^{e}_{\mathrm{XY }} = J^{e}_{\mathrm{YZ }}$ for either $\mu _{\mathrm{X}}=\mu _{\mathrm{Z}}=0$ or $\mu _{\mathrm{X}}=-\mu _{\mathrm{Z}}$ and $T_{\mathrm{X}}=T_{\mathrm{Z}}$. We thus obtain the following expression:

$$\begin{aligned} J^{e}_{\mathrm{YZ }}= & {} \int _{- \infty }^{ \infty }d \omega \,e\left( {\mathcal {T}}^1_{\mathrm{XY}}(\omega )+{\mathcal {T}}^2_{\mathrm{XY}}(\omega ) \right) f_{\mathrm{XZ}} +\int _{- \infty }^{ \infty }d \omega \,e\, ( {\mathcal {T}}^2_{\mathrm{XY}}(\omega )-{\mathcal {T}}^3_{\mathrm{XY}}(\omega )+{\mathcal {T}}^3_{\mathrm{YZ}}(\omega ))f_{{{\bar{X}}}X} \nonumber \\&+\int _{- \infty }^{ \infty }d \omega \,e\, {\mathcal {T}}^3_{\mathrm{YZ}}(\omega )f_{{{\bar{Z}}}Z}\, , \end{aligned}$$

(F3)

where

$$\begin{aligned} f_{\mathrm{XZ}}= & {} f(\omega -\frac{\mu _{\mathrm{X}}}{\hbar },T_{\mathrm{X}}) -f(\omega -\frac{\mu _{\mathrm{Z}}}{\hbar },T_{\mathrm{Z}}) ,\, \,f_{{{\bar{X}}}X}=f(\omega +\frac{\mu _{\mathrm{X}}}{\hbar },T_{\mathrm{X}}) -f(\omega -\frac{\mu _{\mathrm{X}}}{\hbar },T_{\mathrm{X}}) , \\&f_{\bar{\mathrm{Z}}\mathrm{Z}}=f(\omega +\frac{\mu _{\mathrm{Z}}}{\hbar },T_{\mathrm{Z}}) -f(\omega -\frac{\mu _{\mathrm{Z}}}{\hbar },T_{\mathrm{Z}}). \end{aligned}$$

Appendix G: Expressions of Energy Currents $J^{u}_{\mathrm{XY}}$ and $J^{u}_{\mathrm{YZ}}$

The energy currents $J^{u}_{\mathrm{XY}}$ (33) and $J^{u}_{\mathrm{YZ}}$ (34) are written in the generalized basis $\mathbf{a}$ (3) as

$$\begin{aligned} J^{u}_{\mathrm{XY}}= & {} 2\hbar \gamma _{\mathrm{XY}} \left( \gamma _{\mathrm{Y}} \mathrm{Im}[\langle a^\dagger _{2L_{\mathrm{X}}+3}(t) a_{2L_{\mathrm{X}}-1} (t)\rangle ]+ \epsilon _{\mathrm{Y}}\mathrm{Im}[\langle a^\dagger _{2L_{\mathrm{X}}+1}(t) a_{2L_{\mathrm{X}}-1} (t)\rangle ] \right) \nonumber \\&+ 2\gamma _{\mathrm{XY}} \Delta _{\mathrm{Y}}\mathrm{Im}[\langle a^\dagger _{2L_{\mathrm{X}}+3}(t) a_{2L_{\mathrm{X}}} (t)\rangle ], \end{aligned}$$

(G1)

$$\begin{aligned} J^{u}_{\mathrm{YZ}}= & {} 2\hbar \gamma _{\mathrm{YZ}} \left( \gamma _{\mathrm{Y}} \mathrm{Im}[\langle a^\dagger _{2L_{\mathrm{XY}}+1}(t) a_{2L_{\mathrm{XY}}-3} (t)\rangle ]+\epsilon _{\mathrm{Y}}\mathrm{Im}[\langle a^\dagger _{2L_{\mathrm{XY}}+1}(t) a_{2L_{\mathrm{XY}}-1} (t)\rangle ]\right) \nonumber \\&- 2 \gamma _{\mathrm{YZ}} \Delta _{\mathrm{Y}} \mathrm{Im}[\langle a^\dagger _{2L_{\mathrm{XY}}+1}(t) a_{2L_{\mathrm{XY}}-2} (t)\rangle ]. \end{aligned}$$

(G2)

Since the total energy is conserved in the middle Y wire, the values of $J^{u}_{\mathrm{XY}}$ and $J^{u}_{\mathrm{YZ}}$ are equal. These energy currents can also be separated in three parts like the electrical currents in Appendix F. We again need to set the chemical potential $(\mu _{\mathrm{X}},\mu _{\mathrm{Z}})$ to zero for the TS bath(s) of different devices in the following expressions of $J^{u}_{\mathrm{XY}}$ and $J^{u}_{\mathrm{YZ}}$.

$$\begin{aligned} J^{u}_{\mathrm{XY}} =J^{u}_{\mathrm{YZ}}= & {} \int _{- \infty }^{ \infty } d\omega \, \hbar \left( {\mathcal {A}}^1(\omega )+{\mathcal {A}}^2(\omega ) \right) f_{\mathrm{XZ}} +\int _{- \infty }^{ \infty } d \omega \, \hbar \, {\mathcal {A}}^2(\omega ) f_{{{\bar{X}}}X} \nonumber \\&+\int _{- \infty }^{ \infty } d \omega \, \hbar \, {\mathcal {A}}^3(\omega )\, f_{{{\bar{Z}}}Z}, \end{aligned}$$

(G3)

where,

$$\begin{aligned} {\mathcal {A}}^1(\omega )= & {} \left( \frac{2 \gamma _{\mathrm{XY}}^2}{\pi }\right) \sum _{l=2L_{\mathrm{X}}+1}^{2L_{\mathrm{X}}+2} (-1)^{l+1}\text {Im} \left[ \gamma _{\mathrm{Y}} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,l}^*\text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}}-1,l-2}(\omega )] \right. \\&+\epsilon _{\mathrm{Y}} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,l}^*\text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}}-1,l-2}(\omega )] \\&+ \left( \gamma _{\mathrm{Y}}[{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,l}^* \right. \\&- \Delta _{\mathrm{Y}}[{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+2,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,l}^* \\&+ \left. \epsilon _{\mathrm{Y}} [{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,l}^* \right) \\&\times \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}}-1,l-2}(\omega )] \\&+ \left( \gamma _{\mathrm{Y}}[{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,2L_{\mathrm{X}}+1}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+2,l}^* \right. \\&- \Delta _{\mathrm{Y}}[{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+2,2L_{\mathrm{X}}+2} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,2L_{\mathrm{X}}+1}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+2,l}^* \\&+ \left. \epsilon _{\mathrm{Y}} [{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+2,l}^* \right) \\&\left. \times \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{l-2,2L_{\mathrm{X}}-1}(\omega )] \right] , \\ {\mathcal {A}}^2(\omega )= & {} \left( \frac{2 \gamma _{\mathrm{XY}}^2}{\pi }\right) \sum _{l=2L_{\mathrm{X}}+1}^{2L_{\mathrm{X}}+2} (-1)^{l+1} \text {Im} \left[ \Delta _{\mathrm{Y}}[{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,l}^*\text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}},l-2}(\omega )] \right. \\&-\left( \gamma _{\mathrm{Y}}[{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,l}^* \right. \\&-\Delta _{\mathrm{Y}}[{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+2,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,l}^* \\&\left. +\epsilon _{\mathrm{Y}}[{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,l}^*\right) \\&\times \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{2L_{\mathrm{X}},l-2}(\omega )] \\&-\left( \gamma _{\mathrm{Y}}[{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,l}^* \right. \\&-\Delta _{\mathrm{Y}}[{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+2,2L_{\mathrm{X}}+2} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+2,l}^* \\&\left. +\epsilon _{\mathrm{Y}}[{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+2,l}^*\right) \\&\times \text {Im}[{\tilde{G}}^{\mathrm{X +}}_{l-2,2L_{\mathrm{X}}}(\omega )] , \end{aligned}$$

and,

$$\begin{aligned} {\mathcal {A}}^3(\omega )&= \left( -\frac{2 \gamma _{\mathrm{YZ}}^2}{\pi }\right) \sum _{l=2L_{\mathrm{XY}}-1}^{2L_{\mathrm{XY}}} (-1)^{l+1} \text {Im} \left[ \left( \gamma _{\mathrm{Y}} [{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,2L_{\mathrm{XY}}}^* \right. \right. \\&\quad -\Delta _{\mathrm{Y}} [{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+2,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,2L_{\mathrm{XY}}}^* \\&\quad \left. +\epsilon _{\mathrm{Y}} [{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+1} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{XY}}}^* \right) \\&\quad \times [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,l} \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{l-2L_{\mathrm{XY}}+2,2}(\omega )] \\&\quad +\left( \gamma _{\mathrm{Y}} [{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,2L_{\mathrm{XY}}}^* \right. \\&\quad -\Delta _{\mathrm{Y}} [{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+2,2L_{\mathrm{X}}+2} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+3,2L_{\mathrm{XY}}}^* \\&\quad \left. +\epsilon _{\mathrm{Y}} [{\tilde{\Sigma }}_{\mathrm{X }}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{X}}+2} [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+1,2L_{\mathrm{XY}}}^* \right) \\&\quad \times [{\tilde{G}}^+(\omega )]_{2L_{\mathrm{X}}+2,l} \text {Im}[{\tilde{G}}^{\mathrm{Z +}}_{l-2L_{\mathrm{XY}}+2,2}(\omega )] . \end{aligned}$$

Appendix H: Numerical Time-Evolution of Density Matrix

To study the time-evolution dynamics of the devices with the Kitaev chains, we first consider three disjoint wires X, Y and Z, where $L_{\mathrm{X}}, ~L_{\mathrm{Z}}>> L_{\mathrm{Y}}$, so that we can treat X and Z as baths. At $t=t_0$, two opposite ends of Y wire are connected to the X and Z bath. The equilibrium density matrix for sites on the isolated boundary wires kept at zero chemical potential reads

$$\begin{aligned} \langle a^{\dagger }_{l}(t_0) a_{m}(t_0) \rangle = \sum _r {\psi _r^{\alpha }}^*(l)\psi _r^{\alpha }(m) f(\omega ^{\alpha }_r,T_{\alpha }), \end{aligned}$$

(H1)

with $l,m=1,\dots ,2L_{\mathrm{X}}$ for $\alpha =\mathrm{X}$ and $l,m=2L_{ \mathrm XY}+1,\dots ,2L$ for $\alpha =\mathrm{Z}$. Here, $f(\omega ,T)=1/(\exp [\hbar \omega /k_BT]+1)$ is the Fermi function. However, we need to modify the aforesaid expression in case of an N bath kept at a chemical potential $\mu _\alpha$. The Fermi function in the $\langle a^{\dagger }_{2l'-1}(t_0) a_{2m'-1}(t_0) \rangle$ and $\langle a^{\dagger }_{2l'}(t_0) a_{2m'}(t_0) \rangle$ should be modified as $f(\omega ^{\alpha }_r-\mu _\alpha / \hbar ,T_{\alpha })$ and $f(\omega ^{\alpha }_r+\mu _\alpha /\hbar ,T_{\alpha })$, respectively. All other correlations involving odd-even and even-odd indices are zero for an N bath.

Since X, Y, and Z wires are disconnected at $t \le t_0$, their uncorrelated operators at $t=t_0$ satisfy

$$\begin{aligned}&\langle a^{\dagger }_{l}(t_0) a_{p}(t_0) \rangle =\langle a_{l}(t_0) a^{\dagger }_{p}(t_0)\rangle =0 \, ,\nonumber \\&\langle a_{l}(t_0) a_{p}(t_0) \rangle =\langle a_{l}^{\dagger }(t_0) a^{\dagger }_{p}(t_0)\rangle =0 \,, \end{aligned}$$

(H2)

where l and p represent indices corresponding to two different wires of the hybrid device. We further choose some arbitrary initial density matrix for the Y wire, such as,

$$\begin{aligned} \langle a_{2l'}(t_0) a_{2m'-1}(t_0) \rangle = \left\{ \begin{array}{l l} n_{l'} &{} \quad \text{ when }~l'=m' \\ 0 &{} \quad \text{ when }~ l' \ne m' \end{array} \right. \end{aligned}$$

(H3)

for physical sites: $l',m' \in \{L_{\mathrm{X}}+1,\dots , L_{\mathrm{XY}}\}$ and $n_{l'}$ denotes the number of fermions at a site $l'$. We then connect the three wires through the tunneling Hamiltonians at time $t_0$ and investigate the time-evolution of the full device using the Heisenberg equations of motion. The solution of the equations of motion for $t \ge t_0$ is given by

$$\begin{aligned} \mathbf{a}(t)=i {\mathcal {G}}^+(t-t_0)\mathbf{a}(t_0), \end{aligned}$$

(H4)

where ${\mathcal {G}}^+(\tau )=-i e^{-i H^{\mathrm{F}} \tau / \hbar } \theta (\tau )=\sum _{l,m}{\mathcal {G}}_{lm}^+(\tau )a^{\dagger }_la_m$ is the retarded Green’s function of the full device and the matrix elements are ${\mathcal {G}}_{lm}^+ (\tau ) \equiv [-i e^{-i {\mathcal {H}}^{\mathrm{F}} \tau } ]_{lm} \theta (\tau )$, where ${\mathcal {H}}^{\mathrm{F}}_{lm}$ is the matrix Hamiltonian associated with the full Hamiltonian $H^\mathrm{F}$ (21). Suppose, $\Psi _q(m)$ and $\Lambda _q$ denote the eigenvectors and eigenvalues of the full Hamiltonian matrix, ${\mathcal {H}}^{\mathrm{F}}_{lm}$, hence they satisfy

$$\begin{aligned} \sum _{m=1}^{2L} {\mathcal {H}}^{\mathrm{F}}_{lm}\Psi _q(m)=\Lambda _q \Psi _q(l),~~l=1,2,\dots ,2L. \end{aligned}$$

(H5)

Further, the matrix elements of full Green’s function can be expanded in the following form for $t >t_0$:

$$\begin{aligned} {\mathcal {G}}^+_{rs}(t-t_0)=- i \sum _{q=1}^{2L} \Psi _{q}(r) \Psi _{q}^*(s) e^{- i \Lambda _q (t-t_0)} \, , \end{aligned}$$

(H6)

where $r,s \in \{1,\dots , 2L\}$. Finally, the time-evolved density matrix of the full device is expressed as

$$\begin{aligned} \langle a^\dagger _l(t) a_m(t) \rangle = \sum _{r,s=1}^{2L} {\mathcal {G}}^+_{ms}(t-t_0) \langle a^\dagger _{r}(t_0) a_{s}(t_0) \rangle [{\mathcal {G}}^+_{lr}(t-t_0)]^{\dagger }, \end{aligned}$$

(H7)

where we plug the initial density matrix $\langle a^\dagger _{r}(t_0) a_{s}(t_0) \rangle$ from Eqs. H1, H2, H3. Using the above expression for $\langle a^\dagger _l(t) a_m(t) \rangle$, we can numerically calculate time-evolution of both electrical and energy currents at the junctions. An extension of this direct time-evolution method for the devices with the Majorana wire leads is straightforward.

Fig. 12

Comparison of electrical currents at both junctions obtained from the first-principle/direct time-evolution numerics (full lines) and the generalized LEGF method (dashed lines) for a TS-N-TS at TP, where the TS wires are made of Kitaev chains. The initial numbers of spinless electrons at the middle N wire for the direct time-evolution numerics are provided on the heading of the top panel. In both panels, $L_{\mathrm{X}}=L_{\mathrm{Z}}=3000, L_{\mathrm{Y}}=3$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=1$, $\gamma _{\mathrm{Y}}=0.5$, $\Delta _{\mathrm{X}}=\Delta _{\mathrm{Z}}=0.3,\Delta _{\mathrm{Y}}=0$, $\epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Z}}=0$, $\epsilon _{\mathrm{Y}}=0.05$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.25$, $T_{\mathrm{X}}=0.02, T_{\mathrm{Z}}=0.2$ and $\mu _{\mathrm{X}}=\mu _{\mathrm{Z}}=0$. All above parameters except lengths are in units of $\gamma$

Fig. 13

Comparison of electrical currents at both junctions obtained from the first-principle/direct time-evolution numerics (full lines) and the generalized LEGF method (dashed lines) for a TS-TS-N device made of Kitaev chains. The initial numbers of spinless fermions at the middle TS wire for the direct time-evolution numerics are given on the heading of the top panel. In both panels, $L_{\mathrm{X}}=L_{\mathrm{Z}}=2000, L_{\mathrm{Y}}=20$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=1$, $\gamma _{\mathrm{Y}}=0.5$, $\Delta _{\mathrm{X}}=0.3,\Delta _{\mathrm{Y}}=0.1,\Delta _{\mathrm{Z}}=0$, $\epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Y}}=\epsilon _{\mathrm{Z}}=0$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.25$, $T_{\mathrm{X}}=T_{\mathrm{Z}}=0.02$ and $\mu _{\mathrm{X}}=0,\mu _{\mathrm{Z}}=0.5$. All above parameters except lengths are in units of $\gamma$

Fig. 14

Comparison of energy currents at both junctions obtained from the first-principle/direct time-evolution numerics and the generalized LEGF method for an N-TS-N device made of a Kitaev chain. The initial numbers of spinless electrons ($n_{l'}$) at the middle TS wire used for the time-evolution numerics are indicated on the headings. In both panels, $L_{\mathrm{X}}=L_{\mathrm{Z}}=2000, L_{\mathrm{Y}}=10$, $\gamma _{\mathrm{X}}=\gamma _{\mathrm{Z}}=\gamma =1$, $\gamma _{\mathrm{Y}}=0.5$, $\Delta _{\mathrm{X}}=\Delta _{\mathrm{Z}}=0,\Delta _{\mathrm{Y}}=0.15$, $\epsilon _{\mathrm{X}}=\epsilon _{\mathrm{Z}}=0$, $\epsilon _{\mathrm{Y}}=0.01$, $\gamma _{\mathrm{XY}}=\gamma _{\mathrm{YZ}}=0.25$, $T_{\mathrm{X}}=0.02,T_{\mathrm{Z}}=0.02$ and $\mu _{\mathrm{X}}=0.2,\mu _{\mathrm{Z}}=-0.4$. All above parameters except lengths are in units of $\gamma$

We next demonstrate the correctness of our steady-state electrical and energy current formulas in different devices by comparing them with the long-time values of the time-dependent currents obtained using the above first-principle numerics. We show results for steady-state electrical transport in TS-N-Z and TS-TS-Z with Z=N and TS at TP. First, we consider TS-N-TS at TP device in Fig. 12 to verify the validity of the generalized LEGF for such devices. In Fig. 12a, b, we respectively depict $J^e_{\mathrm{XY}}(t)$ and $J^e_{\mathrm{YZ}}(t)$ calculated from the direct time-evolution numerics and the generalized LEGF. These results show an excellent agreement between the two different methods at a long time. However, it takes very long time (equivalently, very long leads) in the direct numerics to achieve the steady state in such TS-N-TS at TP device in comparison to TS-TS-N, TS-N-N. It even requires much longer time to reach the steady state in a TS-TS-TS at TP device in the direct time-dependent numerics. The relaxation time scales for attaining steady-state (as observed in the direct time-evolution numerics) depend on the initial densities of the middle wire, the strength of the tunnel couplings, nature of the sub-gap states, and the bulk-gap of the superconductors [40]. In comparison, it takes much shorter time to numerically evaluate the steady-state values of $J^e_{\mathrm{XY}}(t)$ and $J^e_{\mathrm{YZ}}(t)$ in all devices using the LEGF. This indicates an advantage of using the LEGF method over direct/first-principle numerics for calculating steady-state currents in a hybrid device with relatively longer middle wire. Since the middle wire is made of an N in a TS-N-TS device, we have $J^e_{\mathrm{XY}}= J^e_{\mathrm{YZ}}$, and the steady-state value of the currents for this set of parameters is 0.0002325 (in units of $\gamma e$).

We further show the validity of the generalized LEGF for a TS-TS-N device in Fig. 13 by comparing the steady-state currents with the time-evolution numerics. In the upper (lower) panel, we plot the value of $J^e_{\mathrm{XY}}$ ($J^e_{\mathrm{YZ}}$) found from the generalized LEGF and the direct time-evolution numerics. Again, the LEGF results match accurately with those obtained from the direct time-evolution numerics at a long time. We notice here that $J^e_{\mathrm{XY}} \ne J^e_{\mathrm{YZ}}$ as the middle wire is a TS. We also get such excellent agreements between the LEGF and the direct time-evolution numerics for steady-state electrical currents in a TS-N-N and a TS-TS-TS at TP.

In Fig. 14, we compare the steady-state energy current $J^{u}_{\mathrm{XY}}$ and $J^{u}_{\mathrm{YZ}}$ in (G3) with the long-time values of the time-dependent currents obtained using the first-principle numerics. The values of $J^{u}_{\mathrm{XY}}(t)$ and $J^{u}_{\mathrm{YZ}}(t)$ in Fig. 14 due to a voltage bias are the same at long times for any initialization of the middle wire.