# Neighbor-junction state effect on the fluxon motion in a Josephson stack

## I Introduction

Stacked long Josephson junctions (LJJ’s) have recently received much attention since they show a variety of new physical phenomena  in comparison with single LJJ’s and have potential for applications as a narrow linewidth powerful oscillators for mm and sub-mm wavebands. The naturally layered high-$`T_c`$ superconductors (HTS) can be described as intrinsic stacks of Josephson junctions . Therefore, study of fluxon dynamics in artificial stacks can help to understand the phenomena which take place in HTS.

The inductive coupling model describing the dynamics of Josephson phases in $`N`$ inductively coupled LJJ’s was derived by Sakai et al.. Experimental investigation of stacked junctions became possible after the progress achieved in (Nb-Al-AlO<sub>x</sub>)<sub>N</sub>-Nb technology  which, at the present stage, allows to fabricate stacks with up to about 30 Josephson tunnel junctions having parameter spread between them of less than 10% . Initially, the interest was concentrated on investigation of the simplest symmetric fluxon states since they are promising for oscillator applications. Later on, it was found that it is very interesting to understand the asymmetric states because they show rather nontrivial nonlinear dynamics. Such asymmetric states are also of practical importance, because multilayered oscillators most probably will operate in a regime when only the majority but not all of the junctions oscillate coherently while the other junctions are in resistive or not synchronized flux-flow state.

In a recent work  it has been shown that the dynamic state of one junction in a 2-fold stack affects the static properties of the other junction. As a next step, it is interesting and important to understand how different dynamic states in one LJJ affect the dynamics of fluxons in the other LJJ. In particular, the goal of this work is to study the dynamics of a fluxon in one LJJ when the neighboring LJJ is in the resistive state and compare it with the case when the neighboring LJJ is in Meissner state. We call the resistive state a “phase-whirling” state because the Josephson phase difference rotates very fast and nearly uniformly, in rough approximation. Such a dynamic state often occurs in experiments and, therefore, it is important to understand and describe it adequately. In fact, in the early experiments with stacks it was somewhat naively supposed that the voltage of Flux-Flow Step (FFS) in one LJJ does not depend on the state (Meissner or resistive) of the other LJJ. In fact this is true only for the asymptotic voltage of FFS. Here we show that in the presence of the “phase-whirling” solution in one of the junctions the actual flux-flow voltage across the other LJJ gets lower.

In the section II we present the experimental data which clearly show that in a 2-fold stack the switching of one junction from the Meissner state to the phase-whirling state decreases the velocity of a fluxon moving in the other junction and, therefore, the dc voltage across it. The analytical approach which explains the observed decrease of fluxon velocity in the limit of high fluxon density is developed in section III. The results of numerical simulations confirm analytical results and are shown in section IV. The results of the work are summarized in section V.

## II Experiment

In order to investigate the influence of the phase-whirling state in the neighboring junction on the fluxon dynamics, we have chosen the most clean ring-shape (annular) LJJ stack geometry. Due to magnetic flux quantization in a superconducting ring, the number of fluxons initially trapped in each annular junction of the stack is conserved. The fluxon dynamics can be studied here under periodic boundary conditions which exclude possible complicated interference of the fluxon with the junction edges.

Experiments have been performed with 3 different (Nb-Al-AlO<sub>x</sub>)<sub>2</sub>-Nb annular LJJ stacks prepared in 2 different technological runs (2 samples in one run and the third sample in another run). The sample geometry is shown in Fig. 1. Two annular LJJ’s are stacked one on top of the other, with bias leads attached to the top and bottom electrodes. The physical parameters of all samples, measured at $`T=4.2\mathrm{K}`$, are summarized in Tab. I. The stacks were designed with extra contacts to the middle superconducting electrode  so that the voltages across each LJJ can be measured separately. The inner diameter of all stacks was $`D=122.5\mu \mathrm{m}`$ and the width $`W=10\mu \mathrm{m}`$. Due to technological difficulties of making a stack of identical LJJ’s with contacts to the middle electrode, the two stacked junctions had rather substantial difference in quasiparticle (subgap) resistance $`R_{\mathrm{QP}}`$. The normalized circumference of the ring was $`\pi D/\lambda _J=L/\lambda _J15`$, where $`\lambda _J`$ is the Josephson penetration depth, which was approximately equal in both junctions. Measurements were performed in the temperature range $`4.2`$$`5.8\mathrm{K}`$.

