Metastable soliton necklaces supported by fractional diffraction and competing nonlinearities

We demonstrate that fractional cubic-quintic nonlinear Schr\"odinger equation,characterized by its L\'evy index, maintains ring-shaped soliton clusters ("necklaces") carrying orbital angular momentum. They can be built, in the respective optical setting, as circular chains of fundamental solitons linked by a vortical phase field. We predict semi-analytically that the metastable necklace-shaped clusters persist, corresponding to a local minimum of an effective potential of interaction between adjacent solitons in the cluster. Systematic simulations corroborate that the clusters stay robust over extremely large propagation distances, even in the presence of strong random perturbations.


I. INTRODUCTION
Ring-shaped soliton clusters, constructed as coherent chains of fundamental solitons, are quasistationary multisoliton bound states. These complexes generalize the concept of necklace-ring beams, which were first studied in self-focusing Kerr media [1][2][3], and later extended to the nonlinear Schrödinger equation (NLSE) with saturable nonlinearity [4].
Usually, many-soliton circular structures tend to destroy themselves through expansion or collapse, an experimentally relevant issue being how to maintain their quasi-stationary propagation. To balance the destructive trends, a vortical phase profile can be imposed onto the clusters [5]. In this context, various nonlinearities, including saturable [6], competing quadratic-cubic [7], cubic-quintic [8], and nonlocal [9], were employed to enhance the robustness of the circular soliton clusters. It was also demonstrated that two-dimensional (2D) soliton clusters can be supported in active systems, such as pumped nonlinear optical cavities [10,11]. A possibility of the existence of 3D soliton clusters, composed of spatiotemporal optical solitons, has also been explored [12][13][14]. Creation of necklace beams was experimentally demonstrated in photorefractive [15], nonlocal [16], and Kerr [17] media. More recently, assemblage of metastable necklace clusters of 2D "quantum droplets" was predicted [18] in the framework of the Gross-Pitaevskii equations with the cubic nonlinearity modified by the logarithmic * Electronic address: lpf281888@gmail.com factor, which is induced by the Lee-Huang-Yang correction to the mean-field dynamics of Bose-Einstein condensates (BEC) [19]. Further exploration of physically relevant settings which maintain necklace complexes in a quasi-stable form remains a relevant objective for theoretical and experimental work. In particular, the creation of nearly stable necklace patterns with imprinted overall vorticity in a medium with the cubic-quintic (CQ) nonlinearity (liquid carbon disulfide) was very recently reported in Ref. [20].
In this connection, it is relevant to stress that the vortical soliton clusters and usual axisymmetric vortex solitons [21][22][23]76] are different types of beams carrying the orbital angular momentum (OAM). The clusters carry OAM in the explicit form, being rotating azimuthally modulated necklace-like structures, while the vortex solitons carry OAM in the "hidden" form, being axisymmetric structures, while OAM is stored in their phase circulation. Circular soliton clusters, considered in the present work, are constructed as quasi-steady states, realizing local minima of an effective potential. Vortex solitons are fully stationary modes with a topological charge, s, representing the vorticity. Accordingly, stable vortex solitons keep their shape in the course of the propagation, while unstable ones are broken up by azimuthal perturbations in several fragments, whose number is usually (s+1). On the other hand, the circular clusters ("necklaces" ) should be created differently, as they do not emerge as result of the splitting of unstable vortex solitons. In particular, the number of separate elements ("beads" ) in the necklace is independent of s. Unlike the uniform vortex modes with the "hidden" vorticity, for clusters their steadiness is a major issue.
Theoretical studies of optical solitons in NLSE models were lately extended in two noteworthy directions. On the one hand, the theory has been expanded for non-Hermitian optical systems, which are modeled by NLSEs with parity-time-symmetric complex potentials, see reviews. [24][25][26][27] and an experimental report [28]. On the other hand, much interest has been drawn to the fractional Schrödinger equations (FSEs), which were originally proposed by Laskin [29][30][31] and, in a rigorous mathematical form, by Hu and Kallianpur [32]. In the quantum theory, FSE was derived [29] as a model in which the consideration of Feynman path integrals over Brownian trajectories leads to the standard Schrödinger equation, while path integrals over " skipping" Lévy trajectories lead to the FSE [30]. Experimental schemes have been proposed for realization of FSE in condensed-matter settings [33,34] and in optical media [35], where the effective fractional diffraction may be realized in cavities including optical filters and phase masks [36,37]. Actually, optical cavities offer a testbed to explore the propagation dynamics and eigenmodes of optical beams governed by linear and nonlinear FSEs. Following this approach, many noteworthy properties have been revealed in the framework of FSE [38][39][40][41][42][43][44][45][46][47][48][49][50]. Adding nonlinear terms to FSE [51,52], recent publications have predicted a variety of fractional solitons .
In this work, we aim to construct circular soliton clusters in the model combining fractional diffraction and competing CQ nonlinear terms. The paper is organized as follows. The model, including a discussion of its implementation in optics, and numerically found fundamental solitons, is introduced in Sec. II, which is followed by identification of an effective potential energy of the interaction between adjacent solitons, and construction of necklace-shaped soliton clusters, in a semi-analytical form, in Sec. III. Robust propagation of metastable soliton clusters under the action of perturbations is additionally reported, by means of systematic numerical simulations, in Sec. IV. The paper is concluded by Sec. V.

