Analytical construction of soliton families in one- and two-dimensional nonlinear Schrödinger equations with non-parity-time-symmetric complex potentials
aa r X i v : . [ n li n . PS ] F e b Analytical construction of soliton families in one- and two-dimensional nonlinearSchr¨odinger equations with non-parity-time-symmetric complex potentials
Jianke Yang
Department of Mathematics and Statistics, University of Vermont, Burlington, VT 05405, U.S.A
The existence of soliton families in non-parity-time-symmetric complex potentials remains poorlyunderstood, especially in two spatial dimensions. In this article, we analytically investigate thebifurcation of soliton families from linear modes in one- and two-dimensional nonlinear Schr¨odingerequations with localized Wadati-type non-parity-time-symmetric complex potentials. By utilizingthe conservation law of the underlying non-Hamiltonian wave system, we convert the complex solitonequation into a new real system. For this new real system, we perturbatively construct a continuousfamily of low-amplitude solitons bifurcating from a linear eigenmode to all orders of the small solitonamplitude. Hence, the emergence of soliton families in these non-parity-time-symmetric complexpotentials is analytically explained. We also compare these analytically constructed soliton solutionswith high-accuracy numerical solutions in both one and two dimensions, and the asymptotic accuracyof these perturbation solutions is confirmed.
1. INTRODUCTION
Nonlinear wave phenomena in parity-time ( PT ) symmetric systems have been under intensive studies in the pastdecade (see [1–3] for reviews). Although the concept of PT symmetry originated from non-Hermitian quantum me-chanics [4–6], it was the interpretation of PT symmetry as balanced gain and loss that made it flourish in optics andmany other branches of physics [1–3, 6, 7]. PT symmetric systems are important for at least two reasons. From theintellectual point of view, these systems are the first reported non-Hamiltonian systems that, despite the gain andloss, exhibit many properties of Hamiltonian systems — such as all-real linear spectra and continuous families ofsolitons [1–6]. From the practical point of view, PT symmetry has inspired many interesting applications, such asthe coherent perfect absorber laser [8–10] and single-mode PT lasers [11, 12]. While applications of PT symmetry arestill developing, its peculiar Hamiltonian-like phenomena, such as the existence of all-real spectrum and continuousfamilies of solitons, have already been understood from a mathematical point of view [1, 2, 5, 6]. In particular, thisunderstanding relies entirely on the PT symmetry.In the past few years, it was discovered that certain non- PT -symmetric non-Hamiltonian systems also share prop-erties of Hamiltonian systems. For example, the linear Schr¨odinger operator with certain types of non- PT -symmetriccomplex potentials could still admit all-real spectra [13–16]. In addition, the one- and two-dimensional NLS equationswith Wadati-type non- PT -symmetric complex potentials could still admit continuous families of solitons [15, 17–19]. Furthermore, in the NLS equations with Wadati-type non- PT -symmetric potentials, the linear-stability spectraof solitons still exhibit the quartet eigenvalue symmetry that is typical of Hamiltonian systems [20]. In a genericnon- PT -symmetric non-Hamiltonian system, none of these properties would hold. Thus, why these Hamiltonian-likeproperties arise in certain types of non- PT -symmetric non-Hamiltonian systems is an intriguing theoretical ques-tion. While the all-real spectra of certain non- PT -symmetric complex potentials can be explained by techniques suchas supersymmetry and pseudo-Hermiticity [13–16], analytical explanations for the other properties associated withnonlinear non- PT -symmetric systems remain elusive.This article is concerned with the question of why the NLS equations with Wadati-type non- PT -symmetric com-plex potentials could still admit continuous families of solitons. This phenomenon is peculiar, since these non- PT -symmetric systems are non-Hamiltonian due to the presence of gain and loss, and solitons in non-Hamiltonian systemsare generically isolated and do not exist as continuous families due to the double balancing requirement of nonlinearitywith dispersion and gain with loss [21]. Numerical evidence to support this generic behavior in a non- PT -symmetricsystem can be found in [22], and a more mathematical reason for it can be found in [23]. In view of this genericbehavior and in the absence of PT symmetry, why soliton families could appear in the NLS equations with specialWadati-type non- PT -symmetric potentials is a deep mathematical mystery. It is physically meaningful for us to addthat, unlike PT -symmetric potentials where the spatial gain and loss distributions must be balanced in an exactanti-symmetric way, the Wadati potentials allow the gain and loss distributions to be arbitrary, which could poten-tially accommodate more realistic non-Hamiltonian physical systems in optics and beyond. A physical setup to realizeWadati potentials in a coherent atomic system has been proposed in [24].In the one-dimensional (1D) case, some analytical understanding on this question has been provided in [17, 19].In [17], Konotop and Zezyulin discovered a constant of motion for the underlying soliton equation with Wadatipotentials. Combining this constant of motion with a shooting argument, the authors gave a plausible, but notdefinitive, explanation for these soliton families. In [19], the authors used this constant of motion to convert theoriginal second-order complex soliton equation into a second-order real equation for the amplitude of the soliton.From this real soliton-amplitude equation, it was shown that continuous families of solitons bifurcating from linearmodes could be constructed perturbatively. One drawback of this treatment in [19] is that this real amplitude equationhas some sign ambiguity in front of a square root term, which can cause technical complications. Another drawback,which is more serious, is that this treatment cannot be generalized to two and higher spatial dimensions.In the 2D case, while soliton families in the 2D NLS equations with separable Wadati-type non- PT -symmetricpotentials were briefly mentioned on numerical grounds in [18], there has been absolutely no analytical explanationfor this phenomenon yet, except that a conservation law for the underlying non-Hamiltonian 2D equation was reportedin that same article. Note that in this 2D case, the shooting argument of [17] no longer applies. In addition, the real-amplitude-equation treatment of [19] also fails. Thus, new approaches need to be developed to analytically explainthese 2D soliton families.We would like to mention that continuous families of solitons in the 1D NLS equation perturbed by non- PT -symmetric potentials more general than the Wadati-type were also reported by Kominis et al. [25] through Melnikov’sperturbation method. Since the authors’ analysis was carried out only to the first order of the perturbation series,we suspect that those soliton families in non-Wadati potentials are valid only to the first order of the perturbationtheory, but not to higher orders. If so, then those “soliton families” would be just approximate solutions, but nottrue solitons. This suspicion makes it more imperative to analytically explain the existence of soliton families innon- PT -symmetric Wadati potentials, since such analytical understanding could shed light on the nature of “solitonfamilies” reported in [25] for non-Wadati potentials.In this article, we analytically investigate the bifurcation of soliton families from linear modes in the 1D and2D NLS equations with non- PT -symmetric Wadati-type localized potentials through a new perturbative treatment.Utilizing the constant of motion of the underlying soliton equation, we convert this complex soliton equation into anew real system. The advantage of this new real system is that it allows us to analytically construct low-amplitudesoliton families perturbatively to all orders of the amplitude in both one and two dimensions. Hence, soliton familiesin these 1D and 2D non- PT -symmetric systems are analytically established. The reason this construction can bepursued to all orders is that the linear operator of these perturbation equations possesses two localized functions inits kernel, while the associated adjoint operator contains a single localized or bounded function in its kernel. Thesekernel structures, together with the phase invariance of solitons, ensure that at each order, the Fredholm condition forlocalized perturbation solutions can always be satisfied. Hence, we can construct a low-amplitude soliton solution, as aperturbation series to all orders, at each propagation constant in a continuous interval bordering the linear eigenmodeof the potential. In other words, a soliton family bifurcating from a linear mode is derived in the underlying non- PT -symmetric non-Hamiltonian system. These analytically constructed perturbation-series solutions for the solitonfamilies are also compared to direct numerical solutions, and the asymptotic accuracy of these perturbation seriessolutions is confirmed.
2. CONSTRUCTION OF SOLITON FAMILIES IN THE 1D CASE
We first consider the 1D NLS equation iU t + U xx + V ( x ) U + σ | U | U = 0 (1)with a non- PT -symmetric Wadati potential V ( x ) = g ( x ) + i g ′ ( x ) , (2)where g ( x ) is an asymmetric real function that is differentiable everywhere, the prime represents differentiation, and σ = ± g ( x ) is real and asymmetric, V ∗ ( − x ) = V ( x ), i.e., the complex potential V ( x ) is non- PT -symmetric [1–3, 7]. Potentials of this form appeared in Wadati’s investigation of complex potentialswith real spectra [26], and are thus sometimes referred to as the Wadati potentials in the literature. In the opticalcontext, the complex potential V ( x ) in Eq. (1) corresponds to the complex refractive index of the medium, wherethe imaginary part of V ( x ), i.e., Im( V ), describes the spatial gain and loss distributions, with regions of Im( V ) > V ) < g ( x ) in theWadati potential (2) can be arbitrary, this complex potential then can accommodate optical systems with arbitrarygain and loss distributions. The main constraint of the Wadati potential is that, the real refractive index profile ofthe medium, as described by the real part of the complex potential Re( V ), should be designed accordingly as g ( x ).But this requirement on the real refractive index profile can be readily met given the sophisticated refractive-indexengineering technology that is currently widely available.An important property of the NLS equation (1) with Wadati potentials is that, although this equation is non-Hamiltonian due to the complex potential, it admits a conservation law Q t + J x = 0 , (3)where Q = − U ∗ (i U x − gU ) , J = | U x + i gU | + i U ∗ U t + σ | U | , (4)and the asterisk ‘*’ represents complex conjugation. This conservation law is a special case of the more generalconservation law reported in [18] for the 2D NLS equation with a separable Wadati-type potential.Solitons in Eq. (1) are of the form U ( x, t ) = e i µt u ( x ) , (5)where µ is a real propagation constant, and u ( x ) is a localized function satisfying the soliton equation u xx + ( g + i g ′ ) u − µu + σ | u | u = 0 . (6)Notice that this complex soliton equation is phase-invariant, i.e., if u ( x ) is a solution, so is e i θ u ( x ), where θ is anarbitrary real constant. Substituting the soliton solution (5) into the conservation law (3), we get dJ/dx = 0, where J ( x ) = | u x + i gu | − µ | u | + σ | u | . (7)Since solitons decay to zero as x → ±∞ , we see that J ( x ) = 0, which is a constant of motion for the soliton equation(6). This constant of motion is equivalent to the one reported in Ref. [17] for the same equation (6).Soliton families in Eq. (6), parameterized by the propagation constant µ , for non- PT -symmetric Wadati potentialswere reported numerically in [15], and studied analytically in [17, 19, 23] with limited success. In particular, theperturbative construction of small-amplitude soliton families as proposed in [19, 23] exhibits some difficulties. Theperturbative construction in [23] was based on the complex soliton equation (6). The difficulty with this construction,as explained in [23], is that each order of the perturbation series creates a nontrivial condition which needs to besatisfied, and it is almost impossible to prove that all those infinite number of conditions would hold. The perturbativeconstruction in [19] was based on a real second-order equation for the amplitude | u ( x ) | of the soliton, and this realamplitude equation was derived from the original complex equation (6) with the help of the above constant of motion J ( x ) = 0. This latter construction removed those infinite number of nontrivial conditions of the former, and thusmade the perturbative construction possible, at least in principle. But it does create some technical difficulties.For example, this reduced amplitude equation contains a square root term, whose sign can be ambiguous and causetechnical complications. To remove this ambiguity, some technical assumptions had to be imposed in [19]. A moreserious problem with this latter treatment is that it does not work for the 2D case. In other words, in two (and higher)spatial dimensions, we will not be able to convert the original complex soliton equation into a single real equation forthe amplitude of the soliton.In this section, we will develop a new perturbative construction of low-amplitude soliton families in Eq. (6), whichcan be easily pursued to all orders of the perturbation series. More importantly, this new 1D treatment can be readilygeneralized to the 2D case.For the technical convenience of our perturbative construction, we will assume that the Wadati potential (2) islocalized in space, i.e., the real function g ( x ) in this potential will be assumed to be localized. This assumption oflocality on the potential has two main benefits. One is that such a Wadati potential often admits a discrete realeigenvalue [15, 16], which is the starting point of our perturbative calculation. The other is that under this localityassumption, the eigenfunction associated with this discrete real eigenvalue of the potential features simple and explicitexponential decay at large distances. These explicit decay rates of the eigenfunction facilitate our derivation andunderstanding of the kernels for the linearization operator and its adjoint in the upcoming section 2 2.2. If the Wadatipotential (2) is not localized (for instance, unbounded) but still admits a discrete real eigenvalue, then the analysisof this section can still go through, because the kernel structures of the linearization operator and its adjoint to beestablished in Sec. 2 2.2 would still remain valid. However, if the Wadati potential (2) is periodic, then the situationwould be different. In this case, the periodic potential does not admit any discrete real eigenvalues. Instead, thespectrum of the potential comprises Bloch bands. Low-amplitude solitons, if any, would have to bifurcate out fromedges of these Bloch bands as envelope solitons [27]. The analytical calculation of soliton bifurcation from Bloch-bandedges in a periodic Wadati potential would be very different from the one to be developed in this section, and it willbe left for future studies. It can be checked that the original complex soliton equation (6) is equivalent to two real equations — one is thatthe real part of (6) is zero, and the other is J = 0, where J is given in Eq. (7). The first real equation comes directlyfrom (6), and the second one is the constant of motion discussed below Eq. (7). To see these two real equationscombined could also reproduce the original complex equation (6), we only need to notice that dJ/dx is equal to thereal part of the product between u ∗ x − igu ∗ and the left side of the complex soliton equation (6). Thus, if J = 0 andthe real part of (6) is zero, then the imaginary part of (6) needs to be zero as well.Expressing u ( x ) as u ( x ) = p ( x ) + i q ( x ) , (8)where p ( x ) and q ( x ) are the real and imaginary parts of the complex function u ( x ), these two real equations forsolitons are p xx + ( g − µ ) p − g ′ q + σ ( p + q ) p = 0 , (9)( p x − gq ) + ( q x + gp ) − µ ( p + q ) + σ p + q ) = 0 . (10)This system of two real equations will be the one we use to analytically calculate soliton families. It is important tonotice that this is a third-order real system, which contrasts the original soliton equation (6), which is a fourth-orderreal system when that complex equation is split into two real second-order equations for p ( x ) and q ( x ). This third-order real system also contrasts the second-order real system we derived in Ref. [19] for the amplitude | u ( x ) | of thesoliton.Now, we perturbatively construct a continuous family of low-amplitude solitons bifurcating from a linear discreteeigenmode of a localized Wadati potential. Suppose the Schr¨odinger operator ∂ xx + V ( x ) with a localized Wadatipotential (2) admits a discrete real eigenvalue µ , whose eigenfunction is φ ( x ) + i ψ ( x ), where φ ( x ) and ψ ( x ) arelocalized real functions. Then, (cid:0) ∂ xx + g + i g ′ (cid:1) ( φ + i ψ ) = µ ( φ + i ψ ) . (11)The existence of such a real eigenvalue is common in a Wadati potential. For instance, it was shown in [15] that if g ( x )is a single-humped localized real function, then the spectrum of the corresponding Wadati potential is strictly real. Inthe more general case, it was shown in [16] that eigenvalues in a Wadati potential always come as complex-conjugatepairs and are thus often real. Because this potential is assumed to be localized, its discrete real eigenvalue µ mustbe positive, i.e., µ > µ near µ , and this soliton can be expanded into the following perturbation series, p ( x ; µ ) = ǫ / (cid:2) p ( x ) + ǫp ( x ) + ǫ p ( x ) + · · · (cid:3) , (12) q ( x ; µ ) = ǫ / (cid:2) q ( x ) + ǫq ( x ) + ǫ q ( x ) + · · · (cid:3) , (13)where ǫ = µ − µ and is assumed to be small positive (so that ǫ / is real). This means that we assume the bifurcationis to the right side of µ , i.e., µ > µ . As we will see in later text [see Eq. (34)], this rightward bifurcation can beinduced by a proper choice on the sign of nonlinearity σ . If this sign of nonlinearity is opposite of that choice, thesoliton bifurcation will be to the left side of µ . In that case, we can define ǫ = µ − µ , and the rest of the perturbativecalculation would be very similar.Substituting the above perturbation expansion into the real system (9)-(10), we get a sequence of real equationsfor the functions ( p k , q k ). The equations for ( p , q ) are( ∂ xx + g − µ ) p − g ′ q = 0 , (14)( p x − gq ) + ( q x + gp ) − µ ( p + q ) = 0 . (15)Even though this is a nonlinear system, it is scaling invariant, i.e., if ( p , q ) is a solution, so is ( αp , αq ), where α isan arbitrary real constant. Thus, this system is actually an eigenvalue problem in disguise and is equivalent to thelinear complex eigenvalue problem (11). Its solution then is (cid:20) p q (cid:21) = c (cid:20) φψ (cid:21) , (16)where c is a real constant to be determined. Indeed, since φ + i ψ is a solution to the linear eigenvalue problem (11),the above ( p , q ) then satisfy the original equations (9)-(10) to leading order, which are Eqs. (14)-(15).Utilizing the above ( p , q ) solution, we find that the functions ( p k , q k ) for k ≥ L (cid:20) p k q k (cid:21) = (cid:20) f k g k (cid:21) , (17)where L = (cid:20) ∂ xx + g − µ − g ′ ( φ ′ − gψ ) ∂ x + g ( ψ ′ + gφ ) − µ φ ( ψ ′ + gφ ) ∂ x − g ( φ ′ − gψ ) − µ ψ (cid:21) , (18) (cid:20) f g (cid:21) = c (cid:20) φ − σc ( φ + ψ ) φ ( φ + ψ ) − σc ( φ + ψ ) (cid:21) , (19) (cid:20) f g (cid:21) = (cid:20) (cid:0) − σp − σq (cid:1) p − σp q q c (cid:2) (cid:0) − σp − σq (cid:1) ( p p + q q ) + µ ( p + q ) − ( p x − gq ) − ( q x + gp ) ) (cid:3) (cid:21) , (20) (cid:20) f k g k (cid:21) = (cid:20) M M M M (cid:21) (cid:20) p k − q k − (cid:21) + " N [1] k N [2] k , k ≥ , (21)the matrix elements M ij are k -independent and given by the formulae M = 1 − σp − σq , M = − σp q , M = 1 c (cid:2) p (cid:0) − σp − σq (cid:1) + µ p − ( p x − gq ) ∂ x − g ( q x + gp ) (cid:3) , M = 1 c (cid:2) q (cid:0) − σp − σq (cid:1) + µ q + g ( p x − gq ) − ( q x + gp ) ∂ x (cid:3) , and N [1] k , N [2] k are functions which depend only on k , p , p , . . . , p k − , q , q , . . . , q k − and g ( x ). For example, when k = 3, N [1]3 = − σ (cid:0) p p + 2 p q q + p q (cid:1) , N [2]3 = 12 c (cid:2) ( p + q ) (cid:0) − σp − σq (cid:1) − σ ( p p + q q ) (cid:3) . Next, we will show that we can solve the linear nonhomogeneous equations (17) and obtain localized solutions( p k , q k ) for all k , using the Fredholm alternative method. The key to solving linear nonhomogeneous equations (17) by the Fredholm alternative method is to understand thekernel structures of the linear operator L and its adjoint operator L A . Under the inner product of h F, G i ≡ Z ∞−∞ [ F ( x )] T G ( x ) d x, (22)where the superscript ‘ T ’ represents the transpose of a vector or matrix, the adjoint operator of L is L A = (cid:20) ∂ xx + g − µ − ∂ x ( φ ′ − gψ ) + g ( ψ ′ + gφ ) − µ φ − g ′ − ∂ x ( ψ ′ + gφ ) − g ( φ ′ − gψ ) − µ ψ (cid:21) . (23)First, we consider the kernel structure of operator L . It is easy to check that this kernel contains the following twolocalized functions K ≡ (cid:20) φψ (cid:21) , K ≡ (cid:20) − ψφ (cid:21) , (24)where L K = L K = 0 . (25)Indeed, L K = 0 is equivalent to the complex linear eigenvalue equation (11), and L K = 0 is equivalent to thiscomplex eigenvalue equation with the eigenfunction changing from φ + i ψ to i( φ + i ψ ), which clearly remains aneigenfunction. Another way to understand these kernel functions is that, the first kernel function K is induced bythe scaling invariance of the complex linear eigenvalue equation (11), and the second kernel function K is inducedby the phase invariance of that same equation.It is clear that L is a third-order differential operator. Thus, the system L K = 0 admits one more linearlyindependent solution K in addition to K and K . This third solution is obviously unbounded in space. Indeed,since φ + i ψ is the eigenfunction of the Schr¨odinger operator with a localized potential at the positive eigenvalue µ [see Eq. (11)], both φ ( x ) and ψ ( x ) decay exponentially at the rate of e −√ µ | x | when x → ±∞ . Then, converting thesystem L K = 0 into a system of three first-order equations and using Abel’s formula, we can show that this thirdsolution K ( x ) grows exponentially at the rate of e √ µ | x | when x → ±∞ .Next, we consider the kernel structure of L A . Functions in this kernel can be derived from the functions in thekernel of L . One way to do so is to first rewrite the equation L K = 0 with K ≡ [ K [1] , K [2] ] T as a first-order system ∂ x Y = P ( x ) Y (26)for Y = [ K [1] , K [1] x , K [2] ] T , where P ( x ) is a 3 × Y ( x ) of this first-orderhomogeneous system is given through the three solutions K , K and K of the original system L K = 0 as Y = K [1]1 K [1]2 K [1]3 K [1]1 ,x K [1]2 ,x K [1]3 ,x K [2]1 K [2]2 K [2]3 . (27)The adjoint of the first-order system (26) is − ∂ x Y A = P T Y A , (28)whose fundamental matrix is Y A = ( Y − ) T . Using the large- x asymptotics of the ( K , K , K ) solutions described inthe previous paragraph, together with their Wronskian expression from Abel’s formula, we can readily show that thethird column of Y A is localized with its second component decaying at the rate of e −√ µ | x | at large | x | , while the firstand second columns of Y A are unbounded with their second components growing at the rate of e √ µ | x | at large | x | .The adjoint first-order system (28) has a simple connection with the original adjoint system L A K A = 0. Specifically,if Y A = [ Y A [1] , Y A [2] , Y A [3] ] T , then K A = [ Y A [2] , Y A [3] / ( ψ ′ + gφ )] T . Using this connection, we see that the kernel of L A contains a single localized function, which we denote as K A = (cid:20) φ A ψ A (cid:21) , (29)where L A K A = 0. This K A is obtained from the third column of Y A ; so φ A ( x ) decays at the rate of e −√ µ | x | when x → ±∞ . Regarding the decay rate of ψ A ( x ), using dominant balance on the second equation of the adjoint system L A K A = 0, we can show that ψ A ( x ) decays at the same rate of g ( x ) for large | x | . The other two functions in thekernel of L A are obtained from the first and second columns of Y A and are thus both unbounded. More specifically,their first components grow at the rate of e √ µ | x | , and their second components grow at the rate of e √ µ | x | , when x → ±∞ . Utilizing the above kernel structures of operators L and L A , we can solve the linear nonhomogeneous equations(17) and obtain a localized solution ( p k , q k ) for all k . To do so, we will use the Fredholm solvability condition, whichwill be explained in this subsection.First, we notice that f k on the right side of the nonhomogeneous equations (17) is localized, and its decay rateat large | x | is e −√ µ | x | , multiplied by a certain polynomial function of x . In addition, g k on the right side of theseequations is also localized, and its decay rate at large | x | is e − √ µ | x | , multiplied by another polynomial function of x . The reason for these decay rates of ( f k , g k ) is that p n and q n in the expressions of f k and g k decay at the rate of e −√ µ | x | , multiplied by a polynomial function of x . These decay rates of ( p n , q n ) can be seen from the ǫ expansions(12)-(13) of solitons ( p, q ), which decay at the rate of e −√ µ + ǫ | x | at large | x | . These decay rates of ( p n , q n ) can alsobe seen from the equations (17) which determine them.In view of the decay rates of ( f k , g k ) on the right side of the linear nonhomogeneous equations (17), as well asthe kernel structures of linear operators L and L A delineated in the previous subsection, the Fredholm alternativetheorem says that these nonhomogeneous equations (17) would admit a localized solution ( p k , q k ) if and only if thenonhomogeneous term ( f k , g k ) T is orthogonal to the localized function K A in the kernel of L A , i.e., (cid:28)(cid:20) φ A ψ A (cid:21) , (cid:20) f k g k (cid:21)(cid:29) = Z ∞−∞ (cid:0) φ A f k + ψ A g k (cid:1) d x = 0 . (30)The Fredholm alternative theorem was originally developed for compact operators ([28], page 160), which is restrictive.But this theorem can be generalized to operators with closed range ([28], page 46). In this article, we will not attemptto prove that our operator L has closed range. Instead, we will provide an elementary proof of this Fredholm alternativeresult below.The necessity of the above condition (30) for Eq. (17) to admit a localized solution can be derived quickly by takingthe inner product of this equation with the localized function K A in the kernel of L A . To prove the sufficiency of thiscondition, we can first rewrite Eq. (17) as a first-order system ∂ x Y − P ( x ) Y = F, (31)where Y = [ p k , p k,x , q k ] T , F = [0 , f k , g k / ( ψ ′ + gφ )] T , and P ( x ) is the 3 × Y for the first-order homogeneous system of (31) has been discussed before. Using thisfundamental matrix and variation of parameters, we can derive the general solution to the nonhomogeneous system(31) as Y ( x ) = Y ( x ) (cid:18) c + Z x (cid:2) Y A ( z ) (cid:3) T F ( z ) d z (cid:19) , (32)where c is a constant vector, and Y A = ( Y − ) T is the fundamental matrix of the first-order adjoint system (28). Inview of this explicit solution formula for Eq. (31), as well as the large- x asymptotics of fundamental matrices Y and Y A described earlier, we can readily see that a localized solution Y ( x ) can be obtained, through a proper choice ofthe third element of the c constant, if the following condition is met, Z ∞−∞ (cid:2) Y A ( x ) (cid:3) T F ( x ) d x = 0 , (33)where Y A is the third column of Y A . This third column is connected to the localized function K A through a relationexplained in the last paragraph of the previous subsection. Then, using the expression of F given above, the abovecondition (33) reduces exactly to the Fredholm solvability condition (30). Thus, the sufficiency of this Fredholmcondition to guarantee the existence of a localized solution in Eq. (17) is directly proved. Now, we use the Fredholm solvability condition (30) to determine a soliton solution u ( x ; µ ) through the perturbationseries (12)-(13), to all orders of ǫ ≡ µ − µ , at each µ value near µ . These solutions then constitute a continuousfamily of solitons, parameterized by the propagation constant µ , in the non- PT -symmetric Wadati potential (2).We first consider Eq. (17) for ( p , q ). Substituting the ( f , g ) expressions (19) into the Fredholm solvabilitycondition (30) and simplifying, we see that Eq. (17) admits a localized solution ( p , q ) if and only if the constant c is selected as c = ± vuut R ∞−∞ (cid:2) φφ A + ( φ + ψ ) ψ A (cid:3) d xσ R ∞−∞ (cid:2) ( φ + ψ ) φφ A + ( φ + ψ ) ψ A (cid:3) d x . (34)In order for the quantity under the square root above to be positive, σ must have the same sign as the ratio of thetwo integrals in the above formula. In other words, in order for the soliton bifurcation to appear for µ > µ , thenonlinearity must be of a certain sign. In this case, c has two value choices which differ by a sign. But it is easy tosee that these two sign choices in c would simply lead to two soliton solutions u ( x ) = p ( x ) + i q ( x ) which also differonly by a sign. Since the u ( x ) equation (6) is phase-invariant, solutions differing by a sign are equivalent. Thus, wewill just take the plus sign for c below.When c is selected from the above formula (34), Eq. (17) admits a localized solution for ( p , q ), which we denoteas ( p s , q s ). However, since the kernel of the homogeneous operator L in Eq. (17) contains two localized functions K and K given in Eq. (24), the general localized solution ( p , q ) to the linear nonhomogeneous equations (17) isthen (cid:20) p q (cid:21) = (cid:20) p s q s (cid:21) + c (cid:20) φψ (cid:21) + d (cid:20) − ψφ (cid:21) , (35)where c and d are two real constants.It is important to recognize that the d term above can be removed by phase invariance of the complex solitonsolution u ( x ). To see this more clearly, we put the above perturbation solutions together and get u ( x ) = ǫ / (cid:2) c φ + ǫ ( p s + c φ − d ψ ) + i [c ψ + ǫ (q + c ψ + d φ )] + O( ǫ ) (cid:3) = ǫ / (cid:2) ( c + ǫc + i ǫ d )( φ + i ψ ) + ǫ (p + iq ) + O( ǫ ) (cid:3) = ǫ / e i ǫ d / c (cid:2) ( c + ǫc )( φ + i ψ ) + ǫ (p + iq ) + O( ǫ ) (cid:3) . Notice that the d term only contributes a constant phase of order ǫ to the soliton solution u ( x ). But u ( x ) isphase-invariant. Thus, that d term in (35) can be dropped and we can set (cid:20) p q (cid:21) = (cid:20) p s q s (cid:21) + c (cid:20) φψ (cid:21) (36)without loss of generality.The ( p , q ) solution in the above equation contains an unknown real constant c . This c constant will be determinedfrom the Fredholm solvability condition on the ( p , q ) equations. The equations for ( p , q ) are (17), where thenonhomogeneous terms ( f , g ) are given in Eq. (20). Substituting the ( p , q ) solutions (16) and ( p , q ) solutions(36) into the ( f , g ) expressions (20) and recalling that ( φ, ψ ) satisfy the equation (15), we find that the ( f , g )expressions (20) reduce to (cid:20) f g (cid:21) = (cid:20) f a g a (cid:21) + c (cid:20) f b g b (cid:21) , (37)where f a = (1 − σp − σq ) p s − σp q q s ,f b = (1 − σp − σq ) φ − σp q ψ,g a = 12 c (cid:2) − σp − σq )( p p s + q q s ) + µ ( p s + q s ) − ( p s,x − gq s ) − ( q s,x + gp s ) (cid:3) ,g b = 1 c (cid:2) (1 − σp − σq )( p φ + q ψ ) + µ ( p s φ + q s ψ ) − ( p s,x − gq s )( φ x − gψ ) − ( q s,x + gp s )( ψ x + gφ )] , which are independent of the unknown constant c . Then, the Fredholm solvability condition (30) at k = 2 gives theformula for the constant c as c = − R ∞−∞ (cid:0) φ A f a + ψ A g a (cid:1) d x R ∞−∞ ( φ A f b + ψ A g b ) d x . (38)The rest of the perturbation calculations can then proceed to all orders as follows. When c k − ( k ≥
2) has beenobtained, the ( p k − , q k − ) solutions are completely determined. Meanwhile, the solvability condition (30) for ( p k , q k ) isalso satisfied, and thus there exists a localized solution which we denote as ( p k,s , q k,s ). The general localized solutionsfor ( p k , q k ) can be written as (cid:20) p k q k (cid:21) = (cid:20) p k,s q k,s (cid:21) + c k (cid:20) φψ (cid:21) . (39)The constant c k will be determined from the solvability condition for the ( p k +1 , q k +1 ) equations (17). Specifically,when the above ( p k , q k ) solutions are inserted into the ( f k +1 , g k +1 ) formulae (21), it is easy to see that the solvabilitycondition (30) at k + 1 is a linear equation for c k , which we can easily solve to obtain the value of c k as c k = − *(cid:20) φ A ψ A (cid:21) , (cid:20) M M M M (cid:21) (cid:20) p k,s q k,s (cid:21) + " N [1] k +1 N [2] k +1 φ A ψ A (cid:21) , (cid:20) M M M M (cid:21) (cid:20) φψ (cid:21)(cid:29) , k ≥ . Utilizing the ( p , q ) formula (36) and the fact that ( φ, ψ ) satisfy Eq. (15), we can verify that the denominator in this c k formula is equal to the denominator in the c formula (38). Thus, the above c k formula can be reduced to c k = − *(cid:20) φ A ψ A (cid:21) , (cid:20) M M M M (cid:21) (cid:20) p k,s q k,s (cid:21) + " N [1] k +1 N [2] k +1 ∞−∞ ( φ A f b + ψ A g b ) d x , k ≥ . (40)This process is then repeated to higher orders.The only conditions for the above perturbation calculations to succeed to all orders are that the numerator anddenominator in the c formula (34), as well as the denominator in the c formula (38), are all nonzero. Thus, we onlyhave 3 numbers to check, which can be easily done for each given equation (1) when its Wadati potential V ( x ) isspecified. In this subsection, we compare the above perturbation-series soliton solution (12)-(13) with the high-accuracynumerical solution, for a continuous range of small ǫ values, and confirm the asymptotic accuracy of this analyticalsolution.In our comparison, we choose the non- PT -symmetric Wadati potential (2) as the one with g ( x ) = 0 . x + 2) + 1 . x − . (41)The resulting Wadati potential is shown in Fig. 1(a). This potential admits a discrete real eigenvalue µ ≈ . φ ( x )+i ψ (x) is plotted in Fig. 1(b). Numerically, we find that the adjoint operator L A in Eq. (23) indeed admits a single localized function ( φ A , ψ A ) T in its kernel, and this function is displayed in Fig.1(c). Utilizing these eigenfunctions and adjoint eigenfunctions, the ratio of integrals under the square root in Eq. (34)is found to be positive. Thus, according to our perturbation theory, a continuous family of solitons would bifurcateout for µ > µ under the positive sign of nonlinearity σ = 1 and for µ < µ under the negative sign of nonlinearity σ = − σ = 1, this soliton at µ = µ + 0 . P ( µ ) = Z ∞−∞ | u ( x ; µ ) | d x, (42)is shown in Fig. 1(e). These solitons are computed numerically by the Newton-conjugate-gradient method describedin [18], and their numerical error is below 10 − . Due to their high accuracy, we will call these numerical solutions asexact solutions in the remainder of this subsection.Now, we make a more quantitative comparison between our perturbation-series solution and the exact solution.For this purpose, we first consider the perturbation-series solution (12)-(13) at µ = µ + 0 .
1, i.e., when ǫ = 0 . c , c , c ), ( p , q ), ( p , q ) and ( p , q ) in the previous subsection, and plotted in Fig. 1(d) alongside the exactsolution. As can be seen, this third-order perturbation solution is almost indistinguishable from the exact solution.This is not surprising, since this third-order perturbation solution has relative error of order ǫ , or roughly 0.001 for ǫ = 0 .
1, which is indeed very small.Next, we compare the power function of our perturbation-series solutions (12)-(13) to that of the exact solitonsolutions. For this purpose, we insert the perturbation-series solution (12)-(13) into the power function definition (42)and get P anal ( µ ) = ǫP + ǫ P + ǫ P + · · · , (43)0where µ = µ + ǫ as before, and P = Z ∞−∞ ( p + q ) d x, P = Z ∞−∞ p p + q q ) d x, P = Z ∞−∞ (cid:2) p + q + 2( p p + q q ) (cid:3) d x. Using the ( p , q ), ( p , q ) and ( p , q ) solutions we have numerically obtained, we find that P ≈ . , P ≈ − . , P ≈ − . . (44)Truncating the power-function expansion (43) to the third order, this truncated power function is plotted in Fig. 1(e)alongside the exact power function. Again, the two functions are almost indistinguishable when µ is close to µ .The power series (43) is an asymptotic series. It does not have to be convergent, but it must satisfy the requirementof an asymptotic series, which is that | P ( µ ) − P nk =1 ǫ k P k | = o ( ǫ n ) when ǫ → n [29]. Toverify this asymptotic condition of our power series (43), we examine the difference between the third-order truncatedpower expansion (43) and the exact power function. According to our power expansion, this difference is expected tobe ∆ P ≡ P ( µ ) − ǫP − ǫ P − ǫ P = O ( ǫ ) . (45)If this is indeed true, then the above asymptotic condition for n = 3 would be met. To confirm this ∆ P = O ( ǫ )asymptotics for small ǫ , we show in Fig. 1(f) a log-log plot of ∆ P versus ǫ . Its comparison with the benchmark∆ P = ǫ curve on the same graph shows that this ∆ P is indeed O ( ǫ ) at small ǫ , confirming the asymptotic accuracyof our third-order power expansion. −20 0 2001 x V (a) −20 0 2001 x φ a nd ψ (b) −20 0 2001 x φ A a nd ψ A (c) −20 0 2000.4 x | u | (d) µ P (e) −2 −1 −9 −5 −1 ǫ ∆ P (f) FIG. 1: Comparison of solitons between theory and numerics for the 1D equation (6) with σ = 1 and g ( x ) given by Eq. (41).(a) Wadati potential (2), where solid blue is Re( V ) and dashed red Im( V ). (b) Linear eigenmode φ ( x ) + i ψ (x) of this potential,with solid blue being φ ( x ) and dashed red ψ ( x ). (c) Localized adjoint eigenfunction, with solid blue being φ A ( x ) and dashedred ψ A ( x ). (d) Amplitude profile | u ( x ; µ ) | of the soliton at µ = µ + 0 .
1, where solid blue is from numerical computationand red dots from analytical third-order perturbation series prediction. (e) Power curve of this soliton family, with solid bluefrom numerical computations and red dots from the third-order perturbation expansion (43). (f) Log-log plot of the powerdifference (45) between numerical values and the third-order perturbation expansion versus ǫ = µ − µ . The dashed red line isthe ∆ P = ǫ curve for comparison. In the above numerical example, we chose the focusing nonlinearity (with σ = 1). If the nonlinearity is defocusing,we have found similarly good agreement between perturbation-series solutions and the numerics.1 In the above perturbation calculation, we introduced the tangible small parameter as ǫ = µ − µ . Because of that,we only needed to expand the solutions ( p, q ) into perturbation series. In this treatment, the ( p n , q n ) T solution at eachorder must contain the homogeneous term c n ( φ, ψ ) T , so that c n can be selected judiciously to satisfy the solvabilitycondition of the linear nonhomogeneous ( p n +1 , q n +1 ) T equation.There is an alternative perturbation calculation, where we expand not only the solutions ( p, q ), but also thepropagation constant µ , into perturbation series. In this treatment, the ( p, q ) expansion would still be (12)-(13),while the µ expansion would be µ = µ + µ ǫ + µ ǫ + · · · , (46)where µ , µ , · · · are real constants to be determined. Due to the introduction of these ( µ , µ , · · · ) parameters inthe µ -expansion, we can choose each µ n judiciously to satisfy the solvability condition of the linear nonhomogeneous( p n , q n ) T equation. As a consequence, we do not need to introduce the homogeneous term c n ( φ, ψ ) T in the ( p n , q n ) T solution anymore. In this alternative treatment, we systematically detune the propagation constant µ ; while in theoriginal treatment, we systematically detune the coefficient of the ( φ, ψ ) T term in the ( p, q ) T solution, since thatcoefficient is c + c ǫ + c ǫ + · · · . Algebra-wise, this alternative perturbation calculation turns out to be a littlesimpler, because µ appears in the original two real soliton equations (9)-(10) in a simpler way than p and q , and thusthis µ -detuning introduces less terms in each ( p n , q n ) equation than our present treatment. The slight downside ofthis alternative treatment is that, the “physical” meaning of the small parameter ǫ in it is less clear. Indeed, ǫ inthis alternative treatment is more like a non-tangible arbitrary book-keeping-type small parameter, to which boththe propagation constant µ and the soliton solution ( p, q ) relate in a nontrivial parametric (perturbation-series) way.Overall, these two different perturbation procedures are roughly equivalent, and their choice is largely a personaltaste. Indeed, we have also implemented this alternative perturbation treatment analytically and compared its resultsto the numerics, and found similar agreement as that shown in Fig. 1.
3. CONSTRUCTION OF SOLITON FAMILIES IN THE 2D CASE
Now, we consider the 2D NLS equation iU t + U xx + U yy + V ( x, y ) U + σ | U | U = 0 , (47)where V ( x, y ) is a complex potential, and σ the sign of nonlinearity. It has been shown in [16] that when this potentialis of the form V ( x, y ) = g ( x ) + i g ′ ( x ) + h ( y ) , (48)where g ( x ) and h ( y ) are real functions, then its spectrum can be all-real. This potential is separable, and its x -part isthe 1D Wadati potential (2). So, this 2D potential will also be called Wadati-type in this article. When g ( x ) is even,then this potential admits the partial PT symmetry V ∗ ( x, y ) = V ( − x, y ). In this case, Eq. (47) admits continuousfamilies of solitons, which has been demonstrated numerically and explained analytically in [30]. However, when g ( x )is not even, so that the potential V ( x, y ) is non- PT -symmetric, numerical evidence in [18] indicates that Eq. (47)could still admit continuous families of solitons, which is mysterious in the absence of PT symmetry.In this section, we analytically explain the existence of continuous families of solitons in the 2D NLS equation (47)with a non- PT -symmetric Wadati-type potential (48) by extending the 1D perturbation calculations of the previoussection to the present 2D case. In this potential (48), we require g ( x ) to be localized and differentiable, and h ( y )localized and continuous or piece-wise continuous.Solitons in Eq. (47) are of the form U ( x, y, t ) = e i µt u ( x, y ) , (49)where µ is a real propagation constant, and u ( x, y ) is a localized function satisfying the 2D complex soliton equation (cid:2) ∂ xx + ∂ yy + g ( x ) + i g ′ ( x ) + h ( y ) (cid:3) u − µu + σ | u | u = 0 . (50)2 Similar to the 1D case, an important property of Eq. (47) with the Wadati-type potential (48) is that it admits aconservation law even though it is non-Hamiltonian [18]. Substituting the soliton solution (49) into that conservationlaw, we get a stationary real-valued flux equation ∂J ∂x + ∂J ∂y = 0 , (51)where J = | u x + i g ( x ) u | + [ h ( y ) − µ ] | u | + σ | u | − | u y | , (52)and J = u x u ∗ y + u ∗ x u y − i g ( x )( u y u ∗ − uu ∗ y ) . (53)Following the 1D strategy, instead of working with the complex soliton equation (50), we will work with the real partof that soliton equation, i.e., (cid:2) ∂ xx + ∂ yy + g ( x ) + h ( y ) − µ (cid:3) p − g ′ ( x ) q + σ ( p + q ) p = 0 , (54)where u ≡ p +i q as before [see (8)], together with the real-valued flux equation (51), in our construction of a continuousfamily of 2D solitons. A minor difference from the 1D case is that here, we have to use the flux equation (51), whichis the counterpart of the dJ/dx = 0 equation in the 1D case. This contrasts the 1D case where we used J = 0 directly.This minor difference in the starting equations for solitons will lead to minor differences in the technical constructionsof soliton solutions, as we will see later in this section.The soliton family to be constructed bifurcates from a discrete real eigenvalue µ of the potential. The correspondinglocalized eigenmode ˆ φ ( x, y ) + i ˆ ψ ( x, y ), with real ( ˆ φ, ˆ ψ ), satisfies the linear eigenmode equation obtained by droppingthe nonlinear term in the soliton equation (50), i.e., (cid:2) ∂ xx + ∂ yy + g ( x ) + i g ′ ( x ) + h ( y ) (cid:3) ( ˆ φ + i ˆ ψ ) = µ ( ˆ φ + i ˆ ψ ) . (55)Since the potential in this equation is separable, its linear mode is also separable and can be decomposed asˆ φ ( x, y ) = φ ( x ) ζ ( y ) , ˆ ψ ( x, y ) = ψ ( x ) ζ ( y ) , (56)and µ = µ + µ , where φ ( x ) + i ψ ( x ) is a localized eigenmode of the x -part of the potential (a Wadati potential)with a discrete real eigenvalue µ , i.e., (cid:2) ∂ xx + g ( x ) + i g ′ ( x ) (cid:3) ( φ + i ψ ) = µ ( φ + i ψ ) , (57)and ζ ( y ) is a real localized eigenmode of the y -part of the potential with a discrete real eigenvalue µ , i.e.,[ ∂ yy + h ( y )] ζ = µ ζ. (58)This 2D eigenmode u = ˆ φ ( x, y ) + i ˆ ψ ( x, y ) satisfies the flux equation (51) with the | u | term dropped in J and µ replaced by µ .Bifurcating from this linear eigenmode, we seek a low-amplitude soliton at each real propagation constant value µ near the linear eigenvalue µ , and this soliton is expanded into the following perturbation series, p ( x, y ; µ ) = ǫ / (cid:2) p ( x, y ) + ǫp ( x, y ) + ǫ p ( x, y ) + · · · (cid:3) , (59) q ( x, y ; µ ) = ǫ / (cid:2) q ( x, y ) + ǫq ( x, y ) + ǫ q ( x, y ) + · · · (cid:3) , (60)where ǫ = µ − µ and is assumed to be small positive (so that ǫ / is real). As explained in the 1D case, this positive- ǫ assumption corresponds to a proper sign of nonlinearity σ , and the negative- ǫ case can be treated similarly.Substituting the above perturbation expansion into Eqs. (51) and (54), we get a sequence of real equations for( p k , q k ). The equations for ( p , q ) are just the flux equation (51) with the | u | term dropped in J , and the linearpart of Eq. (54), with µ replaced by µ . Their solutions are obviously (cid:20) p q (cid:21) = c (cid:20) ˆ φ ˆ ψ (cid:21) , (61)3where c is a real constant to be determined. The equations for ( p k , q k ) ( k ≥
1) are the following linear nonhomoge-neous system of equations b L (cid:20) p k q k (cid:21) = (cid:20) ˆ f k ∂ x ˆ g k + ∂ y ˆ g k (cid:21) , (62)where b L is a 2 × b L = ∂ xx + ∂ yy + g + h − µ , b L = − g x , b L = ∂ x h ( ˆ φ x − g ˆ ψ ) ∂ x + g ( ˆ ψ x + g ˆ φ ) + ( h − µ ) ˆ φ − ˆ φ y ∂ y i + ∂ y h ˆ φ y ∂ x + ( ˆ φ x − g ˆ ψ ) ∂ y + g ˆ ψ y i , b L = ∂ x h ( ˆ ψ x + g ˆ φ ) ∂ x − g ( ˆ φ x − g ˆ ψ ) + ( h − µ ) ˆ ψ − ˆ ψ y ∂ y i + ∂ y h ˆ ψ y ∂ x + ( ˆ ψ x + g ˆ φ ) ∂ y − g ˆ φ y i , ˆ f ˆ g ˆ g = c ˆ φ − σc ( ˆ φ + ˆ ψ ) ˆ φ ( ˆ φ + ˆ ψ ) − σc ( ˆ φ + ˆ ψ ) , (63) ˆ f ˆ g ˆ g = (cid:0) − σp − σq (cid:1) p − σp q q c (cid:2) (cid:0) − σp − σq (cid:1) ( p p + q q ) + ( µ − h )( p + q ) + ( p y + q y ) − ( p x − gq ) − ( q x + gp ) ) (cid:3) − c [( p x − gq ) p y + ( q x + gp ) q y ] , (64) ˆ f k ˆ g k ˆ g k = c M c M c M c M c M c M (cid:20) p k − q k − (cid:21) + " b N [1] k b N [2] k , k ≥ , (65)the matrix elements c M ij are k -independent and given by the formulae c M = 1 − σp − σq , c M = − σp q , c M = 1 c (cid:2) p (cid:0) − σp − σq (cid:1) + ( µ − h ) p + p y ∂ y − ( p x − gq ) ∂ x − g ( q x + gp ) (cid:3) , c M = 1 c (cid:2) q (cid:0) − σp − σq (cid:1) + ( µ − h ) q + q y ∂ y + g ( p x − gq ) − ( q x + gp ) ∂ x (cid:3) , c M = ( p x − gq ) ∂ y + p y ∂ x + gq y , c M = ( q x + gp ) ∂ y + q y ∂ x − gp y , and b N [1] k , b N [2] k are functions which depend only on k , p , p , . . . , p k − , q , q , . . . , q k − , g ( x ) and h ( y ). To solve the 2D linear nonhomogeneous equations (62) and obtain localized solutions ( p k , q k ) for all k , we will alsouse the Fredholm alternative method. To do so, we need to understand the kernel structures of the 2D operator b L and its adjoint operator b L A , where elements of the adjoint operator are b L A = ∂ xx + ∂ yy + g + h − µ , b L A = − g x , b L A = h ∂ x ( ˆ φ x − g ˆ ψ ) − g ( ˆ ψ x + g ˆ φ ) − ( h − µ ) ˆ φ − ∂ y ˆ φ y i ∂ x + h ∂ x ˆ φ y + ∂ y ( ˆ φ x − g ˆ ψ ) − g ˆ ψ y i ∂ y , b L A = h ∂ x ( ˆ ψ x + g ˆ φ ) + g ( ˆ φ x − g ˆ ψ ) − ( h − µ ) ˆ ψ − ∂ y ˆ ψ y i ∂ x + h ∂ x ˆ ψ y + ∂ y ( ˆ ψ x + g ˆ φ ) + g ˆ φ y i ∂ y . b L . It is easy to check that this kernel contains two localized functions b K ≡ (cid:20) ˆ φ ˆ ψ (cid:21) , b K ≡ (cid:20) − ˆ ψ ˆ φ (cid:21) , (66)where b L b K = b L b K = 0 , (67)similar to the 1D case and for similar reasons. Since the kernel equation b L b K = 0 is the linearization of the two real“eigenvalue” equations for ( p , q ) [the 2D counterparts of 1D equations (14)-(15)] around the linear mode ( ˆ φ, ˆ ψ ),localized functions in b L ’s kernel can only be induced by amplitude and phase invariances of these ( p , q ) equations,which result in b K and b K above. Thus, there are no other localized functions in b L ’s kernel.Next, we consider the kernel structure of the adjoint 2D operator b L A . Due to the separability of the 2D eigenmode( ˆ φ, ˆ ψ ) in Eq. (56), we can quickly verify that the kernel of b L A contains a bounded function b K A = (cid:20) φ A ( x ) ζ ( y ) − R ψ A ( x )d x (cid:21) , (68)where b L A b K A = 0, and [ φ A ( x ) , ψ A ( x )] T is the unique localized function (29) in the kernel of the 1D adjoint operator L A given in Eq. (23), with µ replaced by µ . One may notice that this kernel function of the 2D adjoint operator doesnot naturally fall back to the 1D adjoint kernel function (29). The reason is twofold. One is that the second columnof the 2D adjoint operator b L A , i.e., [ b L A , b L A ] T given above, contains an additional spatial derivative compared to thesecond column of the 1D adjoint operator L A given in Eq. (23) — a difference caused by our using the divergenceform of the flux equation (51) in 2D instead of its integrated form J ( x ) = 0 in 1D. This difference in the secondcolumn of the adjoint operator explains the integral in the second element of b K A above. The second reason for b K A in 2D not naturally falling back to K A in 1D is that, the second columns of the two adjoint operators contain lineareigenmodes or their derivatives as multiplicative factors, while the first columns of these adjoint operators do not.Thus, in the 2D case, we need to introduce the factor ζ ( y ) from the 2D linear eigenmode (56) into the first elementof the adjoint kernel function b K A in Eq. (68) in order to balance such a term coming from the second column of b L A .We can further show that, if h ( y ) is a slowly varying function, then the above b K A would be the only boundedfunction in the kernel of b L A . To do so, let h ( y ) = ˆ ǫ H ( Y ) be a slowly varying function of Y = ˆ ǫy , where ˆ ǫ is a smallreal parameter. For this h ( y ), its eigenmode from Eq. (58) is ζ ( y ) = ˆ ζ ( Y ), with eigenvalue µ = O (ˆ ǫ ). In this case, b L A can be rewritten as a quadratic function of ˆ ǫ , b L A = b L A ( x, Y ) + ˆ ǫ b L A ( x, Y ) + ˆ ǫ b L A ( x, Y ) , (69)where b L A ( x, Y ) = (cid:20) ∂ xx + g − µ ˆ ζ ( Y ) [ − ∂ x ( φ ′ − gψ ) + g ( ψ ′ + gφ ) − µ φ ] ( − ∂ x ) − g ′ ˆ ζ ( Y ) [ − ∂ x ( ψ ′ + gφ ) − g ( φ ′ − gψ ) − µ ψ ] ( − ∂ x ) (cid:21) = L A (cid:20) − ˆ ζ ( Y ) ∂ x (cid:21) , (70)and L A is the 1D adjoint operator (23) with µ replaced by µ . Since b L A is a function of x , Y and ˆ ǫ , functions b F inits kernel are also functions of these same variables and can be expanded into a perturbation series of ˆ ǫ as b F ( x, y ; ˆ ǫ ) = b F ( x, Y ) + ˆ ǫ b F ( x, Y ) + ˆ ǫ b F ( x, Y ) + . . . . (71)Inserting this expansion and Eq. (69) into b L A b F = 0 and using the kernel structures of the 1D operators L and itsadjoint L A detailed in Sec. 2 2.2, we can sequentially determine b F n ( x, Y ) in the above perturbation expansion andshow that the only bounded function in the kernel of b L A is b F = (cid:20) φ A ( x ) ˆ ζ ( Y ) − R ψ A ( x ) d x (cid:21) , (72)which matches (68) when the eigenmode ζ ( y ) = ˆ ζ ( Y ) is slowly varying. All other functions in the kernel of b L A growexponentially at large | x | or | Y | .5When h ( y ) continuously deforms from slowly varying to the general case of non-slowly varying, the above kernelstructure of b L A generically will not change, i.e., its kernel will generically still contain a single bounded function(68). While we cannot at this time rule out the possibility of additional bounded functions appearing in the kernelof b L A at some special h ( y ) functions during this deformation process, for specific examples of the potentials, we canuse numerics to directly verify this single-bounded-function kernel structure for b L A , so that our analysis below canproceed. With the above kernel structures of b L and b L A in hand, we can now sequentially solve Eq. (62) for localized solutions( p k , q k ) using the Fredholm alternative method. According to this method, if functions ( ˆ f k , ˆ g k , ˆ g k ) on the right sideof the linear nonhomogeneous system (62) are localized (which is the case here), this system would admit a localizedsolution ( p k , q k ) if and only if its right hand side is orthogonal to the bounded function b K A of (68) in the kernel of b L A , i.e., (cid:28)(cid:20) φ A ( x ) ζ ( y ) − R ψ A ( x ) d x (cid:21) , (cid:20) ˆ f k ∂ x ˆ g k + ∂ y ˆ g k (cid:21)(cid:29) = 0 . (73)It is noted that the arbitrary constant out of the indefinite integral R ψ A ( x ) d x gives no contribution to the innerproduct in the above solvability condition. In addition, the above integral is convergent since ( ˆ f k , ˆ g k , ˆ g k ) are alllocalized in space.Our perturbative construction of 2D solitons bifurcating from a linear localized eigenmode of the complex potential(48) proceeds similarly as the 1D case, since the kernel structures in the 2D case resemble those in the 1D case. Wefirst consider Eq. (62) for ( p , q ). Substituting the ( ˆ f , ˆ g , ˆ g ) expressions (63) into the above solvability conditionand performing integration by parts, we get (cid:28)(cid:20) φ A ( x ) ζ ( y ) ψ A ( x ) (cid:21) , (cid:20) ˆ f ˆ g (cid:21)(cid:29) = 0 . (74)Inserting the ( ˆ f , ˆ g ) expressions (63) and ( ˆ φ, ˆ ψ ) formulae (56) into the above equation, we obtain a formula for c as c = ± vuut R ∞−∞ (cid:2) φφ A + ( φ + ψ ) ψ A (cid:3) d xσ R ∞−∞ (cid:2) ( φ + ψ ) φφ A + ( φ + ψ ) ψ A (cid:3) d x R ∞−∞ ζ ( y ) d y R ∞−∞ ζ ( y ) d y . (75)As in the 1D case, the sign of σ must match the sign of the ratio between integrals in the above equation so that thequantity under the square root is positive. In addition, we can choose the plus sign outside the square root withoutloss of generality.When c is selected from the above formula (75), Eq. (62) admits a localized solution for ( p , q ), which we denoteas ( p s , q s ). Since the kernel of the homogeneous operator b L in Eq. (62) contains two localized functions b K and b K given in Eq. (66), the general localized solution ( p , q ) to the linear nonhomogeneous equations (62) is then[ p s , q s ] T + c b K + d b K , where c and d are two real constants. But as in the 1D case, the d term can be removedby phase invariance of the complex soliton solution u ( x, y ). Thus, the ( p , q ) solution can be set as (cid:20) p q (cid:21) = (cid:20) p s q s (cid:21) + c (cid:20) ˆ φ ˆ ψ (cid:21) . (76)The constant c in this solution will be determined from the Fredholm solvability condition on the ( p , q ) equations.The equations for ( p , q ) are (62), where ( ˆ f , ˆ g , ˆ g ) in the nonhomogeneous terms are given in Eq. (64). Substi-tuting the ( p , q ) solutions (61) and ( p , q ) solutions (76) into these nonhomogeneous terms and recalling that theeigenmode u = ˆ φ ( x, y ) + i ˆ ψ ( x, y ) satisfies the flux equation (51) with the | u | term dropped in J and µ replaced by µ , we see that the right side of Eq. (62) for ( p , q ) reduces to (cid:20) ˆ f ∂ x ˆ g + ∂ y ˆ g (cid:21) = (cid:20) ˆ f a ∂ x ˆ g a + ∂ y ˆ g a (cid:21) + c (cid:20) ˆ f b ∂ x ˆ g b + ∂ y ˆ g b (cid:21) , (77)6where ˆ f a = (cid:0) − σp − σq (cid:1) p s − σp q q s , ˆ g a = 12 c (cid:2) (cid:0) − σp − σq (cid:1) ( p p s + q q s ) + ( µ − h )( p s + q s ) + ( p s,y + q s,y ) − ( p s,x − gq s ) − ( q s,x + gp s ) ) (cid:3) , ˆ g a = − c [( p s,x − gq s ) p s,y + ( q s,x + gp s ) q s,y ] , and ˆ f b = (1 − σp − σq ) ˆ φ − σp q ˆ ψ, ˆ g b = 1 c h(cid:0) − σp − σq (cid:1) ( p ˆ φ + q ˆ ψ ) + ( µ − h )( p s ˆ φ + q s ˆ ψ ) + ( p s,y ˆ φ y + q s,y ˆ ψ y ) − ( p s,x − gq s )( ˆ φ x − g ˆ ψ ) − ( q s,x + gp s )( ˆ ψ x + g ˆ φ ) i , ˆ g b = − c h ( p s,x − gq s ) ˆ φ y + p s,y ( ˆ φ x − g ˆ ψ ) + ( q s,x + gp s ) ˆ ψ y + q s,y ( ˆ ψ x + g ˆ φ ) i . Inserting (77) into the Fredholm solvability condition (73) at k = 2, we get a formula for the constant c as c = − (cid:28)(cid:20) φ A ( x ) ζ ( y ) − R ψ A ( x )d x (cid:21) , (cid:20) ˆ f a ∂ x ˆ g a + ∂ y ˆ g a (cid:21)(cid:29)(cid:28)(cid:20) φ A ( x ) ζ ( y ) − R ψ A ( x )d x (cid:21) , (cid:20) ˆ f b ∂ x ˆ g b + ∂ y ˆ g b (cid:21)(cid:29) . (78)When the c value is selected as above, the ( p , q ) solutions (76) are completely determined. In addition, theFredholm solvability condition (73) for the ( p , q ) equations (62) is also satisfied; so these equations admit a localized( p , q ) solution, which we denote as ( p s , q s ). In view of the kernel structure of operator b L and phase invariance ofthe complex soliton solution u ( x, y ), the general localized solutions for ( p , q ) can be written as (cid:20) p q (cid:21) = (cid:20) p s q s (cid:21) + c (cid:20) ˆ φ ˆ ψ (cid:21) , (79)where c is a real constant. This constant c will be determined from the Fredholm solvability condition for the ( p , q )equations. Indeed, inserting this ( p , q ) solution into the right side of Eq. (62) with k = 3, it is easy to see thatthe solvability condition (73) at k = 3 is a linear equation for c , which we can easily solve to obtain the value of c .After this c value is obtained, ( p , q ) is ascertained. In addition, the ( p , q ) equation admits a localized solution,which we denote as ( p s , q s ), and the general ( p , q ) solutions can be written as (79) with the index changed from 2to 3. This process is then repeated to higher orders. Lastly, we compare the above 2D perturbation series soliton solution (59)-(60) with the high-accuracy numericalsolution and confirm the asymptotic accuracy of this 2D analytical solution. In our comparison, we choose thepotential (48) with g ( x ) = 0 . x + 2) + 1 . x − , h ( y ) = 2sech y. (80)Notice that this g ( x ) function is the same as (41) in the 1D example. This potential admits a discrete real eigenvalue µ = µ + µ ≈ . µ ≈ . µ = 1. The corresponding eigenfunction( ˆ φ, ˆ ψ ) is given in Eq. (56), where [ φ ( x ) , ψ ( x )] is as shown in Fig. 1(b), and ζ ( y ) = sech( y ). Numerically, we confirmedthat the 2D adjoint operator b L A indeed admits a single bounded function (68) in its kernel, where [ φ A ( x ) , ψ A ( x )] isthe localized function (29) in the kernel of the 1D adjoint operator L A in Eq. (23), which was plotted in Fig. 2(c).When σ = 1 (focusing nonlinearity), our theory predicts that a continuous family of solitons bifurcates out fromthe above linear discrete eigenmode when µ > µ . This is indeed the case. For the choice of µ = µ + 0 . ǫ = 0 . c , c , p , q , p and q in the previous subsection, and plotted in Fig. 2(a). The high-accuracy numerical solution atthis same µ value is displayed in Fig. 2(b) for comparison. It is seen that these two solutions are visually identical.We have also calculated the difference between these two solutions, and found that the relative error between them isunder 2 . O ( ǫ ) (i.e., order of 0.01) as expected.Next, we compare the power function of our perturbation-series solutions (59)-(60) to that of the exact solitonsolutions. This 2D power function is defined as P ( µ ) = Z ∞−∞ Z ∞−∞ | u ( x, y ; µ ) | d x d y, (81)analogous to the 1D case (42). Inserting the perturbation-series solution (59)-(60) into this power function, we get P anal ( µ ) = ǫP + ǫ P + · · · , (82)where P = Z ∞−∞ Z ∞−∞ ( p + q ) d x d y, P = Z ∞−∞ Z ∞−∞ p p + q q ) d x d y. (83)Using the ( p , q ) and ( p , q ) solutions obtained from Eqs. (61) and (76), we find that P ≈ . , P ≈ − . . Truncating the power-function expansion (82) to these first two terms, this truncated power function is plotted inFig. 2(c) alongside the exact power function. Again, the two functions are almost indistinguishable. To verify theasymptotic accuracy of our perturbation series solutions, we show in Fig. 2(d) a log-log plot of ∆ P ≡ ǫP + ǫ P − P ( µ )versus ǫ . Its comparison with the benchmark ∆ P = ǫ curve on the same graph shows that this ∆ P is O ( ǫ ), whichmatches our asymptotic prediction for this quantity. The above comparison indicates that the true 2D soliton solutionsand our perturbation series solutions are in perfect agreement. xy (a) −10 0 10−10010 xy (b) −10 0 10−100101.4 1.5 1.60123 µ P (c) −3 −2 −1 −9 −5 −1 ǫ ∆ P (d) FIG. 2: Comparison of solitons between theory and numerics for the 2D equation (50) with σ = 1, and g ( x ), h ( y ) givenby Eq. (80). (a) Amplitude profile | u ( x, y ; µ ) | of the second-order perturbation series solution (82) at µ = µ + 0 .
1. (b)Numerically computed soliton | u | at the same µ value of (a). (c) Power curve of this soliton family, with solid blue fromnumerical computations and red dots from the second-order perturbation expansion (82). (d) Log-log plot of the powerdifference between numerical values and the second-order perturbation expansion versus ǫ = µ − µ . Dashed red line is the∆ P = ǫ curve for comparison.
4. SUMMARY AND DISCUSSION
In this article, we have analytically constructed continuous families of low-amplitude solitons bifurcating from linearmodes in one- and two-dimensional NLS equations with localized Wadati-type non- PT -symmetric complex potentials,thus providing an analytical explanation for this counter-intuitive phenomenon of soliton families appearing in thesenon- PT -symmetric non-Hamiltonian systems. Our analytical construction utilized the conservation laws of thesenon- PT -symmetric equations, which allowed us to convert the complex soliton equations into new real systems. Akey advantage of these new real systems is that, during a perturbation expansion of low-amplitude solitons bifurcatingfrom linear modes, the underlying linear operator has two localized functions in its kernel, and the associated adjointoperator has a single localized or bounded function in its kernel. This kernel structure, coupled with the phaseinvariance of the complex soliton, guarantees that at each order of the soliton’s perturbation expansion, the Fredholmsolvability condition can always be satisfied, so that a localized solution at each order of the perturbation series canbe found. As a result, a continuous family of low-amplitude solitons bifurcating from a linear mode is obtained as aperturbation series to all orders of the small soliton amplitude. We have also compared these analytically constructedsoliton solutions to high-accuracy numerical solutions, in both one and two spatial dimensions, and the asymptoticaccuracy of these perturbation solutions is fully confirmed.In this article, the nonlinearity in our 1D and 2D NLS equations (1) and (47) is cubic. But our analytical treatmentfor this cubic nonlinearity can be trivially generalized to other types of nonlinearities of the general form G ( | U | ) U ,where G ( · ) is an arbitrary real function. Indeed, for the 1D and 2D NLS equations (1) and (47) with this more generalform of nonlinearity but the same Wadati-type complex potentials (2) and (48), a conservation law still exists [18, 19].Thus, the analytical treatment of this article still applies.In our perturbative construction of soliton families in the NLS equations (1) and (47) with non- PT -symmetricWadati-type potentials, the conservation laws of those equations played a critical role. If such conservation laws areabsent, such as for non- PT -symmetric complex potentials not of Wadati-type, this construction would not work. Insuch cases, we do not believe true soliton families can still exist. This implies that we do not think the “solitonfamilies” reported in [25] for non-Wadati complex potentials are true soliton solutions.A closely related subject is symmetry breaking of solitons in PT -symmetric Wadati-type potentials (2) and (48),where g ( x ) is an even function. It is known that for generic PT -symmetric potentials, symmetry breaking of solitonsis forbidden. However, for PT -symmetric Wadati-type potentials (2) and (48), symmetry breaking of solitons canoccur, where two branches of non- PT -symmetric solitons bifurcate out from the base branch of PT -symmetric solitonswhen the base branch’s power reaches a certain threshold [18, 31]. So far, there has been no analytical explanationfor these symmetry breakings. For PT -symmetric Wadati-type potentials, the conservation laws (3) and (51) are stillvalid. Then, using our new real system of soliton equations in this article, together with bifurcation conditions forsymmetry breaking, branches of symmetric and asymmetric solitons in these PT -symmetric Wadati-type potentialscould be perturbatively constructed near the symmetry-breaking point. Details of this construction will be left forfuture studies.The analytical construction of soliton solutions is often a precursor of the subsequent linear stability analysis ofthese solitons. Thus, the results of this article could be helpful for the analytical stability investigations of solitons inWadati-type complex potentials. Acknowledgement
This material is based upon work supported by the Air Force Office of Scientific Research under award numberFA9550-18-1-0098 and the National Science Foundation under award number DMS-1910282. [1] V.V. Konotop, J. Yang and D.A. Zezyulin, “Nonlinear waves in PT -symmetric systems”, Rev. Mod. Phys. 88 035002(2016).[2] S.V. Suchkov, A.A. Sukhorukov, J. Huang, S.V. Dmitriev, C. Lee and Y.S. Kivshar, “Nonlinear switching and solitons in PT -symmetric photonic systems”, Laser Photon. Rev. 10, 177 (2016).[3] D. N. Christodoulides and J. Yang (Eds.), Parity-time Symmetry and Its Applications (Springer, 2018).[4] C.M. Bender and S. Boettcher, “Real spectra in non-Hermitian Hamiltonians having PT symmetry”, Phys. Rev. Lett. 80,5243-5246 (1998). [5] A. Mostafazadeh, “Pseudo-Hermitian representation of quantum mechanics”, Int. J. Geom. Methods Mod. Phys. 7, 1191-1306 (2010).[6] C.M. Bender, PT Symmetry in Quantum and Classical Physics (World Scientific, London, 2019).[7] Z. H. Musslimani, K. G. Makris, R. El-Ganainy, and D. N. Christodoulides, “Optical solitons in PT periodic potentials”,Phys. Rev. Lett. 100, 030402 (2008).[8] S. Longhi, “ PT -symmetric laser absorber”, Phys. Rev. A 82, 031801 (2010).[9] Y.D. Chong, L. Ge and A.D. Stone, “PT-Symmetry breaking and laser-absorber modes in optical scattering systems”,Phys. Rev. Lett. 106, 093902 (2011).[10] Z.J. Wong, Y.L. Xu, J. Kim, K. O’Brien, Y. Wang, L. Feng and X. Zhang, “Lasing and anti-lasing in a single cavity”,Nat. Photonics 10, 796-801 (2016).[11] L. Feng, Z.J. Wong, R. Ma, Y. Wang and X. Zhang, “Single-mode laser by parity-time symmetry breaking”, Science 346,972-975 (2014).[12] H. Hodaei, M.A. Miri, M. Heinrich, D.N. Christodoulides and M. Khajavikhan, “Parity-time-symmetric microring lasers”,Science 346 975-978 (2014).[13] F. Cannata, G. Junker and J. Trost, “Schr¨odinger operators with complex potential but real spectrum”, Phys. Lett. A246, 219-226 (1998).[14] M. Miri, M. Heinrich and D.N. Christodoulides, “Supersymmetry-generated complex optical potentials with real spectra”,Phys. Rev. A 87, 043819 (2013).[15] E.N. Tsoy, I.M. Allayarov and F. Kh. Abdullaev, “Stable localized modes in asymmetric waveguides with gain and loss”,Opt. Lett. 39, 4215 (2014).[16] S. Nixon and J. Yang, “All-real spectra in optical systems with arbitrary gain-and-loss distributions”, Phys. Rev. A 93,031802(R) (2016).[17] V.V. Konotop and D.A. Zezyulin, “Families of stationary modes in complex potentials”, Opt. Lett. 39, 5535-5538 (2014).[18] J. Yang, “Symmetry breaking of solitons in two-dimensional complex potentials”, Phys. Rev. E 91, 023201 (2015).[19] S. Nixon and J. Yang, “Bifurcation of soliton families from linear modes in non- PT -symmetric complex potentials”, Stud.Appl. Math. 136, 459-483 (2016).[20] J. Yang and S. Nixon, “Stability of soliton families in nonlinear Schr¨odinger equations with non-parity-time-symmetriccomplex potentials”, Phys. Lett. A 380, 3803-3809 (2016).[21] N. Akhmediev and A. Ankiewicz (Eds), Dissipative Solitons (Springer, Berlin, 2005).[22] J. Cuevas-Maraver, P.G. Kevrekidis, D.J. Frantzeskakis, and Y. Kominis, “Nonlinear beam propagation in a class ofcomplex non- PT -symmetric potentials”, in D. Christodoulides, J. Yang (eds.), Parity-time Symmetry and Its Applications (Springer Tracts in Modern Physics 280, 2018), pp. 557-579.[23] J. Yang, “Necessity of PT symmetry for soliton families in one-dimensional complex potentials”, Phys. Lett. A 378, 367-373(2014).[24] C. Hang, G. Gabadadze and G. Huang, “Realization of non- PT -symmetric optical potentials with all-real spectra in acoherent atomic system”, Phys. Rev. A 95, 023833 (2017).[25] Y. Kominis, J. Cuevas-Maraver, P.G. Kevrekidis, D.J. Frantzeskakis and A. Bountis, “Continuous families of solitary wavesin non-symmetric complex potentials: A Melnikov theory approach”, Chaos, Solitons and Fractals 118, 222-233 (2019).[26] M. Wadati, “Construction of parity-time symmetric potential through the soliton theory”, J. Phys. Soc. Jpn. 77, 074005(2008).[27] D.E. Pelinovsky, Localization in Periodic Potentials: From Schrodinger Operators to the Gross-Pitaevskii Equation (Cam-bridge University Press, 2011).[28] H. Brezis,
Functional Analysis, Sobolev Spaces and Partial Differential Equations (Springer, New York 2011).[29] C.M. Bender and S.A. Orszag,