# Exact solutions for discrete breathers in a forced-damped chain

Oleg Gendelman

► To cite this version:

Oleg Gendelman. Exact solutions for discrete breathers in a forced-damped chain. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2013, 87, pp.062911. 10.1103/PhysRevE.87.062911. hal-01510820

HAL Id: hal-01510820

https://hal.science/hal-01510820

Submitted on 20 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Exact solutions for discrete breathers in a forced-damped chain

O. V. Gendelman

*Faculty of Mechanical Engineering, Technion—Israel Institute of Technology, 32000 Haifa, Israel*

Exact solutions for symmetric on-site discrete breathers (DBs) are obtained in a forced-damped linear chain with on-site vibro-impact constraints. The damping in the system is caused by inelastic impacts; the forcing functions should satisfy conditions of periodicity and antisymmetry. Global conditions for existence and stability of the DBs are established by a combination of analytic and numeric methods. The DB can lose its stability through either pitchfork, or Neimark-Sacker bifurcations. The pitchfork bifurcation is related to the internal dynamics of each individual oscillator. It is revealed that the coupling can suppress this type of instability. To the contrary, the Neimark-Sacker bifurcation occurs for relatively large values of the coupling, presumably due to closeness of the excitation frequency to a boundary of the propagation zone of the chain. Both bifurcation mechanisms seem to be generic for the considered type of forced-damped lattices. Some unusual phenomena, like nonmonotonous dependence of the stability boundary on the forcing amplitude, are revealed analytically for the initial system and illustrated numerically for small periodic lattices.

## I. INTRODUCTION

Discrete breathers (DBs) or intrinsic localized modes (ILMs) are well known in many nonlinear lattices [1,2]. Generally, they are exponentially localized (if a coupling between the neighbors in the lattice has a linear component) and can demonstrate remarkable stability also in two- and three-dimensional lattices [2]. Numerous systems which exhibit the DBs include chains of mechanical oscillators [3], superconducting Josephson junctions [4], nonlinear magnetic metamaterials [5], electrical lattices [6], michromechanical cantilever arrays [7], antiferromagnets [8], and many other physical systems [9].

In the majority of theoretical studies related to the DBs, Hamiltonian models are considered. Still, in many applications the damping cannot be neglected; in order to maintain the DB, one should compensate it by some kind of direct or parametric external forcing [2]. Experimentally, many of the DBs observed in the experiments were in fact created and maintained in the presence of damping and under homogeneous forcing from the external sources.

It is more or less easy to “explain” the nonlinear localization in a chain of forced-damped oscillators on a qualitative level. It is well known that single forced-damped nonlinear oscillators can exhibit stable steady-state responses with different amplitudes, depending on initial conditions [10]. If one or few oscillators in the lattice are excited at high amplitude, and if all the others are at low amplitude, and coupling between the oscillators is weak enough to preserve this structure—one obtains an example of strongly localized excitation in conditions of a homogeneous forcing. However, quantitative description of such breathers remains a major challenge. Lack of Hamiltonian structure changes the properties of the DBs. Instead of a continuous family of localized solutions, one expects to obtain a discrete set of attractors. Consequently, many of the methods devised for computation and analysis of the Hamiltonian DBs are not applicable in

forced-damped systems. Even more important, the DBs in the forced-damped systems can have some qualitative properties absent in their Hamiltonian counterparts. For instance, one can observe stable moving DBs in forced-damped systems [11]. In Hamiltonian lattices, the existence of such moving stable DBs is generally denied [2], with the known exception of the integrable Ablowitz-Ladik model [12]. Also, there exist some special nonintegrable Hamiltonian discrete models where the DBs can be computed exactly [13]. To the best of the author’s knowledge, no such exact solutions are known in forced-damped discrete lattices. Even in continuous nonlinear models with forcing and damping, some exact solutions for the breathers were obtained only for a handful of special cases [14]. This paper is devoted exactly to this problem and suggests a model, which allows one to derive exact solutions for the forced-damped DBs and to study some of their properties.

We are going to demonstrate that the exact solutions can be derived for the DBs in a vibro-impact (VI) chain. Dynamical systems involving impacts have an important peculiarity—they exhibit extreme (actually, the strongest possible) nonlinearity, but the latter reveals itself only at the moments of impacts; between the impacts the system obeys linear equations of motion, if other sources of nonlinearity are absent. Consequently, the VI models can offer a relatively simple description of complicated nonlinear phenomena. Celebrated examples of this sort are problems of bouncing ball [15] and kicked rotor [16], but many other tractable VI models are known and used for description of realistic physical systems, such as Bose-Einstein condensates [17]. One of the few exact solutions for DBs in Hamiltonian systems also uses the VI model [18]; in the current paper, we extend this result for the forced-damped case.

## II. DESCRIPTION OF THE MODEL AND ANALYTIC SOLUTION

The model used here is a homogeneous chain with linear nearest-neighbor interactions; in addition, each particle can move only between inelastic impacts constraints. Such model is topologically equivalent to smooth models with linear
<figure><img src="image_2.png" /><figcaption>FIG. 1. Sketch of the model system.</figcaption></figure>

interaction and on-site potential, widely used in investigations devoted to the DBs [7]. A sketch of the model system is presented in Fig. 1.

Each particle is subject to the same external forcing. The forcing function is considered to be 2π-periodic and antisymmetric, as specified below in the equations of motion. Between the impacts, displacement $u_n$ of the nth particle is described as

$$ 
\begin{aligned}
\ddot{u}_n + \gamma(2u_n - u_{n-1} - u_{n+1}) &= F(t), \\
F(t) &= F(t + 2\pi), \quad (1) \\
F(t + \pi) &= -F(t).
\end{aligned}
 $$

Without affecting the generality, the clearance $\Delta$ and the particle mass are set to unity and the only parameter which characterizes the chain is a coupling coefficient $\gamma$. Equation (1) is valid between successive impacts, i.e., as $|u_n| < 1$ for all particles. If the particle with number $n$ impacts the constraint at time instance $t_{i,n}$, the following conditions are satisfied:

