The Josephson junction as a quantum engine

We treat the Cooper pairs in the superconducting electrodes of a Josephson junction (JJ) as an open system, coupled via Andreev scattering to external baths of electrons. The disequilibrium between the baths generates the direct-current bias applied to the JJ. In the weak-coupling limit we obtain a Markovian master equation that provides a simple dynamical description consistent with the main features of the JJ, including the form of the current–voltage characteristic, its hysteresis, and the appearance under periodic voltage driving of discrete Shapiro steps. For small dissipation, our model also exhibits a self-oscillation of the JJ’s electrical dipole with frequency Ω=2eV/ℏ around mean voltage V. This self-oscillation, associated with ‘hidden attractors’ of the nonlinear equations of motion, explains the observed production of monochromatic radiation with frequency Ω and its harmonics. We argue that this picture of the JJ as a quantum engine resolves open questions about the Josephson effect as an irreversible process and could open new perspectives in quantum thermodynamics and in the theory of dynamical systems.


I. INTRODUCTION
In 1962, Josephson predicted that two superconducting electrodes separated by a thin insulating barrier should exhibit interesting properties associated with the tunneling through that barrier of Cooper pairs (bound states of two electrons, whose condensation at low temperatures accounts for superconductivity) [1].The tunneling current through such a "Josephson junction" (JJ) is usually expressed as where I c is a constant (the "critical current") corresponding to the maximum value of the supercurrent that the JJ can support and ϕ is the difference in the phases of the Ginzburg-Landau wavefunctions for each of the electrodes.If a direct-current (DC) bias is applied, the phase varies in time as where V is voltage jump across the junction.Combining Eqs. ( 1) and ( 2), the constant bias V should generate an alternating current (AC) through the junction of the form Josephson's main predictions were soon confirmed by experiment [2,3] and the JJ has since been an active subject of research in both pure and applied physics; see [4][5][6][7] and references therein.The JJ's conversion of a constant bias V > 0 into an alternating current I J (we refer to such a process as "DC → AC") implies that the JJ should be understood as an electrical self-oscillator [8].The first practical DC→AC converter was built by Hertz in 1887 using a spark gap (two conductors separated by a thin gap filled with ambient air), allowing him to produce electrical oscillations with frequencies ∼ 100 MHz [9].Hertz used this spark-gap oscillator to study the properties of electromagnetic waves.The spark gap then became the basis of radio technology, until it was replaced by the related "singing arc" of Duddell [10], and later by triode self-oscillators (see, e.g., [11,12] and references therein).These developments motivated a sophisticated mathematical theory of classical self-oscillation by Poincaré, Andronov, Hopf, van der Pol, and others.That approach has been based principally on the theory of ordinary differential equations and of nonlinear phenomena. 1eanwhile the subject of DC→AC and of self-oscillation more generally has attracted relatively little attention in theoretical physics.According to Le Corbeiller, there is some danger that the beauty of these kinematic developments could lead the student or reader into forgetting the physical relationships in the actual oscillator, the energy exchanges in particular, as Charles Fabry pointed out to the writer about 1929.Papers [15,16] were inspired by this remark.[17] The present work pursues, in a quantum context, the program suggested long ago by Fabry and Le Corbeiller and advanced more recently by some of the present authors of describing self-oscillators as engines that generate work cyclically and irreversibly, at the expense of an external thermodynamic disequilibrium [8,18].
Motivated by these considerations, in this article we propose a simple model of the JJ as an open quantum system, consisting of the superconducting Cooper pairs, weakly coupled to two external baths of electrons (each bath corresponds to a terminal of the battery that provides the DC bias).Our description is based on the Markovian master equation (MME) for the irreversible dynamics of the open system.Our model modifies the form of the relations between the JJ's tunneling current, voltage, and phase that appear in the literature [Eqs.(1), (2), and (3)], while giving the right properties for the JJ's current-voltage characteristic, including the hysteresis seen when the dissipation is small and the "Shapiro steps" that appear when the JJ is placed in an oscillating electric field.
The nonlinear dynamical equations resulting from our model also exhibit self-oscillations about stable fixed points on the JJ's current-voltage characteristic.For small dissipation and significant bias voltage V , we find that a finite but small perturbation can take the system out of the fixed point's basin of attraction, allowing the JJ's electrical Feynman presented a simple model for JJ in [33], based on a two-component MWF for the Cooper-pair condensate: where n 1,2 are number of Cooper pairs and ϕ 1,2 are phases corresponding to the left-hand and right-hand electrode of the JJ, respectively (see Fig. 1).The dynamics is given by the Schrödinger equation with Hamiltonian where ℏU is the energy difference between a Cooper pair located at the right electrode and one located at the right one.The value of U is related to the voltage bias V across the JJ by the formula while K the rate of tunneling of Cooper pairs across the barrier between the two electrodes.
In this treatment, the tunneling current is obtained from the unitary evolution according to Eqs. ( 5), (6), and (7), but the effect of that current on the charge of the electrodes is ignored.As Feynman explained, this is because the electrodes are connected to the terminals of a battery, which supplies the necessary current to keep the density of Cooper pairs in the superconducting electrodes approximately constant, so that the net charge density of the corresponding material (including the fixed positive charge density of the underlying atomic lattice) remains close to zero [33].This coupling to the terminals of the battery makes the JJ an open system.Describing the JJ using closed-system equations of motion (plus the ad hoc condition of constant Cooper-pair densities in the electrodes) obscures the irreversible dynamics by which the JJ can output work in the form of a self-oscillating [8] electrical dipole that generates non-thermal radiation.

B. RCSJ model
Let I be the external current applied to the JJ, as shown in Fig. 1.In the RCSJ model we have where the first term in the right-hand side of Eq. ( 8) corresponds to the Josephson current of Eq. ( 1).This current is taken to run in parallel ("shunt") with a current C V through capacitor C, and also to a current V /R through resistor R; see, e.g., [34,35].By using Eq. ( 2) to eliminate V in Eq. ( 8), we get the equation of motion In the R → ∞ limit, Eq. ( 9) can be obtained from the Hamiltonian where the charge q separated across the terminals of the junction is taken as a canonically conjugate variable to the phase ϕ.However, Eq. ( 10) is problematic in that the angle ϕ must be taken as a real number (i.e., a non-compact variable) with a "tilted washboard potential" that is unbounded from below.
In this article we will pursue a qualitatively different approach.By treating the JJ as an open system, we obtain a master equation that naturally incorporates a thermodynamic irreversibility in the form of the coupling of the Cooper pairs to two electronic baths that are out of equilibrium with each other.We will see that such a treatment incorporates a resistance R for the JJ that is not related to the phase-slip phenomena often invoked in the literature (see, e.g., [36]).In our model, the DC → AC dynamics of the JJ can be naturally understood as the irreversible operation of an engine, powered by the thermodynamic disequilibrium between the electrons in the two terminals of the external battery.

