Tunable pseudogaps due to non-local coherent transport in voltage-biased three-terminal Josephson junctions

TL;DR

This paper explores how non-local coherent transport in three-terminal Josephson junctions under voltage bias induces tunable pseudogaps and complex structures in the density of states, revealing new proximity effects beyond two-terminal systems.

Contribution

It demonstrates the existence of highly tunable pseudogaps caused by non-local processes in three-terminal junctions, a phenomenon absent in two-terminal systems, and analyzes their dependence on various parameters.

Findings

Presence of deep, tunable pseudogaps in the density of states.

Identification of non-local coherent transport effects absent in two-terminal junctions.

Similarity of multiple Andreev reflection signatures to two-terminal cases.

Abstract

We investigate the proximity effect in junctions between $N=3$ superconductors under commensurate voltage bias. The bias is chosen to highlight the role of transport processes that exchange multiple Cooper pairs coherently between more than two superconductors. Such non-local processes can be studied in the dc response, where local transport processes do not contribute. We focus on the proximity-induced normal density of states that we investigate in a wide parameter space. We reveal the presence of deep and highly tunable pseudogaps and other rich structures. These are due to a static proximity effect that is absent for $N=2$ and is sensitive to an emergent superconducting phase associated to non-local coherent transport. In comparison with results for $N=2$, we find similarities in the signature peaks of multiple Andreev reflections. We discuss the effect of electron-hole decoherence and of various types of junction asymmetries. Our predictions can be investigated experimentally using tunneling spectroscopy.

Equations

$$E_J({\phi }_1,{\phi }_3)=\sum _{m_1,m_3}{E}_{(m_1,m_3)}{e}^{i(m_1{\phi }_1+m_3{\phi }_3)},$$

$$I_j({\phi }_1,{\phi }_3)=\sum _{m_1,m_3}{I}_{j,(m_1,m_3)}{e}^{i(m_1{\phi }_1+m_3{\phi }_3)},$$

$$I_{\mathrm{j},\mathrm{DC}}({\phi }_Q,V)=I_Q=\sum _m{I}_m(V)\exp (im{\phi }_Q),$$

$$\frac{\partial }{\partial \text{x}}\left(\mathcal{D}(\text{x})\check{G}\frac{\partial }{\partial \text{x}}\check{G}\right)-i\left[\check{H},\check{G}\right]=0,$$

$$\displaystyle\check{G}=\left(\begin{array}[]{cc}G^{R}&G^{K}\\ 0&G^{A}\end{array}\right)\ ,\quad\check{G}^{2}=\check{1},\quad\check{H}=\left(\begin{array}[]{cc}\hat{H}&0\\ 0&\hat{H}\end{array}\right);$$

$$\hat{H}=E{\hat{\sigma }}_z+\frac{1}{2}\Delta (\text{x})(i{\hat{\sigma }}_y+{\hat{\sigma }}_x)+\frac{1}{2}{\Delta }^{\ast }(\text{x})(i{\hat{\sigma }}_y-{\hat{\sigma }}_x).$$

$$\frac{\partial }{\partial \text{x}}\check{j}(\text{x})+\frac{i{e}^2\nu }{\hslash }\left[\check{H},\check{G}\right]=0;\check{j}=-\sigma (\text{x})\check{G}\frac{\partial }{\partial \text{x}}\check{G}.$$

$$I_{ic}=\frac{2e^2}{\pi \hslash }\sum _n\frac{T_n^{(i)}[{\check{G}}_i,{\check{G}}_c]}{4+T_n^{(i)}(\{{\check{G}}_i,{\check{G}}_c\}-2)}.$$

$$I_{fc}=$$

$${\check{G}}_f=$$

$$\sum _i{\check{I}}_{ic}+{\check{I}}_{fc}=0.$$

$$\check{M}=\sum _{i,n}{T}_n^{(i)}\left(\frac{{\check{G}}_i}{1+\frac{T_n^{(i)}}{4}(\{{\check{G}}_i,{\check{G}}_c\}-2)}+{\check{G}}_f\right).$$

