Distinguishing Majorana zero modes from impurity states through time-resolved transport

Topological quantum computing [1] is an active field of research based on the key idea to reduce quantum decoherence issues by using topologically protected states [2, 3]. Majorana fermions are their own antiparticles [4], and their condensed-matter analogs, Majorana bound state or Majorana zero mode (MZM), retain this feature [5]. They are thus considered to be promising candidates for technological advances in topological quantum computing [6–8] since their non-abelian statistics allow performing quantum computation protected from environmental perturbations [9]. Even though various experimental signatures of MZM have been reported [10–17], a clear and unambiguous detection and the consequent control of these states has proven difficult so far. For example, other types of bound states [18], surface states 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 [19]. 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 [20–22]. This is especially relevant when nowaday transport measurements are pushing the temporal resolution to sub-picosecond regime [23–30], 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 [31–36] 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 IS, we discover that the dynamics for (i) and (ii) look 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. For case (ii) the corresponding oscillations are suppressed. In addition, we study the transient response of quasi-Majorana states (QMS) in the topologically trivial phase [37–44], 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 could therefore be used to identify the MZM or QMS.

We consider a normal metal–superconductor–normal metal (NSN) junction, see figure 1. The superconducting central region of the junction is a nanowire in proximity to an s-wave bulk SC with order parameter Δ. The nanowire in addition features a strong spin–orbit interaction (e.g. InSb [45, 46] or InAs [13, 47]) which favors aligning the spins along the ±y direction. An external magnetic field parallel to the nanowire breaks time-reversal symmetry and aligns the spins along the ±z direction, introducing a Zeeman splitting V_Z = gμ_BB/2 where g is the Landé factor and μ_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 [5, 10, 48–51]. Specifically the infinite nanowire is in a topologically nontrivial phase for V_Z > Δ > 0 and certain intervals of the chemical potential [48, 49, 51], 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.

Zoom In Zoom Out Reset image size

We write the total Hamiltonian as

$$\begin{eqnarray}&&\hat{H}={\hat{H}}_{e}+{\hat{H}}_{c}+{\hat{H}}_{w},\end{eqnarray} \tag{ 1 }$$

where the individual components for the electrodes and coupling are characterized by the single-particle energy dispersion in the electrodes _kλ and by the coupling matrix elements T_jkλ between the states in the nanowire and the electrodes [33]:

$$\begin{eqnarray}&&{\hat{H}}_{e}=\displaystyle \sum _{k\lambda }{\epsilon }_{k\lambda }{\hat{c}}_{k\lambda }^{\dagger }{\hat{c}}_{k\lambda }\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{\hat{H}}_{c}=\displaystyle \sum _{{jk}\lambda }({T}_{{jk}\lambda }{\hat{c}}_{j}^{\dagger }{\hat{c}}_{k\lambda }+{\rm{h}}.{\rm{c}}.).\end{eqnarray} \tag{ 3 }$$

Here kλ labels the kth basis element in the λth electrode, and j labels the atomic sites on the nanowire. The nanowire, in turn, is characterized by [51, 52]

$$\begin{eqnarray}&&{\hat{H}}_{w}=\displaystyle \sum _{j}\left[-\displaystyle \frac{J}{2}({\hat{c}}_{j}^{\dagger }{\hat{c}}_{j+1}+{\rm{h}}.{\rm{c}}.)-(\mu -J){\hat{c}}_{j}^{\dagger }{\hat{c}}_{j}-\displaystyle \frac{\alpha }{2}({\rm{i}}{\hat{c}}_{j}^{\dagger }{\sigma }_{2}{\hat{c}}_{j+1}+{\rm{h}}.{\rm{c}}.)+{V}_{{\rm{Z}}}{\hat{c}}_{j}^{\dagger }{\sigma }_{3}{\hat{c}}_{j}+{\rm{\Delta }}({\hat{c}}_{j\uparrow }{\hat{c}}_{j\downarrow }+{\rm{h}}.{\rm{c}}.)\right],\end{eqnarray} \tag{ 4 }$$

where J, μ, α, V_Z, and Δ 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 σ_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 (λ ∈ {S, D}, see figure 1) this out-of-equilibrium condition is defined by the source–drain voltage V_SD = V_S − V_D. The transport setup is considered partition-free [53–55] meaning that the whole system is initially contacted in a global thermo-chemical equilibrium at unique chemical potential μ and at inverse temperature $\beta \equiv {\left({k}_{{\rm{B}}}T\right)}^{-1}$.

For a compact notation we introduce Nambu spinors [56–58] ${\underline{\hat{{\rm{\Phi }}}}}_{x}\equiv {\left({\hat{{\rm{\Phi }}}}_{x}^{1},{\hat{{\rm{\Phi }}}}_{x}^{2},{\hat{{\rm{\Phi }}}}_{x}^{3},{\hat{{\rm{\Phi }}}}_{x}^{4}\right)}^{T}\equiv {\left({\hat{c}}_{x\uparrow },{\hat{c}}_{x\downarrow }^{\dagger },{\hat{c}}_{x\downarrow },{\hat{c}}_{x\uparrow }^{\dagger }\right)}^{T}$, and the anticommutation relation is then understood componentwise $\{{\hat{{\rm{\Phi }}}}_{x}^{\mu },{\left({\hat{{\rm{\Phi }}}}_{y}^{\nu }\right)}^{\dagger }\}={\delta }_{{xy}}{\delta }^{\mu \nu }$. Here we denote quantities in the Nambu ⨂ spin space by an underline. This representation allows for writing the Hamiltonian for the nanowire in a Bogoliubov–de Gennes form [59, 60]

$$\begin{eqnarray}&&{\hat{H}}_{w}=\displaystyle \frac{1}{2}\displaystyle \sum _{j}\left[{\underline{\hat{{\rm{\Phi }}}}}_{j}^{\dagger }{\underline{a}}_{j}{\underline{\hat{{\rm{\Phi }}}}}_{j}+({\underline{\hat{{\rm{\Phi }}}}}_{j}^{\dagger }{\underline{b}}_{j}{\underline{\hat{{\rm{\Phi }}}}}_{j+1}+{\rm{h}}.{\rm{c}}.)\right],\end{eqnarray} \tag{ 5 }$$

where we introduced on-site and nearest-neighbor contributions [52]

$$\begin{eqnarray}{\underline{a}}_{j}={\left(\begin{array}{cccc}J-\mu +{V}_{{\rm{Z}}} & -{\rm{\Delta }} & 0 & 0\\ -{\rm{\Delta }} & \mu -J+{V}_{{\rm{Z}}} & 0 & 0\\ 0 & 0 & J-\mu -{V}_{{\rm{Z}}} & {\rm{\Delta }}\\ 0 & 0 & {\rm{\Delta }} & \mu -J-{V}_{{\rm{Z}}}\end{array}\right)}_{j},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}{\underline{b}}_{j}={\left(\begin{array}{cccc}-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{array}\right)}_{j},\end{eqnarray} \tag{ 7 }$$

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

$$\begin{eqnarray}&&{\hat{H}}_{e}=\displaystyle \frac{1}{2}\displaystyle \sum _{k\lambda }{\underline{\hat{{\rm{\Phi }}}}}_{k\lambda }^{\dagger }{\underline{\epsilon }}_{k\lambda }{\underline{\hat{{\rm{\Phi }}}}}_{k\lambda }\end{eqnarray} \tag{ 8 }$$

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

$$\begin{eqnarray}&&{\hat{H}}_{c}=\displaystyle \frac{1}{2}\displaystyle \sum _{{jk}\lambda }({\underline{\hat{{\rm{\Phi }}}}}_{j}^{\dagger }{\underline{T}}_{{jk}\lambda }{\underline{\hat{{\rm{\Phi }}}}}_{k\lambda }+{\rm{h}}.{\rm{c}}.)\end{eqnarray} \tag{ 9 }$$

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

By using the nonequilibrium Green's function approach [32] 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 [57]

$$\begin{eqnarray}&&{\underline{G}}_{{xy}}(z,z^{\prime} )=-{\rm{i}}\langle {{ \mathcal T }}_{\gamma }[{\underline{\hat{{\rm{\Phi }}}}}_{x}(z)\otimes {\underline{\hat{{\rm{\Phi }}}}}_{y}^{\dagger }(z^{\prime} )]\rangle ,\end{eqnarray} \tag{ 10 }$$

where the contour-ordering operator ${{ \mathcal T }}_{\gamma }$ is taken for the variables $z,z^{\prime}$ on the Keldysh contour γ [32]. The form in equation (10) automatically handles both normal and anomalous components of the Green's function [61]. In appendix A we show that the equations of motion for the Green's function are exactly the same as those in [33, 35], 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 -{\rm{i}}{\underline{G}}^{\lt }(t,t)$ from the lesser Green's function (see appendices B–D). In order to obtain a closed solution to the equation of motion we have described the electrodes within wide-band approximation (WBA), 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 Γ_λ.

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 [62–64]. In appendix E we derive the following expression for the bond current between the sites j and j + 1 within the nanowire:

$$\begin{eqnarray}&&{I}_{j,j+1}=2\mathrm{Im}\left[-\displaystyle \frac{J}{2}(\langle {\hat{c}}_{j\uparrow }^{\dagger }{\hat{c}}_{(j+1)\uparrow }\rangle +\langle {\hat{c}}_{j\downarrow }^{\dagger }{\hat{c}}_{(j+1)\downarrow }\rangle )-\displaystyle \frac{\alpha }{2}\left(\langle {\hat{c}}_{j\uparrow }^{\dagger }{\hat{c}}_{(j+1)\downarrow }\rangle -\langle {\hat{c}}_{j\downarrow }^{\dagger }{\hat{c}}_{(j+1)\uparrow }\rangle \right)+2{\rm{\Delta }}\displaystyle \sum _{m=1}^{j}\langle {\hat{c}}_{m\downarrow }{\hat{c}}_{m\uparrow }\rangle \right],\end{eqnarray} \tag{ 11 }$$

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

