Existence, Stability and Dynamics of Nonlinear Modes in a 2D PartiallyPT Symmetric Potential

It is known that multidimensional complex potentials obeying parity-time(PT)symmetry may possess all real spectra and continuous families of solitons. Recently, it was shown that for multi-dimensional systems, these features can persist when the parity symmetry condition is relaxed so that the potential is invariant under reﬂection in only a single spatial direction. We examine the existence, stability and dynamical properties of localized modes within the cubic nonlinear Schrodinger equation in such a scenario of partiallyPT-symmetric potential.


Introduction
The study of P T (parity-time) symmetric systems was initiated through the works of Bender and collaborators [1,2].Originally, it was proposed as an alternative to the standard quantum theory, where the Hamiltonian is postulated to be Hermitian.In these works, it was instead found that Hamiltonians invariant under P T -symmetry, which are not necessarily Hermitian, may still give rise to completely real spectra.Thus, the proposal of Bender and co-authors was that these Hamiltonians are appropriate for the description of physical settings.In the important case of Schrödinger-type Hamiltonians, which include the usual kinetic-energy operator and the potential term, V(x), the P T -invariance is consonant with complex potentials, subject to the constraint that V * (x) = V(−x).
A decade later, it was realized (and since then it has led to a decade of particularly fruitful research efforts) that this idea can find fertile ground for its experimental realization although not in quantum mechanics where it was originally conceived.In this vein, numerous experimental realizations sprang up in the areas of linear and nonlinear optics [3][4][5][6][7][8], electronic circuits [9][10][11], and mechanical systems [12], among others.Very recently, this now mature field of research has been summarized in two comprehensive reviews [7,8].
One of the particularly relevant playgrounds for the exploration of the implications of P T -symmetry is that of nonlinear optics, especially because it can controllably involve the interplay of P T -symmetry and nonlinearity.The relevant efforts have gradually progressed from the simpler setting of coupled waveguides (bearing gain and loss) to entire P T -symmetric lattices and the identification of optical solitons in them.These developments hold considerable promise for the potential realization of more complex settings, such as the one proposed herein.In this context, the propagation of light (in systems such as optical fibers or waveguides [7,8]) is modeled by the nonlinear Schrödinger equation of the form: (1) In the optics notation that we use here, the evolution direction is denoted by z, the propagation distance.Here, we restrict our considerations to two spatial dimensions and assume that the potential U(x, y) is complex valued, representing gain and loss in the optical medium, depending on the sign of the imaginary part (negative for gain, positive for loss) of the potential.In this two-dimensional setting, the condition of full P T -symmetry in two dimensions is that U * (x, y) = U(−x, −y).Potentials with full P T symmetry have been shown to support continuous families of soliton solutions [13][14][15][16][17][18].However, an important recent development was the fact that the condition of (full) P T symmetry can be relaxed.That is, either the condition U * (x, y) = U(−x, y) or U * (x, y) = U(x, −y) of, so-called, partial P T symmetry can be imposed, yet the system will still maintain all real spectra and continuous families of soliton solutions [19].Other subsequent results include stable vortex structures in cubic nonlinear media under the influence of partially P T -symmetric azimuthal potentials [20].
In the original contribution of [19], only the focusing nonlinearity case was considered for two select branches of solutions and the stability of these branches was presented for isolated parametric cases (of the frequency parameter of the solution).Our aim in the present work is to provide a considerably more "spherical" perspective of the problem.The corresponding physical scenario that we have in mind involves a medium featuring a cubic nonlinearity.Given that the potential of [19] was principally featuring four "nodes", we are envisioning an implementation involving a four waveguide system.This is both in the spirit of the pioneering (two waveguide) work of [3] and also in that of subsequent proposals for trimers and quadrimers (i.e., three and four waveguide systems, respectively) [21,22].A more remote possibility is that of a continuous cubic nonlinearity medium bearing a complex index of refraction in accordance with the prescription of [19], although, to the best of our understanding, such a possibility appears less tractable on the basis of current experimental capabilities.
In what follows, we examine the bifurcation of nonlinear modes from all three point spectrum eigenvalues of the underlying linear Schrödinger operator of the partially P T -symmetric potential.Upon presenting the relevant model (Section 2), we perform the relevant continuations (Section 3) unveiling the existence of nonlinear branches both for the focusing and for the defocusing nonlinearity case.We also provide a systematic view towards the stability of the relevant modes (Section 4), by characterizing their principal unstable eigenvalues as a function of the intrinsic frequency parameter of the solution.In Section 5, we complement our existence and stability analysis by virtue of direct numerical simulations that manifest the result of the solutions' dynamical instability when they are found to be unstable.Finally, in Section 6, we summarize our findings and present our conclusions, as well as discussing some possibilities for future studies.

