Solitons supported by singular modulation of the cubic nonlinearity

A model of the optical media with a spatially structured Kerr nonlinearity is introduced. The nonlinearity strength is modulated by a set of singular peaks on top of a self-focusing or defocusing uniform background. The peaks may include a repulsive or attractive linear potential too. We find that a pair of mutually symmetric peaks readily gives rise to the spontaneous symmetry breaking (SSB) of modes pinned to individual peaks, while antisymmetric pinned modes are always unstable, transforming into robust spatially odd breathers. Three- and five-peak structures support symmetric modes, with in-phase or twisted profiles, and do not give rise to asymmetric states. A stability area is found for the twisted states pinned to the triple peaks, while the corresponding in-phase modes are unstable, unless the three modulation peaks are set very close to each other, covered by a single-peak pinned mode. All patterns pinned to five peaks are unstable too. Collisions of moving solitons with the singular-modulation peak are studied too. Slowly moving solitons bounce back from the peak, while the collisions are quasi-elastic for fast solitons. In the intermediate case, the soliton is destroyed by the collision. In a special case, the condition of a resonance of the incident soliton with a trapped mode supported by the peak leads to capture of the soliton.


Introduction and the model
The great potential offered by solitons for fundamental studies and development of applications in optics and other areas of physics is commonly known [1,2,3,4,5]. However, the use of uniform media as a host material for solitons has its limitations, as they admit the existence of few species of solitary modes. In particular, the integrable one-dimensional (1D) setting based on the nonlinear Schrödinger (NLS) equation gives rise to the single stable class of solutions in the form of fundamental solitons [1] (higher-order solutions for breathers are available too [6], but they are not truly stable objects, as their binding energy is zero).
The variety of stable localized objects may be greatly expanded by using pseudopotentials induced by spatial modulation of the local nonlinearity strength [7]. In optics, nonuniform Kerr nonlinearity can be induced by an inhomogeneous distribution of nonlinearity-enhancing dopants [8]. Similar patterns may be realized in composite media assembled of different materials [9,10,11,12]. Still broader flexibility of nonlinearity landscapes is available in Bose-Einstein condensates (BECs), where the landscapes, underlain by the optical Feshbach resonance [13,14,15], can be "painted" by rapidly moving laser beams [16,17]. The magnetic Feshbach resonance can be used too, imposing stationary nonlinearity patterns by means of magnetic lattices [18,19,20].
In particular, structures in the form of a very narrow region with strong local cubic nonlinearity, embedded into a linear host setting, are modeled by the NLS equation with a delta-functional coefficient in front of the nonlinear term, g(x) ∼ δ(x) [21,22,23,24,25] (a discrete version of this setting is represented by a linear dynamical lattice with an embedded nonlinear site [26,27,28,29]). In this limit, the family of solitons pinned to the point-like inhomogeneity (nonlinear defect) is degenerate, in the sense that their norm does not depend on the propagation constant/chemical potential, in terms of optics and BEC, respectively [23]. The norm degeneracy is a distinctive feature of Townes solitons, such as those supported by uniform cubic [30,31,32] and quintic [33] self-focusing in the 2D and 1D geometry, respectively. The fact that solitons of the Townes type are always made unstable by the occurrence of the critical collapse in the same setting makes it necessary to replace the delta-functional nonlinearity profile by more realistic singular ones, approximated by the cusp, with α > 0 [34]. The analysis has produced stable solitons pinned to this nonlinearity peak in the range of 0 < α < 1. In the limit of α = 1, the solitons' amplitude vanishes, and they do not exists at α > 1, when the singularity is too strong ( g(x)dx diverges around x = 0; actually, solitons cannot be pinned to the nonlinearity peak with α ≥ 1 as they are destroyed by the collapse in that case). It is worthy to note that it was also demonstrated in [34] that the singular modulation profile with α < 1 may emulate the spatially uniform attractive cubic nonlinearity in a sub-1D space, with effective fractal dimension A natural continuation of the analysis, which is the objective of the present work, is to study configurations pinned to sets of several nonlinearity peaks. Such patterns can be created in the experiment by means of the same techniques that were proposed for building a single peak. In particular, the possibility of the spontaneous symmetry breaking (SSB) in patterns attached to a pair of identical peaks, which is known in other settings [26,27,28,29], [22,35], opens a way to use the double peaks for power-controlled switching in optical circuitry. Further, Table 1: Different sets of nonlinearity-modulation peaks representing the coefficient in front of the cubic term in Eq. (2).
1 peak sets of several peaks may be considered as a transition to a periodic lattice of nonlinear defects, which is an object of obvious interest too [36,37,38,39]. Thus, we consider the model based on the NLS equation with an x-dependent self-focusing coefficient, g(x) > 0: This equation is written in terms of optics, with z and x being the propagation distance and transverse coordinate in a planar waveguide, and σ representing the nonlinearity coefficient of the uniform background into which the nonlinearity peaks are embedded. In particular, σ > 0 implies competition of the selffocusing peaks with the self-defocusing background, cf. work [12]. The linearpotential term ∼ β takes into account the fact that material perturbations, which induce the local modulation of the nonlinearity, may also give rise to similar local changes of the refractive index. In the application to a cigar-shaped trapping configuration in BEC with longitudinal coordinate x, the evolutional variable, z, in Eq. (2) must be replaced by time t, and the term ∼ β accounts for the local linear potential that may be induced by the same tightly focused laser beams which modify the local strength of the nonlinearity via the Feshbach resonance.
The singular nonlinearity-modulation patterns considered below are specified in Table 1, including two, three, or five peaks. For the completeness' sake, the single peak, which was introduced in [34], is included too. In these patterns, ∆ determines the separation between adjacent peaks (for the double peak, the separation is 2∆), and α is the singularity power in Eq. (1). Coefficients in front of the singular-modulation functions, which are defined in Table 1, are set to be 1 by means of obvious rescaling. Further, the remaining scaling invariance may be used to fix |σ| = 1 in Eq. (2), unless we choose σ = 0 (no uniform-nonlinearity background).
It is relevant to stress the difference of the superposition of two singular peaks, defined as per the second line of Table 1, from the superposition of regularized δ-functions, represented by Gaussians with width a: while the singular peaks remain discernible at any separation 2∆ between them, the two peaks in the superposition of the Gaussians merge into one at a/ (2∆) > 1/ √ 2 [22]. The rest of the paper is organized as follows. Methods used in the analysis are briefly summarized in Section II. These include a quasi-analytical variational approximation (VA) for the mode pinned to a single peak, and a numerical scheme elaborated for the identification of stationary patterns and their stability. The VA presented here generalizes that from [34], which did not include the linear-potential component of the local peak. Basic results for the SSB in patterns pinned to a pair of peaks are reported in Section III, and findings for the set of three and five peaks (the latter being a prototype of the lattice of nonlinear defects) are presented in Section IV. In Section V, a dynamical situation is considered, viz., collisions of free solitons with a single nonlinearity peak (collisions with defects and interfaces of various types were studied before [40,41,42,43,44,45,46], but not for the present one). The paper is concluded by Section VI.
2 The framework of the analysis 2.1 The variational approximation for a single peak Stationary solutions to Eq. (2) with propagation constant k are looked for in the usual form, u (x, z) = e ikz w(x), with real function w satisfying equation Numerical solutions of Eq. (3) with zero boundary conditions at |x| → ∞ were constructed by means of the Newton's method, while simulations of the full NLS equation (2) were carried out by means of an explicit Runge-Kutta algorithm. The stationary solutions are characterized by the total power (norm), The use of the Newton's method is contingent on a possibility to choose a relevant initial guess. It was provided by attaching to each nonlinearity peak a localized mode predicted for the isolated defect by the VA. To this end, the simplest Gaussian ansatz, is substituted in the Lagrangian of Eq. (3) with g(x) taken as per the top line in Table 1, The calculation of integrals in Eq. (6) gives rise to the corresponding effective Lagrangian, where Γ is the Euler's Gamma-function. For given k, parameters of ansatz (5) are determined by the corresponding Euler-Lagrange equations, The use of the VA-predicted profiles for modes pinned to individual peaks as inputs always provide for convergence of the Newton's method, applied to settings with one or several peaks, even in cases when (depending on values of control parameters) the so generated numerical solution might be conspicuously different from the input.

