Families of stable solitons and excitations in the PT-symmetric nonlinear Schrödinger equations with position-dependent effective masses

Since the parity-time-(\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pmb{\mathscr{P}}\pmb{\mathscr{T}}$$\end{document}PT-) symmetric quantum mechanics was put forward, fundamental properties of some linear and nonlinear models with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pmb{\mathscr{P}}\pmb{\mathscr{T}}$$\end{document}PT-symmetric potentials have been investigated. However, previous studies of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pmb{\mathscr{P}}\pmb{\mathscr{T}}$$\end{document}PT-symmetric waves were limited to constant diffraction coefficients in the ambient medium. Here we address effects of variable diffraction coefficient on the beam dynamics in nonlinear media with generalized \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pmb{\mathscr{P}}\pmb{\mathscr{T}}$$\end{document}PT-symmetric Scarf-II potentials. The broken linear \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pmb{\mathscr{P}}\pmb{\mathscr{T}}$$\end{document}PT symmetry phase may enjoy a restoration with the growing diffraction parameter. Continuous families of one- and two-dimensional solitons are found to be stable. Particularly, some stable solitons are analytically found. The existence range and propagation dynamics of the solitons are identified. Transformation of the solitons by means of adiabatically varying parameters, and collisions between solitons are studied too. We also explore the evolution of constant-intensity waves in a model combining the variable diffraction coefficient and complex potentials with globally balanced gain and loss, which are more general than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pmb{\mathscr{P}}\pmb{\mathscr{T}}$$\end{document}PT-symmetric ones, but feature similar properties. Our results may suggest new experiments for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pmb{\mathscr{P}}\pmb{\mathscr{T}}$$\end{document}PT-symmetric nonlinear waves in nonlinear nonuniform optical media.

The Hamiltonians in the quantum mechanics are usually required to be Hermitian, which secures the corresponding spectra to be real 1 . Nevertheless, it had been demonstrated by Bender and Boettcher in 1998 that non-Hermitian Hamiltonians obeying the parity-time (PT ) symmetry may also produce entirely real spectra [2][3][4][5][6][7][8] . The PT symmetry implies that the real and imaginary parts of the complex-valued potential, = + U x V x iW x ( ) ( ) ( ), are, respectively, even and odd functions of the coordinate 2 : . For a given real part of the potential, the spectrum of most PT -symmetric systems remains real, as long as the amplitude of the imaginary component of the potential is kept below a certain critical value (the PT -symmetry-breaking threshold); nevertheless, dynamical models featuring unbreakable PT symmetry are known too 9 . Pioneering theoretical works had predicted a possibility to realize the PT -symmetric wave propagation in optical media with symmetrically placed gain and loss elements [10][11][12][13][14] , which was followed by the experimental implementation in optical and atomic settings, including synthetic photonic lattices, metamaterials, microring lasers, whispering-gallery microcavities, and optically induced atomic lattices [15][16][17][18][19][20][21][22] . In particular, phase transitions between regions of the unbroken and broken PT symmetry have been observed in many experiments.
On the theoretical side, the consideration of PT -symmetric potentials in both one-and multi-dimensional linear and nonlinear Schrödinger (NLS) or Gross-Pitaevskii (GP) equations has revealed many remarkable PT -symmetry-breaking phenomena, including several models that give rise to PT -symmetric solitons  . It is commonly known that the soliton theory has been widely applied to fluid mechanics, plasma physics, Bose-Einstein condensates (BECs), nonlinear optics, and many other fields. In particular, optical solitons can utilize the nonlinearity in optical fibers to balance the group-velocity dispersion, thus stably propagating in long-scale telecommunication links. More recently, stable PT -symmetric solitons were also investigated in the third-order NLS equation 53 , the generalized GP equation with a variable group-velocity coefficient 54 , and the derivative NLS equation 55 . The vast work performed in the field of nonlinear waves in PT -symmetric systems has been summarized in two recent comprehensive reviews 56,57 .
As mentioned above 2 , the usual one-dimensional (1D) PT -symmetric Hamiltonian is , with V(−x) = V(x) and W(−x) = −W(x). However, physics of semiconductors gives rise to a Hamiltonian in which the effective mass of collective excitations 58 , M(x), may be variable (position-dependent) [59][60][61][62][63][64][65][66] : . To the best of our knowledge, Hamiltonians combining a variable effective mass and complex-valued PT -symmetric potentials, in particular, , have not been studied yet. In fact, this model represents more general settings than the above-mentioned one occurring in semiconductors. Indeed, the effective mass for collective excitations propagating in lattice media is determined by local properties of the underlying lattice 67,68 , which may be nonuniform in many situations, thus making the effective mass position-dependent. In the context of optics, essentially the same Hamiltonian governs the light propagation in the spatial domain, where the effective diffraction coefficient, 1/M, can be also made x-dependent in nonuniform photonic lattices 69,70 .
The goal of the present work is to introduce 1D and 2D NLS models with such Hamiltonians, incorporating a particular physically relevant PT -symmetric potential, namely, a generalized Scarf-II potentials (i.e., a hyperbolic version of the quantum-mechanical potential introduced by Scarf 71 ). We reveal characteristic properties of both 1D and 2D linear and nonlinear modes in such models, including solitons in the case of cubic nonlinearity, which are quite generic and can be extended to other physically relevant complex-valued PT -symmetric potentials.