In stacked annular LJJ’s, clean trapping of a single fluxon in a desired junction is rather difficult due to the asymmetry of the required state $`[1|0]`$. In the particular case of the 3 samples mentioned above, the asymmetry in the junction’s resistance allowed to trap the fluxon in the desired $`[1|0]`$ state without many efforts just by applying a small bias current through one of the junctions during cooling the sample below the critical temperature $`T_c`$. After every trapping attempt, the resulting state was checked. The $`I`$$`V`$ characteristic (IVC) of both LJJ’s were traced simultaneously in such a way that the current was applied through the whole structure (through two junctions connected in series) and the voltages were measured individually across each LJJ. The wanted state $`[1|0]`$ with a fluxon in one junction and no fluxon in the other junction was identified by a simultaneous observation of a small critical current $`I_c`$ and fluxon step with the smallest asymptotic voltage $`20\mu \mathrm{V}`$ in LJJ<sup>A</sup>, and a large critical current in LJJ<sup>B</sup>. Both the current amplitude of the fluxon step $`I_{\mathrm{max}}^A(H)`$ and the critical current $`I_c^B(H)`$ are expected to have their maxima at zero applied magnetic field $`H=0`$. To check that we have clean fluxon trapping, i.e. that the fluxon is trapped in a LJJ and not accompanied by the parasitic Abrikosov vortices in the superconducting films surrounding LJJ, we checked the dependences $`I_c^{A,B}(H)`$ and $`I_{\mathrm{max}}(H)`$ after each trapping attempt and repeated it until these dependences were symmetric.

The main experimental result of the paper is shown in Fig. 2. It is IVC’s of both LJJ’s of the sample #2 traced at $`T5\mathrm{K}`$ using rather complex current sweep sequence. Note, that the voltage scales of two IVC’s in Fig. 2 are different and shown on the bottom axis for $`V^A`$ and on the top axis for $`V^B`$. The sweep starts at the bias point A where $`I=0`$ and $`V=0`$ and a fluxon is trapped in LJJ<sup>A</sup> (state $`[1|0]`$). When the current is increased up to about $`I=0.69\mathrm{mA}`$ (point B in Fig. 2), LJJ<sup>A</sup> switches to the fluxon step, while LJJ<sup>B</sup> still remains in the Meissner state. Ideally, the LJJ<sup>A</sup> should switch to the fluxon step at zero bias current since any non zero current applied should drive the fluxon around the stack. In our case the fluxon is pinned, most probably near one of the contacts to the middle electrode, and only current $`I=0.69\mathrm{mA}`$ can tear it away from the pinning center. With the further increase of the bias current the LJJ<sup>A</sup> follows the fluxon step which corresponds to the fluxon rotating in the ring, and the voltage across LJJ<sup>A</sup> is proportional to the fluxon rotation frequency according to the Josephson relation.

In an ideal single annular LJJ, the fluxon step has a relativistic nature and its slope approaches infinity when the fluxon moves with velocity $`u`$ close to the Swihart velocity $`\overline{c}_0`$. In the stack with different $`j_c`$ or with inhomogeneities, the fluxon’s velocity can exceed the Swihart velocity. This results in the emission of the electromagnetic waves travelling behind moving fluxon (Cherenkov radiation) and in a finite slope of the step at any velocity. If the length of the emitted radiation tail is comparable with the circumference of the LJJ, the resonant structures (small steps) can appear on the top of the fluxon step. Such steps are visible on the top of the fluxon step in Fig. 2 and are outlined by the circle. At $`I2.67\mathrm{mA}`$ (point C in Fig. 2) both junctions simultaneously switch to the resistive state (gap voltage) with rapidly whirling Josephson phase. Such a simultaneous switching is called current locking. We studied it in detail for stacks of linear geometry in Ref. . Analyzing the dependence of critical current $`I_c^B(H)`$ and maximum current of the fluxon step $`I_{\mathrm{max}}^A(H)`$ on magnetic field, we conclude that the current locking was driven (initiated) by LJJ<sup>A</sup> (if LJJ<sup>A</sup> is kept in the resistive state, $`I_c^B`$ is substantially higher).

