Nonlocal conductance spectroscopy of Andreev bound states:
Symmetry relations and BCS charges

Tunneling spectroscopy is a well-established tool for studying normal metal-superconductor (NS) hybrid systems. In the context of topological superconductors, tunneling spectroscopy is widely used in attempts to identify Majorana bound states (MBSs) Mourik et al. (2012); Das et al. (2012); Deng et al. (2012, 2016); Nichele et al. (2017); Zhang et al. (2018); Vai , the prediction being that a single isolated Majorana mode should yield a zero-bias peak that is quantized to a conductance of 2$e^{2}/h$ for temperatures much below the scale of the tunneling broadening Sengupta et al. (2001); Law et al. (2009); Flensberg (2010). For overlapping Majorana states, the overlap gives distinct features in the two-probe conductance Flensberg (2010); Hansen et al. (2016) which become very pronounced when probed with a quantum dot Leijnse and Flensberg (2011); Deng et al. (2016); Prada et al. (2017); Clarke (2017); Deng et al. (2018).

While standard two-probe tunneling spectroscopy, with one normal and one grounded superconducting probe, can provide information about the subgap spectrum, it also has a severe limitation in that the interpretation of the data is ambiguous in the context of Majorana wires: Local Andreev states or so-called quasi-MBSs can give signatures that strongly resemble those of a truly topological zero mode Kells et al. (2012); Prada et al. (2012); Fle (2018); Liu et al. (2017); Vui . The reason for this is essentially that one local probe can, on general grounds, not confirm the true nonlocal nature of the MBS.

With the limitations of such a two-probe measurement, it is natural to investigate other types of finite-bias spectroscopy to access the nonlocal properties of the subgap states. One approach is to use Coulomb-blockaded Majorana islands, where two normal probes are connected to the ends of the island and the proximitizing superconductor is floating Albrecht et al. (2016, 2017); Vaitiekenas et al. (2018a); Vai ; Vaitiekenas et al. (2018b); Shen et al. (2018). Linear-response sequential transport is then possible only through states that have support at both ends of the island. The spacing between the peaks in the zero-bias differential conductance as a function of a gate-induced potential offset provides information about the energy of the lowest-energy bound state on the island van Heck et al. (2016); Chiu et al. (2017). Experiments have yielded results that could be consistent with the presence of overlapping exponentially localized MBSs Albrecht et al. (2016); Shen et al. (2018); Vai and also provided information about the quasiparticle dynamics on the island Higginbotham et al. (2015); Albrecht et al. (2017). However, numerical simulations showed that the observed detailed behavior of the bound-state energy could also indicate a significant contribution from trivial Andreev bound states Chiu et al. (2017) and therefore the experimental observations cannot be regarded as conclusive evidence of the presence of MBSs on the islands.

Another approach is to consider a three-terminal setup, with two normal local probes and a grounded superconducting probe (e.g., such as sketched in Fig. 1), and a first step in this direction was already taken in Gramich et al. (2017). Linear-response signals of three-terminal devices were also used in the context of the search for signatures of Cooper-pair splitting Hofstetter et al. (2009); Herrmann et al. (2010); Schindele et al. (2012), following theoretical predictions Recher et al. (2001); Loss and Sukhorukov (2000). Further, it was recently pointed out that for wires close to a topological transition, the nonlocal conductance gives information about the induced gap, the topological gap, as well as the coherence length Rosdahl et al. (2018).

In this paper, we investigate this three-terminal setup in detail. From quasiparticle-current conservation we derive a symmetry relation that dictates that the antisymmetric parts of the local and nonlocal conductance are equal for voltages below the gap. Moreover, we show how the (experimentally accessible) nonlocal conductances contain detailed information about the electron and hole components of the bound states in the superconducting region, or more specifically, about the local BCS charge of the bound states, $|u(z)|^{2}-|v(z)|^{2}$, close to the two leads. When compared with predictions from theoretical models, this could help to differentiate in practice between for example near-topological quasi-MBSs and trivial Andreev bound states. In a parallel paper, these findings are investigated experimentally Men .

We first turn to the calculation of the current in the left and right normal leads, to which voltages $V_{L}$ and $V_{R}$ are applied, respectively. Since we are interested in using the differential conductance to probe the subgap states, we will focus exclusively at voltages below the gap: $V_{L,R}<\Delta$, where $\Delta$ is the gap in the grounded superconducting lead, see Fig. 1. This means that no quasiparticles are entering or leaving through the superconducting lead.

