# Single fluxon in double stacked Josephson junctions: Analytic solution.

## Abstract

We derive an approximate analytic solution for a single fluxon in a double stacked Josephson junctions (SJJ’s) for arbitrary junction parameters and coupling strengths. It is shown that the fluxon in a double SJJ’s can be characterized by two components, with different Swihart velocities and Josephson penetration depths. Using the perturbation theory we find the second order correction to the solution and analyze its accuracy. Comparison with direct numerical simulations shows a quantitative agreement between exact and approximate analytic solutions. It is shown that due to the presence of two components, the fluxon in SJJ’s may have an unusual shape with an inverted magnetic field in the second junction when the velocity of the fluxon is approaching the lower Swihart velocity.

Magnetic field, $`H`$, parallel to layers, penetrates stacked Josephson junctions (SJJ’s) in the form of fluxons. Fluxons in SJJ’s are different both from Abricosov vortices in type-II superconductors, since they don’t have a normal core, and from Josephson vortices in a single Josephson junction (JJ’s), since magnetic field is spread over many JJ’s. Despite of a large body of works concerning vortex related properties of layered superconductors, for review see e.g. Ref. , exact solution for a single fluxon in SJJ’s is still missing. Distribution of magnetic induction, $`B`$, at large distances from a fluxon was studied in Refs.  and , for tightly packed stacks with identical and nonidentical layers, respectively. In Ref.  a solution for weekly coupled double SJJ’s was obtained using perturbation approach. However, this solution is not valid for strongly coupled SJJ’s, which is the most interesting case. The electromagnetic coupling is strong when interlayer space periodicity of the stack, $`s`$, is (much) less than the effective London penetration depth, $`\lambda _{ab}`$, as it is the case for high-$`T_c`$ superconductors.

Recently, a simple approximate analytic solution for double SJJ’s was suggested. The approximate solution was shown to be in a good agreement with numerical simulations in the static case for arbitrary coupling strengths and parameters of the stack. However, the accuracy of the solution was not strictly proven. Another question which is still opened is how the shape of the fluxon is changed in the dynamic case. This is a crucial question since fluxon dynamics determines perpendicular ($`c`$-axis) transport properties of SJJ’s.

In this paper we rigorously derive the approximate analytic solution for a single fluxon in a double SJJ’s, which is valid for arbitrary junction parameters and coupling strengths. Using the perturbation theory we find the second order correction to the solution and analyze its accuracy. Comparison with direct numerical simulations shows a quantitative agreement between exact and approximate analytic solutions. It is shown that the fluxon in a double SJJ’s can be characterized by two components, with different Swihart velocities and Josephson penetration depths. We also studied the transformation of the fluxon shape with increasing propagation velocities. It is shown that due to the presence of two components, the fluxon in SJJ’s in the dynamic case may have unusual shape with an inverted magnetic field in the second junction at high propagation velocities.

We consider a double SJJ’s with the overlap geometry, consisting of JJ’s 1 and 2 with the following parameters: $`J_{ci}`$ -the critical current density, $`t_i`$\- the thickness of the tunnel barrier between the layers, $`d_i`$ and $`\lambda _{Si}`$ \- the thickness and London penetration depth of superconducting layers. Hereafter the subscript $`i`$ on a quantity represents its number. The strength of electromagnetic coupling of junctions is determined by the coupling parameter, $`S`$, which varies from 0 to 0.5 for the stack with identical JJ’s. More details about definitions can be found in Ref. . We consider frictionless fluxon motion with a constant velocity, $`u`$. The fluxon will be placed in junction 1.

The physical properties of SJJ’s are described by coupled sine-Gordon equation (CSGE) for gauge invariant phase differences, $`\phi _i`$. The problem with solving CSGE is the coupling of nonlinear $`sin\left(\phi _i\right)`$ terms. To decouple the variables, we take linear combination of equation in CSGE and rewrite it as:

$$\stackrel{~}{\lambda }_{1,2}^2F_{1,2\xi \xi }^{\prime \prime }sin\left(F_{1,2}\right)=Er_{1,2}\left(\xi \right),$$
(1)   
where $`\xi =xut`$ is the self-coordinate of the fluxon,

