Master equation approach for transport through Majorana zero modes

Let us rewrite the tunnel-coupling Hamiltonian Eq. (2) as

$$H_{\mbox{\tiny SB}}=\sum_{\alpha}\big{(}F^{\dagger}_{\alpha}Q_{\alpha}+{\rm H.c.}\big{)}\,.$$

(3)

Here we denoted $Q_{\rm L}\equiv\gamma_{\rm L}=(f+f^{\dagger})=Q^{\dagger}_{\rm L}$, $Q_{\rm R}\equiv i\gamma_{\rm R}=f-f^{\dagger}=-Q^{\dagger}_{\rm R}$, and $F^{\dagger}_{\alpha}\equiv\sum_{k}t_{\alpha k}c^{\dagger}_{\alpha k}$. For the convenience of latter presentation, let us also introduce the correlation functions $g^{+}_{\alpha}(t-\tau)=\langle F^{\dagger}_{\alpha}(\tau)F_{\alpha}(t)\rangle_{\mbox{\tiny B}}$ and $g^{-}_{\alpha}(t-\tau)=\langle F_{\alpha}(t)F^{\dagger}_{\alpha}(\tau)\rangle_{\mbox{\tiny B}}$, where $\langle\cdots\rangle_{B}$ means the thermal average over the reservoir states of the leads. Straightforwardly, one can more explicitly obtain $g^{\pm}_{\alpha}(t-\tau)=\sum_{k}e^{-i\varepsilon_{\alpha k}(t-\tau)}\,|t_{\alpha k}|^{2}n^{\pm}_{\alpha}(\varepsilon_{\alpha k})$, where $n^{\pm}_{\alpha}$ are the Fermi occupied and empty functions. The sum of these two functions reads as $g_{\alpha}(t-\tau)=\sum_{k}e^{-i\varepsilon_{\alpha k}(t-\tau)}\,|t_{\alpha k}|^{2}$. Through the Fourier transformation

	$\displaystyle g_{\alpha}(t-\tau)$	$\displaystyle=\int^{\infty}_{-\infty}\frac{d\omega}{2\pi}\,e^{-i\omega(t-\tau)}\,\Gamma_{\alpha}(\omega)\,,$		(4a)
	$\displaystyle g^{\pm}_{\alpha}(t-\tau)$	$\displaystyle=\int^{\infty}_{-\infty}\frac{d\omega}{2\pi}\,e^{-i\omega(t-\tau)}\,\Gamma^{\pm}_{\alpha}(\omega)\,,$		(4b)

we also introduce the corresponding rates $\Gamma_{\alpha}(\omega)=2\pi\sum_{k}|t_{\alpha k}|^{2}\delta(\omega-\varepsilon_{\alpha k})$ and $\Gamma^{\pm}_{\alpha}(\omega)=\Gamma_{\alpha}(\omega)n^{\pm}_{\alpha}(\omega)$.

The basic idea of master equation approach to quantum transport is regarding the transport leads as a generalized environment and constructing an equation-of-motion for the reduced state of the central device system, $\rho(t)\equiv{\rm tr}_{\mbox{\tiny B}}\rho_{\mbox{\tiny T}}(t)$, where the trace is over the states of the transport leads and $\rho_{\mbox{\tiny T}}(t)$ is the total density operator of the device-plus-leads. The starting equation is the Schrödinger equation for the total wavefunction, or, equivalently, the Liouvillian equation for $\rho_{\mbox{\tiny T}}(t)$, $\dot{\rho}_{\mbox{\tiny T}}(t)=-i[H_{\mbox{\tiny T}},\rho_{\mbox{\tiny T}}(t)]$. Performing the trace ${\rm tr}_{\mbox{\tiny B}}(\cdots)$ on two sides of the Liouvillian equation, one obtains

$\displaystyle\dot{\rho}(t)$

$\displaystyle\!\!=\!\!-i[H_{\mbox{\tiny S}},\rho(t)]\!-\!i\!\sum_{\alpha}[Q_{\alpha},\rho^{+}_{\alpha}\!(t)]\!-\!i\!\sum_{\alpha}[Q^{\dagger}_{\alpha},\rho^{-}_{\alpha}\!(t)],$

(5)

where the two auxiliary density operators (ADOs) are introduced as

$\displaystyle\rho^{+}_{\alpha}(t)\!\equiv\!\mbox{tr}_{\mbox{\tiny B}}\{\rho_{\mbox{\tiny T}}(t)F^{\dagger}_{\alpha}\}~{}\text{and}~{}\rho^{-}_{\alpha}(t)\!\equiv\!\mbox{tr}_{\mbox{\tiny B}}\{F_{\alpha}\rho_{\mbox{\tiny T}}(t)\}\,.$

(6)

They satisfy $\rho^{+}_{\alpha}(t)=[\rho^{-}_{\alpha}(t)]^{\dagger}$. By means of a path-integral formulation in the fermionic coherent state representation Jin10083013 , it is possible to express the ADOs $\rho^{\pm}_{\alpha}(t)$ in terms of the reduced density operator $\rho(t)$. For the MZMs system under study, the detailed derivation is presented in Appendix A and the final result reads as

	$\displaystyle\dot{\rho}(t)$	$\displaystyle=-i[H_{\mbox{\tiny S}}(t),\rho(t)]+\Gamma_{1}(t){\cal D}[f]\rho(t)+\Gamma_{2}(t){\cal D}[f^{\dagger}]\rho(t)$
		$\displaystyle\quad+\Upsilon(t)\big{(}f\rho f+f^{\dagger}\rho f^{\dagger}\big{)}\,,$		(7)

