A solvable model for solitons pinned to a PT-symmetric dipole

We introduce the simplest one-dimensional nonlinear model with the parity-time (PT) symmetry, which makes it possible to find exact analytical solutions for localized modes ("solitons"). The PT-symmetric element is represented by a point-like (delta-functional) gain-loss dipole {\delta}^{\prime}(x), combined with the usual attractive potential {\delta}(x). The nonlinearity is represented by self-focusing (SF) or self-defocusing (SDF) Kerr terms, both spatially uniform and localized ones. The system can be implemented in planar optical waveguides. For the sake of comparison, also introduced is a model with separated {\delta}-functional gain and loss, embedded into the linear medium and combined with the {\delta}-localized Kerr nonlinearity and attractive potential. Full analytical solutions for pinned modes are found in both models. The exact solutions are compared with numerical counterparts, which are obtained in the gain-loss-dipole model with the {\delta}^{\prime}- and {\delta}- functions replaced by their Lorentzian regularization. With the increase of the dipole's strength, {\gamma}, the single-peak shape of the numerically found mode, supported by the uniform SF nonlinearity, transforms into a double-peak one. This transition coincides with the onset of the escape instability of the pinned soliton. In the case of the SDF uniform nonlinearity, the pinned modes are stable, keeping the single-peak shape.

The optical realizations of the PT symmetry suggest additional interest to nonlinearity in these systems [8], in which stable solitons can be supported by the combination of the Kerr (cubic) nonlinearity and spatially periodic complex potential, whose odd (antisymmetric) imaginary part accounts for the balanced gain and loss, as mentioned above. Stability of such solitons was rigorously analyzed in Ref. [9]. Dark solitons were also investigated in in the framework of models combining the PT symmetry and self-defocusing Kerr nonlinearity [10]. In addition to that, bright solitons were predicted too in PT -symmetric systems with the quadratic (second-harmonic-generating) nonlinearity [11].
Solitons can also be found in linearly-coupled dual-core systems, with the balanced gain and loss acting in the two cores, and the intrinsic Kerr nonlinearity present in each one [12,13].
Further, discrete solitons were predicted in various models based on chains of linear [15] and circular [16] coupled PT -symmetric elements and, more generally, in networks of coupled PT -symmetric oligomers (dimers, quadrimers, etc.) [17]. Parallel to incorporating the usual Kerr nonlinearity into the conservative part of the system, its gain-loss-antisymmetric part can be made nonlinear too, by introducing mutually balanced cubic gain and loss terms [18]. Effects of combined linear and nonlinear PT terms on the existence and stability of optical solitons were studied too [19].
Unlike the usual nonlinear dissipative systems, where solitons exist as isolated solutions (attractors) [20,21], in PT -symmetric settings solitons form continuous families, similar to the generic situation in conservative media. However, the increase of the gain-loss coefficient (γ) in the PTsymmetric nonlinear system leads to shrinkage of existence and stability areas for PT -symmetric solitons, until they vanishes when this coefficient attains a critical value, γ cr .
In all the previously studied models of PT -symmetric nonlinear systems, except for the simplest dual-core model considered in Ref. [12], solitons could only be constructed and investigated in a numerical form (in Ref. [12], the solutions for PT -symmetric solitons were tantamount to those for symmetric solitons in the usual coupler model, the difference being in their stability). The objective of the present work is to propose a solvable one-dimensional nonlinear model with the PT -balanced gain and loss concentrated at a single point, in the form of a PT dipole. A possibility to construct such a tractable model is suggested by recently studied models of dissipative systems (not subject to the condition of the PT symmetry), in which localized gain, competing with spatially uniform loss, was applied at a single [22] or two [23] "hot spots", making it possible to find exact solutions for dissipative solitons pinned to those spots. Similar models, but with hot spots of a finite width, were investigated by means of numerical methods in Refs. [24].
A PT -symmetric model with the uniformly distributed nonlinearity and localized mutually balanced gain and loss, applied at two points in the form of δ-functions, along with the attractive potential, represented by a pair of δ-functions placed at the same points, was elaborated in Ref. [25]. The model dealt with in the present work may be considered as a limit form of the one introduced in Ref. [25], for a vanishingly small separation between the two δ-functions, when the balanced pair of δ-function-shaped gain and loss go over into a single term in the underlying propagation equation, represented by the δ ′ -function (the "PT dipole").
The model with the PT dipole embedded into the Kerr-nonlinear medium admits a full family of analytical solutions for localized modes ("solitons") pinned to the point-like dipole ("defect"), for both signs of the nonlinearity of the host medium, self-focusing (SF) and self-defocusing (SDF).
Unlike what occurs in other PT -symmetric systems, the analytical solutions exist at all values of the gain-loss parameter, γ, without featuring the above-mentioned threshold (critical value), γ cr .
On the other hand, our numerical solution for the model with the δ-and δ ′ -functions replaced by their finite-width regularizations [see Eq. (32) below] demonstrate that, while at γ small enough the numerical solutions are close to their analytical counterparts found for the ideal δ-functions, their stability and existence are always limited by finite γ cr .
The appearance of the threshold with the introduction of the regularization, i.e., finite separation between the gain and loss, can also be explained in an analytical form. To this end, we introduce an additional model, based on the linear host medium with an embedded pair of separated δfunctions, which carry, in addition to the gain and loss with equal coefficients Γ, the local cubic SF nonlinearity, along with the linear δ-functional potentials. This setting also admits a full analytical solution (cf. Refs. [26] and [27], where exact solutions were found for one-and two-component symmetric, antisymmetric, and asymmetric solitons pinned to a pair of points with the localized SF nonlinearity, embedded into the linear medium). Analytical solutions obtained in the system with the separated δ-functions explicitly feature a threshold value, Γ cr , which bounds their existence region.
The solvable models and analytical solutions are introduced in Section II. Their numerically found counterparts, corresponding to the regularized δ ′ and δ-functions, are presented in Section III. The numerical analysis pursues two objectives: to estimate the robustness (structural stability) of the analytical solutions for the localized modes, obtained with the ideal δ ′ -and δ-functions, and to test the dynamical stability of the modes by means of systematic simulations of the perturbed evolution, which is a crucially important issue in the context of PT -symmetric systems. The paper is concluded by Section IV.

