text:
# Demonstration of the stability or instability of multibreathers at low coupling

## 1 Introduction

Discrete breathers are time–periodic, localized oscillations in discrete systems due to a combination of nonlinearity and discreteness. They have become a well understood phenomenon since the publication of the proof of existence in Ref. . This proof is the origin of exact and powerful numerical methods to calculate them and to determine their stability . A deeper insight has been achieved since the introduction of Aubry’s band theory . The last reference, together with Ref.  can be considered as reviews, although a new one is badly needed after the huge development of the subject in the last years.

The stability of one–site breathers is proofed in Refs.  under rather general conditions, the stability of the possible two–sites breathers in analyzed in an unpublished theorem in Ref. , page 69. Other methods use the hypothesis of small amplitude oscillations or the rotating wave approximation  leading to the Nonlinear Schrödinger equation (DNLS), which is the actual equation analyzed, and, there are, of course, the efficient but slow numerical methods mentioned above.

Here we propose a method based in the properties of the bands in Aubry’s bands theory and perturbation theory, to obtain the stability properties of any multibreather at low coupling. In many cases it involves only some modification of the coupling matrix and the knowledge of the hardness/softness of the on–site potential. We will briefly summarize the band theory in the next section, while we introduce the notation and some basic concepts, on which our theory is based. However, we refer to Ref.  for detailed explanations, as it will be long and repetitive to expose it in detail.

In section 3 we develop the method for symmetric on–site potentials or in–phase multibreathers, which is synthesized in a theorem, commented in section 4. Its scope is enlarged to non–symmetric potentials in section 3, and to generalized Klein–Gordon systems in section 6. The method is applied to several interesting examples in sections 7 and 8. In Appendix A we calculate the value of a magnitude $`\gamma `$ used in our theory, and in Appendix B we relate the curvature of the bands with the characteristics of the on site potential.

## 2 Band theory and notation

### 2.1 The Newton operator

We consider Klein–Gordon systems with linear coupling described by dynamical equations of the form:

$$\ddot{u}_n+V^{}(u_n)+\epsilon \underset{m=1}{\overset{N}{}}C_{nm}u_n=\mathrm{\hspace{0.17em}0}n=1,\mathrm{},N$$
(1)
where the variables $`u_n`$ are functions of time $`t`$, $`V(u_n)`$ is an homogeneous on–site potential, $`V^{}`$ its derivative, $`\dot{u}_n`$ denotes derivation with respect to time, $`N`$ is the number of oscillators, $`C`$ is a coupling constant matrix, which can describe nearest neighbor or long–range interaction, and includes the boundary conditions, and $`\epsilon `$ is the coupling parameter. We use a notation similar to quantum mechanics, i.e., $`|u[u_1(t),\mathrm{},u_N(t)]^{}`$ ( meaning the transpose matrix). Defining $`V(u)=[V(u_1),\mathrm{},V(u_N)]^{}`$ and analogously its derivatives, equation (1) can be written as:

$$|\ddot{u}+|V^{}(u)+\epsilon C|u=\mathrm{\hspace{0.17em}0}.$$
(2)
We suppose that the functions $`u_n(t)`$ are time–periodic (with period $`T`$ and frequency $`\omega _\mathrm{b}`$), time–reversible solutions and, therefore, they can be written as cosine Fourier series with real coefficients (some of these assumptions will be relaxed later). The (linear) stability properties of a given solution $`|u`$ depend on the properties of the characteristic equation for the Newton operator $`𝒩_\epsilon `$ given by

$$𝒩_\epsilon (u)|\xi |\ddot{\xi }+V^{\prime \prime }(u)|\xi +\epsilon C|\xi =E|\xi ,$$
(3)
where $``$ product is the list product, i.e., $`f(u)|\xi `$ is the column matrix with elements $`f(u_n(t))\xi _n(t)`$. If $`E=0`$, this equation describes the evolution of small perturbations $`|\xi `$ of $`|u`$.

### 2.2 The Floquet matrix

Any solutions of Eq. (3) can be determined by the column matrix of the initial conditions for positions and momenta $`\mathrm{\Omega }(0)=[\xi _1(0),\mathrm{}\xi _N(0),\pi _1(0),\mathrm{}\pi _N(0)]^{}`$, with $`\pi _n(t)=\dot{\xi }_n(t)`$, and $`\mathrm{\Omega }(t)`$ describes its evolution in the space of coordinates and momenta. A base of solutions is given by the $`2N`$ functions with initial conditions $`\mathrm{\Omega }^\nu (0)`$, $`\nu =1,\mathrm{},2N`$, with $`\mathrm{\Omega }_l^\nu (0)=\delta _{\nu l}`$. We will identify often a given solution $`\xi (t)`$ of the Newton equation with the corresponding matrix of initial conditions $`\mathrm{\Omega }(0)`$. As the Newton operator depends on the $`T`$–periodic solution $`u(t)`$, it is also periodic, and the evolution of the solutions of Eq. (3) can be studied by means of the Floquet operator, which maps $`\xi (t)`$ into $`\xi (t+T)`$. In a finite system this is equivalent to the Floquet matrix $`_E`$ being given by:

$$\left[\begin{array}{c}\{\xi _n(T)\}\\ \{\pi _n(T)\}\end{array}\right]=_E\left[\begin{array}{c}\{\xi _n(0)\}\\ \{\pi _n(0)\}\end{array}\right]$$
(4)
These matrices can be easily constructed numerically, by integrating Eq. (3) $`2N`$ times from $`t=0`$ to $`t=T`$, with initial conditions $`\mathrm{\Omega }^\nu (0)`$. Then, the $`\nu `$ column of $`_E`$ is given by $`\mathrm{\Omega }^\nu (T)`$.

### 2.3 Stability and bifurcations

The eigenvalues $`\{\lambda _l\}_{l=1}^{2N}`$ of $`F_0`$, called the Floquet multipliers, determine the linear stability of the solution $`u(t)`$. If there is any eigenvalue with $`|\lambda _l|>1`$, the corresponding eigenfunction $`\xi ^l(t)`$ grows with time and $`u(t)`$ is unstable; if $`\lambda _l=1`$, $`\xi ^l(t)`$ is $`T`$-periodic; if $`\lambda _l=1`$, $`\xi ^l(t)`$ is $`2T`$–periodic. However, $`F_0`$ (and any $`F_E`$) is a symplectic matrix, as it is derived from the symplectic system Eq. (3), which implies that if $`\lambda _l`$ is a non–zero eigenvalue, so it is $`1/\lambda `$. $`F_E`$ is also real, which in turn implies that $`\lambda _l^{}`$ and $`1/\lambda ^{}`$ are also multipliers. Therefore, the eigenvalues of $`F_E`$ come in groups of four, if they are not real, or in pairs, if they are real or if $`|\lambda _l|=1`$. We can write the Floquet multipliers as $`\{\mathrm{exp}(\mathrm{i}\theta _l)\}`$, with $`\theta _l`$, in general, complex numbers. The complex number $`\mathrm{i}\theta _l`$ are called the Floquet exponents, and $`\theta _l`$ the Floquet arguments.

The only possibility for the stability of $`u(t)`$ is that all the eigenvalues of $`F_0`$ have moduli $`1`$, i.e., they are at the unit circle. Therefore, the condition for linear stability of $`u(t)`$ can be described as all the Floquet arguments for $`F_0`$ being real. Moreover, if $`u(t)`$ is stable and, then, every $`|\lambda _l|=1`$, the Floquet multipliers come in complex conjugate pairs with arguments $`(\pm \theta _l)`$ or as double $`1`$ or $`1`$ ($`\theta =0`$ or $`\theta =\pm \pi `$). If a parameter like the coupling $`\epsilon `$ is changed the multipliers of $`F_0`$ change continuously. Therefore, a bifurcation to an unstable solution can only take place in three different forms: a) Harmonic instability: two complex eigenvalues moving along the unit circle collide at $`\lambda =1`$ ($`\theta =0`$); b) Subharmonic instability: two complex eigenvalues moving along the unit circle collide at $`\lambda =1`$ ($`\theta =\pm \pi `$); c) Oscillatory or Hopf instability: two pairs of complex eigenvalues moving along the unit circle collide at $`\pm \theta 0`$ and abandon the unit circle as a quadruplet $`(\lambda ,1/\lambda ,\lambda ^{},1/\lambda ^{})`$.

A Floquet multiplier of $`F_0`$ is always known. Calculating the derivative of Eq. (2) with respect to time we obtain that $`𝒩\dot{u}=0`$. $`\dot{u}`$ is also $`T`$–periodic, therefore, it is an eigenfunction of $`F_0`$ with multiplier $`1`$, and as they come in pairs, there is always a double $`1`$ multiplier. $`\dot{u}`$ is called the phase mode because its meaning is that if $`\dot{u}(t)`$ is a solution of Eq. (2), $`\dot{u}(t+dt)`$ is also a solution. While the solution $`u(t)`$ exists, this double eigenvalue is always there. Due to the possible forms of the bifurcations, if all the multipliers are isolated except the double $`1`$, the system is structurally stable, i.e., there is a neighborhood in the space of the parameters where $`u(t)`$ is stable.

### 2.4 Aubry’s band theory