We start by using the conservation of probability current, corresponding to unitarity of the scattering matrix, to write the following identities (for more details, see App. A)

where $N_{\alpha}$ is the number of channels in lead $\alpha\in\{L,R\}$, $T_{\bar{\alpha}\alpha}^{e}$ is the total transmission probability of an electron from lead $\alpha$ to the opposite lead $\bar{\alpha}$, $A_{\bar{\alpha}{\alpha}}^{e}$ is the transmission probability of an electron in lead $\alpha$ to a hole in lead $\bar{\alpha}$ (crossed Andreev reflection), and $R^{e}_{\alpha}$ and $A^{e}_{\alpha}$ denote the probability of reflection of an incoming electron in lead $\alpha$ as an electron or hole, respectively. The first equation (1a) expresses the conservation of an incoming electron in lead $\alpha$, while the second equation (1b) expresses that an outgoing electron in lead $\alpha$ must have entered somewhere.

Defining the positive direction of a current always to be into the central scattering region, we can write Takane and Ebisawa (1992)

where $\tilde{f}_{\alpha}(\omega)=f(\omega+eV_{\alpha})-f(\omega)$, with $f(\omega)=1/(e^{\beta\omega}+1)$ the Fermi-Dirac distribution function (where $\beta=1/k_{\rm B}T$ and $\omega$ is measured relative to the equilibrium chemical potential). This equation then allows us to derive the elements of the differential-conductance matrix as

where we will make the important assumption that all transmission and reflection coefficients do not depend on bias voltages, i.e., the voltages only enter through the distribution functions in the leads. In that case, all elements $G_{\alpha\beta}$ only depend on one voltage (corresponding to the second index), so they have the form $G_{\alpha\beta}(V_{\beta})$.

The second conservation law (1b) has the interesting consequence that the antisymmetric parts (in voltage) of $G_{\alpha\alpha}$ and $G_{\alpha\bar{\alpha}}$ are identical. This is easy to see when adding the two, e.g., for $\alpha=L$, and setting $V_{L}=V_{R}=V$

where $f^{\prime}(\omega)$ is the derivative of the Fermi function and

Now inserting $N_{L}-R^{e}_{L}(\omega)$ from Eq. (1b), we can obtain

Due to the general symmetry of the Andreev reflections $A^{e}(\omega)=A^{h}(-\omega)$, we then see that

Defining the symmetric and antisymmetric components of the conductance as

we see that (7) implies that $G^{\rm asym}_{LL}(V)=-G^{\rm asym}_{LR}(V)$, which is one of the main results of the paper, and is illustrated in Fig. 2. We emphasize that $G_{LL}(V)$ and $G_{LR}(V)$ are measured as functions of two different voltages ($V_{L}$ and $V_{R}$, respectively); only when regarded as functions of the same variable $V$ the relation (7) holds. Of course, a similar relation can be derived for $G_{RR}+G_{RL}$, and therefore, if one adds all four elements of the conductance matrix one gets a symmetric function, which is consistent with the fact that a two-terminal NS-junction has a symmetric differential conductance below the gap.

We now study the differential-conductance matrix in more detail, both its symmetric and antisymmetric parts. We do this using a model where the central region has a single Andreev bound state that is coupled to the leads. For metallic wide-bandwidth normal leads, we can use the expression Aleiner et al. (2002) for the $S$-matrix,

where $\nu$ is the density of states in the leads. The central superconducting region is described by a Bogoliubov-de Gennes Hamiltonian for a single level: $H=E_{0}\tau_{z}$, where $\tau_{x,y,z}$ are the Pauli matrices in electron-hole space. The coupling matrix ${\bf W}$ follows as

where $t_{\alpha}u_{\alpha}(-t_{\alpha}v^{*}_{\alpha})$ parameterizes the tunneling coupling between the electron(hole) component of the bound state and an electron state in lead $\alpha$.

The elements of the single-level $S$-matrix can now be inserted into the expressions for the conductance matrix, yielding explicit expressions for the zero-temperature differential-conductance matrix $G_{\alpha\beta}^{0}=(e^{2}/h)g_{\alpha\beta}^{0}$, from which finite-temperature expressions follow straightforwardly by convolution with $-f^{\prime}(\omega)$.