C. Feedback in engines
We use the word engine to denote an open system that can generate positive work over the course of a cycle, at the expense of an external disequilibrium.In the case of a heat engine, this disequilibrium is a difference of temperatures between two external baths, which causes heat to pass through the system.The disequilibrium may also be a difference of chemical potentials between the baths, causing matter (electrons in the case of the JJ) to pass through the open system.In addition to the external disequilibrium, an engine must exhibit feedback between the state of a macroscopic "tool" (such as a piston or turbine) acted upon by the engine's working substance, and the coupling of that working substance to the external baths.It is this feedback that causes an active force that does positive work on the tool over a complete cycle of the engine [18].
The laws of thermodynamics can helps us to understand better the general nature of this feedback.Consider an engine with a homogenous working substance that has internal energy U , quantity of matter N , and chemical potential µ.By the first law of thermodynamics, an infinitesimal transformation of this working substance gives a change of the internal energy where δQ is the infinitesimal heat absorbed by the substance from its surroundings, and δW is the infinitesimal work that the substance does on the tool (the symbol d is used for exact differentials, while δ is used for inexact differentials).Some authors refer to the µ dN term as "chemical work", but here we want to stress the distinction between this term and the actual mechanical work δW that is associated with a force acting on a directional and macroscopic degree of freedom, such as the position of a piston.Although the chemical contribution µ dN carries no entropy, its cyclical conversion into mechanical work is non-trivial: it requires feedback and dissipation, leading to an active force capable of driving the tool; see [18] and references therein.Over a full cycle the substance returns to its initial state, so that the change to U vanishes: The net work is, therefore, By Clausius's theorem, the entropy produced by the engine over one cycle is where T is the temperature of the source from which the working substance receives heat δQ at each instant during the cycle, T is the mean value of the temperature over the cycle's full period, and T d ≡ T − T .Let µ be the instantaneous chemical potential and μ its mean, with µ d ≡ µ − μ.Combining Eq. ( 14) and Eq. ( 15) we obtain This bound on W is achieved by a reversible cycle (Σ = 0).The result of Eq. ( 16) may be generalized to inhomogenous temperatures and chemical potentials by integrating over the maximum work that each part of the working substance may perform.The result of Eq. ( 16) tells us that for a heat engine to generate positive work (W > 0), its working substance must be at a higher temperature (T d > 0) during the part of the cycle during which it is more strongly coupled to the hotter bath (from which it takes heat δQ > 0).Conversely, the working substance must be at a lower temperature (T d < 0) when it is more coupled to the colder bath (into which it rejects heat δQ < 0).Thus, in any heat engine there must be positive feedback between the working substance's instantaneous T d and its coupling to the external baths (which gives the instantaneous δQ) [18].In [18] this general result is called the "Rayleigh-Eddington criterion". 2t is also evident from Eq. ( 16) that a heat engine works most efficiently if heat is absorbed by the working substance only at the maximum T d , while heat is rejected only at the minimum T d .In that case, Eq. ( 16) gives the Carnot bound.If T d = 0 (as in the case of the JJ, where the Cooper pairs are always at the ambient temperature), Eq. ( 16) reduces to Note that Eq. ( 17) applies even if the working substance exchanges heat with its surroundings, as long as the temperature is fixed.Thus, we see that in a chemical engine (i.e., in an engine which runs on a disequilibrium of the external chemical potentials, rather than on a thermal disequilibrium) there must be a positive feedback such that, during the cycle of operation, matter enters the engine's working substance at a higher chemical potential (µ d > 0 when dN > 0) and exits at a lower chemical potential (µ d < 0 when dN < 0).The efficiency of a chemical engine will always be less than unity, even without dissipation, because matter must exit at some non-zero potential µ.But there is no Carnot bound in this case nor (as far as we know!) any other universal bound independent of engine design.In the models for the "putt-putt" steam engine presented in [18,37], as well as in the "electron shuttle" and the "Quincke rotor" described in [18] the engine's dynamics is formulated in terms of two coupled ordinary differential equations: a mechanical equation of motion (second order in time) for a macroscopic degree of freedom (the "tool"), and a kinetic equation (first order in time) for the amount of matter in the working substance.The kinetic equation is explicitly irreversible (i.e., non-invariant under t → −t).A feedback appears between those two equations: the working substance exerts a force on the tool, while the tool modulates the coupling of the working substance to the external baths (and therefore the rate of change of the corresponding amount of matter).This feedback is what makes the force on the tool active, i.e., a non-conservative force that does positive work (W > 0) over a cycle of the tool's motion.Milburn underlines that this kind of dynamical system (which he applies to the electronic Schmitt trigger in sec.2.3 of [43]) is not Hamiltonian.
The relation between voltage V and phase ϕ in the RCSJ model is such that the voltage can be eliminated to give an equation of motion for ϕ alone [Eq.(9)].The absence of a stable equilibrium value for the angle ϕ (which is physically defined modulo 2π) requires reinterpreting ϕ as a non-compact variable with a potential U (ϕ) that is unbounded from below, as in Eq. (11).Our model of the JJ as an engine is qualitatively different.Instead of just the phase ϕ, we will work with the macroscopic coherence z (a complex number whose argument corresponds to ϕ).The equation of motion for z depends on the separation of charge between the two electrodes (or, equivalently, on the instantaneous voltage jump across the junction).A feedback appears in that the rate of change of this separated charge depends on z because of the tunneling between electrodes.The dependence between z and V is such that the dynamics cannot be reduced to an equation of motion for z alone.
When an external current source is applied to the JJ, the feedback can be positive, maintaining |z| against the damping associated with the decay of Cooper pairs in the electrodes into electrons in the external baths.This positive feedback sustains the rotation of the phase ϕ.Thus, in our model the JJ is considered as an engine in which the Cooper pairs in the electrodes constitute the working substance, while the JJ's macroscopic coherence z acts as the tool.The Rayleigh-Eddington criterion of Eq. ( 17) is then met, because more Cooper pairs are formed in the electrode at the higher voltage, while more Cooper pairs are destroyed in the electrode at the lower voltage.The resulting engine dynamics can also produce self-oscillation of the JJ's electric dipole, thus explaining how the device produces monochromatic radiation.
We will consider the precise dynamics of these various processes in the rest of this paper.However, it will be useful to keep in mind how the detailed model fits into this simple thermodynamic picture.At the most basic level, the point is that electrons come from the bath at higher potential, pass through the JJ (in the form of Cooper pairs), and end up in the bath at the lower potential.This flow sustains the JJ's macroscopic coherence z and turns its phase ϕ.This turning phase is a form of work, which can be used to produce non-thermal radiation.It may also be used to pump a current in another circuit, as has been shown experimentally in [30].This work extraction is possible because of a positive feedback between z and the differential rates at which Cooper pairs are produced or destroyed in the electrodes.
In previous theoretical treatments of microscopic engines (in addition to the references already cited see also, e.g., [44][45][46]) the tool degree of freedom is effectively classical.Our model of the JJ is new because the tool is the coherence z, which is a macroscopic but essentially quantum object.We therefore expect, even beyond the specific problem of understanding the Josephson effect as an irreversibly process, that the model that we present here may open new perspectives in the theory of quantum thermodynamics.

III. IRREVERSIBLE DYNAMICS
Our model is an extension of the model proposed by Feynman, in which the MWF will be replaced by the macroscopic density matrix (MDM), which is the appropriate mathematical object to describe the dynamics of the JJ as an open system. 3The coupling to baths of electrons gives a non-unitary evolution of this MDM that can be expressed in terms of a Markovian master equation (MME).From this MME we arrive at a formulation of the dynamics of the JJ in terms of coupled first-order ordinary differential equations for the time-dependence of the JJ's voltage V and its coherence z.

