Numerical investigation of the magnetic flux static distributions in layered josephson junctions

Effective numerical algorithms are worked out for solving the nolinear system of ODE for finding the static distributions of the magnetic flux in N-stacked JJs, as well as the corresponding matrix Sturm-Liouville problem for studying their global stability. The particular case of three stacked JJs is investigated. A correspondence is made between loss of stability of a possible static distribution of the magnetic flux, obtained by solving the static problem, and the switching to dynamic state obtained by solving the dynamic problem. In this work we show by means of numerical simulation that the transient process of switching from static to dynamic state in symmetric three stacked JJs depends on the way of exceeding the external current.


INTRODUCTION
Stacks of long Josephson Junctions (JJs) were intensively studied during the past years.In these systems both nonlinearity and interaction between subsystems play an important role.Such structures make it possible to state and study new physical effects that do not occur in single JJs.One of the most interesting experimental results for two stacked JJs found in resent years is the so-called current locking (CL).The essence of this phenomenon is as follows: there exists a range of the external magnetic field where the different junctions switch to dynamic state simultaneously when the external current exceeds some critical value.It was shown by means of numerical simulation [1] that experimentally found CL for two stacked JJs can be obtained and well explained in the framework of inductive coupling model [2].CL is essentially dynamical phenomenon witch occurs during the complex transient process of magnetic flux penetration into the junctions.However we think that for better understanding of the processes (not only the above mentioned) in real experiments, it is useful to study the possible static distributions of the magnetic flux.For example we can interpret the transitions from static to dynamic state as bifurcations of some stable static solutions under the change of parameters (the applied magnetic field and the external current).

MATHEMATICAL MODEL
A simple scheme of N-stacked JJs (N + 1 superconductor layers divided by N insolating layers) is shown on Figure 1.Black layers are insulators with thickness D and white ones are superconductors with thickness d.Such stack is called symmetric.Main physical quantities that can be measured are the external current, the voltages on the whole system and the voltages of individual junctions.

Numerical investigation of the magnetic flux static distributions in layered josephson junctions
The dynamics of the Josephson phases ϕ (x, t) = (ϕ (x, t), …, ϕ N (x, t)) T (the ith element of ϕ is the phase difference across insulating layer i) in symmetric N-stacked JJs is described by the following system of perturbed sine-Gordon equations [2]: (1) where α is the dissipation coefficient, J = (sin ϕ 1 , …, sin ϕ N ) T is the vector of the Josephson current density, Γ = γ (1, …, 1) T is the vector of the external current density.The matrix L = tridiag(1, s, 1) is the matrix of the inductive interaction between junctions.Here s = −λ L / (D sinh(d/λ L ) + 2λ L cosh(d/λ L )), where λ L is London's penetration depth, 0.5 < s ≤ 0 for arbitrary N.
The system is written in normalized units.Space x is normalized to the Josephson length λ J and time t to the inverse of plasma frequency ω 0 −1 .The definition of dimensionless units can be found in [3].In this work we consider stacks of overlap geometry placed in external magnetic field h e , therefore the system (1.1) should be solved together with the boundary conditions: (2) where H is the vector H = h e (1, …, 1) T and appropriate initial conditions.The existence of Josephson current generates a specific magnetic flux.When the external current is less than some critical value the junctions are in static state and the voltages of all individual junctions are equal to zero.
The static problem corresponding to Eqns (1), ( 2) is To study the global stability of a possible static solution the following Sturm-Liouville problem (SLP) is generated: where .This is equivalent to study the positive definiteness of the second variation of the potential energy of the system.The minimal eigenvalue λ min determines the stability of the distribution under consideration.A minimal eigenvalue equal to zero means a bifurcation caused by change of some parameter, in our case -the external current γ .

NUMERICAL METHOD
In order to solve the nonlinear boundary value problem ( 3), (4) we use an iterative algorithm, based on the continuous analog of Newton's method (CAMN) [4].As initial approximations for the iteration process we take combinations (for the different junctions) of solutions which exist in the one-junction case and h e = 0, γ = 0 : CANM gives a linearized boundary value problem at each iteration step.The linear boundary value problem is solved numerically by means of Galerkin finite element method (FEM) and quadratic finite elements.FEM is used also to reduce the SLP (5), ( 6), (7) to a linear algebraic problem whose few smallest eigenvalues and the corresponding eigenfunctions are found by the subspace iteration method [5].To test the accuracy of the above methods we have used the method of Runge by computing the solutions on sequence of embedded meshes.The numerous experiments made show a super-convergence of order four for both the static problem and the SLP.

