Preventing critical collapse of higher-order solitons by tailoring unconventional optical diffraction and nonlinearities

Self-trapped modes suffer critical collapse in two-dimensional cubic systems. To overcome such a collapse, linear periodic potentials or competing nonlinearities between self-focusing cubic and self-defocusing quintic nonlinear terms are often introduced. Here, we combine both schemes in the context of an unconventional and nonlinear fractional Schrödinger equation with attractive-repulsive cubic–quintic nonlinearity and an optical lattice. We report theoretical results for various two-dimensional trapped solitons, including fundamental gap and vortical solitons as well as the gap-type soliton clusters. The latter soliton family resembles the recently-found gap waves. We uncover that, unlike the conventional case, the fractional model exhibiting fractional diffraction order strongly influences the formation of higher band gaps. Hence, a new route for the study of self-trapped modes in these newly emergent higher band gaps is suggested. Regimes of stability and instability of all the soliton families are obtained with the help of linear-stability analysis and direct simulations. Stabilising localised solitons in higher dimensions is more challenging than their one dimensional counterpart due to the onset of critical collapse. Here, competing nonlinear terms are added to a linearly-modulated optical lattice to predict regimes of stability for different soliton families.

S olitons, alias localized waves, are stable objects balanced by diffraction/dispersion and nonlinearity [1][2][3] . The stable propagation and self-trapping of two-and three-dimensional (2D and 3D) solitons have recently received reinvigorating scientific interests in physics and beyond 4 . A fundamental challenge in this field of research is the stabilization of solitons in multidimensional coordinates, since the 2D and 3D solitons in free space are always unstable and undergo respectively critical and supercritical collapses arising from catastrophic self-focusing nonlinearity 5,6 . It is therefore intriguing to suppress these collapses in the recent soliton research 7,8 .
The stabilization of multidimensional localized states usually relies on nonuniform media. A common way is to introduce linear periodic potentials (alias linear lattices) 7-10 to a material with uniform nonlinearity, e.g., a periodic transverse modulation of the refractive index (photonic crystal and lattice) in optics 11,12 or an optical lattice generated by counter-propagating laser beams in Bose-Einstein condensates [13][14][15][16][17] . The periodic potentials play a paramount stabilizing part for D-dimensional localized modes including fundamental solitons and more complicated states like gap solitons, multihump states, and vortex solitons 4,7,8 . Another potential scheme is assisted by nonlinear pseudopotentials-the well-known nonlinear lattices 7,[18][19][20] . Although nonlinear lattices with smooth variation of nonlinearity may support various species of solitons, to stabilize multidimensional solitons using sinusoidal nonlinear lattices remains an open area of research 7 . Further, the combined linear-nonlinear lattices yield stabilization of various soliton families [21][22][23] .
Conversely, materials which exhibit different orders of nonlinearity, including the well-known cubic-quintic model with competing nonlinearities (self-focusing cubic and self-defocusing quintic nonlinear terms), can create stable multidimensional solitons and suppress the above-mentioned critical-and supercritical collapses 4,7,8,21,22 . Experimentally, fundamental solitons and vortex solitons in 2D and 3D geometries have been observed in CS 2 24,25 and in suspensions of metallic nanoparticles 26 . Theoretical predictions 27,28 and experimental observations 29-32 very recently confirmed the existence of stable multidimensional states of ultracold bosonic gases 8 , appearing as the so-called quantum droplets 27,28 , where the Bose-Bose mixture collapses are suppressed by quantum fluctuations 33 .
Beyond strategies affording the stabilization of multidimensional solitons in standard quantum mechanics, two extensions of which have been elaborated in past decades. One is the non-Hermitian system whose Hamiltonian is parity-time (PT ) symmetric and, due to this unique property, such systemcounterintuitively-exhibits entirely real eigenvalue spectra [34][35][36] . In recent years, the studies of various types of solitons 37 and their propagation dynamics governed by different PT -symmetric systems integrated with PT lattices (periodic potentials with alternating gain and loss regions) are fruitful [38][39][40][41][42] . Another extension, made by Nick Laskin, is space-fractional quantum mechanics (SFQM) 43,44 . The SFQM is achieved if the Lévy flights replaced the Brownian trajectories in Feynman path integrals, leading to a fractional diffraction order characterized by Lévy index α (1 < α ≤ 2) rather than the standard one with α ¼ 2. As the Schrödinger equation roots in the path integral over Brownian paths, the path integral over Lévy trajectories results into fractional Schrödinger equation (FSE) [43][44][45] .
FSE is currently attracting a great deal of attention in different branches of physics, including condensed-matter physics [46][47][48] and quantum physics [49][50][51] . Importantly, a new hallmark of FSE development published in 2015 by Longhi 52 found that a fractional quantum harmonic oscillator can be realized based on transverse laser beam dynamics in aspherical optical cavities, pushing the studies of fractional physical system into optics. Such work triggers a growing prosperity of fractional models for digging deep into the propagation properties of the possible beam solutions in both the linear and nonlinear systems, striking examples include: Gaussian beams propagation 53,54 , conical diffraction of light (PT symmetry) in FSE with a periodic PT -symmetric potential 55 , accessible solitons in FSE with a harmonic potential 56 , solitons 57,58 , and modulational instability 59 in purely nonlinear fractional Schrödinger equation (NLFSE), and diverse soliton families in the NLFSE setting by adding to its linear [60][61][62][63] or nonlinear 20 part with a periodic potential (the aforementioned optical and nonlinear lattices, respectively). Although the NLFSE provides a fertile ground for exploring various types of solitons, their existence and stability property supported by the cubic-quintic nonlinearities and 2D linear periodic potential (optical lattice) are yet to be revealed.
In this work, we propose and demonstrate, theoretically and numerically, a framework of 2D NLFSE-describing light propagation in a nonlinear periodic system with an optical lattice and attractive (self-focusing)-repulsive (self-defocusing) cubic-quintic nonlinearities-which can suppress the critical collapse mentioned before. We reveal that the model produces a variety of stable soliton families, including 2D fundamental gap and vortical solitons as well as gap soliton clusters (solitons are always unstable in the quintic-only model). Noticeably, the Lévy index α and propagation constant b impact deeply upon the profiles of soliton families, the side peaks of which with higher modulations when decreasing α or increasing b. By performing the linearstability analysis and direct numerical simulations, the stability and instability diagrams of all the soliton families are acquired. A detailed insight into the dynamic properties of solitons further shows that the solitons are robustly stable in the middle of the band gaps of the underlying linear Bloch spectrum, while unstable near the edges of the band gaps; and the stability of the solitons is moderately influenced by nonlinear strength. Our work proposes a feasible scheme to stabilize 2D localized modes against critical collapse by considering the fractional diffraction order to light propagation in periodic physical systems with competing selffocusing and self-defocusing nonlinearities in cubic-quintic nonlinear terms, thus offering a new avenue to investigate the existence and dynamic properties of 2D localized modes by managing the diffraction order and tunable band gaps of the systems.