However, it turns out that much information can be obtained by studying also the Floquet arguments for $`E0`$, with is known as Aubry’s band theory . The set of points $`(\theta ,E)`$, with $`\theta `$ a real Floquet argument of $`F_E`$ have a band structure. The fact that the Floquet multipliers come in pairs of complex conjugate pairs brings about that if $`(\theta ,E)`$ belongs to a band, $`(\theta ,E)`$ does it too, i.e., the bands are symmetric with respect to $`\theta `$ and $`E/\theta (0)=0`$. As $`\dot{u}`$ has eigenvalue $`1`$ ($`\theta =0`$), there is always a band tangent to the axis $`E=0`$ at $`\theta =0`$. There are at most $`2N`$ points for a given value of $`E`$ and, therefore, there are at most $`2N`$ bands crossing any horizontal axes in the space of coordinates $`(\theta ,E)`$. The condition for linear stability of $`u(t)`$ is equivalent to the existence of $`2N`$ bands crossing the axis $`E=0`$ (including tangent points with their multiplicity). If a parameter like the coupling $`\epsilon `$ changes, the bands evolve continuously, and they can loss crossing points with $`E=0`$, bringing about the instability.

We know more about the eigenfunctions of the Newton operator. The fact that it is periodic, or, in other words, that it commutes with the operator $`𝒫`$ of translation in time a period, $`𝒫f(t)=f(t+T)`$, implies that its eigenfunctions can be chosen simultaneously as eigenfunctions of $`𝒫`$, which by Bloch theorem, are given by $`\xi (\theta ,t)=\chi (\theta ,t)\mathrm{exp}(\mathrm{i}\theta t/T)`$, $`\chi (\theta ,t)`$ being a column matrix of $`T`$–periodic functions. It is straightforward to check that $`\mathrm{exp}(\mathrm{i}\theta )`$ is the corresponding Floquet multiplier, and $`\theta `$ its Floquet argument ($`F_E`$ is simply the representation of $`𝒫`$ in the base $`\mathrm{\Omega }^l`$).

### 2.5 Bands at the anticontinuous limit

The key concept on the demonstration of breather existence , single breather stability  and the present paper is the anticontinuous limit, i.e., the system with all the oscillators uncoupled, $`\epsilon =0`$, in Eq. (12). At the anticontinuous limit, Eq. (1) reduces to $`N`$ identical equations:

$$\ddot{u}_n+V^{}(u_n)=\mathrm{\hspace{0.17em}0}.$$
(5)
Supposing that we consider time–reversible solutions around a single minimum of $`V`$, there are only three different solutions: a) oscillators at rest $`u_n=\mathrm{\hspace{0.17em}0}`$; b) excited oscillators with identical $`u_n(t)`$, hereafter denoted $`u^0(t)`$; c) excited oscillators with a phase difference of $`\pi `$ with the previous ones, given by $`u_n(t)=u^0(t+T/2)`$. Each site index is given a code $`\sigma _n`$, which takes elements in $`\{0,1,1\}`$, where $`\sigma _n=\mathrm{\hspace{0.17em}0}`$ represents an oscillator at rest, $`u_n=\mathrm{\hspace{0.17em}0}`$; $`\sigma _n=\mathrm{\hspace{0.17em}1}`$, an oscillator with solution $`u^0(t)`$; and $`\sigma _n=1`$, the solution $`u^0(t+T/2)`$. The matrix of codes $`\sigma =[\sigma _1,\mathrm{},\sigma _N]^{}`$ represents the state of the system at the anticontinuous limit.

We suppose that there are $`p`$ oscillators at rest and $`Np`$ excited oscillators. Equation (3) at $`\epsilon =\mathrm{\hspace{0.17em}0}`$ becomes

$$𝒩_0(u)|\xi |\ddot{\xi }+V^{\prime \prime }(u)|\xi =E|\xi ,$$
(6)
or, equivalently, $`Np`$ identical equations:

$$\ddot{\xi }_n+V^{\prime \prime }(u_n)\xi _n=E\xi _n,$$
(7)
with only a periodic eigenfunction for $`E=\mathrm{\hspace{0.17em}0}`$, the phase mode $`\dot{u}_n(t)`$ for the isolated, excited oscillators. The other one is the growth mode, (see Appendix B) which is not bounded and will not be used in this work. They give rise to $`Np`$ bands tangent to the axis $`E=\mathrm{\hspace{0.17em}0}`$, shown in Fig. 1. They have positive curvature at $`(\theta ,E)=(0,0)`$ if the on–site potential is soft ($`H/\omega _b<0`$) and negative if it is hard ($`H/\omega _b>0`$), as it is demonstrated in Appendix B. The possibility of $`H/\omega _b=0`$ is excluded by the conditions of the breather existence theorem , which means that the on-site potential is truly nonlinear for the isolated oscillators with the breather frequency $`\omega _\mathrm{b}`$. If, as we are supposing here, the on–site potentials $`V(u_n)`$ are identical, the solutions $`u_n(t)`$ are also identical except for a change of phase and the bands are superposed, if not, we might have different bands but with the same general shape.

The remaining $`p`$ equations corresponding to the oscillators at rest are of the form:

$$\ddot{\xi }_n+(\omega _0)^2\xi _n=E\xi _n,$$
(8)
with $`\omega _0=\sqrt{V^{\prime \prime }(0)}`$. They have only a $`T`$-periodic solution for $`E=\mathrm{\hspace{0.17em}0}`$ , the null solution, due to the non–resonance condition for breather existence $`n\omega _b\omega _0`$ ($`p`$, i.e., none of the harmonics of the breather resonates with the oscillators at rest. They provide $`p`$ identical bands which are easily calculated solving the equation above. Its solutions are $`\xi =\mathrm{exp}(\pm \mathrm{i}\sqrt{\omega _0^2E}t)`$, with Floquet multipliers $`\mathrm{exp}(\pm \mathrm{i}\sqrt{\omega _0^2E}T)`$ and Floquet arguments $`\theta =\pm \sqrt{\omega _0^2E}T`$. That is, the bands are given by $`E=\omega _0^2\omega _\mathrm{b}^2(\theta /T)^2`$, where $`\theta `$ can be reduced to the first Brillouin zone $`[\pi ,\pi ]`$ by the addition of $`2\pi p`$, $`p`$. If the oscillators have different rest frequencies, the rest bands are not superposed but they have the same characteristics.

The rest bands are also shown in Fig. 1. Note that this figure is only a sketch, for clarity, as very often some excited oscillators bands are very flat and difficult to appreciate at the same scale. This sketch reproduces, however, the basics facts of the band structure.

When the coupling is switched on, the degeneracy of the bands is generically raised. All the tangent bands at $`(\theta ,E)=(0,0)`$ except one, than continues there, move upwards or downwards. If the on–site potential is soft, Fig. 1 (left), and a band moves upwards, a double tangent point with the axis $`E=0`$ is lost, or, in other terms, a pair of Floquet arguments of $`F_0`$ becomes complex and the solution $`u(t)`$ is unstable. If the on site potential is hard, Fig. 1 (right), the same occurs when a band moves downwards.

## 3 Multibreathers stability

In this work use degenerate perturbation theory  to demonstrate the stability or instability of the breathers of any code. Degenerate perturbation theory establishes that if $`𝒩_0`$ is a linear operator with a degenerate eigenvalue $`E_0`$, with eigenvectors $`\{|v_n\}`$, which are ortonormal with respect to a scalar product, i.e., $`v_n|v_m=\delta _{nm}`$, and if $`\epsilon \stackrel{~}{𝒩}`$ is a perturbation of $`𝒩_0`$, with $`\epsilon `$ small; then, to first order in $`\epsilon `$, the eigenvalues of $`𝒩_0+\epsilon \stackrel{~}{𝒩}`$ are $`E_0+\epsilon \lambda _i`$, with $`\lambda _i`$ being the eigenvalues of the perturbation matrix $`Q`$ with elements $`Q_{nm}=v_n|\stackrel{~}{𝒩}|v_m`$. Note that perturbation theory as described in the reference cited, is time independent perturbation theory, but our time variable $`t`$ is their spatial coordinate $`x`$.

Let us consider again Eq. (3) with zero coupling, $`\epsilon =0`$. As explained in the previous section, if there are $`Np`$ excited oscillators, there $`Np`$ zero eigenvalues corresponding to the $`T`$–periodic eigenfunction $`\dot{u}_n`$ (the phase modes of the isolated oscillators), or in other words, a $`Np`$ times degenerate eigenvalue $`E_0=0`$ if we restrict the domain of $`𝒩`$ to periodic functions. What we need to know is the sign of this degenerate eigenvalues when the coupling is switched on.

To apply degenerate theory to multibreather stability we need to identify the perturbation operator, a suitable scalar product and a ortonormal basis of the eigenspace with eigenvalue $`E_0=0`$. The scalar product is defined as:

$$\xi _1|\xi _2=\underset{n=1}{\overset{N}{}}_{T/2}^{T/2}\xi _1^{}(t)\xi _2(t).t.$$
(9)
Let us suppose initially that all the excited oscillator are identical and vibrate in phase. We will denote by $`u^0(t)`$ these identical solutions. In this case, all the phase modes are also identical and will be denoted as $`\dot{u}^0`$. The $`Np`$ elements of the basis are

$$|n=\frac{1}{\mu }[0,\mathrm{},0,\dot{u}^0,0,\mathrm{},0]^{},$$
(10)
with the non-zero element at the position $`n`$, $`n`$ being the index of the excited oscillators, and $`\mu =\sqrt{_{T/2}^{T/2}(\dot{u}^0)^2t}`$. It is straightforward to check that they are ortonormal.

It is enough to know the eigenvalues corresponding to periodic and real solutions, i.e., with Floquet argument $`\theta =\mathrm{\hspace{0.17em}0}`$, as the intersections of the bands with the axis $`\theta =\mathrm{\hspace{0.17em}0}`$ correspond to periodic solutions of (3). Therefore the $`Np`$ $`|n`$’s form the basis of the degenerate eigenvalue $`E_0=0`$ needed to apply perturbation theory.

To obtain the perturbation operator, we expand in Taylor series Eq. (3) at $`\epsilon =\mathrm{\hspace{0.17em}0}`$ and obtain to first order in $`\epsilon `$ :

$`𝒩_\epsilon (u)|\xi `$ $`=`$ $`|\ddot{\xi }+V^{\prime \prime }(u)|\xi +\epsilon \left(V^{\prime \prime \prime }(u)u_\epsilon |\xi +C|\xi \right)=`$ (11)
$`=`$ $`(E_0+\epsilon \lambda _i)|\xi ,`$
with $`u_\epsilon =\left(\frac{u}{\epsilon }\right)_{\epsilon =0}`$ and $`u`$ is also the $`\epsilon =0`$ solution.

Therefore the perturbation operator of $`𝒩_0`$ is

$$\stackrel{~}{𝒩}|\xi =V^{\prime \prime \prime }(u))u_\epsilon |\xi +C|\xi ,$$
(12)
where $`u`$ and $`u_\epsilon `$ are calculated at $`\epsilon =0`$. We do not know $`u_\epsilon `$, but deriving with respect to $`\epsilon `$, at $`\epsilon =\mathrm{\hspace{0.17em}0}`$, the dynamical equations (2) we obtain:

$$|\ddot{u}_\epsilon +V^{\prime \prime }(u)|u_\epsilon +C|u=\mathrm{\hspace{0.17em}0}\text{or}N_0|u_\epsilon =C|u.$$
(13)
The perturbed eigenvalues $`E_i`$ are $`\epsilon \lambda _i`$ , $`\lambda _i`$ being the eigenvalues of the perturbation matrix $`Q`$, with elements $`n|\stackrel{~}{𝒩}|m`$. The matrix $`\stackrel{~}{C}`$ of elements $`n|C|m`$ is simply the matrix $`C`$ without the columns and rows corresponding to the oscillators at rest. The other terms are:

$`n|V^{\prime \prime \prime }(u)u_\epsilon |m=`$
$`{\displaystyle \frac{1}{\mu ^2}}{\displaystyle _{T/2}^{T/2}}[\mathrm{},0,\dot{u}_n^0,0,\mathrm{}]`$ $`[\mathrm{},0,V^{\prime \prime \prime }(u^0)u_{m,\epsilon }\dot{u}_m^0,0,\mathrm{}]^{}t=`$ (14)
$`{\displaystyle \frac{\delta _{nm}}{\mu ^2}}{\displaystyle _{T/2}^{T/2}}\dot{u}^0V^{\prime \prime \prime }(u^0)u_{n,\epsilon }\dot{u}^0t,`$
with $`u_{n,\epsilon }=\left(\frac{u_n}{\epsilon }\right)_{\epsilon =0}`$. Thus, only the diagonal elements $`n|V^{\prime \prime \prime }(u)u_\epsilon |n`$ are non–zero. To calculate the last integral in (14) we will integrate by parts and use that the integral in a period of the derivative of a periodic function is zero. Besides, the functions $`u_{n,\epsilon }`$ are periodic as the coefficients of their Fourier series are given by the derivatives with respect to $`\epsilon `$ of the Fourier coefficients of $`u_n`$. In the deduction below, all the integral limits are $`T/2`$ and $`T/2`$, and the terms between brackets from integration by parts will be zero. The last integral in Eq. (14) becomes:

$`\left[\dot{u}^0u_{n,\epsilon }V^{\prime \prime }(u^0)\right]_{T/2}^{T/2}{\displaystyle V^{\prime \prime }(u^0)\dot{u}^0\dot{u}_{n,\epsilon }t}{\displaystyle V^{\prime \prime }(u^0)\ddot{u}^0u_{n,\epsilon }t}`$ $`=`$
$`[V^{}(u^0)\dot{u}_{n,\epsilon },]_{T/2}^{T/2}+{\displaystyle }V^{}(u^0)\ddot{u}_{\epsilon ,n},t{\displaystyle }V^{\prime \prime }(u^0)\ddot{u}^0u_{n,\epsilon }t`$ $`=`$
$`{\displaystyle \ddot{u}^0(\ddot{u}_{n,\epsilon }+V^{\prime \prime }(u^0)u_{n,\epsilon })t}`$ (15)
The term between parentheses, is the $`n`$ component of the lhs of Eq. (13), i.e., it becomes $`_mC_{nm}u_m^0`$, where $`u_m^0=u^0`$, if the oscillator $`m`$ is excited, and zero otherwise, i.e., it is $`_m\stackrel{~}{C}_{nm}u^0=(_m\stackrel{~}{C}_{nm})u^0`$. Equation (LABEL:eq:larga) becomes:

$`({\displaystyle \underset{m}{}}\stackrel{~}{C}_{nm}){\displaystyle \ddot{u}^0u^0𝑑t}`$ $`=`$
$`({\displaystyle \underset{m}{}}\stackrel{~}{C}_{nm})\left(\left[\dot{u}^0u^0\right]_{T/2}^{T/2}{\displaystyle (\dot{u}^0)^2𝑑t}\right)`$ $`=`$ $`({\displaystyle \underset{m}{}}\stackrel{~}{C}_{nm})\mu ^2.`$ (16)
That is, equation (14), leads to:

$$n|V^{\prime \prime \prime }(u)u_\epsilon |n=\underset{m}{}\stackrel{~}{C}_{nm}\text{at}\epsilon =\mathrm{\hspace{0.17em}0}.$$
(17)
Therefore the diagonal elements of the perturbation matrix $`Q`$ are

$$Q_{nn}=\underset{m}{}\stackrel{~}{C}_{nm}+\stackrel{~}{C}_{nn}=\underset{mn}{}\stackrel{~}{C}_{nm}=\underset{mn}{}Q_{nm}.$$
(18)
To summarize, the perturbation matrix $`Q`$ is given by:

$$Q_{nm}=\stackrel{~}{C}_{nm},nm,Q_{nn}=\underset{mn}{}Q_{nm},$$
(19)
$`\stackrel{~}{C}`$ being the coupling matrix without the $`p`$ rows and columns corresponding to oscillators at rest.

This result can be very easily extended, to the case when there are oscillators out of phase, at least for an even potential $`V(u)`$. In this case $`u^0(t+T/2)=u^0(t)`$ and $`\dot{u}^0(t+T/2)=\dot{u}^0(t)`$. Let us define a new code $`\stackrel{~}{\sigma }`$, which is equal to $`\sigma `$, but with $`+1`$ instead of zeros, and perform the ansatz $`u_n=\stackrel{~}{\sigma }_n\stackrel{~}{u}_n`$. By substitution in the dynamical equations (1), and taking into account that $`V^{}(u_n)=V^{}(\stackrel{~}{\sigma }_n\stackrel{~}{u}_n)=\stackrel{~}{\sigma }_nV^{}(\stackrel{~}{u}_n)`$, because $`V^{}(u)`$ is odd, we obtain:

$$\stackrel{~}{\sigma }_n\ddot{\stackrel{~}{u}}_n+\stackrel{~}{\sigma }_nV^{}(\stackrel{~}{u}_n)+\underset{m}{}C_{nm}\stackrel{~}{\sigma }_m\stackrel{~}{u}_m=\mathrm{\hspace{0.17em}0}.$$
(20)
Multiplying by $`\stackrel{~}{\sigma }_n`$, and using that its square is $`+1`$ we get:

$$\ddot{\stackrel{~}{u}}_n+V^{}(\stackrel{~}{u}_n)+\underset{m}{}\stackrel{~}{\sigma }_nC_{nm}\stackrel{~}{\sigma }_m\stackrel{~}{u}_m=\mathrm{\hspace{0.17em}0}.$$
(21)
As the elements of $`C`$ with an index corresponding to an oscillator at rest will not appear in the calculations, we can define a new coupling matrix, with elements $`\sigma _nC_{nm}\sigma _m`$, i.e., $`\text{Diag}(\sigma )C\text{Diag}(\sigma )`$. The eigenvalues $`E_i`$ will have the right sign, but the eigenfunctions $`\stackrel{~}{\xi }`$ of the Newton operator corresponding to (21) have to be transformed to recover the original ones by $`|\xi =\text{Diag}(\sigma )|\stackrel{~}{\xi }`$

The result can be summarized in the following way:

###### Theorem 1 (Symmetric MST)

Consider a Klein–Gordon system, with homogeneous on–site potential $`V`$, linear coupling matrix $`C`$ and coupling parameter $`\epsilon `$, and a multibreather given by a vector code $`\sigma `$ at zero coupling, with $`p`$ zero elements, and $`Np`$ non–zero elements. The potential $`V`$ must be even if there is any $`1`$ in $`\sigma `$. We define the following matrices: $`S=\text{Diag}(\sigma )`$; $`\stackrel{~}{C}`$ the $`Np`$, squared matrix equal to $`SCS`$ but without the rows and columns corresponding to the zero codes in $`\sigma `$; $`Q`$, the perturbation matrix, with non–diagonal elements $`Q_{nm}=\stackrel{~}{C}_{nm}`$, and diagonal elements $`Q_{nn}=_{mn}\stackrel{~}{C}_{nm}`$. Let $`\{\lambda _i\}_{i=1}^{Np}`$ be the eigenvalues of $`Q`$ and suppose that there is only one zero among them. Then:
a) if $`V(u_n)`$ is hard and there is any negative value in $`\{\epsilon \lambda _i\}`$ the multibreather at low coupling will be unstable, and stable otherwise.
b) if $`V(u_n)`$ is soft and there is any positive value in $`\{\epsilon \lambda _i\}`$ the multibreather at low coupling will be unstable, and stable otherwise.