$$ 
\begin{aligned}
t &= t_{i,n}; u_n(t_{i,n}) = \pm 1; \\
\dot{u}_n(t_{i,n} + 0) &= -k\dot{u}_n(t_{i,n} - 0), \quad 0 \le k \le 1.
\end{aligned}
 $$

Here $k$ is a restitution coefficient. The notation $\pm 0$ is used to denote the particle's velocities immediately after (+) and before (-) the impact. Values of $k<1$ correspond to inelastic impacts; it is the only source of damping in the model. The model is not solvable in a general form. Therefore, we are going to look for partial solutions for the DBs. In this paper, the simplest on-site symmetric DB is considered. Our motivation for this choice is that the DBs of this type are the most explored ones, both experimentally and numerically. Hence, we adopt that responses of all particles are periodic and symmetric, and only a particle with $n=0$ impacts each of its constraints symmetrically with period $2\pi$:

$$ 
\begin{aligned}
& u_n(t) = u_n(t + 2\pi m) = -u_n[t + \pi(2m + 1)], \\
& t_{i,0} = \varphi + \pi m, \quad u_0(\varphi) = 1, \\
& \dot{u}_0(\varphi + 0) = -k\dot{u}_0(\varphi - 0), \quad u_0(\varphi + \pi) = -1, \\
& \dot{u}_0(\varphi + \pi + 0) = -k\dot{u}_0(\varphi + \pi - 0).
\end{aligned}
 $$

Here $\varphi$ is a “phase lag” of the impacts with respect to the external forcing. The effect of each impact is equivalent to a transfer of a certain amount of momentum to the impacting particle. Consequently, the solution we look for

should satisfy the following modification of Eq. (1) (in terms of distributions):

$$ 
\begin{aligned}
& \ddot{u}_n + \gamma(2u_n - u_{n-1} - u_{n+1}) \\
& = F(t) + 2p\delta_{n0} \sum_{j=-\infty}^{\infty} \{\delta[t - \varphi + \pi(2j + 1)] \\
& \qquad - \delta(t - \varphi + 2\pi j)\},
\end{aligned}
 $$

where $\delta(t)$ is the Dirac delta function, $\delta_{nm}$ is the Kronecker symbol. $2p$ is the yet unknown momentum transferred to the “central” particle in the course of each impact. For further analysis, the variables are changed as follows:

$$ u_n = v_n + f(t), \quad \ddot{f}(t) = F(t). \quad (5) $$

Due to the antisymmetry of $F(t)$, it is always possible to find a unique function $f(t)$ satisfying Eq. (5) and the antisymmetry condition $f(t + \pi) = -f(t)$. In terms of the new variables, the external forcing term on the right-hand side of Eq. (4) disappears:

$$ 
\begin{aligned}
& \ddot{v}_n + \gamma(2v_n - v_{n-1} - v_{n+1}) \\
& = 2p\delta_{n0} \sum_{j=-\infty}^{\infty} \{\delta[t - \varphi + \pi(2j + 1)] - \delta(t - \varphi + 2\pi j)\}.
\end{aligned}
 $$

In Eq. (6) the only remaining forcing term is related to the impacts. The stationary regime is possible only if the impact is symmetric in terms of variables $v_n(t)$:

$$ \dot{v}_0(\varphi - 0) = -\dot{v}_0(\varphi + 0) = p. $$

With the account of this relationship, Eq. (2) for the inelastic impact is rewritten as

$$ 
\begin{aligned}
& v_0(\varphi) + f(\varphi) = 1; \\
& \dot{v}_0(\varphi + 0) = -p + \dot{f}(\varphi) \\
& = -k[p + \dot{f}(\varphi)] \Rightarrow \dot{f}(\varphi) = p \frac{1-k}{1+k}.
\end{aligned}
 $$

Equation (7) has three unknowns, $p, \varphi$, and $v_0(\varphi)$. To define the solution completely, an additional relationship between $p, \varphi$, and the value of $v_0$ in time instance of the impact $v_0(\varphi)$ is required. One obtains it by rewriting Eq. (6) with the help of a formal Fourier series:

$$ 
\begin{aligned}
& \ddot{v}_n + \gamma(2v_n - v_{n-1} - v_{n+1}) \\
& = -\frac{4p\delta_{n0}}{\pi} \sum_{l=0}^{\infty} \cos[(2l + 1)(t - \varphi)].
\end{aligned}
 $$

Equation (8) represents a common linear chain with an external forcing applied to the particle $n=0$. The steady-state solution will be spatially localized if all frequencies in the forcing term lie in an attenuation zone of the chain. The latter requirement will be satisfied if $\gamma \in [0, 0.25)$. In this case the steady-state solution of Eq. (8) may be presented in the following form [17]:

$$ v_n(t) = \frac{4p}{\pi} \left(\frac{-1}{2\gamma}\right)^{|n|} \sum_{l=0}^{\infty} \frac{[(2l+1)^2 - 2\gamma - \sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}]^{|n|}}{\sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}} \cos[(2l+1)(t-\varphi)]. \quad (9) $$

For the sake of completeness, the derivation details of Eq. (9) are presented in Appendix A.
Fourier series (9) converges for all $n$ and for $\gamma \in [0,0.25)$, i.e., for all values of coupling, for which the forcing frequency stays in the attenuation zone. An explicit expression for $v_0(\varphi)$ is written as

$$v_0(\varphi) = p\chi(\gamma), \chi(\gamma) = \frac{4}{\pi} \sum_{l=0}^{\infty} \frac{1}{\sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}}. \quad (10)$$

Taking into account Eq. (10), the system (7) can be closed:

$$p\chi(\gamma) + f(\varphi) = 1, \quad \dot{f}(\varphi) = p\frac{1-k}{1+k}. \quad (11)$$

Then, the exact solution for the symmetric one-site DB is expressed as

$$u_n(t) = f(t) + \frac{4p}{\pi} \left(\frac{-1}{2\gamma}\right)^{|n|} \sum_{l=0}^{\infty} \frac{[(2l+1)^2 - 2\gamma - \sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}]^{|n|} \cos[(2l+1)(t-\varphi)]}{\sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}}. \quad (12)$$

Parameters $p$ and $\varphi$ can be determined from Eq. (11) for each specific choice of $F(t)$.