Results
Theoretical model and its linear spectrum. We start from the cubic-quintic NLFSE (with competing focusing-defocusing nonlinearities), which describes a laser beam propagation in a nonlinear medium and can be expressed in scaled form: where E and z represent the field amplitude and propagation distance, respectively, γ > 0 being the coefficient of self-focusing cubic nonlinearity. Note that fractional Laplacian ∇ α ¼ ðÀ∂ 2 x Þ α=2 þ ðÀ∂ 2 y Þ α=2 acts on two transverse coordinates x and y, with α (1 < α ≤ 2) corresponds to Lévy index. In particular, Eq. (1) restores to the conventional nonlinear Schrödinger equation at α ¼ 2. The remaining parameter V is the linear (square) potential trap, which yields where V 0 > 0 stands for the amplitude (modulation depth) of optical lattice.
To explore the existence of 2D nonlinear localized gap modes, we should first give out the band gap structure of guided Bloch modes (linear Bloch waves) stemming from the linearized model of Eq. (1) (see "Methods"). For the given expression of the sinusoidal square optical lattice (2), whose contour plot is portrayed in Fig. 1a, the corresponding band gap structure of Bloch-wave spectrum, as mentioned above, can be described and found by adopting the analogies with the solid energy band theory of crystalline materials and which, in the way as usual, should be based on the lowest (first) Brillouin zone of the lattice potentials. Such reduced Brillouin zone of the square lattice under considered is displayed in Fig. 1b. In the view of these foundations, we have plotted the dispersion diagram for guided Bloch waves with the variation of Lévy index α in Fig. 1c, and shown such dispersion spectrum for a particular value of α (α ¼ 1:3) in Fig. 1d. It is observed from the former panel that both the first and second finite band gaps widen with an increasing α, conforming to the previous reported results in 1D situations 63 . As for the latter panel, choosing α ¼ 1:3 is qualified to guarantee the moderate wideness of the first two band gaps, within which various localized gap modes may be stabilized by the competing cubic-quintic nonlinear terms.
Formation of two-dimensional nonlinear localized gap modes. Having obtained the finite band gaps of the underlying linear spectrum to the model (Eq. (1)), in the next part of this research work we are going to explore the formation of 2D nonlinear localized gap modes (with their corresponding propagation constants lying in such band gaps) supported by the interplay of the given linear lattice and the competing cubic-quintic nonlinearities. Several classes of soliton families are identified, such as fundamental gap solitons and vortical ones as well as the gap-type soliton clusters. The physical mechanism and theoretical analysis are also given out to all the localized modes thus found.
To understand how the Lévy index (α) in the fractional Laplacian (diffraction) term affects the shape of fundamental gap solitons, we plotted the spatial structure of such fundamental localized modes at two different values of α by means of 3D plots, top views and contour profiles in Fig. 2a-f. Comparing the solitons' shape in these panels, we can see that the fundamental gap soliton is in a wavy shape at smaller α (Fig. 2a-c), while such wavy modulation of the gap soliton disappears at relatively large α (Fig. 2d-f); this characteristic resembles those found in analog soliton-based systems 60,63 .
Our systematic simulations by exploiting the undermentioned linear stability and direct perturbed propagation methods (see "Methods" for detailed descriptions) have demonstrated the presence of such fundamental gap solitons dwelling in the first and second band gaps, and identified their stability and instability properties, which are accumulated in the associated P $ b curve in Fig. 3a. From such curve one can clearly observe that the stability region is much broader when preparing the gap solitons in the first band gap than that in the second gap; and intuitively, the unstable gap solitons emerge once the propagation constant is near the edge of the finite band gap, as within where the effective mass of the gap solitons may change the sign 7 . Note that such curves (P $ b) for multipole gap solitons (or gap soliton clusters) Fig. 1 Shape, Brillouin zone, and linear spectrum of the 2D linear (index) potential. a Contour plot of a sinusoidal square lattice in Cartesian space, blue shading corresponds to potential minima. b The first reduced Brillouin zone of such 2D square lattice in reciprocal lattice space; tagged are the high-symmetry points. c Band gap structure for guided Bloch waves with different α. d The dispersion spectrum at Lévy index α ¼ 1:3, with the colored curves denoting Bloch bands. V 0 ¼ 6 and γ ¼ 0:5 are used throughout this paper. with numbers of peaks n ¼ 4 and n ¼ 16 are showed respectively in Fig. 3b, c. Three representative scenarios of such fundamental gap solitons, residing on the first two band gaps, are severally depicted in Fig. 4(a, d), (b, e) and (c, f) (which correspond to the marked circles (A1, A2, A3) in Fig. 3a). When placing in much higher band gap, e.g., the second band gap, the gap soliton becomes more localized than that in the first band gap, as comparing the Fig. 4a, d with 4b, e. It is observed from the Fig. 4c, f that the gap soliton becomes highly confined mode occupying several lattice sites when the propagation constant b is close to the edge of the second band gap, with an emergence of a cusplike behavior, such feature is typical for the gap modes in other physical settings 7,16 . The eigenvalues, calculated from the linear stability equations (see the Eq. (4) in "Methods"), of these three gap modes are displayed in Fig. 4g-i, showing that the former two modes are linearly stable while unstable for the latter mode which lies inside the upper edge of the second band gap.
On the other hand, our research model (Eq. (1)) upholds families of higher-order solitons, which we call gap soliton clusters. Depicted in Fig. 5 are the quadruple-mode solitons of such soliton clusters, with a distance (4s) between isolated gap solitons doubling to the period of the lattice. Because of the existence of repelling interaction between the each adjacent solitons of the soliton clusters, the stability diagram of the quadruple-mode solitons, portrayed as the dependency PðbÞ in Fig. 3b, is much narrower than that for the fundamental gap solitons. The stability and instability properties of such multiple solitons with the number of isolated gap solitons n ¼ 4 are further validated by the corresponding eigenvalues analysis, which are showed in the bottom of Fig. 5.
The successful stabilization of quadruple-mode solitons enlightens us to generate multiple solitons with more number of gap solitons, which can be naturally achieved by the model under studied. For instance, the gap soliton clusters with the  neighboring distance 4s ¼ 2π and n ¼ 16 solitons arranging as 4 4 configuration, are displayed as stable modes localized in the first and second band gaps in Fig. 6a, d, b, e, respectively, and as unstable object when its propagation constant approaches the edge of the band gap in Fig. 6c, f. We stress that the multiple gap solitons or gap soliton clusters reported here are equivalent to the recently found truncated nonlinear Bloch waves or gap waves in nonlinear optical media and ultracold atoms 17,23,[64][65][66][67][68] .
Stable propagations of the localized gap modes thus found (fundamental gap solitons and soliton clusters), located inside the first and second band gaps, are depicted respectively in Fig. 7 which, in practice, demonstrates the stability properties (typical for solitons) of the localized gap modes in a direct and vivid way. Typical propagation properties of unstable gap modes are depicted in Fig. 8, indicating that the unstable modes exhibit a weakly instability (whereas the main peaks can sustain, their side peaks grow quickly during the propagation), which is in an excellent agreement with the relevant linear eigenvalues results as shown in the bottom lines of the Figs. 4-6.
An interesting result is the possibility to find stable localized gap vortices, the gap solitons carrying with topological charge m. A typical gap vortex is usually constructed as quadruple-mode bound states of identical solitons with modulated phase structure. Characteristic examples of such gap vortices with m = 1 and spacing 4s ¼ 2π are shown in Fig. 9. It is shown that the complicate phase structures of these gap vortices are modulated spatially owning to the influence of the lattice potential. Eigenvalues of these vortical solutions, depicted in the bottom of Fig. 9, indicate that there is weak linear instability for the unstable gap vortices, such effect equates to that for the localized gap modes without vorticity (m = 0) (see the panels in the bottom right corner of Figs. 4-6) and is verified as well in their direct perturbed propagation in Fig. 10. For defined values of propagation b and cubic nonlinear coefficient γ, we find that the stability and instability regions of the vortex gap solitons alter with the variation of Lévy index α , as attested by the dependence P(α) in Fig. 10a.

