Small-amplitude nonlinear modes under the combined effect of the parabolic potential, nonlocality and ${\cal PT}$ symmetry

We consider nonlinear modes of the nonlinear Schr\"odinger equation with a nonlocal nonlinearity and an additional PT-symmetric parabolic potential. We show that there exists a set of continuous families of nonlinear modes and study their linear stability in the limit of small nonlinearity. It is demonstrated that either PT symmetry or the nonlocality can be used to manage the stability of the small-amplitude nonlinear modes. The stability properties are also found to depend on the particular shape of the nonlocal kernel. Additional numerical simulations show that the stability results remain valid not only for the infinitesimally small nonlinear modes, but also for the modes of finite amplitude


I. INTRODUCTION
Nonlocal nonlinearities are known to be of relevance for numerous physical applications. In optical applications, nonlocal nonlinearity emerges in description of beam propagation in thermal media [1] or in atomic vapors [2]. Nonlocal nonlinearities have been shown to support optical solitons [3][4][5][6] and to alter their properties, such as, for instance, stability [3], the collapse dynamics [7], or the interactions between solitons [8].
Nonlocality is also relevant in the theory of Bose-Einstein condensates (BECs) [9] in view of the intrinsically nonlocal nature of the inter-atomic interactions [10,11]. Nonlocal nonlinearities are also inherent for plasma problems [12] and for the theory on nematic crystals [13].
In many relevant situations, nonlocal nonlinearity is accompanied by a certain potential which models a confining trap for a BEC or a refractive index modulation in an optical waveguide. Presence of the external potential or nonlocality enriches the variety of possible nonlinear modes allowing to observe the so-called higher (alias multi-pole) nonlinear modes [5,14] apart from the best studied fundamental (alias ground state) nonlinear modes [4,15]. Nonlocal nonlinearities have been considered in combination with the double-well [16] and the periodic potentials [17][18][19]. It looks however that a systematic study of nonlocal nonlinear modes in a parabolic potential has not been performed yet. In our present study, we make a step towards filling this gap and investigate nonlocal nonlinear modes in a complex parabolic PT -symmetric potential of the form V (x) = x 2 − 2iαx, where α is a real coefficient. For α = 0 the potential is reduced to the real parabolic potential whose modes are fairly well studied in the local case [14,20], but have attracted much less attention in the nonlocal context.
For a nonzero α, the potential V (x) is complex and PT symmetric [21], which is reflected by the fact that the real and imaginary parts of V (x) are even and odd functions of x, respectively. PT symmetry can be implemented in diverse physical setting, including optics [22], atomic gases [23], and BECs [24].
Nonlinear modes and solitons have been studied in various PT -symmetric potentials with local [25][26][27] and nonlocal [28] nonlinearities, but the present work has a slightly peculiar focus. We are particularly interested in the examining the combined effect of the parabolic trapping, nonlocality and PT symmetry on the excited (higher-order) nonlinear modes. The focus of our study is additionally shifted towards the existence and stability of small-amplitude nonlinear modes, whose properties can be studied by means of asymptotic expansions describing bifurcations of nonlinear modes from the eigenstates of the underlying linear model. The linear spectrum of the PT -symmetric parabolic potential is available in the analytical, which makes it possible to reduce the stability problem to computing the real zeros of certain analytical expressions. Our analysis shows that either the degree of the nonlocality or the strength of the PT symmetry can be used to manage the stability of the excited modes. On the other hand, we reveal a dependence of the stability properties on the particular shape of the nonlocal kernel. We also touch upon the properties of nonlinear modes of finite amplitude and illustrate numerically that the stability results hold not only for the infinitesimally small nonlinear modes but also remain valid in a finite range of the amplitudes of nonlinear modes.
The paper is organized as follows. In the next section, we present a model to be considered throughout the paper. In Sec. III we describe the linear properties of the model and introduce the families of small-amplitude nonlinear modes. In Sec. IV we describe the methodology of the linear stability analysis. The main results of stability analysis are presented in Sec. V and Sec. VI for the conservative and PT -symmetric cases, respectively. In Sec. VII we demonstrate numerically that the presented asymptotical stability results persist in a finite range of the amplitude of the modes. Concluding Sec. VIII outlines the main results of the work.

