Stabilization of 1D solitons by fractional derivatives in systems with quintic nonlinearity

We study theoretically the properties of a soliton solution of the fractional Schrödinger equation with quintic nonlinearity. Under “fractional” we understand the Schrödinger equation, where ordinary Laplacian (second spatial derivative in 1D) is substituted by its fractional counterpart with Lévy index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}α. We speculate that the latter substitution corresponds to phenomenological account for disorder in a system. Using analytical (variational and perturbative) and numerical arguments, we have shown that while in the case of Schrödinger equation with the ordinary Laplacian (corresponding to Lévy index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha =2$$\end{document}α=2), the soliton is unstable, even infinitesimal difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}α from 2 immediately stabilizes the soliton texture. Our analytical and numerical investigations of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega (N)$$\end{document}ω(N) dependence (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document}ω is soliton frequency and N its mass) show (within the famous Vakhitov–Kolokolov criterion) the stability of our soliton texture in the fractional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha <2$$\end{document}α<2 case. Direct numerical analysis of the linear stability problem of soliton texture also confirms this point. We show analytically and numerically that fractional Schrödinger equation with quintic nonlinearity admits the existence of (stable) soliton textures at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2/3<\alpha <2$$\end{document}2/3


The model
In the spirit of what was said above, here we consider the substitution of the ordinary Laplacian (second derivative in one spatial dimension) in the quintic NLSE of the form Here, ψ = ψ(x, t) is a wave function and parameter χ distinguishes between repulsive (focusing) χ > 0 and attractive (defocusing) χ < 0 nonlinearities. Our aim is to substitute the ordinary second spatial derivative in (2) by the 1D fractional Laplacian (2) i ∂ψ ∂t = − ∂ 2 ψ ∂x 2 + χ|ψ| 4 ψ. www.nature.com/scientificreports/ where 0 < α < 2 is the Lévy index. The fractional Laplacian (3) is Riesz fractional derivative, see Ref. 29 for details. We note that at α = 2 the fractional 1D Laplacian converts into ordinary second spatial derivative. The details of the corresponding transition can be found, for instance in Ref. 30 . The easiest way to check that is via Fourier transformation, which for second derivative gives −k 2 . This implies that the Fourier image of the fractional Laplacian (3) gives −|k| α or explicitly The expression (5) will be used in subsequent calculations. The equation, having soliton solution for χ < 0 , is obtained from (2), (3) by setting where ω is a soliton frequency. This generates following fractional equation for y (x) where A α is defined by (4). To derive the equation (7), we made substitution p = u − x in the integral (3). At α = 2 the equation (7) transforms into which has the exact soliton solution 9 If we calculate the norm of this soliton solution we see that it is independent of soliton frequency ω so that VK stability criterion (necessary condition for soliton texture stability 15 ) dN/dω < 0 does not hold in this case. Rather, here we have dN/dω = 0 , i.e. the texture is marginally stable. Below we shall show numerically that the VK criterion dN/dω < 0 corresponds to the exponential decay of the linear corrections to the soliton texture for α < 2 . In other words, in our 1D fractional case, the above analytical form of the VK criterion correctly reflects the soliton stability.
Moreover, as it had been shown in Ref. 10 , this marginal stability holds only for a single value of the norm N = N 0 (10). At N > N 0 the soliton (9) collapses, while for N < N 0 it decays into the ground state y = 0 10 . This behavior can be well understood qualitatively by the simple calculation of powers in corresponding energy functional. Namely, the energy functional, corresponding to the Eq. (8), has the form If we consider for a while the solution y 0 (x) as a trial function (with inverse localization radius B being variational parameter and fixed norm N = πA 2 /B ) and substitute it to the functional (11). This yields In the expression (12), first term comes from the "kinetic energy" (dispersion term), i.e. first term in the integrand (11). It is seen that the energy E 0 (12) can be represented in the form E 0 = κB 2 − ωN/2 , where κ = |χ |N 3 12π 2 − N 16 . This immediately shows that the energy E 0 has extremum at B = 0 (corresponding to infinite localization radius, i.e. actually to the already decayed state of the soliton), but maximum or minimum depends on the sign of the coefficient κ . Namely, at κ > 0 we have minimum (at B = 0 ) and at κ < 0 -maximum. The condition κ = 0 occurs just at N = N 0 (10), while κ ≷ 0 corresponds to N ≷ N 0 , giving rise to above behavior, i.e. collapse at N > N 0 and decay at N < N 0 . At the same time, the texture (9) corresponds to the case of κ = 0 , i.e. "no parabola", which means neither maximum nor minimum or aforementioned marginal stability. Note that this reasoning is already very rough as extremum at B = 0 already corresponds to the decayed state with infinite localization radius.
It is instructive to compare the above results with the case of cubic nonlinearity, corresponding to the substitution of the nonlinear term |χ|y 5 0 in the equation (8) by gy 3 0 . The contribution from this term into E 0 is BN 2 /(2π) so that the energy has extremum at finite B extr = 4gN/π , i.e. at finite value of the soliton localization radius, proportional to its norm. This shows that the case of quintic nonlinearity has kind of degeneracy as the contributions (4) A α = Ŵ(1 + α) π sin πα 2 , |p| 1+α dp − ωy + |χ|y 5 = 0, www.nature.com/scientificreports/ coming from dispersion and nonlineariry have the same power. At the same time, for cubic nonlinearity this degeneracy is absent. On the other hand, the soliton in this case is also unstable as E 0 has a maximum. This means, that it is harder to stabilize the soliton in the system with pure quintic nonlinearity as here we should first "gain" the finite soliton radius, where its energy has an extremum. Below we show that the fractional derivatives can well cope with this task.