A. Markovian master equation
Let us begin by considering a single superconducting electrode, coupled to an external bath of electrons.The condensate of Cooper pairs occupies the ground state of a certain effective Hamiltonian, which can be accounted for using a single quantum harmonic oscillator with the operator a † a for the number of Cooper pairs.This number can be increased by the fusion of two electrons from the bath into Cooper pair, or decreased by the reverse process of Cooper pair decomposition.The effective Hamiltonian describing such processes takes the form where b † k , b k ′ are creation and annihilation operators for electrons of the bath.The interaction Hamiltonian of Eq. ( 18) corresponds, when restricted to on-shell processes, to the well known Andreev scattering, given by the coupling between electrons and holes in the normal metal and Cooper pairs in the superconductor [32].On the physics of Andreev scattering, see also [47] and references thererin.The results of this paper will not depend on the detailed form of the g kk ′ form factors in Eq. (18).
Under the standard assumption of separation of time scales of the system (Cooper pairs) and baths (electrons), one obtains the following Markovian master equation for the reduced density matrix of Cooper pairs (harmonic oscillator): where E is the renormalized value of ground state energy and γ ↓ > γ ↑ are decomposition and creation rates, respectively.The explicit formulas for those parameters are not relevant at the moment and can be found in the literature (see, e.g., [48]).Using Eq. ( 19) one can find the exact kinetic equation for the number of Cooper pairs The Markovian approximation is justified in this case because of the separation of time scales between the slow dynamics of the system composed of Cooper pairs in the electrodes and the fast internal dynamics of the electrons in the leads.For the Cooper pairs, the relevant time scales are the Josephson frequency of Eq. ( 2), the tunneling rate K, and the damping rate γ (we will define these carefully in Sec.III B).In a common experimental setup, all of these are ∼ 10 −12 s; see, e.g., [49].Meanwhile, the time scale for the internal dynamics of the leads can be estimated from the Drude model of electrical conduction to be ∼ 10 −14 to 10 −15 s.There is, therefore, at least two orders of magnitude of separation between the respective time scales of the system and the baths for the physical setups that we are interested in describing.

B. Josephson Junction model
In order to describe the JJ, we must extend the model above so as to describe two superconducting electrodes separated by a potential barrier, each of them coupled to its own electron reservoir.Consider, therefore, two Cooperpair modes 1 and 2 corresponding to the two electrodes in Fig. 1.The system is then equivalent to two quantum harmonic oscillators, with reduced dynamics given by an extension of Eq. ( 19), to wit: The Hamiltonian H 12 takes the form where E 1,2 are the energies for each of the two electrodes, while K is the tunneling rate between the electrodes, which can be taken positive by absorbing the phase.Note, that we consider so-called local master equation [50].This approach is justified when the perturbation is small compared to the energy scale of the free Hamiltonian [51,52].In our case, we treat the tunneling rate K as a small perturbation compared to the electrode energies E 1,2 in Eq. ( 22), which is consistent with usual approach in theory of superconductors.
The dissipative generators L 1,2 in Eq. ( 21) correspond to the coupling of each of the two electrodes of the JJ to the electronic bath with which that electrode is in physical contact.These generators have the same structure as in Eq. ( 19), with In order to obtain the generalization of Eq. ( 20) corresponding to the two electrodes of the JJ, we need to define a single-particle density matrix.By analogy with the macroscopic wave function of Eq. ( 4), we will call this new object the macroscopic density matrix (MDM).In this case, the MDM is a 2 × 2 positively defined matrix σ, given by the following two-point correlations: A useful parametrization of MDM is given by the populations n j = Tr(ρa † j a j ) and coherence z = Tr(ρa † 1 a 2 ): Notice that, for the pure state corresponding to Eq. ( 4), Eqs. ( 21) -( 25) translate into the following evolution equations for the MDM: where ↑ , and nj ≡ γ ↑ /γ (j) for j = 1, 2. Let us assume that the relaxation rates for both electrodes are equal, i.e., Then, introducing new variables for the sum and difference of the expected number of Cooper pairs on each of the JJ's electrodes, we rewrite Eq. ( 27) as This implies that the average total number of Cooper pairs N relaxes, independently of all the other dynamical variables, to its stationary value (N = N ), which must correspond physically to an electrically neutral JJ.This justifies our having taken equal relaxation rates for both electrodes in Eq. ( 28).
The variable n describes an internal local disequilibrium that causes an electrical double layer to form at the junction, with excess charge Q = en and associated voltage jump V .Accordingly, we introduce the capacitance of the junction C, such that the voltage is given by the formula: We can interpret the energy difference U between the two electrodes as the work needed to transfer a single Cooper pair between the electrodes, against the external voltage.Thus: Note that, by making the voltage across the junction depend on the difference n of the numbers of Cooper pairs, we have introduced a mean-field approximation.This is justified as long as n 1 , n 2 ≫ |n|.The tunneling rate K will also depend on n, but that dependence -unlike the dependence of V on n-is not essential to the JJ dynamics that we will describe here.For reasons of simplicity we will take K as constant.Combining all of the above definitions, we get a dynamical system characterized by the coupled differential equations where we have put: Finally, we introduce a resistance parameter R such that in order to facilitate comparison between our results and those of the RCSJ model.Using this, we may rewrite Eq. ( 35) as where the last term can also be defined, in terms of the Hamiltonian dynamics, as: with q = −2e the charge of a Cooper pair.Comparing Eq. ( 39) to Eq. ( 8) of the RCSJ model, we interpret the quantity I defined by Eq. ( 37) as the fixed current that the external source provides to the JJ. 4 We also see that the term V /R in Eq. ( 39) results from the finite lifetimes of the Cooper pairs (which decay into electrons in the baths), while the tunneling of Cooper pairs contributes a current 4eK|z| sin ϕ.For simplicity, we neglect the contribution to the JJ current from quasiparticles.Note that |z| is not fixed, because the n 1,2 of Eqs. ( 24) and ( 25) are dynamical variables.This is a key difference between our model and that of RCSJ, in which the tunneling current was taken to be I c sin ϕ, for fixed I c .Note also that Eq. (36) implies that the standard relation between voltage and phase as given in Eq. ( 2) and assumed by the RCSJ model, is not exact.Therefore Eq. ( 3) for the AC tunneling current (a quantity which has never been observed directly in experiments), is not applicable to our model.
We see now how the feedback that we discussed in general thermodynamic terms in Sec.II C appears mathematically in our dynamical model: Eq. ( 35) for the real-valued V depends on z via the tunneling term proportional to K Im (z).Meanwhile, Eq. ( 36) for the complex-valued z depends on V , which controls the resonant frequency 2eV /ℏ.It is therefore plausible that, in some regime, this non-Hamiltonian dynamical system will self-oscillate: An oscillation of Im (z) forces V in in Eq. ( 35) into an oscillation with the same frequency.Meanwhile, an oscillation of V about a non-zero average value can excite z via parametric resonance in Eq. ( 36).Thus we expect that the feedback between V and z may, in certain regimes, lead to a self-oscillation.We will investigate this question mathematically in Sec.V.
Note that I > 0 in eq. ( 35) is analogous to "population inversion" in the theory of a laser, in that it corresponds to n1 > n2 [see Eqs. ( 31) and ( 37)], with the energy of the macroscopic quantum state |1⟩ greater than that of |2⟩ by the amount U of Eq. ( 34).The same is true for I < 0, because then n < 0 and (by Eq. ( 34)) we have that the energy of |2⟩ is greater than that of |1⟩.Thus, it is the external current I that makes the JJ an active system, from which work may be extracted.The JJ is therefore analogous to a laser in that the "population inversion" induced by the external I is what, when combined with a positive feedback involving the coherence z (an effect analogous to stimulated emission), can result in a self-oscillation. 5 It will be convenient to transform a set of Eqs. ( 35) and ( 36) into dimensionless units.First, we separate the second equation into real and imaginary part of z obtaining where we have introduced the variable I S = 4eK Re z = 4eK|z| cos ϕ, so that Let us also introduce the characteristic voltage and current parameters and define new dimensionless variables as: The equations of motion for the JJ can now be expressed as 4 A full treatment of the JJ's dynamics would require us to include in the Hamiltonian the electrostatic interaction between the Cooper pairs and also with the background of positive ions in the electrodes, as Feynman suggests in [33].That interaction makes it energetically costly to vary the local density of Cooper pairs from its equilibrium value.By making the relaxation rates equal [Eq.( 28)] (which allows the JJ to relax to state with net zero electric charge, independently of n and z) and by taking the I defined by Eq. ( 37) to be equal to the fixed current provided by the JJ's coupling to the external baths (an assumption justified by the fact that form of the dependence on I of Eq. ( 35) corresponds to an electrostatic capacitor connected to a current source), the main consequences of this interaction are incorporated into our dynamical equations without having to treat it explicitly. 5In this context, it may be worth recalling the lesson of Borenstein and Lamb that the essence of lasing is not intrinsically quantum mechanical [53].Sargent, Scully, and Lamb emphasized the similarity of lasing with classical self-oscillation (which they called "sustained oscillation") in [54].FIG. 2. Current-voltage characteristic [Eq.( 48)], expressed in terms of a dimensionless total applied current itot and voltage v across the junction [for definitions of these dimensionless variables, see Eq. ( 44)]: (a) For α = 6.The critical current ic [Eq.( 50)] and the retrapping current ir [Eq.( 51)] are marked by the horizontal dashed lines.The part of the characteristic with negative slope is drawn as dotted rather than solid (see Sec. IV B for the demonstration that points along this negative slope are unstable equilibria).(b) For α = 2, in which case the characteristic is monotonically increasing and no hysteresis is observed. where Our model is therefore characterized by a single parameter α: the ratio of the squares of the tunneling and dissipation rates.

