Vortex modes supported by spin-orbit coupling in a laser with saturable absorption

We introduce a system of two component two-dimensional (2D) complex Ginzburg-Landau equations (CGLEs) with spin-orbit-coupling (SOC) describing a wide-aperture microcavity laser with saturable gain and absorption. We report families of two-component self-trapped dissipative laser solitons in this system. The SOC terms are represented by the second-order differential operators, which sets the difference, $|\Delta S|=2$, between the vorticities of the two components. We have found stable solitons of two types: vortex-antivortex (VAV) and semi-vortex (SV) bound states, featuring vorticities $\left( -1,+1\right) $ and $\left( 0,2\right) $, respectively. In previous works, 2D localized states of these types were found only in models including a trapping potential, while we are dealing with the self-trapping effect in the latteraly unconfined (free-space) model. The SV states are stable in a narrow interval of values of the gain coefficients. The stability interval is broader for VAV states, and it may be expanded by making SOC stronger (although the system without SOC features a stability interval too). We have found three branches of stationary solutions of both VAV and SV types, two unstable and one stable. The latter one is an attractor, as the unstable states spontaneously transform into the stable one, while retaining vorticities of their components. Unlike previously known 2D localized states, maintained by the combination of the trapping potential and SOC, in the present system the VAV and SV complexes are stable in the absence of diffusion. In contrast with the bright solitons in conservative models, chemical potentials of the dissipative solitons reported here are positive.

