Quasi-energies, parametric resonances, and stability limits in ac-driven \mathcal{PT}-symmetric systems

# Quasi-energies, parametric resonances, and stability limits in ac-driven $\mathcal{PT}$-symmetric systems

## Abstract

We introduce a simple model for implementing the concepts of quasi-energy and parametric resonances (PRs) in systems with the symmetry, i.e., a pair of coupled and mutually balanced gain and loss elements. The parametric (ac) forcing is applied through periodic modulation of the coefficient accounting for the coupling of the two degrees of freedom. The system may be realized in optics as a dual-core waveguide with the gain and loss applied to different cores, and the thickness of the gap between them subject to a periodic modulation. The onset and development of the parametric instability for a small forcing amplitude () is studied in an analytical form. The full dynamical chart of the system is generated by systematic simulations. At sufficiently large values of the forcing frequency, , tongues of the parametric instability originate, with the increase of , as predicted by the analysis. However, the tongues following further increase of feature a pattern drastically different from that in usual (non-) parametrically driven systems: instead of bending down to larger values of the dc coupling constant, , they maintain a direction parallel to the axis. The system of the parallel tongues gets dense with the decrease of , merging into a complex small-scale structure of alternating regions of stability and instability. The cases of and are studied analytically by means of the adiabatic and averaging approximation, respectively. The cubic nonlinearity, if added to the system, alters the picture, destabilizing many originally robust dynamical regimes, and stabilizing some which were unstable.

###### pacs:
11.30.Er; 72.10.Fk; 42.79.Gn; 11.80.Gw

Recently, a great deal of interest has been drawn to systems featuring the  (parity-time) symmetry. Originally, they were introduced as quantum systems with spatially separated and symmetrically placed linear gain and loss. A fundamental property of non-Hermitian -symmetric Hamiltonians is the fact that their spectra remain purely real, like in the case of the usual Hermitian Hamiltonians, provided that the strength of the non-Hermitian part of the Hamiltonian, , does not exceed a certain critical value, , at which the  symmetry breaks down. Such linear systems were recently implemented experimentally in optics, due to the fact that the paraxial propagation equation for electromagnetic waves can emulate the quantum-mechanical Schrödinger equation. In particular, the -symmetric quantum system may be emulated by waveguides with spatially separated mutually balanced gain and loss elements. Additional interest to the optical realizations of the  symmetry has been drawn by the possibility to implement this setting in nonlinear waveguides. Unlike the ordinary models of nonlinear dissipative systems, where stable modes exist as isolated attractors, in -symmetric systems they emerge in continuous families, similar to what is commonly known about conservative nonlinear systems. However, the increase of the gain-loss coefficient, , leads to the shrinkage of existence and stability regions for the -symmetric modes, which completely vanish at some . Especially convenient for the study of such solutions are models of dual-core couplers, with the balanced gain and loss applied to different cores, and the cubic nonlinearity acting in both. The objective of the present work is to extend the well-known concepts of quasi-energy and parametric resonances (PRs) to -symmetric systems. In terms of the coupler model, the parametric drive (with ac amplitude and dc amplitude ) can be implemented by making the thickness of the gap between the two cores in the coupler a periodic function of the evolution variable (propagation distance, in terms of optics), with frequency . Using a combination of analytical and numerical methods, we produce the full dynamical chart for both linear and nonlinear versions of the system in the plane of  for different values of . The charts feature tongues of parametric instability, whose shape has some significant differences from that predicted by the classical theory of the PR in conservative systems. In addition, the cases of small and large  are studied by means of the adiabatic and averaging approximations, respectively.

## I Introduction and the model

A commonly known principle of quantum mechanics is that operators representing physical variables, including the Hamiltonian, must be Hermitian, to guarantee the reality of the respective spectra. Recently, it was found that the class of physically acceptable Hamiltonians may be extended by admitting operators with anti-Hermitian [dissipative and antidissipative (gain)] parts, if they commute with the (parity-time) transformation (1). The spectrum of such a -invariant Hamiltonian may be purely real, provided that the strength of the anti-Hermitian part, , does not exceeds a critical value, , at which complex eigenvalues emerge. In the simplest form, a -invariant Hamiltonian contains a complex potential in which the real part is an even function of coordinates, while the imaginary part is odd.