Perturbation theory
We begin with perturbational treatment of possible soliton solution of the fractional equation (7) as it can be considered to be "exact" in the range of validity of corresponding perturbation theory. In other words, perturbational treatment does not imply any trial function, but has the exact solution (8) as its zeroth approach.
The latter fact implies that the small parameter ε of our perturbation treatment is the "fractionality", i.e. deviation of the Lévy index α from 2: To expand the solution of the equation (7) over small parameter ε (13), we represent it in the form where We have further We see that at α → 2 , the fractional derivative generates logarithmic corrections to the initial "dispersion law" ∼ k 2 . Note, that as at k → 0 k 2 ln n k → 0 , the corresponding integrals at k = 0 are all convergent. The Eqs. (14a)-(14c) can now be recast to the form, suitable for perturbative solution. We have We look for the perturbative solution of (15) in the form where y 0 (x) is defined by (8). In this case, the right hand side of Eq. (15) assumes the form It is seen that, as usually in the perturbation theories, the expressions for higher-order corrections become progressively more cumbersome. Now, in zeroth approach we have, as expected, α = 2 ( ε = 0 ) so that y = y 0 is defined by the expression (9). In the first approach y = y 0 + εy 1 we have while the second approach y = y 0 + εy 1 + ε 2 y 2 yields × εk 2 ln |k| − ε 2 2! k 2 ln 2 |k| + · · · e −ikx dk. www.nature.com/scientificreports/ It turns out that the equations for each perturbative correction y n to the solution y 0 (9) can be solved in quadratures. Namely, it can be shown (this is also seen from the equations (19), (20) for two first corrections), that the general structure of the equation for each perturbation correction has the form where we have already substituted y 0 (9). The properties of operator L in (21) are already known, see, e.g. 31 . This means that we can look for solution of the equations for perturbative corrections (21) in the form of expansion over the eigenfunctions of operator L . But as this operator has also continuous spectrum, latter approach becomes much less profitable than simple method of variations of parameters, known from classical theory of linear differential equations, see, e.g. 32 . To realize the latter method, we need first to find two linearly independent solutions of the homogeneous equation L y = 0 . It is well known that one of the fundamental solutions of the homogeneous equation (21) is always the derivative of the unperturbed soliton solution y 0 (9). It has following explicit form As it is also well-known from the theory of linear differential equations, the second linearly independent solution u 2 (x) can be found as It is seen from expressions (22), (23) that u 1 (x) is odd and decaying function, while u 2 (x) is even and growing function. General solution of the inhomogeneous equation (21) assumes the form Our numerical calculations show that the functions C 1n (x)u 1 (x) are odd, while C 2n (x)u 2 (x) are even. As our soliton texture should be even function, we take the even term only, i.e.
The obtained perturbative solution of our problem permits to calculate the soliton characteristics near α = 2 . As soliton norm plays a major role in the VK soliton stability criterion, we calculate it to obtain where N 0 is defined by (10). In the spirit of VK criterion, as N 0 does not depend on ω , the stability will be determined by the signs of dN 1 /dω and higher order corrections. As dN 1 /dω is the largest term (it is proportional to ε , while subsequent terms are proportional to powers of latter small parameter), the stability would primarily be determined by dN 1 /dω sign. Our numerical calculations (based on inhomogeneous solution (26)) in the first order of perturbation theory show that this is indeed the case, i.e. dN 1 /dω < 0 . This situation is portrayed in Fig. 1. It is clearly seen that small "fractionality" (i.e. small deviation of α from 2, corresponding in this specific case to ε = 0.1 ) indeed stabilizes the soliton, making its norm to be ω dependent. The panel (a) of the Figure www.nature.com/scientificreports/ shows that total norm is positive, while first order correction (even at relatively large ε = 0.1 ) is smaller then N 0 . Figure 1 shows a remarkable property of fractional derivatives: even small "fractionality" stabilizes the initially unstable soliton. This point will be further confirmed below by direct numerical simulation of the soliton texture stability.

