Bright solitons in a PT-symmetric chain of dimers

We study the existence and stability of fundamental bright discrete solitons in a parity-time (PT)-symmetric coupler composed by a chain of dimers, that is modelled by linearly coupled discrete nonlinear Schrodinger equations with gain and loss terms. We use a perturbation theory for small coupling between the lattices to perform the analysis, which is then confirmed by numerical calculations. Such analysis is based on the concept of the so-called anti-continuum limit approach. We consider the fundamental onsite and intersite bright solitons. Each solution has symmetric and antisymmetric configurations between the arms. The stability of the solutions is then determined by solving the corresponding eigenvalue problem. We obtain that both symmetric and antisymmetric onsite mode can be stable for small coupling, on the contrary of the reported continuum limit where the antisymmetric solutions are always unstable. The instability is either due to the internal modes crossing the origin or the appearance of a quartet of complex eigenvalues. In general, the gain-loss term can be considered parasitic as it reduces the stability region of the onsite solitons. Additionally, we analyse the dynamic behaviour of the onsite and intersite solitons when unstable, where typically it is either in the form of travelling solitons or soliton blow-ups.


Introduction
A system of equations is PT-symmetric if it is invariant with respect to combined parity (P) and time-reversal (T) transformations. The symmetry is interesting as it forms a particular class of non-Hermitian Hamiltonians in quantum mechanics [1], which may have a real spectrum up to a critical value of the complex potential parameter, above which the system is in the "broken PT symmetry" phase [2][3][4].
The most basic configuration having PT symmetry is a dimer, that is, a system of two coupled oscillators where one of the oscillators has damping losses and the other one gains energy from external sources. Considerably, dimers are also the most important PT systems as the concept of PT symmetry was first realised experimentally on dimers consisting of two coupled optical waveguides [5,6] (see also the review [7] for PT symmetry in optical applications). The experiments have been rapidly followed by many other observations of PT symmetry in different branches of physics, from mechanical to electrical analogues (see the review [8]).
When nonlinearity is present in a PT system, one may have nontrivial behaviours that do not exist in the linear case, such as the presence of blow-up dynamics in the parameter region of the unbroken phase in the linear counterpart [9][10][11]. When nonlinear dimers are put in arrays where elements with gain and loss are linearly coupled to the elements of the same type belonging to adjacent dimers, one can also obtain a distinctive feature in the form of the existence of solutions localized in space as continuous families of their energy parameter [12]. The system therefore has two arms with each arm described by a discrete nonlinear Schrödinger equation with gain or loss. Here, we study the nonlinear localised solutions, which loosely we also refer to as bright discrete solitons, and their stability analytically and numerically.
In the continuous limit, the coupled equations without gain-loss have been studied in [13][14][15][16], where it has been shown that the system admits symmetric, antisymmetric, and asymmetric solitons between the arms. Unstable asymmetric solutions bifurcate from the symmetric ones through a subcritical symmetry breaking bifurcation, which then 2 Advances in Mathematical Physics become stable after a tangent (saddle-center) bifurcation. When one adds a gain and loss term in each arm, one obtains PT-symmetric couplers, which have been considered in [17][18][19][20][21][22]. In the presence of the linear-gain and loss terms, asymmetric solitons cease to exist, while antisymmetric solitons are always unstable [20], even though those with small amplitudes can live long due to weak underlying instability [17]. Symmetric solitons can be stable in a similar fashion to those in the system without gain-loss [20].
The stability of bright discrete solitons in PT-symmetric couplers was discussed in [12] using variational methods, where it was shown that symmetric onsite solutions can be stable and there is a critical solution amplitude above which the PT symmetry is broken. The case when the polarity of the PT-symmetric dimers is staggered along the chain is considered in [23]. The same equations without gain and loss were considered in [24] where the symmetric soliton loses its stability through the symmetry-breaking bifurcation at a finite value of the energy, similarly to that in the continuous counterpart [13][14][15][16]. Recently, a similar PT chain of dimers with a slightly different nonlinearity was derived [25] to describe coupled chains of parametrically driven pendula as a mechanical analogue of PT-symmetric systems [26]. The stability of bright discrete solitons was established through the applications of the Hamiltonian energy and an index theorem. The nonlinear long-time stability of the discrete solitons was also established using the Lyapunov method in the asymptotic limit of a weak coupling between the pendula [27].
In this work, we determine the eigenvalues of discrete solitons in PT-symmetric couplers analytically using asymptotic expansions. The computation is based on the socalled method of weak coupling or anticontinuum limit. The application of the method in the study of discrete solitons was formulated rigorously in [28] for conservative systems. It was then applied to PT-symmetric networks in [29,30]. However, no explicit expression of the asymptotic series of the eigenvalues for the stability of discrete solitons has been presented before. Here, in addition to the asymptotic limit of weak coupling between the dimers, we also propose to consider expansions in the coefficient of the gain-loss terms. In this case, explicit computations of the asymptotic series of the eigenvalues become possible.
The manuscript is outlined as follows. In Section 2, we present the mathematical model. In Section 3, we use perturbation theory for small coupling to analyse the existence of fundamental localised solutions. Such analysis is based on the concept of the so-called anticontinuum limit approach. The stability of the solitons is then considered analytically in Section 4 by solving a corresponding eigenvalue problem. In this section, in addition to small coupling, the expansion is also performed under the assumption of small coefficient of the gain-loss term due to the nonsimple expression of the eigenvectors of the linearised operator. The findings obtained from the analytical calculations are then compared with the numerical counterparts in Section 5. We also produce stability regions for all the fundamental solitons numerically. In this section, we present the typical dynamics of solitons in the unstable parameter ranges by direct numerical integrations of the governing equation. We present the conclusion in Section 6.