Model, Theoretical Setup and Linear Limit
Motivated by the partially P T -symmetric setting of [19], we consider the complex potential U(x, y) = V(x, y) + iW(x, y) where with real constants β, a = b and c = −d.The potential is chosen with partial P T -symmetry so that U * (x, y) = U(−x, y).That is, the real part is even in the x-direction with V(x, y) = V(−x, y) and the imaginary part is odd in the x-direction with −W(x, y) = W(−x, y).The constants a, b, c, d are chosen such that there is no symmetry in the y direction.
In [19], it is shown that the spectrum of the potential U can be all real as long as |β| is below a threshold value, after which a (P T -) phase transition occurs; this is a standard property of P T -symmetric potentials.For the case of a = 3, b = c = 2, d = 1, the spectrum is real below the threshold value of |β| ≈ 0.214; we focus on β = 0.1, i.e., we operate well below this critical point.Figure 1 shows plots of the potential U. The real part of the potential is shown on the left, while the imaginary part associated with gain-loss is on the right; the gain part of the potential corresponds to W < 0 and occurs for x < 0, while the loss part with W > 0 occurs for x > 0. Figure 2 shows the spectrum of U, i.e., eigenvalues for the underlying linear Schrödinger problem (∇ 2 + U)ψ 0 = µ 0 ψ 0 .The figure also shows the magnitude of the corresponding eigenvectors for the three discrete real eigenvalues µ 0 .It is from these modes that we will seek bifurcations of nonlinear solutions in what follows.It is worthwhile to mention here, in comparison, e.g., with the real four-well potential of [23] that the latter possessed four localized modes, with the fourth antisymmetric, quadrupolar one being absent from the point spectrum in the case of interest herein.

Existence: Nonlinear Modes Bifurcating from the Linear Limit
As is customary, we focus on stationary soliton solutions of (1) of the form Ψ(x, y, z) = ψ(x, y)e iµz .Thus, one obtains the following stationary equation for ψ(x, y).
In [19], it is discussed that a continuous family of solitons bifurcates from each of the linear solutions in the presence of nonlinearity.In order to see this, let µ 0 be a discrete simple real eigenvalue of the potential U (such as one of the positive real eigenvalues in the top left plot of Figure 2).Now, following [19], expand ψ(x, y) in terms of = |µ − µ 0 | << 1 and substitute the expression into Equation (3).This gives the equation for ψ 1 as where ρ = sgn(µ − µ 0 ) and Here, ψ * 0 plays the role of the adjoint solution to ψ 0 .Thus, in order to find solutions of (3) for σ = ±1, we perform a Newton continuation in the parameter µ where the initial guess for ψ is given by the first two terms of (4).The bottom left panel of Figure 3 shows how the (optical) power P(µ) = |ψ| 2 dxdy of the solution grows as a function of increasing µ for σ = 1, or as a function of decreasing µ for σ = −1 (from the linear limit).The first branch begins at the first real eigenvalue of U at µ 0 ≈ 0.286, the second branch at µ 0 ≈ 0.487, and the third branch begins at µ 0 ≈ 0.785.Plots of the solutions and their corresponding time evolution and stability properties are shown in the next section.The curves begin at the lowest power (i.e., at the linear limit) at the discrete real eigenvalues of approximately 0.286 (branch 1), 0.487 (branch 2), 0.785 (branch 3).Each power curve is drawn with its corresponding stability noted: a blue solid curve denotes a stable solution and a red dashed curve denotes an unstable solution.The other three plots track the maximum real part of eigenvalues ν as a function of the continuation parameter µ: the red dashed line represents the max real part of eigenvalues that are real (exponential instability) while the blue dotted line tracks the max real part for eigenvalues that have a nonzero imaginary part (quartets); this case corresponds to oscillatory instabilities.
As a general starting point comment for the properties of the branches, we point out that all the branches populate both the gain and the loss side.In the branch starting from µ 0 ≈ 0.286, all four "wells" of the potential of Figure 1 appear to be populated, with the lower intensity "nodes" being more populated and the higher intensity ones less populated.The second branch starting at µ 0 ≈ 0.487, as highlighted also in [19], possesses an anti-symmetric structure in x (hence the apparent vanishing of the density at the x = 0 line).Both in the second and in the third branch, the higher intensity nodes of the potential appear to bear a higher intensity.

