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

TL;DR

This paper investigates the existence, stability, and dynamics of solitons in PT-symmetric nonlinear Schrödinger equations with position-dependent diffraction, revealing new stable solutions and their behaviors in variable media.

Contribution

It introduces the study of PT-symmetric nonlinear Schrödinger equations with variable diffraction coefficients, finding stable soliton families and analyzing their properties and interactions.

Findings

Stable families of 1D and 2D solitons identified.

Analytical solutions for some stable solitons derived.

Soliton dynamics and collisions studied.

Abstract

Since the parity-time-(PT-) symmetric quantum mechanics was put forward, fundamental properties of some linear and nonlinear models with PT-symmetric potentials have been investigated. However, previous studies of 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 $\mathcal{PT}$-symmetric Scarf-II potentials. The broken linear 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 PT-symmetric ones, but feature similar properties. Our results may suggest new experiments for PT-symmetric nonlinear waves in nonlinear nonuniform optical media.

Equations

$$i\frac{\partial \psi }{\partial z}=\left[-\frac{\partial }{\partial x}\left(m(x)\frac{\partial }{\partial x}\right)+V(x)+iW(x)-g|\psi {|}^2\right]\psi ,$$

$$m(x)=m_{\alpha }(x)=m_0{\mathrm{sech}}^{\alpha }x+1,$$

$$V_{\alpha }(x)$$

$$W_{\alpha }(x)$$

$$\hat{L}\Phi (x)=\lambda \Phi (x),\hat{L}\equiv -\frac{d}{dx}\left(m(x)\frac{d}{dx}\right)+V(x)+iW(x).$$

$$\psi (x,z)=\rho {\mathrm{sech}}^{\frac{\alpha +2}{2}}x\exp \left[i{\theta }_0{\tan }^{-1}(\sinh x)+i\mu z\right],$$