While in the quantum theory this concept remains a rather abstract one, the similarity of the quantum-mechanical Schrödinger equation to the paraxial propagation equation in optics has suggested a possibility to emulate the -symmetry in optical waveguides with symmetrically placed mutually balanced gain and loss elements. This implementation was proposed theoretically (2) and then demonstrated in experiments (3). Recently, another experimental realization of the  symmetry was reported in a system of coupled electronic oscillators (4). There is a possibility to realize such a setting in Bose-Einstein condensates (5) too.

These possibilities have drawn a great deal of interest to the analysis of various -symmetric dynamical models (6); (7). In particular, the natural occurrence of the self-focusing in optical media stimulated the study of nonlinear effects in such systems. The interplay of the paraxial diffraction, Kerr nonlinearity, and a spatially periodic complex potential, which represents the -symmetric part of the system, gives rise to bright solitons (8), whose stability was rigorously analyzed in Ref. (9). Dark solitons have been predicted too in this context, assuming the self-defocusing Kerr nonlinearity (10). In addition, bright solitons have been predicted in -symmetric systems with the quadratic (second-harmonic-generating) nonlinearity (11). Systems of another type, which readily maintain -symmetric solitons, and actually make it possible to obtain such solutions and their stability conditions in an exact analytical form, are provided by models of linearly-coupled dual-core waveguides, with the balanced gain and loss applied to the two cores, and the intrinsic Kerr nonlinearity present in each one (12)-(14).

Further, nonlinear modes were studied in many discrete -symmetric systems with the cubic nonlinearity. These include linear (15) and circular (16) chains of coupled elements and aggregates of coupled -symmetric oligomers (dimers (17), quadrimers, etc.) (18). General features of the dynamics in -symmetric lattices of single-component nonlinear -symmetric lattices have been recently considered in Ref. (19). Related to these models are systems featuring discrete -balanced elements embedded into continuous nonlinear conservative media (20); (21), which, in particular, allow one to find exact solutions for solitons pinned to the -symmetric dipole (21).

The inclusion of the Kerr nonlinearity into the conservative part of the system suggests that its gain-and-loss part can be made nonlinear too, by adopting mutually balanced cubic gain and loss terms  (22); (23); (24). Effects of combined linear and nonlinear  terms on the existence and stability of solitons were studied in (25).

In the usual nonlinear dissipative systems, stable dynamical regimes exist as isolated solutions (attractors), such as limit cycles or strange attractors (26) in discrete settings, or dissipative solitons in dissipative media (27). On the contrary to that, in -symmetric systems stable modes form continuous families, similar to the generic situation in conservative systems. However, the increase of the gain-loss coefficient () in the -symmetric system leads to shrinkage of existence and stability domains for the families, which eventually vanish at .

In many cases, quantum systems are driven by time-periodic (ac) external forces. In that case, the eigenstates are characterized by particular values of the quasi-energy, defined as per the Floquet theorem, see, e.g., Ref. (28). A related concept is the parametric resonance (PR) in linear and nonlinear systems subject to the action of the parametric drive (29). In the context of -symmetric systems, it is interesting to analyze the realization of the quasi-energy and PRs in ac-driven systems, both linear and nonlinear ones, including mutually balanced gain and loss terms. This analysis is the subject of the present work. We perform it in the framework of a basic model with two degrees of freedom (a dimer), represented by complex amplitudes and , and the coupling containing an ac-modulated term, which represents the parametric drive:

 iddtψA = iψA+[V0+2V1cos(ωt)]ψB+χ|ψA|2ψA, iddtψB = −iψB+[V0+2V1cos(ωt)]ψA+χ|ψB|2ψB. (1)

Here the equal gain and loss coefficients in the equations for and are scaled to be , and are strengths of the constant (dc) and ac-modulated couplings, is the modulation frequency, and is the nonlinearity coefficient, which may be taken as or in the linear and nonlinear versions of the system. Note that by the complex conjugation and replacing Eqs. (1) with are transformed into each other, therefore we do not consider . Also, by an obvious transformation , it is possible to reverse the sign of , and shift can independently reverse the sign of , therefore we can fix both and to be positive.

According to the concept of the quasi-energy, quasi-stationary solutions to Eqs. (1) for are sought for as

 ψA(t)=exp(−iEt)ΨA(t), ψB(t)=exp(−iEt)ΨB(t) (2)

where are periodic functions with period , and the quasi-energy belongs to interval . A corollary of Eqs. (1) is the balance equation for the total power in the presence of the mutually symmetric gain and loss (19)

 ddt(|ψA|2+|ψB|2)=2(|ψA|2−|ψB|2). (3)