In order to study the DB described by Eqs. (11) and (12) in more depth, the case of a simple harmonic forcing is considered. In this case, one can solve Eq. (11) explicitly:

$$F(t) = -a \cos t \Rightarrow f(t) = a \cos t; \quad -a \cos \varphi = p\chi(\gamma) - 1, \quad -a \sin \varphi = pq, \quad q = \frac{1-k}{1+k}, \\ p = \frac{\chi(\gamma) \pm \sqrt{\chi^2(\gamma)a^2 - q^2(1-a^2)}}{\chi^2(\gamma) + q^2}, \quad \varphi = -\arcsin \frac{pq}{a}. \quad (13)$$

Here $a$ is the amplitude of the harmonic forcing. A stable response corresponds to a positive sign in expression (13) for $p$. From (13) one obtains the condition for minimum forcing amplitude necessary to support the DB solution, as a function of the restitution coefficient and the coupling:

$$a \ge \frac{q}{\sqrt{\chi^2(\gamma) + q^2}}. \quad (14)$$

To ensure self-consistency of the solution starting from Eq. (4), two conditions should be satisfied:

(1) The particle with $n=0$ should impact the barriers only at the time instances $t = \varphi, \varphi + \pi$;

(2) The particles with $n \ne 0$ should never reach the impact constraints.

Condition (1) is satisfied for solution (12)-(13) for not-too-small values of the restitution coefficient. It is demonstrated in Appendix B that for the values of $k \ge 0.1$ the solution is self-consistent for all possible values of the other parameters. To satisfy Condition (2), it is enough to take care about this condition for particles with $n = \pm 1$, since the amplitudes of all particles with $|n| > 1$ would be even smaller (a proof of this fact and more detailed exploration of the solution self-consistency are also presented in Appendix B). Then, for any time instance, one should obtain

$$|u_1(t)| = \left| a \cos t - \frac{2p}{\gamma \pi} \sum_{l=0}^{\infty} \frac{(2l+1)^2 - 2\gamma - \sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}}{\sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}} \cos[(2l+1)(t-\varphi)] \right| < 1. \quad (15)$$

Conditions (14) and (15) determine the zone of existence of the DB in the space of parameters. This zone is illustrated in Fig. 4, together with the stability thresholds for the DBs.

### III. STABILITY ANALYSIS

The DB solution (12)-(13) is $2\pi$-periodic by construction, and thus its stability may be established by analysis of eigenvalues of a monodromy matrix computed on a single time period [19]. Such matrix cannot be computed for an infinite chain; therefore, it is commonly accepted to compute it for a finite chain and to check whether the stability properties depend on the number of particles [11]. In this paper we adopt the same approach and consider the system with $N$ particles and periodic boundary conditions. For systems with smooth nonlinearity, such monodromy matrices are usually computed numerically, by solving a complete set of $N$ ordinary

differential equations (ODEs) for $2N$ independent sets of initial conditions [11,20]. However, for the vibro-impact chain considered here, in addition to availability of exact solutions, one can compute the monodromy matrix in a much simpler manner. A state vector of the system is defined as

$$\mathbf{V} = \begin{bmatrix} u_0 \\ \vdots \\ u_{N-1} \\ w_0 \\ \vdots \\ w_{N-1} \end{bmatrix}. \quad (16)$$
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 2. (Color online) Zones of existence and stability for the symmetric DB on the (y, a) parameter plane for two different values of the restitution coefficient: (a) k = 0.8; (b) k = 0.4. Crosses denote the pitchfork bifurcation line and diamonds denote the line of Neimark-Sacker bifurcation. Transitions corresponding to illustrations in Figs. 2 and 3 are marked by arrows in (a). “S” and “US” denote zones of stability and instability, respectively.</figcaption></figure>

We adopt $w_k = \dot{u}_k$, $k = 0, \dots, N$. The periodic boundary conditions are imposed as

$$u_0 = u_N, \quad w_0 = w_N. \tag{17}$$

In the exact solution (12)-(13) the impacts occur only at the site $n = 0$ in time instances $t = \varphi, \varphi + \pi$. Consequently, between impacts, the time evolution of the state vector is described by a linear equation. With account of (1) and (16), this equation is written as

$$\dot{\mathbf{V}} = \mathbf{A}\mathbf{V} + \mathbf{F}, \quad \mathbf{A} = \begin{pmatrix} 0_{N \times N} & I_{N \times N} \\ \gamma L_{N \times N} & 0_{N \times N} \end{pmatrix}, \quad \mathbf{F} = F(t) \begin{bmatrix} \mathbf{0}_{N \times 1} \\ \mathbf{I}_{N \times 1} \end{bmatrix} \\ \varphi < t < \varphi + \pi, \quad \varphi + \pi < t < \varphi + 2\pi. \tag{18}$$

Here

$$\mathbf{L}_{N \times N} = \begin{pmatrix} -2 & 1 & 0 & \cdots & 0 & 1 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & \ddots & \ddots & 0 & \vdots \\ \vdots & 0 & \ddots & \ddots & 1 & 0 \\ 0 & \vdots & \ddots & 1 & -2 & 1 \\ 1 & 0 & \cdots & 0 & 1 & -2 \end{pmatrix}$$

is the discrete Laplacian adjacency matrix accounting for the linear coupling.

Due to the linearity of Eq. (18), evolution of small perturbations of the state vector V between the impacts is described by the exponent of matrix A. At the moments of impact, the evolution of small perturbations of the state vector is described by the so-called saltation matrix [21], which takes into account both the perturbation of the state vector at the point of discontinuity and the infinitesimal variation of the impact time instance. For the system under consideration, this matrix at time instance $t = \varphi$ in the lowest-order (linearized) approximation [21] may be written as

$$\mathbf{S} = \begin{pmatrix} -k & 0 & \cdots & 0 \\ 0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & 1 \end{pmatrix} \qquad 0_{N \times N} \\ \begin{pmatrix} \frac{\omega(1+k)}{w_0(\varphi-0)} & 0 & \cdots & 0 & -k & 0 & \cdots & 0 \\ 0 & 0 & \ddots & \vdots & 0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & \vdots & \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & 0 & \cdots & 0 & 1 \end{pmatrix} \tag{19}$$

