Breathers in Hamiltonian ${\cal PT}$-symmetric chains of coupled pendula under a resonant periodic force

We derive a Hamiltonian version of the ${\cal PT}$-symmetric discrete nonlinear Schr\"{o}dinger equation that describes synchronized dynamics of coupled pendula driven by a periodic movement of their common strings. In the limit of weak coupling between the pendula, we classify the existence and spectral stability of breathers (time-periodic solutions localized in the lattice) supported near one pair of coupled pendula. Orbital stability or instability of breathers is proved in a subset of the existence region.


Introduction
Synchronization is a dynamical process where two or more interacting oscillatory systems end up with identical movement. In 1665 Huygens experimented with maritime pendulum clocks and discovered the anti-phase synchronization of two pendulums clocks mounted on the common frame [15]. Since then, synchronization has become a basic concept in nonlinear and complex systems [10]. Such systems include, but are not limited to, musical instruments, electric power systems, and lasers. There are numerous applications in mechanical [22] and electrical [21] engineering. New applications are found in mathematical biology such as synchronous variation of cell nuclei, firing of neurons, forms of cooperative behavior of animals and humans [28].
Recently, Huygen's experiment has been widely discussed and several experimental devices were built [8,20,23]. It was shown that two real mechanical clocks when mounted to a horizontally moving beam can synchronize both in-phase and anti-phase [14]. In all these experiments synchronization was achieved due to energy transfer via the oscillating beam, supporting Huygen's intuition [8].
One of the rapidly developing areas in between physics and mathematics, is the topic of PT -symmetry, which has started as a way to characterize non-Hermitian Hamiltonians in quantum mechanics [6]. The key idea is that a linear Schrödinger operator with a complexvalued potential, which is symmetric with respect to combined parity (P) and time-reversal(T ) transformations, may have a real spectrum up to a certain critical value of the complex potential amplitude. In nonlinear systems, this distinctive feature may lead to existence of breathers (time-periodic solutions localized in space) as continuous families of their energy parameter.
The most basic configuration having PT symmetry is a dimer, which represents a system of two coupled oscillators, one of which has damping losses and the other one gains some energy from external sources. This configuration was studied in numerous laboratory experiments involving electric circuits [31], superconductivity [29], optics [2,30] and microwave cavities [9].
In the context of synchronization of coupled oscillators in a PT -symmetric system, one of the recent experiments was performed by Bender et al. [7]. These authors considered a PTsymmetric Hamiltonian system describing the motion of two coupled pendula whose bases were connected by a horizontal rope which moves periodically in resonance with the pendula. The phase transition phenomenon, which is typical for PT -symmetric systems, happens when some of the real eigenvalues of the complex-valued Hamiltonian become complex. The latter regime is said to have broken PT symmetry.
On the analytical side, dimer equations were found to be completely integrable [1,27]. Integrability of dimers is obtained by using Stokes variables and it is lost when more coupled nonlinear oscillators are added into a PT -symmetric system. Nevertheless, it was understood recently [3,4] that there is a remarkable class of PT -symmetric dimers with cross-gradient Hamiltonian structure, where the real-valued Hamiltonians exist both in finite and infinite chains of coupled nonlinear oscillators. Analysis of synchronization in the infinite chains of coupled oscillators in such class of models is a subject of this work.
In the rest of this section, we describe how this paper is organized. We also describe the main findings obtained in this work.
Section 2 introduces the main model of coupled pendula driven by a resonant periodic movement of their common strings. See Figure 1 for a schematic picture. By using an asymptotic multi-scale method, the oscillatory dynamics of coupled pendula is reduced to a PT -symmetric discrete nonlinear Schrödinger (dNLS) equation with gains and losses. This equation generalizes the dimer equation derived in [4,7]. Section 3 describes symmetries and conserved quantities for the PT -symmetric dNLS equation. In particular, we show that the cross-gradient Hamiltonian structure obtained in [3,4] naturally appears in the asymptotic reduction of the original Hamiltonian structure of Newton's equations of motion for the coupled pendula. Section 4 is devoted to characterization of breathers, which are time-periodic solutions localized in the chain. We show that depending on parameters of the model (such as detuning frequency, coupling constant, driving force amplitude), there are three possible types of breather solutions. For the first type, breathers of small and large amplitudes are connected to each other and do not extend to symmetric synchronized oscillations of coupled pendula. In the second and third types, large-amplitude and small-amplitude breathers are connected to the symmetric synchronized oscillations but are not connected to each other. See Figure 2 with branches (a), (b), and (c), where the symmetric synchronized oscillations correspond to the value E = 0 and the breather amplitude is given by parameter A.
Section 5 contains a routine analysis of linear stability of the zero equilibrium, where the phase transition threshold to the broken PT -symmetry phase is explicitly found. Breathers are only studied for the parameters where the zero equilibrium is linearly stable. Section 6 explores the Hamiltonian structure of the PT -symmetric dNLS equation and characterizes breathers obtained in Section 4 from their energetic point of view. We show that the breathers for large value of parameter E appear to be saddle points of the Hamiltonian function between continuous spectra of positive and negative energy, similar to the standing waves in the Dirac models. Therefore, it is not clear from the energetic point of view if such breathers are linearly or nonlinearly stable. On the other hand, we show that the breathers for smaller values of parameter E appear to be saddle points of the Hamiltonian function with a negative continuous spectrum and finitely many (either three or one) positive eigenvalues. Section 7 is devoted to analysis of spectral and orbital stability of breathers. For spectral stability, we use the limit of small coupling constant between the oscillators (the same limit is also used in Sections 4 and 6) and characterize eigenvalues of the linearized operator. The main analytical results are also confirmed numerically. See Figure 3 for the three types of breathers. Depending on the location of the continuous spectral bands relative to the location of the isolated eigenvalues, we are able to prove nonlinear orbital stability of breathers for branches (b) and (c). We are also able to characterize instabilities of these types of breathers that emerge depending on parameters of the model. Regarding branch (a), nonlinear stability analysis is not available by using the energy method. Our follow-up work [12] develops a new method of analysis to prove the long-time stability of breathers for branch (a).
The summary of our findings is given in the concluding Section 8, where the main results are shown in the form of Table 1.