3.1. Emergence of the MZM

Using equation (11) we calculate the (steady-state) current–voltage characteristics for nanowires of varying lengths. The parameter space of the model in equation (4) is rather broad. Instead of mapping out the individual details of each parameter, we concentrate our discussion on representative points in the parameter (phase) space. Namely, for the nanowire we choose J = 1, α = 0.5, V_Z = 0.25, Δ = 0.1, and μ = 0 [52] using which the nanowire possesses the topological Majorana state for long enough nanowires. We have investigated many parameters in this topological regime and found that they do not result in qualitatively different behavior. This choice also 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 Γ_λ = 0.01. The bias voltage is applied symmetrically for the source and drain electrodes V_S = −V_D ≡ V, and we consider the zero-temperature limit.

In figure 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. The inset of figure 2 shows the exponential localization of the MZM for N_w = 60. These two zero-energy states are far apart and the coupling between them is weak, and we see a single zero-bias peak whose value is exactly one conductance quantum. The coupling between these two zero-energy states is enhanced for shorter nanowires leading to the splitting of the peak.

Zoom In Zoom Out Reset image size

3.2. Transient signature of the MZM

We evaluate transient currents through a N_w = 50 nanowire by considering the two centermost sites in equation (11). In addition to the topological SC with the MZM, we consider an ordinary SC wire (same as figure 2 but for V_Z = 0) and three auxiliary cases of magnetic impurities deposited at the edges of an ordinary SC. In accordance with the Anderson's theorem, non-magnetic impurities cannot reduce the superconducting gap Δ for an ordinary (time-reversal invariant) superconductor, and therefore there cannot be any in-gap states from non-magnetic impurities [65–67]. With this in mind, we introduce an interaction potential containing both magnetic (${ \mathcal M }$) and non-magnetic (${ \mathcal N }$) parts [68]

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{imp}}=\displaystyle \sum _{j}{ \mathcal M }({\hat{c}}_{j\uparrow }^{\dagger }{\hat{c}}_{j\uparrow }-{\hat{c}}_{j\downarrow }^{\dagger }{\hat{c}}_{j\downarrow })+{ \mathcal N }({\hat{c}}_{j\uparrow }^{\dagger }{\hat{c}}_{j\uparrow }+{\hat{c}}_{j\downarrow }^{\dagger }{\hat{c}}_{j\downarrow }),\end{eqnarray} \tag{ 12 }$$

where the index j runs over the host sites in the ordinary superconductor where the impurities are attached to. The parameters ${ \mathcal M }$ and ${ \mathcal N }$ for these IS can be modeled by modified tight-binding parameters [69, 70] for the terminal sites in the nanowire, j = {1, N_w} in equations (6) and (7), $\widetilde{\mu }$ and ${\widetilde{V}}_{{\rm{Z}}}$ (modified parameters signified by a tilde). For the ordinary SC wire supplemented with impurities we also concentrate our discussion on representative points in the nontopological parameter (phase) space. We keep the time-reversal invariant bulk superconductor unchanged by setting $\widetilde{{\rm{\Delta }}}={\rm{\Delta }}=0.1$, although it has also been reported that impurities may suppress or even destroy superconductivity locally [71–74]. For the ordinary SC we have V_Z = 0 but in order to introduce magnetic impurities we set ${\widetilde{V}}_{{\rm{Z}}}\ne 0$, and the energy of these in-gap states may be tuned close to zero by varying ${\widetilde{V}}_{{\rm{Z}}}$. For our purposes the exact formulation of the magnetic impurities is not too important as long as there is a separate state within the gap with different topological character compared to the MZM. In addition to the on-site scattering potential in equation (12), the hybridization of the impurities with the underlying superconductor $\widetilde{J}$ (and $\widetilde{\alpha }$)may be different than the bulk hopping J (and α), resulting in an effective tunnel barrier around the impurity [66, 69, 70, 75, 76].

The different cases considered in our comparative simulations are collected in table 1. We single out the in-gap states 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 in-gap state (zero mode) in the nanowire.

Table 1.  Model parameters for the superconducting nanowire. Modified tight-binding parameters for the terminal sites of the nanowire are signified by tildes.

Case	μ	$\widetilde{\mu }$	J	$\widetilde{J}$	α	$\widetilde{\alpha }$	${V}_{{\rm{Z}}}$	${\widetilde{V}}_{{\rm{Z}}}$	${\rm{\Delta }}$	$\widetilde{{\rm{\Delta }}}$
(1) Ordinary SC	0	0	1	1	0.5	0.5	0	0	0.1	0.1
(2) Topological SC	0	0	1	1	0.5	0.5	0.25	0.25	0.1	0.1
(3) Ordinary SC + IS1	0	0	1	0.2	0.5	0.1	0	0.21	0.1	0.1
(4) Ordinary SC + IS2	0	0	1	0.5	0.5	0.25	0	0.43	0.1	0.1
(5) Ordinary SC + IS3	0	0	1	1	0.5	0.5	0	0.67	0.1	0.1

The steady-state dI/dV signals of cases (2)–(5) look qualitatively similar showing a zero-bias peak (figure 3(a)) as we tuned the magnetic impurity with the parameter for ${\widetilde{V}}_{{\rm{Z}}}$ to appear close to zero energy, like the MZM. The resonances corresponding to the magnetic impurities around zero energy are broadened compared to the MZM for $\widetilde{J}\lt J$ due to an increased effective tunnel barrier [77–81]. The broadening of the MZM resonance is therefore systematically narrower than the one of IS. The MZM emerges completely from the physical mechanism described by the model in equation (4), and the width of the MZM peak is therefore completely specified by the coupling to the leads. While the model parameters of the nanowire are locally modified for the description of the IS, the coupling Hamiltonian in equation (3), importantly, remains independent of the modified central region. In particular, this means that the Hamiltonian of the nanowire and the corresponding eigenstates result in the observed behavior for the same tunneling rate Γ. Due to this broadening, the transient signals in figure 3(b) for these cases can be qualitatively different. First, for the ordinary SC without the MZM, case (1), the current signal is zero on average due to there being no transport channels within the SC gap and the small bias window. Second, for the topological SC with the MZM, case (2), the transient oscillations last for thousands of units of inverse hopping (for hopping energies in the eV scale we have J⁻¹ ∼ 0.658 fs), i.e. up to picoseconds. The MZM at the edges of the nanowire are 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 weak leading to a very narrow resonance and a long lifetime. This generic finding is not limited to our parameter choice for the MZM. Third, for the ordinary SC with magnetic impurities, cases (3)–(5) the transient oscillations may be suppressed (compared to the MZM) due to the broader resonance. For case (3) the transient current rises rapidly but also saturates relatively fast to its stationary value within couple of hundred time units. The IS is directly connected to the electrodes resulting in a strong hybridization and in a relatively fast decay of the transient. When the effective tunnel barrier due to the magnetic impurity is decreased, cases (4) and (5), the transient time scales approach case (2) although they do not exceed it. The decay rate can be approximated by the expectation value of the tunneling rate operator: $\gamma ={\sum }_{j=1}^{2}\langle {\varphi }_{j}| {\rm{\Gamma }}| {\varphi }_{j}\rangle$, where ${\rm{\Gamma }}\equiv {\sum }_{\lambda }{{\rm{\Gamma }}}_{\lambda }$ and $| \varphi \rangle$ are the MZM or IS eigenvectors, see the dashed lines in figure 3(b). For identical wire-electrode coupling, the decay time 1/γ of the MZM transient current is between 1.5 and 5 times the one of the IS.

Zoom In Zoom Out Reset image size

