---
address: 'Max Planck Institute for the Physics of Complex Systems, D-01187, Dresden, Germany'
author:
- 'M. V. Fistul and J. B. Page[@byline]'
title: |
    Penetration of dynamic localized states in DC-driven\
    Josephson junction ladders by discrete jumps
---

Introduction
============

The subject of large-amplitude anharmonic dynamics in lattices has received widespread attention over the past decade. In particular, intense theoretical focus has centered on so-called intrinsic localized modes (ILMs), also known as discrete breathers, with the result that many of their properties are now well understood[@ILMs]. These excitations result from the interplay between nonlinearity and discreteness, and they can be highly localized in perfect lattices, with or without external driving. They can occur in a variety of different lattices: recent experiments have reported vibrational ILMs in a quasi-1D charge-density wave system[@Bishop], spin-wave ILMs in a quasi-1D antiferromagnetic system[@Sievers], and discrete breathers in Josephson junction (JJ) ladders[@Ustinov1; @Orlando].

The latter systems are noteworthy, in that arrays of coupled JJs have served for many years as reliable laboratory systems for studying diverse nonlinear phenomena[@StrogLikh]. The nonlinear dynamics are particularly rich. A single “small” JJ subject to an applied constant DC bias current can be mapped onto the problem of a damped pendulum driven by a constant torque, with the dynamical degree of freedom being the Josephson phase difference[@BandC]. There are thus two qualitatively different states, namely a static (superconducting) state and a dynamic (whirling or resistive) state, with the latter producing a readily measured voltage $V \propto \dot \varphi$ across the whirling junction. When several junctions are assembled to form a regular array, such as the ladder shown in Fig. \[Fig1\], they become inductively coupled. In the coupled system, junctions in the superconducting state can also exhibit steady-state librations, when JJs in the whirling state are present. In view of the mapping onto the pendulum problem, JJ ladders share features with lattices of nonlinearly coupled electric dipole rotors, driven by an external monochromatic AC electric field[@BonPage], but with the important simplification that they can be driven with purely DC bias currents.

Figure \[Fig1\] sketches an anisotropic JJ ladder consisting of small JJs of two types, “horizontal” and “vertical,” which are, respectively, perpendicular and parallel to the applied bias current (arrows). The anisotropy arises from the different areas of the horizontal and vertical junctions and is characterized by the parameter $\eta=A_h/A_v=I_c^h/I_c^v$, where $I_c^h$ and $I_c^v$ are the the critical currents for each type of junction.

References [@Ustinov1; @Orlando; @Ustinov2] report experimental observations of various discrete breathers in ladders driven by a [*homogeneously*]{} applied DC bias current, represented by the dashed arrows in Fig. \[Fig1\]. For these states, the localized voltage patterns have a simple structure involving only two nonzero steady-state voltages (rotational frequencies). The breathers were found to be stable in the limit of small coupling ($\eta \alt 1$) and for bias currents $I_{\rm ext} \alt I_c^v$. For the case of large 2-D JJ arrays subject to a homogeneously applied DC bias, more complicated inhomogeneous states, with [*meandering*]{} voltage patterns, have also been reported[@Misha1].

Here, we study the dynamics of a JJ ladder with an external DC bias current applied at only one edge (solid arrows, Fig. \[Fig1\]). For increasing bias ($I_{\rm ext} \agt I_c^v$) and over a wide range of anisotropies, we find by direct numerical simulations that the dynamic state expands into the ladder one cell at a time, by a sequence of abrupt jumps. This behavior is in marked contrast to the well-known cases of long JJs and JJ parallel arrays ($\eta~=~\infty$), where the entire system abruptly switches to the resistive state at a particular value of the DC bias. It is also different than the breather case, since all of the junctions within the boundary of this localized dynamic state whirl, and at each expansion the number of different frequencies (voltages) grows. The sequence of $I$-$V$ characteristics and threshold currents can be modeled analytically, yielding very good agreement with the numerical results.

Numerical Simulations
=====================