Model
A simple yet universal model widely used to study coupled nonlinear oscillators is the Frenkel-Kontorova (FK) model [19]. It describes a chain of classical particles coupled to their neighbors and subjected to a periodic on-site potential. In the continuum approximation, the FK model reduces to the sine-Gordon equation, which is exactly integrable. The FK model is known to describe a rich variety of important nonlinear phenomena, which find applications in solid-state physics and nonlinear science [11].
We consider here a two-array system of coupled pendula, where each pendulum is connected to the nearest neighbors by linear couplings. Figure 1 shows schematically that each array of pendula is connected in the longitudinal direction by the torsional springs, whereas each pair of pendula is connected in the transverse direction by a common string. Newton's equations of motion are given by ẍ n + sin(x n ) = C (x n+1 − 2x n + x n−1 ) + Dy n , y n + sin(y n ) = C (y n+1 − 2y n + y n−1 ) + Dx n , n ∈ Z, t ∈ R, where (x n , y n ) correspond to the angles of two arrays of pendula, dots denote derivatives of angles with respect to time t, and the positive parameters C and D describe couplings between the two arrays in the longitudinal and transverse directions, respectively. The type of coupling between the two pendula with the angles x n and y n is referred to as the direct coupling between nonlinear oscillators (see Section 8.2 in [28]). We consider oscillatory dynamics of coupled pendula under the following assumptions.
(A1) The coupling parameters C and D are small. Therefore, we can introduce a small parameter µ such that both C and D are proportional to µ 2 .
(A2) A resonant periodic force is applied to the common strings for each pair of coupled pendula. Therefore, D is considered to be proportional to cos(2ωt), where ω is selected near the unit frequency of linear pendula indicating the 1 : 2 parametric resonance between the force and the pendula.
y n x n x n−1 x n+1 y n+1 y n−1 Figure 1: A schematic picture for the chain of coupled pendula connected by torsional springs, where each pair is hung on a common string.
Mathematically, we impose the following representation for parameters C and D(t): where γ, ǫ, Ω are µ-independent parameters, whereas µ is the formal small parameter to characterize the two assumptions (A1) and (A2).
In the formal limit µ → 0, the pendula are uncoupled, and their small-amplitude oscillations can be studied with the asymptotic multi-scale expansion where (A n , B n ) are amplitudes for nearly harmonic oscillations and (X n , Y n ) are remainder terms. In a similar context of single-array coupled nonlinear oscillators, it is shown in [24] how the asymptotic expansions like (3) can be justified. From the conditions that the remainder terms (X n , Y n ) remain bounded as the system evolves, the amplitudes (A n , B n ) are shown to satisfy the discrete nonlinear Schrödinger (dNLS) equations, which bring together all the phenomena affecting the nearly harmonic oscillations (such as cubic nonlinear terms, the detuning frequency, the coupling between the oscillators, and the amplitude of the parametric driving force). A similar derivation for a single pair of coupled pendula is reported in [4].
The system (5) takes the form of coupled parametrically forced dNLS equations (see [5] for references on parametrically forced NLS equations).
There exists an invariant reduction of system (5) given by to the scalar parametrically forced dNLS equation. The reduction (6) corresponds to the symmetric synchronized oscillations of coupled pendula of the model (1) with x n = y n , n ∈ Z.
In what follows, we consider a more general class of synchronized oscillations of coupled pendula. It turns out that the model (5) can be cast to the form of the parity-time reversal (PT ) dNLS equations [4]. Using the variables the system of coupled dNLS equations (5) is rewritten in the equivalent form 2iu n = ǫ (v n+1 − 2v n + v n−1 ) + Ωv n + iγu n + 2 2|u n | 2 + |v n | 2 v n + u 2 nv n , 2iv n = ǫ (u n+1 − 2u n + u n−1 ) + Ωu n − iγv n + 2 |u n | 2 + 2|v n | 2 u n +ū n v 2 n , which is the starting point for our analytical and numerical work. The invariant reduction (6) for system (5) becomes Without loss of generality, one can scale parameters Ω, ǫ, and γ by a factor of two in order to eliminate the numerical factors in the system (9). Also in the context of hard nonlinear oscillators (e.g. in the framework of the φ 4 theory), the cubic nonlinearity may have the opposite sign compared to the one in the system (9). However, given the applied context of the system of coupled pendula, we will stick to the specific form (9) in further analysis.

