Stable patterns in the Lugiato-Lefever equation with a confined vortex pump

We introduce a model of a passive optical cavity based on a novel variety of the two-dimensional Lugiato-Lefever equation, with a localized pump carrying intrinsic vorticity S, and the cubic or cubic-quintic nonlinearity. Up to S = 5, stable confined vortex-ring states (vortex pixels) are produced by means of a variational approximation and in a numerical form. Surprisingly, vast stability areas of the vortex states are found, for both the self-focusing and defocusing signs of the nonlinearity, in the plane of the pump and loss parameters. When the vortex-rings are unstable, they are destroyed by azimuthal perturbations which break the axial symmetry. The results suggest new possibilities for mode manipulations in passive nonlinear photonic media by means of appropriately designed pump beams.


I. INTRODUCTION AND THE MODEL
Optical solitons are a broad class of self-trapped states maintained by the interplay of nonlinearity and dispersion or diffraction in diverse photonic media [1,2].In addition to that, dissipative optical solitons are supported by the equilibrium of loss and gain or pump, concomitant to the nonlinearity-dispersion/diffraction balance [3,4].Dissipative solitons have been studied in detail, theoretically and experimentally, in active setups, with the loss compensated by local gain (essentially, provided by lasing), being modelled by one-and two-dimensional (1D and 2D) equations of the complex Ginzburg-Landau (CGL) type [5,6].
In passive nonlinear optical cavities, the losses are balanced by the pump field supplied by external laser beams, with the appropriate models provided by the Lugiato-Lefever (LL) equations [7].This setting was also studied in the 1D and 2D forms [8][9][10][11].Widely applied in nonlinear optics, equations of the LL type play a crucial role in understanding fundamental phenomena such as the modulation instability (MI) and pattern formation in dissipative environments [8]- [24].The relevance of these models extends to the exploration of complex dynamics of various nonlinear photonic modes, tremendously important applications being the generation of Kerr solitons and frequency combs in passive cavities [12]- [24], as well as the generation of terahertz radiation [25].In addition to rectilinear cavity resonators, circular ones can be used too [26].In many cases, they operate in the whispering-gallery regime [27]- [32].
In most cases, solutions of the one-and two-dimensional LL equations are looked for under the action of the spatially uniform pump, which approximately corresponds to the usual experimental setup.However, the use of localized (focused) pump beams is possible too, which makes it relevant to consider LL equations with the respective shape of the pump terms.In fact, truly localized optical modes in the cavities can be created only in this case, otherwise the uniform pump supports nonzero background of the optical field.In particular, exact analytical solutions of the LL equations with the 1D pump represented by the delta-function, and approximate solutions maintained by the 2D pump in the form of a Gaussian were reported in Ref. [33].In Ref. [26], the LL equation for the ring resonator with localized pump and loss terms produced nonlinear resonances leading to multistability of nonlinear modes and coexisting solitons, that are associated with spectrally distinct frequency combs.
Further, solutions for fully localized robust pixels with zero background were produced by the 2D LL equation, incorporating the spatially uniform pump, self-focusing or defocusing cubic nonlinearity, and a tight confining harmonic-oscillator potential [34].Additionally, this model with a vorticity-carrying pump gives rise to stable vortex pixels.In particular, in the case of the self-defocusing sign of the nonlinear term, the pixels with zero vorticity and ones with vorticity S = 1 were predicted analytically by means of the Thomas-Fermi (TF) approximation.
In this work, we introduce the 2D LL equation for complex amplitude field u (x, y, t) of the light field, with the cubic or cubic-quintic nonlinearity: and a confined pump beam represented by factor f (r).Here i is the imaginary unit, α > 0 is the loss parameter, σ = +1 and −1 corresponds, respectively, to the self-focusing and defocusing Kerr (cubic) nonlinearity, g > 0 or g < 0 represents the self-defocusing or focusing quintic nonlinearity (which often occurs in optical media [35][36][37], in addition to the cubic term), parameter η defines the cavity mismatch, which is ση 2 − gη 4 (the coefficient multiplying the linear term ∼ iu), in terms of the linearized LL equation, and written in terms of polar coordinates (r, θ), corresponds to the confined pump beam with real amplitude f 0 , radial width W, and integer vorticity S ≥ 1. Vortex beams, shaped by the passage of the usual laser beam through an appropriate phase mask, are available in the experiment [37].Equation ( 1) is written in the scaled form.All figures are plotted below in the same notation.In physical units, r = 1 and t = 1 normally correspond to ˜50 ¯m and ∼ 50 ps, respectively.Then, the typical width W = 2, considered below, corresponds to the pump beam with diameter ∼ 100 ¯m, which an experimentally relevant value.Accordingly, the characteristic evolution time in simulations presented below, t ∼ 100, corresponds to times ˜5 ns.
Stationary solutions of Eq. ( 1) are characterized by values of the total power (alias norm), and angular momentum, (with * standing for the complex conjugate), even if the power and angular momentum are not dynamical invariants of the dissipative equation ( 1).In the case of the axisymmetric solutions with vorticity S, i.e., u (x, y) = u(r)e iSθ [38][39][40], the expressions for the power and angular momentum are simplified: Our objective is to construct stable ring-shaped vortex solitons (representing vortex pixels, in terms of plausible applications), as localized solutions of Eq. ( 1) with the same S as in the pump term (2).The stability is a challenging problem, as it is well known from the work with models based on the nonlinear Schrödinger and CGL equations that (in the absence of a tight confining potential) vortex-ring solitons are normally vulnerable to the splitting instability.In the case of a narrow ring shape, the splitting instability may be considered as quasi-one-dimensional MI of the ring against azimuthal perturbations which break its axial symmetry [2,40].The azimuthal MI is driven by the self-focusing nonlinearity, and inhibited by the self-defocusing.
To produce stationary solutions for the vortex solitons in an approximate analytical form (parallel to the numerical solution), we employ a variational approximation (VA).Our results identify regions of the existence and stability of the vortex solitons with 1 ≤ S ≤ 4 in the space of parameters of Eqs.(1) and (2) (in particular, in the plane of ( f 0 , α)) for both signs of the cubic nonlinearity, σ = ±1, while the mismatch parameter is fixed to be η = 1 by dint of scaling.The stability areas are vast, provided that the loss coefficient α is, roughly speaking, not too small.A majority of the results are produced for the pure cubic model, with g = 0, but the effect of the quintic term, with g ̸ = 0, is considered too.Quite surprisingly, a stability area for the vortices with S ≤ 3 is found even in the case of σ = +1, g < 0, when both the cubic and quintic terms are self-focusing, which usually implies strong propensity to the azimuthal instability of the vortex rings [40].
The rest of the paper is structured as follows.The analytical approach, based on the appropriate VA, is presented in Sec.II.An asymptotic expression for the tail of the vortex solitons, decaying at r → ∞, is found too in that section.Systematically produced numerical results for the shape and stability of the vortex solitons are collected (and compared to the VA predictions) in Sec.III.The paper is concluded by Sec.IV.