where the Lindblad superoperator is defined by ${\cal D}[A]\rho=A\rho A^{\dagger}-\frac{1}{2}\{A^{\dagger}A,\rho\}$ and the time-dependent rates are given by

	$\displaystyle\Gamma_{1}(t)$	$\displaystyle\equiv\sum_{\alpha}\big{[}\Gamma^{-}_{\alpha 11}(t)+\Gamma^{+}_{\alpha 22}(t)\big{]}$
		$\displaystyle\quad-\sum_{\alpha}(-)^{\alpha}\big{[}\Gamma^{-}_{\alpha 12}(t)+\Gamma^{+}_{\alpha 21}(t)\big{]}\,,$		(8a)
	$\displaystyle\Gamma_{2}(t)$	$\displaystyle\equiv\sum_{\alpha}\big{[}\Gamma^{+}_{\alpha 11}(t)+\Gamma^{-}_{\alpha 22}(t)\big{]}$
		$\displaystyle\quad-\sum_{\alpha}(-)^{\alpha}\big{[}\Gamma^{+}_{\alpha 12}(t)+\Gamma^{-}_{\alpha 21}(t)\big{]}\,,$		(8b)
	$\displaystyle\Upsilon(t)$	$\displaystyle=\Gamma_{\rm L}(t)-\Gamma_{\rm R}(t)\,.$		(8c)

Assuming a wide-band limit approximation for the transport leads, i.e., $\Gamma_{\alpha}(\omega)\rightarrow\Gamma_{\alpha}$ in Eq. (4), the various individual rates can be more simply expressed as

	$\displaystyle\Gamma^{+}_{\alpha 11}(t)$	$\displaystyle\!=\!2\,{\rm Re}\!\int_{t_{0}}^{t}\!d\tau\,g^{+}_{\alpha}(t-\tau)u^{\dagger}_{11}(t,\tau),$		(9a)
	$\displaystyle\Gamma^{-}_{\alpha 22}(t)$	$\displaystyle\!=\!2\,{\rm Re}\!\int_{t_{0}}^{t}\!d\tau\,g^{-}_{\alpha}(t-\tau)u_{11}(t,\tau),$		(9b)
	$\displaystyle\Gamma^{+}_{\alpha 12}(t)$	$\displaystyle\!=\!2{\rm Re}\int_{t_{0}}^{t}\!d\tau\,g^{+}_{\alpha}(t-\tau)u^{\dagger}_{12}(t,\tau),$		(9c)
	$\displaystyle\Gamma^{-}_{\alpha 21}(t)$	$\displaystyle\!=\!2{\rm Re}\int_{t_{0}}^{t}\!d\tau\,g^{-}_{\alpha}(t-\tau)u_{12}(t,\tau),$		(9d)
	$\displaystyle\Gamma_{\alpha}(t)$	$\displaystyle\!=\!2\,{\rm Re}\!\int_{t_{0}}^{t}\!d\tau\,g_{\alpha}(t-\tau)\,,$		(9e)
while the other rates can be obtained through the following relations
	$\displaystyle\Gamma^{-}_{\alpha 11}(t)$	$\displaystyle=\Gamma_{\alpha}(t)-\Gamma^{+}_{\alpha 11}(t)\,,$		(9f)
	$\displaystyle\Gamma^{+}_{\alpha 22}(t)$	$\displaystyle=\Gamma_{\alpha}(t)-\Gamma^{-}_{\alpha 22}(t)\,,$		(9g)
	$\displaystyle\Gamma^{-}_{\alpha ij}(t)$	$\displaystyle=-\Gamma^{+}_{\alpha ij}(t),(i\neq j)\,.$		(9h)

In the above results, the reduced propagators (RPs) for the central MZMs system $u_{ij}(t,\tau)$, in terms of a matrix form $\bm{u}(t,\tau)$, satisfy the following equation-of-motion

$\displaystyle\frac{d}{dt}{\bm{u}}(t,t_{0})+\!\left(\begin{array}[]{cc}i\varepsilon_{\mbox{\tiny M}}+\Gamma&\Gamma_{\rm d}\\ \Gamma_{\rm d}&-i\varepsilon_{\mbox{\tiny M}}+\Gamma\\ \end{array}\right){\bm{u}}(t,t_{0})\!=\!0\,.$

(12)

Here we also considered the wide-band limit and introduced that $\Gamma=\Gamma_{\rm L}+\Gamma_{\rm R}$ and $\Gamma_{\rm d}=\Gamma_{\rm L}-\Gamma_{\rm R}$. Actually, ${\bm{u}}(t,t_{0})$ is nothing but the superconductor Green’s function (GF) matrix, specified here for the MZMs coupled the transport leads. Therefore, for the normal GFs (the diagonal elements), we have $u_{11}(t_{0},t_{0})=u_{22}(t_{0},t_{0})=1$; and for the anomalous GFs (the off-diagonal elements), we have $u_{12}(t_{0},t_{0})=u_{21}(t_{0},t_{0})=0$. Moreover, we have ${u}^{\dagger}_{11}(t,t_{0})={u}_{22}(t,t_{0})$ and ${u}^{\dagger}_{12}(t,t_{0})={u}_{21}(t,t_{0})$.

Owing to the symmetry properties of the ${\bm{u}}(t,t_{0})$ matrix, from Eq. (9c) and (9d), one can prove that $\Gamma^{+}_{\alpha 12}(t)=-\Gamma^{-}_{\alpha 21}(t)$. Together with Eq. (9h), we then have $\Gamma^{\pm}_{\alpha 12}(t)=\Gamma^{\pm}_{\alpha 21}(t)$. Remarkably, for the master equation, the total rates $\Gamma_{1}(t)$ and $\Gamma_{2}(t)$ of Eq. (8) are determined only by the sum of the diagonal rates $\Gamma^{\pm}_{\alpha 11/22}(t)$. The off-diagonal rates $\Gamma^{\pm}_{\alpha 12/21}(t)$ have no effect on the state evolution. However, as we will see later, the off-diagonal rates have important contribution to the transport currents.

