Role of the quasi-particles in an electric circuit with Josephson junctions

TL;DR

This paper develops a comprehensive computational framework that models the microscopic dynamics of Josephson junctions, including quasi-particles, to better understand their behavior in superconducting circuits and predict non-equilibrium phenomena.

Contribution

It introduces a generalized model that incorporates quasi-particle effects and complex circuit geometries, extending the standard RCJ model for more accurate simulations.

Findings

Out-of-equilibrium current-phase relations show high harmonic AC Josephson effect.

Quasi-particle dynamics significantly alter resonance shapes in RLC circuits.

The approach enables quantitative modeling of superconducting circuits with complex junctions.

Abstract

While Josephson junctions can be viewed as highly non-linear impedances for superconducting quantum technologies, they also possess internal dynamics that may strongly affect their behavior. Here, we construct a computational framework that includes a microscopic description of the junction (full fledged treatment of both the superconducting condensate and the quasi-particles) in presence of a surrounding electrical circuit. Our approach generalizes the standard ResistorCapacitor-Josephson model (RCJ) to arbitrary junctions (including e.g. multi-terminal geometries and/or junctions that embed topological or magnetic elements) and arbitrary electric circuits treated at the classical level. By treating the superconducting condensate and quasi-particles on equal footings, we capture non-equilibrium phenomena such as Multiple Andreev Reflection. We show that the interplay between the quasi-particle dynamics and the electrical environment leads to the emergence of new phenomena. In a RC circuit connected to single channel Josephson junction, we find out-of-equilibrium current-phase relations that are strongly distorted with respect to the (almost sinusoidal) equilibrium one, revealing the presence of high harmonic AC Josephson effect. In an RLC circuit connected to a junction, we find that the shape of the resonance is strongly modified by the quasi-particle dynamics: close to resonance, the current can be smaller than without the resonator. Our approach provides a route for the quantitative modeling of superconducting based circuits.

Figures

Equations

$$i\hslash {\partial }_t\hat{\rho }=[\hat{H}(\phi ,t),\hat{\rho }],$$

$$I(t)=\mathrm{Tr}(\hat{I}\hat{\rho }),$$

$$\frac{d^2\phi }{d{t}^2}=F\left(\phi ,\frac{d\phi }{dt},I(t)\right),$$

$$\hat{H}$$

$$+\sum _{x=-\infty }^{+\infty }\Delta (1-{\delta }_{x,0}){\hat{c}}_{x\uparrow }^{†}{\hat{c}}_{x\downarrow }^{†}+h.c.$$

$$\frac{d^2\phi }{d{t}^2}+\frac{1}{Q}\frac{d\phi }{dt}+\sin (\phi )=\frac{I_0}{I_c},$$

$$I(\phi )=I_c\sin \phi +I_{\mathrm{MAR}}(2e\dot{\phi }/\hslash )$$

Full text

There is currently a huge effort around the world - both within academia and major industrial partners - to promote the superconducting transmon quantum bit Koch et al. (2007); Barends et al. (2013); Paik et al. (2011) from a laboratory object to a viable technology for building a quantum computer. The central element of this approach is a weak normal link between two pieces of superconductors, the Josephson junction. Although tunneling junctions with an insulating (oxide) barrier are the most mature elements, other types of junctions such as atomic contacts Goffman et al. (2000) (with very few propagating channels), semiconducting nanowires Mourik et al. (2012) (with high spin-orbit suitable for stabilizing Majorana bound states) superconducting-ferromagnetic-superconducting Izyumov et al. (2002); Buzdin (2005) (with anomalous current phase relations) or multiterminal devices Riwar et al. (2016) could provide new functionalities to the superconducting toolbox. While the theoretical description of these objects is rather well understood Golubov et al. (2004), many relevant regimes lie outside of what may be treated analytically and the development of numerical methods is important. In fact, the complexity of the circuits that are being created is increasing very rapidly and building predictive numerical tools is a key element for the success of any quantum technology.