Results
1D PT -symmetric nonlinear waves in the effective diffraction. In Kerr-type nonlinear media with the complex-valued PT -symmetric potential and the effective diffraction coefficient defined by the position-dependent mass, ≡ m x Mx ( ) 1/ ( ) (in particular, it represents the above-mentioned variable diffraction coefficient in optics), the scaled 1D modified NLS equation for the wave function ψ x z ( , ) is 2 where z is the propagation distance, x is the transverse coordinate, and g is the strength of the Kerr nonlinearity. Replacing z by time t, Eq. (1) may be used as the GP equation for the BEC wave function, with the effective mass affected by a nonuniform optical lattice 72 . In the spatial-domain optics, real potential V(x) represents the local modification of the refractive index, whereas W(x) stands for the transverse symmetric distribution of the optical gain (W > 0) and loss (W < 0). Under the conditions that V(x) + iW(x) is a PT -symmetric potential and m(x) is an even function of x, it is easy to see that Eq. (1) is invariant under the action of PT -symmetric transformation, where the spatial-reflection operator  and time-reversal operator  are defined as usual 2 ; : , . Equation (1) may be rewritten in the variational form, , with a non-Hermitian but PT -symmetric field Hamiltonian,  , where the asterisk denotes the complex conjugate. Further, the power (norm) of the wave function, Linear spectra of unbroken and broken PT -symmetric phases. We now introduce the diffraction coefficient with a localized spatial modulation, 0 and a physically relevant PT -symmetric ingredient of the model in the form of a generalized Scarf-II potential where m 0 > −1 and α > 0 are both real constants,  θ α α = − + + + (4 6 8) Here m 0 and θ 0 account for the modulation of the local diffraction coefficient m α , even real potential V α , and odd gain-loss distribution W α , respectively. For m 0 = 0, Eq. (2) yields a constant diffraction coefficient, m = 1, hence Eq. (1) reduces to the usual NLS equation. The corresponding complex-valued potential (3) and (4) becomes the usual Scarf-II potential: PT -symmetry breaking when the strength of the imaginary part, 1  , exceeds a threshold  − + 1/4 1 75 . Below, we use the complex potential with m 0 > 0 and α > 0, which may be viewed as a generalization of the basic Scarf-II potential.
First, we consider linear spectra of phases with unbroken and broken PT symmetries in the framework of the linear spectral problem with PT -symmetric operator L that includes variable diffraction coefficient (2) and PT -symmetric potential (3) and (4). The problem is based on the equation for localized eigenfunctions Φ x ( ) and respective eigenvalues λ, In the limit of m 0 = 0 in Eq. (2), L amounts to the standard Hamiltonian operator with the usual PT -symmetric Scarf-II potential, which has been studied in detail by means of analytical and numerical methods 54,75 . For m 0 > 0, we focus on natural values of powers in Eq. (2), α = 1, 2, and 3. By means of the available numerical Fourier spectral algorithm 76, 77 , we find symmetry-breaking boundaries in the (m 0 , θ 0 ) parameter plane, below and above which the PT symmetry is unbroken and broken, respectively, as shown in Figs. 1(a1,b1,c1)). It is seen that, for fixed m 0 , there always exists a critical values of θ 0 , beyond which the symmetry-breaking phase transition makes the spectra complex-valued. On the other hand, for given θ 0 , the phase transition may occur more than once with the increase of m 0 , featuring transitions of the energy spectra first from real to complex, then back to real (restoration of the PT symmetry), and finally again from real to complex values. This scenario is drastically different from what is known in the case of the usual Hamiltonian with the PT -symmetric Scarf-II potential 75 . To illustrate these findings, a few lowest energy levels are displayed in Figs Nonlinear localized modes and their instability. For the given spatial profile of the diffraction coefficient (2) and the generalized PT -symmetric Scarf-II potential (3) and (4), it is possible to find analytically particular exact solutions for stationary nonlinear localized modes (bright solitons) of Eq. (1) for g > 0 in Eq. (1), i.e., the self-focusing sign of the cubic nonlinearity (see Methods): . Without loss of generality, we display subsequent results for the normalization defined by g = 1.
The integral power of the nonlinear localized modes (6) is , which is πρ , whose sign is solely determined by θ 0 for any α > 0. It is clearly seen from Eq. (4) that the signs of gain-loss distribution W α are also determined by the single parameter θ 0 , for m 0 ≥ 0 and α > 0. Thus we conclude that the power always flows from the gain region to one carrying the loss, regardless of the sign of θ 0 .
For different powers α = 1, 2, 3, we aim to study the linear stability of exact nonlinear modes (6) by numerically calculating the largest absolute value of imaginary parts of the linearization eigenvalue δ from Eq. (23), see below, in the modulation-parameter plane (m 0 , θ 0 ). Figures 2(A-C) exhibit the so found stability (blue) and instability (other colors) regions of the localized modes (6) for α = 1, 2, 3, respectively. For these three cases, we, respectively, choose three stable points (i.e., ones with the unbroken PT symmetry): θ = . . . . . . m ( , ) (0 5, 0 1), (0 4, 0 1), (0 2, 0 1) 0 0 , to display the corresponding PT -symmetric potentials in Figs. 2(a-c). The single-well potential becomes deeper, whereas the amplitude of the gain-and-loss distribution slowly decreases, as α increases and m 0 decreases. For these three cases, we show exact solitonic solutions (6) and their numerically found counterparts in Figs , respectively, with a trend to growth. In the existence range of the numerically found solitons, the dependence of the corresponding linear-stability eigenvalues (the largest absolute value of the imaginary part of the perturbation growth rate δ) on the propagation constant −μ is shown in Figs. 2(a3,b3,c3). It is seen that the numerically found soliton solutions in Figs. 2(a1,b1,c1), corresponding to μ = 2.25, 4, 6.25, respectively, are completely stable, in accordance with the stability of the corresponding exact nonlinear modes (see Figs. 2(A-C)).
To validate the linear-stability results, we have simulated the propagation, by taking the input provided by the stationary modes in Figs. 2(a1,b1,c1) with the addition of 2% random perturbations, as is shown in Figs. 2(a4,b4,c4), respectively. In practice, we simulate the beam propagation with the initial input where φ x ( ) is a nonlinear mode, and  is a complex broadband random perturbation. In MATLAB, the 2% white noise  can be realized by utilizing a random matrix such as = .
, where rand(N, 1) is an N-by-1 array of pseudorandom uniform values on the open interval (0, 1) (similarly in other cases, even for the 2D situation). Finally, based on the linear stability results presented in Figs. 2(a3,b3,c3), we choose solutions that are predicted to be unstable (with μ = 10, 5, 10.5), to test the dynamical behavior of the corresponding numerically found soliton solutions. It is found that they are indeed (weakly) unstable, see Figs. 2(a5,b5,c5).
Further, following the results of the linear-stability analysis shown in Figs. 2(A-C), we have tested the propagation dynamics of the exact nonlinear modes (6), varying the parameter m 0 or θ 0 in Eqs. (2)-(4). For the brevity's sake, hereafter we use the words "unbroken" and "broken" to mention if the corresponding linear modes, considered above, do or do not keep their PT symmetry. For fixed α = 1 and m 0 = 0.5, as θ 0 increases continuously from 0.1 to 0.3 (unbroken), nonlinear modes (6) always feature stable propagation dynamics, as shown in Fig. 3(a1). However, as θ 0 increases to slightly larger values, e.g., 0.32 (unbroken), instability sets in, see Fig. 3(a2). Moreover, we have found that a stable nonlinear localized mode (6) with (m 0 , θ 0 ) = (0.75, 0.1) belongs to the region of broken linear PT -symmetry (see Fig. 3(a3)), with the real part of the corresponding complex potential having a single-well shape, see Fig. 3(a4). The latter result implies that the exact PT -symmetric nonlinear modes may be stable while their linear counterparts are not.  (5), with the position-dependent diffraction coefficient (2) and generalized PT -symmetric Scarf-II potentials (3) and (4), is unbroken/broken in the domains below/above the solid blue curves. The results are displayed for α = 1, 2, 3, respectively. (a2,b2,c2) Real and (a3,b3,c3) imaginary parts of first several eigenvalues λ as a function of m 0 for θ 0 = 0.6, 1, 0.7, respectively. The insets in (b2) and (b3,c3) show, severally, the separation and coalescence of the real and imaginary parts of the second and third energy eigenvalues.