II. ANALYTICAL CONSIDERATIONS
This section summarizes the analytical part of the work and results produced by this part.Two directions of the analytical considerations for the present model are possible: the investigation of the decaying "tails" of the localized stationary states, in the framework of the linearized model, and detailed development of VA for the full model, including the nonlinear terms.

A. Asymptotic forms of the vortex solitons
Direct consideration of the linearized version of Eqs. ( 1) and (2) readily produces an explicit result for the soliton's tail decaying at r → ∞: with the power of the pre-Gaussian factor, r S−2 , which is lower than that in the pump term, r S .Due to this feature, the asymptotic expression (6) formally predicts a maximum of local power |u(r)| 2 at S > 2, at A local maximum is indeed observed in numerically found radial profiles of all vortex solitons (see Figs. 2, 4, 7(b), 9(b), 10, and 11 below, for S = 1, 4, 2, 2, 1, and 3, respectively).In fact, for these cases Eq. ( 7) predicts values of r max which are smaller by a factor ≃ 0.6 than the actually observed positions of the maxima.The discrepancy is explained by the fact that the asymptotic expression ( 6) is valid at values of r which are essentially larger than r max .In a looser form, one can try to construct an asymptotic approximation for the solution at moderately large r, by adopting the ansatz which follows the functional form of the pump term (2), viz., where u(r) is a complex slowly varying function, in comparison with those which are explicitly present in ansatz (8).Substituting the ansatz in Eq. ( 1) and omitting derivatives of the slowly varying function, one can develop an approach which is akin to the TF approximation applied to the model with the tight trapping potential in Ref. [33].
The result for the linearized version of Eq. ( 1), which implies a small amplitude of the mode pinned to the pump beam, is where, as said above, η = 1 is substituted.In the limit of r → ∞, Eqs. ( 9) and ( 8) carry over into the asymptotically rigorous expression (6).On the other hand, Eq. ( 9) predicts a maximum of the local power at cf. Eq. ( 7).One may expect that the prediction of the local maximum of the vortex soliton at point ( 10) is valid when it yields values of r 2 max TF which are large enough, i.e., if S and W are relatively large.Indeed, the comparison with the numerically produced profiles of the vortex solitons, displayed below in Figs. 6 and 10, demonstrates that Eq. ( 10) predicts, relatively accurately, (r max ) TF = 4 for S = 5, W = 2, σ = −1, and g = 0.However, the prediction given by Eq. ( 10) is not accurate for S = 1 and 2.
Finally, in the limit of r → 0, the asymptotic form of the solution is simple, u(r, θ) ≈ u 0 r S , but constant u 0 cannot be found explicitly, as it depends on the global structure of the vortex-soliton solution.In particular, the crude TF approximation given by Eq. ( 9) yields .

