Nonlocality and dynamic response of Majorana states in fermionic superfluids

Search for Majorana bound states (MBS) has recently become an active topic in the condensed matter community [1–3]. These exotic states are known to be characterized by the coinciding annihilation and creation operators. This is why it is quite natural to look for such states in superconducting systems where the order parameter Δ is known to mix particles (electrons) and anti-particles (holes) because of the Andreev scattering processes. Standard singlet superconductivity still does not allow the formation of this kind of excitations while the more exotic triplet state can host MBS. Among the available superfluids there exist only a few possible candidates for the triplet pairing such as He-3, ${{\rm{Sr}}}_{2}{{\rm{RuO}}}_{4}$ and heavy fermion compounds [4, 5]. Alternatively, the effective triplet pairing can be induced, e.g., in semiconducting nanowires [6, 7] in the presence of rather strong spin–orbit coupling and external magnetic field. Despite the clear and reliable observation of zero bias peaks (ZBP) in the differential conductance measurements [8, 9] and on the change in the charge periodicity of conductance in Coulomb blockade regime [10] consistent with the existence of MBS it would be extremely important to probe other attributes of these states especially keeping in mind alternative explanations of the ZBP based on Kondo physics [11].

The goal of this paper is to suggest a test revealing the nonlocal dynamic response of the MBS. This issue has recently become a subject of intensive debate in the context of so-called quantum teleportation [12–16]. The Majorana partner states are localized at the length scales of the order of the coherence length ξ and are usually strongly separated provided the distance L between them well exceeds this length ξ (see figure 1). From the standard quantum mechanics one could naively expect that the time ${\tau }_{0}$ of the particle transfer between these localized states should be determined by the inverse tunneling rate roughly proportional to the value ${\rm{\Delta }}{{\rm{e}}}^{-L/\xi }$. Such scenario can be questioned if we remind that two Majorana states form a single fermionic level and, thus, the injected particle should appear simultaneously in both partner states [12–14]. This conclusion is in obvious contradiction with the analysis of the current noise correlations [15, 16]: the latter points towards the existence of a finite charge transfer time between the MBS. Later on the teleportation phenomenon has been argued to be restored due to the nonlocal coupling via the Coulomb blockage [14]. It was concluded that the key omission of the previous studies was related to the treating of the superconducting phase as a constant, and not as a dynamic variable. According to the work [14] the recovering of the nonlocal coupling between the MBS should occur if we consider the phase of the superconducting order parameter as a quantum variable canonically conjugate to the charge of the island.

Zoom In Zoom Out Reset image size

In the present manuscript we show that the previous studies of the nonlocality in the system of the MBS suffer from another key omission, namely they do not take into account the nonequilibrium effects responsible for the mixing of the quasiparticle eigenfunctions with the positive and negative energies in the dynamic processes. In the remaining part of the paper we consider a model describing the corresponding low frequency dynamics of the MBS and make clear predictions for the time-dependent experiment suggested above. Specifically, our analysis demonstrates that the time of the quasiparticle transfer between Majorana states should be of the order of the inverse energy splitting ${\tau }_{0}\sim {\omega }_{0}^{-1}$ caused by their coupling ${\omega }_{0}$. This result imposes restrictions on the time scales of adiabatic manipulation of the Majorana states giving a criterion of their topological protection in time-dependent phenomena. For comparison it is interesting to mention here the work [17] where the dynamics is governed by time of flight of excitations in the normal metal wire coupled to the MBS.

The low frequency dynamics of quasiparticles (QPs) can be described within the time-dependent generalization of the BdG equations (see [18])

$$\begin{eqnarray}{\rm{i}}\displaystyle \frac{\partial }{\partial t}{\hat{g}}_{n}=\left(\begin{array}{cc}{\hat{H}}_{0}-\mu & \hat{{\rm{\Delta }}}\\ {\hat{{\rm{\Delta }}}}^{\dagger } & \mu -{\hat{H}}_{0}^{* }\end{array}\right){\hat{g}}_{n}.\end{eqnarray} \tag{ 1 }$$

Here ${\hat{H}}_{0}$ is the normal state Hamiltonian, μ is the chemical potential, and ${\hat{g}}_{n}({\bf{r}},t)=({u}_{\alpha ,n},{v}_{\alpha ,n})$. The condition of adiabaticity naturally assumes that all the characteristic frequencies are much lower than the superconducting gap Δ, otherwise a full nonequilibrium description of a superconductor should be applied [19]. The coefficients ${u}_{\alpha ,n}$ and ${v}_{\alpha ,n}$ are usually interpreted as electronic- and hole- like parts of the QP wave functions defined by the Bogolubov transformation,

$$\begin{eqnarray}&&{\hat{{\rm{\Psi }}}}_{\alpha }({\bf{r}},t)=\displaystyle \sum _{n}({u}_{\alpha ,n}({\bf{r}},t){\hat{c}}_{n}+{v}_{\alpha ,n}^{* }({\bf{r}},t){\hat{c}}_{n}^{\dagger }),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\hat{{\rm{\Psi }}}}_{\alpha }^{\dagger }({\bf{r}},t)=\displaystyle \sum _{n}({u}_{\alpha ,n}^{* }({\bf{r}},t){\hat{c}}_{n}^{\dagger }+{v}_{\alpha ,n}({\bf{r}},t){\hat{c}}_{n}).\end{eqnarray} \tag{ 3 }$$

Here α is the spin index and ${\hat{c}}_{n}^{\dagger }$, ${\hat{c}}_{n}$ are the fermionic QP creation and annihilation operators, respectively. The index n enumerates the solutions of time-dependent BdG equations for different initial conditions at t = 0 when the expressions (2) take the form of expansion over a certain full set of functions. In equilibrium the time dependence of the wave functions reduces to the standard form ${u}_{\alpha ,n}({\bf{r}},t)={\bar{u}}_{\alpha ,n}({\bf{r}}){{\rm{e}}}^{-{{\rm{i}}{E}}_{n}t}$, ${v}_{\alpha ,n}({\bf{r}},t)={\bar{v}}_{\alpha ,n}({\bf{r}}){{\rm{e}}}^{-{{\rm{i}}{E}}_{n}t}$, where E_n and $({\bar{u}}_{\alpha ,n}({\bf{r}}),{\bar{v}}_{\alpha ,n}({\bf{r}}))$ are the spectrum and eigenfunctions of the stationary BdG equations. Only the states with ${E}_{n}\geqslant 0$ contribute to the equation (2) in this limit while in general the time-dependent solutions ${\hat{g}}_{n}({\bf{r}},t)$ may contain the contributions from all positive and negative levels of the stationary Hamiltonian.

The Majorana-type states in the stationary case can appear provided we have an isolated eigenfunction satisfying the condition ${v}_{\alpha ,0}^{* }={u}_{\alpha ,0}$ corresponding to zero energy. The inverse transformation for this zero energy state can specify only the sum of the fermionic operators

$$\begin{eqnarray}&&\displaystyle \frac{{\hat{c}}_{0}+{\hat{c}}_{0}^{\dagger }}{2}=\displaystyle \sum _{\alpha }\int {\rm{d}}{\bf{r}}({u}_{\alpha ,0}^{* }({\bf{r}}){\hat{{\rm{\Psi }}}}_{\alpha }({\bf{r}})+{u}_{\alpha ,0}({\bf{r}}){\hat{{\rm{\Psi }}}}_{\alpha }^{\dagger }({\bf{r}})).\end{eqnarray} \tag{ 4 }$$