The full results can be found in the App. B, and here we discuss only the limit where the scale of the energy of the bound state $E_{0}$ is much larger than its level broadening. Focusing on the symmetric and antisymmetric components we find close to the resonances, where $\omega\approx\pm E_{0}$,

where $\xi_{\alpha}=\pi\nu|t_{\alpha}|^{2}q_{\alpha}$ and $\gamma_{\alpha}=\pi\nu|t_{\alpha}|^{2}n_{\alpha}$ parameterize the coupling to the leads of the local BCS charges $q_{\alpha}=|u_{\alpha}|^{2}-|v_{\alpha}|^{2}$ and total (local) weights $n_{\alpha}=|u_{\alpha}|^{2}+|v_{\alpha}|^{2}$, respectively. The function

where $\gamma=\gamma_{L}+\gamma_{R}$, is sharply peaked at $\omega=\pm E_{0}$, where $L^{0}(\pm E_{0})=2/\gamma^{2}$, and has a line width of $2\gamma$.

It is interesting to see that $g^{0,\text{sym}}_{LR}$ is proportional to the BCS charge at both terminals, $q_{L}q_{R}$, while $g^{0,\text{asym}}_{LR}$ is proportional only to the BCS charge at the junction where the current is measured, in this case $q_{L}$ (this observation agrees with the rate-equation result derived in Ref. Gramich et al. (2017)). Therefore, the ratio of the peak heights

provides a direct measure for the relative weight of the electron and hole components of the bound state close to end $\alpha$ of the wire. This ratio is closely related to the parameter $\Lambda$ introduced in Ref. Hansen et al., 2018, which can be extracted from the relative heights of neighboring conductance peaks in a Coulomb-blockaded setup Albrecht et al. (2017); Shen et al. (2018) and reveals information about the actual wave function of the bound state.

Our results can easily be extended to include finite temperature, especially when $E_{0}\gg k_{\rm B}T\gg\gamma_{L,R}$. In that case one can approximate the expression given in (15) by $(2\pi/\gamma)\delta(E_{0}-|\omega|)$, and convolution with $-f^{\prime}(\omega)$ then straightforwardly yields

which replaces $L^{0}(\omega)$ in the zero-temperature results. We emphasize that finite temperature affects the line shape of all conductance peaks in the same way, and the relation found in (16) is thus valid at all temperatures, as long as all conductance peaks are well separated.

We thus found a general relation between the local BCS charge of a bound state close to the ends of the wire and the elements of the experimentally accessible differential-conductance matrix. In a sense, these results for an open (strongly coupled) setup complement those of Ref. Hansen et al., 2018 where the wire was treated as a Coulomb-blockaded island with a significant charging energy.

We now investigate the behavior of these BCS charges in more detail for the case where the scattering region is a quasi-one-dimensional semiconducting wire which has, besides proximity-induced superconductivity, strong spin-orbit coupling and a Zeeman splitting that can be made large enough to drive the wire into the topological regime. We thus assume that we can describe the scattering region with the BdG Hamiltonian Oreg et al. (2010); Lutchyn et al. (2010)

acting in the Nambu space $\{u_{\uparrow}(z),u_{\downarrow}(z),v_{\uparrow}(z),v_{\downarrow}(z)\}$, where the Pauli matrices $\sigma$ act in spin space. Here, $m^{*}$ is the effective electronic mass, $\mu$ the chemical potential in the wire, $\alpha$ characterizes the strength of the spin-orbit interaction, $\Delta$ is the proximity-induced pairing potential (assumed real for convenience), and $V_{\rm Z}=\frac{1}{2}g\mu_{\rm B}B$ corresponds to (half) the electronic Zeeman splitting in terms of the $g$-factor $g$ and the applied magnetic field $B$.

We immediately see that, quite generally, the local BCS charge in a bound state described by the Hamiltonian (18) can be related directly to the dependence of the bound-state energy $E_{n}$ on local potentials, which could be controlled through gates such as those sketched in green in Fig. 1. Indeed, if we add a term $V(z,a,b)\tau_{z}$ to the Hamiltonian to account for a local gate potential, where $V(z,a,b)=V_{g}$ for $a<z<b$ and zero otherwise, we obtain straightforwardly

