Nonlinear Beam Propagation in a Class of Complex Non-PT -Symmetric Potentials

Nonlinear Beam Propagation in a Class of Complex Non-PT -Symmetric Potentials J. Cuevas-Maraver, P. G. Kevrekidis, D. J. Frantzeskakis, and Y. Kominis Abstract The subject of PT -symmetry and its areas of application have been blossoming over the past decade. Here, we consider a nonlinear Schr¨odinger model with a complex potential that can be tuned controllably away from being PT - symmetric, as it might be the case in realistic applications. We utilize two parameters: the first one breaks PT -symmetry but retains a proportionality between the imaginary and the derivative of the real part of the potential; the second one, detunes from this latter proportionality. It is shown that the departure of the potential from the PT -symmetric form does not allow for the numerical identification of exact stationary solutions. Nevertheless, it is of crucial importance to consider the dynamical evolution of initial beam profiles. In that light, we define a suitable notion of optimization and find that even for non PT -symmetric cases, the beam dynamics, both in 1D and 2D –although prone to weak growth or decay– suggests that the optimized profiles do not change significantly under propagation for specific parameter regimes. Jes´us Cuevas-Maraver Grupo de F´ısica No Lineal, Universidad de Sevilla, Departamento de F´ısica Aplicada I, Escuela Polit´ecnica Superior. C/ Virgen de ´ Africa, 7, 41011-Sevilla, Spain Instituto de Matem´aticas de la Universidad de Sevilla (IMUS). Edificio Celestino Mutis. Avda. Reina Mercedes s/n, 41012-Sevilla, Spain e-mail: jcuev[email protected] Panayotis G. Kevrekidis Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA 01003- 4515, USA Dimitrios J. Frantzeskakis Department of Physics, National and Kapodistrian University of Athens, Panepistimiopolis, Zografos, Athens 15784, Greece Yannis Kominis School of Applied Mathematical and Physical Science, National Technical University of Athens, Zographou GR-15773, Greece 1 2 J. Cuevas-Maraver, P. G. Kevrekidis, D. J. Frantzeskakis, and Y. Kominis Key words: Solitons, Nonlinear Schr¨odinger Equation, Stability, PT -symmetry, Unbalanced gain and loss, Symmetry breaking 1 Introduction The original suggestion of Bender and collaborators [1, 2] of a new class of systems that respect parity and time-reversal (so-called PT -symmetric systems) was motivated by the consideration of the foundations of quantum mechanics and the examination of the need of Hermitianity within them. The argument of Bender and collaborators was that such systems, even if non-Hermitian and featuring gain and loss, could give rise to real spectra, thus presenting a candidacy for being associated with measurable quantities. This proposal found a fertile ground for its development in areas, arguably, different than where it was originally proposed. In particular, the work of Christodoulides and co-workers in nonlinear optics a decade later, spearheaded an array of experimental realizations of such media (capitalizing on the ubiquitous in optics loss and on controllable gain) [3–8]. Other experiments swiftly followed in areas ranging from electronic circuits [9–11] to mechanical systems [12], bringing about not only experimental accessibility, but also an intense theoretical focus on this theme. These threads of research have now been summarized in two rather comprehensive recent reviews [7, 8]. While PT -symmetric variants of other nonlinear wave models have more recently been proposed, including the PT -symmetric variants of the Dirac equations [13] and of the Klein-Gordon equation [14], the main focus of associated interest has been on models of the nonlinear Schr¨odinger (NLS) type. This is natural given the relevance at the paraxial approximation level of such a model in applications stemming from nonlinear optics and related themes [7, 8]. In this important case, the PT -invariance is consonant with complex external potentials ˜ V, of the form ˜ V=V+iW, subject to the constraint that ˜ V∗(x) = ˜ V(−x). This implies that the real part, V, of the potential needs to be even, while the imaginary part, W, of the potential needs to be odd to ensure PT -symmetry. The expectation, thus, has been that typically Hamiltonian and PT -symmetric systems featuring gain and loss will possess continuous families of soliton solutions; otherwise, the models will possess solutions for isolated values within the parameter space. However, more recent investigations have started to challenge this belief. On the one hand, work on complex, asymmetric so-called Wadati potentials has produced mono-parametric continuous families of stationary solutions [15, 16]. On the other hand, the notion of partial PT -symmetry has been explored, e.g., with models that possess the symmetry in one of the directions but not in another [17, 18]. In fact, in the recent work of [19, 20] that motivated the present study, it was shown that to identify critical points one can localize a soliton 1in a way such that its intensity has 1Below, we use the term “soliton” in a loose sense, without implying complete integrability [21]. Nonlinear Beam Propagation in a Class of Complex Non-PT -Symmetric Potentials 3 a vanishing total overlap with the imaginary part of the potential, assuming that the real part of the potential is proportional to the anti-derivative of the imaginary part (but without making any assumptions on the parity of either). In the present work, we revisit these considerations. In particular, we discuss the results of the important contribution of [22]. This work suggests (and indeed conjectures) that the only complex potentials that could feature continuous families of stationary solutions although non-PT-symmetric are the ones of the Wadati type. In our case, we have considered potentials that depart from this form and either satisfy –or controllably depart from – the simpler mass and momentum balance conditions of [19,20]. We observe that in such settings, waveforms “optimizing” the vector field (which we define as bringing it very –but nor arbitrarily– close to vanishing) may exist, but still are not true solutions, in line with the above conjecture. We develop diagnostics that explore how these optimized beams behave dynamically, and identify their slow growth or decay. We do this for two different broad multi-parametric families of potentials to showcase the generality of our conclusions. We then extend relevant considerations also to 2D settings, showing how symmetry breaking bifurcation scenarios can be traced via our optimized beam approach. Our presentation will be structured as follows. In section 2, we introduce the model, connect our considerations to those of [22] and justify the selection of the complex potential. In section 3, we explore the optimized beams and the associated dynamics of the relevant waveforms numerically. Then, in section 4, we generalize these notions in a two-dimensional setting. In section 5, we proceed to summarize our findings and propose a number of directions for future study. Finally, in the Appendix, details of the numerical method used to optimize the dynamical beams are presented. 2 The one-dimensional potential As explained in the previous section, motivated by the development in the analysis of NLS models with complex potentials, we consider the rather broad setting of the form: i ψ t=− ψ xx + [V(x) + iW(x)] ψ − | ψ |2 ψ ,(1) with subscripts denoting partial derivatives. In the context of optics, ψ (x,t)represents the complex electric field envelope, tis the propagation distance, xcorresponds to the transverse direction, while the variation of the dielectric permittivity plays the role of the external potential, with V(x)and W(x)being its real and imaginary parts, respectively [7,8]. In the recent analysis of [19,20], assuming the existence of bright solitons (as is natural in the focusing nonlinearity setup under consideration), dynamical evolution equations were obtained for the soliton mass and velocity. Here, we use as our motivating point for constructing standing wave structures of Eq. (1) the stationary form of these equations, which read (cf. Eqs. (5)-(6) of Ref. [20]): 4 J. Cuevas-Maraver, P. G. Kevrekidis, D. J. Frantzeskakis, and Y. Kominis Z∞ −∞| ψ (x)|2W(x+x0)dx =Z∞ −∞| ψ (x)|2V′(x+x0)dx =0,(2) where x0denotes the soliton center. The first one among these equations, corresponds to a “power-balance” (or mass balance) condition, implying that the soliton has a transverse profile such that it experiences gain and loss in an overall balanced fashion across its spatial extent. The second equation corresponds to a “momentumbalance” condition, i.e., the total force exerted on the solitary wave vanishes, hence the coherent structure is at an equilibrium. This pair of stationarity conditions in Eq. (2) reduces to a single one, provided that V′=−CW, with Cbeing a constant. In that context, the resulting condition posits the following: if a soliton can be placed relative to the gain/loss profile so that its intensity has an overall vanishing overlap with the imaginary part of the potential, then the existence of a fixed point (and thus a stationary soliton solution) may be expected. However, it should be kept in mind that these conditions are necessary but not sufficient for the existence of a stationary configuration. In particular, a recent ingenious calculation shed some light on this problem for a general potential in the work of [22]. Using a standing wave deomposition ψ =r(x)eiRx θ (x′)dx′ei µ t in Eq. (1), the following ordinary differential equations were derived: rxx − µ r−Vr +r3− θ 2r=0,(3) (r2 θ )x=Wr2.(4) It was then realized that, in the absence of external potential, two quantities, namely J1=r2 θ (the “angular momentum” in the classical mechanical analogy of the problem) and J2=r2 x− µ r2+r4/2+r2 θ 2(the “first integral” or energy in the classical analogue) are conserved, i.e., dJi/dx =0. For J1, Eq. (4) yields its evolution in the presence of the potential while for J2, direct calculation shows: dJ2 dx =V(r2)x+2Wr2 θ =Sx−r2Vx−2(r2 θ )xZWdx,(5) with S=Vr2+2r2 θ RWdx. Combining the last terms, upon substitution of (r2 θ )x from Eq. (4) allows us to infer that this pair of terms will vanish if the coefficient multiplying r2, namely Vx−2WRWdx, vanishes; this occurs if the potential has the form: V+iW =−[g2+ig′(x)] + c, where cis a constant. A shooting argument presented in [22] suggests that there are 3 real constants (2 complex ones, yet one of them can be considered as real due to the phase invariance) in order to “glue” two complex quantities, namely ψ and ψ xat some point within the domain. This can only be done when a conserved quantity exists, which requires the type of potential suggested above, in the form Nonlinear Beam Propagation in a Class of Complex Non-PT -Symmetric Potentials 5 −[g2+ig′(x)]. However, if additional symmetry exists, such as PT-symmetry, the symmetry alone may impose conditions such as Im( ψ (0))=Re( ψ x(0)) = 0, which in turn allows for the shooting to go through (and thus solutions to exist) for a continuous range of µ ’s. Nevertheless, a natural question is: suppose that the potential is not of this rather non-generic form, yet it deviates from the PT-symmetric limit, possibly in ways respecting the above mass and/or momentum balance conditions of Eq. (2); then what is the fate of the system ? Do stationary states perhaps exist or do they not, and what are the dynamical implications of such conditions ? It is this class of questions that we will aim to make some progress towards in what follows. To test relevant ideas, we will use two different potentials ˜ Vj(x) =Vj(x)+iWj(x), with j=1,2. In the first one, Wis of the form: W1(x) = A1k1sech(x−xd− δ 1)tanh(x),(6) where A1,k1,xdand δ 1are constants, with the latter two controlling the breaking of the PT -symmetry. We then use a real potential V1given by the form: V1(x) = −2A1arctantanhxd−x 2coth(xd)−arctantanhx 2csch(xd), (7) which, in the limit xd→0, transforms into V1(x) = −A1sech(x). The motivation behind this selection is that if δ 1=0 in Eq. (6) then V1is proportional to the antiderivative of W1(hence ensures that the pair of conditions of Eq. (2) degenerate to a single one). In addition, for δ 1=0 and in the limit xd→0, the potential is PT - symmetric. In short, the two parameters xdand δ 1both control the departure from PT -symmetry, while the latter affects the departure from proportionality of V′ 1and W1. This selection and these parameters thus allow us to tailor the properties of the potential, controlling its departure from the PT -symmetric limit, but also from the possible degeneracy point of the conditions (2). The second potential is given by W2(x) = A2k2xsech2(x− δ 2−1),(8) and V2(x) = −A2(log[cosh(1−x)]+ xtanh(1−x)),(9) where A2,k2and δ 2are constants. Contrary to the ˜ V1case, this potential does not possess a PT -symmetric limit. For both ˜ V1(x)and ˜ V2(x)potentials, if δ =0 then V′ j(x) = −CjWj(x)and, as shown in Ref. [19], rendering a topic of interest the exploration of the potential existence of stationary solutions in the vicinity of the interface between the lossy and amplifying parts when Eq. (2) applies. In our particular case, the proportionality factor Cjis Cj=1/kj. 6 J. Cuevas-Maraver, P. G. Kevrekidis, D. J. Frantzeskakis, and Y. Kominis 3 Numerical results 3.1 Stationary states We start by seeking stationary localized solutions, of the form ψ (x,t) = ei µ tu(x) with u(x)∈C, which will thus satisfy: F[u]≡ µ u−uxx + [V(x) + iW(x)]u− |u|2u=0.(10) In what follows, we fix A1=0.1 and A2=1, and consider stationary solutions of frequency µ =1. We will make use of periodic boundary conditions. Notice that the potentials of [22] ˜ V(x) = −[g2(x) + ig′(x)] would, in the present notation, necessitate: [V′(x)]2=−4V(x)W2(x).(11) It is important to note that the potentials studied in our chapter do not fulfill this relation for any set of parameters (A1,k1,xd, δ 1)or (A2,k2, δ 2)—except for the “trivial” PT -limit— as it can be easily demonstrated. As a result then, presumably because of the above calculation, the standard fixed point methods that we have utilized fail to converge away from the PT -symmetric limit. For this reason, we make use of minimization algorithms in order to obtain optimized profiles of localized waveforms. With these methods, one can seek for local minima of the norm of F[u]instead of zeros of that function. In our problem, we have made use of the Levenberg–Marquardt algorithm (see Appendix A for more details), which has been successfully used for computing solitary gravity-capillary water waves [23], and established a tolerance of ||F[u]|| <10−3with ||F[u]|| being the L2-norm of F[u]: ||F[u]|| =rZ|F[u(x)]|2dx.(12) In the particular case of potential ˜ V1(x), we have studied the stability of solitons in the PT -symmetric limit xd= δ 1=0 as a function of k1, observing that solitons are stable whenever k1<kc, with kc=8.28. At this point, the soliton experiences a Hopf bifurcation. In order to avoid any connection of the findings below with the effect of such instability, we have fixed in what follows a value of k1far enough from kc. Moreover, since the minimal value attained for ||F[u]|| increases with k1, we have restricted consideration to relatively small values of k1and more specifically will report results in what follows for k1=1/2. Figures 1 and 2 show the potential profile for two different (xd, δ 1) and (k2, δ 2) parameter sets. These figures also show the profile of the waveforms minimizing ||F[u]|| for such potentials, which will be considered further in what follows. These beam profiles will be referred to as “optimized” in the sense of the above minimization. In particular, their real part is nodeless, while their imaginary part features a zero crossing. Naturally, the profiles are asymmetric mirroring the lack of definite parity of the potentials’ real and imaginary part. It is interesting to see that, despite Nonlinear Beam Propagation in a Class of Complex Non-PT -Symmetric Potentials 7 −10 −5 0 5 10 −0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 x V1(x) W1(x), δ1=0 W1(x), δ1=0.05 W1(x), δ1=−0.05 −10 −5 0 5 10 0 0.5 1 1.5 x Re[u(x)] δ1=0 δ1=0.05 δ1=−0.05 −10 −5 0 5 10 −0.02 −0.01 0 0.01 0.02 Im[u(x)] x −10 −5 0 5 10 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 x V1(x) W1(x), δ1=0 W1(x), δ1=0.5 W1(x), δ1=−0.25 −10 −5 0 5 10 0 0.5 1 1.5 x Re[u(x)] δ1=0 δ1=0.5 δ1=−0.25 −10 −5 0 5 10 −0.01 0 0.01 0.02 Im[u(x)] x Fig. 1 Left panels: Real and imaginary part of the potential ˜ V1(x)for A1=0.1, k1=1/2 and xd=0 (top) and xd=1 (bottom); green line corresponds to the imaginary part for δ 1=0, whereas red (black) line corresponds to the imaginary part for δ 1=0.05 ( δ 1=−0.05). Right panels: Beam profiles minimizing ||F[u]|| (real and imaginary parts) for A1=0.1, k1=1/2 and xd=0 (top) and xd=1 (bottom); the blue line corresponds to δ 1=0, and the green (red) line corresponds to δ 1=0.5 ( δ 1=−0.25). −10 −5 0 5 10 −0.5 0 0.5 1 1.5 2 x V2(x) W2(x), δ2=0 W2(x), δ2=0.2 W2(x), δ2=−0.1 −10 −5 0 5 10 0 0.5 1 1.5 x Re[u(x)] δ2=0 δ2=0.2 δ2=−0.1 −10 −5 0 5 10 −10 −5 0 5x 10−3 Im[u(x)] x Fig. 2 Left panel: Real and imaginary part of the potential ˜ V2(x)for A2=1, and k2=0.05; the green line corresponds to the imaginary part for δ 2=0, whereas the red (black) line corresponds to the imaginary part for δ 2=0.2 ( δ 2=−0.1). Right panels: optimized beam profiles (real and imaginary parts) for A2=1, and k2=0.01; the blue line corresponds to δ 2=0, and the green (red) line corresponds to δ 2=0.2 ( δ 2=−0.1). 8 J. Cuevas-Maraver, P. G. Kevrekidis, D. J. Frantzeskakis, and Y. Kominis the breaking of both the PT -symmetry and the violation of conditions such as the one in Eq. (11), there still exist spatially asymmetric structures almost satisfying the equations of motion. This naturally poses the question of the dynamical implications of such profiles in the evolution problem of Eq. (1), as we will see below. 3.2 Dynamics −20 −10 0 10 20 0 0.5 1 1.5 2 x |ψ(x)|2 t=0 t=500 t=1000 −20 −10 0 10 20 −8 −6 −4 −2 0 2 4x 10−5 x F[u(x)] Re(F[u(x)]) Im(F[u(x)]) t x 0 200 400 600 800 1000 −10 −5 0 5 10 0 0.5 1 1.5 0 200 400 600 800 1000 3.2 3.4 3.6 3.8 4 N(t) t Fig. 3 Optimized beam dynamics in the potential ˜ V1(x)for A1=0.1, k1=1/2, xd=1 and δ 1=0. The top left panel shows the density profile at t=0, t=500 and t=1000, while the top left panel shows the real and imaginary part of F[u]. The bottom left panel shows the space-time contour plot of the density evolution, and the bottom right panel shows the evolution of the norm N(t). The values of diagnostic quantities are λ =−3.63 ×10−3and σ =−1.07 ×10−4. We now analyze the dynamics of several case examples for the NLS equation with potential ˜ V1(x), using as initial condition the optimized beam profiles found by the Levenberg-Marquardt algorithm. Figures 3 and 4 show the outcome of the simulations for xd=1 and xd=−1, respectively, when δ 1=0 is fixed; on the other hand, Figs. 5 and 6 correspond, respectively, to δ 1=−0.1 and δ 1=0.1, when xd=1 is fixed. In these figures, we show the density | ψ (x)|2at different time instants (top left), the real and imaginary part of F[u](top right), a space-time contour plot of the evolution of the localized beam density | ψ (x,t)|2(bottom left), and the (squared) L2-norm (power/mass in optics/atomic physics), N(t)(bottom right), defined as Nonlinear Beam Propagation in a Class of Complex Non-PT -Symmetric Potentials 9 −20 −10 0 10 20 0 0.5 1 1.5 2 2.5 x |ψ(x)|2 t=0 t=500 t=1000 −20 −10 0 10 20 −4 −2 0 2 4 6 8x 10−5 x F[u(x)] Re(F[u(x)]) Im(F[u(x)]) t x 0 200 400 600 800 1000 −10 −5 0 5 10 0 0.5 1 1.5 2 0 200 400 600 800 1000 3.6 3.8 4 4.2 4.4 4.6 N(t) t Fig. 4 Same as Fig. 3, but for xd=−1. The values of diagnostic quantities are λ =3.63 ×10−3 and σ =1.07 ×10−4. −20 −10 0 10 20 0 1 2 3 4 5 6 x |ψ(x)|2 t=0 t=500 t=1000 −20 −10 0 10 20 −1 −0.5 0 0.5 1 1.5x 10−3 x F[u(x)] Re(F[u(x)]) Im(F[u(x)]) t x 0 200 400 600 800 1000 −10 −5 0 5 10 0 1 2 3 4 5 0 200 400 600 800 1000 3.5 4 4.5 5 5.5 6 6.5 7 N(t) t Fig. 5 Same as Fig. 3, but for δ 1=0.1. The values of diagnostic quantities are λ =2.07 ×10−2 and σ =6.07 ×10−4. 16 J. Cuevas-Maraver, P. G. Kevrekidis, D. J. Frantzeskakis, and Y. Kominis One can observe that the asymmetric branches exist for µ ≥6.3837. The figure also shows the profile of solitons at µ =7, the same value that was taken in [26]. Notice that all the soliton profiles are symmetric with respect to the y-axis; the symmetric solitons present a couple of humps at (x=±x1,y=0)whereas the asymmetric solitons only possess a single hump at (x=x2,y=0). We have only shown solitons with x2>0 as the solutions with x2<0 are attained simply by making the transform u(x,y)→u(−x,y). 6 6.5 7 7.5 8 0 1 2 3 4 5x 10−3 µ σ 6 6.5 7 7.5 8 0 0.5 1 1.5 2 µ λ Fig. 12 Dependence of diagnostic quantities σ (left) and λ (right) as a function of µ for optimized beams at the 2D potential ˜ V3(x,y)with A3=3 and k3=1. We have also computed the diagnostic quantities λ and σ [see (14) and Eqs. (15), with Sadapted to 2D domains] for the asymmetric soliton and depicted them in Fig. 12. Again, we have considered asymmetric solitons branches centred at x= x2>0. In that case, the norm grows with time, as corresponds to σ >0 and λ >0 whereas the opposite takes place if x2<0. We can observe, as in the 1D case, a clear correlation between both quantities. Finally, we show in Fig. 13 and 14 the dynamics of the asymmetric and PT - symmetric solitons with µ =7. As it was pointed out in [26], the PT -symmetric solitons are unstable past the “bifurcation” point, i.e. when they coexist with the asymmetric branch; as we have shown in Fig. 14, they tend to a state similar to the asymmetric soliton, although displaying some density oscillations. However, it was claimed in the same reference that the asymmetric solitons were stable. For the optimized beam profiles that we have obtained, as shown in Fig. 13, the dynamical evolution does not dramatically alter the shape of the beam, yet it leads to slow growth of N(t). We also considered the stability of the PT -symmetric branch past the relevant bifurcation point. A spectral stability analysis shows that for µ &6.40, the solitons become exponentially unstable as an eigenvalue pair becomes real. Interestingly, although the asymmetric solitons are actually optimized beams (i.e. solutions with minimal ||F(u)|| but not exact solutions), they might be more robust than the exact solutions of the NLS equation corresponding to the PT -symmetric branch, past Nonlinear Beam Propagation in a Class of Complex Non-PT -Symmetric Potentials 17 the corresponding destabilization point; compare the associated dynamics of Fig. 14 with those of Fig. 13. x y −4 −2 0 2 4 −4 −2 0 2 4 0.5 1 1.5 2 2.5 3 x y −4 −2 0 2 4 −4 −2 0 2 4 1 2 3 4 x y −4 −2 0 2 4 −4 −2 0 2 4 −10 −5 0 5 x 10−4 x y −4 −2 0 2 4 −4 −2 0 2 4 1 2 3 x 10−3 0 10 20 30 40 3.4 3.6 3.8 4 4.2 N(t) t Fig. 13 Optimized beam dynamics in the 2D potential ˜ V3(x,y)for A3=3, k3=1 and µ =7. The top panels show the density profile at t=0 (left) and t=40 (right). The bottom left panel shows the real and imaginary part of F[u]and the bottom right panel shows the evolution of the norm N(t). The values of diagnostic quantities are λ =1.39 and σ =4.07 ×10−3. 5 Conclusions & Future Work In the present work, we have revisited a variant of PT -symmetric systems. In particular, we have examined multi-parametric potentials whose parameters control, on the one hand, the potential departure from the PT -symmetric case (such as xd herein), and on the other hand, the potential degeneracy of the conditions (2) for stationary solutions –motivated by the recent works of [19,20]. We have confirmed the results of the important recent contribution of [22], suggesting that in the absence of a special form of the complex potential, no true stationary solutions are found to exist. On the other hand, that being said, we have identified beams that come very close to satisfying the stationary equations. The dynamics of these beams indicate a 18 J. Cuevas-Maraver, P. G. Kevrekidis, D. J. Frantzeskakis, and Y. Kominis t=0 t=10 t=15 t=20 t=25 t=30 t=35 t=40 0 10 20 30 40 4 5 6 7 8 9 N(t) t Fig. 14 Unstable PT -symmetric solitons dynamics in the 2D potential ˜ V3(x,y)for A3=3, k3=1 and µ =7. The top panels show snapshots of the density profile evolution. The bottom shows the evolution of the norm N(t). slow departure from such a configuration. In fact, diagnostics identifying the rate of this growth and connecting it to the proximity of the profiles to a stationary solution (via ||F[u]||) were developed and numerically evaluated, both in 1D and in 2D. Naturally, this work poses a number of questions for the future. One of the most notable such concerns the most general conditions (on, say, a complex potential) under which one may expect to find (or not) families of stationary solutions. PT - symmetry is a sufficient but not a necessary condition for such existence and extending beyond it seems of particular interest. The conjecture of [22] is that the potentials V+iW =−(g2+ig′(x))+crepresent the generic scenario is plausible, but it would be particularly interesting to produce a proof, perhaps revisiting more systematically the relevant shooting argument. It is also important to highlight that such shooting arguments are only valid in one spatial dimension. Hence, examining generalizations of the present setting to higher dimensions is of particular interest in their own right. We have briefly touched upon this aspect here, based on the earlier works of [25,26], but clearly further efforts are necessary to provide a definitive reply in this direction. In particular, while an interesting class of partially-PT -symmetric potentials has been examined herein, it is particularly relevant to attempt to understand Nonlinear Beam Propagation in a Class of Complex Non-PT -Symmetric Potentials 19 the phenomenology in non-PT -symmetric settings, constituting an important extension of our 1D considerations herein. Acknowledgments J.C.-M. thanks financial support from MAT2016-79866-R project (AEI/FEDER, UE). P.G.K. gratefully acknowledges the support of NSF-PHY-1602994, the Alexander von Humboldt Foundation, the Stavros Niarchos Foundation via the Greek Diaspora Fellowship Program, and the ERC under FP7, Marie Curie Actions, People, International Research Staff Exchange Scheme (IRSES-605096). The authors gratefully acknowledge numerous valuable discussions with and input from Professor Jianke Yang during the course of this work. Appendix: The Levenberg–Marquardt algorithm Classical fixed-point methods like Newton–Raphson cannot be used for solving the problem F[u(x)] = 0 in the setting considered in the context of this Chapter, essentially because there might not exist a u(x)that fulfils this relation (to arbitrarily prescribed accuracy). However, it is possible to find a function u(x)that can minimize F[u(x)]. To this aim, an efficient method is the Levenberg–Marquardt algorithm (LMA, for short), which is also known as the damped least-square method. This method is also used to solve nonlinear least squares curve fitting [27, 28]. LMA is implemented as a black box in the Optimization Toolbox of Matlab TM and in MINPACK library for Fortran, and can be considered as an interpolation between the Gauss-Newton algorithm and the steepest-descent method or viewed as a damped Gauss-Newton method using a trust region approach. Notice that LMA can find exact solutions, in case that they exist, as it is the case of the results presented, e.g., in Ref. [23]. Prior to applying LMA, we need to discretize our equation (10). Thus, we take a grid xn=−L/2+nh with n=0,1,2...Mand Lbeing the domain length, and denote un≡u(xn)and Fn≡F[u(xn)]. With this definition uxx can be cast as (un+1+un−1−2un)/h2. In order to simplify the notation in what follows, let us call u≡ {un}M n=1and F(u)≡ {Fn}M n=1. We will also need to define the Jacobian matrix J(u)≡ {Jn,m}M n,m=1with Jn,m= ∂ umFn. In the presently considered optimization framework, F(u)is also knows as the residue vector. Let us recall that fixed point methods typically seek a solution by performing the iteration uj+1=uj+ δ jfrom the seed u0until the residue norm ||F(u)|| is below the prescribed tolerance; here δ jis dubbed as the search direction. In the Newton–Raphson method, the search direction is the solution of the equation system J(uj) δ j=−F(uj). If the Jacobian is non-singular, the equation system can be easily solved (as a linear system); however, if this is not the case, one must look for 20 J. Cuevas-Maraver, P. G. Kevrekidis, D. J. Frantzeskakis, and Y. Kominis alternatives like the linear least square algorithm. It was successfully used for some of the authors for solving the complex Gross–Pitaevskii equation that describes the dynamics of exciton-polariton condensates [29–31]. This technique also allowed us to find optimized beams in the present problem, but presented poor convergence rates, as we were unable to decrease the residue norm controllably below the order of unity. As fixed point methods are unable to give a reasonably small residue norm, we decided to use a trust-region reflective optimization method. Such methods consist of finding the search direction that minimizes the so called merit function m( δ ) = 1 2F(u)TF(u)+ δ TJ(u)TF(u)+ δ TJ(u)TJ(u) δ .(26) In addition, δ must fulfill the relation ||D· δ || < ∆ ,(27) where Dis a scaling matrix and ∆ is the radius of the trust region where the problem is constrained to ensure convergence. There are several trust-region reflective methods, with the LMA being the one that has given us the best results for the problem at hand. This is a relatively simple method for finding the search direction δ by means of a Gauss-Newton algorithm (which is mainly used for nonlinear least squares fitting) with a scalar damping parameter λ >0 according to: (J(uj)TJ(uj)+ λ jD) δ j=−J(uj)TF(uj)(28) with Dbeing the scaling matrix introduced in Eq. (27). There are several possibilities for choosing such matrix. In the present work, we have taken the simplest option, that is D=I(the identity matrix), so (27) simplifies to || δ j|| < ∆ . Notice that for λ j=0, (28) transforms into the Gauss-Newton equation, while for λ j→∞ the equation turns into the steepest descent method. Consequently, the LMA interpolates between the two methods. Notice also the subscript in λ j: this is because the damping parameter must be changed in each iteration, with the choice of a suitable λ jconstituting the main difficulty of the algorithm. The scheme of the LMA is described in a quite easy way in Numerical Recipes book [32, Chapter 15.5.2] and is summarized below: 1. Take a seed u0and compute ||F(u0)|| 2. Choose a value for λ 0. In our particular problem, we have taken λ 0=0.1. 3. Solve the equation system (28) in order to get δ 0and compute ||F(u0+ δ 0)|| 4. •If ||F(u0+ δ 0)|| ≥ ||F(u0)||, then take λ 1=10 λ 0and u1=u0, as with this choice of λ 0the residue norm has not decreased. •If ||F(u0+ δ 0)|| <||F(u0)||, then take λ 1= λ 0/10 and u1=u0+ δ 0, as with this choice of λ 0has succeeded in decreasing the residue norm. 5. Go back to step 3 doing λ 0= λ 1and u0=u1 This algorithm is repeated while ||F(u)|| is above the prescribed tolerance. Nonlinear Beam Propagation in a Class of Complex Non-PT -Symmetric Potentials 21 References 1. Bender, C. M.; Boettcher, S. Real Spectra in Non-Hermitian Hamiltonians Having PT Symmetry, Phys. Rev. Lett. 80, 5243-5246 (1998). 2. Bender, C. M.; Brody, D. C.; Jones, H. F. Complex Extension of Quantum Mechanics, Phys. Rev. Lett. 89, 270401 (2002). 3. Ruter, C. E.; Makris, K. G.; El-Ganainy, R.; Christodoulides, D. N.; Segev, M.; Kip, D. Observation of parity-time symmetry in optics, Nat. Phys., 6, 192-195 (2010). 4. Peng, B.; Ozdemir, S. K.; Lei, F.; Monifi, F.; Gianfreda, M.; Long, G. L.; Fan, S.; Nori, F.; Bender, C. M.; Yang, L., Parity-time-symmetric whispering-gallery microcavities, Nat. Phys. 10, 394-398 (2014). 5. Peng, B.; Ozdemir, S. K.; Rotter, S.; Yilmaz, H.; Liertzer, M.; Monifi, F.; Bender, C. M.; Nori, F.; Yang, L., Loss-induced suppression and revival of lasing, Science 346, 328-332 (2014). 6. Wimmer, M.; Regensburger A.; Miri, M.-A.; Bersch, C.; Christodoulides, D.N.; Peschel, U.; Observation of optical solitons in PT-symmetric lattices, Nature Comms. 6, 7782 (2015). 7. Suchkov, S. V.; Sukhorukov, A. A.; Huang, J.; Dmitriev, S. V.; Lee, C.; Kivshar, Yu. S., Nonlinear switching and solitons in PT-symmetric photonic systems. Laser Photonics Rev. 10, 177-213 (2016). 8. Konotop, V. V.; Yang, J.; Zezyulin, D. A. Nonlinear waves in PT -symmetric systems, Rev. Mod. Phys. 88, 035002 (2016). 9. Schindler, J.; Li, A.; Zheng, M. C.; Ellis, F. M.; Kottos, T. Experimental study of active LRC circuits with PT symmetries. Phys. Rev. A 84, 040101 (2011). 10. Schindler, J.; Lin, Z.; Lee, J. M.; Ramezani, H.; Ellis, F. M.; Kottos, T., PT -symmetric electronics. J. Phys. A: Math. Theor. 45, 444029 (2012). 11. Bender, N.; Factor, S.; Bodyfelt, J. D.; Ramezani, H.; Christodoulides, D. N.; Ellis, F. M.; Kottos, T. Observation of Asymmetric Transport in Structures with Active Nonlinearities, Phys. Rev. Lett. 110, 234101 (2013). 12. Bender, C. M.; Berntson, B.; Parker, D.; Samuel, E. Observation of PT Phase Transition in a Simple Mechanical System. Am. J. Phys. 81, 173-179 (2013). 13. Cuevas-Maraver, J; Kevrekidis, P. G.; Saxena, A.; Cooper, F.; Khare, A.; Comech, A.; Bender, C. M. Solitary Waves of a PT -Symmetric Nonlinear Dirac Equation, IEEE J. Select. Top. Quant. Electron. 22, 5000109 (2016). 14. Demirkaya, A; Frantzeskakis, D. J.; Kevrekidis, P. G.; Saxena, A.; Stefanov, A. Effects of parity-time symmetry in nonlinear Klein-Gordon models and their stationary kinks, Phys. Rev. E 88, 023203 (2013); see also: Demirkaya, A; Kapitula, T.; Kevrekidis, P. G.; Stanislavova M.; Stefanov, A. On the Spectral Stability of Kinks in Some PT-Symmetric Variants of the Classical KleinGordon Field Theories, Stud. Appl. Math. 133, 298–317 (2014). 15. Tsoy, E. N.; Allayarov, I.M.; Abdullaev, F.Kh. Stable localized modes in asymmetric waveguides with gain and loss, Opt. Lett. 39, 4215–4218 (2014). 16. Konotop, V.V.; Zezyulin, D. A. Families of stationary modes in complex potentials, Opt. Lett. 39, 5535–5538 (2014). 17. Yang, J. Partially PT -symmetric optical potentials with all-real spectra and soliton families in multi-dimensions”, Opt. Lett. 39, 1133–1136 (2014). 18. D’Ambroise, J.; Kevrekidis, P. G. Existence, stability and dynamics of nonlinear modes in a 2D partially PT -symmetric potential, Appl. Sci. 7, 223 1–10 (2017). 19. Kominis, Y. Dynamic power balance for nonlinear waves in unbalanced gain and loss landscapes. Phys. Rev. A 92, 063849 (2015) 20. Kominis, Y. Soliton dynamics in symmetric and non-symmetric complex potentials, Opt. Commun. 334, 265–272 (2015). 21. Ablowitz M. J.; Segur, H. Solitons and the Inverse Scattering Transform (SIAM, Philadelphia, 1981). 22. Nixon S. D.; Yang, J. Bifurcation of Soliton Families from Linear Modes in Non-PT - Symmetric Complex Potentials. Stud. App. Maths. 136 (2016) 459. 22 J. Cuevas-Maraver, P. G. Kevrekidis, D. J. Frantzeskakis, and Y. Kominis 23. Dutykh, D.; Clamond, D.; Dur´an, ´ A. Efficient computation of capillary-gravity generalised solitary waves. Wave Motion 65 (2016) 1. 24. Johansson, M.; Aubry, S. Growth and decay of discrete nonlinear Schr¨odinger breathers interacting with internal modes or standing-wave phonons, Phys. Rev. E 61, 5864–5879 (2000). 25. Yang, J. Symmetry breaking of solitons in two-dimensional complex potentials. Phys. Rev. E91 (2015) 023201. 26. Chen, H.; Hu, S. The asymmetric solitons in two-dimensional parity-time symmetric potentials. Phys. Lett. A 380 (2016) 162. 27. Levenberg, K. A Method for the Solution of Certain Non-Linear Problems in Least Squares. Quarterly of Applied Mathematics 2(1944) 164. 28. Marquardt, D. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. SIAM Journal on Applied Mathematics 11 (1963) 431 29. Cuevas, J.; Rodrigues, A.S., Carretero-Gonz´alez, R.; Kevrekidis, P.G.; Frantzeskakis, D.J. Nonlinear excitations, stability inversions, and dissipative dynamics in quasi-one- dimensional polariton condensates. Physical Review B 83 (2011) 245140. 30. Rodrigues, A.S.; Kevrekidis, P.G.; Cuevas, J.; Carretero-Gonz´alez, R.; Frantzeskakis, D. J. Symmetry-breaking effects for polariton condensates in double-well potentials. In B. A. Malomed (Ed.), Spontaneous Symmetry Breaking, Self-Trapping, and Josephson Oscillations. Springer Verlag. 2013, (pp. 509-529) 31. Rodrigues, A.S.; Kevrekidis, P.G.; Carretero-Gonz´alez, R.; Cuevas-Maraver, J.; Frantzeskakis, D. J.; Palmero, F. From nodeless clouds and vortices to gray ring solitons and symmetry-broken states in two-dimensional polariton condensates. Journal of Physics: Condensed Matter 26 (2014) 155801. 32. Press, W.H.; Teukolsky, S.A.; Vetterling, W.T.; Flannery, B.P. Numerical recipes. The art of scientific computing. Third Edition. Cambridge University Press (2007).