where $\omega = \ddot{u}_0(\varphi - 0)$. The complete monodromy matrix at the period $2\pi$ is thus expressed as

$$\mathbf{M} = (\mathbf{Q}\mathbf{S})^2, \quad \mathbf{Q} = \exp(\pi\mathbf{A}). \tag{20}$$

Values of velocity and acceleration of the impacting particle in the moment of impact are taken from the exact solution (12)-(13). Expression (20) still requires the numeric computation of the exponent matrix and the multiplication of the matrices, but all of these tasks are much easier from a computational viewpoint than the solution of N second-order ODEs with 2N independent sets of the initial conditions.

The loss of stability has been recorded numerically in the space of parameters, when the eigenvalue of the monodromy matrix (or couple of the complex conjugate eigenvalues) with the largest absolute value crossed a threshold of 1.005. For all explored examples, it was found that for $N \ge 400$ the critical parameter values do not depend on N within a relative margin of about 0.1%. Thus, it is possible to suggest that these results on stability of the DBs are also relevant for the infinite chain.
<figure><img src="image_5.png" /><figcaption>FIG. 3. Evolution of eigenvalues of the monodromy matrix corresponding to Neimark-Sacker bifurcation. N = 400, a = 0.2, k = 0.8; (a) γ = 0.09, (b) γ = 0.11.</figcaption></figure>

Two generic scenarios of the loss of stability were revealed. The first one corresponds to the transition of two complex conjugate eigenvalues through the unit circle and is related to the Neimark-Sacker bifurcation. The other scenario corresponds to the passage of a single eigenvalue through unity. This scenario corresponds to a pitchfork bifurcation; the latter results in appearance of a pair of stable asymmetric DBs. Typical structure of the eigenvalues of the monodromy matrix for both of these bifurcation scenarios is presented in Figs. 2 and 3.

Therefore, for instance, it is possible to infer from Fig. 2 that for the set of parameters corresponding to case (a) the symmetric DB will be stable, and for case (b) it will be unstable.

Examples for zones of existence and stability for the DBs are presented in Fig. 4. The lower boundary of the zones of existence is described by Eq. (14); the upper one is described by Eq. (15).

One can see that for lower values of the restitution coefficient, the DB exists for a narrower range of amplitudes of the external forcing, but it is stable in a wider range of the values of coupling coefficient. The pitchfork bifurcation in the upper left corner is similar to the loss of symmetry observed in a single vibro-impact oscillator [22]; naturally, this case is equivalent to zero coupling. A similar pitchfork bifurcation, related to a “loss of symmetry,” is also known in the regular forced-damped Duffing oscillator [10]. Thus, this scenario of the loss of stability is caused by the internal dynamics of each individual oscillator. Our computations reveal an interesting

<figure><img src="image_6.png" /><figcaption>FIG. 4. Evolution of eigenvalues of the monodromy matrix corresponding to pitchfork bifurcation. N = 400, a = 0.8, k = 0.8; (a) γ = 0.05, (b) γ = 0.045.</figcaption></figure>

and, it seems, previously unknown fact: The coupling between the oscillators can suppress this scenario of the loss of stability and stabilize the symmetric DBs.

It is also possible to conjecture that the Neimark-Sacker bifurcation in the lower right corner is related to four-wave interaction with the boundary of the propagation zone [7]. One can also see that the size of the interval of existence of the DBs shrinks and approaches zero as the frequency approaches the boundary of the propagation zone. This observation also correlates with numeric and experimental findings presented in Ref. [7], where the amplitude of the DB was seen to decrease as the forcing frequency gets closer to an optical band. These results allow us to conjecture that both bifurcation mechanisms described above are generic for forced lattices with nonlinear on-site potentials, and not only relevant for the VI chain studied here.

In Fig. 4(a) one can see that for relatively low damping (k = 0.8) the stability threshold on the parameter plane is not monotonic. It is not clear whether this property is generic. In the next section, this finding is illustrated by direct numeric simulations.

## IV. NUMERIC SIMULATIONS

Usually, experiments with the DBs are performed with a relatively small number of particles or nodes. Therefore, in order to illustrate the analytic findings, numeric simulations are performed for a system with N = 20 particles and periodic boundary conditions. Additionally, the inelastic impacts were replaced by strongly nonlinear smooth force and damping, which exactly correspond to the Newtonian model of the inelastic impact in the limit case specified below [23]. Also, small on-site linear damping with coefficient μ = 0.005 was included. Thus, the following smooth system has been simulated:

$$ \ddot{u}_n + \gamma(2u_n - u_{n-1} - u_{n+1}) + \mu\dot{u}_n + \lambda(2m+1)u_n^{2m}\dot{u}_n + (2m+1)u_n^{4m+1} = a \cos t, \\ u_0 = u_N. \qquad (21) $$

Here m is a natural number. In Ref. [23] it was proved that the two last terms on the left-hand side of Eq. (21) in the limit m → ∞ faithfully model the inelastic impact with a velocity-independent restitution coefficient. The value of this coefficient is related to parameter λ as follows [23]:

$$ k = \exp\left(\frac{-\pi\lambda}{2\sqrt{1-\lambda^2/4}}\right). \qquad (22) $$

The system has been simulated numerically by the Rosenbrock method suitable for stiff systems. It was observed that for values of m > 5 the simulation results ceased to depend on the exact choice of m.

Convergence to the exact DB solution has been simulated for a = 0.1, k = 0.8, and γ = 0.08. The system has been initialized by initial velocity ˙u₀(0) = 3; remaining initial conditions were set to zero. The simulated response of the system after an initial transient is compared to the exact solution (12)-(13) for the vibro-impact counterpart in Figs. 5(a)–5(c).

One can conclude that the DB for the given set of parameters is a genuine attractor of the dynamics and that
<figure><img src="image_7.png" /><figcaption>FIG. 5. (Color online) Comparison between exact solution for infinite chain (12)-(13) and numeric simulation of system (22) for $N = 20$; periodic boundary conditions, $a = 0.1, k = 0.8, \gamma = 0.08$, and $\mu = 0.005$. The plots present almost indistinguishable exact solution (red, solid line) and numeric result (black, circles) for (a) $u_0(t)$; (b) $u_1(t)$; (c) $u_{10}(t)$.</figcaption></figure>