So far, we have presented the main results of the master equation, i.e., Eq. (II.1) and the various rates associated. This master equation is applicable for arbitrary coupling strengths and bias voltages. In general, the nature of this master equation is non-Markovian, as explicitly reflected by the time-dependent rates in Eq. (9).

Markovian Limit.— As the usual treatment, one may consider the Markovian approximation which is characterized by time independent decohrence/dissipative rates. The Markovian approximation is largely associated with taking a long time limit for the transient rates. As is well known, in most cases, this treatment is very successful. Now, in this work, let us consider the Markovian approximation by taking the long time limit ($t\rightarrow\infty$) for the transient rates given by Eqs. (9). We obtain

	$\displaystyle\Gamma^{\pm}_{\alpha 11}$	$\displaystyle=\Gamma_{\alpha}\!\int\!d\omega\,n^{(\pm)}_{\alpha}(\omega)A_{11}(\omega),$		(13a)
	$\displaystyle\Gamma^{\pm}_{\alpha 22}$	$\displaystyle=\Gamma_{\alpha}\!\int\!d\omega\,n^{(\pm)}_{\alpha}(\omega)A_{22}(\omega),$		(13b)
	$\displaystyle\Gamma^{+}_{\alpha 12}$	$\displaystyle=-\Gamma^{-}_{\alpha 12}=\Gamma_{\alpha}\!\int\!d\omega\,n^{(+)}_{\alpha}(\omega)A_{12}(\omega).$		(13c)

Here, precisely following the Green’s function theory, we defined the spectral function through

$${\bm{A}}(\omega)=\frac{1}{\pi}{\rm Re}[{\bm{u}}(\omega)]\,,$$

(14)

while ${\bm{u}}(\omega)$ is the Fourier transformation of the Green’s function matrix ${\bm{u}}(t)$ in time domain, i.e., ${\bm{u}}(\omega)=\int_{0}^{\infty}\!dte^{i\omega t}{\bm{u}}(t)$. Here the initial condition of ${\bm{u}}(t)$ has been considered, which makes the Fourier and Laplace transformations identical. From Eq. (12), performing a Laplace transformation, we obtain the solution in frequency domain as

$\displaystyle{\bm{u}}(\omega)=\frac{i}{Z}\left(\begin{array}[]{cc}\omega+\varepsilon_{\mbox{\tiny M}}+i\Gamma&-i\Gamma_{\rm d}\\ -i\Gamma_{\rm d}&\omega-\varepsilon_{\mbox{\tiny M}}+i\Gamma\\ \end{array}\right)\,,$

(17)

where $Z=(\omega+\varepsilon_{\mbox{\tiny M}}+i\Gamma)(\omega-\varepsilon_{\mbox{\tiny M}}+i\Gamma)+\Gamma_{\rm d}^{2}$. Then, the spectral functions are obtained as

	$\displaystyle A_{11}(\omega)$	$\displaystyle=\frac{\Gamma\left[(\omega+\varepsilon_{\mbox{\tiny M}})^{2}+4\Gamma_{\rm L}\Gamma_{\rm R}\right]}{\pi|Z|^{2}},$		(18a)
	$\displaystyle A_{22}(\omega)$	$\displaystyle=\frac{\Gamma\left[(\omega-\varepsilon_{\mbox{\tiny M}})^{2}+4\Gamma_{\rm L}\Gamma_{\rm R}\right]}{\pi|Z|^{2}},$		(18b)
	$\displaystyle A_{12}(\omega)$	$\displaystyle=A_{21}(\omega)=\frac{\Gamma_{\rm d}\left(\omega^{2}-\varepsilon_{\mbox{\tiny M}}^{2}-4\Gamma_{\rm L}\Gamma_{\rm R}\right)}{\pi|Z|^{2}}\,.$		(18c)

In general, the spectral functions hold the properties of $\!\int\!d\omega\,A_{11/22}(\omega)=1$ and $\!\int\!d\omega\,A_{12}(\omega)=0$.

In master equation approach, the transport currents should be calculated using the reduced density matrix $\rho(t)$. This approach holds, very naturally, the advantage of rendering the obtained currents time dependent. One may notice that using other approaches, such as the non-equilibrium Green’s function approach or the scattering matrix theory, it is very inconvenient (or impossible) to calculate the time-dependent transport currents. Let us start to consider the average of the current operators $\hat{I}_{\alpha}$. Using the Heisenberg equation, we have $\hat{I}_{\alpha}=-e\frac{d}{dt}\hat{N}_{\alpha}=ie[\hat{N}_{\alpha},H_{\mbox{\tiny T}}]$, where the total electron number operator of lead-$\alpha$ reads as $\hat{N}_{\alpha}=\sum_{k}c^{\dagger}_{\alpha k}c_{\alpha k}$. Then, we obtain

$\displaystyle\hat{I}_{\alpha}=-ie\big{(}Q^{\dagger}_{\alpha}F_{\alpha}-F^{\dagger}_{\alpha}Q_{\alpha}\big{)},$

(19)

The current in lead-$\alpha$ is $I_{\alpha}(t)\equiv\langle\hat{I}_{\alpha}\big{\rangle}_{\rm T}={\rm Tr}_{\rm T}[\hat{I}_{\alpha}\rho_{\mbox{\tiny T}}(t)]$, where the average is over the total states of the whole system-plus-leads, i.e., ${\rm Tr}_{\rm T}\equiv{\rm tr}_{\mbox{\tiny S}}{\rm tr}_{\mbox{\tiny B}}$. Further, we reexpress the current in terms of the ADOs

$\displaystyle I_{\alpha}(t)$

$\displaystyle=-ie\,{\rm tr}_{\mbox{\tiny S}}\big{[}Q^{\dagger}_{\alpha}\rho^{-}_{\alpha}(t)-\rho^{+}_{\alpha}(t)Q_{\alpha}\big{]}\,.$

