# Parity-time symmetric optical coupler with birefringent arms

###### Abstract

In this work, we propose a parity-time (-) symmetric optical coupler whose arms are birefringent waveguides as a realistic physical model which leads to a so-called quadrimer i.e., a four complex field setting. We seek stationary solutions of the resulting linear and nonlinear model, identifying its linear point of symmetry breaking and examining the corresponding nonlinear solutions that persist up to this point, as well as, so-called, ghost states that bifurcate from them. We obtain the relevant symmetry breaking bifurcations between symmetric (circularly polarized) and asymmetric (elliptically polarized) states and numerically follow the associated dynamics which give rise to growth/decay even within the -symmetric phase. Our symmetric stationary nonlinear solutions are found to terminate in saddle-center bifurcations which are analogous to the linear -phase transition. We found that the symmetry significantly changes the stability and dynamical properties of the modes with different polarizations.

###### pacs:

42.65.Jx, 42.65.Tg, 42.65.Wi## I Introduction

An optical coupler with one arm having losses and another one having gain, balanced against each other, recently became a test-bed for many phenomena originating from the interplay of the parity-time () symmetry and nonlinearity. Unidirectional dynamics kot1 (), unversality of the dynamics kot2 (), symmetry breaking properties kip (); pgk (), switching of the beams sukh1 () and of solitons R30add3 (), formation of symmetric and asymmetric bright solitary waves R30add1 (); R30add2 (), breathers baras (), and their stability R30add4 (), dark solitons BKM (), as well as the emergence of ghost states R46a (); R46b (); R46 () and large-scale temporal -symmetric lattices christo_nat2 () are some among the many topics that have been touched upon in very intense recent theoretical and experimental work.

As a direct extension of the previous activity, a large chunk of which has focused on the prototypical setting of the -symmetric dimer, there emerges a problem of effect of symmetry and nonlinearity on the polarization of the electric field. In that regard, the previously proposed settings, to the best of our knowledge, were chiefly focused on effectively scalar models. On the other hand, the “vector” type of problems is natural for experimental settings where the exploited fibers obey birefringence, since the two orthogonal polarizations are to be taken into account Menyuk ().

In this work we consider a -symmetric coupler whose arms are birefringent waveguides. Assuming that the first waveguide is active and the second one is absorbing, we address the problem of a -preserving (in the linear limit) configuration. While being an interesting model from a physical point of view, this setting also offers a different (in comparison to what was studied before) mathematical situation where the nonlinear modes bifurcate from doubly degenerated eigenvalues of the linear problem. This requires the generalization of earlier developed approaches (e.g. like the one reported in ZK ()) for the bifurcation of the nonlinear modes from the linear spectrum. In addition, it presents a rich playground for dynamical systems analysis, due to the emergence of a variety of saddle-center bifurcations (nonlinear analogs of the linear -phase transition), as well as symmetry-breaking (pitchfork) ones. It is these nonlinear states, their emergence, stability, dynamics and the asymptotics of the system that we will focus on hereafter.

The organization of the paper is as follows. In section II, we present the model in its evolution as well as in its stationary form. In section III, we focus on its linear properties. Then, in section IV, we examine the nonlinear modes (and bifurcations). Finally, in section V, we briefly touch upon the dynamical implications of our findings and in section VI, we present some conclusions, as well as some potential directions for future work.

## Ii The model

We specify the problem by imposing that the principal optical axes of the two Kerr-type waveguides are -rotated with respect to each other, as it is schematically represented in Fig. 1. In each arm, labeled by , there are two orthogonal field components of the electric fields which we write down in the form Menyuk () ():

(1) |

Here are the field envelopes depending on the propagation distance , i.e. we consider the stationary – in time – problem, assuming that the carrier wavelength is in the region of the normal group velocity dispersion, thus ruling out a possibility of modulational instability; is a transverse radius vector, and are the positions of the centers of the cores of the coupler. The real parameters are the propagation constants of each of the field components, and are the polarization vectors, which are mutually orthogonal in each arm of the coupler, i.e. . The real functions and describe the transverse distributions of the fields in each waveguide and the normalization coefficient , where is the Kerr coefficient, is introduced for convenience. For and the functions are centered in different points . Also, for the sake of simplicity, we consider (for ), such that the integral

(the integration is performed over the transverse plane) describes the linear coupling between the respective modes. Since in the configuration shown in Fig. 1 we use the single linear coupling coefficient (see also coupling ()).

