text:
# Discrete Breathers in Josephson Ladders

## 1 Introduction

Linear models of crystals have been instrumental in developing a physical understanding of the solid state. Thermodynamic properties such as specific heat, transport properties such as electron relaxation times or superconductivity, and even interactions with radiation can be understood by modeling a crystal as a lattice of atoms with fixed harmonic coupling. This leads to the conventional phonon-like analysis with a basis of plane waves as normal modes. However, certain properties of solids, such as thermal expansion, cannot be understood in this linear model. For example, the elastic constants of the atomic interactions may depend on the temperature or the volume and so make the interaction non-linear. This is termed an anharmonic effect and the usual approach is to use a more generalized Taylor expansion for the lattice coupling that includes more than just the harmonic term.

Until quite recently anharmonic effects were only studied as perturbations to the fully solvable harmonic model. Then it was discovered  and proved  that in classical Hamiltonian systems non-linearity may lead to localized vibrations in the lattice that cannot be analyzed using the standard plane wave approach. These intrinsic localized modes are time periodic and spatially localized solutions and have been termed discrete breathers (DB). Their amplitudes oscillate around a few sites in the lattice and they do not depend on impurities for their localization. The study of DB has been extended to the dynamics of coupled rotor lattices  where the terms “rotating localized modes” or “rotobreathers” were introduced. DB have also been proven to exist in the dynamics of dissipative systems . In this case chaotic localized solutions have been discovered . DB have been found to play an important role in conditions far from equilibrium  and recently have been studied in disordered systems . Excellent reviews on the topic are  and .

Because DB are generic modes in many non-linear lattices, they are the object of great theoretical and numerical attention in many diverse fields like condensed matter physics , mechanical engineering  and biophysics . Only recently have the first experiments been performed which detect intrinsic localized modes in quasi-one-dimensional charge-density-wave compounds , antiferromagnetic anharmonic crystals  and superconducting arrays .

The existence of DB in a Josephson-junction (JJ) array was first proposed by Floría et al.  in the study of the dynamics of an AC-biased anisotropic ladder array. Both, oscillating and rotating localized modes were simulated and studied in this system . Later, rotobreathers were also numerically studied in the dynamics of inductively coupled junctions in DC-biased arrays . Arrays were then designed, fabricated and measured and DB were found .

The Josephson junctions studied in this article are made from a superconducting-insulator-superconductor (SIS) tunneling structure that because of the Josephson effect behave as solid-state non-linear oscillators. In the framework of the resistively and capacitively shunted junction (RCSJ) model a single JJ is modeled by a parallel combination of an ideal junction, a capacitor $`C`$, and a resistance $`R`$. The current of the SIS tunnel junction is then

$$I=C\frac{dV}{dt}+\frac{V}{R}+I_J.$$
(1)
The ideal junction has a constitutive relation of $`I_J=I_c\mathrm{sin}\phi `$ where $`\phi `$ is the gauge-invariant phase difference of the junction and the voltage $`V`$ across a junction is $`V=(\mathrm{\Phi }_0/2\pi )d\phi /dt`$. Thus, the current of SIS junction, $`i`$, in reduced units is given by

$$i=\ddot{\phi }+\mathrm{\Gamma }\dot{\phi }+\mathrm{sin}\phi =𝒩(\phi ).$$
(2)
This current is normalized by the junction’s intrinsic critical current $`I_c`$ and $`\mathrm{\Gamma }`$, the damping, is usually referred to as the Stewart-McCumber parameter $`\beta _c=\mathrm{\Gamma }^2=2\pi I_cCR^2/\mathrm{\Phi }_0`$ ($`\mathrm{\Phi }_0`$ is the flux quantum). Time is normalized by $`\tau =\sqrt{\mathrm{\Phi }_0C/2\pi I_c}`$. The RCSJ model equation is isomorphic to the equation of a driven pendulum. The mass is normalized to one and the viscous damping is $`\mathrm{\Gamma }`$.

Fig. 1 sketches a typical current-voltage (IV) curve of a single junction. The junction is biased by a DC current and the average voltage is measured. The IV curve of a JJ presents two different physical states. The first is a zero-voltage state, the so-called superconducting or quiet state; a current flows through the junction but no voltage difference appears. This state exists for values of the current smaller than the junction critical current $`I_c`$. At this value of the current the junction switches to the superconducting gap voltage, $`V_g`$. The gap voltage results from the breaking of Cooper-pairs and causes the junction’s resistance, and thereby the damping, to change in a complicated nonlinear way. If we increase the applied current further, the junction reaches its normal state and it behaves as a resistor, $`R_n`$. These resistive states are also called rotating or whirling modes. As the current decreases the junction returns to the gap voltage and then to its zero-voltage state at the retrapping current, $`I_r`$. Thus, the resistive state occurs for values of the current above the junction retrapping current and coexists with the superconducting state for currents between the critical and the retrapping values. The amplitude of this hysteretic loop is governed by the value of the damping $`\mathrm{\Gamma }`$. The behavior of the curve (Fig. 1) close to $`V_g`$ can be modeled using the RCSJ model with an appropriate nonlinear voltage dependence for the resistance (see section 5.2) or by using other more sophisticated models .

We can design JJ arrays of different geometries and parameters. Networks of junctions are valuable model systems for the study of coupled non-linear oscillators. For instance, solid-state physical realizations of the Frenkel-Kontorova model for dislocations  and the two-dimensional X-Y model for phase transitions in condensed matter  are two prominent examples. There has also been extensive studies on the soliton dynamics in one-dimensional Josephson arrays . The experimental advantage of Josephson networks is that there is good control over the parameters because they are fabricated microelectronic solid-state circuits. Moreover, these networks can be designed for a wide range of oscillator parameters from the extremely underdamped to overdamped limits.

Here we present an in-depth study of the experiments reported in . We use more complete models to understand the dynamics of the array. In particular we will not assume uniform current bias, the effects of temperature will be discussed, a nonlinear junction resistance is added, and a full-inductance matrix will be used. In the next section we introduce the governing equations and explain intrinsic localization in the ladder. In section 3 we will report on the experimental study of DB in our superconducting arrays. We develop in section 4 a simple circuit model which allows for the understanding of most of the experimental findings. Numerical simulations will be shown section 5. In section 6 we perform a linear analysis that yields resonance frequencies and decay lengths of excitations. Section 7.1 is devoted to a numerical study of single–site DB solutions in the ladder. There we compute regions of existence of DB at different array parameters, and we will analyze typical simulated IV curves for the cases of type-B and type-A DB. We also study the vortex patterns associated to both solutions. The major results are summarized in the final section.

## 2 Josephson Ladders

### 2.1 Ladder equations

Fig. 2 shows the circuit diagram for the Josephson array in ladder geometry and with uniform current injection. The junctions are marked by an “$`\times `$”. Horizontal junctions have critical currents of $`I_{ch}`$ while vertical junctions have a critical current of $`I_{cv}`$. Anisotropic arrays are fabricated by varying the area of the junctions. In our junctions, the critical current and capacitance are proportional to this area. Due to the constant $`I_cR_n`$ product, the normal state resistance is inversely proportional to the junction area. The anisotropy parameter $`h`$ can then be defined as $`h=I_{ch}/I_{cv}=C_h/C_v=R_v/R_h`$.

A ladder is a useful geometry to study DB because vertical junctions can play the role of pendula while the horizontal junction can act as controllable non-linear coupling. Thereby, a ladder can be roughly thought of as a one-dimensional chain of pendula that are coupled by non-linear springs. If the phases of the vertical junctions are interpreted as as particle coordinates, then the ladder is in essence a one-dimensional chain of particles with non-linear on-site potential as well as non-linear interactions since the horizontal junctions are nonlinear Josephson elements. It is known that lattices with non-linear on-site potential or non-linear interactions are needed to support DB excitations.

The junctions in the array are coupled by means of current conservation and fluxoid quantization. Kirchhoff’s current conservation law (KCL) at the top node yields

$$I_{j1}^t+I_{\mathrm{ext}}I_j^tI_j^v=0.$$
(3)
A consequence of the open boundary conditions is that the current on the top horizontal junctions must be equal but opposite to the current in the bottom horizontal junctions. Thus, $`I_j^m=I_j^t=I_j^b`$. We will normalize all the currents by $`I_{cv}`$. Also, we will refer to horizontal junctions by the superscript $`h`$ (which is not to be confused with the anisotropy) when we are dealing with either the top or bottom horizontal junction.

Fluxoid quantization causes the circulation of the gauge invariant phase differences around a loop to be equal to the flux of the total (external plus induced) magnetic field through the loop. When we impose this condition on one of our cells and assume only external and self-induced fields, that is, we neglect mutual inductances between different cells, we find

$$\phi _j^v\phi _{j+1}^v\phi _j^t+\phi _j^b+2\pi f+\frac{1}{\lambda }i_j^m=0.$$
(4)
Here, $`f=\mathrm{\Phi }_{\mathrm{ext}}/\mathrm{\Phi }_0`$ is the normalized applied flux per unit cell and $`\lambda =\mathrm{\Phi }_0/2\pi I_{cv}L_s`$ where $`L_s`$ is the self-inductance of the loop and $`i_j^m`$ is the normalized mesh current so that $`i_j^m/\lambda `$ is the normalized self-induced flux of a cell.

The vorticity, $`n_j`$ is defined through the expression

$$[\phi _j^v][\phi _{j+1}^v][\phi _j^t]+[\phi _j^b]=2\pi (n_jff_j^{\mathrm{ind}})$$
(5)
where $`[\phi ]`$ represents the phases modulus $`2\pi `$, and $`f_j^{\mathrm{ind}}=i_j^m/2\pi \lambda `$. This expression is equivalent to Eq. (4) and thus also referred to as fluxoid quantization.

We let the functional $`𝒩(\phi )=\ddot{\phi }+\mathrm{\Gamma }\dot{\phi }+\mathrm{sin}\phi `$ represents the current through a junction in the RCSJ model. The resulting set of nonlinear coupled equations can be written as

$`𝒩(\phi _j^t)`$ $`=`$ $`{\displaystyle \frac{\lambda }{h}}\{\phi _j^v\phi _{j+1}^v\phi _j^t+\phi _j^b+2\pi f\}`$
$`𝒩(\phi _j^v)`$ $`=`$ $`\lambda \{\phi _{j+1}^v2\phi _j^v+\phi _{j1}^v+\phi _j^t\phi _{j1}^t\phi _j^b+\phi _{j1}^b\}+i_{\mathrm{ext}}`$
$`𝒩(\phi _j^b)`$ $`=`$ $`{\displaystyle \frac{\lambda }{h}}\{\phi _j^v\phi _{j+1}^v\phi _j^t+\phi _j^b+2\pi f\}`$ (6)
We can identify in these equations a discrete Laplacian term $`^2\phi _j^v=\phi _{j+1}^v2\phi _j^v+\phi _{j1}^v`$, which accounts for the interaction between vertical junctions, and discrete first order derivatives $`\delta _x\phi _j^v=\phi _{j+1}^v\phi _j^v`$ and $`\delta _x\phi _{j1}^h=\phi _j^h\phi _{j1}^h`$, which account for the interaction terms between vertical and horizontal junctions respectively.

The external current is normalized as $`i_{\mathrm{ext}}=I_{\mathrm{ext}}/I_{cv}`$. The damping is $`\mathrm{\Gamma }=\sqrt{\mathrm{\Phi }_0/2\pi I_{cv}R_v^2C_v}`$. We note that because the anisotropy in our arrays is caused by varying the junction area, $`\mathrm{\Gamma }`$ is the same for every junction in the array.

In Eq. (6), $`j=1`$ to $`N`$ and at the open boundaries,

$`\phi _0^t`$ $`=`$ $`\phi _N^t=0`$
$`\phi _0^b`$ $`=`$ $`\phi _N^b=0`$
$`\phi _{N+1}^v`$ $`=`$ $`\phi _N^v+2\pi f`$
$`\phi _0^v`$ $`=`$ $`\phi _1^v2\pi f`$ (7)
where the phases at $`j=0`$ and $`j=N+1`$ are for mathematical convenience and do not represent real junctions.

The physical currents through the junctions are $`I_j^t=I_{ch}𝒩(\phi _j^t)=I_j^b=I_{ch}𝒩(\phi _j^b)`$ and $`I_j^v=I_{cv}𝒩(\phi _j^v)`$.

We can use these governing equations to compare the ladder with a JJ parallel array, a system with broad interest. In a parallel array the horizontal junctions of the ladder are replaced by superconducting wires which have a linear current-phase relation. In the ladder, the dynamics of the horizontal junctions can also be described by a linearized constitutive relation under some restricted circumstances. For instance, in the static case, when the horizontal junctions have no time dependence then $`i^h=h\mathrm{sin}\phi ^hh\phi ^h`$. From the normalized Eq. (3),

$$i_j^vi_{\mathrm{ext}}=i_{j1}^ti_j^t$$
(8)
The right hand side is simply $`h(\phi _{j1}^t\phi _j^t)`$ and we can find a similar relation for the bottom horizontal junctions. We can substitute these linearized relations for $`\phi ^t`$ and $`\phi ^b`$ in Eq. (2.1) to get,