II. THE MODELS
We use the dimensionless nonlocal nonlinear Schrödinger equation with an additional parabolic potential and nonlocal nonlinearity: The case α = 0 corresponds to the real parabolic potential, while nonzero α governs the strength of the PT symmetry. The system is subject to gain in the domain x < 0 and is subject to losses for x > 0. In the optical context, equation (1) models beam guidance in a medium whose refractive index n(x) = n r (x) + in i (x) has parabolic modulation of the real part n r (x) = x 2 , and linear modulation of the imaginary part n i = 2αx.
Equation (1) with α = 0 also arises as a Gross-Pitaevskii equation [9] in the mean-field theory the harmonically trapped a quasi-one-dimensional BEC [10,19,29]. In this case, q(x, z) is a quasi-one-dimensional wavefunction [11], and z plays a role of time. Sign σ = ±1 corresponds to attractive (focusing) and repulsive (defocusing) nonlinearities, respectively. Nonlocal properties of the model are described by the kernel K(x) which is a positive function with normalization K(x)dx = 1 (hereafter all the integrals are taken over the whole real axis). Following to the previous studies on the subject [5,6,16,18], we pay the primary attention to the following two kernels: (i) The smooth Gaussian kernel where λ > 0 governs strength of nonlocality.
(ii) Exponential kernel with the singularity at x = 0: (iii) It will be also interesting to compare the obtained results with a slowly (algebraically) decaying kernel. While some studies (for example, [16,19]) deal with the so-called cut-off (CO) kernel which decays as |x| −3 as x → ±∞, in our study we choose the algebraically decaying Lorentzian kernel (the Lorentzian) This choice is motivated by the fact that the Lorentzian kernel (4) allows for the relatively simple analytical evaluation of the convolution integrals that arise in what follows.

III. BIFURCATIONS OF NONLINEAR MODES
The stationary nonlinear modes for Eq. (1) admit the representation q(x, z) = e iβz w(x), where β is a real parameter. It is convenient to introduce the shifted parameter b defined through the relation β = b − α 2 , which allows to rewrite the equation for the nonlinear modes in the following form: subject to the zero boundary conditions at the infinity: lim x→±∞ w(x) = 0. Let us first recall the properties of the underlying linear problem which formally corresponds to σ = 0 in Eq. (5). The resulting equation which is considered as a linear eigenvalue problem for b and w(x), has a purely real and discrete spectrum [27,30,31]. The eigenvaluesb n can be listed asb n = −(2n + 1), n = 0, 1, . . .. The corresponding eigenfunctionsw n (x) read where H n (x) and G n (x) are Hermite and Laguerre polynomials [32], respectively. The choice of the coefficients c n implies normalization |w n | 2 dx = 1. Each Hermite polynomial H n (x) contains only powers of x of the same parity as the number n. This implies that in the conservative limit (α = 0) the eigenfunctions are real-valued functions which have exactly n zeros and have the same parity as the number n. For α = 0, each eigenfunctioñ w n (x) possesses real even part and odd imaginary part for even n and odd real part and even imaginary part for odd n.
Returning now to the nonlinear problem, we look for bifurcations of small-amplitude nonlinear modes from the linear eigenstates. To this end, we introduce the formal asymptotic expansions which describe nonlinear modes w n (x) branching off from the zero solution w(x) ≡ 0 at the eigenvaluesb n [27]: where ε ≪ 1 is a formal small real parameter (for the rigorous analysis of the bifurcations of small-amplitude modes see [33]).
Sincew n (x) are normalized, the squared L 2 -norm of the nonlinear modes U = |w n (x)| 2 dx (which is also associated with the total energy flow in the optical applications of the model) is approximately equal to ε 2 : U ≈ ε 2 for ε ≪ 1. Substituting expansions (8) into Eq. (5) and collecting the terms of the ε 3 -order, one arrives at an inhomogeneous linear differential equation with respect to the function w (3) n . The solvability condition for the latter equation yields the expression for the coefficient b For the conservative case (α = 0) the eigenfunctions w n (x) are real, and hence the coefficients b (2) n are obviously real too for any n. For α = 0, the eigenfunctions w n (x) possess nontrivial real and imaginary parts as described above. However, parities of real and imaginary parts ofw n (x) imply that for a real and even kernel func- n are also real, which shows that continuous families of nonlinear modes can exist in the PT -symmetric model with α > 0. The family with n = 0 represents the fundamental nonlinear modes, and the families n = 1, 2, . . . contain the higher-order modes.