which produce Townes solitons with zero vorticity [32], and their ring-shaped vorticity-carrying extensions [33]- [36]. The instability of the Townes solitons is driven by the critical collapse in the 2D space with the cubic attraction [37,38], while the vortex-ring solitons are subject to a still stronger ring-splitting instability [39,40]. The presence of the SOC terms changes the situation, as they bring a coefficient which fixes a length scale in the system, which, in turn, breaks the scaling invariance that makes norms of the Townes' solitons degenerate (the entire family of these solitons has a single value of the norm). The lifting of the norm degeneracy pushes it below the degenerate value, which determines the collapse-onset threshold, hence the collapse cannot occur anymore. This mechanism secures the stabilization of both the SVs and MMs, and lends them the role of the otherwise missing ground states [30,31].
A difference introduced by the microcavity SOC terms is their second-order derivatives, giving rise to the vorticity difference between two spin (circularly polarized) components to be |∆S| = 2, rather than |∆S| = 1 as in the abovementioned models of atomic BEC. In particular, the analysis of the composite modes in the exciton-polariton system, which includes a harmonic-oscillator trapping potential, was recently elaborated in Ref. [41]. An essential result was the identification of stability regions for MMs and vortex-antivortex (VAV) complexes, with vorticities S = ±1 in the two components (thus complying with the above-mentioned constraint, ∆S = 2). An area of the MM-VAV bistability was identified too, and stable SVs were found when the Zeeman splitting was added to the two-component system.
While trapping potentials are an important ingredient in various systems, see, e.g., [42,43] in addition to the above, it is also relevant to construct stable composite modes supported by SOC in the laterally unbound settings. In particular, laser systems with and without saturable absorbers [20,21,[25][26][27] are known to support self-bound vortex rings [21,27]. A possible practical realisation of a laser model considered below is a wide-aperture semiconductor laser with a saturable absorber section similar to the ones used in [20,25,26], where the SOC effects originate from the above mentioned TE-TM splitting of the cavity resonances. Theoretically, but without the SOC effects, the model used below was elaborated in a series of papers by Rosanov and co-authors, see, e.g., [44][45][46], and references there in. We generalise the Rosanov's approach by including the SOC terms represented by second derivatives mixing the two spin components. The system may also include effects of diffusion, but, unlike the setting explored in Ref. [41] and in many other models [48]- [50], the presence of the diffusion is not necessary for the stability of 2D states in the present case. Below, we identify two species of stable self-bound vortex modes, viz., VAVs and SVs, both built of components obeying the constraint of ∆S = 2. In particular, the VAVs feature two mutually symmetric components with vorticities S = ±1, while the amplitude of the SV's zero-vorticity component is much larger than the amplitude of its vortical counterpart.
It is relevant to mention that VAVs are quite similar to two-component states, produced by systems of 2D [51][52][53] and 3D [54] nonlinearily-coupled GPEs, with opposite vorticities, ±S, and identical density profiles in the components. In the latter context, the bound states of this type were called "hidden-vorticity modes", as, on the contrary to their counterparts with equal vorticities in both components, the total angular momentum of the hidden-vorticity states is zero. In addition to that, conservative systems with the nonlinear (cross-phase-modulation) interaction between two components (and no linear coupling between them) maintain bound states similar to SVs, with zero and nonzero vorticities in the components [55,56].
The subsequent presentation is organized as follows. The system of complex Ginzburg-Landau equations (CGLEs) with the saturable gain, coupled by the SOC terms, is introduced in Section II. The same section also reports some simple analytical results (necessary conditions for the existence of stable 2D dissipative solitons). Numerical results, produced by the systematic investigation of the 2D system, are reported in Section III, for two above-mentioned species of stable modes, viz., VAVs and MMs. The paper is concluded by Section IV.

II. THE MODEL AND ANALYTICAL ESTIMATES
In terms of scaled time t and coordinates (x, y), the system of CGLEs for wave functions ψ ± of the two spin components, with the linear loss, whose coefficient is scaled to be 1, saturable gain and saturable absorption, which are represented, respectively, by coefficients g > 0 and a > 0, and SOC linear mixing with real coefficient β, is written as Here, the coefficient in front of the Laplacian, which represents the diffraction, is also scaled to be 1 (in this notation, physically relevant values of the SOC coefficient, |β|, are definitely smaller than 1), and positive ε < 1 defines a relative saturation strength of the gain and absorption. The dispersive component of the linear loss (diffusion) may be present, with coefficient η > 0, but, in fact, the inclusion of this term, whose physical origin may not be obvious, is not necessary for producing stable 2D modes, therefore we set η = 0 in what follows below. The generic situation may be adequately represented by parameters which implies weak saturation of the gain in comparison with absorption, while the gain and SOC strengths, g and β, will be varied as physically relevant parameters controlling modes generated by the system. Stationary states with real chemical potential µ are sought as where complex stationary wave functions u ± obey the following equations: For the analysis of stability of stationary states, we define perturbed solutions as where δ is a real infinitesimal amplitude of the perturbations, λ is a complex eigenvalue, and complex perturbation eigenmodes, v ± (x, y) and w ± (x, y), satisfy the linearized equations where we have defined a real function, As usual, the stability condition is Re(λ) ≤ 0 for all eigenvalues.
The system of equations (10)-(13) can be written as an eigenvalue problem in the matrix form: where The spectrum of the linear-stability operator L was constructed by means of the Fourier collocation method.
To complete the formulation of the stability-analysis framework, we notice that, being interested in stable dissipative solitons, a necessary condition is the stability of the zero background in Eqs. (1), (2) and ( 3), which obviously amounts to condition g < 1 + a. On the other hand, a necessary condition for the ability of the saturable gain to maintain nontrivial modes is that the largest value of f (n ≡ |ψ + | 2 + |ψ − | 2 ) in Eq. (3), which is attained at density must be positive. The substitution of n 0 in Eq. (3) yields a lower bound for g, which, if combined with the abovementioned upper one, g < 1 + a, defines the interval in which the gain coefficient may take its values: [the compatibility condition for a and ε, following from Eq. ( 23), ]. An additional restriction on g is imposed by the condition that expression (23) must be positive too: Note that for values a = 2 and ε = 0.1 adopted here, interval ( 24) amounts to while interval (25) is much broader and may therefore be disregarded. Our objective is to constructs solutions of the CGLE system of Eqs. (1 ) and (2) in the form of 2D bright solitons morphed as bound states of two components with certain values of integer vorticities, m − 1 and m + 1, so that they comply with the above-mentioned constraint, ∆S = 2. In polar coordinates (r, θ), the relevant solutions with real chemical potential µ can be defined as with complex amplitude functions φ ± (r) satisfying the radial equations: where f is defined by Eq. (3), with |ψ ± | replaced by |φ ± |. The boundary condition for Eq. (28) at r → 0 is that φ ± must be vanishing as r |m∓1| at r → 0, except for the case of m ∓ 1 = 0, when the boundary condition is dφ ± /dr| r=0 = 0. At r → ∞, soliton solutions must feature the exponential decay, with λ r > 0 and, generally, a nonvanishing imaginary part of the decay rate, λ i . The substitution of asymptotic expression (29) in Eq. (28) and the linearization for the exponentially small amplitude functions leads to the relation between the imaginary and real parts, and a quadratic equation for λ 2 r , a larger root of which yields a relevant value, λ 2 r > 0 (here, sign ± is unrelated to the subscript of φ ± ). An essential consequence of Eq. (31) is that, while in the conservative model, in which the linear-loss factor, (1 + a − g), does not appear, λ 2 r > 0 is only possible for a negative chemical potential, µ < 0 (at least, in the physically relevant case of |β| < 1), in the dissipative system µ > 0 is admitted too. Indeed, all the numerical solutions, reported below, have been obtained with µ > 0. Further, it should be stressed that, as in the case of generic dissipative solitons [57]- [59], relevant solutions exist only at isolated (positive) eigenvalues of µ. The solitons are characterized by the total integral power (norm), The analysis is carried out below for two most essential species of the soliton complexes, viz., VAVs corresponding to m = 0 (so called because their components carry opposite vorticities, −1 and +1), and SVs, corresponding to m = 1, with component vorticities 0 and 2. For m ≥ 2, bound states (27) resemble "excited states" of SVs, introduced in the model of the atomic BEC in Ref. [29] , and, as well as in that setting, it is plausible that they are completely unstable.

III. NUMERICAL RESULTS
Numerical solutions for the soliton modes were obtained by splitting complex functions φ ± (r) exp [i(m ∓ 1)θ], defined in Eq. (27), into real and imaginary parts, and solving the resultant equations in the Cartesian coordinates by means of the modified squared-operator method [60]. The respective eigenvalues µ were found simultaneously with the stationary modes. Stability of the stationary states was identified by means of systematic direct simulations of perturbed evolution of the stationary states, governed by Eqs. (1) and (2), and, in parallel, by calculating the stability eigenvalues, as determined by Eqs. (15) and (16).
A. Vortex-antivortex (VAV) complexes, m = 0 VAV modes, with m = 0, were constructed as solutions of the stationary equations, starting with the input adjusted to vorticities ±1 in the two components, with some empirically chosen real parameters φ 0 and α, a typical value being α = 0.05. As said above, parameters η = 0, a = 2, ε = 0.1, were fixed in Eqs. (1), (2) and (28), while the SOC and gain strengths β and g were varied. As a result, three coexisting VAV families were found, two unstable and one stable. Figure 1 represents them by showing their integral power ( 32) and the peak density, vs. g, at fixed β = 0.1. In particular, the pair of unstable branches (continuous and dashed blue lines in Fig. 1) emerge from a bifurcation point at The stable (red) branch is shown only in the relatively narrow window where it remains stable, 2.047 < g < 2.097, which is bounded by vertical dashed lines in Fig. 1. Note that, although being much more narrow than interval (24) determined by the above-mentioned necessary conditions, stability region (36) is located quite close to the center of the broad interval (24). It is relevant to note that, the narrow interval (36), and similar intervals reported below, can be realized in terms of actual physical parameters of laser cavities with the saturable gain and loss, as it follows from Refs. [21,[25][26][27].
As concerns the chemical potential of the stable branch, in the stability interval (36) it decreases nearly linearly, as a function of g, from µ (g = 2.047) = 0.127 to µ (g = 2.097) = 0.081 . We stress that, as mentioned above, these values are positive, on the contrary to necessarily negative chemical potential for solitons in conservative models. As for the two unstable branches, they produce positive and almost constant µ in the same interval (36).
At g < 2.047, the development of the instability of the branches shown by the blue lines in Fig. 1 leads to their decay towards the zero solution, while at g > 2.097 the amplitude features exponential growth (blowup, see Fig. 4). Inside stability interval (36), the stable branch (the red segment in Fig. 1) is an attractor : solitons belonging to either unstable branch spontaneously transform into ones belonging to the stable branch, as shown in Fig. 2. The outcome of the transformation is identical to the VAV mode that may be found as the stationary solution, see an example in Fig. 3. We did not aim to extend the red branch in Fig. 1 into the areas where it is unstable, i.e., g < 2.047 and g > 2.097, see Eq. (36).
As concerns (in)stability eigenvalues, produced by the linearized equations ( 15) and (16), typical examples of stable and unstable spectra are displayed in Figs. 2(b), 3(b), 4(b) and 7(b), respectively. These results are consistent with the conclusions made on the basis of systematic direct simulations of the perturbed evolution of the VAV modes.
The results of the investigation of the VAV modes are summarized in Fig. 5, which displays the stability region in the plane of the varying parameters, g and β (the gain and SOC strengths). It is worthy to note that the instability proceeds via the decay or blowup, without breaking the axial symmetry of the VAVs. This is a drastic difference  from the above-mentioned hidden-vorticity modes, which are also built as bound states of localized components with vorticities S = ±1 in the framework of the conservative system of nonlinearly coupled GPEs, and are partly [51][52][53] or fully [54] unstable against the splitting instability, which breaks the respective vortex ring into fragments. Figure 5 demonstrates that the stability window for the VAV modes remains narrow in comparison with interval (26) defined by the above-mentioned coarse necessary conditions. The window exists too at β = 0: 2.040 < g < 2.090,  Composite soliton modes corresponding to ansatz (27) with m = 1 were found, along with the corresponding eigenvalues µ, as stationary solutions initiated by input  Unlike input (33) which generates the VAV modes, here the amplitude of the vortex component is taken essentially smaller than in the zero-vorticity component, because the vortex component in established SV modes tends to have a relatively small amplitude [29]- [31]. The systematic numerical analysis yields a stability band for the SV modes shown in Fig. 6, which demonstrates a situation generally similar to that displayed for VAV modes in Fig. 1, but with a stability window, 2.095 < g < 2.120, (39) whose width is half of that defined by Eq. (36) for the VAVs. This interval, similar to its counterpart (36), stays close to the center of the broad interval (26), which is defined by the necessary conditions considered above.
In the present case too, three branches of stationary solutions are observed in Fig. 6, viz., two unstable ones, which emerge at the bifurcation point, that virtually exactly coincides with its counterpart for the VAV states, given by Eq. (35), and a stable branch shown by the red segment in Fig. 6. Another similarity to the case of VAV modes is that, outside the stability band (39), unstable SVs decay to zero at g < 2.095, and undergo blowup at g > 2.120. On the other hand, the comparison of Figs. 1 and 6 shows that the integral power of the stable SVs is smaller than the power of the stable VAVs, roughly, by a factor of 5.
The chemical potential of the stable branch decreases nearly linearly, as a function of g, in the stability interval (36), from µ (g = 2.095) = 0.159 to µ (g = 2.120) = 0.106, while µ remains nearly constant in interval (39) for both unstable branches. Note that all these values of the chemical potential are positive, as they are for the VAV branches considered above.
As well as in the case of VAV modes, the stable branch is, within the limits imposed by Eq. (39), an attractor, as initial states corresponding to either of the two unstable families spontaneously transform into the stable counterpart, keeping the SV structure. A typical example of the spontaneous transformation is displayed in Fig. 7.

C. The system without the SOC terms
At β = 0, the vortex component of the SV vanishes, ψ − = 0, while the zero-vorticity one turns into am axisymmetric dissipative soliton produced by the single equation (1):  Fig. 6, with the spontaneous transition between them designated by the vertical arrows.
(recall we actually set η = 0). In the case of β = 0, essentially the same equation, with ψ + replaced byψ + , applies to the two-component system with equal components, of substitution ψ ± = (1/ √ 2)ψ + . The numerical solution of Eq. (40) produces the respective stability window, 2.090 < g < 2.130, (41) which is somewhat broader than its counterpart (39), found for β = 0.1. The relative narrowness of the stability windows (39) and (41) in comparison with the ones given by Eqs. (36) and (37) for VAVs suggests that SV modes are more fragile in comparison with the VAVs. This expectation is confirmed by the fact that, on the contrary to the situation for the VAVs, whose stability area tends to expand with the growth of the SOC strength β in Fig. 3, the increase of β leads to shrinkage of the SV stability window, which closes and does not exist at β ≥ 0.24 (not shown here in detail). Finally, Fig. 8 summarizes the findings in the plane of (a, g) for fixed ε = 0.1. It is clearly seen that the region of stable zero vorticity soliton is relatively broad at large values of a. Typical exaple of the spontaneous transformation of unstable zero vorticity soliton is displayed in Fig. 9.

IV. CONCLUSION
The objective of this work is to construct 2D self-trapped states (dissipative solitons) in the system of CGLEs with saturable gain and absorption for SOC (spin-orbit-coupled) modes of an optical microcavity with saturable gain and loss, in the absence of any geometric trapping. The second-order linear differential operator representing SOC creates complexes with difference ∆S = 2 between vorticities of the two components. Previously, 2D localized dissipative vortex states were predicted, under the action of SOC, in the presence of a trapping potential, but they were not found in the free space. Here, we have identified stable complexes of the VAV (vortex-antivortex) and SV (semi-vortex) types, i.e., ones with vorticities (−1, +1) and (0, 2) in the two components. The 2D solitons of the latter type are quite fragile, being stable in a very narrow (but, nevertheless, existing) window. For the VAVs, the stability interval, defined in terms of the gain coefficient, is rather narrow too, but it may be expanded by applying stronger SOC.
It is worthy to note that, on the contrary to the previously reported 2D localized states of the VAV, MM (mixedmode), and SV types, supported by the trapping potential [41], as well as single-component trapped [50] and selftrapped [48,49] vortices, the stability of the VAVs and SVs in the present system (including its simplified form which does not include SOC) does not require the presence of diffusion terms (dispersive linear losses, whose physical origin may be problematic) in the CGLE model. Due to the absence of the diffusion, the 2D dissipative system introduced in this work keeps the Galilean invariance, which suggests a possibility to introduce moving solitons and simulate collisions between them (cf. Ref. [61]), as well as a possibility of forming bound states of two or several dissipative solitons. These issues should be a subject of a separate work.