Solitons in PT-symmetric ladders of optical waveguides

We consider a PT-symmetric ladder-shaped optical array consisting of a chain of waveguides with gain coupled to a parallel chain of waveguides with loss. All waveguides have the focusing Kerr nonlinearity. The array supports two co-existing solitons, an in-phase and an antiphase one, and each of these can be centred either on a lattice site or midway between two neighbouring sites. We show that both bond-centred (i.e. intersite) solitons are unstable regardless of their amplitudes and parameters of the chain. The site-centred in-phase soliton is stable when its amplitude lies below a threshold that depends on the coupling and gain-loss coefficient. The threshold is lowest when the gain-to-gain and loss-to-loss coupling constant in each chain is close to the interchain gain-to-loss coupling coefficient. The antiphase soliton in the strongly-coupled chain or in a chain close to the PT-symmetry breaking point, is stable when its amplitude lies above a critical value and unstable otherwise. The instability growth rate of solitons with small amplitude is exponentially small in this parameter regime; hence the small-amplitude solitons, though unstable, have exponentially long lifetimes. On the other hand, the antiphase soliton in the weakly or moderately coupled chain and away from the PT-symmetry breaking point, is unstable when its amplitude falls in one or two finite bands. All amplitudes outside those bands are stable.


I. INTRODUCTION
In soliton-bearing optical systems, weak dissipative losses are typically compensated by the application of a resonant pump. The damping and driving terms break the scaling-, phase-and Galilean invariances of the underlying nonlinear Schrödinger equation. As a result, the amplitude, phase, and velocity of the damped-driven solitons are found to be fixed by the driver's strength and dissipation coefficient [1].
Dissipative solitons in systems with competing linear and nonlinear gain and loss -systems modelled by the Ginsburg-Landau equations -are equally inflexible. The balance of nonlinearity and diffraction or dispersion singles out a one-parameter family of solitons, while the competition of gain and loss fixes a particular member of the family [2]. This inflexibility can be a disadvantage in applications where one would like the parameters of the soliton to be determined just by the initial conditions.
FIG. 1. A PT -symmetric ladder of coupled optical waveguides with gain (red) and loss (blue). Each guide is coupled to its left and right nearest neighbours in the horizontal plane. We are occasionally referring to this coupling as "horizontal" in the text. In addition, each guide with gain is coupled to its nearest guide with loss -this will be referred to as the "vertical coupling".
Yet there are ways of supplying energy in a soliton-bearing system with loss, free from the above drawback. The pertinent class of approaches exploits the concept of the parity-time (PT -) symmetry, originally developed in the context of quantum mechanics [3]. One PT -symmetric recipe consists in the manufacturing of two identical replicas of the same system and arranging for the energy supply in one copy and energy drain in the other one. With a judicious choice of coupling, the combined system with gain and loss will contain, as its scalar reduction, the unperturbed nonlinear Schrödinger equation. As a result, this PT -symmetric system will support families of solitons with continuously variable amplitudes, phases, and velocities.
In this paper, we consider an infinite chain of PT -symmetric dimers coupled to their nearest neighbours. The coupling of the neighboring dimers is arranged in such a way that, apart from its "internal" counterpart with gain, each waveguide with loss interacts only with its two lossy neighbours (Fig 1). Likewise, each waveguide with gain is coupled only to its two neighbours with gain -besides its counterpart with loss with who they make up the PT -symmetric entity [40]: iu n + C∆u n + 2|u n | 2 u n + v n = iγu n , iv n + C∆v n + 2|v n | 2 v n + u n = −iγv n . (1.1) Here ∆ is the second difference operator acting on bi-infinite sequences ..., u −2 , u −1 , u 0 , u 1 , u 2 , ... according to the rule ∆u n = u n+1 − 2u n + u n−1 (n = 0, ±1, ±2, ....).
In the basis of infinite-component vectors u = (..., u −2 , u −1 , u 0 , u 1 , u 2 , ...) T , the operator ∆ is represented by a tridiagonal matrix. The u n is the nondimensional complex mode amplitude in the nth waveguide with gain, and v n is the amplitude in the nth waveguide with loss. The overdot indicates the derivative with respect to z, the propagation coordinate. The sign of the cubic terms corresponds to the self-focusing Kerr nonlinearity; the coefficients in front of these terms have been set to 2 by scaling the mode amplitudes. Finally, C > 0 is a dimer-to-dimer coupling constant (the "horizontal" coupling in the language of Fig 1), and γ ≥ 0 is the gain-loss coefficient.
The ladder-shaped chain (1.1) supports four families of soliton excitations [40]. These are distinguishable by the phase difference between u n and v n , and by the centring of the soliton relative to lattice sites. In the continuum limit, their stability has been classified in [45,47]. Ref [51] considered the anticontinuum limit, under an additional assumption of γ → 0. In the present paper we study the discrete solitons in the entire range of their amplitudes, for all values of the coupling C and all physically admissible values of γ.
We will demonstrate that the soliton stability problem can be decomposed into an eigenvalue problem for perturbations belonging to the same symmetric manifold as the stationary soliton itself, and perturbations that take the evolution out of the symmetry. This decomposition alone is sufficient to classify the stability of the bond-centred solitons -which happen to be all unstable against a perturbation in the symmetric manifold. In contrast, the instability of the site-centred solitons can only be nucleated by nonsymmetric perturbations. In the latter case, the advantage of our decomposition is that it reduces the underlying two-component eigenvalue problem to a much simpler one-component problem.
We will show that stability properties of the in-phase and antiphase site-centred soliton are completely determined just by two combinations of C, γ and the soliton's amplitude. The stability domain of the in-phase soliton admits a simple characterisation in terms of a single function relating these two self-similar combinations. As a result, for each C and γ we will be able to determine a critical value of the amplitude above which the soliton becomes unstable.
Stability properties of the antiphase soliton in a strongly-coupled chain or in a weakly or moderately coupled chain close to its PT -symmetry breaking point, follow an opposite pattern. Here, the soliton instability sets in below a critical amplitude. Reducing the inter-dimer coupling or turning down the gain-loss coefficient of the chain, changes the topography of this soliton's parameter space. Namely, the weakly or moderately coupled array away from the PT -symmetry breaking threshold displays one or two bands of unstable amplitudes, with all antiphase solitons outside those bands being stable.
The outline of the paper is as follows. Section II summarises stability properties of two solitons arising in the continuum limit of the array (1.1). In section III, we introduce four discrete solitons: the in-phase and antiphase site-centred soliton, and the in-and out-of-phase bond-centred one. The subsequent section IV lays out the general framework of our stability analysis and draws some general conclusions.
In section V we classify the stability of the in-phase solitons while section VI focusses on the antiphase variety. Some technical results on eigenvalues playing the central role in our analysis have been relegated to eight appendices. Finally, section VII summarises conclusions of this study and contrasts them with the earlier results of [40] and [51].

II. CONTINUUM LIMIT: SUMMARY
If the characteristic wavelength in the chain is much greater than C −1/2 , the PT -symmetric lattice (1.1) reduces to a system of two partial differential equations: This PT -symmetric system of coupled nonlinear Schrödinger equations was considered in [40,45,47,48]. Results of the previous analyses can be summarised as follows.
(1) For each γ < 1, there are two continuous families of coexisting solitons. Each family is characterised by its own value of the phase difference between the uand v-component. Representatives of the two families with equal amplitudes have unequal propagation constants [40]. .) The left panel shows the form of |un| and |vn| for the site-centred soliton (both for the in-phase and antiphase one). The right panel illustrates the magnitude profile of the bond-centred soliton -again, for the in-and out-of-phase one alike.
(2) The soliton stability and internal dynamics are determined by a single self-similar combination η = 2 1 − γ 2 /A 2 of the gain-loss coefficient γ and the soliton's amplitude, A [47].
In this paper, we will extend these considerations to localised solutions of the discrete system (1.1).