Then following the analysis described in details in Menyuk () we end up with the system of equations:

(2a) | |||

(2b) | |||

(2c) | |||

(2d) |

Here describes gain in the first waveguide and dissipation in the second waveguide, with , being the carrier wave frequency, is a properly normalized mismatch between the propagation constants of the orthogonal polarizations and . The asterisk stands for complex conjugation.

We will be interested in the stationary solutions, in particular in their linear stability properties and ensuing nonlinear dynamics which can be found in the two prototypical limiting cases of (i) zero mismatches and (ii) large mismatches when the respective nonlinear terms can be neglected. It is convenient to introduce a parameter which vanishes () in the case (ii) and is unity () in the case (i). Using the standing wave ansatz , where are independent, we obtain the system of algebraic equations:

(3a) | |||

(3b) | |||

(3c) | |||

(3d) |

Below the spectral parameter will be also referred to as the propagation constant.

## Iii Properties of the linear problem

First we address the underlying linear problem [which corresponds to the situation when all cubic terms in Eqs. (3) are negligible]. It can be rewritten in the matrix form where

(hereafter we use tilde in order to distinguish eigenvalues and eigenvectors of the linear problem).

The operator is symmetric, which means that , where is a spatial reversal linear operator

(4) |

and performs element-wise complex conjugation: . The spectrum of operator consists of two double eigenvalues

(5) |

which are real for where will be referred to as a primary critical point: the spontaneous symmetry breaking occurs at above which the eigenvalues are all imaginary. In order to visualize the symmetry of the linear system, following ZK () one can represent it with a graph shown in the right panel of Fig. 1, reminiscent of four linearly coupled waveguides ZK () (notice however the sign difference in the coupling constants) or plaquettes guenther ().

Since the details of our analysis are the same for both eigenvalues and , we drop the subscripts and wherever this does not lead to confusion. Despite having double eigenvalues in its spectrum, is diagonalizable below the -symmetry breaking point. This means that double eigenvalues are semisimple, i.e. for an eigenvalue one can find two linearly independent eigenvectors, i.e. , where . Moreover, each eigenvalue possesses an invariant subspace spanned by and .

Let us also notice the following peculiarity of the case at hand. In a situation where a -symmetric operator has no multiple eigenvalues, the condition of unbroken symmetry (i.e. reality of all the eigenvalues) requires that for each eigenvalue the corresponding eigenvector can be chosen as an eigenstate of the operator, i.e. . However, in the situation at hand arbitrarily chosen linearly independent eigenvectors and may not be eigenstates. However, unbroken symmetry requires that a certain linear combination of and is an eigenstate for the operator. More specifically, it is easy to establish that all the eigenvectors that belong to the invariant subspace of and, at the same time, are the eigenstates for the operator, can be parametrized by a complex parameter as follows:

(6) |

Being interested in linearly independent vectors , it is sufficient to consider only the vectors with

(7) |

Then Eq. (6) yields a monoparametric set of eigenvectors with a real parameter . In particular, setting and one can choose two orthogonal (and therefore linearly independent) eigenvectors:

(8) |

(hereafter we use the standard scalar product ). Any eigenvector corresponding to the eigenvalue can be represented as a linear combination of eigenstates : (of course, this does not mean that any eigenvector is also a eigenstate).

Let us also introduce a Hermitian adjoint operator . Since the matrix is symmetric, one has . As long as symmetry of is unbroken, the spectrum of the adjoint operator also consists of two double eigenvalues which are semi-simple. Any eigenvector corresponding to an eigenvalue of the adjoint operator can be represented as a linear combination of and .

## Iv Nonlinear modes

### iv.1 Bifurcations from the linear limit

Now we develop a perturbation theory for the eigenstates of the linear problem giving rise to monoparametric families of nonlinear modes. We will look for nonlinear modes that are eigenstates of the operator, i.e. . To this end we introduce the expansions

(9) |

Here is a small real parameter, and are the coefficients of the expansions and is to be determined from the symmetry of the solution (see below). We notice that for the expansion to be meaningful the coefficient must be real.

Expansions (9) describe nonlinear modes that bifurcate from the linear limit corresponds to and is given by the eigenvector being a linear combination of such that . Respectively, in the linear limit the parameter is given by the eigenvalue .

Passing from to one has to compute the coefficients and . While the physical sense of the coefficient is clear — it is a deviation of the propagation constant due to small nonlinearity, it turns out that also has a clear geometrical interpretation. Indeed, let us consider the total energy flow through the coupler, which is defined by