A. Statement of the problem
Following to the the standard procedure of the linear stability analysis, we consider a perturbed nonlinear which after linearization with respect to small functions u(x) and v(x) gives the linear stability eigenvalue problem with the linear operator L given as K is the convolution operator, i.e., The nonlinear mode w n (x) is said to be linearly unstable if operator L has at least one eigenvalue ω such that Im ω < 0. Otherwise, w n (x) is linearly stable.

B. Linear stability of small-amplitude modes
Let us now analyze the spectrum of the operator L for small-amplitude nonlinear modes belonging to the nth family and situated in the vicinity of the bifurcation from the linear limit. For nonlinear modes of zero amplitude (i.e., for ε = 0) the operator L acquires the diagonal form where Since the spectrum of the operator L n +b n is known from the discussion in Sec. III, it is easy to deduce that the spectrum of the operatorL n consists of two subsets. Eigenvalues and eigenvectors of the first subset read ω n,n = 0 which remains zero and double when passing from the linear limit ε = 0 to ε > 0 [27]. The operatorL n also has 2n double nonzero eigenvalues: n,2n−k , where k runs from 0 to 2n except for k = n. Each double eigenvalue Ω n,k is semi-simple; the two corresponding linearly independent eigenvectors are given as Passing from the linear limit ε = 0 to small nonzero ε, each double eigenvalue Ω n,k generically splits into two simple eigenvalues which will be either both real or complex conjugated. If each double eigenvalue splits into two real ones, then the small-amplitude modes belonging to the nth family are stable for both attractive and repulsive nonlinearities. Vice versa, if at least one double eigenvalue gives birth to a pair of two complex-conjugated simple eigenvalues, then the small-amplitude nonlinear modes of the nth family are unstable. Notice also that for n = 0 no double eigenvalues Ω n,k exists. Therefore the fundamental solutions from family n = 0 are always stable in the small-amplitude limit.
It is easy to check that the eigenvalue Ω n,k splits in exactly the same way as the opposite double eigenvalue Ω n,2n−k = −Ω n,k , i.e., if Ω n,k splits into two real eigenvalues then Ω n,2n−k also splits into two real eigenvalues and vice versa. Therefore, in order to examine stability of small-amplitude nonlinear modes that belong to the nth family, it is sufficient to analyze only n positive double eigenvalues Ω n,k with k = 0, . . . , n − 1.
In order to examine splitting of the double eigenvalues Ω n,k , we employ Eqs. (8) which yield the following asymptotic expansion for the linear stability operator: Following the standard arguments of the perturbation theory for linear operators [30], in order to address the splitting of the double eigenvalue Ω n,k we consider a 2×2 matrix defined as a(x)dx for any two column vectors a and b (here † means the transposition and complex conjugation). If both the eigenvalues of matrix M n,k are real, then the double eigenvalue Ω n,k splits into two real and simple eigenvalues. If such a situation takes place for all k = 0, 1, . . . , n−1, then the small-amplitude nonlinear modes w n (x) from the nth family are stable. On the other hand, if for some k the matrix M n,k has a complex eigenvalue, then the double eigenvalue Ω n,k splits into a pair of complex-conjugated eigenvalues. This is sufficient to conclude that the small-amplitude nonlinear modes from the nth family are unstable.
Recall that real and imaginary parts ofw n (x) have opposite parities. Using this fact, one can establish that all the entries of the matrices M n,k are real numbers which read Thus in order to check reality of the eigenvalues of the matrix M n,k , it is sufficient to consider the discriminant of its quadratic characteristic equation, i.e., the quantity D n,k = [(M n,k ) 1,1 + (M n,k ) 2,2 ] 2 − 4 det M n,k . For D n,k > 0 the eigenvalues of M n,k are real while a situation with D n,k < 0 corresponds to a pair of complex conjugated eigenvalues. The situation D n,k = 0 requires a more delicate analysis as in this case behavior of the double eigenvalue Ω n,k is determined by the next terms of the asymptotic expansions. Using analytical expressions for the eigenfuctions w n (x), the entries of matrices M n,k and hence the determinant D n,k can be computed analytically with a computer algebra program (at least, for several relevant shapes of the kernel function K(x)). The resulting analytical expressions are too bulky to be presented in the paper. However, for an interested reader to have an idea of how the expressions for the determinants D n,k read, in Appendix A we provide two expressions for D 1,0 and D 2,0 obtained for the Gaussian kernel and α = 0. Complexity of the expressions for D n,k increases drastically if larger numbers n or nonzero values of α are involved. Nevertheless, these expressions are available analytically and can be used for finding the domains of positivity and negativity of the determinants D n,k which are considered as functions of λ and α.