$$\left\{m_0,{\theta }_0\right\}(z)=\{\begin{array}{ll}\left(\left\{m_{02},{\theta }_{02}\right\}-\left\{m_{01},{\theta }_{01}\right\}\right)\sin (\pi z/2000)+\left\{m_{01},{\theta }_{01}\right\},& \text{0≤z<1000},\\ \left\{m_{02},{\theta }_{02}\right\},& \text{z≥1000.}\end{array}$$

$$i\frac{\partial \psi }{\partial z}=\left[-\frac{\partial }{\partial x}m(x,z)\frac{\partial }{\partial x}+V(x,z)+iW(x,z)-g|\psi {|}^2\right]\psi ,$$

$$V(x)=-m(x)v^2(x),W(x)=\frac{d(m(x)v(x))}{dx},$$

$$\psi (x,z)=C\exp \left[i{\int }_0^{x}v(x)dx+ig{C}^{2}z\right].$$

$$i\frac{\partial \psi }{\partial z}=\left\{-\nabla \left[\mathbf{m}(x,y)\nabla \right]+V(x,y)+iW(x,y)-g|\psi {|}^2\right\}\psi ,$$

$$(m_{\alpha }(x)u_x{)}_x+(m_{\alpha }(y)u_y{)}_y-[{\theta }_x^2{m}_{\alpha }(x)+{\theta }_y^2{m}_{\alpha }(y)]u-V(x,y)u-\mu u+g{u}^3=0,$$

$$(m_{\alpha }(x){\theta }_x{u}^2{)}_x+(m_{\alpha }(y){\theta }_y{u}^2{)}_y-W(x,y)u^2=0.$$

$$\psi (x,y,z)=B{e}^{i\theta (x,y)+ig{A}^{2}z},B=\mathrm{const}$$

$$V(x,y)=-\nabla \theta [(\nabla \theta )\mathbf{m}],W(x,y)=\nabla [(\nabla \theta )\mathbf{m}],$$

$$\left[\frac{d}{dx}m(x)\frac{d}{dx}-V(x)-iW(x)+g|\varphi {|}^2\right]\varphi =\mu \varphi .$$

$$v(x)=\frac{1}{m(x)u^2(x)}{\int }_0^{x}W(s)u^2(s)ds.$$

$$\frac{d}{dx}\left[m(x)\frac{du}{dx}\right]=[V(x)+m(x)v^2(x)-g{u}^2+\mu ]u,$$

$$\frac{du(x)}{dx}=\frac{\tilde{u}(x)}{m(x)},$$

$$\frac{d\tilde{u}(x)}{dx}=[V(x)+m(x)v^2(x)-g{u}^2+\mu ]u,$$

$$\hat{w}=\frac{1}{k_x^2+p}{F}_1+\frac{{\lambda }^2}{k_x^2+p}{F}_2,$$

$$\hat{w}=\frac{1}{k_x^2+k_y^2+p}{F}_1+\frac{{\lambda }^2}{k_x^2+k_y^2+p}{F}_2,$$

$$\psi (x,z)=\{\varphi (x)+ϵ[F(x)e^{i\delta z}+G^{\ast }(x)e^{-i{\delta }^{\ast }z}]\}{e}^{i\mu z},$$

$$\left(\begin{array}[]{cc}\hat{L}_{1}&g\phi^{2}\vspace{0.05in}\\ -g\phi^{\ast 2}&-\hat{L}_{1}^{\ast}\end{array}\right)\left(\begin{array}[]{c}F(x)\vspace{0.05in}\\ G(x)\\ \end{array}\right)=\delta\left(\begin{array}[]{c}F(x)\vspace{0.05in}\\ G(x)\\ \end{array}\right),$$

$${\mathbf{m}}_1(x,y)=\mathrm{diag}(m_0\mathrm{sech}x+1,m_0\mathrm{sech}y+1),$$

$$V_1(x,y)=\frac{1}{4}\sum _{\sigma =x,y}\left[15m_0\mathrm{sech}\sigma -(4{\theta }_0^2+15){\mathrm{sech}}^2\sigma -m_0(4{\theta }_0^2+21){\mathrm{sech}}^3\sigma \right]+m_0({\theta }_0^2+\frac{21}{4})(\mathrm{sech}x\mathrm{sech}y{)}^3,$$

$$W_1(x,y)=-{\theta }_0\sum _{\sigma =x,y}\mathrm{sech}\sigma \tanh \sigma (5m_0\mathrm{sech}\sigma +4),$$

$${\varphi }_1(x,y)=\sqrt{m_0(4{\theta }_0^2+21)/(4g)}(\mathrm{sech}x\mathrm{sech}y{)}^{3/2}\exp \left\{i{\theta }_0[\mathrm{ta}{\mathrm{n}}^{-1}(\sinh x)+\mathrm{ta}{\mathrm{n}}^{-1}(\sinh y)]\right\},$$

$${\mathbf{m}}_2(x,y)=\mathrm{diag}(m_0{\mathrm{sech}}^{2}x+1,m_0{\mathrm{sech}}^{2}y+1),$$

$$V_2(x,y)=\sum _{\sigma =x,y}\left[(8m_0-{\theta }_0^2-6){\mathrm{sech}}^2\sigma -m_0({\theta }_0^2+10){\mathrm{sech}}^4\sigma \right]+m_0({\theta }_0^2+10)(\mathrm{sech}x\mathrm{sech}y{)}^4,$$

$$W_2(x,y)=-{\theta }_0\sum _{\sigma =x,y}\mathrm{sech}\sigma \tanh \sigma (7m_0{\mathrm{sech}}^2\sigma +5),$$

$${\varphi }_2(x,y)=\sqrt{m_0({\theta }_0^2+10)/g}(\mathrm{sech}x\mathrm{sech}y{)}^2\exp \left[i{\theta }_0\sum _{\sigma =x,y}\mathrm{ta}{\mathrm{n}}^{-1}(\sinh \sigma )\right],$$

$${\mathbf{m}}_3(x,y)=\mathrm{diag}(m_0{\mathrm{sech}}^{3}x+1,m_0{\mathrm{sech}}^{3}y+1),$$

$$V_3(x,y)=-\frac{1}{4}\sum _{\sigma =x,y}\left[(4{\theta }_0^2+35){\mathrm{sech}}^2\sigma -55m_0{\mathrm{sech}}^3\sigma +m_0(4{\theta }_0^2+65){\mathrm{sech}}^5\sigma \right]$$

$$+\frac{m_0}{4}(4{\theta }_0^2+65)(\mathrm{sech}x\mathrm{sech}y{)}^5,$$

$$W_3(x,y)=-3{\theta }_0\sum _{\sigma =x,y}\mathrm{sech}\sigma \tanh \sigma (3m_0{\mathrm{sech}}^3\sigma +2),$$

Full text

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

Yong Chen1,2, Zhenya Yan*1,2†* † Correspondence and requests for materials should be addressed to Z.Y. (zyyan@mmrc.iss.ac.cn)**, Dumitru Mihalache3 & Boris A. Malomed4,5

1Key Laboratory of Mathematics Mechanization, Institute of Systems Science, AMSS, Chinese Academy of Sciences, Beijing 100190, China

2School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China

3Department of Theoretical Physics, Horia Hulubei National Institute of Physics and Nuclear Engineering, PO Box MG-6, Bucharest, Romania

4Department of Physical Electronics, School of Electrical Engineering, Tel Aviv University, Tel Aviv 59978, Israel

5Laboratory of Nonlinear-Optical Informatics, ITMO University, St. Petersburg 197101, Russia

Since the parity-time-($\mathcal{PT}$-) symmetric quantum mechanics was put forward, fundamental properties of some linear and nonlinear models with $\mathcal{PT}$-symmetric potentials have been investigated. However, previous studies of $\mathcal{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 $\mathcal{PT}$-symmetric Scarf-II potentials. The broken linear $\mathcal{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 $\mathcal{PT}$-symmetric ones, but feature similar properties. Our results may suggest new experiments for $\mathcal{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 ($\mathcal{PT}$) symmetry may also produce entirely real spectra [2, 3, 4, 5, 6, 7, 8]. The $\mathcal{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: $V(x)=V(-x)$, $W(-x)=-W(x)$ [2]. For a given real part of the potential, the spectrum of most $\mathcal{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 $\mathcal{PT}$-symmetry-breaking threshold); nevertheless, dynamical models featuring unbreakable $\mathcal{PT}$ symmetry are known too [9]. Pioneering theoretical works had predicted a possibility to realize the $\mathcal{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 $\mathcal{PT}$ symmetry have been observed in many experiments.

On the theoretical side, the consideration of $\mathcal{PT}$-symmetric potentials in both one- and multi-dimensional linear and nonlinear Schrödinger (NLS) or Gross-Pitaevskii (GP) equations has revealed many remarkable $\mathcal{PT}$-symmetry-breaking phenomena, including several models that give rise to $\mathcal{PT}$-symmetric solitons [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52]. 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 $\mathcal{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 $\mathcal{PT}$-symmetric systems has been summarized in two recent comprehensive reviews [56, 57].

As mentioned above, the usual one-dimensional (1D) $\mathcal{PT}$-symmetric Hamiltonian is $\mathcal{H}=-\partial_{x}^{2}+V(x)+iW(x)$, with $V(-x)=V(x)$ and $W(-x)=-W(x)$ [2]. 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]: $\mathcal{H}_{M}=-\frac{1}{2}\hbar^{2}\partial_{x}\left(\frac{1}{M(x)}\partial_{x}\right)+V(x)$.

To the best of our knowledge, Hamiltonians combining a variable effective mass and complex-valued $\mathcal{PT}$-symmetric potentials, in particular, $\mathcal{H}_{M}^{\left(\mathcal{PT}\right)}=-\frac{1}{2}\hbar^{2}\partial_{x}\left(\frac{1}{M(x)}\partial_{x}\right)+V(x)+iW(x)$, 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 $\mathcal{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 $\mathcal{PT}$-symmetric potentials.

Results

1D $\mathcal{PT}$-symmetric nonlinear waves in the effective diffraction. In Kerr-type nonlinear media with the complex-valued $\mathcal{PT}$-symmetric potential and the effective diffraction coefficient defined by the position-dependent mass, $m(x)\equiv 1/M(x)$ (in particular, it represents the above-mentioned variable diffraction coefficient in optics), the scaled 1D modified NLS equation for the wave function $\psi(x,z)$ is

[TABLE]

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 $\mathcal{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 $\mathcal{PT}$-symmetric transformation, where the spatial-reflection operator $\mathcal{P}$ and time-reversal operator $\mathcal{T}$ are defined as usual [2], $\mathcal{P}:\,x\rightarrow-x;\,\mathcal{T}:i\rightarrow-i,\,z\rightarrow-z$ . Equation (1) may be rewritten in the variational form, $i\psi_{z}=\delta\mathcal{H}(\psi)/\delta\psi^{\ast}$, with a non-Hermitian but $\mathcal{PT}$-symmetric field Hamiltonian, $\mathcal{H}(\psi)=\int_{-\infty}^{+\infty}\{m(x)|\psi_{x}|^{2}+[V(x)+iW(x)]|\psi|^{2}-\frac{g}{2}|\psi|^{4}\}dx$, where the asterisk denotes the complex conjugate. Further, the power (norm) of the wave function, $P(z)=\int_{-\infty}^{+\infty}|\psi(x,z)|^{2}dx$, evolves according to equation $dP/dz=2\int_{-\infty}^{+\infty}W(x)|\psi|^{2}dx$, which can be immediately deduced from Eq. (1).

Linear spectra of unbroken and broken $\mathcal{PT}$-symmetric phases. We now introduce the diffraction coefficient with a localized spatial modulation,

[TABLE]

and a physically relevant $\mathcal{PT}$-symmetric ingredient of the model in the form of a generalized Scarf-II potential

[TABLE]

where $m_{0}>-1$ and $\alpha>0$ are both real constants,   ${\mathcal{V}_{1}}=-\frac{1}{4}(4\theta_{0}^{2}+\alpha^{2}+6\alpha+8)$,  ${\mathcal{V}_{2}}=\frac{m_{0}}{4}(3\alpha^{2}+8\alpha+4)$,  ${\mathcal{W}_{1}}=-\theta_{0}(\alpha+3)$, and ${\mathcal{W}_{2}}=-m_{0}\theta_{0}(2\alpha+3)$. Here $m_{0}$ and $\theta_{0}$ account for the modulation of the local diffraction coefficient $m_{\alpha}$, even real potential $V_{\alpha}$, and odd gain-loss distribution $W_{\alpha}$, 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: $U(x)={\mathcal{V}_{1}}\mathrm{sech}^{2}x+i{\mathcal{W}_{1}}\mathrm{sech}x\,\mathrm{tanh}x$  [73], where the corresponding Hamiltonian $H(x)=-\partial_{x}^{2}+U(x)$ may have two branches of energy eigenvalues and interpreted as the so-called quasi-parity [74]. Moreover, the Hamiltonian $H(x)$ can be shown to exhibit the spontaneous $\mathcal{PT}$-symmetry breaking when the strength of the imaginary part, $|{\mathcal{W}_{1}}|$, exceeds a threshold $-{\mathcal{V}_{1}}+1/4$ [75]. Below, we use the complex potential with $m_{0}>0$ and $\alpha>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 $\mathcal{PT}$ symmetries in the framework of the linear spectral problem with $\mathcal{PT}$-symmetric operator $\hat{L}$ that includes variable diffraction coefficient (2) and $\mathcal{PT}$-symmetric potential (3)-(4). The problem is based on the equation for localized eigenfunctions $\Phi(x)$ and respective eigenvalues $\lambda$,

[TABLE]

In the limit of $m_{0}=0$ in Eq. (2), $\hat{L}$ amounts to the standard Hamiltonian operator with the usual $\mathcal{PT}$-symmetric Scarf-II potential, which has been studied in detail by means of analytical and numerical methods [75, 54]. For $m_{0}>0$, we focus on natural values of powers in Eq. (2), $\alpha=1,2,$ and $3$. By means of the available numerical Fourier spectral algorithm [76, 77], we find symmetry-breaking boundaries in the $(m_{0},\theta_{0})$ parameter plane, below and above which the $\mathcal{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 $\theta_{0}$, beyond which the symmetry-breaking phase transition makes the spectra complex-valued. On the other hand, for given $\theta_{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 $\mathcal{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 $\mathcal{PT}$-symmetric Scarf-II potential [75]. To illustrate these findings, a few lowest energy levels are displayed in Figs. 1(a2, b2, c2) (real parts) and Figs. 1(a3, b3, c3) (imaginary parts).

Nonlinear localized modes and their instability. For the given spatial profile of the diffraction coefficient (2) and the generalized $\mathcal{PT}$-symmetric Scarf-II potential (3)-(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):

[TABLE]

where $\rho=\sqrt{m_{0}(4\theta_{0}^{2}+3\alpha^{2}+10\alpha+8)/(4g)}$, and the propagation constant is $-\mu=-\left(\frac{\alpha}{2}+1\right)^{2}$. 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 $P=\int_{-\infty}^{+\infty}|\psi(x,z)|^{2}dx=\rho^{2}\int_{-\infty}^{+\infty}\mathrm{sech}^{\alpha+2}xdx$, which is $\frac{1}{2}\pi\rho^{2}$, $\frac{4}{3}\rho^{2}$, $\frac{3}{8}\pi\rho^{2}$ for $\alpha=1,2,3$, respectively. It is also relevant to examine the transverse power flow of these modes (alias the Poynting vector), $S(x)=\left(i/2\right)(\phi\phi_{x}^{\ast}-\phi^{\ast}\phi_{x})=\rho^{2}\theta_{0}\mathrm{sech}^{\alpha+3}(x)$, whose sign is solely determined by $\theta_{0}$ for any $\alpha>0$. It is clearly seen from Eq. (4) that the signs of gain-loss distribution $W_{\alpha}$ are also determined by the single parameter $\theta_{0}$, for $m_{0}\geq 0$ and $\alpha>0$. Thus we conclude that the power always flows from the gain region to one carrying the loss, regardless of the sign of $\theta_{0}$.

For different powers $\alpha=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 $\delta$ from Eq. (23), see below, in the modulation-parameter plane $(m_{0},\theta_{0})$. Figures 2(A,B,C) exhibit the so found stability (blue) and instability (other colors) regions of the localized modes (6) for $\alpha=1,2,3$, respectively. For these three cases, we, respectively, choose three stable points (i.e., ones with the unbroken $\mathcal{PT}$ symmetry): $(m_{0},\theta_{0})=(0.5,0.1),(0.4,0.1),(0.2,0.1)$, to display the corresponding $\mathcal{PT}$-symmetric potentials in Figs. 2(a,b,c). The single-well potential becomes deeper, whereas the amplitude of the gain-and-loss distribution slowly decreases, as $\alpha$ increases and $m_{0}$ decreases. For these three cases, we show exact solitonic solutions (6) and their numerically found counterparts in Figs. 2(a1,b1,c1), to corroborate that they are mutually identical. For other values of propagation constant $-\mu$, exact analytical solutions are not available, but we can use the numerical method (validated by the comparison with the exact solutions) to produce fundamental solitons. It follows from Figs. 2(a2,b2,c2) that the solitons’ powers are monotonously increasing functions of $\mu$. While existence ranges of the numerically found solitons have nearly the same upper cutoff $\mu_{\mathrm{up}}\approx 19.8$, the lower cutoffs are different: $\mu_{1\mathrm{low}}=0.6,~{}\mu_{2\mathrm{low}}=1.48,~{}\mu_{3\mathrm{low}}=4.11$, 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 $\delta$) on the propagation constant $-\mu$ 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 $\mu=2.25,4,6.25$, respectively, are completely stable, in accordance with the stability of the corresponding exact nonlinear modes (see Figs. 2(A,B,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 $\psi(x,z=0)=\phi(x)(1+\epsilon)$, where $\phi(x)$ is a nonlinear mode, and $\epsilon$ is a complex broadband random perturbation. In MATLAB, the $2\%$ white noise $\epsilon$ can be realized by utilizing a random matrix such as $\epsilon=0.04(\mathrm{rand}(N,1)-0.5)(1+i)/\sqrt{2}$, where $\mathrm{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 $\mu=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,B,C), we have tested the propagation dynamics of the exact nonlinear modes (6), varying the parameter $m_{0}$ or $\theta_{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 $\mathcal{PT}$ symmetry. For fixed $\alpha=1$ and $m_{0}=0.5$, as $\theta_{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 $\theta_{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},\theta_{0})=(0.75,0.1)$ belongs to the region of broken linear $\mathcal{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 $\mathcal{PT}$-symmetric nonlinear modes may be stable while their linear counterparts are not.

For fixed $\alpha=2$ and $m_{0}=0.4$, as $\theta_{0}$ grows continuously from $0.1$ to $0.5$ (unbroken) and further to $0.55$ (broken), similar results occur, see Figs. 3(b1,b2). Here too, we find a stable nonlinear localized mode, while the $\mathcal{PT}$ symmetry of its linear counterpart is broken, for parameters $(m_{0},\theta_{0})=(0.8,0.3)$, see Fig. 3(b3), even if the real part of the corresponding complex potential is a weak barrier (rather than a well), whose amplitude is smaller than that of the corresponding gain-loss distribution, see Fig. 3(b4). For $\theta_{0}=0.1$, we increased $m_{0}$ continuously from $0.4$ to $0.6$ (unbroken), so that a continuous family of stable solitons could also be readily found.

For $\alpha=3$, when we fix $m_{0}=0.2$ and increase $\theta_{0}$ continuously from $0.1$ to $0.36$ (unbroken), similar results still hold, see Fig. 3(c1). On the other hand, when, for $\theta_{0}=0.1$, $m_{0}$ increases to $0.54$, the nonlinear localized mode may be unstable, as in Fig. 3(c2), even if the $\mathcal{PT}$ symmetry of the linear state remains unbroken. Most interesting, if we increase $\theta_{0}$ to $0.36$ in the same case, the nonlinear localized mode restabilizes itself, while the $\mathcal{PT}$ symmetry of the linear state gets broken, see Fig. 3(c3). Another interesting feature is that the real part of the corresponding complex potential may now exhibit a double-well shape, as shown in Fig. 3(c4).

Interactions between exact bright solitons and sech pulses. To additionally test the robustness of the exact bright solitons (6), we simulated their collisions with boosted (moving) sech-shaped solitary pulses. For $\alpha=1,2,3$, we respectively choose $(m_{0},\theta_{0})=(0.5,0.1),(0.4,0.1),(0.2,0.1)$ with the unbroken $\mathcal{PT}$ symmetry of the linear state and consider initial conditions $\psi_{\alpha}(x,0)=\phi_{\alpha}(x)+\mathrm{sech}(x+20)e^{1.8ix}$ with $\phi_{\alpha}(x)$ given by Eq. (6), and $e^{1.8ix}$ imposing the boost onto the sech pulse with the unitary amplitude and initial position at $x=-20$. Direct simulations demonstrate that the exact nonlinear modes $\phi_{\alpha}(x)$ propagate steadily without any change of shape and velocity after the collision, see Figs. 4(a,c,e). In other simulations, we chose values $(m_{0},\theta_{0})=(0.75,0.1),(0.8,0.3),(0.54,0.36)$, which correspond to the linear modes with broken $\mathcal{PT}$ symmetry, and took the corresponding initial conditions as $\psi_{1}(x,0)=\phi_{1}(x)+\mathrm{sech}(x+40)e^{2.2ix}$ with $\phi_{1}(x)$ given by Eq. (6) for $\alpha=1$, as well as $\psi_{2,3}(x,0)=\phi_{2,3}(x)+\mathrm{sech}(x+40)e^{4ix}$, with $\phi_{2,3}(x)$ given by Eq. (6) for $\alpha=2,3$. In the case of $\alpha=1$, the exact soliton still demonstrates the steady propagation after the collision, while in the cases of $\alpha=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, $m_{0}\rightarrow m_{0}(z)$ or $\theta_{0}\rightarrow\theta_{0}(z)$, cf. Refs. [38, 53, 54, 55]. We choose the following temporal-modulation pattern:

[TABLE]

This implies replacing Eq. (1) by

[TABLE]

where $m(x,z),V(x,z),W(x,z)$ are respectively given by Eqs. [(2), (3), (4)] with $m_{0}\rightarrow m_{0}(z)$ and $\theta_{0}\rightarrow\theta_{0}(z)$.

For the cases of $\alpha=1,2,3$, Figs. 5(a1,b1,c1), respectively, exhibit stable switch of the nonlinear localized modes $\psi(x,z)$ governed by Eq. (8), starting from the initial condition given by Eq. (6). In these figures, we demonstrate the transformation of initially stable nonlinear localized modes for $(m_{01},\theta_{01})=(0.5,0.1),(m_{01},\theta_{01})=(0.4,0.1),(m_{01},\theta_{01})=(0.2,0.1)$, with $\alpha=1,2,3$ to ones corresponding, respectively, to $(m_{02},\theta_{02})=(0.5,0.3),(m_{02},\theta_{02})=(0.4,0.5),(m_{02},\theta_{02})=(0.2,0.36)$, i.e., $m_{0}$ is fixed, while $\theta_{0}$ varies. In these cases, both the initial and final parameter values correspond to linear modes with unbroken $\mathcal{PT}$ symmetry. Likewise, fixing $\theta_{0}$ and setting $m_{0}\rightarrow 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 $\alpha=2$, also ending with parameters corresponding to the linear state with unbroken $\mathcal{PT}$ symmetry.

On the other hand, Fig. 5(a2) reveals unstable transformation of the nonlinear modes (6) with $\alpha=1$, and Fig. 5(c2) exhibits stable transformation of the nonlinear modes (6) with $\alpha=3$, both from initial parameters corresponding to the linear state with unbroken $\mathcal{PT}$ symmetry to final ones corresponding to broken $\mathcal{PT}$ symmetry in the linear state. Finally, the transformation is also implemented by varying both $m_{0}$ and $\theta_{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 $\theta_{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 $\mathcal{PT}$ symmetry of the respective linear state, to another exact nonlinear mode, which corresponds to the broken $\mathcal{PT}$ symmetry in the linear state, varying both parameters. The stable outcome of the simultaneous variation of $m_{0}$ and $\theta_{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, a2, a3). If both separate excitations produce stable outputs, their simultaneous action may give either an unstable output, see Figs. 5(b1,b2,b3), or a stable one, see Figs. 5(c1,c2,c3)).

1D constant-intensity waves. Recently, a class of complex potentials that are more general than the $\mathcal{PT}$-symmetric complex-valued ones has been put forward [39, 78, 79, 80]. Similar to the $\mathcal{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; see Refs. [39, 78, 79, 80]. The real and imaginary parts of these potentials are defined as [81] $V(x)=-v^{2}(x)$, $W(x)=dv(x)/dx$, where $v(x)$ 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

[TABLE]

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 $\mathcal{PT}$-symmetric one. In a more general case, when the potential is not $\mathcal{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)

[TABLE]

where $C$ is a real constant amplitude. The power flow (the Poynting vector) corresponding to the CW solution (10) is $S(x)=C^{2}v(x)$. To examine dynamical properties of the CW solution (10), we consider three different types of the complex potentials (9), in which the diffraction coefficient $m(x)$ is chosen as in Eq. (2), and $v(x)$ is taken as a Hermite-Gauss (HG) function, $H_{n}(x)e^{-\omega x^{2}/2}$ ($H_{n}(x)$ is the Hermite polynomial, and $\omega>0$ is a frequency), or as a simple periodic one, $v_{0}\mathrm{cos}(x)$. Generic examples of this class of complex potentials correspond to $m(x)=0.5\mathrm{sech}(x)+1$, with $v(x)=H_{2}(x)e^{-x^{2}/2}$ (a $\mathcal{PT}$-symmetric HG form), or $v(x)=H_{3}(x)e^{-x^{2}/2}$ (a non-$\mathcal{PT}$-symmetric HG form), or, finally, $v(x)=\mathrm{cos}(x)$ (a $\mathcal{PT}$-symmetric periodic form). The profiles and gain-loss regions of these three complex potentials are displayed in Figs. 6a, 7(a), 8(a), respectively. In the linear limit ($g=0$) in these three cases, if the CW inputs are taken merely as $\psi(x,0)=C$, 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), 7c, and 8(c). It is worthy to note that the beam in the $\mathcal{PT}$-symmetric system with the HG complex potential steadily propagates farther than that in the case of the non-$\mathcal{PT}$ HG complex potential, as the $\mathcal{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 increases too under the action of the $\mathcal{PT}$-symmetric HG potential (see Fig. 6d), but it does not increase in the case of the non-$\mathcal{PT}$-symmetric HG potential (see Fig. 7d).

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. 6e, 7e, and 8e for the self-focusing nonlinearity, and in Figs. 6(f), 7(f), and 8(f) for the self-defocusing nonlinearity.

2D $\mathcal{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 $\mathcal{PT}$-symmetric setting. In this case, the field evolution is governed by the 2D NLS equation with a $\mathcal{PT}$-symmetric potential:

[TABLE]

where $\nabla$ is the 2D gradient $(\partial_{x},\partial_{y})$, and $\mathbf{m}(x,y)$ is a $2\times 2$ matrix function, which we take in the diagonal form, $m_{\alpha}(x,y)=\mathrm{diag}[m_{\alpha}(x),m_{\alpha}(y)]$ with $m_{\alpha}(x)$ given by Eq. (2). We look for stationary 2D modes in the form of $\psi=\phi(x,y)e^{i\mu z}$ with $\phi(x,y)=u(x,y)e^{i\theta(x,y)}$, where $u(x,y)$ and $\theta(x,y)$ are real amplitude and phase, respectively, which satisfy stationary equations:

[TABLE]

2D $\mathcal{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 $\alpha=1,2,$ and $3$, are given below in Eqs. (24a), (25a), and (26a), respectively, along with the exact solutions (see Methods).

In particular, similar to 1D case, it is possible to find stationary constant-intensity (CW) solutions of Eq. (11)

[TABLE]

(cf. Eq. (10)) for a family of 2D complex potentials similar to their 1D counterparts (9):

[TABLE]

where $\theta(x,y)$ is an arbitrary real function. 2D CW solutions (14) will be considered in detail elsewhere.

Comparison of exact 2D solitons and numerical solutions. For $\alpha=1,2,3$, we first choose three parameter sets, $(m_{0},\theta_{0})=(0.5,0.1),(0.4,0.1),(0.2,0.1)$, with the corresponding 2D $\mathcal{PT}$-symmetric complex potentials shown in Figs. 9(a,b), 10(a,b), and 11(a,b), respectively. The corresponding exact soliton solutions, given by Eqs. (24a), (25a), and (26a), with $\mu=9/2,8,$ and $25/2$, are shown in the second row of Figs. 9, 10, 11, respectively.

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, 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 $\mu$ and identify the existence ranges of the numerically found solitons, as shown in Fig. 12(a), where the lower cutoffs are $\mu_{1}\approx 3.4,\mu_{2}\approx 6.3,\mu_{3}\approx 11.1$, and there are no finite upper cutoffs. These numerical soliton solutions are more often than not unstable, especially for larger values of $\mu$. We display their linear-stability spectra in Figs. 12(b,c,d) around $\mu=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, $(\alpha,m_{0},\theta_{0})=(1,0.5,0.1)$, almost all the corresponding numerical and exact solitons are unstable; however, in the second and third cases, viz., $(\alpha,m_{0},\theta_{0})=(2,0.4,0.1)$ and $(3,0.2,0.1)$, there exist some stability regions, which are distributed around $\mu=8$ and $25/2$. To confirm these linear-stability results, in what follows we investigate the dynamics by dint of direct numerical simulations.

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 $\alpha=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\approx 280$ (Fig. 13(c)). However, for $\alpha=2$ and $3$, both Figs. 13(e) and 13(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) and 13(h), which are identical to the corresponding initial profiles in Fig. 10 and Fig. 11.

We have also numerically investigated the stability of the solitons at the lowest-power points $\mu_{1}=3.4,\mu_{2}=6.3,$ and $\mu_{3}=11.1$ for $\alpha=1,2,$ and $3$, respectively. They all turn out to be unstable, although they may seem stable at a shorter propagation distance. Figure 14 specifically exhibits their evolution and the corresponding intermediate and final states. These results indicate a somewhat surprising fact that the localized nonlinear mode with the lowest power is not necessarily stable. On the other hand systematic simulations demonstrate that the numerically found solitons are fully stable in the vicinity of the corresponding exact nonlinear modes, that is, near $\mu_{1}=9/2,\mu_{2}=8,\mu_{3}=25/2$ for $\alpha=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,c,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 $\mathcal{PT}$-symmetric generalized Scarf-II potentials and variable (position dependent) diffraction coefficients. First, in the linear version of the model, linear problem, parameter regions of the unbroken and broken $\mathcal{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 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 $\mathcal{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 $\mathcal{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 $\mathcal{PT}$-symmetric ones, but feature similar properties. These theoretical results suggest new experiments for $\mathcal{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, $\psi(x,z)=\phi(x)e^{i\mu z}$, where $-\mu$ is the propagation constant, and the localized complex wave function satisfies the following ordinary differential equation with the $x$-dependent diffraction coefficient:

[TABLE]

We take the complex wave function as $\phi(x)=u(x)\exp\left(i\int_{0}^{x}v(s)ds\right)$ with real amplitude $u(x)$, and superfluid phase velocity given by

[TABLE]

The amplitude satisfies the following second-order ordinary differential equation:

[TABLE]

which may be transformed by setting $\tilde{u}(x)\equiv m(x)u_{x}$ into a system of coupled first-order equations:

[TABLE]

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 $\phi(x)=\lambda w(x)$, where $\lambda$ 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:

[TABLE]

where $F_{1}=\mathcal{F}\{[m_{x}\mathcal{F}^{-1}(ik_{x}\hat{w})-(V+iW)w-\mu w+pmw]/m\}$, $F_{2}=\mathcal{F}(g|w|^{2}w/m),\lambda^{2}=\int_{-\infty}^{+\infty}[(k_{x}^{2}+p)|\hat{w}|^{2}-F_{1}\hat{w}^{\ast}]\mathrm{d}k_{x}/\int_{-\infty}^{+\infty}F_{2}\hat{w}^{\ast}\mathrm{d}k_{x}$, $\mathcal{F}$ 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 $\phi(x)=\lambda\mathcal{F}^{-1}[\hat{w}(x)]$.

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:

[TABLE]

where $F_{1}=\mathcal{F}\{[m(x)-1]\mathcal{F}^{-1}(-k_{x}^{2}\hat{w})+[m(y)-1]\mathcal{F}^{-1}(-k_{y}^{2}\hat{w})+m_{x}(x)\mathcal{F}^{-1}(ik_{x}\hat{w})+m_{y}(y)\mathcal{F}^{-1}(ik_{y}\hat{w})-(V+iW)w-\mu w+pw\}$, $\mathcal{F}$ 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 $\mathcal{PT}$-symmetric potential $V(x)+iW(x)$, one may solve Eq. (17) (or equivalently, Eqs. (18) and (19)), to obtain stationary soliton solutions $\phi(x)$, by analytical or the above-mentioned methods. Then localized nonlinear modes of Eq. (1) can be found in the stationary form, as $\psi(x,z)=\phi(x)e^{i\mu z}$. To explore the linear stability of the localized modes in the 1D case, we consider a perturbed solution [87, 88],

[TABLE]

where $\epsilon$ is an infinitesimal perturbation amplitude, $F(x)$ and $G(x)$ are the eigenfunctions of the linearized problem, and $-\delta$ is the respective eigenfrequency, the instability taking place if some eigenvalues are not purely real. Inserting the perturbed solution (22) into Eq. (1) and linearizing with respect to $\epsilon$, we obtain the following eigenvalue problem:

[TABLE]

where $\hat{L}_{1}\equiv\partial_{x}\left(m(x)\partial_{x}\right)-[V(x)+iW(x)]+2g|\phi|^{2}-\mu$. The $\mathcal{PT}$-symmetric nonlinear modes are linearly stable provided $\delta$ has no imaginary part, otherwise they are linearly unstable. The whole stability spectrum $\delta$ can be numerically calculated by the Fourier collocation method (see [88]).

For the 2D case (11), the operator $\hat{L}_{1}$ in Eq. (23) is changed to $\hat{L}_{1}=\partial_{x}[m(x)\partial_{x}]+\partial_{y}[m(y)\partial_{y}]-[V(x,y)+iW(x,y)]+2g|\phi|^{2}-\mu$. Other technical details are similar to those in the 1D case. As concerns the numerical computation of the full stability spectrum in 2D, 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},\theta_{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 $\mu$ (see Figs. 12(b, c, 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 2D, the corresponding $\mathcal{PT}$-symmetric linear spectra are not exhibited further. Nevertheless, through repeated numerical tests, we find that it is instructive that the $\mathcal{PT}$-symmetric breaking curves in 1D can provide powerful reference for those in 2D.

2D nonlinear modes for the $\mathcal{PT}$-symmetric potentials with different parameters $\alpha$. The 2D $\mathcal{PT}$-symmetric potentials and exact solutions of Eqs. (12) and (13), which they admit, are listed as follows:

Case 1 ($\alpha=1$):

[TABLE]

Case 2 ($\alpha=2$):

[TABLE]

Case 3 ($\alpha=3$):

[TABLE]