IV. CURRENT-VOLTAGE CHARACTERISTIC
From the dynamical Eqs.(35) and (36) we may deduce the stationary current-voltage (I-V ) relation (or "characteristic") for the JJ.One can easily solve the fixed point equations ( V = ż = 0) by eliminating z to obtain The second term inside the brackets on the right-hand side of Eq. ( 47) is a Lorentzian function centered at V = 0.It comes from the tunneling current's contribution to the DC component of Eq. ( 39), a contribution that is important only near V = 0.For K ̸ = 0, the average power consumed by the JJ is greater than the ohmic P diss = IV .It is this additional power that makes the JJ an active element.In the dimensionless units introduced before (see Eqs. ( 43), (44), and ( 45)), the I-V characteristic of Eq. ( 47) takes the simple form: A. Critical and re-trapping currents For α > 2 (i.e., K > √ 2γ) the function i tot (v) given in Eq. ( 48) is not monotonically increasing.The two positive solutions to i ′ tot (v 0 ) = 0 are such that a negative differential resistance, i.e., i ′ tot (v) < 0, is observed for voltages The positive local maximum i tot (v − ) we identify with the critical Josephson current: The positive local minimum i(v + ) (which we will identify with the retrapping Josephson current) is: Thus, the critical value α = 2 -for which i c = i r = 3 √ 3-discriminates between two regimes: with and without the negative differential resistance of JJ.In Sec.IV B we will show that the points along the part of the characteristic curve that have negative differential resistance are unstable, explaining the hysteresis of the underdamped JJ.
Note that Eq. ( 50) defines a critical current only for α > 2 (the underdamped regime).For α ≤ 2 (the overdamped regime), a natural definition of the critical current is given by the inflection point: Hence, the expression for critical current for full range of parameter α is given by These results are illustrated in Fig. 2, for two different choices of the parameters α.In Fig. 2(a), for α > 2, the characteristic i tot (v) has local maxima and minima, and therefore well-defined values of i c and i r .In Fig. 2(b) the value α ≤ 2 the characteristic i tot (v) is monotonically increasing and therefore has no i r .We will show in Sec.IV B that the first characteristic exhibits hysteresis, whereas the second does not.
Note that the critical current i c is given by the expression: which is a definition of the Stewart-McCumber parameter β c in the RCSJ model.Thus, our model exhibits hysteresis for β c ≥ 3 √ 3, whereas in the RCSJ model, hysteresis is usually said to correspond to β c ≳ 1.In the limit α → ∞ (corresponding to a damping rate γ much smaller than the tunneling rate K), Eq. ( 53) implies that i c → 2α.In that case the physical value of the critical current I c is We will use this limiting value of the critical current in Sec.VI A, where we study the efficiency of the JJ as an engine.

B. Stability and hysteresis
To determine the stability of fixed points we introduce the small deviations δv, δi J , and δi S of the dynamical variables away from their respective stationary values v 0 , i J,0 , and i S,0 , such that Inserting this into Eq.( 45) yields the linearized equation The eigenvalues of this linear set are therefore the three roots {λ i } of the characteristic polynomial Then, according to the Routh-Hurwitz criterion, the equilibrium at δv = δi J = δi S = 0 will be unstable (i.e., at least one root will have Re λ i > 0) if and only if However, let us observe that such that the unstable regime given by condition (59) precisely corresponds with negative differential resistance, i ′ tot (v 0 ) < 0. The instability of points along the current-voltage characteristic that have negative slope implies that, for α > 2 (i.e., in the underdamped regime), there will be a hysteresis loop bounded from above by i c (Eq. ( 50)) and from below by i r (Eq.( 51)).If the applied current i tot is initially zero and is then increased slowly, the resulting voltage will follow the nearly vertical segment of the characteristic (see Fig. 2(a)) until the critical current i c is reached.Since the part of the characteristic with negative slope is unstable, the current will then stay very nearly constant while the voltage v jumps quickly to the value corresponding to the intersection of the approximately ohmic branch of the characteristic with the horizontal line i tot = i c .Further increase of the current will move the voltage up along the ohmic branch.If the voltage is then slowly dialed down, the voltage will decrease smoothly until the value i tot = i r is reached, whereupon the current will stay almost fixed while the voltage jumps quickly to a value close to zero (corresponding to the intersection of the nearly vertical branch of the characteristic with the horizontal line i tot = i r ).
We have verified this hysteretic behavior in numerical integrations of Eq. ( 45) with a time-dependent bias current i tot = i tot (τ ).In particular, we let i tot start at zero, increase linearly up to a value above i c , and then decrease linearly until it returns to zero.The response of the JJ is shown in Fig. 3 for α = 4.The hysteresis loop appears clearly in Fig. 3(a).
We may also study the responses of the Josephson (i.e., tunneling) current i J , the resistive current i res = v, and the capacitive current i cap = dv/dτ , such that i J + i res + i cap = i tot .In Fig. 3(b) we see a sudden jump in i res shortly after i tot surpasses the critical value i c .In Fig. 3(c) we see how i res decreases rapidly shortly after i tot falls below the retrapping i r .These jumps in i res are accompanied by changes in i J in the opposite direction, as well as by brief spikes in i cap .