II. THE MODEL AND ANALYTICAL SOLUTIONS
A. The basic model The underlying equation for the light propagation in a nonlinear planar waveguide, with the PT -symmetric complex potential concentrated near x = 0 (with even real and odd imaginary parts), and the uniform cubic nonlinearity with coefficient σ, is Here z is the propagation distance and x the transverse coordinate, with term u xx accounting for the transverse diffraction in the paraxial approximation, while ε 0 > 0 and γ > 0 are strengths of the real and imaginary parts of the complex potential, and ε 2 represents a possible nonlinear part of the trapping potential [26]. By means of obvious rescaling, one can set |σ| ≡ 1, with σ = +1 and −1 corresponding, respectively, to the SF an SDF spatially uniform nonlinearity. In addition, σ = 0 is possible too, corresponding to the model with the Kerr nonlinearity fully concentrated at the same spot where the PT dipole is set. Rescaling also allows us to fix ε 0 ≡ 1, unless ε 0 = 0, in which case it is possible to fix γ ≡ 1. The sign of ε 2 in Eq. (1) may be either the same as σ = ±1 or opposite to it, the latter case corresponding to the competition between the uniform and localized nonlinearities.
Stationary PT -symmetric localized solutions to Eq. (1) are looked for as with real propagation constant k > 0, where complex function U (x) obeys the following equation: At x = 0, PT -symmetric solutions of Eq. (4) for localized modes are constructed in terms of the commonly known analytical expressions for regular and singular solitons of the nonlinear Schrödinger (NLS) equation with the SF or SDF nonlinearity (σ = +1 and σ = −1, respectively): where θ and ξ are free real parameters. The form of this solution implies that Im (U (x = 0)) = 0, while jumps (∆) of the imaginary part and first derivative of the real part at x = 0 are determined by the integration of the δ-and δ ′ -functions in an infinitesimal vicinity of x = 0: B. The analytical solution for the model with the spatially uniform nonlinearity (ε 2 = 0) Substituting bulk solutions (5) and (6) into boundary conditions (7) and (8) with ε 2 = 0, it is straightforward to obtain the following results, which determine the free constants in the solutions, θ and ξ, as functions of k: which does not depend on k and is the same for σ = ±1, and The total power of the localized mode is As seen from Eq. (10), the solutions exist at As concerns stability of the solutions, it is relevant to mention that expression (11) with σ = +1 and −1 satisfy, respectively, the Vakhitov-Kolokolov (VK) [28] and "anti-VK" [29] criteria, i.e., which are necessary conditions for the stability of localized modes supported, severally, by the SF and SDF nonlinearities, hence in both cases the present solutions have a chance to be stable.
C. Analytical solutions for the model with the inhomogeneous nonlinearity (ε 2 = 0) In the presence of the localized nonlinearity, Eq. (9) remains the same as before, while expression (10), following from Eq. (8), is replaced by a rather cumbersome expression: These solutions may be free of singularities and stable only if they yield ξ > 0. In the case of ε 2 > 0 (the attractive nonlinear potential placed at x = 0), the condition of ξ > 0 for Eq. (14) with σ = +1 holds only with sign + in front of the radical, while Eq. (14) with σ = −1 may give rise to two different solutions, corresponding to both signs + and −. In the case of ε 2 < 0, the situation is opposite: Eq. (14) with σ = −1 makes sense only with sign + in front of the radical, while σ = +1 may generate meaningful solutions for both signs + and −. Thus, two different solutions may exist in the case of the competition between the uniform and localized nonlinearities.
It is also relevant to consider the special case of ε 0 = 0, ε 2 > 0, when the attractive potential of the defect at x = 0 is purely nonlinear. In this situation, solution (14) essentially simplifies, taking the same form for σ = ±1: The corresponding expressions for the total power can be found too, cf. Eq. (11): Note that this expression depends on γ, unlike its counterpart (11).
Solution (15) exists for all values of k > 0, unlike the one given by Eq. (10), whose existence region is limited by condition (12). Further, expression (16) satisfies the VK and anti-VK criteria (13), severally for σ = +1 and σ = −1, hence in both cases solution (15) may be stable. For the nonlinear PT dipole embedded into the linear medium, it is possible to fix ε 2 = ±1, for the SF and SDF localized nonlinearity, respectively. The solution of Eq. (4) for the trapped mode is simple in this case: with the corresponding total power Like the above solution given by Eq. (10), and unlike the one amounting to Eq. (15), the existence of this solution is limited by conditions √ 2k > ε 0 and √ 2k < ε 0 , respectively, in the case of the SF and SDF localized nonlinearity, cf. Eq. (12). Further, as well as the solution families considered above, relation (18) satisfies the VK and anti-VK criteria [see Eq. (13)] for the SF and SDF signs of the localized nonlinearity, i.e., ε 2 = +1 and ε 2 = −1. With ε 0 = 0, Eq. (18) yields the degenerate dependence, dP 0 /dk = 0, which formally implies VK-neutral stability, but in reality the solitons are unstable in this case [30].
E. The model with the separated gain and loss embedded into the linear medium As explained above, it is relevant to supplement the PT -dipole model by a solvable one which features a finite separation, 2l, between the mutually balanced gain and loss δ-like elements with equal strengths Γ. Such a system may be built following the lines of Ref. [25], but replacing the uniform Kerr nonlinearity by its counterpart localized at the same points where the gain and loss are set (otherwise, the system is not analytically solvable): Using obvious rescaling, we can fix here |ε 2 | = 1 for the SF (ε 2 = +1) and SDF (ε 2 = −1) signs of the nonlinearity.
Stationary PT -symmetric localized solutions to Eq. (19) are looked in the same form (2) and (3) as above, with U (x) obeying equation At |x| > l and |x| < l, respectively, PT -symmetric solutions of Eq. (4) for U (x) are constructed as follows, cf. Eq. (17): where A, B and C, D are real amplitudes. The condition of the continuity of solution (21) at x = ±l yields a relation eliminating C and D in favor of A and B: The condition for the jump of the first derivatives induced by the delta-functions at x = ±l, cf.
Eq. (8), can be written as a single cubic complex equation for A and B: Equation (23) may be considered as a system of two homogeneous equations for A and B, a nontrivial solution to which exists when the system's determinant vanishes. After some algebra, an explicit solution of Eq. (23) can be obtained: These solutions exist only in the region where the radical in Eq. (26) is real, i.e., √ 2k Condition (27) holds at k ≤ k max , with k max determined by a transcendental equation, It is easy to see that Eq. (28) has a single physical solution provided that and no solutions at Γ > Γ cr , which is a manifestation of the above-mentioned generic feature of nonlinear PT -symmetric systems: soliton families exist below a certain critical value of the gainloss coefficient. On the other hand, Eq. (29) demonstrates that the critical value diverges in the limit of l → 0, which corresponds to the replacement of the separated gain and loss by the PT dipole in Eqs. (1) and (4). The latter fact helps to understand why the above analytical solutions, found in the PT -dipole models, do not feature the existence threshold.

