text:
# Review of the PQCD approach to exclusive 𝐵 decays

## 1 Introduction

Perturbative QCD (PQCD) seems to be a successful approach to exclusive $`B`$ meson decays. In this talk I will briefly review important aspects of this approach, concentrating on factorization of infrared divergences, gauge invariance of meson light-cone distribution amplitudes, the smearing of end-point singularities, Sudakov suppression from $`k_T`$ and threshold resummations, power counting of various topologies of diagrams for hadronic $`B`$ meson decays, dynamical enhancement of penguin contributions, large strong phases from annihilation diagrams, and CP asymmetries in $`BK\pi `$ and $`\pi \pi `$ decays.

## 2 Factorization theorem

We start with the lowest-order diagram for the $`B\pi `$ form factor in the kinematic region with a fast-recoil pion, which contains a hard gluon exchanged between the $`b`$ quark and the spectator quark. The spectator quark in the $`B`$ meson, forming a soft cloud around the heavy $`b`$ quark, carries momentum of order $`\overline{\mathrm{\Lambda }}=M_Bm_b`$, $`M_B`$ ($`m_b`$) being the $`B`$ meson ($`b`$ quark) mass. The spectator quark on the pion side carries momentum of $`O(M_B)`$ in order to form the fast-moving pion with the $`u`$ quark produced in the $`b`$ quark decay. These dramatic different orders of magnitude in momenta explain why the hard gluon is necessary. Based on the above argument, the hard gluon is off-shell by order of $`\overline{\mathrm{\Lambda }}M_B`$. As explored later, this scale, characterizing heavy-to-light decays, is important for developing the PQCD formalism of exclusive $`B`$ meson decays.

At higher orders, infinitely many gluon exchanges appear. If all these gluons are hard, higher-order contributions are calculable in perturbation theory. Unfortunately, these diagrams generate infrared divergences. There are two types of infrared divergences, soft and collinear. Soft divergences come from the region of a loop momentum $`l`$, where all its components vanish:

$`l^\mu =(l^+,l^{},l_T)(\overline{\mathrm{\Lambda }},\overline{\mathrm{\Lambda }},\overline{\mathrm{\Lambda }}).`$ (1)
Collinear divergences are associated with a massless quark of momentum $`P(M_B,0,0_T)`$. In the collienar region with $`l`$ parallel to $`P`$, the components of $`l`$ behave like

$`l^\mu (M_B,\overline{\mathrm{\Lambda }}^2/M_B,\mathrm{\Lambda }).`$ (2)
In both regions the invariant mass of the radiated gluon diminishes as $`\overline{\mathrm{\Lambda }}^2`$, and the corresponding loop integrand may diverge like $`1/\overline{\mathrm{\Lambda }}^4`$. As the phase space for loop integration vanishes like $`d^4l\overline{\mathrm{\Lambda }}^4`$, logarithmic divergences are generated.

It has been proved $`^\mathrm{?}`$ that the soft divergences in the $`B\pi `$ form factor can be factorized into a light-cone $`B`$ meson distribution amplitude, and that the collinear divergences can be factorized into a pion distribution amplitude. The remaining finite piece is assigned into a hard part. Factorization of infrared divergences is performed in momentum, spin, and color spaces. Factorization in momentum space means that a hard part does not depend on the loop momentum of a soft or collinear gluon, which has been absorbed into a meson distribution amplitude. Factorization in spin and color spaces means that there are separate fermion and color flows between a hard part and a distribution amplitude, respectively. To achieve these, we rely on the eikonal approximation for loop integrals in leading infrared regions, the insertion of the Fierz identity to separate fermion flows, and the Ward identity to sum up diagrams with different color structures. Under the eikonal approximation, a soft or collinear gluon is detached from the lines in a hard part and in other distribution amplitudes. The Fierz identity decomposes a full amplitude into contributions characterized by different twists. The Ward identity is essential for deriving factorization theorem in a nonabelian gauge theory.

I emphasize that the hard scale $`\overline{\mathrm{\Lambda }}M_B`$ is essential for constructing a gauge invariant $`B`$ meson distribution amplitude,

$`\varphi _B(x)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2N_c}M_B}}{\displaystyle \frac{dy^{}}{2\pi }e^{ixP_B^+y^{}}0|\overline{q}(y^{})\frac{\gamma _5\mathit{}}{2}}`$ (3)
$`\times P\mathrm{exp}\left[ig{\displaystyle _0^y^{}}𝑑znA(zn)\right]`$
$`\times b_v(0)|B(P_B),`$
with the $`B`$ meson momentum $`P_B`$, the velocity $`v=P_B/M_B`$, the light-like vector $`n=(0,1,0_T)`$, the number of colors $`N_c`$, the momentum fraction $`x`$ associated with the spectator quark, and the rescaled $`b`$ quark field

$`b_v(w)=\mathrm{exp}(iM_Bvw){\displaystyle \frac{\mathit{}+I}{2}}b(w).`$ (4)
The above nonlocal matrix element is gauge invariant because of the presence of the path-ordered Wilson line integral. A careful investigation shows that the $`O(\alpha _s^2)`$ diagram with the second gluon attaching the hard gluon contributes to this line integral. That is, this diagram contains the soft divergence, which is factorized into $`\varphi _B`$. This is possible, only when the hard gluon is off-shell by the intermediate scale $`\overline{\mathrm{\Lambda }}M_B`$, rather than by $`\overline{\mathrm{\Lambda }}^2`$ or $`M_B^2`$.

A hard part is calculable in perturbation theory. A meson distribution amplitude, though not calculable, is universal, since it absorbs long-distance dynamics, which is insensitive to a specific decay of the $`b`$ quark into light quarks with large energy release. The universality of nonperturbative distribution amplitudes is the fundamental concept of PQCD factorization theorem. Because of this universality, one can extract distribution amplitudes from experimental data, and employ them to make model-independent predictions for other processes.

## 3 End-point singularities

After developing factorization theorem, we calculate the $`B\pi `$ form factor $`F^{B\pi }(q^2)`$ at large recoil, where $`q`$ denotes the lepton pair momentum. A difficulty immediately occurs. The lowest-order diagram for the hard part, which was mentioned above, is proportional to $`1/(x_1x_2^2)`$, $`x_1`$ ($`x_2`$) being the momentum fraction associated with the spectator quark on the $`B`$ meson (pion) side. If the pion distribution amplitude vanishes like $`x_2`$ as $`x_20`$ (in the leading-twist, i.e., twist-2 case), $`F^{B\pi }`$ is logarithmically divergent. If the pion distribution amplitude is a constant as $`x_20`$ (in the next-to-leading-twist, i.e., twist-3 case), $`F^{B\pi }`$ even becomes linearly divergent. These end-point singularities have caused critiques on the perturbative evaluation of the $`B\pi `$ form factor.

Several methods have been proposed to regulate the above end-point singularities. An on-shell $`b`$ quark propagator has been subtracted from the hard part as $`x_20`$ in $`^\mathrm{?}`$. However, this subtraction renders the lepton energy spectrum of the semileptonic decay $`B\pi l\overline{\nu }`$ vanishes as the lepton energy is equal to half of its maximal value. Obviously, this vanishing is unphysical. The subtraction also leads to a value of $`F^{B\pi }(0)`$, which is much smaller than the expected one 0.3. A lower bound of $`x_2`$ of $`O(\overline{\mathrm{\Lambda }}/M_B)`$ has been introduced in $`^\mathrm{?}`$ to make the convolution integral finite. However, the outcomes depend on this cutoff sensitively, and PQCD loses its predictive power.