(20)

As mentioned above and demonstrated in Appendix A, the ADOs $\rho^{\pm}_{\alpha}(t)$ can be expressed in terms of the reduced density matrix $\rho(t)$, see Eq. (74). Then, the current is obtained as

$\displaystyle I_{\alpha}(t)$

$\displaystyle\!=\!I_{\alpha 11}(t)\!+\!I_{\alpha 22}(t)\!+\!(\!-\!)^{\alpha}I_{\alpha 12}(t)\!+\!(\!-\!)^{\alpha}I_{\alpha 21}(t)\,,$

(21)

with the various component currents explicitly given by

	$\displaystyle I_{\alpha 11}(t)$	$\displaystyle=e\Big{[}\Gamma^{+}_{\alpha 11}(t)\rho_{00}(t)-\Gamma^{-}_{\alpha 11}(t)\rho_{11}(t)\Big{]},$		(22a)
	$\displaystyle I_{\alpha 22}(t)$	$\displaystyle=e\Big{[}\Gamma^{+}_{\alpha 22}(t)\rho_{11}(t)-\Gamma^{-}_{\alpha 22}(t)\rho_{00}(t)\Big{]},$		(22b)
	$\displaystyle I_{\alpha 12}(t)$	$\displaystyle=e\left[\Gamma^{+}_{\alpha 12}(t)\rho_{00}(t)-\Gamma^{-}_{\alpha 12}(t)\rho_{11}(t)\right],$		(22c)
	$\displaystyle I_{\alpha 21}(t)$	$\displaystyle=e\left[\Gamma^{+}_{\alpha 21}(t)\rho_{11}(t)-\Gamma^{-}_{\alpha 21}(t)\rho_{00}(t)\right]\,.$		(22d)

In deriving this result, the elements of the density matrix were defined through $\rho_{11}(t)={\rm tr}_{\mbox{\tiny S}}[f^{\dagger}f\rho(t)]$ and $\rho_{00}(t)={\rm tr}_{\mbox{\tiny S}}[ff^{\dagger}\rho(t)]$. In practice, they should be solved from Eq. (II.1). The various rates, $\Gamma^{\pm}_{\alpha ij}(t)$, are those given by Eq. (9). Noting that $\Gamma^{\pm}_{\alpha 12}(t)=\Gamma^{\pm}_{\alpha 21}(t)$, we thus have $I_{\alpha 12}(t)=I_{\alpha 21}(t)$.

In terms of rate process, which applies also to understanding the above master equation, the physical meaning of the various current components can be understood as follow. (i) In the current $I_{\alpha 11}(t)$, the first term is associated with normal tunneling process of an electron from the $\alpha$-lead to the MZMs quasiparticle, while the second term describes the inverse process of the first one. (ii) In the current $I_{\alpha 22}(t)$, the first term describes the Andreev reflection (AR) process of an electron entering the MZMs from the $\alpha$-lead and annihilating the existing quasiparticle, associated with formation of a Cooper pair; and the second term describes the inverse process of the first one, associated with splitting a Cooper pair. (iii) In $I_{\alpha 12}(t)$, the first term describes another contribution of normal tunneling process of an electron from the $\alpha$-lead to the MZMs quasiparticle, while the second term is the inverse process of the first one. (iv) In $I_{\alpha 21}(t)$, again, the first and second terms describe additional contribution of the AR process and its inverse process, associated with formation and splitting of a Cooper pair, respectively.

Importantly, all the rates in $I_{\alpha 11}(t)$ and $I_{\alpha 22}(t)$, i.e., $\Gamma^{\pm}_{\alpha 11}$ and $\Gamma^{\pm}_{\alpha 22}$, involves internal self-energy propagation associated with the normal Green’s functions $u_{11}(t,\tau)$ and $u_{22}(t,\tau)$. By contrast, all the rates in $I_{\alpha 12}(t)$ and $I_{\alpha 21}(t)$, i.e., $\Gamma^{\pm}_{\alpha 12}$ and $\Gamma^{\pm}_{\alpha 21}$, involves internal self-energy propagation associated with the anomalous Green’s functions $u_{12}(t,\tau)$ and $u_{21}(t,\tau)$.

The above formal results of master equation and currents can be applied to transport under arbitrary bias voltages. For the two-lead (three-terminal) transport setup, following Ref. Nil08120403, , we simply consider the equally biased voltage of the leads with respect to the chemical potential of the grounded superconductor (which is set as zero), i.e., $\mu_{\rm L}/e=\mu_{\rm R}/e=V$. For this bias setup, the currents flow back to the leads from the superconductor through the grounded terminal. Between the leads, there is no transmission current from one lead to another.

In Fig. 2(a), we illustrate the transient behavior of the rates. For simplicity, in the numerical studies of this work, we assume the wide-band approximation for the leads. Rather than the usual constant rates, the transient behavior of the rates –even under the wide-band limit–reflects the non-Markovian nature. Using the non-Markovian transient rates, we solve the master equation and calculate the transient current, with the results as shown in Fig. 2(b) and (c) by the solid curves. As an interesting comparison, we also use the asymptotic rates $\Gamma^{\pm}_{\alpha ij}(t\to\infty)$ to solve the master equation and calculate the transient current, obtaining the results as shown in Fig. 2(b) and (c) by the dashed curves. From Fig. 2(b), we find that the occupation probability does not differentiate from each other too much. However, in Fig. 2(c), we find that the total current $I_{\rm L}(t)$ from the Markovian rates does not reveal any transient behavior. A simple explanation for this result is as follows. The total current $I_{\rm L}(t)$ contains two parts: one is from the normal tunneling process, which associates with electron entering the superconductor conditioned on the MZMs unoccupied; another part is from the Andreev reflection process, which is conditioned on the MZMs occupied. Then, the total current becomes free from the occupation of the MZMs.