III. DISCRETE SOLITONS
The dispersion relations for the small-amplitude harmonic waves in the system (1.1) are In what follows, we assume that the trivial solution, u n = v n = 0, is stable: 0 ≤ γ ≤ 1.
For future convenience, we perform the change of variables u n = e iθ a n + b n , v n = e iθ b n − a n , where θ = arcsin γ. This transformation diagonalises the linear part of equations (1.1) [36,48]: where u n and v n denote the linear combinations (3.1). The system (1.1) has two invariant manifolds where the evolution is conservative [47]. Indeed, the equations (3.2) admit a reduction a n = 0, b n = φ n to the scalar discrete nonlinear Schrödinger equation (Note that under this reduction, u n = φ n and v n = e iθ φ n .) A different, independent, scalar reduction arises by letting b n = 0, a n = φ n . The equations (3.2) become iφ n + C∆φ n − ω 0 φ n + 2|φ n | 2 φ n = 0.
Here χ is a phase mismatch taking the initial condition out of the scalar reduction (3.4). In this plot the parameters of the lattice are C = 0.3 and γ = 0.93; the amplitude of the gaussian Ags = 2 and the mismatch χ = 0.3. Note small-amplitude oscillations of the soliton's amplitude -the non-oscillatory soliton would result from the gaussian with χ = 0.
In this paper, we are focussing on localised solutions of (3.3) and (3.4): |φ n | → 0 as |n| → ∞. The simplest solutions of this sort are separable; these are given by φ n = e iβz AR n .
Here A > 0 is an arbitrary amplitude and R n solves 1 h 2 ∆R n − R n + 2R 3 n = 0, (3.5) where h is the amplitude-to-coupling ratio: In what follows, we are referring to the corresponding solutions of (3.2) as solitons. Equation (3.5) determines the profile of the localised solution and h has the meaning of discretisation stepsize in this equation. In the case of the reduction (3.3), the propagation constant β is given by while the equation (3.4) requires We note that equations (3.3)-(3.4) admit no localised separable solutions with complex R n (except for the trivial situation where all R n share the same, n-independent, phase factor) [53]. Therefore, equation (3.5) is the most general equation for localised profiles.
Each localised solution of (3.5) gives rise to two families of separable solutions of the original system (1.1), with C = (A/h) 2 . The first family corresponds to the a n = 0 reduction of Eq. (3.2). When the gain-loss coefficient γ → 0, the u n and v n components of any solution in this family become equal: v n → u n . In contrast, the components of any separable solution corresponding to the b n = 0 reduction become equal in magnitude but opposite in sign: v n → −u n .
Because of this simple and noticeable difference, the localised separable solutions with a n = 0 will be called the in-phase solitons, and those with b n = 0 the antiphase ones. (By continuity, we will use this terminology even in the γ → 1 limit where the difference between the two families becomes negligible.) Note that the soliton solution of Eq.(3.2) with a n = 0 and amplitude A, has a larger propagation constant than the (b n = 0)-soliton of the same amplitude. For this reason, the in-phase separable solution was referred to as the high-frequency soliton in [47], and the one with b n = 0 as the low-frequency one. Equation (3.5) has a variety of localised solutions, in particular two fundamental, i.e., single-hump, nonlinear modes. One of these is centred on the site with n = 0 (the so-called Sievers-Takeno mode) and the other one midway between the sites with n = 0 and n = 1 (the Page mode) [54][55][56][57][58]. We will be referring to the former as the 'site-centred' mode, and the latter as 'bond-centred'. Both solutions exist for all 0 < h < ∞ [59]; both satisfy R n > 0 for all n [58]. We depict these solutions in Fig. 2.
Thus there are four fundamental solitons in the dimer array (1.1): the in-phase and antiphase site-centred soliton, and the in-phase and antiphase bond-centred one. The distribution of the waveguide mode magnitudes in the siteand bond-centred soliton is shown in Fig 2. The phases of u n and v n remain constant as n grows or decreases; hence there is no flow of energy along the array. On the other hand, the phase difference between u n and v n reflects a nonzero energy flux from the gaining to the losing waveguide of each dimer.
Small h correspond to the continuum limit (small soliton amplitudes or strongly coupled sites). As h → 0, the site-centred and bond-centred solutions satisfy R n → sech(hn) and R n → sech[h(n − 1 2 )], respectively. Large h pertain to the anti-continuum limit (large amplitudes or weak dimer-to dimer coupling).
Before proceeding to stability analysis, we note that discrete solitons and their breather-like perturbations emerge from a broad class of initial conditions. Fig 3 shows the formation of a soliton from a gaussian with generic amplitude and phase.

A. The in-phase soliton
Consider the in-phase soliton first. Letting a n = e iβz Re(p (1) n e µz ) + i Re(p (2) n e µz ) , b n = e iβz AR n + Re(q (1) n e µz ) + i Re(q (2) n e µz ) , where R n is a solution of (3.5), β is as in (3.7), and p L q n − ω 0 q n + 4γA 2 R 2 n σ 1 p n = µJ q n , (4.1b) where the operator acts on two-component complex sequences The J in the right-hand sides of (4.1a)-(4.1b) denotes a unimodular skew-symmetric matrix: and σ 1 in (4.1b) is the conventional Pauli matrix. The purpose of the stability analysis is to determine whether the eigenvalue problem (4.1) has eigenvalues µ with nonzero real part. If it has, the soliton is unstable; if it has not, the soliton is classified as linearly stable.
The continuous spectrum of µ's occupies two opposite pairs of bands on the imaginary axis: µ = ±iω 1,2 , where Outside the domain of intersection of the bands (4.2) and (4.3), the problem (4.1) may have discrete eigenvalues. To find the discrete eigenvalues, one does not have to solve the full system (4.1). Indeed, one class of eigenvectors is selected simply by letting p n = 0; this choice selects perturbations belonging to the scalar Schrödinger equation (3.3). The associated q n satisfies L q n − ω 0 q n = µJ q n , On the other hand, the site-centred soliton of the scalar Schrödinger equation (3.3) is stable for all A and C [55][56][57]. The operator J −1 (L − ω 0 I) in this case has only pure imaginary eigenvalues: two zeros, two pairs of opposite nonzero imaginary eigenvalues and two bands of continuous spectrum on the imaginary axis: This does not yet guarantee that the corresponding soliton of the vector equations (1.1) is stable; the stability of the latter will be determined by eigenvectors of (4.1) with p n = 0. These correspond to linear evolutions that do not satisfy the scalar reduction (3.3). The nontrivial sequences p n = 0 can be found as eigenvectors of the eigenvalue problem (4.1a). This is a significant simplification as the dimension of the matrix L + ω 0 I in (4.1a) is half the dimension of the full triangular matrix . Once the p n and the associated eigenvalue µ have been found, Eq.(4.1b) becomes a nonhomogeneous equation with the right-hand side determined by p n : Here N is a nonhermitian operator defined by The two-component eigenvector p n q n (4.7) of the full "triangular" eigenvalue problem (4.1) exists if the equation (4.6) has a solution with Whether the finite-norm solution q n exists or not, depends on whether the kernel space of the adjoint operator is empty, and if not -whether the right-hand side in (4.6) is orthogonal to the kernel. For the given eigenvalue µ in (4.1a), the operator N † has a nonempty kernel space only if there is a sequence y n such that that is, if −µ * coincides with one of the (pure imaginary) eigenvalues of the symplectic operator J −1 (L − ω 0 I), or falls in its continuous spectrum occupying two bands in (4.5). Therefore if µ has a nonzero real part -the situation that concerns us here -the kernel space of N † is empty and the two-component eigenvector (4.7) does exist. Thus the stability analysis of the in-phase soliton reduces to the solution of the eigenvalue problem (4.1a). All unstable eigenvalues of (4.1a) (if exist) are automatically eigenvalues of the full eigenvalue problem (4.1).

B. The antiphase soliton
In the case of the antiphase soliton, we let a n = e iβz AR n + Re(p (1) n e µz ) + i Re(p (2) n e µz ) , where, as before, R n is a solution of (3.5) and p (1,2) n , q (1,2) n , µ are complex. This time, however, β is given by (3.8) rather than (3.7). Linearising equations (3.2) in p (1,2) n and q (1,2) n gives L q n − ω 0 q n = µJ q n , (4.8a) L p n + ω 0 p n − 4γA 2 R 2 n σ 1 q n = µJ p n .  In the case of the site-centred soliton, the spectrum of (4.9) consists of two zeros, two pairs of opposite nonzero imaginary eigenvalues and two bands of continuous spectrum on the imaginary axis: µ = ±iω 1 , where Therefore, in order to detect instabilities of the site-centred soliton one has to turn to eigenvectors with q n = 0.
The eigenvalues pertaining to q n = 0 are determined by solving the "reduced" eigenvalue problem (4.8a). The p n -component of the eigenvector ( p n , q n ) is then recovered from the solution of the nonhomogeneous equation (4.8b). The question that concerns us here is whether there are any eigenvalues with nonzero real part; for those µ, the nonhomogeneous equation (4.8b) has a finite-norm solution (see the argument in the previous subsection). Therefore the full triangular matrix in (4.8) has an unstable eigenvalue as long as the half-size matrix in (4.8a) has one.

