Stable Patterns in the Lugiato-Lefever Equation with a Confined Vortex Pump

1. 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, the LL equations play a crucial role in understanding fundamental phenomena such as the modulation instability (MI) and pattern formation in dissipative environments [8,10,13,15,18,19,20,23,25,27]. Its relevance 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,13,14,15,16,17,18,19,20,21,21,22,23,24], as well 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,28,29,30,31,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, were reported in Ref. [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 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\left(x,y,t\right)$, with the cubic or cubic-quintic nonlinearity:

$$\frac{\partial u}{\partial t}=-\alpha u+\frac{i}{2}{\nabla }^{2}u+i\sigma \left({\left|u\right|}^2-{\eta }^2\right)u-ig\left({\left|u\right|}^4-{\eta }^4\right)u+f\left(r\right)e^{iS\theta },$$

(1)

and a confined pump beam. Here $\alpha >0$ is the loss parameter, $\sigma =+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 $\eta$ defines the cavity mismatch, $\sigma {\eta }^2-g{\eta }^4$, in terms of the linearized LL equation, and

$$f\left(r\right)=f_0{r}^{S}exp\left(-r^2/W^2\right)$$

(2)

represents the confined pump beam with real amplitude $f_0$, radial width W, and integer vorticity $S\ge 1$, written in terms of polar coordinates $\left(r,\theta \right)$. Vortex beams, shaped by the passage of the usual laser beam through an appropriate phase mask, are available in the experiment [37].

Stationary solutions of Eq. (1) are characterized by values of the total power (alias norm),

$$P=\int \int {\left|u\left(x,y\right)\right|}^{2}dxdy\equiv 2\pi {\int }_0^{\infty }{\left|u(r,t)\right|}^{2}rdr,$$

(3)

and angular momentum,

$$M=i\int \int u^{*}\left(y{\partial }_{x}u-x{\partial }_{y}u\right)dxdy$$

(4)

(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\left(x,y\right)=u\left(r\right)e^{iS\theta }$, the expressions for the power and angular momentum are simplified:

$$P=2\pi {\int }_0^{\infty }{\left|u\left(r\right)\right|}^{2}rdr,M=SP.$$

(5)

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). It is a challenging problem, as it is well known from the work with models based on the nonlinear Schrodinger and CGL equations that, in the absence of the tight confining potential, vortex-ring solitons are vulnerable to the splitting instability, which may be considered as MI of the ring against azimuthal perturbations which break its axial symmetry [2,38]. 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\le S\le 4$ in the space of parameters of Eqs. (1) and (2) (in particular, in the plane of $\left(f_0,\alpha \right)$) for both signs of the cubic nonlinearity, $\sigma =\pm 1$, while the mismatch parameter is fixed to be $\eta =1$ by dint of scaling. The stability areas are vast, provided that the loss coefficient $\alpha$ 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\ne 0$, is considered too. Quite surprisingly, a stability area for the vortices with $S\le 3$ is found even in the case of $\sigma =+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 [38].

The rest of the paper is structured as follows. The analytical approach, based on the appropriate VA, is presented in Section 2. An asymptotic expression for the tail of the vortex solitons, decaying at $r\to \infty$, 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 Section 3. The paper is concluded by Section 4.

2. Analytical Considerations

2.1. 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\to \infty$:

$$u(r,\theta )\approx (i/2)W^4{r}^{S-2}exp\left(-r^2/W^2+iS\theta \right),$$

(6)

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 local maximum of local power ${\left|u\left(r\right)\right|}^2$ at $S>2$, at

$$r^2=r_{max}^2\equiv \sqrt{(S/2-1)}W.$$

(7)

A local maximum is indeed observed in numerically found radial profiles of all vortex solitons (see Figure 5, Figure 9, and Figure 10 below, for $S=5$, $S=5$, and $S=3$, respectively). In fact, for these cases Eq. (7) predicts values of $r_{max}$ which are smaller by a factor $\simeq 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.,

$$u(r,\theta )\approx \mathrm{u}\left(r\right)r^{S}exp\left(-r^2/W^2+iS\theta \right),$$

(8)

where $\mathrm{u}\left(r\right)$ is a complex slowly varying function, in comparison with those which explicitly appear 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

$$\mathrm{u}\left(r\right)=f_0{\left[\alpha -i\left(\frac{2r^2}{W^4}-\frac{2(S+1)}{W^2}-\sigma +g\right)\right]}^{-1},$$

(9)

where, as said above, $\eta =1$ is substituted. In the limit of $r\to \infty$, 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

$$r^2={\left(r_{max}^2\right)}_{\mathrm{TF}}=\left(S+1\right)W^2+\frac{\sigma -g}{2}{W}^4,$$

(10)

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 ${\left(r_{max}^2\right)}_{\mathrm{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 Figure 5 and Figure 9, demonstrates that Eq. (10) predicts, relatively accurately, ${\left(r_{max}\right)}_{\mathrm{TF}}=4$ for $S=5$, $W=2$, $\sigma =-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\to 0$, the asymptotic form of the solution is simple, $u(r,\theta )\approx u_0{r}^S$, but constant $u_0$ cannot be found in an explicit form, as it depends on the global structure of the vortex-soliton solution. The crude TF approximation given by Eq. (9) yields $u_0=f_0{\left[\alpha +i\left(\frac{2(S+1)}{W^2}+\sigma -g\right)\right]}^{-1}$.

2.2. 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. [39], which absorbs the dissipative term:

$$u\left(r,\theta ,t\right)=U(r,t)e^{iS\theta -\alpha t},$$

(11)

producing the following time-dependent equation for complex function $U(r,t)$, where, as said above, we set $\eta =1$ by means of scaling:

$$\frac{\partial U}{\partial t}=\frac{i}{2}\left(\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{S^2}{r^2}\right)U+i\sigma \left({\left|U\right|}^2{e}^{-2\alpha t}-1\right)U-ig\left({\left|U\right|}^4{e}^{-4\alpha t}-1\right)U+f\left(r\right)e^{\alpha t}.$$

(12)

The real Lagrangian which precisely produces the time-dependent equation (12) is

$$\begin{array}{c}L={\int }_0^{\infty }\left[\frac{i}{2}\left(\frac{\partial U^{*}}{\partial t}U-U^{*}\frac{\partial U}{\partial t}\right)+\frac{1}{2}{\left|\frac{\partial U}{\partial r}\right|}^2+\left(\frac{S^2}{2r^2}+\sigma -g\right){\left|U\right|}^2\right.\\ -\left.\frac{\sigma }{2}{e}^{-2\alpha t}{\left|U\right|}^4+\frac{g}{3}{e}^{-4\alpha t}{\left|U\right|}^6+if\left(r\right)e^{\alpha t}\left(U^{*}-U\right)\right]rdr.\end{array}$$

(13)

The simplest ansatz which may be used as the basis for VA follows the form of the pump term (2):

$$U(r,t)=U_0{e}^{i\varphi }{r}^{S}exp\left(-\frac{r^2}{W^2}+\alpha t\right),$$

(14)

where variational parameters $U_0$ and $\varphi$ are the real amplitude and phase shift of the solution with respect to the pump. Power (3) for this ansatz is

$$P_S=\pi \Gamma (S+1){\left(\frac{W^2}{2}\right)}^{S+1}{U}_0^2,$$

(15)

where $\Gamma (S+1)\equiv S!$ is the Gamma-function. Note that the local power ${\left|U\left(r\right)\right|}^2$, corresponding to ansatz (14), attains it maximum at $r^2=S{W}^2/2$.

The substitution of ansatz (14) in Lagrangian (13) and straightforward integration yields the respective VA Lagrangian,

$$\begin{array}{ccc}\hfill L_{\mathrm{VA}}& =& U_0{e}^{2\alpha t}\left\{\left[6^{-(3S+2)}\Gamma \left(3S+1\right)g{W}^{2(3S+1)}\right]U_0^5-\left[2^{-4(S+1)}\Gamma \left(2S+1\right)\sigma W^{2(2S+1)}\right]U_0^3\right.\hfill \\ & +& 2^{-(S+2)}\Gamma \left(S+1\right)W^{2S}\left[\left(\frac{d\varphi }{dt}+(\sigma -g)\right)W^2+\left(S+1\right)\right]U_0\hfill \\ & +& \left.2^{-(S+1)}\Gamma \left(S+1\right)W^{2(S+1)}{f}_{0}sin\varphi \right\}.\hfill \end{array}$$

(16)

Then, the Euler-Lagrange equations for $U_0$ and $\varphi$ are obtained as

$$\partial L_{\mathrm{VA}}/\partial U_0=\delta L_{\mathrm{VA}}/\delta \varphi =0,$$

(17)

where $\delta /\delta \varphi$ stands for the variational derivative. Taking into regard that Lagrangian (16) must be substituted into the respective action, $\int L_{\mathrm{VA}}dt$, and then the action must be actually subjected to the variation, one should apply the time differentiation to factor $e^{2\alpha t}$ in Lagrangian (16), while deducing the appropriate form of $\delta L_{\mathrm{VA}}/\delta \varphi$ in Eq. (17). Once the Euler-Lagrange equations (17) have been derived , we consider their stationary (fixed-point) solutions by setting $d{U}_0/dt=d\varphi /dt=0$, which yields

$$U_0\alpha -f_{0}cos\left(\varphi \right)=0,$$

(18)

$$\begin{array}{ccc}\hfill 6^{-(3S+1)}\Gamma \left(3S+1\right)g{W}^{2(2S+1)}{U}_0^5-2^{-2(2S+1)}\Gamma \left(2S+1\right)\sigma W^{2S+2}{U}_0^3& & +\hfill \\ \hfill +2^{-(S+1)}\Gamma \left(S+1\right)\left\{f_0{W}^{2}sin\varphi +\left[(\sigma -g)W^2+\left(S+1\right)\right]U_0\right\}& =& 0.\hfill \end{array}$$

(19)

It is relevant to mention that the evolution equation for the power (3), which follows from Eq. (12), is

$$\frac{dP}{dt}=-2\alpha P+4\pi {\int }_0^{\infty }f\left(r\right)\mathrm{Re}\left\{U(r,t)\right\}{e}^{-\alpha t}rdr.$$

(20)

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,

$$cos\varphi =\frac{\alpha U_0}{f_0},$$

(21)

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 $\varphi =0$ in Eq. (21):

$$U_0\le {\left(U_0\right)}_{max}=f_0/\alpha .$$

(22)

2.3. VA for the Cubic ($g=0$) and Quintic ($g\to \infty$) 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 $\sim f_0$ in Eq. (19), one obtains a relatively simple expression which predicts the squared amplitude of the vortex soliton,

$${\left(U_0^2\right)}_{\mathrm{VA}}^{(g=0)}=\frac{2^{2S+1}}{(2S-1)!!}{W}^{-2S}\left[1+\sigma \left(S+1\right)W^{-2}\right],$$

(23)

the respective power (15) of the underlying ansatz (14) being

$$P_{\mathrm{VA}}^{(g=0)}=\frac{2^{S}S!\pi }{(2S-1)!!}\left[W^2+\sigma \left(S+1\right)\right].$$

(24)

Note that expressions (23) and (24) are always meaningful for the self-focusing sign of the cubic nonlinearity, $\sigma =+1$, while in the case of defocusing, $\sigma =-1$, the expressions are meaningful if they are positive, which imposes a restriction on the fixed width of the Gaussian pump: it must be broad enough, viz.,

$$W^2>S+1.$$

(25)

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 Figure 1(a) and (b) for the self-focusing and defocusing signs of the cubic term, $\sigma =+1$ and $-1$, respectively. Note that in Figure 1(b) for $\sigma =-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\to \pm \infty$, opposite to the pure cubic model considered above, the asymptotic solution of Eq. (19) is

$${\left(U_0^2\right)}_{\mathrm{VA}}^{(g\to \pm \infty )}\approx 2^S\sqrt{3^{3S+1}\frac{S!}{\left(3S\right)!}}{W}^{-2S},$$

(26)

the respective expression for the power (15) being

$$P_{\mathrm{VA}}^{(g\to \pm \infty )}=\frac{3^{(3S+1)/2}\pi {\left(S!\right)}^{3/2}}{2\sqrt{\left(3S\right)!}}{W}^2,$$

(27)

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.

3. Numerical Results

Simulations of Eq. (1) were conducted by means of the split-step pseudo-spectral algorithm. The solution started from the zero input, and was running until convergence to an apparently stable stationary profile. 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 $\alpha$, $\sigma =\pm 1$, g, and $f_0$, W, S (see Eq. (2)). The results are presented below by varying, severally, loss $\alpha$, 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 Figure 11.

3.1. Variation of the Loss Parameter $\alpha$

In Figure 2(a) we display the cross-section (drawn through $y=0$) of the variational and numerical solutions for the stable vortex solitons obtained with $\alpha =0.5$, $1.0$, and $2.0$, while the other parameter are fixed as $\sigma =+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 Figure 2(a) is characterized by the relative power difference from their numerically found counterparts, which is $5.1\%$, $3.1\%$, and $0.5\%$ for

$$\alpha =0.5,1.0,2.0,$$

(28)

respectively. Thus, the VA accuracy improves with the increase of $\alpha$.

Similar results for the self-repulsive nonlinearity, $\sigma =-1$, are presented in Figure 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 Figure 2(a).

There is a critical value of $\alpha$ below which the vortex solitons are unstable. As an example, Figure 3 shows the VA-predicted and numerically produced solutions for $\alpha =0.2$ and $\sigma =+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 Figure 3, the instability boundary is ${\alpha }_{\mathrm{crit}}\approx 0.35$. The stabilizing effect of the loss at $\alpha >{\alpha }_{\mathrm{crit}}$ is a natural feature. On the other hand, the increase of $\alpha$ leads to decrease of the soliton’s amplitude, as seen in Figure 2.

In the case of the self-defocusing ($\sigma =-1$), all the numerically found vortex modes are stable, at least, at $\alpha \ge 0.1$, although the discrepancy in the values of the power between these solutions and their VA counterparts is very large at small $\alpha$, exceeding $75\%$ at $\alpha =0.1$. At still smaller values of $\alpha$, the relaxation of the evolving numerical solution toward the stationary state is very slow, which makes it difficult to identify the stability.

3.2. 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 $\alpha =1$, $f_0=1$, $W=2$ in Eqs. (1) and (2). In the self-focusing case ($\sigma =+1$), the numerically produced solutions are stable for $S=1$ and 2, and unstable for $S\ge 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 an relatively accurate approximation in this case.

In the self-defocusing case ($\sigma =-1$), considering the same values of the other parameters as used above, the numerical solution produces stable vortex solitons at least until $S=5$. As above, the accuracy of VA is much lower for $\sigma =-1$ than for the self-focusing case ($\sigma =+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.

3.3. 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$, $\alpha =1$, $f_0=1$, and $S=1$. In Figure 4(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., $\sigma =+1$ and $\sigma =-1$, respectively. In the former case, the azimuthal MI sets in at $W\ge 2.7$, see an example in Figure 4(b) for $W=2.75$.

In the self-focusing case, $\sigma =+1$, the azimuthal MI for the solitons with higher vorticities, $S=2$, 3, 4, or 5, sets in at $W\ge 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.

3.4. Variation of the Pump’s Strength $f_0$

Effects of the variation of $f_0$ are reported here, fixing other parameters as $g=0$, $\alpha =1$, and $W=2$. In the case of the cubic self-focusing, $\sigma =+1$, the vortex soliton with $S=1$ are subject to MI at $f_0\ge 1.6$. As a typical example, in Figure 5 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\ge 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 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 Figure 6(a). As an example, Figure 6(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\le 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.

Figure 5. (a) The VA-predicted profile of ${\left|u\right|}^2$ in the self-focusing case, with parameters $g=0$, $\sigma =+1$, $\alpha =1$, $f_0=1.7$, $W=2$, and $S=1$. (b) The unstable solution, produced, at $t=100$, by the simulations of Eq. (1) for the same parameters.

Figure 5. (a) The VA-predicted profile of ${\left|u\right|}^2$ in the self-focusing case, with parameters $g=0$, $\sigma =+1$, $\alpha =1$, $f_0=1.7$, $W=2$, and $S=1$. (b) The unstable solution, produced, at $t=100$, by the simulations of Eq. (1) for the same parameters.

Figure 6. (a) The power versus $f_0$ for the confined vortex modes in the self-focusing case ($\sigma =+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$, $\sigma =+1$, $\alpha =1$, $W=2$.

Figure 6. (a) The power versus $f_0$ for the confined vortex modes in the self-focusing case ($\sigma =+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$, $\sigma =+1$, $\alpha =1$, $W=2$.

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\ge S$ of symmetrically placed localized fragments [2,38], while the above examples, displayed in Figure 3(b), Figure 4(b), and Figure 5(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 Figure 7(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 Figure 3(b), Figure 4(b), and Figure 5(b).

In self-defocusing case, $\sigma =-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 Figure 8(a) (cf. Figure 6(a) for $\sigma =+1$). Naturally, the VA-numerical discrepancy increases with the growth of the pump’s strength, $f_0$, see an example in Figure 8(b). Unlike the case of $\sigma =+1$, in the case of the self-defocusing the vortex modes with $S\le 5$ remain stable, at least, up to $f_0=5$ (here, we do not dot consider the case of $S>5$).

3.5. 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 Figure 3, which demonstrated that the vortex soliton with $S=1$, as a solution to Eqs. (1) and (2) with $g=0$, $\sigma =+1$ and $f_0=1$, $W=2$, was unstable if the loss parameter fell below the critical value, $\alpha =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 Figure 3, which was unstable in the absence of the quintic term. This conclusion seems surprising because, in most cases, the inclusion of the self-focusing quintic nonlinearity gives rise to the supercritical collapse in 2D, making all solitons strongly unstable [2,38]. However, as it can be concluded from the stability charts displayed below in Figure 11, the stabilizing effect of the quintic self-focusing occurs only at moderately small powers. In the general case, solitons stability regions shrinks 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, with parameters produced by a numerical solution of Eqs. (18) and (19), is presented in Figure 9. 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 Figure 9, the power-measured discrepancy for $g=-1$ and $+1$ is, respectively, $4\%$ and $64\%$.

Figure 9. Cross-section profiles (drawn through $y=0$) for stable vortex solitns with $S=1$, produced by Eq. (1) with the quintic self-focusing, $g=-1$ (a) or defocusing, $g=+1$ (b) term. In both cases, the cubic self-focusing cubic term, with $\sigma =+1$, is present. The numerically found profiles and their VA-produced counterparts are displayed, respectively, by the solid blue and dashed lines. Other parameter are $\alpha =0.2$, $\eta =1$, and $f_0=1$, $W=2$.

Figure 9. Cross-section profiles (drawn through $y=0$) for stable vortex solitns with $S=1$, produced by Eq. (1) with the quintic self-focusing, $g=-1$ (a) or defocusing, $g=+1$ (b) term. In both cases, the cubic self-focusing cubic term, with $\sigma =+1$, is present. The numerically found profiles and their VA-produced counterparts are displayed, respectively, by the solid blue and dashed lines. Other parameter are $\alpha =0.2$, $\eta =1$, and $f_0=1$, $W=2$.

Another noteworthy finding is the stabilization of higher-vorticity solitons by the quintic term. For instance, it was shown above that, for parameters $\alpha =1$, $\eta =1$, $f_0=1$, $W=2$, and $\sigma =+1$ (the cubic self-focusing) all vortex solitons with $S\ge 3$, produced by Eq. (1), were unstable. Now, we can demonstrate that the soliton with $S=3$ can be stabilized by adding a self-defocusing quintic with a small coefficient, just $g=0.1$, see Figure 10(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 Figure 10(a) (recall that $g=-1$ is sufficient for the stabilization of the vortex soliton with $S=1$ and $\alpha =0.2$ in Figure 9(a)). For the stabilized vortex modes shown in Figure 10(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, $\approx 0.7\%$, while in the presence of the weak quintic self-defocusing in Figure 10(b) the discrepancy is $6.6\%$.

Figure 10. Cross-section profiles (drawn through $y=0$) of vortex solitons with $S=3$, stabilized by the self-focusing (a) or defocusing (b) quintic tterm in Eq. (1), with the respective coefficient $g=-6$ or $g=0.1$, other parameters being $\sigma =-1$ (cubic self-focusing), $\alpha =1$, $\eta =1$, and $f_0=1$, $W=2$ in Eq. ( 2). The numerically found solutions and their VA-predicted counterparts are presented by the solid blue and dashed red lines, respectively.

Figure 10. Cross-section profiles (drawn through $y=0$) of vortex solitons with $S=3$, stabilized by the self-focusing (a) or defocusing (b) quintic tterm in Eq. (1), with the respective coefficient $g=-6$ or $g=0.1$, other parameters being $\sigma =-1$ (cubic self-focusing), $\alpha =1$, $\eta =1$, and $f_0=1$, $W=2$ in Eq. ( 2). The numerically found solutions and their VA-predicted counterparts are presented by the solid blue and dashed red lines, respectively.

3.6. Stability Charts in the Parameter Space

The numerical results produced in this work are summarized in the form of stability areas plotted in Figure 11 in the parameter plane $\left(f_0,\alpha \right)$, 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, $\sigma =+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 ($\sigma =-1$) and quintic focusing ($g=-1$), also for $S=1,2,3$, are plotted in Figure 12.

The choice of the parameter plane $\left(f_0,\alpha \right)$ is relevant, as the strength of the pump beam, $f_0$, and loss the coefficient, $\alpha$, and are amenable to accurate adjustment in the experiment (in particular, $\alpha$ 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 Figure 11 and Figure 12, the increase of $\alpha$ naturally provides effective stabilization of the vortex modes, while none of them may be stable at $\alpha =0$, in agreement with the known properties of vortex-soliton solutions to the 2D nonlinear Schrödinger equation with the cubic and/or cubic-quintic nonlinearity [2,38]. 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 Figure 11 are the general stabilizing/destabilizing effect of the quintic self-defocusing/focusing, and expansion of the splitting-instability area with the increase of S (the latter feature is also exhibited by Figure 12).

Figure 11. Stability areas for families of the vortex solitons with winding numbers $S=1$, $2,$ and 3, in the plain of the loss coefficient ($\alpha$) 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 $\sigma =+1$ (the self-focusing cubic term), $\eta =1$, and $W=2$.

Figure 11. Stability areas for families of the vortex solitons with winding numbers $S=1$, $2,$ and 3, in the plain of the loss coefficient ($\alpha$) 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 $\sigma =+1$ (the self-focusing cubic term), $\eta =1$, and $W=2$.

Figure 12. The same as in Figure 11, but for $\sigma =-1$ and $g=-1$ (the cubic self-defocusing and quintic focusing terms).

Figure 12. The same as in Figure 11, but for $\sigma =-1$ and $g=-1$ (the cubic self-defocusing and quintic focusing terms).

The instability mode which determines the boundary of the stability areas in Figure 11 and Figure 12 is the breaking of the axial symmetry of the vortex rings by azimuthal perturbations, as shown above in Figure 3(b), Figure 4(b), and Figure 5(b). The destabilization through the spontaneous splitting of the rings into symmetric sets of fragments (see Figure 7) occurs deeply inside the instability area, at larger values of $f_0$.

4. Conclusions

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\le 5$. Stable states in the form of vortex solitons (rings) with the same 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 onset of the azimuthal instability breaking their axial symmetry. These findings suggest new possibilities for the design of tightly confined 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.