Solitons under spatially localized cubic-quintic-septimal nonlinearities

TL;DR

This paper investigates the stability of solitons in a nonlinear Schrödinger equation with spatially localized cubic, quintic, and septimal nonlinearities, revealing conditions for stable high-power light beams in optical waveguides.

Contribution

It introduces a model with combined cubic, quintic, and septimal nonlinearities confined in space, providing analytical and numerical stability analysis, and highlights the potential for high-power beam guidance.

Findings

Stability regions align with the Vakhitov-Kolokolov criterion.

Tight confinement increases maximum stable soliton power with defocusing quintic.

Analytical solutions are derived for the limit case of delta-function confinement.

Abstract

We explore stability regions for solitons in the nonlinear Schrodinger equation with a spatially confined region carrying a combination of self-focusing cubic and septimal terms, with a quintic one of either focusing or defocusing sign. This setting can be implemented in optical waveguides based on colloids of nanoparticles. The solitons stability is identified by solving linearized equations for small perturbations, and is found to fully comply with the Vakhitov-Kolokolov criterion. In the limit case of tight confinement of the nonlinearity, results are obtained in an analytical form, approximating the confinement profile by a delta-function. It is found that the confinement greatly increases the largest total power of stable solitons, in the case when the quintic term is defocusing, which suggests a possibility to create tightly confined high-power light beams guided by the spatial modulation of the local nonlinearity strength.

Equations

$$i\frac{\partial \Psi }{\partial z}+\frac{1}{2}\frac{{\partial }^2\Psi }{\partial x^2}+[G_3(x)P|\Psi {|}^2+G_5(x)P^2|\Psi {|}^4+G_7(x)P^3|\Psi {|}^6]\Psi =0,$$

$${\int }_{-\infty }^{+\infty }|\Psi (x,z){|}^{2}dx=1$$

$$G_n(x)=g_n\exp \left(-b^2{x}^2/2\right),$$

$$z$$

$$P$$

$$g_3=g_7=1,$$

$$k\psi =\frac{1}{2}\frac{d^2\psi }{d{x}^2}+\exp \left(-\frac{1}{2}{b}^2{x}^2\right)P{\psi }^3\left(1+g_{5}P{\psi }^2+P^2{\psi }^4\right).$$

$$\left(i\frac{\partial }{\partial z}-k+\frac{1}{2}\frac{{\partial }^2}{\partial x^2}\right)\delta \psi +\exp \left(-\frac{1}{2}{b}^2{x}^2\right)P{\psi }^2(x)\times [(\delta {\psi }^{\ast }+2\delta \psi )+$$

$$g_{5}P{\psi }^2(x)(2\delta {\psi }^{\ast }+3\delta \psi )+P^2{\psi }^4(x)(3\delta {\psi }^{\ast }+4\delta \psi )]=0,$$

$$\frac{1}{2}\frac{d^{2}u}{d{x}^2}+\exp \left(-\frac{1}{2}{b}^2{x}^2\right)P{\psi }^2(x)[(2u+v)$$

$$+g_{5}P{\psi }^2(x)(3u+2v)+P^2{\psi }^4(x)(4u+3v)]-ku=\lambda u,$$

$$\frac{1}{2}\frac{d^{2}v}{d{x}^2}+\exp \left(-\frac{1}{2}{b}^2{x}^2\right)P{\psi }^2(x)[(2v+u)$$

$$+g_{5}P{\psi }^2(x)(3v+2u)+P^2{\psi }^4(x)(4v+3u)]-kv=-\lambda v.$$

$${\psi }_{k\to \infty }(x)\approx \frac{k^{1/4}}{{\left[\cosh \left(3\sqrt{2k}x\right)\right]}^{1/3}},$$

$$k\approx 15.2\cdot P^{-6}$$

$$\exp \left(-b^2{x}^2/2\right)\approx \sqrt{2\pi }{b}^{-1}\delta (x).$$