Interactions between exact bright solitons and sech pulses.
(broken), with the initial solitary pulse + x e sech( 40) ix 4 . Other parameters are α = 1, 2, 3 for the first, second, and third row, respectively. Note that t-axes should be z-axes. try, and took the corresponding initial conditions as given by Eq. (6) for α = 1, as well as given by Eq. (6) for α = 2, 3. In the case of α = 1, the exact soliton still demonstrates the steady propagation after the collision, while in the cases of α = 2, 3 the amplitude of the initial exact soliton rapidly grows after the collision, thus manifesting interaction-induced instability.
Adiabatic transformation of stable nonlinear modes. Here, we elaborate three different scenario of dynamical control of nonlinear localized modes, making use of adiabatically varying parameters of the potential 38,[53][54][55] , m 0 → m 0 (z) or θ 0 → θ 0 (z). We choose the following temporal-modulation pattern: This implies replacing Eq. (1) by , i.e., m 0 is fixed, while θ 0 varies. In these cases, both the initial and final parameter values correspond to linear modes with unbroken PT symmetry. Likewise, fixing θ 0 and setting m 0 → m 0 (z) as per by Eq. (7), we can perform similar transformations of stable exact nonlinear localized modes, as is shown in Figs. 5(a2,b2,c2). Figure 5(b2) displays a typical stable transformation of the nonlinear modes (6) with α = 2, also ending with parameters corresponding to the linear state with unbroken PT symmetry.
On the other hand, Fig. 5(a2) reveals unstable transformation of the nonlinear modes (6) with α = 1, and Fig. 5(c2) exhibits stable transformation of the nonlinear modes (6) with α = 3, both from initial parameters corresponding to the linear state with unbroken PT symmetry to final ones corresponding to broken PT symmetry in the linear state. Finally, the transformation is also implemented by varying both m 0 and θ 0 . Instability of such simultaneous transformation, observed in Fig. 5(a3) may be attributed to the instability of the corresponding single-parameter variation, shown above in Fig. 5(a2). Further, Fig. 5(b3) reveals that the simultaneous variation of m 0 and θ 0 may give rise to an unstable output, while both corresponding the single-component variations are stable, see Figs. 5(b1,b2). Nevertheless, in Fig. 5(c3) we stably transform an initially exact nonlinear mode, with unbroken PT symmetry of the respective linear state, to another exact nonlinear mode, which corresponds to the broken PT symmetry in the linear state, varying both parameters. The stable outcome of the simultaneous variation of m 0 and θ 0 may be possible if the separate variation of each parameter produced a stable output.
Thus, we conclude that the simultaneous variation of the two parameters leads to an unstable output if either of the two corresponding separate variations is unstable, see Figs. 5(a1-a3). If both separate excitations produce stable outputs, their simultaneous action may give either an unstable output, see Figs. 5(b1-b3), or a stable one, see Figs. 5(c1-c3).