The stability test
Stability of stationary patterns was checked by means of the linearization with respect to small perturbations, taking the perturbed solution to Eq. (2) as where ± (x) are components of the perturbation eigenmode, and λ is the instability growth rate. Straightforward linearization reduces the search for eigenvalues λ to the calculation of the spectrum of a linear matrix operator, where C ≡ [σ − g(x)] w 2 (x), and the Sturm-Liouville operator iŝ [in fact, the stability analysis is presented below for solutions with real w(x)].
Dealing with the conservative systems, stable stationary solutions are identified as those for which all eigenvalues λ have zero real parts. The so predicted stability was checked in direct simulations of the evolution of perturbed solutions of Eq. (2).

The system with a symmetric pair of singularmodulation peaks
As said above, the basic issue for the setting with the double peak, which corresponds to Eq. (2) with g(x) taken as per the second line in Table 1, is the possibility of the SSB between solitons pinned to the two peaks, with separation 2∆ between them. The asymmetry is characterized by the normalized intensity-contrast parameter, The symmetry breaking indeed happens with the increase of propagation constant k, as shown in Figs. 1 and 2 for β = 0 [no linear potential in Eq.  Table  1 with ∆ = 0.244. The three panels correspond to different strengths of the background nonlinearity: σ = −1 (a), σ = 0 (b), and σ = +1 (c). Green and red branches represent stable and unstable solutions, respectively. Red segments of the symmetric solutions (Θ = 0) at small k correspond to very weak(actually, formal) instability of modes with the convex shape, see below.
(3)]. These figures clearly demonstrate that the SSB in the present system corresponds to the supercritical pitchfork bifurcation [47], which destabilizes the symmetric states, replacing them by stable asymmetric ones. A qualitatively similar result was obtained for the nonlinearity modulation format represented by a pair of symmetric Gaussians (regularized δ-functions) in [22], which confirms the generic character of this SSB scenario in systems with the double-peak modulation of the local self-focusing nonlinearity. This conclusion is further corroborated by the fact that the presence of the uniform background nonlinearity of either sign does not essentially effect the SSB scenario, although it naturally shifts the total power at the bifurcation point to larger values in the case of the defocusing background, as seen in Fig. 2. Generic examples of stationary symmetric and asymmetric modes are displayed in Figs. 3(c) and 3(a), respectively. The conclusion about the destabilization of the symmetric states at the SSB point, which was produced by the calculation of stability eigenvalues by means of Eq. (9), is confirmed by direct simulations, demonstrating that unstable symmetric states spontaneously transform into their stable asymmetric counterparts, as shown in Fig. 4.
Data concerning the location of the SSB (i.e., bifurcation) points are collected in Fig. 5, which provides a comprehensive map of the parameter space of the double-peak setting. It is seen that, as mentioned above, the dependence of the critical values of the propagation constant (k) and total power (P ) on the strength of the background nonlinearity, σ, is quite weak, while the dependence on the modulation's singularity degree, α [see Table 1], is essential. Naturally, k → ∞ and P → ∞ is observed Fig. 5 at ∆ → 0, as the SSB may happen if the width of the mode pinned to an individual peak, which scales as k −1/2 , has  Table 1, α = 0.5, β = 0 and different values of σ in Eq. (3). Here, and in Fig. 6 below, symbols in the notation box mean: "u." -unstable symmetric solutions, "symm. s." -stable symmetric ones, and "asymm. s." -stable asymmetric states.  the same order of magnitude as ∆.
Unlike the uniform background nonlinearity, the presence of the linear component of the modulation peaks, represented by coefficient β in Eq. (2), produces a strong effect on the SSB scenario. Namely, β > 0 (which correspond to the repulsive linear potential) leads to partial destabilization of the asymmetric states, as shown in Fig. 6. In this case, the pair of the singular-modulation peaks may support two different asymmetric states: one with a relatively small difference between amplitudes of the solution at the two peaks, and a fully asymmetric state, strongly pinned to one of the peaks. The repulsive linear potential destabilizes the weakly asymmetric state, converting it into an asymmetric breather, as shown in Fig. 7(a). On the other hand, the unstable symmetric states are transformed into strongly asymmetric modes pinned to one of the peaks, as shown in Fig. 7(b). This finding may be explained by the fact that the increase of the amplitude of the pinned soliton helps the nonlinear attractive potential to overcome the repulsion induced by the linear potential. On the other hand, symmetric states, unlike the asymmetric ones, are not destabilized by the linear component of the modulation peaks.
In the limit case when the localized modes pinned to adjacent nonlinearity peaks are strongly overlapped, the two modes merge into a nearly flat-top one, with the small deviation of the flatness corresponding to slightly convex or concave shapes, as shown in Figs. 3(d) and 3(e), respectively. Further analysis demonstrates that the concave shapes are completely stable, while the computation of the eigenvalues of operator (9) for convex ones (which appear when the two individual modes become completely overlapping) features a very weak instability, that corresponds to short red segments near the left edge of all the panels in Fig. 1. The eigenvalues that characterize the instability are limited to be 10 −4 . Actually, direct simulations do not reveal any tangible instabil-    Figure 6: The same as in Fig. 2, but for α = 0.5, σ = −1, ∆ = 0.244, and different values of β. It is seen that, in the presence of β > 0, weakly asymmetric stationary states become partially unstable (see, e.g., the solution families for β = 0.05, where a gap between chains of green circles and stars corresponds to the instability). The instability leads to transformation of the weakly asymmetric modes into asymmetric breathers, as shown in Fig. 7(a) . ity of the convex profiles, in accordance with the fact that the corresponding eigenvalues are extremely small. Boundaries between the nearly flat-top convex and concave profiles in the plane of (k, ∆) for different values of α are displayed in Fig. 8. It is also relevant to mention that, in accordance with the results of [34], where it was demonstrated that modes pinned to a single modulation peak, competing with the self-defocusing uniform background (σ = +1), exist if their total power exceeds a certain threshold (minimum) value, the same effect was found here for the symmetric modes, the threshold being virtually exactly equal to twice the one found in [34] (this feature is not shown in detail here). Lastly, the model with two singular-modulation peaks admits solutions for antisymmetric (twisted) pinned modes too. However, both the computation of the stability eigenvalues and direct simulations reveal that all the twisted states are unstable. As shown in Fig. 9, the instability transforms them into robust breathers, which preserve the twisted structure.