C. Bifurcation analysis
Much like in the "putt-putt" engine model of [18,37] and in the treatment of the electronic Schmitt trigger in [43], ours is a non-Hamiltonian dynamical system in which an explicitly irreversible "kinetic equation" [Eq.(35)] is coupled to a "mechanical" equation of motion for an oscillator [Eq.(36)].Positive feedback between these equations can lead to a self-oscillation, which in the models of [18,37,43] appears mathematically as a Hopf bifurcation [55].It is therefore natural to ask whether something similar happens in our model of the JJ.
Of the three roots of the polynomial in Eq. ( 58), one is real (we label it as λ 0 ) and the other two are conjugate: It can be shown that for α > 0 the real part of the complex roots is non-positive, i.e., κ ≤ 0. This implies that a bifurcation appears only for the real λ 0 , which is therefore not a Hopf bifurcation (i.e., the unstable eigenvalue is not associated with an oscillation frequency η > 0).Thus, our model cannot exhibit limit cycles around unstable equilibria, as in the self-oscillating dynamical systems considered in [8].
We will return to the question of self-oscillation in Sec.V.There we will show that our model does exhibit selfoscillation, but that this appears as a hidden attractor around a linearly stable fixed point, rather than as a limit cycle associated with a Hopf bifurcation [19].It will therefore require a small but finite perturbation away from the fixed point to get the system to self-oscillate.

D. Interference
Our model can also be applied to interference phenomena.Consider the superconducting quantum interference device (SQUID) represented in Fig. 4. We modify the tunneling constant K by introducing an external magnetic field (via the vector potential A) and we introduce the two tunneling channels through the insulators a and b, such that where K A,B are the tunneling rates and A,B are the integrals along the paths passing through insulators A and B, respectively.The two terms in the right-hand side of Eq. ( 62) correspond to two tunnelling possible tunneling paths, associated with the two junctions.So far we have considered real tunneling constant.For complex K, the tunneling part of the Hamiltonian ( 22) should be rewritten as We then obtain the following equations: It follows that the I-V characteristics has the same form as in Eq. ( 48), with α = K 2 /γ 2 replaced with |K| 2 /γ 2 .The expressions for the maximal and retrapping currents are also the same as before, with K replaced by |K|, where It we write the magnetic flux through the loop as Φ = A • ds, this becomes where the last term describes a well-known interference effect.The critical current is almost proportional to |K| 2 and therefore also oscillates as a function of the flux Φ.

E. Shapiro steps
Josephson had predicted in [1] that if a JJ were driven by an external oscillating electric field with frequency Ω f , then the I-V characteristic would be modified, with steps appearing at values of the bias voltage V 0 such that for n ∈ Z.This effect was found experimentally by Shapiro [3], and is now used as a metrological standard for the determination of voltages [56].The usual explanation of this effect (see, e.g., [34]) is that a voltage bias causes, by Eqs. ( 1) and ( 2), a tunneling current The sine of the sine can be expanded using Bessel functions, so that The terms of the expansion in Eq. ( 71) all average out to zero over long times unless ω 0 = nΩ f exactly, in which case a DC contribution would appear.Although this argument agrees with the observed locations of the Shapiro steps, it fails to account for their widths [34].Various elaborations of the theory have therefore been proposed, based phenomena such as frequency mixing [57], resistive feedback [58], or non-linear resonance [59].In general, these models assume that the external current I applied to the JJ also oscillates with frequency Ω f even though in experiments (including the original one by Shapiro) what is applied to the JJ is an oscillating electric field.
Unlike in the RCSJ model, in our model it is straightforward to account for the effects of this external AC field.For this, we extend the Hamiltonian of (22) to include the dipole interaction with the field: where ε cos(Ω f t) is a time-dependent external field with the amplitude ε and frequency Ω f , and ed(a † 1 a 1 − a † 2 a 2 ) is an electrical dipole moment of JJ.This gives rise to the following inhomogeneous extension of Eqs.(45): where In Fig. 5 it is presented the I-V characteristics of JJ in a presence of the oscillating electric field, where the voltage is additionally averaged in time.The standard Shapiro's structure is observed with steps appearing at values predicted by Eq. ( 68), which in the introduced variables takes the simple form: v = kω f .

V. SELF-OSCILLATION
We have seen that the simple dynamical model of Eqs. ( 35) and ( 36) can explain the hysteresis seen in the currentvoltage relation for the JJ when the dissipation is small, as well as the frequency locking ("Shapiro steps") response to an external driving by an applied AC electric field.Remarkably, the same model exhibits self-oscillations of the voltage and the electrical dipole around stable points in the I-V characteristic (i.e., along parts of the characteristic with positive differential resistance).Even though those fixed points are stable, a finite but small perturbation can move the system out of the basin of attraction of the fixed point, leading to a small self-oscillation with frequency Ω = 2eV /ℏ, where V is the voltage bias.This corresponds to a "hidden attractor" [19].
In [19] the authors draw a distinction between "self-excited" oscillations, associated with an unstable fixed point, and "hidden" oscillations which are not.An infinitesimal perturbation is enough to take a self-excited system away from the fixed point and towards the limit cycle or strange attractor.The lack of such an unstable equilibrium makes the study and simulation of hidden attractors considerably more difficult.On the other hand, the hidden-attractor oscillations that we consider here for our model of the JJ fall within the scope of Andronov et al.'s conceptualization of a "self-oscillator" as an open system that generates and maintains a periodic variation at the expense of an external source of power that has no corresponding periodicity [8,22].The presence of hidden attractors in our model of the JJ can explain how the JJ acts as an engine, producing work in the form of non-thermal sound and electromagnetic waves at the frequency Ω and its lower harmonics.

A. Hidden attractors
The prediction and analytical characterization of hidden attractors is an open problem in the mathematical theory of dynamical systems.Here we will present a simple argument why the dynamical system described by Eq. ( 45) may exhibit hidden limit cycles.The equations of motion for the system can be expressed as where ζ is a complex variable such that Re ζ = i S and Im ζ = i J , while the dot indicates derivation with respect to the dimensionless time variable τ [see Eqs. ( 43) (44), and ( 45)].In terms of the Fourier expansions and taking i tot = const.,Eqs. ( 76) and ( 77) become We expect the hidden-attractor oscillation to be dominated by the fundamental frequency ω.Taking Eq. ( 79) for n = 0 and Eq. ( 80) for n = 1, and neglecting v n , ζ n for all |n| > 1, we obtain The equilibrium value of ζ is, according to Eq. ( 77), where v 0 is the corresponding equilibrium value for v (the same quantity as in Eq. ( 56)).For |v 0 | ≫ 1 we therefore have Putting Eq. (84) into Eq.( 81) we get i.e., v 0 is approximately equal to the external current bias applied to the JJ.This corresponds to the stable and approximately ohmic branch of the I-V characteristic for the JJ.
In the regime of large applied current (|i tot | ≫ 1) and small damping (4α ≫ 1), Eqs. ( 81) and (82) are consistent with a hidden-attractor self-oscillation (|v 1 | > 0 and |ζ 1 | > 0) with a fundamental frequency Note that in the original variables [see Eqs. ( 43) and ( 44)], this fundamental frequency is This agrees with what we expected from the nature of parametric resonance, as was explained in words in Sec.III B, after we formulated Eqs. ( 35) and (36).In Sec.V B we discuss how such hidden-attractor oscillations can be seen in numerical simulations.There we will measure frequency in terms of an equivalent v 0 .