This relation does not naturally yield the full fermionic operator ${\hat{c}}_{0}={\hat{\gamma }}_{L}+{\rm{i}}{\hat{\gamma }}_{R}$ but only its part ${\hat{\gamma }}_{L}=({\hat{c}}_{0}+{\hat{c}}_{0}^{\dagger })/2$ which indeed meets the Majorana conditions. Another part (${\hat{\gamma }}_{R}$) of the QP operator remains undefined and in this sense the ground state of the superconductor with an isolated zero energy mode appears to be degenerate. The ambiguity of the operator ${\hat{\gamma }}_{R}$ can be resolved by introducing a coupling mechanism of the above isolated state either to the second Majorana-type state or to a fermionic bath [16, 20]. Both these mechanisms destroy the symmetry of the isolated level ${v}_{\alpha ,0}^{* }={u}_{\alpha ,0}$ and shift its energy from zero. Each Majorana pair of states gives one positive and one negative energy level. In equilibrium it is natural to keep only the positive energy level and the corresponding hybridized wave function. Considering the nonequilibrium dynamics at a finite time interval t we can no more disregard the contribution of the negative energy level to the wave function dynamics when the energy uncertainty $\delta E\sim {\hslash }/t$ exceeds the splitting of levels in a Majorana pair. Thus, despite of the obvious fact that both levels correspond to the only fermion the nonequilibrium time-dependent solutions ${\hat{g}}_{n}({\bf{r}},t)$ of the BdG equations contain contributions corresponding to both levels.

To probe the nonlocal dynamics of coupled Majorana states we suggest to study transport through the wire hosting these MBS at its ends modulated by the changes in the coupling of the wire to the external normal metal leads (see figure 1). A natural way to tune this coupling in conditions of the real experiment (see, e.g., [10, 21, 22]) is to apply time-dependent voltages at the gate electrodes controlling the transparencies ${t}_{L,R}(t)$ of the barriers between the wire and the normal lead at the left (L) and right (R) end, respectively. Tuning these transparencies at the two ends of the wire one can easily determine the spatial correlations in the dynamic response of the Majorana partners as well as the scale $1/{\tau }_{0}$ of the frequency dispersion. Considering a possible experimental setup based on a semiconducting nanowire with induced superconductivity one should take this system in a topologically nontrivial state [6, 7] which allows to get the subgap quasiparticle states bound to the wire ends.

Further derivation has been carried out by applying a general approach [23] for the solution of the scattering problem with the quasiparticle waves incoming from the left or right leads at a certain energy $\varepsilon$ and propagating along the one-dimensional p-wave superconducting wire hosting two MBS. We focus here on the case of a weak charging energy of the wire which is different from the situation studied in [14]. The p-wave order parameter is chosen in the form ${\rm{\Delta }}(x)\sim {{\rm{e}}}^{{\rm{i}}{\theta }_{p}}$, where ${\theta }_{p}=0,\pi$ is the trajectory orientation angle. Assuming low energies ($\varepsilon ,{\omega }_{0}\ll {\rm{\Delta }}$) and considering the solution of equation (1) near the left end of the wire one can write it as a superposition

$$\begin{eqnarray}&&g({\bf{r}},t)={{\rm{e}}}^{-{\rm{i}}\varepsilon t+{{\rm{i}}{k}}_{F}s}[{a}_{L}^{+}{w}^{(1)}(s)+{b}_{L}^{+}{w}^{(2)}(s)]+{{\rm{e}}}^{-{\rm{i}}\varepsilon t-{{\rm{i}}{k}}_{F}s}[{a}_{L}^{-}{w}^{(1)}(-s)+{b}_{L}^{-}{w}^{(2)}(-s)]\end{eqnarray} \tag{ 5 }$$

of two independent solutions

$$\begin{eqnarray}&&{w}^{(1)}(s)={{\rm{e}}}^{{\rm{i}}{\hat{\sigma }}_{z}{\theta }_{p}/2}\left[{{\rm{e}}}^{-D(s)/2}\left(\begin{array}{c}1\\ -{\rm{i}}\end{array}\right)\right.+\left.{\rm{i}}\displaystyle \frac{\varepsilon }{\tilde{{\rm{\Delta }}}}\mathrm{sign}(s){{\rm{e}}}^{D(s)/2}\left(\begin{array}{c}1\\ {\rm{i}}\end{array}\right)\right],\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{w}^{(2)}(s)={{\rm{e}}}^{{\rm{i}}{\hat{\sigma }}_{z}{\theta }_{p}/2}{{\rm{e}}}^{D(s)/2}\left(\begin{array}{c}1\\ {\rm{i}}\end{array}\right),\end{eqnarray} \tag{ 7 }$$

found in [24, 25] for the quasiclassical Andreev equations at the trajectory with the coordinate $s=(L/2)\cos {\theta }_{p}+x$. A similar expression can be written near the right end of the wire by changing the subscripts $L\to R$ and the angle ${\theta }_{p}$ from 0 to π, which shifts the origin $x\to x-L$ corresponding to $s({\theta }_{p}=0)\gt 0$ and $s({\theta }_{p}=\pi )\lt 0$.

Here v_F is the Fermi velocity in the wire, ${\tilde{{\rm{\Delta }}}}^{-1}=\tfrac{2}{{v}_{F}}{\int }_{0}^{L/2}{{\rm{e}}}^{-D(s)}{\rm{d}}{s}$, $D(s)=\tfrac{2}{{v}_{F}}\left|{\int }_{0}^{s}{\rm{\Delta }}(s^{\prime} ){\rm{d}}{s}^{\prime} \right|\sim \tfrac{| s| }{\xi }$, and Pauli matrices ${\hat{\sigma }}_{k}$ act in the electron–hole Gor'kov–Nambu space. An appropriate matching of the wavefunctions at the wire ends with the ones in the leads gives us the equations for the coefficients ${a}_{k}^{\pm }={{\rm{e}}}^{\pm {\rm{i}}{\phi }_{k}/2}({A}_{k}\pm {a}_{k})/2$ at the left (k = L) and right (k = R) wire ends (see appendix A for details of calculations)

$$\begin{eqnarray}&&({{\rm{\Gamma }}}_{k}-{\rm{i}}\varepsilon ){A}_{k}={F}_{k},\ (\tilde{{\rm{\Delta }}}-{\rm{i}}\varepsilon {{\rm{\Gamma }}}_{k}/\tilde{{\rm{\Delta }}}){a}_{k}={F}_{k}.\end{eqnarray} \tag{ 8 }$$

Here for simplicity we neglect the MBS coupling ${\omega }_{0}\sim \tilde{{\rm{\Delta }}}{{\rm{e}}}^{-L/\xi }$, ${{\rm{\Gamma }}}_{k}=\tilde{{\rm{\Delta }}}(1-{r}_{k})/(1+{r}_{k})$ is the rate characterizing the coupling of wire states to the kth external lead with ${r}_{k}=\sqrt{1-| {t}_{k}{| }^{2}}$ being the real-valued reflection coefficient of the insulating barrier, ${\phi }_{k}$ are the scattering phases. ${F}_{k}=\tilde{{\rm{\Delta }}}{t}_{k}/(1+{r}_{k})\propto \sqrt{{{\rm{\Gamma }}}_{k}\tilde{{\rm{\Delta }}}}$ are the tunneling sources characterizing the incoming QP flows. Applying the Fourier transform with respect to the energy variable $\varepsilon$ and considering the parameter ${\omega }_{0}/\tilde{{\rm{\Delta }}}\sim {{\rm{e}}}^{-L/\xi }$ pertubatively one can obtain the equations describing the dynamics of a model two-level system in the time frame (see [26, 27]), i.e. the dynamics of the Majorana pair:

$$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{\partial t}+{{\rm{\Gamma }}}_{L}\right){A}_{L}+{\omega }_{0}{A}_{R}={F}_{L}{{\rm{e}}}^{-{\rm{i}}\varepsilon t},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{\partial t}+{{\rm{\Gamma }}}_{R}\right){A}_{R}-{\omega }_{0}{A}_{L}={F}_{R}{{\rm{e}}}^{-{\rm{i}}\varepsilon t}.\end{eqnarray} \tag{ 10 }$$