Two very successful complementary viewpoints are commonly used to describe Josephson junctions circuits. The first one is the RCJ model Likharev (1979); Stewart (1968) (Resistor-Capacitor-Josephson) that views the Josephson junction as a highly non-linear impedance embedded in an electric circuit. In this model, one considers a classical circuit such as the ones shown in the insets of Fig. 1 or Fig. 3 with resistances $V=RI$, capacitances $I=C\partial_{t}V$, inductances $V=L\partial_{t}I$ or other classical elements. The Josephson junction is described by its current-phase relation $I=I_{c}\sin\varphi$ and the Josephson relation $\partial_{t}\varphi=(2e/\hbar)V$. Such a simple model is surprisingly powerful. It captures the hysteresis loops of the $I-V$ curves. Its simple extension, where one adds a Langevin stochastic term to account for finite temperature, accurately describes the noise properties found experimentally including the probability for the junction to switch from the superconducting branch Ben-Jacob et al. (1982). It has also been successfully used for more elaborate circuits that include resonators Cassidy et al. (2017). Its quantum extension provides the model used to design the various sorts of superconducting qubits Xiang et al. (2013) and has been shown to describe very accurately a large corpus of experimental data Devoret et al. (1990). Yet, the model fails dramatically in some very simple limits. For instance, at large voltages, it does not properly reproduce the Ohmic behavior of the circuit, since the latter involves the excitation spectrum of the junction which is not accounted for in the current-phase relation. More importantly, it does not account for some important out-of-equilibrium phenomena such as Multiple Andreev Reflexion Klapwijk et al. (1982); Averin and Bardas (1995) (MAR) processes.

The second model uses a microscopic mean-field description of the junction through the (time-dependent) Bogoliubov-De-Gennes (BdG) equation. BdG models capture most of the salient features of these junctions including those which contain exotic (e.g. topological or magnetic) materials. It naturally describes MAR Averin and Bardas (1995), the interplay with microwaves Cuevas et al. (2002), ac Josephson effects and emergent topological effects in multi-terminal geometries Houzet et al. (2013). Until recently, however, the direct numerical integration of BdG equations has been very limited due to the intrinsic computational complexity Weston and Waintal (2016a) and did not include the electromagnetic environment of the junction.

The present letter builds on recent advances made in time-dependent computational transport Kloss et al. (2018) to construct a numerical method that merges the RCJ model with the BdG equation, thereby providing a fully self-consistent treatment of the Josephson junction and its electromagnetic environment at the BdG level (hereafter called RC-BdG model). The method has an arbitrary precision and is scalable to hundreds of thousands of orbitals, paving the way to the simulations of complex superconducting circuits. It applies to arbitrary BdG equations and classical electromagnetic environments.

Problem formulation. We model our circuits in two parts. First, the junction itself is described with a microscopic BdG Hamiltonian $\hat{H}(\varphi,t)$ that depends explicitly on time $t$ (through e.g. a capacitive gate) and on the phase difference $\varphi(t)$ between the two superconductors (which extends to a vector when more than two superconductors are involved). Note that due to the AC Josephson effect, the problem is intrinsically time-dependent even in the absence of time-dependent perturbations. Integrating the BdG equation provides the density matrix $\hat{\rho}(t)$ from which one can compute the current $I(t)$ that flows through the system. Below, we restrict ourself to the average current, but its quantum fluctuations are also accessible through our formalism Gaury and Waintal (2016). The second part of the model describes the classical circuit, or the electromagnetic environment, that surrounds the junction. The classical equations that describe these circuits take the form of a differential equation for $\varphi$. More complex circuits are described in a similar way with more degrees of freedom describing the classical circuit. The set of equations reads,

[TABLE]