Stability of the Nonlinear Modes: Spectral Analysis
The natural next step is to identify the stability of the solutions.This is monitored by using the linearization ansatz: which yields the order δ linear system where Thus, max(Re(ν)) > 0 corresponds to instability and max(Re(ν)) = 0 corresponds to (neutral) stability.
In the bottom left panel of Figure 3, the power curve is drawn with stability and instability as determined by ν noted by the solid or dashed curve, respectively.The other three plots in Figure 3 show the maximum real part of eigenvalues ν for each of the three branches: the red dashed curve is the max real part of real eigenvalue pairs, and the blue dotted curve is the max real part of eigenvalue quartets with a nonzero imaginary part.The former corresponds to exponential instabilities associated with pure growth, while the latter indicate so-called oscillatory instabilities, where growth is present concurrently with oscillations.In Figure 4, we plot some example eigenvalues in the complex plane for some sample unstable solutions.The dominant unstable eigenvalues within these can be seen to be consonant with the growth rates reported for the respective branches (and for these parameter values) in Figure 3.One can compare the maximal real part with Figure 3.For example, the top left complex plane plot here shows that for branch 1 at µ ≈ 0.71, the eigenvalues with the maximum real part are complex; this agrees with the top left plot of Figure 3 where at µ ≈ 0.71, the blue dotted curve representing complex eigenvalues is bigger.Similarly, one can check whether the other three eigenvalue plots here also agree with what is shown in Figure 3, the top right for branch 1, the bottom left for branch 2 and the bottom right for branch 3.
The overarching conclusions from this stability analysis are as follows.The lowest µ branch, being the ground state in the defocusing case, is always stable in the presence of the self-defocusing nonlinearity.For the parameters considered, generic stability is also prescribed for the third branch under self-defocusing nonlinearity.The middle branch has a narrow interval of stability and then becomes unstable, initially (as shown in the top right of Figure 3) via an oscillatory instability and then through an exponential one.In the focusing case (that was also focused on in [19] for the second and third branch), all three branches appear to be stable immediately upon their bifurcation from the linear limit, yet all three of them subsequently become unstable.Branch 1 (that was not analyzed previously) features a combination of oscillatory and exponential instabilities.Branch 2 features an oscillatory instability which, however, only arises for a finite interval of frequencies µ, and the branch restabilizes.On the other hand, for branch 3, when it becomes unstable, it is through a real pair.Branches 2 and 3 terminate in a saddle-center bifurcation near µ = 1.9.The eigenvalue panels of Figure 4 confirm that the top panels of branch 1 may possess one or two concurrent types of instability (in the focusing case), branch 2 (bottom left) can only be oscillatorily unstable in the focusing case (yet as is shown in Figure 3, it can feature either type of instability in the defocusing case), while for branch 3, when it is unstable in the focusing case, it is via a real eigenvalue pair.