Additionally, the transient oscillations for the MZM, case (2), are more pronounced than the ones for the IS, cases (3)–(5). This difference can be seen by taking the Fourier transforms of the time-dependent signals, see figure 3(c). The Fourier transforms are calculated from an extended temporal window (up to t = 10000J⁻¹) for better frequency resolution, we subtract the steady-state value from the sample points in order to get rid of a divergence at zero frequency, and ultimately we take the absolute value of the result. 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 for the IS are suppressed due to the broader resonances and shorter lifetimes. Before entering the band of all possible transitions outside the SC gap ($\omega \geqslant 2{\rm{\Delta }}=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. The resonances due to the IS are tuned to be close to the zero energy, but they are still strictly speaking nonzero compared to the actual MZM. This is seen as an additional Fourier peak forming around ω = 0 for cases (4) and (5). This could be understood as an artificial intra-level transition between the magnetic IS. Importantly, this peak is not seen for the MZM. Overall the transient features of the MZM can be distinguished from the IS.

3.3. Comparison between MZMs and QMS

Even though we found a 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 ${\rm{\Delta }}\ll {V}_{{\rm{Z}}}\lt \mu$ the resulting in-gap QMS emerge without additional surface or IS but by adding a smooth confining potential [37–40, 42]. These states can also be tuned arbitrarily close to zero energy, thereby mimicking the behavior of the MZM.

We implement a confining potential within the nanowire as a simple function of the lattice coordinate j ∈ [0, N_w) labeling the atomic sites on the nanowire:

$$\begin{eqnarray}{f}_{s}(j)=\left\{\begin{array}{cl}{\sin }^{2}\left[\displaystyle \frac{\pi {sj}}{2({N}_{w}-1)}\right], & j\lt \displaystyle \frac{{N}_{w}}{s},\\ 1, & \displaystyle \frac{{N}_{w}}{s}\leqslant j\lt \displaystyle \frac{(s-1)({N}_{w}-1)}{s},\\ {\cos }^{2}\left[\displaystyle \frac{\pi {sj}}{2({N}_{w}-1)}\right], & j\geqslant \displaystyle \frac{(s-1)({N}_{w}-1)}{s},\end{array}\right.\end{eqnarray} \tag{ 13 }$$

where s controls the smoothness at the edges. We then re-cast the values of the spin–orbit interaction and pair potential accordingly: α → α f_s(j) and Δ → Δf_s(j). Large values of s correspond to an abrupt hard-wall confinement where both α and Δ 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 α = 0.5, V_Z = 1.2, Δ = 0.1, and μ = 2.0 both with abrupt ($s\to \infty$) and smooth (s = {17, 9}) confinement potentials according to equation (13). In figure 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 figures 4(b) and (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. It is also possible that more than one Majorana or quasi-Majorana pair coexist leading to multiple in-gap resonances of both zero and nonzero energies [82, 83]. In these situations we would expect a rich transient signature with oscillations between different Majorana (and quasi-Majorana) states reflecting all the intricacies of the in-gap level structure.

Zoom In Zoom Out Reset image size

We studied the time-dependent features of MZMs and QMS in a superconducting nanowire in contrast with magnetic IS. The transient features related to MZM and QMS were found to be different than the ones resulting from magnetic impurities: 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. Compared to the MZM and QMS, the broadening of the resonances related to the IS and the consequent faster decay times could be attributed to the impurity-induced hopping disorder and the consequent increased effective tunnel barrier around the impurities. This finding could be utilized in possible detection and identification of the MZM or QMS via ultrafast transport measurements [23–30, 84]. In order to estimate a limit for the time scales that one needs to be able to resolve, we refer to the experiment of Mourik et al [10], who found an induced superconducting energy gap of 250 μeV. Using this as an upper limit for the energy scale of the in-gap oscillations, the fastest temporal oscillation period related to such processes would be 16.5 picoseconds (frequency of 60 GHz). While this might be at limits of what is routinely measurable, recent ultrafast transport measurements showed a sub-picosecond time resolution [30] and should definitely be able to resolve the oscillations predicted by us.

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 [40]. Since topological properties of the MZM should be robust against electronic interactions [51], 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 [85–87]. Together with ultrafast optical switching of chiral superconductors [88, 89] or nonequilibrium engineering of topologically nontrivial states of matter [90–98] our findings highlight the great potential of ultrafast techniques for advances towards topological quantum computation.

RT and MAS acknowledge funding by the DFG (Grant No. SE 2558/2-1) through the Emmy Noether program. RVL acknowledges funding from the Academy of Finland (Project No. 317139). EP and GS acknowledge EC funding through the RISE Co-ExAN (Grant No. GA644076). EP 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. GS also acknowledges Tor Vergata University for financial support through the Mission Sustainability Project 2DUTOPI.

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 λ 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

$$\begin{eqnarray}\underline{{\boldsymbol{h}}}=\left(\begin{array}{cccc}{\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{array}\right)\end{eqnarray} \tag{ A.1 }$$

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 (equations (6) and (7) in the main text), 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 equation (A.1) does not refer to the Nambu⨂spin space but it is of dimension (N_e + 1) × (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 equation (10) in the main text (indices x, y belonging either to the electrodes or to the central region) therefore label the transport setup in the same block form

$$\begin{eqnarray}\underline{{\boldsymbol{G}}}=\left(\begin{array}{cccc}{\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{array}\right).\end{eqnarray} \tag{ A.2 }$$

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

$$\begin{eqnarray}&&{\rm{i}}{\partial }_{z}{\underline{G}}_{{xy}}(z,z^{\prime} )={\partial }_{z}\left[\theta (z,z^{\prime} )\langle {\underline{\hat{{\rm{\Phi }}}}}_{x}(z)\otimes {\underline{\hat{{\rm{\Phi }}}}}_{y}^{\dagger }(z^{\prime} )\rangle -\theta (z^{\prime} ,z)\langle {\underline{\hat{{\rm{\Phi }}}}}_{y}{\left(z^{\prime} \right)}^{\dagger }\otimes {\underline{\hat{{\rm{\Phi }}}}}_{x}(z)\rangle \right],\end{eqnarray} \tag{ A.3 }$$

where the step function is defined on the Keldysh contour γ according to the contour-ordering operator ${{ \mathcal T }}_{\gamma }$ [32]. Evaluating the derivative gives

$$\begin{eqnarray}&&{\rm{i}}{\partial }_{z}{\underline{G}}_{{xy}}(z,z^{\prime} )=\delta (z,z^{\prime} )\left\{{\underline{\hat{{\rm{\Phi }}}}}_{x}(z),{\underline{\hat{{\rm{\Phi }}}}}_{y}^{\dagger }(z^{\prime} )\right\}-{\rm{i}}\langle {{ \mathcal T }}_{\gamma }[{\rm{i}}{\partial }_{z}{\underline{\hat{{\rm{\Phi }}}}}_{x}(z)]\otimes {\underline{\hat{{\rm{\Phi }}}}}_{y}^{\dagger }(z^{\prime} )\rangle ,\end{eqnarray} \tag{ A.4 }$$

where the anticommutator gives simply ${\delta }_{{xy}}\underline{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 [32, 33, 35, 99]

$$\begin{eqnarray}&&\left[{\rm{i}}{\partial }_{z}\underline{{\bf{1}}}-\underline{{\boldsymbol{h}}}(z)\right]\underline{{\boldsymbol{G}}}(z,z^{\prime} )=\delta (z,z^{\prime} )\underline{{\bf{1}}},\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&\underline{{\boldsymbol{G}}}(z,z^{\prime} )\left[-{\rm{i}}{\overleftarrow{\partial }}_{z^{\prime} }\underline{{\bf{1}}}-\underline{{\boldsymbol{h}}}(z^{\prime} )\right]=\delta (z,z^{\prime} )\underline{{\bf{1}}},\end{eqnarray} \tag{ A.6 }$$

which the Green's function satisfies being antiperiodic along the contour (Kubo–Martin–Schwinger boundary condition [100, 101]).

We see that the equations of motion are the same as those of [33, 99], hence we may in similar fashion, using the Langreth rules [32, 102], derive an equation for the equal-time lesser Green's function with indices on the central region ${\underline{G}}_{{CC}}^{\lt }$. This is a key quantity as it relates to the time-dependent one-particle reduced density-matrix (TD1RDM) by ${\underline{\rho }}_{{CC}}(t)=-{\rm{i}}{\underline{G}}_{{CC}}^{\lt }(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 [99]

$$\begin{eqnarray}&&{\rm{i}}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\underline{G}}^{\lt }(t,t)-[{\underline{h}}_{{CC}}(t),{\underline{G}}^{\lt }(t,t)]=-\left[{\underline{G}}^{{R}}\cdot {\underline{{\rm{\Sigma }}}}^{\lt }+{\underline{G}}^{\lt }\cdot {\underline{{\rm{\Sigma }}}}^{{A}}+{\underline{G}}^{\rceil }\,\star \,{\underline{{\rm{\Sigma }}}}^{\lceil }\right](t,t)+{\rm{h}}.{\rm{c}}.,\end{eqnarray} \tag{ A.7 }$$

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 }{\rm{d}}\bar{t}f(t,\bar{t})g(\bar{t},t)$ and $[f\,\star \,g](t,t)=-{\rm{i}}{\int }_{0}^{\beta }{\rm{d}}\tau f(t,\tau )g(\tau ,t)$. The superscripts R, A, $\lt ,\rceil ,\lceil$ refer to the retarded, advanced, lesser, right and left Keldysh components, respectively [32, 99]. The embedding self-energy, ${\rm{\Sigma }}$, accounts for the coupling between the central region and the electrodes [33].

We note that the left-hand side of equation (A.7) 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 WBA for the embedding self-energy, equation (A.7) becomes a closed equation for the equal-time lesser Green's function and the TD1RDM can be solved analytically.

In order to close the equation of motion we now describe the electrodes in the framework of 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. [21, 103–105], 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

$$\begin{eqnarray}&&{\underline{{\rm{\Sigma }}}}_{\lambda ,{mn}}^{{\rm{R}}}(\omega )=\displaystyle \sum _{k}{\underline{T}}_{{mk}\lambda }\displaystyle \frac{1}{\omega -{\underline{\epsilon }}_{k\lambda }+{\rm{i}}\eta }{\underline{T}}_{k\lambda n}\approx -{\rm{i}}{\underline{{\rm{\Gamma }}}}_{\lambda ,{mn}}/2.\end{eqnarray} \tag{ B.1 }$$

The advanced component is given simply by conjugating this. The other components of the self-energy ($\lt ,\lceil$) may further be derived from the retarded and advanced components [32, 99]. The time-domain quantities in equation (A.7) are then obtained by Fourier transforming. Looking at equation (A.7) and the earlier work [33, 99] we may use the fact that the same equations have the same solutions, i.e. including the Nambu ⨂ 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 ⨂ spin structure of the matrices.

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

$$\begin{eqnarray}&&\langle {{\rm{\Psi }}}^{{\rm{L}}}| {\underline{h}}_{\mathrm{eff}}=\epsilon \langle {{\rm{\Psi }}}^{{\rm{L}}}| ;\quad {\underline{h}}_{\mathrm{eff}}| {{\rm{\Psi }}}^{{\rm{R}}}\rangle =\epsilon | {{\rm{\Psi }}}^{{\rm{R}}}\rangle ,\end{eqnarray} \tag{ B.2 }$$

where the eigenvectors and eigenvalues correspond to the 4 × 4 Nambu ⨂ spin space. The solution for the TD1RDM expanded in the left eigenbasis takes the explicit form [33]

$$\begin{eqnarray}\begin{array}{rcl}\langle {{\rm{\Psi }}}_{j}^{{\rm{L}}}| \underline{\rho }(t)| {{\rm{\Psi }}}_{k}^{{\rm{L}}}\rangle & = & \displaystyle \sum _{\lambda }\left\{{{\rm{\Gamma }}}_{\lambda ,{jk}}{{\rm{\Lambda }}}_{\lambda ,{jk}}+{V}_{\lambda }{{\rm{\Gamma }}}_{\lambda ,{jk}}\left[{{\rm{\Pi }}}_{\lambda ,{jk}}(t)+{{\rm{\Pi }}}_{\lambda ,{kj}}^{* }(t)\right]\right.\\ & & \left.+{V}_{\lambda }^{2}{{\rm{\Gamma }}}_{\lambda ,{jk}}{{\rm{e}}}^{-{\rm{i}}({\epsilon }_{j}-{\epsilon }_{k}^{* })t}{{\rm{\Omega }}}_{\lambda ,{jk}}\right\},\end{array}\end{eqnarray} \tag{ B.3 }$$

where

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\lambda ,{jk}}=\langle {{\rm{\Psi }}}_{j}^{{\rm{L}}}| {\underline{{\rm{\Gamma }}}}_{\lambda }| {{\rm{\Psi }}}_{k}^{{\rm{L}}}\rangle ,\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\lambda ,{jk}}=\int \displaystyle \frac{{\rm{d}}\omega }{2\pi }\displaystyle \frac{f(\omega -\mu )}{(\omega +{V}_{\lambda }-{\epsilon }_{j})(\omega +{V}_{\lambda }-{\epsilon }_{k}^{* })},\end{eqnarray} \tag{ B.5 }$$

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{\lambda ,{jk}}(t)=\int \displaystyle \frac{{\rm{d}}\omega }{2\pi }\displaystyle \frac{f(\omega -\mu ){{\rm{e}}}^{{\rm{i}}(\omega +{V}_{\lambda }-{\epsilon }_{j})t}}{(\omega -{\epsilon }_{j})(\omega +{V}_{\lambda }-{\epsilon }_{j})(\omega +{V}_{\lambda }-{\epsilon }_{k}^{* })},\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\lambda ,{jk}}=\int \displaystyle \frac{{\rm{d}}\omega }{2\pi }\displaystyle \frac{f(\omega -\mu )}{(\omega -{\epsilon }_{j})(\omega +{V}_{\lambda }-{\epsilon }_{j})(\omega +{V}_{\lambda }-{\epsilon }_{k}^{* })(\omega -{\epsilon }_{k}^{* })}.\end{eqnarray} \tag{ B.7 }$$

Here $f(\omega -\mu )={\left({{\rm{e}}}^{\beta (\omega -\mu )}+1\right)}^{-1}$ is the Fermi function at inverse temperature β and chemical potential μ. 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

$$\begin{eqnarray}&&\langle {\varphi }_{m}| \underline{\rho }(t)| {\varphi }_{n}\rangle =\displaystyle \sum _{{jk}}\displaystyle \frac{\langle {\varphi }_{m}| {{\rm{\Psi }}}_{j}^{{\rm{R}}}\rangle }{\langle {{\rm{\Psi }}}_{j}^{{\rm{L}}}| {{\rm{\Psi }}}_{j}^{{\rm{R}}}\rangle }\displaystyle \frac{\langle {{\rm{\Psi }}}_{k}^{{\rm{R}}}| {\varphi }_{n}\rangle }{\langle {{\rm{\Psi }}}_{k}^{{\rm{R}}}| {{\rm{\Psi }}}_{k}^{{\rm{L}}}\rangle }\langle {{\rm{\Psi }}}_{j}^{{\rm{L}}}| \underline{\rho }(t)| {{\rm{\Psi }}}_{k}^{{\rm{L}}}\rangle ,\end{eqnarray} \tag{ B.8 }$$

which follows from the biorthogonality of the left and right eigenvectors. The TD1RDM is then simply given by evaluating the terms in equations (B.5), (B.6) and (B.7) for all indices j, k and time parameter t, and then inserting into equations (B.3) and (B.8).

The integrands in equations (B.5), (B.6) and (B.7) 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 /(-{\rm{i}}\beta )$. Expressions similar to those in equations (B.5), (B.6), (B.7) have been found, e.g. in [106, 107] and integrated correspondingly using contour integration techniques. In [33] the frequency integrals in equations (B.5), (B.6), (B.7) 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.

Making a change of variables $z=\beta (\omega -\mu )$ in equation (B.5) gives

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\lambda ,{jk}}=\beta {\int }_{-\infty }^{\infty }\displaystyle \frac{{\rm{d}}z}{2\pi }\displaystyle \frac{1}{(z-{z}_{1})(z-{z}_{2})({{\rm{e}}}^{z}+1)},\end{eqnarray} \tag{ C.1 }$$

where we defined z₁ = β(_j − μ_λ) and ${z}_{2}=\beta ({\epsilon }_{k}^{* }-{\mu }_{\lambda })$ with μ_λ = μ + V_λ. This integrand has simple poles at z = z₁, z = z₂ and $z={w}_{n}={\rm{i}}(\pi +2\pi n)$, see figure C1. The spectrum of the complex eigenvalues of the nonhermitian matrix ${\underline{h}}_{\mathrm{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 '${\left(z-{z}_{n}\right)}^{-1}$' contributions the residues are simply one and for the Fermi function we have $\mathrm{Res}\left[{\left({{\rm{e}}}^{z}+1\right)}^{-1},z={w}_{n}\right]=-1$. Then, we can close the integral in equation (C.1) in the UHP as shown in figure C1, and using the residue theorem we get

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\lambda ,{jk}}={\rm{i}}\beta \left[\displaystyle \frac{1}{{z}_{2}-{z}_{1}}\displaystyle \frac{1}{{{\rm{e}}}^{{z}_{2}}+1}-\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{({w}_{n}-{z}_{1})({w}_{n}-{z}_{2})}\right].\end{eqnarray} \tag{ C.2 }$$

The infinite sum can be written as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{\left[{\rm{i}}(\pi +2\pi n)-{z}_{1}\right]\left[{\rm{i}}(\pi +2\pi n)-{z}_{2}\right]} & = & \displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{2\pi {\rm{i}}\left(n+\tfrac{{\rm{i}}{z}_{1}+\pi }{2\pi }\right)2\pi {\rm{i}}\left(n+\tfrac{{\rm{i}}{z}_{2}+\pi }{2\pi }\right)}\\ & = & -\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{(n+a)(n+b)}\\ & = & -\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}\displaystyle \frac{1}{b-a}\left[\psi (b)-\psi (a)\right],\end{array}\end{eqnarray} \tag{ C.3 }$$

where we defined $a=({\rm{i}}{z}_{1}+\pi )/2\pi$, $b=({\rm{i}}{z}_{2}+\pi )/2\pi$, and ψ is the digamma function which is defined as the logarithmic derivative of the gamma function, $\psi (z)=\tfrac{{\rm{d}}}{{\rm{d}}z}\mathrm{log}{\rm{\Gamma }}(z)$ [108]. We can then insert the result of the sum back into equation (C.2) and couple the terms by simplifying

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\lambda ,{jk}}=\displaystyle \frac{{\rm{i}}}{{\epsilon }_{k}^{* }-{\epsilon }_{j}}\left\{\displaystyle \frac{1}{{{\rm{e}}}^{\beta ({\epsilon }_{k}^{* }-{\mu }_{\lambda })}+1}+\displaystyle \frac{1}{2\pi {\rm{i}}}\left[\psi \left(\displaystyle \frac{1}{2}-\displaystyle \frac{\beta ({\epsilon }_{k}^{* }-{\mu }_{\lambda })}{2\pi {\rm{i}}}\right)-\psi \left(\displaystyle \frac{1}{2}-\displaystyle \frac{\beta ({\epsilon }_{j}-{\mu }_{\lambda })}{2\pi {\rm{i}}}\right)\right]\right\},\end{eqnarray} \tag{ C.4 }$$

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 ψ which is broadly known in computational sciences and readily implemented for example in the GNU Scientific Library [109]. Equation (C.4) is our final result for ${{\rm{\Lambda }}}_{\lambda ,{jk}}$ for arbitrary values of β. We note in passing that it would give completely equivalent result if the integral was closed in the LHP.

