Zero-dimensional limit of the two-dimensional Lugiato-Lefever equation

TL;DR

This paper investigates the zero-dimensional limit of the 2D Lugiato-Lefever equation, revealing stable localized optical modes, including vortex solitons, through variational, numerical, and analytical methods.

Contribution

It introduces a detailed analysis of tightly confined optical modes in the 2D LL model, including vortex states, using variational, numerical, and Thomas-Fermi approximations.

Findings

Stable pixel-like localized modes identified.

Vortex solitons with embedded vorticity are stable over broad parameter ranges.

Multi-ring vortex states with spiral phase are also stable.

Abstract

We study effects of tight harmonic-oscillator confinement on the electromagnetic field in a laser cavity by solving the two-dimensional Lugiato-Lefever (2D LL) equation, taking into account self- focusing or defocusing nonlinearity, losses, pump, and the trapping potential. Tightly confined (quasi-zero-dimensional) optical modes (pixels), produced by this model, are analyzed by means of the variational approximation, which provides a qualitative picture of the ensuing phenomena. This is followed by systematic simulations of the time-dependent 2D LL equation, which reveal the shape, stability, and dynamical behavior of the resulting localized patterns. In this way, we produce stability diagrams for the expected pixels. Then, we consider the LL model with the vortical pump, showing that it can produce stable pixels with embedded vorticity (vortex solitons) in remarkably broad sta- bility areas. Alongside confined vortices with the simple single-ring structure, in the latter case the LL model gives rise to stable multi-ring states, with a spiral phase field. In addition to the numeri- cal results, a qualitatively correct description of the vortex solitons is provided by the Thomas-Fermi approximation.

Equations

$$i\left(\gamma +\frac{\partial }{\partial t}\right)\varphi =[-\frac{1}{2}\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)+\Delta$$

$$+\frac{{\Omega }^2}{2}(x^2+y^2)+\sigma |\varphi {|}^2]\varphi +E,$$

$$\varphi (r)\approx -\frac{2E}{{\left(\Omega r\right)}^2}+\frac{4\left(\Delta -i\gamma \right)E}{{\left(\Omega r\right)}^4}.$$

$$\frac{dP}{dt}=-2\gamma P-2\int \int \mathrm{Im}\{E^{\ast }\varphi (x,y,t)\}dxdy,$$

$$P=\int \int |\varphi (x,y,t){|}^{2}dxdy.$$

$$\varphi (x,y,t)\equiv \Phi (x,y,t)\exp \left(-\gamma t\right),$$

$$i\frac{\partial }{\partial t}\Phi$$

$$L$$

$$\Phi \left(x,y,t\right)$$

$$ϵ$$

$$\frac{e^{-2\gamma t}}{\pi {ϵ}^2}{L}_{\mathrm{eff}}$$

$$\frac{1}{2}\left[q_2(\chi )\cos \chi +q_2^{\prime }(\chi )\sin \chi +q_1(\chi )\right]\frac{df}{dt}$$

$$\Delta q_1(\chi )+2f{q}_2(\chi )+\frac{3g}{8}{f}^2{q}_3(\chi )+\gamma q_2(\chi )\sin \chi$$

$$g=\sigma {ϵ}^2\equiv 4\sigma E^2/{\Omega }^4.$$

$$\left({\Omega }^2\sin \chi -\gamma f\right)q_1(\chi )+\left(\Delta f-{\Omega }^2\cos \chi \right)q_1^{\prime }(\chi )$$

$$\Delta q_1(\chi )+2f{q}_2(\chi )+\frac{3g}{8}{f}^2{q}_3(\chi )+\gamma q_2(\chi )\sin \chi =0.$$

$$|\chi |\ll \pi .$$

$$\Delta +\frac{2}{3}f+\frac{g}{2}{f}^2+\frac{\gamma }{3}\chi =0.$$