When both LJJ’s are in the resistive state, we inverse the direction of the sweep, i.e. start reducing the bias current. At $`I=1.08\mathrm{mA}`$ (point D in Fig. 2), LJJ<sup>A</sup> switches from the resistive state to the fluxon step while LJJ<sup>B</sup> still stays in the phase-whirling state. We denote such state of the stack as $`[1|R]`$. The voltage $`V^A`$ in the $`[1|R]`$ is by about $`16\%`$ smaller than $`V^A`$ in the $`[1|0]`$ state at the same bias. In fact, this is one of key observations in our study.

At this point there are two possibilities: first, continue to decrease the bias current down to zero or, second, increase the current and trace up the single fluxon step for $`[1|R]`$ state.

If we continue decreasing the bias current, at $`I=0.916\mathrm{mA}`$ (point E in Fig. 2) LJJ<sup>B</sup> switches from the resistive state (McCumber branch) to the Meissner state and the overall state of the stack becomes $`[1|0]`$. This causes the voltage $`V^A`$ to increase and become equal to voltage of the fluxon step which we traced in the beginning of the bias current sweep. The fact that switching of LJJ<sup>B</sup> caused a voltage jump across LJJ<sup>A</sup> is marked in Fig. 2 by dotted arrow. Thus, we demonstrated experimentally that the change in the state of LJJ<sup>B</sup> affects the velocity of fluxon moving in LJJ<sup>A</sup>. Further decrease of the bias current results in the fluxon pinning at $`I=0.470\mathrm{mA}`$ (point F in Fig. 2) and in the zero voltage across both LJJ’s.

The second possibility is, being in the bias point D, to increase bias current and trace the fluxon step of the $`[1|R]`$ state up. This step has rather peculiar shape as shown in Fig. 2. In addition to the common trend to have smaller voltage than the fluxon step in the $`[1|0]`$, the step in $`[1|R]`$ state bends back and forth which possibly implies some interesting physics behind it. We also noticed that as we trace both IVC’s in $`[R|R]`$ state from bias point C down to bias point D and then in $`[1|R]`$ state from point D up to bias point G, the voltage across the LJJ<sup>B</sup> has a small hysteresis at voltages equal to the sum of the gap voltages of the superconducting electrodes constituiting the LJJ<sup>B</sup>. This small hysteresis is shown magnified in the inset of Fig. 2. The voltage across LJJ<sup>B</sup> is somewhat smaller in the $`[1|R]`$ state than in the $`[R|R]`$ state. As soon as LJJ<sup>A</sup> switches to the resistive state (point G in Fig. 2 and dotted arrow in the inset), this difference vanishes.

We propose the following explanation for the observed back bending. As we increase current starting from point D up to $`I1.5\mathrm{mA}`$, corresponding to nearly vertical slope of the fluxon step, the voltage $`V^B`$ increases from $`55\mu \mathrm{V}`$ up to $`1.9\mathrm{mV}`$, i.e. approaches the gap voltage. The fluxon motion in LJJ<sup>A</sup>, due to the coupling between the junctions, causes oscillations of Josephson phase and, therefore, of electric and magnetic fields in LJJ<sup>B</sup>. This leads to photon-assisted tunneling (PAT) effect  in LJJ<sup>B</sup>. This effect was earlier observed by Giaever  also using a stack of two junctions. The characteristic frequency $`\omega `$ of photons absorbed in LJJ<sup>B</sup> is equal to the fluxon rotation frequency in LJJ<sup>A</sup>. Due to PAT, one expects to observe a step at the gap sum voltage decreased by $`\mathrm{}\omega /e=2V^A`$. In the low bias region shown in the inset of Fig. 2 the gap sum step is not well defined which results in somewhat weaker gap suppression result. In fact, the maximum suppression of the voltage $`V^B`$ we have found is about $`20\mu \mathrm{V}`$ that corresponds to the top of the back bending region at $`I2.6\mathrm{mA}`$. Since the Josephson voltage in LJJ<sup>A</sup> is also about $`20\mu \mathrm{V}`$, its effect on $`V^B`$ is only about 50% of the expected PAT step voltage change. Thus, as we increase current, the gap voltage in LJJ<sup>B</sup> decreases due to photon-assisted tunneling. The resulting decrease of fluxon step voltage of LJJ<sup>A</sup> is associated with the appearance of an additional dissipation channel due to PAT. Since the PAT step on IVC is limited in voltage (by $`2V^A`$) and in current amplitude ($``$ to the amplitude of the first Josephson harmonic which, in the case of fluxon motion, saturates at some bias), the bending to the right caused by Cherenkov radiation appears to be stronger at $`I>2.6\mathrm{mA}`$ so that the fluxon step of LJJ<sup>A</sup> gains the positive slope again. The differential resistance at the top the fluxon step in $`[1|R]`$ state is rather high and no resonances are observed. This picture is typical for fluxon with Cherenkov radiation tail moving in a media with high dissipation.