Zoom In Zoom Out Reset image size

Making the same change of variables in equation (B.6) as in the previous case leads to

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{\lambda ,{jk}}(t)={\beta }^{2}{\int }_{-\infty }^{\infty }\displaystyle \frac{{\rm{d}}z}{2\pi }\displaystyle \frac{{{\rm{e}}}^{\tfrac{{\rm{i}}}{\beta }(z-{z}_{2})t}}{(z-{z}_{1})(z-{z}_{2})(z-{z}_{3})({{\rm{e}}}^{z}+1)},\end{eqnarray} \tag{ C.5 }$$

where we defined z₁ = β(_j − μ), z₂ = β(_j − μ_λ) and ${z}_{3}=\beta ({\epsilon }_{k}^{* }-{\mu }_{\lambda })$. Also in this case we notice poles in the complex plane, similarly as in figure C1. 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

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{\lambda ,{jk}}(t)={\rm{i}}{\beta }^{2}\left[\displaystyle \frac{{{\rm{e}}}^{\tfrac{{\rm{i}}}{\beta }({z}_{3}-{z}_{2})t}}{({z}_{3}-{z}_{1})({z}_{3}-{z}_{2})({{\rm{e}}}^{{z}_{3}}+1)}-\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{{\rm{e}}}^{\tfrac{{\rm{i}}}{\beta }({w}_{n}-{z}_{2})t}}{({w}_{n}-{z}_{1})({w}_{n}-{z}_{2})({w}_{n}-{z}_{3})}\right].\end{eqnarray} \tag{ C.6 }$$