The system with three and five singular-modulation peaks
The transition from the double singular-modulation peaks to a triple one (which corresponds to the third line in Table 1) is of obvious interest, as it makes it possible to understand respective changes in the variety of pinned modes with different symmetries. In this case, asymmetric modes (stable or unstable ones) have not been found, while symmetric states exist in two varieties, viz., with in-phase and out-of-phase (alias twisted) shapes, which are shown in Fig. 10.  The computation of the corresponding eigenvalues, as well as direct simulations, demonstrates that the in-phase modes with k = 1.2 have their stability regions at ∆ < 0.66 and ∆ < 0.61 for σ = +1 and −1, respectively, as shown in Fig.  11(a) (the stability regions have similar shapes at other values of k). In this region, the modes are strongly overlapped, seeming as a single local-power peak, see Fig. 11(b).
Outside of the above-mentioned regions, the in-phase states are always unstable (unless separation ∆ between the peaks is so large in comparison to the width of local modes pinned to individual peaks that they practically do not interact), while twisted modes have a nontrivial stability area in the parameter plane of (∆, α), as shown in Fig. 12. It is worthy to note a conspicuous dependence of the stability boundaries on strength σ of the uniform nonlinearity background in Figs. 11 and 12, unlike a much weaker dependence on σ in the case of the double peak, cf. Fig. 5.
In direct simulations, all unstable modes pinned to the triple peak spontaneously transform into patterns featuring either a state pinned to the central peak, while the edge modes disappear, or two narrow modes pinned to the edge peaks, while the central one disappears, see examples in Fig. 13. In the latter case, the emerging pattern seems to be stable because the narrow modes of which it is built virtually do not interact with each other.
As mentioned above, the consideration of the set of five modulation peaks, which corresponds to the fourth line in Table 1, is relevant too, as it is a prototype of a periodic lattice of nonlinear defects [7]. In this case, no states with broken symmetry could be produced, similar to the above-mentioned situation with the triple peak. As concerns symmetric states, all of them, of both in-phase and twisted types (see examples in Fig. 14), have been found to be unstable,    Table 1 (the twisted modes are unstable on the left-hand side of the continuous boundary lines). In this figure, β = 0 and k = 1.2 are fixed, while the strength of the background uniform nonlinearity σ takes both values σ = +1 and −1 (the self-defocusing and focusing background, respectively). Beyond the dashed boundaries, inphase modes cease being unstable because the interaction between different modes pinned to individual peaks becomes negligible. i.e., the nonlinear lattice of the singular-modulation peaks, unlike the double and triple peaks, is not able to support stable states which attach local modes to all the peaks [unless separation ∆ between adjacent ones is so large in comparison with widths of the modes pinned to individual peaks that they cease to interact, or so small that a single narrow mode can cover all of them, cf. Fig. 11(b)]. In direct simulations (not shown here in detail), unstable states supported by the five-peak set transform themselves into "rarefied" patterns formed by three narrow modes pinned to the central peak and two edge ones, while intermediate peaks (the second and fourth ones) remain empty, quite similar to what is shown in Fig. 13 for unstable states in the case of the triple peak. The emergent pattern seems stable because individual pinned modes do not interact. This instability and its development scenario resemble the modulational instability of extended states pinned to periodic potentials in the NLS equation with the self-attractive nonlinearity [48,49,7].