We consider a ladder with a large but finite number of cells $N$. The ladder’s state is specified by the time-dependent Josephson phases $\{\varphi_i\}$, $\{\psi_i\}$, and $\{\tilde \psi_i\}$ for the vertical, upper horizontal and lower horizontal junctions, respectively, where $i$ denotes the cell. We have found in our simulations that the symmetry condition ${\tilde \psi_i}=-\psi_i$ holds for the phenomena to be discussed here. The ladder dynamics are then determined by the coupled nonlinear equations of motion obtained in Refs.  and : $$\begin{aligned}
\label{GenEq}
\hat L(\varphi_i)&=&\gamma_i+\frac{1}{\beta_L}[\varphi_{i-1}-2\varphi_i+
\varphi_{i+1}+2(\psi_i-\psi_{i-1})], \nonumber \\
&&i=2,\ldots,N-1, \\
\hat L(\psi_i)&=&\frac{1}{\eta \beta_L}(\varphi_i-\varphi_{i+1}-2\psi_i), \qquad
\quad i=1,\ldots,N, \nonumber\end{aligned}$$ where the operator $\hat L(\varphi) \equiv \ddot \varphi + \alpha \dot\varphi
+\sin(\varphi)$. The equations for the vertical junctions at $i=1$ and $N$ are $$\begin{aligned}
\label{BC1}
\hat L(\varphi_1)&=&\gamma_1+\frac{1}{\beta_L}(\varphi_2-\varphi_1+2\psi_1), \\
\hat L(\varphi_N)&=&\gamma_N+
\frac{1}{\beta_L}(\varphi_{N-1}-\varphi_N-2\psi_{N-1}). \nonumber\end{aligned}$$

Equations (\[GenEq\]) and (\[BC1\]) describe each junction within the resistively and capacitively shunted junction (RCSJ) model[@BandC], and the unit of time is the inverse of the plasma frequency $\omega_J \equiv \sqrt{2eI_c/C\hbar}$. Since each junction’s critical current and capacitance scale with the area, $\omega_J$ is independent of the anisotropy parameter $\eta$, as is the effective damping constant $\alpha \equiv 1/(\omega_J RC)$. The normalized bias current $\gamma_i$ is defined as $I_{i,{\rm ext}}/I_c^v$. The inductive coupling between junctions is determined by the parameter $\beta_L \equiv 2\pi L I_c^v/\Phi_0$, where $L$ is the self-inductance of a single cell and $\Phi_0 = hc/2e$ is the elementary flux quantum. Coupling beyond that described by $\beta_L$ is not included.

We performed direct numerical integration of the equations of motion for ladders with $N=20$ cells, using a fifth-order Gear predictor-corrector algorithm[@AandT], for a range of anisotropies: $\eta =$ 0.5, 1.0, 2.0, 3.0, and 5.0. The arrays were underdamped, with $\alpha=0.1$, and we used a moderate value of the coupling parameter $\beta_L=0.5$. The external DC bias was applied at one edge, i.e. $\gamma_1=\gamma$ and all other $\gamma_i=0$. To simulate the $I$-$V$ curves, we started with all phases at zero and gradually increased the external bias $\gamma$ from zero to 50, in increments of 0.005. When junctions were present in the whirling state, the MD time scale was set by the time-average period of the fastest rotating phase. For a given value of gamma, we waited for at least 100 of these reference periods before averaging, in order to avoid transients, following which we computed the time-averages $\langle \dot \varphi_i \rangle$ and $\langle \dot \psi_i\rangle$ over at least 100 additional reference periods. These averages are proportional to the average voltages across the junctions. The current was then incremented, with the initial phase configuration being that from the preceding MD time step. In all runs, the time step was 1/200 of the reference period.