(10) |

Expansions (9) imply that in the vicinity of the linear limit . Therefore, the coefficient governs a slope of the energy curve in the vicinity of the bifurcation point, i.e.

(11) |

For the sake of definiteness, now we concentrate on the case . Then the nonlinear problem (3) is conveniently written in the matrix form

(12) |

where is a diagonal matrix-function describing the nonlinearity: . Substituting (9) into Eq. (12), noticing that , and collecting the terms order of , we obtain

(13) |

Equation (13) implies two possibilities. The first one corresponds to the case when at some the eigenvector of the operator is simultaneously an eigenvector for the matrix . Then the coefficient can be chosen as an eigenvalue of corresponding to the eigenvector (provided that this eigenvalue is real). In this situation, the right hand side of Eq. (13) is zero and it is sufficient to set . Since the matrix is diagonal, its eigenvalues are equal to its diagonal elements and the corresponding eigenvectors are given as columns of the identity matrix.

Let us first assume that all the eigenvalues of are simple. In this case can not be an eigenvector for . The latter fact becomes evident if one notices that has no zero entries for any [see the definition (6)]. Therefore, can be an eigenvector for only if has a multiple eigenvalue. Then could be searched in the form of a linear combination of eigenvectors corresponding to the multiple eigenvalue. However, using the same argument, i.e. the fact that all entries of the vector are nonzero, one can see that even if has a double or an triple eigenvalue, the matrix still can not have among its eigenvectors. Therefore can be an eigenvector of only if all its eigenvalues are equal. Imposing this constraint on the matrix , one obtains and . Noticing that the form of implies, through Eqs. (6) and (7), that we conclude that is an eigenvector of only if the moduli of all the entries of are equal to unity. Then matrix is equal to the identity matrix multiplied by . Requiring the moduli of all the entries of to be equal, one arrives at the equation for whose root is given as

(14) |

Since now the moduli of all the entries of are equal to unity and therefore , and Eq. (11) readily yields that in the vicinity of the slope . Notice that the found value does not depend on or , and thus these modes correspond to their counterpart in pure conservative coupler with birefringent arms with [Eq. (14) is valid in this case since ].

Let us now consider the second possibility to fulfill Eq. (13). If for some the corresponding is not an eigenvector for , then one must satisfy Eq. (13) choosing nonzero . Then the coefficient is to be determined from the solvability condition which requires the right hand side of Eq. (13) to be orthogonal to all the eigenvectors of the invariant subspace of in the spectrum of the adjoint operator . As we have established in Sec. III, any eigenvector of from the invariant subspace of can be represented as a linear combination of and . Requiring the right hand side of Eq. (13) to be orthogonal to an arbitrary linear combination of and , we arrive at the following relations:

(15) |

In Eq. (15) the first equality sign is an equation which is to be solved with respect to . Once a root of the latter equation is found, then is given from the second equality sign. Notice that despite the fact that the vector is complex, the coefficient will be real ZK ().

Substituting the expression for into Eq. (15), one obtains a rather cumbersome equation, which, however can be attacked with a computer algebra program. After some transformations, Eq. (15) yields the following condition:

(16) |

The latter equation has a real root only if , i.e. when the expression under the radical is not positive. This result suggests that there exists a critical value of the gain-loss parameter which we term as the secondary critical point , such that for sufficiently small , namely, , there exists another family bifurcating from the eigenvalue of the linear spectrum. However, this family disappears for , in spite of the fact that the symmetry of the underlying linear problem remains unbroken, i.e. . It is important to point that the newly found family of solutions does not correspond to equal amplitude among the different nodes, and hence pertains to an elliptically (rather than circularly) polarized family of modes.

Computing the corresponding value of for the newly found elliptically polarized family, one finds that the coefficient does depend on and (in contrast to the above considered circularly polarized family, characterized by for any combination of and which does not violate symmetry). In particular, when approaches , the coefficient of the elliptically polarized family tends to , which suggests that at the circularly polarized and elliptically polarized families globally merge.

### iv.2 Algebraic analysis and numerical results

#### iv.2.1 Exact solutions

Having explored nonlinear modes close to the linear limit, where the amplitudes of the modes are small and therefore they can be analyzed by means of perturbation theory, let us now consider nonlinear modes of arbitrary amplitudes (turning again to the general case ). Relying on results of the previous subsection we firstly search for nonlinear modes which have equal intensities in all four waveguides. Making the substitution and imposing the condition , , which is necessary for a nonlinear mode to be an eigenstate of thus leading to the circularly polarized light in each of the coupler arms, system (3) yields the single (complex) algebraic equation