C. Symplectic eigenvalue problem
Since the general instability of the in-and out-of-phase bond-centred solitons has been established, we only need to consider the site-centred ones in what follows.
Starting with the in-phase soliton and forming (infinite-component) vectors p (1) and p (2) out of the elements of the sequences p (1) n and p (2) n , respectively, we rename p (1) = g and p (2) = f for notational convenience. The eigenvalue problem (4.1a) is then written in the matrix form: Here L + and L − are symmetric matrices with elements where ∆ mn = δ m,n+1 + δ m,n−1 − 2δ mn , h is as in (3.6), and R n is the site-centred single-hump solution of Eq. (3.5) (the Sievers-Takeno mode). In (4.10), we have denoted λ = µ/C and introduced a parameter σ = 2ω 0 /C > 0. Turning to the antiphase soliton and treating the sequences q (1) n and q (2) n as infinite-component vectors g and f , respectively, the eigenvalue problem (4.8a) is represented in the same form (4.10) -but with σ = −2ω 0 /C < 0. Therefore, equation (4.10) with σ varying from negative to positive values defines the stability problem both for the in-and the out-of-phase soliton. Here σ = 2ω 0 /C, in-phase soliton; −2ω 0 /C, antiphase soliton. (4.13) Note that the single-hump site-centred solution R n is uniquely specified by the value of the parameter h. This means that the stability properties of the soliton with amplitude A, in the system (1.1) with coupling C and gain-loss coefficient γ, are completely determined by two combinations of those three parameters: σ and h.
The combination σ in (4.13) is a structural parameter. It is fixed by the coupling constant C and the gain/loss coefficient γ -but does not depend on the soliton's amplitude. On the other hand, h in (3.6) is the amplitude of a particular solution under scrutiny normalised by the coupling: h = A/ √ C. Solitons with different amplitudes and in systems with different couplings and gain-loss coefficients have the same stability properties as long as they share the values of σ and h.
It is fitting to note that a self-similarity of this sort was previously encountered in the continuum limit of the system (1.1). There, the stability of solitons was determined by a single similarity parameter η = 2 1 − γ 2 /A 2 [47].
According to Eq.(3.5), the matrix L − has a zero eigenvalue, with the eigenvector R n . Keeping in mind that R n > 0 for all n, equation (3.5) can be cast in the form (4.14) Making use of the identity (4.14), the quadratic form y|L − |y with an arbitrary y = {y n } admits the following representation: This representation implies that the matrix L − does not have any negative eigenvalues.
On the other hand, one can readily check that the matrix L + has a negative eigenvalue. A simple upper bound for it is given by equation (A6) in Appendix A. The negative eigenvalue is single.
In fact, the operators (4.11) and (4.12) are not unfamiliar in the soliton stability literature. These infinitedimensional matrices occurred in the studies of the scalar discrete nonlinear Schrödinger equation -where some of their properties were established. In particular, it was shown that the operator L + corresponding to the sitecentred soliton has a single negative eigenvalue [55,56]. (Had there been two negative eigenvalues, the site-centred scalar soliton would have been unstable. In actual fact, stability is a firmly established property of that solution [55][56][57].) The operator in (4.10), is symplectic, that is, generates a hamiltonian flow. The spectrum of symplectic operators consists of pairs of pure-imaginary values, real pairs, and complex quadruplets. If λ is a real or pure imaginary point of spectrum, then −λ is another one; if a complex λ is in the spectrum, then so are λ * , −λ, and −λ * [60]. The continuous spectrum of λ lies on the imaginary axis and consists of two bands: λ = ±iω, where σ + h 2 ≤ ω ≤ σ + 4 + h 2 . When σ > 0 (the in-phase soliton), the bands are separated by a gap, whereas when σ < 0 (the antiphase soliton), the bands overlap. See Fig 4. When σ = 0, equation (4.10) coincides with the eigenvalue problem for the scalar nonlinear Schrödinger equation, with no dissipative or driving terms. It is important to keep in mind, however, that the problem (4.10) with σ = 0 does not result from the coupled waveguides with γ = 0. The value σ = 0 corresponds to γ = 1 rather than γ = 0. The significance of this value is that as γ is increased through 1, the in-phase and antiphase solitons merge and disappear. As for the "no-gain, no-loss" (γ = 0) pair of waveguides, it is represented by a nonzero σ (σ = ±2/C) and is not special as far as the eigenvalue problem (4.10) is concerned.
The scalar discrete nonlinear Schrödinger equation has two zero eigenvalues in its linearised spectrum, stemming from its U(1) phase invariance and the transformation to the linearly growing phase. Accordingly, the eigenvalue problem (4.10) with σ = 0 has two zero eigenvalues. As σ deviates from zero, the two eigenvalues move out of the origin.
Finally we note that the eigenvalue problem (4.10) with general σ occurred previously in a context unrelated to the PT -symmetry. (Specifically, it appeared in the analysis of the transverse stability of the one-dimensional soliton in the scalar two-dimensional lattice [64].) The authors of [64] have established some useful properties of eigenvalues in the anticontinuum limit h → ∞. We make contact with those results in what follows.

A. Stability domain
In this and the next section the word "soliton" refers to the site-centred soliton (the Sievers-Takeno mode). In the case of the in-phase soliton (σ > 0), the eigenvalue problem (4.10) is amenable to simple analysis. The smallest eigenvalue of the operator L − + σI equals σ; therefore, the matrix L − + σI is positive definite and admits an inverse. Hence the problem (4.10) can be written as a generalised eigenvalue problem for the vector g: The operator on the left in (5.1) is symmetric, and the one on the right is symmetric and positive definite. The lowest eigenvalue −λ 2 is given by the minimum of the Rayleigh quotient Let E 0 = E 0 (h) denote the (single) negative eigenvalue of L + . The minimum in (5.2) is nonnegative (hence the soliton is stable) if the smallest eigenvalue of the operator in the numerator is nonnegative: The eigenvalue E 0 , calculated numerically for a sequence of h values, is plotted in Fig 5 (a). It is clear from the figure that the function σ = −E 0 (h) grows, monotonically, from zero to infinity as h varies over the positive part of the axis. Denoting h = H 0 (σ) its inverse function, the stability condition for the in-phase soliton can be written as The function H 0 (σ), determined numerically, is plotted in Fig.5(c). Equation (A6) of the Appendix A provides a simple upper bound for this function: Using equations (B9) and (D4) in the corresponding Appendices, we obtain its small-and large-σ asymptotic behaviours: and unstable otherwise. Here A 0 is a function of the coupling and gain-loss coefficient: For each fixed value of C we will be referring to the function A 0 (γ, C) as the instability threshold. Equation (5.4) translates into a simple upper bound for the threshold: More accurate estimates for the boundary of the stability domain are available in two opposite limits. In the limit 1 − γ 2 /C → 0, the approximate expression for A 0 stems from equation (5.5): (5.10) (The first term in this expression reproduces the instability threshold for the solitons in the PT -symmetric coupler consisting of two parallel guiding planes [45,47].) In the weak coupling limit (C/ 1 − γ 2 → 0), an approximate expression for the instability threshold follows from equation (5.6): Plotting the curve (5.8) with a variety of C (Fig 6) one observes that the instability threshold is lowest when the coupling takes values near C = 0.5. It is impossible to identify a single "most unstable" value of C because small variations about C = 0.5 raise either the small-γ or large-γ part of the curve (5.8) while moving down the remaining part. Instead of a single value, there is a finite (yet narrow) interval of C encompassing 0.5. The instability thresholds associated with C in this interval form a thin belt underlying instability thresholds associated with C further away from C = 0.5. With this reservation in mind, we can refer to values close to C = 0.5 as "most unstable".
It is interesting to compare the "most unstable" values of the gain-to-gain and loss-to-loss coupling C to the value of the gain-to-loss coupling in the ladder (1.1). To make a fair comparison, we think of each guide with gain as being connected to its counterpart with loss by two bonds, left and right, with each coupling coefficient being equal to 1/2. (Geometrically, this corresponds to seeing each pair of gain and loss sites as a circular necklace consisting just of two elements.) In this symmetric picture, each site of the ladder has four bonds: two horizontal and two vertical. The maximum instability associated with C ≈ 0.5 implies that the soliton's stability domain is narrowest when the coupling coefficients of two horizontal and two vertical bonds are equal.
As C is raised from C ≈ 0.5, the instability threshold is lifted. (See Fig 6.) As C → ∞, it approaches the limit curve

