Solving the hysteresis loop calculation problem for Josephson junction stacks Текст научной статьи по специальности «Медицинские технологии»

UDC 519.6

Solving the Hysteresis Loop Calculation Problem for Josephson Junction Stacks

S. I. Serdyukova

Laboratory of Information Technologies Joint Institute for Nuclear Research 6, Joliot-Curie str., Dubna, Moscow region, Russia, 141980

A detailed investigation of the IVC breakpoint and the breakpoint region width gives important information concerning the peculiarities of stacks with a finite number of intrinsic Josephson junctions. The current-voltage characteristics for a stack of n Josephson junctions is defined from solving the system of n nonlinear differential equations. The current voltage characteristic has the shape of a hysteresis loop. On the back branch of the Hysteresis loop, near the breakpoint lb, voltage V(!) decreases to zero rapidly. The goal of this work is to accelerate the computation of IVC based on numerical solution of the system. A numerical-analytical method was proposed in. This method showed perfect results in IVC calculations for a stack of 9 and 19 intrinsic Josephson junctions and the computation time reduced by five times approximately. The question of choosing a change-over point from "analytical" to numerical calculation was open. In testing computations the change-over point was taken equal to 21b. In the case of periodic boundary conditions an equation, determining the approximate location of lb, was obtained. This moment we succeeded to develop an algorithm determining the approximate value lb in more complicated technically case of non-periodic boundary conditions with 7 =1. All calculations were performed using the REDUCE 3.8 system.

Key words and phrases: stack of Josephson junctions, computation of current-voltage characteristics, hysteresis loop, Cauchy problem for a system of nonlinear differential equations, fourth-order Runge-Kutta method, long-time asymptotic formulas, a numerical-analytical method, computation of formulas using the REDUCE 3.8 system.

1. Introduction

Solving the system

n

Pi = ^ Ai,i' (I — sin(<^') - fjp I'), I = 1,...,n, (1)

i'=1

for different I : I = I0 + kAI < /max; I = Imax — kAI, the current-voltage characteristics of stacks as hysteresis loops are found [1]. For initial value of the current (I = I0) the system (1) is solved with zero initial data on an interval [0, Tmax]. For each next I : I = Ik+1, found already pi(Ik,Tmax),pi(Ik,Tmax) are used as initial data. On the back branch of the Hysteresis loop, near the breakpoint lb, voltage V(I) decreases to zero rapidly. The goal of this search is to accelerate the computation of IVC.

In the case of periodic boundary conditions the A matrix is

/1 + 2« —a 0

0

\ —a

—a 0 1 + 2« —a —a 1 + 2«

0

—a

—a \ 0

—a 1 + 2« —a 0 —a 1 + 2«/

(2)

0

0

Received 27th September, 2013.

/1 + a(1 + 7) —a 0 0 0

—a 1 + 2 a — a 0 0

0 — a 1 + 2 a — a 0

0 0 — a 1 + 2 a —a

V 0 0 0 — a 1 + a(1 + 7)/

square matrix of order n. And in the case of nonperiodic boundary conditions the A matrix is

(3)

where j = s/s0 = s/sn and s, s0, sn are sickness of middle, first and last superconducting layers respectively [1]. The parameter a gives the coupling between junctions, / is the dissipation parameter. The dynamics of phase differences pi(t) had been simulated by solving the equation system (2) using the fourth order Runge-Kutta method [2]. After simulation of the phase differences dynamics the voltages on each junction were calculated as

n

dpi/dt = ^a^VY . (4)

i'=i

The average of the voltage V is given by

T

± max

V* = =-f Vidt. (5)

— max — min J T ■

J- min

Finally the total voltage V of the stack is obtained by summing these averages:

n

V = (6)

i=i

The calculation can be simplified using specific properties of the matrices (2), (3). These matrices are symmetric. They have complete systems of orthonormal eigenvectors Ei with real eigenvalues Xi . The fundamental matrices D (whose columns are Ei) reduce the A-matrices to the diagonal form: D'AD = A = diag (X1, X2,..., Xn).

After changing the variables

n n

pi = , V=

I'=1 I'=1

we get a system:

n

Ipi = -Xi/3ipi + Xi • I • Si - Xi di',i sin(pr),

v=i

where I = 1,... ,n and Si is the sum of Ei elements: Si = d1ti + d2,i + ... + dn,i. Relations (4), (5), (6) result in

d^l _ w. _ ^¡(Tmax) - ^(Tmin)

respectively.

T = Wn Wi = y'((Jm;x) — WW, v = y (7)

dt 1 11 1 h(Tmax -Tmin) ' fe ()

2. Periodic Boundary Conditions

In the case of periodic boundary conditions, the eigenvalue problem of A has the solution:

Xi = 1 + 2a(1 - cos(^)), Pi =