(17) |

Representing , we obtain a bi-quadratic equation for yielding two families of the modes bifurcating from the eigenvalues , given by (5), of the linear spectrum:

(18) |

Respectively, the nonlinear modes have the following form:

(19) |

Using Eqs. (18) one can easily obtain continuous families of nonlinear modes that can be identified as a function of the propagation constant , for given and .

Let us now recall that Eq. (16) predicts that for and , there exist families of elliptically polarized modes having different absolute values of the polarization vectors. Such families were indeed found in our numerics. However, all such modes turned out to be unstable (see Figs. 2–3 and discussion below).

In the case of zero propagation constant mismatch, i.e. when , one also can find families which have different amplitudes of the polarization vectors. The explicit expressions for the families bifurcating from read

(20) |

Remarkably, these modes, which also describe propagation of elliptically polarized light, are stable in a certain range of the parameters.

#### iv.2.2 Families of nonlinear modes

The results of our analysis of the families of nonlinear modes are summarized in Fig. 2. The upper panels show that for each eigenvalue of the linearized problem gives rise to two distinct (circularly and elliptically polarized) families of nonlinear modes. In the case of the slopes of the dependencies are close for the families of both types, and the elliptically polarized families are always unstable while circularly polarized families have both stable and unstable solutions. For one can find stable solutions both for the families with circular and for those with elliptical polarization.

For , the case does not allow for elliptically polarized families [see Fig. 2 C], a feature which is in accordance with the perturbation approach developed above. In this case one can only find circularly polarized modes, which are unstable. On the other hand, for , stable and unstable modes of both types can be found [see Fig. 2 D].

Summarizing at this point, we have identified 4 sets of solutions, two circularly polarized with equal amplitude at the nodes, and two elliptically polarized with unequal such amplitudes. These all degenerate into the two distinct eigenvalues , given by (5), of the linear problem. The circularly polarized solutions are more robust, while the elliptically polarized ones are always unstable for and stable only for small enough amplitudes for . Among the circularly polarized ones, for the more fundamental state (stemming from the negative eigenvalue at the linear limit) is always the stable ground state of the system in continuations over the parameter , while the excited state is only stable for small enough amplitudes.

#### iv.2.3 Continuation over

An alternative and perhaps even more telling way to illustrate the above features stems from fixing some value of , starting from the Hamiltonian limit of and subsequently identifying branches of the nonlinear modes by means of changing , as shown in Fig. 3. It is important to note that this alternative viewpoint affords us the ability to visualize bifurcations which we now explore.

The relevant results for parametric continuations over are given in Figs. 3-4; typical examples of the corresponding linearization spectra for different values of can be found in Fig. 5. Here, it can be seen that a lower amplitude and a higher-amplitude intensity-symmetric (i.e., equal amplitude) branch exist, for fixed , from the Hamiltonian limit of and all the way up to the linear -phase transition point . At that point, the two equal amplitude branches collide and disappear in a saddle-center bifurcation which can be thought of as a nonlinear analog of the linear -phase transition djf (). An additional very interesting feature arises precisely at the point [see (16) and the related discussion], where it can be seen that both branches of equal amplitude between the sites become dynamically unstable for . In fact, it is seen that for the larger amplitude branch (associated with the blue circles), one pair of unstable eigenvalues arises, while for the smaller amplitude (red diamond) branch, two such pairs accompany the symmetry breaking bifurcation occurring at this critical point. A closer inspection reveals that the symmetric branch (blue circles) is destabilized through a subcritical pitchfork bifurcation with its “corresponding” asymmetric state (i.e., the one degenerate with it in the linear limit). In the case of the lower amplitude (excited) state for the same , the situation appears to be more complex. In particular, there exists once again a subcritical pitchfork with the corresponding asymmetric branch, yet this would justify one pair of unstable eigenvalues and we observe two. This is because at the same point, there also exists a supercritical pitchfork, which gives rise to the so-called ghost states, denoted by magenta plus symbols. These states are analogous to the ones to analyzed in R46a (); R46b (); R46 (), but remarkably are not stationary states of the original problem, yet they are pertinent to its dynamical (instability) evolution and for this reason they will be examined in further detail separately in the dynamics section below.