the exact solution (12)-(13) is robust in a sense that it nicely approximates similar excitation in a smooth system with a relatively small number of degrees of freedom.

In order to illustrate the nonmonotonic dependence of stability on the coupling revealed in Fig. 4(a), system (22) was simulated for initial conditions exactly corresponding to solution (12)-(13), the restitution coefficient $k = 0.8$, coupling $\gamma = 0.104$, and three different values of the external forcing $-a = 0.1, 0.2, 0.7$. From Fig. 4(a) it is clear that all these values of the forcing and coupling are rather close to the boundary of stability. The results of this simulation are presented in Fig. 6.

For cases (a) and (c), the system settles in the attractor corresponding to the symmetric DB after some initial transient. For case (b), the system preserves the localization for some time, but then the response becomes delocalized and the

<figure><img src="image_8.png" /><figcaption>FIG. 6. (Color online) Numeric simulation of the DBs in small system with smoothened potentials and damping, $N = 20, k = 0.8, \gamma = 0.104$; (a) $a = 0.1$, (b) $a = 0.2$, (c) $a = 0.7$. Time series for five central particles ($-2 \le n \le 2$) are presented.</figcaption></figure>

DB structure is eventually lost. The results of simulations completely correspond to predictions in Fig. 4(a), including the nonmonotonous behavior of the stability threshold.

One can mention that in the simulation presented in Fig. 6(b) the DB is destroyed. It does not turn out to be quasiperiodic, as one would expect from the generic Neimark-Sacker bifurcation scenario described above. This result suggests more careful study of the DB structure transformation when the stability boundary is crossed. For this sake, the symmetric DB was created according to the exact solution in the region of stability. Then, the amplitude of the forcing
<figure><img src="image_9.png" /></figure>

(a)

<figure><img src="image_10.png" /></figure>

(b)

<figure><img src="image_11.png" /></figure>

(c)

FIG. 7. (Color online) Comparison between exact solution for the symmetric DB in the infinite chain (12)-(13) (red, solid line) and numeric simulation of system (22) (black, circles) for $N = 20$; periodic boundary conditions, $a = 0.82$, $k = 0.8$, $\gamma = 0.04$, and $\mu = 0.005$. The time series are depicted for (a) $u_0(t)$; (b) $u_1(t)$; (c) $u_{10}(t)$.

was varied (at the rate of about $10^{-3}$ with respect to the forcing frequency), in a way that the system has crossed the stability threshold. It was revealed that the localized DBs in such conditions could persist for a rather large variation of the forcing amplitude; however, their structure undergoes significant changes. The simulations were performed for both stability boundaries (corresponding to the pitchfork and the Neimark-Sacker bifurcations). The time series for the final DB shape are presented in Figs. 7 and 8, respectively.

The time series obtained in both simulations correspond to the DBs, but the latter are no more symmetric. Analytic solutions for the unstable symmetric DBs are depicted in these figures by solid lines for the sake of comparison. In Fig. 7 one can observe that the numeric solution is still periodic, but asymmetric. The solution presented in Fig. 8 is no more $2\pi$-periodic, but still spatially localized. One can suppose that this solution might be either quasiperiodic or phase locked with some large period. Both results presented in Figs. 7 and 8 completely conform to the predictions of the stability analysis presented above. A study of precise structure and stability properties of the asymmetric and quasiperiodic DBs

<figure><img src="image_12.png" /></figure>

(a)

<figure><img src="image_13.png" /></figure>

(b)

<figure><img src="image_14.png" /></figure>

(c)

FIG. 8. (Color online) Comparison between exact solution for the symmetric DB in the infinite chain (12)-(13) (red, solid line) and numeric simulation of system (22) (black, circles) for $N = 20$; periodic boundary conditions, $a = 0.15$, $k = 0.8$, $\gamma = 0.104$, and $\mu = 0.005$. The time series are depicted for (a) $u_0(t)$; (b) $u_1(t)$; (c) $u_{10}(t)$.

requires many additional analytic and numeric efforts and will be presented elsewhere.

V. CONCLUDING REMARKS

As was demonstrated above, the vibro-impact model allows derivation of exact analytic solutions for discrete breathers in forced-damped essentially nonlinear chain, without any simplifications or approximations. The zone of existence and stability of the DB solution in the space of parameters has been computed. Two mechanisms of the loss of stability (pitchfork and Neimark-Sacker bifurcations) were revealed. Existence of the asymmetric and quasiperiodic DBs was predicted analytically and verified numerically. It is interesting to mention that the increase of the coupling can both stabilize and destabilize the symmetric DB in different zones on the plane of parameters and for different bifurcation mechanisms.
Also, nonmonotonic behavior of the stability boundary in the space of parameters was revealed analytically for the infinite chain and illustrated numerically for small smooth systems with periodic boundary conditions.

One should note that the numeric simulations were performed for the model with smooth, rather that vibro-impact, on-site potentials. The results were very similar to the theoretical counterparts obtained for the vibro-impact model. This coincidence allows one to suggest that the results concerning the behavior of the on-site DBs (generic bifurcations, “topography” of the existence and stability zones on the plane of parameters) are in fact relevant for a wider class of systems with topologically similar potentials of interaction. This suggestion will be explored in a future work.

## ACKNOWLEDGMENTS

The author is very grateful to Professor Yuli Starosvetsky for useful discussions and suggestions and to the anonymous reviewer for great help in improving the paper.

## APPENDIX A

In this Appendix we present the detailed derivation of Eq. (9). It is natural to look for the steady-state solution of Eq. (8) in the form of the Fourier series:

$$v_n(t) = \sum_{l=0}^{\infty} v_{n,l} \cos[(2l+1)(t-\varphi)]. \quad (A1)$$

Substituting this expression into Eq. (8), one obtains

$$-(2l+1)^2 v_{n,l} + \gamma(2v_{n,l} - v_{n-1,l} - v_{n+1,l}) = -\frac{4p\delta_{n0}}{\pi}. \quad (A2)$$