Discussion
As previously described, the 2D solitons created in pure selffocusing settings (uniform Kerr media) undergo critical collapse, and their vortex counterparts suffer a greater azimuthal instability. To suppress such adverse factors, scientists usually resort to two alternative methods-the cubic-quintic model with competing self-focusing and self-defocusing nonlinearities, and the introduction of external linear potentials such as optical lattices. By combing such two methods, we here numerically explored and with an analysis of the generation of localized modes in the 2D cubic-quintic nonlinear FSE with an optical lattice. We discovered that such model gives rise to a class of 2D stable soliton families, which include 2D fundamental gap solitons and their clusters as well as vortex solitons. Like previous studies in similar but different models, we demonstrated that all the gap soliton families are always unstable physical objects in the quintic-only model (the cubic-quintic model discarding the cubic term). We emphasize that, the two physical parameters, the Lévy index α and propagation constant b, impose a nontrivial effect upon the profiles of soliton families of all kinds, the solitons accompanied by a higher modulated side peaks if we simply   minish α or enlarge b. We identify the stability and instability diagrams of all the soliton families by right of the linear-stability analysis by calculating the eigenvalues of the underlying soliton solutions against linear infinitesimal background perturbations, and the direct numerical simulations. We conclude that the solitons are stable within the band gaps (of the underlying linear system) and unstable close to the band edges.
An obvious conceptional extension of our study is to reveal the underlying physical mechanism for critical collapse related to physical parameters including Lévy index α (diffraction order) and order of a single nonlinearity (e.g., cubic only) as well as the strength of optical lattices. Rich localized modes and nonlinear dynamics are yet to be explored in an incommensurate structure -combined linear and nonlinear lattices [21][22][23] . Another interesting issue is to study the localization and delocalization properties of waves in other periodic lattices 69 with regolabile (fractional) diffraction.