V. RESULTS FOR THE REAL PARABOLIC POTENTIAL
Let us start presentation of our results with an important particular case α = 0, when PT -symmetric contribution is absent, and the model is conservative. First, it is relevant to notice that for α = 0 the eigenvalue problem (10) admits an exact solution with eigenvalue ω = 2 and the respective eigenvector given as Notice however that this exact solution does not hold for nonzero α.
Results of linear stability analysis of small-amplitude nonlinear modes from several first families for the Gaussian kernel (2) are summarized in Table I. For n = 1, 2, . . . 6 and for k = 0, 1, . . . , n − 1, Table I reports sign of the determinant D n,k in the local limit λ = 0 and the approximated values of positive zeros of determinants D n,k being considered as functions of λ [at the same time, those are the values of λ at which functions D n,k = D n,k (λ) change their signs]. From Table I one can recover information on stability or instability of the small-amplitude nonlinear modes bifurcating from the nth linear eigenstate for the given value of the nonlocality parameter λ: the stability takes place if D n,k (λ) > 0 for each k = 0, 1, . . . , n − 1.
It follows from Table I that for n = 1 and k = 0 the respective determinant D 1,0 (λ) is positive for any λ ≥ 0 (see also the explicit expression (A1) for D 1,0 in Appendix A). Therefore, the small-amplitude nonlinear modes bifurcating from the linear eigenstate n = 1 are stable for any λ ≥ 0. The small-amplitude modes with n = 2 are unstable for 0 ≤ λ 1.82 but become stable for λ 1.82. For n > 2 the stability situation might be more complex with the half-axis λ ≥ 0 divided into several alternating intervals of stability and instability. In the local We additionally mention that for each n one has D n,n−1 (λ) > 0 for all λ. This is a consequence of the exact solution (23) which implies that the spectrum of operator L always contains an eigenvalue equal to 2. Therefore, the double eigenvalue Ω n,n−1 = 2 always splits in the pair of real eigenvalues (one of which is equal to 2). As a result, the eigenvalue Ω n,n−1 can not cause instability of small-amplitude nonlinear modes.
Let us now turn to examination of other kernel functions K(x). The results of linear stability analysis for the exponential kernel (3) are reported in Table II which features an essentially different stability picture comparing to that for the Gaussian kernel. Specifically, we observe that for the exponential kernel in the limit of strong nonlocality λ → ∞, only the small-amplitude modes with n = 0, 1, 2 are stable, whereas the modes with n ≥ 3 are unstable both in the local limit λ = 0 and for all positive λ.
To summarize the intermediate results obtained so far, we have encountered different stability situations for two different shapes of the kernel K(x). In the case of the smooth Gaussian kernel, small-amplitude nonlinear modes bifurcating from the nth linear state become stable in the limit of high nonlocality for any n. On contrary, for the exponential kernel the small-amplitude modes are stable in the strongly nonlocal limit only for n ≤ 2 and unstable for n ≥ 3. The found difference in stability situation of strongly nonlocal modes reveals certain correspondence to the previous results on the stability of the high-order nonlocal nonlinear states, such as multipolemode 1D solitons [5] and 2D spatial solitons [6]. Both [5] and [6] report that stable higher-order states can be found in the case of the Gaussian kernel, while the exponential kernel results in instability of all high order solitons, starting from some threshold order. It has been also suggested in [6] that the mentioned difference between Gaussian and exponential kernels is related to the singularity which the exponential kernel has at the origin x = 0. Our approach allows to substantiate the above suggestion. Indeed, for the Gaussian kernel one can establish asymptotic behavior of entries of the matrix M n,k in the limit λ → ∞. Specifically, for λ → ∞ one has Technical details on computing the asymptotics (24)- (25) are presented in Appendix D. Additionally, in Appendices B and C we provide auxiliary relations necessary for the derivation of (24)- (25). Relations (24)-(25) immediately imply that in the limit of large λ for k = 0, 1, . . . , n − 2 the diagonal elements of the matrix M n,k dominate the off-diagonal elements, which guarantees reality of the eigenvalues of M n,k . Since the instability never occurs for k = n − 1 [due to the presence of the exact solution (23), as discussed above], we confirm that for the Gaussian kernel (or, more generally, for an infinitely differentiable kernel) the small-amplitude modes of arbitrarily high order n are stable in the limit of the strong nonlocality.
At the same time, asymptotic results (24)-(25), generically speaking, do not hold for a non-differentiable kernel. For example, for the exponential kernel we have computed that for n = 3 and k = 0 all the entries of the matrix M n,k behave as ∝ λ −2 for λ → ∞, i.e., the diagonal domination in the strongly nonlocal limit is not guaranteed. Table III which reports the stability analysis for the algebraically depending kernel (4) further confirms the stability of strongly nonlocal high-order small-amplitude modes for an infinitely differentiable kernel. Besides of this, we have also performed the same linear stability analysis for several nonphysical but still illustrative nonlocal kernels such as a rectangular kernel with finite sup- a triangular shaped kernel with finite support and a smooth exponentially decaying kernel [the particular choice of the kernels (26)-(28) is motivated by convenience of analytical evaluation of the respective convolution integrals]. Kernels (26) and (28) also feature stabilization of all higher-order strongly nonlocal modes, whereas for the triangular kernel (27) the higher-order strongly nonlocal modes are unstable. This allows to conclude that the main feature responsible for the stability is not the very presence of the singularity of the kernel, but the location of the singularity: the rectangular kernel (26) features two discontinuities at x = ±λ but its behavior is nevertheless similar to that of the Gaussian and Lorentzian. The singularity at x = 0 (the exponential and triangular kernels) does not allow for the stabilization of the higher-order small-amplitude modes.