$$𝒩(\phi _j^v)=\frac{h\lambda }{h+2\lambda }^2\phi _j^v+i_{\mathrm{ext}}$$
(9)
This is the discrete sine-Gordon equation with a renormalized discreteness parameter of $`\lambda _{\mathrm{eff}}=h\lambda /(h+2\lambda )`$ and is equivalent to the governing equations of a Josephson-junction parallel array .

The difference between the parallel array and the ladder is obviously the existence of horizontal junctions in the case of the ladder. The use of $`\lambda _{\mathrm{eff}}`$ to map the ladder to the discrete sine-Gordon model is correct only for the study of a reduced set of possible states of the array: Those for which only the convex part of the inter-phase interaction is relevant. If we can neglect the dynamics of the horizontal junctions (for instance when studying static properties) then the above is a very good approximation as has been rigorously stated in . The rotating DB states studied in this paper are a good example of a set of solutions of the ladder which do not appear in the dynamics of the discrete sine-Gordon model.

### 2.2 Localization in the Josephson ladder

In Fig. 1 we showed that the IV curve of an underdamped JJ has an hysteresis loop for current values between the critical and the retrapping currents. In this current range the zero-voltage ($`V=0`$) and rotating ($`V=V_g`$) attractors coexist, and it is this hysteresis loop that allows for the existence of breathers in the ladder with DC bias current. In the full ladder, the phase space is much more complex. However, when the vertical junctions are weakly coupled then it is possible that each of the junction attractors is essentially independent. A breather is then a localized state where one vertical junction is rotating while the others are in the zero voltage state. Under general considerations of $`N`$ nonlinear oscillators, it is not obvious that the phase space of the coupled system supports attractors of localized solutions. Indeed, in linear systems even when the oscillators are weakly coupled, the phase space does not support stable localized solutions.

The coupling of vertical junctions occurs through horizontal junctions, geometrical inductances, and fluxoid quantization. But the most important contribution is from horizontal junctions and their influence is measured by the parameter $`h`$. If the coupling is too strong, i.e. $`h`$ is too large, then localized solution cannot exist. We numerically found that small $`\lambda `$ and $`h=0.25`$ are adequate values of the coupling parameters for studying DB .

Fig. 3 shows some of the possible configurations for single-site breathers in our ladders. The arrows indicate voltage polarity. Junctions without arrows are in the zero-voltage state. The solutions in Fig. 3 are single-site breathers because only one vertical junction is rotating. We also see that due to Kirchhoff’s voltage law (KVL), there must be at least one other junction rotating in each of the neighboring cells.

The actual number and pattern of rotating horizontal junctions will determine the type of breather. We call breathers that have two rotating horizontal junctions, like Fig. 3(a) and Fig. 3(c), type A breathers. We see that in this case the voltages of all the rotating junctions are the same and so the breather solution is highly constrained. We will use the term asymmetric to refer to this set of type $`A`$ solutions.

In the case of the type $`B`$ solution, the four nearest horizontal junctions rotate as depicted in Fig. 3(b). We will use the term symmetric to refer to type $`B`$ solution. However, it is important to note that symmetric refers only to the fact that top and bottom horizontal junction in the same cell rotate. It does not always correspond to true up-down symmetric solutions. An up-down symmetric solution is a solution for which $`\phi _j^t(t)=\phi _j^b(t)\mathrm{mod}\mathrm{\hspace{0.17em}2}\pi `$. This is a solution of the system as can be inferred from Eqs. (6), where it is clear that $`𝒩(\phi _j^t)=𝒩(\phi _j^b)`$. However, for certain values of the parameters, the underdamped character of the junctions allows for different solutions, such as the type A breathers, that do not obey the up-down symmetry condition. Also, for most type B solutions the magnitude of the voltage of a horizontal rotating junction is one half the voltage of the rotating vertical junction. Again, not all type B solutions obey this condition.

Single-site breather solutions that have 3 horizontal junctions rotating such as those shown in Fig. 3(d) are called hybrids of type A and B. The other possible single-site breathers that are not shown can also be classified as either A, B or hybrid.

In the rest of this paper we will designate librating junctions as shorts and junctions that are rotating by an $`\times `$ as shown at the top of Fig. 4. Each branch represents a junction, but only the ones with a cross are rotating. Fig. 4(a) shows a type A breather and a plot of the numerically calculated DC flux per unit cells for a 9 junction ladder array. Fig. 4(b) shows the decay of the DC flux for a type B breather. The flux decays exponentially so that our breather solutions are exponentially localized. In the case of the hybrid breather shown in Fig. 3(d), the decay of the flux in the right side of the array is the same as for the type B shown in Fig. 4(b) and the decay on the left side of the array is the same as the type A shown in Fig. 4(a). Note that we have plotted the modulus of the flux. The fluxes for all the cells of one side of the vertical rotating junction have opposite sign to the fluxes of the cells in the other side.

Fig. 5 shows different type B breather solutions for which a set of contiguous vertical junctions are in the rotating state. We call them multi-site or $`m`$-site solutions where $`m`$ refers to the number of vertical junctions rotating in the array. In general many different solutions are allowed where each vertical junction can be in the resistive or in the superconducting state while one or both horizontal junctions between vertical ones in different states also rotate.

The above discussion focused on the existence of rotobreathers in DC-biased arrays. In the case of AC-biased arrays the JJ ladder in addition to rotating localized modes also supports oscillating localized modes or oscillobreathers where one vertical junction describes a large amplitude oscillation when comparing with other vertical junctions. Such modes were studied in  for non-inductive arrays. These modes persist when inductances are added to the model and also different type-B and type-A families of AC-biased rotobreathers can be identified .

## 3 Experiments

### 3.1 Experimental observation of breathers

We have designed and measured several anisotropic ladders. The inset of Fig. 6 shows a schematic of the measured arrays. The junctions are fabricated using a Nb-Al<sub>2</sub>O<sub>x</sub>-Nb tri-layer technology with a critical current density of about $`1\mathrm{kA}/\mathrm{cm}^2`$ . The current is injected and extracted through bias resistors in order to distribute it as uniformly as possible in the array. These resistors are large enough so as to minimize any deleterious effects on the dynamics. Our ladders have $`3\times 3\mu \mathrm{m}^2`$ horizontal junctions and $`6\times 6\mu \mathrm{m}^2`$ vertical junctions. Vertical junctions have been designed with four times the area of the horizontal ones. Thus, the anisotropy ratio $`h`$ is approximately $`0.25`$. The bias resistors are $`25\mathrm{\Omega }`$.

There are voltage probes in the fourth, fifth, sixth and ninth vertical junctions to measure $`V_4`$, $`V_5`$, $`V_6`$ and $`V_9`$. The voltage probes can also be used to measure the top horizontal junctions in the middle, which we denote as $`V_{4T}`$ and $`V_{5T}`$, or any other combination of terminals.

From the measured normal state resistance we calculate $`I_{cv}=360\mu \mathrm{A}`$ and $`I_{ch}=90\mu \mathrm{A}`$ at $`T`$ = 0 K. The dimensionless penetration depth $`\lambda `$, which measures the inductive coupling in the array, is defined as $`\mathrm{\Phi }_0/2\pi L_sI_{cv}`$. We estimate the mesh inductance $`L_s=30.2\mathrm{pH}`$ from numerical modeling  so $`\lambda =0.04`$ at $`T`$ = 0 K.

To determine $`\mathrm{\Gamma }`$ we need to measure the subgap resistance. Different approaches can be used. One possibility is to calculate this damping from the measured return current of the junction. The model used for the return current determines the subgap dissipation and resistance . We will estimate the value of $`\mathrm{\Gamma }`$ in our experiments from the retrapping current $`I_r`$ by the relation $`I_r/NI_{cv}=(4/\pi )\mathrm{\Gamma }`$, where $`N`$ is the number of vertical junctions. This expression can be found from a simple energy argument valid for small values of the damping. The energy injected into the system by the applied current must equal the energy to “rotate” the junction one full period . Thus from the experiment shown in Fig. 6 we infer a value of $`\mathrm{\Gamma }0.08`$ ($`\beta _c160`$) at $`T=5.2\mathrm{K}`$.

We show a typical IV curve of a ladder without DB in Fig. 6. We measure the time-averaged voltage of the 9-th vertical junction as a function of the uniformly applied current. The junction is in a zero-voltage state as we increase the applied current from zero. When the applied current reaches the depinning current $`I_{\mathrm{dep}}`$ at about $`2\mathrm{mA}`$ the junction switches from zero-voltage state to the superconducting gap voltage, $`V_g`$, which at this temperature is $`2.5\mathrm{mV}`$. If we increase the applied current further, the junction reaches its normal state and it behaves as a resistor, $`R_n`$, of $`5\mathrm{\Omega }`$. As the current decreases the junction returns to the gap voltage and then to its zero-voltage state at the retrapping current, $`I_r0.2\mathrm{mA}`$.

Sometimes, when we sweep the applied current we find that DB solutions appear spontaneously: They can be thermally excited when the applied current is close to $`I_{\mathrm{dep}}`$. However for our experiments, we have developed a simple reproducible method of exciting a breather: (i) Bias the array uniformly to a current below depinning current; (ii) increase the current injected into the middle vertical junction ($`V_5`$) until its voltage switches to the gap; (iii) reduce this extra current in the middle junction to zero. Other procedures are possible. For instance, we can increase the current for the middle vertical junction first until it rotates and then increase the array bias current as was described in . We could also inject this extra current into a horizontal junction. All of these methods produce breathers in our arrays and thereby hinting at the generic nature of breather solutions in our ladders

Fig. 7 shows the result after we have excited the breather and we have increased the array current. The breather is excited at $`I_a1.4\mathrm{mA}`$ and then the junction voltages are measured as the applied current is increased. We simultaneously measure the voltages of the vertical junctions ($`V_4`$, $`V_5`$ and $`V_6`$) and the top two horizontal junctions, $`V_{4T}`$ and $`V_{5T}`$. The DB is localized in the fifth vertical junction and is a type B breather since the top, $`V_{5T}`$, and bottom, $`V_{5B}`$, horizontal junctions have opposite voltage. We also find that the horizontal voltages are half in magnitude to the rotating vertical junction $`V_5`$. This is the B1 breather state shown in Fig. 5.

The breather exists until a maximum current $`I_+1.95\mathrm{mA}`$ is reached. If the applied current is further increased then Fig. 7 shows that the voltages of the 4th and 6th vertical junctions jump to the gap voltage while those of the horizontal junctions go to zero. In fact all the vertical junction have jumped to the gap voltage, The ladder has now all of the vertical junctions rotating and it is in the “whirling” state. The state is depicted as W in Fig. 5. In section 4 we will present a circuit model that will relate $`I_+`$ to the array depinning current.

In Fig. 8 we measure the voltage as the applied external current is decreased. We excite the breather at $`I_a1.25\mathrm{mA}`$ and then the voltages are measured as the current is decreased. We have measured vertical junctions four, five and six. We have also measured the total voltage of the horizontal junctions five through eight and this sum voltage is referred as $`V_{\mathrm{top}}`$. This allows us to reconstruct the $`m`$-site breather state in terms of the rotating junctions of the array.

Our experiments show that as we decrease the applied current the single-site breather will usually decay into an $`m`$-site breather state. From $`I_a1.25\mathrm{mA}`$ to $`I_a1.05\mathrm{mA}`$ the vertical voltages $`V_4=V_6=0`$. This is the single-site type B breather, the state depicted as B1 in Fig. 5. If we further decrease the applied current from $`I_{}=1.05\mathrm{mA}`$, $`V_4`$ jumps from zero to the gap voltage. This is a two site breather shown schematically in Fig. 5 as B2. As we further decrease the current we can count nine discontinuous curves, each one corresponding to the switching of a vertical junction. At 0.3 mA all of the vertical junctions return to their zero-voltage state via a retrapping mechanism analogous to that of a single pendulum.

From these experiments we conclude that this shifting of the voltage corresponds to at least one vertical junction switching from the zero-voltage state to the rotating state. If we assume that only consecutive junctions switch, then every curve in the measurement can be associated with one $`m`$-site breather. Since we have measured $`V_{\mathrm{top}}`$ we can reconstruct the ladder solutions. Fig. 5 shows a schematic of the states measured in Fig. 8.

The shapes of the IV curves in this multi-site breather regime are influenced by the junction nonlinear resistance and the redistribution of current when each vertical junction switches. This redistribution may also govern the evolution of the system after each transition to one of the other possible breather attractors in the phase space of the array.

Fig. 8 shows a measurement at $`T=6.0`$K and we find 10 different states can be distinguished (8 $`m`$-site breather states, the whirling state, and the zero voltage state). In general we only see three or four $`m`$-site breathers as we decrease the current, and often we see different ones even under similar experimental conditions.

Fig. 9 shows an IV which includes a type A breather. First, a type B breather is excited at $`I_a0.9\mathrm{mA}`$ and the current is decreased. As we decrease the current $`V_4`$ and $`V_6`$ have zero voltage while $`V_5=2V_{5T}`$. At $`I_a0.58\mathrm{mA}`$ the type B breather becomes unstable. In most of our measurements, the array will then suddenly switch to an $`m`$ site breather as shown in Fig. 8, but in this case we find a new switching behavior: The type B breather switches to a type A solution. This state is a type A breather because the voltage $`V_5=V_{5T}`$ and $`V_4`$ and $`V_6`$ are zero. As the current decreases the type A breather disappears at $`I_a\mathrm{\hspace{0.17em}0.53}\mathrm{mA}`$ and the array is in the superconducting state. In our experiments, the ladder always returns to the superconducting state whenever a type A breather reaches its minimum current.