In the non-stationary regime the localized states at the wire ends (being of Majorana nature in the stationary regime) can be described by the wave function amplitudes A_k which are in fact the quantum mechanical amplitudes describing the probability to find the quasiparticle at the kth wire end. The amplitudes a_k correspond to the off-resonant fast-decaying contributions from the states above the gap. The amplitudes A_k and a_k together describe in fact the low frequency dynamics of the function ${\hat{g}}_{n}({\bf{r}},t)$ including contributions from positive and negative levels of the stationary Hamiltonian. Note that in the absence of incoming QP flows, F_k = 0, equations (9), (10) have purely real-valued coefficients corresponding to the Hermitian nature of Majorana operators ${\hat{\gamma }}_{k}$. In this case the average $\langle {{\rm{\Psi }}}_{\alpha }^{\dagger }({\bf{r}},t){{\rm{\Psi }}}_{\alpha }({\bf{r}},t)\rangle$ of the electron number operator is conserved since its change is determined by the sum $| {A}_{L}{| }^{2}+| {A}_{R}{| }^{2}$ of probabilities $| {A}_{k}{| }^{2}$ to find the quasiparticle at the kth wire end. This conservation fixes, in particular, the quasiparticle parity number in the wire by fixing the parameter $| {A}_{L}{| }^{2}+| {A}_{R}{| }^{2}$ even for non-trivial dynamics of $| {A}_{k}{| }^{2}$ themselves. Note that this statement is independent of a strength of Coulomb interaction as the latter only governs the correlations between tunneling rates. The rates ${{\rm{\Gamma }}}_{L,R}$ are determined by the local Andreev reflection processes [16] while the energy splitting of coupled Majorana states ${\omega }_{0}=\tilde{{\rm{\Delta }}}{{\rm{e}}}^{-D(L/2)}\sin \varphi$ is related to the probability of the quasiparticle transfer through the system. Parameters $D(L/2)\sim L/\xi$ and $\varphi ={k}_{F}L+({\phi }_{L}-{\phi }_{R})/2$ depend on the wire length L.

The current flowing from the left and right electrodes can be calculated as [28] (see also appendix B for details of calculations)

$$\begin{eqnarray}&&{I}_{L,R}=e/\pi \int {g}_{L,R}(\varepsilon )({f}_{T}(\varepsilon -{{\rm{eV}}}_{L,R})-{f}_{T}(\varepsilon -{{\rm{eV}}}_{s})){\rm{d}}\varepsilon ,\end{eqnarray} \tag{ 11 }$$

where ${f}_{T}(\varepsilon )={({{\rm{e}}}^{\varepsilon /T}+1)}^{-1}$ is the Fermi–Dirac distribution function with the bath temperature T,

$$\begin{eqnarray}&&{g}_{k}(\varepsilon )\simeq 2{\rm{Re}}[{A}_{k}{a}_{k}^{* }]=2\sqrt{{{\rm{\Gamma }}}_{k}}{\rm{Re}}({A}_{k}{{\rm{e}}}^{{\rm{i}}\varepsilon t}),\end{eqnarray} \tag{ 12 }$$

V_k is the potential of the kth electrode, and V_s is the potential of a superconductor. Generally, the definition of the potential V_s in a nonstationary problem follows from the solution of the equations describing the particular electric circuit [29], e.g., the one in figure 1: ${I}_{L}+{I}_{R}={{C}{\rm{d}}{V}}_{s}/{\rm{d}}{t}+{V}_{s}/R$, where C and R are the capacitance and shunt resistance of the ground connection, respectively. Considering a constant applied bias $V={V}_{L}-{V}_{R}$ and putting ${A}_{L,R}\propto {{\rm{e}}}^{-{\rm{i}}\varepsilon t}$ we obtain a dc differential conductance peak at ${eV}\simeq {\omega }_{0}$ attributed to MBS [8, 9, 30–32].

We now proceed with the analysis of the dynamic response of a pair of Majorana partners and consider two generic examples of the time-dependent transport realized by the modulating tunnel barrier (see figure 1): (i) the phase-shifted sinusoidal driving with ${{\rm{\Gamma }}}_{L}(t)={{\rm{\Gamma }}}_{0}+\tilde{{\rm{\Gamma }}}\cos (\omega t)$ and ${{\rm{\Gamma }}}_{R}(t)={{\rm{\Gamma }}}_{0}+\tilde{{\rm{\Gamma }}}\cos (\omega t+{\phi }_{0});$ (ii) pump-probe driving by ${\rm{\Delta }}t$-broadened delta-functional pulses with different amplitudes G_k^t applied with a time delay τ, i.e., with ${{\rm{\Gamma }}}_{k}(t)={G}_{k}^{0}{\delta }_{{\rm{\Delta }}t}(t)+{G}_{k}^{\tau }{\delta }_{{\rm{\Delta }}t}(t-\tau )$.

To start with, our consideration of the dynamic response of MBS within equations (9), (10) through a single fermionic state formed of a superposition of two partner Majorana states. Indeed, the levels $\pm {\omega }_{0}$ around the zero energy can be introduced as a basis of hybridized states with the amplitudes ${A}_{\pm }={A}_{L}\pm {{\rm{i}}{A}}_{R}$. In equations (9), (10) each of quasiparticle sources F_k excites both amplitudes ${A}_{\pm }$ simultaneously. Due to the coupling to the reservoirs both amplitudes evolve then in time as separate quantities and, thus, cannot be described as an empty and filled state of a single level. As a result, we find beating of the wavefunction between the edge states at the frequency ${\omega }_{0}$. The above arguments concerning the sources of the injected particles should be valid irrespective to the strength of the Coulomb effects and for the sake of simplicity we start our consideration of time-dependent problems from the limit of large capacitance C when these effects can be neglected.

Starting from the case of sinusoidal driving we consider for simplicity the ac amplitude $\tilde{{\rm{\Gamma }}}\ll {{\rm{\Gamma }}}_{0}$ as a perturbation and solve equations (9), (10). For the zero-bias differential conductance we find

$$\begin{eqnarray}&&{\left.\displaystyle \frac{1}{{G}_{{\rm{\Gamma }}}}\displaystyle \frac{{{\rm{d}}{I}}_{L}}{{{\rm{d}}{V}}_{L}}\right|}_{{V}_{L}=0}\simeq 1+\displaystyle \frac{\tilde{{\rm{\Gamma }}}}{2{{\rm{\Gamma }}}_{0}}\left[\cos \omega t+\left.\displaystyle \frac{{\omega }_{0}^{2}-{{\rm{\Gamma }}}_{0}^{2}}{2{{\rm{\Gamma }}}_{0}}\displaystyle \sum _{\eta =\pm 1}{{ \mathcal L }}_{\eta }{F}_{0}^{\eta }-{\omega }_{0}\displaystyle \sum _{\eta =\pm 1}\eta {{ \mathcal L }}_{\eta }{F}_{{\phi }_{0}-\pi /2}^{\eta }\right]\right.,\end{eqnarray} \tag{ 13 }$$

where ${G}_{{\rm{\Gamma }}}=({e}^{2}/\pi )2{{\rm{\Gamma }}}_{0}^{2}/({{\rm{\Gamma }}}_{0}^{2}+{\omega }_{0}^{2})$, ${{ \mathcal L }}_{\eta }={{\rm{\Gamma }}}_{0}/[{(\omega +\eta {\omega }_{0})}^{2}+{{\rm{\Gamma }}}_{0}^{2}]$, ${F}_{\phi }^{\pm }=\cos (\omega t+\phi )+\sin (\omega t+\phi )(\omega \pm {\omega }_{0})/{{\rm{\Gamma }}}_{0}$. One can see that for low-frequencies $\omega \lesssim {\omega }_{0}$ the above expression contains an essential phase ${\phi }_{0}$ dependence, while with increasing ω these contributions decay faster than the other time-dependent terms.