We may manipulate the infinite sum in equation (C.6) as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{{\rm{e}}}^{\tfrac{{\rm{i}}}{\beta }({w}_{n}-{z}_{2})t}}{({w}_{n}-{z}_{1})({w}_{n}-{z}_{2})({w}_{n}-{z}_{3})} & = & \displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{{\rm{e}}}^{\tfrac{{\rm{i}}}{\beta }2\pi {\rm{i}}\left(n+\tfrac{{\rm{i}}{z}_{2}+\pi }{2\pi }\right)t}}{2\pi {\rm{i}}\left(n+\tfrac{{\rm{i}}{z}_{1}+\pi }{2\pi }\right)2\pi {\rm{i}}\left(n+\tfrac{{\rm{i}}{z}_{2}+\pi }{2\pi }\right)2\pi {\rm{i}}\left(n+\tfrac{{\rm{i}}{z}_{3}+\pi }{2\pi }\right)}\\ & = & \displaystyle \frac{{\rm{i}}}{{\left(2\pi \right)}^{3}}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{{\rm{e}}}^{x(n+b)}}{(n+a)(n+b)(n+c)},\end{array}\end{eqnarray} \tag{ C.7 }$$

where we defined $a=({\rm{i}}{z}_{1}+\pi )/2\pi$, $b=({\rm{i}}{z}_{2}+\pi )/2\pi$, $c=({\rm{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}$ [110]:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{i}}}{{\left(2\pi \right)}^{3}}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{{\rm{e}}}^{x(n+b)}}{(n+a)(n+b)(n+c)}=\displaystyle \frac{{\rm{i}}}{{\left(2\pi \right)}^{3}(a-b)(a-c)(b-c)}\\ \ \times \ \left\{{{\rm{e}}}^{{bx}}\left[\displaystyle \frac{b-c}{a}{}_{2}{F}_{1}(1,a,1+a,{{\rm{e}}}^{x})+\displaystyle \frac{c-a}{b}{}_{2}{F}_{1}(1,b,1+b,{{\rm{e}}}^{x})+\displaystyle \frac{a-b}{c}{}_{2}{F}_{1}(1,c,1+c,{{\rm{e}}}^{x})\right]\right\}.\end{array}\end{eqnarray} \tag{ C.8 }$$

The hypergeometric function together with the Pochhammer symbol are defined as [110, 111]

$$\begin{eqnarray}{}_{2}{F}_{1}\left(p,q,r,s\right)=\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{\left(p\right)}_{n}{\left(q\right)}_{n}}{{\left(r\right)}_{n}}\displaystyle \frac{{s}^{n}}{n!},\ \ {\left(p\right)}_{n}=\left\{\begin{array}{cl}1 & n=1,\\ p(p+1)\cdots (p+n-1) & n\gt 0.\end{array}\right.\end{eqnarray} \tag{ C.9 }$$

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

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{{\rm{e}}}^{\tfrac{{\rm{i}}}{\beta }({w}_{n}-{z}_{2})t}}{({w}_{n}-{z}_{1})({w}_{n}-{z}_{2})({w}_{n}-{z}_{3})}=\displaystyle \frac{-{\mathrm{ie}}^{-\pi t/\beta }{{\rm{e}}}^{-{\rm{i}}({\epsilon }_{j}-{\mu }_{\lambda })t}}{{\beta }^{2}({\epsilon }_{k}^{* }-{\epsilon }_{j})({\epsilon }_{k}^{* }-{\epsilon }_{j}-{V}_{\lambda })}\\ \ \times \ \left\{{\mathfrak{F}}({\epsilon }_{k}^{* }-{\mu }_{\lambda },t,\beta )+\displaystyle \frac{{\epsilon }_{k}^{* }-{\epsilon }_{j}-{V}_{\lambda }}{{V}_{\lambda }}{\mathfrak{F}}({\epsilon }_{j}-{\mu }_{\lambda },t,\beta )-\displaystyle \frac{{\epsilon }_{k}^{* }-{\epsilon }_{j}}{{V}_{\lambda }}{\mathfrak{F}}({\epsilon }_{j}-\mu ,t,\beta )\right\},\end{array}\end{eqnarray} \tag{ C.10 }$$

where we defined an auxiliary function

$$\begin{eqnarray}&&{\mathfrak{F}}(z,t,\beta )\equiv \displaystyle \frac{1}{{\rm{i}}\beta z+\pi }\ {}_{2}{F}_{1}\left(1,\displaystyle \frac{1}{2}+\displaystyle \frac{{\rm{i}}\beta z}{2\pi },\displaystyle \frac{3}{2}+\displaystyle \frac{{\rm{i}}\beta z}{2\pi },{{\rm{e}}}^{-2\pi t/\beta }\right).\end{eqnarray} \tag{ C.11 }$$

This calculation was only for the infinite sum in equation (C.6). Inserting the definitions of z's into the first term gives

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{e}}}^{\tfrac{{\rm{i}}}{\beta }({z}_{3}-{z}_{2})t}}{({z}_{3}-{z}_{1})({z}_{3}-{z}_{2})({{\rm{e}}}^{{z}_{3}}+1)}=\displaystyle \frac{{{\rm{e}}}^{-{\rm{i}}({\epsilon }_{j}-{\epsilon }_{k}^{* })t}}{{\beta }^{2}({\epsilon }_{k}^{* }-{\epsilon }_{j})({\epsilon }_{k}^{* }-{\epsilon }_{j}-{V}_{\lambda })}\displaystyle \frac{1}{{{\rm{e}}}^{\beta ({\epsilon }_{k}^{* }-{\mu }_{\lambda })}+1}.\end{eqnarray} \tag{ C.12 }$$

Combining the terms finally gives

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Pi }}}_{\lambda ,{jk}}(t) & = & \displaystyle \frac{{\rm{i}}}{({\epsilon }_{k}^{* }-{\epsilon }_{j})({\epsilon }_{k}^{* }-{\epsilon }_{j}-{V}_{\lambda })}\left\{\displaystyle \frac{{{\rm{e}}}^{-{\rm{i}}({\epsilon }_{j}-{\epsilon }_{k}^{* })t}}{{{\rm{e}}}^{\beta ({\epsilon }_{k}^{* }-{\mu }_{\lambda })}+1}+{\mathrm{ie}}^{-\pi t/\beta }{{\rm{e}}}^{-{\rm{i}}({\epsilon }_{j}-{\mu }_{\lambda })t}\right.\\ & & \left.\times \left[{\mathfrak{F}}({\epsilon }_{k}^{* }-{\mu }_{\lambda },t,\beta )+\displaystyle \frac{{\epsilon }_{k}^{* }-{\epsilon }_{j}-{V}_{\lambda }}{{V}_{\lambda }}{\mathfrak{F}}({\epsilon }_{j}-{\mu }_{\lambda },t,\beta )-\displaystyle \frac{{\epsilon }_{k}^{* }-{\epsilon }_{j}}{{V}_{\lambda }}{\mathfrak{F}}({\epsilon }_{j}-\mu ,t,\beta )\right]\right\}\end{array}\end{eqnarray} \tag{ C.13 }$$

for arbitrary values of β. 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 [112].

In the third case, in equation (B.7), we do the same change of variables as before to get

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\lambda ,{jk}}={\beta }^{3}{\int }_{-\infty }^{\infty }\displaystyle \frac{{\rm{d}}z}{2\pi }\displaystyle \frac{1}{(z-{z}_{1})(z-{z}_{2})(z-{z}_{3})(z-{z}_{4})({{\rm{e}}}^{z}+1)},\end{eqnarray} \tag{ C.14 }$$

where we defined z₁ = β(_j − μ), ${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 figure C1, and we may close also this integral in the UHP. Again, according to the residue theorem we get as a result

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Omega }}}_{\lambda ,{jk}} & = & {\rm{i}}{\beta }^{3}\left\{\displaystyle \frac{1}{({z}_{3}-{z}_{1})({z}_{3}-{z}_{2})({z}_{3}-{z}_{4})({{\rm{e}}}^{{z}_{3}}+1)}+\displaystyle \frac{1}{({z}_{4}-{z}_{1})({z}_{4}-{z}_{2})({z}_{4}-{z}_{3})({{\rm{e}}}^{{z}_{4}}+1)}\right.\\ & & \left.-\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{({w}_{n}-{z}_{1})({w}_{n}-{z}_{2})({w}_{n}-{z}_{3})({w}_{n}-{z}_{4})}\right\}.\end{array}\end{eqnarray} \tag{ C.15 }$$