B. Eigenvalue trajectories
Although equations (5.7)-(5.8) provide the complete solution to the stability problem of the in-phase soliton, they do not give any insight into the motion of the eigenvalues on the complex plane. Yet the knowledge of stable and unstable eigenvalues can address a number of pertinent questions, in particular classify the type of instability and identify modes of internal oscillation. The aim of this subsection is to consider the eigenvalues: variationally, asymptotically and numerically.
The variational principle (5.2) gives rise to simple upper and lower bounds of the lowest eigenvalue −λ 2 of the scalar eigenvalue problem (5.1) (see Appendix E). The two bounds allow one to follow the corresponding pair of opposite symplectic eigenvalues ±λ, as h grows from 0 to ∞.
For h smaller than H 0 , the set of symplectic eigenvalues may only consist of pure imaginary pairs λ = ±iω. We denote ±iω even the pair with the smallest value of ω 2 ; the notation will become clear later in this subsection. Using (E1), the modulus-squared ω 2 even is found to satisfy Recalling that the continuous spectrum of symplectic eigenvalues λ = ±iω in (4.10) occupies the band (σ + h 2 ) 2 ≤ ω 2 ≤ (σ + 4 + h 2 ) 2 , we observe that the right-hand side in (5.12) lies below the bottom edge of the continuous spectrum. Hence This implies that the lowest eigenvalue −λ 2 even cannot bifurcate from the lower edge of the continuum for any finite h. On the other hand, the inequality in (5.13) may become equality when h = 0; hence the bifurcation can occur at this point. Equation (E3) indicates that for small h, the quantity ω 2 even lies above σ(σ − 3h 2 ). Taken together with the inequality (5.12), this fact implies that "the innermost" pair of imaginary eigenvalues (i.e. the pair with the smallest value of ω 2 ) does indeed bifurcate from the continuum at the point h = 0. In fact, there are two pairs of pure imaginary eigenvalues, λ = ±iω even and λ = ±iω odd , bifurcating from the edges of the continuous spectrum as h is increased from zero. According to the asymptotic analysis of Appendix F, one pair has even and the other one odd eigenfunctions; hence the choice of notation. Equation (F6) gives Note that the lowest value of ω 2 is provided by the eigenvalue with the even eigenfunction. This explains the choice of notation in (5.12)-(5.13).
In the previous subsection we have shown that the pure imaginary pair ±iω even (the pair with lowest ω 2 ) converges at the origin and splits into the negative and positive real axis as h is increased through the critical value H 0 , When h grows to infinity, equations (E1) and (E4) give a corridor for the absolute value of the emerging real pair: The anti-continuum limit (Appendix H) provides a more specific result: These variational and asymptotic considerations are corroborated by numerical solutions. Figures 7 show the evolution of the imaginary and real parts of eigenvalues as h grows from 0 to large values, with σ being fixed. Clearly seen is a pair of imaginary eigenvalues which detaches from the continuous spectrum at h = 0, moves on to the real axis, and diverges to infinity at a linear rate, as h → ∞.

VI. ANTIPHASE SOLITON
To classify stability of the antiphase solitons, we consider the eigenvalue problem (4.10) with σ < 0. Unlike equation (4.10) with σ > 0, this eigenvalue problem is amenable to analytical treatment only in two regions on the (σ, h) plane. These are considered in subsections VI A and VI B. For other parameter values we have to rely upon the asymptotic and numerical eigenvalue analysis.
A. Instability against an odd mode When σ < 0, the matrix L − + σI has a negative eigenvalue (equal to σ). Hence, the lowest value of −λ 2 in the generalised eigenvalue problem (5.1) cannot be sought as a minimum of the Rayleigh quotient (5.2) over the entire 2 space of square-summable bi-infinite sequences. However, when σ > −h 2 , we can use the Rayleigh quotient to establish the existence of eigenvalues with odd (antisymmetric) eigenvectors g: g −n = −g n , n = 0, 1, 2, .... Indeed, the spectrum of the infinite matrix L − does not have eigenvalues other than E 0 = 0 and a continuous band h 2 ≤ E ≤ h 2 + 4. (We have verified this by computing numerical eigenvalues of a large finite truncation of L − over a representative sample of h values between 0 and ∞.) The zero eigenvalue corresponds to an even (symmetric) eigenvector z 0 = R. Therefore, the matrix L − + σI with σ > −h 2 , is positive definite on the subspace S of 2 consisting of odd sequences g. Consequently, the lowest eigenvalue of the generalised problem (5.1) associated with an odd eigenvector, is given by Let E 1 denote the second lowest eigenvalue of the matrix L + : The eigenvalue E 1 is positive [55,56] and satisfies 0 < E 1 (h) < h 2 . (See Fig.5(a)). The corresponding eigenvector y 1 is odd and renders the quadratic form g|L + + σI|g minimum in S. Equation (6.1) implies then that −λ 2 is positive respectively negative for E 1 + σ > 0 respectively E 1 + σ < 0.
In the region −h 2 < σ < −E 1 (h), the symplectic eigenvalue problem (4.10) has a pair of real eigenvalues ±λ signifying the instability of the antiphase soliton against an odd mode. As σ approaches −E 1 (h) from below, the two eigenvalues collide at the origin and move onto the imaginary axis. The odd-mode instability is replaced with a mode of internal oscillation.
In terms of the soliton's amplitude, A, the above odd-mode instability region is given by where and h = H 1 (|σ|) is the function inverse to |σ| = E 1 (h). This function is shown in Fig 5(c). The small-and large-|σ| behaviour of H 1 (|σ|) follow from the corresponding asymptotics of the eigenvalue E 1 (h) (see the Appendix A): H 1 = C ln 1/2 (|σ| −1 ) as |σ| → 0, Here C is a positive constant and e.s.t. denotes a positive term that is smaller than |σ| −n with any n > 0. The exponential approach of H 1 to |σ| 1/2 as |σ| → ∞, is clearly visible in Fig. 5(c). The region (6.2) is displayed in Fig.8 for three values of coupling C: C = 1, 4, and 5. It is important to emphasise here that the region (6.2) constitutes only a part of the full odd-mode instability range (namely, the part where this instability can be established by analytical means). While the inequality A < A 1 (C, γ) does demarcate the top Let |σ| > 4. If h 2 < |σ| − 4, or, equivalently, both matrices P + and P − are positive definite. Hence the variational principle (5.2) is valid and the lowest eigenvalue −λ 2 is positive: The upshot of this elementary consideration is that if the chain parameters satisfy the antiphase soliton has a range of stable amplitudes given by (6.6).
It is important to emphasise that the region described by the inequalities (6.8) and (6.6) constitutes only a part of the full stability domain of the antiphase soliton; see Fig 13(a). The significance of this region lies in the fact that stability is guaranteed by an accurate analytical argument here. Unlike the numerical eigenvalue analysis, this argument would have not missed even a weak instability -if there had been any.