Indeed, this statement is clearly visible in the most interesting and representative case ${\omega }_{0}\sim {{\rm{\Gamma }}}_{0}$ in which dc results [30–32] (see also (C.1) in appendix C) are already broadened and inconclusive. In this case to clarify the results we rearrange the functions ${F}_{\phi }^{\pm }=\cos (\omega t+\phi )+\sin (\omega t+\phi )(\omega \pm {\omega }_{0})/{{\rm{\Gamma }}}_{0}={F}_{\phi }^{c}\pm {F}_{\phi }^{s}$ to ${F}_{\phi }^{c}=\cos (\omega t+\phi )+(\omega /{{\rm{\Gamma }}}_{0})\sin (\omega t+\phi )$ and ${F}_{\phi }^{s}=({\omega }_{0}/{{\rm{\Gamma }}}_{0})\sin (\omega t+\phi )$ getting

$$\begin{eqnarray}&&{\left.\displaystyle \frac{1}{{G}_{{\rm{\Gamma }}}}\displaystyle \frac{{{\rm{d}}{I}}_{L}}{{{\rm{d}}{V}}_{L}}\right|}_{{V}_{L}=0}\simeq 1+\displaystyle \frac{\tilde{{\rm{\Gamma }}}}{2{{\rm{\Gamma }}}_{0}}\cos \omega t+\displaystyle \frac{\tilde{{\rm{\Gamma }}}}{{{\rm{\Gamma }}}_{0}}\left[\displaystyle \frac{{\omega }_{0}^{2}-{{\rm{\Gamma }}}_{0}^{2}}{2{\omega }_{0}^{2}}\{{C}_{0}{F}_{0}^{c}-{C}_{1}{F}_{0}^{s}\}\right.-\left.\displaystyle \frac{{{\rm{\Gamma }}}_{0}}{{\omega }_{0}}\{{C}_{0}{F}_{{\phi }_{0}-\pi /2}^{s}-{C}_{1}{F}_{{\phi }_{0}-\pi /2}^{c}\}\right],\end{eqnarray} \tag{ 14 }$$

with the dimensionless coefficients ${C}_{k}=(\pi {\omega }_{0}^{2}/2{{\rm{\Gamma }}}_{0}){\sum }_{\eta =\pm 1}{(-\eta )}^{k}{{ \mathcal L }}_{\eta }$ accumulating the ω-dependence of the prefactors as follows

$$\begin{eqnarray}&&{C}_{0}(\omega )=\displaystyle \frac{({\omega }^{2}+{\omega }_{0}^{2}+{{\rm{\Gamma }}}_{0}^{2}){\omega }_{0}^{2}}{[{(\omega +{\omega }_{0})}^{2}+{{\rm{\Gamma }}}_{0}^{2}][{(\omega -{\omega }_{0})}^{2}+{{\rm{\Gamma }}}_{0}^{2}]},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{C}_{1}(\omega )=\displaystyle \frac{2\omega {\omega }_{0}^{3}}{[{(\omega +{\omega }_{0})}^{2}+{{\rm{\Gamma }}}_{0}^{2}][{(\omega -{\omega }_{0})}^{2}+{{\rm{\Gamma }}}_{0}^{2}]}.\end{eqnarray} \tag{ 16 }$$

Now considering two limits: (a) $\omega \sim {\omega }_{0},{{\rm{\Gamma }}}_{0}$ and (b) $\omega \gg {\omega }_{0},{{\rm{\Gamma }}}_{0}$ illustrated in the corresponding panels of figure 2, one can see that in the first limit (a) $\omega \sim {\omega }_{0},{{\rm{\Gamma }}}_{0}$ the above mentioned coefficients ${C}_{\mathrm{0,1}}\sim 1$ weakly depend on the frequency ω and the phase dependent corrections (the last line in (14)) are of order of the main term $\tfrac{\tilde{{\rm{\Gamma }}}}{2{{\rm{\Gamma }}}_{0}}\cos \omega t$.

Zoom In Zoom Out Reset image size

In the second limit (b) due to the smallness of ${C}_{0}\simeq {\omega }_{0}^{2}/{\omega }^{2}$ and ${C}_{1}\simeq 2{\omega }_{0}^{3}/{\omega }^{3}$ the conductance has relatively small ${\omega }_{0}^{2}/{\omega }^{2}$ phase-dependent corrections to the oscillating terms

$$\begin{eqnarray}&&{\left.\displaystyle \frac{1}{{G}_{{\rm{\Gamma }}}}\displaystyle \frac{{{\rm{d}}{I}}_{L}}{{{\rm{d}}{V}}_{L}}\right|}_{{V}_{L}=0}\simeq \displaystyle \frac{\tilde{{\rm{\Gamma }}}}{2{{\rm{\Gamma }}}_{0}}\left[\cos \omega t+\displaystyle \frac{{\omega }_{0}^{2}-{{\rm{\Gamma }}}_{0}^{2}}{{\omega }^{2}}{F}_{0}^{c}\right.-\left.\displaystyle \frac{2{\omega }_{0}{{\rm{\Gamma }}}_{0}}{{\omega }^{2}}{F}_{{\phi }_{0}-\pi /2}^{s}\right].\end{eqnarray} \tag{ 17 }$$

In figure 2 the $(t,\phi )$-dependence of the differential conductance (14) is plotted for the following parameters ${\omega }_{0}={{\rm{\Gamma }}}_{0}=10\tilde{{\rm{\Gamma }}}$ at (a) $\omega ={\omega }_{0}$ and (b) $\omega =10{\omega }_{0}$ demonstrating the above mentioned arguments.

In the other limit of ${\omega }_{0}\gg {{\rm{\Gamma }}}_{0}$ many beating periods pass before a tunneling event occurs leading to the efficient transport of the charge between the localized states A_k. This can be in some sense viewed as a signature of 'teleportation'. If additionally ${\rm{\Delta }}\omega =\omega -{\omega }_{0}\sim {{\rm{\Gamma }}}_{0}\ll {\omega }_{0}$ one can neglect the contributions from $\eta =+1$ and obtain

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\pi }{{e}^{2}}\displaystyle \frac{{{\rm{d}}{I}}_{L}}{{{\rm{d}}{V}}_{L}}\right|}_{{V}_{L}=0}\simeq \displaystyle \frac{2{{\rm{\Gamma }}}_{0}^{2}}{{\omega }_{0}^{2}}+\displaystyle \frac{\tilde{{\rm{\Gamma }}}{{\rm{\Gamma }}}_{0}}{{\omega }_{0}^{2}}\cos \omega t+\displaystyle \frac{\tilde{{\rm{\Gamma }}}}{2{{\rm{\Gamma }}}_{0}}\displaystyle \frac{{{\rm{\Gamma }}}_{0}^{2}}{{\rm{\Delta }}{\omega }^{2}+{{\rm{\Gamma }}}_{0}^{2}}\left[{F}_{0}^{-}+\displaystyle \frac{2{{\rm{\Gamma }}}_{0}}{{\omega }_{0}}{F}_{{\phi }_{0}-\pi /2}^{-}\right].\end{eqnarray} \tag{ 18 }$$

Clearly this limit describes the sharp peaks at ${\omega }_{0}$ in the frequency dependence of the dynamic response with the amplitude that depends on the phase shift. In the opposite limit of broad peaks the nonlocal correlations in the dynamic response are naturally more difficult to observe since their contributions in the dynamic response become small when ${\omega }_{0}/{{\rm{\Gamma }}}_{0}\ll 1$.

For arbitrary bias and drive amplitudes we should get a multiplication of harmonics and considering the current averaged over the drive period we can expect the appearance of the conductance peaks at voltages ${{eV}}_{L,R}=n\omega \pm {\omega }_{0}$ due to the resonant effect similar to the Shapiro phenomenon in Josephson junctions [33]. Note that the periodic backgate voltage modulation can give another opportunity to observe the resonant features on the current–voltage curve controlling the chemical potential of the wire as a whole. This modulation should cause the change in the energy splitting ${\omega }_{0}$ through its dependence on the Fermi momentum k_F. Assuming ${k}_{F}L=\omega t$ to be linear in time one can obtain resonances at ${{eV}}_{L,R}=n\omega$.