1D constant-intensity waves.
Recently, a class of complex potentials that are more general than the PT -symmetric complex-valued ones has been put forward 39,[78][79][80] . Similar to the PT -symmetric potentials, these more general complex potentials admit the existence of continuous families of stationary states, supported by the balance between gain and loss (unlike isolated stationary solutions found in generic dissipative systems), a part of which may be stable 39,[78][79][80] . The real and imaginary parts of these potentials are defined 81 is an arbitrary real function. Here, we address a generalization of such potentials in the model with the variable diffraction coefficient, m(x), in the form of 2 where v(x) is a known real function of space. If m(x) (e.g., that defined in Eq. (2)) and v(x) are both even functions, then the complex potential given by Eq. (9) is a PT -symmetric one. In a more general case, when the potential is not PT -symmetric, it nevertheless provides for the global balance between the gain and loss, in the case of localized or periodic functions m(x)v(x), because W(x) is defined in Eq. (9) as a full derivative.
For the general potential taken as per Eq. (9), stationary constant-intensity (alias CW, i.e., continuous-wave) solutions of Eq. (1) with any g are found in the form of (see Methods) , without the correct phase given by Eq. (10), the beams grow fast in the center, leading to instability, as shown in Figs. 6(b), 7(b) and 8(b). However, if the inputs are taken as solution (10) with the correct phase, the growth of perturbations is initially suppressed, occurring later as the modulational instability of the CW (Figs. 6(c), 7(c) and 8(c)). It is worthy to note that the beam in the PT -symmetric system with the HG complex potential steadily propagates farther than that in the case of the non-PT HG complex potential, as the PT -symmetric potential naturally provides for a better balance of the gain and loss. When the truncation length of the input CW becomes larger, the stable-propagation distance . Cases (a1,b1,b2,c1) exhibit the stable transformation between parameter sets correspond to stable linear states with unbroken PT symmetry. (c2,c3) Stable transformation between parameter sets corresponding to unbroken and broken PT -symmetry of the linear states. (a2,a3,b3) Unstable transformation of the initially stable nonlinear mode between parameter sets corresponding to unbroken and broken PT -symmetry of the linear states. Other parameters are α = 1, 2, 3 for the first, second, and third rows, respectively. Note that t-axes should be z-axes.  increases too under the action of the PT -symmetric HG potential (see Fig. 6(d)), but it does not increase in the case of the non-PT -symmetric HG potential (see Fig. 7(d)).
We have also investigated the nonlinear evolution of the CW input in the presence of the self-focusing and self-defocusing Kerr nonlinearity (g = +1 and −1, respectively) for the same three complex potentials. As a result, we find that the CW state maintains stable propagation over a smaller distance than that in the corresponding linear model, as can be seen in Figs. 6(e), 7(e) and 8(e) for the self-focusing nonlinearity, and in Figs. 6(f), 7(f) and 8(f) for the self-defocusing nonlinearity.
2D PT -symmetric nonlinear waves. Multidimensional spatial solitons are a subject of great interest to nonlinear optics [82][83][84][85] . We here consider the formation of bright spatially localized solitons in the 2D PT -symmetric setting. In this case, the field evolution is governed by the 2D NLS equation with a PT -symmetric potential: 2 where ∇ is the 2D gradient ∂ ∂ ( , ) x y , and m(x, y) is a 2 × 2 matrix function, which we take in the diagonal form, , where u(x, y) and θ(x, y) are real amplitude and phase, respectively, which satisfy stationary equations: x x y y x y 2 2 3 The 2D PT -symmetric potentials + V x y iW x y ( , ) ( , ) (i.e., V(x, y) = V(−x, −y) and W(−x, −y) = −W(x, y)), which, like in the 1D setting, admit particular exact solutions for α = 1, 2, and 3, are given below in Eqs. (24a), (25a) and (26a), respectively, along with the exact solutions (see Methods).