III. NUMERICAL RESULTS FOR THE REGULARIZED MODEL
Numerical results are presented below, chiefly, for solutions which are counterparts of the analytical ones obtained above in the fully explicit form, i.e., the solutions based on Eqs. (5), (6), (9) and (10) (for the spatially uniform nonlinearity, with ε 2 = 0) or (15), for ε 0 = 0, i.e., the purely nonlinear attractive potential at x = 0.
A. The approximation of the δ-function in numerical solutions As said above, the numerical analysis of the model aims to obtain solutions for the δ-function replaced by its finite-width regularization,δ(x), with the objectives to produce solutions for a modification of the model relevant for the experimental implementation, and also to test the stability of the PT -symmetric modes produced above in the analytical form. In many works,δ(x) was used in the form of a narrow Gaussian, see, e.g., Ref. [26]. However, this is not convenient in the present context, as, replacing the exact solutions in the form of Eqs. (5), (6), and (17) by regularized expressions, it is necessary, inter alia, to replace sgn(x) ≡ −1 + 2 x −∞ δ(x ′ )dx ′ by a continuous function realized as −1 + 2 x −∞δ (x ′ )dx ′ , which would imply using a non-elementary function in the case of the Gaussian. Therefore, we here use the regularization in the form of the Lorentzian, with 0 < a ≪ k −1/2 .
B. The self-focusing uniform nonlinearity (σ = +1, ε 2 = 0) We start the presentation of the results with the case of the uniform SF nonlinearity, fixing ε 0 = 1 in Eqs. (1) and (4) (larger values of ε 0 are used below to report results obtained in the model with the uniform SDF nonlinearity). Stationary solutions were found by solving Eq. (4) with the help of the Newton's method, using the input provided by the analytical solution in the form of Eqs. (5), (9), and (10), with the regularization implemented as per Eq. (32). The stability of the so generated solutions was tested through direct simulations of their perturbed evolution by means of the fourth-order split-step method. It was implemented in domain −10 < x < +10, with periodic boundary conditions (the width of the integration domain is definitely much larger than the size of all the modes considered in this work, see Figs. 1, 2, and 7 below). Sufficient numerical stability and accuracy were achieved with time and space step sizes ∆t = 0.001 and ∆x = 0.039, respectively. Accordingly, values a ≥ 0.02 of the regularization scale were adopted in Eq. (32), as a cannot be essentially smaller than ∆x (in fact, the plots displayed below are generated with a = 0.02, unless it is stated otherwise).
The first result is that, for fixed values of a in Eq. (32), there is a critical value, γ cr , of the PT gain-loss coefficient, such that, at γ < γ cr , the numerical solution features a shape very close to that of the analytical solution corresponding to the ideal δ ′ -and δ-functions, while at γ > γ cr the single-peak shape of the solution transforms into a double-peak one, as shown in Fig. 1(a). In particular, it was found that γ cr (a = 0.02) ≈ 0.24. As shown below, there is the second critical value of γ above which pinned modes do not exist at all, cf. the exact result (29).
The drastic difference between the single-and double-peak modes is that the former ones are completely stable, as confirmed by systematic simulations (not shown here in detail), while all  The unstable evolution (escape) of the double-peak soliton from Fig. 1(b).
the double-peak solutions are unstable. This correlation between the shape and (in)stability of the pinned modes is not surprising: the single-and double-peak structures imply that the pinned mode is feeling, respectively, effective attraction to or repulsion from the local defect. Accordingly, in the latter case the pinned soliton is unstable against spontaneous detachment (escape) from the PT dipole, transforming itself into an ordinary freely moving NLS soliton, see an example in Fig.   2.
The transition between stable single-peak and unstable double-peak pinned modes was earlier reported in Ref. [31], which was dealing with a chain of parametrically driven damped pendulums, with a discrete soliton attached to a local defect in the chain. In that case, the instability development was different, leading not to detachment of the soliton, but rather to π-out-of-phase oscillations of the two lobes of the double-peak structure.
The results for soliton families in the present situation are summarized in Figs. 3(a) and 3(b), in the form of plots for P (k) [cf. Eq. (11)] and P (γ) . The plots also delineate the effective boundary between the stable single-peak modes and unstable double-peak ones. The curves in panel 3(b) terminate at critical points, beyond which no pinned modes are produced by the numerical solution.
Exact analytical results for solitons in the model of the PT -symmetric nonlinear coupler [12] suggest that the termination of the solution branches may be explained by a tangent (saddle-node) bifurcation, i.e., annihilation of the given branch with an additional unstable one (in the coupler model, this is the branch of PT -antisymmetric solitons). However, search for that additional branch in the present model, which is, presumably, a fully unstable one, is a challenging problem.
Note that, at γ 0.1, the numerically found total power almost does not depend on γ in Fig.   3(a), in accordance with analytical result (11). However, P grows with γ at larger values of γ. The analytical curve in Fig. 3(a) terminates at k = 0.5, as predicted by Eq. (10) for ε 0 = 1, but with the growth of γ the cutoff value of k increases.
Finally, Fig. 4 summarizes the findings in the plane of (a, γ) for fixed k = 3.0. It is clearly seen that the region of the unstable double-peak solitons is actually a relatively narrow boundary layer between the broad areas in which the stable single-peak solitons exist, or no solitons exist at all, at large values of γ. Note also that the stability area strongly expands to larger values of γ as the regularized profile (32) becomes smoother, with the increase of a. On the other hand, the stability region remains finite even for very small a. C. The system with the self-focusing uniform nonlinearity and nonlinear pinning potential Another explicit solution produced above, based on Eqs. (5), (9), and (15), pertains to the case of σ = +1, ε 0 = 0, and ε 2 > 0, when the attractive potential of the PT dipole is purely nonlinear.
In this case, stable single-peak modes, close to the aforementioned analytical solution, were found for 0 γ 0.13, while at γ > 0.13 the pinned modes are unstable, featuring a double-peak shape.
Similar to Fig. 4, the existence and stability diagram for the solitons is plotted in the plane of (a, γ)   Fig. 4 (for σ = +1, k = 3.0), but in the model with the purely nonlinear attractive potential, i.e., ε 0 = 0, and ε 2 = 2.0 (a) or ε 2 = 8.0 (b). As before, in the present case stable and unstable solitons feature the single-and double-peak shapes, respectively.
in Fig. 5 for smaller (a) and larger (b) values of ε 2 . The comparison with Fig. 4 demonstrates that, in the case of the nonlinear pinning potential, the stability area is much smaller than it was in the case of the linear attractive potential.
Next, we consider combined linear and nonlinear pinning potentials, restoring ε 0 = 1. With ε 2 > 0, stable single-peak solitons are readily found up to the respective critical value of γ (for instance, at k = 2.0 they are found at γ < 0.15 for any value of ε 2 ). With ε 2 < 0, both stable single-peak solitons and unstable double-peak ones are produced by the numerical solution. For this case, Fig. 6(a) displays P (k) curves with fixed γ = 0.1 and different negative values of ε 2 .
The curves include segments representing both the single-and double-peak modes. Further, the respective stability boundary in the plane of (ε 2 , P ) for fixed γ = 0.1 is shown in Fig. 6(b). The increase of P naturally leads to destabilization of the pinned mode, as the corresponding nonlinear repulsive potential, accounted for by ε 2 < 0, becomes stronger.
D. The self-defocusing nonlinearity (σ = −1) Another basic case corresponds to the exact solutions for the SDF nonlinearity, given by Eqs. (6), (9), and 10) with σ = −1, ε 2 = 0. In this case, the numerical solution reveals solely singlepeak modes. For small γ = 0.01, the comparison between the analytical solutions and their numerically found counterparts in displayed in Fig. 7. Further, soliton families are represented by the respective P (k) curves in Fig. 8. The corresponding analytical dependence, given by Eq.
(11), does not depend on γ, while its numerical counterparts deviate from it with the increase of γ, especially for larger ε 0 , see Fig. 8 Note that, in the framework of Eq. (4), different values of ε 0 can be transformed into ε 0 = 1 by rescaling, but regularization (32) then implies that a will be rescaled by factor ε 0 , hence larger ε implies a farther departure from the model with the ideal δ ′ -and δ-functions. This fact explains the stronger deviation of the numerically found curves from their analytically obtained counterparts in Fig. 8(b) in comparison with 8(a). The same trend is observed below in Fig. 9.
In the model with the uniform SDF nonlinearity, the analytical solution exists in the interval of the propagation constant k < k max = ε 2 0 /2, see Eq. (12). The respective numerically found existence boundaries for the single-peak solitons are displayed in Fig. 9. The observed deviation of the boundary value from k max at γ = 0 is explained by the difference of the regularized δ-function (32) from its ideal counterpart, the existence region further shrinking with the increase of γ.
Finally, the numerical results demonstrate that, in the case of σ = −1 and ε 0 = 0, solitons existing under the action of the nonlinear pinning potential with ε 2 > 0, as predicted by Eqs. (6), (9) and (15), are completely unstable. Recall that, unlike this result, in the model with the SF bulk nonlinearity (σ = +1), a small stability area was found for solitons pinned by the nonlinear potential, see Fig. 5. Finally, we have produced numerical counterparts of the simplest exact solutions given by Eq. (17) for the nonlinear dipole embedded into the linear medium (σ = 0, ε 2 = ±1). The results are summarized in the form of P (k) curves, which are displayed in Figs. 11(a) and 11(b) for ε 2 = +1 and ε 2 = −1, i.e., the SF and SDF signs of the localized nonlinearity, respectively. In the former case, the P (k) dependences obey the VK criterion, dP/dk > 0, at γ < 0.28. Accordingly, the pinned modes are stable in direct simulations in this region, and they are destroyed by an instability at γ > 0.28, when dP/dk becomes negative, see Fig. 11(a).
For the SDF sign of the localized nonlinearity, ε 2 = −1, the P (k) curves displayed in Fig. 11(b) satisfy the anti-VK criterion, dP/dk < 0. In agreement with this condition, the solitons are found to be stable in the direct simulations.