Mathematical Model
The governing equations describing PT-symmetric chains of dimers are of the form [12] The derivative with respect to the evolution variable (i.e., the propagation distance, if we consider their application in fiber optics) is denoted by the overdot, = ( ) and V = V ( ) are complex-valued wave function at site ∈ Z, > 0 is the constant coefficient of the horizontal linear coupling (coupling constant between two adjacent sites), Δ 2 = ( +1 − 2 + −1 ) and Δ 2 V = (V +1 − 2V + V −1 ) are the discrete Laplacian term in one spatial dimension, the gain and loss acting on complex variables and V are represented by the positive coefficient ; that is, > 0. The nonlinearity coefficient is denoted by , which can be scaled to +1 without loss of generality due to the case of focusing nonlinearity that we consider. Bright discrete soliton solutions satisfy the localisation conditions , V → 0 as → ±∞.
The focusing system has static localised solutions that can be obtained from substituting into (1) to yield the equations where and are complex-valued and the propagation constant ∈ R.

Solutions of Weakly Coupled Equations
In the uncoupled limit, that is, when = 0, chain (1) becomes the equations for the dimer. The static equation (3) has been analysed in detail in [29,31], where it was shown that there is a relation between and above which there is no time-independent solution to (3) (see also the analysis below). When is nonzero, but small enough, the existence of solutions emanating from the uncoupled limit can be shown using the Implicit Function Theorem. The existence analysis of [25] can be adopted here despite the slightly different nonlinearity as the Jacobian of our system when uncoupled shares a rather similar invertible structure (see also [29,30] that have the same nonlinearity in the governing equations but different small coupling terms). However, below we will not state the theorem and instead derive the asymptotic series of the solutions.
Using perturbation expansion, solutions of the coupler (3) for small coupling constant can be expressed analytically as By substituting the above expansions into (3) and collecting the terms in successive powers of , one obtains at O(1) and O( ), respectively, the equations It is well-known that there are two natural fundamental solutions representing bright discrete solitons that may exist for any , from the anticontinuum to the continuum limit, that is, an intersite (two-excited-site) and onsite (one-excitedsite) bright discrete mode. Here, we will limit our study to these two fundamental modes.

