Numerical Study of a System of Long Josephson Junctions with Inductive and Capacitive Couplings

The phase dynamics of the stacked long Josephson junctions is investigated taking into account the inductive and capacitive couplings between junctions and the diffusion current. The simulation of the current–voltage characteristics is based on the numerical solution of a system of nonlinear partial differential equations by a fourth order Runge–Kutta method and finite-difference approximation. A parallel implementation is based on the MPI technique. The effectiveness of the MPI/C++ code is confirmed by calculations on the multi-processor cluster CICC (LIT JINR, Dubna). We demonstrate the appearance of the charge traveling wave (CTW) at the boundary of the zero field step. Based on this fact, we conclude that the CTW and the fluxons coexist.


Introduction
The layered high-temperature superconducting materials such as Bi 2 Sr 2 CaCu 2 O 8+δ (BSCCO) can be considered as a stack of coupled intrinsic Josephson junctions (JJ) [1]. This system is one of the promising objects of superconducting electronics [2,3]. Coherent terahertz electromagnetic radiation from this system provides wide possibilities for various applications [4]. However, the mechanism of this radiation is still unclear. The most intense coherent radiation corresponds to the region of the current-voltage characteristic (CVC) where switching from the upper branch to inner branches occurs [3]. In this region of the CVC the Josephson oscillations excite a longitudinal plasma wave due to the parametric resonance [5]. In stacked long Josephson junctions (SLJJ) with the length L larger than the Josephson penetration depth λ J (L > λ J ) the parametric resonance can be realized in the zero field step region [6,7] of the CVC and the fluxons coexist with the longitudinal plasma wave, which can be interpreted as the collective excitation [8]. Other types of the collective states [9] in the stack of coupled JJs are the Josephson plasma [10] and the charge traveling wave (CTW) [ First of all, the investigation of the long JJs stack is an actual problem, because most of the experimental results were obtained for this system. Therefore, the construction of a model that ensures an adequate description of the properties of the SLJJ in the high temperature superconductors is one of the topical tasks of the modern physics of the superconductivity and, also, an actual problem is the construction of effective numerical algorithms for the simulation of the phase dynamics of the SLJJ.
To describe the SLJJ, Sakai, Bodin, and Pedersen [12] proposed a model taking into account the inductive coupling between JJs. Machida and Sakai [13] proposed a generalization of the model to the case of both inductive and capacitive coupling. However, they disregarded the diffusion current [14], the significance of which was emphasized in [15]. Furthermore, the CVC have not yet been studied within the generalized Machida-Sakai model.
In this paper, we investigate the phase dynamics of the SLJJ taking into account the inductive and capacitive couplings [12,13] and the diffusion current [15]. Simulation is based on a numerical solution of a system of nonlinear partial differential equations by a fourth order Runge-Kutta method, finite-difference approximation, and the MPI technique for parallel implementation. The effectiveness of the MPI/C++ code is confirmed by calculations on the multi-processor cluster CICC (LIT JINR, Dubna).

Theoretical model and numerical approach
We use the model with inductive and capacitive couplings [12,13] and with the diffusion current [15]. The SLJJ has a layered structure with N + 1 superconducting and interjacent insulating layers. The x-axis is directed along the length L, the y-axis along the width of the superconducting layers and the z-axis is perpendicular to the superconducting layers. Each superconducting layer is described by the Ginzburg-Landau order parameter Δ l = |Δ 0 | exp (iθ l ). The l-th and the (l − 1)-th superconducting layers form the l-th JJ and it is described by the gauge-invariant phase difference (1) of the Ginsburg-Landau order parameter [13] where θ l is the phase of the order parameter of the l-th superconducting layer, e is the electrical charge, is the Plank constant, c is the speed of the light in vacuum and A z is the vector potential. The system of equations which describes the phase dynamics of the coupled long JJs stack in terms of normalized quantities can be written as follows: where , λ e -the Debye screening length, d I and d s denote the thickness of the insulating and the superconducting layers, s C = −λ e /[d I sinh(d s /λ e )] is the capacitive coupling parameter, V l is the voltage of the l-th JJ normalized to V 0 = ω p /(2e), ω p = 8πd I e j c /( ε) is the plasma frequency of JJ, j c is the critical current of JJ, ε is the dielectric constant of the insulating layer, β = σV 0 /d I j c is the dissipation parameter, σ is the conductance of the JJ, I is the bias current normalized to the critical current j c . Here £ is the matrix of the inductive coupling. The main tasks are to calculate the CVC, the spatio-temporal dependence of the magnetic field in the JJs, and the electric charge in the superconducting layers. First of all, we should solve numerically the system of partial differential equations (2). We introduce the uniform mesh with the stepsize Δx in the coordinate x along the JJ and the stepsize Δt in time. We put Δt = Δx/5 in accordance with the Courant-Friedrichs-Lewy condition in order to provide the stability of the numerical scheme. We denote the discrete coordinate by x i = Δx × (i − 1), where i = 1, . . . , N x and N x = L/Δx + 1 is the number of coordinate nodes. The discrete time is denoted by t j = Δt × j, where j = 0, 1, 2, . . . N t . The x = 0 corresponds to x 1 and x = L -to x N x . In the same way, t = 0 corresponds to t 1 and t = T max to t N t , where T max is the end of the time domain. The standard second order finite difference approximation in the spatial coordinate x is used,

EPJ Web of Conferences
Δx .
The resulting system of ordinary differential equations is solved numerically for a fixed value of the current I by a 4th-order Runge-Kutta (RK) algorithm (here l is the JJ number) in the interval [0, L] along the x-coordinate and [0, T max ] in time and obtain the ϕ l (x, t) and V l (x, t) as functions of x and t. Next, the obtained V l (x, t) is averaged with respect to the coordinate x usinḡ and with respect to the time t using where T min denotes the beginning of the averaging interval. The total voltage of the JJs stack can be calculated using V = N l=1 V l . The integrals (3) and (4) are calculated using the Simpson method and the rectangles method, respectively. Then the bias current value is changed by ΔI and the above procedure is repeated. In the calculations the bias current increases from the start value I = 0.01 to I = I max and then decreases to I = 0. In order to investigate collective excitations in the JJs stack, like the longitudinal plasma wave [5] or charge travelling wave [11], the electric charge-time dependence in the superconducting layers is to be calculated. In this case, the electric charge, normalized to Q 0 = εV 0 /4πd s d I , is calculated as a function of x and t using the expression Q l ( [8]. Then the value of Q l (x, t) is averaged with respect to the coordinate x using the Simpson method. For the external current value corresponding to the fluxon states, the magnetic field B in the JJs is calculated using the expression The magnetic field is normalized to B 0 = c/2eD L λ J . The parallel algorithm is based on a distribution of the calculations in the coordinate nodes x i between the group of P m parallel MPI-processes, where m = 0, 1 . . . , M. At each time-step t j , each process P m calculates the RK coefficients and V l ( At each t j the exchange between neighbor processes is arranged: each process P m (m < M −1) sends the RK coefficients and values of V and ϕ at i = i max −1th point to the P m+1 -process; each P m process (m > 0) sends the RK coefficients and solutions at i = i min to the P m−1 -process. In order to calculate the average value V l , the parallel calculation of the integral (3) is performed at each time-step t j . Each P m -process calculates the partial sum of elements V l (x i , t j ) at each JJ l in accordance to the Simpson quadrature formula. Then the resulting summation is performed in the process P 0 . In the P 0 -process V l is averaged in time and in LJJs, and the resulting value is saved to the file. For some values of I the solutions V l (x i , t j ) and φ l (t i , t j ) are collected in the process P 0 where they are saved to the file together with the respective physical characteristics.

Results and discussions
Let us first of all discuss the effectiveness of the parallel algorithm. The methodical calculations on the CICC multi-processor cluster (LIT JINR) with a different number of parallel MPI-processes are performed. For these calculations we put the number of JJs N = 10, the JJ length L = 5 and L = 10; Δx = 0.05; ΔI = 0.005. Table 1 shows the calculation time (in minutes) of the CVC per number of processes for the stacks with N = 3, N = 10 JJs and with length L = 5, L = 10. One can see that the developed parallel algorithm provides 6-7 times acceleration (depending on the values of N and L) as compared to the sequential simulation.
The next step of the parallel optimization of the long JJs stack simulation code is the parallelization of the calculations at each junction.
Let us now discuss the main features of the single long JJ. Figure 1(a) shows the one loop CVC of the single JJ with L = 10. The calculations have been performed for the β = 0.2, B ext = 0, ΔI = 0.0001 and T max = 200. The stepsize in the coordinate and time were Δx = 0.05 and Δt = Δx/5, respectively. The curve shows the hysteresis and seven step structure at decreasing current. These zero-field steps (ZFS) correspond to the formation of fluxons in the JJs [6,7]. The velocity of the fluxon u is given by the expression [6,7]  of n steps in the CVC [6,7]. In order to directly demonstrate fluxons, we calculate the spatiotemporal dependence of the magnetic field in JJ. Figure 1 Now we discuss the numerical results obtained for the SLJJ. The CVC of the stack with N = 10 JJ and L = 5 is represented in Fig. 2(a), which demonstrates the hysteresis and three ZFS. In order to study in more detail the collective dynamics of the SLJJ, we have calculated the time dependence of the averaged charge with respect to the x-coordinate in the region of the CVC, corresponding to the third ZFS. This region of the CVC is shown by the hollow arrow in Fig. 2(a). The chargetime dependence is shown in Fig. 2(b). The inset in this figure demonstrates the character of the charge oscillation. Analyzing the numerical results we conclude that, at the boundaries of the third ZFS, a charge traveling wave (CTW) is formed. The resonance between the CTW and the Josephson Mathematical Modeling and Computational Physics 2015 02038-p.5 oscillations increases the charge amplitude, which can be seen in these figures. Since the CTW appears in the CVC region corresponding to the ZFS, we may conclude about the coexistence of the CTW and the fluxons.

Conclusions
In this paper, we have investigated numerically the phase dynamics of the SLJJ taking into account the inductive and capacitive couplings between junctions and diffusion current. In our investigation we used the parallel and sequential calculations of CVC. The parallel implementation is based on the MPI technique. We showed that the parallel algorithm results in a seven time acceleration in comparison with the sequential one. We compared the numerical and theoretical results for a single long JJ and obtained a good agreement. For the case of the SLJJ we demonstrated that, at the boundaries of the zero field step, a charge travelling wave is formed. Based on this fact, we concluded that the CTW and the fluxons coexist. This effect is confirmed by a number of calculations with different values of the LJJ length. Due to the resonance between the CTW and the Josephson oscillations, the charge amplitude in the superconducting layers increases.