This system can be implemented as a model of a dual-core optical waveguide operating in the CW (continuous-wave) regime, so that Eqs. (1) do not include terms accounting for the transverse diffraction or group-velocity dispersion, in the spatial-domain and temporal-domain realizations, respectively. Variable is actually the propagation distance, and the ac modulation represents periodic variation of the thickness of the gap separating the cores (30). As suggested in Refs. (12), terms in Eqs. (1) represent the gain and loss applied to the two cores. A similar system, without the ac-modulating terms, but with the addition of cubic gain and loss terms, was proposed as a model of a -symmetric dimer in Ref. (23).

The main issue considered in this work is the effect of the symmetry on PRs in the framework of the linear and nonlinear versions of system (1). First, in Sections II.A-B we develop an analytical approach for the case of moderate , predicting at what values of the ac drive with fixed gives rise to the PR at infinitesimally small , and how tongues of the parametric instability expand with the increase of . Comprehensive dynamical charts of the linearized version of system (1) and the full  nonlinear system, built on the basis of systematic numerical simulations, are presented in Sections II.C-D, respectively. The structure of the charts at large is very different from the picture of PRs in usual (non-) driven dynamical systems. In Section III we present, in a brief form, solutions of the system for slow (small ) and fast (large ) time-periodic modulations, which can be treated, respectively, by means of the adiabatic and averaging approximations. The paper is concluded by Section IV.

## Ii The analysis for moderate driving frequencies ω

### ii.1 Perturbative analysis of the linear system

With a small amplitude of the ac drive, , a perturbative solution to the linearized version of Eqs. (1), with quasi-frequency (or quasi-energy, in terms of quantum -symmetric systems), is looked for in the form of Eq. (2), adopting the lowest-order truncation for the Fourier expansion of the -periodic amplitudes (29),

 ΨA,B(t)=Ψ(0)A,B+Ψ(+)A,Be−iωt+Ψ(−)A,Be+iωt, (4)

with constants . Substituting this ansatz into Eqs. (1) and keeping, in line with the form of truncation (4), the zeroth and first harmonics in the ensuing expansion, gives rise to a system of linear homogeneous equations for the vectorial set of six amplitudes,

 Ψ={Ψ(0)A, Ψ(+)A, Ψ(−)A, Ψ(0)B, Ψ(+)B, Ψ(−)B}. (5)

Nontrivial solutions exist under the condition that the determinant of the system vanishes, which leads to the following equation for :

 ∣∣ ∣ ∣ ∣ ∣ ∣ ∣∣i−E00V0V1V10i−E−ω0V1V0000i−E+ωV10V0V0V1V1−i−E00V1V000−i−E−ω0V10V000−i−E+ω∣∣ ∣ ∣ ∣ ∣ ∣ ∣∣=0 . (6)

Each eigenvalue is associated with eigenvector of matrix (6). Complex solutions of Eq. (6) for imply instability of the linear system due to the spontaneous breaking of the symmetry, or excitation of a PR, or an interplay of both mechanisms (as shown below, all these cases are possible). Real eigenvalues correspond to stable oscillatory solutions.

In the zeroth-order approximation, , Eq. (6) yields an obvious result: the spectrum remains real under the condition

 V20>(V20)min≡1 (7)

[more explicitly, , if scaling for the gain-loss coefficient is not fixed]. In this approximation, the six eigenvalues are

 E=±√V20−1,E=±ω±√V20−1, (8)

where different signs are mutually independent.

Stability and instability areas in the plane of the ac and dc coupling strengths, , at four different values of the driving frequency,

 ω=0.4,0.8,1.5,2√3≈3.46, (9)

as predicted by a numerical solution of algebraic equation (6), is displayed in Fig. 1. The plots also show stability areas (covered by dots) as obtained in Section III by means of systematic direct simulations of the linear version of Eqs. (1), i.e., the one with . Circles on the vertical axis in Fig. 1 corresponds to parametric resonances, as described in the next subsection.

### ii.2 Parametric resonances

In the limit of , when the eigenfrequency of free oscillations in the linearized system is , see Eq. (8), the standard theory (29) produces the PR (i.e., the onset of the instability of the linear system, different from the instability due to breaking of the symmetry) at the following values of the driving frequency:

 ω(n)res=2n√V20−1, (10)