$$k\psi =\frac{1}{2}\frac{d^2\psi }{d{x}^2}+\frac{\sqrt{2\pi }}{b}\delta (x)P{\psi }^3\left(1+g_{5}P{\psi }^2+P^2{\psi }^4\right).$$

$$\psi (x)={\left(2k\right)}^{1/4}\exp \left(-\sqrt{2k}|x|\right).$$

$$\frac{1}{2}\frac{d\psi }{dx}{|}_{x=0}=-\frac{\sqrt{2\pi }}{b}P{\psi }^3\left(1+g_{5}P{\psi }^2+P^2{\psi }^4\right){|}_{x=0}.$$

$$2P^{3}k+g_5{P}^2\sqrt{2k}+\left(P-\frac{b}{\sqrt{2\pi }}\right)=0.$$

$$P_{\max }\left(g_5>0;b\right)=b/\sqrt{2\pi }.$$

$$P_{\min }(g_5$$

$$=$$

Full text

1 Introduction

The importance of spatial solitons, which are maintained by the balance between diffraction and self-focusing, in modern photonics is commonly known [1]-[6]. In particular, the ubiquitous Kerr self-focusing nonlinearity is modeled by the cubic nonlinear Schrödinger equation (NLSE), which predicts stable solitons in optical fibers and planar waveguides, as well as in many other systems [2]. On the other hand, more complex effects, such as the generation of higher-order harmonics [9], filamentation [10], saturation of the Kerr nonlinearity [11], modulational instability in metamaterials [12], creation of quasi-stable optical beams with embedded vorticity [13], and other nonlinear phenomena make it necessary to add higher-order nonlinearities to the NLSE. In particular, the recently reported stable propagation of (2+1)-dimensional fundamental spatial solitons in bulk dielectric media is explained by the action of the cubic-quintic (CQ) [14] or quintic-septimal [15] nonlinearities, in which the sign of the higher term is defocusing. It was also predicted that the CQ nonlinearity is sufficient to support stable two- [16] and three- [17] dimensional (2D and 3D) solitons with embedded vorticity.

Recently, it has been demonstrated that various combinations of optical nonlinearities, including CQ and septimal terms of either sign, can be engineered in colloids composed of metallic nanoparticles [18]. This finding has suggested a detailed analysis of 1D solitons under the action of the combined cubic-quintic-septimal (CQS) nonlinearity [19], the most essential issue being stability of such soliton families for different combinations of signs in front of the CQS terms. In particular, a noteworthy result is that, in spite of the well-known fact that the focusing quintic and septimal terms give rise, respectively, to the critical and supercritical collapse in the 1D NLSE [20, 21], the 1D equation with the focusing signs of all the nonlinear terms in the CQS combination maintains stable soliton families (similar to the previously known fact that the combination of self-attractive cubic and quintic terms also generates stable solitons [22]).

Another direction of the work with nonlinearity, which has also drawn much attention in the course of the last decade, is the use of spatially modulated nonlinear terms for the creation and stabilization of various species of solitons [5]. In particular, a stepwise radial change of the local strength of the cubic self-focusing is sufficient to stabilize Townes solitons in the 2D space [23]. It has also been recently observed that the spatially modulated nonlinearity could also be contributed to earlier onset of Faraday and resonant waves [7, 8]. It was also demonstrated that stable 2D solitons with embedded higher-order vorticity, up to $S=6$, can be supported by localized cubic self-focusing in a combination with the harmonic-oscillator trapping potential [24], while previously no stable states with $S\geq 2$ could be constructed using local nonlinearity and trapping potentials [3, 5].

The above-mentioned findings suggest to seek for possibilities to predict stable solitons supported by the CQS nonlinearity subject to the spatial modulation, which is the objective of the present work. In effectively planar waveguides built of the above-mentioned colloidal materials, which provide the CQS nonlinearity for the creation of spatial solitons, the modulation can be induced by making the waveguide’s thickness a function of the propagation distance. Higher-order nonlinearities may also be provided by resonant dopants (see, e.g., Ref. [25]), in which case the spatial modulation is naturally imposed by inhomogeneity of the dopant density along the propagation distance. We focus on the most interesting “sandwich” setting, when the lowest and highest cubic and septimal terms are self-focusing, while the intermediate quintic one may have either sign. A challenge is to secure stability of the solitons when, in particular, the competition between focusing and defocusing (quintic) terms may support stable solitons, in spite of the possibility of the onset of the supercritical collapse, driven by the septimal self-focusing term. The model, based on an appropriate version of the NLSE, is introduced in Section II. Accurate results of a systematic numerical analysis, focused on the solitons’ stability, are reported in Section III. For the limiting case of tight confinement of the nonlinearity, approximate analytical results, which are found to be in good agreement with the numerical counterparts, and thus help to explain the numerical findings, are presented in Section IV. The paper is concluded in Section V.

2 The model

The 1D NLSE for amplitude $\Psi\left(x,z\right)$ of the electromagnetic wave with spatially-localized CQS nonlinear terms is written in the scaled form for the spatial domain as

[TABLE]

with normalization condition

[TABLE]

and total power $P$. Here $\ z$ and $x$ are the propagation distance and transverse coordinate, and $x$-dependent coefficients in front of the cubic, quintic, and septimal ($n=3,5,7$, respectively) terms are defined as

[TABLE]

where $b$ determines the size of the nonlinearity-bearing region, $|x|~{}\lesssim 1/b$. Further, being interested, as said above, in the case when both the cubic and septimal terms are self-focusing, we use the possibility to rescale the variables and parameters in Eqs. (1) and (3), as

[TABLE]

so as to fix the accordingly transformed cubic and septimal coefficients (primes appearing in Eq. (4) are dropped below),

[TABLE]

while $b$, $g_{5}$ and $P$ are kept as irreducible control parameters, $g_{5}>0$ and $<0$ corresponding, severally, to the focusing and defocusing quintic terms. Of course, the normalization conditions (5) may be chosen in a different from, if that may be more convenient for some purposes.

The dynamical model based on Eq. (1) with $g_{5}>0$ can be also derived as an approximate form of the Gross-Pitaevskii equation (GPE) for a relatively dense atomic Bose-Einstein condensate with attractive inter-atomic interactions. Starting from the 3D GPE with the cubic self-attractive nonlinearity, the reduction of the dimension from 3D to 1D for a condensate tightly trapped in a cigar-shaped potential leads to the NLSE for the 1D mean-field wave function $\Psi\left(x,t\right)$ with nonpolynomial nonlinearity $\sim\left[1-(3/2)g|\Psi|^{2}\right](1-g|\Psi|^{2})^{-1/2}\Psi,$ where $g>0$ is an effective self-attraction coefficient [26]. The truncated expansion of this expression in powers of $|\Psi|^{2}$ up to the third order leads to the NLSE in the form of Eq. (1), with $G_{5}\sim g^{2}>0$, and $z$ replaced by scaled time $t$. The spatial modulation of the nonlinearity coefficient, $g(x)$, can be imposed by means of the experimentally available method based on the Feshbach resonance controlled by a spatially inhomogeneous optical [27] or magnetic [24] field.

Stationary states with propagation constant $k$ are looked for as solutions to Eq. (1) in the form of $\Psi(x,z)=\psi(x)e^{ikz}$, with real function $\psi(x)$ satisfying the ordinary differential equation,

[TABLE]

Stability of the stationary states is addressed by taking perturbed solutions as $\Psi(x,t)=\left[\psi(x)+\delta\psi(x,z)\right]\exp\left(ikz\right)$, where small perturbation $\delta\psi(x,t)$ evolves according to the linearized equation,

[TABLE]

with $\ast$ standing for the complex conjugate. For perturbation eigenfunctions sought for in the usual form [28, 29, 30, 31], $\delta\psi(x,z)=u(x)\exp\left(i\lambda z\right)+v^{\ast}(x)\exp\left(-i\lambda^{\ast}z\right)$, the corresponding eigenvalue problem for $\lambda$ amounts to the system of coupled linear equations:

[TABLE]

As usual, the stability condition is that all eigenvalues $\lambda$ (eigenfrequencies of the small-perturbation modes) must be purely real.

Stationary equation (6) was solved numerically by means of a relaxation algorithm, which is outlined in Ref. [32], starting with spatial stepsize $\Delta x=0.01$, and reducing it for solutions with a shrinking size. For the verification of the results, stationary solutions were also produced, for $g_{5}>0$ and $b>2$, with the help of a shooting Runge-Kutta algorithm, as described in Ref. [33]. The time evolution governed by Eq. (1) was simulated using a Crank-Nicolson algorithm [34], with spatial stepsizes similar to the above-mentioned ones, employed for the realization of the relaxation method, and temporal step $\Delta t=0.001$.

The linear system (8) was solved, using $\psi^{2}(x)$ produced by the relaxation method, and putting the system on a spatial grid composed of up to $1000$ points. For carrying out the solution procedure, the corresponding matrix was diagonalized using a LAPACK routine [35], that eventually provided all the eigenvalues.

3 Numerical results

Typical examples of numerically found solitons, along with their counterparts predicted by the analytical approximation [see Eqs. (13) and (15) below], are displayed in Fig. 1. The shape of the solitons supported by the confined nonlinearity is close to that of peakons, which are usually produced by models with nonlinear dispersion, such as the Camassa-Holm equation [36] and the continual limit of the Salerno model [37] with competing on-site and inter-site cubic terms [38].

Further, in Fig. 2 we present characteristics of numerically found stationary soliton families, provided by the dependence of the propagation constant, $k$, on the total power, $P$, as obtained from the numerical solution of Eq. (6) for $g_{5}=-1$ and $g_{5}=+1$, i.e., defocusing and focusing signs of the quintic term, with different values of parameter $b$ which determines the inverse width of the nonlinearity-localization region in Eq. (3). The stability of the soliton families is also shown in Fig. 2.

At $g_{5}<0$ (the defocusing sign of the quintic term) and each $b$, the stable solitons exist in the interval of $P_{\min}\left(g_{5}<0;b\right)<P<P_{\max}\left(g_{5}<0;b\right)$, with $P_{\min}\left(g_{5}<0;b\right)$ and $P_{\max}\left(g_{5}<0;b\right)$ given (in an approximate form) by Eq. (18) presented below. These smallest and largest values of the total power are designated by bold dots in Fig. 2.

In the limit of $k\rightarrow\infty$, the soliton solution, subject to normalization condition (2), becomes very narrow (hence it does not depend on the value of $b$), with the shape dominated by the septimal term:

[TABLE]

where the limit form of the $k(P)$ relation is

[TABLE]

[numerical factors in Eqs. (9) (where the factor is very close to $1$) and (10) are determined by the value of coefficient $\int_{-\infty}^{+\infty}\left(\mathrm{sech~{}}x\right)^{2/3}dx\approx 4.2$]. The range of (very large) values of $k$ in which the asymptotic solution, given by Eqs. (9) and (10), is valid, is not shown in Fig. 2 because the solutions are completely unstable in that case, as Eq. (10) does not meet the Vakhitov-Kolokolov (VK) criterion, which is discussed in detail below (the instability is also corroborated by numerical results).

It is worthy to note that the self-defocusing sign of the quintic term, $g_{5}<0$, provides for the expansion of the existence and stability limit to much greater values of the total power, $P$. While this result is quite natural in the case of the CQS nonlinearity, it is relevant to stress that a formally similar combination of the focusing quintic and defocusing cubic nonlinearities, in the absence of the septimal term, gives rise to completely unstable solitons [22]. In the present case, the stabilization is provided by the presence of the focusing term with the lowest (cubic) nonlinearity, which is absent in the above-mentioned cubic-quintic setting. As for the situation with the focusing quintic term, $g_{5}>0$, which is displayed in the Fig. 2(b), it is relevant to stress that each $k(P)$ curve features a stability segment.

Furthermore, it is noteworthy too that, taking larger $b$, i.e., a narrower nonlinearity-localization region in Eq. (3), which is the most essential feature of the present model, also helps to extend the stability region of $P$. This effect is quite conspicuous, in the Fig. 2(a), for $g_{5}<0$, while for $g_{5}>0$ the stability interval expands, in the Fig. 2(b), up to some critical value of $b$, which is followed by the transition to completely unstable soliton families. These observations suggests to consider the limit form of the model with the modulation represented by the delta-function, instead of the Gaussian in Eq. (3), which admits an exact analytical solution, as shown in the next section.

In the typical $k(P)$ plots displayed in Fig. 2, each curve has a single boundary between unstable and stable (if any) segments. In a more special case, occurring at $g_{5}>0$ (the fully focusing nonlinearity) near the transition between partly stable and fully unstable soliton families, $k(P)$ curves feature stable segments located between two unstable ones, i.e., in a narrow interval of $P_{\min}<P<P_{\max}$, as can be seen in Fig. 3. In particular, the stability regions displayed in the figure resemble results reported in Ref. [39] for a liquid-gas phase transition.

The largest total power, $P_{\max}$, up to which stable solitons exist, is an obviously important characteristic of the model with finite $b$ and $g_{5}<0$, as suggested by the Fig. 2(a). This dependence is displayed in Fig. 4, along with the analytical prediction obtained in the delta-functional limit.

Results of the systematic numerical analysis are summarized in Fig. 5, which displays $P_{\max}$ as a function of both control parameters, $g_{5}$ and $b$, in a broad area, $\left(0\leq b\leq 10\right)\times\left(-1.5\leq g_{5}\leq+1.5\right)$. This map shows, in the general form, the trend which was stressed above in the connection to Figs. 2 and 4, viz., that, at $g_{5}\leq 0$, the increase of $b$, i.e., narrowing of the spatial localization of the nonlinearity, leads to strong increase of the largest total power of the stable solitons. This effect may find various applications, such as strong confinement of powerful light beams in a narrow guiding channel by means of the spatial modulation of the local nonlinearity, in the context of the general topic of nonlinear light guiding [5].

In addition to the accurate computation of the eigenvalues based on Eq. (8), the stability of the soliton families can be determined by dint of the VK criterion [40, 20], which gives the necessary stability condition in the form of $dk/dN>0$ (in the general case, the VK criterion is not sufficient for the stability, although in relatively simple settings, such as the present ones, it may be sufficient). As is clearly seen in Figs. 2 and 3, in the present model the VK criterion is, indeed, not only necessary but also sufficient for the stability, as it precisely predicts the stability areas, identified by the condition of $\mathrm{Im}\{\lambda\}=0$ [see Eq. (8)].

The so predicted stability/instability of the solitons was tested by direct simulations of their propagation, as shown in Figs. 6 and 7. In particular, at $g_{5}=1.5$, $b=2.7$, and $k=0.04$, $P=1.038$, the soliton with broad tails [their size may be estimated as $\sim k^{-1/2}$, see Eq. (9)], which, according to Fig. 3, is unstable, spontaneously transforms into a robust breather with a small amplitude of internal vibrations, as shown in Figs. 6(a) and 7(a). On the other hand, Figs. 6(c) and 7(c) demonstrate that an unstable soliton with essentially shorter tails, corresponding to $k=0.7$, $P=1.049$, suffers the collapse. Lastly, the predicted stability of a typical soliton with $k=0.25$, $P=1.045$ is confirmed by the simulations displayed in Figs. 6(b) and 7(b).

4 The analytical solution in the limit of the delta-functional

modulation profile

As mentioned above, the model produces noteworthy results, such as the possibility to support the transmission of tightly confined high-power beams, at large values of $b$. In this case, the narrow Gaussian in Eq. (3) may be approximated by the delta-function:

[TABLE]

Then, Eq. (6) is replaced by

[TABLE]

At $x\neq 0$, Eq. (12) amounts to a linear equation, $\psi^{\prime\prime}=-2\mu\psi$, whose solution, satisfying normalization condition (2), is

[TABLE]

The integration of Eq. (12) in an infinitesimal vicinity of $x=0$ gives rise to the following boundary condition at $x=0$, assuming that $\psi(x)$ is an even function:

[TABLE]

The substitution of solution (13) in Eq. (14) yields a quadratic equation for $\sqrt{k}$, which determines it as a function of $P$:

[TABLE]

As said above, an important characteristic of the soliton families is the largest value $P_{\max}$ of the total power, up to which the solitons exist. A straightforward analysis of Eq. (15) yields, for $g_{5}>0$ (i.e., the self-focusing quintic term),

[TABLE]

At $P<P_{\max}\left(g_{5}>0;b\right)$, there is a single soliton family with $dk/dP<0$, which is definitely unstable, according to the VK criterion. Note that $P_{\max}\left(g_{5}>0;b\right)$ given by Eq. (16) does not depend on $g_{5}$ (as long as $g_{5}$ is positive). These analytical findings, including the values of $P_{\max}\left(g_{5}>0;b\right)$ given by Eq. (16), are in good agreement with the numerical counterparts displayed in the Fig. 2(b) for $b\geq 3$.

Further, the analysis of Eq. (15) for $g_{5}<0$ demonstrates that it generates two branches of the $k(P)$ dependence. The branch satisfying the VK criterion, $dk/dP>0$, exists in a finite interval of values of the power, under restriction $\left|g_{5}\right|<2$:

[TABLE]

Equation (18) correctly explains an asymptotically linear dependence $P_{\max}(g_{5}=-1;b)$ at large $b$, which is clearly seen in Fig. 4 . In particular, the slope of the asymptotic dependence is exactly predicted by Eq. (18). Subsequent consideration of Eq. (15) readily shows that $P_{\max}(g_{5}<0)=\infty$ at $g_{5}\leq-2$, which actually implies that there remains a single $k(P)$ branch, with $dk/dP<0$, i.e., a completely unstable one.

5 Conclusion

In this work, we have addressed the question of how the recently demonstrated stability area for solitons supported by the NLSE with the CQS (cubic-quintic-septimal) nonlinearity is affected by the spatial confinement of the nonlinearity, which may be realized in optical waveguides based on colloids of nanoparticles. In particular, it was demonstrated that the localization greatly increases the largest power of the stable guided beams, which opens a path to transmit narrow high-power beams by means of the nonlinear confinement. The stability of the spatial solitons in this model is completely determined by the VK (Vakhitov-Kolokolov) criterion. The results were obtained by means of the systematic numerical analysis, and also in the analytical form, valid in the limit of the tight localization of the nonlinearity.

It may be interesting to extend the analysis for a system of two parallel nonlinearly confined waveguides, thus implementing a nonlinear coupler. A challenging possibility is to explore a two-dimensional generalization of the present model, that would imply a nonlinear core embedded in a linear bulk medium (cf. a similar model with the self-focusing cubic-only nonlinearity, considered in Ref. [23]).

J.B.S. thanks Council of Scientific and Industrial Research (CSIR), India, for the financial support. R.R. acknowledges DST (grant No. SR/S2/HEP-26/2012), Council of Scientific and Industrial Research (CSIR), India (grant 03(1323)/14/EMR-II dated 03.11.2014) and Department of Atomic Energy - National Board of Higher Mathematics (DAE-NBHM), India (grant 2/48(21)/2014 /NBHM(R.P.)/R & D II/15451 ) for financial support in the form of Major Research Projects. H.F. and A.G. thank FAPESP and CNPq (Brazil) for financial support. The work of B.A.M. is supported, in part, by the joint program in physics between the National Science Foundation (U.S.) and Binational Science Foundation (U.S.-Israel), through grant No. 2015616.