IVC Calculation Problem for Josephson Junction Stacks. On Asymptotic Construction near the Breakpoint

Solving the system of 𝑛 essentially nonlinear differential equations for different 𝐼 we obtain the current-voltage characteristics (IVC) for a system of 𝑛 Josephson junctions (JJ) as a hysteresis loop. When the current 𝐼 approaches on the back way the breakpoint 𝐼 𝑏 the voltage 𝑉 ( 𝐼 ) falls sharply to zero. In addition, in numerical modelling (non-periodic boundary conditions (NPBC)) IVC multiple branching is observed near 𝐼 𝑏 . It is interesting to study this phenomenon analytically developing asymptotic methods. There had been developed simple “asymptotic” formulas suitable for calculation of all IVC points except near to 𝐼 𝑏 . A numerical-analytical method allowing to shorten IVC calculation time essentially was proposed. This method showed good results in IVC multiple branching calculation in particular. All calculations were performed using the REDUCE system. We succeeded first to calculate analytically all points of IVC. An approximate solution at the breakpoint region (periodic boundary conditions (PBC)) has been developed using the Bogolyubov–Krylov method.


Introduction
The definition of the singular points of the current voltage characteristics together with the estimation of the width of their influence region provide adequate input for physical experiments aiming at studying the finite JJ stacks [1][2][3].
Here˙is derivative of of , ⊂ [0, max ]. Solving the Cauchy problem for different : = 0 + ∆ max (right way) and = max − ∆ we obtain IVC of Josephson junctions as hysteresis loop. The points ( , ( )), with corresponding to the right way, form the right branch of the loop. And the points ( , ( )), with corresponding to the back way, form the back branch of the loop. ( ) is the total voltage of the stack. For = 0 the system (1) is solved with zero initial data. For each next : = +1 , found ( , max ),˙( , max ) are used as initial data.
The coefficients , ′ of the system (1) are elements of the matrix . In the case of PBC is symmetric square matrix of order : And in the case of NPBC is symmetric square tridiagonal matrix of order : The dynamics of phase differences ( ) had been simulated [4] by solving the equation system (1) using the fourth order Runge-Kutta method. The aim was to shorten the time of IVC calculation, and it was achieved. First the number of the system equations was reduced to one in the case of PBC and halved in the case of NPBC. The longtime "asymptotics" formulas were developed. The mixed numerical-analytical method was suggested: all points of IVC are calculated using the asymptotics except points of the little arc of the back branch of the loop calculated numerically [5]. And Fig. 5 shows that the Runge-Kutta method of the fourth order accuracy can be replaced by the simplest explicit cross-scheme of the second order accuracy. This shortens the time of IVC calculation complimentarily.

On Long-Time "Asymptotics" Construction
In the case of PBC IVC calculation for a stack of intrinsic Josephson junctions is reduced to solving a unique equation [5] ( ) = −˙( ) − sin( ( )) + .
Solving this equation with given initial data (0) = 1 ,˙(0) = 2 is equivalent to solving the following integral equation The simple iteration method starting from zero gives on the second step And the total voltage of the stack of Josephson junctions is calculated as In Fig. 1, Fig. 2, Fig. 3 the solid curve is the IVC back branch calculated using the Runge-Kutta method. In Fig. 1 the pictures of the back way of the hysteresis loop for 9 Josephson junctions are shown. The circles refer to "asymptotic" (using (5)) calculations respectively. Here we put min = 50, max = 1000, = = 0.2, ∆ = 0.05, the step in the Runge-Kutta method ℎ = 0.1. In Fig. 2 the circles refer to calculation performed by the following mixed analyticalnumerical method. The right branch of the hysteresis loop and the back branch on the interval 0.4 < < 1.45 were computed using (5). The rest ones were computed numerically. The REDUCE system [6] was used in all calculations performed.
We proved [7] that in the case of NPBC with = 1 the problem of hysteresis loop calculation reduces to solving the following system of integral equations: The system (6) was solved using simple iterations starting at zero. The results obtained on the second step are regarded as long-time "asymptotics" of the system (6) solution.
We constructed "asymptotics" following to [8]. In the introduction N. Levinson declared: "A method is given for showing that formal "approximate solutions" of nonlinear differential equations are in fact the leading terms in an asymptotic representation of actual solutions".
The Fig. 3, Fig. 4 confirm this declaration. At the same time these pictures show efficiency of the suggested mixed numerical-analytical method. This time ∆ = 0.001. Fig. 3 refers to numerical calculation: all points of IVC were calculated using the Runge-Kutta method. Fig. 4 refers to calculation performed by the mixed numerical-analytical method: the right branch of IVC and the part of the back branch (1.45 > 0.4) were calculated analytically (using the long-time "asymptotics" of (6) system solution). The rest points (0.4 0.2) were calculated numerically using the Runge-Kutta method.   (7), obtained by Bogolyubov-Krylov method The Fig. 5 is similar to Fig. 4, but here the simplest explicit cross-scheme of the second order accuracy instead of the fourth order accuracy Runge-Kutta method was used. Remark that in both cases ∆ = 0.1 was chosen.
Following step by step to [9], we obtain 0 and are determined by given initial data: Let = 0 cos( ), = 0 sin( ), 2 0 = 2 + 2 , the first equation of (5) gives = 1 − 0 , and is determined of the third order polynomial equation This equation has at least one real root. When , are found, 0 = √︀ 2 + 2 and = arcsin( / 0 ). By such a way we could correct Fig. 1. In Fig. 6 the back branches of the hysteresis loop for the case of PBC are presented.
The solid curve is the same as in Fig. 1. The circles on the solid curve were found analytically using mixed analytical method: all points of the hysteresis loop were calculated using (5) -except two points ( = 0.3, 0.25) calculated using (7). Вычисление ВАХ для систем джозефсоновских переходов.