Bright solitons in a symmetric chain of dimers
Abstract
We study the existence and stability of fundamental bright discrete solitons in a paritytime ()symmetric coupler composed by a chain of dimers, that is modelled by linearly coupled discrete nonlinear Schrödinger 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 socalled anticontinuum 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 gainloss 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 blowups.
1 Introduction
A system of equations is symmetric if it is invariant with respect to combined parity () and timereversal () transformations. The symmetry is interesting as it forms a particular class of nonHermitian Hamiltonians in quantum mechanics [1], that may have a real spectrum up to a critical value of the complex potential parameter, above which the system is in the ’broken symmetry’ phase [2, 3, 4].
The most basic configuration having symmetry is a dimer, i.e. 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 systems as the concept of symmetry was first realised experimentally on dimers consisting of two coupled optical waveguides [5, 6] (see also the review [7] for symmetry in optical applications). The experiments have been rapidly followed by many other observations of symmetry in different branches of physics, from mechanical to electrical analogues (see the review [8]).
When nonlinearity is present in a system, nontrivial behaviours may emerge that cease to exist in the linear case, such as the presence of blowup dynamics in the parameter region of the unbroken phase in the linear counterpart [19, 20, 21]. 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 [22]. 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 gainloss have been studied in [9, 10, 11, 12], 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 become stable after a tangent (saddlecenter) bifurcation. When one adds a gain and loss term in each arm, one obtains symmetric couplers, which have been considered in [13, 14, 15, 16, 17, 18]. In the presence of the lineargain and loss terms, asymmetric solitons cease to exist, while antisymmetric solitons are always unstable [16], even though those with small amplitudes can live long due to weak underlying instability [13]. Symmetric solitons can be stable in a similar fashion to those in the system without gainloss [16].
The stability of bright discrete solitons in symmetric couplers was discussed in [22] using variational methods, where it was shown that symmetric onsite solutions can be stable and there is a critical solution amplitude above which the symmetry is broken. The case when the polarity of the 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 symmetrybreaking bifurcation at a finite value of the energy, similarly to that in the continuous counterpart [9, 10, 11, 12]. Recently a similar chain of dimers with a slightly different nonlinearity was derived [25] to describe coupled chains of parametrically driven pendula as a mechanical analogue of symmetric systems [26]. The stability of bright discrete solitons was established through the applications of the Hamiltonian energy and an index theorem. The nonlinear longtime 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 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 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 gainloss 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 socalled 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 gainloss 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.
2 Mathematical model
The governing equations describing symmetric chains of dimers are of the form [22]
(1) 
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, , are complexvalued wave function at site , is the constant coefficient of the horizontal linear coupling (coupling constant between two adjacent sites), and are the discrete Laplacian term in one spatial dimension, the gain and loss acting on complex variables , are represented by the positive coefficient , i.e. . The nonlinearity coefficient is denoted by , which can be scaled to without loss of generality due to the case of focusing nonlinearity that we consider. Bright discrete soliton solutions satisfy the localisation conditions as .
The focusing system has static localised solutions that can be obtained from substituting
(2) 
into (1) to yield the equations
(3) 
where , are complexvalued and the propagation constant .
3 Solutions of weakly coupled equations
In the uncoupled limit, i.e. when , the chain (1) becomes the equations for the dimer. The static equation (3) has been analysed in details in [31, 29], where it was shown that there is a relation between and above which there is no timeindependent 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
(4) 
By substituting the above expansions into Eq. (3) and collecting the terms in successive powers of , one obtains at and , respectively, the equations
(5) 
and
(6) 
It is wellknown that there are two natural fundamental solutions representing bright discrete solitons that may exist for any , from the anticontinuum to the continuum limit, i.e. an intersite (twoexcitedsite) and onsite (oneexcitedsite) bright discrete mode. Here, we will limit our study to these two fundamental modes.
3.1 Intersite soliton
In the uncoupled limit, the mode structure for the intersite soliton is of the form
(7) 
with [31]
(8) 
which is an exact solution of Eq. (5). Note that (8) will have no real solution when . This is the broken region of symmetry. The parameter can be taken as 0, due to the gauge phase invariance of the governing equation (1) and henceforth . The former phase corresponds to the socalled 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 and are the necessary condition for soliton I and II, respectively.
For the first order correction due to the weak coupling, writing , and substituting it into equations (6) will yield
(9) 
Equations (4),(7),(8),(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.
3.2 Onsite soliton
4 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 substitute this into Eq. (1) to obtain the linearised equations at )
(12) 
which have to be solved for the eigenvalue and the corresponding eigenvector . As the stability matrix of the eigenvalue problem (12) is real valued, and are also eigenvalues with corresponding eigenvectors and with , respectively. Therefore, we can conclude that the solution is (linearly) stable only when for all eigenvalues .
4.1 Continuous spectrum
The spectrum of (12) will consist of continuous spectrum and discrete spectrum (eigenvalue). To investigate the former, we consider the limit , introduce the planewave ansatz , , and substitute the ansatz into (12) to obtain
(13) 
where . The equation can be solved analytically to yield the dispersion relation
(14) 
The continuous spectrum is therefore given by and with the spectrum boundaries
(15)  
(16) 
obtained from (14) by setting and in the equation.
4.2 Discrete spectrum
Following the weakcoupling analysis as in Section 3, we will as well use similar asymptotic expansions to solve the eigenvalue problem (12) analytically, i.e. we write
(17) 
with . We then substitute the expansions into the eigenvalue problem (12).
At order , 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 (17) further as
. Hence, we have two small parameters, i.e. and , that are independent of each other. For the sake of presentation, the detailed calculations are shown in the Appendix. Below 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
(18) 
and
(19) 
4.2.2 Intersite soliton II
The intersite soliton II, i.e. the intersite soliton that is antisymmetric between the arms, has three pairs of eigenvalues given by
(20) 
and
(21) 
4.2.3 Onsite soliton I
The onsite soliton has only one eigenvalue for small given asymptotically by
(22) 
4.3 Onsite soliton II
As for the second type of the onsite soliton, we have
(23) 
5 Numerical results
We have solved the steadystate equation (3) numerically using a NewtonRaphson method and analysed the stability of the numerical solution by solving the eigenvalue problem (12). Below we will compare the analytical calculations obtained above with the numerical results.
First, we consider the discrete intersite soliton I. We show in Fig. 1 the spectrum of the soliton as a function of the coupling constant for and . 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 Fig. 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 all the nonzero eigenvalues satisfy (see (1516)), then we will obtain unstable soliton I in the continuum limit . However, when in the anticontinuum limit all the nonzero eigenvalues satisfy , 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 Fig. 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 (the top panels of the figure), we have the condition that all the nonzero eigenvalues satisfy in the anticontinuum limit . The collision between the eigenvalues and the continuous spectrum as the coupling increases creates complex eigenvalues. In the second case using larger (lower panels of the figure), the nonzero eigenvalues satisfy when . 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 Figs. 4 and 5 is the stability of discrete solitons type I and II, respectively.
In Fig. 4, the top left panel shows that for 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 , even though initially in the uncoupled limit the nonzero eigenvalue satisfies , one may obtain an exponential instability (i.e. instability due to a real eigenvalue). The bottom left panel shows the case when the discrete soliton is already unstable even in the uncoupled limit due to the nonzero eigenvalue that is already realvalued.
Fig. 5 shows that the antisymmetric solitons are generally unstable due to a quartet of complex eigenvalues, as shown in the left panels of the figure. When the coupling is increased further, there will be an eigenvalue bifurcating from 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 Fig. 6, we present the (in)stability region of the two types of discrete solitons in the plane for three values of the gainloss parameter . Discrete solitons are unstable above the curves. 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 the left panel generally corresponds to an eigenvalue crossing the origin that becomes realvalued, while the curves in the other panel are due to the appearance of a quartet of complex eigenvalues. In general, we obtain that the gainloss term can be parasitic as it reduces the stability region of the discrete solitons.
Finally, we present in Fig. 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, i.e. in the form of travelling discrete solitons or solution blowups. 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.
6 Conclusion
We have presented a systematic method to determine the stability of discrete solitons in a 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 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 symmetry as future investigations.
Acknowledgement
AAB and HS are grateful to the University of Nottingham for the 2013 Visiting Fellowship Scheme and British Council for the 2015 Indonesia Second City Partnership Travel Grant.
The authors declare that there is no conflict of interest regarding the publication of this paper.
Appendix A Analytical calculation
As mentioned in Section 4.2, to solve the eigenvalue problem (12) analytically we expand the eigenvalue and eigenvector asymptotically as
(24) 
with .
Performing the expansion in , at we obtain the following set of equations
(25) 
where
At and , we obtain
(26)  
(27)  
where
The steps of finding the coefficients of the asymptotic expansions, , are as follows.