Our simulated $I$-$V_i$ curves for an anisotropy of $\eta =2$ are shown in Fig. \[Fig2\]. The most striking finding is the occurrence of extremely sharp voltage jumps. At each of the corresponding threshold currents $\gamma^{thr}_n$, a new cell is added to the ladder’s dynamic state. With the applied current below the first threshold $\gamma^{thr}_1$, all junctions are in a static (zero voltage) state. When $\gamma$ exceeds $\gamma^{thr}_1$, the first vertical junction and its adjacent top and bottom horizontal junctions abruptly switch into the dynamic state, with all other junctions remaining in the zero voltage state—the rotating phases are confined to the first cell. As the bias is increased further, all three average voltages for this 1-cell dynamic state increase linearly until the next threshold current $\gamma^{thr}_2$ is reached, at which point the dynamic state suddenly expands into the second cell, accompanied by sharp changes of the voltages in the first cell. This process continues, yielding successive transitions from $n$-cell dynamic localized states to $(n+1)$-cell dynamic states. The distribution of threshold currents and voltage ratios depends on the ladder’s anisotropy. Over the range $0<\gamma <50$, the ladder reached a 3-cell state for $\eta=0.5$ and 1, a 4-cell state for $\eta=2$, and a 5-cell state for $\eta=3$ and 5. In the following, these states will be termed $n$-cell edge states.

An $n$-cell edge state is in striking contrast to an n-cell breather. The breather occurs away from the ladder’s edge and is homogeneously driven by same DC bias current ($I_{\rm ext} \alt I_c$) applied at every cell, whereas the edge states are driven by a DC bias ($I_{\rm ext} \agt I_c$) applied at just one edge. The edge states have a richer internal structure—[*all*]{} of the junctions within an edge state are in a nonzero voltage state (see Fig. \[Fig3\]), whereas in a breather state, all of the horizontal junctions are in the zero voltage state except for those on the breather’s boundary[@Ustinov1; @Orlando; @Ustinov2]. Moreover, all of a breather’s vertical junction phases rotate at the same average frequency, whereas the $n$-cell edge state exhibits a peculiar distribution of average frequencies. This frequency (voltage) distribution depends on both $n$ and the ladder’s anisotropy. For example, our simulations for the $\eta=2$, 3-cell edge state of Fig. \[Fig3\](d) yield the ratios given in second column of Table \[Table1\]. Comparison with the third column shows that they are in very good agreement with analytic predictions derived below.

  ------------------------- ------- ---------------------------------------
            Ratio             MD                   Predicted
   $\omega_1^v/\omega_2^v$   2.667   $(3\eta^2+8\eta+4)/2\eta(1+\eta)=8/3$
   $\omega_1^v/\omega_3^v$   8.006       $(3\eta^2+8\eta+4)/\eta^2=8$
   $\omega_1^h/\omega_2^h$   2.499    $(\eta^2+6\eta+4)/\eta(2+\eta)=5/2$
   $\omega_1^h/\omega_3^h$   5.017        $(\eta^2+6\eta+4)/\eta^2=5$
   $\omega_3^v/\omega_3^h$   2.005                     2
  ------------------------- ------- ---------------------------------------

  : Average frequency (voltage) ratios for 3-cell edge states. The MD ratios are for the $\eta=2$, 3-cell edge state of Fig. \[Fig3\](d), and the predicted ratios were calculated from Eqs. (\[Vert\]) and (\[Horiz\]).

\[Table1\]

The superconducting state forming ahead of the $n$-cell edge state is also unusual. Fig. \[Fig3\] gives snapshot images of the Josephson phase distribution for several values of the applied DC current bias, for the anisotropy $\eta=2$. In panel (a), the current is just below the first threshold, and one sees a single [*Josephson vortex*]{} in the superconducting part of the ladder. The remaining panels (b)–(e) show the phases just after a new cell is added to the dynamic state. At the threshold currents $\gamma^{thr}_n$, the superconducting state becomes unstable, and the vortex jumps into the next cell as the edge state expands. Our simulations show that in general the superconducting state is sensitive to the anisotropy. Thus for rather small values of $\eta \lesssim 1$, there are no vortices trapped in the superconducting part of the ladder over the range of currents studied. For these cases, the Josephson phases of the vertical junctions in the superconducting portion of the ladder simply decrease with distance from the boundary of the resistive portion, corresponding to the [*Meissner state*]{} of the superconductor. With increasing anisotropy, single Josephson vortices appear in the superconducting portion, as in Fig. \[Fig3\]. For large anisotropies ($\eta \sim 5$) more complex [*vortex trains*]{} are observed, and we also find that the penetration of the edge state changes the nature of the superconducting vortex state, rather than simply pushing it ahead as for $\eta=2$. A detailed discussion of the superconducting state will be given elsewhere (Ref. [@FP]).