In the case of the pump-probe driving the differential conductance of the left electrode contains three contributions

$$\begin{eqnarray}&&\displaystyle \frac{\pi }{{e}^{2}{\rm{\Delta }}t}\displaystyle \frac{{{\rm{d}}{I}}_{L}}{{{\rm{d}}{V}}_{L}}={G}_{L}^{0}{\delta }_{{\rm{\Delta }}t}(t)+{G}_{L}^{\tau }{\delta }_{{\rm{\Delta }}t}(t-\tau )+\sqrt{{G}_{L}^{0}{G}_{L}^{\tau }}{\rm{\cos }}({\omega }_{0}\tau ){\rm{\cos }}({{eV}}_{L}\tau ){\delta }_{{\rm{\Delta }}t}(t-\tau ).\end{eqnarray} \tag{ 19 }$$

Here we impose zero initial conditions on both amplitudes ${A}_{\pm }$. The terms in the first line of equation (19) correspond to the local charging of the single fermionic level, while the term on the second line reflects correlations in the response to two pulses with the time delay τ and shows the non-trivial dynamics of MBS at frequencies $| {{eV}}_{L}\pm {\omega }_{0}|$ (see figure 3). The first pulse excites the quantum beatings between the Majorana edge states at the frequency ${\omega }_{0}$ modifying the response of the system to the second pulse.

Zoom In Zoom Out Reset image size

Taking for the estimate ${\rm{\Delta }}\sim 2.5$ K, $\xi \sim 100\,\mathrm{nm}$ for Al, and $L\gtrsim 1$ μm we find ${\omega }_{0}\lesssim 15\,\mathrm{MHz}$ which gives us a reasonable range of frequencies $\omega \sim {\omega }_{0}$ of the drive and typical time delay $1/{\omega }_{0}\sim 0.06$ μs for the pump-probe setup. The conditions on bias for the observation of the beating phenomenon are less restrictive comparing to the ones of a dc conductance peak (with the restriction of $V\sim {\omega }_{0}\sim 0.01$ μV [32]) as ac measurements are already conclusive at zero bias, see, e.g., equation (13). To get ${{\rm{\Gamma }}}_{L,R}\lesssim {\omega }_{0}$ we should take the barriers with resistances ${R}_{L,R}\sim 0.1\mbox{--}1$ GΩ.

Certainly, the above dynamic response of the MBS will be modified in Coulomb blockade regime. This difference arises from the obvious fact that in the case of Coulomb blockade the charge tunneling processes between the island and left/right electrodes are strongly correlated. The entry and exit of charged particles are always controlled by the overall charge of the island. However, this correlation does not destroy the beating phenomenon and cannot cause the formation of a single eigenstate responsible for the non-local transport through Majorana states (teleportation) for the operating frequencies above the energy splitting of Majorana partners. Let us take the limit of high Coulomb energy and imposing, thus, the restriction on two possible charge states of the island and assume the operating frequencies and energy splitting ${\omega }_{0}$ to be small comparing to the tunneling rates ${{\rm{\Gamma }}}_{L,R}$. The latter limit allows one to consider the charging/discharging processes as instantaneous events changing the fermion parity. On the longer time scales than the injection/ejection rates the fermion parity is fixed due to the fixed electron charge. However, the beating phenomenon as an internal dynamics of Majorana states is present due to the nonequilibrium time-dependent nature of the electron injection and further transformation of the wave function of the injected electron into the Andreev eigenstates both with positive and negative energies. Therefore the current through the system is fully determined by the interplay of two time scales, namely, the inverse beating frequency ${\omega }_{0}^{-1}$ and the delay time τ between the opening of the left/right junctions. The latter is determined either by the operating frequency f and the phase shift ϕ ($\tau =(n+\phi /2\pi )/f$ with an integer n value) for the periodic driving or by the delay time τ for the pump-probe experimental setup. Certainly the above comment on the influence of Coulomb blockade on the beating phenomenon is only qualitative and should be verified by further quantitative analysis based on the use of more elaborated methods taking account of the interaction effects.

To conclude the solution of the above dynamic problems allows us to predict a beating effect at the frequency ${\omega }_{0}$ which is a hallmark of the topologically nontrivial state of the nanowire. We show that due to the exponentially small coupling ${\omega }_{0}$ the MBS are strongly sensitive to any external perturbation. According to our consideration any driving of Majorana states with the typical operating frequency ω exceeding ${\omega }_{0}$ brings the system to the non-equilibrium regime imposing, thus, an important restriction on the operating frequencies of such a device The Majorana nature of these states needed for quantum calculations recovers only in the adiabatic regime $\omega \ll {\omega }_{0}$. On the other hand, the measurement of the characteristic frequency threshold ${\omega }_{0}$ separating the regimes of weak and strong perturbations of the Majorana pairs could be considered as their hallmark characterizing the nonlocality of these pairs. Certainly the beating phenomenon similar to the one discussed in our work should appear in other superconducting systems with subgap Andreev states. To distinguish the beating phenomenon in topological situation from the one caused by the presence of usual Andreev states it may be helpful to study the behavior of the beating frequency as a function of system parameters, gate potentials and magnetic field so that to reveal the features peculiar to the topologically protected levels. The beating phenomenon may also affect non-stationary Josephson-type transport in systems with MBS studied in recent experiments [34, 35].

We are pleased to thank A A Bespalov, Yu G Makhlin, C Marcus, G E Volovik, and A D Zaikin for valuable comments and A J Leggett for correspondence. This work has been supported in part by Microsoft Project Q, by the Nanosciences Foundation, foundation under the aegis of the Joseph Fourier University Foundation (Grenoble, France), by Academy of Finland Projects No. 284594, 272218 (JPP), by the Russian Foundation for Basic Research and German Research Foundation (DFG) Grant No. KH 425/1-1 (IMK), and by the Russian Science Foundation, Grant No. 17-12-01383 (ASM) and Foundation for the advancement of theoretical physics 'BASIS'.

In this section we present the derivation of the equations (9), (10) from the main text for an exemplary system consisting of a one dimensional (1D) p-wave superconducting (S) wire of the length L connected to the left and right one-dimensional normal-metal leads. We choose the x axis along the wire, the origin to be in the middle of the wire and the order parameter in the form ${\rm{\Delta }}\propto {\hat{k}}_{x}+{\rm{i}}{\hat{k}}_{y}$. Such system is known to host the subgap edge states at rather small energies $\pm {\omega }_{0}\sim \pm {\rm{\Delta }}{{\rm{e}}}^{-L/\xi }$. To describe these localized states we start from the quasiclassical version of the Bogolubov-de Gennes equations, i.e., Andreev equations for the envelopes $w=(u,v)$ of the electron and hole waves propagating along the quasiclassical trajectory ${\bf{k}}={k}_{F}(\cos {\theta }_{p},\sin {\theta }_{p})$

$$\begin{eqnarray}&&-{{\rm{i}}{v}}_{F}{\hat{\sigma }}_{z}\displaystyle \frac{\partial }{\partial s}w+{\hat{\sigma }}_{x}{\rm{\Delta }}(s)w=\varepsilon w,\end{eqnarray} \tag{ A.1 }$$

where v_F is the Fermi velocity in the wire, $s=(L/2)\cos {\theta }_{p}+x$ is the coordinate in the wire along the trajectory, and Pauli ${\hat{\sigma }}_{k}$ matrices act in the electron–hole Gor'kov–Nambu space, Considering the p-wave symmetry of superconducting order parameter one can put ${\rm{\Delta }}(x)\sim {{\rm{e}}}^{{\rm{i}}{\theta }_{p}}$. Note that in 1D geometry of the p-wave S wire it is natural to align the trajectory in the positive or negative direction of the x axis which correspond to ${\theta }_{p}=0,\pi$. The phase ${\theta }_{p}$ can be removed from the gap operator Δ by the standard transformation $u(x)\to u(x){{\rm{e}}}^{{\rm{i}}{\theta }_{p}/2}$ and $v(x)\to v(x){{\rm{e}}}^{-{\rm{i}}{\theta }_{p}/2}$.

