Master equation approach for transport through Majorana zero modes

Based on an exact formulation, we present a master equation approach to transport through Majorana zero modes (MZMs). Within the master equation treatment, the occupation dynamics of the regular fermion associated with the MZMs holds a quite different picture from the BdG S-matrix scattering process, in which the"positive"and"negative"energy states are employed, while the master equation treatment does not involve them at all. Via careful analysis for the structure of the rates and the rate processes governed by the master equation, we reveal the intrinsic connection between both approaches. This connection enables us to better understand the confusing issue of teleportation when the Majorana coupling vanishes. We illustrate the behaviors of transient rates, occupation dynamics and currents. Through the bias voltage dependence, we also show the Markovian condition for the rates, which can extremely simplify the applications in practice. As future perspective, the master equation approach developed in this work can be applied to study important time-dependent phenomena such as photon-assisted tunneling through the MZMs and modulation effect of the Majorana coupling energy.

In addition to the steady-state transport study for such as the zero-bias conductance peak, future study may turn more attention to dynamical aspects associated with the MZMs, which may include problems such as motionally moving the MZMs, modulating the coupling between the MZMs, and photon-assisted tunneling through MZMs, etc. In particular, one may consider to insert all these studies into a transport context, which is anticipated to make the unique dynamics of the MZMs measurable through time-dependent transport currents. For isolated Majorana systems, the theoretical tool of * Electronic address: jsjin@hznu.edu.cn † Electronic address: xinqi.li@tju.edu.cn quantum dynamics simulation can be the time-dependent Schrödinger equation. However, for the open system of MZMs coupled to outside environments (e.g. fluctuating background charges and/or transport leads), it is required, usually, to develop a master equation approach such as done in the recent works [21,22], where the dissipative quantum dynamics caused by local noises, transient evolution of the zero-bias conductance peak, and finite bandwidth effect of the leads are numerically investigated.
Indeed, master equation approach is an alternative useful tool for studying transport through mesoscopic systems [23][24][25][26][27][28][29][30][31][32], aside from the well-known Landauer-Büttiker scattering matrix theory [33] and the nonequilibrium Green's function formalism [34]. One of the distinct advantages of the master equation approach is its convenience for studying the various statistical properties of transport currents. In transport through the MZMs, the charge transmission channels contain normal and Andreev processes. In steady state, the Bogoliubov-de Gennes (BdG) S-matrix scattering approach can describe the various channels very clearly. However, in the master equation approach, which describes the occupation dynamics and allows for current calculation, it is unclear how to infer these transport channels. In this work, we carry out explicit results for the master equation and currents of transport through a pair of MZMs. In particular, via careful analysis for the structure of the rates involved in the master equation, we will reveal the intrinsic connection between the master equation and the BdG S-matrix scattering approach. To the best of our knowledge, this important connection is absent in literature. This connection can also help us to better understand the confusing issue of teleportation when the Majorana coupling vanishes.
IG. 1: (Color online) Sketch of the two-lead (three-terminal) setup for transport through a pair of MZMs (γL and γR, with coupling energy εM). The two leads are biased by voltages VL and VR with respect to the Fermi level of the superconductor, while the superconductor is grounded through the third terminal.
The work is organized as follows. In Sec. II we present the main results derived in Appendix A, including the master equation and time-dependent currents. In Sec. III we carry out numerical results for the transient rates and currents, illustrate the bias condition for Markovian approximation, and analyze the steady-state results in detail by a connection with the BdG S-matrix scattering approach. Finally, summary and discussion are presented in Sec. IV.

II. FORMULATION
For the transport setup of a pair of MZMs coupled to transport leads, the total Hamiltonian consists of three parts, The cental TS quantum wire can be described by the low-energy effective Hamiltonian of a pair of MZMs, where ε M is the coupling energy between γ L and γ R . The two Majorana fermions are related to a regular fermion through the transformation of The transport leads can be modeled by the noninteracting electron Hamiltonian as H B = αk ε αk c † αk c αk , with c † αk (c αk ) the creation (annihilation) operators of the electron in the αlead (α = L/R, the left/right lead). The tunnel-coupling between the MZMs and the two leads is described by [7,35,36], Throughout this work, we set = 1, unless otherwise stated.