$$F_{1,2}=\phi _1\kappa _{2,1}\phi _2,$$
(2)   
$`Er_{1,2}`$ $`=`$ $`sin(\phi _1)\kappa _{2,1}sin(\phi _2)sin(\phi _1\kappa _{2,1}\phi _2)`$ (3)   
$``$ $`\kappa _{2,1}\phi _2\left(1cos(\phi _1)\right)+O\left(\phi _2^3\right).`$ (4)   
Here $`\stackrel{~}{\lambda }_{1,2}^2=\lambda _{1,2}^2(1u^2/c_{1,2}^2)`$ and coefficients $`\kappa _{1,2}`$, characteristic Josephson penetration depths, $`\lambda _{1,2}`$, and characteristic Swihart velocities, $`c_{1,2}`$ are given by Eqs.(17,20,21) from Ref.. Such choice for coefficients $`\kappa _{1,2}`$ minimizes $`Er_{1,2}`$ far from the fluxon center. The phase differences should satisfy boundary conditions:

$`\phi _1(\mathrm{})`$ $`=`$ $`0,\phi _1(0)=\pi ,\phi _1(+\mathrm{})=2\pi ;`$ (5)   
$`\phi _2(\pm \mathrm{},0)`$ $`=`$ $`0.`$ (6)   
Therefore, functions $`Er_{1,2}`$ have a form of ripple around zero value and will be considered as perturbation. Solutions of Eq.(1) can be easily found. Solutions of uniform Eq.(1), i.e. with zero r.h.s., are:

$$F_{1,2}=4arctan\left[exp\left(\xi /\stackrel{~}{\lambda }_{1,2}\right)\right],$$
(7)   
From Eq.(2) we obtain the first approximation for phase differences:

$`\phi _1`$ $`=`$ $`{\displaystyle \frac{\kappa _1F_1\kappa _2F_2}{\kappa _1\kappa _2}},`$ ()   
$`\phi _2`$ $`=`$ $`{\displaystyle \frac{F_1F_2}{\kappa _1\kappa _2}},`$ ()   
which coincides with the approximate analytic solution, obtained in Ref. . The approximate solution is asymptotically correct at large distances from the fluxon center and has correct values at $`x=0`$, as follows from Eqs. (3,4).

Next, we look for a solution of nonuniform Eq.(1) in the form $`F_{1,2}=F_{1,2}+\delta F_{1,2}`$, where $`\delta F_{1,2}`$ are corrections due to perturbation terms $`Er_{1,2}`$:

$$\left(\stackrel{~}{\lambda }_{1,2}^2\frac{d^2}{d\xi ^2}1+\frac{2}{cosh^2\left(\xi /\stackrel{~}{\lambda }_{1,2}\right)}\right)\delta F_{1,2}=Er_{1,2},$$
(8)   
with boundary conditions $`\delta F_{1,2}(\pm \mathrm{},0)=0`$. The solution of Eq. (7) is

$$\delta F_{1,2}=a_{1,2}(\xi )f_{1,2}+b_{1,2}(\xi )g_{1,2},$$
(9)   
where

$`f_{1,2}`$ $`=`$ $`1/cosh(\xi /\stackrel{~}{\lambda }_{1,2}),`$ (10)   
$`g_{1,2}`$ $`=`$ $`sinh\left({\displaystyle \frac{\xi }{\stackrel{~}{\lambda }_{1,2}}}\right)+{\displaystyle \frac{\xi }{\stackrel{~}{\lambda }_{1,2}cosh(\xi /\stackrel{~}{\lambda }_{1,2})}},`$ ()   
are partial solutions of the uniform Eq.(7) and