There is always a zero eigenvalue is explained below. If there are more than one the stability is undefined within first order perturbation theory.

## 4 Comments

### 4.1 Global phase mode

The perturbation matrix $`Q`$ has always a zero eigenvalue corresponding to the global phase mode. It can be useful to make it apparent as it reduces the order of the secular equation by an unity. Considering the matrix $`Z=Q\lambda I`$, $`det(Z)=0`$ is the secular equation. Its roots are the same if we change $`Z`$ to $`Z^{}`$ obtained by summing to the first column the rest of them, i.e., $`Z_{n\mathrm{\hspace{0.17em}1}}^{}=_{mn}\stackrel{~}{C}_{n,m}\lambda +_{mn}\stackrel{~}{C}_{n,m}=\lambda `$. Now, subtracting the first row from the rest of them, the first column has only a non zero element, which is equal to $`\lambda `$. Therefore, $`\lambda =0`$ is an eigenvalue of $`Q`$ and the secular equation is reduced by an unity.

The corresponding eigenvector $`\stackrel{~}{x}`$ satisfies $`Q\stackrel{~}{x}=\mathrm{\hspace{0.17em}0}`$, i.e., $`(_{mn}\stackrel{~}{C}_{nm})\stackrel{~}{x}_n+_{mn}\stackrel{~}{C}_{nm}\stackrel{~}{x}_m=\mathrm{\hspace{0.17em}0}`$, $`n`$, with solution $`\stackrel{~}{x}_n=\mathrm{\hspace{0.17em}1}`$, $`n`$. Therefore, $`x=S\stackrel{~}{x}=\sigma `$, which corresponds to the global phase mode.

### 4.2 Extended and reduced forms

Note that the code $`\sigma `$ is the matrix of components of the multibreather at zero coupling in a N–dimensional basis, $`\{|n\}`$, with the only non zero element $`u^0`$ at the position $`n`$, including the indexes of all the oscillators. For numerical calculation and interpretation of the eigenvalues and eigenvectors of $`Q`$ it can be more convenient to define $`\stackrel{~}{C}=SCS`$, i.e, conserving the rows and columns corresponding to the rest oscillators but all their elements set to zero. It will be denoted the extended form while the previous one is the reduced form. The price for using the extended form is to obtain $`p`$ extra zero eigenvalues, which causes no harm, except if there are more than $`p+1`$, as it can make more difficult the interpretation of the eigenvectors.

### 4.3 Diagonal elements

The diagonal elements of $`C`$ play no role, because they do not appear in the coupling matrix. Therefore a coupling of the spring–type as $`\epsilon (\mathrm{\hspace{0.17em}2}u_nu_{n1}u_{n+1})`$ gives identical results as the dipole-type $`\epsilon (u_{n1}+u_{n+1})`$.

### 4.4 2D and 3D lattices

Note that the theorem applies equally to 2D and 3D lattices, provided that we can number the sites with an index $`n`$, or consider $`n`$ as a multi–index, and construct the coupling matrix $`C`$.

## 5 Non–symmetric potentials

The extension to non–symmetric potentials $`V(u_n)`$ is straightforward, but it involves some numerical or analytical calculation. That is the reason for not including it in theorem 1, and we present it here as another theorem instead. If the potential is non–symmetric, there are two time–symmetric, $`T`$-periodic solutions to the equations for the isolated oscillators at zero coupling, Eq. (5): $`u^0(t)`$, with code $`+1`$, and $`u^0(t+T/2)`$, which is no longer $`u^0(t)`$, with code $`1`$. The $`Np`$ elements of the basis, $`|n`$ will be as in Eq. (10) but with $`\dot{u}^0(t+T/2)`$ instead of $`\dot{u}^0(t)`$ if $`\sigma _n=1`$. Let us define the symmetry coefficient, $`\gamma `$, in the following way:

$$\gamma =\frac{_{T/2}^{T/2}\dot{u}^0(t)\dot{u}^0(t+T/2)t}{_{T/2}^{T/2}\dot{u}^0(t)\dot{u}^0(t)t}$$
(22)
If $`\{z_k\}_{k=0}^{k_m}`$ are the Fourier coefficients of a $`k_m`$–truncated Fourier series of $`u^0(t)`$, i.e., $`u^0(t)=z_0+_{k=1}^{k_m}2z_k\mathrm{cos}(k\omega _bt)`$, then:

$$\gamma =\frac{_{k=1}^{k_m}(1)^kk^2z_k^2}{_{k=1}^{k_m}k^2z_k^2}.$$
(23)
Note that $`\gamma `$ depends on the breather frequency through the $`\{z_k\}`$, and it is positive and smaller than $`1`$, being its value $`\gamma =\mathrm{\hspace{0.17em}1}`$ for the symmetric potential case. Following the deduction of theorem 1 it is straightforward to obtain:

###### Theorem 2 (Non–symmetric MST)

The generalization of theorem 1 for non–symmetric potentials consists of constructing the perturbation matrix $`Q`$ in the following way:
a) If $`n`$ or $`m`$ are indexes corresponding to oscillators at rest $`Q_{nm}=\mathrm{\hspace{0.17em}0}`$.
b) If $`n`$ and $`m`$, $`nm`$, are indexes corresponding to excited oscillators in phase, i.e., $`\sigma _n\sigma _m=\mathrm{\hspace{0.17em}1}`$, $`Q_{nm}=C_{nm}`$.
c) If $`n`$ and $`m`$, $`nm`$, are indexes corresponding to excited oscillators out of phase , i.e., $`\sigma _n\sigma _m=\mathrm{\hspace{0.17em}1}`$, $`Q_{nm}=\gamma C_{nm}`$.
d) The diagonal elements are $`Q_{nn}=_{mn}Q_{nm}`$
e) The rest of the procedure is exactly the same as described in theorem 1, getting rid of the rest columns and rows (reduced form) or not (extended form), calculating the eigenvalues $`\{\lambda _i\}_1^{Np}`$ of $`Q`$ and obtaining the stability properties of the multibreather from the signs of the $`\{\epsilon \lambda _i\}`$ and the hardness/softness of the on–site potential.

The coefficient $`\gamma `$ depends on the on–site potential and the frequency $`\omega _\mathrm{b}`$. We demonstrate in the Appendix A that $`\gamma =\omega _\mathrm{b}`$ for the Morse potential. Figure 2 shows some numerically calculated values for different potentials. They tend to $`+1`$ when $`\omega _\mathrm{b}\omega _0=1`$ because at this limit $`\dot{u}^0(t)=A\mathrm{sin}(\omega _\mathrm{b}t)`$ and $`\dot{u}^0(t+T/2)=\dot{u}^0(t)`$ as in the symmetric on–site potential case. Generally speaking, if the frequencies are not very far away from the linear frequency $`\omega _0`$, $`\gamma `$ is usually close enough to the unity and the sign of the eigenvalues do not change with respect to the symmetric on–site potential case $`\gamma =1`$.

## 6 Generalization

Theorems 1 and 2 are equally applicable if the system has inhomogeneities at the coupling matrix but not if they are at the on–site potential. Here we introduce the generalization to a broad class of Klein–Gordon systems, including different masses or inertia momenta, inhomogeneous on-site potentials with several minima, nonlinear coupling, solutions with different phases (i.e, not restricted to time–reversible solutions). Let us consider a Klein–Gordon system with Hamiltonian

$$H=\underset{n}{}\left(\frac{1}{2}m_n\dot{u}_n^2+V_n(u_n)\right)+\epsilon W(u),$$
(24)
and dynamical equations

$$m_n\ddot{u}_n+V_n^{}(u_n)+\epsilon \frac{W(u)}{u_n}=\mathrm{\hspace{0.17em}0},$$
(25)
which at zero coupling become:

$$m_n\ddot{u}_n+V_n^{}(u_n)=\mathrm{\hspace{0.17em}0}n.$$
(26)
Suppose we can label the solutions of these equations, with given frequency $`\omega _\mathrm{b}`$ and period $`T`$, with a code composed of an index for the site $`n`$, another index for the different minima and solutions (biperiodic, etc, different potential wells) $`\sigma _n`$, which can be different for each site. If we considering a periodic solution for the $`n`$–th equation , it can be written as a $`2\pi `$ periodic function $`g_{n\omega _\mathrm{b}\sigma _n}(\omega _\mathrm{b}t+\alpha )`$. Let us suppose that $`g_{n\omega _\mathrm{b}\sigma _n}`$ is chosen such that $`\alpha =0`$ corresponds to the time symmetric solution with $`\dot{g}_{n\omega _\mathrm{b}\sigma _n}(0)=0`$ and $`g_{n\omega _\mathrm{b}\sigma _n}(0)>g_{n\omega _\mathrm{b}\sigma _n}(\pi )`$. Therefore, the solutions of equations (26) are determined by the code $`\{n,\sigma _n,\alpha _n\}`$. If we consider a given solution at the anticontinuous limit we can simply write it as $`u=[u_1^0,\mathrm{},u_N^0]^{}`$, but now, the $`u_n^0`$ are, in general, different, because of the particular multibreather solution $`u`$, whose stability we are interested in.