A QCD-consistent prescription has been proposed in $`^\mathrm{?}`$, where parton transverse moemta $`k_T`$ are retained in internal particle propagators involved in the hard part. The inclusion of $`k_T`$ certainly brings in large double logarithms $`\alpha _s\mathrm{ln}^2(k_T/M_B)`$ through radiative corrections. These large logarithms should be resummed in order to improve the perturbative calculation. The $`k_T`$ resummation $`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$ then sets a distribution of $`k_T`$, such that the average of $`k_T^2`$ is numerically around $`k_T^2\overline{\mathrm{\Lambda }}M_B`$ for $`M_B5`$ GeV. The off-shellness of internal particles then remain of $`O(\overline{\mathrm{\Lambda }}M_B)`$ even in the end-point region, and the end-point singularities are smeared out. This is so-called Sudakov suppression.

Recently, another type of resummation has been observed. It is not difficult to verify that the loop correction to the weak decay vertex produces the double logarithms $`\alpha _s\mathrm{ln}^2x_2`$ $`^\mathrm{?}`$. These double logarithms can be factored out of the hard part systematically, and grouped into an exclusive quark jet function. If the end-point region is important, these large logarithms need to be resummed $`^\mathrm{?}`$ in order to improve the perturbative calculation. The threshold resummation $`^{\mathrm{?},\mathrm{?}}`$ for the jet function results in Sudakov suppression, which decreases faster than any power of $`x_2`$ as $`x_20`$, and removes the end-point singularities $`^\mathrm{?}`$. In conclusion, if the PQCD analysis of the $`B\pi `$ form factor is performed self-consistently, there exist no end-point singularities.

The mechanism of Sudakov suppression can be easily understood by regarding a meson as a color dipole. In the region with vanishing $`k_T`$ and $`x`$, the meson possesses a huge extent in the transverse and longitudinal directions, respectively. That is, the meson carries a large color dipole. At fast recoil, this large color dipole, strongly scattered, tends to emit real gluons. However, these emissions are forbidden in an exclusive process with final-state particles specified. As a consequence, contributions to the $`B\pi l\overline{\nu }`$ decay from the region with vanishing $`k_T`$ and $`x`$ must be highly suppressed.

## 4 Power counting

We then discuss two-body hadronic $`B`$ meson decays, such as $`BPP`$. These modes involve three scales $`^{\mathrm{?},\mathrm{?},\mathrm{?}}`$: the $`W`$ boson mass $`M_W`$, at which the matching conditions of the effective weak Hamiltonian to the full Hamiltonian are defined, the typical scale $`t\sqrt{\overline{\mathrm{\Lambda }}M_B}`$ of a hard part, which reflects the specific dynamics of a decay mode, and the factorization scale $`\overline{\mathrm{\Lambda }}`$, which characterizes infrared divergences. Above the factorization scale, perturbation theory is reliable, and radiative corrections produce two types of large logarithms: $`\mathrm{ln}(M_W/t)`$ and $`\mathrm{ln}(t/\overline{\mathrm{\Lambda }})`$. The former are summed by renormalization-group (RG) equations to give the evolution from $`M_W`$ down to $`t`$ described by the Wilson coefficients $`c(t)`$, while the latter are summed to give the evolution from $`t`$ to $`\overline{\mathrm{\Lambda }}`$. When the argument of a Wilson coefficient $`c(\mu )`$ is set to $`\mu =t`$, we have included the logarithmic piece of the vertex corrections to the four fermion operators to all orders. These vertex corrections have been considered in a modified factorization approach $`^{\mathrm{?},\mathrm{?}}`$, which retains the scale independence of predictions for hadronic $`B`$ meson decays. The finite piece, being of higher orders in the PQCD approach, is dropped.

The hard part for hadronic $`B`$ meson decays contains all possible Feynman diagrams, such as factorizable diagrams, where hard gluons attach the valence quarks in the same meson, and nonfactorizable diagrams, where hard gluons attach the valence quarks in different mesons. The annihilation topology is also included, and classified into factorizable or nonfactorizable one according to the above definitions. Below I will discuss the power counting of these different topologies of diagrams and argue that all of them should be included in the leading-power PQCD analysis.

As explained in Sec. 2, factorizable amplitudes scale like $`1/\overline{\mathrm{\Lambda }}M_B`$, since the end-point singularities do not exist. For a similar reason, each nonfactorizable diagram also scales like $`1/\overline{\mathrm{\Lambda }}M_B`$. However, because of the soft cancellation between a pair of nonfactorizable diagrams, their sum turns out to scale like $`1/M_B^2`$. I emphasize that it is more appropriate to count the power of each individual diagram, instead of the power of sum of diagrams. In some case factorizable contributions are suppressed by a vanishing Wilson coefficient, such that nonfactorizable contributions become dominant. For example, factorizable internal-$`W`$ emisson contributions are strongly suppressed by the Wilson coefficient $`a_2`$ in the $`BJ/\psi K^{()}`$ decays $`^\mathrm{?}`$. In some case, such as the $`BD\pi `$ decays, there is no soft cancellation between a pair of nonfactorizable diagrams, and nonfactorizable contributions also become important $`^\mathrm{?}`$.

A folklore for annihilation contributions is that they are negligible compared to emission contributions. The annihilation conributions from the operators $`O_{1,2,3,4}`$ with the structure $`(VA)(VA)`$ are small because of helicity suppression. Those from the operators $`O_{5,6}`$ with the structure $`(SP)(S+P)`$, though surviving under helicity suppression, are of $`O(1/M_B^2)`$. This assumption is reasonable, because the hard gluon in an annihilaton diagram is off-shell by $`x_2x_3M_B^2`$, where $`x_2O(1)`$ and $`x_3O(1)`$ are the momentum fractions associated with the two fast outgoing light mesons. However, this argument applies only to the real part of annihilation contributions, but not to the imaginary part.

To obtain the imaginary part, the internal quark is required to be on mass shell. Its propagator, proportional to

$`{\displaystyle \frac{1}{x_2M_B^2k_T^2}},`$ (5)
then gives

$`x_2M_B^2=k_T^2\overline{\mathrm{\Lambda }}M_B,x_2{\displaystyle \frac{\overline{\mathrm{\Lambda }}}{M_B}},`$ (6)
after considering Sudakov suppression from $`k_T`$ resummation. This smaller $`x_2`$ implies that the hard gluon off-shell only by $`\overline{\mathrm{\Lambda }}M_B`$ contributes to the imaginary part, and that the imaginary annihilation amplitudes possess a power behavior the same as factorizable emission ones. Note that the suppression from threshold resummation decreases the above estimation by a factor 3 or 4. This is the reason that an annihilation amplitude is almost purely imaginary, and its magnitude is usually few times smaller than the factorizable emission ones.

At last, I argue that two-parton twist-3 distribution amplitudes, though proportional to the ratio $`m_0/M_B`$ with the mass $`m_0`$ related to chiral condensate, should be included. As stated before, the corresponding convolution integral for the $`B\pi `$ form factor is linearly divergent. This integral, regulated in some way with an effective cutoff $`x_c\overline{\mathrm{\Lambda }}/M_B`$, is proportional to the ratio $`M_B/\overline{\mathrm{\Lambda }}`$. Combining the two ratios $`m_0/M_B`$ and $`M_B/\overline{\mathrm{\Lambda }}`$, contributions from two-parton twist-3 distribution amplitudes are in fact not suppressed by a power of $`1/M_B`$:

$`{\displaystyle \frac{m_0}{M_B}}{\displaystyle _{\overline{\mathrm{\Lambda }}/M_B}^1}{\displaystyle \frac{dx_2}{x_2^2}}{\displaystyle \frac{m_0}{\overline{\mathrm{\Lambda }}}},`$ (7)
and should be taken into account. Hence, it is appropriate to claim that the PQCD formalism is complete at leading power.

## 5 Phenomenology

Another important phenomenological consequence related to the special scale $`t\sqrt{\overline{\mathrm{\Lambda }}M_B}`$ is the dynamical enhancement of penguin contributions $`^\mathrm{?}`$. The RG evolution of the Wilson coefficients $`C_{4,6}(t)`$ dramatically increase as $`t<M_B/2`$, while that of $`C_{1,2}(t)`$ almost remain constant $`^\mathrm{?}`$. With this penguin enhancement, the observed branching ratios of the $`BK\pi `$ decays, about three times larger than those of the $`B\pi \pi `$ decays, can be explained for a smaller unitarity angle $`\varphi _3<90^o`$. Note that the former are dominated by penguin contributions and the latter are dominated by tree contributions. In the factorization approach $`^\mathrm{?}`$ and in the QCD factorization approach $`^\mathrm{?}`$, it is assumed that factorizable contributions are not calculable. The leading contribution to a hadronic decay amplitude is then expressed as a convolution of a hard part with a form factor and a meson distribution amplitude. In both approaches the only hard scale is $`M_B`$ and the intermediate scale $`\overline{\mathrm{\Lambda }}M_B`$ can not appear. Therefore, the dynamical enhancement of penguin contributions does not exist. It is then difficult to account for the $`BK\pi `$ data.

To accommodate the $`BK\pi `$ data in the factorization and QCD factorization approaches, one relies on the chiral enhancement by increasng the mass $`m_0`$ to an unreasonablly large value $`m_04`$ GeV $`^\mathrm{?}`$. Whether dynamical enhancement or chiral enhancement is responsible for the large $`BK\pi `$ branching ratios can be tested by measuring the $`B\varphi K`$ modes. In these modes penguin contributions dominate and the mass $`m_0`$ is replaced by the $`\varphi `$ meson mass $`M_\varphi 1`$ GeV. If the branching ratios of the $`B\varphi K`$ decays are around $`4\times 10^6`$ $`^\mathrm{?}`$, the chiral enhancement may be essential for the $`BK\pi `$ decays. Without including annihilation contributions which depend on an end-point cutoff sensitively, similar values of the $`B\varphi K`$ branching ratios have been derived in $`^\mathrm{?}`$ using the QCD factorization approach. If the branching ratios are around $`10\times 10^6`$ as predicted in the PQCD approach $`^\mathrm{?}`$, the dynamical enhancement may be essential.

Because of the large imaginary annihilation contributions in the PQCD formalism, the predicted CP asymmetry ($`30\%`$) in the $`B^0\pi ^\pm \pi ^{}`$ decays $`^\mathrm{?}`$ dominates over that ($`5\%`$) in the QCD factorization approach and over that ($`10\%`$) in the factorization approach $`^\mathrm{?}`$. The CP asymmetries in the $`BK\pi `$ decays can also reach 15% $`^\mathrm{?}`$. Future measurements of CP asymmetries can distinguish these different approaches to hadronic $`B`$ meson decays. For numerical results of all the hadronic modes that have been studied in PQCD, refer to Dr. Y.Y. Keum’s talk in this workshop $`^\mathrm{?}`$.

## 6 Conclusion

In this talk I have briefly reviewed important aspects and phenomenological consequences of the PQCD approach to exclusive $`B`$ meson decays. I have explained that the end-point singularities do not exist in a self-consistent PQCD formalism because of Sudakov suppression from $`k_T`$ and threshold resummations. I have emphasized the special characteristic scale $`\overline{\mathrm{\Lambda }}M_B`$, under which a gauge-invaraint $`B`$ meson light-cone distribution amplitude can be constructed, annihilation diagrams should be included in a leading-power analysis, penguin contributions are dynamically enhanced, and large CP asymmetries in the $`BK\pi `$ and $`\pi \pi `$ decays have been predicted.

More applications of the PQCD approach to other two-body $`B`$ meson decays have been discussed in other talks and posters of this workshop.

## Acknowledgments

I thank I. Bigi, S. Brodsky, H.Y. Cheng, A. Falk, Y.Y. Keum, T. Kuromoto, D. London, B. Melic, T. Mannel, T. Morozumi, A.I. Sanda. and L. Silvestrini for useful discussions and thank Theory Group of SLAC for hospitality during my visit. The work was supported in part by the National Science Council of R.O.C. under the Grant No. NSC-89-2112-M-006-033, and in part by Grant-in Aid for Special Project Research (Physics of CP Violation) and by Grant-in Aid for Scientific Exchange from Ministry of Education, Science and Culture of Japan.

## References


format: text

text:
# Excited electron contribution to the e⁺e⁻→𝛾𝛾 cross-section

## Abstract

The differential cross-section for the process $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ is calculated assuming the presence of excited electrons with a chiral magnetic coupling. This calculation permits constraining the excited electron coupling using the same theoretical framework as the one generally used for direct production searches.

Excited electrons could be produced directly in $`\mathrm{e}^+\mathrm{e}^{}`$ collisions. They are also expected to alter the $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ differential cross-section via the diagrams shown in Fig. 1. Currently, the limits on the size of the $`\mathrm{e}^{}\mathrm{e}\gamma `$ coupling obtained from direct production LEP\_direct\_search  searches are generally calculated using one form of the $`\mathrm{e}^{}\mathrm{e}\gamma `$ interaction whereas limits extracted from the $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ differential cross-section LEP\_indirect\_search  use a different form. By using a common theoretical framework, the $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ process can be used to extend the limits on the $`\mathrm{e}^{}\mathrm{e}\gamma `$ coupling strength for excited electron masses beyond the kinematically allowed region of direct production. This paper presents a brief overview of both forms of the $`\mathrm{e}^{}\mathrm{e}\gamma `$ interaction and the analytic expression for the $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ differential cross-section derived in the framework used by direct production searches. This differential cross-section is compared to the form currently in use and the impact on the limits of the $`\mathrm{e}^{}\mathrm{e}\gamma `$ coupling strength is discussed.

Indirect searches, where the existence of $`\mathrm{e}^{}`$ is inferred from deviations in the $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ differential cross-section, usually express limits on the $`\mathrm{e}^{}\mathrm{e}\gamma `$ coupling strength assuming a general extension of the Standard Model. The interaction between an excited lepton, a lepton and a gauge boson ($`\mathrm{}^{}\mathrm{}V`$) is described by the simplest gauge invariant form of the interaction Lagrangian Litke

$$_\mathrm{}^{}\mathrm{}V=\frac{e}{2}\frac{\kappa }{M_{\mathrm{}^{}}}\overline{\mathrm{}}^{}\sigma _{\mu \nu }\mathrm{}\mathrm{F}^{\mu \nu }+\mathrm{h}.\mathrm{c}.$$
(1)
where $`\mathrm{F}^{\mu \nu }`$ denotes the electromagnetic field tensor, $`\sigma _{\mu \nu }`$ is the covariant bilinear tensor and $`M_{\mathrm{}^{}}`$ is the mass of the excited lepton. The parameter $`\kappa `$ is a measure of the coupling strength. The $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ differential cross-section using this purely magnetic coupling is explicitly calculated in Litke  and given by

$`\left({\displaystyle \frac{d\sigma }{d\mathrm{\Omega }}}\right)_\mathrm{e}^{}=`$ $`\left({\displaystyle \frac{d\sigma }{d\mathrm{\Omega }}}\right)_{\mathrm{QED}}`$ (2)
$`+\alpha ^2\{{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\kappa }{M_\mathrm{e}^{}}}\right)^4(E^2\mathrm{sin}^2\theta +M_\mathrm{e}^{}^2)({\displaystyle \frac{q^4}{\left(q^2M_\mathrm{e}^{}^2\right)^2}}+{\displaystyle \frac{q^4}{\left(q^2M_\mathrm{e}^{}^2\right)^2}})`$
$`+\mathrm{\hspace{0.33em}4}\left({\displaystyle \frac{\kappa }{M_\mathrm{e}^{}}}\right)^4{\displaystyle \frac{M_\mathrm{e}^{}^2E^4\mathrm{sin}^2\theta }{(q^2M_\mathrm{e}^{}^2)(q^2M_\mathrm{e}^{}^2)}}`$
$`+\left({\displaystyle \frac{\kappa }{M_\mathrm{e}^{}}}\right)^2[{\displaystyle \frac{q^2}{q^2M_\mathrm{e}^{}^2}}+{\displaystyle \frac{q^2}{q^2M_\mathrm{e}^{}^2}}+E^2\mathrm{sin}^2\theta ({\displaystyle \frac{1}{q^2M_\mathrm{e}^{}^2}}+{\displaystyle \frac{1}{q^2M_\mathrm{e}^{}^2}})]\}`$
where $`\left(\frac{d\sigma }{d\mathrm{\Omega }}\right)_{\mathrm{QED}}`$ is the Born level Standard Model differential cross-section, $`\theta `$ is the polar angle of one of the photons with respect to the incoming electron, $`E`$ is the beam energy ($`E=\sqrt{s}/2`$), $`q^2=2E^2(1\mathrm{cos}\theta )`$ and $`q^2=2E^2(1+\mathrm{cos}\theta )`$. Since the two outgoing photons are indistinguishable, $`\mathrm{cos}\theta `$ is defined to be positive. Limits on the strength of the $`\mathrm{e}^{}\mathrm{e}\gamma `$ coupling, $`\kappa `$, are expressed as a function of $`M_\mathrm{e}^{}`$ LEP\_indirect\_search .

The interaction Lagrangian of Equation 1 leads to large contributions to the anomalous magnetic moment of electrons and muons and the coupling is therefore severely constrained by existing g-2 precision measurements Renard ; Brodsky . In fact, limits from g-2 measurements are approximately an order of magnitude better than limits from $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ calculated using Equation 2.

Limits on the $`\mathrm{e}^{}\mathrm{e}\gamma `$ coupling strength from the search for directly produced excited leptons are calculated using a different theoretical framework. The effective Lagrangian describing the $`\mathrm{}^{}\mathrm{}V`$ interaction is chosen to have a chiral symmetry which protects Standard Model leptons from acquiring large anomalous magnetic moments. This Lagrangian is generally written as BDK ; Cabibbo ; Kuhn ; Hagiwara ; Boudjema

$$_\mathrm{}^{}\mathrm{}V=\frac{1}{2\mathrm{\Lambda }}\overline{\mathrm{}}^{}\sigma ^{\mu \nu }\left[gf\frac{𝝉}{2}𝑾_{\mu \nu }+g^{}f^{}\frac{Y}{2}B_{\mu \nu }\right]\mathrm{}_L+\mathrm{h}.\mathrm{c}.$$
(3)
where $`g`$ and $`g^{}`$ are the usual Standard Model couplings, the tensors $`𝑾_{\mu \nu }`$ and $`B_{\mu \nu }`$ represent the Standard Model gauge-invariant field tensors, $`𝝉`$ denotes the Pauli matrices and $`Y`$ is the weak hypercharge. The compositeness scale is set by the parameter $`\mathrm{\Lambda }`$ which has units of energy. The size of the $`\mathrm{}^{}\mathrm{}V`$ coupling is governed by the constants $`f`$ and $`f^{}`$ which can be interpreted as weight factors associated to the different gauge groups. In the physical basis, the Lagrangian of Equation 3 leads to the following chiral magnetic vertex BDK

$$\mathrm{\Gamma }_\mu ^\mathrm{}^{}\mathrm{}V=\frac{e}{2\mathrm{\Lambda }}f_Vq^\nu \sigma _{\mu \nu }(1\gamma _5)$$
(4)
where $`q^\nu `$ is the momentum of the gauge boson and $`f_V`$ represents a particular combination of the constants $`f`$ and $`f^{}`$ for a given boson interaction. For an excited electron coupling to photon ($`V=\gamma `$BDK , $`f_\gamma =\frac{1}{2}(f+f^{})`$. In this model, experimental limits on the $`\mathrm{}^{}\mathrm{}V`$ coupling are usually expressed in terms of limits on $`|f|/\mathrm{\Lambda }`$ as a function of the excited lepton mass for the hypothesis 111For the specific $`\mathrm{e}^{}\mathrm{e}\gamma `$ interaction the photon coupling $`f_\gamma `$ is zero for the hypothesis $`f=f^{}`$ but the coupling to $`Z^0`$ and $`W^\pm `$ does not vanish BDK . $`f=f^{}`$ or $`f=f^{}`$LEP\_direct\_search .

Using this chiral magnetic interaction, the coupling is less severely constrained since contributions to the electron and muon anomalous magnetic moment are suppressed. It still however permits observable deviations in the process $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ which are not excluded by g-2 measurements. In addition, limits from indirect searches expressed in this framework could be easily compared and combined with limits coming from direct searches. To achieve this, deviations from the Standard Model $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ differential cross-section must be calculated assuming a chiral magnetic $`\mathrm{e}^{}\mathrm{e}\gamma `$ coupling.

With the existence of an excited electron, the four diagrams shown in Fig. 1, contribute to the Born level production of $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ events. The differential cross-section is calculated using the $`\mathrm{}^{}\mathrm{}V`$ vertex given in Equation 4 and combined with the standard QED interaction $`e\overline{\mathrm{}}\gamma _\mu \mathrm{}\mathrm{A}_\mu `$. The excited electron propagator is taken to be the usual fermion expression with a mass $`M_\mathrm{e}^{}`$. Summing over the outgoing photon polarizations and neglecting the mass of the electron, the resulting differential cross-section is

$`\left({\displaystyle \frac{d\sigma }{d\mathrm{\Omega }}}\right)_\mathrm{e}^{}`$ $`=`$ $`\left({\displaystyle \frac{d\sigma }{d\mathrm{\Omega }}}\right)_{\mathrm{QED}}`$ (5)
$`+{\displaystyle \frac{\alpha ^2}{4}}{\displaystyle \frac{f_\gamma ^4}{\mathrm{\Lambda }^4}}E^2\mathrm{sin}^2\theta \left[{\displaystyle \frac{q^4}{\left(q^2M_\mathrm{e}^{}^2\right)^2}}+{\displaystyle \frac{q^4}{\left(q^2M_\mathrm{e}^{}^2\right)^2}}\right]`$
$`{\displaystyle \frac{\alpha ^2}{8E^2}}{\displaystyle \frac{f_\gamma ^2}{\mathrm{\Lambda }^2}}\left[{\displaystyle \frac{q^4}{\left(q^2M_\mathrm{e}^{}^2\right)}}+{\displaystyle \frac{q^4}{\left(q^2M_\mathrm{e}^{}^2\right)}}\right]`$
where the same notation as for Equation 2 is used. Terms of order $`\left(f_\gamma /\mathrm{\Lambda }\right)^2`$ in Equation 5 describe the interference between the excited electron and Standard Model diagrams.

Figure 2 shows the Standard Model Born level differential cross-section for the process $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ compared with the prediction obtained using a purely magnetic (Equation 2) and chiral magnetic (Equation 5) $`\mathrm{e}^{}\mathrm{e}\gamma `$ coupling. For comparison purposes, the coupling parameters are set to $`\kappa =1`$ in Equation 2, and $`f_\gamma =1`$ and $`\mathrm{\Lambda }=M_\mathrm{e}^{}`$ in Equation 5. The differential cross-section assuming a purely magnetic coupling is larger than the one assuming a chiral magnetic coupling.

In summary, the differential cross-section for the process $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ assuming the existence of an excited electron with chiral magnetic coupling has been calculated. The differential cross-section predictions from chiral magnetic and purely magnetic couplings are compared. Limits on the size of the $`\mathrm{e}^{}\mathrm{e}\gamma `$ coupling determined from the $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ differential cross-section can be combined with limits obtained from direct production searches when the chiral magnetic form is used.

I would like to thank Robert McPherson, Randall Sobie and Matt Dobbs for their useful comments and careful reading of the manuscript. I am also grateful to Kirsten Sachs for stimulating discussions. This work was supported by the Natural Sciences and Engineering Research Council of Canada.


text:
# Breathers in Josephson junction ladders: resonances and electromagnetic waves spectroscopy

## I Introduction

Various nonlinear and discrete systems (nonlinear lattices) have attracted a lot of interest as they display diverse fascinating phenomena.  Well known examples of such phenomena are solitary excitations, propagation of (non)linear waves and appearance of various inhomogeneous structures.

Moreover, the interest in this area was boosted by the prediction, theoretical analysis  and the subsequent observation  of intrinsic dynamic localized excitations (discrete breathers) that are periodic in time and localized in space. Note here, that the origin of such dynamical localization is not the presence of disorder but the interplay between the nonlinearity and discreteness.

These peculiar states have been experimentally verified as vibrational modes in low-dimensional crystals,  localized excitations in spin lattices,  and localized resistive states in Josephson junction arrays . The latter systems are of special interest because they have served for many years as well-controlled laboratory objects to study various nonlinear phenomena.  Moreover, at variance with other systems intrinsic localized modes found in Josephson coupled systems can be excited in the presence of time independent external driving forces.

A well known structure where dynamical localized states appear, is the anisotropic Josephson junction ladder (JJL) . A schematic view of such a ladder is given in Fig.1. The ladder contains small Josephson junctions indicated by crosses in Fig.1, in both longitudinal (vertical junctions) and transverse (horizontal junctions) directions to the dc bias current $`\gamma `$. The anisotropy of the ladder is due to the different sizes of vertical and horizontal junctions and is characterized by the anisotropy parameter $`\eta =\frac{I_{cH}}{I_{cV}}`$, where $`I_{cV}`$ and $`I_{cH}`$ are respectively the critical currents of the vertical and horizontal junctions.

Dynamical localized excitations persist in a JJL due to the intrinsic bistability property of a single small underdamped Josephson junction. One of these stable states is a superconducting one with zero voltage drop across the junction. The other state is a resistive one with a nonzero voltage drop, also called whirling state. A breather state in a JJL is characterized by a few junctions being in the resistive state, while the rest of all junctions are in the superconducting state. The presence of breather states can be verified by measurements of a total dc voltage drop across the ladder, which is used to plot current-voltage ($`IV`$) characteristics. This method does not provide spatially resolved information. It has been successfully combined with snapshots made using low temperature laser microscopy techniques,  which allow for a spatial resolution of the dc voltage drops. Note that both methods provide only with time-averaged voltage drop data, so the internal dynamics are so far not accessible in experiments.

However, the full dynamical picture is much more subtle than the time-averaged picture might suggest. In particular the Josephson junctions in the superconducting state exhibit small librations of the Josephson phase and correspondingly, nonzero ac voltage drops. The amplitude of these librations should decay to zero with increasing distance from the resistive junctions. It is also well known that various systems of coupled Josephson junctions support a delocalized class of excitations, namely, small amplitude electromagnetic waves (EWs).  In JJLs the spectrum of EWs consists of three branches and depends in a complex manner on the anisotropy $`\eta `$ and the dc bias $`\gamma `$. The resonant interaction of these EWs with the homogeneous whirling state (HWS) in the presence of an externally applied magnetic field has been studied in \[Ref. \]. It was shown that such an interaction leads to resonant steps in $`IV`$ curves and the voltage positions of resonant steps can be mapped onto the spectrum of EWs.

Early theoretical studies  dealt with the possibility of resonant interaction of breather states with EWs of the ladder. Due to the intrinsic spatial inhomogeneity of a breather state, EWs can be excited even in the absence of an externally applied magnetic field. The resonant interaction of EWs with the breather state manifests itself through the appearance of resonant steps  and various switchings between different breather states (”voltage jumps”) in $`IV`$ curves. Moreover, the amplitude of the Josephson phase librations at some distance from the breather center increases drastically due to this interaction.

The spatial inhomogeneity of a breather state allows also for the appearance of localized small amplitude electromagnetic excitations (EEs). We will show that the presence of localized EEs and their resonant interaction with the breather is of crucial importance. Our study resolves a long standing puzzle of the nature of breather instability. We show that most of these instabilities are driven by localized EEs. Especially the so-called retrapping (i.e. the switching from a breather to the superconducting state) is due to these instabilities and can not be explained using standard retrapping arguments.

In this paper we present a consistent theoretical study of resonant interactions between the breather states and EEs. We will derive the spectrum of EWs and calculate the ac power of oscillations at some distance from the breather center (”breather tail”). By making use of direct numerical simulations of the dynamics of JJLs in different parameters ranges, we demonstrate that the $`IV`$ curves display a variety of different resonance steps and voltage jumps as the external dc bias, and correspondingly the breather frequency, decrease. All these features are mapped onto various resonances of the breather frequency (or even its second and third harmonics) with EE frequencies, as well as with combinations of frequencies of different EEs. Moreover, by monitoring the power of ac oscillations as a function of the dc bias $`\gamma `$, we develop a spectroscopy of the EWs in the JJL.

The paper is organized as follows: In section II we derive the equations of motion within the framework of the Resistively Shunted Junction (RSJ) model  and obtain the spectrum of linear EWs of a JJL. A symmetry classification of different types of breathers is presented in section III. An improved analysis of the breather tail and the corresponding dependence of the power of ac oscillations on the breather frequency is given in section IV. In section V we will classify different resonances of breathers with EEs. Direct numerical simulations of the dynamics of breather states are presented and discussed in section VI. Section VII is devoted to the interpretation of the obtained resonances and switchings.

## II Dynamics of a JJL and the spectrum of linear electromagnetic waves

The complete dynamics of a JJL is determined by the time dependent Josephson phases of vertical $`\varphi _n^v`$, upper horizontal $`\varphi _n^h`$ and lower horizontal $`\stackrel{~}{\varphi }_n^h`$ junctions. The subscript $`n`$ labels the cell number. By making use of the RSJ model for each junction  we obtain the following set of equations:

$`𝒩(\varphi _n^v)=I_n^v`$ (1)
$`𝒩(\varphi _n^h)={\displaystyle \frac{1}{\eta }}I_n^h`$ (2)
$`𝒩(\stackrel{~}{\varphi }_n^h)={\displaystyle \frac{1}{\eta }}\stackrel{~}{I}_n^h,`$ (3)
where the nonlinear operator $`𝒩`$ is defined as

$`𝒩(y)=\ddot{y}+\alpha \dot{y}+\mathrm{sin}y.`$ (4)
Here, the unit of time is the inverse plasma frequency $`\omega _p^1`$, and the currents $`I_n^v`$, $`I_n^h`$ and $`\stackrel{~}{I}_n^h`$ are measured in units of the critical current of vertical junctions. The dimensionless parameter $`\alpha `$ determines the damping strength in each junction. Note that the positive current direction is chosen to be directed from bottom to top and from left to right. The currents flowing via the Josephson junctions and the Josephson phases are governed by the Kirchhoff laws

$`\gamma =I_n^v+I_n^hI_{n1}^h`$ (5)
$`\gamma =I_n^v\stackrel{~}{I}_n^h+\stackrel{~}{I}_{n1}^h,`$ (6)
and the flux quantization law in each cell:

$`\beta _LI_n^m=\varphi _n^h+\varphi _{n+1}^v\stackrel{~}{\varphi }_n^h\varphi _n^v.`$ (7)
Here, we introduced the mesh currents $`I_n^m`$ and the normalized inductance of the cell, $`\beta _L`$.

The Kirchhoff equations can be subtracted from each other yielding

$`I_n^h+\stackrel{~}{I}_n^h=C,`$ (8)
where $`C`$ is a constant for the whole ladder. This constant corresponds to the net difference between the currents flowing through the upper and lower horizontal junctions. For an open ladder of finite size, $`C`$ is zero. For a ladder of annular geometry with periodic boundary conditions $`C`$ may be nonzero and corresponds to the flux ”trapped” by the ladder ring. In the following we will consider the case of a finite open ladder with $`C=0`$. Then we may always eliminate the currents $`\stackrel{~}{I}_n^h`$ from the set of equations (1) as $`\stackrel{~}{I}_n^h=I_n^h`$. Since the junction width is larger than the London penetration depth, the mesh currents are

$`I_n^m=I_n^h,`$ (9)
and the currents flowing through the vertical junctions are expressed as

$`I_n^v=\gamma I_n^m+I_{n1}^m.`$ (10)
Inserting the relations (7), (9), and (10) into (1) we finally obtain the following set of coupled differential equations  :

$`\ddot{\varphi }_n^v+\alpha \dot{\varphi }_n^v+\mathrm{sin}\varphi _n^v`$ $`=`$ $`\gamma +{\displaystyle \frac{1}{\beta _L}}(\mathrm{}\varphi _n^v+\varphi _{n1}^h\stackrel{~}{\varphi }_{n1}^h)`$ (11)
$`\ddot{\varphi }_n^h+\alpha \dot{\varphi }_n^h+\mathrm{sin}\varphi _n^h`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(\varphi _n^v+\varphi _n^h\stackrel{~}{\varphi }_n^h)`$ (12)
$`\ddot{\stackrel{~}{\varphi }}{}_{n}{}^{h}+\alpha \dot{\stackrel{~}{\varphi }}{}_{n}{}^{h}+\mathrm{sin}\stackrel{~}{\varphi }_n^h`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(\varphi _n^v+\varphi _n^h\stackrel{~}{\varphi }_n^h),`$ (13)
where we use the notations $`\mathrm{}f_nf_{n1}2f_n+f_{n+1}`$ and $`f_nf_{n+1}f_n`$.

Next, we carry out the analysis of the delocalized class of excitations, namely of small amplitude electromagnetic waves (EWs). Note here that the spectrum of EWs $`\omega (q)`$ depends crucially on the state of the system. In the following we will consider a static (superconducting) state, $`\varphi _n^v=\mathrm{arcsin}\gamma `$ and $`\varphi _n^h=\stackrel{~}{\varphi }_n^h=0`$. We decompose the Josephson phases into the particular form:

$`\varphi _n^v=\varphi _n^v+\phi _n^v`$ (14)
$`\varphi _n^h=\varphi _n^h+\phi _n^h`$ (15)
$`\stackrel{~}{\varphi }_n^h=\stackrel{~}{\varphi }_n^h+\stackrel{~}{\phi }_n^h,`$ (16)
where $`\phi _n^v`$, $`\phi _n^h`$ and $`\stackrel{~}{\phi }_n^h`$ describe the small amplitude EWs. Substituting these expressions into system (11) and using the smallness of the amplitude of EWs we obtain

$`\ddot{\phi }_n^v+\alpha \dot{\phi }_n^v+\sqrt{1\gamma ^2}\phi _n^v`$ $`=`$ $`{\displaystyle \frac{1}{\beta _L}}(\mathrm{}\phi _n^v+\phi _{n1}^h\stackrel{~}{\phi }_{n1}^h)`$ (17)
$`\ddot{\phi }_n^h+\alpha \dot{\phi }_n^h+\phi _n^h`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(\phi _n^v+\phi _n^h\stackrel{~}{\phi }_n^h)`$ (18)
$`\ddot{\stackrel{~}{\phi }}{}_{n}{}^{h}+\alpha \dot{\stackrel{~}{\phi }}{}_{n}{}^{h}+\stackrel{~}{\phi }_n^h`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(\phi _n^v+\phi _n^h\stackrel{~}{\phi }_n^h).`$ (19)
In the weakly damped case as the parameter $`\alpha 1`$ we can derive the spectrum $`\omega (q)`$ of EWs neglecting effects of dissipation. By taking the Josephson phases $`\phi _n^v`$, $`\phi _n^h`$ and $`\stackrel{~}{\phi }_n^h`$ in the form

$`\left(\begin{array}{c}\phi _n^v\\ \phi _n^h\\ \stackrel{~}{\phi }_n^h\end{array}\right)=e^{i(qn+\omega t)}\left(\begin{array}{c}\mathrm{\Delta }_v\\ \mathrm{\Delta }_h\\ \stackrel{~}{\mathrm{\Delta }}_h\end{array}\right),`$ (26)
we find that the spectrum consists of three branches. The first is given by

$`\omega _0^2=1,`$ $`\mathrm{\Delta }_v=0,`$ $`\mathrm{\Delta }_h=\stackrel{~}{\mathrm{\Delta }}_h.`$ (27)
This branch is dispersionless and EWs corresponding to this branch are characterized by nonactive vertical junctions and in phase (symmetric) librations of the Josephson phases of upper and lower horizontal junctions.

The two other solutions are generalizations of those discussed in \[Ref. \], namely

$`\omega _\pm ^2`$ $`=`$ $`F\pm \sqrt{F^2G},`$ (28)
$`F`$ $`=`$ $`{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{\beta _L\eta }}+{\displaystyle \frac{1}{2}}\sqrt{1\gamma ^2}+{\displaystyle \frac{1}{\beta _L}}(1\mathrm{cos}q),`$ (29)
$`G`$ $`=`$ $`(1+{\displaystyle \frac{2}{\beta _L\eta }})\sqrt{1\gamma ^2}+{\displaystyle \frac{2}{\beta _L}}(1\mathrm{cos}q).`$ (30)
Both branches have a nonzero dispersion.

The branch $`\omega _+`$ is characterized by $`\mathrm{\Delta }_h=\stackrel{~}{\mathrm{\Delta }}_h`$ for all wave vectors $`q`$, i.e. the upper and lower horizontal phases are antisymmetric. The frequency range of the branch is above the degenerate branch $`\omega _0`$, i.e. $`\omega _+(q)>\omega _0`$ and it depends strongly on $`\beta _L`$. As the parameter $`\beta _L`$ increases, the width of $`\omega _+(q)`$ decreases and the branch approaches the dispersionless one, $`\omega _0`$. In the opposite case of small $`\beta _L`$, the frequencies $`\omega _+(q)`$ increase as $`1/\sqrt{\beta _L}`$. For zero wave number $`q=0`$, the amplitudes of EWs in this branch are characterized by $`\mathrm{\Delta }_v=0`$ and $`\mathrm{\Delta }_h=\stackrel{~}{\mathrm{\Delta }}_h`$, which means that only horizontal junctions are excited.

The branch $`\omega _{}`$ becomes dispersionless for the particular case of $`\gamma =0`$ and it corresponds to the dispersionless band obtained in \[Ref. \]. The frequency range of this branch is located below $`\omega _0`$ i.e. $`\omega _{}(q)<\omega _0(q)`$. For zero wave number $`q=0`$ the horizontal junctions are not active ($`\mathrm{\Delta }_h=\stackrel{~}{\mathrm{\Delta }}_h=0`$) and only vertical junctions are excited.

For a finite size ladder with open boundary conditions and $`N`$ cells, i.e. $`N+1`$ vertical junctions, the spectrum of linear waves is discrete and characterized by the following choice of allowed wave number values:

$$q_l=\frac{l\pi }{N+1},l=0,1,2,\mathrm{},N.$$
(31)
These EWs are cavity modes of the JJL. Odd values of $`l`$ correspond to antisymmetric eigenvectors (with respect to reflections at the center of the ladder), whereas even values correspond to symmetric ones.

Note that the above spectrum of EWs (28) is in general quite different from the EW spectrum of the homogeneous whirling state.  The latter can be obtained by choosing $`\gamma =1`$ in (28). The main difference is the appearance of a gapless (acoustic) lower branch.

## III Symmetries and dc bias dependent frequencies of breathers

In this Section we turn to the analysis of dynamic localized excitations (breathers) in JJLs. As mentioned in the introduction, breather states correspond to a few junctions being in the resistive state with all other junctions being in the superconducting state. The Josephson phases of resistive junctions are unbounded in time and the Josephson phases of superconducting junctions display small amplitude librations with a frequency $`\mathrm{\Omega }`$. This frequency is called the breather frequency.

Experiments  have revealed many different breather structures. All of them can be classified into three symmetry types using the reflection symmetries of the JJL. Some possible realizations are presented schematically in Fig. 2.

Breathers from the first group reveal an ’up-down’ reflection symmetry $`\widehat{S}_{ud}`$ (see, Fig.2a), i.e. they are invariant under exchange of upper and lower horizontal junctions. The second group consists of breathers invariant under a ’left-right’ reflection symmetry $`\widehat{S}_{lr}`$ (see, Fig.2b), i.e. they are invariant under a reflection at a vertical line cutting the ladder (this line is located either in the middle between two vertical junctions, or passes directly through one vertical junction). A third distinct group of breathers possesses an ’inversion’ symmetry $`\widehat{S}_{in}`$ (Fig.2c), i.e. these breathers are invariant under a reflection at a point which is either located on a vertical junction or in the center of a plaquette. A fourth group of breathers has no symmetries at all and does not belong to any of the three listed symmetry types. A particular example of a breather without symmetry is shown in Fig.2d. All of these types of breather excitations have been observed experimentally and numerically.  Each group of breathers can also have a different number of vertical junctions in the resistive state. Note that the particular example in Fig.2a possesses not only $`\widehat{S}_{ud}`$ symmetry, but also $`\widehat{S}_{lr}`$ and $`\widehat{S}_{in}`$ symmetries. However it is also possible to construct more complex breather states which display $`\widehat{S}_{ud}`$ symmetry only.

Next we derive the average voltage drop across a resistive vertical junction $`V=\dot{\varphi }^v`$ for different breather types. For the particular case of a breather with up-down symmetry $`\widehat{S}_{ud}`$, there are $`k`$ rotating vertical junctions in cells $`(i+1)..(i+k)`$ and two rotating horizontal junctions in the $`(i)`$th and $`(i+k)`$th cells respectively. By making use of the flux quantization law (7) we obtain identical voltage drops across the vertical junctions in the resistive state. For the same reasons the voltage drops across the horizontal junctions are two times smaller. Neglecting nonlinear contributions from the time average of $`\mathrm{sin}\varphi `$ on resistive junctions we obtain:

$$\begin{array}{ccc}\alpha \frac{V}{2}=\frac{1}{\eta }I_i^m\hfill & & \\ \alpha V=\gamma +I_{i+1}^mI_i^m\hfill & & \\ \mathrm{}\hfill & & \\ \alpha V=\gamma +I_{i+k}^mI_{i+k1}^m\hfill & & \\ \alpha \frac{V}{2}=\frac{1}{\eta }I_{i+k}^m.\hfill & & \end{array}$$
(32)
Thus, the voltage drop across a resistive vertical junction that corresponds to the experimentally measured voltage drop across the ladder, is given by

$$V=\frac{k\gamma }{\alpha (k+\eta )}.$$
(33)
Similarly we analyze a breather with left-right symmetry $`\widehat{S}_{lr}`$ and with inversion symmetry $`\widehat{S}_{in}`$. Taking into account that in these cases the voltage drops across resistive horizontal and vertical junctions are identical, we find

$$V=\frac{k\gamma }{\alpha (k+2\eta )}.$$
(34)
In a similar manner the result for a breather which has no symmetry (cf. Fig. 2d) reads

$$V=\frac{k\gamma }{\alpha (k+\frac{3}{2}\eta )}.$$
(35)
The above results for the dependence of the average voltage drop on the dc bias may be combined in a single expression

$$V=\frac{k\gamma }{\alpha [k+(3\frac{1}{2}\delta )\eta ]},$$
(36)
where $`k`$ is the number of vertical rotating junctions and $`\delta `$ denotes the number of resistive horizontal junctions. Note that $`\delta =4`$ for breathers with up-down symmetry, $`\delta =2`$ for left-right or inversion symmetry, and $`\delta =3`$ for no symmetry.

In order to analyze the interaction between the breather state and the linear EEs we need to know the frequency $`\mathrm{\Omega }`$ of a breather solution. Noting that the breather frequency is given by the lowest realized voltage drop across a resistive junction, we find

1) up-down symmetry

$$\mathrm{\Omega }=\frac{k\gamma }{2\alpha (k+\eta )},$$
(37)
2) left-right symmetry and inversion symmetry

$$\mathrm{\Omega }=\frac{k\gamma }{\alpha (k+2\eta )},$$
(38)
3) no symmetry

$$\mathrm{\Omega }=\frac{k\gamma }{\alpha (2k+3\eta )}.$$
(39)
## IV Spatial tails of breathers

We consider the spatial dependence of Josephson phases in the presence of a breather state. At some distance from the breather center (”breather tail”) the Josephson phases librate with small amplitudes. In order to analyze the breather dynamics in the tail, we use the linearized system of equations (17). Keeping in mind the time periodicity of the breather solution, the librating Josephson phases take the form (for $`n<0`$, and with the breather center as the origin)

$`\left(\begin{array}{c}\phi _n^v\\ \phi _n^h\\ \stackrel{~}{\phi }_n^h\end{array}\right)=e^{\lambda n+i\mathrm{\Omega }t}\left(\begin{array}{c}\mathrm{\Delta }_v\\ \mathrm{\Delta }_h\\ \stackrel{~}{\mathrm{\Delta }}_h\end{array}\right).`$ (46)
Introducing new sum and difference variables of the amplitudes of Josephson phases of horizontal junctions

$$\mathrm{\Delta }_h^{}=\frac{1}{2}(\mathrm{\Delta }_h\stackrel{~}{\mathrm{\Delta }}_h),\mathrm{\Delta }_h^+=\frac{1}{2}(\mathrm{\Delta }_h+\stackrel{~}{\mathrm{\Delta }}_h),$$
(47)
it follows that

$`(A{\displaystyle \frac{2}{\beta _L}}\mathrm{cosh}\lambda )\mathrm{\Delta }_v{\displaystyle \frac{2}{\beta _L}}(1e^\lambda )\mathrm{\Delta }_h^{}=0`$ (48)
$`{\displaystyle \frac{1}{\beta _L\eta }}(1e^\lambda )\mathrm{\Delta }_v+B\mathrm{\Delta }_h^{}=0`$ (49)
$`(\mathrm{\Omega }^2+i\alpha \mathrm{\Omega }+1)\mathrm{\Delta }_h^+=0,`$ (50)
where the frequency dependent parameters $`A`$ and $`B`$ are

$`A`$ $`=`$ $`\mathrm{\Omega }^2+i\alpha \mathrm{\Omega }+\sqrt{1\gamma ^2}+{\displaystyle \frac{2}{\beta _L}}`$ (51)
$`B`$ $`=`$ $`\mathrm{\Omega }^2+i\alpha \mathrm{\Omega }+1+{\displaystyle \frac{2}{\beta _L\eta }}.`$ (52)
From (48) we immediately find that $`\mathrm{\Delta }_h^+=0`$ and hence $`\mathrm{\Delta }_h^{}=\mathrm{\Delta }_h`$. It follows that $`\phi _n^h=\stackrel{~}{\phi }_n^h`$, i.e. breather tails appear with perfect up-down symmetry. This is at variance with the complex symmetry properties of the resistive breather center (see the Section III). A nontrivial solution to the first two equations in (48) exists if

$$AB\frac{2}{\beta _L}B\mathrm{cosh}\lambda \frac{4}{\beta _L^2\eta }+\frac{4}{\beta _L^2\eta }\mathrm{cosh}\lambda =0.$$
(53)
The dependence of the complex parameter $`\lambda `$ on the breather frequency $`\mathrm{\Omega }`$ is given by

$`\lambda =\mathrm{ln}(z+\sqrt{z^21}),`$ (54)
with

$`z={\displaystyle \frac{4\beta _L^2\eta AB}{42\beta _L\eta B}}.`$ (55)
Note here that this expression can be obtained directly from (28) by assuming that $`q=i\lambda `$ and substituting $`\mathrm{\Omega }^2i\alpha \mathrm{\Omega }`$ instead of $`\omega _\pm ^2`$.

The real $`\mathrm{}(\lambda )`$ and imaginary $`\mathrm{}(\lambda )`$ parts of $`\lambda `$ determine respectively the spatial decay and spatial period of oscillations of Josephson phases in the breather tail. Moreover, $`\mathrm{}(\lambda )`$ and $`\mathrm{}(\lambda )`$ strongly depend on the breather frequency that in turn, can be changed by varying the external dc bias $`\gamma `$. In Fig.3(a,b) we plot the real and imaginary parts of $`\lambda `$ for three breather types (cf. Fig. 2) versus $`\gamma `$. The minima of the real part of $`\lambda `$ correspond to resonances with linear EWs $`\omega _+(q)`$ and $`\omega _{}(q)`$. As one can see, the current positions of these minima shift for different types of breathers. Note that the $`\gamma `$ range in the plots extends from zero to one. In fact the breather will exist only in a narrower current region. This is due to the presence of both finite nonzero retrapping current  (switching to the superconducting state) and a particular current $`\gamma <1`$, where the breather switches to the HWS.

Since the EW frequencies $`\omega _+(q)`$ decrease with increasing inductance of the cell (28), the position of the global minimum also depends strongly on $`\beta _L`$ (see Fig.4(a)).

Although the dynamics of the Josephson phases in the breather tail are completely determined by the dependence of the parameter $`\lambda `$ on the breather frequency, it may be more convenient to measure just the time-average power $`P_{ac}`$ of the libration of a junction at the edge of a JJL. We show below that the monitoring of $`P_{ac}`$ upon decreasing the dc bias allows one to develop a spectroscopy of EWs. The value of $`P_{ac}`$ is determined by the average kinetic energy of the edge vertical junction

$`P_{ac}={\displaystyle \frac{1}{2}}\dot{\varphi }_l^{v^2},l={\displaystyle \frac{N}{2}}`$ (56)
and is derived in Appendix A. The typical dependences of $`P_{ac}`$ on the dc bias and the inductance of the cell for different types of breathers are presented in Figs. 3c, 4b.

## V Breather instabilities and classification of resonances

In early experiments on breathers in JJLs  it was observed that breather states may switch to other breather states or HWSs upon variation of the dc bias $`\gamma `$. These switchings can be either due to the disappearance of the breather state as a solution of the underlying dynamical equations or due to the novel effect of dynamical instability of the breather state. In the latter case the solution continues to exist, but turns unstable, forcing the system to search for another stable attractor state. The cause of such instabilities is the resonant interaction of the breather state with small amplitude (localized or delocalized) EEs.

In order to analyze the breather instabilities in more precise terms, we have to linearize the phase space flow around a given breather solution and study the corresponding Floquet eigenvalue problem (see Appendix B). The outcome of this analysis is a spectrum of Floquet multipliers (i.e. eigenvalues) $`\mu `$. In most cases the multipliers reside on a circle of radius

$`R(\alpha )=e^{\frac{\alpha \pi }{\mathrm{\Omega }}},`$ (57)
in the complex plane, which is less than 1. Thus, the multipliers can be expressed in the general form

$`\mu =R(\alpha )\mathrm{exp}\left(\pm i\omega {\displaystyle \frac{2\pi }{\mathrm{\Omega }}}\right),`$ (58)
where $`\mathrm{\Omega }`$ is the breather frequency and $`\omega `$ is some characteristic EE frequency. The corresponding eigenvectors may be divided into two classes, namely those which are localized on the breather and those which are delocalized. While it is more involved to make analytical predictions for the class of localized Floquet eigenstates, we may immediately proceed with the characterization of the class of delocalized eigenstates. Although the concrete form of the delocalized eigenvectors has to be obtained numerically, the fact that they are delocalized allows one to determine their frequency $`\omega `$. Indeed, since the breather is itself a localized state, delocalized Floquet eigenstates simply correspond to the above discussed linear EWs. Their frequencies have been derived in Section II, and thus we can reconstruct the delocalized part of the Floquet spectrum using them.

Stable breathers are characterized by all Floquet multipliers being located inside or on the unit circle. Instabilities occur after collisions of multipliers on the inner circle with radius (57) and a subsequent detaching from this circle towards larger absolute values. Although in the dissipative case an additional change of the control parameter (dc bias $`\gamma `$) is necessary in order for the corresponding Floquet eigenvalue(s) to cross and escape the unit circle, we may still look for collision conditions and classify possible instability scenarios.

Note that for Hamiltonian systems any collision and detaching from the unit circle leads to an instability. Finite (even though possibly weak) dissipation may drastically change the instability patterns by selecting the strong instabilities (large detachments) over the weak instabilities (small detachments). This is exactly what we observe in our numerical studies.

By applying the general stability analysis of nonlinear discrete systems  we obtain three possibilities for multiplier collisions. The first is realized when the collision takes place on the positive real axis in the complex plane. This implies that a multiplier is colliding with its complex conjugate partner. With (58) it follows that

$`\omega =m\mathrm{\Omega },`$ (59)
for any integer number $`m`$. These are primary resonances of the breather frequency or its higher harmonics with any of the frequencies of the EEs.

The collision of an eigenvalue with its complex conjugated partner can also take place on the negative real axis in the complex plane:

$`\omega =({\displaystyle \frac{1}{2}}+m)\mathrm{\Omega }.`$ (60)
These are parametric resonances of the breather state with EEs .

The third case is realized when the collision takes place away from the real axis. Then a multiplier has to collide with a different one, but not with its own complex conjugate partner. It follows that

$`\omega _1\pm \omega _2=m\mathrm{\Omega }.`$ (61)
This is a combination resonance, as the breather frequency (or its multiple) has to match a sum (difference) of the frequencies of two different EE frequencies.

We stress again that EEs may be localized or delocalized. While we may proceed with an analytical prediction of collisions using delocalized Floquet multipliers, we have to use numerical calculations to observe collisions involving localized ones. Especially the combination resonance may involve either two delocalized, two localized or one delocalized and one localized Floquet eigenvalues. In addition the above mentioned dissipation-induced selection of weak and strong instabilities will result in some possible collisions being harmless (leaving the breather stable) while others will turn out to be important for understanding breather instabilities.

## VI Numerical simulations of breather dynamics

To study the breather dynamics, we performed direct numerical simulations of the set of equations (11). All simulations were carried out for JJLs with $`N=10`$ cells. We impose open boundary conditions and use the 4th order Runge-Kutta method. Time is measured in dimensionless time units. The initial value of the dc bias was $`\gamma =0.8`$. We choose proper initial conditions that lead to the relaxation of the system into a particular breather state of left-right symmetry with one resistive vertical junction, as in Fig.2b . After a waiting time of 500 time units we use the next 500 time units to calculate the time averaged characteristics of the state. We then decrease the dc bias $`\gamma `$ by a tiny step of $`\mathrm{\Delta }\gamma =0.0005`$ and repeat the procedure. We checked that our results do not change upon further increase of the waiting time. We varied the anisotropy $`\eta `$ and the inductance of the cell $`\beta _L`$ while the dissipation $`\alpha =0.1`$ was fixed. We will comment on hysteresis effects due to additional increasing of the current in the conclusion.

There are three different ways to monitor the simulations. The first one is the $`IV`$ characteristics, that is the dependence of the averaged voltage drop across the resistive vertical junction on the dc bias. Furthermore we obtain the power $`P_{ac}`$ (56) of ac oscillations of the vertical junction at the edge of JJLs. Finally we generate time-resolved images (movies) of the full dynamical behavior of the ladder so that we can visually check whether the system still resides in the initially chosen breather state, or switches into another state. Our results are presented in Figs. 5-11, where each figure consists of two parts. In the left hand parts the $`IV`$ characteristics are shown (solid lines) together with the approximate results from section III (dotted lines). The vertically oriented dashed lines indicate the band edges of the linear EWs. In the right hand part of the figures we show the dependence of $`P_{ac}`$ (solid lines) on the dc bias together with our approximate analysis from section IV and Appendix A (dashed lines) where appropriate.

We start with the case of small $`\beta _L`$ values. For $`\beta _L=0.2`$ and $`\eta =1.15`$ (Fig.5) the breather is easily excited, and its frequency is located below $`\omega _+(q)`$. In the $`P_{ac}`$ plot we observe peaks that are due to the resonance of the second and third harmonics of the breather with $`\omega _+(q)`$. These resonances are primary ones as discussed above. Note that their presence is barely seen on the $`IV`$ curve. The series of observed peaks is related to the finite size of the system, and therefore to the resonant interaction of the breather with a discrete set of cavity modes as was discussed at the end of section II. We tested our interpretation by increasing the size of the system and observed the predicted increase in the number of resonance peaks. Close to the lowest possible current (around $`\gamma =0.55`$) we observe a switching to another breather state, which however has the same symmetry and spatial structure. Note that shortly after this switching (upon further lowering of the current) we lose the breather and the system switches to the superconducting ground state.

For lower values of the anisotropy $`\eta =0.35`$ (Fig.6) the resonances are again not detectable in the $`IV`$ curve. However, by monitoring $`P_{ac}`$ we observe the singularities that correspond to the primary resonance $`2\mathrm{\Omega }=\omega _+(q)`$ . Moreover, at the dc bias $`\gamma 0.35`$ we detect a weak third-order primary resonance $`3\mathrm{\Omega }=\omega _+(q)`$ in the breather tail. The dashed line in the right part of the figure is the prediction of $`P_{ac}`$ using our approximate tail analysis. Note that our approximate tail analysis is based on the assumption of a dense spectrum of EWs. Consequently the calculated $`P_{ac}`$ presents an envelope of the numerically observed series of discrete peaks.

Next we increased the inductance of the cell to $`\beta _L=0.5`$. In Fig.7 we show the results for $`\eta =1.15`$. For the initial value of the dc bias $`\gamma =0.8`$ the breather frequency is already located inside the $`\omega _+(q)`$ band of EWs, and this primary resonance is observed in the $`IV`$ curve. Indeed, the slope of the $`IV`$ curve is larger than the prediction (38) which does not take into account resonant interactions with EEs. With decreasing dc bias, the breather frequency is lowered and the above primary resonance disappears. However at lower current values the next primary resonances $`2\mathrm{\Omega }=\omega _+(q)`$ occur and are observable, both in the breather tail and in the $`IV`$ characteristics.

The primary resonance structures ($`\mathrm{\Omega }=\omega _+(q)`$) in the large current domain are also observed for smaller values of the anisotropy parameter $`\eta =0.5`$ (Fig.8). In this case, they manifest themselves through resonant steps in the $`IV`$ curve.  At lower values of the dc bias we again observe primary resonances with $`m=2`$.

For $`\beta _L=1.0`$ and $`\eta =0.5`$ (Fig.9) the breather frequency is located above $`\omega _+(q)`$ for large current values. Upon decreasing the dc bias we observe a peculiar switch to a different breather state with the same spatial structure but a lower frequency located inside $`\omega _+(q)`$. The most interesting feature here is that shortly before the switching the breather frequency is clearly larger and outside of the $`\omega _+(q)`$ region. Upon further lowering of the dc bias we observe primary resonances $`\mathrm{\Omega }=\omega _+(q)`$, and corresponding resonant steps in $`IV`$ curve.

Let us increase the inductance of the cell $`\beta _L`$ even further. For $`\beta _L=3`$ and $`\eta =0.5`$ (Fig.10) we again find that breather frequencies are located above $`\omega _+`$ for large current values. Similar to the previous case we observe a switching when the breather frequency is clearly outside (above) the branch $`\omega _+(q)`$. This brings the system into another breather state with the same spatial structure, but with a frequency again located above $`\omega _+`$. This highly nonlinear state is then lost by switching to the superconducting ground state after further decrease of the dc bias.

Now we come to an interesting observation. Lowering the anisotropy $`\eta =0.35`$ (Fig.11) we again observe the switching at a breather frequency being located above $`\omega _+(q)`$. However the switching increases the voltage drop. The new state is of different internal structure. We remind the reader that all previous numerical results have been obtained for a breather with a structure as in Fig.2b. Here we find that after the switch the new breather state is characterized by three vertical junctions being in the resistive state. At the same time the symmetry is broken. In fact the new state exactly corresponds to the example given in Fig.2d. Note that similar switchings (which lead to an increase of the number of resistive junctions) have been reported in early experimental studies.  The left-right symmetry breaking leads to new interesting features in the breather tails. The no-symmetry breather has two times lower frequency than the voltage drop across the vertical junction. Thus the new breather frequency is inside the upper band $`\omega _+(q)`$, and a primary resonance is clearly observed in the $`P_{ac}`$ dependence on $`\gamma `$.

In order to test the influence of small fluctuations on the $`IV`$ curves we repeated the simulations in the presence of small noise with amplitude $`10^8`$. All obtained results are stable except the switching outcome in Fig.11 . While this switching occurs at the same current value, the new breather state is changed. In particular we observed the left-right symmetry breather with three junctions being in the resistive state. Thus we find extreme sensitiveness of the outcoming breather structure (including its symmetries) to small fluctuations. This implies that the boundaries of the volumes of attraction of different (breather) attractors are entangled in a very peculiar way. We can reliably predict the switching position, but not the outcome of the switching.

We also numerically simulated the breather dynamics in the JJL with an extremely large inductance of the cell $`\beta _L=500`$. We did not find any indication of resonances and instabilities. We argue that the reason for that is the weak dispersion of the linear EWs for such large values of $`\beta _L`$. This implies that interactions along the ladder are weak. The breather is continued to small current values until it switches to the superconducting state at the dc dc bias $`\gamma =0.22`$ for $`\eta =0.35`$. This particular value can be obtained by making use of the simple dc analysis (34) and the standard theory of the retrapping current in a single small Josephson junction

$$\gamma _r=(1+2\eta )\frac{4\alpha }{\pi }.$$
(62)
This equation yields a value of 0.22 for the considered case, in good agreement with the numerical observation. Note that within this theory retrapping occurs purely due to energy considerations, not due to resonances (or instabilities).

It is very important to notice that for cases with small or intermediate values of inductance of the cell (Fig.5-11) the observed currents at which we lose the breather state and switch to the superconducting one exceed the expected retrapping values (62). We will explain this disagreement in the next section.

Motivated by the above findings we investigated the loss of the HWS upon lowering the current. We recall that in this state all vertical junctions are resistive and all horizontal ones are superconducting. Usually it is assumed that the HWS loss is again due to a standard retrapping mechanism. It is important that any numerical simulation of such a process is done with the addition of some weak noise, because the processor will otherwise perform a perfect simulation of a single junction repeated $`N+1`$ times. We chose $`\beta _L=3`$ and $`\eta =0.35`$. The $`IV`$ characteristic is shown by a thick dashed line in Fig.11 . The expected retrapping current $`4\alpha /\pi =0.127`$ is clearly not reached. Instead we observe the loss of the HWS at $`\gamma =0.273`$. At the same time it follows from equations (11) that the HWS exists as a solution down to the retrapping current of a single junction, i.e. down to $`\gamma =0.127`$ ! So in this case we conclude that the numerically observed loss of the HWS at $`\gamma =0.273`$ is due to an instability. The HWS continues to exist as a solution down to the standard retrapping current, but it is an unstable state. This novel result is very important, since very often the current value of HWS loss in the absence of a magnetic field is used to estimate different parameters of the system assuming that the HWS is behaving similar to a single junction. Our results show that this is definitely not the case.

## VII Evaluation of resonances and EWs spectroscopy

This section is devoted to a quantitative explanation of the observed resonances and switchings.

### A Primary resonances

Primary resonances are characterized by $`m\mathrm{\Omega }=\omega `$ where $`\omega `$ is some EE frequency. We detected various primary resonances with extended EWs. The case $`m=1`$, which corresponds to the breather frequency being located inside the $`\omega _+(q)`$ band, shows up with resonant steps in the $`IV`$ curves (see, Figs. 8-9). The finite number of observed resonant steps is due to the discrete spectrum of the excited cavity modes. In addition we observe strong variations of the breather tail amplitudes.

Higher order primary resonances ($`m=2,3`$) are much less pronounced in the $`IV`$ characteristics. They mainly lead to a weaker localization of the breather tail and can be clearly detected in the form of sharp peaks in the $`P_{ac}(\gamma )`$ dependence. Since the breather in our case has left-right symmetry, the only linear cavity modes which can be exited are symmetric (see section II). These modes are characterized by even values of $`k`$ (31). We start our evaluation of these resonances with the case shown in Fig.6. We determine the $`\gamma `$ value of each observed peak in $`P_{ac}(\gamma )`$ and thus obtain the corresponding breather frequency $`\mathrm{\Omega }(\gamma )`$. We then compare its multiples with the discrete spectrum of linear mode frequencies of the $`\omega _+`$ branch. The numbers are listed in Table I. We find that all observed resonances are due to symmetric linear modes (even $`k=2,4,6,8,10`$) as expected.

The same method of analysis allows us to conclude that the three peaks in Fig.5 (in decreasing order of dc bias) are due to the following resonances: $`(m,k)=(2,6);(2,4);(3,8)`$. Similarly the shoulder and the peak in Fig.6 are due to resonances with $`(m,k)=(2,10);(2,8)`$. Finally the resonances in Fig.8 correspond to the values $`(m,k)=(2,8);(2,6);(2,4)`$. Note that in all of these cases the deviations between the theoretical and observed numbers are much less than the frequency difference between adjacent cavity modes.

### B Parametric resonances

So far we did not comment on the nature of the switching from a breather state to the superconducting state for small and intermediate values of $`\beta _L`$. The Floquet analysis results show that all these switchings are due to an instability of the breather. In terms of Floquet multipliers, all of these instabilities are due to a collision of two localized Floquet multipliers on the negative real axis. The breather state continues to exist as a solution to the dynamical equations for lower current values, but it is unstable. Note that the so-called retrapping mechanism instead (as for a single junction) uses the critical current value as a criterion for retrapping. This argument is based purely on energy considerations and does not take into account any resonance mechanism. This is not surprising, as a single junction has no other degrees of freedom it may resonate with. Below the retrapping current, the resistive state disappears in this case. So we may state that the switching from a breather state to the superconducting one as observed in our simulations is usually driven by resonances with localized EEs (frequency matching) and is not due to energy effects (current value matching).

### C Combination resonances

Let us discuss the nature of the switchings of the breather for intermediate $`\beta _L`$ values when the breather frequency is located above the branch $`\omega _+(q)`$. These switchings are again due to an instability. It is characterized by Floquet multipliers colliding away from the real axis. As discussed in Section V this corresponds to a combination resonance. The numerical Floquet analysis shows that one of the two participating multipliers is a localized one (which bifurcates from the lower branch $`\omega _{}(q)`$)., while the second one belongs to the delocalized spectrum of $`\omega _+(q)`$. The Floquet multiplier which finally leaves the unit circle is a localized one. So again the instability of the breather is driven by a localized perturbation.

In Fig.12 the dependence of the arguments and absolute values of all relevant Floquet multipliers is shown for the breather of Fig.11. For convenience we do not plot the complex conjugate multipliers and restrict the arguments to $`0arg(\mu )\pi `$. The narrow band $`\omega _{}(q)`$ and broad band $`\omega _+(q)`$ are nicely observed. The degenerate band $`\omega _0`$ is located slightly above $`\omega _{}(q)`$. This band does not interact with other multipliers when crossing them, as expected from our analytical considerations. The two separated arguments which are located below the $`\omega _{}(q)`$ band have localized eigenvectors.

In the plot of the absolute values we observe the predicted values $`\mu =1`$ and $`e^{\alpha T_b}`$. The multipliers which correspond to lines between these two states generally reside on the circle with radius (57). Many of them depart from this circle due to collisions. At current values of $`0.55<\gamma <0.7`$ we observe parametric resonances $`2\mathrm{\Omega }=\omega _+(q)`$, which belong to the set of weak resonances and do not evolve into a global instability. However it is possible that a slight variation of control parameters (e.g. decreasing the damping $`\alpha `$) might change these resonances into strong ones. Then we would expect sudden instabilities of the breather state at these large current values. In our case the global instability is realized when one of the localized multipliers collides with the $`\omega _+(q)`$ band around $`\gamma =0.48`$. Subsequent lowering of the current leads to a fast escape of this multiplier from the unit circle and to the observed switching.

The importance of localized EEs for the destabilization of a breather is simply due to the localized nature of the latter. It is hard (if not impossible) for a breather to generate a parametric instability through extended EWs alone, as these excitations are damped out far from the breather center. In contrast localized EEs do not travel away from the breather center. These modes can be effectively excited by the breather, leading to an instability of the latter.

To understand the nature of the observed instability of the HWS we show a similar Floquet multiplier plot in Fig.13.

As in the previous plot, we observe weak parametric resonances of the upper EW band at current values $`\gamma 0.38`$ which do not evolve into a global instability. Again the observed instability is driven by a combination resonance at $`\gamma =0.175`$. Since the HWS is an extended state, all Floquet multipliers are also extended. The combination resonance is due to the collision of two Floquet multipliers belonging to the two EW branches $`\omega _\pm (q)`$ with $`q=\pi `$ (we remind the reader that the EW spectrum of the HWS is different from that of the superconducting state and can be obtained by putting $`\gamma =1`$ in (28)). Note that indeed for the present case, the frequency of the HWS at the instability equals $`2.73`$, while the value of the combination $`\omega _{}(\pi )+\omega _+(\pi )=2.69`$. For current values $`\gamma 0.18`$ the HWS is becoming stable again. However around $`\gamma =0.17`$ another even stronger instability due to parametric resonance sets in, which brings the HWS to the next instability well above the expected retrapping current.

## VIII Conclusion

We have presented analytical and numerical studies of breather properties in Josephson junction ladders. Our results confirm and substantially extend early suggestions that breathers may resonate in different ways with localized and extended electromagnetic modes. The numerical studies have been done in a parallel manner to the way experiments are conducted. The variation of the control parameter $`\gamma `$ allows one to continuously change the breather frequency, whereas the linear mode spectrum is not significantly changed in the domain of interest. We observed primary resonances $`m\mathrm{\Omega }=\omega _+(q_k)`$ with extended EWs, parametric resonances $`\mathrm{\Omega }=2\omega `$ with localized EEs, and combination resonances $`m\mathrm{\Omega }=\omega +\omega _+(q_1)`$ with a localized EE and a delocalized EW participating. We also observed a combination resonance which leads to a switching from a small breather (one resistive vertical junction) to a larger one (three resistive vertical junctions) together with a possible symmetry lowering of the breather.

The primary resonances with extended EWs lead to singularities in the breather tails. This allows one to develop a spectroscopy of EWs by monitoring $`P_{ac}`$ versus $`\gamma `$. Such a spectroscopy can be experimentally realized, e.g. with the help of a well known Josephson junction detector technique.  It could be important for obtaining a coherent source of high frequency radiation, since in such a resonance the whole breather tail starts to coherently oscillate with large amplitudes. The resistive breather center serves as a region of energy input via a dc bias.

Our studies show that the main control parameter (in addition to the dc bias) is the self inductance $`\beta _L`$. For small $`\beta _L`$ values the breather frequency is located between the two branches of EWs, $`\omega _\pm (q)`$. One may perform spectroscopy of EWs of the upper branch, or observe parametric instability of a breather due to localized modes. Moreover as $`\beta _L`$ increases, resonant steps in the $`IV`$ characteristics can be observed. For intermediate $`\beta _L>1`$ values, the breather frequency is located above $`\omega _+(q)`$. In this case we observe combination resonances due to localized and delocalized EEs, which may result in a unusual sensitivity of the switching outcome on small fluctuations. Large $`\beta _L`$ values stabilize the breather states, make resonances impossible and lead to a standard retrapping mechanism for breather switching to the superconducting ground state. This is likely the situation for the reported experimental data in \[Ref. \]. However lower $`\beta _L`$ values allow for the appearance of the above listed resonances (and perhaps even other still unobserved resonances). We believe that our findings will help to make the proper parameter choice when designing new ladders for experiments.

Note that throughout our studies we always decrease the dc bias $`\gamma `$. Let us consider a breather state $`Nr.1`$ which becomes unstable upon lowering the dc bias at a certain value $`\gamma _1`$. Let this be a case where the system will switch to another stable breather state $`Nr.2`$. This new breather state is in fact keeping its stability not only upon further lowering of the dc bias, but also upon a reversing (increasing) of the dc bias. Thus we find that there exist dc bias windows in which both the starting breather state $`Nr.1`$ and the new breather state $`Nr.2`$ are stable. Even though their dc spatial structures (cf. Fig.2) may be identical, the average voltage drops (and frequencies) are in general different. Further increasing of the dc bias while staying on the new breather state $`Nr.2`$ will lead to an instability and switching at $`\gamma _2>\gamma _1`$. In case the switching brings the system back to the breather state $`Nr.1`$ we are faced with the well-known two-state hysteresis phenomenon in Josephson junction systems. However we also observed cases when the switching due to an instability of breather state $`Nr.2`$ (upon increasing the dc bias) brings the system to yet another state, which differs from breather state $`Nr.1`$, e.g. simply to the HWS. In such a case the hysteretic behavior is of more complex nature. To keep the discussion of our studies as clear as possible we did not present data for increasing current.

Another observation is that the expected values of the retrapping current based on a pure dc analysis  are too low to match the observed values at which the breather switches to the superconducting ground state. Only for very large $`\beta _L`$ values do we observe agreement. For all other cases the breather switches to the superconducting state via an instability driven by parametric resonance with localized EEs. Moreover the HWS also undergoes an instability which is due to combination resonances with extended EWs.

In this work we always started with a breather configuration as in Fig.2b, at large current values. It may be expected that the results for other starting configurations will show up with similar properties, and perhaps with additional types of resonances as well. This may be due to the fact that the structure of the phase space is very complex, being separated in many different regions of attraction of different attractors. It is this complexity which makes the understanding of breather properties both a fascinating and complicated undertaking.
Acknowledgements
We thank A. Benabdallah, P. Binder, M. Schuster and A. V. Ustinov for valuable discussions and F. V. Kusmartsev and J. J. Mazo for sending us their preprints prior publication. This work was supported by the Deutsche Forschungsgemeinschaft and by the European Union under the RTN project LOCNET HPRN-CT-1999-00163.

## A Power of ac librations at the edge of a JJL

Here we derive the time-average power of ac librations at the edge of a JJL. This characteristics is proportional to the average kinetic energy $`\dot{\varphi }^2/2`$. In order to obtain an expression for the kinetic energy, we have to determine the dynamics of the junction at the JJL edge. For this we write the system of equations (48) in a matrix form:

$`\widehat{A}\stackrel{}{v}=0.`$ (A1)
Here, $`\stackrel{}{v}`$ is an unknown vector

$`\stackrel{}{v}=\left(\begin{array}{c}\mathrm{\Delta }_v\\ \mathrm{\Delta }_h\end{array}\right),`$ (A4)
and $`\widehat{A}`$ is a $`2\times 2`$ matrix

$`\widehat{A}(\lambda )=\left(\begin{array}{cc}a_{11}(\lambda )& a_{12}(\lambda )\\ a_{21}(\lambda )& a_{22}(\lambda )\end{array}\right),`$ (A7)
where $`a_{11}(\lambda )=A\frac{2}{\beta _L}\mathrm{cosh}\lambda `$, $`a_{12}(\lambda )=\frac{2}{\beta _L}(1e^\lambda )`$, $`a_{21}(\lambda )=\frac{1}{\beta _L\eta }(1e^\lambda )`$, and $`a_{22}(\lambda )=B`$. A nonzero solution exists if the determinant of $`\widehat{A}`$ vanishes: $`det\widehat{A}(\lambda _0)=0`$. The parameter $`\lambda _0`$ was determined in section IV (54).

The components of the vector $`\stackrel{}{v}`$ satisfy the condition:

$`\mathrm{\Delta }_v={\displaystyle \frac{a_{22}}{a_{21}}}\mathrm{\Delta }_h.`$ (A8)
To determine the components $`V`$ and $`H`$ separately we have to impose an additional condition at the breather center. This condition is not known exactly due to the complex dynamics in the resistive breather center. Nevertheless the ac librations in the breather tails are weakly depending on it. Here we use the simplest normalization condition

$`|\mathrm{\Delta }_v|^2+|\mathrm{\Delta }_h|^2=1.`$ (A9)
Substituting (A8) into (A5) we obtain

$`\mathrm{\Delta }_h`$ $`=`$ $`{\displaystyle \frac{|a_{21}|^2}{\sqrt{|a_{21}|^2+|a_{22}|^2}}}`$ (A10)
$`\mathrm{\Delta }_v`$ $`=`$ $`{\displaystyle \frac{a_{22}\overline{a}_{21}}{\sqrt{|a_{21}|^4+|a_{22}|^2|a_{21}|^2}}}.`$ (A11)
Due to the up-down symmetry in the breather tail, the dynamics of Josephson phases at the edge of the JJL can be written in the form

$`\phi _n^v`$ $`=`$ $`{\displaystyle \frac{\mathrm{}(a_{22}\overline{a}_{21})e^{\mathrm{}(\lambda _0)n}\mathrm{cos}(\mathrm{}(\lambda _0)n+\mathrm{\Omega }t)}{\sqrt{|a_{21}|^4+|a_{22}|^2|a_{21}|^2}}}`$ (A13)
$`+{\displaystyle \frac{\mathrm{}(a_{22}\overline{a}_{21})e^{\mathrm{}(\lambda _0)n}\mathrm{sin}(\mathrm{}(\lambda _0)n+\mathrm{\Omega }t)}{\sqrt{|a_{21}|^4+|a_{22}|^2|a_{21}|^2}}},`$
$`\phi _n^h`$ $`=`$ $`{\displaystyle \frac{|a_{21}|^2e^{\mathrm{}(\lambda _0)n}\mathrm{cos}(\mathrm{}(\lambda _0)n+\mathrm{\Omega }t)}{\sqrt{|a_{21}|^2+|a_{22}|^2}}}.`$ (A14)
We finally obtain the expression for the average kinetic energy (for $`n<0`$)

$`{\displaystyle \frac{1}{2}}\dot{\phi }_n^{v^2}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }^2|a_{22}\overline{a}_{21}|^2e^{\mathrm{}(\lambda _0)2n}}{4(|a_{21}|^4+|a_{22}|^2|a_{21}|^2)}}.`$ (A15)
## B Linear stability of the breather in JJL.

The stability of periodic motion is analyzed with the help of the Floquet theory.  Linearizing the system (11) around a time periodic breather solution, we obtain

$`\ddot{ϵ}_n^v+\alpha \dot{ϵ}_n^v+A_n^v(t)ϵ_n^v`$ $`=`$ $`{\displaystyle \frac{1}{\beta _L}}(\mathrm{}ϵ_n^v+ϵ_{n1}^h\stackrel{~}{ϵ}_{n1}^h)`$ (B1)
$`\ddot{ϵ}_n^h+\alpha \dot{ϵ}_n^h+A_n^h(t)ϵ_n^h`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(ϵ_n^v+ϵ_n^h\stackrel{~}{ϵ}_n^h)`$ (B2)
$`\ddot{\stackrel{~}{ϵ}}{}_{n}{}^{h}+\alpha \dot{\stackrel{~}{ϵ}}{}_{n}{}^{h}+\stackrel{~}{A}_n^h(t)\stackrel{~}{ϵ}_n^h`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(ϵ_n^v+\delta _n^h\stackrel{~}{ϵ}_n^h),`$ (B3)
where $`A_n(t)`$ are time-periodic coefficients determined by the given breather state.

The substitution

$`\left(\begin{array}{c}ϵ_n^v\\ ϵ_n^h\\ \stackrel{~}{ϵ}_n^h\end{array}\right)=e^{\frac{1}{2}\alpha t}\left(\begin{array}{c}\kappa _n^v\\ \kappa _n^h\\ \stackrel{~}{\kappa }_n^h\end{array}\right)`$ (B10)
allows one to eliminate the dissipation

$`\ddot{\kappa }_n^v+B_n^v(t)\kappa _n^v`$ $`=`$ $`{\displaystyle \frac{1}{\beta _L}}(\mathrm{}\kappa _n^v+\kappa _{n1}^h\stackrel{~}{\kappa }_{n1}^h)`$ (B11)
$`\ddot{\kappa }_n^h+B_n^h(t)\kappa _n^h`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(\kappa _n^v+\kappa _n^h\stackrel{~}{\kappa }_n^h)`$ (B12)
$`\ddot{\stackrel{~}{\kappa }}{}_{n}{}^{h}+\stackrel{~}{B}_n^h(t)\stackrel{~}{\kappa }_n^h`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(\kappa _n^v+\kappa _n^h\stackrel{~}{\kappa }_n^h),`$ (B13)
where $`B_n(t)=\frac{1}{4}\alpha ^2+A_n(t)`$. Introducing new variables

$`z_n^v=\kappa _n^v,z_n^h=\sqrt{\eta }\kappa _n^h,\stackrel{~}{z}_n^h=\sqrt{\eta }\stackrel{~}{\kappa }_n^h,`$ (B15)
we find the system of equations

$`\ddot{z}_n^v+B_n^v(t)z_n^v`$ $`=`$ $`{\displaystyle \frac{1}{\beta _L}}\mathrm{}z_n^v+{\displaystyle \frac{1}{\beta _L\sqrt{\eta }}}(z_{n1}^h\stackrel{~}{z}_{n1}^h)`$ (B16)
$`\ddot{z}_n^h+B_n^h(t)z_n^h`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(z_n^h\stackrel{~}{z}_n^h){\displaystyle \frac{1}{\beta _L\sqrt{\eta }}}z_n^v`$ (B17)
$`\ddot{\stackrel{~}{z}}_n^h+\stackrel{~}{B}_n^h(t)\stackrel{~}{z}_n^h`$ $`=`$ $`{\displaystyle \frac{1}{\eta \beta _L}}(z_n^h\stackrel{~}{z}_n^h)+{\displaystyle \frac{1}{\beta _L\sqrt{\eta }}}z_n^v.`$ (B18)
These equations describe a Hamiltonian system, namely