B. Numerical simulations
As we already alluded to in Sec.IV C, the mathematical appearance of such a self-oscillation is different from the Hopf bifurcations considered in [8,55] and in the "putt-putt" engine model of [18,37].The form of Eqs. ( 76) and ( 77) is closer to the model of the Quincke rotor in [18], but the Quincke rotor's self-sustained rotation appeared as a simple pitchfork bifurcation, which is not the case in our model of the JJ.As we will see, our model of the JJ is more subtle from the point of view of the mathematical theory of dynamical systems, since it exhibits hidden attractors around stable fixed points.
In order to find hidden attractors in our model and characterize some of their properties, we integrate Eqs. ( 76) and ( 77) numerically, for various values of the parameter α and for various choices of initial conditions.The results FIG. 7. Value of the fundamental peak of the frequency spectrum for the self-oscillation of the tunneling current iJ [labelled vspectra and measured in the units defined in Eqs. ( 43), (44)] as a function of the applied voltage v0, for α = 2.2.The spectra are computed from the numerically integrated orbits for 4850 ≤ τ ≤ 5000.
shown in this section were obtained using the LSODA method, as implemented in the ODEPACK library [60].Almost identical results were obtained using the BDF method, as implemented in the CVODE library [61].
As expected from the argument in Sec.V A, we find numerical evidence of hidden attractors around stable equilibrium points for sufficiently large values of α and v.To find them, we displaced the initial conditions slightly away from the equilibrium δv = δi J = δi S = 0 [see Eq. ( 56)] and we then observed the corresponding orbit after long times.Figure 6 shows one such hidden attractor, for an orbit with initial values v = 30, δv = δi J = 0, and δi S = −0.1In general, we find that hidden attractors appear only for values of α large enough that the corresponding I-V characteristic presents hysteresis (see Sec. IV A).Since the equilibrium point ζ 0 inside the hidden attractor is stable, it has a finite basin of attraction.For the same parameters used in the simulation of Fig. 6 but with initial |δi S | ≲ 10 −6 (while keeping the initial δv = δi J = 0), we find that the resulting orbit decays to ζ 0 , rather than giving a persistent oscillation.
Thus, a small but finite initial disturbance is needed to excite the JJ's self-oscillation.This is the defining feature of a "hidden attractor".The power spectrum for i J shows that the hidden attractor corresponds to an oscillation whose fundamental frequency is v 0 (in the dimensionless variables defined in Eq. ( 44)).This relation between the fundamental frequency of the hidden-attractor oscillation and the applied voltage v 0 holds over a range of values of v 0 , as shown in Fig. 7.This result is consistent with the fact that the electromagnetic and phonon radiation from a DC-biased JJ has fundamental frequency Ω = 2eV /ℏ.Note, however, that in our model the simple formula of Eq. ( 3) does not hold.The conversion between the applied DC bias and the observed AC emission (what we called "DC → AC" in Sec.I) appears here as the result of an irreversible engine dynamics, associated with the presence of hidden attractors in the solutions to our equations of motion.
For α < 2 -in which case the I-V characteristic for the JJ shows no hysteresis-we have not found hidden attractors in our numerical simulations.Instead, we find that initial conditions close to the linearly stable equilibrium decay steadily towards that equilibrium.An example of such behavior is shown in Fig. 8.A more precise and extensive characterization of the dynamical system defined by Eq. ( 45) is left for future investigation, as it constitutes and interesting mathematical problem in its own right.Here our goals has been only to establish the presence of hidden attractors whose behavior is consistent with the DC → AC properties of the Josephson effect.
We consider it non-trivial that our simple model of the JJ gives an adequate picture not only of the unusual properties of the I-V characteristic (as explained in Sec.IV B) and of its response to external driving by an AC electric field (as explained in Sec.IV E), but also of its autonomous operation as an engine (self-oscillator) that acts as an irreversible DC → AC converter.We believe that this calls for further investigation not only in terms of the physical theory of open quantum systems, but also of the mathematical theory of dynamical systems.The aim should be to clarify conceptually the relation between these different phenomena, all of which are characteristic of an active electronic device.
The numerical study of hidden attractors is a rather subtle problem that deserves a deeper investigation in the mathematical theory of nonlinear dynamical systems.For our present purposes it suffices to show that self-oscillations persist over very long times: we observed this for τ as large as 5000.We were able to establish this using both the CVODE and the ODEPACK libraries to carry out our numerical integration.Since the ODEs for our model might be stiff (see, e.g., [62]), it is important to establish that the self-oscillations observed are not due to instabilities of the integration method.We therefore solved the equations using various different stable methods, including the CVODE and ODEPACK libraries already mentioned, and compared the results.All of the methods that we used gave long-lived self-oscillations with the correct Josephson frequency.However, these self-oscillations did not appear to be infinitely long lived (i.e., hidden attractors) in simulations with some other commonly used numerical methods.In particular, as shown in Fig. 9, for the same initial conditions as in Fig. 6, the SUNDIALS library [63] (which uses the Radau method) approaches the stable fixed point.
Comparing the LSODA and Radau trajectories shown in Fig. 9, we find that the former (which seems to show a hidden attractor) is much smoother than the latter (which does not).We are therefore confident that the former is more likely to reflect the correct long-term behavior of our dynamical system of interest.The trouble with the Radau solution of Fig. 9 seems to be related to accumulation of round-off error and to the time steps being too large.This method's implementation varies adaptively between 5th and 13th order.The sensibility at different time steps may be misestimated, leading to a wrong interpolation.
Another interesting problem raised by our model concerns the possibility of chaos (strange attractors), both in the homogeneous system of Eq. ( 45) and in the inhomogeneous system with an external periodic driving of Eq. ( 74).Although we must leave this for future investigation, here we make some general remarks about this problem.
In van der Pol's simple model and related formulations of self-oscillation (see [8] and references therein), the total phase-space for the homogeneous system has only two dimensions and (by the Poincaré-Bendixson theorem) strange attractors are not possible; see, e.g., [64].This is also true of the RCSJ model as described by Eq. ( 9).However, as we stressed in Sec.II C, these models are not physically realistic because they invoke a negative linear damping, rather than describing the active dynamics of the self-oscillators in terms a feedback between the oscillator and another degree of of freedom.If that degree of freedom is added, the Poincaré-Bendixson theorem becomes inapplicable.Thus, e.g., the presence of both limit cycles and strange attractors has been shown numerically for the homogenous "leaking elastic capacitor" model, which describes a classical engine based on feedback between a mechanical and an electrical degree of freedom [65].
On the other hand, one of the earliest results in what would later be called "chaos theory" was the work of Cartwright and Littlewood on the forced (inhomogenous) van der Pol equation (see [66,67] and references therein).The appearance of chaos is associated with the presence of a quasiperiodic regime in which the ratio of the frequency of the external driving and the natural frequency of the self-oscillator is irrational.An interesting open question is whether something similar occurs for the inhomogeneous Eq. ( 74) and whether this might have something to do with the structure of the Shapiro steps that we obtained in Sec.IV E.
So far we have not identified any chaotic regimes for Eqs.(45) or (74).We expect that the investigation of that problem will require a better handle on the numerical methods suited to the detailed study of hidden attractors in general and of our dynamical model of the Josephson effect in particular.Identifying attractors in numerical simulations of open systems is a non-trivial problem that may be of considerable physical interest (see, e.g.[68]).