B. The variational approximation (VA)
A consistent global analytical fit for the vortex solitons may be provided by VA, based on the Lagrangian of the underlying equation [2].While Eq. ( 1), which includes the linear dissipative term, does not have a Lagrangian structure, it can be converted into an appropriate form by the substitution suggested by Ref. [41], which absorbs the dissipative term: producing the following time-dependent equation for complex function U(r, t), where, as said above, we set η = 1 by means of scaling: The real Lagrangian which precisely produces the time-dependent equation ( 12) is The simplest ansatz which may be used as the basis for VA follows the form of the pump term (2): where variational parameters U 0 and ϕ are the real amplitude and phase shift of the solution with respect to the pump.Power (3) for this ansatz is where Γ(S + 1) ≡ S! is the Gamma-function, and the time-dependent factors, exp (±αt), mutually cancel when relations ( 11) and ( 14) are substituted in expression (3).Note that the local power |U(r)| 2 , corresponding to ansatz (14), attains it maximum at r 2 = SW 2 /2.The substitution of ansatz (14) in Lagrangian ( 13) and straightforward integration yields the respective VA Lagrangian, Then, the Euler-Lagrange equations for U 0 and ϕ are obtained as where δ/δϕ stands for the variational derivative.Taking into regard that Lagrangian ( 16) must be substituted in the respective action, L VA dt, and then the action must be actually subjected to the variation, one should apply the time differentiation to factor e 2αt in Lagrangian (16), while deducing the appropriate form of δL VA /δϕ in Eq. (17).Once the Euler-Lagrange equations (17) have been derived , we consider their stationary (fixed-point) solutions by setting dU 0 /dt = dϕ/dt = 0, which yields It is relevant to mention that the evolution equation for the power (3), which follows from Eq. ( 12), is The stationary states must satisfy the balance condition, dP/dt = 0.Then, the substitution of ansatz ( 14) and expression (2) for the pump in this condition yields a simple relation, which is identical to Eq. ( 18).In particular, Eq. ( 21) implies that, for the fixed pump's amplitude f 0 , the amplitude of the established localized pattern cannot exceed the maximum value, which corresponds to ϕ = 0 in Eq. ( 21): FIG.
1.The power of the VA solution, in the case of g = 0 (no quintic nonlinearity), vs. the pump's squared width, as given by Eq. ( 24), which neglects weak effects of the pump and loss (small f 0 and α), for the self-focusing (σ = +1) in (a) and defocusing (σ = −1) in (b) signs of the cubic term.The black curve with circles, the red one with squares, and the blue one with triangles pertain to vorticities S = 1, 2, and 3, respectively.
C. VA for the cubic (g = 0) and quintic (g → ∞) models First, we aim to predict stationary states, as solutions of Eqs. ( 18) and ( 19), for the model with cubic-only nonlinearity, i.e., g = 0, under the assumption that the loss and pump terms in Eq. ( 1) may be considered as small perturbations.In the lowest approximation, i.e., dropping the small term ∼ f 0 in Eq. ( 19), one obtains a relatively simple expression which predicts the squared amplitude of the vortex soliton, the respective power (15) of the underlying ansatz (14) being Note that expressions (23) and ( 24) are always meaningful for the self-focusing sign of the cubic nonlinearity, σ = +1, while in the case of defocusing, σ = −1, the expressions are meaningful if they are positive, which imposes a restriction on the width of the Gaussian pump: it must be broad enough, viz., The dependence of the power given by Eq. ( 24) on the pump's squared width W 2 for three different values of the vorticity, S = 1, 2, 3, is plotted in Figs.1(a) and (b) for the self-focusing and defocusing signs of the cubic term, i.e., σ = +1 and −1, respectively.Note that in Fig. 1(b) for σ = −1, there is no solution in the region in which condition (25) does not hold, hence the VA solution does not exist.
In the limit of the dominant quintic nonlinearity, i.e., g → ±∞, opposite to the pure cubic model considered above, the asymptotic solution of Eq. ( 19) is (in which the large coefficient g cancels out), the respective expression for power (15) being cf. Eqs. ( 23) and ( 24).
In the following section, we report results of numerical solution of Eq. ( 1), comparing them to solutions of the full system of the VA equations ( 18) and (19), which include effects of the pump and loss terms.