Comparison of exact 2D solitons and numerical solutions.
To verify the analytical soliton solutions, we have numerically found the corresponding stationary fundamental solitons of Eq. (11), which are displayed in the third row of Figs. 9, 10 and 11, respectively. As is shown in Figs. 9(c), 10(c) and 11(c), the difference between the exact solutions and their numerical counterparts falls below 10 −9 , hence the exact solutions are correct. . Furthermore, using the numerical method, we calculate the integral power at other values of the soliton parameter μ and identify the existence ranges of the numerically found solitons, as shown in Fig. 12(a), where the lower cutoffs are µ ≈ .
, and there are no finite upper cutoffs. These numerical soliton solutions are more often than not unstable, especially for larger values of μ. We display their linear-stability spectra in Figs. 12(b-d) around μ = 9/2, 8, and 25/2, respectively, i.e., around the values at which the corresponding exact nonlinear modes exist. As a result, we find that in first case, α θ = .
. Dynamical behavior of 2D nonlinear modes. To display the evolution of the 2D exact or numerically found nonlinear modes, we plot the corresponding intensity isosurfaces. For α = 1, as shown in Figs. 13(a,d), instability is produced by simulations of the long-time evolution, although at shorter times, such as t = 200, the solution seems as a stable one, see Fig. 13(b). The instability sets in at ≈ t 280 (Fig. 13(c)). However, for α = 2 and 3, both Figs. 13(e,g) exhibit fully robust evolution of the initial-state solitons taken from Figs. 10 and 11, respectively. To confirm their stability, we display the corresponding final-state soliton profiles in Figs. 13(f,h), which are identical to the corresponding initial profiles in Figs. 10 and 11. We have also numerically investigated the stability of the solitons at the lowest-power points μ 1 = 3.4, μ 2 = 6.3, and μ 3 = 11.1 for α = 1, 2, and 3, respectively. They all turn out to be unstable, although they may seem stable at a shorter propagation distance. Figure 14  for α = 1, 2, 3. For this reason, the exact soliton solutions are especially important ones, as they help to spot stability areas for broad soliton families.
It is worthy to note that the dynamical-stability results are in good agreement with the predictions produced above by the linear-stability analysis (see Figs. 12(b-d)). Thus, the latter analysis, in the combination with systematic direct simulations, make it possible to identify soliton-stability regions in a reliable form.