$$\chi =\frac{30\gamma f}{10{\Omega }^2+f\left[10\Delta +8f+3g{f}^2\right]},$$

$$f$$

$$\chi$$

$$f$$

$$\chi$$

$$P$$

$$\frac{dP}{dt}$$

$$\sin (\chi )=\frac{2\gamma }{{\Omega }^2}f.$$

$$f\le {\Omega }^2/(2\gamma ).$$

$$\frac{g{\Omega }^4{\chi }^3}{80{\gamma }^3}+\frac{{\Omega }^2{\chi }^2}{15{\gamma }^2}+\frac{\Delta \chi }{6\gamma }-\frac{1}{6}$$

$$\frac{g{\Omega }^4{\chi }^2}{8{\gamma }^2}+\frac{1}{3\gamma }\left({\Omega }^2+{\gamma }^2\right)\chi +\Delta$$

$$\left(\Delta -i\gamma +\frac{{\Omega }^2}{2}{r}^2+\sigma |\varphi {|}^2\right)\varphi =-E.$$

$$\varphi (r)\approx -E{\left(\Delta +\frac{{\Omega }^2}{2}{r}^2\right)}^{-1}.$$

$$\varrho (x)\equiv {\int }_{-\infty }^{+\infty }|\varphi (x,y){|}^{2}dy,$$

$${\varrho }_{\mathrm{TFA}}(x)=\frac{2\pi E^2}{{\Omega }^4}{\left(\frac{2\Delta }{{\Omega }^2}+x^2\right)}^{-3/2}.$$

Full text

I Introduction

The Lugiato-Lefever (LL) equation LL in one and two dimensions (1D and 2D) is a fundamental model governing the dynamics of optical fields in pumped lossy laser cavities with the intrinsic Kerr nonlinearity, which may have self-focusing or defocusing sign. This equation is well known as an important tool for the analysis of pattern formation, with various applications in nonlinear optics patterns-early ; pixel . The progress in theoretical and experimental studies has recently drawn a great deal of renewed interest to the use of the LL equation in diverse settings patterns-2 -NJP . A natural extension of these studies is incorporation of external potentials into the LL equation, which can be easily fabricated in laser cavities as transverse landscapes of the refractive-index inhomogeneity, and may be used as an efficient means for the control of optical fields Kestas .

One of essential applications of the LL equation is its use for modeling well-localized pixels (i.e., sharply bounded bright spots) in the cavity pixel . In most cases, pixels are considered as anti-dark solitons, i.e., bright objects created on top of a uniformly pumped background field. In this work, we aim to demonstrate a possibility to create completely localized robust pixels (i.e., bright solitons with zero background), by adding to the model a confining potential corresponding to an isotropic 2D harmonic oscillator. Furthermore, we demonstrate that the same setting makes it possible to create stable vortex pixels, by applying a vortically structured pump. The consideration reported below combines an analytical approach, chiefly based on the variational and Thomas-Fermi approximations (VA and TFA), and systematic direct simulations, in imaginary and real time alike, with the purpose to create confined modes and test their stability.

The paper is organized as follows. The model, based on the 2D LL equation with the harmonic-oscillator trapping potential, is formulated in Section II. Analytical treatment, which makes use of the VA, power-balance equation, and TFA, is presented in Section III. Numerical results for the existence and stability of the fundamental (zero-vorticity) and vortical trapped modes are reported in Sections IV and V, respectively. The latter section also reports simple analytical results for the vortex states, obtained by means of the TFA. The paper is concluded by Section VI.

II The model

The 2D LL equation for the amplitude $\phi(x,y,t)$ of the electromagnetic field in a pumped lossy laser cavity is (see, e.g., Ref. Kestas )

[TABLE]