The infinite sum may again be manipulated as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{\left[{\rm{i}}(\pi +2\pi n)-{z}_{1}\right]\left[{\rm{i}}(\pi +2\pi n)-{z}_{2}\right]\left[{\rm{i}}(\pi +2\pi n)-{z}_{3}\right]\left[{\rm{i}}(\pi +2\pi n)-{z}_{4}\right]}\\ =\ \displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{2\pi {\rm{i}}\left(n+\tfrac{{\rm{i}}{z}_{1}+\pi }{2\pi }\right)2\pi {\rm{i}}\left(n+\tfrac{{\rm{i}}{z}_{2}+\pi }{2\pi }\right)2\pi {\rm{i}}\left(n+\tfrac{{\rm{i}}{z}_{3}+\pi }{2\pi }\right)2\pi {\rm{i}}\left(n+\tfrac{{\rm{i}}{z}_{4}+\pi }{2\pi }\right)}\\ =\ \displaystyle \frac{1}{{\left(2\pi \right)}^{4}}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{(n+a)(n+b)(n+c)(n+d)},\end{array}\end{eqnarray} \tag{ C.16 }$$

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

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{{\left(2\pi \right)}^{4}}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{(n+a)(n+b)(n+c)(n+d)}=\displaystyle \frac{1}{{\left(2\pi \right)}^{4}}\left[\displaystyle \frac{\psi (a)}{(a-b)(a-c)(a-d)}+\displaystyle \frac{\psi (b)}{(b-a)(b-c)(b-d)}\right.\\ \ \ \left.+\ \displaystyle \frac{\psi (c)}{(c-a)(c-b)(c-d)}+\displaystyle \frac{\psi (d)}{(d-a)(d-b)(d-c)}\right].\end{array}\end{eqnarray} \tag{ C.17 }$$

Inserting the expressions for a, b, c and d, and then further the expressions for z₁, z₂, z₃ and z₄ leads to

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{1}{({w}_{n}-{z}_{1})({w}_{n}-{z}_{2})({w}_{n}-{z}_{3})({w}_{n}-{z}_{4})}\\ =\ \displaystyle \frac{{\rm{i}}}{2\pi {\beta }^{3}}\left\{\displaystyle \frac{1}{({\epsilon }_{k}^{* }-{\epsilon }_{j})({\epsilon }_{k}^{* }-{\epsilon }_{j}-{V}_{\lambda }){V}_{\lambda }}\left[\psi \left(\displaystyle \frac{1}{2}+\displaystyle \frac{{\rm{i}}\beta ({\epsilon }_{j}-\mu )}{2\pi }\right)-\psi \left(\displaystyle \frac{1}{2}+\displaystyle \frac{{\rm{i}}\beta ({\epsilon }_{k}^{* }-{\mu }_{\lambda })}{2\pi }\right)\right]\right.\\ \left.-\ \displaystyle \frac{1}{({\epsilon }_{k}^{* }-{\epsilon }_{j})({\epsilon }_{k}^{* }-{\epsilon }_{j}+{V}_{\lambda }){V}_{\lambda }}\left[\psi \left(\displaystyle \frac{1}{2}+\displaystyle \frac{{\rm{i}}\beta ({\epsilon }_{j}-{\mu }_{\lambda })}{2\pi }\right)-\psi \left(\displaystyle \frac{1}{2}+\displaystyle \frac{{\rm{i}}\beta ({\epsilon }_{k}^{* }-\mu )}{2\pi }\right)\right]\right\}.\end{array}\end{eqnarray} \tag{ C.18 }$$

Combining the terms in equation (C.15) gives as the final result

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Omega }}}_{\lambda ,{jk}} & = & \displaystyle \frac{\tfrac{{\rm{i}}}{{{\rm{e}}}^{\beta ({\epsilon }_{k}^{* }-\mu )}+1}-\tfrac{1}{2\pi }\left[\psi \left(\tfrac{1}{2}-\tfrac{\beta ({\epsilon }_{j}-{\mu }_{\lambda })}{2\pi {\rm{i}}}\right)-\psi \left(\tfrac{1}{2}-\tfrac{\beta ({\epsilon }_{k}^{* }-\mu )}{2\pi {\rm{i}}}\right)\right]}{({\epsilon }_{k}^{* }-{\epsilon }_{j})({\epsilon }_{k}^{* }-{\epsilon }_{j}+{V}_{\lambda }){V}_{\lambda }}\\ & & -\displaystyle \frac{\tfrac{{\rm{i}}}{{{\rm{e}}}^{\beta ({\epsilon }_{k}^{* }-{\mu }_{\lambda })}+1}-\tfrac{1}{2\pi }\left[\psi \left(\tfrac{1}{2}-\tfrac{\beta ({\epsilon }_{j}-\mu )}{2\pi {\rm{i}}}\right)-\psi \left(\tfrac{1}{2}-\tfrac{\beta ({\epsilon }_{k}^{* }-{\mu }_{\lambda })}{2\pi {\rm{i}}}\right)\right]}{({\epsilon }_{k}^{* }-{\epsilon }_{j})({\epsilon }_{k}^{* }-{\epsilon }_{j}-{V}_{\lambda }){V}_{\lambda }}\end{array}\end{eqnarray} \tag{ C.19 }$$

for arbitrary values of β.