Discussion
We have reported analytical and numerical results for new classes of 1D and 2D stable spatial solitons in cubic nonlinear media with the PT -symmetric generalized Scarf-II potentials and variable (position dependent) diffraction coefficients. First, in the linear version of the model, parameter regions of the unbroken and broken PT symmetry have been numerically delineated. Then, in the presence of the Kerr nonlinearity, particular exact solutions for nonlinear localized modes with real eigenvalues have been obtained in the analytical form, and verified numerically. These solitons are shown to be stable through the linear-stability analysis and by means of Figure 11. 2D PT -symmetric nonlinear waves for α = 3. The same as in Fig. 9, but for α = 3, and μ = 25/2 in panel (c). Parameters are θ = .
. m ( , ) (0 2, 0 1) 0 0 . direct simulations, in wide ranges of the governing parameters. It is worthy to note that the addition of the Kerr nonlinearity can fix the broken PT symmetry of the linear system, transforming complex eigenvalues into real ones. In addition, stable bright solitons have been found in parameter regions where the PT symmetry of the linear states is broken, for various shapes of the underlying real part of the potential, such as single-and double-well forms. Finally, interactions and adiabatic transformations of the exact solitons have been studied in detail, and the existence range and propagation dynamics of numerically found solitons have been examined too. We also study the evolution of constant-intensity waves in a model combining the variable diffraction coefficient and complex potentials with globally balanced gain and loss, which are more general than PT -symmetric ones, but feature similar properties. These theoretical results suggest new experiments for PT -symmetric nonlinear waves in nonlinear and nonuniform optical media, and provide useful theoretical guidance for studies in related fields, such as BECs.

