Nonlinear Beam Propagation in a Class of Complex Non-PT -Symmetric Potentials
J. Cuevas-Maraver P. G. Kevrekidis D. J. Frantzeskakis Y. Kominis
nlin.PS


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




Jesús Cuevas-Maraver
Grupo de Física No Lineal, Universidad de
Sevilla, Departamento de Física Aplicada I,
Escuela Politécnica Superior.
C/ Virgen de África, 7, 41011-Sevilla, Spain 

Instituto de Matemáticas de la Universidad de Sevilla (IMUS). Edificio Celestino Mutis. Avda. Reina Mercedes s/n, 41012-Sevilla, Spain
jcuevas@us.es
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


*
    

J. Cuevas-Maraver
  
P. G. Kevrekidis

D. J. Frantzeskakis

Y. Kominis

    December 30, 2023
=============================================================================



The subject of 𝒫𝒯-symmetry and its areas
  of application have been blossoming over the past decade. Here, we
  consider a nonlinear Schrödinger model with a complex potential that
  can be tuned controllably away from being 𝒫𝒯-symmetric, as it
  might be the case in realistic applications.
  
  We utilize two parameters: the first one breaks 𝒫𝒯-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 𝒫𝒯-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 𝒫𝒯-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.






§ INTRODUCTION


The original suggestion of Bender and collaborators <cit.>
of a new class of systems that respect parity and time-reversal
(so-called 𝒫𝒯-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) <cit.>.
Other experiments swiftly followed in areas ranging from
electronic circuits <cit.> to
mechanical systems <cit.>, 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 <cit.>.

While 𝒫𝒯-symmetric variants of
other nonlinear wave models have more recently been proposed, including
the 𝒫𝒯-symmetric variants of the Dirac equations <cit.>
and of the Klein-Gordon equation <cit.>, the main focus of
associated interest has been on models of the nonlinear Schrödinger (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 <cit.>. In this
important case, the 𝒫𝒯-invariance is consonant with
complex external potentials Ṽ, of the form Ṽ=V+i W,
subject to the constraint that Ṽ^∗(x)=Ṽ(-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 𝒫𝒯-symmetry.
The expectation, thus, has been that typically Hamiltonian
and 𝒫𝒯-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 <cit.>.
On the other hand, the notion of partial 𝒫𝒯-symmetry
has been explored, e.g., with models that possess the symmetry
in one of the directions but not in another <cit.>.
In fact, in the recent work of <cit.> that motivated
the present study, it was shown that to identify critical points
one can localize a soliton
[Below, we use the term “soliton” in a loose sense, without
implying complete integrability <cit.>.]
in a way such that its intensity has
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 <cit.>. This work suggests (and indeed
conjectures) that the only complex potentials that could
feature continuous families of stationary solutions
although non-𝒫T-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 <cit.>.
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 <cit.> 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.




§ 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)+i W(x)]ψ-|ψ|^2ψ,


with subscripts denoting partial derivatives. In the context
of optics, ψ(x,t) represents the complex electric field
envelope, t is the propagation distance, x corresponds 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
<cit.>. In the recent analysis of <cit.>,
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. (<ref>)
the stationary form of
these equations, which read
(cf. Eqs. (5)-(6) of Ref. <cit.>):


    ∫_-∞^∞ |ψ(x)|^2 W(x+x_0) dx=
      ∫_-∞^∞ |ψ(x)|^2 V'(x+x_0) dx=0,