Dynamics of Unstable Solutions
Figures 5-7 show the time evolution of three unstable solutions, one on each branch, in the focusing case for the value of σ = 1.That is, they each correspond to a µ-value that is bigger than the initial discrete value µ 0 and pertain to the focusing case.The time evolution figures show a similar feature for the unstable solutions, namely that over time the magnitude of the solutions will increase on the left side of the spatial grid.This agrees with what is expected from P T -symmetry since the left side of the spatial grid corresponds to the gain side of the potential U. Importantly, also, the nature of the instabilities varies from case to case, and is consonant with our stability expectations based on the results of the previous section.
In Figure 5, branch 1 (for the relevant value of the parameter µ) features an oscillatory instability (but with a small imaginary part).In line with this, we observe a growth that is principally exponential (cf.also the top panel for the power of the solution), yet features also some oscillation in the amplitude of the individual peaks.It should be noted here that although two peaks result in growth and two in decay (as expected by the nature of W in this case), one of them clearly dominates between the relevant amplitudes.In Figure 6, it can be seen that branch 2, when unstable in the focusing case, is subject to an oscillatory instability (with a fairly significant imaginary part).Hence the growth is not pure, but is accompanied by oscillations as is clearly visible in the top left panel.In this case, among the two principal peaks of the solution of branch 2, only the left one (associated with the gain side) is populated after the evolution shown.Lastly, in branch 3, the evolution (up to z = 42) manifests the existence of an exponential instability.The latter leads, once again, to the indefinite growth of the gain part of the solution, resulting in one of the associated peaks growing while the other (for x > 0 on the lossy side) features decay.It is worthwhile to mention that in the case of branch 2-the only branch that was found (via our eigenvalue calculations) to be unstable in the self-defocusing case-we also attempted to perform dynamical simulations for σ = −1.Nevertheless, in all the cases considered, it was found that, fueled by the defocusing nature of the nonlinearity, a rapid spreading of the solution would take place (as z increased), resulting in the interference of the wave pattern with the domain boundaries.In Figure 8, we show an example where the maximum real part of the eigenvalues is ≈ 0.016.This weak instability is of exponential type with eigenvalues having a nonzero real part lying on the real axis similar to the eigenvalues in the bottom right panel of Figure 4 (but very close to the origin).Figure 8 shows that the solution initially seems to follow the expected pattern of dominating on the gain side (left) and there is a power growth similar to the focusing cases, but accompanied by a quick increase of the spatial extent of the solution leading to interference with the computational domain boundary.5 (with an evolution up to distance z = 300).Here, the plots correspond to the evolution of the branch 2 solution for the value µ ≈ 0.244.In the top right plot, the colors follow the same pattern as in Figure 5 representing the four peaks.While the gain side of the solution seems to dominate, we also observe a quick increase in the spatial footprint of the solution.Shortly after the plotted time interval, the solution interferes with the boundary of the computed spatial grid.

Conclusions and Future Challenges
In the present work, we have revisited the partially P T -symmetric setting originally proposed in [19] and have attempted to provide a systematic analysis of the existence, stability and evolutionary dynamics of the nonlinear modes that arise in the presence of such a potential for both self-focusing and self-defocusing nonlinearities.It was found that all three linear modes generate nonlinear counterparts.Generally, the defocusing case was found to be more robustly stable than the focusing one.In the former, two of the branches were stable for all the values of the frequency considered, while in the focusing case, all three branches developed instabilities sufficiently far from the linear limit (although all of them were spectrally stable close to it).The instabilities could be of different types, both oscillatory (as for branch 2) and exponential (as for branch 3) or even of mixed type (as for branch 1).The resulting oscillatorily or exponentially unstable dynamics, respectively, led to the gain overwhelming the dynamics and leading to indefinite growth in the one or two of the gain peaks of our four-peak potential.
Naturally, there are numerous directions that merit additional investigation.For instance, and although this would be of less direct relevance in optics, partial P T symmetry could be extended to three dimensions.There, it would be relevant to appreciate the differences between potentials that are partially P T symmetric in one direction vs. those partially P T symmetric in two directions.Another relevant case to explore in the context of the present mode is that where a P T phase transition has already occurred through the collision of the second and third linear eigenmode considered herein.Exploring the nonlinear modes and the associated stability in that case would be an interesting task in its own right.Such studies are presently under consideration and will be reported in future publications.

