CC BY
31
UDC 519.642.2
Numerical Method for Computation of Sliding Velocities for Vortices in Nonlocal Josephson Electrodynamics
E. V. Medvedeva
Department of Higher Mathematics-1 National Research University of Electronic Technology 5, Passage 4806, Zelenograd, Moscow, Russia, 124498
In this paper, a model of infinite Josephson layered structure is considered. The structure consists of alternating superconducting and tunnel layers and it is assumed that (i) the electrodynamics of the structure is nonlocal and (ii) the current-phase relation is presented by sum of Fourier harmonics instead of one sinusoidal harmonic for the case of the sine-Gordon equation. The governing equation is a nonlocal generalization of the nonlinear Klein-Gordon equation with periodic nonlinearity that depends on external parameter of nonlocality A. The velocity of vortices (2^-kinks) in models of such kind are not arbitrary, but belong to some discrete set. The paper presents a method for computation of these velocities (called also "sliding velocities") and the shapes of kinks. The estimation of error of the method is given. The results of computations are the families of 2^-kinks parametrized by A. It is observed that the 2^-kinks corresponding to different families for the same A have nearly the same central part but differ in asymptotics of the tails. The numerical algorithm has been incorporated into a program complex "Kink solutions" in MatLab environment. The complex enables to compute the shapes and velocities of 2^-kinks for nonlinearities represented by sums of up to ten Fourier harmonics, as well as to model the propagation of these kinks.
Key words and phrases: Josephson junction, nonlocal Josephson electrodynamics, embedded solitons, sliding velocities, nonsinusoidal nonlinearity.
1. Introduction
In 90-es, in series of papers (see the survey [1]) it has been shown that in some situations the electrodynamics of distributed Josephson junction becomes nonlocal. The basic equation to describe the vortex dynamics is no longer the sine-Gordon one, and it should be replaced by its nonlocal generalization (nonlocal sine-Gordon equation, NSGE). It was found [2] that contrary to the sine-Gordon equation NSGE, does not support travelling simplest Josephson vortices (called also 2k-kinks or fluxons). At the same time it describes fast 4^-, 6^- etc kinks and each of them can travel only with its own velocity that depends on its shape.
Later, nonlocal Josephson electrodynamics of layered structures also was developed [3,4]. It was shown that the equation for vortex dynamics in this case is also NSGE. In the recent paper [5] an extension of the approach of [3] to the case of non-sinusoidal current-phase relation (CPR) was done. Complex CPRs are relevant to many situations, in particular to SIS or SNINS and SFIFS junctions [6]. In [5] it was assumed that CPR is described by two Fourier harmonics, therefore the governing equation was the nonlocal double sine-Gordon equation (NDSGE). It was found that NDSGE does support 2^-kinks, however the velocity of 2^-kink cannot be arbitrary, contrary to the local case. There exists a discrete set of so-called sliding velocities and each of them is associated with the shape of the corresponding vortex.
The study of sliding velocities for NSGE and related equations is quite complex problem even from the numerical viewpoint. The difficulty is that the velocity and the shape of the vortex should be found simultaneously. Therefore the methods that are applicable in the local case (see e.g. [7]) hardly can be applied to the nonlocal problems. So, special technique should be developed, see [8].
Author is grateful to Prof G.L. Alfimov for reading of the manuscript and for many useful comments and to Prof. A.S. Malishevskii for useful instructive discussions.
In this paper I describe a numerical algorithm which allows to find numerically the sliding velocities for nonlocal generalizations of the sine-Gordon equations. This algorithm is a part of computer application "Kink solutions". This computer application allows to find the set of sliding velocities for the case of arbitrary CPR and then to study the vortex propagation. The algorithm is based on statements from the dynamical system theory (see [8]). Some numerical results and the estimation of numerical error are presented.
2. The Model
Consider a Josephson structure consisting of alternating superconducting and tunnel layers (see Fig. 1) assuming that (a) superconducting layers S are identical; (b) tunnel layers I are identical; (c) vortex formation is symmetric, i.e. the phase difference and the magnetic field in all the JJs are identical; (d) the superconducting layers S are thin when compared with the penetration depth of magnetic field into superconductors; (e) the structure is assumed to include infinitely many layers, in order to disregard boundary effects on the edges.
Figure 1. Layered structure: S — superconducting layers, I — tunnel layers. A vortex formation moves from left to right
Assume that the CPR is of the form
J(<p) = Jc sin + J2 sin + J3 sin + ..., (1)
where Jc is the critical Josephson current. The straightforward generalization of the model of [5] to infinitely many harmonics in CPR yields
F (?) + ufvtt = A ^ J dx' e_|x_x' l/Ae« (x', t). (2)
Here = <p(x, t) is the phase difference of the order parameters across each tunnel layer. The parameters are: uj is the Josephson plasma frequency and
A = (A- + d)X* , Aeff . XL^V
2\lV(L + № * L + d
where Xj is the Josephson length, XL is the London penetration depth, 2d is the thickness of the each tunnel layer, 2L is the thickness of the each superconducting layer. Nonlinearity in Eq. (2) comes from (1),
F(<p) = sin + A2 sin 2^ + A3 sin 3^ + ...,, (3)
where A2 = J2/Jc, A3 = J3/Jc, .... In dimensionless variables,
l~L + d x
c "V x^Td X], T = ujt>
Eq. (2) takes the form
F (?) +<Prr = ~ y^'e-K-O/V (C, r), (4)
where ^ = (p(Ç, t) and the parameter of nonlocality A is
Aeff I L + d VLAl
A
A
Aj V Al + d \j\J\l + d' For the traveling wave solutions, p(z) = p(£ — vr), Eq. (4) reads
F(p) + v2pzz = —/ e~lZ~zVX^'(*') dz>. (5)
A single Josephson vortex of one flux quantum corresponds to the 2^-kink solution of Eq. (5). It obeys the boundary conditions
lim ip(z) = 0, lim ip(z) = 2ir. (6)
Eq. (5) with the boundary conditions (6) is the main object of the study in this paper.
3. Numerical Method 3.1. The Idea of the Method and the Numerical Algorithm
In order to study 2-^-kink solutions of Eq. (5) the approach of [8] was adopted. Following this approach Eq. (5) should be rewritten in the form of the system
v2pZz = qz + F (p), —A2 qzz + q = pz. (7)
The first integral of (7) is
1 )2 A 2 1
1 = yp2 + At — 112 + u (p), (8)
where U(p) is a primitive of F(p). The system (7) has Hamiltonian structure [2]. Reduction of the nonlocal equation to a system of ODE drastically simplifies both numerical and analytical studies of the problem.
Consider 4D phase space for the system (7) with coordinates p,pz,q, qz. Then 2-^-kink solutions of the system (7) correspond to heteroclinic orbits which connect equilibrium states Oo(0,0,0, 0) and 02li(2tt, 0,0,0). Singular points O0 and 02^ are equilibrium states of saddle-center type. This means that for O0 there exist a pair of trajectories 7+2(O0) which approach O0 as z ^ —<x (outgoing trajectories) and one pair of trajectories 7-2(O0) which approach O0 as z ^ (incoming trajectories). The same situation takes place for O2li, the corresponding trajectories are 7+2(O2^) (incoming) and 7-2(O2^) (outgoing). The existence of the heteroclinic orbit means merging of the trajectories: one of the outgoing from O0 and one of incoming in O2-K. In the 4D phase space of the Hamiltonian system (7) this merging is not a general case, but this may occur for some values v, when A is fixed. In this case heteroclinic orbit has to pass through the plane of symmetry (p = k, qz = 0) (this follows from the statement 4 in [8]).
Numerical algorithm for finding the values of sliding velocity v consists in "adjusting" the parameter v in such a way that one of the outgoing trajectories from O0
(7+(Oo), for definiteness) passes through the plane (p = K,qz = 0). Technically, it can be implemented by finding zeroes of the function R(v,X), which for given parameters A and v is defined as
R(v,X) = qz(zo; X,v),
where (p(z; X, v), >pz(z; X, v), q(z; X, v), qz(z; X, v)) is an outgoing trajectory for O0 and z0 satisfies the condition p(z0; X,v) = The zeroes can be found by means of dichotomy method.
Summarizing, the method can be algorithmized as follows.
1. Set initial conditions S = (^(0; X,v),pz (0; X,v),q(0; X,v),qz (0; X,v)) using the linearization of (7) in the vicinity of the point O0 (i.e. to find the point S with a controllable accuracy on the trajectory 7+(00) near the equilibrium point 00).
2. Solve numerically the Cauchy problem with initial conditions at the point S until the moment when the component p on the trajectory takes the value ■k. Fix R(v, X) = qz(z0; X, v) at this point.
3. Changing the parameter v find the values of v at which R(v,X) = 0. The component p on the corresponding trajectory yields the 2^-kink solution.
3.2. Error Estimation
There are two main sources of numerical error when calculating the values v with this algorithm. They are: (a) the error of the numerical solution of the Cauchy problem; (b) the error of choosing initial data S. Let us describe methods for control the accuracy for the points (a) and (b) consecutively.
a. Error of the numerical solution of the Cauchy problem. The Cauchy problem was solved numerically by 4-th order Runge-Kutta method. The error estimation in this case can be done in a standard way using the Richardson rule, see [9]. Consider a sequence of uniform grids with steps h, h/2, h/4, etc. The nodes of every grid coincide with even nodes of the next grid. Let (z) be a solution for some of these grids and let <p2(z) be a solution for the next, half-step grid. The estimation of the error at the point z is defined as A(z) = [p2(z) — (z)]/(2p — 1), where p = 4 is theoretical order of accuracy for a numerical method. Consider the norm of the error
||A||C = max |A(zn)|, (9)
where N is the total number of points in the grid for which the norm was calculated. Fig. 2, panel A, shows the values of the norm for the error of the 4-th order Runge-Kutta method versus a number of nodes N ~ 1/h. Since the asymptotics of the error obeys power law, double logarithmic scale was used: log10 N along the axis of abscissa and log10 ||A|| along the axis of ordinate. As N increases the value log10 ||A|| follow a straight line with a slope tan a = —4, however, this takes place until the roundoff errors become comparable with the error of the method. According to the calculations, asymptotical estimation of error, ~ Ch4, is valid for values of h £ [0.0004, 0.04].
b. Error of initial condition. Initial data
5 = (p(0; A, v), pz(0; A, v), q(0; X, v),qz(0; A, v)) was found by linearization of the system (7) at the point O0. The linearized system is
v2pzz = qz + —X2qzz + q = <pz,
where Q = Y1 k ^-kk and k = k(v, X) is a single positive root of the equation
v2K2 = ^ ~ „ + Q.
K
2
B
-4 -3
'°9ln(S)
Figure 2. Error behavior in logarithmic scale (all calculations were fulfilled
with parameters A = 0.5, v = 0.5, A2 = 1/8, Am = 0, m = 3, 4,..., the points corresponding to the calculated values are connected by bold lines, thin line corresponds to the theoretical slope): (A) norm of the error for 4-th order Runge-Kutta method versus a number of nodes N of the grid (p = 4, S = 10-3, the length of the interval of integration is equal to « 13); (B) the value of &s
versus S (a = 10, p =3)
Then for the outgoing trajectory the asymptotics takes the form
X,v) ~ eKZ, (z; X,v) ~ xeKZ,
q(z; X,v)
eKZ, qz (z; X,v)
1 - A2k2
Thereby, the point S = S(5) can be chosen as
k5
1- X2k2
S(S) = 8,
(V k5,
k25
0 •
where S is sufficiently small. Consider the Cauchy problem with initial conditions S(S). Introduce the function
Rs(v,X) = qz(zo; X,v),
where (<p(z; X, v), (z; X, v), q(z; X, v), qz(z; X, v)) is a trajectory with initial data S(5) and z0 is such that >^(z0; X,v) = k. Due to the fact that the nonlinearity is odd, the value ARs(v,X) = (v,X) — R(v,X)l depends on 5 not quadratically, but cubically, i.e. ARs(v,X) = 0(53). This fact was verified numerically. For fixed values of v and A the value Rs(v,X) was calculated for the set of 5, 5/a, S/a2, etc with the highest possible accuracy of the Runge-Kutta method (see the point a). Fig. 2, panel B, shows the dependence of the value
6 (v,X)= R (v,X) — (v,X)
on 5/an in the double logarithmic scale, a = 10. It follows from Fig. 2 that this dependence is well approximated by a straight line with the slope coefficient p = 3. Therefore for 5 ^ 0 one has Rs(v,X) = R(v,X) + C53 + o (53), This implies the estimation for the error
6. (v,X) ARs_ (v,X) "
a3- 1
2
K
1 - A2K
2
a
It turns out that for the majority of calculations it is sufficient to take h = 0.0004 and 5 = 10-3, that ensures accuracy ~ 10-9.
4. Numerical Results
Numerical method described in Section 3 was implemented in the program complex "Kink solutions" fulfilled in Matlab environment. This complex allows to calculate the values of sliding velocities vn, n = 0,1,... for 2^-kinks if A is fixed. Also the "Kink solutions" permits to simulate the propagation of the kink being found. The nonlinearity (i.e. CPR law) should be set by the amplitudes of Fourier harmonics A2 - A10, see Eq. (3).
The typical result of the computations is presented in Fig. 3. If the parameter of nonlocality A is fixed, each of found 2^-kinks corresponds to its own sliding velocity vn, n = 1, 2,.... The curves w„(A) corresponding to the three highest velocities v0, v\ and v2 are shown in the left panel of Fig. 3. Two profiles of the 2^-kinks corresponding to points A and B (A = 0.2) are shown in the right panel of Fig. 3. All the curves w„(A) in the left panel of Fig. 3 are originated at the point A = 0, v = 1. The analysis in vicinity of this point can be done using "weak nonlocality" limit, see [5]. One should note that the kinks corresponding to different branches have nearly the same core but differ in asymptotics of the tails.
1
0.9 0.8 0.7 0.6 V 0.5 0.4 0.3 0.2 0.1 0
V\a
........\ \ ...... \V0
.......... <P
V2
0 1 2 3 4 5 6 7
Figure 3. Kink solutions of Eq. (5) for Ag = 1/8, Am = 0, m = 3, 4,.... Values of v„ versus A are shown for the first three solutions. Profiles of the first and third kinks for A = 0.2 are shown in the inserts
5. Conclusion
In this paper, the numerical method for finding 2^-kink solutions of problems of nonlocal Josephson electrodynamics is presented. The model is described by nonlocal equation (5). This equation depends on the parameters Am, m = 2, 3 ..., corresponding to amplitudes of harmonics in current-phase-relation (3) and A (the "strength" of nonlocality). The main difficulty is that the shape of 2^-kink and its velocity v should be computed simultaneously.
The method is based on the ideas of the theory of dynamical systems. In the paper detailed description of the algorithm is presented. The estimation of error of the algorithm is given.
The algorithm has been implemented as a part of the program complex "Kink solutions". The program complex is suitable for finding the set of parameters vn, n = 0,1,..., when A is fixed. The results of numerical calculations with this complex are reported.
References
1. A. A. Abdumalikov, G. L. Alfimov, A. S. Malishevskii, Nonlocal Electrodynamics of Josephson Vortices in Superconducting Circuits, Superconductor Science and Technology 22 (2), art 023001.
2. G. L. Alfimov, V. M. Eleonsky, L. M. Lerman, Solitary Wave Solutions of Nonlocal Sine-Gordon Equations, Chaos 8 (1) (1998) 257-271.
3. Y. M. Aliev, K. N. Ovchinnikov, V. P. Silin, Nonlocal Josephson Electrodynamics of Layered Structures, Journal of Experimental and Theoretical Physics 80 (3) (1995) 551-559.
4. S. Savel'ev, V. A. Yampol'skii, A. L. Rakhmanov, F. Nori, Terahertz Josephson Plasma Waves in Layered Superconductors: Spectrum, Generation, Nonlinear and Quantum Phenomena, Reports on Progress in Physics 73 (2), art 026501.
5. G. L. Alfimov, A. V. Malishevskii, E. V. Medvedeva, Discrete Set of Kink Velocities in Josephson Structures: the Nonlocal Double Sine-Gordon Model, Physica D 282 (2014) 16-26.
6. A. A. Golubov, M. Y. Kupriyanov, E. Il'ichev, The Current-Phase Relation in Josephson Junctions, Reviews of Modern Physics 76 (2) (2004) 411-469.
7. P. K. Atanasova, T. L. Boyadjiev, Y. M. Shukrinov, E. V. Zemlyanaya, Numerical Modeling of Long Josephson Junctions in the Frame of Double Sine-Gordon Equation, Mathematical Models and Computer Simulations 3 (2011) 389-398.
8. G. L. Alfimov, E. V. Medvedeva, Moving Nonradiating Kinks in Nonlocal 4>4 and
— fi6 Models, Physical Review E 84 (5), art 056606.
9. N. N. Kalitkin, P. V. Koryakin, Numerical Methods, Book 2, Methods of Mathematical Physics, Academia, Moscow, 2013, in Russian.
УДК 519.642.2
Численный метод нахождения скоростей скольжения вихрей в нелокальной джозефсоновской электродинамике
Э. В. Медведева
Кафедра высшей математики - 1 Национальный исследовательский университет «МИЭТ» проезд 4806, дом 5, г. Зеленоград, Россия, 124498
В статье рассматривается модель бесконечной джозефсоновской слоистой структуры. Структура состоит из чередующихся сверхпроводящих и туннельных слоёв, при этом предполагается, что (i) электродинамика структуры является нелокальной и (ii) ток-фазовая зависимость представлена в виде суммы гармоник Фурье вместо одной синусоидальной гармоники для случая уравнения синус-Гордона. Основным уравнением модели является нелокальное обобщение нелинейного уравнения Клейна—Гордона с периодической нелинейностью, которое зависит от внешнего параметра нелокальности А. Скорости вихрей (решения типа 2^-кинков) в моделях такого рода не являются произвольными, а принадлежат некоторому дискретному множеству. В работе предложен метод для вычисления таких скоростей (называемых также «скоростями скольжения») и формы кинков. Приводится оценка погрешности этого метода. Результатами вычислений являются семейства решений типа 2^-кинка, параметризуемые значением А. Из результатов численного счёта вытекает, что центральные части профилей 2^-кинков, соответствующих различным семействам при одном и том же значении А, схожи между собой. Отличие наблюдается в асимптотике «хвостов» этих решений. Численный алгоритм был использован в комплексе программ «Kink solutions», написанный в среде MatLab. Комплекс позволяет вычислять формы и скорости решений типа 2^-кинка для нелинейностей, представленных суммой до десяти гармоник Фурье, а также моделировать распространение этих кинков.
Ключевые слова: джозефсоновский переход, нелокальная джозефсоновская электродинамика, вложенные солитоны, скорости скольжения, несинусоидальная нелинейность.
Литература
1. Abdumalikov A. A., Alfimov G. L., Malishevskii A. S. Nonlocal Electrodynamics of Josephson Vortices in Superconducting Circuits // Superconductor Science and Technology. — 2009. — Vol. 22, No 2. — Art 023001.
2. Alfimov G. L., Eleonsky V. M., Lerman L. M. Solitary Wave Solutions of Nonlocal Sine-Gordon Equations // Chaos. — 1998. — Vol. 8, No 1. — Pp. 257-271.
3. Aliev Y. M., Ovchinnikov K. N., Silin V. P. Nonlocal Josephson Electrodynamics of Layered Structures // Journal of Experimental and Theoretical Physics. — 1995. — Vol. 80, No 3. — Pp. 551-559.
4. Terahertz Josephson Plasma Waves in Layered Superconductors: Spectrum, Generation, Nonlinear and Quantum Phenomena / S. Savel'ev, V. A. Yampol'skii, A. L. Rakhmanov, F. Nori // Reports on Progress in Physics. — 2010. — Vol. 73, No 2. — Art 026501.
5. Alfimov G. L., Malishevskii A. V., Medvedeva E. V. Discrete Set of Kink Velocities in Josephson Structures: the Nonlocal Double Sine-Gordon Model // Physica D. — 2014. — Vol. 282. — Pp. 16-26.
6. Golubov A. A., Kupriyanov M. Y., Il'ichev E. The Current-Phase Relation in Josephson Junctions // Reviews of Modern Physics. — 2004. — Vol. 76, No 2. — Pp. 411-469.
7. Numerical Modeling of Long Josephson Junctions in the Frame of Double Sine-Gordon Equation / P. K. Atanasova, T. L. Boyadjiev, Y. M. Shukrinov, E. V. Zemlyanaya // Mathematical Models and Computer Simulations. — 2011. — Vol. 3, No 3. — Pp. 389-398.
8. Alfimov G. L., Medvedeva E. V. Moving Nonradiating Kinks in Nonlocal ф4 and ф4 - ф6 Models // Physical Review E. — 2011. — Vol. 84, No 5. — Art 056606.
9. Калиткин Н. Н., Корякин П. В. Численные методы, Книга 2, Методы математической физики. — Москва: Академия, 2013.