For $n \neq 0$ Eq. (A2) is homogeneous. If one requires vanishing of the solution as $n \to \pm\infty$, then the appropriate solution of (A2) is expressed as

$$v_{n,l} = v_{0,l} \xi^{-|n|}, \\ \xi = - \frac{(2l+1)^2 - 2\gamma - \sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}}{2\gamma}. \quad (A3)$$

For $n = 0$ one substitutes Eq. (A3) into Eq. (A2) and obtains the following expression for $v_{0,l}$:

$$v_{0,l} = \frac{4p}{\pi\sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}}. \quad (A4)$$

Substituting Eqs. (A3) and (A4) into Eq. (A1), one obtains Eq. (9).

## APPENDIX B

In this Appendix the proof of the self-consistency of analytic solution (12)-(13) for not-too-small values of the restitution coefficient $k$ is presented for the case of simple harmonic forcing. The proof is based on two statements, which will be analyzed separately:

(1) The particle with $n = 0$ undergoes impacts with the constraints at time instances $t = \varphi, \varphi + \pi$ and does not impact the constraints at any other time instance;

(2) The upper boundary for existence of the single-site DB is limited by the curve corresponding to the condition

$$\max_t [u_1(t)] = 1. \quad (B1)$$

This condition implies

$$\max_t [u_n(t)] \le 1 \quad (B2)$$

for all values of $|n| > 1$.

If Statements 1 and 2 are correct, they guarantee that no particle besides the one with $n = 0$ is engaged in the impacts with the constraints, and the central particle meets the constraints only at $t = \varphi, \varphi + \pi$. Thus, the transition from Eqs. (1) and (2) to Eq. (4) and all analytic treatment that follows would be justified.

Let us first consider Statement 1. It is equivalent to the statement that

$$|u_0(t)| < 1 \text{ for } t \neq \varphi + \pi j, j \in \mathbb{Z}. \quad (B3)$$

Due to the symmetry of the considered solution, it is enough to prove (B3) for $\varphi < t < \varphi + \pi$. If one recalls that $u_0(\varphi) = 1, u_0(\varphi + \pi) = -1$, it is clear that it is sufficient to demonstrate a monotonous decrease of $u_0(t)$ for $\varphi < t < \varphi + \pi$.

From Eq. (12), velocity of the particle with $n = 0$ is written as

$$V_0(t) = \dot{u}_0(t) = -a \sin t - \frac{4p}{\pi} \sum_{l=0}^{\infty} \frac{\sin[(2l+1)(t-\varphi)]}{\sqrt{(2l+1)^2 - 4\gamma}}. \quad (B4)$$

Due to a discontinuity of the velocity at the impacts, series (B4) does not converge uniformly. In order to obtain suitable estimations of the velocity, this series should be presented in somewhat modified shape:

$$\begin{align*} & \frac{4p}{\pi} \sum_{l=0}^{\infty} \frac{\sin[(2l+1)(t-\varphi)]}{\sqrt{(2l+1)^2 - 4\gamma}} \\ &= \frac{4p}{\pi} \left[ \sum_{l=0}^{\infty} \frac{\sin[(2l+1)(t-\varphi)]}{2l+1} + \sum_{l=0}^{\infty} \sin[(2l+1)(t-\varphi)] \left\{ \frac{1}{\sqrt{(2l+1)^2 - 4\gamma}} - \frac{1}{2l+1} \right\} \right] \\ &= p + \frac{16\gamma p}{\pi} \sum_{l=0}^{\infty} \frac{\sin[(2l+1)(t-\varphi)]}{(2l+1)(2l+1 + \sqrt{(2l+1)^2 - 4\gamma})(2l+1)^2 - 4\gamma} \end{align*}$$
$$
\begin{align*}
&= p + \frac{16\gamma p}{\pi} \left( \frac{\sin(t - \varphi)}{\sqrt{1 - 4\gamma}[1 + \sqrt{1 - 4\gamma}]} + \sum_{l=1}^{\infty} \frac{\sin[(2l + 1)(t - \varphi)]}{(2l + 1)[2l + 1 + \sqrt{(2l + 1)^2 - 4\gamma}]\sqrt{(2l + 1)^2 - 4\gamma}} \right) \\
&\ge p + \frac{16\gamma p \sin(t - \varphi)}{\pi} \left( \frac{1}{\sqrt{1 - 4\gamma}(1 + \sqrt{1 - 4\gamma})} - \sum_{l=1}^{\infty} \frac{1}{[2l + 1 + \sqrt{(2l + 1)^2 - 4\gamma}]\sqrt{(2l + 1)^2 - 4\gamma}} \right) \\
&\ge p + \frac{16\gamma p \sin(t - \varphi)}{\pi} \left( \frac{1}{2} - \frac{1}{8} \sum_{l=1}^{\infty} \frac{1}{l^2} \right) = p + \frac{16\gamma p \sin(t - \varphi)}{\pi} \left( \frac{1}{2} - \frac{\pi^2}{48} \right). \tag{B5}
\end{align*}
$$

The first inequality in (B5) explores the fact that $\sin[(m+1)x] \le (m+1)\sin x$ for all non-negative whole numbers $m$ and for $x \in [0, \pi]$. This fact easily follows from recurrence relationships for the second-order Chebyshev polynomials [24]. In the last transformation we took into account that $\gamma \in [0, 0.25)$. Using estimation (B5), one obtains from Eq. (B4)

$$V_0(t) = -a \sin t - \frac{4p}{\pi} \sum_{l=0}^{\infty} \frac{\sin[(2l+1)(t-\varphi)]}{\sqrt{(2l+1)^2 - 4\gamma}} \le -a \sin t - p - \frac{16\gamma p \sin(t-\varphi)}{\pi} \left(\frac{1}{2} - \frac{\pi^2}{48}\right) \le -a \sin t - p. \quad (B6)$$