In Fig. 2(e)-(f), we show the transient behavior of the component currents. As expected, even using the asymptotic rates, the switching-on transient behavior of $I_{{\rm L}11}(t)$ and $I_{{\rm L}22}(t)$ is obvious, in contrast to the total current $I_{\rm L}(t)$. We notice that, at the very beginning stage, $I_{{\rm L}22}(t)$ is negative. This is because, under the finite bias voltage $V$ (not large enough), the empty-occupation of the MZMs allows for happening of the inverse AR process, which splits a Cooper pair and results in occupation of the MZMs and another electron entering the lead from the superconductor, causing thus the negative $I_{{\rm L}22}(t)$ as observed. In Fig. 2(f), we find that the small off-diagonal current components are negative. This can be understood using the negative and very small off-diagonal rates shown in Fig. 2(a). Actually, as to be further analyzed in the following, at large $V$ limit, the off-diagonal rates will vanish. Thus, the off-diagonal currents will be zero under such limit. Finally, we remark that the results displayed by the solid curves in Fig. 2(e)-(f) are from using the non-Markovian transient rates as shown in Fig. 2(a). The more complicated transient behaviors are owing to the gradual formation of the asymptotic rates, during the rate process dynamics.

Markovian approximation can considerably simplify the master equation approach, making it extremely convenient for studying, not only the transport currents, but also current correlations and full-counting statistics. In this subsection, we analyze the conditions for Markovian approximation to the MME.

As briefly shown in Sec. II A, the Markov approximation is largely associated with taking a long-time limit for the transient rates, in present work, which are given by Eqs. (9). Hereafter, let us term the long-time asymptotic rates as Markovian rates. In Fig. 3, we display the behaviors of the Markovian rates, via their dependence on the bias voltage and the coupling asymmetry parameter $\eta=\Gamma_{\rm L}/\Gamma_{R}$. Owing to the a few symmetry properties of the rates $\Gamma^{\pm}_{\alpha ij}$, here we only plot the results of $\Gamma^{+}_{{\rm L}11/22}$ and $\Gamma^{+}_{{\rm L}12}$.

In Fig. 3(a), for the diagonal rates $\Gamma^{+}_{{\rm L}11/22}$, we find that they increase with the bias voltage $V$ (from negative to positive), while at the large $\pm V$ limits they approach, respectively, saturated values and zero. This bias-voltage dependence can be understood based on the lead-electron-occupation, which determines the integrated range of contribution to the rates. In Fig. 3(b), as an example of the off-diagonal rates, we find an antisymmetric dependence on $\pm V$. This behavior is determined by the following properties of the off-diagonal element of the spectral function matrix from Eq. (18c), i.e., $A_{12}(\omega)=A_{12}(-\omega)$ and $\int^{\infty}_{-\infty}d\omega A_{12}(\omega)=0$. Using them, we can easily prove that $\Gamma^{+}_{{\rm L}12}(V)+\Gamma^{+}_{{\rm L}12}(-V)=0$. In Fig. 3(c) and (d), we show the total rates $\Gamma_{1}$ and $\Gamma_{2}$ which characterize, respectively, the annihilation and creation rates of the MZMs-quasiparticle. For $\varepsilon_{\rm M}=0$, we find that the both total rates are bias ($V$) independent; while for $\varepsilon_{\rm M}\neq 0$, they reveal a symmetric dependence on $\pm V$ and take maximum and minimum at $V=0$ for $\Gamma_{1}$ and $\Gamma_{2}$, respectively. The reason for this behavior is as follows. First, based on Eq. (18c), we know $A_{12}(\omega)=A_{21}(\omega)$. Then, we can prove that $\Gamma^{-}_{{\rm L}12}(V)+\Gamma^{+}_{{\rm L}12}(V)=\Gamma_{\rm L}\int^{\infty}_{-\infty}d\omega A_{12}(\omega)=0$. Second, for $\varepsilon_{M}\neq 0$, we know that the peaks of $A_{11/22}(\omega)$ are centered at $\pm\varepsilon_{M}$. Then, we know that $\Gamma^{+}_{{\rm L}11}(V)$ increases with $V$ later than $\Gamma^{+}_{{\rm L}22}(V)$, when crossing the region around zero bias (from negative to positive). Because of the same reason, $\Gamma^{-}_{{\rm L}11}(V)$ decreases with $V$ earlier than $\Gamma^{-}_{{\rm L}22}(V)$, when crossing the region around zero bias. As a result, the total rates, $\Gamma_{1}(V)=\Gamma^{-}_{{\rm L}11}(V)+\Gamma^{+}_{{\rm L}22}(V)$ and $\Gamma_{2}(V)=\Gamma^{+}_{{\rm L}11}(V)+\Gamma^{-}_{{\rm L}22}(V)$, exhibit the behavior as shown in Fig. 3(c) and (d).

In Fig. 4 we show further the transient behaviors of the rates under different bias voltages. We see that for small bias voltage, the transient behavior (non-Markovian feature) is remarkable, while with increase of the bias voltage the rates more rapidly approach the (Markovian) asymptotic results. In practice, when the bias voltage is considerably larger than the broadening width, i.e., in the sequential tunneling regime, the Markovian master equation approach should work very well. In Fig. 4(a) and (b), we show the transient behaviors of diagonal and off-diagonal rates for single-lead coupling, while in (c) and (d) for a two-lead coupling configuration. We find qualitatively similar behaviors for both coupling configurations. This is because the rates $\Gamma^{\pm}_{\alpha ij}$ are dominantly affected by the Fermi occupation of the $\alpha$-lead electrons, while influence of the opposite lead is weakly mediated through the propagating function ${\bm{u}}(\omega)$ via the self-energy effect.