Variational treatment. Lévy index limits for soliton solution existence
In the case of "ordinary" soliton with α = 2 , the variational functional has the well-known form (11), where the square of 1D gradient (first spatial derivative) being varied, gives the second spatial derivative in the resulting differential equation (8). Our analysis shows that the fractional analog of the first term in the integrand (11) has the form of the square of the fractional gradient The most convenient way to define the fractional gradient is through its Fourier transform, namely With respect to the above, the variational functional for fractional case can be written in the form It turns out that many general properties of soliton texture (like minimal Lévy index α min , at which the soliton solution of the Eq. (7) seaces to exist) can be obtained without knowledge of explicit functional form of the trial function (although this will be done subsequently for the sake of quantitaive comparison with perturbative and numerical results), but rather of more general "scaling" form of it. Namely, we consider trial function of the form where A and B are variational parameters. The contribution of the first term in (31) can be calculated most conveniently via the Fourier image y(k) of the function y(x). We have The "fractional soliton" norm, calculated in the first order of perturbation theory at several nonlinearity parameters χ and α = 1.9 ( ε = 0.1 , legend in panel (a)). As first order corrections (panel (b)) are negative at ω > 1 , panel (a) shows the overall norm, which turns out to be positive as it should be. Panel (c) shows the derivatives dN/dω ≡ dN 1 /dω , which are negative. At large ω the corrections at any χ sease to depend on ω so that dN 1 /dω → 0 at ω → ∞ . The larger |χ| , the more pronounced is latter behavior.  (33) into the integrand of the first term in (31) generates factor e −ix(k+k ′ ) in it. After integration over x this yields 2πδ(k + k ′ ) , where δ(z) is Dirac δ function. After integration over k ′ we than arrive at the following expression for the first term Here we used the fact that Fourier image (33) is the even function of its argument. The final expression for variational energy then reads where �f n � = ∞ −∞ f n (x)dx . In this case soliton norm As the functions f 2 and f 6 are all positive and independent of A and B, we can find the extremum of E α (35) with respect to these parameters. After lengthy calculations we arrive at following result As A 4 0 and B α 0 should be positive, the function E α has extremum at 3α − 2 > 0 or α > 2/3 . Thus we arrive at the following criterion of soliton solution existence This criterion coincides with that obtained from rigorous mathematical treatment of the localized solutions of the fractional differential equations with 2p + 1 degree nonlinearity in d dimensional space 28 . In the paper 28 this criterion relates the degree of nonlinearity p with space dimensionality d and critical (lower) Lévy index α cr . In our notations this criterion reads α cr = pd/(1 + p) , which for p = 2 and d = 1 gives α cr = 2/3 . Below we shall demonstrate the correctness of criterion (38) by direct numerical solution of Eq. (7).
Note that the expression (36) with respect to (37) permits to obtain the explicit expression for N(ω) in fractional case. Substitution of (37) into (36) yields The expression (39) shows immediately, that as soon as α becomes less then 2 (fractional case), N acquires the dependence on ω , thus giving nonzero derivative dN/dω .
This shows that our variational (actually scaling as no explicit form of trial function has been used) approach predicts the fractional soliton stability according to VK criterion. This is one more chief result of the present consideration, which will be confirmed below by direct numerical simulations.
The above shows that simple scaling arguments, put in the context of a variational method, permit to capture correctly the soliton "phase diagram", i.e. the range of Lévy indices, where the localized solution exists. Moreover, our analysis shows that the extremum (37) corresponds to the minimum, which permits to hope that our soliton texture is stable within the range (38) of admissible α values. Our explicit numerical solution of the fractional equation (7) along with simulations of the soliton texture stability will show that this is indeed the case.
To further study the character of extremum in the fractional case α < 2 , we, similar to the case α = 2 investigate the dependence E α (B) at fixed N. For that we eliminate A from the expression (35) with the help of the relation (36) to get The structure of the expression (40) reveals the fact that as the contributions coming from dispersion ( ∼ B α ) and nonlinearity ( ∼ B 2 ) have different exponents, this immediately implies that contrary to the case α = 2 (12), .

