Solitonic dynamics and excitations of the nonlinear Schrödinger equation with third-order dispersion in non-Hermitian PT-symmetric potentials

Solitons are of the important significant in many fields of nonlinear science such as nonlinear optics, Bose-Einstein condensates, plamas physics, biology, fluid mechanics, and etc. The stable solitons have been captured not only theoretically and experimentally in both linear and nonlinear Schrödinger (NLS) equations in the presence of non-Hermitian potentials since the concept of the parity-time -symmetry was introduced in 1998. In this paper, we present novel bright solitons of the NLS equation with third-order dispersion in some complex -symmetric potentials (e.g., physically relevant -symmetric Scarff-II-like and harmonic-Gaussian potentials). We find stable nonlinear modes even if the respective linear -symmetric phases are broken. Moreover, we also use the adiabatic changes of the control parameters to excite the initial modes related to exact solitons to reach stable nonlinear modes. The elastic interactions of two solitons are exhibited in the third-order NLS equation with -symmetric potentials. Our results predict the dynamical phenomena of soliton equations in the presence of third-order dispersion and -symmetric potentials arising in nonlinear fiber optics and other physically relevant fields.

The representative nonlinear Schrödinger (NLS) equation can be used to describe distinguishing wave phenomena arising in many nonlinear physical fields such as nonlinear optics, Bose-Einstein condensates (alias Gross-Pitaevskii equation), the deep ocean, DNA, plasmas physics, and even financial market, etc. [1][2][3][4][5][6][7][8][9] . The dynamics of fundamental bright and dark solitons (also vortices and light bullets in higher-dimensional cases) of the NLS model and self-similar modes in its generalized forms have been addressed (see refs 10-17 and references therein). Recently, inspired by the non-Hermitian PT -symmetric potentials first suggested by Bender and Boettcher 18,19 in the classical Hamiltonian operators, Musslimani et al. 20 first introduced the complex PT -symmetric potentials in the NLS model such that some novel phenomena with stable modes were found. After that, a variety of distinguishing PT -symmetric or non-PT -symmetric potentials were introduced to the continuous or discrete NLS equations to explore dynamical behaviors of PT -symmetric nonlinear modes (see, e.g., refs 21-42 and references therein). Meanwhile, more and more physical experiments have also been designed to observe new wave phenomena in the sense of non-Hermitian PT -symmetric potentials [43][44][45][46][47][48] . Here the parity  and temporal  operators are defined as 19 : . The one-dimensional complex potential U(x) is PT -symmetric provided that the sufficient (not necessary) conditions U R (x) = U R (− x) and U I (− x) = − U I (x) hold 19 , where U(x) is also called the refractive-index in optical fibre.
In the study of ultra-short (e.g., 100 fs 1 ) optical pulse propagation, the higher-order dispersive and nonlinear effects become significant such as third-order dispersion (TOD), self-steepening (SS), and the self-frequency shift (SFS) arising from the stimulated Raman scattering. The third-order NLS equation was introduced from the Maxwell equation 49,50 . The generalized inhomogeneous third-order NLS equation with modulating coefficients in where Raman effect, nonlinear dispersion terms (e.g., self-steepening term and self-frequency shift effect), and higher-order dispersion terms are neglected 49,50,55,56 , ψ ≡ ψ(x, z) is a complex wave function of x, z, z denotes the propagation distance, the real parameter β stands for the coefficient of TOD, the PT -symmetric potential requires that V(x) = V(− x) and W(x) = − W(− x) describing the real-valued external potential and gain-and-loss distribution, respectively, and g > 0 (or < 0) is real-valued inhomogeneous self-focusing (or defocusing) nonlinearity. The power of Eq. (1) is given by 2 and one can readily know that . Equation (1) becomes the usual higher-order NLS equation in the absence of the gain-and-loss distribution 53 . Equation (1) with β = 0 becomes the PT -symmetric nonlinear model, which has been studied [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40] . In the following we consider the case in the presence of TOD term (β ≠ 0) and gain-and-loss distribution. Here our following results are also suitable for the case x → t in Eq. (1).
Linear spectrum problem with PT-symmetric potential. We start to study the physically interesting potential in Eq. (1) as the PT -symmetric Scarff-II-like potential where the real constant V 0 and TOD parameter β can be used to modulate the amplitudes of the reflectionless potential V(x) 57 and gain-and-loss distribution W(x), respectively. Moreover, V(x) and W(x) are both bounded (i.e., 0 and V(x), W(x) → 0 as |x| → ∞ (see Fig. 1a). It is easy to see that the gain-and-loss distribution is always balanced since ∫ = −∞ +∞ W x dx ( ) 0 and has only the limit effect on linear and nonlinear modes since W(x) ~ 0 as |x| > M > 0. The sole difference between the potential (2) and the usual Scraff-II potential 58 is that the gain-and-loss distribution in Eq. (2) will more quickly approaches to zero than one (i.e., β sec h x tan h x) in Scarff-II potential for the same amplitudes.
We firstly consider the linear spectrum problem (i.e., Eq. (1) with g = 0) in the PT -symmetric Scarff-II-like potential (2) and use the stationary solution transformation ψ(x, z) = Φ (x)e −iλz to yield where λ and Φ (x) are the corresponding eigenvalue and eigenfunction, respectively, and lim |x|→∞ Φ (x) = 0. Since the discrete spectrum of a complex PT -symmetric potential is either real or appears in complex conjugated pairs, thus we may find some proper parameters V 0 , β for which the complex PT -symmetric potential keep unbroken.
Here we consider V 0 < 0 such that the shape of potential V(x) seems to be V-shaped with zero boundary conditions (see Fig. 1a). We numerically study the discrete spectra of the operator L (see Methods). Figure 1b exhibits the regions of broken and unbroken PT -symmetric phases on the (V 0 , β) space. Two almost parallel straight lines (β ≈ ± 0.12) separate the limited space {(V 0 , β)|− 0.02 ≤ V 0 ≤ − 3, |β| ≤ 0.5}. The regions of broken and unbroken PT -symmetric phases are outside and inside between two lines, respectively. For the given TOD parameter β = 0.1 and varying V 0 , the spontaneous symmetry breaking does not occur from the two lowest states since the maximum absolute value of imaginary parts of λ is less than 6 × 10 −14 and they can be regarded as zero (see Fig. 1(c,d)). However, for the given V 0 = − 2 and varying β, the spontaneous symmetry breaking occurs from two lowest states starting from some β = 0.12 (see Fig. 1(e,f)).

Nonlinear localized modes and stability.
For the given PT -symmetric Scarff-II-like potential (2), based on some transformations, we can find the unified analytical bright solitons of Eq. (1) for both the self-focusing and defocusing cases (see Methods) with v = ± 1, the potential is μ = (3κ 2 + 3βκ − βκ 3 − 3)/6, and the existent condition g(V 0 − κβ + 1) > 0 is required. For the signs of parameter v and nonlinearity g, we find the following four cases for the existent conditions of bright solitons (4) (let In the following we numerically 59 study the robustness (linear stability) of nonlinear localized modes (4) for both self-focusing and defocusing cases (g = ± 1) via the direct propagation of the initially stationary state (4) with a noise perturbation of order about 2% in Fig. 2 (see Methods).  Fig. 1b], the stable nonlinear mode is generated (see Fig. 2(a2)). For the fixed V 0 = − 0.8, if we change β = 1.1 corresponding to the domain of the broken linear PT -symmetric phase, then a stable nonlinear mode is found too, that is, the focusing nonlinear term can modulate the unstable linear modes (broken PT -symmetric phase) to stable nonlinear modes (see Fig. 2(a3)). But if β increases a little bit (e.g., β = 1.5), then the nonlinear mode becomes unstable (see Fig. 2(a4)). In particular, for V 0 = − 1.1, β = 0.7, which does not satisfy the required existent condition of solution V 0 > − α, that is, the expression (4) for this case, ψ = . − .
is not an analytical solution of Eq. (1), but we still use it as the initial solution with a noise perturbation of order 2% to make numerical simulations such that we find the initial mode φ 0 (x, 0) can be excited to a stable and weakly oscillatory (breather-like behavior) situation (see Fig. 2(a5)). For Case (ii) v = g = − 1, we also have the similar results (see Fig. 2(b3-b5)). For the last two Cases (iii) and (iv), we fix V 0 > 0, in which the potential is similar to a Gaussian-like profile and the linear problem (3) has no discrete spectra, but we still find the stable nonlinear modes (see Fig. 2(c3-c5), (d3,d5)) and unstable (see Fig. 2(c2,d2,d4)) nonlinear modes. For the above-obtained nonlinear modes (4), we have the corresponding transverse power-flow or Poynting vector given by . We here consider the only case β > 0 (the case β < 0 can also be considered). For v = 1 (or − 1), we have S > 0 (or < 0), which implies that the pamaeter v change the directions of power flows from gain to loss regions. The power of the solutions (4) is , which is conserved.

Exciting stable nonlinear localized modes in equation (1).
Nowadays we turn to the excitation of nonlinear modes by means of a slow change of the control TOD parameter β → β(z) in Eq. (1) which is regarded as a function of propagation distance z, that is, we focus on simultaneous adiabatic switching on the TOD and gain-and-loss distribution, modeled by  ( /2000) , 0 1000, , 1000 with β 1,2 being real constants. It is easy to see that the solutions (4) with β → β(z) do not satisfy Eq. (5), but for z = 0 and z ≥ 1000 solutions (4) indeed satisfy Eq. (5). Figure 4 displays the wave propagation of the solutions ψ(x, z) of Eq. (5) subject to the initial condition given by Eq. (4) with β → β(z) given by Eq. (6). For v = 1 and different parameters g, V 0 , β 1,2 , Fig. 4(a-c) exhibit the stable modes in which the initial states |ψ(x, 0)| 2 given by Eq. (4) with z = 0, β = β 1 are all of the higher amplitudes and then the amplitudes decrease slowly as z increases such that they reach the alternative stable sates beginning from about z = 1000. For v = − 1 and different parameters g, V 0 , β 1,2 , Fig. 4(e,f) also exhibit the stable modes in which the initial states are all of the lower amplitudes and then the amplitudes grow step and step as z increases such that they reach the stable and weakly oscillatory (breather-like behavior) situations beginning from about z = 1000, but Fig. 4d shows that the stable mode keeps from z = 0 to z = 1100 and then the wave slowly increases a little bit to reach another stable and weakly oscillatory (breather-like behavior) feature. In particular, in Fig. 4(b,f), we can excite the initial states subject to inexact solitons (4) of Eq. Nonlinear model with the spatially varying TOD. We here consider nonlinear modes of the generalized form of Eq. (1) with x-spatially varying TOD coefficient β → β(x), that is, Scientific RepoRts | 6:23478 | DOI: 10.1038/srep23478 x 0 2 with β 0 ≠ 0 being a real amplitude of the Gaussian profile and another complex PT -symmetric harmonic-Gaussian potential βγ γ β γ βγ in which β = β(x) is given by Eq. (8) and we have introduced the function γ σ = − e x /2 2 with σ ≠ 0 being a real constant. In fact, we can also consider the general case β(x). We know that the potential V(x) approaches to the harmonic potential x 2 /2 (which is related to the usual physical experiments) and W(x) → 0 as |x| → ∞ (see Fig. 6(b,d)).
For the given PT -symmetric harmonic-Gaussian potential (9), we studied the discrete spectra of linear problem (3) such that we give the separated regions for the broken and unbroken PT -symmetric phase (see Fig. 5(a)). When |σ| is greater than about 0.16, we only find the discrete spectra for β 0 very approaching to zero (here we consider |β 0 | ≤ 0.3). For the given amplitude β 0 = 0.1 of TOD coefficient, we give the first six lowest eigenvalues (see Fig. 5(b,c)), where the spontaneous symmetry breaking occurs from the two lowest states.
We find the exactly analytical solutions of Eq. (7) with the Gaussian TOD coefficient in the PT -symmetric harmonic-Gaussian potential (9) in the form where g > 0 and erf(·) is an error function. In the following we take g = 1 without loss of generality.
In what follows we numerically investigate the robustness (linear stability) of nonlinear localized modes (10) for self-focusing case (g = 1) via the direct propagation of the initially stationary state (10) with a noise perturbation of order about 2% in Fig. 6 (see Methods). The linear stability of the solutions (10) is displayed in Fig. 6a) such that we find the solutions (10) are possibly stable nearby σ = 0. For σ = β 0 = 0.1, in which the potential becomes almost a harmonic potential x 2 /2 (see Fig. 6b), we find the stable nonlinear mode (see Fig. 6c). But for σ = 0.2, β 0 = 0.35, in which the potential is a double-well potential (see Fig. 6d), we also obtain the stable nonlinear mode (see Fig. 6e).
For the above-obtained nonlinear modes (10), we have the corresponding transverse power-flow or Poynting vector given by , which implies that sgn(S) = sgn(σ). Since the  gain-and-loss distribution W(x) given by Eq. (9) depends on TOD parameter β 0 and σ, which generate that there are more one intervals for gain (or loss) distribution, thus though the power flows along the positive (negative) direction for the x axis as σ > 0 (< 0), it is of the complicated structures from the gain-and-loss view.
Moreover, we also study the interactions of two bright solitons in the PT -symmetric potential. which makes the potential V(x) and gain-or-loss distribution W(x) given by Eq. (9), and solutions (10) change, but it does not change the TOD coefficient given by Eq. (8). We consider the wave evolution of the solution ψ(x, z) satisfied by Eq. (5) with Eqs (8), (9) and (12) subject to the initial condition given by Eq. (10) such that we find the stable nonlinear modes by using the excitation, that is, we can smoothly excite one stable mode (see Fig. 6(c)) to another stable mode (see Fig. 8(b)). Moreover, if we simultaneously consider the effect of varying β 0 (z) and σ(z) given by Eqs (11) and (12) then we also find the similar result (see Fig. 8(c)).   12), (c) β 0 (z) and σ(z) given by Eqs (11) and (12).