A. Exact Formulation of master equation
Let us rewrite the tunnel-coupling Hamiltonian Eq. (2) as Here we denoted For the convenience of latter presentation, let us also introduce the correlation functions g + where · · · B means the thermal average over the reservoir states of the leads. Straightforwardly, one can more explicitly obtain g ± , where n ± α are the Fermi occupied and empty functions. The sum of these two functions reads as g α (t − τ ) = k e −iε αk (t−τ ) |t αk | 2 . Through the Fourier transformation we also introduce the corresponding rates Γ α (ω) = 2π k |t αk | 2 δ(ω − ε αk ) and Γ ± α (ω) = Γ α (ω)n ± α (ω). The basic idea of master equation approach to quantum transport is regarding the transport leads as a generalized environment and constructing an equation-ofmotion for the reduced state of the central device system, ρ(t) ≡ tr B ρ T (t), where the trace is over the states of the transport leads and ρ 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 ρ T (t), ρ T (t) = −i[H T , ρ T (t)]. Performing the trace tr B (· · · ) on two sides of the Liouvillian equation, one obtainṡ where the two auxiliary density operators (ADOs) are introduced as By means of a pathintegral formulation in the fermionic coherent state representation [32], it is possible to express the ADOs ρ ± α (t) in terms of the reduced density operator ρ(t). For the MZMs system under study, the detailed derivation is presented in Appendix A and the final result reads aṡ where the Lindblad superoperator is defined by D[A]ρ = AρA † − 1 2 {A † A, ρ} and the time-dependent rates are given by Assuming a wide-band limit approximation for the transport leads, i.e., Γ α (ω) → Γ α in Eq. (4), the various individual rates can be more simply expressed as while the other rates can be obtained through the following relations In the above results, the reduced propagators (RPs) for the central MZMs system u ij (t, τ ), in terms of a matrix form u(t, τ ), satisfy the following equation-of-motion Here we also considered the wide-band limit and introduced that Γ = Γ L + Γ R and Γ d = Γ L − Γ R . Actually, 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 † 11 (t, t 0 ) = u 22 (t, t 0 ) and u † 12 (t, t 0 ) = u 21 (t, t 0 ). Owing to the symmetry properties of the u(t, t 0 ) matrix, from Eq. (9c) and (9d), one can prove that Γ + α12 (t) = −Γ − α21 (t). Together with Eq. (9h), we then have Γ ± α12 (t) = Γ ± α21 (t). Remarkably, for the master equation, the total rates Γ 1 (t) and Γ 2 (t) of Eq. (8) are determined only by the sum of the diagonal rates Γ ± α11/22 (t).
The off-diagonal rates Γ ± α12/21 (t) have no effect on the state evolution. However, as we will see later, the offdiagonal rates have important contribution to the transport currents.
So far, we have presented the main results of the master equation, i.e., Eq. (7) 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 → ∞) for the transient rates given by Eqs. (9). We obtain Here, precisely following the Green's function theory, we defined the spectral function through while u(ω) is the Fourier transformation of the Green's function matrix u(t) in time domain, i.e., u(ω) = ∞ 0 dte iωt u(t). Here the initial condition of u(t) has been considered, which makes the Fourier and Laplace transformations identical. From Eq. (10), performing a Laplace transformation, we obtain the solution in frequency domain as Then, the spectral functions are obtained as In general, the spectral functions hold the properties of dω A 11/22 (ω) = 1 and dω A 12 (ω) = 0.