III. NUMERICAL RESULTS
Simulations of Eq. ( 1) were conducted by means of the split-step pseudo-spectral algorithm.The solution procedure started from the zero input, and was running until convergence to an apparently stable stationary profile (if this outcome of the evolution was possible).This profile was then compared to its VA counterpart, produced by a numerical solution of Eqs. ( 18) and (19) with the same values of parameters α, σ = ±1, g, and f 0 , W, S (see Eq. ( 2)).The results are presented below by varying, severally, loss α, vorticity S, the pump's width W and strength f 0 , and, eventually, the quintic coefficient g.The findings are eventually summarized in the form of stability charts plotted in Fig. 12.

A. Variation of the loss parameter α
In Fig. 2(a) we display the cross-section (drawn through y = 0) of the variational and numerical solutions for the stable vortex solitons obtained with α = 0.5, 1.0, and 2.0, while the other parameter are fixed as σ = +1 (the self-focusing cubic nonlinearity), g = 0, f 0 = 1, W = 2, and S = 1.The accuracy of the VA-predicted solutions presented in Fig. 2(a) is characterized by the relative power difference from their numerically found counterparts, which is 5.1%, 3.1%, and 0.5% for α = 0.5, 1.0, 2.0, (28) respectively.Thus, the VA accuracy improves with the increase of α.
Similar results for the self-defocusing nonlinearity, σ = −1, are presented in Fig. 2(b), which shows an essentially larger discrepancy between the VA and numerical solutions, viz., 18.9%, 17.4%, and 3.4% for the same set (28) of values of the loss parameter, other coefficients being the same as in Fig. 2(a).The larger discrepancy is explained by the fact that localized (bright-soliton) modes are not naturally maintained by the self-defocusing, hence the ansatz (14), which is natural for the self-trapped solitons in the case of the self-focusing, is not accurate enough for σ = −1.
In the same vein, it is natural that, in the latter case, the discrepancy is more salient for stronger nonlinearity, i.e., smaller α, which makes the respective amplitude higher.
There is a critical value of α below which the vortex solitons are unstable.As an example, Fig. 3 shows the VApredicted and numerically produced solutions for α = 0.2 and σ = +1 (the self-focusing nonlinearity).The observed picture may be understood as a result of the above-mentioned azimuthal MI which breaks the axial symmetry of the vortex soliton.More examples of the instability of this type are displayed below.For the values of other parameters fixed as in Fig. 3, the instability boundary is α crit ≈ 0.35.The stabilizing effect of the loss at α > α crit is a natural feature.On the other hand, the increase of α leads to decrease of the soliton's amplitude, as seen in Fig. 2.
In the case of the self-defocusing (σ = −1), all the numerically found vortex modes are stable, at least, at α ≥ 0.1, although the discrepancy in the values of the power between these solutions and their VA counterparts is very large at small α, exceeding 75% at α = 0.1.As mentioned above, the growing discrepancy is explained by the increase of FIG. 3. Profiles of |u| 2 produced by VA (a) and numerical solution at t = 100 (b) in the case of the self-focusing (σ = +1) for α = 0.2.Other parameter are the same as in Fig. 2(a).the soliton's amplitude with the decrease of α.At still smaller values of α, the relaxation of the evolving numerical solution toward the stationary state is very slow, which makes it difficult to identify the stability.