Discussion
In conclusions, we have introduced some non-Hermitian (e.g., complex PT -symmetric) potentials in the nonlinear Schrödinger equation with third-order dispersion. For the chosen physically interesting PT -symmetric Scarff-II-like and harmonic-Gaussian potentials, we found exact analytical bright solitons of this equation. In the presence of these PT -symmetric potentials, we study the broken and unbroken of PT -symmetric phases of the corresponding linear problem (third-order linear operator with complex potentials) for TOD and potential parameters such that we find the TOD parameter has a strong effect on the spectra (it only admit a few discrete spectra). We have studied the linear stability of exact bright solitons. In particular, we find the stable nonlinear modes for some control parameters for which even if the corresponding linear PT -symmetric phase is broken. Moreover, we also investigate the problems of nonlinear modes excitations, which can excite initial nonlinear modes to reach stable modes. The method in this paper can also be extended to explore other higher-order or/and higher-dimensional NLS equations in the presence of non-Hermitian potentials and may open a new window to investigate similar problems. Our results may be useful to provide theoretical researchers and experimental scientists with more new data about the PT -symmetric nonlinear modes in higher-order nonlinear wave models.

Methods
Linear spectrum problem. For Eq. (1) in the absence of nonlinear term (g = 0), we assume that ψ(x, z) = Φ (x)e −iλz , then we have Eq. (3), which with the PT -symmetric potential (2), as |x| → ∞ , reduces to , whose characteristic equation , whose roots, in general, are complex numbers and complicated. For example, if λ = 1/(3β 2 ) and β > 0 (without loss of generality), we have its three roots , which probably lead to the result that the corresponding eigenfunctions should satisfy periodic boundary conditions. Additionally, if β depends on the space x, e.g., β(x) = β 0 exp(− x 2 ), and V(x), W(x) are given by Eq. (9) , which generally corresponds to zero boundary conditions and discrete spectra. Therefore, in order to verify these results, we use the Fourier collocation method [59][60][61] to numerically study the above-mentioned linear spectrum problems and obtain the agreeable conclusions as ones by the theoretical analysis.
Nonlinear stationary modes. We consider the stationary solutions of Eq. (1) in the form ψ(x, z) = φ(x) e −iμz , where φ(x) is a complex field function and μ the corresponding propagation constant. We have For the given PT -symmetric potential (2) we can find the exact bright solitons (4) of Eq. (1). Similarly, we can also find the solutions of Eq. (5) in the PT -symmetric potential (9).
Linear stability of nonlinear stationary modes. To further study the linear stability of the above-obtained nonlinear stationary solutions ψ(x, z) = φ(x)e −iμz of Eq. (1), we considered a perturbed solution 59,62  , δ and F(x) and G(x) are the eigenvalue and eigenfunctions of the linearized eigenvalue problem. We substitute the expression into Eq. (1) and linearize it with respect to ϵ to yield the following linear eigenvalue problem . It is easy to see that the PT -symmetric nonlinear modes are linearly stable if δ has no imaginary component, otherwise they are linearly unstable.