Symmetries and conserved quantities
The system of coupled dNLS equations (9) is referred to as the PT -symmetric dNLS equation because the solutions remain invariant with respect to the action of the parity P and timereversal T operators given by The parameter γ introduces the gain-loss coefficient in each pair of coupled oscillators due to the resonant periodic force. In the absence of all other effects, the γ-term of the first equation of system (9) induces the exponential growth of amplitude u n , whereas the γ-term of the second equation induces the exponential decay of amplitude v n , if γ > 0. The system (9) truncated at a single site (say n = 0) is called the PT -symmetric dimer. In the work of Barashenkov et al. [4], it was shown that all PT -symmetric dimers with physically relevant cubic nonlinearities represent Hamiltonian systems in appropriately introduced canonical variables. However, the PT -symmetric dNLS equation on a lattice does not typically have a Hamiltonian form if γ = 0.
Nevertheless, the particular nonlinear functions arising in the system (9) correspond to the PT -symmetric dimers with a cross-gradient Hamiltonian structure [4], where variables (u n ,v n ) are canonically conjugate. As a result, the system (9) on the chain Z has additional conserved quantities. This fact looked like a mystery in the recent works [3,4].
Here we clarify the mystery in the context of the derivation of the PT -symmetric dNLS equation (9) from the original system (1). Indeed, the system (1) of classical Newton particles has a standard Hamiltonian structure with the energy function Since the periodic movement of common strings for each pair of pendula result in the timeperiodic coefficient D(t), the energy H x,y (t) is a periodic function of time t. In addition, no other conserved quantities such as momenta exist typically in lattice differential systems such as the system (1) due to broken continuous translational symmetry.
After the system (1) is reduced to the coupled dNLS equations (5) with the asymptotic expansion (3), we can write the evolution problem (5) in the Hamiltonian form with the standard straight-gradient symplectic structure where the time variable t stands now for the slow time µ 2 t and the energy function is The energy function H A,B is conserved in the time evolution of the Hamiltonian system (13). In addition, there exists another conserved quantity which is related to the gauge symmetry (A, B) → (Ae iα , Be iα ) with α ∈ R for solutions to the system (5). When the transformation of variables (8) is used, the PT -symmetric dNLS equation (9) is cast to the Hamiltonian form with the cross-gradient symplectic structure where the energy function is The gauge-related function is written in the form The functions H u,v and Q u,v are conserved in the time evolution of the system (9). These functions follow from (14) and (15) after the transformation (8) is used. Thus, the crossgradient Hamiltonian structure of the PT -symmetric dNLS equation (9) is inherited from the Hamiltonian structure of the coupled oscillator model (1).