Conclusion
The objective of this work is to introduce the solvable model of the nonlinear PT -symmetric medium, in which the gain-loss combination is represented by the point-like dipole, ∼ δ ′ (x), which is embedded into the uniform Kerr-nonlinear SF (self-focusing) or SDF (self-defocusing) medium, in the combination with the linear and/or nonlinear potential pinning the wave field to the PT dipole. The host medium may be linear too. The full set of analytical solutions for pinned modes has been found for this model, along with the solutions for the system of separated PT -symmetric point-like gain and loss sites with the localized Kerr nonlinearity, embedded into the linear medium (the solution for the latter variety of the solvable system makes it possible to explicitly demonstrate the nonexistence of the PT -symmetric modes above the critical value of the gain-loss coefficient).
The analytical solutions were compared with numerical ones, obtained in the model with the ideal δ ′ -and δ-functions replaced by their Lorentzian regularization. It has been concluded that, with the increase of the gain-loss-dipole strength, γ, the shape of the pinned mode supported by the SF bulk nonlinearity gradually deviates from the analytical limit, changing from the single-peak form into the double-peak one, which coincides with the destabilization of the pinned soliton against escape. On the contrary, all the pinned modes found in the model with the SDF sign of the bulk nonlinearity are stable, featuring the single-peak shape.
The corresponding stationary equation [cf. Eq. (4)] is, in principle, solvable, although the respective algebra turns out to be cumbersome. On the other hand, it may be interesting to introduce a two-dimensional version of the system, with a gain-loss quadrupole emulating the corresponding singular expression, δ ′ (x)δ ′ (y).

IV. ACKNOWLEDGMENT
The work of T.M. was supported by the Thailand Research Fund through grant No. RMU5380005.