One- and two-dimensional modes in the complex Ginzburg-Landau equation with a trapping potential

We propose a new mechanism for stabilization of confined modes in lasers and semiconductor microcavities holding exciton-polariton condensates, with spatially uniform linear gain, cubic loss, and cubic self-focusing or defocusing nonlinearity. We demonstrated that the commonly known background instability driven by the linear gain can be suppressed by a combination of a harmonic-oscillator trapping potential and effective diffusion. Systematic numerical analysis of one- and two-dimensional (1D and 2D) versions of the model reveals a variety of stable modes, including stationary ones, breathers, and quasi-regular patterns filling the trapping area in the 1D case. In 2D, the analysis produces stationary modes, breathers, axisymmetric and rotating crescent-shaped vortices, stably rotating complexes built of up to $8$ individual vortices, and, in addition, patterns featuring vortex turbulence. Existence boundaries for both 1D and 2D stationary modes are found in an exact analytical form, and an analytical approximation is developed for the full stationary states.


I. INTRODUCTION
if the local coefficient of the cubic losses grows from the center to periphery faster than |x|. While this setting seems somewhat "exotic", it is interesting to test a possibility to build a more straightforward one, by combining the uniform linear gain and the HO trapping potential, ∼ x 2 , which may approximate deep guiding channels in various photonic settings [23,24]. Indeed, although the linear gain amplifies perturbations far from the localized mode, the HO potential will push the perturbations towards the core area, where the localized mode may have a chance to suppress them with the help of the usual cubic loss. Recently, this possibility was tested in the two-dimensional (2D) model of a spinor (two-component) exciton-polariton condensate, with its components linearly mixed by the spin-orbit coupling (SOC) [25]. It was found that the interplay of the linear gain, HO trap, complex cubic nonlinearity, and SOC gives rise to families of stable 2D dissipative solitons in the form of mixed modes (so called because they combine zero-vorticity and vortex terms in each component), vortex-antivortex complexes, and semi-vortices (which feature vorticity in one component, provided that the Zeeman splitting between the components is present).
However, the stability is only possible if the model also includes an imaginary part (η) of the diffraction coefficient, i.e., spatial filtering [26]. It may represent ionization of the medium [27], or relatively poor finesse of the laser cavity [28]- [30], or diffusion of excitons in semiconductor microcavities [31]- [35]). The stability of the above-mentioned 2D confined modes is secured if η exceeds a certain threshold value.
These results suggest a possibility to consider a setting in a still more fundamental form of the single-component cubic CGLE with the linear gain, spatial filtering (effective diffusion), and the HO potential, in 1D and 2D geometries alike. The model may apply both to the optical laser cavities and microcavities populated by the exciton-polariton condensate. This analysis is the objective of the present work and it will be chiefly performed by means of numerical methods. Nevertheless, the existence boundary for dissipative solitons is obtained in an exact analytical form, in the 1D and 2D cases alike, due to the fact that the linearized form of both the 1D and 2D CGLE, including the spatial filtering and HO potential, admits (novel) exact solutions.
The paper is organized as follows. The models, both 1D and 2D ones, are introduced in Section II, where we also produce exact analytical solutions of the linearized CGLEs, which determine the existence boundaries of the 1D and 2D confined states, and some additional approximate analytical results for nonlinear stationary modes. Numerical results for the 1D and 2D models are reported, severally, in Sections III and IV. It is found that the parameter space of the 1D CGLE is populated by stable stationary modes or breathers (the latter occupying a narrow stripe) and persistent multi-peak quasi-regular states, which fill out the entire spatial domain confined by the HO potential. The 2D model also gives rise to stable stationary modes and breathers, as well as to vortices (both isotropic ones and deformed rotating "crescents") and rotating robust multi-vortex bound states, with the net vorticity up to S = 8. The parameter space of the 2D model also features an area of vortex turbulence.
The paper is concluded by Section V.