Collisions of free solitons with the nonlinearity-modulation peak
Equation Eq. (2) with the self-focusing uniform nonlinearity (σ = −1) obviously supports free NLS solitons far from the location of the peaks: where P is the soliton's total power, and v is its velocity [in terms of the spatial optical solitons, it is the tilt of beam in the (x, z) plane]. In this case, it is natural to consider collisions of freely moving solitons with peaks. Here, we limit the consideration to the basic case of the collision with a single peak, as this case was not studied previously either.  Figure 14: Examples of (a) fully in-phase and (b) combined in-phase-twisted unstable states attached to the set of five singular-nonlinearity-modulation peaks, which corresponds to the fourth line in Table 1. The parameters are: σ = −1, α = 0.5, ∆ = 1.95, k = 2 and β = 0.
An essential result of systematic simulations is that, in most cases, the incident soliton bounces back from the singular-modulation peak if its velocity (tilt) is smaller than a certain critical value, v < v cr . The soliton is destroyed by the collision into a radiation wave train at v v cr , and, finally, the collision is quasi-elastic at v v cr , see typical examples in Figs. 15. The dependence of the trapped, reflected, and transferred shares of the total power on the incidence velocity is shown in Fig. 16. The dependence of v cr on power P of the impinging soliton for different values of the singular-modulation power α is shown in Fig.  17.
The fact that the incident soliton readily bounces back from the singularmodulation peak, which represents an attractive nonlinear defect, is remarkable by itself, as rebound is normally expected from repulsive defects. Previously, rebound of solitons impinging on an (attractive) linear potential well was reported in [50,51].
Capture of the incident soliton by the peak was observed too, but only under a specific condition, which may be considered as a resonant one: a soliton with total power P gets captured, as shown in Fig. 15(a), if its propagation constant in the quiescent state, i.e., k sol (v = 0) = P 2 /8, see Eq. (12), is close to the value of k for the mode with the same value of P trapped by the self-focusing peak, while the velocity of the incoming soliton is relatively small. This special case is illustrated by dependences of the of the trapped, reflected, and transferred shares of the total power on v shown in Fig. 16(c). The resonant situation can be approximately predicted by equating k sol (v = 0) to the propagation constant predicted, for the same value of P , by the VA, i.e., by the Euler-Lagrange equations following from the effective Lagrangian (7) (not shown here in detail).  Figure 15: Collisions of an incident free soliton (12) with the nonlinearitymodulation peak, corresponding to the first line in Table 1 Figure 16: Shares of the total power of the incident free soliton with (a) P sol = 2, which corresponds to the propagation constant k = 0.5 of the free soliton with zero velocity (tilt) and (b) P sol = 4, which corresponds to propagation constant k = 2 of the free soliton with zero velocity (tilt), that are trapped (the red curve) , reflected (the green curve) and keep the original direction of the motion ("transmitted power", depicted by the blue curve), versus the incidence velocity, v. The solitons impinge on the nonlinearity peak with σ = −1, β = 0, and α = 0.4. The respective critical velocity, at which the soliton is completely destroyed by the collision is v cr ≈ 6 for the soliton with P sol = 2, and v cr ≈ 13 for the one with P sol = 4. Panel (c) illustrates the special case of the resonant capture of the incident soliton, an example of which is displayed in Fig. 15(a). In this case, the parameters are σ = −1, β = 0, α = 0.6, and P sol = 8.  Figure 17: v cr as a function of power P of the moving soliton, and singularmodulation parameter α, at σ = −1 and β = 0.

Conclusions
We have introduced the 1D model with the cubic self-focusing nonlinearity locally modulated by a set of singular profiles, while the background uniform nonlinearity may be focusing or defocusing. The possibility of the modulation peaks including the linear-potential components is considered too. Numerical analysis has demonstrated that the setting with a pair of symmetric peaks readily gives rise to the SSB (spontaneous symmetry breaking) of the modes pinned to the individual peaks, via the supercritical bifurcation. The antisymmetric states were found too, turning out to be unstable and transforming into antisymmetric breathers. Sets of three and five peaks do not give rise to asymmetric states, rather supporting symmetric ones, with in-phase or twisted profiles. The twisted states pinned to the triple peaks have their stability area, while the inphase modes are unstable, except for the strongly overlapping ones. The sets of five peaks fails to support any stable nontrivial modes.