Modeling tunneling for the unconventional superconducting proximity effect

In the proximity effect, superconductivity can be produced locally in a normal material placed in close contact with a superconductor [1–4]. Proximity-induced superconductivity in semiconductors has been gaining increasing attention lately [5–15], due to the potential to generate new topologically protected [16–18] excitations and practical applications [19–23]. A few groups have turned their attention to the high critical temperature (T_c) superconducting proximity effect as it may generate novel pairing symmetries, [24] as well as its technical advantages [25–29]. Indeed beyond the large gaps and T_c, the cuprates have a small density of states and linear relation between energy and momentum, making it easier to match boundary conditions with topological insulators [30, 31]. Furthermore, the spatial extent of the proximity effect is set by either the coherence length or the ratio of the density of states and diffusion coefficients of the two materials [32, 33]. Thus the small density of states and the nearly two-dimensional electronic structure of the high T_c superconductors suggests these systems are ready-made for proximity with topological insulators.

Indeed, numerous groups have demonstrated a high T_c proximity effect with insulating cuprates [34–38], Au quantum dots [39], and in LCMO/YBCO thin films [40–43]. Recently we were able to achieve a proximity effect between Bi₂Sr₂CaCu₂O${}_{8+\delta }$ [44] and two doped topological insulators (Bi₂Se₃, Bi₂Te₃), which are known to superconduct with heavy doping [45, 46] or pressure [47]. This observation was enabled by an interpretation of the tunneling spectra, consistent with previous work on s-wave proximity junctions in superconductors. However these results are now controversial since ARPES experiments on thin film heterostructures conflict about the observation of proximity [48–50]. Part of the difficulty lies in the different means by which the junctions were fabricated. In addition there is a need for a framework to model the tunneling spectra of anisotropic superconductors when the proximity effect is present. Recently an approach to junctions with nanoscale disorder at the surface of a d-wave superconductor has emerged [51]. Here we employ a phenomenological approach that combines one used successfully in s-wave proximity devices and the standard approach for calculating tunneling into anisotropic superconductors. Nonetheless, we note there is a need for a more general theory that self-consistently calculates the induced superconductivity.

Blonder, Tinkham, and Klapwijk (BTK) [31] showed that the Bogoliubov equations are very suitable to describe the reflection and transmission of quasiparticles at a normal material/(s-wave) superconductor interface. They assumed that, at the interface, the superconducting gap (${\rm{\Delta }}(x)$) increases instantaneously from zero to a constant value in the superconductor (${{\rm{\Delta }}}_{0}$). The solutions of the Bogoliubov equations in the superconductor and the normal material are found by applying the appropriate boundary conditions. Ultimately the conductance is calculated using these wavefunctions to determine the transmission and reflection probabilities. Kashiwaya et al later modified the BTK theory to take into account the momentum ($\vec{k}$) dependence of the superconducting order parameter often found in unconventional superconductors [52]. Here we extend the calculations of Kashiwaya et al to study the effect of a gradual variation of the magnitude of the pair potential (${\rm{\Delta }}(x)$) near the superconducting-normal interface on the reflection and transmission coefficients of the quasiparticles. Ideally, ${\rm{\Delta }}(x)$ should be calculated self-consistently, including the various possible superconducting symmetries that may arise on the normal material side. Interesting the exact form of superconductivity expected from the proximity effect with topological insulators remains controversial, especially for a d-wave superconductor [53, 54]. Thus we have chosen to focus on the case where the nature of the order parameter does not change across the interface. We note that our approach could be applied to other combinations, by changing the form of the angular dependence of ${\rm{\Delta }}(x)$, however to better explore the role of various parameters, we limit ourselves to a phenomenlogical ${\rm{\Delta }}(x)$ that has worked well for s-wave proximity junctions [1]. We find good agreement with existing differential conductance data obtained from an optimally doped Bi₂Sr₂CaCu₂O${}_{8+x}$/Bi₂Se₃ proximity device (figure 1). The resulting calculation demonstrates an induced and a reduced superconducting gap due to the proximity effect, consistent with previous observations of induced superconductivity via proximity to unconventional superconductors [44]. Furthermore, our model of ${\rm{\Delta }}(x)$ allows us to vary different parameters to see their connection to the spectra observed, enabling us to reveal the origin of various features in the spectra. However we note that this model assumes the normal material already has a pairing potential that would enable the proximity effect, and that the induced superconductivity retains the same momentum dependence as the original superconductor.

Zoom In Zoom Out Reset image size