$$\displaystyle{G}^{R}_{i}=\frac{1}{\xi}\left(\begin{array}[]{cc}\epsilon&\Delta_{i}\\ -\Delta_{i}^{*}&-\epsilon\end{array}\right)\ ;\quad{G}^{A}_{i}=-\frac{1}{\xi^{*}}\left(\begin{array}[]{cc}\epsilon^{*}&\Delta_{i}\\ -\Delta_{i}^{*}&-\epsilon^{*}\end{array}\right)\ ,$$

$$G^K=(G^R-G^A)\tanh (\beta E/2),$$

$$\check{G_i}(t_1,t_2)=e^{i{\hat{\sigma }}_{z}e{V}_i{t}_1}{\check{G}}_i(t_1-t_2)e^{-i{\hat{\sigma }}_{z}e{V}_i{t}_2}.$$

$$I_{t}=\frac{e^{2}}{2\pi\hbar}T_{\text{tun}}\Tr\Big(\hat{\sigma}_{z}[\check{G}_{t},\check{G}_{c}]^{K}\Big{missing}),$$

Full text

I Introduction

Quantum transport in Josephson junctions has been the focus of extensive research, predominantly studying junctions between $N=2$ superconductors. Short junctions exhibit a strong proximity effect that manifests in equilibrium as an induced minigap in the density of states. A finite minigap is accompanied by a non-dissipative superconducting current. On the contrary, out of equilibrium dynamics due to voltage bias leads to entirely dissipative quasiparticle transport in two-terminal junctions. When the bias voltage $V$ is below the superconducting gap of the leads $\Delta$, the dissipative quasiparticle motion is described by multiple Andreev reflections (MAR) MAR. Electrons and holes cross the structure, being Andreev-reflected at each junction interface. Each crossing provides the energy $eV$, giving rise to features in the $I(V)$ curveScheer at integer fractions of $2\Delta/e$. In the regime dominated by MAR, the density of states no longer manifests a clear minigap, Pierre; Bardas instead exhibiting peaks located at energy intervals separated by $eV$.

Recently, the study of junctions between $N\geq 3$ superconductors has attracted considerable interest, both theoretical Cuevas; Houzet; Freyn; Jonckheere and experimental Pfeffer; Giazotto; Ronen; JeanEude. Unique features have been revealed, that do not manifest in the $N=2$ junctions. In equilibrium, mapping the subgap Andreev spectrum has revealed level crossings at zero energy for non-trivial phase values Akhmerov; Padurariu; Riwar; Giazotto. For $N\geq 4$, the crossing point was shown to have analogous topological properties to Weyl points in topological semi-metals Riwar; Yulinew.

Voltage bias further emphasizes the complex phenomenology of $N\geq 3$ junctions. When the voltages are chosen such that the ratio of any two is a rational number (commensurate bias), the transport is no longer entirely dissipative as is the case in $N=2$ junctions. Previous works have shown that a non-dissipative dc current componentCuevas; Freyn; Jonckheere arises in the junction due to coherent exchange of multiple Cooper pairs non-locally between three or more superconductors. The non-local current is sensitive to bias, as well as an emerging stationary phase that is obtained by combining the phases of multiple superconductors.

The simplest setup consists of the three-terminal Josephson junction (TTJ) where the non-dissipative current is expected to be largest when the two independent phases are affected by opposite voltage bias, $V$ and $-V$, as shown in Fig. 1. Under these conditions the elementary non-local transport process has been termed the quartet process. It corresponds to the exchange of two Cooper pairs, four quasiparticles, between the three superconductors, as shown in Fig. 2. The situation has been recently investigated in Refs. Pfeffer, JeanEude, and Ronen. Two of the first experiments study a diffusive TTJ, where robust transport anomaliesPfeffer and Shapiro stepsJeanEude were observed as a function of two applied voltages $V_{1,3}$, that have been interpreted in terms of three quartet modes. The third experiment studies a phase-coherent TTJ realized in a semiconducting nanowireRonen, showing positive current cross-correlation that are interpreted as evidence of the non-local quartet processes.