The characteristic equations for the Newton operator $`𝒩_\epsilon `$ are given by:

$$(𝒩_\epsilon \xi )_nm_n\ddot{\xi }_n+V_n^{\prime \prime }(u_n)\xi _n+\epsilon \underset{m}{}\frac{^2W(u)}{u_nu_m}\xi _m=E\xi _n$$
(27)
which to first order in $`\epsilon `$ calculated at $`\epsilon =0`$ lead to

$$(𝒩_0\xi )_n+\epsilon (\stackrel{~}{𝒩}\xi )_n=(E_0+\epsilon \lambda )\xi _n,$$
(28)
where:

$`(𝒩_0\xi )_n`$ $``$ $`m_n\ddot{\xi }_n+V_n^{\prime \prime }(u_n^0)\xi _n=E_0\xi _n,`$ (29)
$`(\stackrel{~}{𝒩}\xi )_n`$ $``$ $`V^{\prime \prime \prime }(u_n^0)u_{n,\epsilon }\xi _n+{\displaystyle \underset{m}{}}{\displaystyle \frac{^2W(u^0)}{u_nu_m}}\xi _m=\lambda \xi _n.`$ (30)
with $`u_{n,\epsilon }=\left(\frac{u_n}{\epsilon }\right)_{\epsilon =0}`$. The first set of equations is the characteristic equation for the Newton operator at zero coupling, and the second one for the perturbation operator $`\stackrel{~}{𝒩}`$. The only periodic solution of each equation (29) with $`E_0=0`$ for a given solution $`u_0`$ of (26) , is $`\dot{u}_n^0`$. Thus, we can construct a ortonormal basis of $`\mathrm{ker}(N_0)`$, restricted to periodic solutions (and, therefore, that can be chosen real) of given period $`T`$, as $`|n=\frac{1}{\mu _n}[0,\mathrm{},0,\dot{u}_n^0,0,\mathrm{},0]^{}`$, with the only nonzero element at the $`n`$ position, and $`\mu _n=\sqrt{_{T/2}^{T/2}(\dot{u}_n^0)^2t}`$.

Deriving Eq. (25) with respect to $`\epsilon `$ at $`\epsilon =0`$ we get:

$$\left(𝒩_0\frac{u}{\epsilon }\right)_nm_n\ddot{u}_{n,\epsilon }+V_n^{\prime \prime }(u_n^0)u_{n,\epsilon }=\underset{m}{}\frac{^2W(u^0)}{u_nu_m}$$
(31)
Following exactly the same procedure as in section 3 we get the perturbation matrix $`Q`$, in extended form, with non diagonal elements for indexes corresponding to excited oscillators:

$$Q_{nm}=n|_{uu}^2W(u^0)|m=\frac{1}{\mu _n\mu _m}_{T/2}^{T/2}\dot{u}_n\frac{^2W(u^0)}{u_nu_m}\dot{u}_mt,nm$$
(32)
and diagonal elements

$$Q_{nn}=\underset{mn}{}\frac{\mu _m}{\mu _n}Q_{nm}.$$
(33)
If there are $`p`$ oscillators at rest and $`Np`$ excited oscillators, the extended matrix $`Q`$ has $`p`$ zero columns and rows corresponding to the indexes of the oscillators at rest. The reduced form of $`Q`$ is the $`Np`$ squared matrix, equal to $`Q`$ but stripped of those zero rows and columns. As a result we can establish the theorem:

###### Theorem 3 (Generalized MST)

Given a generalized Klein–Gordon system (24), a specific multibreather solution at zero coupling $`\{u_n^0\}`$, determined by a suitable set of codes, $`u`$ the corresponding solution at low coupling, $`\{\lambda \}_{i=1}^{Np}`$ the eigenvalues of the reduced matrix $`Q`$, with only one zero, and being the wells corresponding to the $`Np`$ specific $`u_n^0`$, of the same type, hard or soft, then:
The solution $`u`$ is stable if:
a) The on–site potentials are soft and there is not any positive value in $`\{\epsilon \lambda _i\}_{i=1}^{Np}`$.
b)The on–site potentials are hard and there is not any negative value in $`\{\epsilon \lambda _i\}_{i=1}^{Np}`$

This theorem includes both theorems 1 and 2. Certainly, it is of not so straightforward applicability, although for certain cases there are only a few matrix elements $`Q_{nm}`$ to calculate. As an example, if there are just two types of sites, the on–site potentials are symmetric, the coupling is linear and we consider only time–reversible multibreathers, there are only two different matrix elements to calculate.

## 7 Applications

### 7.1 Two–site breathers

A trivial example is the theorem about the stability of the two–site breather found in Ref.  page 69. Here we generalize it to non–symmetric on–site potentials. The coupling matrix corresponding to periodic boundary conditions and of the standard (spring–type) attractive form, has elements $`C_{nm}=\mathrm{\hspace{0.17em}2}\delta _{nm}\delta _{|nm|,1}`$ and $`C_{N\mathrm{\hspace{0.17em}1}}=C_{1N}=1`$, that is:

$$C=\left[\begin{array}{cccccc}\hfill \mathrm{\hspace{0.33em}2}& \hfill 1& \hfill 0& \hfill \mathrm{}& \hfill \mathrm{\hspace{0.33em}0}& \hfill 1\\ \hfill 1& \hfill 2& \hfill 1& \hfill 0& \hfill \mathrm{}& \hfill 0\\ \hfill \mathrm{}& \hfill \mathrm{}& \hfill \mathrm{}& \hfill \mathrm{}& \hfill \mathrm{}& \hfill \mathrm{}\\ \hfill 1& \hfill 0& \hfill \mathrm{}& \hfill 0& \hfill 1& \hfill 2\end{array}\right].$$
(34)
If the two oscillator are in phase, the perturbation matrix is

$$Q=\left[\begin{array}{cc}\hfill \mathrm{\hspace{0.33em}1}& \hfill 1\\ \hfill 1& \hfill 1\end{array}\right],$$

with eigenvalues $`(\lambda _1,\lambda _2)=(0,2)`$, the first one corresponds to the phase mode, the second to the antisymmetric eigenmode. Thus, one of the two bands tangent to the axis $`E=0`$ at zero coupling will remain there, and the other will move upwards. If the on–site potential is soft, i.e., with positive curvature, a intersection/tangent point will be lost, and the two site breather is unstable. For the hard on–site potential, the curvature of the bands is negative, and when the band moved upwards, the intersection points will be kept, and the system is be stable. If $`\epsilon <0`$ the signs of the eigenvalues $`\epsilon E_i`$ are reversed and so are the conclusions.

For the out–of–phase, two-site breather, with code $`\sigma =[1,1]`$, the perturbation matrix is:

$$Q=\left[\begin{array}{cc}\hfill \gamma & \hfill \gamma \\ \hfill \gamma & \hfill \gamma \end{array}\right],$$

with eigenvalues $`(\lambda _1,\lambda _2)=(0,2\gamma )`$. The conclusions above are reversed, being the instability mode, when appropriate, the symmetric one. Note that the boundary conditions in $`C`$ have no effect.

###### Theorem 4 (Two–site breathers)

Consider time–reversible, two–site breathers in an homogeneous Klein-Gordon system with attractive linear coupling, then:
a) The two-site breathers with codes $`\pm [1,1]`$ and soft on–site potential are unstable and with $`\pm [1,1]`$ are stable.
b) The two-site breathers with code $`\pm [1,1]`$ and hard on–site potential are stable and with $`\pm [1,1]`$ are unstable.
c) If the coupling is repulsive, the conclusions are reversed.

### 7.2 Three–site breathers

Let us renumber the sites $`1,2,3`$ and suppose that the code of the site in the middle is $`\sigma _2,=+1`$ as a reference, i.e., $`\sigma =[\sigma _1,1,\sigma _3]`$. We define $`\alpha _i=1`$ if $`\sigma _i=1`$ and $`\alpha _i=\gamma `$ if $`\sigma _i=1`$, for $`i=1,2`$. Using $`C`$ from Eq. (34) the perturbation matrix is:

$$Q=\left[\begin{array}{ccc}\hfill \alpha _1& \hfill \alpha _1& \hfill \mathrm{\hspace{0.33em}0}\\ \hfill \alpha _1& \hfill \alpha _1+\alpha _3& \hfill \alpha _3\\ \hfill 0& \hfill \alpha _3& \hfill \alpha _3\end{array}\right].$$
(35)
Its nonzero eigenvalues are $`\lambda _\pm =\alpha _1+\alpha _3\pm \sqrt{(\alpha _1+\alpha _3)^2\mathrm{\hspace{0.17em}3}\alpha _1\alpha _3}`$. In order to have both eigenvalues with the same sign, we need that $`\alpha _1\alpha _3>0`$, i.e., either $`\alpha _3=\alpha _1=\mathrm{\hspace{0.17em}1}`$, or $`\alpha _3=\alpha _1=\gamma `$. That is, if $`\sigma _1=\sigma _3=1`$, $`\lambda _\pm =1,3>0`$, if $`\sigma _1=s_3=\gamma `$, $`\lambda _\pm =\gamma ,3\gamma <0`$, and if $`\sigma _1\sigma _2=1`$, then $`\lambda _\pm `$ have different signs. Therefore:

###### Theorem 5 (Three–site breathers)