VI. RESULTS FOR THE PT -SYMMETRIC PARABOLIC POTENTIAL
The linear stability analysis conducted above can be naturally extended on the case of nonzero PT -symmetric modulation α > 0. In this case, the stability results can be conveniently represented in the plane (α, λ). Then the axis α = 0 recovers the results of the previous section, whereas the horizontal axis λ = 0 corresponds to the case of local PT -symmetric parabolic potential studied in [27]. The origin, i.e., the point λ = α = 0 is "the most standard" case of the real potential with the local nonlinearity [20].
Similar to the previous section, we focus our attention on the Gaussian, exponential and Lorentzian kernels. According to results of our analysis, for n = 1 the smallamplitude nonlinear modes are stable for any values of λ and α from the considered range. For the families with n ≥ 2, however, the stability situation becomes more complex and the (α, λ)-diagrams turn out to be divided into the domains of stability and instability. Such stability diagrams for the families with n = 2 and n = 3 are presented in Fig. 1.
Turning to the local limit (λ = 0), one can observe that for sufficiently large α the small-amplitude modes become stable. Notice also that in the limit λ → 0 all the considered kernels tend to the delta-function δ(x) and therefore display the identical behavior. However, in accordance with results of the previous subsection, the considered kernels feature different behavior in the limit α = 0. As one can see in the panels of the first and the third rows, for the Gaussian and Loretzian kernels, in the limit α = 0 large values of λ lead to stabilization of the modes both for n = 2 and n = 3. For the exponential kernel (panels in the second row) the stabilization occurs for n = 2 only, while for n = 3 the modes remain unstable for any arbitrarily large λ. However, this situation changes drastically when a nonzero PT symmetry α > 0 is brought into consideration: say, already for α = 0.5 the nonlinear modes for the exponential kernel and n = 3 become stable for any λ (from the considered range). Overall, we can conclude that the diagrams presented in Fig. 1 feature nontrivial and fairly interesting structure. In particular, one can observe that in some cases increase of α can lead either to stabilization or destabilization of the modes. The same remains true if one considers increase or decrease of λ. The stability diagrams allow to conjecture that for any given λ the small-amplitude modes eventually become stable in the limit α → ∞. However, in order to verify this conjecture, an additional study must be undertaken.
Notice also that for the case α > 0 allows for instability caused by the eigenstate with k = n − 1 (i.e., k = 1 for n = 2 and k = 2 for n = 3). No such instabilities is possible in the case α = 0 due to the exact solution (23).