Breathers (time-periodic solutions)
We characterize the existence of breathers supported by the PT -symmetric dNLS equation (9). In particular, breather solutions are continued for small values of coupling constant ǫ from solutions of the dimer equation arising at a single site, say the central site at n = 0. We shall work in a sequence space ℓ 2 (Z) of square integrable complex-valued sequences. Time-periodic solutions to the PT -symmetric dNLS equation (9) are given in the form [18,26]: where the parameter E is considered to be real, the factor 1/2 is introduced for convenience, and the sequence ( (19) is considered to be PTsymmetric with respect to the operators in (11) if V =Ū . The reduction (10) for symmetric synchronized oscillations is satisfied if The time-periodic breathers (19) with E = 0 generalize the class of symmetric synchronized oscillations (20). The time-independent sequence (U, V ) ∈ ℓ 2 (Z) can be found from the stationary PTsymmetric dNLS equation: The PT -symmetric breathers with V =Ū satisfy the following scalar difference equation Note that the reduction (20) is compatible with equation (22) in the sense that if E = 0 and U n = R n e −iπ/4 , then R satisfies a real-valued difference equation. Let us set ǫ = 0 for now and consider solutions to the dimer equation at the central site n = 0: The parameters γ and Ω are considered to be fixed, and the breather parameter E is thought to parameterize continuous branches of solutions to the nonlinear algebraic equation (23). The solution branches depicted on Figure 2 are given in the following lemma.  (a) Ω > |γ| -two symmetric unbounded branches exist for ±E > E 0 , Proof. Substituting the decomposition U 0 = Ae iθ with A > 0 and θ ∈ [−π, π) into the algebraic equation (23), we obtain Excluding θ by using the fundamental trigonometric identity, we obtain the explicit parametrization of the solutions to the algebraic equation (23) by the amplitude parameter A: The zero-amplitude limit If |Ω| < |γ| , the solution branches (if they exist) are bounded away from the zero solution. Now we analyze the three cases of parameters γ and Ω formulated in the lemma.
Remark 1. The reduction (20) corresponds to the choice: Every solution of Lemma 1 can be extended to a breather on the chain Z which satisfies the spatial symmetry condition in addition to the PT symmetry: In order to prove the existence of the symmetric breather solution to the difference equation (22), we use the following implicit function theorem.
Implicit Function Theorem (Theorem 4.E in [32]). Let X, Y and Z be Banach spaces and let F (x, y) : There are r > 0 and σ > 0 such that for each y with With two applications of the implicit function theorem, we prove the following main result of this section.
There exists ǫ 0 > 0 sufficiently small and C 0 > 0 such that for every ǫ ∈ (−ǫ 0 , ǫ 0 ), there exists a unique solution U ∈ l 2 (Z) to the difference equation (22) satisfying the symmetry (30) and the bound where A and θ are defined in Lemma 1. Moreover, the solution U is smooth in ǫ.
Proof. In the first application of the implicit function theorem, we consider the following system of algebraic equations where U 0 ∈ C is given, in addition to parameters γ, Ω, and E. Let x = {U n } n∈N , X = ℓ 2 (N), y = ǫ, Y = R, and Z = ℓ 2 (N). Then, we have F (0, 0) = 0 and the Jacobian operator D x F (0, 0) is given by identical copies of the matrix with the eigenvalues λ ± := E ± Ω 2 − γ 2 . By the assumption of the lemma, λ ± = 0, so that the Jacobian operator D x F (0, 0) is one-to-one and onto. By the implicit function theorem, for every U 0 ∈ C and every ǫ = 0 sufficiently small, there exists a unique small solution U ∈ ℓ 2 (N) of the system (32) such that where the positive constant C 1 is independent from ǫ and U 0 .
Thanks to the symmetry of the difference equation (22), we find that U −n = U n , n ∈ N satisfy the same system (32) with −n ∈ N, with the same unique solution.
In the second application of the implicit function theorem, we consider the following algebraic equation where U 1 ∈ C depends on U 0 , γ, Ω, and E, satisfies the bound (33), and is uniquely defined by the previous result. Let x = U 0 , X = C, y = ǫ, Y = R, and Z = C. Then, we have F (Ae iθ , 0) = 0, where A and θ are defined in Lemma 1. The Jacobian operator D x F (Ae iθ , 0) is given by the matrix We show in Lemma 2 below that the matrix given by (35) is invertible under the conditions γ = 0 and Ω = −2|γ|. By the implicit function theorem, for every ǫ = 0 sufficiently small, there exists a unique solution U 0 ∈ C to the algebraic equation (34) near Ae iθ such that where the positive constant C 2 is independent from ǫ. The bound (31) holds thanks to the bounds (33) and (36). Since both equations (32) and (34) are smooth in ǫ, the solution U is smooth in ǫ.
In the following result, we show that the matrix given by (35) is invertible for every branch of Lemma 1 with an exception of a single point E = 0 on branch (c) for Ω = −2|γ|. Proof. The matrix given by (35) has zero eigenvalue if and only if its determinant is zero, which happens at Eliminating E 2 by using parametrization (25) and simplifying the algebraic equation for nonzero A 2 , we reduce it to the form 2(Ω + 4A 2 ) 3 = Ωγ 2 .
We now check if this constraint can be satisfied for the three branches of Lemma 1.
(a) If Ω > |γ|, the constraint (37) is not satisfied because the left-hand side exceeds the right-hand side Ωγ 2 .
Remark 2. In the asymptotic limit E 2 = 64A 4 + O(A 2 ) as A → ∞, see Lemma 1, the matrix (35) is expanded asymptotically as with the two eigenvalues λ 1 = E and λ 2 = −2E. Thus, the matrix given by (38) is invertible for every branch extending to sufficiently large values of E.

Stability of zero equilibrium
Here we discuss the linear stability of the zero equilibrium in the PT -symmetric dNLS equation (9). The following proposition yields a simple result.
Proof. Truncating the PT -symmetric dNLS equation (9) at the linear terms and using the Fourier transform u n (t) = 1 2π we obtain the linear homogeneous system The determinant ofD(k) is zero if and only if ω(k) is found from the quadratic equation For any |γ| < γ 0 , where γ 0 is given by (39), the two branches ±ω(k) found from the quadratic equation (41) are real-valued and non-degenerate for every k ∈ [−π, π]. Therefore, the zero equilibrium is linearly stable. On the other hand, for any |γ| > γ 0 , the values of ω(k) are purely imaginary either near k = ±π if Ω > 0 or near k = 0 if Ω < 0. Therefore, the zero equilibrium is linearly unstable.
If ǫ = 0, the zero equilibrium is only linearly stable for |γ| < |Ω|. Since the localized breathers cannot be stable when the zero background is unstable, we shall study stability of breathers only for the case when |γ| < |Ω|, that is, in the regime of unbroken PT -symmetry.

Variational characterization of breathers
It follows from Theorem 1 that each interior point on the solution branches shown on Figure  2 generates a fundamental breather of the PT -symmetric dNLS equation (9). We shall now characterize these breathers as relative equilibria of the energy function.
Thanks to the cross-gradient symplectic structure (16), the stationary PT -symmetric dNLS equation (21) can be written in the gradient form Keeping in mind the additional conserved quantity Q u,v given by (18), we conclude that the stationary solution (U, V ) is a critical point of the combined energy function given by If we want to apply the Lyapunov method in order to study nonlinear stability of stationary solutions in Hamiltonian systems, we shall investigate convexity of the second variation of the combined energy functional H E at (U, V ). Using the expansion u = U + u, v = V + v and introducing extended variables Φ and φ with the blocks we can expand the smooth function H E up to the quadratic terms in φ: where H ′′ E is the self-adjoint (Hessian) operator defined on ℓ 2 (Z) and the scalar product was used in the following form: x, y l 2 = k∈Z x kȳk .
Using (17) and (18), the Hessian operator can be computed explicitly as follows where blocks of L at each lattice node n ∈ Z are given by and ∆ is the discrete Laplacian operator applied to blocks of φ at each lattice node n ∈ Z: In the expression for L n , we have used the PT -symmetry condition V =Ū for the given stationary solution (U, V ).
We study convexity of the combined energy functional H E at (U, V ). Since the zero equilibrium is linearly stable only for |γ| < |Ω| (if ǫ = 0), we only consider breathers of Theorem 1 for |γ| < |Ω|. In order to study eigenvalues of H ′′ E for small values of ǫ, we use the following perturbation theory.
Perturbation Theory for Linear Operators (Theorem VII.1.7 in [17]). Let T (ǫ) be a family of bounded operators from Banach space X to itself, which depends analytically on the small parameter ǫ. If the spectrum of T (0) is separated into two parts, the subspaces of X corresponding to the separated parts also depend analytically on ǫ. In particular, the spectrum of T (ǫ) is separated into two parts for any ǫ = 0 sufficiently small.
With an application of the perturbation theory for linear operators, we prove the following main result of this section.
• If |E| < E 0 and Ω < −|γ|, the spectrum of H ′′ E in ℓ 2 (Z) includes an infinite-dimensional negative part and either three or one simple positive eigenvalues for branches (b) and (c) of Lemma 1 respectively.
Proof. If ǫ = 0, the breather solution of Theorem 1 is given by U n = 0 for every n = 0 and U 0 = Ae iθ , where A and θ are defined by Lemma 1. In this case, the linear operator H ′′ E = L decouples into 4-by-4 blocks for each lattice node n ∈ Z.
For n = 0, the 4-by-4 block of the linear operator L is given by Using relations (24) and (25), as well as symbolic computations with MAPLE, we found that the 4-by-4 matrix block L 0 admits a simple zero eigenvalue and three nonzero eigenvalues µ 1 , µ 2 , and µ 3 given by For each branch of Lemma 1 with γ = 0 and E = ±E 0 , we have 4A 2 + Ω = 0, so that µ 1 = 0. Furthermore, either µ 2 = 0 or µ 3 = 0 if and only if Expanding this equation for nonzero A yields constraint (37). With the exception of a single point E = 0 at Ω = −2|γ|, we showed in Lemma 2 that the constraint (37) does not hold for any of the branches of Lemma 1. Therefore, µ 2 = 0 and µ 3 = 0 along each branch of Lemma 1 and the signs of µ 1 , µ 2 , and µ 3 for each branch of Lemma 1 can be obtained in the limit A → ∞ for branches (a) and (b) or A → 0 for branch (c). By means of these asymptotic computations as A → ∞ or A → 0, we obtain the following results for the three branches shown on Figure 2: (c) µ 1 < 0, µ 2 > 0, and µ 3 < 0.
For n ∈ Z\{0}, the 4-by-4 block of the linear operator L is given by Each block has two double eigenvalues µ + and µ − given by Since there are infinitely many nodes with n = 0, the points µ + and µ − have infinite multiplicity in the spectrum of the linear operator L. Furthermore, we can sort up the signs of µ + and µ − for each point on the three branches shown on Figure 2: (1),(3) If |E| > E 0 := Ω 2 − γ 2 , then µ + > 0 and µ − < 0.
By using the perturbation theory for linear operators, we argue as follows: • Since H ′′ E is Hermitian on ℓ 2 (Z), its spectrum is a subset of the real line for every ǫ = 0.
• The zero eigenvalue persists with respect to ǫ = 0 at zero because the eigenvector (47) belongs to the kernel of H ′′ E due to the gauge invariance for every ǫ = 0.
• The other eigenvalues of L are isolated away from zero. The spectrum of H ′′ E is continuous with respect to ǫ and includes infinite-dimensional parts near points µ + and µ − for small ǫ > 0 (which may include continuous spectrum and isolated eigenvalues) as well as simple eigenvalues near µ 1,2,3 (if µ 1,2,3 are different from µ ± ).
The statement of the theorem follows from the perturbation theory and the count of signs of µ 1,2,3 and µ ± above. Remark 4. In the asymptotic limit E 2 = 64A 4 +O(A 2 ) as A → ∞, we can sort out eigenvalues of H ′′ E asymptotically as: where the remainder terms are O(1) as |E| → ∞. The values µ 1 , µ 3 , and µ + are close to each other as |E| → ∞.

Spectral and orbital stability of breathers
Spectral stability of breathers can be studied for small values of coupling constant ǫ by using the perturbation theory [26]. First, we linearize the PT -symmetric dNLS equation (9) at the breather (19) by using the expansion where (u, v) is a small perturbation satisfying the linearized equations The spectral stability problem arises from the linearized equations (52) after the separation of variables: where φ := (ϕ, ψ, χ, ν) is the eigenvector corresponding to the spectral parameter λ. Note that (ϕ, ψ) and (χ, ν) are no longer complex conjugate to each other if λ has a nonzero imaginary part. The spectral problem can be written in the explicit form where we have used the condition V =Ū for the PT -symmetric breathers. Recalling definition of the Hessian operator H ′′ E in (46), we can rewrite the spectral problem (53) in the Hamiltonian form: where S is a symmetric matrix with the blocks at each lattice node n ∈ Z given by We note the Hamiltonian symmetry of the eigenvalues of the spectral problem (54).

Proposition 2.
Eigenvalues of the spectral problem (54) occur either as real or imaginary pairs or as quadruplets in the complex plane.
If Ω < −|γ| and |E| < E 0 := Ω 2 − γ 2 (points 2 and 4 shown on Figure 2), Theorem 2 implies that the self-adjoint operator H ′′ E in ℓ 2 (Z) is negative-definite with the exception of either three (point 2) or one (point 4) simple positive eigenvalues. In this case, we can apply the following Hamilton-Krein index theorem in order to characterize the spectrum of SH ′′ E .
Hamilton-Krein Index Theorem (Theorem 3.3 in [16]). Let L be a self-adjoint operator in ℓ 2 with finitely many negative eigenvalues n(L), a simple zero eigenvalue with eigenfunction v 0 , and the rest of its spectrum is bounded from below by a positive number. Let J be a bounded invertible skew-symmetric operator in ℓ 2 . Let k r be a number of positive real eigenvalues of JL, k c be a number of quadruplets {±λ, ±λ} that are neither in R nor in iR, and k − i be a number of purely imaginary pairs of eigenvalues of JL whose invariant subspaces lie in the negative subspace of L. Let D = L −1 J −1 v 0 , J −1 v 0 ℓ 2 be finite and nonzero. Then, Lemma 3. Fix γ = 0, Ω < −|γ|, and 0 < |E| < E 0 , where E 0 := Ω 2 − γ 2 > 0. For every ǫ > 0 sufficiently small, K HAM = 2 for branch (b) of Lemma 1 and K HAM = 0 for branch (c) of Lemma 1 with Ω < −2 √ 2|γ|. For branch (c) with Ω ∈ (−2 √ 2|γ|, −|γ|), there exists a value E s ∈ (0, E 0 ) such that K HAM = 1 for 0 < |E| < E s and K HAM = 0 for E s < |E| < E 0 .
Proof. If γ = 0, Ω < −|γ|, |E| < E 0 , and ǫ > 0 is sufficiently small, Theorem 2 implies that the spectrum of H ′′ E in ℓ 2 (Z) has finitely many positive eigenvalues and a simple zero eigenvalue with eigenvector σΦ. Therefore, the Hamilton-Krein index theorem is applied in ℓ 2 (Z) for L = −H ′′ E , J = iS, and v 0 = σΦ. We shall verify that where σΦ is given by (47) and ∂ E Φ denotes derivative of Φ with respect to parameter E. The first equation H ′′ E (σΦ) = 0 follows by Theorem 2. By differentiating equations (21) in E, we obtain H ′′ E (∂ E Φ) = SσΦ for every E, for which the solution Φ is differentiable in E. For ǫ = 0, the limiting solution of Lemma 1 is differentiable in E for every E = 0 and E = ±E 0 . Due to smoothness of the continuation in ǫ by Theorem 2, this property holds for every ǫ > 0 sufficiently small. By using (57) with S −1 = S, we obtain where we have used the definition of Q u,v in (18). We compute the slope condition at ǫ = 0: where relations (24) and (25) have been used. For branch (b) of Lemma 1 with Ω < −|γ|, we have dA 2 /dE 2 > 0 and |Ω| < 4A 2 , so that dQ u,v /dE > 0. By continuity, dQ u,v /dE remains strictly positive for small ǫ > 0. Thus, D < 0 and K HAM = 2 by the Hamilton-Krein index theorem.
If K HAM = 0 and D = 0, orbital stability of a critical point of H E in space ℓ 2 (Z) can be proved from the Hamilton-Krein theorem (see [16] and references therein). Orbital stability of breathers is understood in the following sense. Definition 1. We say that the breather solution (19) is orbitally stable in ℓ 2 (Z) if for every ν > 0 sufficiently small, there exists δ > 0 such that if ψ(0) ∈ ℓ 2 (Z) satisfies ψ(0) − Φ ℓ 2 ≤ δ, then the unique global solution ψ(t) ∈ ℓ 2 (Z), t ∈ R to the PT -symmetric dNLS equation (9) satisfies the bound inf The definition of instability of breathers is given by negating Definition 1. The following result gives orbital stability or instability for branch (c) shown on Figure 2.
Orbital stability of breathers for branches (a) and (b) of Lemma 1 does not follow from the standard theory because K HAM = ∞ for |E| > E 0 and K HAM = 2 > 0 for branch (b) with |E| < E 0 . Nevertheless, by using smallness of parameter ǫ and the construction of the breather (U, V ) in Theorem 1, spectral stability of breathers can be considered directly. Spectral stability and instability of breathers is understood in the following sense.
Definition 2. We say that the breather solution (19) is spectrally stable if λ ∈ iR for every bounded solution of the spectral problem (54). On the other hand, if the spectral problem (54) admits an eigenvalue λ / ∈ iR with an eigenvector in ℓ 2 (R), we say that the breather solution (19) is spectrally unstable.
The following theorem gives spectral stability of breathers for branches (a) and (b) shown on Figure 2.
where the eigenvector σΦ is given by (47) and the generalized eigenvector ∂ E Φ denotes derivative of Φ with respect to parameter E. For every E such that the following non-degeneracy condition is satisfied, the breather (U, V ) is spectrally stable.
Proof. If ǫ = 0, the breather solution of Theorem 1 is given by U n = 0 for every n = 0 and U 0 = Ae iθ , where A and θ are defined by Lemma 1. In this case, the spectral problem (53) decouples into 4-by-4 blocks for each lattice node n ∈ Z. Recall that H ′′ E = L at ǫ = 0. For n = 0, eigenvalues λ are determined by the 4-by-4 matrix block −iSL 0 . Using relations (24) and (25), as well as symbolic computations with MAPLE, we found that the 4-by-4 matrix block −iSL 0 has a double zero eigenvalue and a pair of simple eigenvalues at λ = ±λ 0 , where For n ∈ Z\{0}, eigenvalues λ are determined by the 4-by-4 matrix block −iSL n , where L n is given by (50). If |γ| < |Ω|, E = 0, and E = ±E 0 , where E 0 := Ω 2 − γ 2 , each block has four simple eigenvalues ±λ + and ±λ − , where so that λ ± ∈ iR. Since there are infinitely many nodes with n = 0, the four eigenvalues are semi-simple and have infinite multiplicity. If ǫ > 0 is sufficiently small, we use perturbation theory for linear operators from Section 6.
• The double zero eigenvalue persists with respect to ǫ = 0 at zero because of the gauge invariance of the breather (U, V ) (with respect to rotation of the complex phase). Indeed, H ′′ E (σΦ) = 0 follows from the result of Theorem 2. The generalized eigenvector is defined by equation SH ′′ E Ψ = σΦ, which is equivalent to equation H ′′ E Ψ = (V,V , U,Ū ) T . Differentiating equations (21) in E, we obtain Ψ = ∂ E Φ. Since dim[Ker(H ′′ E )] = 1 and the second generalized eigenvectorΨ ∈ ℓ 2 (Z) exists as a solution of equation SH ′′ EΨ = ∂ E Φ if and only if dQ u,v /dE = 0. It follows from the explicit computation (59) that if ǫ = 0, then dQ u,v /dE = 0 for every E along branches (a) and (b) of Lemma 1. By continuity, dQ u,v /dE = 0 for small ǫ > 0. Therefore, the zero eigenvalue of the operator −iSH ′′ E is exactly double for small ǫ > 0.
• Using the same computation (59), it is clear that λ 0 ∈ iR for every E along branches (a) and (b) of Lemma 1. Assume that λ 0 = ±λ + and λ 0 = ±λ − , which is expressed by the non-degeneracy condition (62). Then, the pair ±λ 0 is isolated from the rest of the spectrum of the operator −iSH ′′ E at ǫ = 0. Since the eigenvalues λ = ±λ 0 are simple and purely imaginary, they persist on the imaginary axis for ǫ = 0 because they cannot leave the imaginary axis by the Hamiltonian symmetry of Proposition 2.
• If |γ| < |Ω|, E = 0, and E = ±E 0 , the semi-simple eigenvalues ±λ + and ±λ − of infinite multiplicity are nonzero and located at the imaginary axis at different points for ǫ = 0. They persist on the imaginary axis for ǫ = 0 according to the following perturbation argument. First, for the central site n = 0, the spectral problem (53) can be written in the following abstract form where L 0 (ǫ) denotes a continuation of L 0 in ǫ. Thanks to the non-degeneracy condition (62) as well as the condition λ ± = 0, the matrix SL 0 − iλ ± I is invertible. By continuity, the matrix SL 0 (ǫ) − iλI is invertible for every ǫ and λ near ǫ = 0 and λ = λ ± . Therefore, there is a unique φ 0 given by which satisfies |φ 0 | ≤ Cǫ(|φ 1 | + |φ −1 |) near ǫ = 0 and λ = λ ± , where C is a positive ǫand λ-independent constant. Next, for either n ∈ N or −n ∈ N, the spectral problem (53) can be represented in the form where L n (ǫ) denotes a continuation of L n given by (50) in ǫ, whereas the operator ∆ is applied with zero end-point condition at n = 0. We have L n (ǫ) = L n + O(ǫ 2 ) and ǫSφ 0 = O(ǫ 2 ) near ǫ = 0 and λ = λ ± . Therefore, up to the first order of the perturbation theory, the spectral parameter λ near λ ± is defined from the truncated eigenvalue problem which is solved with the discrete Fourier transform (40). In order to satisfy the Dirichlet end-point condition at n = 0, the sine-Fourier transform must be used, which does not affect the characteristic equation for the purely continuous spectrum of the spectral problem (66). By means of routine computations, we obtain the characteristic equation in the form, see also equation (41), where k ∈ [−π, π] is the parameter of the discrete Fourier transform (40). Solving the characteristic equation (67), we obtain four branches of the continuous spectrum where the two sign choices are independent from each other. If |Ω| > |γ| is fixed and ǫ > 0 is small, the four branches of the continuous spectrum are located on the imaginary axis near the points ±λ + and ±λ − given by (64).
In addition to the continuous spectrum given by (67), there may exist isolated eigenvalues near ±λ + and ±λ − , which are found from the second-order perturbation theory [25]. Under the condition E = 0 and E = ±E 0 , these eigenvalues are purely imaginary. Therefore, the infinite-dimensional part of the spectrum of the operator −iSH ′′ E persists on the imaginary axis for ǫ = 0 near the points ±λ + and ±λ − of infinite algebraic multiplicity.
The statement of the lemma follows from the perturbation theory and the fact that all isolated eigenvalues and the continuous spectrum of −iSH ′′ E are purely imaginary.

Remark 9.
Observe in the proof of Theorem 4 that λ ± / ∈ iR if |Ω| < |γ|. In this case, branch (b) of Lemma 1 is spectrally unstable. This instability corresponds to the instability of the zero equilibrium for |Ω| < |γ|, in agreement with the result of Proposition 1.
Before presenting numerical approximations of eigenvalues of the spectral problem (54), we compute the Krein signature of wave continuum. This helps to interpret instabilities and resonances that arise when isolated eigenvalues ±λ 0 cross the continuous bands near points ±λ + and ±λ − . The Krein signature of simple isolated eigenvalues is defined as follows.
Definition 3. Let φ ∈ ℓ 2 (Z) be an eigenvector of the spectral problem (54) for an isolated simple eigenvalue λ 0 ∈ iR. Then, the energy quadratic form H ′′ E φ, φ ℓ 2 is nonzero and its sign is called the Krein signature of the eigenvalue λ 0 .
Definition 3 is used to simplify the presentation. Similarly, one can define the Krein signature of isolated multiple eigenvalues and the Krein signature of the continuous spectral bands in the spectral problem (54) [16]. The following lemma characterizes Krein signatures of the spectral points arising in the proof of Theorem 4. (c) all subspaces of −iSH ′′ E in ℓ 2 (Z) near ±λ + , ±λ − , and ±λ 0 (if λ 0 ∈ iR) have negative Krein signature.
The signs of all eigenvalues are nonzero and continuous with respect to parameter ǫ. Therefore, the count above extends to the case of small nonzero ǫ.
The spectrum of −iSH ′′ E is shown at Figure 3. Panels (a), (b) and (c) correspond to branches shown at Figure 2. (a) We can see on panel (a) of Figure 3 that λ 0 , λ ± do not intersect for every E > E 0 and are located within fixed distance O(1), as |E| → ∞. Note that the upper-most λ 0 and λ + have positive Krein signature, whereas the lowest λ − has negative Krein signature, as is given by Lemma 4.
(b) We observe on panel (b) of Figure 3 that λ + intersects λ 0 , creating a small bubble of instability in the spectrum. The insert shows that the bubble shrinks as ǫ → 0, in agreement with Theorem 4. There is also an intersection between λ − and λ 0 , which does not create instability. These results are explained by the Krein signature computations in Lemma 4. Instability is induced by opposite Krein signatures between λ + and λ 0 , whereas crossing of λ − and λ 0 with the same Krein signatures is safe of instabilities.
Note that for small E, the isolated eigenvalue λ 0 is located above both the spectral bands near λ + and λ − . The gap in the numerical data near E = E 0 indicates failure to continue the breather solution numerically in ǫ, in agreement with the proof of Theorem 1.
(c) We observe from panel (c) of Figure 3 that λ 0 and −λ − intersect but do not create instabilities, since all parts of the spectrum have the same signature, as is given by Lemma 4. In fact, the branch is both spectrally and orbitally stable as long as λ 0 ∈ iR, in agreement with Theorem 3. On the other hand, there is E s ∈ (0, E 0 ), if Ω ∈ (−2 √ 2|γ|, −|γ|), such that λ 0 ∈ R for E ∈ (0, E s ), which indicates instability of branch (c), again, in agreement with Theorem 3.
As we see on panel (b) of Figure 3, λ 0 intersects λ + for some E = E * > E 0 . In the remainder of this section, we study whether this crossing point is always located on the right of E 0 . In fact, the answer to this question is negative. We shall prove for branch (b) that the intersection of λ 0 with either λ + or −λ − occur either for E * > E 0 or for E * < E 0 , depending on parameters γ and Ω. where Moreover, if Ω < −5|γ|, there exists a resonance λ 0 = −λ − at E = E * with E * ∈ (0, E 0 ).

Summary
We have reduced Newton's equation of motion for coupled pendula shown on Figure 1 under a resonant periodic force to the PT -symmetric dNLS equation (9). We have shown that this system is Hamiltonian with conserved energy (17) and an additional constant of motion (18). We have studied breather solutions of this model, which generalize symmetric synchronized oscillations of coupled pendula that arise if E = 0. We showed existence of three branches of breathers shown on Figure 2. We also investigated their spectral stability analytically and numerically. The spectral information on each branch of solutions is shown on Figure 3. For branch (c), we were also able to prove orbital stability and instability from the energy method. The technical results of this paper are summarized in Table 1 and described as follows. Table 1: A summary of results on breather solutions for small ǫ. Here, IB is a narrow instability bubble seen on panel (b) of Figure 3.
|E| > E 0 |E| < E 0 Parameter intervals Ω > |γ| Ω < −|γ| Ω < −|γ| Ω < −|γ| For branch (a), we found that it is disconnected from the symmetric synchronized oscillations at E = 0. Along this branch, breathers of small amplitudes A are connected to breathers of large amplitudes A. Every point on the branch corresponds to the saddle point of the energy function between two wave continua of positive and negative energies. Every breather along the branch is spectrally stable and is free of resonance between isolated eigenvalues and continuous spectrum. In the follow-up work [12], we will prove long-time orbital stability of breathers along this branch.
For branch (b), we found that the large-amplitude breathers as E → ∞ are connected to the symmetric synchronized oscillations at E = 0, which have the smallest (but nonzero) amplitude A = A + . Breathers along the branch are spectrally stable except for a narrow instability bubble, where the isolated eigenvalue λ 0 is in resonance with the continuous spectrum. The instability bubble can occur either for E > E 0 , where the breather is a saddle point of the energy function between two wave continua of opposite energies or for E < E 0 , where the breather is a saddle point between the two negative-definite wave continua and directions of positive energy. When the isolated eigenvalue of positive energy λ 0 is above the continuous spectrum near λ + and ±λ − , orbital stability of breathers can be proved by using the technique in [13], which was developed for the dNLS equation.
Finally, for branch (c), we found that the small-amplitude breathers at E → E 0 are connected to the symmetric synchronized oscillations at E = 0, which have the largest amplitude A = A − . Breathers are either spectrally stable near E = E 0 or unstable near E = 0, depending on the detuning frequency Ω and the amplitude of the periodic resonant force γ. When breathers are spectrally stable, they are also orbitally stable for infinitely long times.