The last transformation holds since $\sin(t - \varphi)$ is non-negative for $\varphi \le t \le \varphi + \pi$. It follows from (B6) that $u_0(t)$ will monotonously decrease in the interval $\varphi < t < \varphi + \pi$ for all values of $\gamma \in [0, 0.25)$ if it monotonously decreases for $\gamma = 0$. It is clear from (B6) that the latter condition will be fulfilled if $p > a$ in the interval $1 \ge a \ge \frac{q}{\sqrt{\chi^2(0)+q^2}}$; from expression (13), it is easy to derive that the inequality holds for all $k \ge 0.1$. So, one can conclude that Statement 1 is correct for not-too-small values of the restitution coefficient. It should be mentioned that for very small values of $k$ Statement 1 really turns out to be wrong—that happens, for instance, for the set of parameters $k = 0.005$, $a = 0.99$, $\gamma = 0$. Therefore, although it might be possible to tighten the estimations used in inequalities (B5) and (B6), still the consistency of solution (12)-(13) will be violated for small values of the restitution coefficient $k$.

The proof of Statement 2 essentially relies on the fact that the harmonic component with unit frequency is dominant in all sums $u_n(t)$ for $n \ne 0$.

Solution (12)-(13) may be transformed as follows:

$$
\begin{align}
u_n(t) &= a \cos t + \frac{4p}{\pi} \left(\frac{-1}{2\gamma}\right)^{|n|} \sum_{l=0}^{\infty} \frac{((2l+1)^2 - 2\gamma - \sqrt{(2l+1)^4 - 4\gamma(2l+1)^2})^{|n|} \cos[(2l+1)(t-\varphi)]}{\sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}} \nonumber \\
&= A_n(t) + R_n(t). \tag{B7} \\
A_n(t) &= a \cos t + \frac{(-1)^n 4p \xi_0^n}{\pi \sqrt{1-4\gamma}} \cos(t-\varphi). \nonumber \\
R_n(t) &= \frac{(-1)^n 4p}{\pi} \sum_{l=1}^{\infty} \frac{\xi_l^n \cos[(2l+1)(t-\varphi)]}{\sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}}. \nonumber
\end{align}
$$

In this expression $A_n(t)$ is the component with the unit frequency, and $R_n(t)$ denotes the sum of all other harmonics,

$$
\begin{align}
\xi_l &= \frac{(2l + 1)^2 - 2\gamma - \sqrt{(2l + 1)^4 - 4\gamma(2l + 1)^2}}{2\gamma} \nonumber \\
&= \frac{2\gamma}{(2l + 1)^2 - 2\gamma + \sqrt{(2l + 1)^4 - 4\gamma(2l + 1)^2}} = \frac{2l + 1 - \sqrt{(2l + 1)^2 - 4\gamma}}{2l + 1 + \sqrt{(2l + 1)^2 - 4\gamma}}. \tag{B8}
\end{align}
$$

Note that we consider only $n>0$; the opposite case is similar due to the symmetry of the exact solution. Using the former notations, Eq. (15) can be transformed in the following way:

$$\max_t [u_1(t)] \le 1 \Leftrightarrow \max_t [A_1(t) + R_1(t)] \le 1 \Rightarrow \max_t [A_1(t)] \le 1 + \max_t [R_1(t)]. \quad (B9)$$

Further, the following notations will be used:

$$\hat{A}_n = \max_t A_n(t), \quad \hat{R}_n = \max_t R_n(t).$$

As for $\hat{R}_n$, the following estimation may be obtained:

$$
\begin{align}
\hat{R}_n = \max_t [R_n(t)] &= \max_t \left\{ \frac{(-1)^n 4p}{\pi} \sum_{l=1}^{\infty} \frac{\xi_l^n \cos[(2l+1)(t-\varphi)]}{\sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}} \right\} \nonumber \\
&\le \frac{4p}{\pi} \sum_{l=1}^{\infty} \frac{\xi_l^n}{\sqrt{(2l+1)^4 - 4\gamma(2l+1)^2}} \le \frac{P}{\pi} \left(\frac{\gamma}{4}\right)^n \zeta(2n+2) = p\gamma^n\mu_n, \quad \mu_n = \frac{\zeta(2n+2)}{4^n\pi}. \tag{B10}
\end{align}
$$
To derive estimation (B11), we used expression (B8) and the fact that $\gamma \in [0,0.25)$. $\zeta(x)$ denotes the Riemann zeta function. For values of $\hat{A}_n$ one can obtain exact expressions from definition (B7) with the help of conditions (13) as follows:

$$
\begin{align*}
A_n(t) &= a \cos t + \frac{(-1)^n 4 p \xi_0^n}{\pi \sqrt{1-4\gamma}} \cos(t-\varphi) = a \cos \varphi \cos(t-\varphi) - a \sin \varphi \sin(t-\varphi) \\
&\quad + \frac{(-1)^n 4 p \xi_0^n}{\pi \sqrt{1-4\gamma}} \cos(t-\varphi) = \left[1 - p\chi(\gamma) + \frac{(-1)^n 4 p \xi_0^n}{\pi}\right] \cos(t-\varphi) + pq \sin(t-\varphi) \\
&\Rightarrow \hat{A}_n = \max_t [A_n(t)] = \sqrt{(1-pr_n)^2 + p^2 q^2}, \\
r_n &= \chi(\gamma) - \frac{(-1)^n 4 \xi_0^n}{\pi \sqrt{1-4\gamma}}. \tag{B11}
\end{align*}
$$

From inequality (B9) it follows that

$$
\hat{A}_1^2 \le (1 + \hat{R}_1)^2 \le 1 + 2p\gamma\mu_1 + p^2\gamma^2\mu_1^2 \Leftrightarrow p \le \frac{2(r_1 + \gamma\mu_1)}{r_1^2 + q^2 - \gamma^2\mu_1^2}. \quad (B12)
$$

If one considers inequality (B2) in a similar manner, it is possible to derive the following estimation:

$$
\max_t [u_n(t)] = \max_t [A_n(t) + R_n(t)] \le \hat{A}_n + \hat{R}_n \le 1 \Leftrightarrow p \le \frac{2(r_n - \gamma^n \mu_n)}{r_n^2 + q^2 - \gamma^{2n} \mu_n^2}. \quad (B13)
$$

In order to prove the latter inequality [and, consequently, inequality (B2)] it is enough to demonstrate that

$$
\frac{2(r_1 + \gamma \mu_1)}{r_1^2 + q^2 - \gamma^2 \mu_1^2} \le \frac{2(r_n - \gamma^n \mu_n)}{r_n^2 + q^2 - \gamma^{2n} \mu_n^2}. \quad (B14)
$$