where the function $F\left(\varphi,d\varphi/dt,I(t)\right)$ describes the dynamics of the classical circuit (RC or RLC equation in the examples below) with $I(t)$ as an external source term. We numerically solve the BdG equation within the Keldysh formalism using the approach developed in Gaury et al. (2014); Weston and Waintal (2016b) where the problem is unfolded onto a set of Schrödinger equations for scattering wave functions. An efficient algorithm has been constructed in Weston and Waintal (2016a) to integrate the corresponding equations. The corresponding software, “t-Kwant” relies on the Kwant package Groth et al. (2014) and will be released as open source in a near future. Eq.(1a) alone amounts to solving a few hundred time-dependent Schrödinger equations (the actual number depending on the required energy resolution). The self-consistent condition Eq.(1b) and (1c) makes the problem significantly more challenging since it creates non-linear couplings between these Schrödinger equations. Following Kloss et al. (2018), we address this non-linear coupling by taking advantage of the separation of time scales in the problem between the microscopic time scales of the BdG equation (which imposes discretized time steps of lengths much smaller than $\hbar/E_{F}$ with $E_{F}$ the Fermi energy) and the evolution of the electromagnetic variables $I(t)$ and $\varphi(t)$ that takes place on much slower time scales (typically GHz frequencies as compared with PHz for the Fermi energy in actual devices). Hence, we use a doubly adaptive predictor-corrector approach for $\varphi(t)$ as explained in Kloss et al. (2018): Eq.(1a) is integrated with a “predicted” function $\varphi(t)$ and Eq.(1b) and (1c) are used on the larger time scale to construct this prediction. A straightforward time adaptive fourth order Runge-Kutta Hairer et al. (1993) is used for the integration of Eq.(1c). We used the algorithm of Istas et al. (2018) to calculate the (Andreev) bound states of the model with precision. The method is general to any Hamiltonian $\hat{H}$ that is quadratic in creation/destruction operators. Hence, it may handle electron-electron interaction effects at the mean field or Random Phase Approximation levelKloss et al. (2018) but does not capture correlation effects. Static or dynamic disorder may be added directlyGaury et al. (2014). The full code used for generating the data of this article can be found in Rossignol et al. (2019)

To be specific, we now turn to a particular BdG Hamiltonian that describes a single channel junction. The BdG Hamiltonian describes a one dimensional systems with the two superconductor corresponding to $x<0$ and $x>0$ while the normal region is formed by a single site at $x=0$ placing the system in the short junction limit,

[TABLE]

Here $\varphi(t)=(e/\hbar)\int^{t}V(t^{\prime})dt^{\prime}$ where $V(t)$ is the voltage difference across the junction, $\Delta$ is the superconducting gap inside the superconductors, $U$ is a potential barrier used to tune the transmission probability $D$ of the junction. In the calculations below we used $E_{F}=2$, $\Delta=0.1$ (which will be used as our unit of energy), and $U=2$ which corresponds to a junction with intermediate transmission $D=0.5$. For this value, the equilibrium current-phase relation has small deviations with respect to a sinusoidal shape but the $I-V$ characteristics of the isolated junction exhibits distinct cusps at voltages $eV=\Delta/n$, $n\in\{1,2,3\dots\}$ (MAR) Averin and Bardas (1995). The precise relation $I(\varphi)=I_{c}\sin(\varphi)$ is recovered in the tunneling regime $D\ll 1$ and $eV<\Delta$ with $I_{c}=2e\Delta D/h$.

Results for the RC-BdG model. The first electromagnetic environment we consider is a simple RC circuit as sketched in Fig. 1. The capacitance $C$ typically accounts for the electron-electron interaction in the junction itself while the resistance $R$ accounts for the finite residual resistance in the whole circuit. This $RC$ circuit is the minimum electromagnetic environment that must be considered. The RCJ model for this circuit (where the BdG equation is replaced by the current-phase relation) reads,

[TABLE]