Numerical calculations
To check how both perturbational and variational methods work, it is instructive to compare them with the direct numerical solution of the Eq. (7). For that, we choose the expression (9) as a trial function with A and B being variational parameters. This gives �f 2 � = π/2 , �f 6 � = π/4 and α (34) to be calculated numerically for each α . The substitution permits to obtain for z(x) the equation, which does not contain χ . In other words, the equation for z(x) is similar to Eq. (7) for |χ| = 1 . This means that in our numerical calculations, we can well use the latter equation. The solutions for other χ 's can be easily obtained from those for z(x) with the help of the substitution (42).
The solutions of the equation (7) for |χ| = 1 and ω = 1 are portrayed in Fig. 2a for different α . The shape of solutions at other ω is qualitatively similar to those at ω = 1 . It is seen that as α approaches the critical value α cr = 2/3 ≈ 0.667 , the solution becomes progressively more peaked and at α → α cr , the peak height goes to infinity. In other words, at α = α cr = 2/3 we have "fractionality collapse" of the soliton. We may speculate here, that as Lévy index α can be viewed as a measure of disorder in a system (see 3,23-25 and references therein; the more α deviates from 2, the stronger is the disorder), the value α cr = 2/3 may be regarded as some threshold disorder strength, at which the solitons cease to exist in a system. Really, at α < α cr , our iterative process of the numerical solution diverges so that no localized solutions of the equation (7) could be found.
Latter outcome is in accord with the results of variational treatment (38). This shows that our variational treatment captures well the "phase diagram" (range of admissible Lévy index values) of the fractional soliton. Below we will see that it gives also very good (indistinguishable in the scale of plot) approximation to the numerical solution of the fractional equation (7). This shows that this method could be considered as an efficient analytical tool to obtain and study the solution of the soliton equations with fractional derivatives. Note, that as no explicit expression for trial function has been introduced (see (32)), the above variational method is not sensitive to the trial function details and thus could be of general nature. This means that it can be well applied for more complex models of "fractional" solitons, like those with mixed or saturated nonlinearities.
We now compare the variational, perturbational (in the range of validity, i.e. at α close to 2) and numerical results for both soliton structure and its norm. Such comparison is reported in Fig. 3. The main panel shows the comparison of variational and numerical soliton structures y(x) at some selected values of α . It is seen that the closer is α to 2, the better is a coincidence. On the other hand, for α close to its critical value 2/3, the overall coincidence worsens (although it remains within limits of 5%) especially at the "tails" of the y(x) function. This is because our simple trial function (7) has exponential asymptotics at infinities, while it is well-known (see, e.g., [2][3][4]28,33 ) that the asymptotics of the fractional differential equations solutions is usually power-law (1). Latter power-law asymptotics starts to manifest itself for α 's somewhere close to 2/3. This means that a more sophisticated, multi-parametric, trial function would improve coincidence with numerics near threshold α value.
Right inset in Fig. 3 reports a comparison of numerical and variational dependences N(ω) , obtained from the corresponding solutions at different α . This shows that the variational solution gives a pretty good approximation to the numerical one also at a wide range of the frequencies 0 < ω < 10 . The coincidence of N(ω) curves even at α close to 2/3 (e.g. for α = 0.8 ) shows that for integral characteristics, like N(ω) , the details of solution behaviour is not that important. Thus, for soliton characteristics, involving integrals of y(x), even the simplest possible trial function (9) gives very good coincidence with the numerical curve, which requires a substantial amount of time and computer resources (to calculate N(ω) , we should first solve (7) numerically for each ω ) to be obtained.
Left inset in Fig. 3 compares numerical, variational and perturbational (up to the first order in small parameter ε (13)) solutions y(x). The solutions, obtained at α = 0.9 , are coded by colors and it is seen that they are indistinguishable in the scale of the plot. Our analysis shows that maximal error here is around 0.1%. We choose α = 0.9 ( ε = 0.1 ) for the first order of perturbation theory to be valid. It is obvious, that if we consider also the higher orders in the small parameter ε , we can descent further in α . As consideration of the higher orders of perturbation theory is usually associated with the growth of complexity of calculations, the variational solution in this sense is much more profitable. This is because it gives the approximate (as we saw above, the accuracy of the approximation is generally very good) analytical solution of the problem for all admissible values of Lévy index α . On the other hand, the perturbational solution can be regarded as "exact" (it does not use any variational ansätse) in the range of its validity. This means that perturbational soliton may be considered as a complementary to variational one near α = 2 . This is especially true in view that in the present problem the solution for α = 2 (zeroth approximation in perturbation theory) (9) is exact.
We are now in a position to study the soliton stability numerically. This is because the direct use of VK result dN/dω < 0 can sometimes adulterate the real (i.e. numerical) stability criterion. This is especially true for the fractional case as the theorem of Vakhitov and Kolokolov 15 had been stated and proved for the solitons having ordinary second spatial derivative as a dispersion term. To accomplish the above task, we consider the linear stability problem for the soliton solution of the initial fractional NLSE (2) with fractional Laplacian (3). Namely, we consider following substitution Here y(x) is undisturbed (by the small perturbations p and q ) soliton texture (i.e. the above numerical solution of the Eq. (7)) and asterisk means complex conjugation. Substitution of (43) into fractional NLSE (2), (3) with its subsequent linearization over p and q generates following eigenvalue problem (43) ψ(x, t) = y(x) + p(x)e t + q * (x)e * t e iωt , p, q << y. www.nature.com/scientificreports/ where −|�| α/2 is defined by (3) and function z(x) by (42). By solving numerically the spectral problem (44), it can be shown that our texture y(x) is stable if all real parts r of the eigenvalues are zero. It turns out that this is the case for the entire domain of α (38) and ω values. Our analysis shows that the spectrum of imaginary part, i , consists of negative values only, i.e. i i = −p n , where n enumerates the eigenvalues. This means that the corrections ∼ e t = e −p n t (43) to the soliton solution y(x) decay exponentially so that our soliton texture is stable. At the same time for α = 2 , the imaginary part i equals zero so that soliton is marginally stable. Similar to our preceding results, the soliton gains stability at α < 2 only. Also, the modulus of the smallest (i.e. "most negative") eigenvalue grows as α → 2/3 . For instance, for ω = 1 our numerical calculations show that p min (α = 1.9) = −0.618 ( p min is the smallest eigenvalue so that the correction decays primarily like e −0.618t ), p min (α = 1.7) = −1.294 , p min (α = 1.3) = −3.435 , p min (α = 1.0) = −9.162 and p min (α = 0.74) = −210.1 the lowest α , accessible in our numerical calculations. Most probably, at α = 2/3 p min → −∞ . This means that in fractional case our soliton texture may be prone to collapse at α = 2/3 since the collapsed configuration can be considered as "absolutely stable". However, as in real physical setups one cannot change α (say, it is related to the degree of disorder in a sample), for each specific Lévy index value, the soliton texture is stable.
To demonstrate the considered soliton stability, it is instructive to check it by the direct numerical simulations of the soliton texture time evolution. In the spirit of linear stability ansats (43), for our numerical simulations, we consider the following substitution where y(x) is a static texture (real function), obtained by the numerical solution (42) of the equation (7) for each α . In this case, the soliton stability would mean the decay of δψ(x, t) as time goes to infinity. In other words, as time passes, for a stable case, the initial texture tends to y(x). As function δψ(x, t) is complex, for our pictorial demonstration of the soliton stability, we simulate the dynamics of the modulus |ψ(x, t)| = √ ψ(x, t)ψ * (x, t) . Note, that if we choose too large amplitudes of the real and imaginary parts of the initial perturbation δψ(x, 0) , the soliton will lose its stability. Namely, it will either collapse or "smear down" into the ground state ψ = 0 . Also, in this case, the numerical method instability may come into play. The criterion of our choice of perturbation δψ(x, 0) was that for the "regular" case α = 2 the texture stays in its initial state for a sufficiently long time t > 500 . Also, our observations show, that the choice of δψ(x, 0) in the form, which yields the negative contribution into ψ(x, t) , may "stabilize" the soliton so that slightly larger amplitudes of δψ(x, 0) are admissible. The results are reported in Fig. 4 for three representative values of Lévy index α . In Fig. 4, three upper panels correspond to positive initial perturbations and three lower -to the negative one. To demonstrate, how relaxation to the static case "works", we intentionally take α = 1.999 rather than α = 2 . It is seen, that for positive perturbation, the relaxation to the static texture occurs from above, while for negative one from below. The main observation, however, is that the relaxation to the static texture y(x) occurs quicker at smaller α in accordance with the above stability analysis. Namely, while at α = 1.999 , the relaxation to (both for positive and negative perturbations) y(x) "finishes" at t = 100 , for α = 1.3 this time is 100 times smaller, i.e. for t > 1 , we already have the static texture. As we approach the critical value α = 2/3 , this time tends to zero.