A.1. Low energy modes inside the wire

Considering the low energy modes with $\varepsilon \ll \tilde{{\rm{\Delta }}}\sim {\rm{\Delta }}$ inside the wire one can take the sum of two independent solutions ${w}^{(\mathrm{1,2})}(s)$ of Andreev equation (A.1) found in [24, 25]

$$\begin{eqnarray}&&{w}^{(1)}(s)={{\rm{e}}}^{{\rm{i}}{\hat{\sigma }}_{z}{\theta }_{p}/2}\left[{{\rm{e}}}^{-D(s)/2}\left(\begin{array}{c}1\\ -{\rm{i}}\end{array}\right)\right.+\left.{\rm{i}}\displaystyle \frac{\varepsilon }{\tilde{{\rm{\Delta }}}}\mathrm{sign}(s){{\rm{e}}}^{D(s)/2}\left(\begin{array}{c}1\\ {\rm{i}}\end{array}\right)\right],\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&{w}^{(2)}(s)={{\rm{e}}}^{{\rm{i}}{\hat{\sigma }}_{z}{\theta }_{p}/2}{{\rm{e}}}^{D(s)/2}\left(\begin{array}{c}1\\ {\rm{i}}\end{array}\right),\end{eqnarray} \tag{ A.3 }$$

where $\varepsilon$ is the energy variable, ${\tilde{{\rm{\Delta }}}}^{-1}=2{\int }_{0}^{L/2}{{\rm{e}}}^{-D(s)}{\rm{d}}{s}/{v}_{F}$ and

$$\begin{eqnarray}&&D(s)=\displaystyle \frac{2}{{v}_{F}}\left|{\int }_{0}^{s}{\rm{\Delta }}(s^{\prime} ){\rm{d}}{s}^{\prime} \right|\sim \displaystyle \frac{| s| }{\xi }.\end{eqnarray} \tag{ A.4 }$$

The full wave function near the left end of the wire being an eigenfunction of the stationary version of Bogolubov-de Gennes equations (1) in the main text can be written as a combination of the above envelopes with the corresponding oscillating factors ${{\rm{e}}}^{-{\rm{i}}\varepsilon t\pm {\rm{i}}{\bf{k}}\cdot {\bf{r}}}$ for the left and right movers with certain coefficients ${a}_{k}^{\pm }$ and ${b}_{k}^{\pm }$

$$\begin{eqnarray}&&g({\bf{r}},t)={{\rm{e}}}^{-{\rm{i}}\varepsilon t+{{\rm{i}}{k}}_{F}s}[{a}_{L}^{+}{w}^{(1)}(s)+{b}_{L}^{+}{w}^{(2)}(s)]+{{\rm{e}}}^{-{\rm{i}}\varepsilon t-{{\rm{i}}{k}}_{F}s}[{a}_{L}^{-}{w}^{(1)}(-s)+{b}_{L}^{-}{w}^{(2)}(-s)].\end{eqnarray} \tag{ A.5 }$$

A similar expression can be also written near the right end of the wire by changing the subscripts $L\to R$ and the angle ${\theta }_{p}$ from 0 to π, which shifts the origin $x\to x-L$ corresponding to $s({\theta }_{p}=0)\gt 0$ and $s({\theta }_{p}=\pi )\lt 0$. Matching the wave functions of the left and right movers one can find the equations for the coefficients

$$\begin{eqnarray}&&{{\rm{i}}{b}}_{R}^{\pm }={a}_{L}^{\pm }{{\rm{e}}}^{-D(L/2)\pm {{\rm{i}}{k}}_{F}L}\mp \displaystyle \frac{\varepsilon }{\tilde{{\rm{\Delta }}}}{a}_{R}^{\pm },\end{eqnarray} \tag{ A.6 }$$

$$\begin{eqnarray}&&-{{\rm{i}}{b}}_{L}^{\pm }={a}_{R}^{\pm }{{\rm{e}}}^{-D(L/2)\mp {{\rm{i}}{k}}_{F}L}\mp \displaystyle \frac{\varepsilon }{\tilde{{\rm{\Delta }}}}{a}_{L}^{\pm }.\end{eqnarray} \tag{ A.7 }$$

Here for simplicity we assume the following symmetry $D(x+L/2)-D(L/2)=D(L/2)-D(x-L/2)$ originated from the assumption of a symmetric order parameter ${\rm{\Delta }}(L-s)={\rm{\Delta }}(s)$. As a result, we get a smooth function describing the solution within the interval $| x| \lt L/2$

$$\begin{eqnarray}&&g(x)=\displaystyle \sum _{\eta =\pm 1}{{\rm{e}}}^{-{\rm{i}}\varepsilon t+{\rm{i}}\eta {k}_{F}x}\left[{a}_{L}^{\eta }{{\rm{e}}}^{{\rm{i}}\eta {k}_{F}L/2}{{\rm{e}}}^{-D(x+L/2)/2}\left(\begin{array}{c}1\\ -{\rm{i}}\end{array}\right)\right.+\left.{{\rm{i}}{a}}_{R}^{\eta }{{\rm{e}}}^{-{\rm{i}}\eta {k}_{F}L/2}{{\rm{e}}}^{-D(x-L/2)/2}\left(\begin{array}{c}1\\ {\rm{i}}\end{array}\right)\right].\end{eqnarray} \tag{ A.8 }$$

At the ends of the wire we should put

$$\begin{eqnarray}&&g(-L/2)=\displaystyle \sum _{\eta =\pm 1}\left(\displaystyle \genfrac{}{}{0em}{}{{a}_{L}^{\eta }+{b}_{L}^{\eta }}{-{{\rm{i}}{a}}_{L}^{\eta }+{{\rm{i}}{b}}_{L}^{\eta }}\right){{\rm{e}}}^{-{\rm{i}}\varepsilon t+{\rm{i}}\eta 0},\end{eqnarray} \tag{ A.9 }$$

$$\begin{eqnarray}&&g(L/2)={\rm{i}}\displaystyle \sum _{\eta =\pm 1}\left(\displaystyle \genfrac{}{}{0em}{}{{a}_{R}^{\eta }+{b}_{R}^{\eta }}{{{\rm{i}}{a}}_{R}^{\eta }-{{\rm{i}}{b}}_{R}^{\eta }}\right){{\rm{e}}}^{-{\rm{i}}\varepsilon t+{\rm{i}}\eta 0},\end{eqnarray} \tag{ A.10 }$$

where we marked the left (right) movers by the exponents ${{\rm{e}}}^{\pm {\rm{i}}0}$. One can see that in the vicinity of the wire ends the wave function exhibits a 'jump' which occurs at the length scale of the coherence length ξ [24, 25].

A.2. Scattering problem

As a next step we use the scattering matrix approach to get the solution of a scattering problem for an electron plane wave ${\alpha }_{L(R)}{{\rm{e}}}^{\pm {{\rm{i}}{k}}_{F}x}$ incident from the left or right normal electrode. Note that it is enough to consider only incoming electrons, but not holes, if one integrates over the whole energy interval of the Fermi distribution to calculate the current. Moreover all the sources should be considered separately by putting only one of them to be non-zero at the same time and summing over all contributions in the observable to avoid any fake interference effects. Assuming the absence of the electron-hole conversion in the barriers and using the electron–hole symmetry in a superconductor one can separate complex conjugate electron and hole blocks in the total scattering matrix of the kth barrier ${\hat{S}}_{k}=\left(\genfrac{}{}{0em}{}{{s}_{k}\quad 0}{0\quad {s}_{k}^{* }}\right)$. We take a standard representation of the unitary matrix ${s}_{k}=\left(\genfrac{}{}{0em}{}{{R}_{k}\,{T}_{k}}{{T}_{k}\,-{R}_{k}^{* }{T}_{k}/{T}_{k}^{* }}\right)$ which transforms the incoming electron plane waves from the superconductor (e.g., ${a}_{L}^{-}+{b}_{L}^{-}$ for k = L) and from the normal reservoir (${\alpha }_{L}$) to the outgoing ones (${a}_{L}^{+}+{b}_{L}^{+}$ and ${u}_{L}=-{\alpha }_{L}{R}_{L}{T}_{L}/{T}_{L}^{* }+{T}_{L}({a}_{L}^{-}+{b}_{L}^{-})$ for k = L) at both interfaces (see figure A1 for all notations). Here ${R}_{L}={r}_{L}{{\rm{e}}}^{{\rm{i}}{\phi }_{L}}$, ${R}_{R}={r}_{R}{{\rm{e}}}^{-{\rm{i}}{\phi }_{R}}$ and T_k are the reflection and transmission matrix coefficients, ${r}_{k}=\sqrt{1-| {T}_{k}{| }^{2}}$ and ${\phi }_{k}$ are reflection amplitude and phase.