i.e., the slope of the bound-state energy as a function of $V_{g}$ corresponds to the total BCS charge locally at the position of the gate. Similarly, the dependence on a global gate voltage, which effectively controls $\mu$, is connected to the total integrated BCS charge. We thus understand in very simple terms why the BCS charge, as deduced from the differential-conductance matrix, often changes sign at extrema of the conductance-versus-gate voltages traces, such as observed in Men . This of course assumes that the gate voltage only couples to the charge density and does not change other effective parameters such as the spin-orbit coupling strength.

Another common parameter to sweep in experiment is the magnetic field $B$. Assuming that the field predominantly affects the Zeeman splitting as in (18), i.e., neglecting any orbital contributions, we find

which provides a connection between $E_{n}$ and the spin polarization of the BCS charge. We note that in the large-field limit, where the system is strongly spin-polarized, $dE_{n}/dV_{\rm Z}$ approaches the total BCS charge as well.

We now illustrate these findings with numerical examples. First, we calculate the zero-temperature nonlocal conductance $g_{LR}^{0}$ using Eqs. (2,3,9) with a discretized version of the Majorana-wire Hamiltonian (18) to describe the scattering region. We used $N=800$ lattice sites and we set the local coherence factors in (10) equal to the numerically found values for $u$ and $v$ at site 1 (for $L$) and site $N$ (for $R$). We further used $\alpha=0.28$ eVÅ, $\Delta=180~\mu$eV, $m^{*}=0.023\,m_{e}$, and we set the length of the wire to $L=1500$ nm. The resulting intersite hopping energy is $t=471~$meV and we chose $\pi\nu|t_{L}|^{2}=\pi\nu|t_{R}|^{2}=t/2$.

In Figs. 3(a–d) we show the resulting antisymmetric component of the nonlocal conductance, both as a function of $\mu$ (a,c) and $V_{\rm Z}$ (b,d). When we zoom in on the conductance associated with transport through the lowest-energy state (c,d), we already note a correlation between the magnitude and sign of $g_{LR}^{0,{\rm asym}}$ and the apparent slope of the energy of the bound state that is probed, as indeed predicted qualitatively by Eqs. (12,19,20).

To investigate these relations in more detail, we show in Figs. 3(e,f) the local BCS charge $Q_{R}$ as calculated with Eq. (16) using the conductance peak values extracted from the numerical data (green crosses; we omitted regions where $g_{LR}^{0,{\rm asym}}$ is very small, leading to quick divergences due to small numerical inaccuracies). We also calculate the “actual” BCS charges $q_{R}/n_{R}$ of the bound-state wave functions, in an unconnected wire, simply by diagonalizing the wire Hamiltonian (18), and we plot the resulting $Q_{R}$ (black solid lines). There are no relative scaling factors involved, and the agreement with the $Q_{R}$ found from the conductance measurements is very clear. For comparison, we also show the slope of the bound-state energy (in the unconnected wire) and, as expected based on Eqs. (19,20), we observe a stronger correlation with $Q_{R}$ in Fig. 3(e) than in Fig. 3(f), where the spin polarization of the bound state also plays a role. The correlation in Fig. 3(e) is of course also not expected to be perfect since the slope of the energy is related to the integrated BCS charges whereas the conductance measurements probe the local BCS charges at the ends of the wire. In that sense, a comparison of $Q_{R}$ (as found from conductance measurements) and $dE_{0}/d\mu$ provides information about the degree of localization of the bound state close to the ends of the wire: a larger discrepancy implies more relative weight of the bound state in the center of the wire. This could be a valuable tool in the search for and characterization of (quasi-)MBSs.

In conclusion, we have studied a three-terminal device with one of the terminals being a superconducting lead and the two other normal leads. From the general scattering matrix of this system (for quadratic Hamiltonians), we showed that there is a correspondence between the antisymmetric local and nonlocal differential conductances below the superconducting gap. For a single (Andreev bound state) level in the central region there is a furthermore a correspondence between nonlocal conductance and the BCS charges of the bound state (at the terminals). This allows for a study of the electron-hole texture of in-gap bound states and, in particular, it offers a way to test for the signatures characteristic for Majorana bound states with a small overlap, namely that the BCS charge and the energy splitting oscillates out of phase Ben-Shach et al. (2015); Hansen et al. (2018).