PT-symmetric dimer of coupled nonlinear oscillators

PHYSICAL REVIEW A 88, 032108 (2013) PT -symmetric dimer of coupled nonlinear oscillators Jes´ us Cuevas Nonlinear Physics Group, Departamento de F´ ısica Aplicada I, Universidad de Sevilla. Escuela Polit´ ecnica Superior, C/ Virgen de ´ Africa, 7, 41011-Sevilla, Spain Panayotis G. Kevrekidis Department of Mathematics and Statistics, University of Massachusetts, Amherst, Massachusetts 01003-9305, USA Avadh Saxena Center for Nonlinear Studies and Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA Avinash Khare Indian Institute of Science Education and Research (IISER), Pune 411008, India (Received 23 July 2013; published 16 September 2013) We provide a systematic analysis of a prototypical nonlinear oscillator system respecting PT symmetry i.e., one of them has gain and the other an equal and opposite amount of loss. Starting from the linear limit of the system, we extend considerations to the nonlinear case for both soft and hard cubic nonlinearities identifying symmetric and antisymmetric breather solutions, as well as symmetry-breaking variants thereof. We propose a reduction of the system to a Schr¨ odinger-type PT -symmetric dimer, whose detailed earlier understanding can explain many of the phenomena observed herein, including the PT phase transition. Nevertheless, there are also significant parametric as well as phenomenological potential differences between the two models and we discuss where these arise and where they are most pronounced. Finally, we also provide examples of the evolution dynamics of the different states in their regimes of instability. DOI: 10.1103/PhysRevA.88.032108 PACS number(s): 11.30.Er, 63.20.Pw, 63.20.Ry I. INTRODUCTION The topic of parity-time (PT ) symmetry and its relevance to physical applications, on the one hand, as well as its mathematical structure, on the other, have drawn considerable attention from both the physics and the mathematics communities. Originally, this theme was proposed by C. Bender and coworkers as an additional possibility for operators associated with real measurable quantities within linear quantum mechanics [1–3]. However, one of the major milestones (and a principal thrust of recent activity) regarding the physical and experimental realizability of the corresponding Hamiltonians stemmed from progress in optics both at the theoretical [4,5] and at the experimental [6,7] levels. In particular, the realization that, in optics, the ubiquitous loss can be counteracted by an overwhelming gain in order to create a PT -symmetric setup, e.g., in a waveguide dimer [7], paved the way for numerous developments, especially so at the level of nonlinear systems, as several researchers studied nonlinear stationary states, stability, and dynamics of few site configurations [8–16]as well as of infinite lattices [17–19]. Interestingly, most of this nonlinear activity has been centered around Schr¨ odinger-type systems and for good reason, since the original proposal by Bender involved quantummechanical settings, where this is natural and, in addition, the optics proposal was placed chiefly on a similar footing (i.e., the Schr¨ odinger model as paraxial approximation to the Maxwell equations). Nevertheless, there have been a few notable exceptions where nonlinear oscillator models (involving second-order differential equations in time) have been considered. Perhaps the most relevant example, also for the considerations presented herein, has involved the realization of PT -symmetric dimers in the context of electrical circuits; for the first work in this context, see Ref. [20], while that and follow-up activity has recently been summarized in a review [21]. Chiefly, the experimental considerations of these works focused on the linear variant of the gain-loss oscillator system. More recently, nonlinear variants of PT -symmetric dimers in the form of a chain have been proposed in the context of magnetic metamaterials and in particular for systems consisting of split-ring resonators [22]. The latter setting, while nonlinear, is also far more complex (involving external drive and nonlinear couplings between adjacent sites) and, hence, was tackled for the nonlinear model chiefly at the level of direct numerical simulations. An alternative physical possibility has been highlighted in a series of publications (see, in particular, Ref. [23]; see also Ref. [24] for a relevant topical review). More specifically, a system with a suitable non-Hermitian Hamiltonian can be rewritten as a Schr¨ odinger equation with a Hermitian Hamiltonian and a nonlinear source term. This nonlinear term becomes most significant at the so-called exceptional points, where resonance states overlap. In a different vein, which nevertheless bears considerable similarities to the above direction, it was recently proposed that a nonautonomous Hermitian system with a higher number of degrees of freedom can be reduced to a lower-dimensional PT -symmetric system, provided that suitable conditions are satisfied for the timedependent parameters of the higher-dimensional problem [25]. In what follows, however, we will consider cases that differ from the ones above in that our starting point will involve a non-Hermitian nonlinear system of a second order in time oscillator (rather than Schr¨ odinger, see below) type, in which 032108-1 1050-2947/2013/88(3)/032108(11) ©2013 American Physical Society CUEVAS, KEVREKIDIS, SAXENA, AND KHARE PHYSICAL REVIEW A 88, 032108 (2013) the nonlinearity will be intrinsic rather than the result of a reduction. Our aim herein is to provide a simple, prototypical nonlinear model whose linear analog is effectively the one used in the experimental investigations of Ref. [20]. Yet the nonlinear structure is such that it allows us to obtain detailed numerical and even considerable analytical insights into the phenomenology of such a nonlinear PT -symmetric oscillator dimer. In particular, after formulating and briefly analyzing the linear PT -symmetric coupled oscillator model, we incorporate into it a local cubic nonlinearity (which can, in general, be of soft or hard form, i.e., bearing a prefactor of potentially either sign). This type of potential, especially in its bistable form, is well known to be a canonical example of relevance to numerous physical settings, including (but not limited to) phase transitions, superconductivity, and field theories, as well as high-energy and particle physics; see, e.g., Ref. [26] and the review [27], as well as references therein. For the resulting PT -symmetric, coupled nonlinear oscillator model, we provide a detailed analysis of the existence and stability of breathing (i.e., time-periodic) states in the system. We focus particularly on symmetric and antisymmetric states that arise from the linear limit of the problem. We observe that symmetry-breaking-type bifurcations can arise for both the symmetric and antisymmetric branches and eventually highlight a nonlinear analog of PT phase transition whereby the two branches terminate hand in hand in a saddle-center bifurcation. To provide an analytical insight into the above results, we use the rotating-wave approximation (RWA), which approximates the system by the corresponding nonlinear Schr¨ odingertype PT -symmetric dimer for which everything can be solved analytically, including the stationary states, the symmetry breaking bifurcations, and even the full dynamics [9,11,28,29]. A direct comparison of the RWA-derived Schr¨ odinger dimer reveals natural similarities but also significant differences between the two models. For instance, in the way of similarities, both models bear symmetric and antisymmetric branches of solutions, both models bear principal symmetry-breaking [of the symmetric branch in the soft (attractive) nonlinearity and of the antisymmetric one for the hard (repulsive) nonlinearity] bifurcations, and both have PT phase transitions, involving the collision and disappearance of these two branches. On the other hand, in terms of substantial differences, it is naturally expected that for soft nonlinearities, the oscillator model can escape the potential well (and, hence, have collapse features) even when the Schr¨ odinger model cannot; see for a detailed recent discussion of such features in the Hamiltonian limit theworkofRef.[30]. More importantly for our purposes, another significant difference is that it turns out that both branches, namely both the symmetric and the antisymmetric one, have destabilizing symmetry-breaking bifurcations, even though in the Schr¨ odinger reduction, only one of the two branches (the symmetric for soft and the antisymmetric for hard, as mentioned above) sustains such bifurcations. It should also be mentioned that after completing the examination of the prototypical time-periodic states and their Floquet theory–based stability, we corroborate our bifurcation results by means of direct numerical simulations in order to explore the different dynamical evolution possibilities that arise in this system. These include, among others, the indefinite growth for the hard potential and the finite-time blowup for the soft potential. Our presentation of the two models and their similarities and differences proceeds as follows. In Sec. II, we briefly discuss the underlying linear model and the prototypical nonlinear extension thereof. In Sec. III, we discuss the numerical setup for analyzing the model and its solutions, while in Sec. IV, we provide a means of theoretical analysis in the form of the rotating-wave approximation. In Sec. V, we present the numerical results, separating the cases of the soft and hard potentials. Finally, in Sec. VI, we summarize our findings and present our conclusions, as well as some directions for future study. II. MODEL EQUATIONS AND LINEAR ANALYSIS We consider the system motivated by recent experimental realizations in electrical circuits of the form ¨ u=−ω2 0u+sv +γ˙ u, (1) ¨v=−ω2 0v+su −γ˙v. (2) Here ω0characterizes the internal oscillator at each mode; in the case of the electrical circuit model this is the oscillation of each of the charges within the dimer [20]. The term “proportional to s” reflects the coupling between the two elements in the dimer, while γis proportional to the amplification and resistance within the system. One can try to identify the eigenvalues of the system by using u=Aeit and v=Beit, but one obtains in this case a quadratic pencil for the relevant eigenvalue problem. It is thus easier to formulate this as a 4 ×4 first-order (linear) dynamical system according to ˙ u=p, (3) ˙ p=−u+γp +sv, (4) ˙v=q, (5) ˙ q=su −v−γq, (6) [where ω0has been rescaled without loss of generality to unity, and other quantities such as sand γand time have been rescaled by ω2 0(the first) and ω0(the latter two), respectively]. We then can seek solutions of the form u=Aeλt ,p=Beλt , v=Ceλt , and q=Deλt , to obtain a first-order eigenvalue problem which yields the following eigenvalues: λ=±−2+γ2±4s2−4γ2+γ4 √2.(7) These two pairs of imaginary (for small γ) eigenvalues will collide and give rise to a quartet for γ>γ PT, where γPT satisfies the condition γ4−4γ2+4s2=0.(8) Hence, Eq. (8) will define, at the linear level, the point of the so-called [1–3]PT phase transition and of bifurcation into the complex plane. Now our main interest in what follows will be to examine a prototypical nonlinear variant of the problem, which will be formulated as follows. In particular, we set up the form of the 032108-2 PT -SYMMETRIC DIMER OF COUPLED ... PHYSICAL REVIEW A 88, 032108 (2013) equations as follows: ¨ u=−u+sv +γ˙ u+u3,(9) ¨v=−v+su −γ˙v+v3.(10) Here, in parallel to what is done in the PT -symmetric Schr¨ odinger dimer typically [7,9,11], we have added a cubic onsite nonlinearity on each one of the nodes. For >0, this nonlinearity is soft, imposing a finite (maximal energy) height type of potential, enabling the possibility of indefinite growth by means of the escape scenario considered earlier, e.g., in Ref. [30] (see also references therein). On the other hand, for <0, the nonlinearity is hard, and the potential is monostable, bearing only the ground state at 0 and no possibility for such finite-time collapse (in the Hamiltonian analog of the model, only the potential for oscillations around the 0 state exists in this case of the hard potential). We now discuss the setup and numerical methods, as well as the type of diagnostics that we use for this system. In our description below, we follow an approach reminiscent of that in Ref. [31]. III. SETUP, DIAGNOSTICS, AND NUMERICAL METHODS A. Existence of periodic orbit solutions In order to calculate periodic orbits in the PT nonlinear oscillator dimer, we make use of a Fourier space implementation of the dynamical equations and continuations in frequency or gain (loss) parameter are performed via a pathfollowing (Newton-Raphson) method. Fourier space methods are based on the fact that the solutions are Tbperiodic; for a detailed explanation of these methods, the reader is referred to Refs. [32–34]. 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)= km  k=−km ykexp(ikωbt), (11) v(t)= km  k=−km zkexp(ikωbt), with kmbeing the maximum of the absolute value of the running index kin our Galerkin truncation of the full Fourier series solution. In the numerics, kmhas been chosen as 21. After the introduction of (11), the dynamical equations (9) and (10) yield a set of 2 ×(2km+1) nonlinear, coupled algebraic equations, Fk,1≡−ω2 bk2yk−iγωbkyk+Fk[V(u)] −szk=0,(12) Fk,2≡−ω2 bk2zk+iγωbkzk+Fk[V(v)] −syk=0,(13) with V(u)=u−u3.HereFkdenotes the discrete Fourier transform, Fk[V(u)] =1 N km  n=−km Vkm  p=−km ypexp i2πpn N ×exp −i2πkn N,(14) with N=2km+1. The procedure for Fk(v) is similar to the previous one. As u(t) and v(t) must be real functions, it implies that y−k=y∗ k,z −k=z∗ k. An important diagnostic quantity for probing the dependence of the solutions on parameters such as the gain or loss strength γ, or the oscillation frequency ωbis the averaged over a period energy, defined as H=1 TT 0 H(t)dt, (15) with the Hamiltonian (of the case without gain or loss) being H=˙ u2+˙v2+u2+v2 2− 4(u4+v4)−suv (16) and constituting a conserved quantity of the dynamics in the Hamiltonian limit of γ=0. B. Linear stability equations 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. (9) and (10) according to u=u0+ξ1,v= v0+ξ2. Then the equations satisfied to first order in ξnread ¨ ξ1+V(u0)ξ1−γ˙ ξ1−sξ2=0,(17) ¨ ξ2+V(v0)ξ2+γ˙ ξ2−sξ1=0,(18) 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 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)}.(19) The 4 ×4 monodromy eigenvalues =exp(iθ)are dubbed the Floquet multipliers (FMs) 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 [33,34]) and that the eigenvalues come in pairs {,∗}. As a consequence, due to the “simplicity” of the FM structure (one pair always at 1 and 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 1or−1) on the real line leading to one multiplier (in absolute value) larger than 1 and one smaller than 1. We will explore a scenario of this kind of instability in what follows. Having set up the existence and stability problem, we now complete our theoretical analysis by exploring the outcome of the RWA. 032108-3 CUEVAS, KEVREKIDIS, SAXENA, AND KHARE PHYSICAL REVIEW A 88, 032108 (2013) IV. AN ANALYTICAL APPROACH: THE ROTATING WAVE APPROXIMATION The RWA provides a means of connection with the extensively analyzed PT -symmetric Schr¨ odinger dimer [9,11,35]. This link follows a path similar to what has been earlier proposed, e.g., in Refs. [36–38]. In particular, the following ansatz is used to approximate the solution of the periodic orbit problem as a roughly monochromatic wave packet 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), (20) v(t)≈φ2(t)exp(iωbt)+φ∗ 2(t) exp(−iωbt). By supposing that ˙ φnωbφnand ¨ φnωb˙ φn(i.e., that φ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 (9) and (10) transform into a set of coupled Schr¨ odinger-type equations, 2iωb˙ φ1=ω2 b−1+3|φ1|2+iωbγφ1+sφ2, (21) 2iωb˙ φ2=ω2 b−1+3|φ2|2−iωbγφ2+sφ1, 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. (11) 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=κz1+|y1|2y1+iy1,(22) Ez1=κy1+|z1|2z1−iz1,(23) with E=1−ω2 b 3,κ=s 3,=γω b 3.(24) Recast in this form, Eqs. (22) and (23) are identical to Ref. [35, Eq. (6)]. We express y1and z1in polar form as follows: y1=Aexp(iθ1),z 1=Bexp(iθ2),ϕ=θ2−θ1,(25) and then rewrite the stationary equations as EA =κB cos(ϕ)+A3,(26) EB =κAcos(ϕ)+B3,(27) sin(ϕ)=−A/(κB)=−B/(κA).(28) In the Hamiltonian case γ=0, and, consequently, sin ϕ= 0. Three different solutions may exist therein, namely the symmetric, antisymmetric, and asymmetric solutions, given by A=B, (29) A2=E−κ=1−ω2 b−s 3(symmetric solution), A=−B, (30) A2=E+κ=1−ω2 b+s 3(antisymmetric solution), B=κ/A =s/(3A), A2=(E±E2−4κ2)/2(31) =1−ω2 b±1−ω2 b2−4s2/(6) (asymmetric solution). The symmetric solution derives from the linear mode located at ωS=√1−s, whereas the antisymmetric solution bifurcates from the mode at ωA=√1+s. Straightforwardly (by examining the quantity under the radical in its profile), the asymmetric solution exists for ωb⩽√1−2s, bifurcating via a symmetry-breaking pitchfork bifurcation from the symmetric solution if the potential is soft (>0). On the contrary, if the potential is hard, the asymmetric solution bifurcates from the antisymmetric solution and exists for ωb⩾√1+2s. The emerging (“daughter”) asymmetric solutions inherit the stability of their symmetric or antisymmetric “parent” and are therefore stable, whereas the respective parent branches become destabilized past the bifurcation point. For γ= 0, the asymmetric solution is no longer a stationary solution and only symmetric and antisymmetric solutions exist as exact stationary states in the PT -symmetric Schr¨ odinger dimer [as is directly evident, e.g., from Eq. (28)]. These solutions have A=Btaking the following values: A=1−ω2 b−s2−γ2ω2 b(3)1/2 (symmetric solution),(32) A=1−ω2 b+s2−γ2ω2 b(3)1/2 (antisymmetric solution),(33) and sin ϕ=−/κ.As|sin ϕ|⩽1, the solutions must fulfill γ⩽γSC, with γSC =s ωbi.e., there is a saddle-center bifurcation (namely the RWA-predicted nonlinear analog of the PT phase transition) taking place at γ=γSC. The average energy for both the symmetric and antisymmetric solutions is given by the same expression, H≈4ω2 bA2+3A4.(34) As both solutions coincide at the PT bifurcation critical point, their energy will also be the same therein. We now turn to the linear stability of the 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 (21) and then solving the ensuing [to O(δ)] eigenvalue problem, φ1(t)=y1+δ(a1e−iθt/T +b∗ 1eiθ∗t/T ), (35) φ2(t)=z1+δ(a2e−iθt/T +b∗ 2eiθ∗t/T ), with T=2π/ωbbeing the orbit’s period and θbeing the Floquet exponent (FE). The nonzero FEs are given by θ=±2π ω2 b2s2−γ2ω2 b−1−ω2 bs2−γ2ω2 b1/2 (symmetric solution),(36) 032108-4 PT -SYMMETRIC DIMER OF COUPLED ... PHYSICAL REVIEW A 88, 032108 (2013) θ=±2π ω2 b2s2−γ2ω2 b+1−ω2 bs2−γ2ω2 b1/2 (antisymmetric solution).(37) As instability is marked by an imaginary value of θ, the above expression implies that there is a stability change when the square-root argument becomes zero, i.e., γ=γstab, with γstab =4s2−1−ω2 b2 2ωb .(38) A straightforward analysis shows, in addition, that the symmetric (antisymmetric) solution experiences the change of stability bifurcation when ωbis smaller (larger) than 1. Thus, the symmetric (antisymmetric) solution is always stable when the potential is hard (soft) and stable if γ<γ stab when the potential is soft (hard). As we will see below, it is precisely this prediction of the RWA that will be “violated” from the full numerics of the PT -symmetric oscillator model. In particular, it will be found that in the latter model, both branches in each (soft or hard) case can become unstable through such symmetry-breaking bifurcations within suitable parametric regimes. As a final theoretical remark, it is relevant to point out that while no additional stationary solutions have been argued to exist in the Schr¨ odinger dimer, a “reconciliation” with the expected picture of a pitchfork bifurcation has been offered, e.g., in Ref. [35] (see also references therein) through the notion of the so-called ghost states. These are solutions for which the parameter Ebecomes complex (and the pitchfork bifurcation resurfaces in diagnostics such as the imaginary part of E). Nevertheless, and especially because complex eigenvalue parameters are of lesser apparent physical relevance in models such as the oscillator one considered herein, we will not further pursue an analogy to such ghost states here but will instead restrain our considerations hereafter to the symmetric and antisymmetric branches of time-periodic solutions. V. NUMERICAL RESULTS AND COMPARISON WITH THE ROTATING-WAVE APPROXIMATION In this section, we identify the relevant, previously discussed, periodic orbits by numerically solving in the Fourier space the dynamical equations set (9) and (10).Wehave considered the two cases of =±1, with =1 corresponding to the soft potential, while =−1 corresponds to the hard potential case. We analyze the properties of phase-symmetric and phaseanti-symmetric solutions, characterized respectively by the following properties: u(0) =v(0),˙ u(0) =−˙v(0) (symmetric); (39) u(0) =−v(0),˙ u(0) =˙v(0) (antisymmetric), or, in terms of the Fourier coefficients, yk=z∗ k(symmetric),y k=−z∗ k(antisymmetric).(40) It is worth remarking that, for γ=0, the solutions are time reversible and, consequently, we select ˙ u(0) =˙v(0) =0. We recall that in the Hamiltonian limit, the RWA predicts that the symmetric (antisymmetric) solution becomes unstable at ωb=ωS≡√1−2s(ωb=ωA≡√1+2s) in the soft (hard) case, giving rise to an asymmetric solution. For the particular case of s=√63/32, the results for the soft case are shown (for γ=0) in the left panel of Fig. 1, while those for the hard potential case are shown in the right panel of the figure (for s= √15/8). The solid lines of the direct numerical computation appear to have good agreement with the dashed lines of the RWA, as regards the predicted amplitudes of the asymmetric node equilibrium values, at least near the bifurcation point. Interestingly, this agreement is considerably better in the hard case than in the soft nonlinearity case. This will be a continuous theme within the results that follow, i.e., we will see that the hard case is generally very accurately described by the RWA (even in the presence of gain or loss), while the soft case is only well approximated sufficiently close to the linear limit. It should be noted that a fundamental difference between the nonlinear Schr¨ odinger-type dimer and the φ4oscillator one is expected as the amplitude of the solution increases (and 0.5 0.55 0.6 0.65 0.7 0.4 0.5 0.6 0.7 0.8 0.9 1 u(0), v(0) ωb 1.4 1.45 1.5 1.55 1.6 1.65 1.7 0.4 0.6 0.8 1 1.2 1.4 1.6 u(0), −v(0) ωb FIG. 1. (Color online) Amplitude of the couplers for the asymmetric solution in the Hamiltonian case of γ=0andfors=√63/32, =1 (left panel) and s=√15/8, =−1 (right panel). The dashed lines correspond to the predictions of the RWA theory (see the discussion in the text). It is worth noticing that in the left panel, asymmetric solutions are unstable towards finite-time blowup for frequencies smaller than those shown in the figure. On the contrary, in the right panel, solutions are always stable even for frequencies higher than those shown in the picture. 032108-5 CUEVAS, KEVREKIDIS, SAXENA, AND KHARE PHYSICAL REVIEW A 88, 032108 (2013) 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 1.2 γ ωb I II III IV Curve 1 Curve 2 Curve 3 Curve 4 V FIG. 2. (Color online) Plane with curves and regions of solutions that share the same properties when =1ands=√15/8(see text). In the nonlabeled region of the top right, neither symmetric nor antisymmetric solutions exist. Here Curve 1 corresponds to the linear limit, Curve 2 denotes the PT phase transition curve, Curve 3 indicates the destabilization of the A branch, while Curve 4 is associated with the destabilization of the symmetric branch. The detailed description of the different regions enclosed by the curves is offered in the text. Dashed lines correspond to the RWA predictions. Notice that the colors of the dashed lines are inverted with respect to those of the numerical results for a better visualization, i.e., the blue dashed line is the theoretical prediction corresponding to the red solid one, while the red dashed line represents the prediction corresponding to the blue solid one. This inversion pattern is followed also in all figures comparing theory and numerical computations from here onward. Dots mark the parameters for which simulations are performed in Fig. 5. so does the deviation from the symmetry-breaking point). In particular, the former model due to its norm conservation does not feature finite-time collapse (or any type of infinite growth for that matter when γ=0). On the contrary, the latter model has the potential for finite-time collapse when the amplitude of the nodes exceeds the unit height of the potential (see a detailed analysis of this “escape” phenomenology in the recent work of Ref. [30] and references therein). It is thus rather natural that the two models should significantly deviate from each other as this parameter range is approached. A. Soft potential We analyze, first, the soft case (=+1) in the presence now of the gain or loss term proportional to γ. Figure 2shows aγ-ωbfull two-parameter plane summarizing the existence properties of the solutions and separating the different regimes thereof. (i) Curve 1 corresponds to the linear modes, obtained in Sec. II. The increasing part of the curve fulfills ω1= 1−γ2/2−s2−γ2+γ4/4 and corresponds to symmetric linear modes at γ=0; the decreasing part ω2= 1−γ2/2+s2−γ2+γ4/4 holds for the branch stemming from the antisymmetric linear modes of γ=0[cf. with Eq. (7)]. These two classes of linear modes collide and disappear hand in hand at the value of γpredicted by Eq. (8). For this soft case, solutions are expected to exist (in analogy with the Schr¨ odinger case) for ωb<ω 1for the symmetric branch and for ωb<ω 2for the antisymmetric branch (for a given value of γ). (ii) Curve 2 corresponds to the PT phase transition, i.e., at the nonlinear level it corresponds to the saddle-center bifurcation leading to the termination of both the antisymmetric and symmetric branches. Above this curve, there do not exist any periodic orbits and the system dynamics generically leads to indefinite growth. This curve overlaps with Curve 1 for high ωb. (iii) Curves 3 and 4 separate stable and unstable solutions of the antisymmetric and symmetric branches, respectively. In particular, they represent the threshold for the symmetrybreaking bifurcation of the corresponding branches. The regions limited by the above four curves are the following ones: (a) Region I: Both symmetric and antisymmetric solutions are unstable, as they have both crossed the instability inducing curves 3 and 4. (b) Region II: Symmetric solutions are stable (as they have not crossed curve 4), whereas antisymmetric solutions are unstable (since after bifurcating from the decreasing part of curve 1, they have crossed the instability threshold of curve 3). (c) Region III: Symmetric solutions do not exist (as such solutions only exist to the left of the increasing part of curve 1) and antisymmetric solutions are unstable (as they have crossed the instability threshold of curve 3). (d) Region IV: Symmetric solutions do not exist (for the same reason as in III) and antisymmetric solutions are stable, i.e., they are stable between their bifurcation point (the decreasing part of curve 1) and instability threshold (curve 3). (e) Region V: Stable antisymmetric and symmetric solutions coexist in this narrow region prior to their termination in the saddle-center bifurcation occurring on curve 2. In Fig. 3, some typical examples of monoparametric continuations of the relevant solutions are given. The top panels illustrate continuations as a function of the gain or loss parameter γfor a fixed value of the frequency ωb=0.45, while the bottom ones illustrate a continuation as a function of ωbfor agivenvalueofγ=0.4. The comparison of the numerically obtained symmetric and antisymmetric solutions with the ones obtained analytically by virtue of the RWA (in reverse colors, see the figure) is also offered. It can be inferred that, generally, the RWA offers a reasonable qualitative match to the numerically exact, up to a prescribed tolerance, solutions, although clearly quantitative comparison is less good, at least for the low (i.e., far from the linear limit) value of ωbin the top panels. This deficiency of the method (explained also previously) as one departs far from the linear limit is more clearly illustrated in the ωbdependence. Close to the limit, the RWA does an excellent job of capturing both branches, but things become progressively worse as ωbdecreases. Furthermore, as discussed above, the right panels showcase a stability change as occurring for both branches, while the RWA predicts a destabilization solely of the symmetric branch for this soft nonlinearity case. In the above two-parametric diagram, we have only varied the frequency of the breathers ωband the strength of the gain or loss γ. To illustrate how the results vary as the final (coupling) parameter of the system varies, we have shown the same 032108-6 PT -SYMMETRIC DIMER OF COUPLED ... PHYSICAL REVIEW A 88, 032108 (2013) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 <H> γ 00.2 0.4 0.6 0.8 1 100 102 104 106 108 1010 |Λ| γ 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 <H> ωb 0.4 0.6 0.8 1 100 105 1010 1015 |Λ| ωb FIG. 3. (Color online) Averaged energy (top left) and modulus of the Floquet multipliers (top right) [only the multipliers with moduli higher than one are shown] as a function of the gain or loss parameter γfor =1, ωb=0.45, and s=√15/8. Notice the logarithmic scale in the yaxis of the latter graph. In the left panel, the blue [lower] (red [upper]) solid line corresponds to the symmetric (antisymmetric) solution, while the red (blue), i.e., reversed colors, dashed lines correspond to the symmetric (antisymmetric) branch RWA predictions. One can observe a reasonable agreement between the numerical energy and its theoretical prediction counterpart, although some discrepancy is clearly visible for this low (i.e., far from the linear limit) value of the frequency. The bottom panel shows similar comparisons but now for the dependence on ωbfor fixed γ=0.4. Here it is evident that the agreement is good close to the linear limit of larger values of ωband progressively worsens as one departs from this limit (by lowering the breather frequency). features as in Figs. 2and 3for roughly half the coupling strength in Fig. 4. We observe that the region of stability of the different solutions (and especially of the symmetric one) has nontrivially changed on the considered parametric variation. Nevertheless, sufficiently close to the linear limit of emergence of the two solutions, the RWA remains a reasonable description of their existence and stability, as well as of the saddle-center bifurcation leading to their disappearance. On the other hand, as one further deviates from this limit towards lower frequencies, the RWA fails to capture the observed phenomenology by deviating from the critical point for the saddle-center bifurcation, missing the complex boundary of stability of the symmetric branch and missing altogether the destabilization of the antisymmetric branch. Some examples of the evolution of unstable solutions for s=√15/8areshownatFig.5. In these cases, the perturbation was induced solely from numerical truncation errors. In most cases, the instabilities lead to finite-time blowup, whereas in some cases a switching between the oscillators is observed, i.e., a modulated variant of the time-periodic solution arises as a result of the instability. However, generally speaking, if the perturbation is above a threshold and/or the value of the gain or loss parameter γis sufficiently high and/or the solution frequency sufficiently low, the instability manifestation will typically result in a finite-time blowup. This is strongly related to the escape dynamics considered in Ref. [30]. On the contrary, we want to highlight that this is different than the “worst case scenario” of the Schr¨ odinger dimer of the RWA. As illustrated in Ref. [28], in the latter, at worst an exponential (indefinite) growth of the amplitude may arise (but no finitetime blowup). Our numerical computations indicate that the switching behavior reported above is only possible provided that the growth rate (i.e., the FE) of the periodic solution is small enough. For symmetric (antisymmetric) solutions, this can be achieved close to curve 4 (3). This condition is, however, not sufficient, as shown in the top right panel of Fig. 5.Wehave also analyzed the outcome for solutions with γ>γ SC, taking as initial condition a solution for γ<γ SC; in that case, we have observed that, although the generic scenario corresponds to a blowup, isolated cases of modulated dynamics may arise when such a profile is used as initial condition for the simulation. B. Hard potential We now briefly complement these results with ones of the far more accurately approximated (by the RWA) hard potential case. This disparity in the much higher level of adequacy of the theoretical approximation here is clearly induced by 032108-7 CUEVAS, KEVREKIDIS, SAXENA, AND KHARE PHYSICAL REVIEW A 88, 032108 (2013) 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 γ ωb Curve 1 Curve 2 Curve 3 Curve 4 0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 <H> γ 0 0.05 0.1 0.15 0.2 0.25 0.3 100 101 102 |Λ| γ FIG. 4. (Color online) Same as Fig. 2(top panel) and as the top panel of Fig. 3(bottom panels) but for the soft nonlinearity case of lower coupling s=√63/32. Again, the full numerical results are offered by the solid lines for the two branches [blue for symmetric (lower) and red for antisymmetric (upper)], while the dashed lines with reverse colors correspond to the RWA. 0 500 1000 1500 2000 2500 3000 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 (a) u(t), v(t) t0100 200 300 400 −5 −4 −3 −2 −1 0 1 (b) u(t), v(t) t 0 5 10 15 20 25 −10 −8 −6 −4 −2 0(c) u(t), v(t) t0500 1000 1500 2000 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 (d) u(t), v(t) t FIG. 5. (Color online) Evolution of unstable solutions for the soft potential; the instability is driven only by the numerical truncation errors. (Top left panel) Symmetric solution with ωb=0.45 and γ=0.01. (Top right panel) Symmetric solution with ωb=0.45 and γ=0.45. (Bottom left panel) The antisymmetric solution with ωb=0.70 and γ=0.66. (Bottom right panel) The antisymmetric solution at γ=0.16 and ωb=1.2 is taken as initial condition for γ=0.51; notice that for this value of ωb, the saddle-center bifurcation takes place at γSC =0.168. 032108-8 PT -SYMMETRIC DIMER OF COUPLED ... PHYSICAL REVIEW A 88, 032108 (2013) 0.8 0.9 11.1 1.2 1.3 1.4 0 0.1 0.2 0.3 0.4 0.5 γ ωb I II III Curve 1 Curve 2 Curve 3 Curve 4 FIG. 6. (Color online) Same as in Fig. 2but in the hard potential case, i.e., when =−1. The only change with respect to that figure consists of an interchange of the color and thickness between Curves 3 and 4. Dots mark the parameters for which simulations are performed in Fig. 8. the absence of finite-time collapse in the latter model in consonance (in this case) with its RWA analog. In this case (see Fig. 6), the branches exist to the right of the linear curve (for a given γ), as is expected in the case of a hard and defocusing nonlinearity. Curve 1 illustrates the linear limit; once again its growing part (ω1) is associated with the symmetric solutions, while its decreasing part (ω2) is associated with the antisymmetric solutions, with their collision point representing the linear PT transition point. In fact, from that point, emanates curve 2 which is the nonlinear PT phase transition curve, i.e., the locus of points where the symmetric and antisymmetric solutions collide and disappear for the nonlinear problem. Notice the very good comparison of this curve with the theoretical prediction of the RWA for γSC. Curve 3, on the other hand, denotes the point of destabilization of the antisymmetric branch, which is, in fact, expected also from the RWA, whose prediction is once again (red dashed line) in remarkable agreement with the full numerical result. The final curve in the graph, i.e., the green solid line of curve 4 denotes a narrow parametric window beyond which (or, more appropriately, between which and curve 2) even the symmetric branch of time-periodic solutions is destabilized. This, as indicated previously, is a feature that is not captured by the RWA but is particular to the oscillator system (as is correspondingly the destabilization of the antisymmetric branch in the soft potential case). We now discuss the existence and stability of the branches in each of the regions between the different curves. (a) Region I: Only symmetric solutions exist (i.e., to the right of the increasing part of curve 1, but to the left of its decreasing part). They bifurcate from the linear modes when ωb⩾ω1. The symmetric modes are stable in this regime. (b) Region II: Symmetric and antisymmetric solutions exist and are both stable. The antisymmetric solutions have now bifurcated for ωb>ω 2. (c) Region III: Symmetric solutions are stable, whereas antisymmetric solutions are unstable. Here, the symmetry breaking bifurcation destabilizing the antisymmetric solutions 00.1 0.2 0.3 0 0.5 1 1.5 2 <H> γ 0 0.1 0.2 0.3 1 1.5 2 2.5 3 3.5 |Λ| γ 0.8 0.9 11.1 1.2 1.3 1.4 0 0.5 1 1.5 2 2.5 3 3.5 4 <H> ωb 0.8 0.9 1 1.1 1.2 1.3 1.4 1 1.5 2 2.5 3 3.5 4 |Λ| ωb FIG. 7. (Color online) Similar to Fig. 3but for ωb=1.25 and =−1inthetoppanelsandforγ=0.3and=−1 in the bottom panels. Both the existence (left) and stability (right) diagrams of symmetric [blue (upper in left panels) solid] and antisymmetric [red (lower in left panels) solid] branches and their rotating-wave approximations (the latter with reverse colors and dashed lines) are shown. 032108-9