$`\dot{\stackrel{}{z}}_n`$ $`=`$ $`{\displaystyle \frac{}{\stackrel{}{p}_n}}`$ (B19)
$`\dot{\stackrel{}{p}}_n`$ $`=`$ $`{\displaystyle \frac{}{\stackrel{}{z}_n}},`$ (B20)
where $`\stackrel{}{z}_n=(z_n^v,z_n^h,\stackrel{~}{z}_n^h),\stackrel{}{p}_n=(p_n^v,p_n^h,\stackrel{~}{p}_n^h)`$ and the Hamiltonian $`(\stackrel{}{z}_n,\stackrel{}{p}_n,t)`$ is

$``$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{n}{}}[p_n^{v^2}+p_n^{h^2}+\stackrel{~}{p}_n^{h^2}]`$ (B25)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{n}{}}[B_n^vz_n^{v^2}+B_n^hz_n^{h^2}+\stackrel{~}{B}_n^h\stackrel{~}{z}_n^{h^2}]`$
$`+{\displaystyle \frac{1}{2\beta _L}}{\displaystyle \underset{n}{}}(z_n^vz_{n1}^v)^2`$
$`+{\displaystyle \frac{1}{2\beta _L\eta }}{\displaystyle \underset{n}{}}(z_n^h\stackrel{~}{z}_n^h)^2`$
$`+{\displaystyle \frac{1}{\beta \sqrt{\eta }}}{\displaystyle \underset{n}{}}z_n^v[z_{n1}^h\stackrel{~}{z}_{n1}^hz_n^h+\stackrel{~}{z}_n^h].`$
Since the particular Hamiltonian can be represented in a general quadratic form, the symplectic product of two different trajectories $`\{\stackrel{}{p}_n(t),\stackrel{}{z}_n(t)\}`$ and $`\{\stackrel{}{p}_{n}^{}{}_{}{}^{}(t),\stackrel{}{z}_n^{}(t)\}`$ does not change in time