B. Variation of the pump's vorticity S
To analyze effects of the winding number (vorticity) S, we here fix g = 0 (the pure cubic nonlinearity) and set α = 1, f 0 = 1, W = 2 in Eqs. ( 1) and ( 2).In the self-focusing case (σ = +1), the numerically produced solutions are stable for S = 1 and 2, and unstable for S ≥ 3.In the former case, the power difference between the VA and numerical solutions is 3.1% and 4.9% for S = 1 and 2, respectively, i.e., the VA remains a relatively accurate approximation in this case.
In the self-defocusing case (σ = −1), considering the same values of the other parameters as used above, the numerical solution produces stable vortex solitons at least until S = 5.For the same reason as mentioned above, the accuracy of VA is much lower for σ = −1 than for the self-focusing case (σ = +1), with the respective discrepancies in the power values being 17.4%, 6.7%, 24.1%, 29.0%, and 29.6% for S = 1, 2, 3, 4, and 5, respectively.
For the cogent verification of the stability of the localized vortices in the case of the self-defocusing, we have also checked it for smallest value of the loss parameter considered in this work, viz., α = 0.1, again for S = 1, 2, 3, 4, and 5, and the above-mentioned values of the other coefficients, i.e., f 0 = 1 and W = 2. Naturally, the discrepancy between the VA and numerical findings is still higher in this case, being 75%, 72.5%, 65.7%, 55.2%, and 41.1% for S = 1, 2, 3, 4, and 5, respectively.The result is illustrated by Fig. 4 for a relatively large vorticity, S = 4.In particular, the pattern of |u(x, y)| 2 and the corresponding cross section, displayed in Figs.4(b) and (c), respectively, exhibit the established vortex structure and background "garbage" produced by the evolution.

C. Variation of the pump's width W
To address the effects of the variation of parameter W in Eq. ( 2), we here fix g = 0, α = 1, f 0 = 1, and S = 1.In Fig. 5(a), the power of the VA-predicted and numerically found stable vortex-soliton solutions is plotted as a function of W for both the self-focusing and defocusing cases, i.e., σ = +1 and σ = −1, respectively.In the former case, the azimuthal MI sets in at W ≥ 2.7, see an example in Fig. 5(b) for W = 2.75.
In the self-focusing case, σ = +1, the azimuthal MI for the solitons with higher vorticities, S = 2, 3, 4, or 5, sets in at W ≥ 2.1, 1.7, 1.6, and 1.4, respectively.In the self-defocusing case, no existence/stability boundary was found for the vortex modes with S = 1, S = 2, and S = 3 (at least, up to W = 5).At higher values of the vorticity, the localized vortices do not exist, in the defocusing case, at W > 3.5 and W > 2.5, for S = 4, and S = 5, respectively.