In the case of , only one pair of unstable eigenvalues emerges for the lower amplitude branch at the secondary critical point of (while the larger amplitude branch remains stable throughout the continuation in ). Hence, in this case, once again a saddle-center bifurcation will mark the nonlinear -phase transition, yet the number of unstable eigendirections of each symmetric branch (fundamental and excited) is decreased by one (0 and 1 real pairs instead of 1 and 2, respectively, for ). In this case, in fact, both asymmetric branches persist up to the linear -phase transition (rather than terminate in a subcritical pitchfork as above), and collide and disappear with each other. Interestingly all 3 branches (the lower amplitude, excited symmetric one and the two asymmetric ones) become unstable at the secondary critical point , which again points to the existence of corresponding ghost states. For the lower amplitude symmetric branch, the bifurcating ghost states are again identified by the magenta plus symbols in Fig. 4.

## V Dynamics of the polarization

To examine the potentially symmetry breaking (and more generally instability driven) nature of the dynamical evolution past the critical points identified above, we have also performed direct numerical simulations which are illustrated in Fig. 6; see also Fig. 7. Here, it can be seen that although the relevant parameters are below the critical point for the linear -phase transition , nevertheless, symmetry breaking phenomena are observed due to the dynamical instability of the relevant states (the ones denoted by dashed lines in Fig. 2). This dynamics may, in principle, be associated with the so-called ghost states of complex propagation constant that have recently been proposed as relevant for the dynamical evolution in R46 (). To substantiate this claim, we note that it is observed in the left panel of Fig. 6 that the relative phase of the two gain sites that lock into an equal growing amplitude, is , as is those of the decaying amplitude lossy sites. In light of this, we seek ghost states with precisely this phase difference and are able to explicitly identify them via the ansatz , setting for . For these branches, the propagation constant is complex. This highlights the potential growth or decay of such states. Importantly also, note that these states are “ghosts” because they may be solving the stationary problem of Eqs. (3), but the U invariance of the original model does not permit them to be a solution of the dynamical Eqs. (2).

The algebraic conditions that this family of solutions satisfies are

(21) | ||||

(22) | ||||

(23) | ||||

(24) | ||||

(25) | ||||

(26) |

Notice that the imaginary part of the propagation constant is proportional to the difference . Hence, prior to the symmetry breaking, the relevant solutions bear a real propagation constant. Past the bifurcation point one (unstable) branch has , while the stable branch has . The relevant ghost state branches and their bifurcation from the equal amplitude ones are explored in Fig. 3-4. Given that these are only ghost solutions of the original dynamical problem, the interpretation of their linearization spectrum (shown for completeness in Fig. 5) is still an open problem.

These ghost states appear, in fact, to exhibit very similar evolution dynamics to regular unstable states. To illustrate this, we observed the particular behavior of the unstable modes and how it depends on the form of the initial perturbation. A typical example in which the gain sites lead to growth and the lossy sites to decay is shown in the left panel of Fig. 6. It is interesting that the evolution appears to be very proximal to that of the ghost states identified above. This is clearly showcased in Fig. 7, through the comparison of the growth pattern in the relevant sites (and the decay pattern in the lossy sites) with the exact, shifted in the propagation distance to fit the onset of growth, ghost state solution for the same parameters. On the other hand, in the right panel of Fig. 6 a different scenario of evolution is illustrated. Instead of the gain nodes growing and the lossy ones decaying, a breathing oscillation settles between the two pairs. These two scenarios, illustrated in Fig. 6, are the prototypical instability evolution ones that we have obtained in this system.

## Vi Conclusion

In conclusion, in the present work, we have proposed a novel, physically realistic variant of a symmetric dimer where the effect of birefringence has been taken into consideration. The existence of polarization of the electric field within the coupler yields two complex dynamical equations for each of the fibers, providing a physical realization of a plaquette model with both linear and nonlinear coupling between the elements. The stationary states of the model were identified and both linear and nonlinear -phase transitions were obtained. The degenerate nature of the linear limit complicated the problem in comparison to other ones studied earlier in this context. Furthermore, the emergence of symmetry breaking phenomena and associated (subcritical or supercritical) pitchfork bifurcations, as well as their dynamical implications in leading to indefinite growth and decay (of the corresponding waveguide amplitudes) were elucidated. A connection was also given to ghost states.