where is the order of the resonance. Thus, for given driving frequency , Eq. (10) predicts the resonance at

 (V0)(n)res=√1+(nω/2)2. (11)

In Fig. 1, circles on the vertical axis in each plot mark positions of the PR as predicted by Eq. (11).

Note that in the case of , Eq. (11) yields values

 (V0)(1)res=√1.04≈1.02, (V0)(2)res=√1.16≈1.08, (12)

which are close to the above-mentioned critical value, , induced by the balanced combination of the gain and loss terms, see Eq. (7). The wide and narrow instability stripes, observed in the panel of Fig. 1 corresponding to , originate from two critical values (12) at small values of . This situation actually implies a merger of the -breaking and PR-induced instabilities for relatively small values of .

At larger , the and PR instabilities are clearly separated as . For instance, at , Eq. (11) yields for the fundamental parametric resonance (). In accordance with this, the respective panel of Fig. 1 features a separate instability tongue, originating, at small , from point , and expanding with the increase of . The additional tongue observed in the same panel, originating from at , is predicted by Eq. (11) with . These tongues are quite similar to those known in the usual (non-) model of the parametric instability.

### ii.3 Direct simulations of the linear system on the parameter grid

To produce an accurate stability chart for the linear version of the underlying system, i.e., Eqs. (1) with , we ran systematic simulations of the equations on a square grid in the plane of with limits , and steps . We show results at the four fixed values of the driving frequency from set (9). Initial values for each run of the simulations were calculated as per Eqs. (2) and (4), with set (5) calculated as an eigenvector of matrix (6). The simulations were performed by means of a standard fourth-order Runge-Kutta code. To check the accuracy of the calculations, we verified that the numerical solutions would satisfy balance equation (3).

To distinguish between stable and unstable solutions in the present system, we checked whether the maximum of or would double after a sufficiently long period of time, in comparison to the modulation period, . It has been found that the lack of the doubling is an appropriate indicator of the stability of a given solution (in other words, the solutions which could be identified as stable ones would never double their largest absolute value). Each parameter set categorized as stable is marked by a dot in Fig. 1. Unstable solutions were identified as those breaking the above-mentioned no-doubling condition, which was followed by visual checks, which corroborate that the so identified unstable solutions indeed suffer blowup. In Fig. 1, unstable solutions are represented by missing dots on the selected grid.

The stability of -periodic solutions can also be numerically analyzed by computing the eigenvalues of the monodromy matrix for small perturbations around the given solution (31) (in the case of the linear problem, this effectively amounts to the stability of the zero solution). Writing (1) with as a system of four real equations with Jacobian implies that , the Jacobian of the Poincaré map, satisfies the variational equation

 ddt(JP)=JF(JP) (13)

with initial condition . Solving (13) and computing the eigenvalues of gives the Floquet multipliers which indicate instability if any of the multipliers have modulus greater than one. In Fig. 2 we show that this Floquet multiplier calculation predicts the exact same (in)stability intervals as that which we observed using the doubling condition of our direct integration as explained above.

Figure 1 clearly demonstrates both similarities and nontrivial differences between the present stability charts and those well known in the PR theory. Indeed, different numerical stability tongues originating at do not bend down, as in the usual theory, but keep the direction parallel to the axis of . This difference suggests that the truncation adopted in Eq. (4) becomes irrelevant at large , as a larger number of harmonics get involved into the dynamics, as corroborated by spectra of particular solutions displayed below. With the decrease of , the number of tongues becomes large, while their widths and gaps between adjacent ones get smaller. At small , the stability chart features a complex small-scale structure, which has no counterpart in the ordinary PR theory. It is important to point out that, while the quasi-energy-based analysis presented in Section II.A is expected (and indeed found in Fig. 1) to be valid for small , the PR theory of section II.B is relevant in a different limit. In particular, it accurately captures the dynamics at large , where the PR points are no longer clustered, and their tongues shrink accordingly. It is clearly observed in Fig. 1 that the horizontal tongues for the cases of, e.g., and naturally emanate from the PR points of . An additional relevant point to mention in that light is that higher-order resonances appear to be initiated “deeper” along the axis, although this observation may be affected by a finite extension of our simulations, and of a weak growth rate of these instabilities at small .

Figure 3 shows characteristic examples of the solutions for select combinations of parameter values, both stable (quasi-periodic solutions) and unstable (growing ones). Due to the presence of the gain in the system, unstable solutions grow indefinitely, with (at most) an exponential rate (19). As a result, at some point in time, after the temporal interval displayed in respective panels of Fig. 3, each amplitude increases rapidly (eventually the evolution accrues numerical errors).