D. Variation of the pump's strength f 0
Effects of the variation of f 0 are reported here, fixing other parameters as g = 0, α = 1, and W = 2.In the case of the cubic self-focusing, σ = +1, the vortex soliton with S = 1 are subject to MI at f 0 ≥ 1.6.As a typical example, in Fig. 6 we display the VA-predicted solution alongside the result of the numerical simulations for f 0 = 1.7.For higher vorticities, S = 2, 3, 4, and 5, the azimuthal instability sets in at f 0 ≥ 1.1, 0.6, 0.3, and 0.08, respectively.Naturally, the narrow vortex rings with large values of S are much more vulnerable to the quasi-one-dimensional azimuthal MI.
The power of the vortex solitons with S = 1 and 2, as produced by the VA and numerical solution, is plotted vs. the pump amplitude f 0 in Fig. 7(a).As an example, Fig. 7(b) showcases an example of the cross-section profile of the vortex soliton with S = 2, demonstrating the reliability of the VA prediction.In the range of f 0 ≤ 1, the highest relative difference in the power between the numerical and variational solutions cases is 5.5% and 7.5% for S = 1 and S = 2, respectively.
It is relevant to mention that the "traditional" azimuthal instability of vortex-ring solitons with winding number S demonstrates fission of the original axially symmetric shape into a set consisting of a large number N ≥ S of symmetrically placed localized fragments [2,40]  In self-defocusing case, σ = −1, a summary of the results produced by the VA and numerical solution for the stable vortex solitons with S = 1 and 2, in the form of the dependence of their power on f 0 , is produced in Fig. 9(a) (cf.Fig. 7(a) for σ = +1).Naturally, the VA-numerical discrepancy increases with the growth of the pump's strength, f 0 , see an example in Fig. 9(b).Unlike the case of σ = +1, in the case of the self-defocusing the vortex modes with S ≤ 5 remain stable, at least, up to f 0 = 5 (here, we do not consider the case of S > 5).

E. Influence of the quintic coefficient g
In the above analysis, the quintic term was dropped in the LL equation (1), setting g = 0. To examine the impact of this term, we first address the case shown above in Fig. 3, which demonstrated that the vortex soliton with S = 1, as a solution to Eqs. ( 1) and ( 2) with g = 0, σ = +1 and f 0 = 1, W = 2, was unstable if the loss parameter fell below the critical value, α = 0.35.We have found that, adding to Eq. ( 1) the quintic term with either g = −1 or g = +1 (the self-focusing or defocusing quintic nonlinearity, respectively) leads to the stabilization of the vortex mode displayed in Fig. 3, which was unstable in the absence of the quintic term.The stabilization of the soliton by the quintic self-defocusing is a natural fact.More surprising is the possibility to provide the stabilization by the self-focusing quintic nonlinearity because, in most cases, the inclusion of such a term gives rise to the supercritical collapse in 2D, making all solitons strongly unstable [2,40].However, it is concluded from the stability charts displayed below in Fig. 12 that the stabilizing effect of the quintic self-focusing occurs only at moderately small powers, for which the the quintic term is not a clearly dominant one.In the general case, solitons stability regions naturally shrink in Figs.
12 under the action of the quintic self-focusing quintic term, with g < 0. In the presence of the quintic term, the comparison of the numerically found stabilized vortex-soliton profiles with their VA counterparts, whose parameters are produced by a numerical solution of Eqs. ( 18) and ( 19), is presented in Fig. 10.Similar to the results for the LL equation with the cubic-only nonlinearity, the VA is essentially more accurate in the case of the self-focusing sign of the quintic term (g < 0) than in the opposite case, g > 0. In particular, in the case shown in Fig. 10, the power-measured discrepancy for g = −1 and +1 is, respectively, 4% and 64%.An explanation for this observation is provided by the fact that the soliton's amplitude is much higher in the latter case.
Another noteworthy finding is the stabilization of higher-vorticity solitons by the quintic term.For instance, it was shown above that, for parameters α = 1, η = 1, f 0 = 1, W = 2, and σ = +1 (the cubic self-focusing) all vortex solitons with S ≥ 3, produced by Eq. (1), were unstable.Now, we demonstrate that the soliton with S = 3 is stabilized by adding the self-defocusing quintic term with a small coefficient, just g = 0.1, see Fig. 11(b).As a counter-intuitive effect, the stabilization of the same soliton by the self-focusing quintic term is possible too, but the necessary coefficient is large, g = −6, see Fig. 11(a) (recall that g = −1 is sufficient for the stabilization of the vortex soliton with S = 1 and α = 0.2 in Fig. 10(a)).Nevertheless, similar to what is said above, the set of Figs.12(c,f,i) demonstrates the natural shrinkage of the stability area under the action of the quintic self-focusing.For the stabilized vortex modes shown in Fig. 11(a), in the case when both the cubic and quintic terms are self-focusing, the relative power-measured discrepancy between the numerical and VA-predicted solutions is very small, ≈ 0.7%, while in the presence of the weak quintic self-defocusing in Fig. 11(b) the discrepancy is 6.6%.