2n(l - 1)

n

I = 1,..., ns,

where ns is integer part of (n + 1)/2, ns = [(n + 1)/2], "H r sin(^)

E^-V = \l n

sin(2<^ )

.sin(npi).

E21-1 = A/ —

'cos(pi) ■ cos(2^)

.cos(npi).

I = 2,..., ns.

When n is even, n = 2 ■ ns, matrix A has additional eigenvalue Ans+1 = 1 + 4a, the corresponding pns+1 = n and corresponding eigenvector is

En = 1,1, —1,1,..., —1,1]*.

So every of Xi, 2 < I < ns, has a pair of eigenvectors. Let me remind that Si is a sum of Ei elements. This time we have S1 = y/n, Si = 0,1 = 2,... ,n.

As a result, the hysteresis calculation problem is reduced [3] to solving the unique equation

V(t) = £1 + ^^(1 - e-^) + ut - 11(1 - e-^(i-s)) sin(V(s))ds.

0

(8)

Solving this equation we find ^1(t) = y/n^(t). The rest components (t), j = 2,... ,n are equal zeros. The equation (8) is solved by the simple iterations method. Starting from rq0 = 0. we obtained at the third iteration step

Vs(t)

u —

cos($)

+

sin($)

t + A -

cos(^ + ■&)

2u(u2 + p2) 2fi(u2 + (32))" ' " f3w cos(2^ + &) sin(wt + A + & + atg)

+

u

u2

+ 0(œ-4 + e-^). (9)