Consider time–reversible, three–site breathers in an homogeneous K-G system with nearest–neighbor, linear coupling, attractive for $`\epsilon >0`$. Then:
a) The three–site breathers with codes $`\sigma =\pm [1,1,1]`$ are stable if $`\epsilon >0`$ ($`\epsilon <0`$)and the on–site potential is soft (hard).
b) The three–site breathers with code $`\sigma =\pm [1,1,1]`$ are stable if $`\epsilon >0`$ $`(\epsilon <0)`$ and the on–site potential is hard (soft).
c) All other three-site breathers are unstable for any sign of $`\epsilon `$ and type of on–site potential.

We can compare the predictions of the values given by the theory with the numerically calculated values of the intersection points of the bands with the axis $`\theta =0`$. As an example, the numerically calculated eigenvalues and the predictions given by theorems 1 and 2 for the three–site breather with code $`[1,1,1]`$ are plotted in Figure 3. The very good accuracy is evident.

### 7.3 A model with long–range interaction

An example of system with long range interaction is the twist model used in Ref. . The dynamical equations are:

$$\ddot{u}_n+V^{}(u_n)+\epsilon \underset{m=nN/2}{\overset{n+N/2}{}}\frac{\mathrm{cos}[\theta _{tw}(nm)]}{|nm|^3}u_m=0$$
(36)
$`V`$ being the Morse potential, and $`\theta _{tw}`$ the angle between two dipole moments corresponding to neighboring base pairs in a simplified model of DNA. The coupling matrix elements are:

$$C_{nm}=\frac{\mathrm{cos}[\theta _{tw}(nm)]}{|nm|^3},$$
(37)
Some of the results, also checked numerically and with good accuracy with respect to the eigenvalues $`E`$, given by theorem 1, in spite of $`V`$ being non–symmetric are:

| Code | $`\theta _{tw}=0`$ | $`\theta _{tw}=\pi `$ |
| --- | --- | --- |
| 11 | Stable | Unstable |
| 1 -1 | Unstable | Stable |
| 101 | Stable | Stable |
| 111 | Stable | Unstable |
| 1-1 1 | Unstable | Stable |
| 11-1 | Unstable | Stable |

As an example of the effect of the symmetry coefficient $`\gamma `$ in theorem 2, we can consider the code $`[1,1,1]`$ and plot the Newton eigenvalues $`E`$ as a function of $`\theta _{tw}`$. The perturbation matrix is:

$$Q=\left[\begin{array}{ccc}\gamma ab& \gamma a& b\\ \gamma a& \mathrm{\hspace{0.33em}\hspace{0.33em}2}\gamma a& \gamma a\\ b& \gamma a& \gamma ab\end{array}\right],$$
(38)
with $`a=\mathrm{cos}(\theta _{tw})`$ and $`b=\mathrm{cos}(2\theta _{tw})/8`$. Figure 4 shows the dependence of the eigenvalues $`\lambda `$ on the twist angle $`\theta _{tw}`$ for three values of the symmetry coefficient $`\gamma `$. It can be seen that only values of $`\gamma `$ very far from $`1`$ lead to a change on the stability prediction.

## 8 Multibreathers, phonobreathers and dark breathers. Parity instabilities.

For larger, time–reversible multibreathers in homogeneous Klein–Gordon systems with nearest–neighbor, linear coupling, it is trivial to obtain the stability properties of the two different solutions with wave number $`q=0`$, all the oscillators in phase, or $`q=\pi `$, neighboring oscillators out of phase. Consider any number $`N^{}`$ of contiguous oscillators in phase in a system coupled by the matrix $`C`$ in Eq. (34). The $`N^{}\times N^{}`$ perturbation matrix is given by:

$$Q=\left[\begin{array}{cccccc}\hfill \mathrm{\hspace{0.33em}1}& \hfill 1& \hfill \mathrm{\hspace{0.33em}0}& \hfill \mathrm{}& \hfill \mathrm{\hspace{0.33em}0}& \hfill \mathrm{\hspace{0.33em}0}\\ \hfill 1& \hfill 2& \hfill 1& \hfill 0& \hfill \mathrm{}& \hfill 0\\ \hfill \mathrm{}& \hfill \mathrm{}& \hfill \mathrm{}& \hfill \mathrm{}& \hfill \mathrm{}& \hfill \mathrm{}\\ \hfill 0& \hfill 0& \hfill \mathrm{}& \hfill 1& \hfill 2& \hfill 1\\ \hfill 0& \hfill 0& \hfill \mathrm{}& \hfill 0& \hfill 1& \hfill 1\end{array}\right].$$
(39)
The characteristic equation for $`Qx=\lambda x`$ corresponds to the following equations:

$`x_1x_2`$ $`=`$ $`\lambda x_1,`$
$`2x_nx_{n+1}x_{n1}`$ $`=`$ $`\lambda x_n,n\mathrm{\hspace{0.17em}1},2,`$
$`x_N^{}x_{N^{}1}`$ $`=`$ $`\lambda x_N^{}.`$ (40)
Substitution of the trial solution $`x_n=\mathrm{cos}(qn\varphi _q)`$, leads to $`\varphi _q=q/2`$, $`q=m\pi /N^{}`$, with $`m=0,\mathrm{},N^{}1`$, and $`\lambda _m=4\mathrm{sin}^2(m\pi /(2N^{}))`$, i.e., positive eigenvalues without a degenerate zero. In case all the oscillators of the group are out of phase, $`Q\gamma Q`$ and $`\lambda _m\gamma \lambda _m`$. Therefore:

###### Theorem 6

If $`\epsilon >0`$, the time–reversible, in-phase, multibreathers are unstable (stable) with soft (hard) on–site potential. For any out–of–phase multibreathers, and $`\epsilon <0`$ the conclusions are reversed.

Although, we have no mathematical proof, according to numerical calculations of the eigenvalues of the perturbation matrices corresponding to groups with different codes, the numbers of negative and positive eigenvalues are equal to the numbers of $`1`$ and $`+1`$ in $`\{\sigma _n\sigma _{n+1}\}_{n=1}^{N^{}1}`$. Therefore being unstable for any on–site potential.

Note that if there are several non–contiguous multibreathers, the perturbation matrix if composed of independent blocks of perturbation submatrices of the same type, each one with a zero eigenvalue by construction as seen in subsection 4.1. The stability of such a system cannot be determined by the MST. The only possible conclusion is that if there are eigenvalues of different sign this system would be unstable.

Phonobreathers are multibreathers without oscillators at rest. For time–reversible phonobreathers in a system either with free–ends, or fixed-ends boundary conditions, the perturbation matrix $`Q`$ is the same as in Eq. (39). If the boundary conditions are periodic, the only change is that $`Q_{\mathrm{1\hspace{0.17em}1}}=Q_{N,N}=2`$, $`Q_{N\mathrm{\hspace{0.17em}1}}=Q_{1,N}=1`$, and the eigenvalues are the same. That is, the mode $`q=0`$ with attractive coupling and hard potential is stable. Changing the hardness, the type of coupling, or the mode, changes each time the stability. In this way, we reproduce the results in Refs.  as a consequence of modulational instability and using the DNLS approximation, respectively, and in theorem 9 in Ref. , based in the properties of the action.

Parity instabilities are a interesting effect for phonobreathers with periodic boundary conditions. It is a consequence of the parity of the total number of sites $`N`$. If, for example, $`V`$ is soft, the coupling is attractive, and $`N`$ is even, all the eigenvalues of the perturbation matrix are negative, i.e., the system is stable. But if $`N`$ is odd, there is an isolated positive eigenvalue, bringing about the instability. There is a clear physical origin of this phenomenon, as the boundaries are equivalent to an inhomogeneity. According to our theory it is also obvious: we have, in fact two oscillators at the ends coupled and in phase, which gives rise to a positive eigenvalue.

Dark breathers are multibreathers with only one or a few oscillators at rest. Their stability depends on the characteristics of the contiguous groups of excited oscillators, as a consequence of theorem 6. We have to be careful about the number of particles, for example the $`\pi `$ mode with code around the dark site $`[\mathrm{},1,1,0,1,1,\mathrm{}]`$, needs an odd number of sites to avoid having two in–phase oscillators at the ends, and a parity instability. The one with $`[\mathrm{},1,1,0,1,1,\mathrm{}]`$ needs an even number of sites. If the system has free of fixed periodic conditions, the perturbation matrix is decoupled in two and the stability is undefined as commented above.

We have compared the predictions of the instabilities and the values of the intersection points of the bands with the axis $`\theta =0`$, in many cases, as for example in Ref. , with very good results. The question of up to which values of the coupling parameter $`\epsilon `$ the theory is valid, is however unsolved. Generically speaking values up to $`\epsilon =0.1`$ (compared to the rest frequency $`\omega _0=1`$), are safe, but there are some exceptions. In the same reference, for the dark breather hard on–site potential and attractive coupling, a mode not predicted by the MST gives rise to an instability as soon as $`\epsilon =0.022`$.

As the width of the phonon band can be calculated, and the width of the background is given by the eigenvalues of the perturbation matrices, it is also possible to predict the values of the $`\epsilon `$, for which oscillatory and subharmonic instabilities occur, but we do not extend further here.

## 9 Summary and conclusions

We have developed using degenerate perturbation theory a method for obtaining the stability properties of multibreathers of any code at low coupling. The method is synthesized in three different versions of a theorem, referred to as the multibreather stability theorem or MST for short. The two simplest versions correspond to time–reversible multibreathers, with linear coupling and homogeneous on–site potentials. If the on–site potential is symmetric or we consider only in–phase multibreathers, it involves only some simple modifications of the coupling matrix. If none of those conditions are fulfilled then it involves the analytical or numerical calculations of a magnitude $`\gamma `$, called the symmetry coefficient, as its value is $`1`$ for symmetric on–site potentials. For soft potentials and values of the frequency until about one half of the rest frequency, a good approximation is the multibreather frequency, result which we demonstrate is exact for the Morse potential. For hard potentials its value its close to unity. The generalized version of the MST is of not so straightforward applicability, although its complexity depends on the characteristics of the system and multibreather considered.