Intersite Soliton.
In the uncoupled limit, the mode structures (0) and (0) for the intersite soliton are of the form with [31]̂0 which is an exact solution of (5). Note that (8) will have no real solution when | | > 1. This is the broken region of PT-symmetry. The parameter can be taken as 0, due to the gauge phase invariance of the governing equation (1) and henceforth = arcsin and − arcsin . The former phase corresponds to the so-called symmetric configuration between the arms, while the latter is called antisymmetric one. Herein, we also refer to the symmetric and antisymmetric soliton as soliton I and II, respectively. Equation (8) informs us that > √1 − 2 > 0 and > −√1 − 2 are the necessary conditions for solitons I and II, respectively.
For the first-order correction due to the weak coupling, writing (1) =̃, 1 and (1) =̃, 1 and substituting them into (6) will yield Equations (4), (7), (8), and (9) are the asymptotic expansion of the intersite solitons. One can continue the same calculation to obtain higher order corrections. Here, we limit ourselves to the first-order correction only, which is sufficient to determine the leading order behaviour of the eigenvalues later.

Onsite Soliton.
For the onsite soliton, that is, a oneexcited-site discrete mode, one can perform the same computations to obtain the mode structure of the form with (8). After writing (1) =̃, 1 and (1) =̃, 1 , the first order correction from (6) is given bỹ

Stability Analysis
After we find discrete solitons, their linear stability is then determined by solving a corresponding linear eigenvalue problem. To do so, we introduce the linearisation ansatz = ( +̃( + ) ) and V = ( +̃( + ) ) , |̃| ≪ 1, and substitute this into (1) to obtain the linearised equations at O(̃): Advances in Mathematical Physics

Discrete Spectrum.
Following the weak-coupling analysis as in Section 3, we will as well use similar asymptotic expansions to solve the eigenvalue problem (12) analytically; that is, we write with ◻ = , , , , . We then substitute the expansions into the eigenvalue problem (12). At order O(1), one will obtain the stability equation for the dimer, which has been discussed for a general value of in [31]. The expression of the eigenvalues is simple, but the expression of the corresponding eigenvectors is not, which makes the result of [31] rather impractical to use. Therefore, here we limit ourselves to the case of small | | and expand (16) further as = 0, 1, 2, . . .. Hence, we have two small parameters, that is, and , that are independent of each other. For the sake of presentation, the detailed calculations are shown in the Appendix. Here we will only cite the final results.

4.2.1.
Intersite Soliton I. The intersite soliton I (i.e., the symmetric intersite soliton) has three pairs of eigenvalues for small and . One pair bifurcate from the zero eigenvalue. They are asymptotically given by

Intersite Soliton II.
The intersite soliton II, that is, the intersite soliton that is antisymmetric between the arms, has three pairs of eigenvalues given by

Onsite Soliton I.
The onsite soliton has only one eigenvalue for small given asymptotically by

