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

###### Abstract

Josephson junctions provide highly non-linear impedances at the root of many applications such as quantum limited parametric amplifiers Macklin2015; Bergeal2010 or superconducting qubits Vion2002; Makhlin2001. These junctions are often described by a sinusoidal relation which relates the current to the integral over time of the voltage across the junction . This relation properly captures the contribution of the superconducting condensate but not the quasi-particles that appear when the system is driven out-of-equilibrium Averin1995. Here, we construct a unifying framework that includes a microscopic description of the junction (full fledged treatment of the time-dependent Bogoliubov-De-Gennes equation) in presence of a classical electronic circuit. Our approach generalizes the standard Resistor-Capacitor-Josephson model (RCJ) to arbitrary junctions (including e.g. multi-terminal geometries and/or junctions that embed topological or magnetic elements) and classical circuits. We apply our technique to two situations. First, a RC circuit connected to single channel Josephson junction that exhibits Multiple Andreev Reflection (MAR) phenomena. We show that the theory properly describes both MAR and the hysteresis loops due to the electromagnetic environment. We show that out-of-equilibrium, the current-phase relation of the junction becomes strongly distorted from the simple sinusoidal form. Second, we embed the junction into a RLC circuit and show that the out-of-equilibrium non-sinusoidal current phase relation leads to a strong change of the shape of the resonance.

There is currently a huge effort around the world - both within academia and major industrial partners - to promote the superconducting transmon quantum bit Koch2007; Barends2013; Paik2011 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 (oxyde) barrier are the most mature elements, other types of junctions such as atomic contacts Grimes1966 (with very few propagating channels), semi-conducting nanowires Mourik2012 (with high spin-orbit suitable for stabilizing Majoranas bound states) superconducting-ferromagnetic-superconducting Izyumov2002; Buzdin2005 (with anomalous current phase relations) or multiterminal devices Riwar2016 could provide new functionalities to the superconducting toolbox. While the theoretical description of these objects is rather well understood Golubov2004, 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 models are commonly used to describe Josephson junctions circuits. The first one is the RCJ model Likharev1979; Stewart1968 (Resistor-Capacitor-Josephson) and its extensions to other classical circuits. In this model, one considers a classical circuit such as the ones shown in the insets of Fig. 1 or Fig. 3 with resistances , capacitances , inductances or other classical elements. The Josephson junction is described by its current-phase relation and the Josephson relation . Such a simple model is surprisingly powerful. It captures the hysteresis loops of the curves. Its simple extension, where one add 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 BenJacob1982. It has also been successfully used for more elaborate circuits that include resonators Cassidy2017. Its quantum extension provides the model used to design the various sorts of superconducting qubits Xiang2013 and has been shown to describe very accurately a large corpus of experimental data Devoret1990. 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 later 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 Klapwijk1982 (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 Averin1995, the interplay with microwaves Cuevas2002, ac Josephson effects and emergent topological effects in multi-terminal geometries Houzet2013. Until recently however, the direct numerical integration of BdG equations have been very limited due to the intrinsic computational complexity Weston2016 and did not include the electromagnetic environment of the junction.

The present letter builds on recent advances made in time-dependent computational transport Kloss2018 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 equation and classical electromagnetic environment.

Problem formulation. We model our circuits in two parts. First, the junction itself is described with a microscopic BdG Hamiltonian that depends explicitly on time (through e.g. a capacitive gate) and on the phase difference 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 from which one can compute the current that flows through the system. Below, we restrict ourself to the average current, but its quantum fluctuations are also accessible through our formalism Gaury2016. The second part of the model describe the classical circuit, or electromagnetic environment, that surrounds the junction. The classical equations that describe these circuits take the form of a differential equation for . More complex circuits are described in a similar way with more degrees of freedom describing the classical circuit. The set of equations reads,

(1a) | |||

(1b) | |||

(1c) |

where the function describes the dynamics of the classical circuit (RC or RLC equation in the examples below) with as an external source term. We numerically solve the BdG equation within the Keldysh formalism using the approach developed in Gaury2014; Weston2016b where the problem is unfolded onto a set of Schrödinger equations for scattering wave functions. An efficient algorithm has been constructed in Weston2016 to integrate the corresponding equations. The corresponding software, “t-Kwant” relies on the Kwant package Groth2014 and will be released as open source in a near future. Eq.(1a) alone amounts to solving a few hundreds 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 Kloss2018, 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 with the Fermi energy) and the evolution of the electromagnetic variables and 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 as explained in Kloss2018: Eq.(1a) is integrated with a “predicted” function and Eq.(1b) and (1c) are used on the larger time scale to construct this prediction. A straightforward time adaptive fourth order Range-Kutta Hairer1993 is used for the integration of Eq.(1c). We used the algorithm of Istas2018 to calculate the (Andreev) bound states of the model with precision.