II. THE MODEL AND FUNDAMENTAL SOLITONS
We consider beam propagation along the z-axis in a nonlinear isotropic medium with the CQ nonlinear correction to the refractive index, n nonlin (I) = n 2 I − n 4 I 2 , where I is the light intensity, n 2 and n 4 being coefficients accounting, respectively, for the cubic self-focusing and quintic defocusing. The respective fractional NLSE, with propagation distance z and transverse coordinates (x, y), is where A(z, x, y) is the local amplitude of the optical field, the intensity being I ≡ |A| 2 , k 0 = 2πn 0 /λ is the wavenumber, which is determined by carrier wavelength λ, and n 0 is the background refractive index. The fractional-diffraction operator in Eq. (1), with the corresponding Lévy index α, belonging to interval 1 ≤ α ≤ 2, is defined as [77][78][79], where F and F −1 are operators of the direct and inverse Fourier transform, k x,y being wavenumbers conjugate to transverse coordinates (x, y). By means of rescaling, Ψ(ζ, ξ, η) = n 4 /n 2 A(z, x, y), ζ = k 0 n 2 2 /n 0 n 4 z, and (ξ, η) = 2k 2 0 n 2 2 /n 0 n 4 1/α (x, y), Eq. (1) is cast in the normalized form, Equation (1) with α ≤ 1 and cubic-only self-focusing gives rise to the wave collapse [71], which destabilizes solitons; however, the quintic self-defocusing term suppresses the instability and makes it possible to construct stable soliton states for α ≤ 1 as well [72].
In the linear limit, Eq. (3) models the free transmission of optical beams under the action of fractional diffraction. Experimental setups which realize this propagation regime were proposed in Refs. [35] and [39]. In particular, the former setting is designed as a Fabry-Perot resonator, with two convex lenses and two phase masks inserted into it. The first lens converts the input beam into the Fourier (dual ) space, then a central phase mask, whose position defines the position of the system's Fourier plane, performs the transformation of the beam in the dual space, as per Eq. (2), and, eventually, the second lens converts the output beam back from the dual domain into the real (coordinate) one. Thus, the fractional diffraction is effectively executed in the Fourier representation of the field amplitude, and the linear version of Eq. (3) governs the circulation of the beam in the cavity, averaged over many cycles. As a result, the linearized version of Eq. (3) gives rise to conical diffraction of 1D and 2D input Gaussian beams [39,44].
The usual cubic nonlinearity can be incorporated in the setup by inserting a piece of a Kerr material in the cavity. To maintain the effectively local form of the nonlinearity, the material should be inserted between an edge mirror and the lens closest to it, where the beam's propagation takes place in the spatial domain (rather than in the dual one). The CQ terms may be similarly induced by a sample of an optical material combining the self-focusing and defocusing cubic and quintic nonlinearities, whose absorption is low enough to be neglected. This may be, as mentioned above, carbon disulfide in the liquid form, which was used for the creation of stable (2+1)-dimensional spatial solitons [80], quasi-stable (2+1)-dimensional optical solitons with embedded vorticity [81], and, very recently, of quasi-stable necklace arrays of optical beams with imprinted vorticity [20]. Another appropriate CQ material is a colloidal suspension of metallic nanoparticles, in which strengths of cubic and quintic nonlinearities may be accurately adjusted by selecting the size and concentration of the particles [82].
Proceeding to the consideration of the full nonlinear equation (3), we first address stationary solutions with a real propagation constant, β: where complex function ψ(ξ, η) obeys the equation Further, the complex function is substituted in the Madelung's form, ψ(ξ, η) = U (ξ, η)e iφ(ξ,η) , with real amplitude U (ξ, η) and phase φ(ξ, η). In polar coordinates (r, θ), related to the underlying Cartesian coordinates as usual, ξ = r cos θ, η = r sin θ, soliton solutions to Eq. (3) with integer embedded vorticity s are looked for as [76] The self-trapped fundamental states, with s = 0, are known as Townes solitons of the conventional 2D NLSE with the cubic-only term. In 1D, exact fundamental solitons of the equation with the quintic-only self-focusing term also play the role of solitons of the Townes type [83,84]. These solitons are unstable because the same equations admit the critical collapse [85]. In the case of the CQ nonlinearity, the fundamental solitons, in 2D and 1D settings alike, realize the system's ground state (which is missing in the presence of the collapse).
The self-trapping of the solitons implies that U (r) exponentially decays at r → ∞. Using the same argument as recently applied to the FSE with operator −∂ 2 /∂x 2 α/2 [71], viz., the consideration of the variation of the integral expression, which defines the nonlocal operator (2), with respect to rescaling of coordinates x and y, it is straightforward to conclude that a general form of the exponentially decaying tails of the soliton is where the positive constant C α is not universal, except for C α=2 = 1. Note also that the proportionality multiplier, dropped in Eq. (7), includes a power-law preexponential factor ∼ r −γα with some constant γ α > 0 (obviously, γ α=2 = 1/2). The power (alias integral norm) and z-component of the angular momentum of the vortex are where ρ = (ξ, η) is the vector of the integration coordinates, and e z is a unit vector in the propagation direction. The Hamiltonian of Eq. (3) is (10) Note that, with regard to the definition given by Eq. (2), expression (10) produces a real Hamiltonian [29]. Families of fundamental solitons with different values of the Lévy index α were numerically generated by dint of the Newton-conjugate-gradient method [87,88]. We display dependences of their propagation constant β on power P for different values of the Lévy index in Fig.  1. It is seen that two branches of the β(P ) curves merge at points with the vertical tangent, dβ/dP = 0. Top (dβ/dP > 0) and bottom (dβ/dP < 0) branches are, severally, stable and unstable, as established below by the linear-stability analysis and direct simulations. These conclusions are correctly predicted by the celebrated Vakhitov-Kolokolov criterion [85,89]. Note that this criterion is irrelevant for the stability of vortex solitons with s ≥ 1 against azimuthal instabilities that may split the ring-shaped vortex modes [76]. The point with dβ/dP = 0 corresponds to a threshold (minimum) value of the total power, P (s=0) th , necessary for the existence of the fundamental solitons. In the cases shown in Fig. 1 the numerically found threshold values for the fundamental solitons are P is common for models with competing nonlinearities, the solitons develop a flat-top shape in the limit of large powers.
The goal of this paper is to construct soliton clusters composed of fundamental solitons in free space (without an external potential). Before constructing the clusters, the stability of the fundamental solitons needs to be established, which can be done in the framework of the linearization of Eq. (3) for small perturbations. To this end, we searched for the perturbed solutions of Eq.
where ψ is the stationary solution with propagation constant β defined as per Eq. (4), the asterisk stands for the complex conjugation, while u (ξ, η) and v (ξ, η) are complex eigenmodes of perturbations, with an infinitesimal amplitude ǫ [91]. Further, δ ≡ δ R + iδ I is a complex eigenvalue of Eq. (12), δ R being the linear-instability growth rate. Obviously, the stability condition amounts to δ R = 0, for all eigenvalues. The substitution of expression (11) in Eq. (3) and subsequent linearization leads to the following eigenvalue problem: with matrix elements L 11 = − −∇ 2 ⊥ α/2 + 2 |ψ| 2 − 3 |ψ| 4 −β and L 12 = ψ 2 +2 |ψ| 2 ψ 2 . In the context of BEC, linearized equations for small perturbations are usually called Bogoliubov-de Gennes equations [86]. Note that the choice of the ansatz for the perturbed solution in the form of Eq. (11), with the combination of amplitudes u and v * , is convenient because it leads to Eq. (12) for two-component eigenfunctions which does not include complex conjugation.
The linear problem based on Eq. (12) was solved by means of the Newton-conjugate-gradient method [87,88]. Differently from the conventional NLSE with the CQ nonlinearity, corresponding to α = 2 in Eq. (1), in which the fundamental solitons are completely stable [22], the fundamental solitons are found to be unstable in intervals 0.006 < β ≤ 0.016 (for α = 1), see Fig. 2(a), 0.006 < β ≤ 0.030 (for α = 1.5) in Fig. 2(b), and 0.006 < β ≤ 0.011 (for α = 1.9) in Fig. 2(c). These results corroborate that the fundamental-soliton branches with a negative slope, dβ/dP < 0, are linearly unstable, while the positive-slope ones are stable, i.e., the stability of the fundamental solitons fully complies with the Vakhitov-Kolokolov criterion. In Figs. 2(d) and 2(e), the results of direct simulations demonstrate that unstable fundamental solitons at first keep their integrity, and then start to spread out, see also movies Visualization 1 and Visualization 2 in Supplemental Material. In Fig.  2(f), a linearly-stable fundamental soliton robustly propagates, under the action of random perturbations with a 5% relative amplitude, over the distance of ζ = 2000, which is tantamount to ∼ 50 diffraction lengths of this soliton.