where x_0 denotes 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 “momentum-balance” 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. (<ref>) reduces to a single
one, provided that V'=-C W, with C being 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 <cit.>. Using a standing wave deomposition

    ψ=r(x)e^i ∫^x θ(x') dx' e^i μ t
 in Eq. (<ref>), the following
ordinary differential equations were derived:

    r_xx - μ r - V r + r^3 - θ^2 r = 0,
    
          (r^2 θ)_x = W r^2.

It was then realized that, in the absence of external potential,
two quantities, namely J_1=r^2 θ (the “angular momentum”
in the classical mechanical analogy of the problem) and
J_2=r_x^2 - μ r^2 + r^4/2 + r^2 θ^2 (the “first integral”
or energy in the classical analogue) are conserved, i.e., 
d J_i/dx=0. For J_1, Eq. (<ref>) yields its evolution
in the presence of the potential while for J_2, direct calculation
shows:

    dJ_2/dx=V (r^2)_x + 2 W r^2 θ = S_x - r^2 V_x-2 (r^2 θ)_x
      ∫ W dx,

with S=V r^2 +2 r^2 θ∫ W dx. Combining the last terms,
upon substitution of (r^2 θ)_x from Eq. (<ref>)
allows us to infer that this pair of terms will vanish if the coefficient
multiplying r^2, namely V_x -2 W ∫ W dx, vanishes; this occurs
if the potential has the form:

    V+i W=- [g^2 + i g'(x)] + c,

where c is a constant.
A shooting argument presented in <cit.> 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 ψ_x
at 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 -[g^2 + i g'(x)]. However,
if additional symmetry exists, such as 𝒫T-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 𝒫T-symmetric limit, possibly in ways respecting
the above mass and/or momentum balance conditions of Eq. (<ref>);
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
 Ṽ_j(x)=V_j(x)+iW_j(x), with j=1,2. In the first one, W is of the form:


    W_1(x)=A_1k_1sech(x-x_d-δ_1)tanh(x),


where A_1, k_1, x_d and δ_1 are constants, with the
latter two controlling the breaking of the 𝒫𝒯-symmetry.
We then use a real potential V_1 given by the form:


    V_1(x)=-2A_1[arctan(tanh(x_d-x/2))(x_d)-
    arctan(tanh(x/2))csch(x_d)],


which, in the limit x_d→0, transforms into V_1(x)=-A_1sech(x).
The motivation behind this selection is that
if δ_1=0 in Eq. (<ref>) then
V_1 is proportional to the anti-derivative of W_1
(hence ensures that the pair of conditions of Eq. (<ref>)
degenerate to a single one). In addition, for δ_1=0 and in the limit x_d→0, the potential is 𝒫𝒯-symmetric. In short,
the two parameters x_d and δ_1 both control the departure from
𝒫𝒯-symmetry, while the latter affects the departure from
proportionality of V'_1 and W_1.
This selection and these parameters
thus allow us to tailor the properties of the potential,
controlling its departure from the 𝒫𝒯-symmetric limit,
but also from the possible degeneracy point of the conditions (<ref>).

The second potential is given by


    W_2(x)=A_2k_2xsech^2(x-δ_2-1),


and


    V_2(x)=-A_2(log[cosh(1-x)]+xtanh(1-x)),


where A_2, k_2 and δ_2 are constants. Contrary to the Ṽ_1 case, this potential does not possess a 𝒫𝒯-symmetric limit.

For both Ṽ_1(x) and Ṽ_2(x) potentials, if δ=0 then V'_j(x)=-C_jW_j(x) and, as shown in Ref. <cit.>,
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. (<ref>) applies.
In our particular case, the proportionality factor
C_j is C_j=1/k_j.



§ NUMERICAL RESULTS




 §.§ Stationary states


We start by seeking
stationary localized solutions, of the form ψ(x,t)=e^i μ tu(x) with u(x)∈ℂ,
which will thus satisfy:


    F[u]≡μ u-u_xx+[V(x)+i W(x)]u-|u|^2u=0.


In what follows, we fix A_1=0.1 and A_2=1, and consider stationary solutions of frequency μ=1. We will make use of periodic boundary conditions.

Notice that the potentials of <cit.>
Ṽ(x)=-[g^2(x)+ig'(x)] would, in the present notation,
necessitate:

    [V'(x)]^2=-4V(x)W^2(x).

It is important to note that the potentials studied in our chapter do not
fulfill this relation
for any set of parameters (A_1, k_1, x_d, δ_1) or (A_2, k_2, δ_2) —except for the “trivial” 𝒫𝒯-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 𝒫𝒯-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 <cit.>, and established a tolerance of ||F[u]||<10^-3 with ||F[u]|| being the L^2-norm of F[u]:


    ||F[u]||=√(∫ |F[u(x)]|^2dx).


In the particular case of potential Ṽ_1(x), we have studied the stability of solitons in the 𝒫𝒯-symmetric limit x_d=δ_1=0 as a function of k_1, observing that solitons are stable whenever k_1<k_c, with k_c=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 k_1 far enough from k_c.
Moreover, since the minimal value attained for
||F[u]|| increases with k_1, we have restricted
consideration to relatively small values of k_1
and more specifically will report results in what follows for    k_1=1/2.

Figures <ref> and <ref> show the potential profile for two different (x_d,δ_1) and (k_2,δ_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 the breaking of both the 𝒫𝒯-symmetry
and the violation of conditions such as the one in Eq. (<ref>),
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. (<ref>), as we will see below.









 §.§ Dynamics


We now analyze the dynamics of several case examples for the NLS equation with potential Ṽ_1(x), using as initial condition the optimized beam profiles found by the Levenberg-Marquardt algorithm. Figures <ref> and <ref> show the outcome
of the simulations for x_d=1 and x_d=-1, respectively, when δ_1=0 is fixed;
on the other hand, Figs. <ref> and <ref> correspond, respectively,
to δ_1=-0.1 and δ_1=0.1, when x_d=1 is fixed.
In these figures, we show the density |ψ(x)|^2 at 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) L^2-norm (power/mass in optics/atomic physics), N(t) (bottom right),
defined as


    N(t)=∫ |ψ(x,t)|^2 dx.






One can observe a clear correlation between the qualitative shape of Im{F[u]} and the growing/decaying character of the dynamics. In other words, in the growing case, this quantity is predominantly positive, whereas for the decaying case, it is predominantly negative.









Moreover, it seems that a larger growth rate (i.e., a faster
increase or decrease of N) is associated to a larger ||F[u]||. In order to
showcase this fact, we have depicted in Fig. <ref> the dependence of
diagnostic quantities λ and σ, that we have
accordingly defined as


    λ=dN/dt|_t=0.



    σ=S||F[u]||,


with


    S=sgn{∫ Im{F[u(x)]}dx}



The quantity σ  takes into account both the (minimized) norm of ||F[u]|| and the form of Im{F[u(x)]} through S –that is, if the imaginary part of F[u(x)] is chiefly positive or negative.




Notice that the blank regions
    correspond to solutions for which ||F[u]|| is higher than the prescribed tolerance of 10^-3. On the other hand, λ characterizes the
    rate of “departure” from the optimized beam profile obtained from
    this minimization procedure.

Figure <ref> shows a clear correlation between σ and λ.
Notice the symmetry between the outcomes when the transformation
(x_d,δ_1)→(-x_d,-δ_1) is applied, which is also manifested
in the values of λ and σ displayed in the captions of
Figs. <ref> and <ref>. From this
figure it is also clear that, roughly speaking,
when x_d δ_1<0, N(t) grows with time, whereas the opposite takes place
when x_d δ_1>0. This is not always true, as there is a critical value
δ_1c (close to zero) separating the growing (λ>0) and decaying (λ<0) dynamics, which is tantamount to the separation of the regions
with σ>0 and σ<0. The dependence of δ_1c versus x_d
is also depicted in Fig. <ref>; having in mind the continuous dependence of σ and λ with x_d and δ_1, it is clear that σ=0 just at the curve δ_1c(x_d)=0, so one can find stationary soliton solutions.
This is manifested in Fig. <ref>, where, for a set of parameters
very close to the curve δ_1c(x_d)=0 (in particular, x_d=1 and δ_1=0.014038),
the decay is very slow (with λ≲10^-7), but not identically zero,
as ||F[u]||∼ 10^-8. Interestingly, as shown in the bottom left panel of the figure,
the relation (<ref>) is not fulfilled. Consequently, there is a range of
parameter values for which states with a very small value of ||F[u]||
can be obtained even if the
potential is not of the form -(g^2+ig').







In the case of the NLS equation with potential Ṽ_2(x), we only focus
on the dependence of λ and σ with respect to parameters k_2 and δ_2,
as the outcome of simulations is essentially the same as in the previous case.
Namely, for non-vanishing values of λ and σ,
a growth or decay of the solutions is identified for typical values
of δ_2, as shown in Fig. <ref>.
However, this growth or decay is quite slow, as achieved by the
optimization of the beam via the Levenberg–Marquardt algorithm.
Notice there is an anti-symmetry in the outcome when
the
transformation k_2→-k_2
is applied. In addition, both σ and λ are equal to
zero at k_2=0 as at that point the potential is real and the solutions are stationary.
Once again, the nearly parabolic curve in the
(δ_2,k_2) plane where λ=σ=0
enables us to identify parameter values in the vicinity of which
states with particularly small ||F[u]|| appear to exist.



§ SYMMETRY BREAKING IN TWO-DIMENSIONAL POTENTIALS


It is of particular interest to extend the above one-dimensional
considerations
towards the emergence of asymmetric optimized beam families
in the 2D version of Eq. (<ref>) that reads:


    iψ_t=-(ψ_xx+ψ_yy)+[V(x,y)+i W(x,y)]ψ-|ψ|^2ψ.


In this case, stationary solutions, ψ(x,y,t)=e^i μ tu(x,y) with u(x,y)∈ℂ,
will satisfy:


    F[u]≡μ u-(u_xx+u_yy)+[V(x,y)+i W(x,y)]u-|u|^2u=0.


In Ref. <cit.>, it is shown that not only symmetric solitons exist but also symmetry breaking is possible if the potential Ṽ(x,y)=V(x,y)+i W(x,y) is of the form


    Ṽ(x,y)=-[g^2(x)+α g(x)+ig'(x)+h(y)]


with g(x) being a spatially even real function, h(y) being a real function and α a real constant. Notice that this potential is partially-𝒫𝒯-symmetric (denoted also as P𝒫𝒯-symmetric), i.e.,


    Ṽ^*(x,y)=Ṽ(-x,y)


The linear spectrum of this potential can be purely real. In this case, a family of 𝒫𝒯-symmetric solitons can emerge from the edge of the continuous spectrum; two degenerate branches of asymmetric solitons, which do not respect the P𝒫𝒯 symmetry, bifurcate from the symmetric soliton branch through a pitchfork bifurcation.

The symmetry breaking bifurcation can also be observed either if the potential possesses double P𝒫𝒯 symmetry


    Ṽ^*(x,y)=Ṽ(-x,y)    and Ṽ^*(x,y)=Ṽ(x,-y)


or 𝒫𝒯- and one P𝒫𝒯-symmetry simultaneously


    Ṽ^*(x,y)=Ṽ(-x,-y)    and    Ṽ^*(x,y)=Ṽ(-x,y) or Ṽ^*(x,y)=Ṽ(x,-y).


In such cases of double symmetries, there is no need for the potential to have
a special form as in Eq. (<ref>). In addition, the soliton branch that emerges from the spectrum edge possesses both symmetries whereas the bifurcating branch loses one of the symmetries although it retains the other.

A later work <cit.> reports the existence of the same branching behaviour in a 𝒫𝒯-symmetric potential which
also features a partial 𝒫𝒯-symmetry along the x-direction.

More specifically, the potential used in <cit.> is given by


    V_3(x,y)=-[G^2(x,y)+G(x,y)],    W_3(x,y)=k_3∂_x G(x,y)


with


    G(x,y)=A_3e^-y^2(e^-(x-1)^2+e^-(x+1)^2).

Notice that the symmetries mentioned above are applicable as
a result of the even nature of the G(x,y).

To give an associated example
of the resulting symmetry breaking, we use, as in <cit.>, A_3=3 and k_3=1. The resulting profile of the potential is shown in Fig. <ref>. 𝒫𝒯-symmetric solitons are calculated by means of
the Newton–Raphson method and the corresponding branch emerges from μ=5.810; asymmetric solitons (actually, optimized beams) are attained by using the Levenberg-Marquardt algorithm, with a tolerance of ||F[u]||<10^-2. Now, the L^2-norm is defined as


    ||F[u]||=√(∫ |F[u(x)]|^2dxdy).







Fig. <ref> represents P≡ N(t=0) versus μ for the symmetric and asymmetric soliton branches; notice that N(t) is now defined as


    N(t)=∫ |ψ(x,y,t)|^2 dxdy.


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 <cit.>. Notice that all the soliton profiles are symmetric with respect to the y-axis; the symmetric solitons present a couple of humps at (x=± x_1,y=0) whereas the asymmetric solitons only possess a single hump at (x=x_2,y=0). We have only shown solitons with x_2>0 as the solutions with x_2<0 are attained simply by making the transform u(x,y)→ u(-x,y).



We have also computed the diagnostic quantities λ and σ [see (<ref>) and Eqs. (<ref>), with S adapted to 2D domains] for the asymmetric soliton and depicted them in Fig. <ref>. Again, we have considered asymmetric solitons branches centred at x=x_2>0. In that case, the norm grows with time, as corresponds to σ>0 and λ>0 whereas the opposite takes place if x_2<0. We can observe, as in the 1D case, a clear correlation between both quantities.

Finally, we show in Fig. <ref> and <ref> the dynamics of the asymmetric and 𝒫𝒯-symmetric solitons with μ=7. As it was pointed out in <cit.>, the 𝒫𝒯-symmetric solitons are unstable past the “bifurcation” point, i.e. when they coexist with the asymmetric branch; as we have shown in Fig. <ref>, 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. <ref>, 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 𝒫𝒯-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 𝒫𝒯-symmetric branch, past the corresponding destabilization
point; compare the associated dynamics of Fig. <ref> with those
of Fig. <ref>.








§ CONCLUSIONS & FUTURE WORK


In the present work, we have revisited a variant of 𝒫𝒯-symmetric
systems. In particular, we have examined multi-parametric potentials
whose parameters control, on the one hand, the potential departure
from the 𝒫𝒯-symmetric case (such as x_d herein), and on
the other hand, the potential degeneracy of the conditions
(<ref>) for stationary solutions –motivated by the recent
works of <cit.>. We have confirmed the
results of the important recent contribution of <cit.>,
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 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. 𝒫𝒯-symmetry
is a sufficient but not a necessary condition for such existence
and extending beyond it seems of particular interest.
The conjecture of <cit.> is that the potentials
V+i W=- (g^2 + i g'(x)) + c represent 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 <cit.>,
but clearly further efforts are necessary to provide
a definitive reply in this direction.










§ 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 <cit.>. 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. <cit.>.

Prior to applying LMA, we need to discretize our equation (<ref>). Thus, we take a grid x_n=-L/2+nh with n=0,1,2… M and L being the domain length, and denote u_n≡ u(x_n) and F_n≡ F[u(x_n)]. With this definition u_xx can be cast as (u_n+1+u_n-1-2u_n)/h^2. In order to simplify the notation in what follows, let us call 𝐮≡{u_n}_n=1^M and 𝐅(𝐮)≡{F_n}_n=1^M. We will also need to define the Jacobian matrix 𝐉(𝐮)≡{J_n,m}_n,m=1^M with J_n,m=∂_u_mF_n. In the presently considered
optimization framework, 𝐅(𝐮) is also knows as the
residue vector.

Let us recall that fixed point methods typically seek a solution by performing the iteration 𝐮_j+1=𝐮_j+δ_j from the seed 𝐮_0 until the residue norm ||𝐅(𝐮)|| is below the prescribed tolerance; here δ_j is dubbed as the search  direction. In the Newton–Raphson method, the search direction is the solution of the equation system 𝐉(𝐮_j)δ_j=-𝐅(𝐮_j). 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 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 <cit.>. 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/2𝐅(𝐮)^T𝐅(𝐮)+
         δ^T𝐉(𝐮)^T𝐅(𝐮)+
         δ^T𝐉(𝐮)^T𝐉(𝐮)δ.


In addition, δ must fulfill the relation


    ||𝐃·δ||<Δ,


where 𝐃 is 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)^T𝐉(𝐮_j)+λ_j𝐃)δ_j=
        -𝐉(𝐮_j)^T𝐅(𝐮_j)


with 𝐃 being the scaling matrix introduced in Eq. (<ref>). There are several possibilities for choosing such matrix. In the present work, we have taken the simplest option, that is 𝐃=𝐈 (the identity matrix), so (<ref>) simplifies to ||δ_j||<Δ. Notice that for λ_j=0, (<ref>) 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 λ_j
constituting the main difficulty of the algorithm.

The scheme of the LMA is described in a quite easy way in Numerical Recipes book <cit.> and is summarized below:




  * 

Take a seed 𝐮_0 and compute ||F(𝐮_0)||


  * 

Choose a value for λ_0. In our particular problem, we have taken λ_0=0.1.


  * 

Solve the equation system (<ref>) in order to get δ_0 and compute ||𝐅(𝐮_0+ δ_0)||


  * 




  * 

If ||F(𝐮_0+δ_0)||≥||𝐅(𝐮_0)||, then take λ_1=10λ_0 and 𝐮_1=𝐮_0, as with this choice of λ_0 the residue norm has not decreased.


  * 

If ||F(𝐮_0+δ_0)||<||𝐅(𝐮_0)||, then take λ_1=λ_0/10 and 𝐮_1=𝐮_0+δ_0, as with this choice of λ_0 has succeeded in decreasing the residue norm.




  * 

Go back to step 3 doing λ_0=λ_1 and 𝐮_0=𝐮_1



This algorithm is repeated while ||𝐅(𝐮)|| is above the prescribed tolerance.



100

  Bender1 Bender, C. M.; Boettcher, S.  Real Spectra in Non-Hermitian Hamiltonians Having 𝒫𝒯 Symmetry, Phys. Rev. Lett. 80, 5243-5246 (1998).

Bender2 Bender, C. M.; Brody, D. C.; Jones, H. F.  Complex Extension of Quantum Mechanics,
    Phys. Rev. Lett. 89, 270401 (2002).

Ruter 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).

Peng2014 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).

peng2014b 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).


ncomms2015 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).

RevPT 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).

Konotop  Konotop, V. V.; Yang, J.; Zezyulin, D. A. Nonlinear waves in 𝒫𝒯-symmetric systems,  
  Rev. Mod. Phys. 88, 035002 (2016).

Schindler1  Schindler, J.;  Li, A.;  Zheng, M. C.;  Ellis, F. M.; Kottos,  T.  Experimental study of active LRC circuits with 𝒫𝒯 symmetries. Phys. Rev. A 84, 040101 (2011).

Schindler2 Schindler, J.;  Lin, Z.;  Lee, J. M.; Ramezani, H.;  Ellis, F. M.; Kottos, T., 𝒫𝒯-symmetric electronics. J. Phys. A: Math. Theor. 45, 444029 (2012).

Factor  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).

