The Josephson junction as a quantum engine

Feynman presented a simple model for JJ in [32], based on a two-component MWF for the Cooper-pair condensate:

$$\begin{equation} \Psi = \begin{pmatrix} \psi_1 \\ \psi_2 \\ \end{pmatrix} = \begin{pmatrix} \sqrt{n_1} e^{i\phi_1} \\ \sqrt{n_2} e^{i\phi_2} \\ \end{pmatrix}, \end{equation} \tag{ 4 }$$

where $n_{1,2}$ are Cooper-pair numbers and $\phi_{1,2}$ are phases corresponding to the left-hand and right-hand electrode of the JJ, respectively (see figure 1). The dynamics is given by the Schrödinger equation

$$\begin{equation} \frac{d}{dt}\Psi\left(t\right) = -\frac{i}{\hbar} H \Psi\left(t\right) \end{equation} \tag{ 5 }$$

with Hamiltonian

$$\begin{equation} H = \hbar \begin{pmatrix} U/2 & K \\ K & - U/2 \\ \end{pmatrix}, \end{equation} \tag{ 6 }$$

where $\hbar U$ is the energy difference between a Cooper pair located at the right electrode and one located at the left electrode. The value of U is related to the voltage bias V across the JJ by the formula

$$\begin{equation} U = \frac{2e V}{\hbar}, \end{equation} \tag{ 7 }$$

while K is the rate of tunneling of Cooper pairs across the barrier between the two electrodes.

Zoom In Zoom Out Reset image size