Theoretical Analysis
====================

The unusual voltage distributions in the $n$-cell edge states can be explained analytically by making use of Kirchhoff’s laws, applied to the time-average currents (normalized to $I^v_c$) and corresponding dimensionless voltages in each cell. The key assumption is the coexistence of the resistive and superconducting states in different portions of the ladder. For a cell $i$ within an $n$-cell edge state, current conservation gives $I^v_i+I^h_i=I^h_{i-1}$, while the voltage condition is $I^v_i-\frac{2}{\eta}I^h_i- I^v_{i+1}=0$. Combining these yields an equation for just the horizontal currents: $$\label{DifEq}
I^h_{i+1}+I^h_{i-1}-2\left(1+\frac{1}{\eta}\right)I^h_i=0.$$ This equation and the two from which it was derived apply to all cells $i$ in $1 \le i \le n$, provided we define $I^h_0 \equiv \gamma,\; 
I^v_{n+1} \equiv 0$, and $I^h_{n+1} \equiv I^h_n$, in order to take proper account of the $n$-cell edge state’s boundaries.

Equation (\[DifEq\]) is readily solved by substituting $I^h_i = \lambda^i$, which yields two roots, namely $\lambda \equiv 
1+\frac{1}{\eta} + \sqrt{(1+\frac{1}{\eta})^2-1}$ and $1/\lambda$. Hence the general solution of Eq. (\[DifEq\]) is $I^h_i = c_1 \lambda^i+ c_2 \lambda^{-i}$, where the constants $c_1$ and $c_2$ are obtained from the above definitions of $I^h_0$ and $I^h_{n+1}$. With the horizontal currents thus determined, the vertical currents may be computed from $I^v_i=I^h_{i-1} - I^h_i$. The currents are then converted into the average voltages via $V^v_i = I^v_i/\alpha$ and $V^h_i=I^h_i/(\alpha \eta)$. The resulting voltage distribution within an $n$-cell edge state is ($1\le i\le n$): $$\label{Vert}
V^v_i=\frac{\gamma(1-\lambda)(\lambda^{i-1} -\lambda^{2n+1-i})}
{\alpha(\lambda^{2n+1}+1)},$$ $$\label{Horiz}
V^h_i=\frac{\gamma (\lambda^{i} +\lambda^{2n+1-i})}
{\alpha \eta (\lambda^{2n+1}+1)}.$$

Equations (\[Vert\]) and (\[Horiz\]) give the predicted voltage ratios in Table \[Table1\], which for $\eta=2$ are seen to be in very good agreement with our MD simulations. Indeed, we find that for all of the values of $\eta$ studied, the predicted $I$-$V$ curves are in very good agreement with the MD curves, such as those of Fig. \[Fig2\]. Only the values of the current thresholds for the jumps are left undetermined by these equations.

We can also predict the distribution $\{\gamma^{thr}_n\}$ of threshold currents for each $\eta$ by assuming that the superconducting state associated with the $(n-1)$-cell edge state becomes unstable and converts to the $n$-cell edge state when the current $I^h_n$ reaches a depinning current $I_{dp}$, which we take to be independent of $n$. This yields an expression for the threshold currents $$\label{Ithr}
\gamma^{thr}_n =
I_{dp}\frac{\cosh{[(n-\frac{1}{2}) \ln \lambda]}}
{\cosh{(\frac{1}{2} \ln \lambda)}}.$$ The ratios $\gamma^{thr}_n/I_{dp}$ predicted by Eq. (\[Ithr\]) for $n=$ 2, 3, 4 and 5 are plotted versus $\eta$ in Fig. \[Fig4\]. To compare with the MD results, we fit $I_{dp}$ to the first observed MD threshold current for each $\eta$, namely $\gamma^{thr}_1 =$ 1.295, 1.510, 2.040, 2.525, and 4.010, for $\eta =$ 0.5, 1.0, 2.0, 3.0, and 5.0, respectively. The resulting MD ratios are shown by the circles in Fig. \[Fig4\] and agree well with the predictions.