Figure 1 .
Figure 1.The plots show the spatial distribution of real (V, left panel) and imaginary (W, right panel) parts of the potential U with x 0 = y 0 = 1.5.

Figure 2 .
Figure 2. The top left plot shows the spectrum of Schrödinger operator associated with the potential U in the complex plane (see also the text).Plots of the magnitude of the normalized eigenvectors for the three discrete eigenvalues µ 0 are shown in the other three plots.

Figure 3 .
Figure 3.The bottom left plot shows the power of the solution ψ plotted in terms of the continuation parameter µ.The curves begin at the lowest power (i.e., at the linear limit) at the discrete real eigenvalues of approximately 0.286 (branch 1), 0.487 (branch 2), 0.785 (branch 3).Each power curve is drawn with its corresponding stability noted: a blue solid curve denotes a stable solution and a red dashed curve denotes an unstable solution.The other three plots track the maximum real part of eigenvalues ν as a function of the continuation parameter µ: the red dashed line represents the max real part of eigenvalues that are real (exponential instability) while the blue dotted line tracks the max real part for eigenvalues that have a nonzero imaginary part (quartets); this case corresponds to oscillatory instabilities.

Figure 4 .
Figure 4. Eigenvalues are plotted in the complex plane (Re(ν), Im(ν)) for a few representative solutions.One can compare the maximal real part with Figure3.For example, the top left complex plane plot here shows that for branch 1 at µ ≈ 0.71, the eigenvalues with the maximum real part are complex; this agrees with the top left plot of Figure3where at µ ≈ 0.71, the blue dotted curve representing complex eigenvalues is bigger.Similarly, one can check whether the other three eigenvalue plots here also agree with what is shown in Figure3, the top right for branch 1, the bottom left for branch 2 and the bottom right for branch 3.

Figure 5 .
Figure 5.This figure shows the time evolution of the branch 1 solution for the value µ ≈ 0.71.The bottom left plot shows the magnitude of the solution |Ψ| at z = 0. Observe that this solution has four peaks in its magnitude over the two-dimensional spatial grid.The bottom right plot shows the solution at z = 23.Observe that the magnitudes of the peaks on the left side have increased.The top left plot shows the time evolution of the power of the solution as a function of the evolution variable z.The top right plot here shows the evolution of the four peaks in the magnitude of the solution as a function of z (blue = bottom left peak, red = top left peak, green = bottom right peak, cyan = top right peak).

Figure 6 .
Figure 6.This figure is similar to Figure 5 (the final evolution distance however is about z = 114).Here, the plots correspond to the time evolution of the branch 2 solution for the value µ ≈ 0.94.In the top right plot, the blue curve corresponds to the left peak of the magnitude of the solution over z and the green corresponds to the right peak of the magnitude.

Figure 7 .
Figure 7.This figure is similar to Figure5(with an evolution up to distance z = 42).Here, the plots correspond to the evolution of the branch 3 solution for the value µ ≈ 1.0.In the top right plot, the blue curve corresponds to the left peak of the magnitude of the solution over z and green corresponds to the right peak of the magnitude.Clearly, once again, the gain side of the solution eventually dominates.

Figure 8 .
Figure 8.This figure is similar to Figure5(with an evolution up to distance z = 300).Here, the plots correspond to the evolution of the branch 2 solution for the value µ ≈ 0.244.In the top right plot, the colors follow the same pattern as in Figure5representing the four peaks.While the gain side of the solution seems to dominate, we also observe a quick increase in the spatial footprint of the solution.Shortly after the plotted time interval, the solution interferes with the boundary of the computed spatial grid.