To be specific, we now turn to a particular BdG Hamiltonian that describes a single channel junction. We emphasize that more complex Hamiltonians hosting, e.g. Majorana bound states Weston2015 or more channels can be studied with the same technique. The BdG Hamiltonian describes a one dimensional systems with the two superconductor corresponding to and while the normal region is formed by a single site at placing the system in the short junction limit,

(2) |

Here where is the voltage difference across the junction, is the superconducting gap inside the superconductors, is a potential barrier used to tune the transmission probability of the junction. In the calculations below we used , (which will be used as our unit of energy), and which corresponds to a junction with intermediate transmission . For this value, the equilibrium current-phase relation has small deviations with respect to a sinusoidal shape but the characteristics of the isolated junction exhibits distinct cusps at voltages , (MAR) Averin1995. The precise relation is recovered in the tunneling regime and with .

Results for the RC-BdG model. The first electromagnetic environment we consider is a simple RC circuit as sketched in Fig. 1. The capacitance typically accounts for the electron-electron interaction in the junction itself while the resistance accounts for the finite residual resistance in the whole circuit. This 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,

(3) |

where the time has been rescaled as . is the intrinsic frequency of the circuit for small oscillating amplitudes and is the corresponding quality factor. The physics of this model is well understood Stewart1968: for all the current passes through the junction (super-current branch) while for the equilibrium solution is unstable and a voltage develops across the junction. Interestingly, this model is hysteretic for underdamped circuits and a dynamical solution with exists for some values of . At high bias current , most of the current is dissipated by the resistor and the RCJ model predicts (where is the average voltage difference seen by the junction. This prediction misses an important contribution from the junction, its intrinsic resistance in the normal state. Indeed, at large bias current, one expects

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 ramp of . 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. We have checked that the ramp is slow enough to be considered as quasi-static, so that the entire characteristics of the device can be extracted from the data of a single simulation. From this data, we calculate the voltage across the junction, averaged over a small time window (to get rid of the ac Josephson signal). Fig. 2 show the resulting plot of versus (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 electromagnetic environment. The pure BdG model displays the usual kinks characteristics of the opening of a new MAR channel Cuevas2002. The RC-BdG simulations reconciles the two limits: the pure MAR curve at high bias and the supercurrent branch of the RCJ model at small bias. It also correctly captures for the first time (to the best of our knowledge) the second part of the hysteresis loop (upon decreasing the current). The most interesting features of the system show up in its dynamics. Recording the phase difference across the junction and the current that flows through it, the dynamics is properly captured by the corresponding out-of-equilibrium current phase relation obtained from the corresponding parametric plot. the result is shown in the upper panel of Fig. 2. Such out-of-equilibrium could be reconstructed from a high frequency measurement of the different harmonic of . As a reference, Fig. 2 includes the equilibrium 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 relation contains small deviations to the sinusoidal form. However, the 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 resonator as sketched in the inset of Fig. 3. The electromagnetic circuit is slightly more complex than the previous 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 and a resonance pulsation . The corresponding impedance takes the form and filters frequency around . Such an environment has been studied in a series of recent experiments using tunnel junctions Westig2017; Simmonds2004; Reagor2016; Zaretskey2013.

The RLC circuit provides a direct probe of the AC signal present in the system. We expect a main resonance for 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 . Likewise, the non linearities imply that the RLC circuit can also be driven parametrically at leading to features at .

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,

(4) |

where is the DC non-linear characteristic of the junction in the absence of electromagnetic environment (dashed line of Fig. 4). The numerical results for the average current versus voltage are shown on Fig. 4 for four different RLC circuits with different frequencies . We concentrate on the main features around 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 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 characteristic, the two models are strikingly different and the improved RLCJ model no longer applicable (green and red lines). In these regimes, we find that for , the current is reduced with respect to 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 there is more in it than the simple juxtaposition of the two physics, as evidenced by the failure of the improved RLCJ model. Our approach provides a practical route to study the engineering of electromagnetic environments in presence of junctions that go beyond simple tunneling devices. Beside the example studied in this letter (a single channel junction with arbitrary transparency), other systems such as Josephson FETs Akazaki1996, Majorana devices Sarma2015, multi-terminal devices McCaughan2014 that are being developed by the community could be studied with the same technique.

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