Distinguishing Majorana Zero Modes from Impurity States through Time-Resolved Transport

TL;DR

This paper demonstrates that time-resolved charge transport measurements can distinguish Majorana zero modes from impurity states in superconducting nanowires, revealing transient oscillations unique to Majoranas and highlighting the potential for topological quantum computing.

Contribution

The study introduces a method using time-dependent Landauer-Büttiker theory to differentiate Majorana zero modes from trivial impurity states through transient transport signatures.

Findings

Transient oscillations distinguish Majorana zero modes from impurity states.

Quasi-Majorana states can produce similar zero-bias peaks as true Majoranas.

Time-resolved transport provides a practical probe for topological quantum computing.

Abstract

We study time-resolved charge transport in a superconducting nanowire using time-dependent Landauer-B{\"u}ttiker theory. We find that the steady-state Majorana zero-bias conductance peak emerges transiently accompanied by characteristic oscillations after a bias-voltage quench. These oscillations are absent for a trivial impurity state that otherwise shows a very similar steady-state signal as the Majorana zero mode. In addition, we find that Andreev bound states or quasi-Majorana states in the topologically trivial bulk phase can give rise to a zero-bias conductance peak, also retaining the transient properties of the Majorana zero mode. Our results imply that (1) time-resolved transport may be used as a probe to distinguish between the Majorana zero mode and interfacial impurity states; and (2) the quasi-Majorana states mimic the transient signatures of the topological Majorana zero modes, indicating of the possibility of utilizing also the quasi-Majorana states in topological quantum computing as anyons with non-Abelian braiding statistics.

Figures

Equations

$$\hat{H}={\hat{H}}_e+{\hat{H}}_c+{\hat{H}}_w,$$

$${\hat{H}}_e=\sum _{k\lambda }{ϵ}_{k\lambda }{\hat{c}}_{k\lambda }^{†}{\hat{c}}_{k\lambda }$$

$${\hat{H}}_c=\sum _{jk\lambda }(T_{jk\lambda }{\hat{c}}_j^{†}{\hat{c}}_{k\lambda }+\text{h.c.}).$$

$${\hat{H}}_w$$

$$\hat{H}_{w}=\frac{1}{2}\sum_{j}{\Big{[}}\underline{\hat{\varPhi}}_{j}^{\dagger}\underline{a}_{j}\underline{\hat{\varPhi}}_{j}+{\large(}\underline{\hat{\varPhi}}_{j}^{\dagger}\underline{b}_{j}\underline{\hat{\varPhi}}_{j+1}+\text{h.c.}{\large)}{\Big{]}},$$

$${\underline{a}}_j=$$

$$\displaystyle\begin{pmatrix}J-\mu+V_{\text{Z}}&-\varDelta&0&0\\ -\varDelta&\mu-J+V_{\text{Z}}&0&0\\ 0&0&J-\mu-V_{\text{Z}}&\varDelta\\ 0&0&\varDelta&\mu-J-V_{\text{Z}}\end{pmatrix}_{j}\ ,\$$

$$\underline{b}_{j}=\begin{pmatrix}-J/2&0&-\alpha/2&0\\ 0&J/2&0&-\alpha/2\\ \alpha/2&0&-J/2&0\\ 0&\alpha/2&0&J/2\end{pmatrix}_{j}\ ,\$$

$${\hat{H}}_e=\frac{1}{2}\sum _{k\lambda }{\underline{\hat{\Phi }}}_{k\lambda }^{†}{\underline{ϵ}}_{k\lambda }{\underline{\hat{\Phi }}}_{k\lambda }$$

$${\hat{H}}_c=\frac{1}{2}\sum _{jk\lambda }({\underline{\hat{\Phi }}}_j^{†}{\underline{T}}_{jk\lambda }{\underline{\hat{\Phi }}}_{k\lambda }+\text{h.c.})$$

$${\underline{G}}_{xy}(z,z^{\prime })=-\mathrm{i}⟨{\mathcal{T}}_{\gamma }[{\underline{\hat{\Phi }}}_x(z)\otimes {\underline{\hat{\Phi }}}_y^{†}(z^{\prime })]⟩$$

$$I_{j,j+1}$$

$$f_s(j)=\{\begin{array}{ll}{\sin }^2\left[\frac{\pi sj}{2(N_w-1)}\right],& j<\frac{N_w}{s},\\ 1,& \frac{N_w}{s}\le j<\frac{(s-1)(N_w-1)}{s},\\ {\cos }^2\left[\frac{\pi sj}{2(N_w-1)}\right],& j\ge \frac{(s-1)(N_w-1)}{s},\end{array}$$

$$\underline{\mbox{\boldmath$h$}}=\begin{pmatrix}\underline{h}_{11}&0&\cdots&\underline{h}_{1C}\\ 0&\underline{h}_{22}&\cdots&\underline{h}_{2C}\\ \vdots&\vdots&\ddots&\vdots\\ \underline{h}_{C1}&\underline{h}_{C2}&\cdots&\underline{h}_{CC}\end{pmatrix}$$

$$\underline{\mbox{\boldmath$G$}}=\begin{pmatrix}\underline{G}_{11}&\underline{G}_{12}&\cdots&\underline{G}_{1C}\\ \underline{G}_{21}&\underline{G}_{22}&\cdots&\underline{G}_{2C}\\ \vdots&\vdots&\ddots&\vdots\\ \underline{G}_{C1}&\underline{G}_{C2}&\cdots&\underline{G}_{CC}\end{pmatrix}.$$

$$\mathrm{i}{\partial }_z{\underline{G}}_{xy}(z,z^{\prime })={\partial }_z\left[\theta (z,z^{\prime })⟨{\underline{\hat{\Phi }}}_x(z)\otimes {\underline{\hat{\Phi }}}_y^{†}(z^{\prime })⟩-\theta (z^{\prime },z)⟨{\underline{\hat{\Phi }}}_y(z^{\prime }{)}^{†}\otimes {\underline{\hat{\Phi }}}_x(z)⟩\right]$$

$$\mathrm{i}{\partial }_z{\underline{G}}_{xy}(z,z^{\prime })=\delta (z,z^{\prime })\left\{{\underline{\hat{\Phi }}}_x(z),{\underline{\hat{\Phi }}}_y^{†}(z^{\prime })\right\}-\mathrm{i}⟨{\mathcal{T}}_{\gamma }[\mathrm{i}{\partial }_z{\underline{\hat{\Phi }}}_x(z)]\otimes {\underline{\hat{\Phi }}}_y^{†}(z^{\prime })⟩$$

$$\displaystyle\left[{\mathrm{i}}\partial_{z}\underline{\mbox{$\mathbf{1}$}}-\underline{\mbox{\boldmath$h$}}(z)\right]\underline{\mbox{\boldmath$G$}}(z,z^{\prime})$$

$$\displaystyle\underline{\mbox{\boldmath$G$}}(z,z^{\prime})\left[-{\mathrm{i}}\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{z^{\prime}}\underline{\mbox{$\mathbf{1}$}}-\underline{\mbox{\boldmath$h$}}(z^{\prime})\right]$$

$$\mathrm{i}\frac{\mathrm{d}}{\mathrm{d}t}{\underline{G}}^{<}(t,t)-[{\underline{h}}_{CC}(t),{\underline{G}}^{<}(t,t)]=-\left[{\underline{G}}^{\text{R}}\cdot {\underline{\Sigma }}^{<}+{\underline{G}}^{<}\cdot {\underline{\Sigma }}^{\text{A}}+{\underline{G}}^{\rceil }\star {\underline{\Sigma }}^{\lceil }\right](t,t)+\text{h.c.}$$

$${\underline{\Sigma }}_{\lambda ,mn}^{\text{R}}(\omega )=\sum _k{\underline{T}}_{mk\lambda }\frac{1}{\omega -{\underline{ϵ}}_{k\lambda }+\mathrm{i}\eta }{\underline{T}}_{k\lambda n}\approx -\mathrm{i}{\underline{\Gamma }}_{\lambda ,mn}/2.$$

$$⟨{\Psi }^{\text{L}}|{\underline{h}}_{\text{eff}}=ϵ⟨{\Psi }^{\text{L}}|;{\underline{h}}_{\text{eff}}|{\Psi }^{\text{R}}⟩=ϵ|{\Psi }^{\text{R}}⟩,$$

$$⟨{\Psi }_j^{\text{L}}|\underline{\rho }(t)|{\Psi }_k^{\text{L}}⟩=\sum _{\lambda }\left\{{\Gamma }_{\lambda ,jk}{\Lambda }_{\lambda ,jk}+V_{\lambda }{\Gamma }_{\lambda ,jk}\left[{\Pi }_{\lambda ,jk}(t)+{\Pi }_{\lambda ,kj}^{\ast }(t)\right]+V_{\lambda }^2{\Gamma }_{\lambda ,jk}{\mathrm{e}}^{-\mathrm{i}({ϵ}_j-{ϵ}_k^{\ast })t}{\Omega }_{\lambda ,jk}\right\},$$

$${\Gamma }_{\lambda ,jk}$$

$${\Lambda }_{\lambda ,jk}$$

$${\Pi }_{\lambda ,jk}(t)$$

$${\Omega }_{\lambda ,jk}$$

$$⟨{\phi }_m|\underline{\rho }(t)|{\phi }_n⟩=\sum _{jk}\frac{⟨{\phi }_m|{\Psi }_j^{\text{R}}⟩}{⟨{\Psi }_j^{\text{L}}|{\Psi }_j^{\text{R}}⟩}\frac{⟨{\Psi }_k^{\text{R}}|{\phi }_n⟩}{⟨{\Psi }_k^{\text{R}}|{\Psi }_k^{\text{L}}⟩}⟨{\Psi }_j^{\text{L}}|\underline{\rho }(t)|{\Psi }_k^{\text{L}}⟩,$$

$${\Lambda }_{\lambda ,jk}=\beta {\int }_{-\infty }^{\infty }\frac{\mathrm{d}z}{2\pi }\frac{1}{(z-z_1)(z-z_2)({\mathrm{e}}^z+1)}$$

$${\Lambda }_{\lambda ,jk}=\mathrm{i}\beta \left[\frac{1}{z_2-z_1}\frac{1}{{\mathrm{e}}^{z_2}+1}-\sum _{n=0}^{\infty }\frac{1}{(w_n-z_1)(w_n-z_2)}\right].$$

$$\sum _{n=0}^{\infty }\frac{1}{\left[\mathrm{i}(\pi +2\pi n)-z_1\right]\left[\mathrm{i}(\pi +2\pi n)-z_2\right]}$$

$${\Lambda }_{\lambda ,jk}$$

$${\Pi }_{\lambda ,jk}(t)={\beta }^2{\int }_{-\infty }^{\infty }\frac{\mathrm{d}z}{2\pi }\frac{{\mathrm{e}}^{\frac{\mathrm{i}}{\beta }(z-z_2)t}}{(z-z_1)(z-z_2)(z-z_3)({\mathrm{e}}^z+1)},$$

Full text

I Introduction

Topological quantum computing Nayak et al. (2008) is an active field of research based on the key idea to reduce quantum decoherence issues by using topologically protected states Hasan and Kane (2010); Qi and Zhang (2011). Majorana fermions are their own antiparticles Majorana (1937), and their condensed-matter analogues, Majorana bound state or Majorana zero mode (MZM), retain this feature Kitaev (2001). They are thus considered to be promising candidates for technological advances in topological quantum computing Kim et al. (2018); Lutchyn et al. (2018) since their non-abelian statistics allow performing quantum computation protected from environmental perturbations Ivanov (2001). Even though various experimental signatures of MZM have been reported Mourik et al. (2012); Finck et al. (2013); Nadj-Perge et al. (2014); Albrecht et al. (2016); Suominen et al. (2017); Gül et al. (2018); Wang et al. (2018); Liu et al. (2018), a clear and unambiguous detection and the consequent control of these states has proven difficult so far. For example, other types of bound states Ruby et al. (2015) or interfacial impurity states (IS) also give rise to in-gap states that contribute to transport or scanning tunneling spectroscopy signals. Therefore probes that unambiguously distinguish between MZM and IS are highly desirable.

Time-resolved spectroscopies allow for studying the dynamics of various processes such as charge transport Ochoa et al. (2015). For instance, in a transport setup exhibiting the MZM, there is no guarantee of instantly relaxing to a steady-state configuration once the junction has been “switched on” by, e.g., applying an external perturbation. In contrast, the nonequilibrium problems are often much richer and more interesting than equilibrium properties Kundu and Seradjeh (2013); Covito et al. (2018); Bondyopadhaya and Roy (2018). This is especially relevant when nowaday transport measurements are pushing the temporal resolution to sub-picosecond regime Prechtel et al. (2012); Cocker et al. (2013); Hunter et al. (2015); Rashidi et al. (2016); Cocker et al. (2016); Marguerite et al. (2017); Jelic et al. (2017); McIver et al. (2018), and these ultrafast processes can be observed in real time.

In this paper we propose time-resolved transport as a probe in order to reveal the difference between topological MZM and ordinary IS. We simulate the transient dynamics in a quantum wire coupled to metallic electrodes using the time-dependent Landauer–Büttiker formalism Perfetto et al. (2008); Stefanucci and van Leeuwen (2013a); Tuovinen et al. (2014a); Ridley et al. (2015); Tuovinen et al. (2016a); Ridley and Tuovinen (2017) extended to include superconducting states in a Nambu spinor representation. By comparing the time-dependent build-up of a steady-state current after a sudden quench of the bias voltage between (i) a topological state with MZM and (ii) a non-topological state with trivial impurity end states, we discover that the dynamics for (i) and (ii) look significantly different. For case (i) the time-resolved current shows pronounced oscillations that shift with the applied bias voltage and correspond to transitions between the biased electrodes and the MZM. By contrast, for case (ii) no such oscillations are observed. In addition, we study the transient response of quasi-Majorana states (QMS) in the topologically trivial phase Kells et al. (2012); Liu et al. (2017); Vuik et al. (2018); Moore et al. (2018); Yavilberg et al. (2019), and we find that the QMS mimic the signatures of the MZM both in the stationary and transient regimes. The resulting Fourier spectra of the time-resolved current can therefore be used to identify the MZM or QMS.

II Model and method

We consider a normal metal–superconductor–normal metal (NSN) junction, see Fig. 1. The superconducting central region of the junction is a nanowire in proximity to an $s$-wave bulk SC with order parameter $\varDelta$. The nanowire in addition features a strong spin-orbit interaction (e.g., InSb Plissard et al. (2012); van Weperen et al. (2015) or InAs Fasth et al. (2007); Albrecht et al. (2016)) which favors aligning the spins along the $\pm y$ direction. An external magnetic field parallel to the nanowire breaks time-reversal symmetry and aligns the spins along the $\pm z$ direction, introducing a Zeeman splitting $V_{\text{Z}}=g\mu_{\text{B}}B/2$ where $g$ is the Landé factor and $\mu_{\text{B}}$ the Bohr magneton. A suitable combination of these effects has been shown to host a MZM in the nanowire, exponentially localized at the edges Kitaev (2001); Lutchyn et al. (2010); Oreg et al. (2010); Alicea (2010); Stoudenmire et al. (2011); Mourik et al. (2012). Specifically the infinite nanowire is in a topologically nontrivial phase for $V_{\text{Z}}>\varDelta>0$  Lutchyn et al. (2010); Oreg et al. (2010); Stoudenmire et al. (2011), from which MZM emerge in the case of a finite wire. For the present study, the specific structure of the electrodes, other than being a normal metal with relatively broad bandwidth (e.g., Au, Ag or Cu), is unimportant as we concentrate on the effects within the nanowire.

We write the total Hamiltonian as

[TABLE]

where the individual components for the electrodes and coupling are characterized by the single-particle energy dispersion in the electrodes $\epsilon_{k\lambda}$ and by the coupling matrix elements $T_{jk\lambda}$ between the states in the nanowire and the electrodes Tuovinen et al. (2014a):

[TABLE]

and

[TABLE]

Here $k\lambda$ labels the $k$-th basis element in the $\lambda$-th electrode, and $j$ labels the atomic sites on the nanowire. The nanowire, in turn, is characterized by Stoudenmire et al. (2011); Wu and Cao (2012)

[TABLE]

where $J,\mu,\alpha,V_{\text{Z}}$, and $\varDelta$ are parameters for hopping, chemical potential, spin-orbit coupling, Zeeman splitting, and pairing potential, respectively. The operators $\hat{c}_{xs}^{(\dagger)}$ annihilate (create) electrons with spin $s\in\{\uparrow,\downarrow\}$ in a region specified by $x$. The spin indices are summed when suppressed and $\sigma_{2,3}$ are Pauli matrices. For indices $x,y$ belonging either to the electrodes or to the nanowire, the creation and annihilation operators satisfy the anticommutation relations $\{\hat{c}_{xs},\hat{c}_{ys^{\prime}}^{\dagger}\}=\delta_{xy}\delta_{ss^{\prime}}$.

At times $t>0$ the electrode energy levels are suddenly shifted, corresponding to a quench of the bias voltage, $\epsilon_{k\lambda}\to\epsilon_{k\lambda}+eV_{\lambda}$. For a two-terminal device ($\lambda\in\{\text{S},\text{D}\}$, see Fig. 1) this out-of-equilibrium condition is defined by the source-drain voltage $V_{\text{SD}}=V_{\text{S}}-V_{\text{D}}$. The transport setup is considered partition-free Cini (1980); Stefanucci and Almbladh (2004); Ridley and Tuovinen (2018) meaning that the whole system is initially contacted in a global thermo-chemical equilibrium at unique chemical potential $\mu$ and at inverse temperature $\beta\equiv(k_{\text{B}}T)^{-1}$.

For a compact notation we introduce Nambu spinors Nambu (1960); Zeng et al. (2003); Xing et al. (2007) $\underline{\hat{\varPhi}}_{x}\equiv(\hat{\varPhi}_{x}^{1},\hat{\varPhi}_{x}^{2},\hat{\varPhi}_{x}^{3},\hat{\varPhi}_{x}^{4})^{T}\equiv(\hat{c}_{x\uparrow},\hat{c}_{x\downarrow}^{\dagger},\hat{c}_{x\downarrow},\hat{c}_{x\uparrow}^{\dagger})^{T}$, and the anticommutation relation is then understood componentwise $\{\hat{\varPhi}_{x}^{\mu},(\hat{\varPhi}_{y}^{\nu})^{\dagger}\}=\delta_{xy}\delta^{\mu\nu}$. Here we denote quantities in the Nambu$\otimes$spin space by an underline. This representation allows for writing the Hamiltonian for the nanowire in a Bogoliubov–de Gennes form Bogoliubov (1958); de Gennes (1966)

[TABLE]

where we introduced on-site and nearest-neighbor contributions Wu and Cao (2012)

[TABLE]

[TABLE]

respectively. The electrode and coupling parts of the Hamiltonian are then also expanded in the Nambu$\otimes$spin basis although they do not involve the SC pairing potential:

[TABLE]

with $\underline{\epsilon}_{k\lambda}=\epsilon_{k\lambda}\text{diag}(1,-1,1,-1)$ and

[TABLE]

with $\underline{T}_{jk\lambda}=T_{jk\lambda}\text{diag}(1,-1,1,-1)$.

By using the nonequilibrium Green’s function approach Stefanucci and van Leeuwen (2013a) we conveniently access both transient and steady-state responses in the setup above. The one-electron Green’s function is defined as a contour-ordered tensor product of the spinor field operators Zeng et al. (2003)

[TABLE]

where the contour-ordering operator $\mathcal{T}_{\gamma}$ is taken for the variables $z,z^{\prime}$ on the Keldysh contour $\gamma$ Stefanucci and van Leeuwen (2013a). The form in Eq. (10) automatically handles both normal and anomalous components of the Green’s function Stefanucci et al. (2010a). In the Supplemental Material sm we show that the equations of motion for the Green’s function are exactly the same as those in Refs. Tuovinen et al. (2014a, 2016a), and hence we derive in a similar fashion a closed expression for the time-dependent one-particle reduced density-matrix (TD1RDM) within the nanowire, $\underline{\rho}(t)\equiv-{\mathrm{i}}\underline{G}^{<}(t,t)$ from the lesser Green’s function. In order to obtain a closed solution to the equation of motion we have described the electrodes within wide-band approximation, where the electronic levels of the nanowire are in a narrow range compared to the bandwidth of the electrodes. The coupling strength between the nanowire and the electrodes is characterized by the frequency-independent tunneling rate $\varGamma_{\lambda}$.

As the TD1RDM gives us full information on the local charge and current densities within the nanowire, we calculate the total current through the nanowire by considering a bond current between two atomic sites. In addition, the traditional bond-current operator has to be adapted to include the contribution from the spin-orbit coupling and from the SC pair potential Solomon et al. (2010); Stefanucci et al. (2010b); Winkler et al. (2017). In the Supplemental Material sm we derive the following expression for the bond current between the sites $j$ and $j+1$ within the nanowire:

[TABLE]

where $\langle\cdot\rangle$ denotes elements of the TD1RDM.

III Results

III.1 Emergence of the Majorana zero mode

Using Eq. (11) we calculate the (steady-state) current-voltage characteristics for nanowires of varying lengths. For the nanowire we choose the parameters $J=1$, $\alpha=0.5$, $V_{\text{Z}}=0.25$, $\varDelta=0.1$, and $\mu=0$ Wu and Cao (2012). This fixes the units to the hopping energy; if the values of this quantity are in the eV regime, then times are measured in the units of inverse hoppings which is on the order of femtoseconds. The coupling strength from the terminal sites of the nanowire to the electrodes is chosen such that the tunneling rate $\varGamma_{\lambda}=0.01$. The bias voltage is applied symmetrically for the source and drain electrodes $V_{S}=-V_{D}\equiv V$, and we consider the zero-temperature limit.

In Fig. 2 we show the differential conductance against the applied bias voltage (around a low voltage window). We observe clearly how the MZM behaves as a “half a fermion” on both terminals of the nanowire leading to two peaks of half the conductance quantum. When the coupling between the MZM becomes weaker by elongating the nanowire, the two peaks merge into one at exactly zero bias voltage resulting in one conductance quantum. The inset of Fig. 2 shows the exponential localization of the MZM for $N_{w}=60$.

III.2 Transient signature of the Majorana zero mode

We evaluate transient currents through a $N_{w}=50$ nanowire by considering the two centermost sites in Eq. (11). We single out the MZM by applying a small bias window so that the oscillations in the time-resolved signal are only due to virtual transitions from the biased Fermi level of the electrode to the zero-energy mode in the nanowire. In Fig. 3(a) we show the differential conductance for a $N_{w}=50$ nanowire for three different cases: (1) ordinary SC wire (same as Fig. 2 but for $V_{\text{Z}}=0$), (2) topological SC wire corresponding to Fig. 2, and (3) an ordinary SC wire with an impurity state localized at its edges. We model the impurity states by modified tight-binding parameters Robinson et al. (2008); Wehling et al. (2010) for the terminal sites in the nanowire, $j=\{1,N_{w}\}$ in Eqs. (II) and (7). More specifically, we use (modified parameters denoted by a tilde) $\widetilde{\mu}=\widetilde{J}=0.1J$, $\widetilde{\alpha}=0.1\alpha$, and $\widetilde{\varDelta}=\widetilde{V}_{\text{Z}}=0$. For our purposes the exact formulation is not too important as long as there is a separate state within the gap with different topological character compared to the MZM.

Importantly, while the steady-state $dI/dV$ signals of cases (2) and (3) look qualitatively similar (Fig. 3(a)), the transient signals in Fig. 3(b) for the three cases is qualitatively different. For the pristine wire without the MZM the current signal is zero on average due to there being no transport channels within the SC gap and the small bias window. When the IS is present as an in-gap state, the transient current rises rapidly but also saturates relatively fast to its stationary value within couple of hundred units of inverse hopping. (For hopping energies in the eV scale we have $J^{-1}\sim 0.658$ fs.) The IS is directly connected to the electrodes resulting in a strong hybridization and in a relatively fast decay of the transient. In contrast, the MZM at the edges of the nanowire have a different topological character being weakly coupled to each other although they are far apart, and even though the MZM is also directly connected to the electrodes, the hybridization of the MZM is weaker resulting in transient oscillations for thousands of time units, i.e., up to picoseconds. The decay rate can be approximated by the expectation value of the tunneling rate operator: $\gamma=\sum_{j=1}^{2}\langle\varphi_{j}|\varGamma|\varphi_{j}\rangle$, where $\varGamma\equiv\sum_{\lambda}\varGamma_{\lambda}$ and $|\varphi\rangle$ are the IS or MZM eigenvectors, see the dashed lines in Fig. 3(b). For identical wire-electrode coupling, the decay time $1/\gamma$ of the MZM transient current is roughly $5$ times the one of the IS.

Crucially, the MZM additonally shows transient current oscillations, unlike the IS. This striking difference between the MZM and IS cases is clearly seen by taking the Fourier transforms of the time-dependent signals, see Fig. 3(c). The low-frequency regime shows pronounced peaks for the MZM case, and the frequency of the first peak exactly corresponds to the difference between the biased Fermi level of the electrode and the MZM (indicated by (i) in the figure). The analogous peaks in the case of the IS are strongly diminished. Importantly, we have also checked that by artificially increasing the decay time for the IS case by decreasing the wire-electrode coupling, there are still no pronounced transient oscillations in the IS case. Before entering the band of all possible transitions outside the SC gap (${\omega}\geq 2\varDelta=0.2$, indicated by (iii) in the figure) we observe additional transitions between the MZM and states close to the gap edge (indicated by (ii) in the figure). These resonances remain independent of the applied voltage confirming that they result from intra-level transitions within the nanowire. Overall the transient features of the MZM are distinctly different from the IS.

III.3 Comparison between Majorana zero modes and quasi-Majorana states

Even though we found a clear distinction between trivial IS and the topological MZM, one may still wonder whether and how other in-gap states deeply in the topologically trivial regime for the same model of the nanowire would contribute to the time-resolved signal. Recently, it has been studied that in the parameter regime $\varDelta\ll V_{\text{Z}}<\mu$ the resulting in-gap QMS emerge without any additional impurities but by adding a smooth confining potential Kells et al. (2012); Liu et al. (2017); Vuik et al. (2018); Moore et al. (2018). These states can also be tuned arbitrarily close to zero energy, thereby mimicking the behaviour of the MZM.

We implement a confining potential within the nanowire as a simple function of the lattice coordinate $j\in[0,N_{w})$ labeling the atomic sites on the nanowire:

[TABLE]

where $s$ controls the smoothness at the edges. We then re-cast the values of the spin-orbit interaction and pair potential accordingly: $\alpha\to\alpha f_{s}(j)$ and $\varDelta\to\varDelta f_{s}(j)$. Large values of $s$ correspond to an abrupt hard-wall confinement where both $\alpha$ and $\varDelta$ remain constant (nonzero) throughout the nanowire (cf. previous subsections). For smaller values of $s$ the spin-orbit coupling and the induced superconductivity go to zero smoothly at the edges of the nanowire.

To study the topologically trivial parameter regime we focus our discussion on three additional cases: For an $N_{w}=50$ nanowire we set $\alpha=0.5$, $V_{Z}=1.2$, $\Delta=0.1$, and $\mu=2.0$ both with abrupt ($s\to\infty$) and smooth ($s=\{17,9\}$) confinement potentials according to Eq. (12). In Fig. 4(a) we show how the QMS is brought to a MZM-like state (peak at zero bias) by making the confining potential smoother. We have checked that other shapes for the potential profile do not modify the results qualitatively. The transient signature, see Figs. 4(b) and 4(c), of these states is also similar to the MZM: (1) The current oscillates with a dominant frequency corresponding to the lead-nanowire transition, and (2) the lifetime of the oscillations is similar or even longer compared to the MZM case. However, unless the QMS appears exactly at zero energy, the transient oscillations are suppressed, and the Fourier peak corresponding to the ‘smooth enough’ case is considerably more pronounced.

IV Conclusion

We studied the time-dependent features of Majorana zero modes and quasi-Majorana states in a superconducting nanowire in contrast with trivial impurity bound states. The transient features related to MZM and QMS were found to be completely different than the ones resulting from a simple impurity model: The MZM and QMS transients were found to decay very slowly with a pronounced oscillation frequency due to a weaker hybridization of the MZM and QMS with the electrode states compared to the IS. This finding could be utilized in possible detection and identification of the MZM or QMS via ultrafast transport measurements Prechtel et al. (2012); Cocker et al. (2013); Hunter et al. (2015); Rashidi et al. (2016); Cocker et al. (2016); Marguerite et al. (2017); Jelic et al. (2017); McIver et al. (2018); Sato et al. (2019).

We also found that even though the QMS are only protected by the smoothness of the confining potential (in contrast to the topological protection of the MZM), the QMS may still mimic the transient signature of the MZM. This effect could also be utilized by employing braiding schemes for the QMS in topological quantum computation Vuik et al. (2018). Since topological properties of the MZM should be robust against electronic interactions Stoudenmire et al. (2011), it would be a promising direction for future work to understand this effect for the QMS and how it might be manipulated and controlled.

In practice the sudden switch of the bias voltage employed here could be replaced by a short light pulse in the THz regime to excite the system away from its thermal equilibrium. In the case of an ultrashort laser excitation the current response of MZM or QMS could initially be suppressed and then recover transiently with the oscillations as a characteristic signature, similarly to the amplitude mode oscillations of laser-driven ordered phases Kemper et al. (2015); Sentef et al. (2016); Kemper et al. (2017). Together with ultrafast optical switching of chiral superconductors Dehghani and Mitra (2017); Claassen et al. (2018) or nonequilibrium engineering of topologically nontrivial states of matter Oka and Aoki (2009); Kitagawa et al. (2011); Lindner et al. (2011); Jiang et al. (2011); Sentef et al. (2015); Claassen et al. (2016); Hübener et al. (2017); Thakurathi et al. (2017); Topp et al. (2018) our findings highlight the great potential of ultrafast techniques for advances towards topological quantum computation.

Acknowledgements.

R.T. and M.A.S. acknowledge funding by the DFG (Grant No. SE 2558/2-1) through the Emmy Noether program. R.v.L. acknowledges funding from the Academy of Finland (Project No. 317139). E.P. and G.S. acknowledge EC funding through the RISE Co-ExAN (Grant No. GA644076). E.P. also acknowledges funding from the European Union project MaX Materials design at the eXascale H2020-EINFRA-2015-1, Grant Agreement No. 676598 and Nanoscience Foundries and Fine Analysis-Europe H2020-INFRAIA-2014-2015, Grant Agreement No. 654360. G.S. also acknowledges Tor Vergata University for financial support through the Mission Sustainability Project 2DUTOPI.

I Transport setup and partitioning the Green’s function

Even though in the main text we considered a two-terminal device, the description readily allows for a more general treatment, and we now label by $\lambda$ an arbitrary number of electrodes. The central region $C$, for which we had the superconducting nanowire in the main text, can also take a more arbitrary shape. We only assume there to be no direct connection between any of the electrodes but the coupling is always through the central region. Then, the Hamiltonian for the full transport setup may be partitioned accordingly

[TABLE]

with $(\underline{h}_{\lambda\lambda^{\prime}})_{kk^{\prime}}=\underline{\epsilon}_{k\lambda}\delta_{\lambda\lambda^{\prime}}\delta_{kk^{\prime}}$ for the electrodes, and $(\underline{h}_{C\lambda})_{mk\lambda}=\underline{T}_{mk\lambda}$ for the couplings. For the central region, $\underline{h}_{CC}$, we may use the “on-site” and “nearest-neighbor” contributions [Eqs. (II) and (7)], or consider some other arbitrary structure. We further denote the matrices for the full transport setup as boldface symbols. It is important to notice how the electrode blocks, $\underline{h}_{\lambda\lambda}=\underline{h}_{\lambda\lambda}(z)$, are different for the vertical and horizontal branches of the Keldysh contour due to the shift in energy levels at $t>0$. Also, we stress here that the block structure in Eq. (S1) does not refer to the Nambu$\otimes$spin space but it is of dimension $(N_{e}+1)\times(N_{e}+1)$ where $N_{e}$ is the number of electrodes. Each block then accounts for the individual dimension of the corresponding partition. The matrix elements in the Green’s function in Eq. (10) (indices $x,y$ belonging either to the electrodes or to the central region) therefore label the transport setup in the same block form

[TABLE]

We may derive the equation of motion for the Green’s function by

[TABLE]

where the step function is defined on the Keldysh contour $\gamma$ according to the contour-ordering operator $\mathcal{T}_{\gamma}$ Stefanucci and van Leeuwen (2013b). Evaluating the derivative gives

[TABLE]

where the anticommutator gives simply \delta_{xy}\underline{\mbox{1}} and the evolution of the spinor operator can further be derived from its equation of motion. Depending on which region the index $x$ belongs to (and the corresponding structure of the Hamiltonian in that region), the time-evolution of the field operator is completely specified. The equations of motion for the whole transport setup then take the matrix form Stefanucci and van Leeuwen (2013b); Tuovinen et al. (2013, 2014b, 2016b)

[TABLE]

which the Green’s function satisfies being antiperiodic along the contour (Kubo–Martin–Schwinger boundary condition Kubo (1957); Martin and Schwinger (1959)).

We see that the equations of motion are the same as those of Ref. Tuovinen et al. (2013, 2014b), hence we may in similar fashion, using the Langreth rules Langreth (1976); Stefanucci and van Leeuwen (2013b), derive an equation for the equal-time lesser Green’s function with indices on the central region $\underline{G}_{CC}^{<}$. This is a key quantity as it relates to the time-dependent one-particle reduced density-matrix (TD1RDM) by $\underline{\rho}_{CC}(t)=-{\mathrm{i}}\underline{G}_{CC}^{<}(t,t)$. From now on we will only discuss quantities in the subspace of the central region, so we will drop the subscript ‘$CC$’. The lesser Green’s function at the equal-time limit is given by Tuovinen et al. (2013)

[TABLE]

where the time-convolutions on the horizontal and vertical branches of the Keldysh contour are defined as $[f\cdot g](t,t)=\int_{0}^{\infty}\mathrm{d}\bar{t}f(t,\bar{t})g(\bar{t},t)$ and $[f\star g](t,t)=-{\mathrm{i}}\int_{0}^{\beta}\mathrm{d}\tau f(t,\tau)g(\tau,t)$. The superscripts R,A,$<,\rceil,\lceil$ refer to the retarded, advanced, lesser, right and left Keldysh components, respectively Stefanucci and van Leeuwen (2013b); Tuovinen et al. (2013). The embedding self-energy, $\varSigma$, accounts for the coupling between the central region and the electrodes Tuovinen et al. (2014b).

We note that the left-hand side of Eq. (S7) corresponds to a Liouville-type of equation for the density matrix of an isolated central region whereas the right-hand side gives rise to an open transport setup as in connection to the electrode environment. The time-convolutions on the right-hand side can further be identified as source and drain terms, and the ones including the imaginary track of the Keldysh contour to include the initial contacting of the separate regions. Importantly, within the so-called wide-band approximation (WBA) for the embedding self-energy, Eq. (S7) becomes a closed equation for the equal-time lesser Green’s function and the TD1RDM can be solved analytically.

II Solution to the equation of motion

In order to close the equation of motion we now describe the electrodes in the framework of wide-band approximation (WBA), where the electronic levels of the central region are in a narrow range compared to the electrode bandwidth. The validity of WBA has been discussed in, e.g., Refs. Zhu et al. (2005); Maciejko et al. (2006); Verzijl et al. (2013); Covito et al. (2018), and for the purpose of the present work (weak coupling of the central region to electrodes of large bandwidth), this is a well-justfied approximation. In frequency space the retarded Keldysh component of the embedding self-energy can then be written as

[TABLE]

The advanced component is given simply by conjugating this. The other components of the self-energy ($<,\lceil$) may further be derived from the retarded and advanced components Stefanucci and van Leeuwen (2013b); Tuovinen et al. (2013). The time-domain quantities in Eq. (S7) are then obtained by Fourier transforming. Looking at Eq. (S7) and the earlier work in Refs. Tuovinen et al. (2013, 2014b) we may use the fact that the same equations have the same solutions, i.e., including the Nambu$\otimes$spin structure in the Hamiltonian of the central region (e.g., spin-orbit coupling, Zeeman splitting and pairing field) adds no extra complication to the evolution of the Green’s function. The only difference is in the Nambu$\otimes$spin structure of the matrices.

It is useful to introduce a nonhermitian effective Hamiltonian $\underline{h}_{\text{eff}}=\underline{h}_{CC}-{\mathrm{i}}\underline{\varGamma}/2$ for which the left and right eigenvalue equations are

[TABLE]

where the eigenvectors and eigenvalues correspond to the $4\times 4$ Nambu$\otimes$spin space. The solution for the TD1RDM expanded in the left eigenbasis takes the explicit form Tuovinen et al. (2014b)

[TABLE]

where

[TABLE]

Here $f({\omega}-\mu)=({\mathrm{e}}^{\beta({\omega}-\mu)}+1)^{-1}$ is the Fermi function at inverse temperature $\beta$ and chemical potential $\mu$. Evaluating the TD1RDM in a physically relevant basis, e.g., the localized site basis of the central region $\{|\varphi\rangle\}$, is then readily done as a basis transformation from the left eigenbasis to the desired one

[TABLE]

which follows from the biorthogonality of the left and right eigenvectors. The TD1RDM is then simply given by evaluating the terms in Eqs. (S12), (S13) and (LABEL:eq:omega) for all indices $j,k$ and time parameter $t$, and then inserting into Eqs. (S10) and (S15).

The integrands in Eqs. (S12), (S13) and (LABEL:eq:omega) have a fairly simple analytic structure: The “$1/({\omega}-z)$” type of terms have simple poles at ${\omega}=z$ whereas the Fermi function has simple poles at the Matsubara frequencies given by ${\omega}_{n}=(2n+1)\pi/(-{\mathrm{i}}\beta)$. Expressions similar to those in Eqs. (S12), (S13), (LABEL:eq:omega) have been found, e.g., in Refs. Jauho et al. (1994); Galperin and Nitzan (2003) and integrated correspondingly using contour integration techniques. In Ref. Tuovinen et al. (2014b) the frequency integrals in Eqs. (S12), (S13), (LABEL:eq:omega) were evaluated analytically in the zero-temperature limit to obtain a result for the TD1RDM in terms of logarithms and exponential integral functions. Here we evaluate these integrals analytically at arbitrary (inverse) temperature in the Fermi functions, and we will detail these steps next.

III Details of the Fermi integrals

Making a change of variables $z=\beta({\omega}-\mu)$ in Eq. (S12) gives

[TABLE]

where we defined $z_{1}=\beta({\epsilon}_{j}-\mu_{\lambda})$ and $z_{2}=\beta({\epsilon}_{k}^{*}-\mu_{\lambda})$ with $\mu_{\lambda}=\mu+V_{\lambda}$. This integrand has simple poles at $z=z_{1}$, $z=z_{2}$ and $z=w_{n}={\mathrm{i}}(\pi+2\pi n)$, see Fig. S1. The spectrum of the complex eigenvalues of the nonhermitian matrix $\underline{h}_{\text{eff}}$ is such that the eigenvalues, ${\epsilon}_{j}$, lie in the lower-half plane (LHP) whereas the complex conjugated ones, ${\epsilon}_{k}^{*}$, lie in the upper-half plane (UHP).

For the “$(z-z_{n})^{-1}$” contributions the residues are simply one and for the Fermi function we have $\text{Res}\left[({\mathrm{e}}^{z}+1)^{-1}\ ,\ z=w_{n}\right]=-1$. Then, we can close the integral in Eq. (S16) in the UHP as shown in Fig. S1, and using the residue theorem we get

[TABLE]

The infinite sum can be written as

[TABLE]

where we defined $a=({\mathrm{i}}z_{1}+\pi)/2\pi$, $b=({\mathrm{i}}z_{2}+\pi)/2\pi$, and $\psi$ is the digamma function which is defined as the logarithmic derivative of the gamma function, $\psi(z)=\frac{\mathrm{d}}{\mathrm{d}z}\log\Gamma(z)$ Weisstein (a). We can then insert the result of the sum back into Eq. (S17) and couple the terms by simplifying

[TABLE]

where we also inserted back the definitions of $z$’s. It is important to notice that we did not do anything but manipulations after using the residue theorem; the infinite sum was rewritten in terms of a special function $\psi$ which is broadly known in computational sciences and readily implemented for example in the GNU Scientific Library Galassi et al. . Eq. (S19) is our final result for ${\varLambda}_{\lambda,jk}$ for arbitrary values of $\beta$. We note in passing that it would give completely equivalent result if the integral was closed in the LHP.

Making the same change of variables in Eq. (S13) as in the previous case leads to

[TABLE]

where we defined $z_{1}=\beta({\epsilon}_{j}-\mu)$, $z_{2}=\beta({\epsilon}_{j}-\mu_{\lambda})$ and $z_{3}=\beta({\epsilon}_{k}^{*}-\mu_{\lambda})$. Also in this case we notice poles in the complex plane, similarly as in Fig. S1. In this case, however, we may close the integral only in the UHP due to the exponential in the numerator, and we get according to the residue theorem

[TABLE]

We may manipulate the infinite sum in Eq. (S21) as

[TABLE]

where we defined $a=({\mathrm{i}}z_{1}+\pi)/2\pi$, $b=({\mathrm{i}}z_{2}+\pi)/2\pi$, $c=({\mathrm{i}}z_{3}+\pi)/2\pi$ and $x=-2\pi t/\beta$. In this case the infinite sum will give another type of special function, the hypergeometric function ${}_{2}{F}_{1}$ Weisstein (b):

[TABLE]

The hypergeometric function together with the Pochhammer symbol are defined as Weisstein (b, c)

[TABLE]

Inserting the definitions for $a$, $b$, $c$ and $x$ (and also the previously introduced variables $z$) leads to

[TABLE]

where we defined an auxiliary function

[TABLE]

This calculation was only for the infinite sum in Eq. (S21). Inserting the definitions of $z$’s into the first term gives

[TABLE]

Combining the terms finally gives

[TABLE]

for arbitrary values of $\beta$. Similarly here, after using the residue theorem, we only manipulated the expressions so that we could identify a known function ${}_{2}{F}_{1}$. Conveniently, the hypergeometric function is also widely used in computational sciences, and both fast and accurate implementations of it are available Michel and Stoitsov (2008).

In the third case, in Eq. (LABEL:eq:omega), we do the same change of variables as before to get

[TABLE]

where we defined $z_{1}=\beta({\epsilon}_{j}-\mu)$, $z_{2}=\beta({\epsilon}_{j}-\mu_{\lambda})$, $z_{3}=\beta({\epsilon}_{k}^{*}-\mu_{\lambda})$ and $z_{4}=\beta({\epsilon}_{k}^{*}-\mu)$. The pole structure is again similar to the one shown in Fig. S1, and we may close also this integral in the UHP. Again, according to the residue theorem we get as a result

[TABLE]

The infinite sum may again be manipulated as

[TABLE]

where we defined $a=({\mathrm{i}}z_{1}+\pi)/2\pi$, $b=({\mathrm{i}}z_{2}+\pi)/2\pi$, $c=({\mathrm{i}}z_{3}+\pi)/2\pi$ and $d=({\mathrm{i}}z_{4}+\pi)/2\pi$. Also this sum has an expression in terms of the digamma function

[TABLE]

Inserting the expressions for $a$, $b$, $c$ and $d$, and then further the expressions for $z_{1}$, $z_{2}$, $z_{3}$ and $z_{4}$ leads to

[TABLE]

Combining the terms in Eq. (S30) gives as the final result

[TABLE]

for arbitrary values of $\beta$.

Finally, inserting Eqs. (S19), (S28) and (III) into Eq. (S10) gives then the TD1RDM at arbitrary temperature. When the asymptotic behaviour of the digamma and hypergeometric function is studied, the results in Eqs. (S19), (S28) and (III) can be shown to reduce to those in Ref. Tuovinen et al. (2014b) at the zero-temperature limit ($\beta\to\infty$) Tuovinen (2016). We also note that congruent results involving equivalent special functions have recently been reported in Refs. Ridley et al. (2016a, b, 2017).

IV Inclusion of sudden electromagnetic fields in the central region

It is also possible to include a sudden switch-on of an electromagnetic field in the Hamiltonian of the central region. For example, this includes the possibility for a static potential profile (e.g. a gate voltage) $u_{mn}$, between basis states $m,n$ of the central region, to be added to the “on-site” contribution $\underline{a}$ [Eq. (II)]. Also, for the “nearest-neighbor” contribution $\underline{b}$ [Eq. (7)], it is possible to consider a Peierls phase $\gamma_{mn}=-\gamma_{nm}$ accounting for a magnetic field (normalized to the flux quantum $\phi_{0}=h/2e$) when traversed along a closed loop of states $m,n$. For a general description, we simply consider a perturbed Hamiltonian $\widetilde{\underline{h}}_{CC}$ out of equilibrium (signified by a tilde), and use the unperturbed Hamiltonian $\underline{h}_{CC}$ in equilibrium. Then, a formula for the TD1RDM similar to Eq. (S10) can be derived as Tuovinen et al. (2014b)

[TABLE]

where the introduced terms $\widetilde{\varGamma}$, $\widetilde{\varPi}$ and $\widetilde{\varOmega}$ take a slightly more intricate form compared to those in Eq. (S10) as the eigenbases of the unperturbed and perturbed Hamiltonians, in general, do not need to be the same. Therefore, we need to take the corresponding overlaps into account

[TABLE]

where the tildes signify that the corresponding quantities are calculated from the perturbed Hamiltonian $\widetilde{\underline{h}}_{CC}$, and we explicitly defined a “bias-voltage matrix” \widetilde{\underline{V}}_{\lambda}\equiv V_{\lambda}\underline{\mbox{1}}-(\widetilde{\underline{h}}_{CC}-\underline{h}_{CC}). The eigenvalues $\{{\epsilon},\widetilde{{\epsilon}}\}$ and eigenvectors $\{{\varPsi}^{\text{L/R}},\widetilde{{\varPsi}}^{\text{L/R}}\}$ refer to the complex eigenvalues and to the left/right eigenvectors of $\underline{h}_{\mathrm{eff}}$ and $\widetilde{\underline{h}}_{\mathrm{eff}}=\widetilde{\underline{h}}_{CC}-{\mathrm{i}}\underline{\varGamma}/2$, respectively. In the limit $\widetilde{\underline{h}}_{CC}\to\underline{h}_{CC}$ the result in Eq. (S35) can be checked to reduce to Eq. (S10) Tuovinen (2016).

Similarly, for the TD1RDM with sudden electromagnetic fields in the central region in Eq. (S35), we can take the integrals in Eqs. (S37), (S38) and (S39) and evaluate them in the same manner. This time the pole structure is only a little more intricate due to different eigenvalues for the unperturbed and perturbed Hamiltonians but it can be handled exactly in the same way as above. For perturbed central regions at arbitrary $\beta$ the explicit results are

[TABLE]

Again, inserting Eqs. (S40), (S41) and (S42) into Eq. (S35) gives then the TD1RDM for a perturbed central region at arbitrary temperature. Also here, the zero-temperature limit ($\beta\to\infty$) presented in Ref. Tuovinen et al. (2014b), is recovered by the asymptotics of the digamma and hypergeometric functions in Eqs. (S40), (S41) and (S42) Tuovinen (2016). By careful inspection of Eqs. (S40), (S41) and (S42) in the limit of unperturbed central region ($\widetilde{{\varPsi}}\to{\varPsi}$ and $\widetilde{{\epsilon}}\to{\epsilon}$) it can be verified that they reduce to Eqs. (S19), (S28) and (III) Tuovinen (2016).

V Derivation of the bond current

We define the bond current flowing between site $j$ and $j+1$ in the nanowire (central device) by the rate of change of the number of particles in the region comprising the left electrode and the first $j$ sites in the nanowire:

[TABLE]

where $k$ and $m$ respectively label the basis elements in the left electrode and the sites in the nanowire, and $s$ is a spin index. The bond current between sites $j$ and $j+1$ is then defined by

[TABLE]

The temporal change in the number of particles can be derived from

[TABLE]

where we now separate the ‘normal’ and ‘superconducting’ contributions as $\hat{H}_{\text{total}}=\hat{H}_{\text{normal}}+\hat{H}_{\varDelta}$ with

[TABLE]

where we separated the “on-site” and “nearest-neighbor” contributions in the spin-dependent matrix elements of $\epsilon$ and $J$. Also, the nanowire is coupled to the electrodes only via the terminal sites ($1$ and $N_{w}$), so the coupling Hamiltonian takes the form

[TABLE]

With the normal part of the Hamiltonian, the commutator in Eq. (S45) is nonzero only for the following terms

[TABLE]

As the coupling Hamiltonian $\hat{H}_{c}$ has one creation (annihilation) operator in the nanowire and one annihilation (creation) operator in the electrode, so in principle the first two terms above can give a nonzero commutator, but it turns out they cancel each other out. We are then left with the term on the second line which can be simplified to give

[TABLE]

The remaining calculation is the commutator with the ‘superconducting’ part where the nonzero contribution comes from

[TABLE]

We may then insert Eqs. (S50) and (S51) into Eq. (S45) and further into Eq. (S44), and use the model parameters for the nanowire (we have assumed a real pairing field) to obtain:

[TABLE]

This may further be simplified as in Eq. (11) in the main text.