Numerical Results
We have solved the steady-state equation (3) numerically using a Newton-Raphson method and analysed the stability of the numerical solution by solving the eigenvalue problem (12). Here we will compare the analytical calculations obtained above with the numerical results.
First, we consider the discrete intersite soliton I. We show in Figure 1 the spectrum of the soliton as a function of the coupling constant for = 1.2 and = 0, 0.5. On the real axis, one can observe that there is only one unstable eigenvalue that bifurcates from the origin. As the coupling increases, the bifurcating eigenvalue enters the origin again when → ∞. Hence, in that limit we obtain a stable soliton I (i.e., a stable symmetric soliton). The dynamics of the nonzero eigenvalues as a function of the coupling constant is shown in the right panels of the figure, where one can see that the eigenvalues are on the imaginary axis and simply enter the continuous spectrum as increases.
In Figure 2, we plot the eigenvalues for large enough. Here, in the uncoupled limit, all the three pairs of eigenvalues are on the real axis. As the coupling increases, two pairs go back toward the origin, while one pair remains on the real axis (not shown here). In the continuum limit → ∞, we therefore obtain an unstable soliton I (i.e., an unstable symmetric soliton).
In both figures, we also plot the approximate eigenvalues in solid (blue) curves, where good agreement is obtained for small .
From numerical computations, we conjecture that if in the limit → 0 all the nonzero eigenvalues satisfy 2 > 2 1− (see (15)), then we will obtain unstable soliton I in the continuum limit → ∞. However, when in the anticontinuum limit → 0 all the nonzero eigenvalues satisfy 2 1+ < 2 < 2 2− , we may either obtain a stable or an unstable soliton I in the continuum limit. Next, we consider intersite solitons II (i.e., antisymmetric intersite solitons). Shown in Figure 3 is the spectrum of the discrete solitons for two values of . In both cases, there is an eigenvalue bifurcating from the origin. For the smaller value of (Figures 3(a)-3(c)), we have the condition that all the nonzero eigenvalues satisfy 2 < 2 2− in the anticontinuum limit → 0. The collision between the eigenvalues and the continuous spectrum as the coupling increases creates complex eigenvalues. In the second case using larger (Figures 3(d)-3(f)), the nonzero eigenvalues satisfy 2 > 2 1− when = 0. Even though not seen in the figure, the collision between one of the nonzero eigenvalues and the continuous spectrum also creates a pair of complex eigenvalues. Additionally, in the continuum limit both values of as well as the other values of the parameter that we computed for this type of discrete solitons yield unstable solutions.
We also study onsite solitons. Shown in Figures 4 and 5 is the stability of discrete solitons types I and II, respectively. Figure 4(a) shows that for ( − √1 − 2 ) small enough we will obtain stable discrete solitons. For coupling constant small, we indeed show it through our analysis depicted as the blue solid line. Numerically, we obtain that this soliton is also stable in the continuum limit → ∞. However, when is large enough compared to √1 − 2 , even though initially in the uncoupled limit the nonzero eigenvalue satisfies 2 < 2 2− , one may obtain an exponential instability (i.e., instability due to a real eigenvalue). Figure 4(c) shows the case when the discrete soliton is already unstable even in the uncoupled limit due to the nonzero eigenvalue that is already real-valued. Figure 5 shows that the antisymmetric solitons are generally unstable due to a quartet of complex eigenvalues, as shown in Figures 5(a) and 5(c). When the coupling is increased further, there will be an eigenvalue bifurcating from ± 1− that will move towards the origin and later becomes a pair of real eigenvalues. These solitons are also unstable in the continuum limit.
Unlike intersite discrete solitons that are always unstable, onsite discrete solitons may be stable. In Figure 6, we present the (in)stability region of the two types of discrete solitons in the ( , )-plane for three values of the gain-loss parameter . Discrete solitons are unstable above the curves in panel (a) and between the curves in panel (b). Indeed as we mentioned before, for soliton I there is a critical that depends on below which the soliton is stable in the continuum limit, while soliton II is always unstable in that limit. Another difference between the two figures is that the stability curves in Figure 6(a) generally correspond to an eigenvalue crossing the origin that becomes real-valued, while the curves in the other panel are due to the appearance of a quartet of complex eigenvalues. In general, we obtain that the gain-loss term can be parasitic as it reduces the stability region of the discrete solitons.
Finally, we present in Figure 7 the time dynamics of some of the unstable solutions shown in the previous figures. What we obtain is that typically there are two kinds of dynamics, that is, in the form of travelling discrete solitons or solution blow-ups. The first type was the typical dynamics of the intersite soliton I. The second dynamics is typical for the other types of unstable discrete solitons.

Conclusion
We have presented a systematic method to determine the stability of discrete solitons in a PT-symmetric coupler by computing the eigenvalues of the corresponding linear eigenvalue problem using asymptotic expansions. We have compared the analytical results that we obtained with numerical computations, where good agreement is obtained. From the numerics, we have also established the mechanism of instability as well as the stability region of the discrete solitons. The application of the method in higher dimensional PT-symmetric couplers (see, e.g., [32]) is a natural extension of the problem that is addressed for future work. Additionally, we also address the computation of eigenvalues of discrete solitons in the neighbourhood of broken PTsymmetry as future investigations.
(ii) Determine (1) by taking the vector inner product of both sides of (A.4) with the null-space of the Hermitian transpose of the block matrix that consists of ( 0 − (0) 4 ) along the diagonal, where 4 is the 4×4 identity matrix. (iv) Determine (2) by taking the vector inner product of both sides of (A.5) with the null-space of the Hermitian transpose of the block matrix that consists of ( 0 − (0) 4 ). The procedure repeats if one would like to calculate the higher order terms.