Motivated by these recent experiments, in this paper we describe the proximity effect in a short, metallic TTJ under voltage bias, $V$ and $-V$, as shown in Fig. 1. We argue that driving a dc current in terminal $2$, that is assumed at zero voltage, enables the control of the static non-local phase governing the quartet process, $\varphi_{Q}$. We calculate the normal density of states (NDOS) in a wide parameter regime by employing the quantum circuit formulation of the quasiclassical Usadel equation Yuli; Yulibook; Vanevic. For comparison we study the NDOS in the biased two-terminal junction. We reveal the characteristic rich structure of the NDOS originating from MAR, that is similar between two- and three-terminal junctions. We additionally reveal features characteristic only to the three-terminal junction. The most striking of these are the pseudogaps appearing in the NDOS in the regime where coherent non-local processes give rise to bound states. pseudogaps differ from the proximity-induced minigap in that their edges are not as sharp, they do not in all regimes resemble the edges of the bulk gap, and may be less pronounced. What makes pseudogaps unique is the combination of properties: i. they are tunable by the quartet phase, and ii. they depend strongly on voltage bias.

Our study includes the importance of electron-hole decoherence, introduced phenomenologically using the quasiparticle dwell time in the normal region, $\tau_{d}$. Despite describing a short junction on the scale of the coherence length, the dwell time can become appreciable compared to $\hbar/\Delta$ if the contact resistance at the $SN$ interfaces, $R_{b}$, is much larger than the intrinsic resistance of the junction $G_{N}^{-1}$. The Thouless energy Thouless is proportional to the inverse dwell time and can be decreased by a factor $R_{b}G_{N}\gg 1$. For this reason the Thouless energy can become comparable to or smaller than the superconducting gap, $\Delta$, even in short junctions. The magnitude of the proximity-induced minigap in the NDOS is drastically modified by decoherence effects in a large variety of Josephson junctions Dubos; proximity_expt; Gueron; minigap.

We begin our presentation in Section II with a phenomenological description of dynamics in a TTJ under voltage bias. The theoretical method and equations of quasiclassical circuit theory are presented in Section III. Section LABEL:twoterminals discusses the NDOS of a voltage-biased two-terminal Josephson junction, with peaks interpreted in terms of MAR processes. Section LABEL:biasedTTJ discusses the NDOS of a biased TTJ, revealing the signature of MAR processes as well as pseudogaps originating from non-local processes. Section LABEL:conclusions presents our conclusions.

II Phenomenological description

II.1 Local and non-local Josephson effect

The Josephson effect in an $N$-terminal Josephson junction is governed by $N-1$ independent superconducting phase differences. Due to $2\pi$-periodicity, the phase-dependent part of the junction energy can be expanded in harmonics. For $N=3$ we choose the gauge $\varphi_{2}=0$ and express $E_{J}$ as a Fourier series in $\varphi_{1}$ and $\varphi_{3}$,

[TABLE]

where $m_{1}$ and $m_{3}$ are integers running along the entire real axis, and the Fourier coefficients $E_{(m_{1},m_{3})}$ are generally complex energies chosen such that $E_{J}$ is real.

We explore non-local transport by choosing to evaluate the current flowing from terminal $2$ into terminals $j=\{1,3\}$, given by $I_{j}=(2e/\hbar)\partial E_{J}(\varphi_{1},\varphi_{3})/\partial\varphi_{j}$,

[TABLE]

where $I_{j,(m_{1},m_{3})}=(2e/\hbar)im_{j}E_{(m_{1},m_{3})}$. The total current flowing into terminal 2 is obtained from current conservation, $I_{1}+I_{2}+I_{3}=0$. Any possible current flowing from terminal $1$ into terminal $3$ does not modify the discussion.