Finally, inserting equations (C.4), (C.13) and (C.19) into (B.3) gives then the TD1RDM at arbitrary temperature. When the asymptotic behavior of the digamma and hypergeometric function is studied, the results in equations (C.4), (C.13) and (C.19) can be shown to reduce to those in [33] at the zero-temperature limit ($\beta \to \infty$) [113]. We also note that congruent results involving equivalent special functions have recently been reported in [114–116].

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}$ (equation (6) in the main text). Also, for the 'nearest-neighbor' contribution $\underline{b}$ (equation (7) in the main text), it is possible to consider a Peierls phase γ_mn = −γ_nm accounting for a magnetic field (normalized to the flux quantum ϕ₀ = 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 equation (B.3) can be derived as [33]

$$\begin{eqnarray}&&\langle {\widetilde{{\rm{\Psi }}}}_{j}^{{\rm{L}}}| \widetilde{\underline{\rho }}(t)| {\widetilde{{\rm{\Psi }}}}_{k}^{{\rm{L}}}\rangle =\displaystyle \sum _{\lambda }\left[{\widetilde{{\rm{\Gamma }}}}_{\lambda ,{jk}}{\widetilde{{\rm{\Lambda }}}}_{\lambda ,{jk}}+{\widetilde{{\rm{\Pi }}}}_{\lambda ,{jk}}(t)+{\widetilde{{\rm{\Pi }}}}_{\lambda ,{kj}}^{* }(t)+{\widetilde{{\rm{\Omega }}}}_{\lambda ,{jk}}(t)\right],\end{eqnarray} \tag{ D.1 }$$

where the introduced terms $\widetilde{{\rm{\Gamma }}}$, $\widetilde{{\rm{\Pi }}}$ and $\widetilde{{\rm{\Omega }}}$ take a slightly more intricate form compared to those in equation (B.3) 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

$$\begin{eqnarray}&&{\widetilde{{\rm{\Gamma }}}}_{\lambda ,{jk}}=\langle {\widetilde{{\rm{\Psi }}}}_{j}^{{\rm{L}}}| {\underline{{\rm{\Gamma }}}}_{\lambda }| {\widetilde{{\rm{\Psi }}}}_{k}^{{\rm{L}}}\rangle ,\end{eqnarray} \tag{ D.2 }$$

$$\begin{eqnarray}&&{\widetilde{{\rm{\Lambda }}}}_{\lambda ,{jk}}=\int \displaystyle \frac{{\rm{d}}\omega }{2\pi }\displaystyle \frac{f(\omega -\mu )}{(\omega +{V}_{\lambda }-{\widetilde{\epsilon }}_{j})(\omega +{V}_{\lambda }-{\widetilde{\epsilon }}_{k}^{* })},\end{eqnarray} \tag{ D.3 }$$

$$\begin{eqnarray}&&{\widetilde{{\rm{\Pi }}}}_{\lambda ,{jk}}(t)=\displaystyle \sum _{{mn}}\displaystyle \frac{\langle {\widetilde{{\rm{\Psi }}}}_{j}^{{\rm{L}}}| {{\rm{\Psi }}}_{m}^{{\rm{R}}}\rangle \langle {{\rm{\Psi }}}_{m}^{{\rm{L}}}| {\widetilde{\underline{V}}}_{\lambda }| {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{R}}}\rangle {\widetilde{{\rm{\Gamma }}}}_{\lambda ,{nk}}}{\langle {{\rm{\Psi }}}_{m}^{{\rm{L}}}| {{\rm{\Psi }}}_{m}^{{\rm{R}}}\rangle \langle {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{L}}}| {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{R}}}\rangle }\int \displaystyle \frac{{\rm{d}}\omega }{2\pi }\displaystyle \frac{f(\omega -\mu ){{\rm{e}}}^{{\rm{i}}(\omega +{V}_{\lambda }-{\widetilde{\epsilon }}_{j})t}}{(\omega -{\epsilon }_{m})(\omega +{V}_{\lambda }-{\widetilde{\epsilon }}_{n})(\omega +{V}_{\lambda }-{\widetilde{\epsilon }}_{k}^{* })},\end{eqnarray} \tag{ D.4 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\widetilde{{\rm{\Omega }}}}_{\lambda ,{jk}}(t) & = & \displaystyle \sum _{{mnpq}}\displaystyle \frac{\langle {\widetilde{{\rm{\Psi }}}}_{j}^{{\rm{L}}}| {{\rm{\Psi }}}_{m}^{{\rm{R}}}\rangle \langle {{\rm{\Psi }}}_{m}^{{\rm{L}}}| {\widetilde{\underline{V}}}_{\lambda }| {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{R}}}\rangle {\widetilde{{\rm{\Gamma }}}}_{\lambda ,{np}}\langle {\widetilde{{\rm{\Psi }}}}_{p}^{{\rm{R}}}| {\widetilde{\underline{V}}}_{\lambda }^{\dagger }| {{\rm{\Psi }}}_{q}^{{\rm{L}}}\rangle \langle {{\rm{\Psi }}}_{q}^{{\rm{R}}}| {\widetilde{{\rm{\Psi }}}}_{k}^{{\rm{L}}}\rangle }{\langle {{\rm{\Psi }}}_{m}^{{\rm{L}}}| {{\rm{\Psi }}}_{m}^{{\rm{R}}}\rangle \langle {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{L}}}| {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{R}}}\rangle \langle {\widetilde{{\rm{\Psi }}}}_{p}^{{\rm{R}}}| {\widetilde{{\rm{\Psi }}}}_{p}^{{\rm{L}}}\rangle \langle {{\rm{\Psi }}}_{q}^{{\rm{R}}}| {{\rm{\Psi }}}_{q}^{{\rm{L}}}\rangle }\\ & & \times {{\rm{e}}}^{-{\rm{i}}({\widetilde{\epsilon }}_{j}-{\widetilde{\epsilon }}_{k}^{* })t}\displaystyle \int \displaystyle \frac{{\rm{d}}\omega }{2\pi }\displaystyle \frac{f(\omega -\mu )}{(\omega -{\epsilon }_{m})(\omega +{V}_{\lambda }-{\widetilde{\epsilon }}_{n})(\omega +{V}_{\lambda }-{\widetilde{\epsilon }}_{p}^{* })(\omega -{\epsilon }_{q}^{* })},\end{array}\end{eqnarray} \tag{ D.5 }$$

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{1}-({\widetilde{\underline{h}}}_{{CC}}-{\underline{h}}_{{CC}})$. The eigenvalues $\{\epsilon ,\widetilde{\epsilon }\}$ and eigenvectors $\{{{\rm{\Psi }}}^{{\rm{L}}/{\rm{R}}},{\widetilde{{\rm{\Psi }}}}^{{\rm{L}}/{\rm{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}}-{\rm{i}}\underline{{\rm{\Gamma }}}/2$, respectively. In the limit ${\widetilde{\underline{h}}}_{{CC}}\to {\underline{h}}_{{CC}}$ the result in equation (D.1) can be checked to reduce to equation (B.3) [113].

Similarly, for the TD1RDM with sudden electromagnetic fields in the central region in equation (D.1), we can take the integrals in equations (D.3), (D.4) and (D.5) 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 β the explicit results are

$$\begin{eqnarray}&&{\widetilde{{\rm{\Lambda }}}}_{\lambda ,{jk}}=\displaystyle \frac{{\rm{i}}}{{\widetilde{\epsilon }}_{k}^{* }-{\widetilde{\epsilon }}_{j}}\left\{\displaystyle \frac{1}{{{\rm{e}}}^{\beta ({\widetilde{\epsilon }}_{k}^{* }-{\mu }_{\lambda })}+1}+\displaystyle \frac{1}{2\pi {\rm{i}}}\left[\psi \left(\displaystyle \frac{1}{2}-\displaystyle \frac{\beta ({\widetilde{\epsilon }}_{k}^{* }-{\mu }_{\lambda })}{2\pi {\rm{i}}}\right)-\psi \left(\displaystyle \frac{1}{2}-\displaystyle \frac{\beta ({\widetilde{\epsilon }}_{j}-{\mu }_{\lambda })}{2\pi {\rm{i}}}\right)\right]\right\},\end{eqnarray} \tag{ D.6 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\widetilde{{\rm{\Pi }}}}_{\lambda ,{jk}}(t)=\displaystyle \sum _{{mn}}\displaystyle \frac{\langle {\widetilde{{\rm{\Psi }}}}_{j}^{{\rm{L}}}| {{\rm{\Psi }}}_{m}^{{\rm{R}}}\rangle \langle {{\rm{\Psi }}}_{m}^{{\rm{L}}}| {\widetilde{\underline{V}}}_{\lambda }| {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{R}}}\rangle {\widetilde{{\rm{\Gamma }}}}_{\lambda ,{nk}}}{\langle {{\rm{\Psi }}}_{m}^{{\rm{L}}}| {{\rm{\Psi }}}_{m}^{{\rm{R}}}\rangle \langle {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{L}}}| {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{R}}}\rangle }\\ \ \times \ \displaystyle \frac{{\rm{i}}}{({\widetilde{\epsilon }}_{k}^{* }-{\widetilde{\epsilon }}_{n})({\widetilde{\epsilon }}_{k}^{* }-{\epsilon }_{m}-{V}_{\lambda })}\left\{\displaystyle \frac{{{\rm{e}}}^{-{\rm{i}}({\widetilde{\epsilon }}_{j}-{\widetilde{\epsilon }}_{k}^{* })t}}{{{\rm{e}}}^{\beta ({\widetilde{\epsilon }}_{k}^{* }-{\mu }_{\lambda })}+1}+{\mathrm{ie}}^{-\pi t/\beta }{{\rm{e}}}^{-{\rm{i}}({\widetilde{\epsilon }}_{j}-{\mu }_{\lambda })t}\right.\\ \left.\times \ \left[{\mathfrak{F}}({\widetilde{\epsilon }}_{k}^{* }-{\mu }_{\lambda },t,\beta )-\displaystyle \frac{{\widetilde{\epsilon }}_{k}^{* }-{\epsilon }_{m}-{V}_{\lambda }}{{\widetilde{\epsilon }}_{n}-{\epsilon }_{m}-{V}_{\lambda }}{\mathfrak{F}}({\widetilde{\epsilon }}_{n}-{\mu }_{\lambda },t,\beta )+\displaystyle \frac{{\widetilde{\epsilon }}_{k}^{* }-{\widetilde{\epsilon }}_{n}}{{\widetilde{\epsilon }}_{n}-{\epsilon }_{m}-{V}_{\lambda }}{\mathfrak{F}}({\epsilon }_{m}-\mu ,t,\beta )\right]\right\},\end{array}\end{eqnarray} \tag{ D.7 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\widetilde{{\rm{\Omega }}}}_{\lambda ,{jk}}=\displaystyle \sum _{{mnpq}}\displaystyle \frac{\langle {\widetilde{{\rm{\Psi }}}}_{j}^{{\rm{L}}}| {{\rm{\Psi }}}_{m}^{{\rm{R}}}\rangle \langle {{\rm{\Psi }}}_{m}^{{\rm{L}}}| {\widetilde{\underline{V}}}_{\lambda }| {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{R}}}\rangle {\widetilde{{\rm{\Gamma }}}}_{\lambda ,{np}}\langle {\widetilde{{\rm{\Psi }}}}_{p}^{{\rm{R}}}| {\widetilde{\underline{V}}}_{\lambda }^{\dagger }| {{\rm{\Psi }}}_{q}^{{\rm{L}}}\rangle \langle {{\rm{\Psi }}}_{q}^{{\rm{R}}}| {\widetilde{{\rm{\Psi }}}}_{k}^{{\rm{L}}}\rangle }{\langle {{\rm{\Psi }}}_{m}^{{\rm{L}}}| {{\rm{\Psi }}}_{m}^{{\rm{R}}}\rangle \langle {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{L}}}| {\widetilde{{\rm{\Psi }}}}_{n}^{{\rm{R}}}\rangle \langle {\widetilde{{\rm{\Psi }}}}_{p}^{{\rm{R}}}| {\widetilde{{\rm{\Psi }}}}_{p}^{{\rm{L}}}\rangle \langle {{\rm{\Psi }}}_{q}^{{\rm{R}}}| {{\rm{\Psi }}}_{q}^{{\rm{L}}}\rangle }\ {{\rm{e}}}^{-{\rm{i}}({\widetilde{\epsilon }}_{j}-{\widetilde{\epsilon }}_{k}^{* })t}\\ \ \ \times \ \left\{\displaystyle \frac{1}{({\epsilon }_{m}-{\widetilde{\epsilon }}_{n}+{V}_{\lambda })({\epsilon }_{m}-{\widetilde{\epsilon }}_{p}^{* }+{V}_{\lambda })({\epsilon }_{m}-{\epsilon }_{q}^{* })}\displaystyle \frac{1}{2\pi }\psi \left(\displaystyle \frac{1}{2}-\displaystyle \frac{\beta ({\epsilon }_{m}-\mu )}{2\pi {\rm{i}}}\right)\right.\\ \ \ +\ \displaystyle \frac{1}{({\widetilde{\epsilon }}_{n}-{\epsilon }_{m}-{V}_{\lambda })({\widetilde{\epsilon }}_{n}-{\widetilde{\epsilon }}_{p}^{* })({\widetilde{\epsilon }}_{n}-{\epsilon }_{q}^{* }-{V}_{\lambda })}\displaystyle \frac{1}{2\pi }\psi \left(\displaystyle \frac{1}{2}-\displaystyle \frac{\beta ({\widetilde{\epsilon }}_{n}-{\mu }_{\lambda })}{2\pi {\rm{i}}}\right)\\ \ \ +\ \displaystyle \frac{1}{({\epsilon }_{q}^{* }-{\epsilon }_{m})({\epsilon }_{q}^{* }-{\widetilde{\epsilon }}_{n}+{V}_{\lambda })({\epsilon }_{q}^{* }-{\widetilde{\epsilon }}_{p}^{* }+{V}_{\lambda })}\left[\displaystyle \frac{{\rm{i}}}{{{\rm{e}}}^{\beta ({\epsilon }_{q}^{* }-\mu )}+1}+\displaystyle \frac{1}{2\pi }\psi \left(\displaystyle \frac{1}{2}-\displaystyle \frac{\beta ({\epsilon }_{q}^{* }-\mu )}{2\pi {\rm{i}}}\right)\right]\\ \ \ \left.+\ \displaystyle \frac{1}{({\widetilde{\epsilon }}_{p}^{* }-{\epsilon }_{m}-{V}_{\lambda })({\widetilde{\epsilon }}_{p}^{* }-{\widetilde{\epsilon }}_{n})({\widetilde{\epsilon }}_{p}^{* }-{\epsilon }_{q}^{* }-{V}_{\lambda })}\left[\displaystyle \frac{{\rm{i}}}{{{\rm{e}}}^{\beta ({\widetilde{\epsilon }}_{p}^{* }-{\mu }_{\lambda })}+1}+\displaystyle \frac{1}{2\pi }\psi \left(\displaystyle \frac{1}{2}-\displaystyle \frac{\beta ({\widetilde{\epsilon }}_{p}^{* }-{\mu }_{\lambda })}{2\pi {\rm{i}}}\right)\right]\right\}.\end{array}\end{eqnarray} \tag{ D.8 }$$