Outlook
The chief message of the present consideration is that the "fractionality" (deviation of Lévy index α from 2), stabilizes the soliton texture. The physical origins of the "fractionality" may be different. Although the explicit discussions of the underlying physical mechanisms are scarce, we can distinguish two main implementations. One is possible solid-state realization, called Lévy crystal 34 . This model deals essentially with the discretization of the fractional Schrödinger equation to obtain the effective tight-binding lattice model, where the hoppings between lattice cites obey Lévy statistics. To generate latter statistics experimentally, hoppings in the lattice must be carefully adjusted, which may be a prohibitively complex task. Another implementation is related to the light dynamics in an optical cavity 35 . Although this is readily realizable physically, the reason for the usage of fractional derivatives in corresponding linear Schrödinger equation was phenomenological, i.e. the specific dispersion law, generated by the fractional derivative, was needed to realize the physical setup. We mention also one more solid-state model, dealing with polariton condensates 36 . In this study, however, the question of the physical origin of "fractionality" had not been considered at all. On the other hand, recently we proposed [23][24][25] the possible physical origin (to some extent justification of usage) of "fractionality" coming from the strong (i.e. non-Gaussian) disorder in a system. The main point here 25 is that in a disordered system we usually have a broad, non-Gaussian distribution of its physical quantities. In quantum mechanics, such strong disorder usually leads to localization of initially (before the introduction of a disorder by, say, doping by some impurities) itinerant states. This is the essence of the famous Anderson localization phenomenon 37 . In the context of our solitonic problem, this means "more localization" of the soliton (for instance its collapse at α = 2/3 ) so that it becomes "more stable" than in the ordinary non-fractional case, corresponding to α = 2 . We may call our effect the soliton stabilization by disorder.
Although here we consider the simple model with quintic nonlinearity in one spatial dimension, our preliminary consideration shows that this is the case for more complicated models, containing more intricate nonlinearities (like competing cubic-quintic, higher order like septic and/or saturable nonlinearities) as well as external space and time dependent potentials. A noteworthy result is the variational analysis of the problem, based on the scaling argument only. In other words, we take the trial function in the scaling form (32) without www.nature.com/scientificreports/ any suppositions on its exact shape. The only supposition was that the integrals f 2 and f 6 should be convergent. Such scaling analysis permits us not only to establish (known from the mathematical literature, see 28 and references therein) the range of α's, where the soliton solution exist but to show that already at an infinitesimal deviation of α from 2, the soliton acquires stability according to VK criterion, see Eq. (39). As the conclusion of soliton stability has been confirmed by the numerical analysis, we can assert, that our simple "scaling-variational" procedure is suitable to study also the properties of more complicated soliton textures. One more important result is the construction of perturbation theory near the ordinary case α = 2 . The point is that such perturbational expansion is the only analytical tool to investigate the soliton textures in models of any complexity when variational ansats is hard to conceive for some reasons. Although perturbational calculations (especially in higher orders) may be very cumbersome, they may serve as guidance for numerical simulations as the nonlocal character of fractional derivatives makes even the numerical calculations tricky. Also, the stability effects (for instance using VK criterion) can be investigated analytically within the latter approach.
An interesting generalization of our fractional soliton problem is a substitution of time derivative in Eq. (2) by its fractional counterpart. Our preliminary analysis shows that such substitution has implications on soliton dynamics and stability. We postpone the discussion of the latter interesting dynamical effects for future publications. The results reported here are relevant for disordered multicomponent and nonlinear optical systems both in 1D and higher dimensions. They also may be used to control (e.g. by varying Lévy index α ) the properties of Bose-Einstein condensates.

Methods
The details of our theoretical methodology and those of working with fractional derivatives and fractional Laplacians, in particular, have been described in the sections "The model", as well as in the two following sections "Perturbation theory" and "Variational treatment. Lévy index limits for soliton solution existence". The numerical simulations of soliton structure and stability have been conducted using the commercial Mathematica software package as well as C++ routines, partially written ad hoc and partially taken from the standard libraries like LAPACK.  www.nature.com/scientificreports/