We give some examples of application of the method and compare it with numerical results. The systems considered are two and three–site breathers, with nearest neighbor or long–range interactions, multibreathers, phonobreathers and dark breathers. A parity instability can appear in finite systems depending on the parity of the number of oscillators. All these examples are interesting in themselves, but also illustrate the use of the MST, and show that it provides a powerful method of easy applicability to determine multibreathers stability.

Some other applications under study are multibreathers in 2D and 3D systems, prediction of the oscillatory and subharmonic instabilities and the peculiarities of some disordered systems

This work has been supported by the European Union under the RTN project, LOCNET, HPRN–CT–1999–00163. J Cuevas acknowledges an FPDI grant from ‘La Junta de Andalucía’.

## Appendix A Calculation of the symmetry parameter $`\gamma `$ for a Morse potential

In this appendix, we calculate the value of the symmetry parameter $`\gamma `$ in theorem 2 for the particular case of a Morse potential. The choice of this kind of potential relies in the fact that the expressions of the orbits are easy to manage. We will proof that, for this potential, $`\gamma =\omega _\mathrm{b}`$.

Let us suppose an isolated oscillator submitted to a Morse potential. The energy of the system is:

$$=\frac{1}{2}\dot{x}^2+\frac{1}{2}(\mathrm{exp}(x)1)^2.$$
(41)
We look for solutions with period $`T=2\pi /\omega _\mathrm{b}`$ that are time–reversible, i.e, with $`\dot{x}(0)=0`$, or, in other words, at $`t=0`$ the system is at a turning point. There are two of them, obtained from the equation above: $`x_1=\mathrm{log}(1+\sqrt{2})<0`$ and $`x_2=\mathrm{log}(1\sqrt{2})>0`$. If we consider the solution with $`x(0)=x_1`$, t(x) is given by:

$$t(x)=_{x_1}^x\frac{\mathrm{d}x}{\sqrt{2((\mathrm{exp}(x)1)^2/2)}},$$
(42)
which leads to:

$$t(x)=\frac{1}{\sqrt{12}}\left(\frac{\pi }{2}\mathrm{arcsin}\left[\frac{1}{\sqrt{2}}\left(1+\frac{21}{\mathrm{exp}(x)}\right)\right]\right).$$
(43)
By substitution of $`t=T/2=\pi /\omega _\mathrm{b}`$ and $`x=x_2`$ we obtain $`\omega _\mathrm{b}=\sqrt{12}`$.

The inversion of Eq. (42) leads to:

$$x(t)=\mathrm{log}\frac{1\sqrt{1\omega _\mathrm{b}^2}\mathrm{cos}\omega _\mathrm{b}t}{\omega _\mathrm{b}^2}.$$
(44)
To calculate $`\gamma `$, we need the Fourier coefficients of $`x(t)`$. If the solution $`x(t)`$ is expressed as:

$$x(t)=z_0+2\underset{k=1}{\overset{\mathrm{}}{}}z_k\mathrm{cos}(k\omega _\mathrm{b}t),$$
(45)
the coefficients $`\{z_k\}_{k=0}^{\mathrm{}}`$ are given by:

$$z_k=\frac{1}{T}_{T/2}^{T/2}x(t)\mathrm{cos}(k\omega _\mathrm{b}t)t.$$
(46)
The calculation of these integrals is straightforward and leads to:

$$z_0=\mathrm{log}\frac{1+\omega _\mathrm{b}}{2\omega _\mathrm{b}^2};z_k=\frac{1}{n}\left(\frac{1\omega _\mathrm{b}}{\sqrt{1\omega _\mathrm{b}^2}}\right)^k.$$
(47)
The parameter $`\gamma `$ is given by:

$$\gamma =\frac{_{k1}(1)^kk^2z_k^2}{_{k1}k^2z_k^2}.$$
(48)
By substitution of the Fourier coefficients in Eq. (47), we obtain

$$\gamma =\frac{_{k1}(r)^k}{_{k1}r^k}\text{with}r=\frac{1\omega _\mathrm{b}}{1+\omega _\mathrm{b}}.$$
(49)
As $`|r|<1`$, $`\gamma `$ can be easily calculated using the expressions for the sum of an infinite geometric progression:

$$\underset{k1}{}r^k=\frac{r}{1r},\underset{k1}{}(r)^k=\frac{r}{1+r}$$
(50)
leading to $`\gamma =\omega _\mathrm{b}`$, as we wanted to proof.

## Appendix B Curvature of the bands

In this section we demonstrate a key point used in our theory: the fact that the band corresponding to an isolated excited oscillators has negative curvature if the on site potential is hard and positive if it is soft. Unfortunately, the demonstration is quite long if we want it to be self–consistent. A similar one can be found in Ref.  and different ones in Ref. .

Let us consider the dynamical equation of an isolated oscillator with Hamiltonian $`H=1/2p^2+V(u)`$, with $`u`$ a single–value, real function of time and $`p=\dot{u}`$:

$$\ddot{u}+V^{}(u)=0$$
(51)
The characteristic equation for the Newton operator $`𝒩`$ corresponding to a given periodic solution of Eq. (51)is

$$𝒩\xi \ddot{\xi }+V^{\prime \prime }(u)\xi =E\xi ,$$
(52)
where $`\xi `$ is a $`𝒞^2`$ function of $`t`$. For each eigenvalue $`E`$ there are only 2 independent solutions, determined, for example, by the values of $`\xi (0)`$ and $`\dot{\xi }(0)`$. We can obtain them for $`E=0`$, by deriving Eq. (51) with respect to $`t`$ and with respect to $`\omega _\mathrm{b}`$:

$$\stackrel{\dot{}\dot{}\dot{}}{u}+V^{\prime \prime }(u)\dot{u}=0𝒩\dot{u}=0;\ddot{u}_w+V^{\prime \prime }(u)u_w=0𝒩u_w=0.$$
(53)
As we have a single oscillator we can suppose that $`u(t)`$ is time symmetric with a suitable origin of time, i.e., $`\dot{u}`$ is a time-antisymmetric function of $`t`$ with the same period $`T=2\pi /\omega _\mathrm{b}`$ as $`u(t)`$: $`\dot{u}(t+T)=\dot{u}(t)`$, $`\dot{u}(0)=\dot{u}(T)=0`$. The properties of $`u_w=du/d\omega _\mathrm{b}`$, can be also obtained. We can write $`u`$ as a cosine Fourier series $`u(t)=_{k=0}^{\mathrm{}}c_k\mathrm{cos}(\omega _\mathrm{b}t)`$, with the coefficients $`c_k`$ depending on $`\omega _\mathrm{b}`$. Therefore, $`du/d\omega _\mathrm{b}=_{k=0}^{\mathrm{}}(dc_k/d\omega _\mathrm{b})\mathrm{cos}(\omega _\mathrm{b}t)+_{k=0}^{\mathrm{}}kt\mathrm{sin}(\omega _\mathrm{b}t)`$, or:

$$u_w\frac{du}{d\omega _\mathrm{b}}=\gamma (t)+\frac{t}{\omega _\mathrm{b}}\dot{u}(t)$$
(54)
$`\gamma (t)`$ being a time symmetric, periodic function of time with the same period $`T`$ as $`u`$. This functions constitute a base of the eigenspace corresponding to $`E=0`$.

On the other hand, $`𝒩`$ (52) is a periodic operator because it depends on a periodic function $`u`$. Its eigenfunctions are given by the Bloch theorem:

$$\xi (\theta ,t)=\chi (\mathrm{i}\theta ,t)\mathrm{exp}(\theta t/T),$$
(55)
with $`\chi `$ a T–periodic function of time and $`\theta =\theta (E)`$. If we consider only values of $`E`$ for which $`\theta `$ is real, the set of points $`(\theta ,E(\theta )`$ constitutes a band. If $`\xi `$ is a solution of $`𝒩\xi =E\xi `$, the fact that $`𝒩`$ is a real operator implies that $`\xi ^{}`$ is another solution, i.e., $`(\theta ,E)`$ also belongs to the band and this is symmetric with respect to $`\theta `$ and, therefore, $`E/\theta (0)=0`$. As $`𝒩`$ is time–symmetric $`\xi (\theta ,t)`$ is another, and has to be proportional to $`\xi ^{}`$ as there can be only two independent solutions.

The Floquet matrix $`F_E`$ is defined as:

$$\left[\begin{array}{c}\xi (T)\\ \dot{\xi }(T)\end{array}\right]=F_E\left[\begin{array}{c}\xi (0)\\ \dot{\xi }(0)\end{array}\right].$$
(56)
The eigenvalues $`\{\lambda _l\}`$ of $`F_E`$ are the Floquet multipliers, if we write them as $`\lambda _l=\mathrm{exp}\mathrm{i}\theta _l`$, with $`\theta `$ real or complex, $`\mathrm{i}\theta _l`$ are the Floquet exponents, and $`\theta _l`$, the Floquet arguments. An eigenfunction with $`\theta _l`$ real is bounded, if $`\theta _l=1`$ the eigenfunction has period $`T`$. The condition for linear stability of the solution $`u`$ is that all the Floquet arguments of $`F_0`$ are real.