Bender3  Bender, C. M.;  Berntson, B.; Parker, D.; Samuel, E.  Observation of 𝒫𝒯 Phase Transition in a Simple Mechanical System. Am. J. Phys.
81, 173-179 (2013).

ptnlde Cuevas-Maraver, J; Kevrekidis, P. G.; Saxena, A.; Cooper, F.;
  Khare, A.; Comech, A.; Bender, C. M.
  Solitary Waves of a 𝒫𝒯-Symmetric Nonlinear Dirac Equation,
  IEEE J. Select. Top. Quant. Electron. 22, 5000109 (2016).

ptkg 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 Klein–Gordon Field Theories,
  Stud. Appl. Math. 133, 298–317 (2014).

abdu  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).

kono Konotop, V.V.; Zezyulin, D. A.
  Families of stationary modes in complex potentials,
  Opt. Lett. 39, 5535–5538 (2014).

jian Yang, J.
  Partially 𝒫𝒯-symmetric optical potentials with
  all-real spectra and soliton families in multi-dimensions",
  Opt. Lett. 39, 1133–1136 (2014).

jenn D'Ambroise, J.; Kevrekidis, P. G.
  Existence, stability and dynamics of nonlinear modes
  in a 2D partially 𝒫𝒯-symmetric potential,
  Appl. Sci. 7, 223 1–10 (2017).