As the current decreases the type A breather is only accessible for a small current range of $`0.05\mathrm{mA}`$. However, it is possible to bias on the type A breather and increase the current to trace out the hysteresis loop. The tracing of the type A breather voltage step is also shown in Fig. 9. We see that the type A breather exists up to a current of $`0.7\mathrm{mA}`$. Once it becomes unstable the array dynamics usually jumps back to the B1 breather.

### 3.2 Temperature and magnetic field dependence

By sweeping the temperature and magnetic field we can study how the current range in which our breather exists is affected by a change of the array parameters. We define and study the evolution of four current values of importance: the current when the array returns to the zero-voltage, $`I_r`$; the maximum zero-voltage state current, $`I_{\mathrm{dep}}`$; and the maximum and minimum current for a breather state, $`I_+`$ and $`I_{}`$.

Fig. 10 shows the typical dependence of $`I_+`$ and $`I_{}`$ for a B1 breather as a function of applied magnetic field. The IV curves were measured by applying a perpendicular magnetic field of 0 to 300 mG using a magnetic coil that is mounted on the radiation shield of our probe. There is some $`f`$ dependence. This can be expected since $`I_+`$ is related to the array depinning current. But at this temperature $`\lambda =0.06`$ so we expect to have a large Meissner current and a correspondingly relatively flat $`I_c`$ dependence vs $`f`$. If $`\lambda `$ were larger then the breather dynamics and the $`I_c`$ should show a stronger magnetic field dependence.

We can study further the existence region of the breather by changing the temperature of the sample. In this way, we can vary the parameters to a certain degree. The temperature causes the $`I_{cv}`$ of the junction to change and hence change $`\mathrm{\Gamma }`$ and $`\lambda `$. For our arrays, the junction parameters can range from $`0.031<\mathrm{\Gamma }<0.61`$ and $`0.04<\lambda <0.43`$ as the temperature varies from 4.2 K to 9.2 K.

Fig. 11 shows how the maximum and minimum current of both type A and type B breathers are affected as we vary $`\mathrm{\Gamma }`$. In Fig. 11 $`\mathrm{\Gamma }<0.2`$ corresponds to $`T<6.7\mathrm{K}`$ and $`\lambda 0.05`$. At these low temperatures, $`I_{cv}`$ essentially remains constant so $`\lambda `$ does not vary. However, the subgap resistance varies substantially as can be seen from the retrapping current measurements in Fig. 11. Therefore, there is a larger variation in $`\mathrm{\Gamma }`$. Figure also shows a nice agreement in between the experiments (points) and the theory (lines) except for the case of the maximum current of the type A solutions. Also, experimentally we did not find breathers for $`\mathrm{\Gamma }>0.2`$.

## 4 Circuit model

### 4.1 Introduction

In this section, we will develop a simplified DC circuit model of our array in order to understand the region of existence of the breather solutions in our experiments. Also, this model will allow for evaluating the effect of the bias resistors in the dynamics of the array.

In the experiment we apply the external current through bias resistors as shown in the inset of Fig. 6 in order to distribute it uniformly. If the resistors are very small, the horizontal junctions are effectively shunted by a small resistance and no DB solutions can be excited. If the resistors are large so that they dominate over every other impedance, then applied current will be almost uniform throughout the array. However, one drawback of using large bias resistors is that they will create local heating of the sample and affect the measurements.

Thus, when studying DB, there are at least two complications with the bias resistors. The first is that when we excite a breather state, the applied current is not completely uniform. So questions arise about how this non-uniformity affects localization. Also as we decrease the current we see transitions between different $`m`$-site breather states. In each of these transitions vertical junctions switch from the zero voltage state to the rotating state. This switching causes a change in the junction impedance and thereby affects the current distribution. This redistribution might be important to understand the dominant drive of the pattern selection process between $`m`$-site breathers. However, we will show below that the effect of the bias resistors only adds a small correction to the calculation of the existence region of breather states.

To get some physical intuition we will use a simple model where rotating vertical junctions have a resistance of $`R_v`$ and rotating horizontal junctions have a resistance of $`R_h`$. Librating junction will be modeled as shorts. We will reduce the array to a simple network of resistors and calculate DC properties. The equivalent resistor network for a single-site symmetric breather located on junction 5 in our 9 junction array is shown in Fig. 12.

When the array is in the superconducting or in the whirling states, the applied current distributes uniformly across the vertical junctions. This will not be the case when we have a breather since some junctions are in the resistive state while others are superconducting. Moreover, when we have a breather and we measure the voltage of the fifth vertical junction (Fig. 8), we find that its voltage shifts to a higher value when new vertical junctions switch to the rotating state. These shifts are identified as jumps in the effective current biasing the junction due to the redistribution of currents in the array. Roughly speaking, most of the applied current wants to flow through the superconducting junctions. Whenever a vertical junction switches from a superconducting to a resistive state there is some extra current that is distributed throughout the array and thereby the effective bias of the fifth junction becomes a little larger. This extra bias results in a jump of the measured voltage to a larger value. The size of the voltage jump is dependent on the non-linear subgap resistance and is usually large even for small changes of the effective bias current.

### 4.2 IV curves and current distribution

In this section we will focus on interpreting this redistribution and shifting of the voltage by using a simple DC model of the ladder. We use $`I_a`$ to represent the total applied current in the array. This is the current of the experimental current source. The current applied to a particular junction through the bias resistors, i.e. the current through the bias resistor, will be designated $`I_j`$. We will use $`I_l`$ to denote the current from node a to node b. This current $`I_l`$ is the sum of the currents $`I_1`$ through $`I_4`$. Also, because all of our bias resistors are the same and due to KVL $`I_1=I_2=I_3=I_4`$. So, $`I_l=4I_4`$ for instance. The current through node a to node c is $`I_5`$ since the breather is located in junction 5.

For the type B breather we measured, the voltage is twice as large as the voltage of the horizontal junction

$$2V_h=V_v.$$
(10)
Therefore

$$I_v=2\frac{R_h}{R_v}I_h,$$
(11)
and we will assume that $`R_v=hR_h`$.

The circuit has left-right symmetry so the current through the right branches must equal the current through the left. KCL at node $`a`$ yields $`I_a=2I_l+I_5`$ where $`I_a`$ is the total current applied, and at node $`c`$ yields $`I_5=2I_h+I_v`$. KVL on the top left mesh gives $`I_lR_b/4I_hR_hI_5R_b=0`$. Combining the equations results in

$$\frac{I_a}{N}=\left\{2+\frac{2}{h}+\left(1\frac{1}{N}\right)\frac{R_h}{R_b}\right\}I_h.$$
(12)
To generalize to type A breathers we just note that the horizontal junction voltage is the same as the vertical junction voltage $`V_h=V_v`$. We can write $`sV_h=V_v`$ where $`s=1`$ for type A breathers and $`s=2`$ for type B breathers. We can also generalize for $`m`$-site breathers by modifying the equivalent DC circuit accordingly. For instance, when $`m=2`$ both rotating vertical junctions can be lumped into an equivalent impedance of $`R_v/2`$. Generalizing our circuit results in <sup>1</sup><sup>1</sup>1These equations are valid for type A and type B solutions with $`m`$ consecutive rotating vertical junctions and with junction 1 and $`N`$ in the no rotating state. In a similar way, we can compute the equations for the case when junction 1 or $`N`$ is rotating, for hybrid breathers, or for multi-site breathers for which the rotating junctions are not consecutive.

$$\frac{I_a}{N}=\left\{\frac{2}{m}+\frac{s}{h}+\left(1\frac{m}{N}\right)\frac{R_h}{R_b}\right\}I_h.$$
(13)
We can also calculate the IV curve of a vertical rotating junction by substituting $`I_h=(h/s)I_v=(h/s)V_v/R_v`$ in Eq. (13). The result is

$$\frac{I_a}{N}=\left\{1+\frac{2h}{sm}+\left(1\frac{m}{N}\right)\frac{hR_h}{sR_b}\right\}\frac{V_v}{R_v}.$$
(14)
Another important variable is the amount of effective current biasing a vertical rotating junction, say $`I_5`$,

$$\frac{I_a}{N}=\left\{1+\frac{mh}{2h+sm}\left(1\frac{m}{N}\right)\frac{R_h}{R_b}\right\}I_5.$$
(15)
From our experiments $`R_h/R_b0.8`$, $`h=0.25`$ and $`N=9`$. Thus for the case of a type-B single-site breather ($`m=1`$, $`s=2`$) we get $`I_5=0.934I_a/9`$ and $`I_j=1.008I_a/9`$ ($`j5`$). Thus the vertical rotating junction is biased by a smaller DC current than the quiet vertical junctions. Also this non uniformity disappears as the bias resistance is made larger.

### 4.3 Estimation of $`I_{}`$

To calculate $`I_{}`$, we assume that the instability that destroys an $`m`$-site breather is due to a junction retrapping mechanism. Even if it appears that some other mechanism is important, like a resonance, this instability occurs in the subgap region of the junction where the voltage varies rapidly while the current remains relatively constant. Therefore, a retrapping current should give a good estimate of $`I_{}`$ regardless of the physical mechanism.

From simple energy consideration of an isolated junction, the retrapping current can be estimated as $`4\mathrm{\Gamma }I_c/\pi `$. Then, when the horizontal junction reaches its retrapping current $`I_h=4\mathrm{\Gamma }I_{ch}/\pi `$. Thus, $`I_v=4\mathrm{\Gamma }I_{ch}s/\pi h=4\mathrm{\Gamma }I_{cv}s/\pi `$, and the vertical junction is at $`s`$ times the junction retrapping current. Conversely, if the vertical junction is at the retrapping current, then the horizontal junction is at $`1/s`$ of the junction retrapping current. Therefore, as the applied current decreases, the horizontal junction always reaches the retrapping current first for a type B breather and both the horizontal and vertical junctions reach the retrapping at the same time for a type A breather. This is the reason why type B breathers can decay into $`m`$-site breathers while type A breathers apparently decay into the superconducting state. When a type A breather reaches $`I_{}`$ all of the rotating junctions retrap and the resulting state is more likely to be the superconducting state. On the other hand, when a type B breather reaches $`I_{}`$ only the horizontal junction retraps while the vertical junctions remain whirling.

From Eq. (13),

$$\frac{I_{}}{N}=\left\{\frac{2}{m}+\frac{s}{h}+\left(1\frac{m}{N}\right)\frac{R_h}{R_b}\right\}\frac{4}{\pi }\mathrm{\Gamma }I_{ch}$$
(16)
when $`I_h`$ reaches its retrapping value.

To use this formula we need the horizontal junction parameters. We can estimate $`C_h`$ from the specific capacitance of the tri-layer and the junction area. We find $`C_h=300\mathrm{fF}`$. From the constant $`I_cR_n`$ product, we find that $`I_{ch}=90\mu \mathrm{A}`$. However, the remaining parameter $`R_h`$ is more difficult to estimate since it depends highly on the nonlinear subgap region. Instead of trying to calculate $`R_h`$ directly, we will fit Eq. (16) to our measured $`I_{}`$ using $`R_h`$ as our fitting parameter. We include the effect of $`R_h`$ through the definition of $`\mathrm{\Gamma }=\sqrt{\mathrm{\Phi }_0/2\pi I_{ch}R_h^2C_h}`$. Then from Eq. (16),

$$\frac{I_{}}{N}=\left\{\frac{2}{m}+\frac{s}{h}+\left(1\frac{m}{N}\right)\frac{R_h}{R_b}\right\}\frac{4}{\pi }\sqrt{\frac{\mathrm{\Phi }_0}{2\pi R_h^2C_h}I_{ch}}$$
(17)
and $`R_h`$ as the only free parameter. This is how we estimate $`R_h`$, or equivalently $`\mathrm{\Gamma }`$, for a given temperature.

The inset of Fig. 8 shows the fitted values for $`R_h`$ and the dashed lines in Fig. 8 show the resulting $`I_{}`$. We note that for solutions B7, B8, and $`W`$ we used a different equivalent circuit based on the schematic of Fig. 5. We see that $`R_h`$ is approximately $`50\mathrm{\Omega }`$ for all $`m`$. This value is not totally unexpected. Roughly speaking, the horizontal junctions are shunted by two bias resistor of $`25\mathrm{\Omega }`$ each and the retrapping current depends strongly on the equivalent junction impedance. Since the subgap resistance for our junctions can be several $`\mathrm{k}\mathrm{\Omega }`$ the equivalent horizontal junction impedance will be dominated by the shunts which add up to $`50\mathrm{\Omega }`$.

We can also easily understand why $`I_{}`$ decreases as $`m`$ increases. For a constant array bias at $`m=1`$ some fraction of the current will flow through both horizontal junctions. At the same bias current for $`m=2`$, there is a horizontal junction that is not rotating in between the rotating horizontal junctions. By symmetry considerations there is no current flowing through this quiet horizontal junction. Since the applied bias current is the same, this implies that a larger fraction of the applied current must flow through the rotating horizontal junctions when $`m=2`$ than when $`m=1`$. Thereby, $`I_{}`$ is smaller for $`m=2`$.