It is easy to check some properties. The Bloch functions Eq. (55) are eigenfunctions of $`F_E`$ with Floquet exponent $`\theta `$, therefore $`F_E`$ is diagonalizable (over $`𝒞`$) with Floquet multipliers $`\mathrm{exp}(\pm \mathrm{i}\theta )`$ and Floquet arguments $`\pm \theta `$. However, at $`E=0`$, there is only one eigenfunction of $`F_0`$, $`\dot{u}`$, the other independent function of the subspace $`𝒩\xi =0`$, $`u_w`$, it is not an $`F_0`$–eigenfunction. That is $`F_0`$ has a degenerate Floquet multiplier $`1`$, or Floquet argument $`0`$ and can be transformed into a Jordan block. Another consequence of the latter is that $`(\theta ,E)=(0,0)`$ belongs to the band. We need to calculate the curvature of the band at this point.

Let us come back to the Bloch functions $`\xi (\theta ,t)`$ in Eq. (55) and $`\xi ^{}(\theta ,t)`$ corresponding to a pair of symmetric points $`(\pm \theta ,E)`$ of the band, and suppose that $`\theta 0`$ (and then $`E0`$). At $`\theta =0`$ the two functions collide in $`\xi (0,t)=\chi (0,t)=\chi ^{}(0,t)`$, therefore $`\xi (0,t)`$ is a real, $`T`$–periodic function. But we already now that at $`E=0`$, there is a $`T`$–periodic function solution of $`𝒩\xi =0`$, $`\dot{u}`$. Therefore, $`\xi (0,t)=\dot{u}`$ ($`𝒩`$ is lineal, so we can adjust the norm of $`\xi `$ in order that the latter equation is fulfilled).

We can reobtain the missing function that spans $`\mathrm{ker}(𝒩)`$. If we derive Eq. (52) with respect to $`\theta `$, we obtain $`𝒩\xi _\theta =E_\theta \xi +E\xi _\theta `$ and at $`E=0`$ ($`E_\theta (0)=0`$) we obtain: $`𝒩\xi _\theta (0,t)=0`$. Deriving $`\xi (\theta ,t)`$ in Eq. (55) we get:

$$\xi _\theta =\chi _\theta (\theta ,t)\mathrm{exp}(\mathrm{i}\theta /T)+\mathrm{i}\frac{t}{T}\chi (\theta ,t)\mathrm{exp}(\mathrm{i}\theta /T)=\chi _\theta \mathrm{exp}(\mathrm{i}\theta /T)+\mathrm{i}\frac{t}{T}\xi .$$
(57)
At $`\theta =0`$, $`\xi =\dot{u}`$ and $`\xi _\theta `$ has to be proportional at $`u_w`$, which we know is the other independent function. Therefore

$`\xi _\theta (0,t)=\chi _\theta +\mathrm{i}{\displaystyle \frac{t}{T}}\dot{u}=\mathrm{i}{\displaystyle \frac{\omega _\mathrm{b}}{T}}\left({\displaystyle \frac{T}{\mathrm{i}\omega _\mathrm{b}}}\chi _\theta +{\displaystyle \frac{t}{\omega _\mathrm{b}}}\dot{u}\right)=\mathrm{i}{\displaystyle \frac{\omega _\mathrm{b}^2}{2\pi }}\left({\displaystyle \frac{2\pi }{\omega _\mathrm{b}^2}}\chi _\theta +{\displaystyle \frac{t}{\omega _\mathrm{b}}}\dot{u}\right).`$
Comparing with Eq. (54) we obtain

$$\xi _\theta (0,t)=\mathrm{i}\frac{\omega _\mathrm{b}^2}{2\pi }u_w.$$
(58)
The following step is to express the curvature $`E_{\theta \theta }`$ of the band $`E(\theta )`$ at $`(\theta ,E)=(0,0)`$ in terms of $`u`$ and its derivatives. Let us derive Eq. (52) twice with respect to $`\theta `$, we obtain:

$$𝒩\xi _{\theta \theta }=E_{\theta \theta }\xi +2E_\theta \xi _\theta +E\xi _{\theta \theta }𝒩\xi _{\theta \theta }=E_{\theta \theta }\dot{u}\mathrm{at}\theta =0.$$
(59)
Multiplying to the right by $`\dot{u}`$ and integrating a period we get

$$E_{\theta \theta }_{t=0}^T\dot{u}^2t=_{t=0}^T\dot{u}𝒩\xi _{\theta \theta }t.$$
(60)
For any pair of functions $`f(t)`$ and $`g(t)`$

$`{\displaystyle _{t=0}^T}f𝒩g`$ $`=`$ $`{\displaystyle _{t=0}^T}(f\ddot{g}+fV^{\prime \prime }g)t=[f\dot{g}]_{t=0}^T[\dot{f}g]_{t=0}^T+{\displaystyle _{t=0}^T}(g\ddot{f}+gV^{\prime \prime }g)𝑑t`$
$`=`$ $`[f\dot{g}]_{t=0}^T[\dot{f}g]_{t=0}^T+{\displaystyle _{t=0}^T}g𝒩ft.`$
Applying the latter equation:

$`{\displaystyle _{t=0}^T}\dot{u}𝒩\xi _{\theta \theta }t=[\dot{u}\dot{\xi }_{\theta \theta }]_{t=0}^T[\ddot{u}\xi _{\theta \theta }]_{t=0}^T+{\displaystyle _{t=0}^T}(\xi _{\theta \theta }𝒩\dot{u})t.`$ (61)
Using that $`𝒩\dot{u}=0`$, $`\dot{u}(0)=\dot{u}(T)=0`$ and $`\ddot{u}(0)=\ddot{u}(T)`$ the previous equation leads to

$$_{t=0}^T\dot{u}𝒩\xi _{\theta \theta }t=\ddot{u}(T)\left(\xi _{\theta \theta }(0,T)\xi _{\theta \theta }(0,0)\right).$$
(62)
To calculate the last difference, let us derive Eq. (57) with respect to $`\theta `$. We obtain

$`\xi _{\theta \theta }`$ $`=`$ $`\chi _{\theta \theta }\mathrm{exp}(\mathrm{i}\theta t/T)+\mathrm{i}{\displaystyle \frac{t}{T}}\chi _\theta \mathrm{exp}(\mathrm{i}\theta t/T)+\mathrm{i}{\displaystyle \frac{t}{T}}\xi _\theta `$
$`=`$ $`\chi _{\theta \theta }\mathrm{exp}(\mathrm{i}\theta t/T)+\mathrm{i}{\displaystyle \frac{t}{T}}(\xi _\theta \mathrm{i}{\displaystyle \frac{t}{T}}\xi )+\mathrm{i}{\displaystyle \frac{t}{T}}\xi _\theta `$
$`=`$ $`\chi _{\theta \theta }\mathrm{exp}(\mathrm{i}\theta t/T)+2\mathrm{i}{\displaystyle \frac{t}{T}}\xi _\theta +{\displaystyle \frac{t^2}{T^2}}\xi `$
$``$ $`\xi _{\theta \theta }=\chi _{\theta \theta }+2\mathrm{i}{\displaystyle \frac{t}{T}}\xi _\theta +{\displaystyle \frac{t^2}{T^2}}\dot{u}\mathrm{at}\theta =0.`$
$`\chi `$ and therefore $`\chi _{\theta \theta }`$ are $`T`$–periodic functions, therefore $`\xi _{\theta \theta }(0,0)=\chi _{\theta \theta }(0,0)=\chi _{\theta \theta }(0,T)`$ and $`\xi _{\theta \theta }(0,T)\xi _{\theta \theta }(0,0)=2\mathrm{i}\xi _\theta (0,T).`$ Then

$$_{t=0}^T\dot{u}𝒩\xi _{\theta \theta }t=\ddot{u}(2\mathrm{i}\xi _\theta (0,T))=\frac{H}{u(T)}\mathrm{\hspace{0.17em}2}\mathrm{i}(\mathrm{i}\frac{\omega _\mathrm{b}^2}{2\pi }u_w(T))=\frac{\omega _\mathrm{b}^2}{\pi }\frac{H}{\omega _\mathrm{b}}.$$

Substituting in Eq. (60), and naming the action $`I=_{t=0}^T\dot{u}^2t`$ we finally obtain at $`\theta =0`$:

$$E_{\theta \theta }(0)=\frac{\omega _\mathrm{b}^2}{\pi I}\frac{H}{\omega _\mathrm{b}}.$$
(63)
$`I`$ is always positive, therefore, if the potential is hard ($`H/\omega _\mathrm{b}>0`$), the curvature $`E_{\theta \theta }`$ is negative and the band is tangent to the axis $`E=0`$ from below, if it is soft ($`H/\omega _\mathrm{b}<0`$) the curvature is positive and the band is tangent form above, as we wanted to demonstrate.

Note that this demonstration is also valid for the whole coupled system with the obvious changes: $`u[u_1,\mathrm{},u_N]^{}`$; $`V(u)V(u)+\epsilon W(u)`$, $`W(u)`$ being the coupling potential; $`V^{}(u)[V/u_1,\mathrm{},V/u_N]^{}`$; $`V^{\prime \prime }(u)`$ the matrix with elements $`^2V/u_nu_m`$ and analogously. We also need the additional constraint that the solution $`u(t)`$ is chosen time–symmetric and that the pair of eigenvalues with $`\theta =0`$ is isolated. This means that the band tangent at $`(\theta ,E)=(0,0)`$ continues having the same curvature properties depending on the hardness/softness of the potential