Finally, we carry out the analytic results under wide-band and large-bias limits. At large $V$ limit, $\mu_{\alpha}\to\infty$, we have $n^{+}_{\alpha}=1$ and $n^{-}_{\alpha}=0$. Then, straightforwardly, one can prove that $g^{+}_{\alpha}(t-\tau)\to\Gamma_{\alpha}\delta(t-\tau)$, which yields $\Gamma^{+}_{\alpha 11}=\Gamma^{+}_{\alpha 22}=\Gamma_{\alpha}$. Owing to $n^{-}_{\alpha}=0$, we simply have $\Gamma^{-}_{\alpha ij}=0$. Substituting these results into Eq. (II.1), we obtain a very simple Lindblad-type master equation which, however, contains the central physics of Andreev process associated with the MZMs, and makes the results different from transport through the conventional single-level quantum dot.

Let us introduce the more conventional (retarded) superconductor Green’s functions

	$\displaystyle G^{r}_{11}(t,\tau)$	$\displaystyle\!=\!-i\Theta(t-\tau)\langle\{f(t),f^{\dagger}(\tau)\}\rangle\,,$		(23a)
	$\displaystyle G^{r}_{22}(t,\tau)$	$\displaystyle\!=\!-i\Theta(t-\tau)\langle\{f^{\dagger}(t),f(\tau)\}\rangle\,,$		(23b)
	$\displaystyle G^{r}_{12}(t,\tau)$	$\displaystyle\!=\!-i\Theta(t-\tau)\langle\{f(t),f(\tau)\}\rangle\,,$		(23c)
	$\displaystyle G^{r}_{21}(t,\tau)$	$\displaystyle\!=\!-i\Theta(t-\tau)\langle\{f^{\dagger}(t),f^{\dagger}(\tau)\}\rangle\,.$		(23d)

From the equation-of-motion of ${\bm{u}}(t,\tau)$, i.e., Eq(..), we know that ${\bm{u}}(t,\tau)=i{\bm{G}}^{r}(t,\tau)$. Then, the spectral function matrix can expressed as

$\displaystyle{\bm{A}}(\omega)=\frac{1}{\pi}{\rm Re}[{\bm{u}}(\omega)]=-\frac{1}{\pi}{\rm Im}[{\bm{G}}^{r}(\omega)]$

(24)

From the Dyson equation, we have

	$\displaystyle{\rm Im}[{\bm{G}}^{r}(\omega)]=\frac{1}{2i}[{\bm{G}}^{r}(\omega)-{\bm{G}}^{a}(\omega)]$
	$\displaystyle~{}~{}~{}~{}=-i[{\bm{G}}^{r}(\omega)\widetilde{\bm{\Sigma}}^{r}(\omega){\bm{G}}^{a}(\omega)]$		(25)

The retarded self-energy matrix is given by

$\displaystyle\widetilde{\bm{\Sigma}}^{r}(\omega)=-\frac{i}{2}\sum_{\alpha={\rm L,R}}\left[{\bm{\Gamma}_{\alpha}}(\omega)+{\bm{\Gamma}_{\alpha}}(-\omega)\right]\,,$

(26)

while the rate matrices read as

	$\displaystyle{\bm{\Gamma}}_{\rm L}(\omega)$	$\displaystyle={\Gamma}_{\rm L}(\omega)\left(\begin{array}[]{cc}1&1\\ 1&1\\ \end{array}\right),$		(27c)
	$\displaystyle{\bm{\Gamma}}_{\rm R}(\omega)$	$\displaystyle={\Gamma}_{\rm R}(\omega)\left(\begin{array}[]{cc}1&-1\\ -1&1\\ \end{array}\right).$		(27f)

As to be clear in the following, we may denote ${\bm{\Gamma}}_{\alpha}(\omega)\equiv{\bm{\Gamma}}^{e}_{\alpha}$ and ${\bm{\Gamma}}_{\alpha}(-\omega)\equiv{\bm{\Gamma}}^{h}_{\alpha}$.

Based on Eq. (22), i.e., the time-dependent current formula of the master equation approach, the non-equilibrium Green’s function results of steady-state currents can be easily obtained. As an example, let us consider the steady-state current $I_{\rm L11}=\Gamma^{+}_{\rm L11}-\Gamma_{\rm L}\bar{\rho}_{11}$. We have

	$\displaystyle\Gamma^{+}_{{\rm L}11}$	$\displaystyle=$	$\displaystyle-2\Gamma_{\rm L}\int\frac{d\omega}{2\pi}n^{+}_{\rm L}(\omega){\rm Im}[{\bm{G}}^{r}(\omega)]_{11}\,,$
	$\displaystyle\bar{\rho}_{11}$	$\displaystyle=$	$\displaystyle\int\frac{d\omega}{2\pi}{\rm Im}[G_{11}^{<}(\omega)]\,.$		(28)

Here we have applied the relationship between $\rho_{11}(t)={\rm tr}_{\rm s}[f^{\dagger}f\rho(t)]$ and $G_{11}^{<}(t,\tau)=i\langle f^{\dagger}(\tau)f(t)\rangle$. After similar treatment to all the other component currents and applying Eq. (III.3) and the well known Keldysh equation for ${\bm{G}}^{<}(\omega)$, combination (together also with a procedure of symmetrization) of the component currents yields

	$\displaystyle I_{\rm L}$	$\displaystyle=\frac{e}{2h}\int\!d\omega\big{\{}{\cal T}^{ee}_{\rm LR}(\omega)\big{[}n^{+}_{\rm L}(\omega)-n^{+}_{\rm R}(\omega)\big{]}$
		$\displaystyle\quad\quad\quad\quad+{\cal T}^{hh}_{\rm LR}(\omega)\big{[}n^{-}_{\rm R}(-\omega)-n^{-}_{\rm L}(-\omega)\big{]}\big{\}}$
		$\displaystyle+\frac{e}{2h}\int\!d\omega\big{[}{\cal T}^{eh}_{\rm LLA}(\omega)+{\cal T}^{he}_{\rm LLA}(\omega)\big{]}\big{[}n^{+}_{\rm L}(\omega)-n^{-}_{\rm L}(-\omega)\big{]}$
		$\displaystyle+\frac{e}{2h}\int\!d\omega\big{\{}{\cal T}^{eh}_{\rm LRA}(\omega)\big{[}n^{+}_{\rm L}(\omega)-n^{-}_{\rm R}(-\omega)\big{]}$
		$\displaystyle\quad\quad\quad\quad+{\cal T}^{he}_{\rm LRA}(\omega)\big{[}n^{+}_{\rm R}(\omega)-n^{-}_{\rm L}(-\omega)\big{]}\big{\}}\,.$		(29)

