PT-Symmetric dimer in a generalized model of coupled nonlinear oscillators

PT-Symmetric dimer in a generalized model of coupled nonlinear oscillators Jes´ us Cuevas–Maraver Nonlinear Physics Group, Departamento de F´ ısica Aplicada I, Universidad de Sevilla, Escuela Polit´ ecnica Superior, C/ Virgen de ´ Africa, 7, 41011-Sevilla, Spain Instituto de Matem´ aticas de la Universidad de Sevilla (IMUS). Edificio Celestino Mutis. Avda. Reina Mercedes s/n, 41012-Sevilla, Spain Avinash Khare Indian Institute of Science Education and Research (IISER), Pune 411008, India Panayotis G. Kevrekidis and Haitao Xu Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA 01003-9305, USA Avadh Saxena Center for Nonlinear Studies and Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA In the present work, we explore the case of a general PT -symmetric dimer in the context of two both linearly and nonlinearly coupled cubic oscillators. To obtain an analytical handle on the system, we first explore the rotating wave approximation converting it into a discrete nonlinear Schr¨ odinger type dimer. In the latter context, the stationary solutions and their stability are identified numerically but also wherever possible analytically. Solutions stemming from both symmetric and anti-symmetric special limits are identified. A number of special cases are explored regarding the ratio of coefficients of nonlinearity between oscillators over the intrinsic one of each oscillator. Finally, the considerations are extended to the original oscillator model, where periodic orbits and their stability are obtained. When the solutions are found to be unstable their dynamics is monitored by means of direct numerical simulations. I. INTRODUCTION The notion of parity-time (PT ) symmetry has recently been receiving increasing attention over a wide variety of settings; see e.g. [1–3]. The original proposal involved a non-Hermitian variant of quantum mechanics which might still produce real eigenvalues (and hence be associated with measurable quantities). However, it was instead the analogy of this model with the paraxial approximation in optics which led to the proposal that such mathematical constructs can be realized in optical settings [4, 5], and which eventually led to their experimental realization [6]. This series of developments, in turn, prompted researchers towards a more detailed understanding of the stationary states of such PT-systems (and how they differ from their Hamiltonian analogues), an effort to appreciate their stability properties and finally an attempt to quantify their nonlinear dynamics. This effort emerged both at the level of few-site configurations [7–15] (which were chiefly experimentally accessible), as well as at that of infinite-size lattices [16–18]. Although the quantum-mechanical and paraxial-optical focal points of this activity have provided an emphasis on the study of Schr¨ odinger type settings, a number of recent studies, especially on the experimental side, have led to an increased interest in oscillator systems (which one can think of as oligomers -few site settingsof the Klein-Gordon type). More specifically, a mechanical system realizing PT -symmetry has been proposed in [19], while a major thrust of research has focused on the context of electronic circuits; see e.g. the original realization of [20] and the more recent review of this activity in [21]. As an aside, we note that additional intriguing realizations of PT -symmetry have also emerged e.g. in the realm of whispering-gallery microcavities [22]. Mostly, the efforts on this oscillator realm have been limited to the study of linear systems, yet recently a number of nonlinear variants have been explored both theoretically/numerically and even experimentally. As notable such examples, we mention the split-ring resonator chain in the context of magnetic metamaterials proposed in the work of [23], as well as the experimental realization of a PT-symmetric dimer of Van-der-Pol oscillators that arose in the work of [24]. On the theoretical side, some of these studies raised a number of intriguing theoretical questions. For instance, the theoretical modeling of the linear PT -symmetric analogue of the system [22] led to the realization that such linear oscillator pairs may be Hamiltonian although one of them has gain and the other has loss [25]. This, in turn, led the authors of [26] to appreciate that this feature (the Hamiltonian nature of a PT-symmetric oscillator system) can be extended to the nonlinear case, if the nonlinearity contains both selfand crossinteractions between the oscillators and if these interactions have an appropriate ratio (the ratio utilized between crossand self-interactions in that work was 3). Importantly, the latter work also extended consideration of that model to the Schr¨ odinger variant thereof (through a multiple scales expansion), finding that nonlinearity may, in that context, “soften” the PT -symmetric phase transition. That is, it may enable the existence of stable periodic and quasi-periodic states at any value of the gain-loss parameter γ. Our aim in the present work is to revisit this context of two coupled nonlinear oscillators, one of which bears gain and the 2 other loss. We will consider the nonlinear case (almost exclusively, briefly touching upon the linear case as a special limit). Importantly, we will also explore the ratio of crossto self-interaction of the oscillators as a free parameter. Interestingly, this will enable us to identify a series of special cases, including the integrable one recently explored in [26]. For all values of this nonlinear parameter and as a function of the loss/gain parameter (γ) and of the frequency parameter ωb, we will study systematically both dimer systems. That is, we will first derive and analyze the discrete nonlinear Schr¨ odinger (DNLS) dimer, to obtain a simplified understanding of the existence, stability and dynamics properties. Then, in a way reminiscent of our earlier work (involving no cross interactions) [27] –and complementing the earlier work of [26] who did not focus on the periodic orbit solutions of the original nonlinear oscillator dimer–, we will return to the oscillator system and explore its own solutions, in terms of their existence, stability and dynamical properties. When the solutions are identified as unstable, a brief discussion will also be given of their dynamical evolution. This paper is organized as follows. In the next section (II) we provide the model equations, discuss their symmetries and potential Hamiltonian structure and indicate where corresponding exact solutions for the generalized coupled nonlinear oscillators can be obtained. In Sec. III we invoke the rotating wave approximation (RWA) and provide the stability equations and analytical results as well as perform numerical analysis of both the symmetric and asymmetric solutions for the resulting generalized DNLS dimer. Section IV contains the corresponding analysis of the Klein-Gordon dimer. A discussion of the dynamics of unstable solutions is given in Sec. V. Our main results and conclusions are summarized in Sec. VI, where a number of directions for future study are also highlighted. Details of the numerical analysis are relegated to Appendix A. II. THE MODEL As per the above discussion, we consider the system of coupled oscillators given by: ¨u=−u+kv +γ˙u+ϵu3+δuv2, ¨v=−v+ku −γ˙v+ϵv3+δvu2.(1) This model is an extension of that in [27], which can be obtained by taking δ= 0. Additionally, it is an extension of the specific case of δ= 3ϵconsidered in [26]. In the linear limit, δ=ϵ= 0, there are two branches of solution eigenfrequencies given by: ω±=√1−γ2/2±√k2−γ2+γ4/4(2) with ω+(ω−) corresponding to symmetric (anti-symmetric) linear modes at γ= 0. Here, we proceed with the understanding that ±ω±are of relevance but we will focus our attention on the positive frequencies hereafter. The two pairs of real (for small γ) eigenfrequencies will collide and give rise to a frequency quartet for γ > γP T,L, where γP T,L satisfies the condition: γ4 P T,L −4γ2 P T,L + 4k2= 0.(3) Thus, for fixed k, the lowest value of γP T,L corresponds to ω=4 √1−k2. Additionally, to this linear analysis, we observe that the nonlinear dynamical equations (1) possess several symmetries that leave them invariant: •u→ −u,v→ −v, •u→ −u,k→ −k,v→v, •u→u,k→ −k,v→ −v, •t→ −t,γ→ −γ. •u→αu,v→αv,ϵ→ϵ/α2,δ→δ/α2. In the limit γ= 0, (1) is a Hamiltonian system, with Hgiven by H=˙u2+ ˙v2+u2+v2 2−ϵ 4(u4+v4)−kuv −δ 2u2v2,(4) and, for the case δ= 3ϵ, dynamical equations (1) are Hamiltonian for any value of γ[26], with a Hamiltonian of the form: 3 H2=pupv+γ 2(upu−vpv) + (1 −γ2 4)uv −k 2(u2+v2)−ϵ(u3v+v3u),(5) In this case, pu= ˙v+γv/2, pv= ˙u−γu/2. The aim of this paper is to identify periodic orbits of frequency ωbof the model (1) (and to compare them also to the results of the DNLS approximation). Toward achieving this aim, Fourier space techniques have been utilized in order to expand the solution in time and to obtain its numerically exact form (up to a prescribed numerical tolerance). Finally, Floquet theory has been used to explore the stability of the pertinent configurations. More details about the numerical methods have been given in Appendix A. An important diagnostic quantity for probing the dependence of the solutions on parameters such as the gain/loss strength γ, or the oscillation frequency ωb, is the energy averaged over a period, defined as: < H >=1 Tb∫Tb 0 H(t) dt, (6) with the Hamiltonian (of the case without gain/loss) given by (4) and Tb= 2π/ωbbeing the oscillation period. In what follows, we will restrict to the values of δ/ϵ ={1,3/2,3}, for which as will be seen below, the solutions and/or dynamical equations possess special properties. In addition, we restrict to γ≥0,0< k < 1and |ϵ|= 1. Unless stated otherwise k=√15/8≈0.48 has been fixed; this value implies γP T,L = 0.5. III. THE ROTATING WAVE APPROXIMATION A. The DNLS dimer: the model, stability equations and analytical results The RWA provides a means of connecting with the extensively analyzed PT-symmetric Schr¨ odinger dimer [8, 10, 15, 28–31]. This link follows a path similar to what has been earlier proposed e.g. in [33, 34]. In particular, the following ansatz is used to approximate the solution of the periodic orbit problem as a roughly monochromatic wavepacket of frequency ωb(for ϕ1,2in what follows we will seek stationary states). u(t)≈ϕ1(t) exp(iωbt) + ϕ∗ 1(t) exp(−iωbt), v(t)≈ϕ2(t) exp(iωbt) + ϕ∗ 2(t) exp(−iωbt).(7) By supposing that ˙ ϕn≪ωbϕnand ¨ ϕn≪ωb˙ ϕn(i.e., ϕvaries slowly on the scale of the oscillation of the actual exact time periodic state), discarding the terms multiplying exp(±3iωbt), the dynamical equations (1) transform into a set of coupled Schr¨ odinger type equations: 2iωb˙ ϕ1= [(ω2 b−1) + 3ϵ|ϕ1|2+ 2δ|ϕ2|2+ iωbγ]ϕ1+ [k+δϕ∗ 1ϕ2]ϕ2, 2iωb˙ ϕ2= [(ω2 b−1) + 3ϵ|ϕ2|2+ 2δ|ϕ1|2−iωbγ]ϕ2+ [k+δϕ∗ 2ϕ1]ϕ1,(8) i.e., forming, under these approximations, a PT-symmetric Schr¨ odinger dimer. The stationary solutions of this dimer can then be used in order to reconstruct via Eq. (A1) the solutions of the RWA to the original PT-symmetric oscillator dimer. These stationary solutions for ϕ1(t)≡y1and ϕ2(t)≡z1satisfy the algebraic conditions Ey1= (p+qz1y∗ 1)z1+ (|y1|2+ 2q|z1|2)y1+ iΓy1, Ez1= (p+qy1z∗ 1)y1+ (|z1|2+ 2q|y1|2)z1−iΓz1,(9) with E=1−ω2 b 3ϵ, p =k 3ϵ, q =δ 3ϵ,Γ = γωb 3ϵ.(10) Notice that when q= 1/2, i.e. δ/ϵ = 3/2, coupling in Eq. (8) resembles that in the Manakov model [32]. If we express y1and z1in polar form: y1=Aexp(iθ1), z1=Bexp(iθ2), φ =θ2−θ1,(11) 4 the stationary equations can be rewritten as EA =pB cos φ+qAB2cos 2φ+A(A2+ 2qB2),(12) EB =pA cos φ+qBA2cos 2φ+B(B2+ 2qA2),(13) −ΓB=Asin φ(p+ 2qAB cos φ),(14) −ΓA=Bsin φ(p+ 2qAB cos φ).(15) In the case γ= 0, there can be symmetric or anti-symmetric solutions, fulfilling A2=B2. Contrary to the δ= 0 setting, where sin φ= 0 only, here we have, apart from this case, the possibility of a phase different than 0 or π, i.e. cos φ= −p/(2qAB). Consequently, we have two pairs of symmetric / anti-symmetric solutions with A=Bat the Hamiltonian limit: A2=E−p 1+3q=1−ω2 b−k 3(ϵ+δ), φ = 0 S0solution (16) A2=E+p 1+3q=1−ω2 b+k 3(ϵ+δ), φ =πA0solution (17) A2=E 1 + q=1−ω2 b 3ϵ+δ, φ = cos−1[−p(1+q) 2qE ]= cos−1[−k(δ+3ϵ) 2δ(1−ω2 b)]Sϕsolution (18) A2=E 1 + q=1−ω2 b 3ϵ+δ, φ =π+ cos−1[p(1+q) 2qE ]=π+ cos−1[k(δ+3ϵ) 2δ(1−ω2 b)]Aϕsolution (19) Recall that A=Bin all the previous cases, i.e. the sign of the anti-symmetric solutions has been introduced into the phase. Apart from the previous solutions, there is an asymmetric solution (AS) whose properties strongly depend on δ/ϵ. This solution is given by: A2= (1 −ω2 b)±√(1 −ω2 b)2−4k2 (1−δ/ϵ)2 6ϵ, B =±k 3(ϵ−δ)A, φ = 0 (π)AS solution.(20) Note that the asymmetric solution exists only if δ=ϵ. When they are equal it is easily checked from RWA equations that there is no asymmetric solution. It is easy to show that at γ= 0 and ϵ > 0,S0solutions exist for ωb< ωS=√1−k,A0solutions exist for ωb< ωA= √1 + kand both Sϕand Aϕsolutions only exist when ωb≤ωϕ+=√1−k(1 + 3ϵ/δ)/2; for ϵ < 0,S0solutions exist for ωb> ωS=√1−k,A0solutions for ωb> ωA=√1 + kand both Sϕand Aϕsolutions only exist when ωb≥ωϕ−= √1 + k(1 + 3ϵ/δ)/2. In addition, asymmetric solutions only exist for ωb< ωAS+ =√1+2k/(1 −δ/ϵ)if ωb<1and for ωb> ωAS−=√1−2k/(1 −δ/ϵ)if ωb>1. Using the identifications ϕ1(t)≡y1and ϕ2(t)≡z1introduced after (8), the averaged energy within the RWA can be written as: < H >= (1 + ω2 b)(|y1|2+|z1|2)−2kRe(y1z∗ 1)−3ϵ 2(|y1|4+|z1|4)−δ[Re(y2 1z∗2 1)+2|y1|2|z1|2](21) and, by making use of (11), the average energy for each of the previous solutions at γ= 0 is given by the following expressions: < H > =ω4 S+ 2ω2 Sω2 b−3ω4 b 3(ϵ+δ),S0solution < H > =ω4 A+ 2ω2 Aω2 b−3ω4 b 3(ϵ+δ),A0solution < H > =1+2ω2 b−3ω4 b 3ϵ+δ+k2 2δ,Sϕand Aϕsolutions < H > =1+2ω2 b−3ω4 b 6ϵ+k2 3(δ−ϵ).AS solution (22) 5 Notice that the average energy of both Sϕand Aϕare the same for every δand that also coincide with that of the AS solution for δ= 3ϵ. When γ= 0 only symmetric and anti-symmetric solutions can exist and Eqs. (12)-(15) can be simplified as a quartic equation for A2: 4 ∑ j=0 PjA2j= 0 (23) with P0= (Γ2+E2)(Γ2+E2−p2), P1= 2E[(1 + q)p2−2(1 + 2q)(Γ2+E2)], P2= 4(1 + 2q)2E2+ 2(1 + 3q)(1 + q)(Γ2+E2)−(1 + q)2p2, P3=−4E(1 + q)(1 + 2q)(1 + 3q), P4= (1 + 3q)2(1 + q)2, whereas the phase fulfills the equation: tan φ=−Γ E−(1 + q)A2.(24) Just as one could give an expression for Awithout involving ϕ, similarly by eliminating A, one finds that ϕmust satisfy the constraint Eq sin(2ϕ)±p(1 + q) sin(ϕ) + Γ[1 + q+ 2qcos2(ϕ)] = 0 , B =±A . (25) Notice that there is a phase degeneracy that must be removed by applying, e.g., Eq. (15) together with the previous one. We now turn to the linear stability of different solutions within the RWA. The spectral analysis of the symmetric and antisymmetric solutions can be obtained by considering small perturbations [of order O(ε), with 0< ε ≪1] of the stationary solutions. The stability can be determined by substituting the ansatz below into (8) and then solving the ensuing [to O(ε)] eigenvalue problem: ϕ1(t) = y1+ε(a1e−iθt/Tb+b∗ 1eiθ∗t/Tb), ϕ2(t) = z1+ε(a2e−iθt/Tb+b∗ 2eiθ∗t/Tb),(26) with Tb= 2π/ωbbeing the orbit’s period and θbeing the Floquet exponent (FE). The FEs can be expressed as: θ=π ω2 b iΩ (27) with Ωbeing the eigenfrequencies of the stability matrix M, which is defined as Ω(a1, a2, b∗ 1, b∗ 2)T=M(a1, a2, b∗ 1, b∗ 2)T. In the case of symmetric and anti-symmetric solutions, the matrix can be written as: M=   M1M2M3M4 M2M∗ 1M4M∗ 3 −M∗ 3−M4−M∗ 1−M2 −M4−M3−M2−M1   (28) with the elements being M1= (ω2 b−1) + 2(3ϵ+δ)A2+ iωbγ, (29) M2= 4δA2cos φ+k, (30) M3= [3ϵexp(−iφ) + δexp(iφ)]A2,(31) M4= 2δA2.(32) Thus, the non-zero eigenvalues λcan be expressed in terms of A2and φ, which must be determined by solving Eqs. (23)-(24): Ω2/2 = [δ2(1−16 cos2φ)−6ϵδ(5−2 cos2φ)−27ϵ2]A4−[8kδ cos φ−4(ω2 b−1)(3ϵ+δ)]A2−[(ω2 b−1)2+k2−γ2ω2 b],(33) 6 B. Numerical analysis of symmetric and anti-symmetric solutions We show below the properties of the A0,S0,Aϕand Sϕsolutions in the cases δ=ϵ,δ= 3ϵ/2and δ= 3ϵfor both soft (ϵ= +1) and hard (ϵ=−1) potentials. A summary of the existence and stability regions is displayed in Fig. 1, where the panels depict the γ-ωbplanes. Notice that although A0,S0,Aϕand Sϕsolutions are, strictly speaking, defined only at γ= 0, we will use this notation for solutions at γ= 0 that are obtained by continuation from the Hamiltonian (γ= 0) limit. Prior to starting the analysis for arbitrary γ, we will briefly show the properties of the asymmetric solutions at γ= 0. As explained above, we will choose k=√15/8. Notice that for this parameter value, ω2 AS+ <0if δ= 3ϵ/2>0and consequently, there is no asymmetric solution for this regime. However, there are asymmetric solutions if δ= 3ϵ/2<0and ωb> ωAS−≈1.7136. In fact, at ωb=ωAS−there is a pitchfork bifurcation, as the S0solution is unstable for ωb< ωAS−and becomes stable past the bifurcation point, where a pair of branches corresponding to unstable asymmetric solutions emerge. For δ= 3ϵ, the situation is similar in the hard case in what regards the existence of solutions (now ωAS−=ωA≈1.2182); for the soft case, the asymmetric solution exists for ωb≤ωAS+ =ωS≈0.7182 and bifurcates from the A0solution. Notice that the A0(for the soft case) and the S0(for the hard case) are all stable; furthermore, the AS solution appears to be marginally stable and highly degenerate as all the eigenfrequencies Ωare equal to zero; recall also the special, completely integrable nature of this special limit. We analyze now the properties of the soft potential when γ= 0. In the δ=ϵcase, there are two main regions: in region I, only A0solutions exist, as S0solutions bifurcate from the left arm of the γL(ωb)curve (i.e. ω+), which corresponds to the symmetric linear modes; at the right of region I, no solutions are found because A0solutions bifurcate from the right arm (i.e. ω−) of the linear dispersion relation. In this soft case of ϵ > 0, the bifurcations occur to the left of γL(ωb), whereas in the hard case of ϵ < 0, they arise to the right of γL(ωb). It is easy to show that from (2), γLis given by: γL(ωb) = √(k2−1) + 2ω2 b−ω4 b ωb .(34) Consequently, region I is bounded between ωb=ωS≈0.7182,ωb=ωA≈1.2182 and γ=γP T,L = 0.5. In region II, both A0and S0solutions exist, and experience the PT phase transition at the curve designated as γPT(ωb). Notice also that all the solutions existing in both regions I and II are stable. In addition, Aϕand Sϕsolutions can only be found for ωb< ωϕ+, but their existence range is quite small as ωϕ+≈0.1782 and only exist for γ < 0.1. For δ= 3ϵ/2, both regions I and II have the same properties as before. In addition, region III is included, which is below the curve γ1(ωb). In that region, all four solutions exist and are stable except for A0. This solution becomes stable only nearby i.e. between the curves γ2(ωb)and γ1(ωb). This small stability region can be observed between red and green curves of the inset of the corresponding panel (notice that this phenomenon was also observed in the δ=ϵcase, but was not showcased due to the very small range of existence of Sϕand Aϕsolutions therein). In region II only two solutions exist; for ωb< ωϕ+≈0.5233, i.e. above the curve γ1(ωb),S0and Sϕsolutions coexist and collide/disappear at γPT(ωb), whereas for ωb> ωϕ+, the coexisting solutions are A0and S0. As a side comment, the reason for the existence of Sϕfor ωb< ωϕ+and of A0for ωb> ωϕ+within Region II has to do with the fact that these solutions effectively “morph” from one to the other (smoothly) as this frequency is crossed. For δ= 3ϵ, the scenario is similar to the last one, except for two points: first, the γ1curve finishes at ωb=ωSand encompasses an accordingly broader region III; and second, the solutions in region III, A0and Aϕ, are stable for any value of γand ωb. We focus now on the hard potential (i.e. ϵ=−1) properties. In all the cases, we can find the region I, with the same properties as in the soft case (although now it is the S0solution that exists and the A0that bifurcates into existence beyond the boundary of the region). In addition, region II is present in every case, enclosed between curves γPT(ωb)and γ1(ωb); there are two solutions therein: the A0solution and the S0for ωb< ωϕ−and the A0and Aϕfor ωb> ωϕ−. Under the curve γ1(ωb), whose minimum value takes place at ωb=ωϕ−(so that ωϕ−=ωAfor δ= 3ϵ), the four kinds of solutions coexist, so that S0and Sϕcollide and disappear at this line. As in the soft case, there is a “morphing” from the S0to Aϕsolutions when the frequency ωϕ−is crossed. Thus, the most significant difference between the three considered regimes lies in the existence of region IV and curve γ2(ωb). Region III is characterized by the fact that all the solutions exist (as mentioned above) and are stable. However, below the curve γ2(ωb)(i.e. in region IV), solution S0becomes unstable. Notice that for δ=ϵthis region exists for every ωb> ωϕ−≈1.4029. However, if δ= 3ϵ/2, region IV is shrunk to the range ωϕ−≈1.3138 < ωb.1.71. Finally, region IV has totally vanished at δ= 3ϵ. IV. ANALYSIS OF THE OSCILLATOR DIMER In this section, we complete the description of the system by returning to the original oscillator system and analyzing its exact periodic orbits (that up to now we had only approximated using the RWA). This is done by numerically solving in the Fourier space the dynamical equations set (1) [cf. Appendix A]. That is, we express the solution in the form: 7 ϵ= 1,δ=ϵ ϵ =−1,δ=ϵ 0.2 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 2 2.5 γ ωb I II γPT γL 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 γ ωb I II IV γPT γL γ1 γ2 III ϵ= 1,δ= 3ϵ/2ϵ=−1,δ= 3ϵ/2 0.2 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 2 2.5 3 γ ωb I II γPT γL γ1 0.2 0.3 0.4 0.5 0 0.05 0.1 γ1 γ2 III 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 γ ωb I II III IV γPT γLγ1 γ2 ϵ= 1,δ= 3ϵ ϵ =−1,δ= 3ϵ 0.2 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 2 2.5 3 γ ωb I II III γPT γL γ1 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 γ ωb I II III γPT γL γ1 FIG. 1: γ-ωbplane for k=√15/8. Details on the meaning of each curve and region can be found in the text. The linear limit of the oscillator system is denoted by γL, while the upper PT -symmetric threshold of solution existence is denoted by γPT . An additional delimiter of the existence of further solutions Aϕand Sϕis also given by γ1. The existence regions of the different solutions are encompassed by these curves both in the soft ϵ= 1 case (left panels) and in the hard ϵ=−1case (right panels). u(t) = ∑ n ynexp(inωbt), v(t) = ∑ n znexp(inωbt).(35) We have considered the same cases as in Section III, namely, δ/ϵ equal to 1, 3/2and 3, with ϵ=±1. Prior to showing the results, we want to remark that the Fourier coefficients of S0and A0solutions (due to their symmetry) have the following property: yn=z∗ n(S0), yn=−z∗ n(A0).(36) In what follows, we will first show the properties of the solutions at the Hamiltonian limit γ= 0. Afterwards, we will be focusing in the different cases of δ/ϵ > 0for γ= 0. In most cases, results will be compared with the previously found results for the RWA. 8 A. Solutions for γ= 0 We start by analyzing the modes that can be expressed analytically at the γ= 0 limit. In fact, these can be expressed in terms of Jacobi elliptic functions. If ϵ > 0(soft potential), the solutions are of the form: u(t) = Asn[β t;m], v(t) = ±Asn[β t;m],(37) with A=β√2m ϵ+δ, β2=1∓k 1 + m, ωb=πβ 2K(m).(38) with the upper (lower) sign corresponding to the S0(A0) solution, K(m)the complete elliptic integral of the first kind with modulus m[37], and 0< m < 1. As K(m)> π/2, it is easy to deduce that ωb< ωS=√1−kfor the S0solution and ωb< ωA=√1 + kfor the A0solution, as within the RWA. If ϵ < 0(hard potential), these modes can be expressed as: u(t) = Acn[β t;m], v(t) = ±Acn[β t;m],(39) with A=β√−2m ϵ+δ, β2=1∓k 1−2m, ωb=πβ 2K(m),(40) where 0< m < 1/2. Similar to the soft case, ωb> ωSfor the S0solution and ωb> ωAfor the A0solution. Notice that for these solutions to exist δ < −ϵ. At δ=ϵand for a hard potential, the Sϕand Aϕsolutions are given by: u(t) = Asn(βt;m) + B√mcn(βt;m), v(t) = ±[A√msn(βt;m)−B√mcn(βt;m)] ,(41) provided A=√(3m−2)β2+ 2 4ϵ, B =√1−(2m+ 1)β2 2ϵ, β2= 2k/m, ωb=πβ 2K(m).(42) The AS solution can be analytically expressed whenever δ= 3ϵ. If ϵ > 0, it is given by: u(t) = A+sn[β+t;m+] + A−sn[β−t;m−], v(t) = A+sn[β+t;m+]−A−sn[β−t;m−],(43) with A±=β±√2m± ϵ, β2 ±=1∓k 1 + m± , ωb=πβ+ 2K(m+)=πβ− 2K(m−),(44) whereas if ϵ < 0, the AS solution is: u(t) = A+cn[β+t;m+] + A−cn[β−t;m−], v(t) = A+cn[β+t;m+]−A−cn[β−t;m−],(45) with A±=β±√−2m± ϵ, β2 ±=1∓k 1−2m± , ωb=πβ+ 2K(m+)=πβ− 2K(m−), . (46) It can be numerically observed that A0,S0and AS solutions exist for every δwhereas Aϕand Sϕdo not exist for δ≤3ϵ/2in the ϵ= +1 case. On the contrary, a new solution denoted as A3exists for the soft potential and all of the considered values of δ; this new solution, which was not found in the δ= 0 case, is characterized by a high increase of the third harmonic in the Fourier 9 series, and, consequently, cannot be predicted by the RWA. The existence of this new solution can be caused by the hybridization of the S0mode with frequency ωband the A0mode with frequency 3ωb; this symmetry breaking effect could happen whenever ωb< ωA/3≈0.4061. Notice that the A3mode bifurcates from the S0mode at ω3, which exactly coincides with ωA/3when δ= 3ϵ, but is smaller than this when δ < 3ϵ(e.g. for δ= 3ϵ/2,ω3≈0.384 whereas ω3≈0.365 for δ=ϵ.). There is no stability change at this bifurcation. The asymmetric (AS) solution preserves the properties of the RWA. That is, it bifurcates from the A0solution in soft potentials and from the S0solution for hard potentials. The AS solution does not exist for ϵ=δand for δ= 3ϵ/2>0. Besides, all the Floquet exponents are θ= 0 (or, equivalently, the Floquet multipliers are +1) for δ= 3ϵ. For δ= 3ϵ/2<0, the AS solution is unstable, as in the RWA. In addition, for δ= 3ϵ, the AS, Sϕand Aϕsolutions bifurcate from the A0solution at ωb=ωSif ϵ= +1, with the Sϕand Aϕsolutions being stable and the A0stable (unstable) for ωb> ωS(ωb< ωS); if ϵ=−1, the AS, Sϕand Aϕsolutions bifurcate from the S0mode at ωb=ωA, with the S0solution being stable and the Sϕunstable, whereas the Aϕis marginally stable as are the AS solutions (all the Floquet exponents are zero). In the δ= 3ϵ/2<0case, the Sϕand Aϕsolutions, which are stable, bifurcate from the S0solution at ωb≈1.306 which is close to ωϕ−; the S0solution is stable (unstable) for ωbsmaller (higher) than the bifurcation point. This latter bifurcation also occurs for δ=ϵ=−1, taking place in this case at ωb≈1.386. In addition, the AS solution (which is unstable) bifurcates from the S0solution (which changes its stability) at ωb≈1.708, a value which is close to (but not exactly at) ωAS−. We must also mention that for the case analyzed in [27], i.e. δ= 0, the AS solution bifurcates from the S0(A0) solution in the soft (hard) potential. This situation is reversed in the present observations for sufficiently large δ= 0, which suggests the existence of a critical point. Finally, as within the RWA, the energy coincides for the AS, Sϕand Aϕsolutions when δ= 3ϵ. All the previous properties are summarized in Fig. 2 where the Hamiltonian energy is depicted versus ωband compared with the averaged Hamiltonian for the RWA. This figure is complemented by Figs. 3 and 4 where the time evolution of the different solutions are displayed. Importantly, we should point out here that it is evident that the approximations involved in the RWA become demonstrably less accurate especially in the soft nonlinearity case and particularly as the frequency ωdecreases away from the linear limit (and hence nonlinear terms become more significant). Nevertheless, the qualitative agreement of the features of Fig. 2 is still fairly satisfactory for the regime of parameters considered herein. On the other hand, for the hard nonlinearity case, the agreement seems to be even quantitatively accurate for the frequency range considered. B. δ=ϵcase: existence of exact solutions One of the main features of this case is the existence of two exact periodic solutions to (1): u(t) = Asin(ωbt), v(t) = ±Acos(ωbt)(47) fulfilling that: k=∓γωb, A =√1−ω2 b ϵ,(48) with the upper sign corresponding to the symmetric solution and the lower one to the anti-symmetric solution. It is important to note that for a given k, the frequency is proportional to 1/γ. Thus, the two solutions collide as γ→ ∞, when ωb→0. That is, contrary to the “standard” model of δ= 0, since for the case considered herein there exist nonlinear solutions for all values of γ that are not subject to the relevant transition [38]. This solution can actually be cast as yn=zn= 0 ∀|n|>1and ϕ=±π/2. If we fix the value of ϵ, it is clear from Eq. (48) that the properties of the solutions only depend on two parameters, as k=k(ωb, γ). That is, contrary to what we have discussed so far, here we do not fix kand vary γand ωb, but rather than varying γand ωb, we fix a value of kassociated with them through Eq. (48). Thus, we will consider the effect on the stability of varying parameters ωband γin the case ϵ= 1 (soft potential) and ϵ=−1(hard potential). Notice also that given the restrictions formulated in (48) and the symmetry properties of the dynamical equations, the Floquet spectrum for a given set of parameters is the same for both solutions. Fig. 5 shows the stability/instability regions for these solutions. Shaded areas correspond to stable solutions. The black line therein indicates the locus in the γ-ωbplane where k=√15/8(i.e., the value used for other results in the present work). From this line, it can be deduced that solutions with ϕ=±π/2when δ=ϵare stable in the range ωb∈[0.8535,1] if ϵ= 1 and in ωb∈[1,1.029] ∪[1.2206,∞)if ϵ=−1. The averaged energy is, for both solutions, < H >=1+2ω2−3ω4 3ϵwhich, for ϵ= 1 has a maximum at ωb= 3−1/2≈0.5774; for ϵ=−1, this function is monotonically decreasing. It is worth mentioning that for γ= 0, where the averaged energy coincides with the Hamiltonian, there is a stability change at ω= 3−1/2, the value at which ∂H/∂ω changes its slope. This correlation between energy maximum and stability changes, which resembles the Vakhitov-Kolokolov criterion for NLS systems, 16 δ= 3ϵ.ϵ= +1 δ= 3ϵ.ϵ=−1 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 2 γ ωb γPT γL γa γs γ3 γ1 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 γ ωb γPT γLγ1 γ3 γ2 FIG. 11: Planes with curves separating regions of solutions that share the same properties when δ= 3ϵand s=√15/8(see text). δ= 3ϵ,ϵ= 1,ωb= 0.6δ= 3ϵ,ϵ= 1,ωb= 0.3δ= 3ϵ,ϵ=−1,ωb= 2 0 0.2 0.4 0.6 0.8 1 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 S0 A0 Sφ Aφ H2 γ 0 0.5 1 1.5 2 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 S0 A0 Sφ Aφ A3 − A3 + H2 γ 0 0.2 0.4 0.6 0.8 −4 −3 −2 −1 0 1 2 3 4 S0 A0 Sφ Aφ H2 γ FIG. 12: Dependence of H2on γfor each mode at different frequencies and ϵfor the case δ= 3ϵ. work of [26], only a specific value of the cross interaction was considered (δ= 3ϵ), revealing remarkably the Hamiltonian nature of the model, and then restricting consideration to its DNLS analogue. Here, we have extended considerations to three relevant cases, namely δ=ϵ,δ= 3ϵ/2and δ= 3ϵ, exploring how the existence, nonlinear bifurcation and even dynamical trends develop as we move from weaker to stronger cross-interaction between the nonlinear oscillators. Importantly, the relevant pictures were developed not only for the rotating wave approximation model of the DNLS form, but also for the full model of the coupled oscillators. Generally, the two cases, namely the monochromatic approximation and the full system were very similar, except for the highly nonlinear regime, especially in the soft nonlinearity case. Numerous important features were identified along the way including, e.g., new families of solutions at relative phase angles other than 0and π(introduced by the cross-interaction between oscillators), as well as solutions tractable solely in a numerical form from the four principal families explored. Yet another feature was the existence in the oscillator system of families of solutions not only in the γ= 0 but even in the γ= 0 case in explicit form; one such pair of families appears to “defy” the PT phase transition (in the δ=ϵcase), existing for all values of the gain/loss parameter γ. Finally, the instabilities identified in the analysis were monitored in the full dynamics of the system, revealing the possibility of either indefinite growth or that of bounded quasi-periodic oscillations, as the pertinent dynamical outcome. There are numerous questions that naturally emerge as a result of the present work. Among the most immediate ones, it is worthwhile to extend considerations to the case of, e.g., three oscillators and perhaps even to that of four such, forming effectively a two-dimensional plaquette and a building block for the consideration of higher dimensional systems, in the spirit also of [9]. Furthermore, here only the case of cubic nonlinearities has been explored, but it might be also of interest, as another prototypical nonlinear system to examine the case of quadratic nonlinearities and how their nonlinear states are “deformed” in the presence of gain and loss. At a perhaps deeper level, however, there are also some intriguing questions that we feel are raised. For one, an apparently PT -symmetric and viewed as a gain-loss bearing system at δ= 3ϵis found to be Hamiltonian. This raises the natural yet difficult question: can we discern such a potential Hamiltonian nature and classify a system as Hamiltonian (and not PT ) possibly through an appropriate (to be identified) transformation? If so, what is the relevant criterion and how can we exclude the presence of a yet-unknown transform that may convert a system classified as PT into one which is genuinely 17 δ=ϵ,ϵ= 1,ωb= 0.5,γ= 0.6(S0)δ=ϵ,ϵ= 1,ωb= 0.8,γ= 0.5(A0) 0 20 40 60 80 100 −1 0 1 2 3 4 5 6 u(t), v(t) t 0 5 10 15 20 25 30 −6 −5 −4 −3 −2 −1 0 1 u(t), v(t) t δ= 3ϵ,ϵ= 1,ωb= 0.5,γ= 0.05 (A0)δ= 3ϵ,ϵ= 1,ωb= 0.7,γ= 0.02 (Aϕ) 0 20 40 60 80 100 −20 −15 −10 −5 0 5 10 15 20 u(t), v(t) t 0 1000 2000 3000 4000 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 u(t), v(t) t FIG. 13: Evolution of unstable solutions for the soft potential. Three examples provide the different combination examples where the oscillator amplitudes may grow indefinitely, while the fourth example presents a bounded apparently quasi-periodic scenario. Hamiltonian in a different set of variables ? Potential progress along these veins will be reported in future work. P.G.K. acknowledges support from the National Science Foundation under grants CMMI-1000337, DMS-1312856, from FP7-People under grant IRSES-606096, from the US-AFOSR under grant FA9550-12-10332 and from the Binational (USIsrael) Science Foundation through grant 2010239. This work was supported in part by the U.S. Department of Energy. A.K. acknowledges financial support from Dept. of Atomic Energy, Govt. of India through a Raja Ramanna Fellowship. P.G.K. also acknowledges useful discussions with Igor Barashenkov. APPENDIX A: NUMERICAL ANALYSIS OF PERIODIC ORBITS In order to calculate periodic orbits, we make use of a Fourier space implementation of the dynamical equations and continuations in frequency or gain/loss parameter are performed via a path-following (Newton-Raphson) method. Fourier space methods are based on the fact that the solutions are Tb-periodic; for a detailed explanation of these methods, the reader is referred to Refs. [35, 36]. The method has the advantage, among others, of providing an explicit, analytical form of the Jacobian. Thus, the solution for the two nodes can be expressed in terms of a truncated Fourier series expansion: u(t) = nm ∑ n=−nm ynexp(inωbt), v(t) = nm ∑ n=−nm znexp(inωbt),(A1) with nmbeing the maximum of the absolute value of the running index kin our Galerkin truncation of the full Fourier series solution. In the numerics, nmhas been chosen as 21. After the introduction of (A1), the dynamical equations (1) yield a set of 2×(2nm+ 1) nonlinear, coupled algebraic equations: 18 δ=ϵ,ϵ=−1,ωb= 1,γ= 0.16 (S0)δ=ϵ,ϵ=−1,ωb= 2,γ= 0.46 (S0) 0 500 1000 1500 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 u(t), v(t) t 0 200 400 600 800 1000 −3 −2 −1 0 1 2 3 u(t), v(t) t δ= 3ϵ,ϵ= 1,ωb= 2,γ= 0.2(Aϕ)δ= 3ϵ,ϵ= 1,ωb= 2,γ= 0.3(Sϕ) 0 100 200 300 −15 −10 −5 0 5 10 15 u(t), v(t) t 0 20 40 60 80 100 120 −20 −15 −10 −5 0 5 10 15 20 u(t), v(t) t FIG. 14: Evolution of unstable solutions for the hard potential. The top panels feature examples of quasi-periodic oscillations, while the bottom panels illustrate indefinite growth of one of the oscillators coupled with a decaying oscillation of the other (possibly very slowly, as in the case of the bottom left panel). Fn,1≡ −ω2 bn2yn−iγωbnyn+Fn[V′(u, v)] −kzn= 0,(A2) Fn,2≡ −ω2 bn2zn+ iγωbnzn+Fn[V′(v, u)] −kyn= 0,(A3) with V′(u1, u2) = u1−ϵu3 1−δu1u2 2. Here, Fndenotes the Discrete Fourier Transform: Fn[V′(u)] = 1 N nm ∑ q=−nm V′(nm ∑ p=−nm ypexp [i2πpq N])exp [−i2πnq N],(A4) with N= 2nm+ 1. The procedure for Fn(v)is similar to the previous case. As u(t)and v(t)must be real functions, it implies that y−n=y∗ n, z−n=z∗ n. In order to study the spectral stability of periodic orbits, we introduce a small perturbation {ξ1, ξ2}to a given solution {u0, v0} of Eqs. (1) according to u=u0+ξ1,v=v0+ξ2. Then, the equations satisfied to first order in ξnread: ¨ ξ1= (3ϵu2 0+δv2 0−1)ξ1+γ˙ ξ1+ (k+ 2δu0v0)ξ2, ¨ ξ2= (3ϵv2 0+δu2 0−1)ξ2−γ˙ ξ2+ (k+ 2δu0v0)ξ1,(A5) or, in a more compact form: N({u(t), v(t)})ξ= 0 , where N({u(t), v(t)})is the relevant linearization operator. In order to study the spectral (linear) stability analysis of the relevant solution, a Floquet analysis can be performed if there exists Tb∈Rso that the map {u(0), v(0)}→{u(Tb), v(Tb)}has a fixed point (which constitutes a periodic orbit of the original system). Then, the stability properties are given by the spectrum of the Floquet operator M(whose matrix representation is the monodromy) defined as: ({ξn(Tb)} {˙ ξn(Tb)})=M({ξn(0)} {˙ ξn(0)}).(A6) 19 The 4×4monodromy eigenvalues Λ = exp(iθ)are dubbed the Floquet multipliers and θare denoted as Floquet exponents (FEs). This operator is real, which implies that there is always a pair of multipliers at 1(corresponding to the so-called phase and growth modes) and that the eigenvalues come in pairs {Λ,Λ∗}. As a consequence, due to the “simplicity” of the FE structure (one pair always at 1and one additional pair) there cannot exist Hopf bifurcations in the dimer, as such bifurcations would imply the collision of two pairs of multipliers and the consequent formation of a quadruplet of eigenvalues which is impossible here. Nevertheless, in the present problem, the motion of the pair of multipliers can lead to an instability through exiting (through 1 or −1) on the real line leading to one multiplier (in absolute value) larger than 1and one smaller than 1. [1] C. M. Bender, Rep. Prog. Phys. 70, 947 (2007). [2] See special issues: H. Geyer, D. Heiss, and M. Znojil, Eds., J. Phys. A: Math. Gen. 39,Special Issue Dedicated to the Physics of NonHermitian Operators (PHHQP IV) (University of Stellenbosch, South Africa, 2005) (2006); A. Fring, H. Jones, and M. Znojil, Eds., J. Math. Phys. A: Math Theor. 41,Papers Dedicated to the Subject of the 6th International Workshop on Pseudo-Hermitian Hamiltonians in Quantum Physics (PHHQPVI) (City University London, UK, 2007) (2008); C.M. Bender, A. Fring, U. G¨ unther, and H. Jones, Eds., Special Issue: Quantum Physics with non-Hermitian Operators, J. Math. Phys. A: Math Theor. 41, No. 44 (2012). [3] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, PT symmetric periodic optical potentials, Int. J. Theor. Phys. 50, 1019 (2011). [4] A. Ruschhaupt, F. Delgado, and J. G. Muga, J. Phys. A: Math. Gen. 38, L171 (2005). [5] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, Phys. Rev. Lett. 100, 103904 (2008); S. Klaiman, U. G¨ unther, and N. Moiseyev, ibid.101, 080402 (2008); O. Bendix, R. Fleischmann, T. Kottos, and B. Shapiro, ibid.103, 030402 (2009); S. Longhi, ibid.103, 123601 (2009); Phys. Rev. B 80, 235102 (2009); Phys. Rev. A 81, 022102 (2010). [6] A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides, Phys. Rev. Lett. 103, 093902 (2009); C. E. R¨ uter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, Nature Phys. 6, 192 (2010); A. Regensburger, C. Bersch, M.-A. Miri, G. Onishchukov, D. N. Christodoulides, and U. Peschel, Nature 488, 167 (2012). [7] S.V. Dmitriev, A.A. Sukhorukov, and Yu.S. Kivshar, Opt. Lett. 35, 2976 (2010). [8] K. Li and P.G. Kevrekidis, Phys. Rev. E 83, 066608 (2011). [9] K. Li, P.G. Kevrekidis, B.A. Malomed, and U. G¨ unther, J. Phys. A Math. Theor. 45, 444021 (2012). [10] H. Ramezani, T. Kottos, R. El-Ganainy, and D.N. Christodoulides, Phys. Rev. A 82, 043803 (2010). [11] S.V. Suchkov, B.A. Malomed, S.V. Dmitriev, and Yu.S. Kivshar, Phys. Rev. E 84, 046609 (2011). [12] A.A. Sukhorukov, Z. Xu, and Yu.S. Kivshar, Phys. Rev. A 82, 043818 (2010). [13] D.A. Zezyulin and V.V. Konotop, Phys. Rev. Lett. 108, 213906 (2012). [14] A.E. Miroshnichenko, B.A. Malomed, and Yu.S. Kivshar, Phys. Rev. A 84, 012123 (2011) [15] H. Cartarius and G. Wunner, Phys. Rev. A 86, 013612 (2012) [16] V.V. Konotop, D.E. Pelinovsky, and D.A. Zezyulin, EPL 100, 56006 (2012). [17] A.A. Sukhorukov, S.V. Dmitriev, S.V. Suchkov, and Yu.S. Kivshar, Opt. Lett. 37, 2148 (2012). [18] M.C. Zheng, D.N. Christodoulides, R. Fleischmann, and T. Kottos, Phys. Rev. A 82, 010103(R) (2010). [19] C. M. Bender, B. Berntson, D. Parker, and E. Samuel Am. J. Phys. 81, 173 (2013). [20] J. Schindler, A. Li, M.C. Zheng, F.M. Ellis, and T. Kottos, Phys. Rev. A 84, 040101 (2011). [21] J. Schindler, Z. Lin, J. M. Lee, H. Ramezani, F. M. Ellis, and T. Kottos, J. Phys. A: Math. Theor. 45, 444029 (2012). [22] B. Peng, S.K. Ozdemir, F. Lei, F. Monifi, M. Gianfreda, G.L. Long, S. Fan, F. Nori, C.M. Bender, L. Yang, Nature Physics 10 (2014) 394. [23] N. Lazarides, G.P. Tsironis, Phys. Rev. Lett. 110, 053901 (2013). G.P. Tsironis, N. Lazarides Appl. Phys. A 115, 449 (2014). [24] N. Bender, S. Factor, J. D. Bodyfelt, H. Ramezani, D. N. Christodoulides, F. M. Ellis, and T. Kottos Phys. Rev. Lett. 110, 234101 (2013). [25] C.M. Bender, M. Gianfreda, S.K. ¨ Ozdemir, B. Peng and L. Yang, Phys. Rev. A 88, 062111 (2013). [26] I.V. Barashenkov and M. Gianfreda, J. Phys. A: Math. Theor. 47 282001 (2014). [27] J. Cuevas, P.G. Kevrekidis, A. Saxena, and A. Khare, Phys. Rev. A 88, 032108 (2013). [28] E.-M. Graefe, J. Phys. A: Math. Theor. 45, 444015 (2012). [29] J. Pickton, H. Susanto, Phys. Rev. A 88, 063840 (2013). [30] A.S. Rodrigues, K. Li, V. Achilleos, P.G. Kevrekidis, D.J. Frantzeskakis, and C.M. Bender, Romanian Rep. Phys. 65, 5 (2013). [31] I.V. Barashenkov, G.S. Jackson and S. Flach, Phys. Rev. A 88, 053817 (2013). [32] S.V. Manakov, Sov. Phys. JETP 38, 248 (1974). [33] Yu.S. Kivshar and M. Peyrard. Phys. Rev. A 46, 3198 (1992). [34] K.W. Sandusky, J.B. Page, and K.E. Schmidt. Phys. Rev. B 46, 6161 (1992). [35] J.F.R. Archilla, R.S. MacKay, and J.L. Mar´ ın, Physica D134, 406 (1999). [36] J. Cuevas, J.F.R. Archilla, and F.R. Romero, J. Phys. A: Math. Theor. 44, 035102 (2011). [37] Here, we use the definition K(m) = ∫1 0 dx √(1−x2)(1−mx2). [38] We acknowledge here that Igor Barashenkov in his recent talk at the SIAM conference on Nonlinear Waves and Coherent Structures (Cambridge, August 2014) reported an apparently similar feature as part of ongoing work with Dimitry Pelinovsky [39] Notice that, contrary to the Hamiltonian H, it does not need to be averaged because H2is a constant of motion for δ= 3ϵ