The negative bias part of the IVC reproduces all the features described above for the positive half except for the small hysteresis between bias point B and F. Very similar IVC’s were found for other 2 measured samples. The minor difference was in the particular values of the bias current in the points B, C, D, E, F and G which were also dependent on $`T`$. All samples showed that the voltage $`V^A`$ of the fluxon step in the $`[1|0]`$ state is somewhat higher than the voltage of the same step at the same bias in the $`[1|R]`$ state. The fluxon step in the $`[1|R]`$ state showed back and forth bending for all samples. The hysteresis between points B and F was intersecting with the hysteresis between points D and E for some samples and temperatures, so we had to use even more complex sweep sequence in order to trace out all possible dynamical states.

## III Theory

The main objective of this section is to analyze the origin of the decrease in the fluxon velocity in one LJJ due to the switching of the neighboring junction into the resistive state. Here we use the standard RSJ model which does not take into account the dependence of the dissipation on voltage at $`VV_g2.4\mathrm{mV}`$. Thus the gap related effects like PAT discussed above are neglected.

The fluxon dynamics in the system under investigation can be described in the framework of the inductive coupling model which for the case of two coupled junctions takes the form:

$`{\displaystyle \frac{\varphi _{xx}}{1S^2}}\varphi _{tt}\mathrm{sin}\varphi {\displaystyle \frac{S\psi _{xx}}{1S^2}}`$ $`=`$ $`\alpha \varphi _t\gamma ;`$ (1)   
$`{\displaystyle \frac{\psi _{xx}}{1S^2}}\psi _{tt}\mathrm{sin}\psi {\displaystyle \frac{S\varphi _{xx}}{1S^2}}`$ $`=`$ $`\alpha \psi _t\gamma ,`$ (2)   
where $`\varphi `$ and $`\psi `$ are the Josephson phases across the LJJ<sup>A,B</sup>, respectively, $`1<S<0`$ is a dimensionless coupling parameter, $`\alpha 1`$ is the damping coefficient describing the dissipation in the system due to quasiparticle tunneling, and $`\gamma =j/j_c`$ is the normalized the density of the bias current flowing through the stack. The coordinate $`x`$ and time $`t`$ are measured, respectively, in units of the Josephson length $`\lambda _J`$ and inverse plasma frequency $`\omega _p^1`$ of single-layer LJJ. Most of relevant parameters of the junctions, such as effective magnetic thicknesses and specific capacitances, for the sake of simplicity are assumed to be equal in both LJJ’s.

To understand the origin of an additional friction force we start from unperturbed (without r.h.s.) Eqs. (1) and (2) and use the force balance equations to derive the shape of the IVC. We do not directly solve Eqs. (1) and (2), but, rather, use trial solutions for $`\varphi (x,t)`$ and $`\psi (x,t)`$, which approximate the Josephson phase profiles in the states $`[1|0]`$ or $`[1|R]`$. The choice of the trial functions is suggested by the results of numerical simulations presented in the following section. To simplify the mathematics and concentrate attention on the physical sense, the dense fluxon chain approximation is used. For this case we adopt the following trial solutions:

$`\varphi (x,t)=H(xut)+A_{r,m}\mathrm{sin}\left[H(xut)\right];`$ (3)   
$`\psi ^m(x,t)=B_m\mathrm{sin}\left[H(xut)\right];`$ (4)   
$`\psi ^r(x,t)=\omega t{\displaystyle \frac{1}{\omega ^2}}\mathrm{sin}(\omega t)+B_r\mathrm{sin}\left[H(xut)\right],`$ (5)   
where $`\psi ^m`$ and $`\psi ^r`$ are phases in LJJ<sup>B</sup> in the Meissner state and resistive state, accordingly (these are the two cases which we are going to compare); $`u`$ is the velocity of the fluxon chain; $`A_{r,m}`$, $`B_{r,m}1`$ are the constants which we are going to determine for resistive and Meissner states, respectively; $`H`$ is the average normalized magnetic field in the LJJ<sup>A</sup>:

$$H=\frac{2\pi N}{\mathrm{}},$$
(6)   
where $`N`$ is the number of fluxons trapped in annular LJJ<sup>A</sup> (for $`[1|0]`$ and $`[1|R]`$ states $`N=1`$), and $`\mathrm{}=L/\lambda _J`$ is the normalized length of the junctions. Dense fluxon chain approximation implies that $`H1`$.

Substituting Eqs. (3)–(5) into Eqs. (1) and (2), and using the following approximations (which are justified in our case)

$`\mathrm{sin}\varphi \mathrm{sin}\left[H(xut)\right];`$ (7)   
$`\mathrm{sin}\psi ^mB\mathrm{sin}\left[H(xut)\right];`$ (8)   
$`\mathrm{sin}\psi ^r\mathrm{sin}(\omega t),`$ (9)   
we arrive to the equations from which we can determine $`A_{r,m}`$ and $`B_{r,m}`$. The final result for the $`[N|R]`$ (resistive) state is

$`A_r`$ $`=`$ $`{\displaystyle \frac{D}{H^2Q}};`$ (10)   
$`B_r`$ $`=`$ $`{\displaystyle \frac{S}{H^2Q}},`$ (11)   
where we introduced notations

$`D`$ $`=`$ $`1u^2(1S^2);`$ (12)   
$`Q`$ $`=`$ $`\left(1u^2\right)^2u^4S^2.`$ (13)   
Note, that both $`D>0`$ and $`Q>0`$ for $`u<\overline{c}_{}`$.

The result for the $`[N|0]`$ (Meissner) state is

$`A_m`$ $`=`$ $`{\displaystyle \frac{DH^2+1S^2}{H^2\left(QH^2+D\right)}};`$ (14)   
$`B_m`$ $`=`$ $`{\displaystyle \frac{S}{QH^2+D}}.`$ (15)   
To calculate IVC, $`\gamma (u)`$, we write the force balance equation

$$2\pi N\gamma =F_\alpha ^A+F_\alpha ^B.$$
(16)   
Here $`F_\alpha ^{A,B}`$ are the friction forces which develop in LJJ<sup>A,B</sup>. The expression for the friction force is well known from the perturbation theory 

$$F_\alpha =\alpha _0^{\mathrm{}}\varphi _x\varphi _t𝑑x,$$
(17)   
where we will use $`\mathrm{}=2\pi N/H`$ following from (6). Since we are interested in the average friction force to get an IVC, we have to perform friction force averaging in time:

$$\overline{F}_\alpha =\frac{1}{T}_0^TF_\alpha (t)𝑑t.$$
(18)   
Since there are two characteristic frequencies in the system, fluxon (Josephson) frequency and the phase-whirling frequency of the resistive state, we have to choose the averaging interval $`T`$ in Eq. (18) so that it will contain an integer number of periods of each frequency i.e. $`T=2\pi k/\omega =2\pi m/Hu`$, and $`k`$, $`m`$ being the integer constants. After averaging, we get the following expressions for friction forces

$`\overline{F}_\alpha ^A`$ $`=`$ $`\pi N\alpha Hu\left(A_{r,m}^2+2\right);`$ (19)   
$`\overline{F}_\alpha ^B`$ $`=`$ $`\pi N\alpha HuB_{r,m}^2.`$ (20)   
All information about the actual state is contained in $`A_{r,m}`$ and $`B_{r,m}`$ calculated above for the Meissner and resistive state of the LJJ<sup>B</sup>.

Finally, we insert (19) and (20) into the force balance equation (16) and get IVC’s

$$\gamma _{r,m}(u)=\frac{\alpha Hu}{2}\left(A_{r,m}^2+B_{r,m}^2+2\right).$$
(21)   
Now we can prove that IVC for the $`[N|R]`$ state is shifted to the region of lower velocities in comparison with the IVC for $`[N|0]`$ state, i.e. that

$$\delta (u)=\gamma _r(u)\gamma _m(u)>0\text{ for all }|u|<\overline{c}_{}.$$
(22)   
Substituting $`A_{r,m}`$ and $`B_{r,m}`$ from Eqs. (10), (11), (14), (15) into the expression (21) and using the obtained expressions for $`\gamma _{r,m}(u)`$ in (22) we get

$$\delta (u)=\frac{\alpha u\left(X_1H^2+X_2\right)}{2H^3Q^2\left(H^2Q+D\right)^2},$$
(23)   
where $`X_1`$ and $`X_2`$ are defined as