At the interface between a superconductor and a normal material, the Cooper pairs can leak into the normal material and create a smaller superconducting gap (induced gap ${{\rm{\Delta }}}_{{\rm{i}}}$) locally. As shown in figure 2, this also results in the superconducting gap at the surface being suppressed (reduced gap ${{\rm{\Delta }}}_{{\rm{r}}}$). Some discontinuity in the gap function is allowed at the interface due to the presence of a barrier between the two materials. To simulate a realistic heterostructure, a delta function barrier potential of the form Hδ(x) is assumed at the interface (as is common to do in tunneling calculations [31, 52]). To simplify the formulas, we define the dimensionless parameter, $Z({\rm{\Theta }})=\tfrac{{mH}}{{{\hslash }}^{2}{k}_{{\rm{F}}}{\rm{\cos }}({\rm{\Theta }})}$, where Θ is the angle between the incident quasiparticle and the interface normal ($\hat{x}$). The δ function barrier takes into account the effect of any oxide barrier that may be present at the interface, as well as an effective barrier that might arise from the fermi velocity mismatch between the normal metal and the superconductor [31]. We define the x-direction as the axis perpendicular to the superconductor-normal interface (figure 2), and assume only the magnitude of the pair potential varies in the x-direction. The angular dependence does not change as we go along x, since we expect the pair potential to be uniform in the y–z plane (perpendicular to the x-direction). We note that a key assumption made here is that the in-plane momentum is conserved upon tunneling, and no other orders emerge at the interface/proximity region. While this may not be realistic in interfaces between systems of different symmetry, it is a good first approximation. Furthermore, we anticipate the loss of in-plane momentum conservation to primarily broaden features, which could be included through adding lifetime effects to the model. However the subject of inelastic scattering in such structures is quite complex, [55, 56] though we do not anticipate significant changes in our results. Ultimately, we assume a parabolic x dependence for the shape of the pair potential with x_N and x_S being the length of the induced superconducting region and the length of the reduced superconducting region respectively (figure 2). This approach has been shown to work well for s-wave superconductors, even when other functional forms of ${\rm{\Delta }}(x)$ with gradual variations are used, so long as x_N, x_S, and the values of the induced and reduced gap (${{\rm{\Delta }}}_{{\rm{i}}}$ and ${{\rm{\Delta }}}_{{\rm{i}}}$) are the same [1]. We also note that this phenomenological model could be applied to an S/I/S' junction of two d-wave superconductors connected along the c-axis. Such a situation would require one to have ${T}_{{\rm{c}}}\gt T\gt {T}_{{\rm{c}}}^{\prime }$, ${x}_{N}=\infty ,{x}_{S}=0$, ${{\rm{\Delta }}}_{{\rm{r}}}={{\rm{\Delta }}}_{0}\gt {{\rm{\Delta }}}_{{\rm{i}}}$, however since we are interested in proximity junctions, we do not address this situation here.

Zoom In Zoom Out Reset image size

We use the Bogoliubov equations to describe the quasiparticle excitations in a superconductor [1, 31], determine the wavefunctions throughout the structure and ultimately to determine the probabilities of Andreev (A(E)) and normal reflection (B(E)) which govern the conductance. We note that in principle the problem could be solved analytically. In particular one would begin by considering only incoming electrons and outgoing holes (resulting from Andreev reflection) in the normal region. Then in the induced and reduced regions one would have pairs of electrons and holes, both incoming and outgoing, with their wavefunctions determined by the Bogoliubov equations with a parabolic ${\rm{\Delta }}(x)$. Finally in the superconducting region only outgoing quasiparticles would need to appear. For $x\lt {x}_{N}$ one could use standard Bloch states, and for the $x\gt {x}_{S}$ the solutions found by Kashiwaya et al [52]. Ultimately one would need to apply the continuity conditions at each interface to determine the transmission and reflection probabilities and thus the conductance. However even in a purely s-wave case, such an approach is quite cumbersome, and for the d-wave case even more so. We also note that were such an approach to be undertaken, it would be best to ultimately solve the entire problem self-consistently to determine the gap. Since our goal is to investigate the role of various parameters in the tunneling spectra and provide further evidence for the proximity effect, we focus on a numerical solution described in detail below.

The elementary excitations of a superconductor are made from the creation of an electron with wavevector k and annihilation of an electron with wavevector –k (or creation of a hole with wavevector k). In this formalism, the excitation is represented by a two-element vector ψ [1, 31, 52]: $\psi (\vec{r})=(f(\vec{r}),g(\vec{r}))$. For an anisotropic superconductor the Bogoliubov equations are written as [52]:

$$\begin{eqnarray}{Ef}(r,\hat{k}) & = & \left[-\displaystyle \frac{{{\hslash }}^{2}{{{\rm{\nabla }}}_{r}}^{2}}{2m}-\mu +V(\vec{r})\right]\,f(r,\hat{k})\\ & & +{\rm{\Delta }}(r,\hat{k})g(r,\hat{k})\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}{Eg}(r,\hat{k}) & = & -\left[-\displaystyle \frac{{{\hslash }}^{2}{{{\rm{\nabla }}}_{r}}^{2}}{2m}-\mu +V(\vec{r})\right]g(r,\hat{k})\\ & & +{\rm{\Delta }}(r,\hat{k})f(r,\hat{k}),\end{eqnarray} \tag{ 2 }$$

where $\hat{k}=\vec{k}/| \vec{k}|$, $V(\vec{r})$, $\mu (\vec{r})$, and ${\rm{\Delta }}(\vec{r})$, are the electrical, chemical, and pair potentials experienced by the quasi-particles. For simplicity, we assume the chemical potential is a constant across the interface and take the usual assumption that $H\delta (X)$ accounts for the interface (as described above). In the normal material, far from the interface, ${\rm{\Delta }}=0$ and as expected the equations (1) and (2) reduce to the Schrodinger equation for electrons(holes). To allow for a d-wave superconductor we assume ${\rm{\Delta }}(\hat{k})$ has a momentum-dependence described by ${\rm{\Delta }}(\hat{k})={{\rm{\Delta }}}_{\infty }\cos 2\alpha$, where α is defined as the angle between $\hat{k}$ and the ${\hat{k}}_{y}$-axis in the y, z plane (see figure 2(B)). To simplify the numerical integration, we divide out this oscillation with a ${{\rm{e}}}^{{\rm{i}}\vec{{k}_{{\rm{F}}}}\cdot \vec{r}}$ term and take the trial solutions for f and g: $f(\vec{r},\hat{k})=u(\vec{r},\hat{k}){{\rm{e}}}^{{\rm{i}}{\vec{k}}_{{\rm{F}}}\cdot \vec{r}}$ and $g(\vec{r},\hat{k})=v(\vec{r},\hat{k}){{\rm{e}}}^{{\rm{i}}{\vec{k}}_{{\rm{F}}}\cdot \vec{r}}$. We obtain

$$\begin{eqnarray}&&-{{\rm{i}}{v}}_{{\rm{F}}}\hat{k}\cdot \triangledown u(\vec{r},\hat{k})={Eu}(\vec{r},\hat{k})-{\rm{\Delta }}(\vec{r},\hat{k})v(\vec{r},\hat{k})\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{i}}{v}}_{{\rm{F}}}\hat{k}\cdot \triangledown v(\vec{r},\hat{k})={Ev}(\vec{r},\hat{k})-{\rm{\Delta }}(\vec{r},\hat{k})u(\vec{r},\hat{k}),\end{eqnarray} \tag{ 4 }$$

Ultimately we are interested in calculating the probability of a quasiparticle incident on the interface traversing across, with no incoming quasi-particles on the superconducting side. For convenience it is easier to start with outgoing waves on the superconducting side, and integrating equations (3) and (4) through the intermediate regions (${x}_{N}\lt x\lt {x}_{S}$). Since we are working with an anisotropic superconductor, we have to account for the pair potentials difference for the electron (${k}_{S}^{+}$) and hole (${k}_{S}^{-}$) like quasiparticles: ${{\rm{\Delta }}}_{\pm }={\rm{\Delta }}(\pm {\hat{k}}_{S}^{\pm })=| {{\rm{\Delta }}}_{\pm }| {\rm{\exp }}({\rm{i}}{\phi }_{\pm })$, with ${{\rm{e}}}^{{\rm{i}}{\phi }_{\pm }}=\tfrac{{{\rm{\Delta }}}_{\pm }}{| {{\rm{\Delta }}}_{\pm }| }$ the phases of the pair potentials and ${{\rm{\Omega }}}_{\pm }=\sqrt{{E}^{2}-| {{\rm{\Delta }}}_{\pm }{| }^{2}}$. For a general x we expect the wavefunction to be of the form: ${\psi }_{j}(\vec{r},\hat{k})=({u}_{{aj}}(\vec{r},\hat{k})$, ${v}_{{aj}}(\vec{r},\hat{k})){\rm{\exp }}(-{\rm{i}}{\vec{k}}_{S}^{-}\cdot \vec{r})+({v}_{{bj}}(\vec{r}$, $\hat{k}),{u}_{{bj}}(\vec{r},\hat{k})){\rm{\exp }}({\rm{i}}{\vec{k}}_{S}^{+}\cdot \vec{r})$, where the subscripts a/b account for the left/right moving waves and j = (1, 2) refers to solutions to equations (3) and (4), respectively. In addition, due to the anisotropic nature, we will need to solve the equations for each value of the energy of the quasiparticle (as in BTK), but also for every incident ($\hat{k}\cdot \hat{x}={\rm{\cos }}({\rm{\Theta }})$) and azimuthal (α) angle. Below we first describe how we solve for a given (${\rm{\Theta }},\alpha$), then in the next section we show how to put this together to determine the conductance.

In the superconductor ($x\gt {x}_{S}$) we only have the outgoing waves and thus the two independent solutions:

$$\begin{eqnarray}{\psi }_{1}^{S} & = & ({u}_{0}^{+}(E),{{\rm{e}}}^{-{\rm{i}}{\phi }_{+}}{v}_{0}^{+}(E)){{\rm{e}}}^{{\rm{i}}{\vec{k}}_{S}^{+}\cdot \vec{r}}\\ {\psi }_{2}^{S} & = & ({v}_{0}^{-}(E),{{\rm{e}}}^{-{\rm{i}}{\phi }_{-}}{u}_{0}^{-}(E)){{\rm{e}}}^{-{\rm{i}}{\vec{k}}_{S}^{-}\cdot \vec{r}},\end{eqnarray} \tag{ 5 }$$

where:

$$\begin{eqnarray}&&{u}_{0}^{\pm }(E)=\sqrt{\displaystyle \frac{E+{{\rm{\Omega }}}_{\pm }}{2E}}\,\,{\rm{a}}{\rm{n}}{\rm{d}}\,\,{v}_{0}^{\pm }(E)=\sqrt{\displaystyle \frac{E-{{\rm{\Omega }}}_{\pm }}{2E}\mathrm{.}}\end{eqnarray} \tag{ 6 }$$

Thus equations (5) and (6) form the initial condition for each energy, and set of angles. In the intermediate region, ($-{x}_{N}\lt x\lt {x}_{S}$), the wavefunctions have the general form:

$$\begin{eqnarray}{\psi }_{j}^{l} & = & ({u}_{j}^{+}(E,\vec{{r}}),{{\rm{e}}}^{-{\rm{i}}{\phi }_{+}(\vec{{r}})}{v}_{j}^{+}(E,\vec{{r}})){{\rm{e}}}^{{\rm{i}}{\vec{k}}_{S}^{+}\cdot \vec{{r}}}\\ & & +({v}_{j}^{-}(E,\vec{{r}}),{{\rm{e}}}^{-{\rm{i}}{\phi }_{-}(\vec{{r}})}{u}_{j}^{-}(E,\vec{{r}})){{\rm{e}}}^{-{\rm{i}}{\vec{k}}_{S}^{-}\cdot \vec{{r}}},\end{eqnarray} \tag{ 7 }$$

where $l={\rm{RS}},{\rm{IS}}$ for the reduced and induced superconducting regions respectively, (${u}_{j}^{\pm }(E,\vec{r}),{v}_{j}^{\pm }(E,\vec{r})$) are solutions to equations (3) and (4). At $x={x}_{S}$, we set ${\psi }^{S}({x}_{S})={\psi }^{{\rm{RS}}}({x}_{S})$ and $\tfrac{{\rm{d}}{\psi }^{S}}{{\rm{d}}{x}}{\left|{}_{{x}_{S}}=\tfrac{{\rm{d}}{\psi }^{{\rm{RS}}}}{{\rm{d}}{x}}\right|}_{{x}_{S}}$, neglecting the terms proportional to $\tfrac{{{\rm{d}}{u}}_{j}^{\pm }}{{\rm{d}}{x}}$, $\tfrac{{{\rm{d}}{u}}_{j}^{\pm }}{{\rm{d}}{x}}$, ${\phi }_{\pm }(x)$ and ${\vec{k}}_{S}^{\pm }$. Next we use these conditions to numerically integrate equations (3) and (4) from $x={x}_{S}$ to x = 0. At the boundary between the materials (x = 0), we again match the boundary conditions, however accounting for the potential barrier of height (Z) (i.e. $\tfrac{{\rm{d}}{\psi }^{{\rm{RS}}}}{{\rm{d}}{x}}{\left|{}_{0}-\tfrac{{\rm{d}}{\psi }^{{\rm{IS}}}}{{\rm{d}}{x}}\right|}_{0}=2Z({\rm{\Theta }}){\psi }^{{\rm{RS}}}{| }_{0}$). These then become the initial conditions for the numerical integration of (3) and (4) from the interface (x = 0) to the normal region ($x=-{x}_{N}$). Finally we reach the normal region where the wavefunction has the form:

$$\begin{eqnarray}{\psi }_{j}^{N} & = & (1,0){{\rm{e}}}^{{\rm{i}}{\vec{k}}_{N}^{+}\cdot \vec{r}}+{a}^{+}(E,{\rm{\Theta }},\alpha )(0,1){{\rm{e}}}^{-{\rm{i}}{\vec{k}}_{N}^{-}\cdot \vec{r}}\\ & & +{b}^{+}(E,{\rm{\Theta }},\alpha )(1,0){{\rm{e}}}^{-{\rm{i}}{\vec{k}}_{N}^{+}\cdot \vec{r}}\\ & & +(0,1){{\rm{e}}}^{-{\rm{i}}{\vec{k}}_{N}^{-}\cdot \vec{r}}+{a}^{-}(1,0)(E,{\rm{\Theta }},\alpha ){{\rm{e}}}^{-{\rm{i}}{\vec{k}}_{N}^{+}\cdot \vec{r}}\\ & & +{b}^{-}(E,{\rm{\Theta }},\alpha )(0,1){{\rm{e}}}^{-{\rm{i}}{\vec{k}}_{N}^{-}\cdot \vec{r}},\end{eqnarray} \tag{ 8 }$$

where the coefficients ${a}^{\pm }(E,{\rm{\Theta }},\alpha ),{b}^{\pm }(E,{\rm{\Theta }},\alpha )$ are the amplitudes of the Andreev and ordinarily reflected wave for an incident electron/hole. By applying the continuity conditions at $x={x}_{N}$, we ultimately determine ${a}^{\pm },{b}^{\pm }$.

When an incoming electron is incident on a normal-superconductor interface, it can either ordinarily reflect at the interface, transmit through the interface with a wavevector on the same side of the fermi surface, transmit with a wavevector on the opposite side of the fermi surface, or Andreev-reflect [31]. Andreev reflection involves two incoming electrons crossing the interface in low-barrier superconducting junctions and forming a Cooper pair, typically resulting in a doubling of conductance. Using the conservation of probability, the reflection and transmission coefficients can be related to one another. This leads to a simplified expression for the current at the normal-superconductor interface (I_NS), which depends solely on the probability of Andreev reflection ($A(E,{\rm{\Theta }},\alpha )=| {a}^{+}(E,{\rm{\Theta }},\alpha ){| }^{2}$) and the probability of ordinary reflection ($B(E,{\rm{\Theta }},\alpha )=| {b}^{+}(E,{\rm{\Theta }},\alpha ){| }^{2}$): [31].

$$\begin{eqnarray}{I}_{{\rm{NS}}} & \propto & {\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{-\pi /2}^{\pi /2}{\displaystyle \int }_{-\infty }^{\infty }[f(E-{eV})-f(E)]\\ & & \times [1+A(E,{\rm{\Theta }},\alpha )+B(E,{\rm{\Theta }},\alpha )],\end{eqnarray} \tag{ 9 }$$

where f(E) is the fermi function  $\left(f(E)=\tfrac{1}{1+{{\rm{e}}}^{E/{K}_{B}T}}\right)$, V is the applied voltage on the junction. Ultimately we are interested in the conductance across the interface (${\sigma }_{{\rm{NS}}}$), which can be determined from the current via:

$$\begin{eqnarray}{\sigma }_{{\rm{NS}}}(E,{\rm{\Theta }},\alpha ) & = & \displaystyle \frac{{\rm{d}}{I}}{{\rm{d}}{V}}(E,{\rm{\Theta }},\alpha )=\displaystyle \frac{{{\rm{d}}{I}}_{{\rm{NS}}}}{{\rm{d}}{V}}\\ \sigma (E) & = & \displaystyle \frac{\displaystyle \int {\rm{d}}\alpha {\rm{d}}{\rm{\Theta }}\,{\rm{\cos }}\,{\rm{\Theta }}{\sigma }_{{\rm{NS}}}(E,{\rm{\Theta }},\alpha )}{\displaystyle \int {\rm{d}}\alpha {\rm{d}}{\rm{\Theta }}\,{\rm{\cos }}\,{\rm{\Theta }}{\sigma }_{N}(E,{\rm{\Theta }},\alpha )}\mathrm{.}\end{eqnarray} \tag{ 10 }$$

The normalized conductance ($\sigma (E)$) is calculated numerically. Special care must be taken in dividing up the range for integration as the resulting $a(E,{\rm{\Theta }},\alpha ){\rm{a}}{\rm{n}}{\rm{d}}\,b(E,{\rm{\Theta }},\alpha )$ oscillate strongly for certain values of ${\rm{\Theta }},\alpha$, depending on the initial parameters of $E,Z,{x}_{N},{x}_{S}$.

To check the validity of the calculation we first compared the resulting dI/dV to data from an optimally doped Bi₂Sr₂CaCu₂O${}_{8+x}$ (Bi-2212)/Bi₂Se₃ proximity device (described previously) [44]. The red curve in figure 1 shows a typical differential conductance from these proximity devices at 10 K, revealing a zero bias peak with a width much below the full gap (40 meV) and a second coherence peak also at a value below the Bi-2212 gap. These features were previously attributed to Andreev reflection arising from the induced superconductivity in the normal material, and the resulting reduced gap in the Bi-2212. The overall features are well reproduced by the calculations outlined here as shown by the red line where we used the parameters ${{\rm{\Delta }}}_{{\rm{induced}}}=10\,\mathrm{mV}$, ${{\rm{\Delta }}}_{{\rm{reduced}}}=24\,\mathrm{mV}$, ${{\rm{\Delta }}}_{0}=40\,\mathrm{mV}$, Z = 0.5, x_S = $2.6{\epsilon }_{0}$, x_N = $1.3{\epsilon }_{0}$ (where ${\epsilon }_{0}$ = ${\hslash }{v}_{F}/(\pi {{\rm{\Delta }}}_{0})$ is the BCS coherence length). The mis-match in the fermi velocities are essentially included in the choice of Z, and only the ratio of the sizes of reduced and induced superconducting regions to the coherence length (x_S/${\epsilon }_{0}$ and x_N/${\epsilon }_{0}$) affect the calculation.

Since our interface is perpendicular to the c-axis of the Bi-2212, the hole-like and the electron-like quasiparticles transmitted into the superconductor experience the same effective pair potentials, which have similar dependence on the azimuthal angle α in the ab-plane ${\rm{\Delta }}+={\rm{\Delta }}-={{\rm{\Delta }}}_{0}{\rm{\cos }}(2\alpha )$. The calculated dI/dV is in good agreement with the experimental data (figure 1). We note that small discrepancies likely emerge from the fact that lifetime effects are not included in the calculation, nor are realistic band structures of the two materials, or strong-correlations. This would likely explain the apparent additional broadening in the spectra as well as the electron–hole asymmetry. Nonetheless the calculation clearly captures the qualitative features of the spectra and thus provides a good estimation of the spatial extent and shape of the induced superconducting gap.

In order to study the dependence of the conductance shape on the various parameters employed, we changed only one parameter while holding the others constant. The results are shown in figure 3. First, we start with the simple case of no proximity effect (i.e. no leaking of the order parameter exists between the two materials). Therefore ${{\rm{\Delta }}}_{{\rm{reduced}}}$ = ${{\rm{\Delta }}}_{0}$ = 40 mV and ${{\rm{\Delta }}}_{{\rm{induced}}}$ = 0. Figure 3(a) shows the calculated dI/dV spectra for x_S = $2.6{\epsilon }_{0}$ and x_N = $1.7{\epsilon }_{0}$ at various barrier heights. As expected, in the high barrier limit we obtain a coherence peak originating from the density of states of the superconductor (maximum at ${{\rm{\Delta }}}_{0}$). Upon reducing the barrier, Andreev reflection begins to appear as an enhancement in zero bias conductance which eventually leads to a zero bias peak whose width reflects the full superconducting gap. Beyond confirming our naive expectations, these spectra are similar to the analytical calculation of [52] (figure 2(A) in [52]).

Zoom In Zoom Out Reset image size

Next, we study the Z-dependence of the dI/dV spectra for a proximity-induced gap. Figure 3(B) shows the calculated dI/dV spectra for ${{\rm{\Delta }}}_{{\rm{induced}}}$ = 10 mV, ${{\rm{\Delta }}}_{{\rm{reduced}}}$ = 32 mV, ${{\rm{\Delta }}}_{0}=40\,\mathrm{mV}$, x_S = ${\epsilon }_{0}$, x_N = $2{\epsilon }_{0}$. The barrier height Z is varied from 0.3 to 1. As Z is increased, the peak structures are gradually enhanced, similar to the calculation of conductance spectra in high-barrier tunnel junctions for anisotropic superconductors by Kashiwaya et al (figure 2 in [52]). However here we transfer from a regime where too low of a barrier leads to an Andreev reflection dominated by the original superconductor (i.e zero bias peak extending to the full gap), to a regime where a large barrier enables us to separate the proximity induced Andreev reflection (zero bias peak in blue curve extending to $\approx {{\rm{\Delta }}}_{{\rm{induced}}}$) from tunneling into the full gap of the superconductor (peak in the blue curve at ${{\rm{\Delta }}}_{0}$).

One may wonder why the reduced gap does not reveal any signatures in the spectra, and suggests we may not be capturing its effect on the tunneling. To investigate this further we varied the length of the induced (x_N) and reduced (x_S) superconducting regions. The results are shown in in figures 3(c) and (d), where all parameters are the same as in figure 3(b), however Z is fixed at 3. We find that the detailed spectra are quite sensitive to the choice of the induced and reduced regions, which is consistent with the sensitivity of tunneling to the shape of the wavefunction. For example in figure 3(c), we generally always observe the peak at the reduced gap, however the relative spectral weight of this and the zero bias peak are both quite sensitive to the spatial extent of the induced gap region. We believe this may result from a larger induced region allowing for a higher probability of Andreev tunneling into the induced gap. Interestingly, so far we have not observed a separate peak due to the full gap of Bi-2212 as well as the reduced gap. This is at odds with some of the data seen in junctions of Bi-2212 and Bi₂Te₃ or Bi₂Se₃ [44].

An explanation this discrepancy emerges when the size of the reduced superconducting region is varied below a coherence length. In particular as shown in figure 3(d), when x_S is reduced below ${\epsilon }_{0}$, an additional peak shows up in the calculated spectra (the blue curve in figure 3(d)), corresponding to the bulk superconducting gap. This peak appears due to an additional Andreev reflection process when the energy of the incoming quasiparticles is between the reduced and the bulk superconducting gap (${{\rm{\Delta }}}_{{\rm{reduced}}}\lt E\lt {{\rm{\Delta }}}_{0}$). At larger x_S, the variation of the pair potential at this region is very gradual and no peak is seen at ${{\rm{\Delta }}}_{0}$. Perhaps not surprisingly, varying x_S has a negligible, if any, effect on the width of the zero-bias peak as it is primarily governed by ${{\rm{\Delta }}}_{{\rm{induced}}}$ and Andreev tunneling into the proximity region. We note that the blue spectra in figure 3(d) corresponding to a reduced region half the width of the coherence length is consistent with most of our mechanically bonded Bi-2212/normal material proximity devices [44]. Such a small region over which the superconducting order is reduced is perhaps not surprising given the short coherence length in the Cuprates and the fact that the tunneling occurred along the c-axis. Finally, for a complete consistency check, we have also varied the size of ${{\rm{\Delta }}}_{{\rm{induced}}}$ or ${{\rm{\Delta }}}_{{\rm{reduced}}}$, while keeping all other parameters fixed. As shown in figures 3(e) and (f), the resulting spectra match our expectations. Specifically, lowering the induced gap results in a narrowing of the zero-bias peak, while lowering the reduced gap shifts the high energy peak in the spectrum to lower energies.

Interestingly we have found our approach can capture another subtle aspect of superconductor-normal metal heterostructures, namely coherent geometric resonances. These commonly observed features [2, 57] result from an interference of successive ordinary and Andreev reflections in thin layers, and are termed Tomasch and McMillan-Rowell oscillations, respectively [2, 58, 59]. An example is shown in figure 4(a), where we plot the conductance for ${{\rm{\Delta }}}_{{\rm{induced}}}$ = 12 mV, ${{\rm{\Delta }}}_{{\rm{reduced}}}$ = 28 mV, ${{\rm{\Delta }}}_{0}=40\,\mathrm{mV}$, x_S = $26{\epsilon }_{0}$, x_N = x_S/2, and Z = 0.6. As discussed earlier, when the region over which the superconductivity is reduced is much larger then the coherence length, we expect to only observe features from the induced and reduced gaps. While both features are clearly seen, we also observe additional peaks between the reduced gap and the bulk superconducting gap values. To understand the origin of these features we recall the details of McMillan-Rowell and Tomasch oscillations.

Zoom In Zoom Out Reset image size

McMillan-Rowell oscillations occur due to interference effects in the normal material. Specifically, an electron in the normal material propagating towards the interface at an energy $E\lt {{\rm{\Delta }}}_{S}$ is Andreev-reflected, sending a hole back. The reflected hole has a chance to ordinarily reflect when it hits the normal material edge (x = d_N), return to the superconducting interface and Andreev-reflect again (figure 4(b)). This leads to a series of reflections that interfere with each other, and appear as oscillatory features in the conductance spectra of this tunneling device. Similar to the described bound states in the normal layer, interference effects can occur in the potential well which is formed by the insulating barrier or the superconducting energy gap (Tomasch oscillations) [60]. The oscillations we observe appear between ${{\rm{\Delta }}}_{{\rm{reduced}}}$ and ${{\rm{\Delta }}}_{0}$ (figure 4(a)) and thus are likely coming from an interference effect in the reduced superconducting region. We note that higher values of Z (e.g. see figure 3(b)) did not produce any bound states that can emerge as seen in s-wave superconducting junctions [60]. However this is likely due to the subtle difference between d-wave and s-wave superconductors. Namely the node in the gap means bound-states are typically appearing at or very close to zero-energy [57]. Furthermore such states typically require a change in sign of the superconductor, and thus are seen in ab-plane tunneling into the [110] surface of the cuprates [57, 61]. Indeed previous experiments on high-T_c junctions have only observed McMillan-Rowell and not Tomasch oscillations [62, 63]. However since we are studying c-axis heterostructures we do not expect to see such effects. We note that in superconductor/normal material heterostructures one may also expect to see resonances arising from Andreev Bound states, which are not accounted for in our approach. Such bound states are unlikely to be observed given the large innelastic scattering expected in such junctions. Nonetheless, the observation of geometrical resonances in our dI/dV calculations demonstrates the robustness of this approach. Our calculations not only capture the overall spectral features typically measured in unconventional superconducting proximity devices, but subtle bound states that can emerge as a result of interference effects.

In this paper, we have used a phenomenological approach to compute the differential conductance spectra of a superconductor-normal proximity device for the case of an anisotropic superconductor by solving the Bogoliubov equations. We combine a successful approach in modeling s-wave proximity devices with a standard method for calculating the conductance of interfaces with d-wave superconductors. A parabolic, d-wave pair potential ${\rm{\Delta }}(k,x)$ was assumed for the induced and reduced superconducting regions. A δ function was assumed to include any physical barrier as well any fermi velocity mismatch between the two materials. Under similar conditions the numerical calculation reproduces previous analytic solutions for d-wave superconductors. Furthermore in the proximity regime the calculated conductance spectra reveal a zero-bias peak whose width corresponds to the size of the induced superconducting gap in the normal material. Generally we also observe a peak corresponding to Andreev tunneling into the superconductor's reduced gap at the interface. Upon reducing the size of the reduced gap region to less then a coherence length, we also observe a peak at the full gap of the superconductor, consistent with the expectation that Andreev reflection now reaches the fully gapped region. The calculation is also able to reproduce geometric resonances that are well known to appear in thin heterostructures of normal metal-superconductors. Ultimately, we find that our results are consistent with experimental measurements of high-T_c Bi-2212/Bi₂Se₃ devices and can be used to provide insights in interpreting tunneling into unconventional, superconducting proximity heterostructures. Nonetheless, this phenomenological model does not truly capture the real situation where the order parameter may be different in the induced region. Thus it would be helpful to extend this work by either self-consistently solving the BdG equations, properly modeling the wavefunction in the normal material (for example the Dirac states in topological materials), as well as accounting for a change in the phase of the order parameter across the boundary.

The work at the University of Toronto was supported by the Natural Sciences and Engineering Research Council of Canada, the Canadian Foundation for Innovation, and the Ontario Ministry for Innovation. KSB acknowledges support from the National Science Foundation (grant DMR-1410846). The work at Brookhaven National Laboratory was supported by the Office of Basic Energy Sciences, U S Department of Energy, under Contract No. DE-SC00112704. Z J X was supported by the Center for Emergent Superconductivity, an Energy Frontier Research Center, headquartered at Brookhaven National Laboratory and funded by the U S Department of Energy, under DE-2009-BNL-PM015.