Zoom In Zoom Out Reset image size

The scattering matrices impose the following boundary conditions on the plane wave amplitudes

$$\begin{eqnarray}&&{T}_{L}{\alpha }_{L}+{R}_{L}({a}_{L}^{-}+{b}_{L}^{-})={a}_{L}^{+}+{b}_{L}^{+},\end{eqnarray} \tag{ A.11 }$$

$$\begin{eqnarray}&&{R}_{L}^{* }({a}_{L}^{+}-{b}_{L}^{+})={a}_{L}^{-}-{b}_{L}^{-},\end{eqnarray} \tag{ A.12 }$$

$$\begin{eqnarray}&&{T}_{R}{\alpha }_{R}+{R}_{R}({a}_{R}^{+}+{b}_{R}^{+})={a}_{R}^{-}+{b}_{R}^{-},\end{eqnarray} \tag{ A.13 }$$

$$\begin{eqnarray}&&{R}_{R}^{* }({a}_{R}^{-}-{b}_{R}^{-})={a}_{R}^{+}-{b}_{R}^{+}.\end{eqnarray} \tag{ A.14 }$$

Substituting equations (A.6), (A.7) and introducing the notations ${a}_{k}^{\pm }={{\rm{e}}}^{\pm {\rm{i}}{\phi }_{k}/2}({A}_{k}\pm {a}_{k})/2$ one can obtain the following set of equations

$$\begin{eqnarray}\left(\begin{array}{cccc}{{\rm{\Gamma }}}_{L}-{\rm{i}}\varepsilon & {\omega }_{0} & 0 & {\rm{i}}{\tilde{\omega }}_{0}\\ -{\omega }_{0} & {{\rm{\Gamma }}}_{R}-{\rm{i}}\varepsilon & {\rm{i}}{\tilde{\omega }}_{0} & 0\\ 0 & {\rm{i}}{\tilde{\omega }}_{0} & \tfrac{{\tilde{{\rm{\Delta }}}}^{2}}{{{\rm{\Gamma }}}_{L}}-{\rm{i}}\varepsilon & {\omega }_{0}\\ {\rm{i}}{\tilde{\omega }}_{0} & 0 & -{\omega }_{0} & \tfrac{{\tilde{{\rm{\Delta }}}}^{2}}{{{\rm{\Gamma }}}_{R}}-{\rm{i}}\varepsilon \end{array}\right)\left(\begin{array}{c}{A}_{L}\\ {A}_{R}\\ {a}_{L}\\ {a}_{R}\end{array}\right)=\tilde{{\rm{\Delta }}}\left(\begin{array}{c}{\rho }_{L}{\alpha }_{L}\\ {\rho }_{R}{\alpha }_{R}\\ \tfrac{1}{{\rho }_{L}}{\alpha }_{L}\\ \tfrac{1}{{\rho }_{R}}{\alpha }_{R}\end{array}\right),\end{eqnarray} \tag{ A.15 }$$

where ${{\rm{\Gamma }}}_{k}=\tilde{{\rm{\Delta }}}{\rho }_{k}^{2}$, ${\rho }_{k}^{2}=(1-{r}_{k})/(1+{r}_{k})$, ${\omega }_{0}=\tilde{{\rm{\Delta }}}{{\rm{e}}}^{-D(L/2)}\sin \varphi$, ${\tilde{\omega }}_{0}=\tilde{{\rm{\Delta }}}{{\rm{e}}}^{-D(L/2)}\cos \varphi$, $\varphi ={k}_{F}L+({\phi }_{L}\,-{\phi }_{R})/2$. The phase ${\chi }_{k}$ of the transmission coefficients ${T}_{k}=| {T}_{k}| {{\rm{e}}}^{{\rm{i}}{\chi }_{k}}$ does not affect any measurable quantity, therefore we choose it equal to ${\chi }_{k}={\phi }_{k}/2$ for the sake of simplicity.

A standard recipe to describe the low-frequency (ω) dynamics is to replace the energy $\varepsilon$ by the time derivative ${\rm{i}}\partial /\partial t$. In the isolated wire, ${\rho }_{k}\to 0$, the fast decaying modes ${a}_{k}\sim {\alpha }_{k}{\rho }_{k}$ disappear as they correspond to the states of the continuous spectrum in the wire and do not satisfy the boundary conditions. Resulting equations in the closed wire give two energy levels $\varepsilon =\pm {\omega }_{0}$ and correspond to the beating between A_L and A_R in the time domain (see equations (A.16), (A.17) below). Assuming naturally that ${{\rm{\Gamma }}}_{k}\omega \ll {\tilde{{\rm{\Delta }}}}^{2}$ one can find that fast decaying modes ${a}_{k}\approx {\alpha }_{k}{\rho }_{k}-{\rm{i}}{\tilde{\omega }}_{0}{\rho }_{k}^{2}{A}_{k^{\prime} }\mp {\omega }_{0}{\rho }_{k}^{2}{\rho }_{k^{\prime} }{\alpha }_{k^{\prime} }$ corresponding to the continuous spectrum contributions give small corrections in ${{\rm{e}}}^{-L/\xi }$ to the equations for the low-decaying ones A_k. Here and further $k^{\prime} =R(L)$ for $k=L(R)$. Indeed, this leads to a relative renormalization of the decay rates ${{\rm{\Gamma }}}_{k}$ and sources ${\rho }_{k}$ by a small values $\sim {\omega }_{0}^{2}/{\tilde{{\rm{\Delta }}}}^{2}$ and to the addition of the ${\alpha }_{R(L)}$ source proportional to ${\tilde{\omega }}_{0}/\tilde{{\rm{\Delta }}}\sim {{\rm{e}}}^{-L/\xi }$ to the equation for ${A}_{L(R)}$. All these terms corresponds to a direct tunneling of electron(s) from the lead to the opposite end of the wire. Further we neglect these contributions taking into account only a local tunneling from the kth leads to the kth end of the wire and considering therefore only first two equations for the amplitudes ${A}_{L,R}$ of Majorana states in (A.15) without a_k.

Transforming the equations to the Schrödinger representation one can obtain equations (9), (10) from the main text

$$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{\partial t}+{{\rm{\Gamma }}}_{L}\right){A}_{L}+{\omega }_{0}{A}_{R}={F}_{L}{{\rm{e}}}^{-{\rm{i}}\varepsilon t},\end{eqnarray} \tag{ A.16 }$$

$$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{\partial t}+{{\rm{\Gamma }}}_{R}\right){A}_{R}-{\omega }_{0}{A}_{L}={F}_{R}{{\rm{e}}}^{-{\rm{i}}\varepsilon t},\end{eqnarray} \tag{ A.17 }$$