F. Stability charts in the parameter space
The numerical results produced in this work are summarized in the form of stability areas plotted in Fig. 12 in the parameter plane ( f 0 , α), for the vortex-soliton families with winding numbers S = 1, 2, 3, and three values of the quintic coefficient, g = −1, 0, +1, while the cubic term is self-focusing, σ = +1, and the width of the pump beam is fixed, W = 2.In addition to that, stability charts corresponding to the combination of the cubic self-defocusing (σ = −1) and quintic focusing (g = −1), also for S = 1, 2, 3, are plotted in Fig. 13.
The choice of the parameter plane ( f 0 , α) in the stability diagrams displayed in Fig. 12 is relevant, as the strength of the pump beam, f 0 , and loss coefficient, α,are amenable to accurate adjustment in the experiment (in particular, α may be tuned by partially compensating the background loss of the optical cavity by a spatially uniform pump, taken separately from the confined pump beam).As seen in all panels of Figs. 12 and 13, the increase of α naturally provides effective stabilization of the vortex modes, while none of them may be stable at α = 0, in agreement with the known properties of vortex-soliton solutions of the 2D nonlinear Schrödinger equation with the cubic and/or cubicquintic nonlinearity [2,40].The apparent destabilization of the vortices with the increase of the pump's amplitude f 0 is explained by the ensuing enhancement of the destabilizing nonlinearity.Other natural features exhibited by Figs.
12 are the general stabilizing/destabilizing effect of the quintic self-defocusing/focusing (as discussed above), and expansion of the splitting-instability area with the increase of S (the latter feature is also exhibited by Fig. 13).The latter finding is natural too, as larger S makes the ring-shaped mode closer to the quasi-1D shape (see, in particular, Fig. 8), which facilitates the onset of the above-mentioned azimuthal MI (modulational instability).
Thus, the inference is that the instability mode which determines the boundary of the stability areas in Figs. 12 and 13 is the breaking of the axial symmetry of the vortex rings by azimuthal perturbations, as shown above, in particular, in Figs.3(b), 5(b), and 6(b).The destabilization through the spontaneous splitting of the rings into symmetric sets of fragments (see Fig. 8) occurs deeply inside the instability area, i.e., at larger values of f 0 .