C. Surprising stability of large-amplitude solitons
The two regions (6.2) and (6.8),(6.6) exhaust all parameter values where the stability of the antiphase solitons can be classified using exact analytical criteria. In the complementary part of the parameter space, one has to resort to numerical and asymptotic techniques.
Results of our numerical study of the eigenvalue problem (4.10) with negative σ are summarised in Figs 9 and 10. A striking difference between these two figures on the one hand -and stability properties of the corresponding continuous system [47] on the other -is in the discrete soliton stabilisation as its amplitude becomes large enough. (We remind the reader that once the coupling constant C has been fixed, the amplitude can be represented by the quantity h = A/ √ C.) Indeed, all we see in Figs 9 and 10 with a sufficiently large h is a pair of opposite pure imaginary eigenvalues. These eigenvalues are given by the power expansion (H5) with coefficients as in (H7): As h is decreased from large values, each of the two imaginary eigenvalues seems to collide with the corresponding branch of the continuous spectrum and acquire a real part. A scrutiny of the collision area reveals that the imaginary eigenvalue approaching the continuum collides not with the edge of the continuous-spectrum band itself but with another pure imaginary eigenvalue that bifurcates from the edge just before the collision. The distance in h between the point of bifurcation and point of collision is so short that this second eigenvalue is hardly discernible for large |σ| (Fig 9). However it is clearly visible when |σ| is chosen smaller than 1 (Fig 10).
The collision is nothing but the Hamiltonian Hopf bifurcation. It results in a quadruplet of complex eigenvalues ±λ, ±λ * . In what follows, the bifurcation value of h is denoted by H 3 (|σ|). Recalling that h = A/ √ C and σ = −2ω 0 /C, Left column: σ = −1; right column: σ = −0.1. Two light-green curves trace a complex quadruplet with even eigenvectors while the dark-green arc marks complex eigenvalues with odd eigenvectors. As the imaginary parts of the "even" quadruplet emerge from the continuous spectrum, the quadruplet dissociates into two pairs of imaginary eigenvalues described by the black and red curves. The "odd" quadruplet also dissociates as h is increased; the products include two pairs of real eigenvalues shown in light-blue and deep-blue, respectively. The "light-blue" eigenvalues converge on the origin and disappear in the gapless continuous spectrum. The "deep-blue" pair moves away from the origin but then reverses and crosses to the imaginary axis as well. By the time the "deep-blue" doublet reaches the origin, a gap will have opened in the continuous spectrum. The doublet reaches the boundaries of the gap and immerses in the continuum. The pink oval in (a) represents a pair of real eigenvalues with even eigenvectors. this gives the upper boundary of the soliton instability domain: where the function H 3 (|σ|) admits a simple asymptote [64]:

The asymptotic behaviour for the upper boundary of the soliton instability domain is then
D. Small amplitude solitons: stable or long-lived The limit h → 0, with small or finite |σ| (σ < 0), is also amenable to asymptotic analysis. See Appendix G and F, respectively. As h grows from zero, either two pairs of pure imaginary or two quadruplets of complex eigenvalues bifurcate from the endpoints of the continuous spectrum (more specifically, from the Im λ = ±σ endpoints). Whether the bifurcating eigenvalues have real parts or not, their imaginary parts are given by equations (F6) and (G3): Numerically, the bifurcating imaginary eigenvalues are clearly visible when |σ| > 2. See e.g. Fig 9(b) for σ = −5 and Fig 9(d) for σ = −3. When |σ| < 2, by contrast, the point Im λ = σ -the bottom of the upper band of the continuous spectrum -is embedded in the lower band, which covers the interval −σ − 4 ≤ Im λ ≤ −σ. (Compare the panel (b) of Fig 4 with panels (c) and (d).) In this case, the bifurcating eigenvalues with small or zero real parts are concealed in the continuum bands and are much harder to track.
The asymptotic analysis of the problem (4.10) with a negative σ of order one, provides no clue on whether the eigenvalues bifurcating from the edge of the continuum have real parts or not. The only conclusion one can draw from it is that if an eigenvalue has a nonzero real part, this real part will have to lie beyond all orders of h n . (The proof of this fact is in Appendix F.) The asymptotic argument in Appendix G is of more use here as it guarantees the existence of an (exponentially small) real part when |σ| is small enough. Also useful is the variational principle of section VI B; this principle ensures that all eigenvalues of the operator (4.15) with small h are pure imaginary if |σ| > 4. Taken together, these two considerations suggest that there is a critical value σ c (|σ c | ≤ 4) such that the bifurcating eigenvalues acquire an exponentially small real part when |σ| < |σ c | but remain pure imaginary when |σ| ≥ |σ c |.
Our numerics support this conjecture and indicate that σ c is close to −1.96. For all |σ| ≥ 1.96 that we have examined, the two pairs of bifurcating eigenvalues were found to be pure imaginary. They only develop nonzero real parts once their imaginary parts collide with the continuous spectrum bands at some finite h. (Note the square-root dependence of Re λ on h in the vicinity of the collision point in Fig 9 (a) and (c).) When |σ| < 1.96, by contrast, eigenvalues corresponding to small h exhibit nonzero real parts. In this case, the shape of the Re λ(h) curve is indicative of an exponential decay as h → 0. See Fig 10 (a,c) and note the difference in the behaviour of the curve from the square-root law near the collision point in Fig 9 (a,c).
The proximity of the critical value σ c to −2 sheds some light on the source of the weak oscillatory instability. Unlike the case |σ| > 2, the imaginary parts of the bifurcating eigenvalues of the operator (4.15) with |σ| ≤ 2 are embedded in the continuous spectrum right from the moment they were born, that is, as early as at h = 0. Consequently, the emergence of a nonzero real part is due to the resonance between the frequency corresponding to the imaginary part of the bifurcating eigenvalue, and frequencies of the continuous spectrum.

E. Weakly coupled chain away from PT -symmetry breaking: disjoint stability domain
When the structural parameter |σ| is greater than 4 (that is, when C < 1/2 and 0 ≤ γ ≤ √ 1 − 4C 2 ), the complex quadruplet associated with the continuous-band crossing is the only set of unstable eigenvalues that arise as h varies from 0 to ∞. See Fig 9 (a,b). In this case, as h is decreased from H 3 (the point of collision with the "inner" edge of the continuum), the imaginary part of the eigenvalue travels through the continuous band and crosses through its "outer" edge. At the point where the imaginary part emerges from the continuum, the real part of the eigenvalue appears to vanish. We denote this point by h = H 2 (|σ|) and the corresponding value of the amplitude by A 2 : The function H 2 admits a simple asymptotic representation [64]: This gives an asymptotic expression for the lower boundary of the instability region: As with the inner-edge crossing, the numerical scrutiny confirms that this process is more complex than it seems. It involves two pairs of complex eigenvalues which converge on two opposite points on the imaginary axis and dissociate into two pure imaginary pairs. After that, one pair of opposite imaginary eigenvalues immerses in the continuum whereas the other pair moves away from it. The interval of h where two imaginary pairs coexist is so tiny that the complex eigenvalues appear to be stripped of their real parts as soon as the imaginary parts have emerged out of the continuum.
While considering diagrams pertaining to large |σ|, a natural question concerns the odd-mode instability of section VI A. The domain of the odd-mode instability includes (though is not limited to) the interval |σ| 1/2 < h < H 1 (|σ|). Why is this instability not manifest in the λ(h) dependencies with large |σ| (Fig.9)?
The answer is that since the function H 1 (|σ|) approaches |σ| 1/2 exponentially fast as |σ| grows, the above interval of unstable h shrinks and becomes indiscernible (see Fig 5(c)). Furthermore, the numerical analysis indicates that the full odd-mode instability domain also shrinks as |σ| grows. Therefore, although the solitons with large negative σ are prone to the odd-mode instability, the interval of "unstable" h is minuscule.
The shrinking of the odd-mode instability band (6.2) as the coupling C is reduced, is clearly visible in Fig 8. In the weakly coupled chain of waveguides, the numerical detection of this band is only feasible if γ is close to 1 (so that the argument of H 1 in (6.3) is small enough).
Our conclusions can be summarised as follows. The stability domain in a chain with C < 1/2 and 0 ≤ γ ≤ √ 1 − 4C 2 consists of two disjoint parts. The soliton is unstable if its amplitude falls between A 2 and A 3 , and stable when A lies outside this band. Here A 2 is as in (6.12) and A 3 as in (6.10), with the corresponding asymptotic behaviours given by (6.13) and (6.11), respectively. (There is also a narrow band of the odd-mode instability immersed in the stability domain but it is of little significance due to its exponentially small width.) The disjoint stability domain on the (γ, A)-plane is illustrated in Fig 13(a). The dashed curve in this panel demarcates the region where the soliton's stability has been established using the variational principle (6.7). It is obvious that the analytical argument secures only a small portion of the full stability domain.