Then, (B13) would follow from (B12). Inequality (B13) can be rewritten as

$$
(r_1 - r_n)(r_1 r_n - q^2) - (\gamma^n \mu_n r_1^2 + \gamma \mu_1 r_n^2) + \gamma^{2n} \mu_1^2 r_1 - \gamma^2 \mu_1^2 r_n + \gamma^{n+2} \mu_1^2 \mu_n + \gamma^{2n+1} \mu_n^2 \mu_1 \ge 0. \quad (B15)
$$

From (B11) one can derive the following estimations:

$$
\frac{4}{\pi} \left[ \frac{1 - (-1)^n \xi_0^n}{\sqrt{1 - 4\gamma}} + \frac{\pi^2}{8} - 1 \right] \le r_n < \frac{4}{\pi} \left[ \frac{1 - (-1)^n \xi_0^n}{\sqrt{1 - 4\gamma}} + \frac{\pi^2}{24} \right]. \quad (B16)
$$

Using these estimations, one can easily demonstrate the
correctness of inequality (B15) for all $\gamma \in [0,0.25)$, $q \in$
$[0,1]$. It should be mentioned that the equality in (B10)

and (B15) holds only for $\gamma = 0$ and in (B13), and,
consequently, in (B2)—only for a degenerate case with
$a = 1, \gamma = 0$.

[1] A. A. Ovchinnikov, Zh. Exp. Theor. Phys. **57**, 263 (1969) [Sov. Phys. JETP-USSR **30**, 147 (1970)]; S. Flach and C. R. Willis, Phys. Rep. **295**, 181 (1998); S. Aubry, Physica D **216**, 1 (2006); R. S. MacKay and S. Aubry, Nonlinearity **7**, 1623 (1994).

[2] S. Flach and A. Gorbach, Phys. Rep. **467**, 1 (2008).

[3] J. Cuevas, L. Q. English, P. G. Kevrekidis, and M. Anderson, Phys. Rev. Lett. **102**, 224101 (2009).

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

[5] N. Lazarides, M. Eleftheriou, and G. P. Tsironis, Phys. Rev. Lett. **97**, 157406 (2006).

[6] L. Q. English, F. Palmero, P. Candiani, J. Cuevas, R. Carretero-Gonzalez, P. G. Kevrekidis, and A. J. Sievers, Phys. Rev. Lett. **108**, 084101 (2012).

[7] M. Sato, B. E. Hubbard, and A. J. Sievers, Rev. Mod. Phys. **78**, 137 (2006); M. Sato, S. Imai, N. Fujita, S. Nishimura, Y. Takao, Y. Sada, B. E. Hubbard, B. Ilic, and A. J. Sievers, Phys. Rev. Lett. **107**, 234101 (2011).

[8] U. T. Schwarz, L. Q. English, and A. J. Sievers, Phys. Rev. Lett.

83, 223 (1999).

[9] A. J. Sievers and S. Takeno, Phys. Rev. Lett. **61**, 970 (1988); S. Takeno and A. J. Sievers, Solid State Commun. **67**, 1023 (1988); E. Kenig, B. A. Malomed, M. C. Cross, and R. Lifshitz, Phys. Rev. E **80**, 046202 (2009).

[10] A. H. Nayfeh and D. T. Mook, *Nonlinear Oscillations* (Wiley, New York, 1996).

[11] J. L. Marin, F. Falo, P. J. Martinez, and L. M. Floria, Phys. Rev. E **63**, 066603 (2001).

[12] M. J. Ablowitz and J. F. Ladik, J. Math. Phys. **17**, 1011 (1976).

[13] A. A. Ovchinnikov and S. Flach, Phys. Rev. Lett. **83**, 248 (1999).

[14] I. V. Barashenkov and Yu. S. Smirnov, Phys. Rev. E. **54**, 5707 (1996); I. V. Barashenkov, M. M. Bogdan, and V. I. Korobov, Europhys. Lett. **15**, 113 (1991).

[15] N. B. Tufillaro, T. Abbott, and J. Reilly, *An Experimental Approach to Nonlinear Dynamics and Chaos* (Addison-Wesley, Reading, MA, 2002).

[16] G. Casati, B. V. Chirikov, F. M. Izrailev, and J. Ford, in *Stochastic Behavior in Classical and Quantum Hamiltonian*
*Systems*, edited by G. Casati and J. Ford, Lecture Notes in Physics Vol. 93 (Springer, New York, 1979).

[17] S. W. Shaw and P. J. Holmes, J. Appl. Mech. **50**, 849 (1983); J. B. Page, Phys. Rev. B **41**, 7835 (1990); K. W. Sandusky, J. B. Page, and K. E. Schmidt, *ibid.* **46**, 6161 (1992); G. Casati, J. Ford, F. Vivaldi, and W. M. Visscher, Phys. Rev. Lett. **52**, 1861 (1984); O. V. Gendelman and A. V. Savin, *ibid.* **92**, 074301 (2004); A. V. Savin, G. P. Tsironis, and A. V. Zolotaryuk, *ibid.* **88**, 154301 (2002); V. Rom-Kedar and D. Turaev, Chaos **22**, 026102 (2012).

[18] O. V. Gendelman and L. I. Manevitch, Phys. Rev. E **78**, 026609 (2008).

[19] L. D. Landau and E. M. Lifshitz, *Mechanics* (Pergamon Press, Oxford, 1960).

[20] R. Seydel, *Practical Bifurcation and Stability Analysis* (Springer, New York, 2010).

[21] M. H. Fredriksson and A. B. Nordmark, Proc. R. Soc. London, Ser. A **456**, 315 (2000); M. di Bernardo, C. J. Budd, A. R. Champneys, and P. Kowalczyk, *Piecewise-smooth Dynamical Systems. Theory and Applications* (Springer, New York, 2008).

[22] S. W. Shaw and P. Holmes, Phys. Rev. Lett. **51**, 623 (1983).

[23] O. V. Gendelman, Chaos Solitons Fractals **28**, 522 (2006).

[24] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical Functions* (Dover, New York, 1965).