Hysteresis
==========

Despite their rich structure of frequency ratios, the $n$-cell edge states are found to be highly stable in our simulations, for the case of increasing current. However, we also find notable hysteresis effects when the simulations are started with an $n$-cell edge state for large $n$ and the applied DC current is gradually [*decreased*]{} to zero. Figure \[Fig5\] is representative of the hysteretic behavior encountered. The threshold currents and stability properties for the sequence of down-conversions $\{n \rightarrow n-1\}$ are very different than for the increasing-current case. In particular, we observe resonant steps, switching processes, and nonlinear regions in the $I$-$V$ curves. We believe that all of these features arise from the resonant interaction between the n-cell edge states and other excitations, both localized and delocalized, as will be discussed elsewhere[@FP].

Summary
=======

In summary, our numerical simulations have revealed unusual localized dynamic states in anisotropic JJ ladders subject to a DC bias current at one edge. Increasing the bias causes these states to expand by adding single cells in a sequence of sudden jumps, giving rise to a diverse set of voltage distributions and sharp changes in the $I$-$V$ curves. This behavior occurs for a wide range of parameters and should be observable through the $I$-$V$ characteristics or by direct visualization using low temperature scanning laser microscopy techniques[@Ustinov1; @Ustinov2; @Misha1].

We thank A. V. Ustinov and S. Flach for useful discussions. J. B. Page gratefully acknowledges the Max Planck Institute for the Physics of Complex Systems, Dresden, for their support and hospitality.

Permanent address: Department of Physics and Astronomy, Arizona State University, Tempe, AZ 85287-1504.

See for instance, S. Flach and C. R. Willis, Phys. Rep. [**295**]{}, 181 (1998); A. J. Sievers and J. B. Page, in [*Dynamical Properties of Solids VII Phonon Physics*]{}, edited by G. K. Horton and A. A. Maradudin (Elsevier, Amsterdam, 1995); and references therein.

B. I. Swanson, J. A. Brozik, S. P. Love, G. F. Strouse, A. P. Shreve, A. R. Bishop, W.-Z. Wang, and M. I. Salkola, Phys. Rev. Lett. [**82**]{}, 3288 (1999).

U. T. Schwarz, L. Q. English, and A. J. Sievers, Phys. Rev. Lett. [**83**]{}, 223 (1999).

P. Binder, D. Abraimov, A. V. Ustinov, S. Flach, and Y. Zolotaryuk, Phys. Rev. Lett. [**84**]{}, 745 (2000).

E. Trias, J. J. Mazo, and T. P. Orlando, Phys. Rev. Lett. [**84**]{}, 741 (2000).

S. H. Strogatz, [*Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering*]{} (Addison-Wesley, Reading, MA, 1994); K .K. Likharev, Rev. Mod. Phys. [**51**]{}, 101 (1979).

A. Barone and G. Paternò, [*Physics and Applications of the Josephson Effect*]{} (Wiley, New York, 1982).

D. Bonart and J. B. Page, Phys. Rev. E [**60**]{}, R1134 (1999).

P. Binder, D. Abraimov, and A. V. Ustinov, Phys. Rev. E [**62**]{}, 2858 (2000).

D. Abraimov, P. Caputo, G. Filatrella, M. V. Fistul, G. Yu. Logvenov, and A. V. Ustinov, Phys. Rev. Lett. [**83**]{}, 5354 (1999).

P. Caputo, M. V. Fistul, A. V. Ustinov, B. A. Malomed, and S. Flach, Phys. Rev. B [**59**]{}, 14050 (1999).

G. Grimaldi, G. Filatrella, S. Pace, and U. Gambardella, Phys. Lett. A [**223**]{}, 463 (1996).

See, for instance, M. P. Allen and D. J. Tildesly, [ *Computer Simulations of Liquids*]{} (Clarendon, Oxford, 1987).

M. Fistul and J. B. Page, unpublished.