Again, inserting equations (D.6), (D.7) and (D.8) into (D.1) gives then the TD1RDM for a perturbed central region at arbitrary temperature. Also here, the zero-temperature limit ($\beta \to \infty$) presented in [33], is recovered by the asymptotics of the digamma and hypergeometric functions in equations (D.6), (D.7) and (D.8) [113]. By careful inspection of equations (D.6), (D.7) and (D.8) in the limit of unperturbed central region ($\widetilde{{\rm{\Psi }}}\to {\rm{\Psi }}$ and $\widetilde{\epsilon }\to \epsilon$) it can be verified that they reduce to equations (C.4), (C.13) and (C.19) [113].

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:

$$\begin{eqnarray}&&{\hat{N}}_{j}=\displaystyle \sum _{{ks}}{\hat{c}}_{{kLs}}^{\dagger }{\hat{c}}_{{kLs}}+\displaystyle \sum _{m=1}^{j}\displaystyle \sum _{s}{\hat{c}}_{{ms}}^{\dagger }{\hat{c}}_{{ms}},\end{eqnarray} \tag{ E.1 }$$

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

$$\begin{eqnarray}&&{I}_{j,j+1}\equiv \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle {\hat{N}}_{j}\rangle .\end{eqnarray} \tag{ E.2 }$$

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

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle {\hat{N}}_{j}\rangle =-{\rm{i}}\left\langle \left[{\hat{N}}_{j},{\hat{H}}_{\mathrm{tot}}\right]\right\rangle ,\end{eqnarray} \tag{ E.3 }$$

where we now separate the 'normal' and 'superconducting' contributions as ${\hat{H}}_{\mathrm{tot}}={\hat{H}}_{\mathrm{nor}}+{\hat{H}}_{{\rm{\Delta }}}$ with

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{nor}}={\hat{H}}_{L}+{\hat{H}}_{R}+{\hat{H}}_{c}+\displaystyle \sum _{n=1}^{{N}_{w}}\displaystyle \sum _{{ss}^{\prime} }\left({\epsilon }^{{ss}^{\prime} }{\hat{c}}_{{ns}}^{\dagger }{\hat{c}}_{{ns}^{\prime} }+{\rm{h}}.{\rm{c}}.\right)+\displaystyle \sum _{n=1}^{{N}_{w}-1}\displaystyle \sum _{{ss}^{\prime} }\left({J}^{{ss}^{\prime} }{\hat{c}}_{{ns}}^{\dagger }{\hat{c}}_{(n+1)s^{\prime} }+{\rm{h}}.{\rm{c}}.\right),\end{eqnarray} \tag{ E.4 }$$

$$\begin{eqnarray}&&{\hat{H}}_{{\rm{\Delta }}}={\rm{\Delta }}\displaystyle \sum _{n=1}^{{N}_{w}}{\hat{c}}_{n\uparrow }{\hat{c}}_{n\downarrow }+{\rm{h}}.{\rm{c}}.,\end{eqnarray} \tag{ E.5 }$$

where we separated the 'on-site' and 'nearest-neighbor' contributions in the spin-dependent matrix elements of 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

$$\begin{eqnarray}&&{\hat{H}}_{c}=\displaystyle \sum _{{kss}^{\prime} }\left[\left({T}_{{kL}1}^{{ss}^{\prime} }{\hat{c}}_{{kLs}}^{\dagger }{\hat{c}}_{1s^{\prime} }+{\rm{h}}.{\rm{c}}.\right)+\left({T}_{{{kRN}}_{w}}^{{ss}^{\prime} }{\hat{c}}_{{kRs}}^{\dagger }{\hat{c}}_{{N}_{w}s^{\prime} }+{\rm{h}}.{\rm{c}}.\right)\right].\end{eqnarray} \tag{ E.6 }$$

With the normal part of the Hamiltonian, the commutator in equation (E.3) is nonzero only for the following terms

$$\begin{eqnarray}\begin{array}{rcl}\left[{\hat{N}}_{j},{\hat{H}}_{\mathrm{nor}}\right] & = & \left[\displaystyle \sum _{{ks}}{\hat{c}}_{{kLs}}^{\dagger }{\hat{c}}_{{kLs}},{\hat{H}}_{c}\right]+\left[\displaystyle \sum _{m=1}^{j}\displaystyle \sum _{s}{\hat{c}}_{{ms}}^{\dagger }{\hat{c}}_{{ms}},{\hat{H}}_{c}\right]\\ & & +\left[\displaystyle \sum _{m=1}^{j}\displaystyle \sum _{s}{\hat{c}}_{{ms}}^{\dagger }{\hat{c}}_{{ms}},\displaystyle \sum _{n=1}^{{N}_{w}-1}\displaystyle \sum _{s^{\prime} s^{\prime\prime} }\left({J}^{s^{\prime} s^{\prime\prime} }{\hat{c}}_{{ns}^{\prime} }^{\dagger }{\hat{c}}_{(n+1)s^{\prime\prime} }+{J}^{s^{\prime\prime} s^{\prime} }{\hat{c}}_{(n+1)s^{\prime\prime} }^{\dagger }{\hat{c}}_{{ns}^{\prime} }\right)\right].\end{array}\end{eqnarray} \tag{ E.7 }$$

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

$$\begin{eqnarray}&&\left[{\hat{N}}_{j},{\hat{H}}_{\mathrm{nor}}\right]=\displaystyle \sum _{{ss}^{\prime} }{J}^{{ss}^{\prime} }\left({\hat{c}}_{{js}}^{\dagger }{\hat{c}}_{(j+1)s^{\prime} }-{\hat{c}}_{(j+1)s}^{\dagger }{\hat{c}}_{{js}^{\prime} }\right).\end{eqnarray} \tag{ E.8 }$$

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

$$\begin{eqnarray}\begin{array}{rcl}\left[{\hat{N}}_{j},{\hat{H}}_{{\rm{\Delta }}}\right] & = & \left[\displaystyle \sum _{m=1}^{j}\displaystyle \sum _{s}{\hat{c}}_{{ms}}^{\dagger }{\hat{c}}_{{ms}},{\rm{\Delta }}\displaystyle \sum _{n=1}^{{N}_{w}}{\hat{c}}_{n\uparrow }{\hat{c}}_{n\downarrow }+{\rm{h}}.{\rm{c}}.\right]\\ & = & \displaystyle \sum _{m=1}^{j}\left({\rm{\Delta }}{\hat{c}}_{m\downarrow }{\hat{c}}_{m\uparrow }-{\rm{\Delta }}{\hat{c}}_{m\uparrow }{\hat{c}}_{m\downarrow }-{{\rm{\Delta }}}^{* }{\hat{c}}_{m\uparrow }^{\dagger }{\hat{c}}_{m\downarrow }^{\dagger }+{{\rm{\Delta }}}^{* }{\hat{c}}_{m\downarrow }^{\dagger }{\hat{c}}_{m\uparrow }^{\dagger }\right).\end{array}\end{eqnarray} \tag{ E.9 }$$

We may then insert equations (E.8) and (E.9) into (E.3) and further into equation (E.2), and use the model parameters for the nanowire (we have assumed a real pairing field) to obtain:

$$\begin{eqnarray}\begin{array}{rcl}{I}_{j,j+1} & = & -{\rm{i}}\left\langle -\displaystyle \frac{J}{2}\left[{\hat{c}}_{j\uparrow }^{\dagger }{\hat{c}}_{(j+1)\uparrow }+{\hat{c}}_{j\downarrow }^{\dagger }{\hat{c}}_{(j+1)\downarrow }-\left({\hat{c}}_{(j+1)\uparrow }^{\dagger }{\hat{c}}_{j\uparrow }+{\hat{c}}_{(j+1)\downarrow }^{\dagger }{\hat{c}}_{j\downarrow }\right)\right]\right.\\ & & -\displaystyle \frac{\alpha }{2}\left[{\hat{c}}_{j\uparrow }^{\dagger }{\hat{c}}_{(j+1)\downarrow }-{\hat{c}}_{j\downarrow }^{\dagger }{\hat{c}}_{(j+1)\uparrow }-\left({\hat{c}}_{(j+1)\downarrow }^{\dagger }{\hat{c}}_{j\uparrow }-{\hat{c}}_{(j+1)\uparrow }^{\dagger }{\hat{c}}_{j\downarrow }\right)\right]\\ & & \left.+\displaystyle \sum _{m=1}^{j}{\rm{\Delta }}\left({\hat{c}}_{m\downarrow }{\hat{c}}_{m\uparrow }-{\hat{c}}_{m\uparrow }{\hat{c}}_{m\downarrow }-{\hat{c}}_{m\uparrow }^{\dagger }{\hat{c}}_{m\downarrow }^{\dagger }+{\hat{c}}_{m\downarrow }^{\dagger }{\hat{c}}_{m\uparrow }^{\dagger }\right)\right\rangle .\end{array}\end{eqnarray} \tag{ E.10 }$$

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