Methods
Additional details on the numerical procedure. To obtain the stationary solution U, we first rewrite the field amplitudes as E = U exp(ibz) with the propagation constant b and plug them into the Eqs. (1), then we get the following stationary equation: An essential physical quantity is the soliton power P, which is defined as P ¼ R R Uðx; yÞ j j 2 dxdy. For another, the linear stability analysis of the thus-found gap solitons is a critical issue. To this, we express the perturbed field amplitude as E ¼ ½Uðx; yÞ þ pðx; yÞexpðλzÞ þ q Ã ðx; yÞexpðλ Ã zÞexpðibzÞ, here Uðx; yÞ represents undisturbed field amplitude, p(x, y) and q Ã ðx; yÞ being small perturbations at eigenvalue λ. Substituting such perturbed solution into the Eq. (1), we readily  obtain the following eigenvalue problem: iλp ¼ þ 1 2 ∇ α p þ ½b þ Vðx; yÞp À γUð2U Ã p þ UqÞ þ U 2 U Ã ð3U Ã p þ 2UqÞ; iλq ¼ À 1 2 ∇ α q À ½b þ Vðx; yÞq þ γU Ã ð2Uq þ U Ã pÞ À ðU Ã Þ 2 Uð3Uq þ 2U Ã pÞ: We stress that from the above eigenvalue Eq. (4), the perturbed solutions are stable as long as all the real parts ðλ R Þ are zero (λ R ¼ 0). The numerical recipes we used follow the ways: the calculations of stationary solutions of Eq. (3) and numerical propagation simulations of Eq. (1) under small initial perturbations were, respectively, performed based on the modified squared-operator method and splitstep Fourier method 70 .

Data availability
The data that supports the results of this work are available from the corresponding author upon reasonable request.