VI. ELECTROMAGNETIC RADIATION
Self-oscillation of the JJ produces both electromagnetic and acoustic waves, with fundamental frequency Ω = eV /ℏ and weaker harmonics.In the electromagnetic case, the corresponding wavelength is much longer than the JJ's linear dimension, suggesting a collective character of the emission process.Moreover, since Cooper pairs are bosons, their FIG.9. Integrated orbit found using the LSODA (in blue) and Radau (dotted green).The green dot denotes the initial point, while the orange one denotes the fixed point ζ0.This curve has the same parameters as the Fig. 6, but with initial condition δiS = −1 × 10 −4 .The hidden attractor is only observed in the LSODA integration wave function is symmetric under permutations.This makes our model of the JJ well suited for implementing Dicke's superradiance mechanism [69].Such radiation is highly monochromatic and coherent, and when the JJ is placed inside a resonant cavity the emission is considerably enhanced, as shown experimentally in [25,27].

A. Superradiance dynamics
Here we will briefly discuss the phenomenon of superradiance for an idealized system of N two-level "atoms" interacting with an electromagnetic field at zero temperature.We assume that the size of the sample is small compared to the wavelength of emitted radiation, leading to collective emission of radiation with intensity proportional to N 2 .If the initial state of the atomic ensemble is a product of N identical 2 × 2 density matrices, one can show that in the limit of large N the product structure is preserved and that the time evolution of the system may be described by the non-linear master equation for the corresponding MDM, normalized to N : For details on this formulation, see [48].
Here, ω A is an atomic frequency and γ e is the spontaneous emission rate for the single atom.Denoting by |1⟩, |2⟩ the excited and ground state respectively, we define the "spin" matrices by Inserting, the same parametrization of the MDM in terms of n 1,2 and z as in Eq. ( 25) into Eq.(88), we obtain Because the total number of atoms N = n 1 +n 2 is preserved, we can rewrite this in terms of the variable n Note that, due to the mathematical structure of Eq. (88), the initially pure MDM remains pure (i.e., |z| 2 = n 1 (N −n 1 )) and the evolution of excited state population satisfies the well known equation This reproduces the N 2 proportionality and the bell-shaped time dependence of the intensity of radiated energy that are seen in superradiance.The radiated power is given by Note that it is not easy to implement this superradiance effect with actual atoms (see, e.g., [70]).This requires that the distance between the atoms be much smaller than the radiation's wavelength and also that the atoms occupy a symmetric configuration in space.Only under those two conditions would it be possible to express the interaction of the atoms with the electromagnetic field as the sum of dipoles multiplied by a single field operator.However, in the case of the JJ the two levels correspond to the left-hand and right-hand electrodes, while the transition between levels corresponds to tunneling of a Cooper pair across the junction.In this case, the two conditions for superradiance listed above are satisfied automatically.First, the linear size of the junction is of the order of micrometers, while the observed radiation is typically submillimeter.Second, Cooper pairs are (to a very good approximation) non-interacting bosons whose state is symmetric with respect to exchanging the two electrodes.Although the superradiance model was introduced in a very different context, it is therefore ideally suited to describe the production of non-thermal radiation by the JJ.
To calculate the efficiency with which the JJ transforms the applied external power into radiation, we may add to the JJ equations of motion for V and z [Eqs.(35) and (36)] the superradiant contributions resulting of Eq. ( 91), recalling that V = en/C [Eq.(33)].The modified JJ dynamics is then governed by The correction −2eγ e |z| 2 to the external current I in Eq. ( 94) corresponds, when multiplied by the voltage V , to the radiated power.We can therefore introduce a parameter characterizing the JJ's radiation efficiency under stationary conditions: where z 0 is the stationary value of z in Eqs. ( 94) and (95).Note that the radiated power 2eV γ e |z 0 | 2 computed from Eq. ( 94) agrees with the expression of Eq. (93) if we take an "atomic frequency" ω A = 2eV /ℏ.Rather than using the dimensionless variables that we introduced in Sec.III [see Eq. ( 44)], here it will be more convenient to work with the angular frequency and the Cooper-pair number current Introducing the frequency parameter Eqs. ( 94) and (95) take the form The stationary solution z = z 0 = const.,Ω = Ω 0 = const.corresponds to In the weakly damped regime (which, by the results of Sec.V, is that one in which we expect a self-oscillating dipole) we have that where I c is the physical value of the critical Josephson current [see Eq. ( 55)].Therefore, combining Eqs. ( 96) and (102), we find that the radiation efficiency of the weakly damped JJ is:

B. Radiation into open space
For a single JJ that radiates microwaves into open space, we can use the expression for γ e computed for a two-level atom with an electric dipole d.This is given by the formula where now and where the dipole moment of a "JJ atom" is given by the product of the distance ℓ between the superconducting electrodes and the charge 2e for a single Cooper pair.This leads to the final formula for η rad , valid in the weak damping regime: Typically, radiation is observed from JJ's in which I ≃ I c ≃ 10 −3 A, C ≃ 3 × 10 −13 F, V ≃ 10 −3 V, and ℓ = 10 −9 m, which gives This efficiency can be enhanced by placing the JJ inside a high-quality optical cavity, as we will discuss in Sec.VI C. The resulting estimate of the efficiency is reasonable, but a question remains because the coherence length ξ 0 of the Cooper pairs may be much larger than the separation ℓ between the electrodes [71].The proper interpretation of the dipole moment of the "JJ atom" may therefore require more careful consideration, an issue that we must leave for future research.

C. Cavity enhancement
In the experimental setup of [25,27], the JJ is placed inside a resonant cavity.This leads to much stronger and highly coherent radiation, which the authors of [27] call a "Josephson junction laser".To estimate the efficiency in this case, we use a well-known result of cavity QED that describes enhanced spontaneous emission by an atom in resonance (the "Purcell effect" [72]): where Q is a quality factor of the cavity, L is a cavity linear dimension, and λ A = 2πc/ω A is the wavelength of the emitted radiation.Typically L ≃ λ A , so that the enhancement of the efficiency factor is largely determined by Q. Therefore we expect that which is consistent with the observed enhancement by a few orders of magnitude, relative to Eq. (107).For further treatment of the Purcell effect in cavity QED, see [73] and references therein.