VII. NONLINEAR MODES OF FINITE AMPLITUDE
The stability results of the previous sections apply to nonlinear modes of infinitesimally small amplitude. One can check persistence of the stability results by constructing numerically the families of nonlinear modes of finite amplitude and computing the eigenvalues of the corresponding linear stability problem. The families of nonlinear modes can be visualized in the form of continuous curves on the plane σU vs. b where U = |w n (x)| 2 dx is the total energy flow. Then the small-amplitude modes are situated in the vicinity of the points b n = −(2n + 1), σU = 0.
Let us first recall the results related to the local case (λ = 0). The corresponding diagrams for the conservative case [20] and for the PT -symmetric case [27] are presented in the two upper panels of Fig. 2. The panels show several families that bifurcate from the lowest linear eigenvalues n = 0, 1, 2, 3. In the conservative case, only the two lowest families (n = 0, 1) are stable in FIG. 1: Linear limit stability diagrams for n = 2 and n = 3 in the plane (α, λ). Darker domains correspond to instability. Labels "k = 0", "k = 1", etc. at the unstable domains show the number of the double eigenvalue Ω n,k responsible for the instability.
the small-amplitude limit, whereas the small-amplitude modes with n = 2, 3 are unstable, In the PT -symmetric case, all the four shown families are stable in the linear limit, but at least three of the four families loose stability as the strength of nonlinearity becomes sufficiently large.
The diagrams for the nonlocal case λ = 4 (with the Gaussian kernel) are shown in the lower panels of Fig. 2. According to Table I and Fig. 1, at the chosen value of the nonlocality parameter, the nonlinear modes with n = 0, 1, 2, 3 are stable in the small-amplitude limit. The numerical results in Fig. 2 indicate that the stability persists at least for small and moderate strengths of the nonlinearity (|σU | 40). Notice also that in the panel with (α = 1, λ = 4) the nonlocality cancels the merging between the subsequent families which can be observed in the panel with (α = 1, λ = 0). We however emphasize that an accurate description of the limit of strong nonlinearity |σU | → ∞ should be elaborated by means of a proper asymptotic technique. While in the local case such studies have been already initiated [34,35], analysis of the nonlocal PT -symmetric case in the limit of strong nonlinearity remains the open issue at the moment.

VIII. CONCLUSION
In this study, we have demonstrated that the combined effect of the parabolic potential, nonlocality and PT symmetry results in a nontrivial stability picture for high-order small-amplitude nonlinear modes. The specific results of our work can be outlined as follows: • A set of continuous families of nonlinear modes exists in the nonlocal PT -symmetric nonlinear Schrödinger equation. The spectrum of the corresponding linear eigenvalue problem is available in the analytical form, and thus the stability problem for the small-amplitude nonlinear modes can be reduced to the searching real roots of certain analytical expressions.
• In the conservative case, small-amplitude nonlinear modes of arbitrary order become stable for sufficiently strong nonlocality, provided that the kernel is a sufficiently smooth function in the vicinity of the origin. If the kernel feature a singularity at x = 0, no stabilization of high-order nonlinear modes is observed.
• Both the degree of nonlocality and the strength of the PT symmetry can be used to manage stability of the small-amplitude modes.
• The above stability conclusions remain valid for nonlinear modes of finite amplitude.
Then the Plancherel theorem says that where g 1,2 (s) = F {f 1,2 }; and from the convolution theorem one has where Another well-known relation expresses the pth central moment of a function via the Fourier transform (p = 0, 1, . . .): .
The terms decaying as λ −1 cancel each other.
(D7) From definition of functions F n,k (s) and F n,2n−k (s) it follows that in vicinity of s = 0 they behave as F n,k (s) = i n−k (n − k)! n! k! 2 n−k s n−k + o(s n−k ) s→0 , F 2n−k,k (s) = i n−k (n − k)! (2n − k)! n! 2 n−k s n−k + o(s n−k ) s→0 .