{'fraction_non_alphanumeric': 0.02963769430835299, 'fraction_numerical': 0.009707857013246234, 'mean_word_length': 4.528674134760074, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Breathers in a single plaquette of Josephson junctions: existence, stability and resonances

We present a theoretical study of inhomogeneous dynamic (resistive) states in a single plaquette consisting of three Josephson junctions. These breather states are found to exist in a large range of control parameters (dc bias, anisotropy and self inductance). We perform an analytical analysis of these states, and their stability and resonance behaviour. Numerical calculations confirm our results and show that instabilities arise due to resonant interactions between the breather and electromagnetic oscillations mediated by the breather state. These instabilities manifest themselves by various nonlinearities, resonant steps and voltage jumps in the current-voltage characteristics.


I. INTRODUCTION
In recent years we witnessed a considerable advancing in the area of study of localized excitations in nonlinear lattices. 1-3 These excitations are coined discrete breathers and are time-periodic and spatially localized solutions of the underlying dynamical equations of motion. They generically appear in spatially homogeneous systems due to the interplay between discreteness and nonlinearity. Discrete breathers were detected experimentally in different systems, e.g. in weakly coupled optical waveguides, 4 in the lattice dynamics of solids, 5 in antiferromagnets 6 and in systems of interacting Josephson junctions. [7][8][9][10][11] Especially the latter example served for many years as a well controlled laboratory object to study various nonlinear phenomena. Josephson junction systems allow for an adequate theoretical description of the underlying dynamics. The originally proposed and realized ladder geometry of small interacting Josephson junctions with an external dc bias became a system of intense theoretical and experimental studies. [12][13][14][15] The dynamic localized states in Josephson junction ladders persist due to a well known property of a single Josephson junction, namely the presence of two different states. These are a superconducting state (zero voltage drop across the junction) or a resistive state (nonzero voltage drop across the junction). If the junction is underdamped, both states may coexist in a wide range of an externally applied dc bias. A well known consequence of that is a hysteretic currentvoltage (I − V ) characteristics.
A breather state in such a ladder is characterized by a few junctions being in the resistive state, while the rest of all junctions reside in the superconducting one. These inhomogeneous dynamic states lead to specific total dc voltage drops across the ladder which are used to plot I − V characteristics. This method has been successfully combined with snapshots made using the low temperature scanning laser microscopy technique (LTSM) which allows for a spatial resolution of dc voltage drops. 8,9 Note that the experimentally accessible information provides only with time averaged data, so the internal dynamics of breather states is so far not directly observable in experiments. The full dynamical picture is however more subtle. E.g. the superconducting junctions display nonzero ac voltage drops. The amplitude of these ac voltages decays to zero with increasing spatial distance from the resistive junctions. Moreover, a breather state may be tuned into resonance with other dynamical modes of the system which leads to an increasing complexity of the breather dynamics. Experimental studies revealed a whole "zoo" of different breather states, instabilities, resonances and switching scenarios depending on the parameters of the ladder. [7][8][9][10][11] Similar features have been found in numerical simulations.
To clarify the role of control parameters and the cause of the primary (generic) instabilities, resonances and switchings we carry out a consistent study of a single plaquette model containing three Josephson junctions. Whereas the occurrence of breathers in an extended ladder may be interpreted as a dynamical breaking of discrete translational invariance, in the single plaquette model the analogy of the translational invariance becomes a permutational invariance. It is this permutational invariance which will be violated by dynamical states coined breathers in analogy to the ladder case. The approach of studying inhomogeneous dynamic states in reduced systems 15,16 has been successfully used in the past to get deeper understanding for different breather properties as e.g. energy radiation and quantization. 2 Such a single plaquette model is the simplest system which supports permutationally broken states in the presence of a homogeneous dc bias γ (Fig.1). It consists of two vertical Josephson junctions parallel to the bias direction and a horizontal one that is transversal to the bias. The Josephson junctions are indicated by crosses. In a time-averaged picture we expect only four different states to be realized. The first one is characterized by all junctions being in the superconducting state (SS), with zero current flowing through the horizontal junc-tion. The second one is the homogeneous whirling state (HWS) when both vertical junctions are in the resistive state. Due to permutational invariance of this state the current flowing through the superconducting horizontal junction is still zero. Note that both states are invariant under permutation, i.e. the exchange of two vertical junctions. The states of interest are the remaining two breather states. Each of these states is characterized by only one vertical junction and in addition the horizontal junction being resistive. The second vertical junction is in a superconducting state, although there is a nonzero current flowing through it. Obviously both breather states are not invariant under permutation.
In the following we will systematically study the dependence of breather state properties on the control parameters of the system (dc bias, anisotropy, self-inductance). Our paper is organized in the following way. In section II we present the model, introduce the equations of motion and provide with an approximative analytical treatment. In section III we introduce a number of numerical methods for computing breathers, obtaining their stability and following their evolution. With their help we find the domains of existence, stability and resonances of breather states in the space of control parameters and relate them to computed I − V characteristics.
φ Im Sketch of the plaquette with three Josephson junctions (marked by crosses) in the presence of uniform dc bias. Arrows indicate the direction of external current flow (dc bias γ).

A. Equations of motion
The dynamics of a single plaquette of three Josephson junctions (see Fig. 1) is determined by the time dependent Josephson phases of vertical junctions φ v 1,2 , and the horizontal junction φ h . We will consider an anisotropic system and the parameter of anisotropy is η = I ch /I cv , where I ch and I cv are respectively the critical currents of the horizontal and vertical junctions. By making use of the Resistive Shunted Junction model for each junction 17 we obtain the following set of equations The nonlinear operator N is defined as (2) Here, the unit of time is the inverse plasma frequency, the dc bias γ is normalized to the critical current of the vertical junctions I cv , and the dimensionless parameter α determines the damping in the system. In (1) we introduce the mesh current I m that in turn, is determined by the flux quantization law 18 where β L = 2πI cv L/Φ 0 is the normalized selfinductance L of the plaquette. Note here, that the voltage drop across a junction is given byφ.

B. Breather solutions
We are interested in the properties of inhomogeneous breather states. In the limit η = 0 we may exactly solve the equations of motion and implicitly obtain the breather solution in terms of the resistive state of a single vertical junction. In this case the mesh current vanishes I m = 0 and the two vertical junctions decouple, while the horizontal one follows the dynamics of the vertical junctions (below we will consider a particular breather case as the left vertical junction is in the superconducting state): Here, the function f (α, γ, t) describes the evolution of a single vertical junction in the resistive state with bias γ and damping α. It is well known that this state exists as 17 where γ r is the retrapping current of a single Josephson junction. On the other side the superconducting state of the first vertical junction (see Eq. (4)) will persist only if γ < 1. Consequently we obtain that the breather solution (4) exists for η = 0 provided γ r ≤ γ ≤ 1.
For nonzero values of η = 0 exact analytical solutions are not available. In the following we obtain approximate solutions using a dc analysis, i. e. neglecting the Josephson phase oscillations. We write the Josephson phases in the form: where Ω is the breather frequency. Substituting these expressions into (1) and (3) the mesh current reads Correspondingly the dc bias dependent breather frequency and the Josephson phase shifts are c v 1 = asin Here, we may choose c h = 0 without loss of generality. As a result we find the I −V characteristics of the Josephson junction plaquette in the presence of a breather state: (11) differs from that of the HWS which is γ = V α. The breather is characterized by a smaller voltage drop across the resistive vertical junction for a given current value γ.
As the dc bias current decreases the breather state can disappear as a solution of the dynamic equations and the system undergoes a transition to the superconducting state. By making use of the dc analysis and a standard theory of the retrapping current in a single small Josephson junction 17 we expect that the breather state disappears as the current flowing through the horizontal (resistive) junction is equal γ r . This happens for the following value of the external dc bias γ b,r : Remarkably this value of the retrapping current does not depend on the self inductance parameter β L . The dependence of γ b,r on η is shown in Fig.2.
For large values of the dc bias the breather state switches to the HWS. This switching occurs because the left vertical junction can not anymore support the superconducting state. From (9) we obtain that this switch should take place at The crucial assumption of the above carried out analysis is the neglecting of the Josephson phase oscillations. Indeed, the inhomogeneous dynamic breather state induces the electromagnetic oscillations (EOs) in the plaquette and in turn, these EOs change the simple breather solution (6). However, if the breather frequency Ω is large with respect to the characteristic frequencies of EOs, the influence of such a nonresonant interaction between the breather and EOs is small. Thus, in the nonresonant case the amplitude A of the breather induced EOs is (2η+1)Ω 2 and correspondingly, it leads to a corrected prediction of the critical dc bias value γ b,HW S where the breather switches to the HWS: The result is shown in Fig.2. Note that similar existence regions for the breathers on extended ladders have been reported in. 10 The main outcome of the presented analytical treatment is a triangular shaped region in the space of control parameters γ and η where we expect to find breather states. Notably the prediction is that for small values of damping (α ≪ 1) the breather should exist for rather large values of the anisotropy η, e.g. for α = 0.05 the anisotropy can be as large as η ∼ 5. Such parameter regions are the opposite of what has been initially anticipated, when the breather states were found by making use of the continuation arguments for the solution (4) at η = 0 to small nonzero η values. 19 In this case the horizontal junctions serve as weak links. Instead we predict the existence of breather states for cases of large η values where the vertical junctions serve as weak links, and the horizontal one is a strong link!

C. EO frequencies and resonances
Next we turn to the analysis of the resonant interaction between the breather state and small amplitude EOs. In order to obtain the characteristic frequencies of EOs we linearize the dynamic equations around the breather state (6): Here δ v 1,2 (t), δ h (t) correspond to small amplitude EOs excited in the presence of the breather state. In the weakly damped case (α ≪ 1) we may neglect the damping term and skip the time-periodic terms in (16) and (17). As a result we derive two characteristic frequencies of EOs: These frequencies depend on the anisotropy η and the bias current through the phase shift c v 1 (Eq. (9)). For large values of β L the lower frequency ω − ∼ 1 √ βL and the large frequency ω + ∼ cos(c v 1 ). In the opposite limit of small Since the breather is a time-periodic solution and the EOs describe small amplitude deviations from the breather state, we can expect that the following scenario holds. In the case of a proper frequency matching the breather may drive EOs into resonance. In turn the increase in the EO amplitudes will proceed on the expense of the breather state. In other words the possibility of a breather mediated resonance for EOs turns into a resonant interaction between the breather and the EOs. In such a case we expect in the full dynamical system various instabilities of the breather state and corresponding switchings (the voltage jumps) or resonant steps in the I − V characteristics to occur.
Let us classify the possible expected resonance cases. The first class of primary resonances may appear if where m is an integer. These primary resonances correspond to a matching of the breather frequency Ω or its multiples (higher harmonics of the breather state) with the EOs frequencies.
A second class consists of parametric resonances which are characterized by Finally we may also expect the appearance of combination resonances occurring as In such a case, which is at variance with the first two types of resonances, the breather mediates a resonant interaction between two different EO frequencies.
All listed scenarios of resonance involve the breather frequency and the EO frequencies. While the approximate expression for the breather frequency (8) is independent on β L , the EO frequencies depend sensitively on the self inductance (18). Consequently the expected positions of resonances vary with β L . Since the existence window for breathers (cf. Fig.2) is expected to not depend on β L , certain resonances may simply be shifted outside of the breather existence region and thus, become irrelevant. E.g. for small values of β L the EO frequency ω + becomes very large and we expect then only primary and parametric resonances involving ω − . If however intermediate values for β L are chosen, additional primary and parametric resonances involving ω + appear. But most important we expect combination resonances to appear, which can be realized only if the two EO frequencies, their sum or difference are of the same order as the breather frequency.

A. Description of numerical methods
Let us describe the numerical methods of monitoring and characterizing breather states. In all cases a Runge-Kutta integration method of the differential equations (1) is chosen. The first way is to fix all parameters including the bias γ, to choose some proper initial conditions and to integrate the equations for sufficiently long time until the system relaxes onto a breather attractor. Once the attractor is reached (with sufficient accuracy) all characteristics of this state may be obtained, including the average voltage drop V . After that we may perform a small step in γ and repeat the procedure. Then we find the I-V characteristics. This method is especially good when approaching an instability, since we obtain explicitly the switching to another state. We also generate time-resolved images (movies) of the full dynamical behaviour in order to visually check the realized state or monitor the process of switching. The drawback of this method is that it is not capable of explaining the nature of an observed instability. In addition it may become time consuming if the relaxation times become large.
Another way of studying the system states is to use a Newton scheme that allows to iteratively find timeperiodic states, e.g. breathers. This method is very fast and possesses a high precision. Once a solution is found, it may be characterized similarly to the way described above. The drawback here is that we can not predict the outcome of a switching if a breather state disappears or turns unstable. The method is even insensitive to the stability property of a state under consideration, so we may even miss the switching.
In order to obtain reliable information about the stability properties of breather states we perform in addition a Floquet analysis of small amplitude deviations from a breather state which we generated with the help of a Newton scheme. We consider the time evolution of small perturbations {ǫ i (0),ǫ i (0)} (i = 1, 2, 3) at t = 0 of a breather solution (φ v 1,2 , φ h ) over one period T b of the breather state. The outcome is a map of the space of small perturbations onto itself, which is characterized by a Floquet matrix M (of rank 6): This matrix is generated numerically. The breather will be linearly stable if and only if the matrix M has no eigenvalues with modulus larger than one. The quasisymplectic properties of this map impose the following conditions on its eigenvalues: if λ is an eigenvalue, then λ * , e −αT b /λ and e −αT b /λ * are also eigenvalues (here a * denotes complex conjugation of a). 20 Furthermore there is always one eigenvalue λ t = 1 which corresponds to perturbations tangent to the trajectory of the breather solution. Consequently there always exists a second eigenvalueλ t = e −αT b . These two eigenvalues will not be of further interest. The remaining four eigenvalues can be grouped into two pairs each consisting of an eigenvalue and its complex conjugated partner. These two groups correspond to the two EO frequencies (18). In the normal case of a stable breather two eigenvalues are simply characterized by Thus the four relevant eigenvalues will be located on an inner circle with radius e −αT b /2 < 1 in the complex plane. The variation of parameters system leads to a change (moving) of Floquet eigenvalues. Resonances correspond to collisions of an eigenvalue with its complex conjugated partner on the positive real axis (primary resonance), on the negative real axis (parametric resonance) or to the collision of two different eigenvalues away from the real axis (combination resonance). In all cases colliding eigenvalues will generally leave the inner circle, one moving inside and one outside. Upon further variation of some parameter the outer one may cross the unit circle, which marks the point of instability of the breather state. For small values of α the two events -collision and crossing the unit circle -take place close to each other, so that a clear connection between instabilities and resonances may be drawn. The Floquet analysis thus allows in addition to characterize the type of observed instability (resonance).
A combination of all three methods -direct numerical integration, Newton scheme and Floquet analysisturns out to be a powerful tool which provides with good insight of the dynamical processes under consideration.

B. Domains of breather existence and instability
windows.
In order to obtain the breather solution for different values of anisotropy η and correspondingly, the domains of breather existence in the parameter space (γ, η), we used the following procedure. First, we fixed the value of the external dc bias γ and found a breather solution in the form (4). After that we increased the parameter η by a small step ∆η = 10 −2 , and the previously found solution was used as an initial guess for a Newton method. The procedure was repeated for other values of dc bias γ. We varied the self inductance parameter β L while in most of our simulations the damping parameter α = 0.1 was fixed.
In the particular case of zero self inductance β L the equations (1) are simplified as φ h = φ v 1 − φ v 2 , and the system is characterized by just two degrees of freedom. We found that the domains of breather existence in the parameter space γ-η have a triangular form, and it is shown in Fig. 3.
One can also see that the existence domain expands in the region of large η as the damping parameter α decreases ( Fig. 3a and Fig. 3b). The maximum value of the anisotropy η max still allowing to obtain the breather state is estimated as η max ≃ 1 4α (see inset in Fig. 3b).
Note here, as the anisotropy is not large (η ≤ 1.5 for the system with damping α = 0.1) the dc analysis developed in the Section II (see, Fig. 2 and Eqs. (12), (13) ) gives a good approximation for numerically simulated domain of breather existence (Fig. 3). As the anisotropy increases, η > 1, the region of dc bias γ where the breather state was found, shrinks with respect to the one obtained by making use of the dc analysis.  1 (Fig. 3a) and 0.05 (Fig. 3b). The inset in Fig. 3b shows the dependence of the maximum value of anisotropy ηmax on the damping α.
In order to study the interaction of the breather state with EOs we carried out direct numerical simulations of the system of equations (1) and (3) for different values of anisotropy. Moreover, we monitored the current-voltage (I-V ) characteristics of the system decreasing (increasing) the dc bias by a step of ∆γ = 0.0005. 21 The initial values of the dc bias γ corresponded approximately to the center of the breather existence domain for the particular value of the anisotropy η.
We start with the case of a small anisotropy η = 0.1 (Fig. 4). The breather was excited at γ = 0.5 and it persists in a wide range of the dc bias. The breather state switches to the HWS as γ increases beyond some critical value, and smoothly transits to the superconducting state as the dc bias approaches the retrapping current (Eq. (12)). Note here, that for such small values of anisotropy the dc analysis (Eq. (11) and dashed line in Fig. 4) gives a good approximation practically for a whole I-V curve and it means that the resonant interaction of the breather with EOs is very weak. As we increase the value of anisotropy to η = 1, the I-V curve displays peculiar features (Fig. 5). Indeed, as the dc bias decreases the I-V curve deviates from the linear behaviour ((Eq. (11) and dashed line in Fig. 5), displays a cusp at γ = 0.34 and eventually switches to the HWS state. It is interesting that this scenario is realized if we start from rather large values of the dc bias γ ≥ 0.33. For lower dc bias values we observe a second disconnected branch (the initial value of dc bias was 0.27) which displays a back bending and switching to the HWS in the I-V curve as the dc bias was increased. With decreasing current this breather state transits to the superconducting state. All these features occurring in the current region 0.28 ≤ γ ≤ 0.33, can be explained as a parametric resonant interaction between the breather state and excited EOs. First, we found a good agreement between the breather frequency Ω = γ α(1+η) = 1.6 and the frequency of EOs 2ω − = 1.52 ( for γ = 0.32 ). To study this resonant interaction more precisely we carried out the Floquet analysis of the breather state. Fig. 6 illustrates the evolution of the modulus and the argument of the eigenvalues of the Floquet matrix M with γ. The arguments of the eigenvalues can be restricted as 0 ≤ Arg(λ) ≤ π. As the dc bias decreases we observed that for the particular dc bias γ = 0.34 an eigenvalue of the Floquet matrix crosses the unit circle at −1 (Fig. 6). Thus the original breather state becomes unstable (although it continues to exist as a solution of the equations of motion) and this instability is driven by the parametric resonance (Fig. 6b) with the EO of the frequency ω − . Note, that as the current decreases further this eigenvalue returns into the unit circle, the breather state is again stable and thus an instability window of the breather is found (the region A 1 in Fig.3(a,b)). The position of this instability window in the parameter space γ-η can be predicted using Eq. (22) or in the form (thin solid lines in Fig. 3): The appearance of parametric resonance is often accompanied by a period doubling of the state mediating the resonance, and thus we checked the time dependence of the Josephson phase φ v 1 (t). The results for two different dc bias values are presented in Fig. 7. One can see that for the dc bias γ = 0.2965 the period of the Josephson phase φ v 1 (t) is two times less than the period of φ v 1 (t) for the slightly larger value of γ = 0.2970. Thus, the "back bending" (and in same manner the"cusp") of the I-V curve (Fig. 5) is related to this period doubling process. As the anisotropy η increases and approaches the tip of the domain of breather existence (Fig. 3) the I-V characteristics start to display different features. First, for such a large value of anisotropy (η ≥ 2) we observed a large deviation from the dc analysis ( Fig. 8-10). It indicates a strong increase of the amplitudes of EOs. Presumably this is also the reason why the switching to the HWS occurs for considerably smaller γ values as compared with the dc analysis (see also Ref. [ 23]). This circumstance leads to a shrinking of the domain of breather existence at the tip.
We also found that the switching outcome in this region of anisotropy upon lowering the dc bias γ crucially depends on η. Indeed, for η = 2.1 we observed the switching to the superconducting state (Fig. 8) but the breather state switches to the HWS (Fig. 9) as the parameter η becomes slightly above the value 2.11. Moreover, for η = 2.3 (Fig. 10) we observed that the switching to the HWS occurs via a peculiar intermediate state. This state is not periodic anymore (either quasiperiodic or chaotic in time), yet it shows up with typical breather features when averaged over time. It was stable over the whole simulation time of 10 6 . We carried out the Floquet analysis of the breather state and found the instability region (region A 2 in Fig  3) that is also due to the parametric resonance between the breather state and EOs. At variance with the considered above instability window A 1 , this instability region adjoins the lower boundary of the breather existence domain.

Finite self-inductance case (βL = 0)
The breather properties depend crucially on the selfinductance parameter β L . First, we found that the domain of the breather existence shrinks for small and moderate values of β L (see Fig. 11). Moreover, the upper boundary of the breather existence domain displays a "saw tooth" form for moderate values of the anisotropy (η > 0.5, see Fig. 11). The reason for such a peculiar shape of the breather existence domain is the resonant interaction between the breather state and EOs (thin solid lines in Fig. 11) leading to large amplitudes of EOs and corresponding ac induced instabilities of the superconducting state of the left vertical junction. As the parameter β L increases to large values (see Fig. 12) the resonant interaction between the breather state and EOs becomes weak again. In this case the domain of breather existence practically coincides with the analytical predictions in Section II (see Fig. 2). Also, similar to the case of β L = 0 we found various instabilities of the breather state driven by the resonant interaction with EOs. These instability windows for a particular case of β L = 1 are shown in Fig. 13. For large values of the dc bias (γ ≃ 0.6) a parametric resonance of the breather state with high frequency EOs Ω = 2ω + leads to a resonant step in the I-V curve. Moreover, the breather state biased to this resonance displays a period doubling bifurcation which manifests through a sharp corner in the I-V curve at the bottom of the resonant step (A+ region). As the dc bias decreases to the value γ ≃ 0.4 we observe the combination resonance Ω = ω + + ω − . This resonance is characterized by a large back bending of the I − V curve and the breather state dynamics displays a non-periodic behaviour (see inset in Fig. 14 and C region).
As the dc bias decreases further to the value of γ ≃ 0.37, a large resonant step in the I-V curve appears. It corresponds to the primary resonance with Ω = ω + . This resonant breather state switches to the superconducting state at the value of γ = 0.27 as the result of an instability driven by the parametric resonance with low frequency EOs Ω = 2ω − . These features are clearly seen in the dependence of Floquet multipliers on the dc bias (see Fig. 15). The type of resonance can be easily detected from the behaviour of the arguments of the Floquet multipliers. While unstable solutions (with absolute values of multipliers being larger than one) which corresponding to parametric resonances are characterized by 'frozen' arguments at value π, combination resonances are characterized by a merging of arguments over a particular interval of dc bias.
For even smaller values of the dc bias, namely 0.21 < γ < 0.26) the breather state becomes again stable. For γ ≃ 0.26 this breather state shows up with a non-periodic behaviour (see inset in Fig. 14 (b)). Most important we observe three new resonances: a combination resonance Ω = ω + − ω − , a parametric resonance Ω = 2 3 ω + and yet another combination resonance 2Ω = ω + + ω − . The breather finally switches to the superconducting state at γ = 0.21.

IV. CONCLUSION
We have presented an analytical and numerical study of inhomogeneous dynamic states (breathers) in a single plaquette consisting of three small underdamped Josephson junctions. We found that as the damping α is small the breather state as a solution of dynamic equations exists in a large region of parameters: dc bias γ and anisotropy η (Fig. 2, 3, 12 and 13). For small and moderate values of anisotropy η (η ≤ 1) the dc analysis that completely neglects the breather induced EOs (Eqs. (6)-(13)), results in a good estimation of the boundary of breather existence. However, we found that the domain of the breather existence shrinks as the anisotropy η increases, and moreover, for moderate values of the self inductance β L the upper part of the boundary displays a "saw tooth" feature. These effects are due to the excitation of EOs with large amplitudes.
The second major result of the presence of EOs is that the domain of the breather existence contains a large number of instability regions. The physical origin of these instabilities is the resonant interaction between the breather state and the EOs. It is important that due to the presence of two different characteristic frequencies of EOs (for the dependence of ω ± on the self inductance and anisotropy, see Eqs. (18)), a rich set of resonances (primary, parametric or combination ones) occurs. Such resonant interactions can completely destroy the breather state and correspondingly lead to significant shifts of the switchings of the breather state to the HWS or the superconducting state. These switchings appear through voltage jumps in the I-V curves. Note here, that in most cases we observed the switching to the superconducting state as a result of breather instability (parametric or combination) but not due to the standard retrapping mechanism (Eq. (12)). Similar findings have been also reported for large Josephson junction ladders. 22 The most interesting observation is that the breather state can survive the resonant instabilities by increasing its dynamic complexity. This dynamic complexity manifests itself by resonant steps of various types in the I-V curves. Indeed, we observed a resonant type breather as a result of a primary resonance with EOs (B + region in Fig. 14a), a period doubling process resulting from the parametric resonant interaction (Fig. 5, regions A + , A − in Fig. 14) and even a non-periodic (quasi periodic or chaotic) behaviour in parameter regions where combination resonances appear (C region and the region of small current in Fig. 14).
Finally, all observed effects (switching outcomes, instabilities) are very sensitive to the damping parameter α, and different combinations of the resonances can appear as the damping parameter α decreases.
What can we learn about breather states in extended ladder geometries? Our analysis shows that stable breathers can be induced for large self inductance. If on the other hand one needs switchings between different breather states and resonances, intermediate β L values should be chosen. In addition variations of the damping (e.g. via temperature control) may increase the complexity of expected switching scenarios. At the same time it should be clear that if already such a small system (with three degrees of freedom) shows up with a high complexity of the phase space structure, the situation for large numbers of degrees of freedom as for Josephson ladders will be even more complex.