### 4.4 Estimation of $`I_+`$

To calculate $`I_+`$ we look for the librating junction that supports the maximum current. When this junction reaches $`I_c`$, the breather has reached its maximum current. It is straightforward to find that the critical junction is the first vertical junction that is not rotating (the one nearest to the rotating ones). Let $`I^{}`$ be the current through the junction. KCL at node b of the array yields $`I^{}=I_h+I_4`$. Here we have assumed that there is not current in the horizontal quiet junctions. Since $`I_l=4I_4`$ we can substitute for the currents to solve for $`I^{}`$ in terms of $`I_a`$. We can also generalize to $`m`$-site breathers. The result is

$$\frac{I_a}{N}=\frac{2h/m+s+h\left(1m/N\right)R_h/R_b}{h(1+2/m)+s+hR_h/R_b}I^{}.$$
(18)
The maximum applied current $`I_+`$ occurs when $`I^{}=I_{cv}`$. In the limit $`R_bR_h`$,

$$\frac{I_+}{N}=\frac{2h+sm}{h(m+2)+sm}I_{cv}.$$
(19)
This current will underestimate the actual value for $`I_+`$. This is because we have not taken into account any of the horizontal junction currents that are in the quiet state.

The effect of fluxoid quantization is to redistribute the currents of the quiet junctions. The currents in the bias resistors and the rotating junctions will remain unaffected. If we consider the effect of the next nearest mesh to the breather, then KCL at node b will be $`I_4+I_h=I^{}+I_{ch}\mathrm{sin}(k_h)`$. Here, $`I_{ch}\mathrm{sin}(k_h)`$ is the current in the next horizontal junction. Fluxoid quantization in this quiet mesh yields $`k^{}k_v2k_h=0`$ when $`f=0`$ and we neglect the self induced field. Here $`k^{}=\mathrm{sin}^1(I^{}/I_{cv})`$ and $`k_v`$ is the phase of the next quite vertical junction. Adding self fields will tend to decrease $`I_{ch}\mathrm{sin}(k_h)`$ because smaller screening currents are needed as the inductance becomes larger. This correction will then tend to overestimate $`I_+`$. We note that to calculate $`I_+`$ we set $`I^{}=I_{cv}`$ so $`k^{}=\pi /2`$.

With fluxoid quantization $`k_h=\pi /4k_v/2`$ at $`I_a=I_+`$. To first order we expect the current in the quiet vertical junction to be $`I_a/N`$ and $`k_v=\mathrm{sin}^1(I_a/NI_{cv})`$. For consistency, we apply KCL so $`I_3+I_{ch}\mathrm{sin}(k_h)=I_{cv}\mathrm{sin}(k_v)`$. Again, we neglect the current of the next quiet horizontal junctions. We can solve for $`k_v`$ in terms of $`k_h`$ and substitute back into the fluxoid quantization condition. Using $`\mathrm{sin}xx`$ and $`\mathrm{cos}x1x^2/2`$ we can get a closed expression

$$k_h=\frac{h+\sqrt{h^2+8(1I_3/I_{cv})}}{4}.$$
(20)
This expression is only valid when $`I_aI_+`$.

The current $`I_3`$ equals $`I_4=I_l/4`$ and can be calculated from Fig. 12 and is simply

$$I_3=\frac{2/m+s/h+R_h/R_b}{2/m+s/h+(1m/N)R_h/R_b}\frac{I_a}{N}.$$
(21)
In the limit $`R_bR_h`$ $`I_3=I_a/N`$.

The maximum current will now be Eq. (18) when $`I^{}=I_{cv}+I_{ch}\mathrm{sin}(k_h)`$ with $`k_h`$ defined in Eq. (20) and Eq. (21). The resulting equation is transcendental. But we know the correction due to $`k_h`$ should be small so we can linearize the sine term. The equation then has dependencies of $`I_a`$ on both sides but can be solved by isolating the square root and squaring. This results in a quadratic equation in $`I_a`$ with easily extractable roots. The equation can also be solved iteratively. One iteration usually results in a good approximation. Hence, if $`I_+^0`$ is the uncorrected value from Eq. (18) when $`I^{}=I_{cv}`$, then the first correction $`I_+^1`$ can be calculated when $`I^{}=I_{cv}+I_{ch}\mathrm{sin}(k_h)`$ with $`I_a=I_+^0`$ substituted in in Eq. (21).

### 4.5 Comparison with the experiments

The result of the simplified DC circuit model is compared to the measurements in Fig. 8 and Fig. 11. We used the measurement in Fig. 8 to estimate the values of $`R_h`$ and we use Eq. (17) to calculate $`I_{}`$ for the type A and B breathers. In Fig. 11, we see that the result agrees quite well for $`I_{}`$. We have also calculated $`I_+`$ from our uncorrected expression Eq. (18) and with the correction due to Meissner currents. The measured and expected results agree for type B breathers. However, the measured $`I_+`$ for the type A breathers is much lower than what our circuit model predicts.

We also note that the effect of the bias resistors and $`R_h`$ is essentially to provide an offset to the quasi-linear slope in Fig. 11. In the limit $`R_bR_h`$, the predicted values for $`I_{}`$ would intersect ($`0,0`$). The bias resistor, then, only provide a small correction by allowing a better fit to the measurements.

Finally, there appears to be a maximum amount of damping where breathers can exist in the ladder. This maximum $`\mathrm{\Gamma }`$ occurs when $`I_+`$ coincides with $`I_{}`$ which occurs at $`0.2`$. This can be expected from the fact that the DB discussed in this paper require hysteresis in the IV curves.

## 5 Simulations at the experimental values

### 5.1 Introduction

In section 2.1 we derived the standard model for the dynamics of the array. This model assumes the RCSJ for the junction dynamics, uniform bias condition, and the effect of mesh self-inductances. We will use numerical simulations of the governing equations to provide for detailed analysis of the DB solutions. We numerically integrate Eq. (6) using a 4th order Runge-Kutta method. The initial condition is found via a similar procedure as in the experiments. We first bias the array near the depinning current. Then we increase the current of the middle vertical junction until it starts to rotate and finally we decrease this extra current to zero. Usually the resulting solution is periodic and we verify that it is stable by calculating Floquet multipliers. If the solution is not stable we integrate the equations of motion for a long time with an added small current noise source in every junction. In this way we perturb the solution and sample nearby trajectories in phase space. Usually the final result is a new periodic stable solution.

Sometimes, and usually close to the destabilization points, when we integrate the equations we find aperiodic solutions that appear stable. We again add some noise to check the stability of the solution against fluctuations and use standard methods like Poincaré section analysis or study of the Liapunov exponents of the system to gain more information about the behavior of the solution.

### 5.2 Numerically integrated IV curves with a nonlinear resistance

We will present numerical simulations of the dynamics of the ladder. In our samples $`h=0.25`$ and the damping and the penetration depth vary with the temperature. Thus when the temperature varies from 4.2 to 9.2 K, $`\mathrm{\Gamma }`$ varies from 0.03 to 0.6 and $`\lambda `$ from 0.04 to 0.4. However, experimentally we only find breathers for $`T<6.7`$K and that corresponds to $`0.03<\mathrm{\Gamma }<0.2`$ and $`0.04<\lambda <0.05`$. We will then mostly present simulations at $`h=0.25`$, $`\lambda =0.05`$, however the dynamics of the array is very rich and multiple transitions between different attractors occurs when changing the parameters.

Fig. 13 shows simulated IV curves of single-site type B breather solutions. These simulations were done with a fixed $`\mathrm{\Gamma }=0.18`$ and $`\lambda =0.05`$. We have included a subgap resistance in this simulation in order to compare with the experimental measurements more closely. We use the usual approach to extend the RCSJ model where now $`R`$ depends on the junction instantaneous voltage $`V(t)`$. We define a conductance such that