We define the non-local component of the current flowing from terminal $2$ into terminals $j=\{1,3\}$ by $I_{\rm j,NL}=\partial^{2}I_{j}/\partial\varphi_{1}\partial\varphi_{3}$. The harmonic structure of the Josephson current permits identification of local terms, giving $I_{\rm j,NL}=0$, and non-local terms, giving rise to a finite $I_{\rm j,NL}$. Three contributions correspond to the local Josephson effect between terminals: $1$ and $2$ given by harmonics $(m_{1},0)$; $2$ and $3$ given by harmonics $(0,m_{3})$, and $1$ and $3$ given by harmonics $(-m,m)$, with $-m_{1}=m_{3}=m$. All other pairs of harmonics correspond to the non-local Josephson effect.

The non-local Josephson term lowest in the order of harmonics corresponds to $(m_{1},m_{3})=(1,1)$. It has been named the quartet term, as it implies a coherent exchange of two Cooper pairs, four quasiparticles, between the superconductors as shown in Fig. 2. In the following we show how the quartet term can be filtered from terms corresponding to the rest of the harmonics when driving the junction under commensurate voltage bias, $V_{1}=-V_{3}=V$.

II.2 Out-of-equilibrium dynamics

Under commensurate voltage bias, $V_{1}=-V_{3}=V$, $V_{2}=0$, the phases are given by: $\varphi_{2}=0$, $\varphi_{1}=\varphi_{10}+\omega t$, and $\varphi_{3}=\varphi_{30}-\omega t$, where $\omega=2eV/\hbar$ is the Josephson frequency. The effect of biasing is to separate the harmonics of the Josephson energy in frequency space.

Under these biasing conditions, the quartet term and its harmonics $(m,m)$ give rise to dc current in terminal $j=\{1,3\}$,

[TABLE]

where $\varphi_{Q}=\varphi_{10}+\varphi_{30}$ is the quartet phase and $I_{m}=(2e/\hbar)imE_{(m,m)}$. A detailed discussion of the coefficients $I_{m}$ and their dependence on the bias voltage will be presented elsewhere.

The quartet phase can be tuned independently of the bias voltage by imposing an external current in terminal $2$. Current conservation leads to a current-phase dependence, $I_{DC}(\varphi_{Q})$, similar to the dc Josephson effect, $I_{DC}=-I_{2}=2I_{Q}(V,\varphi_{Q})$. The indirect control of the quartet phase by current bias is similar to the control of the phase drop in a two-terminal Josephson junction by dc current bias. In analogy, the dc current is $2\pi$-periodic in the quartet phase. If $I_{DC}$ surpasses a certain critical value, depending on the details of the junction, the dc behavior of the junction becomes resistive. This situation, together with a discussion of the current flowing between terminals $1$ and $3$, will be presented in detail elsewhere. For discussing the proximity-induced normal density of states (NDOS) in the junction, we will use $V$ and $\varphi_{Q}$ as independent control parameters.

III Microscopic model

We describe transport in a metallic TTJ using quasiclassical equations of non-equilibrium superconductivity. These take the form of a diffusive equation for the quasiclassical Keldysh-Nambu Green’s function, Larkin also known as the Usadel equation (see also Ref. Yulibook)

[TABLE]

In addition to Keldysh-Nambu space (denoted with a check hat, $\check{G}$), the quasiclassical Green’s function generally depends on two times (or energies $\bf E$) and on spatial coordinates $\check{G}(\bf E,\textbf{x})$. The Pauli matrices are defined in Nambu space (denoted with a hat) $\vec{\hat{\sigma}}={\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z}}$, and ${\cal D}(\textbf{x})$ denotes the diffusion coefficient. Matrix products in the Usadel equation are understood as convolutions of the quantities in the double time (or energy) representation, as detailed in the Appendix.

The Usadel equation applies to the most common experimental situation where the junction dimensions are larger than the elastic mean-free path. It is a conservation equation for the Keldysh-Nambu current density, $\check{j}(\textbf{x})$,

[TABLE]

Here, $\nu$ is the electronic NDOS and $\sigma(\textbf{x})$ is the conductivity. The two quantities are related by $\sigma=e^{2}{\cal D}\nu$.