F. Numerical eigenvalue trajectories
While the complex quadruplet crossing the continuous spectrum band is the only set of unstable eigenvalues for |σ| > 4, chains with smaller values of |σ| feature several more instabilities. This subsection summarises results of our numerical study of the problem (4.10) with 0 < |σ| < 4. As before in this section, we keep σ < 0.
The complex quadruplet exhibited by the operator (4.15) with |σ| > 4 persists for |σ| < 4. In Figs 9 and 10 the trajectories of complex eigenvalues constituting the quadruplet in question, are delineated in light green. Similarly to the |σ| > 4 scenario, the quadruplet is born in a Hamiltonian Hopf bifurcation as h is decreased through the point H 3 (|σ|). A further decrease of h with a fixed |σ| ≥ 1.96 brings about an inverse Hopf bifurcation, at the point h = H 2 (|σ|), where the quadruplet dissociates into two opposite pairs of imaginary eigenvalues. When |σ| < 1.96, by contrast, the quadruplet persists over the entire interval 0 < h ≤ H 3 . (Its asymptotic behaviour as h → 0 was discussed in section VI D. ) We note that the eigenvectors associated with the "light-green" quadruplet, are symmetric (that is, even in n). The left-right symmetry of the eigenvectors distinguishes this quadruplet from another quadruplet of complex eigenvalues that is shown in dark green in Fig 10. The eigenvectors of the "dark-green" eigenvalues are antisymmetric (odd in n).
The left-right symmetry of the eigenvectors gives rise to a symmetric instability growth. Fig 11 illustrates two possible outcomes of the decay of the antiphase soliton unstable against the "light-green" complex quadruplet. In  Fig 11 (a) the growing oscillations culminate in the exponential blow-up of the soliton. In Fig 11 (b), the oscillatory instability transforms the soliton into a breather -a spatially localised structure with a steadily oscillating amplitude and width.
The pink loop in Fig 9(c) and 10(a) describes a pair of opposite real eigenvalues (also with even eigenvectors). These hail from the imaginary axis: as h is decreased, a pair of pure imaginary eigenvalues bifurcates from the continuous spectrum, converges at the origin and moves on to the real axis. As h is decreased further, the eigenvalues reach their maximum absolute value and then return to the origin.
The "pink loop" is present in all diagrams with |σ| between 4 and approximately 0.4. When |σ| < 1.96, these real eigenvalues coexist with the complex quadruplet and do not produce any further reduction of the soliton stability domain. The stability domain has a simple structure in this case: the soliton is stable when h > H 3 (|σ|) and unstable otherwise.
Unlike the oscillatory instability discussed above, the monotonic growth associated with the real mode cannot give rise to any breathing structures. A typical evolution is illustrated in Fig 12 (a). The unstable antiphase soliton sheds some power and transforms into a stable soliton of the same variety.
The inequality 1.96 ≤ |σ| < 4 selects two more regions on the (C, γ)-plane: Each of these two regions features two non-overlapping intervals of unstable amplitudes, A 4 < A < A 5 and A 2 < A < A 3 , where A n = √ CH n 2 1 − γ 2 /C , n = 2, 3, 4, 5. (6.16) Two corridors of instability are clearly visible in the left part of Fig 13(b) and middle section of Fig 13(a).
Another instability that has a sizeable domain in chains with |σ| < 2, is due to a pair of real eigenvalues with odd eigenvectors. The inequality |σ| 1/2 < h < H 1 (|σ|) of section VI A defines the upper part of the odd-mode instability domain; the lower part has to be determined numerically. The trajectory of the odd-mode eigenvalues on the (h, Re λ)-plane is described by the deep-blue curve in Fig 10 (a,c). Like the even real mode discussed in the previous paragraphs, this pair of real eigenvalues is born as an imaginary pair which bifurcates from the continuous spectrum. [See the deep-blue curves in Fig 10 (b,d).] As h is decreased to the value of H 1 (|σ|), the two imaginary eigenvalues collide and move to the real axis -in agreement with the variational argument of section VI A.
What the variational principle is powerless to determine though, is the trajectory of these eigenvalues on the complex plane as h is decreased below |σ| 1/2 . The numerical analysis, on the other hand, reveals the arrival of yet another pair of real eigenvalues from the imaginary axis via the origin (shown in light blue in Fig 10 (a,c)). Like the "dark-blue" eigenvalues, the newly born pair has odd eigenvectors. At some h < |σ| 1/2 , the "light-blue" and "deep-blue" eigenvalues collide, pairwise, and emerge into the complex plane.
The evolution of the resulting complex quadruplet depends on whether |σ| is greater or smaller than 2. When |σ| < 2, the quadruplet persists all the way to h = 0, with the real parts of the complex eigenvalues decaying to zero. (See the stiletto-shaped dark-green curves in Fig 10(a)). According to section VI D, the decay is exponentially fast as h → 0. When 2 < |σ| < 4, the complex eigenvalues converge, pairwise, at two opposite points on the imaginary axis outside the continuous spectrum band. After that one pair of imaginary eigenvalues immerses in the continuum while the other one remains outside the band, approaching the band edges only as h → 0. (See the dark-green arcs in Fig 9(d).) Regardless of the value of σ, the point h = H 1 (|σ|) where the odd-mode instability sets in, is to the left of H 3 (|σ|), the Hamiltonian Hopf bifurcation point where the even-mode ("light-green") complex quadruplet is born. When |σ| < 2, the even-mode quadruplet persists as h is decreased from H 3 to 0. Accordingly, in this σ range the "deepblue" real eigenvalues and their descendent "dark-green" complex quadruplet coexist with the complex quadruplet shown in light green. Therefore although the odd-mode instability is detectable in direct computer simulations of the antiphase soliton, its presence does not affect the soliton stability domain on the (γ, A)-plane.
Due to the coexistence of the real eigenvalues and the complex quadruplet, the odd-mode monotonic instability and the even-mode oscillatory one, develop simultaneously. A typical outcome of the soliton decay in the parameter range with |σ| < 2 is depicted in Fig 12 (b). The antiphase soliton breaks into an asymmetric assembly of two or more travelling breathers.
When 2 < |σ| < 4, the point h = H 1 (|σ|) is to the left of H 2 (|σ|), the lower boundary of the "light-green" quadruplet domain of existence. When |σ| falls in this range, the intervals of the odd-and even-mode instability do not overlap. However the odd-mode instability band is so narrow here that this instability would be hard to observe in computer simulations of the soliton. The odd instability band is marked as a white line in Fig 13(a).
We close this subsection by acknowledging an earlier numerical study of the eigenvalue problem (4.10). The authors of Ref [64] computed eigenvalues for three fixed values of h (h = 0.50, 0.71 and 3.16) and varied σ. This is obviously not equivalent to the approach of the present paper where we fix σ and vary h.

