Confinement into a state with persistent current by thermal
quenching of loop of Josephson junctions
Jorge Berger
Physics Department, Ort Braude College,
P. O. Box 78, 21982 Karmiel, Israel∗
Abstract
We study a loop of Josephson junctions that is quenched through its critical temperature. For
three or more junctions, symmetry breaking states can be achieved without thermal activation, in
spite of the fact that the relaxation time is practically constant when the critical temperature is
approached from above. The probability for these states decreases with quenching time, but the
dependence is not allometric. For large number of junctions, cooling does not have to be fast. For
this case, we evaluate the standard deviation of the induced flux. Our results are consistent with
the available experimental data.
PACS numbers: 74.40.+k, 74.81.Fa, 05.70.Fh, 11.15.Ex
http://arxiv.org/abs/0704.0561v1
We consider a process in which a superconducting loop that contains n identical Joseph-
son junctions is cooled through its critical temperature, in the absence of applied fields, and
monitor the spontaneous generation of metastable states with persistent current. From the
theoretical point of view, this process is enlightening, because it provides and additional
example of a phase transition that is dominated by the time evolution of the system param-
eters, rather than by thermal equilibrium; this is a subject that is still far from being closed,
and is thought to be relevant both to condensed matter physics and to cosmology. From
the practical point of view, this process has significant importance, since flux trapping is a
major obstacle for reproducible functioning of large scale ultra-high-speed superconductiv-
ity digital applications [1] and we would like to comprehend how it depends on the system
parameters.
The best known theory for the description of dynamic phase transitions is the Kibble–
Zurek mechanism (KZM) [2, 3, 4]. It states that in these transitions the equilibrium critical
scalings predict various aspects of the nonequilibrium dynamics of symmetry breaking, in-
cluding the density of residual topological defects. Several numeric simulations have tested
the predictions of the KZM, particularly the dependence of the density of defects on the
quenching time. Some of these simulations lead to refinements of the KZM [5] and others
disagree with it [6]. The KZM has also been tested in several experiments; among them we
will mainly be interested on those performed in superconducting loops [7, 8, 9].
The study of a loop of Josephson junctions is appealing, because the “rules of the game”
are particularly simple [10] and the system parameters can be tailored practically at will.
This system may be neater and qualitatively different from other systems, because it can be
trully divided into n identical subsystems, whereas in other systems the division depends on
a continuously varying coherence length. Moreover, experimental results are already avilable
The supercurrent through Josephson junction i is given by
IJi = Ic(t) sin γi , (1)
where γi is the gauge-invariant phase difference across the junction. Ic vanishes above the
critical temperature and increases when the temperature is lowered. We will consider a
quenching process in which Ic(t) grows from 0 to Ic0. This supercurrent may be interpreted
as arising from a potential energy term −
i Ic cos γi. This potential energy gives rise to
metastable asymmetric states in which the system can be trapped.
There are two typical processes for the formation of topological defects in the KZM. In one
case defects become confined when the order parameter becomes unable to follow the change
of the parameters of the system. Another scenario is activation due to thermal fluctuations
close to the Ginzburg temperature. We shall see that the present system does not fit in
either of these cases; the relaxation time does not diverge at the critical temperature and
no activation energy is required in order to enter a metastable state. More precisely, let
R be the resistance of each junction, L the self-inductance of the loop, let us assume that
the capacitance is small and the resistance of the loop itself (above the critical temperature
Tc) is much smaller than that of the junctions. Then there are two relaxation times in the
problem: one of them is ~/2eRIc, which is infinite regardless of the temperature above Tc,
and the other is L/nR, which remains constant and refers to a process that does not attempt
to align the order parameters of neighboring segments into the same phase. In a sense, our
problem is similar to that of decompression of He4 from the λ-line [11].
Let Ic0 be the maximal superconducting current through the junctions at low temperature
and let us take R, Ic0, ~ and 2e as units. Accordingly, the units of voltage, energy, inductance
and time will be RIc0, 2eRIc0, ~/2eIc0 and ~/2eRIc0.
The state of the loop will be described by the set of values {γi}. The sum of these values
can be interpreted as minus the magnetic flux enclosed by the loop. Since we assume that
there is no applied magnetic flux,
γi = −LI , (2)
where I is the current around the loop. Our goal is to find the probabilities for the possible
values of
i=1 γi after the loop has been cooled.
The rules for the evolution of {γi} are stated in several textbooks [10]. The ac Josephson
relation is
dγi/dt = Vi , (3)
where Vi is the voltage across junction i. The total current through junction i is
I = IJi + Vi + CdVi/dt+ INi , (4)
where IJi is the supercurrent, given by Eq. (1), C is the capacitance of each junction, and
INi is the Johnson current. We will assume that during a period of time τ1 the system is
kept above the critical temperature and Ic = 0, then, during a period of time τ2 the system
is quenched and Ic grows up to Ic0 and finally, during a period of time τ3, Ic = Ic0.
In most of our calculations we will assume that the capacitance is negligible. In this case,
from Eqs. (2), (3) and (4),
= −IJi − INi −
γi . (5)
The case C 6= 0 will be discussed below. We integrate Eq. (5) by Euler iterations. For
this purpose we divide the process into short periods of time of duration ∆t. The Johnson
current averaged over a single period is given by
INi = ηgi/
∆t , (6)
where gi is a random number with zero average and variance 1 and η = (2kBT/R)
1/2, with
kB Boltzmann’s constant and T the temperature. We assume that the temperature remains
close to the critical temperature during the entire process; accordingly, η will be taken as
constant.
We first consider the case in which the loop is cooled instantaneously, i.e., τ2 = 0. If we
ignore the Johnson current, Eq. (5) is equivalent to viscosity-dominated motion of a particle
in n-dimensional space that feels a potential energy (1/2L)(
i=1 γi)
2 − Ic
i=1 cos γi. In
this situation, {γi} evolves towards a local minimum of the potential energy. During the
first stage, Ic(t) = 0 and the only local minimum is the plane
i=1 γi = 0. During the last
stage, Ic(t) = 1. At this stage the absolute minimum is located at the origin, γi = 0, but
for sufficiently large values of L additional local minima may also exist. We ask whether
the values {γi} could wander in the plane
i=1 γi = 0 and then, when the temperature is
lowered, flow into a local minimum different from γi = 0. For our present purpose, two
states such that their respective values of γi differ by integer multiples of 2π and
i=1 γi is
the same for both, will be considered equivalent.
Figure 1 shows that this situation is possible for the case n = 3. The figure shows
five evolution curves that start at the plane
i=1 γi = 0 and, in the absence of thermal
fluctuations, flow to a local minimum of the potential energy. Note that although the
temperature is assumed to change instantaneously,
i=1 γi builds up during a lapse of time
of the order of ~/eRIc0.
FIG. 1: Several evolution lines for the case of three junctions, for L = 4 and in the absence of
thermal fluctuations. All these curves start at the plane γ1+γ2+γ3 = 0 and converge to the point
γ1 = γ2 = −0.98, γ3 = 2π − 0.98.
We studied the evolution including thermal fluctuations, for n = 3 and n = 4, for several
values of η. In all cases, we started from values of {γi} randomly located in the interval
−π < γi < π. In order to achieve initially an equilibrium distribution, evolution was followed
during a period τ1 with Ic = 0. After that, evolution was followed during a period τ3 with
Ic = Ic0. In order to decide what is the “final” state, we should average over an additional
period of time, in order to filter out thermal fluctuations. We found it easier to turn off at
this stage the fluctuations and let the state converge to the nearest local minimum. For each
set of values, this process was reapeated 1000 times and the probability for confinement in a
given state was evaluated as the number of times this state was obtained, divided by 1000.
In most cases the final state was the ground state γi = 0, but the first excited state was also
reached. Due to the symmetry of the problem, these states are degenerate, i.e., the γi’s can
be permuted and all the signs can be inverted. For the parameters we considered, we did
not find cases with higher excited states.
The results are shown in Fig. 2. We avoided values of η that might be too small to enable
thermalization during the period τ1. The probabilities shown in the graph correspond to the
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
FIG. 2: Probability for a current-carrying metastable state as a function of the size of thermal
fluctuations. Parameters used: τ1 = 40000, τ2 = 0, τ3 = 10000, ∆t = 0.1, C = 0. ♦ n = 3, L = 4;
⋆ n = 4, L = 4; � n = 4, L = 2.
total probability of reaching any of the (degenerate) excited states. As a general trend, we see
that the probability of ending at an excited state increases with the number of junctions and
with the normalized self-inductance. We also see that this probability is fairly independent
of the size of thermal fluctuations, until a sufficiently large value of η is reached. Beyond
this value, there is a fast decrease of this probability.
It is reasonable to expect that the probability for the metastable state will decrease
when the thermal energy kBT becomes comparable to the energy barrier that confines this
state, i.e., the difference between the energy at the saddle-point and the energy at the local
minimum. For n = 3, L = 4, a local minimum is at γ1 = γ2 = −0.98, γ3 = 2π − 0.98, the
corresponding saddle point is at γ1 = γ2 = −0.67, γ3 = π + 0.67 and the energy difference
is 0.25; for n = 4, L = 4, a local minimum is at γ1 = γ2 = γ3 = −0.83, γ4 = 2π − 0.83,
the saddle point is at γ1 = γ2 = γ3 = −0.54, γ4 = π + 0.54 and the energy difference is
0.42; similarly, for n = 4, L = 2, the energy barrier is 0.13. In all cases we find that the
probability for the metastable state decreases to about half its maximum value when the
thermal energy is about one eighth of the barrier energy. Clearly, the precise value depends
on τ3; in principle, for τ3 → ∞, the metastable state should always decay. If the thermal
energy becomes of the order of the energy difference between the excited and the ground
state, then the probability for the excited state will increase with temperature (equilibrium
probability), but we are not interested in this regime.
Let us now study the influence of the cooling time on the probability for a metastable
0 10 20 30 40 50 60
FIG. 3: Probability for a metastable state as a function of the cooling time for loops with 3 or 4
junctions. The empty (filled) symbols correspond to Ic(t) proportional to t (to t
2) and are fitted
by dashed (continuous) lines. For visibility, the line for n = 3, L = 4 and Ic(t) ∝ t2 has been
lowered by 0.025. Parameters used: τ1 = 40000, τ2 + τ3 = 10000, ∆t = 0.1, C = 0. Unless stated
otherwise, η = 0.1. ♦ n = 3, L = 4; � n = 4, L = 2; △ n = 4, L = 4, η = 0.2.
0 500 1000 1500 2000
0 5000
FIG. 4: Like Fig. 3, for n = 5 and n = 6. The inset shows P (τ2) in the range 0 ≤ τ2 ≤ 104. For
visibility, the line for n = 5, L = 2 and Ic(t) ∝ t has been raised by 0.025.♦ n = 5, L = 2; △ n = 5,
L = 1; ◦ n = 6, L = 2.
state. At the moment that Ic(t) becomes different from zero, there will be an incentive for
leaving this state; on the other hand, as long as Ic(t) is small, the confining barrier will also
be small and the way out will be easy. We therefore expect that the trapping probability
will decrease with τ2. In most of our calculations we assumed that Ic(t) is proportional to
the temperature below Tc and therefore increases linearly with t, but we also considered
the case Ic(t) ∝ t2, which is more realistic for strong coupling [9]. For the present purpose,
simulations were repeated 104 times. Figures 3 and 4 show our results for several values of
n, L and η.
The topological charge for a given state may be defined as the sum of the topological
charges of all junctions, where the topological charge of junction i is the closest integer to
γi/2π. For all the cases considered in Figs. 3–4, the topological charge was 0 or ±1.
The probabilities shown in Figs. 3–4 are also the expectations of the absolute value of
the topological charge. It is therefore tempting to identify this probability with the density
of defects, and anticipate that it will decrease as a power of τ2. However, the arguments
that lead to the time dependence of the defect density in [4] seem to be irrelevant in the
present case; there is no obvious way to associate the presence of topological charge to some
primordial coherence length and, indeed, our results cannot be fitted by a power dependence.
Denoting the probability by P , most of the curves in our results (typically for small n and
L) can be fitted by the empiric form P (τ2) ∝ exp[−(τ2/τ0)
2/n] in the case Ic(t) ∝ t and
by P (τ2) ∝ exp[−(τ2/τ0)
2/(n−1)] in the case Ic(t) ∝ t2. The characteristic time τ0 depends
very strongly on the number of junctions and only weakly on the size of the energy barrier
or on the temperature. For n = 3, 16 . τ0 . 17; for n = 4, 26 . τ0 . 39; for n = 5,
350 . τ0 . 540.
Our empirical fits suggest that for n ≫ 1 the probabilities for metastable states decrease
very slowly with the cooling time. Indeed, in the experiment that involved 214 junctions
[7], the distribution of permanent currents was found to be independent of the cooling time
(up to the order of a minute).
Part of the probabilities shown in Figs. 3–4 do not decrease at a uniform rate. Instead,
they seem to decay in two stages. A possible explanation might be that the region in phase
space that in the absence of thermal fluctuations would flow into a metastable state can be
divided into two subregions, such that escape from one subregion is much easier than escape
from the other.
Let us now consider large values of n and L, as were encountered in the experiment.
In this case many different final metastable states are possible, and the most significant
experimental quantity will be the variance of the induced flux. Our results are shown in
Fig. 5. As in the case of small values of n and L, the general trend is increase of the typical
flux with increase of either n or L. However, these individual increases appear to saturate.
For instance, for L . 200, the standard deviation of the flux actually decreases with n in
50 100 200 500
2eLIc0�Ñ
50 100 200
FIG. 5: Standard deviation of the induced flux, Φ/Φ0 = −
i γi, as a function of the self-inductance
and the number of junctions. Two curves are for fixed n, and L is shown in the lower abscissa; the
other two are for fixed L, and n is shown in the upper abscissa. The abscissas are in logarithmic
scale. The symbols have been joined for visibility. The parameters are as in Fig. 2 and η = 0.1.
Each simulation was repeated 400 times. � n = 100; △ n = 214; ◦ L = 100; ⋆ L = 600.
the range 100 . n . 200.
In the experiment [7], n = 214 and 〈Φ2〉1/2/Φ0 = 7.4 ± 0.7, where Φ is the induced flux
and Φ0 the quantum of flux. Comparison with Fig. 5 indicates that 2eIc0L/~ should be in
the range between ∼ 100 and ∼ 150. The estimates of Ref. [7] suggest that 2eIc0L/~ ∼ 600.
Since the experimental estimate was not based on a measurement, but rather on a plausibility
argument for the size of Ic0, and the junctions are not really all identical, the agreement is
reasonable.
Let us finally consider the case C 6= 0. In this case we integrated Eqs. (3) and (4) as a
system of differential equations. In digital applications, a preferred value is C = ~/2eIc0R
which provides for fast switching without oscillations. We have repeated our calculations
for this case and for several representative values of the other parameters. We found that a
capacity of this size has no qualitative effect.
In summary, we have performed simulations that describe the formation of symmetry
breaking states when a loop of Josephson junctions is quenched. Among the typical systems
in which symmetry breaking occurs in a dynamics-dominated process, the present system
constitutes a class of its own. Our results agree with the experiment in the case of large n
and raise predictions for the case of small n.
This work has been supported by the Israel Science Foundation under grant 4/03-11.7. I
am grateful to Alan Kadin for useful comments.
∗ Electronic address: phr76jb@tx.technion.ac.il
[1] M. Jeffery, T. Van Duzer, J. R. Kirtley, and M. B. Ketchen, Appl. Phys. Lett. 67, 1769 (1995).
[2] T.W. B. Kibble, J. Phys. A 9, 1387 (1976); Phys. Rep. 67, 183 (1980).
[3] W. H. Zurek, Nature (London) 317, 505 (1985); Acta Phys. Pol. B 24, 1301 (1993).
[4] W. H. Zurek, Phys. Rep. 276, 177 (1996).
[5] P. Laguna and W. H. Zurek, Phys. Rev. Lett. 78, 2519 (1997); A. Yates and W. H. Zurek,
Phys. Rev. Lett. 80, 5477 (1998).
[6] M. Hindmarsh and A. Rajantie, Phys. Rev. Lett. 85, 4660 (2000).
[7] R. Carmi, E. Polturak, and G. Koren, Phys. Rev. Lett. 84, 4966 (2000).
[8] J. R. Kirtley, C.C. Tsuei, and F. Tafuri, Phys. Rev. Lett. 90, 257001 (2003).
[9] R. Monaco, J. Mygind, M. Aaroe, R. J. Rivers, and V. P. Koshelets, Phys. Rev. Lett. 96,
180604 (2006).
[10] A. M. Kadin, Introduction to Superconducting Circuits (Wiley, New York, 1999); M. Tinkham,
Introduction to Superconductivity (McGraw-Hill, New York, 1996); K. K. Likharev, Dynamics
of Josephson Junctions and Circuits (Gordon and Breach, New York, 1986); A. Barone and
G. Paternò, Physics and Applications of the Josephson Effect (Wiley, New York, 1982).
[11] P.C. Hendry, N.S. Lawson, R.A.M. Lee, P.V.E. McClintock, and C.H.D. Williams, in: Forma-
tion and Interactions of Topological Defects, ed. A.C. Davis and R.N. Brandenberger (Plenum,
New York,1995).
mailto:phr76jb@tx.technion.ac.il
	References

  We study a loop of Josephson junctions that is quenched through its critical
temperature. For three or more junctions, symmetry breaking states can be
achieved without thermal activation, in spite of the fact that the relaxation
time is practically constant when the critical temperature is approached from
above. The probability for these states decreases with quenching time, but the
dependence is not allometric. For large number of junctions, cooling does not
have to be fast. For this case, we evaluate the standard deviation of the
induced flux. Our results are consistent with the available experimental data.

Introduction to Superconducting Circuits (Wiley, New York, 1999); M. Tinkham,
Introduction to Superconductivity (McGraw-Hill, New York, 1996); K. K. Likharev, Dynamics
of Josephson Junctions and Circuits (Gordon and Breach, New York, 1986); A. Barone and
G. Paternò, Physics and Applications of the Josephson Effect (Wiley, New York, 1982).
[11] P.C. Hendry, N.S. Lawson, R.A.M. Lee, P.V.E. McClintock, and C.H.D. Williams, in: Forma-
tion and Interactions of Topological Defects, ed. A.C. Davis and R.N. Brandenberger (Plenum,
New York,1995).
mailto:phr76jb@tx.technion.ac.il
	References