Hereafter we employ a discretized version of the Usadel equation that describes the system in terms of finite quantum circuit elements Yuli; Yulibook. The bulk superconducting terminals $S_{i}$ are described by coordinate-independent Keldysh-Nambu Green’s functions $\check{G}_{i}$. The junction area is represented by a single node described by the unknown Green’s function $\check{G}_{c}$. The node is separated from each terminal $S_{i}$ by a connector that models the transparency of the contact via transmission coefficients $T^{(i)}_{n}$ corresponding to channel $n$ in contact $i$. The Keldysh-Nambu matrix current flowing between terminal $S_{i}$ and the node takes the compact form,Yuli

[TABLE]

The fraction notation for matrix inversion is justified since $\check{G}_{i}$ and $\check{G}_{c}$ commute with $\{\check{G}_{i},\check{G}_{c}\}$.

Decoherence between electrons and holes is accounted for phenomenologically by connecting the node to a fictitious terminal Yuli. In contrast to the other three terminals, that correspond to the superconductors, the Keldysh-Nambu current flowing in the fictitious terminal does not contain particle or energy currents. The Green’s function of the fictitious terminal is chosen such that the corresponding Keldysh-Nambu current describes only the leakage of electron-hole coherence. The Keldysh-Nambu current matrix to the fictitious terminal is given by,

[TABLE]

where $\tau_{d}$ is the dwell time of quasiparticles in the junction, including the connectors. By including $I_{fc}$, the transport equation can be written as a conservation of the current of coherences,

[TABLE]

Since each of the currents are given by a commutation relation between the unknown Green’s function of the central node $\check{G}_{c}$ and a matrix defined by Eqs. (11) and (12), it is convenient to rewrite the current conservation as a commutation relation $[\check{G}_{c},\check{M}]=0$, where the matrix denoted by $\check{M}$ adds up the terms corresponding to the four currents,

[TABLE]

It is important to note that matrix $\check{M}$ depends non-linearly on the unknown Green’s function of the central node $\check{G}_{c}$, as well as on the known Green’s functions of the terminals. The relation $[\check{G}_{c},\check{M}(\check{G}_{c})]=0$ is a non-linear equation to be solved numerically for $\check{G}_{c}$.

III.1 Green’s functions of superconducting terminals

In equilibrium, transport is stationary and the Green’s functions depend on a single energy (or, in time representation, on the difference of the two times and independent of their sum). As a function of energy, the Green’s functions of the superconducting terminals are given by,

[TABLE]

where complex energies have been introduced $\epsilon=E+i0^{+}$ and $\xi=\sqrt{\epsilon+|\Delta|}\sqrt{\epsilon-|\Delta|}$. Here, $\Delta_{i}=|\Delta|e^{i\varphi_{i}}$. The positive, vanishing imaginary part of $\epsilon$ specifies the position with respect to the branch cut of the square root function in the complex plane.

The advanced and retarded Green’s functions are related by $G^{A}=-\hat{\sigma}_{z}\left(G^{R}\right)^{\dagger}\hat{\sigma}_{z}$ and the Keldysh Green’s function $G^{K}$ is obtained from:

[TABLE]

where $\beta=(k_{B}T_{e})^{-1}$ ($T_{e}$ is the temperature).

We consider voltage-biased terminals, $V_{1}=-V_{3}=V,V_{2}=0$. According to the second Josephson relation, $\dot{\varphi}_{i}=2eV_{i}/\hbar$, constant voltage bias gives rise to time-dependent superconducting phase differences that in general give rise to non-stationary transport. As a result, Green’s functions acquire a non-trivial dependence on both energies, or equivalently in time representation, on both the difference, $(t_{1}-t_{2})$, as well as the sum, $(t_{1}+t_{2})$ of the two times. We relate the out-of-equilibrium Green’s function of terminal $i$ to its equilibrium value by the following gauge transformation,

[TABLE]

III.2 Numerical implementation