$`={\displaystyle \underset{n}{}}[\stackrel{}{p}_n^{}(t)\stackrel{}{z}_n(t)\stackrel{}{p}_n(t)\stackrel{}{z}_n^{}(t)].`$ (B26)
Rewriting our set of equations in the form

$`\delta \dot{\stackrel{}{z}}_n`$ $`=`$ $`{\displaystyle \frac{^2}{\stackrel{}{p}_n^2}}\delta \stackrel{}{p}_n+{\displaystyle \frac{^2}{\stackrel{}{p}_n\stackrel{}{z}_n}}\delta \stackrel{}{z}_n`$ (B27)
$`\delta \dot{\stackrel{}{p}_n}`$ $`=`$ $`{\displaystyle \frac{^2}{\stackrel{}{z}_n\stackrel{}{p}_n}}\delta \stackrel{}{p}_n+{\displaystyle \frac{^2}{\stackrel{}{z}_n^2}}\delta \stackrel{}{z}_n,`$ (B28)
and using the notation

$`𝒥=\left(\begin{array}{cc}0& E\\ E& 0\end{array}\right),`$ (B31)
where $`E`$ is the identity matrix, we obtain

$`\left(\begin{array}{c}\delta \dot{\stackrel{}{p}}_n\\ \delta \dot{\stackrel{}{z}_n}\end{array}\right)=𝒥^1^2\left(\begin{array}{c}\delta \stackrel{}{p}_n\\ \delta \stackrel{}{z}_n\end{array}\right),`$ (B36)
where $`^2`$ is the Hessian of $``$.