For a few stable oscillatory solutions, Fig. 4 shows the power spectra generated by means of the DFFT (discrete fast Fourier transform),

 power A,B(fk)=∣∣ ∣∣1NN∑j=1ψA,B(tj)exp(2πi(j−1)kN)∣∣ ∣∣2, (14)

where , is the total integration time, and for even . All the spectra show a prominent peak at the driving frequency, . One also observes a set of weaker higher-frequency peaks corresponding to the frequency multiplication by the parametric drive.

To resolve the higher frequencies (smaller periods), we have performed the FFT on smaller windows. For example, within the first few seconds of the simulations for small , we would observe a dominant frequency in the spectrum consistent with the prediction, , given by Eq. (8) with . As and evolve in time, the high frequencies lead to modulations. For parameter values which are buried deepest in stability regions in Fig. 1, the solutions feature the least amount of such fluctuations at high frequencies. These fluctuations account for the small-amplitude pattern in the high-frequency part of the spectra in Fig. 4, which appears when the FFT is performed on a large interval, such as . All of the frequencies which appear in Fig. 4 are linear combinations of those predicted at and the driving frequency, .

### ii.4 The nonlinear propagation

As indicated above, the initial values of , used in the simulations of the linear version of Eqs. (1), were generated as per Eqs. (2) and (4), into which eigenvectors (5) of the matrix from Eq. (6) were substituted. We now use the same inputs for the systematic simulations of the full nonlinear form of Eqs. (1) with .

Figure 5 shows the so-generated nonlinear counterpart of Fig. 1, on the same grid of parameter values. In the upper two predicted-as-stable (gray) regions of Fig. 5, the tests generally show that the solutions are more likely to be unstable, in comparison with the linear system. In the inner portions of the lower gray regions, many solutions stay stable in the presence of the nonlinearity, while they are more likely to be unstable at the boundary of the region, in comparison with the linear case.

Some stable solutions of the linear system remain apparently robust when the nonlinearity is added, but with different amplitudes, in comparison to their linear counterparts. Profiles shown in Figs. 6 and 7 are nonlinear extensions of those displayed above in Fig. 3. As in the case of the linear system, the balance equation (3) remains valid for stable solutions, while the unstable solutions grow indefinitely with the (typically) exponential growth rate predicted in (19).