A. Conclusions
In this paper we have considered stability of the site-and bond-centred solitons in the chain of PT -symmetric dimers with gain and loss. Both classes of solitons come in two varieties: the in-phase and the antiphase solitons.
Our principal results are as follows.
1. We have shown that small perturbations about the soliton can be decomposed into a part tangent to the conservative invariant manifold containing the soliton, and a part that is transversal to that manifold. This decomposition proves instability of both varieties of the bond-centred solitons, for all amplitudes and all values of the chain parameters.
2. We have demonstrated that stability properties of a site-centred soliton with amplitude A in the chain with the inter-dimer coupling C and gain-loss coefficient γ, are completely determined just by two combinations of A, C and γ. These are the structural parameter σ = ±2 1 − γ 2 /C (where the top and bottom sign correspond to the in-phase and antiphase soliton, respectively) and the normalised amplitude h = A/ √ C. 3. We have proved that the in-phase site-centred soliton is stable if its amplitude lies below a threshold value (5.8) and unstable otherwise. A simple upper bound for the threshold, equation (5.9), has been established and its explicit expressions in the small-coupling limit (C/ 1 − γ 2 → 0) and in the limit where 1 − γ 2 /C → 0 have been provided. (See equations (5.11) and (5.10), respectively.) The lowest threshold occurs in chains where the value of the dimer-to-dimer coupling C is close to the value of the gain-to-loss coupling within each dimer.
4. With regard to the stability of the site-centered antiphase soliton, we have identified three characteristic coupling ranges: (a) C < 1/2, (b) 1/2 < C < 1.02, and (c) C > 1.02. 4(a). In the weak-coupling range (C < 1/2), the structure of the stability domain depends on whether the gainloss coefficient (i) is smaller than , the soliton is unstable if its amplitude satisfies A 2 < A < A 3 and stable otherwise. As C/ 1 − γ 2 → 0, the functions A 2 (C, γ) and A 3 (C, γ) are given by the asymptotic expressions (6.13) and (6.11), respectively. Otherwise A 2 and A 3 are as in (6.12) and (6.10), with H 2 and H 3 determined numerically. (ii) When γ falls between √ 1 − 4C 2 and √ 1 − 0.96C 2 , there are two instability bands: A 2 < A < A 3 and A 4 < A < A 5 , with A n as in (6.16) and H n determined numerically. (iii) When γ is greater than √ 1 − 0.96C 2 , the amplitudes above A 3 are stable and those below A 3 unstable.
In addition to the above one or two instability bands, we need to note an extremely thin corridor of odd-mode instability that crosses the soliton stability domain when C < 1/2 and γ < √ 1 − 0.96C 2 . The upper part of this corridor is described by the inequalities in (6.2), with the corresponding asymptotic expressions as in (6.4). 4(b). In the middle range (1/2 < C < 1.02), the stability criteria depend on whether the gain-loss coefficient is greater or smaller than √ 1 − 0.96C 2 . When γ < √ 1 − 0.96C 2 , the soliton is unstable if its amplitude falls in either of two bands, A 2 < A < A 3 or A 4 < A < A 5 -and stable otherwise. When γ > √ 1 − 0.96C 2 , the amplitudes above the threshold value A 3 are stable and those below the threshold are unstable. The boundaries A n (C, γ) (n = 2, ..., 5) are as in (6.16), with the functions H n being determined numerically. 4(c). In the strong coupling range (C > 1.02), solitons with amplitudes above A 3 are stable and those below A 3 unstable. The threshold amplitude is given by equation (6.10) where the function H 3 (|σ|) is determined numerically. 4(d) In the parameter regimes where the small-amplitude solitons are unstable (in the strong-coupling range C > 1.02 and in the range C < 1.02 with γ > √ 1 − 0.96C 2 ), the instability growth rate becomes exponentially small as the soliton's amplitude tends to zero. Hence unstable solitons with small amplitudes can live long enough to be regarded as practically stable.

B. Comparison with earlier work
It is instructive to compare our approach and conclusions with those of earlier authors who performed asymptotic analysis in the C, γ → 0 limit [51] and studied soliton stability numerically [40,51].
What we see as an advantage of our approach, is a particular choice of variables that leads to an eigenvalue problem in triangular form (see equations (4.1) and (4.8)). The eigenvectors of the triangular matrices admit a natural decomposition into two disjoint classes: the eigenvectors that belong to the scalar reduction of the two-component nonlinear Schrödinger equation (1.1) and eigenvectors that lie outside the reduction manifold.
This decomposition alone was sufficient for us to prove the instability of all bond-centred solitons, regardless of their amplitude and irrespective of the C and γ parameters of the chain. In contrast, the earlier approaches could only establish the instability of the bond-centred soliton in the C → 0, γ → 0 limit [51] -and for several isolated parameter triplets that were examined numerically [40,51]. Numerical methods require particular caution as they may not be accurate enough to discern small real parts of eigenvalues. This may mislead one into thinking that a weakly unstable soliton is stable (as in the case of the in-phase bond-centred soliton with large C [51].) In the case of the site-centred soliton, the above decomposition reduces the eigenvalue problem to the form familiar from the theory of the scalar nonlinear Schrödinger equation. This reduction allowed us to use the standard analytical methods for the (discrete) Schrödinger operators. The advantage of halving the dimension of the vector space extends , stable amplitudes lie on each side of the band A2 < A < A3. This corridor continues into the interval where it is joined by the second instability band, A4 < A < A5. In the interval √ 1 − 0.96C 2 < γ < 1, stable amplitudes lie only above A3. What looks like a white line through the stability domain is a very narrow corridor of the odd-mode instability. The upper boundary of this band is given by A1(γ); the lower boundary is determined numerically. The black dashed curve is given by (6.6); it bounds the region where the soliton is stable due to a rigorous analytical argument. Panel (b): When γ ≤ √ 1 − 0.96C 2 , there are two instability bands: A4 < A < A5 (lower) and A2 < A < A3 (upper band). In the interval √ 1 − 0.96C 2 < γ < 1, stable amplitudes lie only above A3. Panel (c): stable amplitudes lie only above A3(γ).
beyond the availability of analytics though. In particular, it simplifies power expansions in the anticontinuum limit where, unlike [51], we did not have to assume that γ is small.
Another difference between our approach and previous work is that we recognise that the stability properties of the soliton with amplitude A in the chain with coupling C and gain-loss coefficient γ, are completely determined just by two combinations of these three parameters. Letting h = A/ √ C and |σ| = 2 1 − γ 2 /C vary over their respective domains, we cover the entire (A, C, γ) space.
Our numerical results on the site-centred solitons in the limit C → 0, γ → 0, are in agreement with those of [51]. On the other hand, the eigenvalue trajectories of the site-centred antiphase soliton with larger C or γ close to 1the situation unexplored in [51] -are found to be much more complex than in the C, γ → 0 limit.