A. Complex Ginzburg-Landau equations
All the ingredients of the setting under the consideration, which was outlined in the introduction, are included in the 1D CGLE, with spatial coordinate x, which governs the evolution of local amplitude ψ of the electromagnetic wave: Here η ≥ 0 is the above-mentioned effective diffusivity (in particular, for excitons in a semiconductor microcavity), γ > 0 and Γ > 0 are, respectively, the linear gain and cubic loss, the nonlinearity coefficient takes one of three normalized values, with +1 and −1 corresponding, severally, to the self-focusing and defocusing, and coefficient Ω 2 determines the strength of the trapping HO potential. Variable t in Eq. (1) is the temporal one in the model of the exciton-polariton condensate, or propagation coordinate in the model of the planar laser cavity in optics. On the other hand, a confined 1D laser cavity may be described by Eq. (1) with t remaining the temporal variable, similar to models of the Lugiato-Lefever type [36]. It is relevant to mention that more fundamental models of laser cavities are based on coupled equations of the Maxwell-Bloch type; however, in many case, the atomic degrees of freedom may be adiabatically eliminated, reducing the model to the field-evolution equation of the CGLE type [37].
It is also relevant to mention that, depending on a particular setting, the instability of the zero state in laser cavities may emerge at a finite wavenumber, rather than at the infinitely small one (the short-wave instability, instead of its long-wave counterpart). In that case, the basic nonlinear model is represented not by the CGLE, but by the complex Swift-Hohenberg equation [38]. The latter model may be a subject for a separate work.
Keeping normalization condition (2), one can still perform rescaling of Eq.
Stationary states in the polariton condensate with real chemical potential µ (alias propagation constant, −µ, in the models of planar laser cavities, with t playing the role of the propagation coordinate) are looked for as solutions to Eq. (1) in the form of with complex function u(x) satisfying the equation Stationary states are characterized by the integral power (alias the norm), even if its not a dynamical invariant of the present dissipative model.
The stability analysis for stationary states (4) is based on introducing a perturbed solution, where ε determines the smallness of the perturbation, v(x) and w(x) are complex components of the perturbation eigenmode, and complex λ is the respective stability eigenvalue. The substitution of ansatz (7) in Eq.
(1) and linearization with respect to ε leads to the system of stationary linear equations: Instability is accounted for by eigenvalues with Re(λ) > 0.
The 2D model is based on the corresponding version of Eq. (1), which is written, in terms of the polar coordinates, (r, θ), as with coefficients subject to conditions (2) and (3). Its stationary solutions with integer vorticity S = 0, 1, 2, ... are looked for as where complex function u(r) satisfies the radial equation: with boundary condition u(r) ∼ r S at r → 0. The integral power (norm) of the 2D state is For the stability test, perturbed 2D solutions were looked for as [cf. Eq. (7)] where integer m is an independent angular index of the perturbation mode, and the radial equations for perturbation amplitudes are This eigenvalue problems defined by Eqs. (8) and (14) were solved by means of numerical methods.

B. Analytical solutions
The amplitude of stationary states vanishes at their existence boundary, where, accordingly, one-dimensional equation (5) is replaced by its linearized version, Straightforward consideration of Eq. (15) demonstrates that it gives rise to an exact groundstate (GS) solution, with arbitrary amplitude A 0 , provided that diffusivity η is related to linear gain γ as follows: This exact solution follows the commonly known pattern of the GS wave function of the one-dimensional HO potential in quantum mechanics. However, the complex coefficient in front of the second derivative in Eq. (15) introduces an essential difference, in the form of the spatial chirp of the wave function (a phase term ∼ x 2 ).
Equation (18) Similarly, the linearized version of the 2D stationary equation (11), gives rise to exact solutions for both the GS (S = 0) and vortices (S ≥ 1): the corresponding threshold being These exact solutions follow the pattern of the GS and eigenstates carrying the orbital angular momentum for the two-dimensional HO potential, the difference from the quantummechanical wave functions being the presence of the radial chirp, i.e., a phase term ∼ r 2 .
At γ > γ (2D) thr , the parameter space of the 2D nonlinear model is populated by several types of robust static and dynamical modes, as shown below, while, in the absence of the nonlinearity (σ = Γ = 0), an exponentially growing solution to Eq.
For the 1D nonlinear model (1), an approximate analytical solution can be found in the case of σ = 0, Γ > 0 (i.e., if the nonlinearity is represented solely by the cubic loss), treating the gain and loss terms as small perturbations. In the zero-order approximation, the stationary solution is represented by the GS of the HO potential, the first-order approximation, by the power-balance condition: The substitution of the zero-order wave function in Eq. (24) predicts the squared amplitude, The norm (6) of the solution with amplitude (25) is In particular, Eq. (26) predicts that nonzero states exist at γ > γ thr = Ωη/4. For small η, this result agrees with the exact threshold given by Eq. (18).
A similar prediction can be elaborated for 2D vortex modes produced by Eq. (11) with σ = 0 and Γ > 0. In this case, the zero-order wave function is u The substitution of u which corresponds to the 2D integral power (12) cf. Eq. (26). The threshold condition, following from Eq. (29), γ > γ thr (1 + S) /2, in the limit of small η is consistent with the exact 2D threshold condition (23).
The analytical results, both exact and analytical ones, are compared to their numerical counterparts below.

