Dynamics of hybrid junctions of Majorana wires
Nilanjan Bondyopadhaya Dibyendu Roy












=227mm


    [ ∂̣|̌ł⟨⟩̊δ1 -1mm l^1 Integrated Science Education and Research Centre,  Visva-Bharati University, Santiniketan, WB 731235,  India^2 Raman Research Institute, Bangalore 560080, India
        We investigate the dynamics of hybrid junctions made of the topological superconductor (TS) and normal metal (N) wires. We consider an X-Y-Z configuration for the junctions where X, Y, Z = TS, N. We assume the wires X and Z being semi-infinite and in thermal equilibrium. We connect the wires X and Z through the short Y wire at some time, and numerically study time-evolution of the full device. For TS-N-TS device, we find a persistent, oscillating electrical current at both junctions even when there is no thermal or voltage bias, and the TS wires are identical. The amplitude and period of the oscillating current depend on the initial conditions of the middle N wire indicating the absence of thermalization. This zero-bias current vanishes at a long time for any of X and Z being an N wire or a TS wire near a topological phase transition. Using properties of different bound states within the superconducting gap, we develop a clear understanding of the oscillating currents. 
        
         
        
        
        
        Dynamics of hybrid junctions of Majorana wires
            Nilanjan Bondyopadhaya^1 and Dibyendu Roy^2
            December 30, 2023
        ===============================================
        
        
        
        Recent progress in search of Majorana zero modes (MZMs) in condensed matter systems has generated massive interest for better understanding, control, and engineering of such systems featuring these exotic quasiparticles <cit.>. The MZMs have been theoretically proposed to emerge at the edges of a topological superconductor (TS) wire and to exhibit non-Abelian particle statistics <cit.>. Several experiments have observed significant evidence of MZMs in electrical transport measurements with engineered TS wires <cit.>. These TS wires hosting Majorana quasiparticles are expected to be an essential component of future quantum devices such as fault-tolerant quantum computers. Therefore, it is necessary to investigate various properties including thermalization, transport and braiding operations in hybrid junctions of TS and normal metal (N) wires.
        
        In the past ten years, there are many theoretical studies for electrical transport in different junctions of TS and N wires <cit.>. Most of these studies discuss steady-state transport to derive zero-bias conductance <cit.>, current-voltage characteristics <cit.>, equilibrium current-phase relation for fractional Josephson effect <cit.> and robust Majorana conductance peaks <cit.>. There are also some time-evolution study after a quantum quench <cit.> such as when an N wire is suddenly connected to a TS wire. Interestingly, a systematic study for the dynamics of tunneling current in these hybrid devices is missing. In this Letter, we study transient and steady-state quantum transport in hybrid junctions of X-Y-Z configuration where X, Y, Z = TS, N. We primarily address route to equilibration in these junctions by investigating the dynamics of tunneling current. Ideally, we consider the wires X and Z to be semi-infinite and in thermal equilibrium. However, we can only investigate L_X,Z/L_Y >>1  in our numerics for time-evolution where L_X,Z and L_Y are the length of the wire X or Z and Y respectively. We connect the wires X and Z through the middle Y wire at some time t_0 and explore time-evolution of the full system for various initialization of the middle wire. We probe electrical current through the X-Y and Y-Z junctions at every stage of time-evolution. Here we do not consider a voltage bias which can be incorporated through time-dependent phases in the tunnel couplings across the junctions <cit.>. We also ignore any effect which can arise due to phase difference across the junctions of TS wires; therefore we set the phase of the TS wires to be same (or zero for simplicity) everywhere. 
        
        We find a non-decaying, oscillating electrical current at both junctions of a TS-N-TS device even in the absence of any thermal or voltage bias and for identical TS wires (no phase difference). The amplitude and period of the oscillating current depend on the energy eigenstates of the wires which are localized around the junctions and within the superconducting gap of the TS wires, and on the initial conditions of the N wire. Since the wave functions of the energy eigenstates of the N wire in the superconducting gap decay exponentially deep inside the TS wire, they act as bound states. The presence of such bound states prevent equilibration (thermalization) of the middle N wire with the boundary wires <cit.>, and this results into the non-uniqueness (initial-condition dependence) of the persistent current. There is no zero-bias  oscillating current at long time in an N-TS-N and an X-N-Z device when any of the X and Z being an N wire or a TS wire near a topological phase transition. In these latter cases, the energy spectrum of one of the boundary wire is gapless. Thus, there is no more bound state from the middle N wire, and the middle wire gets equilibrated with the boundary wire(s). Using properties of different bound states within the superconducting gap in our devices, we develop a theory to reproduce the simulated oscillating current with little numerics. 
        
        
        We now introduce the Hamiltonian for various components of the hybrid junctions. We model the N wire as a noninteracting tight-binding chain of spinless fermions and the TS wire by the Kitaev chain of a spinless p-wave superconductor <cit.>. We write below a general Hamiltonian H_α which can be used to represent N and TS wires by tuning the parameters:  
        
        H_α   =   -γ_α∑_l'=1^L_α-1(c^_l'c_l'+1+c^_l'+1c_l')+ϵ_α∑_l'=1^L_αc^_l'c_l'
           +   Δ_α∑_l'=1^L_α-1(c^_l'c^_l'+1+c_l'+1c_l'),
        where c_l' (c^_l') indicates annihilation (creation) operator of a spinless fermion at site l' of the wire segment α=X,Y,Z. Here, the parameters γ_α,ϵ_α and Δ_α denote respectively hopping, on-site energy and superconducting pairing energy, and we take them to be real. The Hamiltonian H_α indicates a N wire in the absence of pairing (Δ_α=0). In the presence of pairing Δ_α 0, the superconducting wire undergoes a topological phase transition as ϵ_α is increased across 2γ_α. The wire is in a topological phase for ϵ_α<2γ_α and the TS wire hosts two MZMs at the opposite ends of the wire for a relatively long wire. The wire is in a topologically trivial phase without the MZMs for ϵ_α>2γ_α. The topological phase transition near ϵ_α=2γ_α is also accompanied by a bulk-gap closing. The superconducting wire has a large bulk-gap in its spectrum both in the topologically non-trivial and trivial phases, and the gap vanishes at the topological phase transition around ϵ_α=2γ_α.  
        
        
        The tunneling Hamiltonian for the X-Y and Y-Z junction is independent of nature of X, Y, Z, and we take it of the following form:
        
        H_αβ=-γ_αβ(c_l'^c_l'+1+c_l'+1^c_l'),
        
        where αβ= XY, l'=L_X, and αβ= YZ, l'=L_X+L_Y respectively for the X-Y and Y-Z junction. For simplicity, we assume here the tunneling rate γ_αβ (with γ_αβ≪γ_α, γ_β) to be same for both the junctions, i.e., γ_XY=γ_YZ=γ'. The full Hamiltonian of the hybrid device is H_ F=H_ X+H_ Y+H_ Z+H_ XY+H_ YZ. The full device consists of L=L_X+L_Y+L_Z number of fermionic sites. Using conservation of electrical charges across the junctions, we define the electrical current at the junctions as
        
        J_αβ=-iγ_αβ⟨ (c_l'^c_l'+1-c_l'+1^c_l')⟩,
        where again αβ= XY, l'=L_X, and αβ= YZ, l'=L_X+L_Y respectively for the X-Y and Y-Z junction. The expectation ⟨ .. ⟩ defines averaging over the initial states of the wires.
        
        
        Due to the pairing term Δ_α in the superconducting parts of the hybrid device, it is convenient to use a basis a≡(a_1, a_2, …, a_2L-1, a_2L)^T=(c_1, c^_1, …, c_L, c^_L)^T to write the quadratic Hamiltonian H_ F= a^ℋ^ F a=∑_l,mℋ_lm^ Fa^_la_m in terms of the matrix ℋ_lm^ F. Thus, a_2l=a^_2l-1. Clearly the index l in a_l (or a_l^†) does not represent actual physical site of the wire. 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. We consider the wires X and Z are in thermal equilibrium at temperatures T_ X,Z and chemical potentials μ_ X,Z before we connect them through the middle wire Y at time t_0. For isolated X and Z  wires, we have H_α=∑_l,mℋ_lm^αa^_la_m with l,m=1,2,…,2L_ X for α= X, and l,m=2(L_X+L_Y)+1,2(L_X+L_Y)+2,…,2L for α= Z. Therefore, we can write
        ∑_mℋ_lm^αψ^α_q(m)=λ^α_qψ^α_q(l),
        
        where ψ^α_q(m) and λ^α_q represent the eigenvectors and eigenvalues of the wire α=X,Z. Thus, the equilibrium density matrix for sites on the isolated boundary wires 
        ⟨ a^†_l(t_0) a_m(t_0) ⟩ = ∑_q ψ_q^α^*(l)ψ_q^α(m) f(λ^α_q,μ_α,T_α),
        with l,m=1,2,…,2L_ X for α= X and l,m=2(L_X+L_Y)+1,2(L_X+L_Y)+2,…,2L for α= Z. Here, f(λ^α_q,μ_α,T_α)=1/(e^(λ^α_q- μ_α)/k_BT_α+1) is the Fermi function.
        
        We also assume that the operators from different X, Y, Z wires are uncorrelated when they are disconnected at t_0. Therefore, we have 
              ⟨ a^†_l(t_0) a_m(t_0) ⟩ =⟨ a_l(t_0) a^†_m(t_0)⟩
              =⟨ a_l(t_0) a_m(t_0) ⟩=⟨ a_l^†(t_0) a^†_m(t_0)⟩=0  ,
        where l and m represent indices corresponding to two different wires of the hybrid device. We here do not take a thermal distribution for the middle wire Y at t_0. Instead, we choose some arbitrary initial density matrix of wire Y, such as, 
        ⟨ a_2l'(t_0) a_2m'-1(t_0) ⟩ = {[   n_l'  l'=m';      0 l'  m' ]. 
        for physical sites: l',m' ∈{L_ X+1,…, L_ X+L_ Y}. We wish to check whether the density matrix of the full device in the long-time limit t→∞ becomes independent of the initial density matrix of the middle wire. Such independence would indicate thermalization of the Y wire by the X and Z wires. We carry out this job by calculating time-evolution of electrical current at the junctions. So, we connect the wires by the tunneling Hamiltonians at time t_0 and study the time-evolution of the full device using the Heisenberg equations of motion. The solution of these equations of motion is given by
         a(t)=i 𝒢^+(t-t_0) a(t_0),
        
        where 𝒢^+(τ)=-i e^-i H_ Fτθ (τ) is the Green's function of the full device. Here, θ (τ) is the Heaviside step function. Suppose, Ψ_q(m) and Λ_q denote the eigenvectors and eigenvalues of the full Hamiltonian H_ F. Therefore, 
        ∑_m=1^2Lℋ^ F_lmΨ_q(m)=Λ_q Ψ_q(l),  l=1,2,…,2L.
        
        We apply the above relations to write the full Green's function in the following form for t >t_0:
        𝒢^+_r,s(t-t_0)=- i ∑_q=1^2LΨ_q(r) Ψ_q^*(s) e^- i Λ_q (t-t_0)  ,
        
        where r,s ∈{1,…, 2L}. Now, we can write the time-evolved density matrix of the full device as
        ⟨ a^†_l(t) a_m(t) ⟩= ∑_r,s=1^2L𝒢^+_m,s(t-t_0) ⟨ a^†_r(t_0) a_s(t_0) ⟩ [𝒢^+_l,r(t-t_0)]^†,
        
        where we plug the initial density matrix ⟨ a^†_r(t_0) a_s(t_0) ⟩ from Eqs. <ref>,<ref>,<ref>. Using Eq. <ref> in Eq. <ref>, we evaluate the time evolution of the electrical current at the junctions.
        
        In the following numerical analysis, we fix γ_ X=γ_ Z=γ >γ_ Y; so that the bands of the boundary wires are broader than that of the middle. An N-N-N device in the above set-up has been studied in details in Ref. <cit.>, and it has been shown that there is a unique nonequilibrium steady-state (independent of initial values of n_l') in the device when there is no single particle bound state from the middle wire (the band of wire Y lies within that of X and/or Z).
        
        First, we investigate dynamics of an N-TS-N device whose steady-state transport characteristics are extensively studied both theoretically <cit.> as well as experimentally <cit.> for detection of MZMs [An N-TS junction is only considered in many of these studies.]. We here prepare the decoupled N wires of our N-TS-N device initially in thermal equilibrium. We find from our numerics <cit.> that the nonequilibrium steady-state transport in the  N-TS-N device seems to be independent of the initial conditions of the finite TS wire when the band of the N wires is wider than that of the TS wire. Our present finding of unique nonequilibrium steady-state in the N-TS-N device validates all those previous steady-state transport analyses <cit.>  in this system.             
        
        
        
        Next, we consider a TS-N-TS device which has been investigated for the Josephson effect in topological systems <cit.>. We take the temperature of both the TS wires to be the same and very low, and set their chemical potential to zero. Then, one would naively expect zero electrical current in such junctions of identical TS wires at a long time. Surprisingly, we find a persistent and oscillating electrical current at both the junctions of the TS-N-TS device. In Fig. <ref>, we show time evolution of J_ XY(t) from t_0=0 for different initial density n_l' of the middle N wire. The amplitude and period of the current oscillation in Fig. <ref>(a,b) depend on the initial density of the N wire. For example, there is only a single oscillation period in Fig. <ref>(a) for initial density, n_1=n_2=n_3=0, while there are two oscillation periods (a short time and a long time) in Fig. <ref>(b) for n_1=0,n_2=n_3=1. Thus, our TS-N-TS device does not equilibrate. 
        
         The absence of unique long-time steady-state (equilibration) in a TS-N-TS device is due to the presence of bound states near the junctions. The wave functions of these bound states decay exponentially deep inside the TS wires. The energy of these bound states lies within the superconducting gap of the TS wires. Apart from the Majorana bound states of the TS wires near the intersections, there can be such bound states originating from the N wire. We observe that the amplitude and period of the current oscillation in the TS-N-TS device also depend on the total number of such bound states near the junctions <cit.>. For example, there are two Majorana bound states with energy ± 0.12236 at the edges of TS wires near the intersections as well as two energy eigenstates (at energy ± 0.16795) from the middle N wire inside the superconducting gap for the TS-N-TS device in Fig. <ref>. The energy of MZMs near the junctions of the TS wires is non-zero due to an overlap of these modes through the short middle wire.  From these bound state energies, we derive the periods of short and long time oscillations respectively as 2π/(0.16795+0.12236)≈ 21.6 and 2π/(0.16795-0.12236)≈ 137.8 which are in good agreement with the simulated periods in Fig. <ref>.      
        
        
        
        
        To further illustrate the role of these bound states in the absence of equilibration, we study time-evolution in a TS-N-Z device where we either apply an N wire for Z or tune the on-site energy of a superconducting Z wire to sweep through a topological phase transition. We do not find any persistent, oscillating electrical current in a TS-N-N device in the absence of bias. The band of a semi-infinite N wire in our model is continuum without a bulk-gap. Therefore, the Majorana bound state of the TS wire near the junctions as well as the energy levels of the middle N wire do equilibrate with the boundary N wire; thus there is no oscillating current <cit.>. 
        
        
        We find persistent, oscillating currents (e.g., check Fig. <ref>(a)) in a TS-N-Z device in the absence of voltage or thermal bias when ϵ_ Z<2γ_ Z, such that there are Majorana bound states at the edges of the topological Z wire, and the bulk-gap (∼ |ϵ_ Z-2γ_ Z|) in the energy spectrum of the Z wire is much larger than γ_ YZ^2– dissipation induced by the coupling between Y and Z wires. Again, we can separately estimate the periods of oscillating currents in  Fig. <ref>(a) using the bound state energies, and they are 2π/(0.15455+0.11557)≈ 23.3 and 2π/(0.15455-0.11557)≈ 161.2. Regardless of the initial density of the middle N wire, there is no persistent, oscillating electrical current at the long time in a TS-N-Z device for ϵ_ Z≈ 2γ_ Z when there is either no bulk-gap in the energy spectrum or a bulk-gap which is smaller or same order of γ_ YZ^2. In Fig. <ref>(b), we show rapid decays of J_ XY(t) and J_ YZ(t) to zero with time from initial time t_0=0 when ϵ_ Z=2γ_ Z. We also notice in Fig. <ref>(b) that the amplitude of J_ YZ(t) is much smaller than J_ XY(t) in the time duration when they are non-zero. The last fact indicates that such decays of J_ XY(t) and J_ YZ(t) are due to equilibration by the thermal Z wire which affects the YZ junction more than the XY junction. Finally, we show oscillating J_ XY(t) and J_ YZ(t) in Fig. <ref>(c) when the Z wire is in a topologically trivial phase with a large bulk-gap for ϵ_ Z>2γ_ Z. However, the amplitude of J_ YZ(t) is almost two order smaller than J_ XY(t) due to the large bulk-gap and the absence of MZM in the Z wire making the effective coupling between the Y and Z wires much smaller than that between the Y wire and the Majorana bound state in the TS wire. The pattern of oscillation in Fig. <ref>(c) differs from that in  Fig. <ref>(a) due to the presence of an extra Majorana bound state in the topological Z wire near the junction for Fig. <ref>(a). The energy of MZM near the junction of the left TS wire remains almost zero due to negligible hybridization when the Z wire is topologically trivial. The period of current oscillation in Fig. <ref>(c) is then estimated solely by the bound state energy from the middle N wire as 2π/0.14833 ≈ 42.4.   
        
        
        We learn from our above discussion for Figs. <ref>,<ref> that the persistent current oscillation in different hybrid devices in the absence of bias is due to the bound states and their initial preparation. We here show that we can reproduce the long-time persistent current oscillations in Figs. <ref>,<ref> by solely keeping the contributions of the  bound states in the current operator in Eq. <ref>. To perform such calculation, we apply truncation both in real space as well as in energy eigenstates in our computation which becomes numerically inexpensive. While we keep only energy eigenstates (bound states) which are localized near the junctions, we constraint the span of real space depending on the localization lengths of the bound states around the junctions. We use these truncation rules to Eqs. <ref>,<ref>, and employ those truncated quantities along with Eqs. <ref>,<ref> in Eq. <ref> to evaluate Eq. <ref>. The calculated currents from such truncated computation with bound states are shown in Fig. <ref> to compare them with the currents from the full simulation. They show an excellent match at long-time.     
        
        The fractional Josephson effect with a characteristic 4π periodic current-phase relation in topological superconductors has been investigated for an unambiguous detection of MZMs <cit.>. There are also many tunneling spectroscopy studies with N and superconducting tip to probe various magnetic and non-magnetic bound states in conventional and unconventional superconductors <cit.>. Our present dynamical study can in principle be tested in both these above set-ups. Our finding of lack of equilibration in the presence of bound states in a device with superconducting boundary wires with a bulk-gap poses challenge for direct detection of dc or ac fractional Josephson effect as well as electrical current measurements in such device. We conclude that the necessary condition for equilibration (or a unique steady-state) in these hybrid junctions is either the absence of bound states near the junctions or the band of one or both the boundary wires being continuum without a bulk-gap.    
        
        \begin* Acknowledgments.
        We are grateful to A. Dhar and D. Sen for discussion. NB acknowledges funding from DST-FIST programme. DR acknowledges funding from the DST, India via the Ramanujan Fellowship. 
        \end*
        
        
        
        
        
        
        
        
        
        
        
        
         Supplemental Material for Dynamics of hybrid junctions of Majorana wires
        
        Nilanjan Bondyopadhaya and Dibyendu Roy
        
        
        
        § N-TS-N DEVICE
        
        In the main text, we have mentioned that the nonequilibrium steady-state transport in an N-TS-N device seems to be independent of the initial conditions of the finite TS wire when the band of the N wires is wider than that of the TS wire. In Fig. <ref>, we present plots to support our claim. In Ref. <cit.>, the time-evolution of the many-electron wave function of an N-TS junction is studied after a quantum quench in which the N and TS wires of equal and extended lengths are connected at some initial time. The decoupled N and TS wires are initially prepared in their respective ground states, and the overlap between the time-evolved wave function with the initial state (the Loschmidt echo) is found to decay with a universal power law in time for large times after the quench. We here instead prepare the decoupled N wires of our N-TS-N device initially in thermal equilibrium. In Fig. <ref>, we show that the electrical currents at a long time  do not depend on the initial density of the middle TS wire. We find from our numerics that the currents at both the junctions reach a constant non-zero value when there is a voltage bias across the TS wire due to the chemical potential difference in the boundary N wires. The simulated steady-state current in Fig. <ref> from time-evolution of the Heisenberg equations of motion matches with that from the Fourier transform solution at steady-state <cit.>.
        
        
        
        § LONGER LENGTH OF N WIRE IN TS-N-TS DEVICE
        
        In Fig. <ref>, we show how the amplitude and period of the current oscillation in a TS-N-TS device change with an increasing length of the middle N wire. The number of bound states from the middle N wire can increase with increasing wire length. The current oscillations in the TS-N-TS device seem to survive for relatively long N wire length. We have seen finite persistent current oscillation in the junctions for L_ Y=50. However, the pattern of current oscillation, e.g., oscillation period, becomes a bit complicated in the presence of multiple bound states for a longer L_ Y. We present in Fig. <ref> plots for L_ Y=4,6 for zero initial density of the middle N wire. 
        
        
        
        § TS-N-N DEVICE
        
        We have got unique nonequilibrium steady-state in a TS-N-N device, and we do not see here any persistent electrical current oscillation in the absence of voltage or temperature bias. In Fig. <ref>, we plot electrical currents at both the junctions of a TS-N-N device for two different initial density in the middle N wire. We also apply a voltage bias across the middle N wire by using a non-zero chemical potential for the boundary N wire. We find that the electrical currents at both junctions are independent of the initial density in the middle N wire. The nonequilibrium steady-state currents are also the same at both intersections. Therefore, we conclude that the middle N wire in a TS-N-N device equilibrates with the boundary wires due to a continuous band of the boundary N wire.
        
        
        
        § ROLE OF TUNNEL COUPLING
        
        In the main text, we have discussed equilibration in a TS-N-Z device when the Z wire transits through a topological phase transition. The closing of superconducting bulk-gap near the topological phase transition is the reason for equilibration. In Fig. <ref>, we further show that the electrical current in a TS-N-Z device decays to zero at a long time when there is a bulk-gap in the Z wire which is the smaller or same order of γ_ YZ^2. The tunnel coupling γ_ YZ induces broadening of the energy of the bound states. Thus, the energy-broadened bound states near the junctions can overlap with the band of the Z wire when the bulk-gap is small. Therefore, we expect equilibration of bound states and decay of electrical currents when the ratio between the bulk-gap and γ_ YZ^2 is low. For example, this ratio is 0.16, 0.8, 5 respectively for three plots in Fig. <ref>(a,b,c). While the current J_ XY(t) in Fig. <ref>(a,b) show clear decay with time, the decay of J_ XY(t) is relatively slow in Fig. <ref>(c). Therefore, a stronger tunnel coupling can effectively reduce the bulk-gap in superconducting boundary wires, and hence can help in equilibration. Nevertheless, there can also be other mechanisms of equilibration such as electron-phonon interaction at the middle wire or irradiation of radiation  on the Josephson junctions for the measurement of Shapiro steps in current-voltage curves. The bound states can equilibrate when the phonon or the radiation at the junctions is wide-band such that the band of phonon or radiation closes the superconducting gaps.  
        
        
        
        § TRUNCATED NUMERICS WITH BOUND STATES
        
        We have seen in the main text that the bound states play a crucial role in the persistent electrical current oscillation in different hybrid junctions in the absence of bias. These bound states may be formed inside the middle wire as well as at the edges of the boundary TS wires. To study the contribution of these bound states to the time evolution of density matrix, let us first rewrite the full Green's function 𝒢^+_r,s(t-t_0) in the following way <cit.>:
              𝒢^+_r,s(t-t_0)= 𝒢^b+_r,s(t-t_0)+𝒢^c+_r,s(t-t_0) 
           =   -i∑_bΨ_b(r) Ψ_b^*(s) e^- i Λ_b (t-t_0)
        -i ∫ dξρ^c_rs(ξ)e^-i ξ (t-t_0)  ,
        
        where Ψ_b and Λ_b are the full system's eigenvectors and eigenvalues corresponding to the bound states localized around the junctions, and the density matrix ρ_rs^c is given by a sum over extended (continuum) states Ψ_q of the system ρ^c_rs (ξ)=∑_q^cΨ_q(r) Ψ_q^*(s)δ(ξ-Λ_q). Clearly, 𝒢^b+(t-t_0) includes all the contribution from the bound states of the full system, while 𝒢^c+(t-t_0) incorporates the contribution of the continuum states. If the system does not support bound states, the first term becomes zero and the second term controls the steady-state behavior of the system. However, the scenario changes drastically in the presence of bound states in the system. 
        
        Since the bound states are localized, the magnitude of the elements of 𝒢^b+(t-t_0) becomes small if we move away from the spatial regions where the bound states are localized. This property of 𝒢^b+(t-t_0) is very useful for our computational purpose. Suppose, we know that all the bound states are localized within the region ℛ_0. 
        Then it is legitimate to choose the elements of 𝒢^b(t-t_0) to be zero outside the region ℛ_0, and we calculate the elements 𝒢^b+(t-t_0) in the region ℛ_0 only. Although the boundary of ℛ_0 is not unique, one can extract the essential features of 𝒢^b+(t-t_0) even with a choice of region ℛ_0 which merely encompasses the spatial extent of all the overlapping bound states.
        
        In our numerics for studying persistent current oscillation, we consider a TS-N-TS device where both the boundary TS wires have the same chemical potential and temperature. Moreover, we also kept both the superconducting wires in the topologically non-trivial phase away from the topological phase transition. For such a system, there exist two Majorana bound states at the ends of each TS wires. If the TS wires are taken to be sufficiently long, the overlap of the Majorana modes from opposite ends of an individual TS wire becomes almost negligible. However, the overlap between the Majorana modes from the right end of left boundary TS wire and the left end of right boundary TS wire remains finite as long as the length of the middle N wire is not very large. Further, there may also exist bound states within the middle N wire, and it can have non-zero overlap with the Majorana modes localized at the inner edges of the boundary TS wires. In such a situation, essential behavior of the currents in the junctions can be studied with the bound state dependent part of the Green's function, i.e., 𝒢^b+(t-t_0). As it has been argued already, this 𝒢^b+(t-t_0) is chosen to be zero for the region outside ℛ_0 where ℛ_0 is extended from L_X-δ to L_X+L_Y+δ with L_Y is the length of the N wire, L_X is the length of left boundary TS wire and δ (<< L_X,  L_Y,  L_Z) is the number of sites inside the boundary wires. The value of δ solely depends on the spatial extent of bound state wave-functions inside the boundary wires. For example, if the wave-functions of the bound states around the junctions have non-zero amplitude up to the L_δ number of sites inside the boundary wires from the junction edge, then δ≈ L_δ. For our simulation, we find that it is sufficient to choose δ∽ 10 for the boundary wires with size L_X(Z)∽ 10^3 and the middle wire length L_Y ∽ 10.
        
        Given the above discussion, it may be possible to study the time evolution of the bound state dependent part of the density matrix using the following expression:
        
        ⟨ a^†_l(t) a_m(t) ⟩_b    =   ∑_r, s = 2(L_X-δ) ^2 L_X𝒢^b+_m,s (t-t_0) ⟨ a^†_r(t_0) a_s(t_0) ⟩ [𝒢^b+_l,r(t-t_0)]^†+ ∑_r,s=2L_X+1^2(L_X+L_Y)𝒢^b+_m,s(t-t_0) ⟨ a^†_r(t_0) a_s(t_0) ⟩ [𝒢^b+_l,r(t-t_0)]^†
              +∑_r, s=2(L_X+L_Y)+1^2(L_X+L_Y+ δ )𝒢^b+_m,s(t-t_0) ⟨ a^†_r(t_0) a_s(t_0) ⟩ [𝒢^b+_l,r(t-t_0)]^†,
        
        where ⟨ a^†_r(t_0) a_s(t_0) ⟩ in the first and third term denotes the initial equilibrium density matrix of the left and right boundary wire respectively, and ⟨ a^†_r(t_0) a_s(t_0) ⟩ in the second term represents the initial density matrix of the middle wire. This simplified analytic expression is quite useful in studying the characteristic features of the full density matrix, η_lm(t)=⟨ a^†_l(t) a_m(t) ⟩, provided the system hosts bound states with non-zero overlap among themselves. However, if the system does not possess bound states, or the bound states have negligible overlap, the expression (<ref>) may fail to depict the actual dynamics of η(t) as the continuum part of the full Green's function (𝒢^c+(t-t_0) ) starts to play a dominant role. If  Eq. <ref> fails, then we need to restore to the original expression in Eq. <ref> of the main text for studying the time evolution; even though it is numerically expensive.
        
        Presence of Majorana bound states at the edges of the boundary TS wires helps us to further simplify the analytical expression for ⟨ a^†_l(t) a_m(t) ⟩_b in some cases. To this end, let us write the initial equilibrium density matrix of the boundary wires as a summation of two contributions,
        ⟨ a^†_r(t_0) a_s(t_0) ⟩   =   ⟨ a^†_r(t_0) a_s(t_0) ⟩_b'+⟨ a^†_r(t_0) a_s(t_0) ⟩_c 
               =∑_b'ψ_b'^α^* (r) ψ_b'^α (s) f(λ_b'^α,μ_α, T_α)+∑_cψ_c^α^* (r) ψ_c^α(s) f(λ_c^α,μ_α, T_α),
        
        where ψ_q^α and λ^α_q are the eigenvectors and eigenvalues of the boundary wires α=X,  Z (see Eq. <ref> of the main text), and r,s=2(L_ X-δ),…, 2L_ X for α= X, and r,s=2(L_ X+L_ Y)+1, …, 2(L_ X+L_ Y+δ) for α= Z. In the second line of Eq. <ref>, the first term represents the contribution of Majorana bound states localized at the edges of boundary TS wires in topologically non-trivial phase, and the second term denotes the participation of the continuum states of the boundary wires. Here b' and c subscripts designate the bound states and the continuum states of the boundary TS wires respectively. 
        
        Substituting <ref> in <ref>, we get 
              ⟨ a^†_l(t) a_m(t) ⟩_b = ∑_r, s = 2(L_X-δ) ^2 L_X𝒢^b+_m,s (t-t_0) ⟨ a^†_r(t_0) a_s(t_0) ⟩_b'  [𝒢^b+_l,r(t-t_0)]^†+ ∑_r,s=2L_X+1^2(L_X+L_Y)𝒢^b+_m,s(t-t_0) ⟨ a^†_r(t_0) a_s(t_0) ⟩ [𝒢^b+_l,r(t-t_0)]^†
              +∑_r,s=2(L_X+L_Y)+1^2(L_X+L_Y+ δ )𝒢^b+_m,s(t-t_0) ⟨ a^†_r(t_0) a_s(t_0) ⟩_b'  [𝒢^b+_l,r(t-t_0)]^†+∑_r,s = 2(L_X-δ) ^2 L_X𝒢^b+_m,s (t-t_0) ⟨ a^†_r(t_0) a_s(t_0) ⟩_c  [𝒢^b+_l,r(t-t_0)]^†
              +∑_r,s=2(L_X+L_Y)+1^2(L_X+L_Y+ δ )𝒢^b+_m,s(t-t_0) ⟨ a^†_r(t_0) a_s(t_0) ⟩_c  [𝒢^b+_l,r(t-t_0)]^†,
        
        where the initial density matrix in the first and third term depends on the bound states of the boundary wires, that of the second term relies entirely on the initial conditions at the middle wire, whereas the initial density matrix in the fourth and final term depends on the continuum states of the boundary wires. 
        For the sake of brevity, we can further rewrite the Eq. <ref> in the following form 
        ⟨ a^†_l(t) a_m(t) ⟩_b=⟨ a^†_l(t) a_m(t) ⟩_b^B+⟨ a^†_l(t) a_m(t) ⟩_b^c, 
        
        where ⟨ a^†_l(t) a_m(t) ⟩_b^B represents first three terms at the right hand side of Eq. <ref>, and ⟨ a^†_l(t) a_m(t) ⟩_b^c represents rest of the terms at the right hand side of the same equation.
        
        
        Interestingly, we observe that the temporal evolution of the zero-bias current in a TS-N-TS device is controlled by the initial density of the middle N wire, the density of the Majorana bounds states formed at the inner edges of the boundary wires and the bound states of the full system. Therefore, the junction currents in Eq. <ref> can be approximated as 
        
        J_ XY(t)   ≈   -2γ_ XY Im[⟨ a^_2L_X+1(t)a_2L_X-1(t) ⟩_b^B],
        
        J_ YZ(t)   ≈   -2γ_ YZ Im[⟨ a^_2(L_X+L_Y)+1(t)a_2(L_X+L_Y)-1(t) ⟩_b^B].
        
        It can be noted that the above approximation remains valid as long as the boundary wires are in topologically non-trivial phase away from the phase boundary. However, near the topological phase transition point, continuum states of TS wires contribute substantially to the junction currents; thus the above approximation fails.
        
         
        
        Strong initial density dependence of junction currents in a TS-N-TS device can be explained with the help of Eq. <ref>, <ref>. For a uniform initial density n_l' of the middle N wire, the amplitudes of the oscillating modes solely depend on the wave functions of the bounds states of the full system. Thus, the initial density of the N wire can be treated as a constant factor. However, the amplitudes of the oscillating modes are also affected by the initial densities for a  non-uniform initial density at the middle wire. As the wave functions of different bound states are also spatially non-uniform, the combined effect of wave functions of bound states and inhomogeneous initial density results in suppression of the amplitudes of some selected frequency modes and (or) enhancement of others in the persistent current oscillations.
        So, the overall time-evolution dynamics of the density matrix behaves very differently than that for a uniform initial density. ]