$$G(V(t))=\{\begin{array}{cc}\frac{1}{R_{sg}(T)}\hfill & \text{if }|V(t)|<V_g\hfill \\ \frac{1}{R_n}\hfill & \text{otherwise}\hfill \end{array}$$
(22)
Here $`R_{sg}(T)`$ is taken to be only a function of temperature so for a given set of parameters it is constant. Our functional $`𝒩(\phi )`$ now becomes $`\ddot{\phi }+\mathrm{\Gamma }g(v)\dot{\phi }+\mathrm{sin}\phi `$, where $`g(v)=R_n/R_{sg}(T)`$ and $`\mathrm{\Gamma }`$ is calculated from $`R_n`$.

A simple approach is to model $`g`$ as a continuous hyperbolic tangent. We will use

$$g(v)=g_{sg}+(1g_{gsg})\frac{1\mathrm{tanh}K(1v)}{2}.$$
(23)
In our simulation, we take the value of $`K=100`$. Thereby Eq. (23) approaches a piecewise linear function with a conductance $`g_{sg}`$ at $`v<1`$ and 1 when $`v>1`$.

We find that the simulated curves are similar to the experimentally measured arrays. The inset compares $`I_{}`$ of the simulation to the prediction of the circuit model Eq. (16) when $`R_bR_h`$. The deviation is due to inductance effects. Our simple circuit model neglects the effects of the inductances. We have also verified that $`I_+`$ in the simulations is within 99% of the adjusted $`I_+`$ of our model Eq. (18) and Eq. (21). We have found similar results for numerically computed type A breathers.

There is little effect of the subgap resistance in our simulations besides changing the shape of the IV curve. Instead of using a subgap resistance, we can redefine $`\mathrm{\Gamma }`$ so that it includes the effects of $`g`$. Then, $`\mathrm{\Gamma }`$ is calculated from $`R_{sg}`$ and $`𝒩(\phi )=\ddot{\phi }+\mathrm{\Gamma }\dot{\phi }+\mathrm{sin}\phi `$ as before.

### 5.3 Breather instabilities and switching mechanism

Another aspect of Fig. 13 is the resulting state of the array once the single-site type B breather becomes unstable. Some insight can be gained by first studying how the single-site breather destabilizes. Fig. 14(b) shows the distribution of the Floquet multipliers at several values of the applied current close to $`I_{}`$ for the type B breather when $`g_{sg}=0.5`$ in Fig. 13. We use $`I_{\mathrm{ext}}=0.36,0.35975,\mathrm{},0.3425`$. The voltage of the vertical junction $`V_5`$ is approximately $`0.55V_g`$ in this small range of currents. Correspondingly, the voltage of the top horizontal junction is $`V_{5T}=V_5/2=0.275V_g`$. In our normalization, the fundamental frequency of the type B breather is then $`\omega =V_{5T}/\mathrm{\Gamma }=1.53`$. Most multipliers lie on a circle of radius $`e^{\pi \mathrm{\Gamma }g_{sg}/\omega }0.83`$ and this can be verified in Fig. 14(b).

In Fig. 14(d) we decrease the current and show the value of the voltage and the modulus of the Floquet multipliers (only for periodic solutions). This is the typically observed bifurcation scenario for small $`\lambda `$. It seems that for small $`\lambda `$ and underdamped junctions, this instability introduces more frequencies in the solution. When the periodic type B breather losses stability at $`I=0.3425`$, the solution becomes a quasi-periodic type B breather similar to the one shown in Fig. 23. This quasi-periodic type B solution persists up to a current value of 0.33775 when the array jumps to a periodic type A solution.

For large $`\lambda `$, however, we usually observe a period-doubling bifurcation where a multiplier crosses the unit circle at $`1`$, though the behavior also depends on the damping.

Our simulations in Fig. 13 show that the array sometimes switches to the superconducting state when the type B solution ceases to exist. The bifurcation scenario for the different subgap resistances is similar to Fig. 14(b).

At $`g_{sg}=0.5`$, on the other hand, when the single-site type B breather becomes unstable, the array switches to a single-site A breather. This type A breather exists for a range of currents. Fig. 14(a) shows the Floquet multipliers for a type A breather at current values: $`0.27,0.26975,\mathrm{},0.2495`$. The voltage of the vertical junction $`V_5`$ is approximately $`0.33V_g`$ in this range of current. In our normalization, the fundamental frequency of the type A breather is then $`\omega =V_5/\mathrm{\Gamma }=1.81`$. Most multipliers lie on a circle of radius $`e^{\pi \mathrm{\Gamma }g_{sg}/\omega }0.85`$ and this can be verified in Fig. 14(a).

In Fig. 14(c) we decrease the current and show the value of the voltage and the modulus of the Floquet multipliers (only for periodic solutions). Below $`I=0.2495`$ the DB is unstable and the solution switches to the superconducting state.

The selection rules between different $`m`$-site type B solutions shown by the experiments is an important feature that can not be explained in the framework of the DC circuit model and also is not well predicted by our simulations of the array. Once the breather solution reaches a critical value of the current the system chooses between many different attractors that have complex boundaries. In this process randomness and thermal fluctuations play a role as it is shown by the different switching patterns under similar experimental conditions. In our experiments we usually observe transitions from a $`m`$-site type B solution to a $`m+n`$-site type B. However sometimes we find jumps from type B to type A in the experiments (Fig. 9) The simulations (like in Fig. 13) show more frequently transitions from a type B solution to a type A solution and then, by decreasing the current, to the superconducting state.

To try to understand this switching process of the cascade of $`m`$-site type B breathers, we have introduced more elaborated models to (i) include the bias circuit in our simulations, (ii) study the effect of external fields, (iii) introduce a full-inductance matrix formalism, (iv) take into account the nonlinear character of the junction resistance, (v) include disorder randomizing the junctions critical currents and (vi) include thermal effects by adding a noise term to the junction currents. However, none of these simulations reliably predicts the observed cascade of $`m`$-site type B breathers found in our experiments as the current is decreased.

The bias resistors can be added to the model by rewriting the KVL and KCL for the new circuit. We have also included an external magnetic field. However at small values of $`\lambda `$ the field is expelled from the array because of strong inductive effects and does not affect the IV curves. This was also found experimentally in Fig. 10.

It also relatively straightforward to use a full-inductance matrix. We find that the additions of extra coupling changes significantly the decay of the fields within the array. For instance, due to the non-local coupling the breathers are no longer exponentially localized, and the field decay is algebraic away from the breather. Nonetheless, the IV remains relatively unaltered.

Finally, thermal effects may play an important role, especially if the attractors have complex boundaries. We use the standard Langevin approach and replace the resistor of the RCSJ model by a noiseless resistor in parallel with a Johnson current noise source,

$$C_j\dot{v}_j+\frac{v_j}{R_j}+I_c^j\mathrm{sin}\phi _j=I_{\mathrm{ext}}+I_j^N$$
(24)
where $`<I_j^N(t)I_k^N(t^{})>=(2kT/R_j)\delta (tt^{})\delta _{jk}`$. This results in the usual current noise spectrum density $`S_j=2kT/R_j`$. We normalize our equations as in Sec. 2.1 and $`𝒩(\phi _j)=\ddot{\phi }_j+\mathrm{\Gamma }\dot{\phi }_j+\mathrm{sin}\phi _j+i_j^n`$ and the spectrum of $`i_j^n`$ is $`S_j=2kTh_j\mathrm{\Gamma }/E_J`$ where the Josephson energy $`E_J=(\mathrm{\Phi }_0/2\pi )I_{cv}`$ and $`h_j=I_c^j/I_{cv}`$.

Our dimensionless temperature is then $`\stackrel{~}{T}=kT/E_J`$. At $`4.2\mathrm{K}`$ our vertical junctions have $`I_{cv}=345\mu \mathrm{A}`$ so that $`E_J=8.2\times 10^3\mathrm{K}`$. With this normalization $`4.2\mathrm{K}`$ equals $`\stackrel{~}{T}5.1\times 10^4`$. Similarly, for the maximum temperature we observe breathers, $`6.7\mathrm{K}`$, $`\stackrel{~}{T}1.3\times 10^3`$.

Our simulations were done at $`h=0.25`$ and $`\lambda =0.05`$. At these values of the parameters, we generally find in our simulations with Johnson noise that when a type B breather becomes unstable it jumps to a type A breather. Sometimes when a type B breather becomes unstable an $`m`$-site B breather will form as in our experiments, but this is much less common. However, we have found that at other values of the parameters (typically larger $`h`$ and $`\lambda `$) different switching scenarios occur. We have run several simulation that include thermal noise, the full inductance matrix, bias resistors, and possible stray fields. All these effects result in a very detailed model for the dynamics of the array and different breather solutions are found. Also, adding inhomogeneities, like a distribution of junction critical currents, or a magnetic field does not seem to change the fact that in the array simulations when most type B breather solutions become unstable a type A breather solution is formed.

To investigate the experimentally observed cascade of $`m`$-site type B solutions we study a toy model that allows for only type-B breathers. The simplest model with this characteristic is a variation of the standard model where an up-down symmetry for the superconducting phases is assumed, thus $`\phi _j^t(t)=\phi _j^b(t)`$ (this was also assumed in  and ). We remove the bottom branches of the full ladder, since their dynamic is by definition identical to the top horizontal junctions. Also, fluxoid quantization, Eq. (4), is modified since $`\phi _j^t=\phi _j^b`$,

$$\phi _j^v\phi _{j+1}^v2\phi _j^t+2\pi f+\frac{1}{\lambda }i_j^m=0.$$
(25)
This along with KCL yields a nonlinear coupled system that can be obtained from Eq. (6) defining $`\phi _j^b=0`$ and changing $`\phi _j^t2\phi _j^t`$ on the right hand sides.

Fig. 15 shows a simulation of the up-down symmetric ladder with $`\mathrm{\Gamma }=0.07`$ and $`\lambda =0.05`$. The curves are the average voltages of each vertical junctions. We start at $`I_{\mathrm{ext}}=0.7`$ with a single rotating junction. As we decrease the current, the single-site breather becomes unstable at $`I_{\mathrm{ext}}=0.32`$ and a 4-site breather is created. This solution persists until $`I_{\mathrm{ext}}=0.23`$ when a 6-site breather is formed and finally at $`I_{\mathrm{ext}}=0.15`$ the arrays switches back to the zero-voltage state. These current values are estimated well from our circuit model Eq. (16) and the discrepancies are due to the effect of the inductances.

The up-down symmetric ladder only allows up-down symmetric solutions such as a type B breather. When the single-site breathers becomes unstable the array can no longer jump to an A type breather. It appears that this constraint is enough to allow the formation of $`m`$-site breather when the single-site breather becomes unstable. This toy model, thus allows for the study of the switching seen in the experiments, though it is not clear the physical relation of this model to the experiments.

## 6 Linear Analysis

### 6.1 Resonances in a Ladder

Before we embark on an analytical study of the different localized solutions in the ladder we first need to understand the basic linearized excitations that can occur.

An important characteristic frequency of the ladder occurs when the frequency of a junction resonates with the lattice eigenmodes. To calculate the resonant frequencies, we linearize Eq. (6) around a solution. Every breather solution is approximately up-down symmetric far from the rotating junctions. Therefore we make the approximation $`\phi _j^t=\phi _j^b`$ and let $`\phi _j^t=\phi _0^t+\delta \phi _j^t`$ and $`\phi _j^v=\phi _0^v+\delta \phi _j^v`$. The resulting linear equations are

$`h(\delta \ddot{\phi }_j^t+\mathrm{\Gamma }\delta \dot{\phi }_j^t+\mathrm{cos}(\phi _0^t)\delta \phi _j^t)`$ $`=`$ $`\lambda (\delta \phi _{j+1}^v+\delta \phi _j^v2\delta \phi _j^t)`$ (26)
$`\delta \ddot{\phi }_j^v+\mathrm{\Gamma }\delta \dot{\phi }_j^v+\mathrm{cos}(\phi _0^v)\delta \phi _j^v`$ $`=`$ $`\lambda (\delta \phi _{j+1}^v2\delta \phi _j^v+\delta \phi _{j1}^v+2\delta \phi _j^t2\delta \phi _{j1}^t)`$
If the ladder is in the uniform whirling state then the approximate solution is $`\phi _0^t=0`$ and $`\phi _0^v=\omega t+zj`$ where $`z`$ is the wavelength of the vortex train. To calculate the dispersion relation we let $`\delta \phi _j^t=ϵ^te^{i(zj+\omega t)}`$ and $`\delta \phi _j^v=ϵ^ve^{i(zj+\omega t)}`$. We substitute in Eq. (26) and take the $`\mathrm{\Gamma }=0`$ limit since our junctions are underdamped and it can be shown that for $`\mathrm{\Gamma }<1`$, the damping terms are only small corrections to the frequency. The $`\mathrm{cos}(\phi _0^v)\delta \phi _j^v`$ results in terms of with coefficients of $`e^{i2\omega t}`$ and $`e^0`$ and all other terms have coefficients of $`e^{i\omega t}`$. We only keep terms of $`e^{i\omega t}`$. The resulting matrix equation is

$$\left[\begin{array}{cc}\omega ^2+1+\frac{2\lambda }{h}& \frac{\lambda }{h}(e^{iz}1)\\ 2\lambda (e^{iz}1)& \omega ^2+4\lambda \mathrm{sin}^2(\frac{z}{2})\end{array}\right]\left[\begin{array}{cc}ϵ_t& \\ ϵ_v& \end{array}\right]=\left[\begin{array}{cc}0& \\ 0& \end{array}\right].$$
(27)
The solution to Eq. (27) is $`\omega ^2=F\pm \sqrt{F^2G}`$ where $`F=(1+2\lambda /h+4\lambda \mathrm{sin}(z/2)^2)/2`$ and $`G=4\lambda \mathrm{sin}(z/2)^2`$. From physical grounds we expect the wavelength to be well approximated by $`z=2\pi f`$, i.e. the average distribution of vortices in the array. The resulting resonant frequencies are important when studying properties of moving vortices in the ladder. Then, there is a traveling wave of vortices with density $`f`$ that can resonate with the lattice modes.

Our breathers solutions are clearly not uniform. So the above result, while instructive, is of limited value. We can think instead of a solution that is valid far from the localized breather core. Far from the breather core, the solution to first order is

$`\phi _0^t`$ $`=`$ $`\pi f`$
$`\phi _0^v`$ $`=`$ $`\mathrm{sin}^1i_{\mathrm{ext}}`$ (28)
and this satisfies the ladder equations, Eq. (6). Since the core junctions in the breather are rotating, they will induce librations in all of the junctions in the array . After substituting the first order solution into Eq. (26) and expanding the perturbations as plane waves, we are left with the following matrix that must have a zero determinant to allow for nontrivial solutions:

$$\left[\begin{array}{cc}\omega ^2+i\mathrm{\Gamma }\omega +\mathrm{cos}(\pi f)+\frac{2\lambda }{h}& \frac{\lambda }{h}(e^{iz}1)\\ 2\lambda (e^{iz}1)& \omega ^2+i\mathrm{\Gamma }\omega +p+4\lambda \mathrm{sin}^2(\frac{z}{2})\end{array}\right]\left[\begin{array}{cc}ϵ_t& \\ ϵ_v& \end{array}\right]=\left[\begin{array}{cc}0& \\ 0& \end{array}\right]$$
(29)
where $`p=\sqrt{1(i_{\mathrm{ext}})^2}`$. We set $`f=0`$ as in the simulations and $`z`$ represents the wavelength of the perturbations. This dispersion relation describes the linearized frequencies that can resonate with a breather. In the case where $`p`$ is approximately 1 and $`\mathrm{\Gamma }=0`$ the eigenvalues have a particular simple form

$`\omega _+=`$ $`\sqrt{1+{\displaystyle \frac{2\lambda }{h}}+4\lambda \mathrm{sin}^2\left(z/2\right)}`$
$`\omega _{}=`$ $`1`$ (30)
we note that $`\omega _+>\omega _{}`$ and that $`\omega _{}`$ is associated with the $`L_JC`$ resonance of the ladder and $`\omega _+`$ with the $`L_sC`$ resonance. We will show that this is a good approximation for our simulations.

A more generalized resonance condition can be calculated by using a harmonic balance technique. For the added complexity of the harmonic balance analysis, it only yields a result that is very similar to the linearized calculations. However, the harmonic balance approach can also be used to calculate the IV curve.

### 6.2 Decay length

To calculate a resonance frequency we substitute into (30) the appropriate wavelength. For instance, the resonance of largest frequency occurs when $`z=\pi `$. However, whenever our frequency or wavelength falls outside of the dispersion relation then any linearized excitations must decay. In general we should let our wavelength be a complex number $`z=k+i\xi `$ where $`\xi `$ represents the decay length.

We first define $`\mathrm{sin}^2(z/2)=x+iy`$ and substitute into the determinant of matrix Eq. (29). We expand terms and solve for $`x`$ and $`y`$,

$$\left[\begin{array}{cc}4\lambda (\omega ^21)& 4\mathrm{\Gamma }\lambda \omega \\ 4\mathrm{\Gamma }\lambda \omega & 4\lambda (\omega ^21)\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}\omega ^4(1+p+2\frac{\lambda }{h}+\mathrm{\Gamma }^2)\omega ^2+2\frac{\lambda }{h}p+p\\ 2\mathrm{\Gamma }\omega ^3\mathrm{\Gamma }(1+p+2\frac{\lambda }{h})\omega \end{array}\right].$$
(31)
Here, the first row are the real components and the second row are the imaginary ones. This equation can be solved for $`x`$ and $`y`$,

$`x`$ $`=`$ $`{\displaystyle \frac{h\omega ^6+[\mathrm{\Gamma }^2h2\lambda (2+p)h]\omega ^4}{4h\lambda [\omega ^4+(\mathrm{\Gamma }^22)\omega ^2+1]}}`$ (32)
$`+{\displaystyle \frac{[h(1+2p)+2\lambda (1+p)\mathrm{\Gamma }^2hp2\mathrm{\Gamma }^2\lambda ]\omega ^2[2\lambda +h]p}{4h\lambda [\omega ^4+(\mathrm{\Gamma }^22)\omega ^2+1]}}`$
and

$$y=\frac{\mathrm{\Gamma }\omega }{4\lambda }+\frac{\mathrm{\Gamma }\omega (p1)}{2h[\omega ^4+(\mathrm{\Gamma }^22)\omega ^2+1]}.$$
(33)
The solution simplifies drastically when $`p=1`$,