Methods
Nonlinear stationary modes. Stationary solution of Eq. (1) are looked for in the usual form, where −μ is the propagation constant, and the localized complex wave function satisfies the following ordinary differential equation with the x-dependent diffraction coefficient: 2 We take the complex wave function as with real amplitude u(x), and superfluid phase velocity given by The amplitude satisfies the following second-order ordinary differential equation: which may be transformed by setting ≡  u x m x u ( ) ( ) x into a system of coupled first-order equations: solvable as a boundary-value problem by means of standard shooting methods 60 . To achieve a higher precision and computation speed, we actually used the spectral renormalization method 86 with some necessary modifications. The method is spectrally efficient and relatively easy to implement not only in the 1D case but also in the higher-dimensional settings.

Implementation of the numerical solution.
To construct 1D localized solutions, one first needs to develop a convergent iteration, to guarantee that the amplitude neither blows up nor decays to zero. This may be realized by setting φ λ = x wx ( ) ( ), where λ is a constant to be determined. Using the Fourier transform and the modified spectral renormalization method 86 , we thus arrive as the following iteration scheme: ,  denotes the 1D Fourier transform, and p is an appropriate positive constant (here p = 10 is taken for the 1D case). A Gaussian or sech can be taken as an input, which eventually leads to absolute errors <10 −9 , in both the convergence criterion and the numerically obtained solution satisfying Eq. (17). In general, one may restrict the number of iterations to <100; however, for the sake of high precision, we admitted up to 1000 iterative steps. Once the above-mentioned conditions for two absolute errors are satisfied simultaneously, the desired numerical soliton is obtained as In the 2D case, we needed to accordingly change F 1 in Eq. (20), although the iteration scheme ran similar to its counterpart in the 1D case:  Fig. 10, and (f) its final shape. (g) The stable evolution of the soliton from Fig. 11, and (h) its final shape. Note that t-axes should be z-axes in (a,e) and "t=" should be "z=" in (b,c,d,f,g,h).
Scientific RepoRts | 7: 1257 | DOI:10.1038/s41598-017-01401-3 ) },  denotes, the 2D Fourier transformation, and p is an appropriate positive constant (here p = 100 is taken for the 2D case). Other settings and procedures are similar to those in the 1D case.
The linear-stability analysis. For the given position-dependent function m(x) and complex-valued PT -symmetric potential V(x) + iW(x), one may solve Eq. (17) (or equivalently, Eqs. (18) and (19)), to obtain stationary soliton solutions φ x ( ), by analytical or the above-mentioned methods. Then localized nonlinear modes of Eq. (1) can be found in the stationary form, as ψ φ = µ x z x e ( , ) ( ) i z . To explore the linear stability of the localized modes in the 1D case, we consider a perturbed solution 87,88 , where  is an infinitesimal perturbation amplitude, F(x) and G(x) are the eigenfunctions of the linearized problem, and −δ is the respective eigenfrequency, the instability taking place if some eigenvalues are not purely real.
. The PT -symmetric nonlinear modes are linearly stable provided that δ has no imaginary part, otherwise they are linearly unstable. The whole stability spectrum δ can be numerically calculated by the Fourier collocation method (see ref. 88 . Other technical details are similar to those in the 1D case. As concerns the numerical computation of the full stability spectrum in the 2D case, the Fourier-collocation method is usually of low precision for a small number of Fourier modes. If one increases the number of the modes for a higher accuracy, the necessary size of the dense matrix corresponding to the eigenvalue problem may become prohibitively large 88 . Therefore the full 2D linear-stability spectrum in the (m 0 , θ 0 ) space is not displayed any more, because of the necessary large space size and number of Fourier modes. But we will roughly depict the 2D linear-stability spectrum with respect to the soliton parameter μ (see Figs. 12(b-d)), in contrast to the dynamical stability of long-time wave propagation. By the way, on the account of the same reason that large spatial domains and number of the Fourier modes are necessary for a high accuracy in the 2D case, the corresponding PT -symmetric linear spectra are not exhibited further. Nevertheless, through repeated numerical tests, we find that it is instructive that the PT -symmetric breaking curves in the 1D case can provide powerful reference for those in the 2D case. where μ = 25/2.