The theoretical framework outlined so far is sufficiently general to describe out-of-equilibrium transport for arbitrary bias. However, the non-linear equations that determine the unknown Green’s function of the node, $\check{G}_{c}(t_{1},t_{2})$, are very difficult to solve in general. The dependence on two times (or two energies) must be solved on a discrete grid, where each grid point corresponds to an entry of the unknown matrix $\check{G}_{c}(t_{1},t_{2})$ (keeping in mind that each entry is a $4\times 4$ matrix in Keldysh-Nambu space). In the general case the size of matrices involved grows quickly giving rise to an overwhelming computational problem.

To proceed, we use the properties of commensurate bias. In general, transport is governed by two Josephson frequencies corresponding to the two independent voltage differences. For commensurate bias, the two Josephson frequencies are harmonics of a single frequency $\omega_{0}$, the greatest common divisor. For the specific bias $V_{1}=-V_{3}=V,V_{2}=0$, the greatest common divisor is the Josephson frequency $\omega_{0}=2eV/\hbar$. We take advantage of this property by performing a double-time Fourier transform, (detailed in the Appendix) previously used in a different context in Ref. Jonckheere2009. In the transformed representation the Green’s functions depend on a single energy (as in equilibrium) and on the harmonics of $\omega_{0}$ counted by two indices, $\check{G}(E,n,m)$. The definition contains redundancy in the indices, $\check{G}(E,n,m)=\check{G}(E-p\omega_{0},n+p,m+p)$, therefore the Green’s functions are determined by the value in the energy interval $[-\omega_{0}/2,\omega_{0}/2]$. (here we have set $\hbar=1$) An alternative representation with only one harmonic index has been used in Ref. Bezuglyi2 for a two-terminal Josephson junction in the tunnel limit.

The practical numerical implementation involves the truncation of the harmonics by a value $N_{m}$, whereby the Green’s functions are square matrices of dimension $4(2N_{m}+1)$ defined on a one-dimensional grid in the energy interval $[-\omega_{0}/2,\omega_{0}/2]$. The matrix entries decay quickly at large harmonics. It is sufficient to truncate the harmonic expansion at $N_{m}=(2\Delta/eV)$.

We take the following steps to solve for the unknown Green’s function of the central node, $\check{G}_{c}$, iteratively at each energy: (i) we start with a guess value for $\check{G}_{c,n}$, (ii) obtain and diagonalize $\check{M}_{n}(\check{G}_{c,n})$, (iii) obtain a new value, $\check{G}_{c,\text{new}}$ that commutes with $\check{M}_{n}$ and has as eigenvalues the signs of the real part of the eigenvalues of $\check{M}_{n}$, (iv) we check if $\check{G}_{c,\text{new}}$ is within a certain tolerance of $\check{G}_{c,n}$, to ascertain convergence, and finally, (v) if convergence was not achieved, we use a modified matrix as the guess of the next iteration step, $\check{G}_{c,n+1}=||\check{G}_{c,n}+\alpha\check{G}_{c,\text{new}}||$, where $\alpha<1$ is a convergence parameter that must be reduced at energies where the transport depends sharply on energy, and $||...||$ denotes the normalization that ensures $\check{G}_{c,n+1}^{2}=1$. In the calculation, a finite imaginary part is added to the energy, $\varepsilon=E+i\eta$, with $\eta/\Delta\ll 1$, to generate numerically smooth transport resonances. The parameter $\eta$ may be understood as a phenomenological description of weak inelastic effects. Convergence to the $\eta=0+$ limit is especially slow for all superconducting multi-terminal calculations REGIS, requiring small convergence parameters of the order $\alpha\simeq\eta/\Delta$. In the numerical calculation we have used $\eta=0.01$.

III.3 The density of states (NDOS)

The NDOS can be measured using a tunnel probe, as has been already realized for a three terminal junction in equilibrium Giazotto. We model the tunnel probe by adding a normal terminal tunnel coupled to the junction. The current to the tunnel probe is given by,

[TABLE]

where $K$ denotes the Keldysh part of the matrix and $T_{\text{tun}}\ll 1$ describes the coefficient of the tunnel contact. Given that $G_{t}^{R,A}=\pm\hat{\sigma}_{z}$ and $G_{t}^{K}=4$