Finally, it is relevant to stress that the addition of the nonlinearity does not necessarily lead to a partial destabilization of the solutions, as the effect may be the opposite (a possibility of a partial stabilization of the symmetry under the action of nonlinearity was also recently demonstrated, in a different context, in Ref. (32). In particular, there are a few dots in Fig. 5 that do not appear in Fig. 1. Profiles of a few solutions in the linear and nonlinear systems, corresponding to three of these dots, are displayed in Fig. 8.

## Iii The cases of slow and fast modulations (adiabatic and averaging approximations)

In addition to the analysis presented above, some simple results can be obtained, like in other dynamical models (33), by means of the adiabatic approximation, which makes it possible to study some specific instabilities of the corresponding nonlinear systems (34). In both the adiabatic approximation and the averaging method for fast modulations given below, we will attempt to reveal the dependence of the solution on relevant parameters such as , , and .

The adiabatic approximation is valid for very small frequencies , when eigenvalues (8) are nearly degenerate, and

 ~V0(t)≡V0+2V1cos(ωt) (15)

is a slowly varying function of . Treating coefficient (15) as a constant, solutions of nonlinear system (1) with then follow from known results (18): setting

 ψA=e−iEtC,ψB=e−iEtC∗, (16)

one finds two solutions,

 C1 = √E−√~V20−1exp[i2arcsin(1~V0)], (17) C2 = √E+√~V20−1iexp[−i2arcsin(1~V0)]. (18)

These solutions are valid as long as conditions and hold at all , with the upper and lower signs corresponding to solutions (17) and (18), respectively. Combining these conditions for solution produces a validity range,

 1+2V1

which exists as long as . For the solution, the above conditions amount to

 V0>{2V1+1,  for  E>0,2V1+√1+E2,  for  E<0. (20)

The above solutions correspond to those reported in Ref. (18) by setting and . In the limit of , solution (18) is known to be stable, which is corroborated by our simulations of the system at small , with parameters from range (20). On the other hand, it is known too (18) that, at , solution (17) is stable only for satisfying an additional condition,

 ~V0≥√1+E2/4. (21)

In direct simulations performed at , the actually observed stability boundary rises to values exceeding the one (21). For example, if and , the lower existence bound for solution (17) in Eq. (19) is , whereas Eq. (21) yields , which is immaterial, as the existence threshold is higher. Nevertheless, even for as small as direct simulations reveal a higher instability boundary, at . Only at extremely small values of the numerically found the instability boundary drops down to its value (21) predicted for the static system.

The opposite case of fast modulations in system (1) implies the consideration of large , when the averaging method may be applied (35). The corresponding solution is easily obtained in the first approximation with respect to small :

 {ψAψB}={ψ(0)Aψ(0)B}−2iV1ωsin(ωt){ψ(0)Bψ(0)A}, (22)

where the unperturbed solution is given by Eqs. (16) and (17), (18) with (provided that condition holds). Thus, the respective final results are

 Missing or unrecognized delimiter for \left (23)
 {ψ(0)Aψ(0)B}=±ie−iEt√E+√V20−1⎧⎪⎨⎪⎩exp[i2arcsin(1V0)]−exp[−i2arcsin(1V0)]⎫⎪⎬⎪⎭. (24)

## Iv Conclusions

We have introduced the simplest model which opens the way to studying the concept of quasi-energy and PRs (parametric resonances), with ensuing instability, in -symmetric systems, that incorporate spatially separated and mutually balanced gain and loss elements. The parametric drive is introduced in the form of the periodic (ac) modulation of the coefficient of the coupling between the two degrees of freedom. The main results are reported for moderate values of driving frequency . The onset and development of the parametric instability at an infinitesimal amplitude of the ac drive was predicted analytically. The full picture was obtained by means of systematic simulations of both linear and nonlinear versions of the system. At large values of , the tongues of the parametric instability originate, with the increase of the driving amplitude (), at values of the coupling constant () predicted by the analytical consideration of the resonance points. However, the continuation of the PR tongues with the increase of in the -symmetric system follows a scenario which is very different from the usual parametrically driven system: instead of bending down to larger values of , the tongues continue in directions parallel to the axis. The decrease of makes the system of parallel tongues denser, eventually replacing the relatively simple system of parallel instability tongues by a complex small-scape structure. While the parametric resonance picture is an accurate predictor of the tongue formation at large , the quasi-energy method yields reasonably accurate results at small . The inclusion of the nonlinearity gradually deforms the picture, destabilizing a part of originally stable dynamical regimes, and stabilizing a few originally unstable ones. Nevertheless, as it might be expected, the chief effect of the nonlinearity is the expansion of the instability area. We have also presented a short account solutions obtained by means of the adiabatic and averaging approximations, for very slow and fast modulations, respectively.

This analysis may be extended in various ways, adding more terms to the system. In particular, a straightforward generalization may include nonlinear gain and loss terms, in addition to the linear pair of such terms, cf. Refs. (22); (23); (24). Extensions beyond the dimer case to oligomers and eventually to lattices would also be particularly interesting to consider (15); (16); (18); (19). In particular, discrete -symmetric solitons may be expected in the latter case. These themes are currently under consideration and will be presented in future publications.

## Acknowledgment

This work was supported, in a part, by the Binational (US-Israel) Science Foundation through grant No. 2010239. PGK gratefully also acknowledges support from NSF-DMS-1312856, NSF-CMMI-1000337 and AFOSR-FA9550-12-1-0332.

### References

1. C. M. Bender, Rep. Prog. Phys. 70, 947 (2007).
2. A. Ruschhaupt, F. Delgado, and J. G. Muga, J. Phys. A: Math. Gen. 38, L171 (2005).
3. 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üter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, Nature Phys. 6, 192 (2010).
4. N. Bender, S. Factor, J. D. Bodyfelt, H. Ramezani, D. N. Christodoulides, F. M. Ellis, and T. Kottos, Phys. Rev. Lett. 110, 234101 (2013); see also J. Schindler, A. Li, M. C. Zheng, F. M. Ellis and T. Kottos, Phys. Rev. A 84, 040101 (2011).
5. E. M. Graefe, U. Günther, H. J. Korsch, and A. E. Niederle, J. Phys. A: Math. Theor. 41, 255206 (2008); H. Cartarius and G. Wunner, Phys. Rev. A 86, 013612 (2012).
6. See special issues: H. Geyer, D. Heiss, and M. Znojil, Eds., J. Phys. A: Math. Gen. 39, Special Issue Dedicated to the Physics of Non-Hermitian 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. Bender, A. Fring, U. Günther, and H. Jones, Eds., Special Issue: Quantum Physics with non-Hermitian Operators, J. Math. Phys. A: Math Theor. 41, No. 44 (2012).
7. K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, Int. J. Theor. Phys. 50, 1019 (2011).
8. Z. H. Musslimani, K. G. Makris, R. El-Ganainy, and D. N. Christodoulides, Phys. Rev. Lett. 100, 030402 (2008); Z. Lin, H. Ramezani, T. Eichelkraut, T. Kottos, H. Cao, and D. N. Christodoulides, Phys. Rev. Lett. 106, 213901 (2011); X. Zhu, H. Wang, L.-X. Zheng, H. Li, and Y.-J. He, Opt. Lett. 36, 2680 (2011); C. Li, H. Liu, and L. Dong, Opt. Exp. 20, 16823 (2012); C. M. Huang, C. Y. Li, and L. W. Dong, ibid. 21, 3917 (2013).
9. S. Nixon, L. Ge, and J. Yang, Phys. Rev. A 85, 023822 (2012).
10. H. G. Li, Z. W. Shi, X. J. Jiang, and X. Zhu, Opt. Lett. 36, 3290 (2011); V. Achilleos, P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-González, Phys. Rev. A 86, 013808 (2012); see also V. Achilleos, P. G. Kevrekidis, D. J. Frantzeskakis, R. Carretero-González, Localized Excitations in Nonlinear Complex Systems, R. Carretero-Gonzlez, J. Cuevas-Maraver, D.ÊFrantzeskakis, N.ÊKarachalios, P.ÊKevrekidis, F.ÊPalmero-Acebedo, (Eds.) Springer-Verlag (Heidelberg, 2013).
11. F. C. Moreira, F. K. Abdullaev, V. V. Konotop, and A. V. Yulin, Phys. Rev. A 86, 053815 (2012); F. C. Moreira, V. V. Konotop, and B. A. Malomed, ibid. 87, 013832 (2013).
12. R. Driben and B. A. Malomed, Opt. Lett. 36, 4323 (2011); N. V. Alexeeva, I. V. Barashenkov, A. A. Sukhorukov, and Y. S. Kivshar, Phys. Rev. A 85, 063837 (2012); I. V. Barashenkov, S. V. Suchkov, A. A. Sukhorukov, S. V. Dmitriev, and Y. S. Kivshar, ibid. 86, 053809 (2012).
13. F. K. Abdullaev, V. V. Konotop, M. Ögren, and M. P. Sørensen, Opt. Lett. 36, 4566 (2011); R. Driben and B. A. Malomed, EPL 96, 51001 (2011); ibid. 99, 54001 (2012).
14. Yu. V. Bludov, V. V. Konotop, and B. A. Malomed, Phys. Rev. A 87, 013816 (2013).
15. S. V. Dmitriev, A. A. Sukhorukov, and Y. S. Kivshar, Opt. Lett. 35, 2976 (2010); S. V. Suchkov, B. A. Malomed, S. V. Dmitriev, and Y. S. Kivshar, Phys. Rev. E 84, 046609 (2011); S. V. Suchkov, A. A. Sukhorukov, S. V. Dmitriev, and Y. S. Kivshar, EPL 100, 54003 (2012); D. A. Zezyulin and V. V. Konotop, Phys. Rev. Lett. 108, 213906 (2012).
16. D. Leykam, V. V. Konotop, and A. S. Desyatnikov, Opt. Lett. 38, 371 (2013); I. V. Barashenkov, L. Baker, and N. V. Alexeeva, Phys. Rev. A 87, 033819 (2013).
17. J. Pickton and, H. Susanto, Phys. Rev. A 88, 063840 (2013).
18. K. Li and P. G. Kevrekidis, Phys. Rev. E 83, 066608 (2011); V. V. Konotop, D. E. Pelinovsky, and D. A. Zezyulin, EPL 100, 56006 (2012); J. D’Ambroise, P. G. Kevrekidis, and S. Lepri, J. Phys. A: Math. Theor. 45, 444012 (2012); K. Li, P. G. Kevrekidis, B. A. Malomed, and U. Günther, ibid. 45, 444021 (2012); M. Kreibich, J. Main, H. Cartarius, and G. Wunner, Phys. Rev. A 87, 051601(R) (2013).
19. P. G. Kevrekidis, D. E. Pelinovsky, D. Y. Tyugin, SIAM J. Appl. Dyn. Syst. 12 1210-36 (2013); see also P.G. Kevrekidis, D.E. Pelinovsky, and D.Y.Tyugin,ÊÊJ. Phys. A: Math. Theor.Ê46,Ê365201ÊÊ(2013).
20. H. Cartarius, D. Haag, D. Dast, and G. Wunner, J. Phys. A: Math. Theor. 45, 444008 (2012).
21. T. Mayteevarunyoo, B. A. Malomed, and A. Reoksabutr, Phys. Rev. E 88, 022919 (2013).
22. F. Kh. Abdullaev, Y. V. Kartashov, V. V. Konotop, and D. A. Zezyulin, Phys. Rev. A 83, 041805(R) (2011); D. A. Zezyulin, Y. V. Kartashov, V. V. Konotop, Europhys. Lett. 96, 64003 (2011).
23. A. E. Miroshnichenko, B. A. Malomed, and Y. S. Kivshar, Phys. Rev. A 84, 012123 (2011).
24. M. Duanmu, K. Li, R. L. Horne, P. G. Kevrekidis and N. Whitaker, Phil. Trans. R. Soc. A 371, 20120171 (2013).
25. Y. He, X. Zhu, D. Mihalache, J. Liu, and Z. Chen, Phys. Rev. A 85, 013831 (2012); Opt. Commun. 285, 3320 (2012).
26. E. Ott, Chaos in Dynamical Systems (Cambridge University Press: Cambridge, 1993).
27. B. A. Malomed, Physica D 29, 155 (1987); S. Fauve, O. Thual, Phys. Rev. Lett. 64, 282 (1990); W. van Saarloos and P. C. Hohenberg, Phys. Rev. Lett. 64, 749 (1990); V. Hakim, P. Jakobsen, and Y. Pomeau, Europhys. Lett. 11, 19 (1990); B. A. Malomed and A. A. Nepomnyashchy, Phys. Rev. A 42, 6009 (1990); P. Marcq, H. Chaté, and R. Conte, Physica D 73, 305 (1994); T. Kapitula and B. Sandstede, J. Opt. Soc. Am. B 15, 2757 (1998); A. Komarov, H. Leblond, and F. Sanchez, Phys. Rev. E 72, 025604 (2005); J. N. Kutz, SIAM Rev. 48, 629 (2006).
28. H. P. Breuer and M. Holthaus, Ann. Phys. 211, 249 (1991).
29. L. D.  Landau and E. M. Lifshitz, Mechanics (Nauka Publishers: Moscow, 1973).
30. P. L. Chu, B. A. Malomed, G. D. Peng, and I. Skinner, Phys. Rev. E 49, 5763 (1994).
31. T. Tel and M. Gruiz, Chaotic Dynamics (Cambridge University Press: New York, 2006).
32. Y. Lumer, Y. Plotnik, M. C. Rechtsman, and M. Segev, Phys. Rev. Lett. 111, 263901 (2013).
33. A. P. Itin, R. de la Llave, A. I. Neishtadt and A. A. Vasiliev, Chaos 12, 1043 (2001); Leoncini, A. Neishtadt, and A. Vasiliev, Phys. Rev. E 79, 026213 (2009).
34. A. P. Itin, S. Watanabe, and V. V. Konotop, Phys. Rev. A 77, 043610 (2008); A. P. Itin and S. Watanabe, Phys. Rev. Lett. 99, 223903 (2007).
35. V. Zharnitsky, I. Mitkov, and N. Grønbech-Jensen, Phys. Rev. E 58, R52 (1998); F. K. Abdullaev and R. A. Kraenkel, Phys. Rev. A 62, 023613 (2000); A. P. Itin, Phys. Rev. E 63, 028601 (2001); F. K. Abdullaev, J. G. Caputo, R. A. Kraenkel, and B. A. Malomed, Phys. Rev. A 67, 013605 (2003); F. K. Abdullaev and R. Galimzyanov, J. Phys. B: At. Mol. Opt. Phys. 36, 1099 (2003); A. P. Itin and A. I. Neishtadt, Phys. Rev. E 86, 016206 (2012).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters