I'm trying to write a program to plot the I-V charateristics the josephson junction using a python program for some given values of resistance $R$, capacitance $C$, and critical current $I_0$, using the resistively and capacitively shunted junction (RCSJ) model. The josephson equations are \begin{align} I + I_N(t) &= I_0 \sin(\delta) + \frac{U}{R} + C \dot{U} \\ \dot \delta &= \frac{2\pi}{\Phi_0} U \end{align}

where $I$ is the current through the junction, $U$ is the voltage across the junction, and $\Phi_0$ is the flux quantum. The plot given below is given in the material I was referencing (pdf). I was able to write a code and time average it for an interval to solve the RCSJ model, but I'm looking for something to match experimental results which account for thermal rounding (the curvature of the graph changes in experimental results as compared to the model). I'm also looking to encompass the fact that the rate of change of Voltage with Current is not as high as predicted by the model.

enter image description here

As visible in the above picture given in the referenced slides, the curvature entirely changes when we count for the thermal smearing. I'm looking for some help regarding matching the experimental results.

  • $\begingroup$ It seems to me that you have completely changed the question. Initially you were asking about the number of free parameters in the RCSJ model, and now you are asking how to generalize the RCSJ model to include thermal fluctuations. I would suggest asking a separate question for the latter. $\endgroup$
    – J. Murray
    Commented May 20 at 18:55
  • $\begingroup$ @J.Murray well, I wanted a real life simulation of it which definitely needs to include the thermal fluctuations and the RCSJ model I had mentioned contains the In(t) term attributed to the fluctuations, so I've basically just added more detail about what I was able to achieve. $\endgroup$
    – L lawliet
    Commented May 20 at 18:59
  • $\begingroup$ OK I updated my answer to include an elementary (and very physicist) explanation of solving stochastic differential equations, to handle the case of thermal noise. $\endgroup$
    – user34722
    Commented May 22 at 5:33
  • $\begingroup$ Because you are offering a 200 rep bounty, I guess you really want an answer. It would help attract good answers (at least from me) if the question were clear. What is $I_N$? Where is the circuit diagram of the experiment you're trying to simulate? You're also changing the requirements, which discourages people (at least me) from writing an answer. I could answer the original question, but now you're bringing in issues of thermal noise... how do we know you won't change the requirements again later, after we spend an hour writing an answer? Clean up your question and you'll get a good answer. $\endgroup$
    – DanielSank
    Commented May 22 at 6:14
  • $\begingroup$ What you're trying to do here comes down to numerically solving a Lengevin equation. $\endgroup$
    – DanielSank
    Commented May 22 at 6:16

2 Answers 2


First, think about the point $U=0$. Any constant value of $\delta$ is a solution to the second equation. This means the current can take any value between $-I_0$ and $I_0$, as you see at $V=0$ in your plot. This is a manifestation of superconductivity; current can flow through the junction even in the absence of an external voltage. In a sense, the extra variable $\delta$ keeps track of the history of the applied voltage to determine the current at any given time.

This 'extra variable problem' is present in the rest of the plot as well. For non-zero $U$, the current oscillates with time. You could get a single value for the current by time-averaging, but that's not actually what's being plotted here.

However, the plot actually fixes $I$ and then computes the time average of $U$, which is the more interesting case. To do so, they use your second equation to replace $U$ with $\delta$ to get this differential equation for $\delta$.

$$I = I_0 \sin \delta + \frac{\Phi_0}{2\pi R}\dot\delta + \frac{\Phi_0 C}{2\pi}\ddot{\delta}$$

This is in the slides, except I replaced $I+I_N(t)$ with $I$. Fixing $I$ gives a differential equation for $\delta$, which you could solve with something like scipy.integrate.odeint. You can then use your second equation to compute $U$ from $\dot{\delta}$. $U$ is time-dependent, actually oscillatory, so they average over a timescale larger than the oscillation period

$$V \equiv \langle U\rangle = \frac{1}{t_\text{av}}\int_0^{t_\text{av}}U(t)dt$$

Now we have a procedure to compute $V$ from $I$. This is all described on the second slide of section IV, but it's pretty dense. Hopefully, this clarifies things.

Edit based on the question of thermal noise

The slides mention the importance of Johnson Noise, which is the intrinsic noise generated by a resistor at finite temperatures. Johnson noise is white noise, which makes it especially easy to model. White noise means that the power spectrum is flat, and the value of the noise source is independent and uncorrelated. This means you can simulate it by generating a new random value at each time step.

The new answer below is fine, but if you want to simulate it yourself (e.g. for learning's sake), it may be helpful to go back to working with two first-order differential equations, rather than one second-order differential equation.

$$dU = dt (I + I_\text{noise} - I_0 \sin\delta - U/R) / C$$

$$d\delta = dt 2\pi U / \Phi_0$$

I wrote this as differentials instead of derivatives since it's natural to solve these equations by finite difference. Informally, I multiplied both sides of the equation by $dt$. Numerically, $dt$, $dU$ and $d\delta$ will just be numbers, and you'll fix $dt$ to some small value that balances simulation speed with numerical accuracy.

We want to write down a noisy variable that has the right statistics for $I_\text{noise}$ independent of $dt$, since $dt$ is chosen arbitrarily. It's convenient to define the noise increment $dW(t)$ as a Gaussian random variable with mean $0$ and variance $dt$, which you can generate numerically (one random value per time step). To motivate this definition of $dW$, note that it is self-consistent; if we compute the integral of $dW$ over a longer time increment $\delta t$, the variance of $\delta W(t) \equiv \int_t^{t+\delta t} dW(\tau)$ is (on average)

$$E(\delta W^2) = E\left(\int_t^{t+\delta t}\int_t^{t+\delta t} dW(\tau) dW(\tau')\right) = \int_t^{t+\delta t} E(dW(\tau)^2) = \int_t^{t+\delta t} dt = \delta t$$

in line with the variance of $dW$ itself. Formalizing this is the business of Ito calculus, which is worth reading about on Wikipedia or something (especially Ito's lemma, which is useful for deriving and solving differential equations).

Now we can replace $dt I_\text{noise} / I_0$ in the first equation with $\sqrt{4\Gamma} dW$ as per LPZ's answer and solve via finite difference. You will get better results if you use a specialized stochastic differential equation solver, but hopefully, the above shows that there is no magic here.

  • $\begingroup$ I have one question, for a fixed value of I, I can de definitely solve the differential to solve for the phase but what about the initial conditions ? I'm pretty confused about that because how will I obtain U without that. $\endgroup$
    – L lawliet
    Commented May 17 at 19:19
  • $\begingroup$ Also, what would be a good way to incorporate thermal noise ? $\endgroup$
    – L lawliet
    Commented May 17 at 20:57
  • $\begingroup$ I think because you're time averaging, the initial value of $\delta$ won't actually matter. I suggest looking at the solution for different initial conditions and then choosing the integration window so that you capture the steady-state behavior $\endgroup$
    – user34722
    Commented May 18 at 18:03
  • $\begingroup$ Would need more clarification/context about what aspect of thermal noise you want to capture. I would suggest a separate question for that. $\endgroup$
    – user34722
    Commented May 18 at 18:05
  • $\begingroup$ Can you please check out the discussion on chat.stackexchange.com/rooms/info/153224/…. $\endgroup$
    – L lawliet
    Commented May 22 at 7:42

For your simulation, it’s best to start from the dimensionless equation (slide 29): $$ \beta_c\ddot\delta+\dot\delta+\sin\delta=i+i_N \tag{1} $$ In your problem, $i$ is imposed (DC, additional alternating etc.). For your graphs, $i$ is a fixed DC current. What you are actually computing is the effective conductance. Indeed, (1) allows you to solve for $\delta$ from which you deduce $V=\langle U\rangle$ using the Josephson equation (introduced the dimensionless voltage): $$ u:= \frac U{I_0R} =\dot\delta\\ v:=\langle u\rangle $$ The classical analogy for (1) is a damped simple pendulum, with $i$ an external torque (imagine a rocket attached to the pendulum) and $\beta_c$ the moment of inertia. Even without $i_N$ you’ll need to resort to a numerical solver.

With $i_N$, you’ll need a stochastic solver. By the fluctuation dissipation theorem, the fluctuating current is given by the usual Johnson noise formula for the current source spectrum density: $$ S_f=\frac{4kT}R $$ Thus your SDE can be written in the Langevin form (physicists’ notation): $$ \beta_c\ddot\delta+\dot\delta+\sin\delta=i+\sqrt{4\Gamma}\eta \tag{2} $$ with $\Gamma$ the dimensionless temperature: $$ \Gamma = \frac{2\pi kT}{I_0\phi_0} $$ and $\eta$ a gaussian white noise: $$ \langle\eta(t_1)\eta(t_2)\rangle=\delta(t_1-t_2) $$ For mathematicians, this translates to ($W$ being a Wiener process): $$ \beta_cdu=(-u-\sin\delta +i)dt+\sqrt{4\Gamma}dW\\ d\delta=udt $$ You can now plug the SDE in your favorite numerical solver (if I recall correctly, Julia has a better library than Python for SDE solvers). The most elementary method is the Euler-Maruyama method if you want to code it yourself. This gives you $u$ and upon averaging you recover the linked graphs.

Your thermal smearing/rounding can be understood qualitatively. For $1>i>0$ your system normally settles down in the potential minima of: $$ u_J=-\cos\delta-i\delta $$ When $\Gamma>0$, thermal fluctuations allow for the system to hop over the energy barrier. It is the thermal analogue of quantum tunneling (just add in a couple of $i$’s). You can make this quantitative in the weak noise limit $\Gamma\to0$ where you should recover Arrhenius’ law.

Btw, (2) does not hold anymore if you want a $1/f$ noise as in slide 37. You’ll have to change the noise term accordingly.

  • $\begingroup$ Hi, thanks for an amazing answer. Is there some portal we can maybe converse over, as I would like to show you some of the generated plots and their differences from the experimentally obtained data. Then I can award the bounty. $\endgroup$
    – L lawliet
    Commented May 21 at 2:58
  • $\begingroup$ Thanks! Yes there is the chat platform (supports mathjax). I created a room if you want to exchange more: chat.stackexchange.com/rooms/info/153224/… $\endgroup$
    – LPZ
    Commented May 21 at 16:22
  • $\begingroup$ okay I've added more details on it. $\endgroup$
    – L lawliet
    Commented May 21 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.