Eq.(A1) implies
Next, noting that and that R n is the null eigenvector of the matrix L − , we have n,m Using the inequality (A3), this gives n,m The lowest eigenvalue of the matrix L + is given by the minimum of the Rayleigh quotient Letting here y n = R n and using (A5), we obtain an upper bound on the lowest eigenvalue: Appendix B: The negative eigenvalue of L+ in the continuum limit In this and the next appendix, we establish the asymptotic behaviour of the discrete eigenvalues E 0 and E 1 as h → 0. Although (4.11) is exactly the matrix arising in the linearisation of the scalar discrete nonlinear Schrödinger soliton, these results are not in the existing literature on the subject.
In the limit h → 0, the site-centred soliton solution of equation (3.5) can be sought as R n = R(X n ), where R(X) is a continuous function of its argument, and X n = nh. The function R(X) satisfies and can be easily constructed as a power expansion Here e.s.t. stands for exponentially small terms. The first two terms in the expansion are and To determine the eigenvalues of the matrix L + in this limit we let y n = Y (X n ), where Y (X) is a continuous function. The system of linear algebraic equations is transformed into a differential eigenvalue problem and using (B2), we equate coefficients of like powers of h 2 in (B4). The coefficient of h 2α+2 gives us where we have introduced the Sturm-Liouville operator In particular, setting α = 0 we obtain from (B6): Eigenvectors y n decaying to zero as |n| → ∞ correspond to eigenfunctions Y (X) decaying as |X| → ∞; that is, discrete eigenvalues of (B3) correspond to discrete eigenvalues of (B8). The eigenvector associated with E 0 , the lowest eigenvalue of the matrix L + , is even in n: y −n = y n . We normalise it by y 0 = 1. The corresponding eigenfunction Y (0) should also be even in its argument: Y (0) (−X) = Y (0) (X). The eigenvalue of the operator L + , associated with an even eigenfunction, is E (0) = −3, and the eigenfunction is where the operator L + − E (0) in the left-hand side has a zero eigenvalue. A bounded solution of the above equation exists only if the right-hand side is orthogonal to the corresponding eigenvector, Y (0) . This requirement fixes the coefficient E (1) : E (1) = −1/5. Consequently, the asymptotic behaviour for the negative eigenvalue of L + as h → 0, is Higher order corrections to (B9) can be evaluated recursively, by considering equations (B6) with larger α.
Appendix C: The positive eigenvalue of L+ in the continuum limit This appendix is a continuation of the previous one. Here, we consider E 1 , the second lowest eigenvalue of the matrix L + . (This eigenvalue is positive [55,56].) We show that as h → 0, E 1 becomes exponentially small in h 2 .
First, differentiating (B1) with respect to X we observe that the derivative R X is an eigenfunction of the operator with the eigenvalue E = 0. Accordingly, equations (B6) with Y (α) = R (α) X and E (α) = 0 become a string of identities: Multiplying both sides of (C1) with the eigenfunction R (0) X of L + and integrating, these give Assume that Y (α) = R (α) X and E (α) = 0 for all 0 ≤ α ≤ α 0 − 1 with some α 0 , and consider equation (B6) with α = α 0 : This equation admits a bounded solution only if the right-hand side is orthogonal to the eigenfunction of L associated with its zero eigenvalue -that is, orthogonal to R (0) X . Making use of (C2) this solvability condition reduces to which implies E (α) = 0. Comparing the right-hand side of (C3) to the equation (C1), we obtain Y (α) with α = α 0 : X , where C is an arbitrary constant. Without loss of generality, we can choose C = 0; this is equivalent to dividing the eigenfunction Y (X) by the normalisation factor 1 + Ch 2α . This completes the proof.
The bottom line of our analysis is that all E (α) in the expansion (B5) are equal to zero and so E 1 is smaller than any power of h 2 in the h → 0 limit. Intuitively, this result is consistent with the exponential smallness of the obstacles to the translation motion of the discrete soliton (such as the value of the Melnikov function [56] or Peierls-Nabarro barrier [61]).
The function E 1 (h), computed numerically, is shown in Fig 5(b). The e −C/h 2 -type of decay of the eigenvalue as h → 0, is clearly visible in the figure.
Appendix D: Discrete eigenvalues of L+ in the anticontinuum limit In this appendix, we establish the asymptotic behaviours of the eigenvalues E 0 and E 1 as h → ∞.
To this end, we use the power-series solution of the scalar equation (3.5). The fundamental mode (the single-hump solution even in n) centred on the zero site, has the following asymptotic expansion as h → ∞: and The eigenvalue E 0 can be sought in the form and the corresponding eigenvector can also be written as a power series: Substituting (D1), (D2) and (D3) in (B3) and equating coefficients of like powers of h, we obtain the required asymptotic expansion Turning to the eigenvalue E 1 , associated with an odd eigenvector, we wish to show that where e.s.t. stands for the "exponentially small terms" as h → ∞. In (D5), the leading term h 2 corresponds to the lower edge of the continuous spectrum band.
It is convenient to introduce the distance of the eigenvalue to the continuum, E = E 1 −h 2 . The solution of equation (B3) with E = 0, decaying as n → ∞, is y n = e −κn . Here κ > 0 is defined by Letting y n = e −κn z n , the equation (B3) with n > 0, is cast in the form −e −κ z n+1 − e κ z n−1 + 2z n − 6 −1 R 2 n z n = E z n , n > 0, where we have denoted = h −2 . As n → ∞, the values z n approach a nonzero constant or grow in proportion to n. We supplement (D7) with the boundary condition corresponding to odd y n : The eigenvalue and each component of the eigenvector in (D6)-(D8) are functions of : E = E ( ), z n = z n ( ). When = 0, the eigenvector is z n (0) = n, and the eigenvalue E = 0. When = 0, we assume that an eigenvalue and eigenvector exist such that E ( ) → 0 − and z n ( ) → n as → 0. We wish to show that such E ( ) has to be nonanalytic in .
To this end, we assume the opposite; that is, we assume that E ( ) admits a Taylor expansion with some integer α > 0. When n is much greater than (α + 1)/2, the potential in (D7) becomes much smaller than E and can be disregarded: −e −κ z n+1 − e κ z n−1 + 2z n = −4 sinh 2 (κ/2)z n .
The only solution of (D10) that grows slower than exponentially as n → ∞, is z n ( ) = 1. (The other, linearly independent, solution is e 2κn .) However z n ( ) = 1 is not a continuous perturbation of z n (0) = n, and this fact contradicts our assumptions.
This leads us to conclude that E ( ) cannot have the Taylor expansion (D9). Instead, E has to be smaller than any power of . If that is the case, the potential in (D7) will remain greater than E even if n → ∞ so that the asymptotic behaviour of z n will not be governed by equation (D10).
The exponential approach of E 1 (h) to the value h 2 , as h → ∞, is clearly visible in Fig 5(b).

Appendix E: Stability eigenvalues for the in-phase soliton -upper and lower bounds
In this Appendix, we establish an upper and lower bounds on the lowest eigenvalue of the generalised eigenvalue problem (5.1). Here, σ > 0.
First, we evaluate the Rayleigh quotient (5.2) using g n = R n as a test function. Equation (A4) gives where S 2 and S 4 are as in (A2). Using (A3) this gives an upper bound on the ground-state eigenvalue −λ 2 : To bound −λ 2 from below, we note two inequalities, valid for an arbitrary g: Substituting these in (5.2) we obtain the lower bound: −λ 2 ≥ σ(E 0 + σ).
Making use of (B9), the lower bound becomes, in the limit h → 0, In order to obtain the lower bound in the opposite limit h → ∞, we use (D4) instead. This gives Appendix F: Symplectic eigenvalues with h → 0 and finite σ The subject of this Appendix is the eigenvalue problem (4.10) in the limit h → 0, with finite σ = 0. This limit corresponds to solitons of small amplitude (A → 0) in the chain with finite coupling C and gain-loss coefficient γ. We show that there are only two pairs of eigenvalues, with each eigenvalue being pure imaginary to all orders in h.
The derivation simplifies if we transform to new variables [47] A n = f n − ig n , B n = −(f n + ig n ).
(The alternative choice is ν 0 = −σ and B (0) = 0; this simply leads to the other member of the ±λ-pair.) The order h 2α , with α = 1, 2, ..., gives then a linear system where R β = 2 β µ=0 R (µ) R (β−µ) , β = 0, 1, .... Equation (F4) with α = 1 has the form of an eigenvalue problem where L 0 is an operator of the Pöschl-Teller variety: The Pöschl-Teller operator (F5) has only two discrete eigenvalues [47], ν is odd and has one zero crossing. (That operator L 0 does not have any other eigenvalues follows from the fact that its eigenvalue associated with the eigenfunction with n zero crossings cannot lie below the corresponding eigenvalue of the operator L + , equation (B7). The operator L + is well known to have only two discrete eigenvalues; hence the operator (F5) cannot have more than two.) Once B (0) (X) has been chosen, equation (F3) gives A (1) (X): The conclusion of our analysis is that the symplectic operator (4.15) has only two pairs of eigenvalues λ in the continuum limit h → 0 with σ = O(1). These are given by the following two-term asymptotic expansions: The notation reflects the fact that the first pair corresponds to even and the second one to odd eigenfunctions. A simple recursion argument proves that the eigenvalues λ even and λ odd remain pure imaginary to all orders in h. Indeed, consider some α ≥ 2 and assume we know A (β) , ν β with β = 0, 1, ..., α − 1 and B (β) with β = 0, 1, ..., α − 2. Assume all these quantities are real. The operator L 0 − ν 1 I in the left-hand side of equation (F4) is singular and the equation admits a bounded solution only if the right-hand side is orthogonal to B (0) , the function spanning its kernel space. This solvability condition determines ν α , a real value. Once the orthogonality condition has been satisfied, we can solve (F4) for B (α−1) . After that equation (F3) gives us A (α) . Both B (α−1) and A (α) come out real.
It is fitting to note here that our conclusion does not rule out the existence of an exponentially small real part of λ (Re λ ∼ e −κ/h , κ = const > 0). In particular, the exponentially small real part arises in the situation of small negative σ -see the next appendix.
When η → −∞, there are two complex quadruplets of eigenvalues. The imaginary parts are given [47] by the same equations (G3) whereas the real parts undergo an exponential decay (at equal rates) [52]: .

(G4)
The upshot of this analysis is that the symplectic problem (4.10) with σ < 0, where h and |σ| are both asymptotically small, has two quadruplets of complex eigenvalues. As h 2 /σ → 0, the real parts of these eigenvalues tend to zero faster than any power of h 2 /σ.
The first few coefficients in (H6) are then straightforward from (H2)-(H3): The outcome of our analysis is that the anti-continuum limit is characterised by a pair of opposite eigenvalues, with |λ| growing in proportion to h as h → ∞. The eigenvalues are real if σ > 0 and pure imaginary if σ < 0.