VII. DISCUSSION
Despite having been amply studied and applied since it was first proposed in 1962, the Josephson effect has remained something of a conundrum because of the failure of theorists to treat it consistently as the dynamics of an open system.The JJ's ability to convert DC→AC should not be regarded as a spontaneous breaking of time-translation invariance (in the manner of a "time crystal" [74]), but rather as the dynamics of an engine that cyclically extracts work from an external disequilibrium.This DC→AC conversion is a self-oscillation, and as such a thermodynamically irreversible process resulting from a feedback between the state of a "tool" (corresponding, in this case, to the coherence z in the macroscopic density matrix of Eq. ( 24)) and the coupling of the engine's working substance (here the Cooper pairs in the JJ) to the external baths.In our model, this feedback is captured by the coupling of the differential equations for the evolution of z and the voltage V of the junction [Eqs.(35) and (36)].
We have shown that the dynamical system characterized by the coupled Eqs. ( 35) and ( 36) has a rich structure, consistent with the key features of the JJ.One of these features is the voltage-current relation of Eq. ( 47), which exhibits ohmic behavior for large voltage bias, as well as negative differential resistance and hysteresis when the damping of the system is small enough.In Sec.IV we recovered the main qualitative results of the RCSJ model, but without invoking a "tilted washboard potential" unbounded from below, or introducing any ad hoc dissipation mechanism extraneous to the JJ's dynamics.Note also that the current-voltage characteristic for the JJ, as shown in Fig. 2(a) displays a "pinched" hysteresis loop (i.e., a closed hysteresis curve that passes through the origin, at zero current and zero voltage), meeting Chua's updated definition of a "memristor" [75].Our model might therefore be useful in clarifying the debate about whether memristors are passive or active devices [76,77].
Our model of the JJ can describe the well known quantum interference used in SQUIDs for the precise measurement of magnetic fluxes, as we showed in Sec.IV D. Another important feature of our treatment of the JJ as an open system is that it accounts for the frequency locking ("Shapiro steps") observed when the JJ is subjected to an external AC voltage driving, as follows directly from the electromagnetic irradiation used in most experimental setups.As we discussed in Sec.IV E, this effect emerges naturally in our description, without having to invoke any accompanying modulation of the external current (as must be done in the RCSJ model).
In Sec.V we found by numerical integration that this same dynamical system exhibits self-oscillations for sufficiently large external DC bias, in the form of "hidden attractors" that can be excited by small but finite perturbations about a linearly stable equilibrium.This self-oscillation of z and V has a fundamental frequency equal to the Josephson frequency Ω = 2eV /ℏ, and can therefore account for the emission of monochromatic photon and phonon radiation observed in many experiments (a radiation which, in the case of photons, may be enhanced by the coupling of the JJ to a resonant cavity).This monochromatic radiation can be understood as a form of work extracted by the JJ engine from the disequilibrium between the two external electronic baths to which it is coupled.This is consistent with the picture of JJ radiation presented in Sec.VI, based on Dicke's theory of superradiance in quantum optics.However, the precise relation between the hidden attractors of Sec.V and the superradiant mechanism discussed in Sec.VI calls for further investigation.
The model that we have presented here is the simplest dynamical description consistent with these key features of the JJ.Further work is needed to improve and clarify the detailed modelling of the JJ, including such effects as the dependence of the tunneling rate K on the junction voltage V and an explicit treatment of the electrostatic interactions of the Cooper pairs.The methods that we have introduced here might also help to clarify the dynamics of the Esaki diode, an active device used to build electrical self-oscillators that, like the JJ, exhibits a negative differential resistance associated with electron tunneling [78].
Beyond the concrete problem of understanding the Josephson effect as an irreversible process, we consider that the present work is also of interest more broadly to the physical theory of quantum thermodynamics as well as to the mathematical theory of dynamical systems.In previous models of quantum engines, the work is extracted by a degree of freedom that is effectively classical or semi-classical [45].In our model of the JJ work is extracted by a "tool" (or "piston") that is a macroscopic but intrinsically quantum object: the coherence z in the density matrix of Eq. (24).This could make it important to developing a better understanding the physics of practical qubits and how the can be maintained against decoherence.
In the theory of dynamical systems, most mathematical models of autonomous engines have been based on Hopf bifurcations and the associated limit cycles.In our model, the JJ's engine dynamics is associated with hidden attractors [19].The precise characterization of these attractors and their relation to other features of our model deserves further study.Following the suggestion made long ago by Fabry and Le Corbeiller [16,17], we believe that a fruitful dialogue between thermodynamics and the theory of differential equations may be the key to a better understanding of active devices [79], including the JJ and other electrical self-oscillators.

FIG. 1 .
FIG. 1.(a) Josephson junction (JJ), composed of two superconducting electrodes separated by a thin insulator.An external current I is applied to the JJ by an external source.(b) In the RCSJ model, this JJ is treated as an idealized junction J obeying the current-phase relation of Eq. (1), in parallel ("shunt") with a resistor R and a capacitor C.

FIG. 3 .
FIG. 3. Hysteretic response of the JJ model of Eq. (45), with α = 4 and a time-dependent total current itot(τ ).(a) A hysteresis loop is seen when itot is first increased and then decreased linearly, with ditot/dτ = ±0.01.The data points obtained by numerical simulation are shown in black.The solid red curve shows the stationary characteristic of Eq. (48).(b) Responses of the Josephson (tunneling) current iJ , the resistive current ires = v, and the capacitive current icap = dv/dτ (such that iJ + ires + icap = itot), during the rising phase (ditot/dτ = 0.01).The value of the critical ic [Eq.(50)] is marked by a horizontal dashed line.(c) Responses of iJ , ires, and icap during the falling phase (ditot/dτ = −0.01).The value of the retrapping ir [Eq.(51)] is marked by a horizontal dashed line.

FIG. 4 .
FIG. 4. Scheme of a superconducting quantum interference device (SQUID).Current between the two superconducting electrodes (1 and 2) may tunnel either through the thin insulator labelled A or the thin insulator labeled B. An external magnetic flux Φ is applied through the loop formed by the two dotted paths.

FIG. 5 .
FIG. 5. Current-voltage characteristic of the Josephson junction in a presence of the oscillating driving field with frequency ω f = 20 and the field voltage v f = 300 for α = 3.The I-V characteristics reveals the common Shapiro's staircase structure, where the grey horizontal lines correspond to v = kω f for integer values of k.

FIG. 6 .
FIG.6.Hidden-attractor self-oscillation around a stable fixed point ζ0 (with Re ζ0 = īS and Im ζ0 = īJ ) of the dynamical system defined by Eqs.(76) and (77), for α = 2.2.(a) The initial point (in green), marked with reference to the characteristic curve (in blue) of Eq. (48).(b) Hidden attractor in the ζ plane.The orbit is plotted for times 500 ≤ τ ≤ 550.The equilibrium point ζ0 is marked in orange.(c) Power spectrum for the oscillation of the tunneling current iJ , computed from the numerical integration for 3500 ≤ τ ≤ 3650.

FIG. 8 .
FIG.8.Approach to a stable fixed point of the dynamical system defined by Eq. (45) for α = 0.8.(a) The initial point, marked in green with reference to the characteristic curve (in blue) of Eq. (48).(b) Orbit in the complex ζ plane for 0 ≤ τ ≤ 550.