where $E$ is the pump field, $\gamma>0$ the dissipation rate, $\Delta\gtrless 0$ detuning of the pump with respect to the cavity, and $\Omega^{2}$ the strength of the confining potential, while $\sigma=-1$ and $+1$ correspond to the self-focusing and defocusing nonlinearity, respectively. By means of rescaling, one may fix $\gamma=1$, although it may be convenient to keep $\gamma$ as a free parameter, as shown below.

Stationary solutions to Eq. (1) have a simple asymptotic form at $r\equiv\sqrt{x^{2}+y^{2}}\rightarrow\infty$:

[TABLE]

We also note that the following exact power-balance equation ensues from Eq. (1):

[TABLE]

where power $P$ (alias norm) of the solitary wave is defined as

[TABLE]

The objective is to reduce the 2D LL equation (1) to a quasi-zero-dimensional limit (a dynamical system for a pixel, similar to those realized by theoretically predicted pixel and experimentally created NJP spatial solitons) in the case of tight confinement, represented by large $\Omega^{2}$. First, we do it by means of the VA, defining

[TABLE]

and thus casting Eq. (1) in the form of

[TABLE]

Unlike the original LL equation (1), the transformed one (6) can be directly derived from a real time-dependent Lagrangian,

[TABLE]

III Analytical considerations

III.1 The variational approximation

For the 1D LL equation without trapping potentials, the VA was developed in Ref. VA-LL . To derive this approximation in a form appropriate for the present model, we note that, in the lowest approximation, Eq. (6) gives rise to the following asymptotic form of solutions at $r\rightarrow\infty$: $\Phi=-2E\left(\Omega r\right)^{-2}e^{\gamma t}$, cf. Eq. (2). This form suggests us to adopt an ansatz based on the fractional expression, with real variables $f(t)$ and $\chi(t)$, which may be combined into a complex one, $F(t)=f(t)\exp\left(i\chi(t)\right)$:

[TABLE]

The insertion of ansatz (8) in Eq. (7) and subsequent integration gives rise to an effective Lagrangian,

[TABLE]