From the physical point of view, we emphasize that the use of symmetry significantly changes the stability of the modes with different polarizations. For example, in the case of large propagation constant mismatches, we saw that the symmetric (circularly polarized) states become destabilized in the presence of gain/loss and may even cease to exist past a certain critical point of the relevant parameter. Instead of them, the dynamics may lead to a breathing exchange of “mass” between the gain and the lossy waveguides, or most commonly an indefinite growth of the former at the expense of the latter. This is in line with what is observed also in the evolution of the emerging ghost states of the system. The instabilities and associated dynamics should be observable in suitable generalizations of existing experiments such as kip (). Such properties are relevant also to possibilities of solitonic (waveguide-array) generalizations of the coupler system as well as towards a possible use of symmetric coupler for measurement techniques based on the use of several nonlinear modes.

From the mathematical point of view, we believe that these studies may pave the way for considering multi-component, as well as multi-dimensional (generalizing the plaquettes considered here or those of guenther ()) lattice models of -symmetric form. In generalizing to multi-plaquette configurations, it would be especially interesting to examine which of the symmetry-breaking and nonlinear -phase transition phenomena examined herein are preserved and what new phenomena may arise as additional degrees of freedom are added. Such studies will be deferred to future publications.

###### Acknowledgements.

VVK and DAZ acknowledge support of the FCT grants PTDC/FIS/112624/2009, PEst-OE/FIS/UI0618/2011, and SFRH/BPD/64835/2009. PGK gratefully acknowledges support from the Alexander von Humboldt Foundation, and the Alexander S. Onassis Public Benefit Foundation, as well as from the US-NSF (grants DMS-0806762, CMMI-1000337) and from the US-AFOSR.## References

- (1) H. Ramezani, T. Kottos, R. El-Ganainy and D. N. Christodoulides, Phys. Rev. A 82, 043803 (2010).
- (2) M. C. Zheng, D. N. Christodoulides, R. Fleischmann and T. Kottos, Phys. Rev. A 82, 010103(R) (2010).
- (3) C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, D. Kip, Nature Phys. 6, 192 (2010).
- (4) K. Li and P. G. Kevrekidis, Phys. Rev. E 83, 066608 (2011).
- (5) A. A. Sukhorukov, Z. Xu and Yu. S. Kivshar, Phys. Rev. A 82, 043818 (2010).
- (6) F. Kh. Abdullaev, V.V. Konotop, M. Ögren and M. P. Sørensen, Opt. Lett. 36, 4566 (2011).
- (7) R. Driben and B. A. Malomed, Opt. Lett. 36, 4323 (2011).
- (8) R. Driben and B. A. Malomed, Europhys. Lett. 96, 51001 (2011).
- (9) I. V. Barashenkov, S. V. Suchkov, A. A. Sukhorukov, S. V. Dmitriev, and Yu. S. Kivshar, Phys. Rev. A 86, 053809 (2012).
- (10) N. V. Alexeeva, I. V. Barashenkov, A. A. Sukhorukov, and Yu. S. Kivshar, Phys. Rev. A 85, 063837 (2012).
- (11) Yu. V. Bludov, V. V. Konotop, and B. A. Malomed, Phys. Rev. A 87, 013816 (2013).
- (12) H. Cartarius and G. Wunner, Phys. Rev. A 86, 013612 (2012); H. Cartarius, D. Haag, D. Dast and G. Wunner, J. Phys. A 45, 444008 (2012);
- (13) E.-M. Graefe, J. Phys. A 45, 444015 (2012).
- (14) A. S. Rodrigues, K. Li, V. Achilleos, P. G. Kevrekidis, D. J. Frantzeskakis, C. M. Bender, arXiv:1207.1066.
- (15) A. Regensburger, C. Bersch, M.-A. Miri, G. Onishchukov, D. N. Christodoulides and U. Peschel, Nature 488, 167 (2012).
- (16) C. R. Menyuk IEEE J. Quant. Eelectron., 23, 174-176 (1987)
- (17) D. A. Zezyulin and V. V. Konotop, Phys. Rev. Lett. 108 213906 (2012).
- (18) A. W. Snyder and Y. Chen, Opt. Lett. 14, 517 (1989); Y. Chen, A. W. Snyder, and D. N. Payne, IEEE J. Quant. Electro. 28, 239 (1992); M. Romagnoli S. Trillo, S. Wabnitz, Opt. and Quant. Electron. 24, S1237 (1992).
- (19) K. Li, P.G. Kevrekidis, B. A. Malomed and U. Guenther, J. Phys. A: Math. Theor. 45 444021 (2012).
- (20) V. Achilleos, P.G. Kevrekidis, D.J. Frantzeskakis, and R. Carretero-González Phys. Rev. A 86, 013808 (2012).