Diode effect in Josephson junctions with a single magnetic atom – Nature

Experimental details

The Pb(111) crystal was cleaned by several cycles of Ne⁺ sputtering and subsequent annealing under ultrahigh-vacuum conditions. Using an electron-beam evaporator, magnetic adatoms (chromium and manganese) were deposited on the clean substrate held at 30 K. The as-prepared sample was then investigated in a CreaTec STM at 1.3 K. The tungsten tip was coated with a sufficiently thick layer of Pb by dipping it into the crystal surface until a full superconducting gap is observed (Δ_tip = Δ_sample). Differential-conductance spectra at large junction resistance show the quality of the superconductor–superconductor junction by a superconducting gap of size 2Δ_tip + 2Δ_sample = 4Δ around the Fermi level, flanked by a pair of coherence peaks (grey spectra in Fig. 3a,b).

As Josephson spectroscopy is performed at junction conductances of 20 μS or higher, exceptional tip stability is required to withstand the forces acting at these conductances. Smaller indentations are performed to improve the stability and sharpness of the tip. Individual Pb atoms from the tip apex were deposited by controlled approaches to the flat surface. Measurements were then done on individual Pb, Cr or Mn adatoms on the Pb(111) surface. The Cr and Mn atoms were pulled out from the initial adsorption site by advancing with the STM tip³⁸.

Josephson spectroscopy was performed by increasing the current set point at a constant bias voltage of 10 mV until reaching the desired junction conductance. After tip stabilization, a large series resistor R_series = 1 MΩ was introduced into the bias line. This resistance is sufficiently large compared with the junction resistances, so that the junction is effectively current biased. Current-biased Josephson spectroscopy was then performed by sweeping the current bias back and forth between positive and negative values at a rate of 100 nA s⁻¹ to 320 nA s⁻¹. Datasets with the same tip were recorded at the same ramp rate for direct comparison of magnetic and non-magnetic (Pb) adatoms. Small variations in the ramp rate do not lead to notable changes of the switching and retrapping current. This is in agreement with their logarithmic dependence on the ramp rate. Positive current corresponds to tunnelling of electrons from tip to sample. For statistical analysis, we perform between 500 and 2,000 sweeps in each direction. The STM feedback was turned off during measurement.

We analysed the data using a dedicated Python program. Switching and retrapping events were extracted by taking the derivative of the individual V–I curves, which were previously smoothed by a standard Gaussian routine. We also determined G_PD from the slope of the V–I curve in the trapped state. In analysing the data, we account for several instrumental effects. (1) A slow creep of the piezoelectric elements causes the tip to drift towards the surface, effectively changing the junction conductance. We continuously monitor G_PD to characterize the junction and plot all switching and retrapping currents versus G_PD. (2) The differential amplifier used during the Josephson measurements introduces a slowly shifting voltage offset, which we subtract from the individual V–I curves. (3) The voltage/current source has a small offset. For this reason, we correct the entire dataset, including the data measured on the magnetic adatoms, by the mean offset for all recorded data on the pristine Pb–Pb junction under identical measurement conditions, that is, identical tip and identical tip locations. (4) At high junction conductances, the voltage drop across the series resistance of the external circuit becomes non-negligible in the voltage-biased measurements. We correct for this by calibrating the voltage to the superconducting gap size of the Pb–Pb junction.

Statistical analysis of switching and retrapping currents

As described in the main text, we create a Josephson junction by advancing the STM tip to the surface at a bias voltage (10 mV) far above the superconducting energy gap until the desired normal-state junction conductance (a few tens of μS) is reached. We effectively current bias the junction by inserting a large series resistor (R_series = 1 MΩ) into the bias line and sweep the current (a few nA) in both directions. The transition from the resistive to the low-resistance state (I_re) is seen as a sudden drop in the voltage, whereas switching from the low-dissipation to the dissipative branch (I_sw) occurs as a sudden increase in the voltage. Both events are stochastic in nature, owing to Johnson–Nyquist noise. For this reason, we complement single sweeps by histograms of switching and retrapping currents extracted from a larger set of V–I curves. Non-reciprocity of the switching and retrapping current is then seen as asymmetries between the histograms for positive and negative bias. Extended Data Fig. 1 shows corresponding histograms extracted from 500 to 2,000 sweeps recorded on Pb, Cr and Mn junctions with G_N equal to 50 μS. For the Pb junction, the histograms of the switching currents |I_sw,+| and |I_sw,−| exhibit broad Gaussian-like distributions, with the same average ((5.9 ± 0.4) nA) for both bias directions. The histograms of the retrapping currents |I_re,+| and |I_re,−| are narrower ((1.8 ± 0.1) nA) but also independent of bias direction. The histograms for Cr and Mn junctions are qualitatively different. The histograms of the retrapping currents exhibit a clear relative shift between the two bias directions, leading to different absolute values of the averages of I_re,+ ((1.9 ± 0.2) nA for Cr and (1.86 ± 0.04) nA for Mn) and I_re,− ((−1.4 ± 0.2) nA for Cr and (−2.18 ± 0.06) nA for Mn). The histograms of the switching current exhibit a noticeable but much weaker dependence on the bias direction.

Analysis of switching and retrapping currents as a function of G
_PD

The histograms in Extended Data Fig. 1 reflect the stochastic nature of the switching and retrapping processes but are further broadened by piezoelectric creep over the course of the measurement. The creep effectively increases the junction conductance (as quantified by the phase-diffusion conductance G_PD) with time. For each, we minimize the creep-induced broadening by using 100 consecutive sweeps for separate histograms with an associated average G_PD. Extended Data Figs. 2–4 illustrate this analysis. The histogram for the earliest 100 sweeps are shown at the bottom of each panel. Histograms obtained from subsequent batches of 100 sweeps correspond to larger junction conductances G_PD, as indicated in the figures. This increase is accompanied by an increase in |I_sw| and |I_re| as seen by a shift of the corresponding histograms. This scheme is the basis for Fig. 2, which collects the average retrapping currents, along with the standard deviations of all of these histograms.

Comparison of switching currents

Extended Data Fig. 5 shows the switching currents of Cr and Mn junctions as a function of G_PD, in both cases compared with Pb junctions measured with the same tip. For identical tips, the switching currents |I_sw| show almost the same linear dependence on G_PD, justifying the use of G_PD as a suitable measure of junction conductance.

Influence of the STM tip

The STM tip is an integral part of our atomic-scale Josephson junctions. To ensure that the main findings remain valid independent of details of the tip apex, we investigated several tips obtained through reshaping by large tip indentations into the Pb substrate. Extended Data Fig. 6 shows the non-reciprocity of the retrapping current as a function of G_PD for junctions including Cr and Mn atoms but measured with different tips. All tips show a positive value of the asymmetry ΔI_re = |I_re,+| − |I_re,−| in case of Cr, a negative value for Mn and no asymmetry for Pb adatoms. Although these qualitative observations are robust for all tips, there are small differences in the magnitude of the asymmetry at the same value of G_PD. We tentatively ascribe these variations to tip-dependent Josephson coupling energies and quasiparticle currents, as well as noise levels.

YSR states and particle–hole symmetry

We trace the asymmetric quasiparticle current to YSR states, which the magnetic impurity induces within the superconducting gap of the substrate superconductor. Here we briefly explain the origin of this asymmetry and its relation to breaking of particle–hole symmetry. The interaction of the magnetic impurity (impurity spin S) with the conduction electrons of the substrate (field operators ψ_k,σ with wavevector k and spin σ) takes the form

Here ({bf{s}}=frac{1}{2}{boldsymbol{sigma }}) with the vector of Pauli matrices σ, J denotes the strength of the exchange coupling and K the strength of potential scattering. Focusing for simplicity on ({rm{spin-}}frac{1}{2}) impurities, this interaction can be obtained from the Anderson impurity model (impurity level with energy ϵ < 0, on-site interaction U > 0, hybridization t) by a Schrieffer–Wolff transformation⁴³. This yields

In the particle–hole symmetric situation ({epsilon }=-,frac{U}{2}), the excitation energies are identical for the empty configuration (|ϵ|) and the doubly occupied configuration (ϵ + U). In this case, the potential scattering K vanishes. Potential scattering becomes non-zero when the empty and doubly occupied configurations have different excitation energies, reflecting broken particle–hole symmetry. The sign of K depends on which of the configurations has the higher excitation energy.

A standard calculation^44,45,46,47 shows that the YSR state induced by the magnetic adatom induces a pair of subgap states whenever the exchange coupling J is non-zero. For K = 0, the electron wavefunction u and the hole wavefunction v of the bound state are equal to each other, as expected in a particle–hole symmetric situation. An asymmetry between the electron and hole wavefunctions appears when potential scattering is non-zero. The sign of the asymmetry depends on the sign of K.

This asymmetry between the electron and hole wavefunctions explains the asymmetric current–voltage characteristic of junctions including a magnetic impurity. Single-electron tunnelling at positive (negative) bias is proportional to |u|² (|v|²), in which the wavefunctions are evaluated at the tip position (Fermi’s golden rule). Similarly, several Andreev reflections will also involve these factors when the final state involves an excitation of the subgap YSR state. As particle–hole symmetry requires fine tuning of the impurity parameters, one generically expects that the current–voltage characteristics of junctions involving a magnetic adatom are asymmetric. The direction of the asymmetry depends on the details of the atomic physics of the adatom, consistent with our observation of opposite asymmetries for Mn and Cr.

We note that these considerations are independent of whether time-reversal symmetry is broken or not. The asymmetry of the YSR wavefunctions is controlled by potential scattering—and thus only by breaking of particle–hole symmetry—even when the adatom spin is polarized and time-reversal symmetry is explicitly broken^44,45,46,47. Of course, broken time-reversal symmetry may lead to asymmetries in the current–phase relation as well, which might induce a coexisting non-reciprocity of the switching current.

In the absence of thermal fluctuations, the switching and retrapping currents can be obtained from the junction dynamics as follows. In the absence of fluctuations, the junction switches at the critical current, that is, at the current bias at which the tilt of the washboard potential eliminates the minima. In this limit, an asymmetry in the switching current requires an asymmetric washboard potential or, equivalently, an asymmetric current–phase relation. Fluctuations will then reduce the switching current below the critical current but, in the limit of weak damping (pronounced hysteresis), the asymmetry of the switching current is largely inherited from the asymmetry in the critical current.

On the other hand, the retrapping current is the result of different physics. In the absence of fluctuations and at weak damping, the junction retraps, once the energy gain owing to the bias current (that is, owing to the tilt in the language of the washboard potential) becomes smaller than the frictional energy loss during the motion. The energy gain depends on the tilt but not on the shape (or the asymmetry) of the washboard potential. Thus, an asymmetry can only arise from asymmetries in the frictional energy loss, which is associated with the quasiparticle current at the microscopic level. Fluctuations tend to increase the retrapping current but the asymmetry is essentially inherited from the asymmetry in the retrapping currents of the junction in the absence of fluctuations.

Simulations of the RCSJ model

Our theoretical simulations underlying Fig. 4 are based on the RCSJ model for a current-biased junction^28,29,

Here I_bias is the current bias, V the voltage drop at the junction, C the junction capacitance and φ the phase difference across the junction. We assume a symmetric and sinusoidal current–phase relation ({I}_{{rm{s}}}(phi )={I}_{{rm{c}}}sin phi ). We allow for a general nonlinear dissipative current I_d(V), with associated Nyquist noise δI with correlator (leftlangle delta I(t)delta I({t}^{{prime} })rightrangle propto delta (t-{t}^{{prime} })) (see below). When combined with the Josephson relation V = (ħ/2e)dφ/dt, equation (4) gives a Langevin equation for the phase difference across the junction. We solve the Langevin equation by Monte Carlo integration, accounting for the current sweep, to obtain the results shown in Fig. 4 (with further details in Extended Data Fig. 7), as well as in Extended Data Figs. 8 and 9.

The dissipative current I_d(V) includes the quasiparticle current I_qp(V), which we extract from experimental I–V traces (see below for details). To account for the observed phase diffusion in the trapped state, we also incorporate frequency-dependent friction. Following Kautz and Martinis³⁰, we shunt the junction by an extra RC element with ohmic resistor (widetilde{R}) and capacitor (widetilde{C}) to model dissipation induced by the electromagnetic environment. The total dissipative current is then the sum of the quasiparticle current and the current flowing by means of the RC element,

in which (widetilde{V}) is the voltage drop across the capacitor, which satisfies the equation

The RC element is inconsequential at low frequencies (running state), so that damping is dominated by the quasiparticle current. By contrast, it dominates friction at high frequencies (trapped state), allowing for phase diffusion. We assume (V,/widetilde{R}gg {I}_{{rm{qp}}}(V)), so that the quasiparticle current is effectively shorted at high frequencies, ({I}_{{rm{d}}}(V)simeq V,/widetilde{R}). The Nyquist noise associated with the quasiparticle current has correlator (leftlangle delta {I}_{{rm{qp}}}(t)delta {I}_{{rm{qp}}}({t}^{{prime} })rightrangle =2Tleft[{I}_{{rm{qp}}}(V)/Vright]delta (t-{t}^{{prime} })), whereas the Nyquist noise associated with the resistor (widetilde{R}) has correlator (leftlangle delta {I}_{widetilde{R}}(t)delta {I}_{widetilde{R}},({t}^{{prime} })rightrangle =2T{widetilde{R}}^{-1}delta (t-{t}^{{prime} })).

Measuring time in units of the inverse plasma frequency, τ = ω_pt with ({omega }_{{rm{p}}}={left[2e{I}_{{rm{c}}}/hbar Cright]}^{1/2}), and currents in units of the critical current, i = I/I_c, the resulting RCSJ equations become