Solve the eigenvalue problem (25), which is a system of equations, for and .

Determine by taking the vector inner product of both sides of (26) with the nullspace of the Hermitian transpose of the block matrix that consists of along the diagonal, where is the identity matrix.

Solve (26) for .

Determine by taking the vector inner product of both sides of (27) with the nullspace of the Hermitian transpose of the block matrix that consists of .
The procedure repeats if one would like to calculate the higher order terms.
The leading order eigenvalue of (25) has been solved in [31]. However, the expression of the corresponding eigenvector was very lengthy, that makes it almost impractical to be used to determine the higher order corrections of . Therefore, in every equation at order obtained from (12), we also expand the variables in , i.e.
where again , and obtain equations at order . The steps to determine and are the same as mentioned above.
References
References
 [1] N. Moiseyev, NonHermitian Quantum Mechanics (Cambridge University Press, Cambridge, 2011).
 [2] C. M. Bender and S. Boettcher, Phys. Rev. Lett. 80, 5243 (1998).
 [3] C. M. Bender, S. Boettcher, and P. N. Meisinger, J. Math. Phys. 40, 2201 (1999).
 [4] C. M. Bender, Rep. Prog. Phys. 70, 947 (2007).
 [5] A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. VolatierRavat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides, Phys. Rev. Lett. 103, 093902 (2009).
 [6] C.E. Rüter, K. G. Makris, R. ElGanainy, D. N. Christodoulides, M. Segev, and D. Kip, Nat. Phys. 6, 192 (2010)
 [7] S.V. Suchkov, A. A. Sukhorukov, J. Huang, S. V. Dmitriev, C. Lee, and Yu. S. Kivshar, Laser & Photon. Rev. 10, 177 (2016).
 [8] V.V. Konotop, J. Yang, D.A. Zezyulin, Rev. Mod. Phys. 88, 035002 (2016).
 [9] P. L. Chu, B. A. Malomed, and G. D. Peng, J. Opt. Soc. Am. B 10, 1379 (1993).
 [10] E. M. Wright, G. I. Stegeman, and S. Wabnitz, Phys. Rev. A 40, 4455 (1989).
 [11] N. Akhmediev and J. M. SotoCrespo, Phys. Rev. E 49, 4519 (1994).
 [12] B. A. Malomed, I. M. Skinner, P. L. Chu, and G. D. Peng, Phys. Rev. E 53, 4084 (1996).
 [13] R. Driben and B. A. Malomed, Opt. Lett. 36, 4323 (2011).
 [14] F. K. Abdullaev, V. V. Konotop, M. Ogren, and M. P. Sørensen, Opt. Lett. 36, 4566 (2011).
 [15] R. Driben and B. A. Malomed, Europhys. Lett. 96, 51001 (2011).
 [16] N. V. Alexeeva, I. V. Barashenkov, A. A. Sukhorukov, and Y. S. Kivshar, Phys. Rev. A 85, 063837 (2012).
 [17] I. V. Barashenkov, S. V. Suchkov, A. A. Sukhorukov, S. V. Dmitriev, and Y. S. Kivshar, Europhys. Lett. 99, 54001 (2012).
 [18] P. Li, L. Li, and B.A. Malomed, Phys. Rev. E 89, 062926 (2014)
 [19] J. Pickton and H. Susanto, Phys. Rev. A 88, 063840 (2013)
 [20] P. G. Kevrekidis, D. E. Pelinovsky, and D. Y. Tyugin, J. Phys. A 46, 365201 (2013).
 [21] I. V. Barashenkov, G. S. Jackson, and S. Flach, Phys Rev A 88, 053817 (2013).
 [22] S V Suchkov, B A Malomed, S V Dmitriev, Yu. S Kivshar, Phys. Rev. E 84, 046609 (2011).
 [23] J. D’Ambroise, P.G. Kevrekidis, and B.A. Malomed, Phys. Rev. E 91, 033207 (2015).
 [24] G. Herring, P. G. Kevrekidis, B. A. Malomed, R. CarreteroGonzález, and D. J. Frantzeskakis, Phys. Rev. E 76, 066606 (2007).
 [25] A Chernyavsky and D E Pelinovsky, Symmetry 8, 59 (2016).
 [26] C.M. Bender, B. Berntson, D. Parker and E. Samuel, Am. J. Phys. 81, 173 (2013).
 [27] A Chernyavsky and D E Pelinovsky, arXiv:1606.02333
 [28] R.S. MacKay and S. Aubry, Nonlinearity 7, 1623 (1994)
 [29] V.V. Konotop, D.E. Pelinovsky, and D.A. Zezyulin, EPL 100, 56006 (2012)
 [30] D.E. Pelinovsky, D.A. Zezyulin and V.V. Konotop, J. Phys. A: Math. Theor. 47, 085204 (2014)
 [31] K. Li and P. G. Kevrekidis, Phys. Rev. E 83, 066608 (2011)
 [32] Z. Chen, J. Liu, S. Fu, Y. Li, and B.A. Malomed, Optics Express 22, 2967929692 (2014).