III. THE INTERACTION ENERGY AND CONSTRUCTION OF SOLITON CLUSTERS
The necklace-shaped soliton cluster of radius R is constructed as a superposition of N identical fundamental solitons with a superimposed phase distribution carrying the integer global vorticity M . Thus, the initial field distribution is set as where U 0 is a stable fundamental soliton, and r n = {R cos(2πn/N ), R sin(2πn/N )} is the position of the center of the n-th soliton in the cluster. The evolution of the pattern is governed by interaction forces between adjacent solitons, which, in turn, depend on the phase shift and separation between them. Note that ansatz (13) includes a linear phase distribution with respect to the angular coordinate, θ, rather than a staircase phase structure, because the former one is more likely to form a quasi-stable structure [13].
As concerns the input necessary for the creation of soliton clusters introduced in the form of Eq. (13), in the experiment it can be put together as a set of several parallel beams placed along the ring, while a vortex phase plate may be used to imprint the phase pattern, which represents the overall vorticity [M in Eq. (13)], onto the multi-beam cluster.
Here, we emphasize the difference of physical properties between the vortex soliton and necklace-shaped soliton cluster. Vortex solitons in Ref. [76] are stationary solutions of Eq. (3), which may be stable or unstable. However, necklace-shaped soliton clusters, discussed in this paper, are not stationary solutions of Eq. (3). They are composed of stable fundamental solitons of Eq. (3), placed along a ring. While Eq. (3) produces no completely stable necklace-shaped soliton clusters, they can be found in a quasi-stable form. Their existence is predicted by minimization of the effective potential of interaction between adjacent solitons in the cluster, cf. Ref. [18].
Since the necklace may be considered as an azimuthally modulated vortex soliton, one may expect that necklace clusters can be maintained by Eq. (3) if their integral power exceeds the known threshold value, P (s=1) th [76], above which this equation gives rise to stable vortex solitons. This circumstance is elaborated below, see Eq. (16). The force of the interaction between fundamental solitons with phase difference χ, separated by large distance L, is determined by an effective interaction potential, which is determined by the integral of overlapping between an exponentially decaying tail of each soliton and the core of its mate. Taking into regard the general asymptotic form of the tails, given by Eq. (7), we conclude that the potential of the soliton-soliton interaction may be estimated, as per Ref. [90], as where W 0 is a positive factor, which may be a slowly decaying function of L (for α = 2, W 0 = const · L −1/2 ), while factor cos χ is a universal one [90]. For the ansatz defined by Eq. (13), one has L = 2R sin (π/N ), χ = 2πM/N . For the whole necklace cluster, the scaled potential, normalized to the energy of the same set of infinitely separated (non-interacting) solitons, may be defined as where H is the full Hamiltonian (10). It can be computed numerically as a function of R for given vorticity M and number of individual " bids in the necklace" , i.e., individual fundamental solitons, N , by substituting ansatz (13) in integral expression (10). First, we address the case of α = 1.5. In this case, Eq. (3) maintains stable solitons with vorticity s = 1 if the integral power exceeds the above-mentioned threshold value, P (s=1) th ≈ 383, see Table 1 in Ref. [76]. Therefore, a set of fundamental solitons, each taken (for instance) with P sol = 84.1 and β = 0.12, may have a chance to build robust clusters with M = 1 and which actually implies N ≥ 5. In Fig. 3, the dependence of the effective potential (15) on the cluster's radius R is presented for different integer values of M and N . These results reveal that the necklace clusters with vorticities M = 0 and 1 exhibit visible local minima of the interaction potential at a specific value of the cluster's radius, R min . These numerically identified values for M = 0 and 1 are summarized in Table I for different numbers of solitons in the cluster, N . The results demonstrate that R min increases with the increase of N . For the clusters with M = 2, local energy minima are observed in Figs. 3(e,f) for N = 9 and 10, but, according to data from Table 1 in Ref. [76], the respective values of the total power of initial cluster (13) are lower than the corresponding stability-threshold value for the axisymmetric vortex soliton, cf. Eq. (16). Therefore, the minima of the potential energy for M = 2 cannot produce stable necklace patterns. Lastly, Fig. 3 does not produce local energy minima for the clusters with M = 3.
To check the predictions produced by the effective potential, we numerically simulated the propagation of the necklace clusters built of N = 5 fundamental solitons [see Fig. 3(a)] by means of the split-step Fourier method. In the course of the simulations, the evolution of the cluster's mean radius R (ζ), mean angular velocity ω (ζ), and depth of the azimuthal modulation Θ (ζ) was monitored. These quantities were defined as follows: and where |Ψ| 2 max and |Ψ| 2 min are the maximum and minimum values of |Ψ| 2 , as functions of θ, along the circumference at which the solution has a radial maximum in each cluster's configuration.
Results produced by the simulations are presented in Fig. 4. For M = 0, interactions between adjacent fundamental solitons are attractive, according to Eq. (14), therefore the multi-soliton cluster fuses into a fundamental state, as shown in the top row of Fig. 4, even for the initial radius taken as R = R min . On the other hand, for M = 2, the soliton-soliton interactions, with χ = π, are repulsive, thus giving rise to a resulting repulsive force acting on each individual soliton in the radial direction. For this reason, the effective potential energy (15) has no local minimum in this case, see Fig. 3(a). As a result, the cluster expands and rotates, as seen in Fig. 4(b). For M = 1, i.e., χ = 3π/5, the interaction potential (14) is weakly attractive, but the net angular momentum of the circular cluster prevents fusion into an axisymmetric vortex soliton. In this case, the cluster with initial radius R > R min performs quasiperiodic cycles of expansion and shrinkage (in combination with rotation), as shown in Fig. 4(c). The most essential finding is that the cluster with initial radius R ≈ R min gives rise a quasi-stationary state, which performs persistent rotation with minimal radial oscillations, see Fig. 4(d). This mechanism of the effective stabilization of the rotating necklace cluster is somewhat similar to that recently reported in the framework of the usual 2D NLSE with the nonlinearity corresponding to quantum droplets [18]. In Fig. 5, the evolution of the cluster's mean radius, cluster's mean angular velocity, and depth of the azimuthal modulation is monitored, to display further details of the variation of the cluster's structure and rotational speed in the course of the propagations. The first panel in the middle row of Fig. 5 shows the case of the cluster's mean angular velocity equal to zero, therefore it does not rotate. In the second panel in the middle row of Fig. 5, the cluster exhibits rotation with a decreasing angular velocity. In the last two panels in the middle row of Fig.  5, the cluster's mean angular velocity fluctuates around non-decaying mean values, which indicates that the re-spective clusters perform persistent rotation with small radial oscillations in the transverse plane. To further elucidate the robustness of the soliton clusters with the properly chosen initial radius, R = R min , we added random perturbation to the initial profiles of the soliton cluster in Eq. (13), multiplying it by where ρ 1 and ρ 2 are uniformly distributed random numbers in the interval of [−0.5, 0.5]. In Fig. 6, the simulations demonstrate that the rotating necklace clusters, composed of different numbers of fundamental solitons (N = 5, 6, 7, and 8), with initial radius R = R min , can survive, as metastable patterns, in the course of long evolution, even under the action of relatively strong random perturbations. In the Supplemental Material, we additionally display the evolution of the clusters with vortic- ity M = 1 and initial radius R = R min , for larger numbers of the fundamental solitons, viz., N = 9 and 10, see movies Visualization 3 and Visualization 4, respectively.
The results indicate that the rotating necklace clusters remain robust against the action of random perturbation, although exhibiting radial oscillations. In the course of the propagations, the evolution of the cluster's mean radius and mean angular velocity quasi-periodically oscillate with a small amplitude, see the top and middle rows in Fig. 7. The depth of the azimuthal modulation is also monitored, see the bottom row in Fig. 7, which indicates that the clusters do not fully fuse into the usual solitons in the course of the long-distance propagation. Next, we address the necklace-shaped clusters composed of the fundamental solitons with Lévy index α = 1 [recall it is the critical value for Eq. (3) with the cubiconly self-attraction, therefore it is relevant to address this case too]. It was recently found that the axisymmetric vortex soliton (not a cluster) with s = 1 and α = 1 is stable if its integral power exceeds the respective threshold value, P > P (s=1) th ≈ 1553 (see Table 1 in Ref. [76]). Here, we construct the necklace patterns with vorticity M = 1 as circular sets of N stable fundamental solitons, each one with power P sol = 332.35 and propagation constant β = 0.105. Thus, they have a chance to be quasi-stable for N > P (s=1) th /P sol ≈ 4.7 [cf. Eq. (16) for α = 1.5], i.e., as a matter of fact, for N ≥ 5. In Fig. 8, we show the dependence of the effective potential energy of the necklace clusters with vorticity M = 1, computed as per Eq. (15), on the cluster's radius for α = 1. The energy has a local minimum at specific values of the radius, viz., R min (N = 5) = 19.8, R min (N = 6) = 23.55, R min (N = 7) = 27.6, and R min (N = 8) = 31.6, respectively. Robust propagation of the rotating clusters at parameters corresponding to the energy minima identified in Fig. 8 is confirmed by direct simulations performed in the presence of perturbations produced by the same factor (20), as used above. As shown in Fig. 9, the perturbed propagation is quasi-stable, with small radial oscillations (cf. Fig. 6). In Fig. 10, we also display the cluster's mean radius, the mean angular velocity, and the depth of the azimuthal modulation as functions of ζ in the course of the propagation.

V. CONCLUSION
We have investigated ring-shaped soliton clusters (necklaces) in the framework of the NLSE, which includes the fractional diffraction, characterized by the respective Lévy index, 1 ≤ α < 2, and the competing cubic-quintic nonlinearity. Fundamental and axisymmetric vortical solitons have their stability and instability regions in this model, the stability of the fundamental solitons being exactly predicted by the Vakhitov-Kolokolov criterion. Computing the effective interaction potential for the multi-soliton clusters, the analysis predicts its minima at particular values of the necklace's radius, at which attractive interactions between adjacent fundamental solitons forming the cluster are balanced by the net angular momentum. Such rotating metastable soliton clusters withstand the action of relatively strong random perturbations, featuring long-scale quasi-stable propagation.
These findings suggest other interesting scenarios for the generation of complex nonlinear states in fractional dimensions. In particular, it is relevant to extend the analysis for vector solitons in two-and many-component settings, as well as for a PT -symmetric system generalizing the present one.

VII. DISCLOSURES
The authors declare no conflicts of interest.