with the choice of sources ${F}_{k}=\tilde{{\rm{\Delta }}}{\rho }_{k}{\alpha }_{k}$ appropriate to the replacement $\varepsilon \to {\rm{i}}\partial /\partial t$. Beyond the stationary regime one can consider the parameters ${\omega }_{0}$, ${{\rm{\Gamma }}}_{k}$ and ${\rho }_{k}$ to be time-dependent keeping the equations (A.16), (A.17) intact for the typical frequency ω of the drive small compared to the gap $\omega \ll \tilde{{\rm{\Delta }}}$.

Note that the equations of motion for Majorana amplitudes A_k in the Heisenberg representation (see the first two lines in equation (A.15)) correspond to the scattering matrix through a scatterer with an internal structure described in [36] and applied for the p-wave superconducting wire, e.g., in the [16, 27].

In this section we consider for simplicity only the case of the non-zero left source ${\alpha }_{L}$, since the results for the right source can be derived using the symmetry $L\leftrightarrow R$. According to [28] the energy resolved contribution to the differential conductance g_k of the kth interface can be written as a sum of the quasiparticle flows of the left and right moving electrons and holes with the corresponding signs

$$\begin{eqnarray}&&{g}_{k}(\varepsilon )=1-{R}_{k}^{e}+{R}_{k}^{h}=| {a}_{k}^{+}+{b}_{k}^{+}{| }^{2}+| {a}_{k}^{+}-{b}_{k}^{+}{| }^{2}-| {a}_{k}^{-}+{b}_{k}^{-}{| }^{2}-| {a}_{k}^{-}-{b}_{k}^{-}{| }^{2}.\end{eqnarray} \tag{ B.1 }$$

We used here the conservation of the quasiparticle flow at the interface which results in the unitarity of the scattering matrix. Substituting the expressions for ${b}_{k}^{\pm }$ (A.6), (A.7) and for ${a}_{k}^{\pm }={{\rm{e}}}^{\pm {\rm{i}}{\phi }_{k}/2}({A}_{k}\pm {a}_{k})/2$ through the amplitudes A_k and a_k into the equation (B.1) one can obtain

$$\begin{eqnarray}{g}_{k}(\varepsilon ) & = & \displaystyle \frac{2}{{\tilde{{\rm{\Delta }}}}^{2}}\mathrm{Re}[{A}_{k}{a}_{k}^{* }({\tilde{{\rm{\Delta }}}}^{2}+{\varepsilon }^{2})-{A}_{k}{A}_{k^{\prime} }^{* }\varepsilon {\tilde{\omega }}_{0}-{{\rm{i}}{A}}_{k^{\prime} }{a}_{k}^{* }\varepsilon {\omega }_{0}\\ & & -{a}_{k^{\prime} }{a}_{k}^{* }\varepsilon {\tilde{\omega }}_{0}-{{\rm{i}}{A}}_{k}{a}_{k^{\prime} }^{* }\varepsilon {\omega }_{0}+{A}_{k^{\prime} }{a}_{k^{\prime} }^{* }({\omega }_{0}^{2}+{\tilde{\omega }}_{0}^{2})].\end{eqnarray} \tag{ B.2 }$$

Omitting the terms which are small in the parameters ${\omega }_{0}/\tilde{{\rm{\Delta }}},{\tilde{\omega }}_{0}/\tilde{{\rm{\Delta }}}\sim {{\rm{e}}}^{-L/\xi }$ (see the previous section) one can keep only the first term in the latter equation

$$\begin{eqnarray}&&{g}_{k}(\varepsilon )=2\mathrm{Re}[{A}_{k}{a}_{k}^{* }]+O\left({A}_{k}^{2}\displaystyle \frac{\varepsilon }{\tilde{{\rm{\Delta }}}}{{\rm{e}}}^{-L/\xi }\right).\end{eqnarray} \tag{ B.3 }$$

In the stationary regime one can express both A_k and a_k from equations (A.15)

$$\begin{eqnarray}&&{A}_{k}=\tilde{{\rm{\Delta }}}\displaystyle \frac{{\rho }_{k}{\alpha }_{k}({{\rm{\Gamma }}}_{k^{\prime} }-{\rm{i}}\varepsilon )\mp \omega {\rho }_{k^{\prime} }{\alpha }_{k^{\prime} }}{({{\rm{\Gamma }}}_{L}-{\rm{i}}\varepsilon )({{\rm{\Gamma }}}_{R}-{\rm{i}}\varepsilon )+{\omega }_{0}^{2}},\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&{a}_{k}\approx {\alpha }_{k}{\rho }_{k}\end{eqnarray} \tag{ B.5 }$$

and show that equation (B.3) transforms into equation (C.1) of the next section, due to incoherence of the left and right sources (${\alpha }_{L}{\alpha }_{R}^{* }\to 0$). Note that we neglect here all the direct tunneling processes in the wire which give the exponentially small corrections to the equation (B.3) in the parameter $L/\xi$.

In general to calculate the zero temperature differential conductance g_k(eV) of the kth interface of the wire in the system with time-dependent parameters one should solve equations (A.16), (A.17) and substitute the solutions A_k into the equation (B.3) together with the expression (B.5) for a_k.

Here we consider the dc transport for a constant applied bias $V={V}_{L}-{V}_{R}$ using the formalism of the previous section and putting ${A}_{L,R}\propto {{\rm{e}}}^{-{\rm{i}}\varepsilon t}$. As a result we obtain

$$\begin{eqnarray}&&{g}_{L,R}(\varepsilon )=\displaystyle \frac{2{{\rm{\Gamma }}}_{L,R}({{\rm{\Gamma }}}_{R,L}{\omega }_{0}^{2}+{{\rm{\Gamma }}}_{L,R}({\varepsilon }^{2}+{{\rm{\Gamma }}}_{R,L}^{2}))}{{({\varepsilon }^{2}-{\omega }_{0}^{2}-{{\rm{\Gamma }}}_{L,R}{{\rm{\Gamma }}}_{R,L})}^{2}+{\varepsilon }^{2}{({{\rm{\Gamma }}}_{L,R}+{{\rm{\Gamma }}}_{R,L})}^{2}}.\end{eqnarray} \tag{ C.1 }$$

It is convenient to discuss separately the limits $R\to \infty$ and $R\to 0$. For the first limit the zero-temperature differential conductance of the device in the symmetric case ${{\rm{\Gamma }}}_{L}={{\rm{\Gamma }}}_{R}$ and ${V}_{R}=-{V}_{L}$ takes the form ${\rm{d}}{I}/{\rm{d}}{V}={e}^{2}{g}_{L}({eV}/2)/2\pi$ with $I={I}_{L}=-{I}_{R}$ and ${g}_{L}(\varepsilon )={g}_{R}(\varepsilon )$. In the opposite limit of $R\to 0$ similar formulas for the differential conductances hold for each interface separately, i.e., ${{\rm{d}}{I}}_{L,R}/{{\rm{d}}{V}}_{L,R}={e}^{2}{g}_{L,R}({{eV}}_{L,R})/\pi$. Thus, in both limits we obtain the conductance peak near the zero bias at ${eV}\sim {\omega }_{0}$. It is this peak which is usually considered [8, 9] as an experimental evidence for the Majorana states in semiconducting wires with the induced superconducting order. The nonlocal nature of the Majorana pair reveals itself in the zero bias dip which is of course smeared due to finite rates ${{\rm{\Gamma }}}_{L,R}$ of tunneling to the fermionic baths. As a result, for the exponentially small splitting ${\omega }_{0}$ the dip completely disappears for ${{\rm{\Gamma }}}_{L}\sim {{\rm{\Gamma }}}_{R}\sim {\omega }_{0}$ and can survive only in a rather exotic limit of strong asymmetry between the couplings to the left and right reservoirs. The latter situation can be realized, in particular, for the dip in the curve ${{\rm{d}}{I}}_{L}/{{\rm{d}}{V}}_{L}$ for $R\to 0$ and ${{\rm{\Gamma }}}_{R}=0$ [30–32]. A more realistic case with both nonzero tunneling rates and the peak broadening due to the finite temperature and inelastic effects can make the experimental observation of the ${\omega }_{0}$ scale in dc transport difficult.