$`x`$ $`=`$ $`{\displaystyle \frac{h\omega ^22\lambda h}{4\lambda h}}`$
$`y`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }\omega }{4\lambda }}.`$ (34)
With these formula we can calculate the decay length. We use the fact that $`\mathrm{cos}(k+i\xi )=\mathrm{cos}k\mathrm{cosh}\xi i\mathrm{sin}k\mathrm{sinh}\xi `$. Then,

$`x=`$ $`\mathrm{}\left\{\mathrm{sin}^2(z/2)\right\}=`$ $`{\displaystyle \frac{1}{2}}(1\mathrm{cos}k\mathrm{cosh}\xi )`$
$`y=`$ $`\mathrm{}\left\{\mathrm{sin}^2(z/2)\right\}=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}k\mathrm{sinh}\xi .`$ (35)
We can now solve for the $`k`$ and $`\xi `$. The solution is not simple in the general case but when $`y=0`$, it simplifies considerably. When $`k`$ is zero,

$$\xi =\mathrm{cosh}^1|12x|,$$
(36)
and the solutions decay exponentially. When $`\xi =0`$ we recover the normal modes of the system and the frequencies are given by the dispersion relation Eq. (30).

When $`y0`$ the solutions decay exponentially for all frequencies. This is due to the damping and is crucial for the existence of the resonant breather and chaotic breather solutions found in Sec. 7.2 and Sec. 7.3.

Fortunately, our ladders are underdamped so the $`\mathrm{\Gamma }0`$ limit is appropriate. Then, from Eq. (32) and Eq. (36),

$$\xi =\mathrm{cosh}^1\left|\frac{h\omega ^4+(2\lambda h+h(p+1)+2\lambda )\omega ^22\lambda h2\lambda pph}{2\lambda h(\omega ^21)}\right|$$
(37)
As $`p1`$, this simplifies to

$$\xi =\mathrm{cosh}^1\left|\frac{2\lambda (h+1)+h}{2\lambda h}\frac{\omega ^2}{2\lambda }\right|.$$
(38)
In the opposite limit of $`\mathrm{\Gamma }\mathrm{}`$, $`1/\xi `$ approaches $`0`$.

Eq. (37) gives the result for the decay length of linearized excitations in the ladder when $`\mathrm{\Gamma }=0`$. Another limit where the equations simplify occurs when $`\omega 0`$.

In the limit of $`\omega 0`$ we get

$$\xi =\mathrm{cosh}^1\left(\frac{2\lambda h+p(2\lambda +h)}{2\lambda h}\right),$$
(39)
and it can be verified that this result is independent of $`\mathrm{\Gamma }`$.

We have taken the $`\omega 0`$ and the $`\mathrm{\Gamma }0`$ limits. There is one other important limit when inductances can be neglected ($`\lambda \mathrm{}`$). Then, with $`\lambda \mathrm{}`$ and $`p1`$, the decay length is

$$\xi =\mathrm{cosh}^1\left(\frac{h+1}{h}\right)$$
(40)
and it can also be verified that this is independent of $`\mathrm{\Gamma }`$. Also as $`\lambda 0`$, $`1/\xi `$ approaches $`0`$.

The decay length given by Eq. (39) describes the decay of the DC flux in the array. For instance, using the parameters in Fig. 4 we calculate that $`1/\xi =0.32`$ and this agrees with the simulations. This results also implies that the decay is exponential and that this exponential localization has an upper bound. The decay length $`1/\xi `$ is always less than $`1/\mathrm{cosh}^1\{(h+p)/h\}`$ if $`\lambda h`$ or $`1/\mathrm{cosh}^1\{(2\lambda +p)/2\lambda \}`$ if $`h\lambda `$.

## 7 Numerical and analytical study of single-site DB solutions

### 7.1 Regions of existence of single-site DB solutions

In section 5 we numerically simulated the behavior of the DB solutions found in the experiments. Such experiments were done at moderate to small values of the damping, anisotropy and penetration depth ($`\mathrm{\Gamma }0.1`$, $`h0.25`$ and $`\lambda 0.05`$). It is also important to study the existence and behavior of the localized solutions at other values of the parameters and to estimate the critical values at which DB solutions destabilize. Varying the temperature we studied experimentally the behavior of the solutions at different values of the damping $`\mathrm{\Gamma }`$. The results were presented in Fig. 11, where $`h0.25`$ and $`\lambda 0.05`$. For these values of the parameters good agreement was found with the theoretical predictions.

The equations found in section 4 allow for a calculation of the IV curves and the maximum and minimum values of external currents supporting DB. For $`R_bR_h`$ and single-site breather solutions we find

$$\frac{I_a}{N}=\left(1+\frac{2h}{s}\right)\frac{V_v}{R_v}$$
(41)
for the IV curves and

$`{\displaystyle \frac{I_{}}{NI_{cv}}}`$ $`=`$ $`\left(2h+s\right){\displaystyle \frac{4}{\pi }}\mathrm{\Gamma },`$
$`{\displaystyle \frac{I_+}{NI_{cv}}}`$ $`=`$ $`{\displaystyle \frac{2h+s}{3h+s}}`$ (42)
for the currents.

Figs. 16 and 17 show the predictions given by the circuit model. The size of the existence regions decrease rapidly when the damping or the anisotropy increase. On the other hand, if the damping is small enough we should find localized solutions at large values of $`h`$. Also, the existence region is larger for type A solutions.

This simple model, however, does not account for any dependence of the curves with the parameter $`\lambda `$. This is an important limitation of the model and we have been unable to develop a more complete, yet still simplified, approach which incorporates $`\lambda `$. We have confirmed in the numerical simulations that $`\lambda `$ affects our predictions in two important ways. First, it affects the value of the array retrapping current. The value used in our circuit models has been calculated from a single junction and should be corrected by $`\lambda `$ in the case of the array. Indeed, some of the curves we will show below can be fitted assuming a simple linear dependence of the retrapping current with $`\lambda `$. Second, as studied in the previous section, it governs the values of the voltage at which resonances between the breather and the normal modes of the array play an important role. Roughly speaking, the resonances split the diagrams in two different regions: The small and the large $`\lambda `$ regions. When $`\lambda `$ is small, the resonance frequency is smaller than the DB frequency, and when $`\lambda `$ is large the resonance frequency is larger. Thereby, complications of damped resonances between the DB and the lattice eigenmodes are avoided in these limits.

Far from the resonance values the effect of $`\lambda `$ is a small correction to our IV curves. This is shown by the numerical simulations. See for instance Figs. 21 and 26, where IV curves numerically integrated agree quite well with the predictions of the model, Eq. (41).

We have done numerical simulations based on Eqs. (6) with $`f=0`$ in order to study the $`\lambda `$ dependence of the breather existence region. The results are presented in Figs. 18 through 20. In these diagrams we show the maximum and minimum values of the parameters for which a localized solution has been numerically found. As we will see, in some cases the characterization of the solutions inside the existence regions is quite complex in which several resonances and transitions between periodic and aperiodic states appear. In the figures we have also marked the values of the parameters at which the experiments were done, all far from these problematic areas and belonging to a region of the diagram where both type A and type B breathers are predicted to exist.

The data were calculated by integrating the governing equations with a small quantity of noise. We start with a type B rotobreather and $`h0.001`$. As we increase $`h`$, type B solutions become unstable and the solution evolves to a type A rotobreather. As we further increase $`h`$ this rotobreather becomes unstable and the system usually jumps to either a superconducting or a whirling state. To verify that our method is accurate, we have calculated Floquet multipliers for periodic rotobreather solutions and found results consistent with those shown.

Fig. 18 shows the existence regions in the anisotropy versus damping plane when $`\lambda =0.04`$ and $`I=0.6`$. When $`h`$ is large the agreement with the predictions of Eq. (42) is good. There is an abrupt deviation of the curve for type B solutions at small values of $`h`$ and a region at larger $`\mathrm{\Gamma }`$ where new type-B solutions appear. Vertical lines correspond $`\mathrm{\Gamma }=0.08`$ and $`\mathrm{\Gamma }=0.2`$, studied in Fig. 20. The asterisk corresponds to the experimental parameters.

Fig. 19 shows the existence regions in the current versus $`\lambda `$ plane when $`\mathrm{\Gamma }=0.08`$ and $`h=0.25`$. These are the values of $`h`$ and $`\mathrm{\Gamma }`$ in our experiments. We can see that Eq. (42) gives a good estimation of the maximum and minimum values of the current for localized solutions, except for the case of the minimum current with a moderate $`\lambda `$. In this case, resonances and other dynamical effects cause a substantial deviation from our simple model.

Fig. 20 shows the diagram in the anisotropy versus $`\lambda `$ plane for $`I=0.6`$ and $`\mathrm{\Gamma }=0.2`$ (left) and $`0.08`$ (right). As expected from Fig. 16 localized solutions exist at larger values of $`h`$ when $`\mathrm{\Gamma }`$ is smaller. Unexpectedly, at $`\mathrm{\Gamma }=0.2`$ and small values of $`\lambda `$ type-B solutions exist only for small values of $`h`$. This behavior is also described in Fig. 18. The asterisk in the $`\mathrm{\Gamma }=0.08`$ figure approximately correspond to the value of the parameter where our experiments were done.

Clearly, the existence diagrams are complex and further research is necessary to fully understand them. However, these diagrams show that DB solutions might be understood within our simple model at limiting values of small and large $`\lambda `$. We will use these limits in the numerical and analytical analysis of type B and type A solutions in the following sections.

### 7.2 Type B breathers

#### 7.2.1 Simulations

Fig. 21 shows typical simulated IV’s at different $`\lambda `$’s for single-site type B breathers. The horizontal dashed line is the minimum retrapping current expected from Eq. (16) when $`R_bR_h`$ in the circuit model. The vertical dashed lines are the two resonant voltages of the ladder as approximated with Eq. (30) when $`z=\pi `$ (the largest lattice frequency). As can be seen in the IV’s, $`\omega _+`$ gives a good approximation of the location of the resonant structure. Also, $`\omega _+`$ is always larger than $`\omega _{}`$. The diagonal dashed line is the expected voltage of the fifth vertical junction calculated from the DC circuit model, Eq. (14) with $`R_bR_h`$. The voltages of the fifth vertical junction and the fifth top horizontal junction are plotted and are usually related by $`V_5=2V_{5T}`$. Here, we have used the simple RCSJ model without a subgap resistance and we have set $`\mathrm{\Gamma }=0.1`$, the experimental value measured from the subgap resistance.

Fig. 21(a) shows an IV when $`\lambda =5`$. It corresponds to the large $`\lambda `$ regime. At these parameter values the upper resonance is above any of the junction voltages. However, we still see some resonant behavior at $`I0.9`$. The solution when $`I=0.8`$ is shown in Fig. 22(left). This type B breather is fully up-down symmetric in that $`\phi _j^t(t)=\phi _j^b(t)`$. Also, the averaged voltage of the horizontal junction is always half of the voltage of the vertical junction. We will label this type of solution as B1$`\alpha `$. As we decrease the current, the breather becomes unstable at $`I=0.32`$ as predicted by the retrapping model, Eq. (16), and a type A breather forms which itself becomes unstable at $`I=0.26`$.

Fig. 21(b) shows an IV when $`\lambda =2`$. We see that the $`V_5`$ branch is now separated into two parts by the $`\omega _+`$ resonance while $`V_{5T}`$ is still below the resonance. The solution when all the voltages are below the resonance is still B1$`\alpha `$. We will label as B1$`\beta `$ the breather solution where $`V_{5T}`$ is below and $`V_5`$ is above the resonance. There is also a hysteresis loop that forms at the resonance that is not shown in the figure. We find that, as the B1$`\alpha `$, the B1$`\beta `$ solution is also up-down symmetric, but in the case of the B1$`\beta `$ there exists a phase difference in between the velocity of the vertical and the horizontal junctions.

When $`\lambda =0.8`$ we see the IV shown in Fig. 21(c). We find that we can interpret the $`\beta `$ solution as a resonant type B breather. Below the resonance, there is a small remnant of the $`\alpha `$ breather. When the $`\beta `$ solution becomes unstable but $`V_{5T}`$ is still below the resonance, there is another type B breathers which we have labeled B1$`\gamma `$. An unusual signature of this breather is that $`V_52V_{5T}`$. Further analysis shows that the $`\gamma `$ solution is an aperiodic type B breather. There is also a hysteresis loop at the resonance that, depending on the parameters, might surround the aperiodic breather. This hysteresis would make it difficult to experimentally access this attractor by only using the applied current.

We can continue to decrease $`\lambda `$. Fig. 21(d) shows an IV when $`\lambda =0.5`$. The upper resonance has now divided the voltage of $`V_{5T}`$ into two branches. Below the resonance we get the solutions described above. We see that the $`\beta `$ solution is getting “squeezed” by the $`\omega _+`$ resonance and the retrapping current. Indeed, there is a critical value of $`\lambda `$ where the type $`\beta `$ solution disappears.

Above the resonance, we find a solution labeled B1. This small $`\lambda `$ regime is the same as in the experiments. The single-site breathers measured in our ladders are of this B1 type while the other breathers described above were only found numerically. In Fig. 22(right) we see that this B1 solution is not up-down symmetric although $`V_{5T}=V_5/2=V_{5B}`$.

We have done a Poincaré section analysis of the aperiodic solution B1$`\gamma `$ at the same parameter values as Fig. 21(d) In Fig. 23 we show the value of $`\dot{\phi }`$ versus $`\phi `$ at times equal to $`t_0+n\tau `$ where $`\tau =2\pi /V_5`$. The simplest periodic type B solutions have a period $`T=2\tau `$ since horizontal junction is half the voltage of the vertical junction. The solution shown in Fig. 23 seems to be quasiperiodic with two incommensurate frequencies, one of which seems to have a period equal to $`4T`$.

There are many other type B breathers that were found numerically at other values of the parameters but are not discussed here. For instance, there is a family of solutions which is not left-right symmetric. However, we shall focus our discussion to left-right symmetric solutions that have the above characteristics.

At other values of the parameters, the IV curves show chaotic localized solutions. Fig. 24 shows Poincaré section for vertical phases in the case of a chaotic localized type B solution. The values of the parameters for this solution are $`h=0.15`$, $`\lambda =0.2`$, $`\mathrm{\Gamma }=0.2`$ and $`I=0.6`$. Such chaotic region can be located in the central part of Fig. 20.

#### 7.2.2 Analysis

In this section, we will use a harmonic balance technique to derive some analytical descriptions of our DB. When studying periodic solutions it is almost always easier to work in Fourier space. Harmonic balance is a technique where the variables are decomposed into their Fourier components and substituted back into the governing equations. This creates a large set of coupled algebraic equations. If the governing equations are linear then each resulting algebraic equation is independent and the full system is easily solved. However, nonlinearities tend to mix harmonic components. If the mixing effect is large then the resulting set of algebraic equations is usually intractable. However, in our underdamped Josephson arrays the capacitances act as filters that allow the transmission of only a few frequencies. Typically it is only one frequency. In this case, we can truncate the set of algebraic equations and a harmonic balance technique can provide useful approximations.

From simulations almost all breather solutions have left-right symmetry (and this is always true for average voltages). We can then make a transformation from the full ladder to what we call a half-ladder. A half-ladder is ladder where the breather is now on the first junction and we assume left-right symmetry.

To make a mapping of the equations from a half-ladder to a full ladder we first consider placing a breather in junction 5 of a 9-junction ladder. Current conservation at that node yields,

$$I_5^v=I_5^tI_4^t+I^e.$$
(43)
Mirror symmetry implies that $`I_5^t=I_4^t`$ thereby current conservation becomes

$$I_5^v=2I_5^t+I^e$$
(44)
If we have a half ladder, current conservation at the first node (labeled 5 for comparison) yields

$$I_5^v=I_5^t+I^e$$
(45)
To make a mapping between the circuit equations, we need to multiply $`I_5^t`$ by 2. Since $`I^t=h𝒩(\phi )`$, this is equivalent to setting $`h_{\mathrm{half}\mathrm{ladder}}=2h_{\mathrm{ladder}}`$. However, we want to maintain the same flux pattern in the half-ladder so fluxoid quantization must remain unchanged. This is simply done by setting $`\lambda _{\mathrm{half}\mathrm{ladder}}=2\lambda _{\mathrm{ladder}}`$.

Since the breather solution is highly localized we can, as a first approximation, neglect every cell except the first one. The resulting reduced system of equations can be written as

$`𝒩(\phi ^t)`$ $`=`$ $`𝒩(\phi ^b)`$ (46)
$`h_c𝒩(\phi ^t)+𝒩(\phi ^l)`$ $`=`$ $`i_{\mathrm{ext}}`$ (47)
$`h_c𝒩(\phi ^t)+𝒩(\phi ^r)`$ $`=`$ $`i_{\mathrm{ext}}`$ (48)
$`\phi ^t+\phi ^b+\phi ^l\phi ^r`$ $`=`$ $`{\displaystyle \frac{h_c}{\lambda _c}}𝒩(\phi ^t)`$ (49)
where $`h_c=2h_{\mathrm{ladder}}`$ and $`\lambda _c=2\lambda _{\mathrm{ladder}}`$.

Most type B breathers have two voltages in the system corresponding to two frequencies. Our simulations indicate that the voltage of the horizontal rotating junction can be approximately decomposed as

$$v^{\{t,b\}}=\pm \frac{\omega }{2}+a^{\{t,b\}}\mathrm{cos}\left(\frac{\omega }{2}t+\theta _a^{\{t,b\}}\right)+b^{\{t,b\}}\mathrm{cos}\left(\omega t+\theta _b^{\{t,b\}}\right)$$
(50)
Similarly for the rotating vertical junction,

$$v^l=\omega +a^l\mathrm{cos}\left(\frac{\omega }{2}t+\theta _a^l\right)+b^l\mathrm{cos}\left(\omega t+\theta _b^l\right).$$
(51)
For the librating junction

$$v^r=a^r\mathrm{cos}\left(\frac{\omega }{2}t+\theta _a^r\right)+b^r\mathrm{cos}\left(\omega t+\theta _b^r\right).$$
(52)
We integrate to get the phases and use exponential notation. Our phases are then

$`\phi ^t`$ $`=`$ $`\mathrm{}\left\{{\displaystyle \frac{\omega }{2}}t+k^t+{\displaystyle \frac{\pi }{2}}iA^te^{i\frac{\omega }{2}t}iB^te^{i\omega t}\right\}`$
$`\phi ^b`$ $`=`$ $`\mathrm{}\left\{{\displaystyle \frac{\omega }{2}}t+k^b+{\displaystyle \frac{\pi }{2}}iA^be^{i\frac{\omega }{2}t}iB^be^{i\omega t}\right\}`$
$`\phi ^l`$ $`=`$ $`\mathrm{}\left\{\omega t+k^l+{\displaystyle \frac{\pi }{2}}iA^le^{i\frac{\omega }{2}t}iB^le^{i\omega t}\right\}`$
$`\phi ^r`$ $`=`$ $`\mathrm{}\left\{k^riA^re^{i\frac{\omega }{2}t}iB^re^{i\omega t}\right\}`$ (53)
where $`A^x=(2a^x/\omega )e^{i\theta _a^x}`$ and $`B^x=(b^x/\omega )e^{i\theta _b^x}`$ with $`x=t,b,l`$ and $`r`$. By convention we will take the real parts for the actual variables. Also, the integration constant has been taken to be $`k^{\{t,b,l\}}+\pi /2`$ for the rotating junctions.

To substitute into our governing equations we must linearize the sine term. Substituting our phase ansatz into the sine term would yield an infinite Fourier series whose coefficients are Bessel functions of the amplitudes. In this sense the sine terms mix all of the harmonics. In the present case, our amplitudes are small and we can linearize $`\mathrm{cos}x=1`$ and $`\mathrm{sin}x=x`$. The sine for the top rotating phase is approximated as

$$\mathrm{sin}\phi ^t=\mathrm{}\{e^{i(\frac{\omega }{2}t+k^t)}\frac{A^t}{2}e^{ik^t}+\frac{A^t}{2}e^{i(\omega t+k^t)}\frac{B^t}{2}e^{i(\frac{\omega }{2}tk^t)}\}$$
(54)
and for the bottom phase

$$\mathrm{sin}\phi ^b=\mathrm{}\{e^{i(\frac{\omega }{2}tk^b)}+\frac{A^b}{2}e^{ik^b}\frac{A^b}{2}e^{i(\omega tk^b)}+\frac{B^b}{2}e^{i(\frac{\omega }{2}t+k^b)}\}.$$
(55)
The sine of the vertical rotating phase is approximated as

$$\mathrm{sin}\phi ^l=\mathrm{}\{e^{i(\omega t+k^l)}\frac{A^l}{2}e^{i(\frac{\omega t}{2}+k^l)}\frac{B^l}{2}e^{ik^l}\}.$$
(56)
For the librating phase

$$\mathrm{sin}\phi ^r=\mathrm{}\{ie^{ik^r}i(A^re^{i\frac{\omega }{2}t}+B^re^{i\omega t})\mathrm{cos}k^r\}.$$
(57)
We have neglected all terms of frequencies not equal to $`\pm \omega `$ or $`\pm \omega /2`$.

If we substitute our ansatz into our governing equations and expand the sine terms as indicated above, we transform the original 4 differential equations into a linear system of 20 equations with 21 unknowns. This extra degree of freedom is associated with the time translational invariance of the equations.

In principle it is possible to solve the full algebraic system by fixing one of the unknowns, but here we will just estimate the amplitude oscillations of the breather. To make headway, we will use some trends found in the simulations to reduce the number of unknowns and take the limit $`\mathrm{\Gamma }=0`$.

We first calculate $`|A^t|`$ and approximate the phases from the numerical simulations. From simulations, we find that as $`\lambda \mathrm{}`$, the phase difference of the $`\omega /2`$ harmonic between the waveforms of $`\ddot{\phi }^t`$ and $`\mathrm{sin}\phi ^t`$ is $`\pi `$. We also find that the phase difference of the $`\omega /2`$ harmonic between the waveforms of $`\phi ^t`$ and $`\phi ^l`$ is zero and between the waveforms of $`\phi ^l`$ and $`\phi ^r`$ is $`\pi `$. These phase relations can be used to reduce the number of unknowns in the algebraic system.

One solution to Eq. (46) obeys the up-down symmetry: $`B^t=B^b`$ and $`A^t=A^b`$. We can combine the $`\omega /2`$ harmonic parts of Eq. (47), Eq. (48), and Eq. (49) and solve for $`|A^t|`$. This yields

$$|A^t|=\left|\frac{2h_c\lambda _ch_c(\omega /2)^2}{(\omega /2)^2(2\lambda _c+2h_c\lambda _ch_c(\omega /2)^2)}\right|$$
(58)
Using the same technique to approximate phases using simulations for the $`\omega `$ harmonics, we find

$$|B^r|=\left|\frac{1}{\omega ^2(2+2/h_c\omega ^2/\lambda _c)}\right|$$
(59)
with $`h|B^t|=|B^r|`$. We can then find all of the amplitudes of the harmonics in terms of $`\omega `$.

As $`\lambda `$ decreases the breather enters a resonant regime. In this regime the phase relationships change from the above limit. We now find the phase difference of the $`\omega /2`$ harmonic of the waveform of $`\ddot{\phi }^t`$ and $`\mathrm{sin}\phi ^t`$ is zero, the phase difference between the waveforms of $`\phi ^t`$ and $`\phi ^l`$ is $`\pi `$, and between the waveforms of $`\phi ^l`$ and $`\phi ^r`$ is also $`\pi `$. Substituting the new ansatz in the governing equations results in the same equation as Eq. (58) but with an overall negative sign. The same occurs with the $`B^r`$ amplitude. In the resonant regime Eq. (59) just changes sign.

In the $`\lambda 0`$, the equations simplify. From Eq. (49) we see that $`𝒩(\phi ^t)`$ must tend to zero. This implies that

$$|A^t|=\frac{1}{(\omega /2)^2}.$$
(60)
This also puts strong constraints on $`A^r`$ and $`A^l`$. One possible solution is $`A^r=0`$ and $`A^l=0`$. Then to satisfy Eq. (49), $`A^t=A^b`$ contrary to the expected up-down symmetry. However, this solution is consistent with Eq. (46) and it is the stable solution we find numerically for small $`\lambda `$.

Fig. 25 shows a summary of the analytical results and how they compare to the simulations. The top graph is a simulation for a cell while the bottom graphs shows results for our ladder. Both system have $`\mathrm{\Gamma }=0.1`$. For each system, we have excited a type B breather at $`i_{\mathrm{ext}}=0.5`$ and $`\lambda =10`$. We plot the value of $`|A^t|`$ and $`\omega `$ as we decrease $`\lambda `$.

As expected, both graphs are very similar when $`\lambda _{\mathrm{cell}}=2\lambda _{\mathrm{ladder}}`$. The top solid line is from the simulated frequency of the rotating vertical junction. The bottom solid line is from the simulated frequency of the horizontal junction which is precisely half of the vertical frequency. The dashed line shows the $`\omega _+`$ resonance. Because there are two frequencies, the resonance divides the solution space into three regions that result in the three different solution found in Fig. 21 and studied analytically using harmonic balance.

The solid circles in Fig. 25 are the simulated $`A^t`$ and the dashed lines are the analytical estimate. For small $`\lambda `$ we have the B1 solution whose amplitude is given by Eq. (60) and the predicted amplitude lies on top of the simulated values. For the resonance regime of solution B1$`\beta `$ the amplitude is given by Eq. (58), and for the large $`\lambda `$ regime the amplitude is also given by Eq. (58).

For the simulations of a single cell, the harmonic balance gives a very good approximation. This is not too surprising, since the analysis was performed for a single cell. However, the bottom of Fig. 25 shows that the analysis is also valid for the full ladder as long as we transform $`\lambda `$ and $`h`$.

### 7.3 Type A breathers

Fig. 26 shows typical simulated IV’s as a function of $`\lambda `$ for type A breathers. As in the previous discussion, the horizontal dashed line is the expected retrapping current while the vertical dashed lines are the two resonant voltages from Eq. (30). Here the voltage of the fifth vertical junction and the fifth horizontal junction are the same and this makes the analysis of the solutions much simpler. There are basically just two solutions.

A harmonic balance technique can be used to calculate the amplitudes of the oscillating parts of the breather. We use a similar approach to Sec. 7.2.2 but approximate the phases with one frequency instead of two. Our phases are then

$`\phi ^t`$ $`=`$ $`\mathrm{}\left\{\omega t+k^t+{\displaystyle \frac{\pi }{2}}iA^te^{i\omega t}\right\}`$
$`\phi ^b`$ $`=`$ $`\mathrm{}\left\{k^biA^be^{i\omega t}\right\}`$
$`\phi ^l`$ $`=`$ $`\mathrm{}\left\{\omega t+k^l+{\displaystyle \frac{\pi }{2}}iA^le^{i\omega t}\right\}`$
$`\phi ^r`$ $`=`$ $`\mathrm{}\left\{k^riA^re^{i\omega t}\right\}.`$ (61)
When $`\omega _+`$ is above the junction voltages, as in Fig. 26(a) the solution is up-down symmetric with respect to the amplitudes: $`A^t=A^b`$, and $`A^rA^t`$. This means that all of the core junctions of the breather have a small oscillating amplitude. As $`\lambda `$ decreases the breather can resonates with the lattice eigenmodes as shown in Fig. 26(b). When $`\lambda `$ is so small that $`\omega _+`$ is below the breather frequencies we find a solution with $`A^b0`$ and $`A^r0`$.

The analysis of the amplitudes is straight forward and similar to what is done in Sec. 7.2.2. We substitute the ansatz in the governing equations and keep only DC and $`\omega `$ harmonics. We then use simulations to approximate some of the phase relations and thereby reduce the number of equations.

In the $`\lambda =0`$ limit, we have from Eq. (49) that $`𝒩(\phi ^t)=0`$. If we set $`\mathrm{\Gamma }=0`$, this implies that

$$|A^t|=\frac{1}{\omega ^2}$$
(62)
In the opposite limit, we can again use Eq. (49). We substitute the ansatz found in the simulations of $`A^t=A^b`$ and approximate $`A^r`$ as $`A^t`$, then

$$|A^t|=\left|\frac{h/\lambda }{3h\omega ^2/\lambda }\right|$$
(63)
Fig. 27 shows a comparison the analytical results and how they compare to the simulations. The top graph is a simulation for a cell while the bottom graphs shows results for the ladder. Both system have $`\mathrm{\Gamma }=0.1`$. For each system, we have excited a type A breather at $`i_{\mathrm{ext}}=0.4`$ and $`\lambda =10`$. We plot the value of $`|A^t|`$ and $`\omega `$ as we decrease $`\lambda `$.

As expected, both graphs are very similar when $`\lambda _{\mathrm{cell}}=2\lambda _{\mathrm{ladder}}`$. The solid line is from the simulated frequency of the rotating vertical junction. The dashed line shows the $`\omega _+`$ resonance and we see that there are two types of solutions: Breathers with frequencies above $`\omega _+`$, and breathers with frequencies below $`\omega _+`$.

The solid circles in Fig. 27 are the simulated $`A^t`$ and the dashed lines are the analytical estimate. For small $`\lambda `$ the amplitude is given by Eq. (62) and the predicted amplitude match the simulated values. For the large $`\lambda `$ regime the amplitude is given by Eq. (63) and the approximation breaks down at the resonance. In any case, we see that comparing a breather in one cell to the breather in a ladder works quite well.

### 7.4 Vortex patterns in breathers

The DB that we observe in the ladder can also be understood as bound vortex states. In Eq. (5) we define the vorticity $`n_j`$ in a cell as

$$n_j=\frac{1}{2\pi }\{[\phi _j^v][\phi _{j+1}^v][\phi _j^t]+[\phi _j^b]\}+f_j^{\mathrm{ind}}$$
(64)
where $`[\phi ]`$ represents the phases modulus $`2\pi `$, and $`f_j^{\mathrm{ind}}=i_j^m/2\pi \lambda `$. As we will show, $`n_j`$ changes in time differently for the three analyzed solutions. However, we want to stress that the magnetic field flux, which is due to the mesh currents in the ladder, librate around a mean value as shown in Fig. 4 and do not circulate. Only the fluxoid, or vorticity, which is the quantized quantity, moves from cell to cell. This behavior is different for the cases of type B or type A solutions and in the case of type B solutions for small or large values of $`\lambda `$.

#### 7.4.1 Large $`\lambda `$ type B solutions

In Fig. 28 we show a sequence of time snapshots of the ladder for a type $`B`$ breather in the large $`\lambda `$ limit when $`f=0`$. There is one vertical junction rotating and four horizontal rotating junctions. The voltage of the vertical junction is twice of the horizontal junction. For every full rotation of a horizontal junction, the vertical does two full rotations. The solid circle is a vortex and the “$`\times `$” is an anti-vortex. In this large $`\lambda `$ case, the solution is completely up-down symmetric (see Fig. 22). Thus $`\phi _j^b(t)=\phi _j^t(t)`$.

At the initial condition (a) there are no vortices in the ladder. As soon as the vertical junction goes over $`\pi `$ then a vortex-antivortex pair is created as shown in (b). This pair remains pinned in the ladder even though the applied current is large enough to depin isolated vortices. When the horizontal junctions go over $`\pi `$ then they create two vortices or antivortices as shown in (d). The result is as if the vortex and antivortex pair change cells. This solutions remains until the vertical junction rotates through $`\pi `$. Then, another vortex-antivortex pair is created that annihilates the pair of opposite polarity that was in the ladder. Now, there are no vortices in the ladder and the sequence repeats itself. The period double nature of the solution is evident. The vertical junction both creates a vortex-antivortex pair, and also annihilates the vortex-antivortex pair created by the horizontal junctions.

#### 7.4.2 Small $`\lambda `$ type B solutions

Fig. 29 shows a type B breather in the small $`\lambda `$ limit at $`f=0`$. This solution is similar to Fig. 28 but we no longer have up-down symmetry. Instead, the top and bottom rotating horizontal junctions have a $`\pi `$ phase difference (see Fig. 22). We again start when the ladder has no vortices in (a). The vertical junction goes over $`\pi `$ and creates a vortex-antivortex pair. Then in (c) the bottom horizontal junctions go over $`\pi `$ and annihilate the vortex-antivortex pair. After some time, the vertical junction again creates a vortex-antivortex pair that disappears when the top horizontal junctions go over $`\pi `$ and annihilate the vortex-antivortex pair.

In the small $`\lambda `$ limit vortex-vortex interactions are large. Therefore, vortices prefer not to enter through the top and bottom horizontal junctions at the same time. In order to minimize this repulsive interaction, vortices and anti-vortices enter the cells at different times when $`\lambda `$ is small and the up-down symmetry of the solution is broken.

#### 7.4.3 Type A solutions

Fig. 30 shows a type A breather for $`f=0`$. The solution in this case is easier to interpret that the type B solutions since every junction rotates at the same voltage. First, the vertical junction goes over $`\pi `$ and creates a vortex-antivortex. The horizontal junctions then rotate over $`\pi `$ and annihilate the vortex-antivortex pair as shown in (c).

If we apply a field, then a phase difference develops between the horizontal junctions. For instance, in the appropriate field the right horizontal junction can go over $`\pi `$ first, thereby injecting a vortex to the right cell of the rotating vertical junction. Then, the vertical junction goes over $`\pi `$ creating a vortex-antivortex pair. The antivortex is created in the right cell and the vortex on the left. Since there was already a vortex on the right cell, the anti-vortex annihilates this vortex. The resulting state is as if the vortex moved from the right cell to the left cell. Then, the left horizontal junction goes over $`\pi `$ creating an anti-vortex in the left cell and thereby annihilating the vortex. This results in the ladder having no vortices. The full process can be interpreted as a rotating vortex around the common node of the rotating junctions. Therefore, type A breathers in a magnetic field are rotating vortices.

## 8 Discussion and Conclusions

Nonlinearity can cause localization in otherwise perfectly periodic systems. We have fabricated an experimental system to study these types of localized solutions called discrete breathers (DB). We have experimentally detected different one-site and multi-site DB states in a superconducting Josephson ladder network. DB can be excited and annihilated in a controlled way in the ladder by using a local bias current. Then, by varying the external current and temperature we have studied the domain of existence and the instability mechanisms of these localized solutions. Experimentally we find two families of DB: type B and type A solutions. In the case of type B, four of the horizontal junctions in the array are in the resistive state, while for type A only two horizontal junctions are in the resistive state. We find that some of the type B and all the type A breathers do not obey the natural up-down symmetry of the ladder equations. The existence of type A solutions was predicted in  and recently has also been reported measurements showing a diversity (type B, type A and hybrid) of DB in an open JJ ladder array .

We have developed a series of models and done theoretical analysis and simulations to explain much of the behavior of the experimental data. At the simplest level, we have developed a DC circuit model for the system. In this model the junctions in the resistive state are studied as linear resistors and the junctions in the superconducting state as shorts. This model allows for evaluating the effect of the bias circuit and computing IV curves. We can estimate the minimum current for which both type A and B breathers can exist by using this model and assuming that the minimum applied current is determined by a single junction retrapping mechanism. We can also predict the maximum current for which type B breathers can exist by using the junction’s critical current. All these predictions agree with the experimental results. However, the type A maximum current found in our experiments is lower than predicted and is not yet understood.

Using numerical simulations we have found the existence of many different DB attractors. We find periodic, quasiperiodic and chaotic type B solutions. Such solutions are identified with different regions in the numerically computed IV curves, and when we vary the parameters we find many different scenarios for the behavior of the solutions. The DB also excite waves in the ladder. We have calculated the frequencies of these small waves by using a linear analysis of the system equations. A resonance can occur when the DB frequency matches a frequency of the linear wave. We have numerically verified that these resonances can cause instabilities in the localized solutions, but in most of the cases the resonances are not strong enough for destroying the localization of the breather. We have also performed a harmonic balance approximation which allows the characterization of the amplitudes of type A and type B breathers.

Since DB do not exist for all the values of the model parameters, we have studied the regions of existence of localized solutions when these parameters vary. This study showed good agreement with the predictions of the simple model in many cases. In brief, DB exist at small values of the anisotropy, but the underdamped character of the junctions is crucial. In general, a smaller $`\mathrm{\Gamma }`$ results in a larger existence region. Thus, for small values of the damping DB are predicted to exist even at high values ($`>1`$) of the anisotropy parameter. However, we find a complicated dependence of the DB instabilities with respect to $`\lambda `$ which is only partly explained by resonances between DB frequencies and the lattice linear waves.

When varying $`\lambda `$ we find different physical behaviors for large and small values. At large values of $`\lambda `$, type B solutions are up-down symmetric, and at small values of $`\lambda `$ this up-down symmetry is broken. In the intermediate $`\lambda `$ regions the resonances between the breather frequency and the linear modes are important and the behavior is much more complex. It is in this region where we found chaotic localized solutions.

One important conclusion of these simulations at different values of $`\lambda `$ is the observation of resonances in the IV curves which do not destroy the DB solution. This can be understood after the analysis presented in section 6.2. For the existence of DB, in the case of a forced and damped system the non-resonance condition is not necessary, because for any frequency the decay length is different from zero.

The induced magnetic fluxes in the cells of the array librate around non-zero DC values and have different signs in opposite regions of the array with respect to the localization of the breather. To complete the physical picture of the dynamics of the DB we also analyzed the type B and type A solutions in terms of the vortex-antivortex behavior in the array. Briefly, vortex-antivortex pairs are created in the center of the array and destroyed when new vortices and antivortices enter the array through the horizontal junctions. In the case of a type A solution with an applied external magnetic field, this vortex picture corresponds to the presence of a rotating vortex.

Many of the aspects of the behavior of the solutions, including ranges of existence, were understood using simple models. However transitions between different states when decreasing the external current is not fully understood. Although different behaviors have been found, a typical experimental result is that type B breathers become unstable through a cascade of switching events to multi-site type B solutions. To understand this phenomenon we have done more detailed simulations of the dynamics of the array. These simulations are based in the RCSJ model and we have numerically studied the effect of the bias circuit, the nonlinear resistance, full inductance matrix, thermal fluctuations and inhomogeneities in the junction parameters. Although sometimes a switching to a multi-site state was found, the typical numerical result for small $`\lambda `$ and $`h=0.25`$ shows destabilization from type B to type A solutions. Another aspect is that all these transitions occurs close to the region of the gap voltage where it might be necessary to used models more sophisticated than the RCSJ one.

We have also studied a reduced model where up-down symmetry in the system is imposed by setting $`\phi _j^t(t)=\phi _j^b(t)`$. Then, transitions to multi-site states are typically found. However this model is unable to explain the existence of type A solutions and the transitions from type B to type A solutions that have been reported here and in . This approach using a reduced model has been used in  and  to study the dynamics of the array.

DB and vortices are two examples of localized excitations in a Josephson array. Future work will explore the interplay between DB and vortices (kinks or solitons) solutions which can coexist in the array . Another interesting direction is a experimental study DB solutions in other geometries such as triangular arrays, and the excitation of breathers by AC driving currents.

After this article was accepted we learnt of  which reports on numerical work of switching mechanisms in rotobreathers at $`h=0.44`$ and $`\lambda 0.21`$.

This work was supported by NSF (DMR-9610042 and DMR-9988832), DGES (PB98-1592), EU (HPRN-CT-1999-00163, LOCNET) and and the Commission for Cultural, Educational and Scientific Exchange between the United States of America and Spain. JJM also thanks Terry Orlando and his group for his hospitality at MIT.