HOME

Natural oscillation period [frames]
Bias period [frames]
RC time [frames]

We begin with the resistively and capacitively shunted Josephson junction with resistance R in ohms, capacitance C in farads, and critical current in amperes $I_0$. The current bias is I(t), also in amperes, and is some general function of time which comes from the input signal. The first Josephson equation is:

$$ I(t) = I_0\sin{\delta}, $$

where $\delta$ is the difference in phase between the wave function on each side of the tunnel barrier, and it is a function of time. The second Josephson equation is that the voltage across the junction is

$$ V = \frac{\Phi_0}{2\pi}\dot{\delta} $$

$\Phi_0$ is the flux quantum, which is $h/2e$, or about $2.07\times 10^{-15}$ webers, and $\dot{\delta}$ is the time derivative of the phase difference in radians.

The full equation of motion with a resistor and capacitor is

$$ C\dot{V} + \frac{V}{R} + I_0\sin{\delta} = I(t) $$

From this we can subsitute the second Josephson equation above to get

$$ \frac{\Phi_0}{2\pi}C\ddot{\delta} + \frac{\Phi_0}{2\pi}\frac{1}{R}\dot{\delta} = I(t) - I_0\sin{\delta} $$

now we define the Josephson inductance to be

$$ L_J \equiv \frac{\Phi_0}{2\pi}\frac{1}{I_0}, $$

Define the normalized current bias to be $i(t) = I(t)/I_0$, and divide each side of the equation by the critical current to get

$$ L_JC\ddot{\delta} + \frac{L_J}{R}\dot{\delta} = i(t) - \sin{\delta} $$

Now we define the plasma frequency to be $$ \omega_p = \frac{1}{\sqrt{L_J C}} $$ and multiply both sides of the equation by $\omega_p^2$ to get:

$$ \ddot{\delta} + \frac{1}{RC}\dot{\delta} = \omega_p^2(i(t) - \sin{\delta}) $$

And now we note that the right side of this equation can be seen as position derivative of an effective potential $U(\delta)$ which is

$$ U(\delta) = E_J(1 - \cos{\delta} - i(t)\delta) $$

This is enough information to plot the oscillating potential energy surface. The normalized current bias is an AC bias with period of a few frames in our animation. Our animations are generally 10 to 30 frames, so the period should be a few frames. The period in frames is $T_{bias}$

$$ i(t) = A\cos{\left(\frac{2\pi \textrm{[frameIndex]}}{T_{bias}}\right)} $$

We vary the position x from 0 on the left side of the square plot to $[\textrm{squaresize}]$ pixels. The phase $\delta$ varies across the horizontal axis with a period defined by $T_{phase}$ in units of pixels. We want to center the plot on the center of the square, at the x value of $[\textrm{squresize}]/2$, and so the value of $\delta$ across the square as a function of x is $$ \delta = \frac{2\pi}{T_{phase}}\left(x - \frac{[\textrm{squaresize}]}{2}\right) $$ and the y value will be $$ y = \frac{[\textrm{squaresize}]}{2} - E_J(1 - \cos{\delta} - i(t)\delta), $$ where $E_J$ is the Josephson energy in units of pixels. Note that $\delta$ is always in units of radians.

position mouse to set the amplitude of oscillation(X) and the josephson energy(Y). hit the s key to save a plot.

We will now use the fourth order Runge Kutta method to turn this into a difference equation.

The second order differential equation can be re-written as a first order differential equation on a vector consisting of the phase space variables $\delta$ and $\dot{\delta}$.

$$ \frac{d}{dt} \begin{bmatrix}\delta \\ \dot{\delta}\end{bmatrix} = \begin{bmatrix}\dot{\delta} \\ \omega_P^2\left(i(t) - \sin{\delta}\right) - \frac{1}{RC}\dot{\delta} \end{bmatrix} $$ $$ \vec{y} \equiv \begin{bmatrix}\delta \\ \dot{\delta}\end{bmatrix} $$

$$ \frac{d\vec{y}(t)}{dt} = f(\vec{y}(t),t), y(t_0) = y_0 $$

$$ \vec{y_{n+1}} = \vec{y_n} + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4), $$

$$ t_{n+1} = t_n + h $$

h is the time increment which is always one frame by default so that we are computing positions frame by frame. All time units are in frames in our p5js code. The coefficients are also vectors and are evaluated one after another.

$$ \vec{k_1} = f(t_n,\vec{y_n}) $$ $$ \vec{k_2} = f(t_n + \frac{h}{2},\vec{y_n} + h\frac{\vec{k_1}}{2}) $$ $$ \vec{k_3} = f(t_n + \frac{h}{2},\vec{y_n} + h\frac{\vec{k_2}}{2}) $$ $$ \vec{k_4} = f(t_n + h,\vec{y_n} + h\vec{k_3}) $$

this is how this is simulated.