IV. CONCLUSION
We have introduced the two-dimensional LL (Lugiato-Lefever) equation including the self-focusing or defocusing cubic or cubic-quintic nonlinearity and the confined pump with embedded vorticity (winding number), S ≤ 5. Stable states in the form of vortex solitons (rings) for these values of S are obtained, in parallel, in the semi-analytical form by means of the VA (variational approximation) and numerically, by means of systematic simulations of the LL equation starting from the zero input.The VA provides much more accurate results in the case of the self-focusing nonlinearity than for the defocusing system.Stability areas of the vortex solitons with S = 1, 2, 3 are identified in the plane of experimentally relevant parameters, viz., the pump amplitude and loss coefficient, for the self-focusing and defocusing signs of the cubic and quintic terms.Stability boundaries for the vortex rings are determined by the FIG.12. Stability areas for families of the vortex solitons with winding numbers S = 1, 2, and 3, in the plain of the loss coefficient (α) and amplitude of the pump beam ( f 0 ), for three different values of the quintic coefficient, g = −1, 0, +1 (recall that g < 0 and g > 0 correspond, respectively, to the self-focusing and defocusing).Other parameters of Eqs. ( 1) and ( 2) are σ = +1 (the self-focusing cubic term), η = 1, and W = 2. FIG. 13.The same as in Figs.12(a-c), but for σ = −1 and g = −1 (the cubic self-defocusing and quintic focusing terms).onset of the azimuthal instability which breaks their axial symmetry.These findings suggest new possibilities for the design of tightly confined robust optical modes, such as vortex pixels.
As an extension of this work, it may be interesting to construct solutions pinned to a symmetric pair of pump beams with or without intrinsic vorticity.In this context, it is possible to consider the beam pair with identical or opposite vorticities.In the case of the self-focusing sign of the nonlinearity, one may expect onset of spontaneous breaking of the symmetry in the dual-pump configuration.Results for this setup will be reported elsewhere.

FIG. 4 .
FIG. 4. (a) The VA-predicted pattern and (b) the corresponding result of the direct simulation of Eq. (1) with σ = −1 and g = 0 (cubic self-defocusing) at t = 100, initiated by the zero input at t = 0, for α = 0.1, f 0 = 1, W = 2, and vorticity S = 4 in the pump term (2).(c) The respective cross sections drawn through y = 0. (d) The evolution of the total power P (see Eq. (3) of the numerical solution in the course of the simulation.
, while the above examples, displayed in Figs.3(b), 5(b), and 6(b), demonstrate the appearance of a single bright fragment and a "garbage cloud" distributed along the original ring.At larger values of f 0 , our simulations produce examples of the "clean" fragmentation, viz., with N = 4 produced by the unstable vortex rings with S = 1 in Fig. 8(a), N = 5 produced by S = 2 and 3 in (b) and (d), N = 7 by S = 2 and 3 in (c) and (e), and N = 8 by S = 3 in (f).These outcomes of the instability development are observed at the same evolution time t = 100 as in Figs.3(b), 5(b), and 6(b).The gradual increase of the number of the fragments on S is explained by the dependence of the azimuthal index of the fastest growing eigenmode of the breaking instability on the underlying winding number S, which is a generic property of vortex solitons[2,40].

FIG. 7 .
FIG. 7. (a)The power versus f 0 for the confined vortex modes in the self-focusing case (σ = +1).The VA solutions for S = 1 and 2 are shown by solid blue and dashed black lines, respectively.The corresponding numerical solutions are represented by red circles and green squares, respectively.Recall that the numerical solutions are stable, in this case, at f 0 < 1.6 and f 0 < 1.1, for S = 1 and 2, respectively.(b) The VA-predicted and numerically obtained (the dashed red and solid blue lines, respectively) profiles of the stable solution with S = 2 and f 0 = 0.4, drawn as cross sections through y = 0 .The other parameters are g = 0, σ = +1, α = 1, W = 2.

FIG. 9 .
FIG. 9. (a) The power versus f 0 for the vortex modes in the self-defocusing case (σ = −1).The VA solutions for S = 1 and 2 are shown by solid blue and dashed black lines, respectively.The corresponding numerical solutions are presented by red circles and green squares, respectively.(b) The VA-predicted and numerically obtained (the dashed red and solid blue lines, respectively) profiles of the solution with S = 2 and f 0 = 2, drawn as the cross sections through y = 0 .The other parameters are g = 0, σ = −1, α = 1, W = 2.