B. Time-Dependent Currents
In master equation approach, the transport currents should be calculated using the reduced density matrix ρ(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 nonequilibrium 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Î α . Using the Heisenberg equation, we haveÎ α = −e d dtN α = ie[N α , H T ], where the total electron number operator of lead-α reads asN α = k c † αk c αk . Then, we obtain The current in lead-α is where the average is over the total states of the whole system-plus-leads, i.e., Tr T ≡ tr S tr B . Further, we reexpress the current in terms of the ADOs As mentioned above and demonstrated in Appendix A, the ADOs ρ ± α (t) can be expressed in terms of the reduced density matrix ρ(t), see Eq. (A24). Then, the current is obtained as (17) with the various component currents explicitly given by In deriving this result, the elements of the density matrix were defined through ρ . In practice, they should be solved from Eq. (7). The various rates, Γ ± αij (t), are those given by Eq. (9).
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 α11 (t), the first term is associated with normal tunneling process of an electron from the α-lead to the MZMs quasiparticle, while the second term describes the inverse process of the first one. (ii) In the current I α22 (t), the first term describes the Andreev reflection (AR) process of an electron entering the MZMs from the α-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 α12 (t), the first term describes another contribution of normal tunneling process of an electron from the α-lead to the MZMs quasiparticle, while the second term is the inverse process of the first one. (iv) In I α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.

A. Transient Rates and Currents
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. 13, 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., µ L /e = µ 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 Γ ± αij (t → ∞) 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 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 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 component currents. As expected, even using the asymptotic rates, the switching-on transient behavior of I L11 (t) and I L22 (t) is obvious, in contrast to the total current I L (t). We notice that, at the very beginning stage, I L22 (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 L22 (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. Parameters used are the same as in Fig. 2.

B. Markovian Limit
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 η = Γ L /Γ R . Owing to the a few symmetry properties of the rates Γ ± αij , here we only plot the results of Γ +

L11/22
and Γ + L12 . In Fig. 3(a), for the diagonal rates Γ + L11/22 , we find that they increase with the bias voltage V (from negative to positive), while at the large ±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 ±V . This behavior is determined by the following properties of the off-diagonal element of the spectral function matrix from Eq. (14c), i.e., A 12 (ω) = A 12 (−ω) and them, we can easily prove that Γ + L12 (V ) + Γ + L12 (−V ) = 0. In Fig. 3(c) and (d), we show the total rates Γ 1 and Γ 2 which characterize, respectively, the annihilation and creation rates of the MZMs-quasiparticle. For ε M = 0, we find that the both total rates are bias (V ) independent; while for ε M = 0, they reveal a symmetric dependence on ±V and take maximum and minimum at V = 0 for Γ 1 and Γ 2 , respectively. The reason for this behavior is as follows. First, based on Eq. (14c), we know A 12 (ω) = A 21 (ω). Then, we can prove that Γ − L12 (V ) + Γ + L12 (V ) = Γ L ∞ −∞ dωA 12 (ω) = 0. Second, for ε M = 0, we know that the peaks of A 11/22 (ω) are centered at ±ε M . Then, we know that Γ + L11 (V ) increases with V later than Γ + L22 (V ), when crossing the region around zero bias (from negative to positive). Because of the same reason, Γ − L11 (V ) decreases with V earlier than Γ − L22 (V ), when crossing the region around zero bias. As a result, the total rates, , 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 con-figurations. This is because the rates Γ ± αij are dominantly affected by the Fermi occupation of the α-lead electrons, while influence of the opposite lead is weakly mediated through the propagating function u(ω) via the self-energy effect.

C. Steady-State Currents: Connection with Other Approaches
Let us introduce the more conventional (retarded) superconductor Green's functions From the equation-of-motion of u(t, τ ), i.e., Eq(..), we know that u(t, τ ) = iG r (t, τ ). Then, the spectral function matrix can expressed as From the Dyson equation, we have The retarded self-energy matrix is given by while the rate matrices read as As to be clear in the following, we may denote Γ α (ω) ≡ Γ e α and Γ α (−ω) ≡ Γ h α .
Based on Eq. (18), i.e., the time-dependent current formula of the master equation approach, the nonequilibrium Green's function results of steady-state currents can be easily obtained. As an example, let us consider the steady-state current I L11 = Γ + L11 − Γ Lρ11 . We have Here we have applied the relationship between ρ 11 (t) = tr After similar treatment to all the other component currents and applying Eq. (21) and the well known Keldysh equation for G < (ω), combination (together also with a procedure of symmetrization) of the component currents yields In this result, the various transmission coefficients are given by The meaning of these transmission coefficients is clear: T eh(he) 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; T ee(hh) LR describes the normal transport process of an electron (or a hole) from the left to the right lead; and T eh(he) 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 L = i,j=1,2 I Lij , we have performed the following combination of elements from the four component currents Here we introduced the 2 × 2 matrix M e R = G r Γ e R G a . It is also worth noting that the individual elements Γ e L M e Rij in the above combination are from the rate components Γ ± Lij in the master equation, i.e., Γ ± Lij = Γ e L dωn ± L (ω)A ij (ω). Finally, we remark that the electron-and hole-type rates, Γ e/h α (±ω), are originated from the BdG-type consideration for the original electrons in the leads. Owing to the particle-hole symmetry, the energy replacement of ω → −ω should be made when considering the transformation from an electron to a hole. Therefore, from the original rates, we redefined Γ α (ω) = Γ e α and Γ α (−ω) = Γ h α . Moreover, associated with this conversion, the electron and hole occupation function are different, i.e., the former is n + α (ω) and the latter is n − α (−ω), as one can find in Eq. (25).
In this context, based on Eq. (26) and in particular its construction as exemplified by Eq. (27) (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 ε M → 0. To see this more clearly, let us reexpress the rate matrices in Eq. (23) in terms of the projection operators as where |Φ L/R = (|1 ± |2 )/ √ 2, with |1 and |2 the basis states which support the expression of the matrices in Eq. (23). In terms of the BdG formalism language, they correspond to the positive and negative energy states |E 0 and |−E 0 , respectively, with E 0 = ε M the energy of the regular fermion (the f -particle) associated with the MZMs. Accordingly, the states |Φ L and Φ R correspond to the left and right Majorana modes.
Based on Eq. (28), one can examine that Γ L G r (ω)Γ R → 0 at the limit ε M → 0, which is valid for all the ee, hh, eh and he transmission components. Actually, using the explicit solution of Eq. (13), we have Φ L |u|Φ R = u 11 − u 22 − u 12 + u 21 = 2 i ε M /Z, while noting that G r = −iu. 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 Andreevreflection process through a pair of MZMs vanish at the limit ε M → 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 ε M → 0. For instance, let us look at It does not vanish as ε M → 0. Therefor, the result that an electron cannot transmit through a pair of MZMs at the limit ε M → 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. (26) can be connected with the stationary S-matrix scattering approach as follows. Let us introduce the S-matrix (operator) as where the two coupling operators between the left/right Majorana modes and the left/right leads are give by 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 All the other coefficients (for the local Andreev reflection) in Eq. (26) can be similarly formulated by means of the S-matrix scattering approach.

IV. SUMMARY AND DISCUSSION
We have presented a master equation approach for transport through MZMs. The master equation approach will be convenient for studying the various statistical properties of transport currents and important time-dependent phenomena such as photon-assisted tunneling through the MZMs and modulation effect of the Majorana coupling energy. In this initial work, we carried out explicit results for the master equation and currents. We illustrated the behaviors of transient rates, occupation dynamics and currents and exploited the Markovian condition for the rates, which can extremely simplify the applications in practice. Via careful analysis for the structure of the rates, we also revealed the intrinsic connection between the master equation and the BdG S-matrix scattering approach, by noting that the charge transfer picture and occupation dynamics in both formulations are quite different. This connection is absent in literature and can help us to better understand the issue of Majorana teleportation.
About the possible application to time-depend transports, the simplest case is to consider modulating the central system under Markovian approximation for the leads. In this case, in the master equation, only the central system Hamiltonian in the coherent term becomes time dependent, and the rates in the dissipative terms are unaffected by the time-dependent modulation of the central system. Beyond Markovian approximation, more rigorous treatment for the rates needs to consider the effect of modulation on the reduced propagating function u(t, τ ), for the central system of MZMs, whose equationof-motion is given by Eq. (10). If we consider a modulation to the coupling energy between the MZMs, we only need to insert the time-dependent ε M (t) into Eq. (10) to solve u(t, τ ) and calculate the time-dependent rates. Then, solving the master equation allows us to obtain the time-dependent currents.
An important example of time-dependent transport is to consider adding an ac voltage to the dc bias. In this case, rather than treating time-dependent chemical potentials for the leads, one can alternatively keep the occupation of the lead electrons (labelled by "αk") unchanged (determined by the dc bias), but at the same time let the energy levels of the lead electrons be time dependent, i.e., ε αk (t) = ε αk + ∆ α (t). In this treatment, we only need to perform a replacement for the correlation functions of the lead electrons (g α , g ± α ) → (g α ,g ± α ), while the latter read as where g α (t − τ ) and g ± α (t − τ ) are given by Eq. (4). Inserting these extended functions into the expressions of the rates, one can straightforwardly solve the master equation and calculate the time-dependent currents. Obviously, this formulation works for arbitrary timedependent voltages, such as step-like modulation and rectangular bias pulse.
Noting that the path integral in Eq. (A3) is of the quadratic type, it is thus possible to carry out the exact result. Applying the variational calculus, one can obtain the "classical trajectory" (i.e. the extremal trajectory); and the exponential factor of the path integral result is given by the action determined by the "classical trajectory". We have where N (t) is the normalization factor and ξ(t), ξ * (t 0 ), ξ ′ (t 0 ), ξ ′ * (t) are determined by the "classical trajectory" which obeys the equations-of-motion