in which we defined the dimensionless voltages v = 2eV/ħω_p and (widetilde{v}=2ewidetilde{V}/hbar {omega }_{{rm{p}}}), the dimensionless currents i_b = I_bias/I_c, ({i}_{{rm{s}}}(phi )={I}_{{rm{s}}}(phi )/{I}_{{rm{c}}}=sin phi ) and i_qp(v) = I_qp(ħω_pv/2e)/I_c, the effective quality factor (widetilde{Q}=widetilde{R}C{omega }_{{rm{p}}}) at large frequencies, as well as the reduced temperatures θ = T/E_J and (widetilde{theta }=widetilde{T}/{E}_{{rm{J}}}). (Here E_J = ħI_c/2e is the Josephson energy and (widetilde{T}) is the temperature of the resistor (widetilde{R})). We also defined dimensionless Langevin currents ξ₁ and ξ₂ with normalized correlations (langle {xi }_{i}(tau ){xi }_{j}({tau }^{{prime} })rangle ={delta }_{ij}delta (tau -{tau }^{{prime} })) corresponding to δI_qp and (delta {I}_{widetilde{R}}), respectively. We estimate the experimental parameters as R_N ≈ 20 kΩ, Δ ≈ 1.5 meV, T ≈ 0.1 meV and C ≈ 10⁻¹⁵ F. This gives I_c ≈ 100 nA, E_J ≈ 0.2 meV and ħω_p ≈ 0.3 meV. The reduced temperature is thus θ = 0.5. For the RC element, we choose parameters (widetilde{Q}=10), (widetilde{tau }=1,000) and (widetilde{theta }=theta ). We sweep the bias current with a rate dI_bias/dt = 10⁻⁷I_cω_p ≈ 1 nA μs⁻¹. The experimental sweep rate is smaller by a factor of about 10⁻³ but this would make the numerical simulations forbidding. Along with the simplified current–phase relation and the order-of-magnitude estimates of experimental parameters, this implies that one can only expect qualitative but not quantitative agreement between simulations and experiment.

Theoretical simulations of single traces are shown in Fig. 4. The corresponding close-up view of the absolute values of the currents in Extended Data Fig. 7 brings out the large asymmetry of the retrapping current and the smaller (and, according to the histograms, largely statistical) asymmetry of the switching current. In Extended Data Fig. 8, we show the histograms of the absolute values of switching and retrapping currents extracted from 100 sweeps in each current direction. Note that the panels only differ in the precise form of I_qp(V), which is extracted from the I–V curves of Pb, Cr and Mn, respectively. The simulations based on the I_qp(V) of Pb do not show asymmetry in the switching or the retrapping currents. The simulations based on the I_qp(V) of Cr and Mn exhibit weak asymmetry in the switching currents and strong asymmetry in the retrapping currents, correctly reproducing the qualitative features of the experimental histograms in Extended Data Figs. 2–4.

Asymmetric current–phase relation

To rule out the possibility that the observed asymmetry stems from the current–phase relation I_s(φ) rather than from the dissipative quasiparticle current, we now demonstrate that an asymmetric current–phase relation leads to strong asymmetry in the switching currents and weak asymmetry in the retrapping currents, in contrast to our experimental observations. To this end, we simulate equation (7) using Pb I–V data for I_qp(V) together with an asymmetric current–phase relation

We choose φ₀ = 0.5 = b and fix I₀ ≃ 54.2 nA by requiring that the current entering the definition of the plasma frequency, that is, the slope of I_s around the stable minimum, is still 100 nA (which we continue to use as the unit of current). The critical current now depends on direction, with I_c,+ ≃ 53.3 nA and I_c,− ≃ 80.0 nA. Histograms of switching and retrapping currents obtained by simulating equation (7) with the current–phase relation given in equation (8) are presented in Extended Data Fig. 9. The asymmetry of the switching currents is clearly much greater than that of the retrapping currents. Thus, a symmetric dissipative current together with an asymmetric current–phase relation cannot explain the phenomenology of strongly asymmetric retrapping currents and weakly asymmetric switching currents observed for the Cr and Mn Josephson junctions.

Extraction of quasiparticle current

We extract the quasiparticle contribution to the dissipative current I_qp(V) from voltage-biased measurements of Pb, Cr and Mn junctions at the normal-state conductance of G_N = 50 μS (see Fig. 3f). As well as the quasiparticle current, these traces include a Josephson peak originating from incoherent Cooper-pair tunnelling. We remove the Josephson contribution I_J(V) by fitting to the phenomenological expressions⁴⁸

over a voltage range e|V| ≪ Δ, which contains the Josephson peak. (We chose e|V| < 0.32 meV). We also account for offsets in the measured voltage and current through the parameters V_offset and I_offset. The fit parameters are collected in Extended Data Table 1. We then subtract the Josephson contribution as well as the offsets from the measured data to isolate the quasiparticle contribution. To reduce the fluctuations at small V associated with the Josephson contribution, a Gaussian filter (width σ = 5 data points ≃ 0.55 mV) is applied to the isolated quasiparticle current data. Finally, I_qp(V) is obtained by interpolation using a linear splining procedure, enforcing I_qp(0) = 0.