$`a_{1,2}`$ $`=`$ $`{\displaystyle \frac{1}{2\stackrel{~}{\lambda }_{1,2}}}{\displaystyle _0^\xi }Er_{1,2}\left(x^{}\right)g_{1,2}\left(x^{}\right)𝑑x^{},`$ (11)   
$`b_{1,2}`$ $`=`$ $`{\displaystyle \frac{1}{2\stackrel{~}{\lambda }_{1,2}}}{\displaystyle _{\mathrm{}}^\xi }Er_{1,2}\left(x^{}\right)f_{1,2}\left(x^{}\right)𝑑x^{}.`$ ()   
The perturbation corrections to the approximate fluxon solution, Eq.(6) are

$`\delta \phi _1`$ $`=`$ $`{\displaystyle \frac{\kappa _1\delta F_1\kappa _2\delta F_2}{\kappa _1\kappa _2}},`$ ()   
$`\delta \phi _2`$ $`=`$ $`{\displaystyle \frac{\delta F_1\delta F_2}{\kappa _1\kappa _2}}.`$ ()   
The corrections $`\delta \phi _{1,2}`$ are used to improve the approximate solution Eq.(6) and to estimate it’s accuracy.

Lets observe that from Eq. (3), $`\kappa _1Er_1\kappa _2Er_2`$. From Eqs.(9-11) it can be seen that for $`\stackrel{~}{\lambda }_2=\stackrel{~}{\lambda }_1`$, $`\kappa _1\delta F_1=\kappa _2\delta F_2`$ and $`\delta \phi _10`$ for arbitrary $`\phi _20`$. Moreover, taking into account Eq.(6 b), it can be shown that $`d(\delta \phi _1)/d\stackrel{~}{\lambda }_{2(\stackrel{~}{\lambda }_2=\stackrel{~}{\lambda }_1)}=0`$.

In Fig. 1, central part of a fluxon is shown for a stack of two identical JJ’s with strong coupling, $`S`$=0.495, for the static case, $`u`$=0. Phase distribution in the full scale is shown in Fig. 3. Parameters of the stack are: $`d_{13}=t_{1,2}=0.01\lambda _{J1}`$, $`\lambda _{S13}=0.1\lambda _{J1}`$, where $`\lambda _{J1}`$ is the Josephson penetration depth of the single junction 1. Solid lines represent results of direct numerical integration of CSGE, Eq. (1), dashed lines show the approximate analytic solution, Eq.(6), and dotted line shows the corrected analytic solution, $`\phi _{2a}+\delta \phi _{2a}`$, Eqs.(6,11), in junction 2. Solid and dashed curves, marked as $`\delta \phi _{1,2}`$, represent the overall discrepancy between numerical and approximate analytic solutions and the perturbation correction, Eq.(11), respectively.

From Fig. 1 it is seen that correction to the fluxon image, $`\phi _2`$, in the second junction vanishes far from and in the fluxon center, while the accuracy decreases at distances $`\lambda _{J1}`$ from the center. Such behavior is expected from the shape of perturbation functions $`Er_{1,2}`$ in Eq. (1). On the other hand, for junction 1, correction $`\delta \phi _1`$ is small in the whole space region and analytic solution gives an excellent fit to the ”exact” numerical solution, in agreement with discussion above. The most crucial test for the approximate solution is the accuracy of derivative at $`x=0`$

$$\frac{\delta \phi _1^{}(0)}{\phi _1^{}(0)}=\frac{\kappa _1b_1(0)/\stackrel{~}{\lambda }_1\kappa _2b_2(0)/\stackrel{~}{\lambda }_2}{\kappa _1/\stackrel{~}{\lambda }_1\kappa _2/\stackrel{~}{\lambda }_2}.$$
(12)   
An estimation for $`u=0`$ yields