In this treatment, the tunneling current is obtained from the unitary evolution according to equations (5)–(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 [32]. 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.

Let I be the external current applied to the JJ, as shown in figure 1. In the RCSJ model we have

$$\begin{equation} I = I_c \sin \phi + C \dot V + \frac{V}{R}, \end{equation} \tag{ 8 }$$

where the first term in the right-hand side of equation (8) corresponds to the Josephson current of equation (1). This current is taken to run in parallel ('shunt') with a current $C \dot V$ through capacitor C, and also to a current $V/R$ through resistor R; see, e.g. [33, 34]. By using equation (2) to eliminate V in equation (8), we get the equation of motion

$$\begin{equation} \ddot \phi + \frac{\dot \phi}{R C} + \frac{2 e I_c}{\hbar C} \sin \phi = \frac{2 e I}{\hbar C} ~. \end{equation} \tag{ 9 }$$

In the $R \to \infty$ limit, equation (9) can be obtained from the Hamiltonian

$$\begin{equation} H\left(q,\phi\right) = \frac{e}{\hbar C} \cdot q^2 - I \phi - I_c \cos \phi, \end{equation} \tag{ 10 }$$

where the charge q separated across the terminals of the junction is taken as a canonically conjugate variable to the phase φ. However, equation (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'

$$\begin{equation} U\left(\phi\right) = - I \phi - I_c \cos \phi \end{equation} \tag{ 11 }$$

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. [35]). 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.

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 help 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

$$\begin{equation} \mathrm{d} U = \delta Q - \delta W + \mu \, \mathrm{d} N, \end{equation} \tag{ 12 }$$

where $\delta Q$ is the infinitesimal heat absorbed by the substance from its surroundings, and $\delta 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 $\mu \, \mathrm{d} N$ term as 'chemical work', but here we want to stress the distinction between this term and the actual mechanical work $\delta 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 $\mu \, \mathrm{d} N$ 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:

$$\begin{equation} \Delta U = \oint \left( \delta Q - \delta W + \mu \, \mathrm{d} N \right) = 0 . \end{equation} \tag{ 13 }$$

The net work is, therefore,

$$\begin{equation} W = \oint \delta W = \oint \left( \delta Q + \mu \, \mathrm{d} N\right) . \end{equation} \tag{ 14 }$$

By Clausius's theorem, the entropy produced by the engine over one cycle is

$$\begin{equation} \Sigma = - \oint \frac{\delta Q}{T} = - \oint \frac{\delta Q}{\bar T \left( 1 + T_d / \bar T \right)} \unicode{x2A7E} 0 ~, \end{equation} \tag{ 15 }$$

where T is the temperature of the source from which the working substance receives heat $\delta Q$ at each instant during the cycle, $\bar T$ is the mean value of the temperature over the cycle's full period, and $T_d \equiv T - \bar T$. Let µ be the instantaneous chemical potential and $\bar \mu$ its mean, with $\mu_d \equiv \mu - \bar \mu$. Combining equations (14) and (15) we obtain

$$\begin{equation} W \unicode{x2A7D} \oint \delta Q \left( 1 - \frac{1}{1 + T_d / \bar T} \right) + \oint \mu \, \mathrm{d} N = \oint \frac{\delta Q \cdot T_d}{\bar T + T_d} + \oint \mu_d \, \mathrm{d} N ~. \end{equation} \tag{ 16 }$$

This bound on W is achieved by a reversible cycle ($\Sigma = 0$). The result of equation (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 equation (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 \gt 0$) during the part of the cycle during which it is more strongly coupled to the hotter bath (from which it takes heat $\delta Q \gt 0$). Conversely, the working substance must be at a lower temperature ($T_d \lt 0$) when it is more coupled to the colder bath (into which it rejects heat $\delta Q \lt 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 $\delta Q$) [18]. In [18] this general result is called the 'Rayleigh-Eddington criterion' ⁴ .

It is also evident from equation (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, equation (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), equation (16) reduces to

$$\begin{equation} W \unicode{x2A7D} \oint \mu_d \, \mathrm{d}N . \end{equation} \tag{ 17 }$$

Note that equation (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 ($\mu_d \gt 0$ when dN > 0) and exits at a lower chemical potential ($\mu_d \lt 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, 36], 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 \to -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 section 2.3 of [42]) 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 (equation (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(\phi)$ that is unbounded from below, as in equation (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 equation (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 [27]. 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. [43–45]) 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.

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^\dagger a$ for the number of Cooper pairs. This number can be increased by the fusion of two electrons from the bath into a Cooper pair, or decreased by the reverse process of Cooper-pair decomposition. The effective Hamiltonian describing such processes takes the form

$$\begin{equation} H = E_0 a^\dagger a + \sum_k \epsilon_k b_k^\dagger b_k + \sum_{k,k^{^{\prime}}} \left( g_{kk^{^{\prime}}} a b_k^\dagger b_{k^{^{\prime}}}^\dagger + \bar g_{kk^{^{\prime}}} a^\dagger b_{k^{^{\prime}}} b_k \right) \end{equation} \tag{ 18 }$$

where $b_k^\dagger$ and $b_k$ are, respectively, creation and annihilation operators for electrons in the bath. The interaction Hamiltonian of equation (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 [31]. On the physics of Andreev scattering, see also [46] and references therein. The results of this paper will not depend on the detailed form of the $g_{kk^{^{\prime}}}$ form factors in equation (18).

Under the standard assumption of separation of time scales of the system (Cooper pairs) and baths (electrons), one obtains the following MME for the reduced density matrix of Cooper pairs (harmonic oscillator):

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} \rho = - \frac{i}{\hbar} E \left[a^\dagger a, \rho \right] + \frac{1}{2} \gamma_\downarrow \left( \left[ a, \rho a^\dagger \right] + \left[ a \rho, a^\dagger \right] \right) + \frac{1}{2} \gamma_\uparrow \left( \left[ a^\dagger, \rho a \right] + \left[ a^\dagger \rho, a \right] \right) \end{equation} \tag{ 19 }$$

where E is the renormalized value of ground state energy and $\gamma_{\downarrow} \gt \gamma_{\uparrow}$ 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. [47]). Using equation (19) one can find the exact kinetic equation for the number of Cooper pairs $n(t) = \mathrm{Tr}[\rho(t) a^{\dagger}a]$:

$$\begin{equation} \dot n = - \gamma \left(n - \bar n\right), \quad \gamma \equiv \gamma_{\downarrow} - \gamma_{\uparrow}, \quad \bar n \equiv \frac{\gamma_{\uparrow}}{\gamma}. \end{equation} \tag{ 20 }$$

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 equation (2), the tunneling rate K, and the damping rate γ (we will define these carefully in section 3.2). In a common experimental setup, all of these are ${\sim}10^{-12}$ s; see, e.g. [48]. Meanwhile, the time scale for the internal dynamics of the leads can be estimated from the Drude model of electrical conduction to be ${\sim}10^{-14}$ to 10⁻¹⁵ 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.

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 Cooper-pair modes 1 and 2 corresponding to the two electrodes in figure 1. The system is then equivalent to two quantum harmonic oscillators, with reduced dynamics given by an extension of equation (19), to wit:

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} \rho = - \frac{i}{\hbar} \left[ H_{12}, \rho\right] + {\cal L}_1 \rho + {\cal L}_2 \rho . \end{equation} \tag{ 21 }$$

The Hamiltonian H₁₂ takes the form

$$\begin{equation} H_{12} = E_1 a_1^\dagger a_1 + E_2 a_2^\dagger a_2 + \hbar K \left( a_1^\dagger a_2 + a_2^\dagger a_1 \right), \end{equation} \tag{ 22 }$$

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 the so-called local master equation [49]. This approach is justified when the perturbation is small compared to the energy scale of the free Hamiltonian [50, 51]. In our case, we treat the tunneling rate K as a small perturbation compared to the electrode energies $E_{1,2}$ in equation (22), which is consistent with usual approach in theory of superconductors.

The dissipative generators ${\cal L}_{1,2}$ in equation (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 equation (19), with

$$\begin{equation} {\cal L}_j \rho = \frac{1}{2} \gamma^{\left(j\right)}_{\downarrow} \left( \left[ a_j, \rho a_j^\dagger \right] + \left[ a_j \rho, a_j^\dagger \right] \right) + \frac{1}{2} \gamma^{\left(j\right)}_{\uparrow} \left( \left[ a_j^\dagger, \rho a_j \right] + \left[ a_j^\dagger \rho, a_j \right] \right) . \end{equation} \tag{ 23 }$$

In order to obtain the generalization of equation (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 equation (4), we will call this new object the MDM. In this case, the MDM is a $2 \times 2$ positively defined matrix σ, given by the following two-point correlations:

$$\begin{equation} \sigma \equiv \left[\sigma_{jk}\right], \quad \sigma_{jk} = \textrm{Tr}\left(\rho a_j^\dagger a_k\right), \quad j,k = 1,2 . \end{equation} \tag{ 24 }$$

A useful parametrization of MDM is given by the populations $n_j = \textrm{Tr}(\rho a_j^\dagger a_j)$ and coherence $z = \textrm{Tr}(\rho a_1^\dagger a_2)$:

$$\begin{equation} \sigma = \begin{pmatrix} n_1 & z \\ z^\ast & n_2 \\ \end{pmatrix}, \quad |z|^2 \unicode{x2A7D} n_1 n_2 . \end{equation} \tag{ 25 }$$

Notice that, for the pure state corresponding to equation (4),

$$\begin{equation} \sigma = |\Psi\rangle\langle \Psi | \quad \Rightarrow \quad z = \sqrt{n_1 n_2} e^{i\phi}, \quad \quad \phi = \phi_1 - \phi_2. \end{equation} \tag{ 26 }$$

Equations (21)–(25) translate into the following evolution equations for the MDM:

$$\begin{align} \dot n_1 & = - \gamma^{\left(1\right)} \left(n_1 - \bar n_1\right) - i K \left(z - z^\ast\right) \nonumber \\ \dot n_2 & = - \gamma^{\left(2\right)} \left(n_2 - \bar n_2\right) + i K \left(z - z^\ast\right)\\ \dot z & = \left[ \frac{i}{\hbar} U - \frac{1}{2} \left(\gamma^{\left(1\right)} +\gamma^{\left(2\right)}\right) \right] z - i K \left(n_1 - n_2\right) \nonumber \end{align} \tag{ 27 }$$

where $U = E_1 - E_2$, $\gamma^{(j)} \equiv \gamma^{(j)}_{\downarrow}- \gamma^{(j)}_{\uparrow}$, and $\bar n_j \equiv \gamma^{(j)}_{\uparrow} /\gamma^{(j)}$ for $j = 1,2$.

Let us assume that the relaxation rates for both electrodes are equal, i.e.

$$\begin{equation} \gamma^{\left(1\right)} = \gamma^{\left(2\right)} = \gamma . \end{equation} \tag{ 28 }$$

Then, introducing new variables for the sum and difference of the expected number of Cooper pairs on each of the JJ's electrodes,

$$\begin{equation} N = n_1 + n_2, \quad n = n_1 - n_2, \end{equation} \tag{ 29 }$$

we rewrite equation (27) as

$$\begin{align} \dot N & = - \gamma \left(N - \bar N\right), \quad \bar N = \bar n_1 + \bar n_2, \end{align} \tag{ 30 }$$

$$\begin{align} \dot n & = - \gamma\left(n - \bar n\right) - i 2K \left(z - z^\ast\right), \quad \bar n = \bar n_1 - \bar n_2, \end{align} \tag{ 31 }$$

$$\begin{align} \dot z & = \left( \frac{i}{\hbar} U - \gamma \right) z - i K n . \end{align} \tag{ 32 }$$

This implies that the average total number of Cooper pairs N relaxes, independently of all the other dynamical variables, to its stationary value ($N = \bar N$), which must correspond physically to an electrically neutral JJ. This justifies our having taken equal relaxation rates for both electrodes in equation (28).

The variable n describes an internal local disequilibrium that causes an electrical double layer to form at the junction, with excess charge $Q = e n$ and associated voltage jump V. Accordingly, we introduce the capacitance of the junction C, such that the voltage is given by the formula:

$$\begin{equation} V = \frac{Q}{C} = - \frac{e n}{C}. \end{equation} \tag{ 33 }$$

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:

$$\begin{equation} U = - 2 e V = \frac{2 e^2 n}{C}. \end{equation} \tag{ 34 }$$

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 \gg |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

$$\begin{eqnarray} \dot V & = & - \gamma V + i \frac{2Ke}{C} \left(z - z^\ast \right) + \frac{I}{C}, \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} \dot z & = & - \left( i \frac{2eV}{\hbar} + \gamma \right) z + i \frac{KC}{e} V, \end{eqnarray} \tag{ 36 }$$

where we have put:

$$\begin{equation} I = - e \gamma \bar n = - e \left( \gamma^{\left(1\right)}_{\uparrow} - \gamma^{\left(2\right)}_{\uparrow} \right) . \end{equation} \tag{ 37 }$$

Finally, we introduce a resistance parameter R such that

$$\begin{equation} \gamma = \frac{1}{RC}, \end{equation} \tag{ 38 }$$

in order to facilitate comparison between our results and those of the RCSJ model. Using this, we may rewrite equation (35) as

$$\begin{equation} I = C \dot V + \frac{V}{R} + I_J, \quad \hbox{where} \quad I_J = 4 e K \, \hbox{Im}~z = 4 e K |z| \sin \phi, \end{equation} \tag{ 39 }$$

where the last term can also be defined, in terms of the Hamiltonian dynamics, as:

$$\begin{equation} I_J = q \dot n_1 = \frac{i q}{\hbar} \textrm{Tr} \left( \left[H_{12}, a_1^\dagger a_1\right] \rho \right) = - 2 q K |z| \sin \phi, \end{equation} \tag{ 40 }$$

with $q = -2e$ the charge of a Cooper pair.

Comparing equation (39) to equation (8) of the RCSJ model, we interpret the quantity I defined by equation (37) as the fixed current that the external source provides to the JJ ⁶ . We also see that the term $V /R$ in equation (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 $4 e K |z| \sin \phi$. For simplicity, we neglect the contribution to the JJ current from quasiparticles. Note that $|z|$ is not fixed, because the $n_{1,2}$ of equations (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 \phi$, for fixed I_c . Note also that equation (36) implies that the standard relation between voltage and phase as given in equation (2) and assumed by the RCSJ model, is not exact. Therefore equation (3) for the AC tunneling current (a quantity that has never been measured directly in experiments), is not applicable to our model.

We see now how the feedback that we discussed in general thermodynamic terms in section 2.3 appears mathematically in our dynamical model: equation (35) for the real-valued V depends on z via the tunneling term proportional to $K \, \hbox{Im}~(z)$. Meanwhile, equation (36) for the complex-valued z depends on V, which controls the resonant frequency $2eV/\hbar$. It is therefore plausible that, in some regime, this non-Hamiltonian dynamical system will self-oscillate: an oscillation of $\hbox{Im}~(z)$ forces V in in equation (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 equation (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 section 5.

Note that I > 0 in equation (35) is analogous to 'population inversion' in the theory of a laser, in that it corresponds to $\bar n_1 \gt \bar n_2$ (see equations (30) and (37)), with the energy of the macroscopic quantum state $|1 \rangle$ greater than that of $|2 \rangle$ by the amount U of equation (34). The same is true for I < 0, because then n < 0 and (by equation (34)) we have that the energy of $|2 \rangle$ is greater than that of $|1 \rangle$. 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 ⁷ .

It will be convenient to transform a set of equations (35) and (36) into dimensionless units. First, we separate the second equation into real and imaginary part of z obtaining

$$\begin{align} \begin{aligned} &C \dot V + \frac{1}{R} V + I_J = I, \hfill \\ &\dot I_J + \frac{1}{RC} I_J + \left( \frac{2 e}{\hbar} I_S - 4 K^2C \right) V = 0,\hfill \\ &\dot I_S + \frac{1}{RC} I_S - \frac{2 e V}{\hbar} I_J = 0, \hfill\\ \end{aligned} \end{align} \tag{ 41 }$$

where we have introduced the variable $I_S = 4 e K \, \textrm{Re}~{z} = 4 e K |z| \cos \phi$, so that

$$\begin{equation} z = \frac{1}{4 e K} \left(I_S + i I_J\right) . \end{equation} \tag{ 42 }$$

Let us also introduce the characteristic voltage and current parameters

$$\begin{equation} \tilde V = \frac{\hbar}{2 e RC}, \quad \tilde I = \frac{\tilde V}{R}, \end{equation} \tag{ 43 }$$

and define new dimensionless variables as:

$$\begin{equation} \tau = \gamma t, \quad v = V/\tilde V, \quad i_\textrm{tot} = I/\tilde I, \quad i_J = I_J/\tilde I, \quad i_S = I_S/\tilde I. \end{equation} \tag{ 44 }$$

The equations of motion for the JJ can now be expressed as

$$\begin{align} \begin{aligned} &\frac{\mathrm{d}}{\mathrm{d} \tau} v + v + i_J = i_\textrm{tot}, \hfill \\ &\frac{\mathrm{d}}{\mathrm{d} \tau} i_J + i_J + \left(i_S -4 \alpha\right) v = 0,\hfill \\ &\frac{\mathrm{d}}{\mathrm{d} \tau} i_S + i_S - i_J v = 0, \hfill\\ \end{aligned} \end{align} \tag{ 45 }$$

where

$$\begin{equation} \alpha \equiv \frac{K^2}{\gamma^2} . \end{equation} \tag{ 46 }$$

Our model is therefore characterized by a single parameter α: the ratio of the squares of the tunneling and dissipation rates.

For α > 2 (i.e. $K \gt \sqrt 2 \gamma$) the function $i_\textrm{tot}(v)$ given in equation (48) is not monotonically increasing. The two positive solutions to $i_\textrm{tot}^{^{\prime}}(v_0) = 0$ are

$$\begin{equation} v_{\pm} = \sqrt{2 \alpha -1 \pm 2 \sqrt{\left(\alpha-2\right) \alpha}}, \end{equation} \tag{ 49 }$$

such that a negative differential resistance, i.e. $i_\textrm{tot}^{^{\prime}}(v) \lt 0$, is observed for voltages $v \in [v_-, v_+]$ or $v \in [-v_+, -v_-]$ (provided that α > 2). The positive local maximum $i_\textrm{tot}(v_-)$ we identify with the critical Josephson current:

$$\begin{equation} i_c \equiv i_\textrm{tot}\left(v_-\right) = \left(1 + \alpha + \sqrt{\left(\alpha-2\right)\alpha}\right) \sqrt{2\alpha -1 - 2 \sqrt{\alpha\left(\alpha-2\right)}}. \end{equation} \tag{ 50 }$$

The positive local minimum $i(v_+)$ (which we will identify with the retrapping Josephson current) is:

$$\begin{equation} i_r \equiv i_\textrm{tot}\left(v_+\right) = \left(1+\alpha - \sqrt{\left(\alpha-2\right)\alpha}\right) \sqrt{2\alpha -1 + 2 \sqrt{\alpha\left(\alpha-2\right)}}. \end{equation} \tag{ 51 }$$

Thus, the critical value α = 2 —for which $i_c = i_r = 3\sqrt{3}$— discriminates between two regimes: with and without the negative differential resistance of JJ. In section 4.2 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 equation (50) defines a critical current only for α > 2 (the underdamped regime). For $\alpha \unicode{x2A7D} 2$ (the overdamped regime), a natural definition of the critical current is given by the inflection point: $i_c = i_\mathrm{tot}(v_0)$ when $i_\mathrm{tot}^{^{\prime\prime}}(v_0) = 0$. This corresponds to $v_0 = \sqrt 3$, so that

$$\begin{equation} i_c = \sqrt{3}\left(1+\alpha\right), \quad \text{for} \quad \alpha \unicode{x2A7D}2. \end{equation} \tag{ 52 }$$

Hence, the expression for critical current for full range of parameter α is given by

$$\begin{align} i_c = \left\{ \begin{array}{ll} \sqrt{3}\left(1+\alpha\right) & \quad \text{for} \quad 0 \unicode{x2A7D} \alpha \unicode{x2A7D} 2 \\ \left(1 + \alpha + \sqrt{\left(\alpha-2\right) \alpha}\right) \sqrt{2 \alpha - 1- 2 \sqrt{\alpha\left(\alpha-2\right)}} & \quad \text{for} \quad \alpha > 2 .\end{array} \right. \end{align} \tag{ 53 }$$

These results are illustrated in figure 2, for two different choices of the parameters α. In figure 2(a), for α > 2, the characteristic $i_\textrm{tot}(v)$ has local maxima and minima, and therefore well-defined values of i_c and i_r . In figure 2(b) the value $\alpha \unicode{x2A7D} 2$ the characteristic $i_\textrm{tot}(v)$ is monotonically increasing and therefore has no i_r . We will show in section 4.2 that the first characteristic exhibits hysteresis, whereas the second does not.

Zoom In Zoom Out Reset image size

Note that the critical current i_c is given by the expression:

$$\begin{equation} i_c = \frac{I_c}{\tilde I} = \frac{2 e R^2 C I_c}{\hbar} \equiv \beta_c, \end{equation} \tag{ 54 }$$

which agrees with the definition of the Stewart–McCumber parameter β_c in the RCSJ model. Thus, our model exhibits hysteresis for $\beta_c \unicode{x2A7E} 3\sqrt{3}$, whereas in the RCSJ model, hysteresis is usually said to correspond to $\beta_c \gtrsim 1$.

In the limit $\alpha \to \infty$ (corresponding to a damping rate γ much smaller than the tunneling rate K), equation (53) implies that $i_c \to 2 \alpha$. In that case the physical value of the critical current I_c is

$$\begin{equation} I_c = \frac{\hbar}{2e R^2 C} 2 \alpha = \frac{\hbar C}{e} K^2 . \end{equation} \tag{ 55 }$$

We will use this limiting value of the critical current in section 6.1, where we study the efficiency of the JJ as an engine.

To determine the stability of fixed points we introduce the small deviations $\delta v$, $\delta i_J$, and $\delta i_S$ of the dynamical variables away from their respective stationary values $v_0, i_{J,0}$, and $i_{S,0}$, such that

$$\begin{equation} v = v_0 + \delta v, \quad i_J = i_{J,0} + \delta i_J, \quad i_S = i_{S,0} + \delta i_S, \quad \text{where} \quad i_{J,0} = \frac{4 \alpha v_0}{1+ v_0^2}, \quad i_{S,0} = \frac{4 \alpha v_0^2}{1+ v_0^2}. \end{equation} \tag{ 56 }$$

Inserting this into equation (45) yields the linearized equation

$$\begin{eqnarray} &&\frac{d}{d\tau} \begin{pmatrix} \delta v \\ &&\delta i_J \\ \delta i_S\end{pmatrix} = \begin{pmatrix} - 1 & - 1 & 0 \\ \frac{4 \alpha}{1+v_0^2} & - 1 & - v_0 \\ \frac{4 \alpha v_0}{1+v_0^2} & v_0 & -1 \end{pmatrix} \begin{pmatrix} \delta v \\ &&\delta i_J \\ &&\delta i_S\end{pmatrix} . \end{eqnarray} \tag{ 57 }$$

The eigenvalues of this linear set are therefore the three roots $\{\lambda_i\}$ of the characteristic polynomial

$$\begin{equation} \lambda^3 + 3 \lambda^2 + \left(v_0^2 + \frac{4 \alpha}{1+ v_0^2} + 3 \right) \lambda + v_0^2 + 4 \alpha \ \frac{1- v_0^2}{1+ v_0^2} + 1 . \end{equation} \tag{ 58 }$$

According to the Routh–Hurwitz criterion, the equilibrium at $\delta v = \delta i_J = \delta i_S = 0$ will be unstable (i.e. at least one root will have $\textrm{Re}~{\lambda_i} \gt 0$) if and only if

$$\begin{equation} v_0^2 + 4 \alpha \ \frac{1- v_0^2}{1+ v_0^2} + 1 < 0 . \end{equation} \tag{ 59 }$$

Note that

$$\begin{equation} v_0^2 + 4 \alpha \, \frac{1-v_0^2}{1+ v_0^2} + 1 = \left( 1+ v_0^2 \right) \frac{d}{d v} \left[v \left( 1 + \frac{4 \alpha}{1 + v^2} \right) \right]_{v = v_0} = \left(1+v_0^2 \right) \, i^{^{\prime}}_\textrm{tot}\left(v_0\right) , \end{equation} \tag{ 60 }$$

such that the unstable regime given by condition (59) precisely corresponds with negative differential resistance, $i_\textrm{tot}^{^{\prime}}(v_0) \lt 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 (equation (50)) and from below by i_r (equation (51)). If the applied current $i_\textrm{tot}$ is initially zero and is then increased slowly, the resulting voltage will follow the nearly vertical segment of the characteristic (see figure 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 as 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_\textrm{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_\textrm{tot} = i_r$ is reached, whereupon the current will stay almost fixed as 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_\textrm{tot} = i_r$).

We have verified this hysteretic behavior in numerical integrations of equation (45) with a time-dependent bias current $i_\textrm{tot} = i_\textrm{tot} (\tau)$. In particular, we let $i_\textrm{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 figure 3 for α = 4. The hysteresis loop appears clearly in figure 3(a).

Zoom In Zoom Out Reset image size

We may also study the responses of the Josephson (i.e. tunneling) current i_J , the resistive current $i_\textrm{res} = v$, and the capacitive current $i_\textrm{cap} = \mathrm{d}v/\mathrm{d}\tau$, such that $i_J + i_\textrm{res} + i_\textrm{cap} = i_\textrm{tot}$. In figure 3(b) we see a sudden jump in $i_\textrm{res}$ shortly after $i_\textrm{tot}$ surpasses the critical value i_c . In figure 3(c) we see how $i_\textrm{res}$ decreases rapidly shortly after $i_\textrm{tot}$ falls below the retrapping i_r . These jumps in $i_\textrm{res}$ are accompanied by changes in i_J in the opposite direction, as well as by brief spikes in $i_\textrm{cap}$.

Much like in the 'putt-putt' engine model of [18, 36] and in the treatment of the electronic Schmitt trigger in [42], ours is a non-Hamiltonian dynamical system in which an explicitly irreversible 'kinetic equation' (equation (35)) is coupled to a 'mechanical' equation of motion for an oscillator (equation (36)). Positive feedback between these equations can lead to a self-oscillation, which in the models of [18, 36, 42] appears mathematically as a Hopf bifurcation [54]. It is therefore natural to ask whether something similar happens in our model of the JJ.

Of the three roots of the polynomial in equation (58), one is real (we label it as λ₀) and the other two are conjugate:

$$\begin{equation} \lambda_{\pm} = \kappa \pm i \eta, \quad \kappa, \eta \in \mathbb{R}. \end{equation} \tag{ 61 }$$

It can be shown that for α > 0 the real part of the complex roots is non-positive, i.e. $\kappa \unicode{x2A7D} 0$. This implies that a bifurcation appears only for the real λ₀, 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 section 5. There we will show that our model does exhibit self-oscillation, 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.

Our model can also be applied to interference phenomena. Consider the superconducting quantum interference device (SQUID) represented in figure 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

$$\begin{equation} K = K_A \exp \left( -i\frac{2e}{\hbar} \int_A \textbf{A} \cdot \mathrm{d}\textbf{s} \right) +K_B \exp \left( -i\frac{2e}{\hbar} \int_B \textbf{A} \cdot \mathrm{d}\textbf{s} \right), \end{equation} \tag{ 62 }$$

where $K_{A,B}$ are the tunneling rates and $\int_{A,B}$ are the integrals along the paths passing through insulators A and B, respectively. The two terms in the right-hand side of equation (62) correspond to two possible tunneling paths, associated with the two junctions.

Zoom In Zoom Out Reset image size

So far we have considered a real tunneling constant. For complex K, the tunneling part of the Hamiltonian (22) should be rewritten as

$$\begin{equation} \hbar K\,^\ast a_1^\dagger a_2 + \hbar K a_2^\dagger a_1 . \end{equation} \tag{ 63 }$$

We then obtain the following equations:

$$\begin{align} \dot n & = - \gamma\left(n - \bar n\right) - i 2K\,^\ast z + i2K z^\ast, \quad \bar n = \bar n_1 - \bar n_2, \end{align} \tag{ 64 }$$

$$\begin{align} \dot z & = \left( \frac{i}{\hbar} U - \gamma \right) z - i K n . \end{align} \tag{ 65 }$$

It follows that the I–V characteristics has the same form as in equation (48), with $\alpha = K^2 / \gamma^2$ replaced with $|K|^2 / \gamma^2$. The expressions for the maximal and retrapping currents are also the same as before, with K replaced by $|K|$, where

$$\begin{eqnarray} &&|K| = \left| K_A \exp \left( -i\frac{2e}{\hbar} \int_A \textbf{A} \cdot \mathrm{d}\textbf{s} \right) +K_B \exp \left( -i\frac{2e}{\hbar} \int_B \textbf{A} \cdot \mathrm{d}\textbf{s} \right) \right| = \left| K_A +K_B \exp \left( -i\frac{2e}{\hbar} \oint \textbf{A} \cdot \mathrm{d}\textbf{s} \right) \right| . \end{eqnarray} \tag{ 66 }$$

It we write the magnetic flux through the loop as $\Phi = \oint \textbf{A} \cdot \mathrm{d}\textbf{s}$, this becomes

$$\begin{equation} |K|^2 = |K_A|^2 + |K_B|^2 + 2 |K_A K_B| \cos \left(\frac{2e}{\hbar} \Phi \right), \end{equation} \tag{ 67 }$$

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 Φ.

Josephson had predicted in [1] that if a JJ were driven by an external oscillating electric field with frequency $\Omega_f$, then the I–V characteristic would be modified, with steps appearing at values of the bias voltage V₀ such that

$$\begin{equation} \omega_0 \equiv \frac{2 e V_0}{\hbar} = n \Omega_f \end{equation} \tag{ 68 }$$

for $n \in \mathbb{Z}$. This effect was found experimentally by Shapiro [3], and is now used as a metrological standard for the determination of voltages [55].

The usual explanation of this effect (see, e.g. [33]) is that a voltage bias

$$\begin{equation} V = V_0 + V_1 \cos \Omega_f t \end{equation} \tag{ 69 }$$

causes, by equations (1) and (2), a tunneling current

$$\begin{equation} I_J \left(t\right) = I_c \sin \left( \phi_0 + \omega_0 t + \frac{2 e V_1}{\hbar \Omega_f} \sin \Omega_f t \right) . \end{equation} \tag{ 70 }$$

The sine of the sine can be expanded using Bessel functions, so that

$$\begin{equation} I_J \left(t\right) = I_c \sum_{n = -\infty}^\infty \left(-1\right)^n J_n \left( \frac{2 e V_1}{\hbar \Omega_f} \right) \sin \left( \phi_0 + \omega_0 t - n \Omega_f t\right) . \end{equation} \tag{ 71 }$$

The terms of the expansion in equation (71) all average out to zero over long times unless $\omega_0 = n \Omega_f$ exactly, in which case a DC contribution

$$\begin{equation} \bar I_J = \left(-1\right)^n J_n \left( \frac{2 e V_1}{\hbar \Omega_f} \right) \sin \phi_0 \end{equation} \tag{ 72 }$$

would appear. Although this argument agrees with the observed locations of the Shapiro steps, it fails to account for their widths [33]. Various elaborations of the theory have therefore been proposed, based on phenomena such as frequency mixing [56], resistive feedback [57], or non-linear resonance [58]. In general, these models assume that the external current I applied to the JJ also oscillates with frequency $\Omega_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:

$$\begin{equation} H_{12} = E_1 a_1^\dagger a_1 + E _2 a_2^\dagger a_2 + \hbar K \left( a_1^\dagger a_2 + a_2^\dagger a_1 \right) + e d \varepsilon \cos\left(\Omega_f t\right) \left( a_1^\dagger a_1 - a_2^\dagger a_2 \right). \end{equation} \tag{ 73 }$$

where $\varepsilon \cos(\Omega_f t)$ is a time-dependent external field with the amplitude ε and frequency $\Omega_f$, while $e d (a_1^\dagger a_1 - a_2^\dagger a_2)$ is an electrical dipole moment of JJ. This gives rise to the following inhomogeneous extension of equation (45):

$$\begin{align} \begin{aligned} &\frac{\mathrm{d}}{\mathrm{d} \tau} v + v + i_J = i_\textrm{tot}, \hfill \\ &\frac{\mathrm{d}}{\mathrm{d} \tau} i_J + i_J + i_S \left[v + v_f \cos\left(\omega_f \tau\right)\right] - 4 \alpha v = 0,\hfill \\ &\frac{\mathrm{d}}{\mathrm{d} \tau} i_S + i_S - i_J \left[v + v_f \cos\left(\omega_f \tau\right)\right] = 0, \hfill\\ \end{aligned} \end{align} \tag{ 74 }$$

where

$$\begin{equation} v_f = \frac{d \varepsilon}{\tilde V} \quad \text{and} \quad \omega_f = \frac{\Omega_f}{\gamma}. \end{equation} \tag{ 75 }$$

Figure 5 shows the I–V characteristic of the JJ in the presence of an oscillating electric field, with the voltage averaged over many periods of the applied field. The standard Shapiro-step structure is observed, with steps appearing at values given by equation (68), which in terms of the dimensionless variables introduced above is simply $v = k \omega_f$, for integer $k$.

Zoom In Zoom Out Reset image size

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 equation (45) may exhibit hidden limit cycles. The equations of motion for the system can be expressed as

$$\begin{eqnarray} & \dot v + v + \frac{1}{2i}\left(\zeta - \zeta^*\right) = i_\textrm{tot}, \end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray} & \dot \zeta + \zeta + i \left(\zeta - 4\alpha\right) v = 0 . \end{eqnarray} \tag{ 77 }$$

where ζ is a complex variable such that $\textrm{Re}~ \zeta = i_S$ and $\textrm{Im}~\zeta = i_J$, while the dot indicates derivation with respect to the dimensionless time variable τ (see equations (43)–(45)). In terms of the Fourier expansions

$$\begin{eqnarray} v\left(\tau\right) & = & \sum_n v_n e^{-in\omega \tau}, \quad v_{-n} = v_n^* \nonumber \\ \zeta \left(\tau\right) & = &\sum_n \zeta_n e^{-in\omega \tau}, \end{eqnarray} \tag{ 78 }$$

and taking $i_\textrm{tot}$ = const., equations (76) and (77) become

$$\begin{align} &\left(1 - i n \omega\right) v_n + \frac{1}{2i} \left(\zeta_n - \zeta_{-n}^*\right) - i_\textrm{tot} \delta_{n,0} = 0, \end{align} \tag{ 79 }$$

$$\begin{align} &\left(1 - i n \omega\right) \zeta_n + i \sum_m \zeta_{n-m} v_m - 4i \alpha v_n = 0. \end{align} \tag{ 80 }$$

We expect the hidden-attractor oscillation to be dominated by the fundamental frequency ω. Taking equation (79) for n = 0 and equation (80) for n = 1, and neglecting $v_n, \zeta_n$ for all $|n| \gt 1$, we obtain

$$\begin{align} & v_0 + \frac{1}{2i} \left( \zeta_0 - \zeta_0^\ast \right) \simeq i_\textrm{tot}, \end{align} \tag{ 81 }$$

$$\begin{align} &\left[1 + i \left(v_0 - \omega\right) \right] \zeta_1 \simeq i \left(4 \alpha - \zeta_0 \right) v_1 . \end{align} \tag{ 82 }$$

The equilibrium value of ζ is, according to equation (77),

$$\begin{equation} \zeta_0 = \frac{4 \alpha}{1 - \frac{i}{v_0}}, \end{equation} \tag{ 83 }$$

where v₀ is the corresponding equilibrium value for v (the same quantity as in equation (56)). For $|v_0| \gg 1$ we therefore have

$$\begin{equation} \zeta_0 \simeq 4 \alpha . \end{equation} \tag{ 84 }$$

Putting equation (84) into equation (81) we get

$$\begin{equation} v_0 \simeq i_\textrm{tot}, \end{equation} \tag{ 85 }$$

i.e. v₀ 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_\textrm{tot}| \gg 1$) and small damping ($4 \alpha \gg 1$), equations (81) and (82) are consistent with a hidden-attractor self-oscillation ($|v_1| \gt 0$ and $|\zeta_1| \gt 0$) with a fundamental frequency

$$\begin{equation} \omega \simeq v_0 \simeq i_\textrm{tot} ~. \end{equation} \tag{ 86 }$$

Note that in the original variables (see equations (43) and (44)), this fundamental frequency is

$$\begin{eqnarray} &&\Omega = \gamma \omega = \frac{2e V_0}{\hbar} = \frac{2e R I}{\hbar}. \end{eqnarray} \tag{ 87 }$$

This agrees with what we expected from the nature of parametric resonance, as was explained in words in section 3.2, after we formulated equations (35) and (36). In section 5.2 we discuss how such hidden-attractor oscillations can be seen in numerical simulations. There we will measure frequency in terms of an equivalent v₀.

As we already alluded to in section 4.3, the mathematical appearance of such a self-oscillation is different from the Hopf bifurcations considered in [8, 54] and in the 'putt-putt' engine model of [18, 36]. The form of equations (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 equations (76) and (77) numerically, for various values of the parameter α and for various choices of initial conditions. The results 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 section 5.1, we find numerical evidence of hidden attractors around stable equilibrium points for sufficiently large values of α and $\bar v$. To find them, we displaced the initial conditions slightly away from the equilibrium $\delta v = \delta i_J = \delta i_S = 0$ (see equation (56)) and we then observed the corresponding orbit after long times. Figure 6 shows one such hidden attractor, for an orbit with initial values $\bar v = 30$, $\delta v = \delta i_J = 0$, and $\delta i_S = -0.1$ In general, we find that hidden attractors appear only for values of α large enough that the corresponding I-V characteristic exhibits hysteresis (see section 4.1). Since the equilibrium point ζ₀ inside the hidden attractor is stable, it has a finite basin of attraction. For the same parameters used in the simulation of figure 6 but with initial $|\delta i_S| \lesssim 10^{-6}$ (while keeping the initial $\delta v = \delta i_J = 0$), we find that the resulting orbit decays to ζ₀, rather than giving a persistent oscillation.

Zoom In Zoom Out Reset image size

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₀ (in the dimensionless variables defined in equation (44)). This relation between the fundamental frequency of the hidden-attractor oscillation and the applied voltage v₀ holds over a range of values of v₀, as shown in figure 7. This result is consistent with the fact that the electromagnetic and phonon radiation from a DC-biased JJ has fundamental frequency $\Omega = 2eV/\hbar$. Note, however, that in our model the simple formula of equation (3) does not hold. The conversion between the applied DC bias and the observed AC emission (what we called 'DC → AC' in section 1) 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.

Zoom In Zoom Out Reset image size

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 figure 8. A more precise and extensive characterization of the dynamical system defined by equation (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.

Zoom In Zoom Out Reset image size

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 section 4.2) and of its response to external driving by an AC electric field (as explained in section 4.5), 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 figure 9, for the same initial conditions as in figure 6, the SUNDIALS library [63] (which uses the Radau method) approaches the stable fixed point.

Zoom In Zoom Out Reset image size

Comparing the LSODA and Radau trajectories shown in figure 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 figure 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 equation (45) and in the inhomogeneous system with an external periodic driving of equation (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 equation (9). However, as we stressed in section 2.3, 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 equation (74) and whether this might have something to do with the structure of the Shapiro steps that we obtained in section 4.5.

So far we have not identified any chaotic regimes for equation (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]).

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². If the initial state of the atomic ensemble is a product of N identical $2 \times 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:

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\sigma = - i \omega_A \left[S_3, \sigma\right] + \frac{1}{2} \gamma_e \left[ \textrm{Tr} \left(\sigma S^+\right) S^- - \textrm{Tr} \left(\sigma S^-\right) S^+, \sigma\right] . \end{equation} \tag{ 88 }$$

For details on this formulation, see [47].

Here, ω_A is an atomic frequency and γ_e is the spontaneous emission rate for the single atom. Denoting by $|1\rangle, |2\rangle$ the excited and ground state respectively, we define the 'spin' matrices by

$$\begin{equation} S_3 |1\rangle = \frac{1}{2} |1 \rangle, \quad S_3|2 \rangle = -\frac{1}{2} |2 \rangle, \quad S^+ |2 \rangle = |1 \rangle, \quad \hbox{and} \quad S^- |1\rangle = |2 \rangle . \end{equation} \tag{ 89 }$$

Inserting, the same parametrization of the MDM in terms of $n_{1,2}$ and z as in equation (25) into equation (88), we obtain

$$\begin{eqnarray} \dot n_1 & = & - \gamma_e |z|^2, \quad \dot n_2 = \gamma_e |z|^2, \nonumber \\ \dot z & = & - i \omega_A z + \frac{1}{2} \gamma_e\left(n_1- n_2\right) z . \end{eqnarray} \tag{ 90 }$$

Because the total number of atoms $N = n_1 + n_2$ is preserved, we can rewrite this in terms of the variable $n = n_1 - n_2$ as

$$\begin{eqnarray} \dot n & = & - 2\gamma_e |z|^2, \nonumber \\ \dot z & = & - i \omega_A z - \frac{1}{2} \gamma_e n z . \end{eqnarray} \tag{ 91 }$$

Note that, due to the mathematical structure of equation (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

$$\begin{equation} \dot n_1 = - \gamma_e \left(N- n_1\right) n_1 . \end{equation} \tag{ 92 }$$

This reproduces the N² proportionality and the bell-shaped time dependence of the intensity of radiated energy seen in superradiance. The radiated power is given by

$$\begin{equation} {\cal P}\left(t\right) = -\hbar \omega_A \dot n_1 = \hbar \gamma_e \omega_A |z\left(t\right)|^2 . \end{equation} \tag{ 93 }$$

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 for describing 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 (equations (35) and (36)) the superradiant contributions resulting of equation (91), recalling that $V = e n/C$ (equation (33)). The modified JJ dynamics is then governed by

$$\begin{eqnarray} \dot V & = &- \gamma V + i \frac{2Ke}{C} \left(z - z^\ast\right) + \frac{1}{C} \left(I - 2e\gamma_e |z|^2 \right), \end{eqnarray} \tag{ 94 }$$

$$\begin{eqnarray} \dot z & = & -\left( \gamma + \gamma_e \frac{CV}{2e} \right) z - i \frac{2 e V}{\hbar} z + i K \frac{CV}{e} . \end{eqnarray} \tag{ 95 }$$

The correction $-2e\gamma_e |z|^2$ to the external current I in equation (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:

$$\begin{equation} \eta_\textrm{rad} = \frac{\mathcal{P}_\textrm{rad}}{\mathcal{P}_\textrm{in}} = \frac{2e V \gamma_e |z_0|^2}{VI} = \frac{2e \gamma_e |z_0|^2}{I}, \end{equation} \tag{ 96 }$$

where z₀ is the stationary value of z in equations (94) and (95). Note that the radiated power $2eV \gamma_e |z_0|^2$ computed from equation (94) agrees with the expression of equation (93) if we take an 'atomic frequency' $\omega_A = 2 eV / \hbar$.

Rather than using the dimensionless variables that we introduced in section 3 (see equation (44)), here it will be more convenient to work with the angular frequency

$$\begin{equation} \Omega = \frac{2eV}{\hbar} \end{equation} \tag{ 97 }$$

and the Cooper-pair number current

$$\begin{equation} J \equiv \frac{I}{2e} . \end{equation} \tag{ 98 }$$

Introducing the frequency parameter

$$\begin{equation} \omega_c \equiv \frac{4e^2}{\hbar C}, \end{equation} \tag{ 99 }$$

Equations (94) and (95) take the form

$$\begin{eqnarray} \dot \Omega & = & - \gamma \Omega + i {K \omega_c} \left(z - z^\ast\right) + \omega_c \left( J - \gamma_e |z|^2 \right) \nonumber \\ \dot z & = & - \left( i \Omega + \gamma + \gamma_e \frac{\Omega}{\omega_c} \right) z +i\frac{2K}{\omega_c} \Omega. \end{eqnarray} \tag{ 100 }$$

The stationary solution $z = z_0 =$ const., $\Omega = \Omega_0 =$ const. corresponds to

$$\begin{equation} z_0 = \frac{2K}{\omega_c} \left[ \frac{\Omega_0^2 + i \left(\gamma+ \gamma_e \Omega_0 / \omega_c \right) \Omega_0}{\Omega_0^2 + \left( \gamma + \gamma_e \Omega_0 / \omega_c \right)^2} \right] . \end{equation} \tag{ 101 }$$

In the weakly damped regime (which, by the results of section 5, is that one in which we expect a self-oscillating dipole) we have that

$$\begin{equation} |z_0|^2 \simeq \frac{4 K^2}{\omega_c^2} \simeq \frac{I_c}{e \omega_c}, \end{equation} \tag{ 102 }$$

where I_c is the physical value of the critical Josephson current (see equation (55)). Therefore, combining equations (96) and (102), we find that the radiation efficiency of the weakly damped JJ is:

$$\begin{equation} \eta_\textrm{rad} = \frac{2 \gamma_e}{\omega_c} \frac{I_c}{I} = \gamma_e \frac{\hbar C}{2 e^2} \frac{I_c}{I} . \end{equation} \tag{ 103 }$$

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

$$\begin{equation} \gamma_e = \frac{4\pi\mu_0}{3\hbar c}\mathbf{d}^2 \omega_A^3 . \end{equation} \tag{ 104 }$$

where now

$$\begin{equation} \omega_A = \frac{2eV}{\hbar}, \end{equation} \tag{ 105 }$$

and where the dipole moment of a 'JJ atom' is given by the product of the distance $\ell$ 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:

$$\begin{equation} \eta_\textrm{rad} = \frac{I_c}{I} \frac{64 \pi \mu_0 e^3}{3 \hbar^3 c} C V^3 \ell^2 . \end{equation} \tag{ 106 }$$

Typically, radiation is observed from JJ's in which $I \simeq I_c \simeq 10^{-3}$ A, $C \simeq 3 \times 10^{-13}$ F, $V \simeq 10^{-3}$ V, and $\ell = 10^{-9}$ m, which gives

$$\begin{equation} \eta_\textrm{rad} \simeq 10^{-7} . \end{equation} \tag{ 107 }$$

This efficiency can be enhanced by placing the JJ inside a high-quality optical cavity, as we will discuss in section 6.3. The resulting estimate of the efficiency is reasonable, but a question remains because the coherence length ξ₀ of the Cooper pairs may be much larger than the separation $\ell$ between the electrodes [71]. The proper interpretation of the dipole moment of the 'JJ atom' could therefore require more careful consideration, an issue that we must leave for future research.

In the experimental setup of [22, 24], the JJ is placed inside a resonant cavity. This leads to much stronger and highly coherent radiation, which the authors of [24] call a 'Josephson junction laser'. To estimate the efficiency in this case, we use a well-known result of cavity quantum electrodynamics (QED) that describes enhanced spontaneous emission by an atom in resonance (the 'Purcell effect' [72]):

$$\begin{equation} \gamma_e^\textrm{cav} = \gamma_e \frac{3Q}{4\pi^2} \left( \frac{\lambda_A}{L}\right)^3, \end{equation} \tag{ 108 }$$

where Q is a quality factor of the cavity, L is a cavity linear dimension, and $\lambda_A = 2\pi c/ \omega_A$ is the wavelength of the emitted radiation. Typically $L \simeq\lambda_A$, so that the enhancement of the efficiency factor is largely determined by Q. Therefore we expect that

$$\begin{equation} \eta_\textrm{rad}^\textrm{cav} \simeq 0.1 \, \eta_\textrm{rad} Q, \end{equation} \tag{ 109 }$$

which is consistent with the observed enhancement by a few orders of magnitude, relative to equation (107). For further treatment of the Purcell effect in cavity QED, see [73] and references therein.