In this result, the various transmission coefficients are given by

	$\displaystyle{\cal T}^{eh}_{\rm LLA}(\omega)$	$\displaystyle=$	$\displaystyle{\rm Tr}\big{[}{\bm{\Gamma}}^{e}_{\rm L}{\bm{G}}^{r}(\omega){\bm{\Gamma}}^{h}_{\rm L}{\bm{G}}^{a}(\omega)\big{]}\,,$
	$\displaystyle{\cal T}^{he}_{\rm LLA}(\omega)$	$\displaystyle=$	$\displaystyle{\rm Tr}\big{[}{\bm{\Gamma}}^{h}_{\rm L}{\bm{G}}^{r}(\omega){\bm{\Gamma}}^{e}_{\rm L}{\bm{G}}^{a}(\omega)\big{]}\,,$
	$\displaystyle{\cal T}^{ee}_{\rm LR}(\omega)$	$\displaystyle=$	$\displaystyle{\rm Tr}\big{[}{\bm{\Gamma}}^{e}_{\rm L}{\bm{G}}^{r}(\omega){\bm{\Gamma}}^{e}_{\rm R}{\bm{G}}^{a}(\omega)\big{]}\,,$
	$\displaystyle{\cal T}^{hh}_{\rm LR}(\omega)$	$\displaystyle=$	$\displaystyle{\rm Tr}\big{[}{\bm{\Gamma}}^{h}_{\rm L}{\bm{G}}^{r}(\omega){\bm{\Gamma}}^{h}_{\rm R}{\bm{G}}^{a}(\omega)\big{]}\,,$
	$\displaystyle{\cal T}^{eh}_{\rm LRA}(\omega)$	$\displaystyle=$	$\displaystyle{\rm Tr}\big{[}{\bm{\Gamma}}^{e}_{\rm L}{\bm{G}}^{r}(\omega){\bm{\Gamma}}^{h}_{\rm R}{\bm{G}}^{a}(\omega)\big{]}\,,$
	$\displaystyle{\cal T}^{he}_{\rm LRA}(\omega)$	$\displaystyle=$	$\displaystyle{\rm Tr}\big{[}{\bm{\Gamma}}^{h}_{\rm L}{\bm{G}}^{r}(\omega){\bm{\Gamma}}^{e}_{\rm R}{\bm{G}}^{a}(\omega)\big{]}\,.$		(30)

The meaning of these transmission coefficients is clear: ${\cal T}^{eh(he)}_{\rm LLA}$ describes the local Andreev-reflection process associated with incidence of an electron (a hole) and reflection of a hole (an electron) in the same (left) lead; ${\cal T}^{ee(hh)}_{\rm LR}$ describes the normal transport process of an electron (or a hole) from the left to the right lead; and ${\cal T}^{eh(he)}_{\rm LRA}$ describes the crossed Andreev-reflection process with incidence of an electron (a hole) from the left lead and reflection of a hole (an electron) in the right lead.

We may remark that, when reaching the compact form of the above results, we have properly reorganized the individual contributions in the various component currents. For instance, based on $I_{\rm L}=\sum_{i,j=1,2}I_{{\rm L}ij}$, we have performed the following combination of elements from the four component currents

	$\displaystyle\sum_{i,j=1,2}\Gamma^{e}_{\rm L}({\bm{G}}^{r}{\bm{\Gamma}}^{e}_{\rm R}{\bm{G}}^{a})_{ij}$
	$\displaystyle=\Gamma^{e}_{\rm L}\left(M^{e}_{{\rm R}11}+M^{e}_{{\rm R}21}+M^{e}_{{\rm R}12}+M^{e}_{{\rm R}22}\right)$
	$\displaystyle={\rm Tr}({\bm{\Gamma}}^{e}_{\rm L}{\bm{M}}^{e}_{\rm R})$		(31)

Here we introduced the $2\times 2$ matrix ${\bm{M}}^{e}_{\rm R}={\bm{G}}^{r}{\bm{\Gamma}}^{e}_{\rm R}{\bm{G}}^{a}$. It is also worth noting that the individual elements $\Gamma^{e}_{\rm L}M^{e}_{{\rm R}ij}$ in the above combination are from the rate components $\Gamma^{\pm}_{{\rm L}ij}$ in the master equation, i.e., $\Gamma^{\pm}_{{\rm L}ij}=\Gamma^{e}_{\rm L}\int d\omega n^{\pm}_{{\rm L}}(\omega)A_{ij}(\omega)$.

Finally, we remark that the electron- and hole-type rates, ${\bm{\Gamma}}^{e/h}_{\alpha}(\pm\omega)$, are originated from the BdG-type consideration for the original electrons in the leads. Owing to the particle-hole symmetry, the energy replacement of $\omega\to-\omega$ should be made when considering the transformation from an electron to a hole. Therefore, from the original rates, we redefined ${\bm{\Gamma}}_{\alpha}(\omega)={\bm{\Gamma}}^{e}_{\alpha}$ and ${\bm{\Gamma}}_{\alpha}(-\omega)={\bm{\Gamma}}^{h}_{\alpha}$. Moreover, associated with this conversion, the electron and hole occupation function are different, i.e., the former is $n^{+}_{\alpha}(\omega)$ and the latter is $n^{-}_{\alpha}(-\omega)$, as one can find in Eq. (III.3).