$$\frac{\delta \phi _1^{}(0)}{\phi _1^{}(0)}\frac{2\alpha \kappa _1\kappa _2\left(\stackrel{~}{\lambda }_1\stackrel{~}{\lambda }_2\right)^2\left(\stackrel{~}{\lambda }_1^1+\stackrel{~}{\lambda }_2^1+\frac{\lambda _0}{\stackrel{~}{\lambda }_1\stackrel{~}{\lambda }_2}\right)}{\left(\kappa _1\kappa _2\right)^2\left(\stackrel{~}{\lambda }_1+\stackrel{~}{\lambda }_2+\frac{2\stackrel{~}{\lambda }_1\stackrel{~}{\lambda }_2}{\lambda _0}\right)\left(1+\frac{\stackrel{~}{\lambda }_1}{\lambda _0}\right)\left(1+\frac{\stackrel{~}{\lambda }_2}{\lambda _0}\right)},$$
(13)   
where $`\alpha `$ is a factor of the order of unity and $`\lambda _0`$ given by Eq.(26) from Ref.. From Eq. (13) it is seen that both $`\delta \phi _1^{}(0)/\phi _1^{}(0)`$ and it’s derivative with respect to $`\stackrel{~}{\lambda }_2`$ goes to zero at $`\stackrel{~}{\lambda }_2=\stackrel{~}{\lambda }_1`$, as discussed above.

In Fig. 2, the maximum of $`\delta \phi _1\left(x\right)`$ (top panel) and the relative correction to derivative at $`x=0`$, $`\delta \phi _1^{}(0)/\phi _1^{}(0)`$, (bottom panel) for $`u`$=0 are shown as a function of $`J_{c2}/J_{c1}`$ for four different coupling parameters $`S`$=0.495 (solid lines), $`S`$=0.433 (dashed lines), $`S`$=0.312 (dotted lines) and $`S`$=0.127 (dashed-dotted lines). For $`S`$=0.495 parameters of the stack are the same as in Fig.1; for $`S`$=0.433 $`d_i=0.5\lambda _{Si}`$; for $`S`$=0.312 $`d_i=\lambda _{Si}`$; and for $`S`$=0.127 $`d_i=2\lambda _{Si}`$. From Fig. 2 it is seen that the accuracy of solution improves with decreasing $`S`$. This is naturally explained by a decrease of $`\phi _2`$, see Eq.(3) and decrease of splitting between $`\lambda _{1,2}`$, see Eq. (20) from Ref.. However, even for strongly coupled case, the analytic solution, Eq.(6), gives quantitatively good approximation not only for the value, but also for the derivative of $`\phi _1`$ for arbitrary parameters of the stack. The gray solid line in the bottom panel of Fig. 2 shows the estimation for $`\delta \phi _1^{}(0)/\phi _1^{}(0)`$, calculated from Eq.(13) for $`S`$=0.495. It is seen, that estimation gives qualitatively correct result. Namely, $`\delta \phi _1^{}(0)/\phi _1^{}(0)`$ vanishes both for $`J_{c2}/J_{c1}0`$, as $`J_{c2}/J_{c1}`$, and for $`J_{c2}/J_{c1}\mathrm{}`$, as $`\sqrt{J_{c1}/J_{c2}}`$. Note, that for $`J_{c2}/J_{c1}0`$, $`\delta \phi _1`$ vanishes even though the splitting of $`\lambda _{1,2}`$ becomes extremely large.

Sofar we have considered the static case, $`u`$=0. On the other hand, according to Eq. (6), radical changes should take place in the dynamic state, and the shape of the fluxon in SJJ’s may become qualitatively different from that in the single Josephson junction. Indeed, as the velocity approaches the lower characteristic velocity, $`uc_1`$, $`\stackrel{~}{\lambda }_10`$, i.e. the $`F_1`$ component of the fluxon contracts, while contraction of the $`F_2`$ component remain marginal. This implies, that at $`uc_1`$, the fluxon in SJJ’s consists of a contracted core and uncontracted ”tails” decaying at distances much larger than the core size. Such behavior is clearly different from that in a single Josephson junction. We note that characteristic velocities, $`c_{1,2}`$ may depend on $`u`$, therefore, contraction of each component, $`F_{1,2}`$, may be different from Lorentz contraction. Another interesting consequence of the approximate analytic solution, Eq.(6), is that with increasing $`u`$, the magnetic field in the second junction may change the sign with respect to that in junction 1. Such behavior was predicted in Ref.  from the approximate analytic solution, Eq. (6), and it was suggested, that this will result in attractive fluxon interaction in SJJ’s. The magnetic induction in SJJ’s is equal to