NUMERICAL RESULTS
The simplest generalizable model of stacked JJs is the case of three stacked JJs because it takes into account the difference in the behavior of the interior and exterior junctions.The numerical results presented here are for the particular case of three stacked JJs, but the method of investigation and its program realization are developed for the general N-junction case (N > 1).We briefly discuss the numerical results.We investigate numerically the possible solutions of problem (3), (4) and seek for critical values of the external current γ at given applied magnetic field h e .Changing the value of γ for given h e when the geometrical parameters are fixed we get the critical currents for the solution under consideration.The results for some special solutions are shown on Figure 2.
In principal, for a given magnetic field, there may be several allowed static distributions of the magnetic flux differing in number of vortices they contain in the different junctions.Each of these solutions has its own critical current.We don't answer the question when the system can be at a given static state in the real physical experiment.However at low fields all the junctions could be in Meissner state and the results for the transient process of switching from static to dynamic state in this case correspond to the real physical situation.In order to verify the critical current γ cr found by solving the SLP ( 5), ( 6), ( 7) we excite all static solutions shown on Figure 2. by exceeding the external current below and above γ cr by using the dynamical model.Solving the static problem gave us a possibility to use precise initial values for the dynamic problem (1), (2).We have expected that for excitations below γ cr the system will remain in static state and for excitations above γ cr it will switch to dynamic state.The numerical results confirmed our expectations.At external field h e = 0.5 we found interesting result.To analyze the transient process of switching to dynamic state we excited Meissner (MMM) solution in four different ways (Figure 3.).For (MMM) solution the components ϕ 1 , ϕ 3 for the exterior junctions are of the same type.Dew to the symmetry of the stack we have ϕ 1 = ϕ 3 = ϕ ex , ϕ 2 = ϕ in .

262
Numerical investigation of the magnetic flux static distributions in layered josephson junctions 0.72 0.71 0 100 Time -t Excitation γ 1 (t) is under γ cr .Excitations γ 3 (t) and γ 4 (t) are above γ cr , they have the same final value, but they differ in slope.Excitation γ 2 (t) is above γ cr , it has the same slope as γ 4 (t), but lower final value.As we have expected, for γ 1 (t) all the junctions remain in Meissner state (Figures 4, 5).
The numerical experiments confirmed also our expectation that if γ cr is some how exceeded, at least one of the junctions will switch to dynamic state.For excitation γ 4 (t) all three junctions switch to resistive state (Figures 6, 7).For γ 2 (t) and γ 3 (t) only the interior junction switches to resistive state (Figures 8, 9).In all these cases the transient process starts with penetration of fluxons in the interior junction, i.e., the interior junction drives the      transient process.In the case γ 4 (t) the switching of the interior junction to resistive state triggers the switching of the exterior ones, while in the cases γ 2 (t) and γ 3 (t) it does not.

CONCLUSIONS
Numerical technique is proposed and realized for solving the nonlinear stationary problem and the corresponding matrix Sturm-Liouville problem for stack of N inductively coupled LJJs.The particular case of three stacked JJs is investigated.A perfect agreement between the results found     by solving the Sturm-Liouville problem and those found by solving the dynamic problem is established.The numerical simulation shows that the switching from static to dynamic state in symmetric three stacked JJs depends on the way of increasing the external current.

-
Meissner solutions (denoted further by M) of the form ϕ (x) = 2 kπ, k = 0, ± 1, … -fluxon (antifluxon) solutions, for which there are exact analytical expressions in the case of infinite junctions (l = ∞).The single fluxon/antifluxon solution has the well known form ϕ (x) = 4 arctan(exp(±x)) +2 kπ, k = 0, ±1, … Further for n-fluxon distributions we use the simple notation F n , n = ±1, ±2, … For junctions of finite length objects of type F n are not fluxons in a strong sense, but by analogy the same terminology is used.

Figure 3 Figure 2
Figure 3 Increasing the external current in different ways in vicinity of the critical current of Meissner type solution at h e = 0.5 for 2l = 10, α = 0.1, s = −0.3.