with $q_{1}(\chi)\equiv\chi/\sin\chi$, $q_{2}(\chi)\equiv[\left(\sin\chi-\chi\cos\chi\right)/\sin^{3}\chi$, and $q_{3}(\chi)\equiv[2\chi-\sin\left(2\chi\right)]/\sin^{3}\chi$. The last term in Eq. (10) is cast in the integral form as a result of “renormalization”: the respective term in the original Lagrangian formally diverges logarithmically at $R\rightarrow\infty$, but the diverging part actually does not depend on $f$ and $\chi$, and it may be cancelled by means of the differentiation with respect to $\chi$ and subsequent integration, also with respect to $\chi$.

The Euler-Lagrange equations following from Lagrangian (10) are (taking into account that the Lagrangian must be substituted into the action, $\int Ldt$, and then the action must be subjected to the variation; this makes it necessary to apply the time differentiation to $e^{2\gamma t}$):

[TABLE]

[TABLE]

where a renormalized nonlinearity coefficient is [see Eq. (9)]

[TABLE]

Note that, although it may seem that Eqs. (11) and (12) are singular at $\chi=0$, in reality all the singularities cancel. A singularity is instead possible at $\chi=\pi$.

We consider stationary (fixed-point) solutions of Eqs. (11) and (12) by setting $df/dt=d\chi/dt=0$, which yields

[TABLE]

[TABLE]

Further, it is possible to find approximate solutions of Eqs. (14) and (15), assuming that they have

[TABLE]

In this case, Eq. (15), in the first approximation, assumes the form of

[TABLE]

Similarly, in the lowest approximation Eq. (14) yields an expression for $\chi$:

[TABLE]

The assumption (16) may be then secured by a natural assumption of the strong confinement, i.e., considering large values of $\Omega$. In this case, Eqs. (17) and (18) can be further simplified to

[TABLE]

Obviously, Eqs. (19) and (20) produce a physically relevant result under condition $g\Delta<2/9$.

One can construct another approximate solution for large detuning $\Delta$:

[TABLE]

In the general case, stationary solutions of Eqs. (14) and (15), where, as said above, we may fix $\gamma=1$, depend on three parameters: $\Delta\gtrless 0$, $g\gtrless 0$ [see Eq. (13)], and $\Omega^{2}>0$.

In addition to the consideration of the stationary solutions (fixed points), the full dynamical version of the VA, based on Eqs. (11) and (12), can be also used to analyze their stability, as well as evolution of unstable solutions. However, in practical terms such a dynamical analysis turns out to be quite cumbersome, direct numerical simulations being actually more efficient, as shown below.

III.2 The power-balance condition

The substitution of ansatz (8) in the definition of power (4) and power-balance equation (3) yields

[TABLE]

(in these expressions, $f>0$ is implied). Equation (24) predicts the equilibrium condition, $dP/dt=0$, at

[TABLE]

Note that $E$ drops out from Eq. (25), and condition $\sin(\chi)\leq 1$, following from Eq. (25), imposes a restriction on $f$,

[TABLE]

Finally, for $|\chi|\ll\pi$ Eq. (25) simplifies to $\chi\approx\left(2\gamma/\Omega^{2}\right)f$. Using this to eliminate $f$ in favor of $\chi$, Eqs. (17) and (18) give rise to the following system of equations:

[TABLE]

Of course, the system of two equations (27) and (28) for the single unknown $\chi$ is overdetermined, and a solution of this system may exist only if a special restriction is imposed on parameters, as shown in Fig. 1, in the plane of ($\Delta$, $g$), for $\gamma=1$ and three different fixed values of the confinement strength, $\Omega=2$, $\Omega=4$, and $\Omega=10$. Note that these curves do not depend on the pumping strength, $E$. Indeed, this parameter is related only to the power of the solution, see Eq. (23) and Fig. 2. The meaning of the overdetermined system is that, realizing the VA and power-balance condition simultaneously, its solution has a better chance to produce an accurate approximation. This expectation is qualitatively corroborated below, see Fig. 9 and related text in the next section.

Generic properties of the modes predicted by ansatz (8) are characterized by the corresponding dependence of power $P$ on the pumping strength, $E$. Using, for this purpose, the simplified approximation given by Eqs. (19) and (20), we display the dependences for the defocusing nonlinearity ($g=1$) in Fig. 2(a), at fixed values of the detuning, $\Delta=-1$, $-4$, and $-10$. Figure 2(b) displays the same dependences in the case of the self-focusing nonlinearity ($g=-1$), for $\Delta=1$, $4$, and $10$. Note that power $P$ is not symmetric with respect to the reversal of the signs of nonlinearity $g$ and detuning $\Delta$.

In Fig. 2(c) we compare the VA results for the self-defocusing ($g=1$ and $\Delta=-10$) and focusing ($g=-1$ and $\Delta=10$) cases. In addition, Fig. 2(c) includes full numerical results (for details see the next Section). It is seen that the simplified VA produces a qualitatively correct prediction, which is not quite accurate quantitatively. Below, we demonstrate that the VA is completely accurate only in small black regions shown in Fig. 9.

III.3 The Thomas-Fermi approximation (TFA)

In the case of the self-defocusing sign of the nonlinearity, and positive mismatch, $\Delta>0$, the ground state, corresponding to a sufficiently smooth stationary solution of Eq. (1), $\phi=\phi(r)$, may be produced by the TFA, which neglects derivatives in the radial equation TFA :

[TABLE]

In particular, the TFA is relevant if $\Delta$ is large enough.

The TFA-produced equation (29) for the ground-states’s configuration is not easy to solve analytically, as it is a cubic algebraic equation with complex coefficients. The situation greatly simplifies in the limit case of a very large mismatch, viz., $\Delta\gg\gamma$ and $\Delta^{3}\gg\sigma E^{2}$. Then, both the imaginary and nonlinear terms may be neglected in Eq. (29), to yield

[TABLE]

This simple approximation, which may be considered as a limit form of ansatz (8), can be used to produce estimates for various characteristics of the ground state (see, in particular, Fig. 5 below). In fact, the TFA will be the most relevant tool in Section V, as an analytical approximation for trapped vortex modes, for which the use of the VA, even in its stationary form, is too cumbersome.

The TFA cannot be applied to nonstationary solutions, hence it does not provide direct predictions for stability of stationary modes. However, it usually tends to produce ground states, thus predicting stable solutions. This expectation is corroborated by results produced below.

IV Numerical results for fundamental modes

IV.1 Stationary trapped modes

To obtain accurate results, and verify the validity of the VA predictions which are presented in the previous section, we here report results obtained as numerical solutions of Eq. (1). First, we aim to find ground-state localized states by means of imaginary-time propagation. In the framework of this method, one numerically integrates Eq. (1), replacing $t$ by $-it$ and normalizing the solution at each step of the time integration to maintain a fixed total power IT . For testing stability of stationary states, Eq. (1) was then simulated in real time, by means of the fourth-order split-step method implemented in a *GNU Octave * program Eaton_Octave (for a details concerning the method and its implementations in MATLAB, see Ref. Yang_10 ).

In Fig. 3 we show 1D integrated intensity profiles $\varrho(x)$, defined as

[TABLE]

and obtained from the imaginary-time solution of Eq. (1), along with their analytical counterparts produced by the VA based on Eq. (8), for three different values of detuning $\Delta$, viz., (a) $\Delta=-1$, (b) $\Delta=-4$, and (c) $\Delta=-10$, for $g=1$ (the self-defocusing nonlinearity) and $\Omega=10$ and $E=10$. We used Eqs. (19) and (20) to produce values of the parameters $f$ and $\chi$ in ansatz (8), which was then used as the initial guess in direct numerical simulations.

In Fig. 4 we display results similar to those shown in Fig. 3, but for the self-focusing nonlinearity ($g=-1$) and three different (positive) values of $\Delta$, viz., (a) $\Delta=1$, (b) $\Delta=4$, and (c) $\Delta=10$, with fixed $\Omega=10$ and $E=10$. In both cases of $g=\pm 1$, the VA profiles show good match to the numerical ones, although the accuracy slightly deteriorates with the increase of $|\Delta|$.

Note that the results displayed in Fig. 4, for the situations to which the TFA does not apply, because the nonlinearity is self-focusing in this case, demonstrate  the growth of the maximum value, $\varrho(x=0)$, with the increase of mismatch $\Delta$. In the case of self-defocusing it is natural to expect decay of $\varrho_{\mathrm{TFA}}(x=0)$ with the increase of $\Delta$. As shown in Fig. 5, this expectation is confirmed by the numerical results and the TFA alike. In particular, for the integrated intensity profile defined by Eq. (31), the simplest version of the TFA, produced by Eq. (30), easily gives

[TABLE]

Figure 5 also corroborates that the TFA, even in its simplest form, becomes quite accurate for sufficiently large values of $\Delta>0$.

IV.2 Stability of the stationary modes

The stability of the trapped configurations predicted by ansatz (8) was tested in real-time simulations of Eq. (1), adding $5\%$ random noise to the input. We display the results, showing the evolution of the solution’s norm (total power) in the case of the defocusing nonlinearity ($g=1$), for $\Delta=-1$ and $-10$, in Figs. 6(a) and (b), respectively. The insets show the asymptotic behavior at large times. In the case of the defocusing nonlinearity, the solution quickly relaxes to a numerically exact stationary form, and remains completely stable at $t>10$ (in fact, real-time simulations always quickly converge to stable solutions at all values of the parameters). However, in the case of the self-focusing with $\Omega=10$, the solutions are unstable, suffering rapid fragmentation, as seen in Fig. 8. This behavior is also exemplified in results shown in Figs. 6(c) and 6(d) for the temporal evolution of the solution’s total power in the case of the self-focusing nonlinearity ($g=-1$), for $\Delta=1$ and for $\Delta=10$, respectively. The instability of the fundamental modes in this case is a natural manifestation of the modulational instability in the LL equation MI . Note that the large size of local amplitudes in small spots, which is attained in the course of the development of the instability observed in Fig. 8, implies the trend to the onset of the 2D collapse driven by the self-focusing cubic nonlinearity collapse .

In Figs. 7 and 8 we display the time evolution of density profiles $|\phi|^{2}$ produced by the simulations of Eq. (1) with the self-defocusing and focusing nonlinearity, respectively. The input profiles are again taken as per the VA ansatz (8) with the addition of $5\%$ random noise. In Figs. 7(a) and 7(c) we show the perturbed input profiles in the case of self-defocusing, for $\Delta=-1$ and $\Delta=-10$, respectively, while the corresponding profiles at $t=1000$ are displayed in Figs. 7(b) and 7(d). Note that the agreement between the variational and numerical profiles tends to deteriorate with the increase of $|\Delta|$ [the same trend as observed in Fig. 2(c)].

Further, in Figs. 8(a) and 8(c) we display the perturbed input profiles in the case of the self-focusing nonlinearity for $\Delta=1$ and $\Delta=10$, respectively, with the corresponding profiles at $t=1000$ displayed in Figs. 8(b) and 8(d). These results clearly confirm the instability of the perturbed solutions, as suggested by the evolution of the total power depicted in Figs. 6(c) and 6(d). Strong instability is observed for all values of $g<0$, which corresponds to the self-focusing.

The findings for the existence and stability of the localized pixels are summarized by diagrams displayed in Fig. 9. To produce them, we analyzed the temporal evolution of total power (4), parallel to monitoring the spatial profile of each solution at large times ($t=100$ and $t=1000$). In Fig. 9, we address three different values of strength $\Omega$ of the trapping potential: (a) $\Omega=2$, (b) $\Omega=4$, and (c) $\Omega=10$. The stability area is represented by gray and white boxes, which correspond, respectively, to robust static outputs and those which feature small residual oscillations, while the parameter area not covered by boxes corresponds to unstable solutions. This includes the area of $g\leq 0$ (self-focusing), where the modes suffer strong instability observed in Figs. 8(b) and 8(d) at $\Omega=10$, but may be stable at $\Omega=2$ and $4$ [in the latter case, the stability domain for $g>0$ is very small, as seen in Fig. 9(b)]. On the other hand, at $g>5$ and $\Omega=10$, the solution undergoes fragmentation under the action of the strong self-defocusing nonlinearity, for all values of $-10\leq\Delta\leq+10$. An example of that is displayed in Fig. 10 for $g=10$ and two extreme values of the mismatch, $\Delta=-10$ and $+10$.

In the stability area, black spots highlight values of the parameters at which the output profiles of the static solutions, observed at $t=1000$, are very close to the respective input profiles, i.e., the VA provides very accurate predictions. Generally, the shape of the stability area in the form of the vertical stripe, observed in Fig. 9(c), roughly follows the vertical direction of the dotted black line in Fig. 1, which pertains to the same value of $\Omega=10$. On the other hand, the expansion of the stability area in the horizontal direction for $\Omega=2$ and $\Omega=4$, which is observed in Figs. 9(a,b), qualitatively complies with the strong change of the curves in Fig. 1 for the same values of $\Omega$. Looking at Fig. 9, one can also conclude that large positive values of $\Delta$ help to additionally expand the stability region.

We stress that the results shown in Fig. 9 are extremely robust: real-time simulations lead to them, even starting with zero input. The input provided by the VA ansatz (8) is used above to explore the accuracy of the VA, which is relevant, as similar approximations can be applied to similar models, incorporating the pump, linear loss, and Kerr nonlinearity (self-defocusing or focusing).

V Vortex solitons

V.1 Analytical considerations: the Thomas-Fermi approximation

In the previous sections, we considered uniform pump field $E$, which generates fundamental modes without vorticity. Here we explore the confined LL model with space-dependent pump carrying the vorticity. It is represented by the driving term

[TABLE]

in Eq. (1), where $\theta$ is the angular coordinate and $E_{0}=\mathrm{const}$. This term naturally corresponds to the pump supplied by a vortex laser beam (with vorticity $1$) vortex beam . In the case of multiple vorticity $m>1$ (which will be considered elsewhere), Eq. (33) is replaced by $E=E_{0}r^{m}e^{im\theta}$.

Patterns supported by the vortex pump correspond to factorized solutions of the stationary version of Eq. (1), taken as

[TABLE]

with complex amplitude $A$ satisfying the following radial equation:

[TABLE]

As an analytical approximation, the TFA for vortex solitons may be applied here, cf. Ref. TFA-vortex . In the general case, the TFA implies dropping the derivatives in the radial equation, which leads to a complex cubic equation for $A$, cf. Eq. (29), under the conditions $\sigma>0$ (self-defocusing) and $\Delta>0$ (positive mismatch):

[TABLE]

Equation (36), as well as its counterpart (29) for the zero-vorticity states, strongly simplifies in the limit of large $\Delta>0$, when both the imaginary and and nonlinear terms may be neglected:

[TABLE]

In particular, the simplest approximation provided by Eq. (37) makes it possible to easily predict the radial location of maximal intensity in the ring-shaped vortex mode:

[TABLE]

Comparison of values given by Eq. (38) with their counterparts extracted from numerically found vortex-ring shapes, which are displayed below in Figs. 12(a) and (b) for $\Delta\geq 0$, demonstrates that the analytically predicted values are smaller than the numerical counterparts by $11\%$ for $\Delta=0$, and by $6\%$ for $\Delta=10$. Naturally, the TFA provides better accuracy for large $\Delta,$ but even for $\Delta=0$ the prediction is reasonable. Furthermore, Eq. (37) predicts a virtually exact largest intensity, $\left|A(r=r_{\max})\right|^{2}$, for the small-amplitude mode displayed in Fig. 12(b).

V.2 Numerical results

Equation (1) with vortex pump profile (33) was numerically solved with zero input. This simulation scenario is appropriate, as vortex states, when they are stable, are sufficiently strong attractors to draw solutions developing from the zero input.

The results, produced by systematic real-time simulations, are summarized in Fig. 13 for $\Omega=2$ in Eq. (1) and $E_{0}=1$ or $2$ in Eq. (33). The figure displays stability areas for the vortex modes in the plane of free control parameters ($\Delta$, $\sigma)$ (the mismatch and nonlinearity strength). It is worthy to note that the stability domain for the self-focusing nonlinearity ($\sigma<0$) is essentially larger than in the diagram for the fundamental (zero-vorticity) modes, which is displayed, also for $\Omega=2$, in Fig. 9. This fact may be naturally explained by the fact that the vanishing of the vortex drive (33) at $r\rightarrow 0$, in the combination with the intrinsic structure of the vortex states, makes the central area of the pattern nearly “empty”, thus preventing the onset of the modulational instability in it.

In the gray areas in Fig. 11, the stable vortex modes have a simple ring-shaped structure, with typical radial profiles shown in Fig. 12. In the case of zero mismatch, $\Delta=0$ [Fig. 11(a)], the vortex state naturally acquires a higher amplitude under the action of the self-focusing. On the other hand, in the case of large positive mismatch [Fig. 11(b)], the small amplitude is virtually the same under the action of the focusing and defocusing, which is explained, as mentioned above, by the TFA that reduces to Eq. (37).

In unstable (white) areas in Fig. 11, direct simulations lead to quick fragmentation of vortically driven patterns into small spots, that feature a trend to developing the above-mentioned critical collapse collapse . A typical example of the unstable evolution is displayed in Fig. 13.

More sophisticated stable vortex profiles are observed in yellow areas in Fig. 13. They are characterized by a multi-ring radial structure, and a spiral shape of the vorticity-carrying phase distribution, as shown in Fig. 14. The yellow areas are defined as those in which the spiral phase field performs a full turn by $360$ degrees, as can be seen in Fig. 14(b). Note that this area exists for both the focusing and defocusing signs of the nonlinearity in Fig. 13(a), and solely for zero nonlinearity in Fig. 13(b), which corresponds to the stronger pump.

The spiral shape of the phase pattern is explained by the fact that radial amplitude $A(r)$ in solution (34) is a complex function, as is explicitly shown, in particular, by Eqs. (2) and (36). The spirality of vortices is a well-known feature of 2D complex Ginzburg-Landau equations GL . However, unlike the present situation, the spirality is not usually related to a multi-ring radial structure. Patterns with multi-ring shapes usually exist as excited states, on top of stable ground states in the same models, being unstable to azimuthal perturbations Javid . For this reason, the stability of complex modes, like the one displayed in Fig. 14, is a noteworthy finding.

Lastly, a typical example of a stable vortex at the boundary between the simple (non-spiral) and complex (spiral-shaped) ones is presented in Fig. 15. It features emerging spirality in the phase field, but the radial structure keeps the single-ring shape.

VI Conclusion

We have introduced the 2D model based on the LL (Lugiato-Lefever) equation with confinement imposed by the harmonic-oscillator trap. In spite of the action of the uniform pump, the confinement creates well localized patterns, which may be used for the creation of robust small-area pixels in applications. The VA (variational approximation), based on a novel fractional ansatz, as well as a simple TFA (Thomas-Fermi approximation), were elaborated to describe the fundamental (zero-vorticity) confined modes. The VA effectively reduces the 2D LL equation to the zero-dimensional version. The VA is additionally enhanced by taking into regard the balance condition for the integral power. The comparison with the full numerical analysis has demonstrated that the VA provides qualitatively accurate predictions, which are also quantitatively accurate, in some areas of the parameter space. The systematic numerical analysis has produced overall stability areas for the confined pattern in the underlying parameter space, which demonstrate that the patterns tend to be less stable and more stable under the action of the self-focusing and defocusing nonlinearity, respectively (although very strong self-defocusing causes fragmentation of the patterns). The increase of the confinement strength leads to shrinkage of the stability area, although it does not make all the states unstable. On the other hand, large positive values of the cavity’s detuning tends to expand the region of the stability in the parameter space.

We have also explored vortex solitons (which may be used to realize vortical pixels in microcavities) supported by the pump with embedded vorticity. In this case, the simple TFA provides a qualitatively correct description, and systematically collected numerical results reveal a remarkably large stability area in the parameter space, for both the self-defocusing and focusing signs of the nonlinearity. In addition to simple vortices, stable complex ones, featuring the multi-ring radial structure and the spiral phase field, have been found too. As an extension of the present work, a challenging issue is to look for confined states with multiple embedded vorticity.

A summary of authors’ contributions to the work: the numerical part has been carried out by W.B.C. Analytical considerations were chiefly performed by B.A.M. and L.S. All the authors have contributed to drafting the text of the paper.

Acknowledgements.

WBC acknowledges the financial support from the Brazilian agencies CNPq (grant #458889/2014-8) and the National Institute of Science and Technology (INCT) for Quantum Information. LS acknowledges for partial support the 2016 BIRD project “Superfluid properties of Fermi gases in optical potentials” of the University of Padova. The work of B.A.M. is supported, in part, by grant No. 2015616 from the joint program in physics between NSF (US) and Binational (US-Israel) Science Foundation.