Kom15
Kominis, Y. Dynamic power balance for nonlinear waves in unbalanced gain and loss landscapes.
Phys. Rev. A 92, 063849 (2015)

kom15b Kominis, Y.
  Soliton dynamics in symmetric and non-symmetric complex potentials,
  Opt. Commun. 334, 265–272 (2015).


mja Ablowitz M. J.; Segur, H. Solitons and the Inverse Scattering
Transform (SIAM, Philadelphia, 1981).

ny16 Nixon S. D.; Yang, J. Bifurcation of Soliton Families from Linear Modes in Non-𝒫𝒯-Symmetric Complex Potentials.
Stud. App. Maths. 136 (2016) 459.

dcd16 Dutykh, D.; Clamond, D.; Durán, Á. Efficient computation of capillary-gravity generalised solitary waves. Wave Motion 65 (2016) 1.

aubry Johansson, M.; Aubry, S.
  Growth and decay of discrete nonlinear Schrödinger breathers interacting with internal modes or standing-wave phonons,
  Phys. Rev. E 61, 5864–5879 (2000).

yang2d Yang, J. Symmetry breaking of solitons in two-dimensional complex potentials. Phys. Rev. E 91 (2015) 023201.

ch16 Chen, H.; Hu, S. The asymmetric solitons in two-dimensional parity-time symmetric potentials. Phys. Lett. A 380 (2016) 162.


leven Levenberg, K. A Method for the Solution of Certain Non-Linear Problems in Least Squares. Quarterly of Applied Mathematics 2 (1944) 164.

marq Marquardt, D. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. SIAM Journal on Applied Mathematics 11 (1963) 431

pola1 Cuevas, J.; Rodrigues, A.S., Carretero-González, 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.

pola2 Rodrigues, A.S.; Kevrekidis, P.G.; Cuevas, J.; Carretero-González, 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)

pola3 Rodrigues, A.S.; Kevrekidis, P.G.; Carretero-González, 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.

numrec 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).