$`X_1`$ $`=`$ $`2QD\left[S^2Q\left(1S^2\right)+D^2\right];`$ (24)   
$`X_2`$ $`=`$ $`D^2\left(D^2+S^2\right)\left(1S^2\right)^2Q^2,`$ (25)   
Obviously, Eq. (23) is positive when both $`X_1`$ and $`X_2`$ are positive. To prove the latter, we express $`D^2`$ as a function of $`Q`$ using Eqs. (12) and (13)

$$D^2=Q\left(1S^2\right)+S^2.$$
(26)   
Substituting (26) into (24) and (25) we get

$`X_1`$ $`=`$ $`4QDS^2>0;`$ (27)   
$`X_2`$ $`=`$ $`S^2\left[3Q\left(1S^2\right)+2S^2\right]>0.`$ (28)   
Thus (22) is proved and $`\gamma _r(u)>\gamma _m(u)`$ for any $`u<\overline{c}_{}`$.

This result is in agreement with our experiment and simulation (see the following section). From physical point of view, the origin of the effect lays in the difference between Eqs. (8) and (9) where the main term depends on the state of LJJ<sup>B</sup> resulting in different phase profile and different friction force for $`[1|0]`$ and $`[1|R]`$ states. The IVC’s of fluxon steps $`u^A(\gamma )`$ for the states $`[1|0]`$ and $`[1|R]`$ calculated using Eq. (21) are shown in Fig. 3. According to the calculations presented above the difference $`\delta (u)`$ diverges as $`u`$ approaches $`\overline{c}_{}`$. Actually, for this case of a single fluxon in our relatively long junction, the dense fluxon chain approximation is not fully valid so one needs to perform more exact analysis or numerical simulations. The region of validity of our approximation is $`A_r<1`$ and $`B_r<1`$. Using Eqs. (10) and (11) for the same parameters as in Fig. 3 we get that our approximation is valid up to $`u0.809`$ while $`\overline{c}_{}0.816`$.

## IV Simulation

To check the limitations of the analysis presented above we performed direct numerical simulations. Our simulations show that the effect observed in experiment and explained in the framework of high fluxon density approximation exists for any fluxon densities and any velocity even very close to $`\overline{c}_{}`$. Another advantage of the simulation over the analytical approach is that the simulation fully reproduce the dynamics of the inductive coupling model and, therefore, all the effects possible in its framework.

The numerical procedure works as follows. For a given set of LJJ’s parameters we simulate the IVC of the system, i.e. calculate $`\overline{V}^A(\gamma )`$ and $`\overline{V}^B(\gamma )`$ while increasing $`\gamma `$ from zero up to $`1`$. To calculate the voltages $`\overline{V}^A(\gamma )`$ and $`\overline{V}^B(\gamma )`$ for each value of $`\gamma `$, we simulate the dynamics of the phases $`\varphi ^{A,B}(x,t)`$ by solving the Eqs. (1) and (2) with the periodic boundary conditions:

$`\varphi ^{A,B}(0,\stackrel{~}{t})`$ $`=`$ $`\varphi ^{A,B}(\mathrm{},\stackrel{~}{t})+2\pi N^{A,B}`$ (29)   
$`\varphi _{\stackrel{~}{x}}^{A,B}(0,\stackrel{~}{t})`$ $`=`$ $`\varphi _{\stackrel{~}{x}}^{A,B}(\mathrm{},\stackrel{~}{t}),`$ (30)   
numerically using an explicit method \[expressing $`\varphi ^{A,B}(t+\mathrm{\Delta }t)`$ as a function of $`\varphi ^{A,B}(t)`$ and $`\varphi ^{A,B}(t\mathrm{\Delta }t)`$\], and treating $`\varphi _{xx}`$ with a five point, $`\varphi _{tt}`$ and $`\varphi _t`$ with a three point symmetric finite difference scheme. Numerical stability was checked by doubling the spatial and temporal discretization steps $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }t`$ and checking its influence on the fluxon profiles and on the IVC. The discretization values used for simulation were $`\mathrm{\Delta }x=0.01`$, $`\mathrm{\Delta }t=0.0025`$. After simulation of the phase dynamics for $`T_0=20`$ time units we calculate the average dc voltages $`\overline{V}^{A,B}`$ during this time interval as

$$\overline{V}^{A,B}=\frac{1}{T}_0^T\varphi _t^{A,B}(t)𝑑t=\frac{\varphi ^{A,B}(T)\varphi ^{A,B}(0)}{T}.$$
(31)   
For faster convergence, we use the fact that $`\overline{V}^{A,B}`$ does not depend on $`x`$ and, therefore, we also take advantage of the spacial averaging of the phases $`\varphi ^{A,B}`$ in (31).

When the values of $`\overline{V}^{A,B}`$ are found from (31), the dynamics of the phases $`\varphi ^{A,B}(x,t)`$ is simulated further during $`1.2T_0`$ time units, the dc voltages $`\overline{V}^{A,B}`$ are calculated for this new time interval and are compared with the previously calculated values. We repeat such iterations further increasing the time interval by a factor 1.2 until the difference in dc voltages $`|\overline{V}(1.2^{n+1}T)\overline{V}(1.2^nT)|`$ obtained in two subsequent iterations becomes less than a given accuracy $`\delta V=10^3`$. The particular value of the factor $`1.2`$ was found to be quite optimal to provide fast convergence as well as more effective averaging of low harmonics on subsequent steps. Very small value of this factor, e.g. $`1.01`$ can result in very slow convergence in the case when $`\varphi (t)`$ contains harmonics with the period comparable or larger than $`T`$. Large values of the factor, e.g. 2 or higher, will consume a lot of CPU time already during the second or third iteration even when the convergence is good. After the voltage averaging for current $`\gamma `$ is complete, $`\gamma `$ is changed by a small amount $`\delta \gamma `$ to calculate the voltages in the next point of the IVC. As initial conditions here we use a distribution of phases (and their derivatives) achieved in the previous point of the IVC.

An example of calculated IVC is shown in Fig. 4. To trace both the Meissner and resistive states we use the following sweep sequence: $`\gamma `$ increases from 0 up to 1 with a step $`\delta \gamma =0.01`$, then decreases down to 0.5 with a step $`\delta \gamma =0.01`$, and further down with a step $`\delta \gamma =0.002`$ until the state $`[N|R]`$ is reached as shown in Fig. 4. From this point, we either sweep further down to $`\gamma =0`$ or up to $`\gamma =1`$ until both junctions switch to the resistive state. The decrease of voltage in $`[1|R]`$ state in comparison with $`[1|0]`$ state is very clearly seen in Fig. 4. We performed this kind of simulations for the wide range of the parameters i.e. for $`|S|=0.1`$, $`0.2\mathrm{}0.5`$, and $`N=1`$, $`2\mathrm{}5`$ (in total 25 pairs of IVC’s) and found similar IVC’s in all cases. For relatively dense fluxon chain $`N/L=1`$ we also compared the IVC’s obtained by means of numerical simulation with IVC’s derived analytically and found a good agreement as shown in Fig. 3. Small difference in the slope of analytical and numerical IVC’s is related to final density of fluxons in simulation $`H=2\pi `$ while theoretical curve corresponds to $`H1`$.

The limitations of our model due to voltage independent loss term $`\alpha `$ prevented proper calculation of the upper part of the fluxon step in $`[1|R]`$ state. Here numerical curve shows a series of small voltage jumps (see Fig. 4) in LJJ<sup>A</sup> that are not observed in experiment. These jumps are related to excitation of fluxon-antifluxon states in LJJ<sup>B</sup> at high voltages. Such states were not found in experiment, most probably due to the increased dumping in LJJ<sup>B</sup> at gap voltage.

## V Conclusion

We investigated experimentally and theoretically the motion of a fluxon in one of two magnetically coupled long Josephson junctions. Two different cases are studied: (1) when the neighboring junction (one which does not contain any fluxon) is in the Meissner state, and (2) when the neighboring junction is in a phase-whirling (resistive) state. We found that the phase-whirling state in LJJ<sup>B</sup> slows down the fluxon motion in LJJ<sup>A</sup> and results in a shift of the fluxon step to lower voltage. This effects is detected experimentally, reproduced in simulations based on the inductive coupling model, and derived analytically in the high fluxon density approximation. In addition, experiment shows quite peculiar back and forth bending of the fluxon step in $`[1|R]`$ state which we explain as a result of increased dumping due to photon assisted tunneling effect in LJJ<sup>B</sup>. The results of our study are also relevant for characterization of stacked Josephson junctions with large number of layers.

###### Acknowledgements.

We thank Norbert Thyssen for the sample fabrication. Partial support of this work by Deutsche Forschungsgemindschaft (DFG) is also acknowledged.