$`B_1`$ $`=`$ $`{\displaystyle \frac{H_0\lambda _{J1}}{2\left(1S^2\right)}}\left[\phi _1^{}+S\sqrt{{\displaystyle \frac{\mathrm{\Lambda }_2}{\mathrm{\Lambda }_1}}}\phi _2^{}\right],`$ ()   
$`B_2`$ $`=`$ $`{\displaystyle \frac{H_0\lambda _{J1}}{2\left(1S^2\right)}}\left[S\sqrt{{\displaystyle \frac{\mathrm{\Lambda }_1}{\mathrm{\Lambda }_2}}}\phi _1^{}+{\displaystyle \frac{\mathrm{\Lambda }_1}{\mathrm{\Lambda }_2}}\phi _2^{}\right],`$ ()   
where $`H_0=\frac{\mathrm{\Phi }_0}{\pi \lambda _{J1}\mathrm{\Lambda }_1}`$ and $`\mathrm{\Lambda }_{1,2}`$ are defined in Ref. .

From Eq.(14) it is seen that estimation of magnetic induction in SJJ’s requires the accuracy of derivatives in both junctions, while Eq.(6) is not valid with the accuracy of derivative at $`x=0`$ for $`\phi _2`$. Therefore, more elaborate analysis is needed for the study of magnetic field distributions in the fluxon.

In Fig. 3, we show a) phase distributions and b) magnetic field distributions in the fluxon for different fluxon velocities. Parameters of the stack are the same as in Fig. 1. Solid lines in Fig. 3 a) represent results of direct numerical simulations of CSGE, Eq.(1), and dotted lines show the approximate analytic solution, Eq.(6). It is seen that quantitative agreement between ”exact” and approximate solutions sustain up to $`c_1`$. For the case of identical junctions, considered here, exactly one half of the fluxon belongs to each of the components, $`F_{1,2}`$. Indeed, from Fig. 3 a) it is seen that for $`uc_1`$ there is a contracted core at $`x=0`$ with a one $`\pi `$ step in $`\phi _1`$. On both sides of the core, there are two $`\pi /2`$ tails, which are slowly decaying at distances $`\stackrel{~}{\lambda }_2\stackrel{~}{\lambda }_1`$. Solid and dashed curves in Fig. 3 b) represent numerically simulated profiles, $`B_{1,2}\left(x\right)`$, in junctions 1 and 2, respectively. From Fig. 3 b) it is clearly seen, that with increasing fluxon velocity, a dip in $`B_2`$ develops in the center of the fluxon. At velocities close to $`c_1`$, $`B_2\left(0\right)`$ changes sign and finally at $`u=c_1`$, $`B_2\left(0\right)=B_1\left(0\right)`$. Such behavior is in agreement with predictions of Ref. . From numerical simulations we have found, that the fluxon shape in double SJJ’s is well described by Eq.(6) up to at least $`u0.98c_1`$ for any reasonable parameters of the stack, although the accuracy of the approximate solution may decrease with increasing $`u`$. The decrease of accuracy is caused by the increase of $`\phi _2`$ as is seen from Fig. 3. In this case, perturbation correction, Eq. (11), should be taken into account.

In conclusion, a simple approximate analytic solution for a single fluxon in a double stacked Josephson junctions for arbitrary junction parameters and coupling strengths is derived. It is shown that the fluxon in a double SJJs can be characterized by two components, with different Swihart velocities and Josephson penetration depths. Using the perturbation theory we find the second order correction to the solution and analyze its accuracy. Comparison with direct numerical simulations shows a quantitative agreement between exact and approximate analytic solutions for all studied parameters of the stack and fluxon velocities up to at least 0.98$`c_1`$. It is shown that due to the presence of two components, the fluxon in SJJ’s may have an unusual shape with an inverted magnetic field in the second junction at large propagation velocities. This may lead to attractive fluxon interaction in the dynamic state of SJJ’s

Discussions with D.Winkler are gratefully acknowledged. The work was supported by the Russian Foundation for Basic Research under Grant No. 96-02-19319.