where the time $t$ has been rescaled as $t\rightarrow\omega_{0}t$. $\omega_{0}=\sqrt{\hbar I_{c}/(2eC)}$ is the intrinsic frequency of the circuit for small oscillating amplitudes and $Q=RC\omega_{0}$ is the corresponding quality factor. The physics of this model is well understood Stewart (1968): for $I_{0}<I_{c}$ all the current passes through the junction (super-current branch) while for $I_{0}>I_{c}$ the equilibrium solution is unstable and a voltage develops across the junction. Interestingly, this model is hysteretic for underdamped circuits $Q>1$ and a dynamical solution with $\dot{\varphi}\neq 0$ exists for some values of $I_{0}<I_{c}$. At high bias current $I_{0}\gg I_{c}$, most of the current is dissipated by the resistor $R$ and the RCJ model predicts $I_{0}=R\bar{V}$ (where $\bar{V}$ is the average voltage difference seen by the junction. This prediction misses an important contribution from the junction, its intrinsic resistance $R_{J}=h/(2e^{2}D)$ in the normal state. Indeed, at large bias current, one expects $I_{0}=(1/R+1/R_{J})\bar{V}$

We now turn to the full simulation of the RC-BdG model. The bare simulation data is shown in Fig. 1 where the dashed line corresponds to a slow (quasi-static) ramp of $I_{0}$ so that the entire $I-V$ characteristics of the device can be extracted from a single simulation. We ramp the current first up and then down to zero in order to capture the hysteresis loop of the junction. The blue line corresponds to the voltage across the junction as a function of time. As shown in the inset, the blue line contains an important oscillating part that corresponds to the AC Josephson effect. From this data, we calculate the voltage $\bar{V}$ across the junction, averaged over a small time window (to get rid of the AC Josephson signal). Fig. 2 show the resulting plot of $I_{0}$ versus $\bar{V}$ (blue line). The dotted line corresponds to the various asymptotic of the RCJ model discussed above while the dashed line corresponds to the pure BdG model in the absence of the electromagnetic environment. The pure BdG model displays the usual kinks characteristics of the opening of a new MAR channel Cuevas et al. (2002). The RC-BdG simulations reconcile the two limits: the pure MAR curve at high bias and the supercurrent branch of the RCJ model at small bias. In the crossover between these two extreme limits, it provides the minimum model that captures all important physical contributions, hence quantitatively predicts the full hysteresis loop including the retrapping current. The most interesting features of the system show up in its dynamics. Recording the phase difference $\varphi(t)$ across the junction and the current $I(t)$ that flows through it, the dynamics is properly captured by the corresponding out-of-equilibrium current phase $I-\varphi$ relation obtained from the corresponding parametric plot. the result is shown in the upper panel of Fig. 2. Such out-of-equilibrium $I-\varphi$ could be reconstructed from a high-frequency measurement of the different harmonic of $V(t)$. As a reference, Fig. 2 includes the equilibrium $I-\varphi$ characteristics of the junction (dotted line) obtained by taking all contributions into account (i.e. both Andreev bound states and the small contribution from the continuous part of the spectrum). This equilibrium $I-\varphi$ relation contains small deviations to the sinusoidal form. However, out-of-equilibrium relations can be strongly different from the simple sinusoidal shape. This is true in particular in the returning part of the hysteresis loop (red line, square, visible component of the second harmonic) and close to the MAR cusps (yellow line, triangle, strongly non-sinusoidal). In these regimes, the excursions in voltage across the junction are wide (as can be seen directly from Fig. 1) and the junction effectively highly non-linear.

Results for the RLC-BdG model. We now turn to a second circuit where the junction is put in series with a classical $RLC$ resonator as sketched in the inset of Fig. 3. The electromagnetic circuit is slightly more complex than the previous $RC$ model, but in return, the highly non-linear behavior shown in the previous example manifests itself already on DC observables. The resonator has a quality factor $Q=R\sqrt{C/L}$ and a resonance pulsation $\omega_{0}=1/\sqrt{LC}$. The corresponding impedance $Z(\omega)$ takes the form $R/Z(\omega)=1+iQ[\omega/\omega_{0}-\omega_{0}/\omega]$ and filters frequency around $\omega_{0}$. Such an environment has been studied in a series of recent experiments using tunnel junctions Westig et al. (2017); Simmonds et al. (2004); Reagor et al. (2016); Zaretskey et al. (2013).

The RLC circuit provides a direct probe of the AC signal present in the system. We expect a main resonance for $2eV_{J}/\hbar=\omega_{0}$ when the AC Josephson effect drives the RLC circuit. Due to the non-linear character of the junction, higher harmonic of the AC Josephson effects are generated, so that additional features are expected at $2eV_{J}/\hbar=\omega_{0}/n$. Likewise, the non linearities imply that the RLC circuit can also be driven parametrically at $2\omega_{0}$ leading to features at $eV_{J}/\hbar=\omega_{0}/n$.

We compare the RLC-BdG calculations with an improved RLCJ model. The improved RLCJ model captures the super-current branch and the MAR non-linear I-V curve as,

[TABLE]

where $I_{\rm MAR}(V)$ is the DC non-linear $I-V$ characteristic of the junction in the absence of electromagnetic environment (dashed line of Fig. 4). The numerical results for the average current $\bar{I}$ versus voltage are shown in Fig. 4 for four different RLC circuits with different frequencies $\omega_{0}$. We concentrate on the main features around $2eV_{J}/\hbar=\omega_{0}$ and disregard the smaller peaks associated with higher harmonics and/or parametric pumping. The improved RLCJ model (dotted line) presents a Lorentzian like resonance at $2eV_{J}/\hbar=\omega_{0}$ for all four RLC circuits. When the resonance lies in the tunneling regime of the junction (blue line), there is a very good agreement between the improved RLCJ model and the full RLC-BdG simulations. The agreement is also qualitatively (but not quantitatively) good when the resonance corresponds to high voltages in the almost “Ohmic” regime of the junction (yellow line). However, for the two circuits where the resonance lies in the vicinity of a kink of the $I_{\rm MAR}(V)$ characteristic, the two models are strikingly different and the improved RLCJ model no longer applicable (green and red lines): the improved RLCJ model is typically off by $\pm 50\%$ including in the linewidth. In these regimes, we find that for $2eV_{J}/\hbar>\omega_{0}$, the current is reduced with respect to $I_{\rm MAR}(V)$ instead of the Lorentzian increase observed in the improved RLCJ model. This reduction of the current is a direct manifestation of the non-linear AC physics happening in the device. This DC prediction is the counterpart of the highly non-sinusoidal non-equilibrium current-phase relations discussed above for the RC-BdG case. However, the fact that the observable is in DC makes this prediction more easily accessible to an experimental test.

Conclusion. The Environment-BdG model presented in this manuscript unifies simple RCJ like models with microscopic models that include the quasi-particle spectrum of the junctions as well as its dynamics out-of-equilibrium. We have shown that the interplay between the two physics strongly modifies the behavior of the system and lead to new phenomena such as the voltage induced non-sinusoidal current phase relations. Our approach provides a practical route to study the engineering of electromagnetic environments in the presence of junctions that go beyond simple tunneling devices. Besides the example studied in this letter (a single channel junction with arbitrary transparency), other systems such as Josephson FETs Akazaki et al. (1996), Majorana devices Sarma et al. (2015), multi-terminal devices McCaughan and Berggren (2014) that are being developed by the community could be studied with the same technique Weston et al. (2015). On the technical level, the approach could be extended to include electron-electron interaction and/or a self-consistent calculation of the superconducting gap at the time-dependent mean-field level Kloss et al. (2018).

Acknowledgment. This work was supported by the ANR Fully Quantum (ANR-16-CE30-0015-02), ANR QTERA (ANR-15-CE24-0007-02), the NSF-ANR partnership PIRE: HYBRID and U.S. Office of Naval Research. We thank Manuel Houzet and Julia Meyer for useful discussions.