Here w = 1/(3, A = ^ + (£2 — , atg = arctg(^/w), § = — cos(A)/(u/3). Remark that V(I,n) = -nW1(I) (see (7)) and

W1(I) = -E(v(I,Tmax) — V(I,Tmin))/(Tmax — ^min).

The approximate break point location w can be found [4] from (9) as a solution of the equation F(w) = 0, where

t-,/ \ sin(^ - atg)

F (w) = u +-v , 67 +

2Pu^/u2 + /32

T

2 sin(w(Tmax - Tmin)/2) cos(w(Tmax + Tm[n)/2 + A + ê + atg)

u

VP2 + U2 (Tmax - Tmin)

1

1

t

Here A = -u/f, according to = £2 = 0. The polynomial P(x) = 4f2u4(cu2 + f2) -1 has the unique positive root xt = 1.35232. We find that F(xt) = 1.447... and F(1) = -1.434.... After this no difficulty is to calculate the approximate break point location using the interval bisection method, Ib = 0.210248 .... Roughly spiking the jump to numerical calculations must be done at 2Ib. In our calculations we put Tmin = 50, Tmax = 1000, AI = 0.05. The step in the Runge-Kutta method was h = 0.1. All calculations were performed by using the REDUCE 3.8 system [5] .

3. Nonperiodic Boundary Conditions

In the case of nonperiodic with 7 =1 boundary conditions, the A matrix has following eigenvalues and eigenvectors: Xj = 1 + 2a(1 — cos(jtf)), $ = n/(n +1), j = 1,... ,n;

Ej = cn[sin(j§), sin(2j$),..., sin(nj$)] , j = 1,... ,n, cn = V2/( n + 1). For even j, j = 2,4,..., 2 k < n, Sj = 0. And for odd j

Sj = cn ctg(j$/2), j = 1, 3,..., 2k - 1 <n.

We proved [3] that in the case of non-periodic boundary conditions with 7 = 1 the problem of Hysteresis loop calculation reduces to solving the following system of n integral equations:

^21-1 = U2i-ii + Hi(2l - 1) + 6(2 -,1) -U2-1 (1 - exp(-fX2i-it))

fX21— 1

1 t n / ns \

-1 (1 - exp (-fX2i-i(t - s)))^2dm,21-1 ds, (10)

f 0 m=1 Vk=1 J

where 1 = 1, 2,..., ns, ns is integer part of (n + 1)/2 and u2i-1 = S2-11/f.

For each I and given initial data ^(2 I - 1), £2(2 I - 1) the system (10) was solved using simple iterations starting at zero. The results obtained after three iterations are regarded as "asymptotics" of the solution for large t. These "asymptotics" were used in [3].

In result of a number of useless attempts to find approximate breakpoint location by analogy with [4], we concluded that this can be done as follows. It is sufficient to calculate V(I) for different 1,1 = 0.5 - 0.05 ■ j, solving the system (10) with zero initial data "analytically", until /0, satisfying V(/0) ■ V(I0 + 0.05) < 0 was obtained. Found I0 is taken for the approximate breakpoint location. In the case of the non-periodic boundary conditions we succeeded to realize this algorithm only after refusing a number of algebraic manipulations "eating" the time and leading to computing errors accumulation as well. The calculating time decreased more than 9 times (against 5 declared in [3]). Below we present result of calculating It, = 0.3 for the stack 0f 19 Josephson junctions.

vb :={{0.5,33.7167375273},{0.45,40.235725455},{0.4,63.5053115073}, {0.35,62.9218862448}, {0.3, - 0.0454510614144}}

In Fig. 1 the pictures of the back way of the hysteresis loop are shown. The solid and dotted lines refer to numerical and "analytical" calculations respectively. In Fig. 2 the solid line is the same as in Fig. 1, while the circles on this line refer to calculation performed by the following mixed numerical-analytical method. The right way of the hysteresis loop and the back way on the interval 1.45 > I > 0.45 = 1.5 ■ lb are computed using the "asymptotic" formulas. The rest points of the hysteresis loop are computed numerically. The calculations were performed for a = 0.2, f = 0.2 using

the REDUCE 3.8 system [5] with Tmin = 50,Tmax = 1000 and the step h = 0.1 was chosen in the numerical calculations.

Figure 1. The solid line refers to the back branch of the hysteresis loop for n = 19, calculated numerically using the fourth order Runge-Kutta method. The dotted line refers to the back branch of the hysteresis loop, calculated "analytically" using the "asymptotic" formulas

Figure 2. The solid line refers to the back branch of the hysteresis loop calculated numerically for n =19. The circles on this line refer to calculation performed by the mixed analytical- numerical method: the whole right branch of the hysteresis loop, together with the back branch at 1.45 > I > 0.45 = 1.5 • lb,

have been computed using the "asymptotic" formulas. The points at 0.45 ^ I ^ 0.2 were computed numerically. The point (0.5, 43.561 ...), marked in Fig. 2 by the arrow, is the last point of the hysteresis loop calculated

"analytically"

References

1. Shukrinov Y. M., Mahfousi F., Pedersen N. F. Investigation of the Breakpoint Region in Stacks with a Finite Number of Intrinsic Josephson Junctions // Phys. Rev. B. — 2007. — Vol. 75. — P. 104508.

2. Бахвалов Н. С., Жидков Н. П., М. К. Г. Численные методы. — М.: Мир, 1977. [Bakhvalov N. S, Zhidkov N.P., Kobelkov G.M. Numerical Methods. — Moscow: Mir, 1977. — (in russian). ]

3. Serdyukova S. I. Numerical-Analytical Method for Computing the Current-Voltage Characteristics for a Stack of Josephson Junctions // Computational Mathematics and Mathematical Physics. — 2012. — Vol. 52, No 11. — Pp. 1590-1596.

4. Serdyukova S. I. Determination of IVC Breakpoint for Josephson Junction Stack. Periodic and Non-Periodic with 7 = 0 Boundary Conditions // Particles and Nuclei, Letters. — 2013. — No 3. — Pp. 269-272.

5. Neun W. REDUCE User's Guide for Unix Systems. Version 3.8. — Berlin: ZIB, 2004.

УДК 519.6

Решение проблемы вычисления петли гистерезиса для систем джозефсоновских переходов

С. И. Сердюкова

Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, Дубна, Московская область, Россия, 141980

Детальное исследование критических точек ВАХ и оценка области их влияния представляют большой интерес для изучения свойств систем с конечным числом внутренних джозефсоновских переходов. Вольт-амперная характеристика для системы п внутренних джозефсоновских переходов определялась по решению системы п нелинейных дифференциальных уравнений. Вольт-амперная характеристика (сокращённо ВАХ) имеет вид петли гистерезиса. На обратной ветви петли гистерезиса, при подходе к точке излома 1ь, напряжение V(I) резко спадает к нулю. Цель этой работы — ускорить процесс вычисления ВАХ, основанный на численном решении системы. Был предложен смешанный численно-аналитический алгоритм. Этот метод показал прекрасные результаты при вычислении ВАХ для систем 9 и 19 внутренних джозефсоновских переходов. При этом время счета по смешанному методу сократилось приблизительно в пять раз. Оставался открытым вопрос выбора точки перехода от «аналитического» счета к численному. При тестовых расчётах точка перехода принималась равной 2Д. В случае периодических граничных условий было получено уравнение, определяющее приближенное значение 1ь. В настоящий момент удалось разработать алгоритм, определяющий приближенное значение 1ь в более сложном технически случае непериодических граничных условий. Все вычисления производились с использованием системы REDUCE 3.8.

Ключевые слова: система джозефсоновских переходов, вычисление вольт-амперных характеристик, петля гистерезиса, задача Коши для систем нелинейных дифференциальных уравнений, метод Рунге—Кутта четвёртого порядка точности, асимптотические формулы решения задачи Коши при больших t, численно-аналитический метод вычисления ВАХ, вывод асимптотических формул, используя систему REDUCE 3.8.