In this context, based on Eq. (III.3) and in particular its construction as exemplified by Eq. (III.3) (i.e., the combination of multiple channels), we briefly discuss the issue of teleportation, in an attempt to shed light on the reason of the vanishing of the teleportation channel when $\varepsilon_{\mbox{\tiny M}}\to 0$. To see this more clearly, let us reexpress the rate matrices in Eq. (27) in terms of the projection operators as

	$\displaystyle{\bm{\Gamma}}_{\rm L}$	$\displaystyle=$	$\displaystyle 2\Gamma_{\rm L}|\Phi_{\rm L}\rangle\langle\Phi_{\rm L}|\,,$
	$\displaystyle{\bm{\Gamma}}_{\rm R}$	$\displaystyle=$	$\displaystyle 2\Gamma_{\rm R}|\Phi_{\rm R}\rangle\langle\Phi_{\rm R}|\,,$		(32)

where $|\Phi_{\rm L/R}\rangle=(|1\rangle\pm|2\rangle)/\sqrt{2}$, with $|1\rangle$ and $|2\rangle$ the basis states which support the expression of the matrices in Eq. (27). In terms of the BdG formalism language, they correspond to the positive and negative energy states $|E_{0}\rangle$ and $|-E_{0}\rangle$, respectively, with $E_{0}=\varepsilon_{\mbox{\tiny M}}$ the energy of the regular fermion (the $f$-particle) associated with the MZMs. Accordingly, the states $|\Phi_{\rm L}\rangle$ and $\Phi_{\rm R}\rangle$ correspond to the left and right Majorana modes.

Based on Eq. (III.3), one can examine that ${\bm{\Gamma}}_{\rm L}{\bm{G}}^{r}(\omega){\bm{\Gamma}}_{\rm R}\to 0$ at the limit $\varepsilon_{\mbox{\tiny M}}\to 0$, which is valid for all the $ee$, $hh$, $eh$ and $he$ transmission components. Actually, using the explicit solution of Eq. (17), we have $\langle\Phi_{\rm L}|{\bm{u}}|\Phi_{\rm R}\rangle=u_{11}-u_{22}-u_{12}+u_{21}=2\,i\,\varepsilon_{\mbox{\tiny M}}/Z$, while noting that ${\bm{G}}^{r}=-i{\bm{u}}$. This is the remarkable issue of teleportation vanishing when the Majorana coupling energy approaches to zero, which implies that all the electron transmission and the crossed Andreev-reflection process through a pair of MZMs vanish at the limit $\varepsilon_{\mbox{\tiny M}}\to 0$.

Here, we would like to point out that this result is a consequence of the combination of multiple transmission channels. It does not involve only a single particle transmission, but involves several different charge configurations. For single transmission process, the vanishing of teleportation does not take place when $\varepsilon_{\mbox{\tiny M}}\to 0$. For instance, let us look at

	$\displaystyle\Gamma^{e}_{\rm L}M^{e}_{{\rm R}11}=\Gamma^{e}_{\rm L}\langle 1|{\bm{G}}^{r}{\bm{\Gamma}}^{e}_{\rm R}{\bm{G}}^{a}|1\rangle$
	$\displaystyle=2\Gamma^{e}_{\rm L}\Gamma^{e}_{\rm R}\langle 1|{\bm{G}}^{r}|\Phi_{\rm R}\rangle\langle\Phi_{\rm R}|{\bm{G}}^{a}|1\rangle\,.$		(33)

It does not vanish as $\varepsilon_{\mbox{\tiny M}}\to 0$. Therefor, the result that an electron cannot transmit through a pair of MZMs at the limit $\varepsilon_{\mbox{\tiny M}}\to 0$ is an effective picture owing to combination of multiple processes. Under large bias condition (i.e., with the bias voltage larger than the zero-energy-level broadening), the multiple processes cannot take place simultaneously. Then, one should not expect the vanishing phenomenon of teleportation.

The result of transmission coefficients given by Eq. (III.3) can be connected with the stationary $S$-matrix scattering approach as follows. Let us introduce the $S$-matrix (operator) as

$\displaystyle{\bm{S}}_{\rm LR}(\omega)={\bm{W}}_{\rm L}{\bm{G}}^{r}(\omega){\bm{W}}^{\dagger}_{\rm R}\,,$

(34)

where the two coupling operators between the left/right Majorana modes and the left/right leads are give by

	$\displaystyle{\bm{W}}_{\rm L}$	$\displaystyle=$	$\displaystyle\sum_{p=e,h}\sqrt{\Gamma_{\rm L}}\,|p_{\rm L}\rangle\langle\Phi_{\rm L}|\,,$
	$\displaystyle{\bm{W}}^{\dagger}_{\rm R}$	$\displaystyle=$	$\displaystyle\sum_{q=e,h}\sqrt{\Gamma_{\rm R}}\,|\Phi_{\rm R}\rangle\langle q_{\rm R}|\,.$		(35)

Here, $p$ and $q$ are introduced to denote the electron and hole states in the left and right leads, respectively. Then, the coefficients of transmission from the left to right leads can be reexpressed in terms of the $S$-matrix elements as

$\displaystyle{\cal T}^{pq}_{\rm LR}(\omega)=|\langle p_{\rm L}|{\bm{S}}_{\rm LR}(\omega)|q_{\rm R}\rangle|^{2}\,.$

(36)

All the other coefficients (for the local Andreev reflection) in Eq. (III.3) can be similarly formulated by means of the $S$-matrix scattering approach.