Let us consider the following map by integrating the equations (B36) over one period $`T_b`$ of the initial solution

$`\left(\begin{array}{c}\delta \stackrel{}{p}_n(T_b)\\ \delta \stackrel{}{z}_n(T_b)\end{array}\right)=U(T_b)\left(\begin{array}{c}\delta \stackrel{}{p}_n(0)\\ \delta \stackrel{}{z}_n(0)\end{array}\right).`$ (B41)
Since the form $``$ is symplectic (B26) , $`U(T_b)`$ is symplectic too. As a result, we find that the eigenvalues of $`U(T_b)`$ have to fulfill the condition that if $`\nu `$ is an eigenvalue then $`\frac{1}{\nu }`$, $`\nu ^{}`$ and $`\frac{1}{\nu ^{}}`$ are also eigenvalues. Note that for a marginally stable periodic motion of a Hamiltonian system the Floquet eigenvalues $`\nu `$ are located on the unit circle. Switching to an unstable state is realized by collisions of eigenvalues on the unit circle and departing from it.

All of the obtained relations for eigenvalues (and eigenvectors) can be rewritten for the original $`\delta _n(t)`$ variables (call the corresponding Floquet eigenvalues $`\mu `$). First of all, we expect that most of the eigenvalues will be located on a circle of radius (57). Further it follows that if $`\mu `$ is an eigenvalue then $`\frac{e^{\alpha T_b}}{\mu }`$, $`\mu ^{}`$ and $`\frac{e^{\alpha T_b}}{\mu ^{}}`$ are eigenvalues too. So, all scenarios of collisions of the eigenvalues are similar to the Hamiltonian case. For a periodic motion to be stable we need all eigenvalues $`\mu `$ to have absolute values less or equal to one, i.e. the complex numbers $`\mu `$ should reside inside the unit circle in the complex plane. Since the radius $`R(\alpha )1`$ it needs further finite variation of the control parameters (after a collision) to enforce an eigenvalue to traverse the finite distance to the unit circle and to exit it. Thus not every collision will lead to an instability. However, the decrease of the dissipation parameter $`\alpha `$ may tune the system closer to the Hamiltonian case, making its states more sensitive to any occurring collisions. Note that due to the periodicity of the breather state there is always one eigenvalue $`\mu =1`$, whose eigenvector is tangent to the breather orbit. Consequently we always find another eigenvalue $`\mu =e^{\alpha T_b}`$ which is located on the positive real axis inside the unit circle. During all numerical computations of Floquet eigenvalues and eigenvectors the above properties were tested and complete agreement was found.