III. BASIC RESULTS: THE 1D SETTING
Numerical solutions of stationary equations (5) and (11) were constructed by means of the squared-operator iterative method [39]. The evolution governed by Eqs.
(1) and (9) was simulated using the fourth-order split-step Fourier-transform algorithm. The integration domain with dimensions −8 < x, y < +8 was sufficient to accommodate all confined modes investigated herein. In most cases, the numerical grid of size 128 × 128 was sufficient to produce the modes with necessary accuracy. The stability and instability of all stationary modes, predicted by the numerical solution of linearized equations (8), was accurately corroborated by direct simulations. As said above, control parameters of the model are γ = Γ [see Eq. (3)], η, and Ω 2 , as well as the discrete nonlinearity coefficient, σ = −1, 0, +1 in Eq. (2).
First, in the absence of losses, γ = Γ = η = 0, it is easy to construct a family of real GS solutions of Eq. (5) for either sign of σ, which correspond to stable solutions of Eq.
(1). For σ = +1 and −1 (the self-focusing/defocusing cubic term), the shape of the GS is close, respectively, to that of the usual bright soliton placed at the bottom of the HO potential trap, or to the GS predicted by the Thomas-Fermi (TF) approximation [40]. In  Proceeding to the 1D model with γ = Γ = 1, we note, first, that all the stationary modes are unstable in the absence of the effective diffusion, η = 0, similar to what was reported in Ref. [25]. This instability transforms static modes into quasi-turbulent states which fill the entire space admitted by the HO trapping potential, as shown in Fig. 2. quasi-regular multi-peak periodically oscillating states, developing from unstable stationary modes at small η [ Fig. 4]. In the limit of η → 0, the latter pattern goes over into the quasi-turbulent one displayed in Fig. 2. The increase of η naturally leads to simplification of the quasi-regular patterns via the decrease of the number of peaks, from seven in Fig.   4(a) to five in 4(b), to three 4(c) and, eventually, to one in Fig. 4(d). The latter dynamical mode is close to a regular breather, cf. Fig. 3(b).
The results obtained in the case of self-focusing, σ = +1, are summarized in Fig. 5, where panel (a) shows the numerically found dependence of the integral power, N on the crucially important control parameter, η, and panel (b) presents the findings in the form of stability map displayed in the plane of γ = Γ and η. In particular, we stress that a numerically found existence boundary for stable GSs, shown by the black dashed line in (b), precisely coincides with the analytical prediction given by Eq. (18), which is shown by the red dashed line.
Further, we observe that vast areas populated by stable GSs and persistent quasi-regular states are separated by a narrow sliver supporting stable breathers. We also stress that no  Other parameters are the same as in Fig. 3. by the cubic loss, Γ > 0. These conclusions are illustrated by Fig. 6, which summarizes the results for σ = 0 and Ω 2 = 6, in the same way as it is done in Fig. 5 for σ = +1 and Ω 2 = 2.
In particular, the magenta line in panel 5 demonstrates that the analytical approximation, elaborated for σ = 0 in the form of Eq. 26), is an accurate one. Overlapping black and red dashed curves display, severally, the numerically found and analytically predicted [see Eq.

IV. BASIC RESULTS: THE 2D SETTING
In the 2D model, the systematic numerical analysis has identified several types of stationary and dynamical states. First, similar to the 1D case, the largest stability area is found for stationary GS states, which are trapped isotropic modes, see an example in Fig.   7. Next, persistent axisymmetric breathers, developing from unstable stationary states, are found too, although in a small parameter area (see below). An example of the breather is displayed in Fig. 8 The decrease of η leads, as in the 1D model, to destabilization of the stationary modes.
The evolution of unstable modes generates a variety of different dynamical states featuring vortex structures, cf. Refs. [41,42], where the 2D CGLE with the cubic-quintic nonlinearity and a radially localized linear gain (but without a trapping potential and diffusion) also gave rise to a large number of different vortical modes; stable complexes built of one, two, or three vortices were reported too in the study of the cubic-quintic CGLE model with the uniform linear loss and trapping HO potential [14]. Among such structures, we first identify stable gives rise to the above-mentioned vortex-turbulent state, an example of which is presented in Fig. 11.
Numerically measured characteristics of the rotating complexes displayed in Fig. 10, viz., the values of their radius, defined as the distance of maxima of the local intensity from the rotation pivot, and the angular velocity of the rotation, are collected in Table I In all panels, spatial scales are the same as in Fig. 9(a).    Table I. As concerns the angular velocity, ω, of the rotating crescent, we note that, for a narrow ring of radius ρ, the linear Schrödinger equation produces an elementary result, ω = S/ρ 2 . In reality, the crescent, displayed in the top row of Fig. 10, is not quite narrow, and the comparison of this formula with the respective angular velocity given in Table I suggests that ω is determined by the inner layer of the crescent, which has ρ 0.5. It is relevant to note that the rotation velocity of the crescent-shaped vortex with S = 1 is essentially higher than in all the complexes with S ≥ 2.
Note that the sequence of multi-vortex complexes, displayed in Fig. 10 for σ = +1 (the self-focusing sign of the nonlinearity), does not include a bound state of five vortices.
Actually, such a stable state can be found too, but for σ = −1, as shown in Fig. 12.
Results produced by systematic numerical analysis of the 2D model are summarized in the plots are quite similar to the ones displayed here. Furthermore, for σ = −1, the plots are also similar to what is shown in Fig. 13(b) for σ = 0. We stress that the analytical existence boundary for the stationary mode, predicted by Eq. (23), is completely identical to its numerically found counterpart. Comparing these stability maps with their counterparts reported above for the 1D model, cf. Figs. 6 for σ = 0 and 5 for σ = +1, we conclude  that the axisymmetric vortices with S = 1 and the vortex-turbulent states play, roughly speaking, the same role as, respectively, breathers and quasi-regular patterns in 1D, while the multi-vortex bound states do not have their 1D counterparts. FIG. 14: The same as in Fig. 13, but for σ = +1.

V. CONCLUSION
We have demonstrated a new mechanism for the stabilization of confined modes in the 1D and 2D CGLEs (complex Ginzburg-Landau equations) with the cubic-only nonlinearity and the spatially uniform linear gain. Namely, we have found that the background instability is suppressed by the combination of the HO (harmonic-oscillator) trapping potential and effective diffusion. Analytical solutions for the linearized version of the CGLE provide existence boundaries for trapped modes, which are exactly confirmed by numerical results. Our numerical analysis has produced a set of stable stationary modes, persistent breathers, and more complex patterns, viz., quasi-regular multi-peak structures in 1D, and vortices (axisymmetric and rotating deformed crescent-shaped ones), and rotating multi-vortex complexes, with the net topological charge up to S = 8, in 2D. The models addressed in this work apply to laser cavities and semiconductor microcavities with exciton-polariton condensates.