Spin-orbit coupling and nonlinear modes of the polariton condensate in a harmonic trap

We consider a model of the exciton-polariton condensate based on a system of two Gross-Pitaevskii equations coupled by the second-order differential operator, which represents the spin-orbit coupling (SOC) in the system. Also included are the linear gain, effective diffusion, nonlinear loss, and the standard harmonic-oscillator trapping potential, as well as the Zeeman splitting. By means of combined analytical and numerical methods, we identify stable two-dimensional modes supported by the nonlinear system. In the absence of the Zeeman splitting, these are mixed modes, which combine zero and nonzero vorticities in each of the two spinor components, and vortex-antivortex complexes. We have also found a range of parameters where the mixed-mode and vortex-antivortex states coexist and are simultaneously stable. Sufficiently strong Zeeman splitting creates stable semi-vortex states, with vorticities 0 in one component and 2 in the other.

the first-order derivatives, in the leading approximation [3,19]. Thus, depending on physics of the SOC system under the consideration, vorticities in the two spinor components differ by either ∆S = 2 or 1, for the second-and first-order-derivative coupling terms, respectively.
Adding trapping or lattice potentials to the nonlinear SOC systems is a very active direction of the current research, see, e.g., Refs. [21][22][23]. In the polariton microcavities, which are the subject of the present work, a variety of practical methods have been developed to secure high quality of the engineered trapping landscapes [23]. As concerns polaritons in a simple single potential well, there are only a few studies that addressed nonlinear modes of these traps, largely disregarding SOC, see, e.g., Refs. [24][25][26][27]. Nevertheless, a recent experimental work on open polaritonic resonators supplies a system with strong SOC effects [28].
In this work, we address an effectively 2D polariton condensate modeled by a system of two GPEs with the SOC terms represented by the above-mentioned second spatial derivatives. Usual ingredients of microcavity models, viz., the linear gain, viscosity (diffusion of the effective fields), nonlinear cubic loss, and an isotropic harmonic-oscillator trapping potential, are included. In the general case, the Zeeman splitting between the two fields is present too. Our objective is to identify stable 2D states in this system, which turn out to be mixed modes and vortex-antivortex bound states, that tend to be stable, severally, under the action of weak and strong SOC, with a small bistability area in the parameter space. These results for the symmetric system, which does not include Zeeman splitting, are obtained by means of analytical and numerical methods in Section II. In addition to that, in Section III we demonstrate that sufficiently strong Zeeman splitting creates stable semi-vortex modes, with vorticities 0 and 2 in its components. The latter states are often referred to as the half-vortices in the polariton context [20,30,31], and are also known in models without SOC [30], or with the spin-only coupling (direct linear mixing of the Rabi-coupling type) [27]; see also Ref. [32], as concerns vortex lattices in a scalar model, which does not include spin effects. The paper is completed by a summary in Section IV.

A. Basic equations
The model that we consider in this work consists of coupled GPEs for the two-component polariton spinor wave function (Ψ + , Ψ − ) [20,24]. We assume that the polariton condensate is pumped by means of a nonresonant optical scheme, which creates a density of coherent and incoherent excitons [29,38]. Adiabatic elimination of the latter yields the following system of coupled equations [38]: Here, gain ε > 0 is determined by the rate of the pumping of the polariton density, σ > 0 is the gain-saturation coefficient that prevents an unphysical blow-up and allows for the time-independent finite density solutions to exist in the system, and η is an effective polariton diffusion coefficient (viscosity), linked to the carrier diffusion [38,39]. Further, the strength of the repulsive intra-component interactions is normalized to be 1, while α = −0.05 accounts for the inter-component attraction, and β is the strength of the polariton SOC [20], while the strength of the harmonicoscillator potential is scaled to be 1. The potential corresponds to the shift of the cavity detuning, and originates from the trapping of the photonic component of exciton-polariton condensate. The potential can be created, for example, if either the above-mentioned open cavity with a parabolic top mirror is used, or one of the many available patterning techniques is applied to the top mirror in the traditional closed-resonator setting with Bragg reflectors [23]. Lastly, Ω characterizes the Zeeman-splitting energy, which is proportional to the applied magnetic field. Accessible physical values of the parameters in the linear Hamiltonian part of Eqs. (1), (2), and, in particular, the SOC strength in the presence of the harmonic-oscillator potential have been discussed in detail in Ref. [28]. That work demonstrates that SOC energy may reach 1meV, which is comparable to the characteristic harmonic-oscillator energy, thereby making values of β between 0 up to 1 physically relevant. Achieving large polarization energy splitting has been important, since the beginning of the studies of microcavity polaritons, in the context of the polariton spin Hall effect [33,34], see also Ref. [4] for a detailed review of the previous work. Currently, the design of strong SOC is getting crucially important for the work towards experimental observation and applications of polariton topological insulators [5,[35][36][37]. Stationary solutions to Eqs. (1) and (2), with chemical potential µ, are sought in the form of where ψ ± are complex functions of polar coordinates (r, θ). The substitution of expressions (3) in Eqs. (1) and (2) leads to the following stationary equations: which are used below as the basic system.

B. Linear and quasi-linear modes: vortex-antivortex states
Before tackling the full nonlinear system, it is relevant to address the simplest version of stationary equations (4) and (5), which neglects the nonlinearity, gain, and diffusion, as well as assumes the symmetry between the components, i.e., Ω = 0 (no Zeeman splitting): It is easy to see that Eq. (6) gives rise to two exact eigenmodes of the linearized dissipation-free system, corresponding to independent signs ± in Eq. (7): where A max is an arbitrary constant, and superscript (0) attached to µ ± refers to the linear approximation. These exact eigenmodes, which, as a matter of fact, are found following the pattern of standard solutions for the 2D isotropic harmonic oscillator in quantum mechanics, may be naturally called vortex-antivortex states. Note that the symmetric vortex-antivortex state with the lower eigenvalue of the energy, µ exists if the SOC is not too strong, viz., at β ≤ 0.5, while its antisymmetric counterpart with higher energy, µ 1 + 2β, exists at all values of β. In fact, the disappearance of the eigenmode (7) corresponding to the top sign in ∓ at β = 0.5, where the radial size of the eigenmode shrinks to zero as (1 − 2β) 1/4 , may be considered as a quantum phase transition in the polariton condensate, cf. Ref. [40].
Both types of the vortex-antivortex states, antisymmetric and symmetric ones, were experimentally observed, in a quasi-linear regime, in Ref. [28], where they were named, respectively, modes of types (i) and (iii). In terms of the present notation, the respective SOC coupling, estimated in Ref. [28], was β ≈ 0.04.
In a sense, the exact vortex-antivortex modes are similar to the above-mentioned semi-vortices found in the nonlinear dissipation-free 2D SOC model of the atomic condensate considered in Ref. [7], which feature zero vorticity in one component and vorticity ±1 in the other. Indeed, the difference of the vorticities in the components of the semivortex is ∆S = 1 because (as mentioned above too) SOC in the atomic condensates is represented by the linear differential operator of the first order [3], while in the present system the difference between the components of the vortex-antivortex is ∆S = 2, as the SOC is represented by the second-order operator in Eqs. (1) and (2).
If the nondissipative nonlinear terms are kept in Eqs. (4) and (5) as small perturbations, it is easy to find the first nonlinear correction to the chemical potential, µ (1) ± , using the commonly known method in quantum mechanics, that generates the first correction as the spatial average of the perturbation potential [18,41]. A simple calculation, which makes use of the unperturbed wave functions (7), yields where A 2 max should be small enough, to keep correction (9) small in comparison with the result (8), obtained in the linear limit.
Furthermore, if the linear gain, diffusion, and nonlinear loss are also kept as small perturbations in Eqs. (4) and (5), then amplitude A max of the linear mode (7), which is treated in Eqs. (7) and (9) as a free parameter, is uniquely selected by the obvious balance condition for the total norm (it is defined below by Eq. (14)), written as for the mode |ψ + (r)| = |ψ − (r)| taken as per Eq. (7). The result of the simple calculation is which makes sense if the pump rate ε is small enough. A more sophisticated analytical approximation for the vortexantivortex in the full system is developed below, see Eqs. (16)- (18).

C. Mixed-mode states
In addition to the two species of vortex-antivortex complexes, built as per Eqs. (7), the linearized coupled GPEs (6) give rise to another type of eigenstates, which, following Ref. [7], may be called mixed modes, as they contain terms with zero vorticity and vorticities ±2 in each component. Like the symmetric vortex-antivortex states (7), the mixed modes obey constraint ψ + = ψ * − , with * standing for the complex conjugate. Linear mixed modes cannot be represented by an exact solution (unlike Eq. (7) for the vortex-antivortex complexes), but they can be easily constructed in an approximate form for β 1, starting from an obvious ground state of the harmonic-oscillator in each component at β = 0: At larger β, the mixed mode eigenmodes were constructed numerically by means of the imaginary-time method, applied to the time-dependent version of linearized equations (6), starting with the ground-state input, Ψ + = Ψ − = exp(−r 2 /2). The numerically found mixed-mode shapes are displayed in Fig. 1, and the comparison with the analytical approximation (12), for a relevant small value β = 0.1, is presented in Fig. 1(b). Further, eigenvalues corresponding to the mixed modes were compared to the analytically found eigenvalues for the vortex-antivortex states, given by Eq. (8), with the aim to identify the ground state, which must have the smallest eigenvalue of energy, µ, for given β. As a result, it was found that, at β < 0.46, the mixed modes realize the ground state (for small β, this is evident from the comparison of Eqs. (8) and (12)), but in a narrow interval of 0.46 < β < 0.5, the symmetric vortex-antivortex complex, with the upper sign in Eq. (7), produces a slightly smaller value of µ, which is given by Eq. (8), as shown in Fig. 1(c); recall that vortex-antivortex state with the upper sign does not exist at µ > 0.5. The same figure 1(c) shows that the mixed mode does not exists either at β > 0.5 (which is another essential manifestation of the above-mentioned quantum phase transition), hence the only eigenmode which survives at β > 0.5 is the antisymmetric vortex-antivortex complex with the lower sign in Eq. (7). It is worthy to mention here that, as shown below, solely vortex-antivortex modes with opposite signs of their components are found as stable states in the full system of Eqs. (1) and (2), while vortex-antivortex bound states with identical signs of its components cannot be generated by the full system, even if they may realize the ground state of the linearized system at β close to 0.5. The next step is to consider Eqs. (1) and (2) in which nonlinear terms are kept, but the dissipative ones are still ignored: Solving these equations by means of the imaginary-time method makes it possible to produce stationary solutions with a fixed norm, and eventually identify nonlinear ground states, which realize the minimum of the Hamiltonian, (here c.c. stands for the complex-conjugate expression), for the given norm. Figure 2(a) shows the so found profiles of |Ψ + | and |Ψ − | for the ground state at β = 0.48, 0.45, 0.42. In agreement with the above-mentioned transition of the ground state from the mixed mode to the symmetric vortex-antivortex in the linearized system, as the SOC strength changes from β < 0.46 to β > 0.46, Fig. 2 demonstrates that, in the dissipation-free nonlinear system, the same transition takes place as the SOC strength increases from β = 0.45 to β = 0.48.

D. Stable modes in the full nonlinear dissipative system
The systematic numerical analysis demonstrates that the full system of Eqs. (1) and (2) gives rise to stable stationary states only if the diffusion term is present, η > 0 (at η = 0, chaotic solutions appear, which are not shown here in detail). For β not two large, the system with small values of η > 0 generates stable patterns of the mixed-mode type, which thus plays the role of the fundamental state in the full system, see an example in Fig. 3(a).
At larger values of the SOC strength, β, and η > 0, the full system readily produces the fundamental state in the form of stable vortex-antivortex complexes with the opposite sign of the components, which, in the linear limit, correspond to the lower sign (antisymmetric bound state) in Eq. (7). A typical example is presented in Fig. 4, where panel (a) shows real parts of Ψ + and Ψ − (with the opposite signs, as said here), and panel (b) shows |Ψ ± (x)| (the solid line) in the cross section of y = 0. We stress that the numerical solution of the full system, including the nonlinear and dissipative terms, has not produced any stable vortex-antivortex complex with identical signs of its two component, in spite of the fact that such complexes play the role of the ground state in the dissipation-free system close to β = 0.5, as shown above.
The vortex-antivortex states produced by the full system can be quite efficiently approximated in an analytical form. To this end, the following ansatz is adopted for their two components: with free parameters A max and ρ, cf. Eq. (7). First, ignoring the dissipative terms, it is possible to predict values of A max and ρ by minimizing the system's energy (15) at a fixed value of norm N (this time, we do not set N = 1). The results is Further, restoring the dissipative terms ∼ ε, η, and σ, the norm, which appears as a free parameter in Eq. (17), can be predicted by the balance equation (10). The final result is The comparison of the vortex-antivortex modes predicted by Eqs. (16)-(18) with their numerically found counterparts is presented in Figs. 4(b) and (c), showing good accuracy of the analytical approximation. Note that Fig. 4(c) demonstrates that this family of the vortex-antivortex states does not extend to β < 0.4. Results of the systematic numerical analysis are summarized in Fig. 5(a), which displays the existence areas (between the solid and dashed lines) for stable mixed-mode and vortex-antivortex states, in the plane of the most essential control parameters, viz., the SOC strength, β, and diffusion coefficient, η . As mentioned above, the mixed modes and antisymmetric vortex-antivortex complexes tend to be stable, as ground states, at smaller and larger values of β, respectively, with a small bistability region revealed by the figure. Above the dashed lines, the solutions decay to zero, while below the solid lines, the evolution makes them chaotic.

III. STABLE SEMI-VORTEX STATES SUPPORTED BY THE ZEEMAN SPLITTING
The Zeeman splitting, represented by Ω = 0 in Eqs. (1) and (2), can strongly impact properties of nonlinear modes in the present system, as is known to happen in other polariton models, see, e.g., Ref. [42]. We here aim to outline general trends imposed by the Zeeman splitting, without producing full details. As suggested by the analysis of the spin-orbit-coupled two-component atomic condensates [43], the enhancement of the Zeeman splitting tends to replace mixed modes by the semi-vortices, in which, in the present case, one component has zero vorticity, and the other carries vorticity 2. The present system does not give rise to semi-vortices at Ω = 0; however, they appear with the increase of Ω.
The case of very large Ω can be considered by means of an analytical approximation. In this case, the chemical potential contains a large term, −Ω, suggesting one to substitute where the time dependence in Φ ± is assumed to be slow in comparison to e iΩt . Then, using an approximation similar to that developed in Ref. [43], we use Eq. (1) to eliminate Φ + in favor of Φ : Finally, the substitution of this approximation in Eq. (2) reduces it to a single fourth-order equation for the wave function Φ − : In particular, in the case of η = ε = σ = 0, we arrive at a linear or nonlinear Schrödinger equation with a kinetic super-energy term, expressed by the fourth-order operator ∼ β 2 /Ω: It is obvious that the present approximation gives rise to semi-vortices, with a large zero-vorticity component produced by isotropic equation (21), and a small vortex one generated from it by Eq. (20) An example of a stable semi-vortex is presented in Fig. 3(b), which shows profiles of absolute values of the zerovorticity component (the solid curve), |Ψ + |, and its counterpart with vorticity 2 (|Ψ − |, the dashed line). Results of the systematic study of the semi-vortices are summarized in Fig. 5(b), in the parameter plane of (β, Ω). Stable semivortices exist above the boundary shown in the figure, provided that Ω exceeds a finite threshold value, Ω min ≈ 0.16.

IV. SUMMARY
The aim of this work is to identify species of robust 2D localized modes which play the role of fundamental states in two-component (spinor) exciton-polariton condensates, subject to the action of SOC (spin-orbit coupling), represented by the second-order linear differential operator which mixes the two components. The system includes ordinary terms which are common to semiconductor-microcavity models, such as the linear gain and effective diffusion, nonlinear loss, self-repulsive nonlinearity in each component, and the isotropic harmonic-oscillator trapping potential. Starting from the analysis of the linearized dissipation-free system, and then proceeding to its full form, we have found, by means of numerical and analytical methods, that basic states supported by the system are antisymmetric vortex-antivortex complexes and mixed modes. The latter ones combine zero-vorticity and vorticity-carrying terms in each component. The mixed modes and vortex-antivortex complexes tend to realize the fundamental states (similar to the ground states of conservative systems) when SOC is, respectively, weak or strong, with a small region of bistability. The presence of the effective diffusion is a necessary condition for the stability of modes of both types in the full system. We have also found that the addition of the Zeeman splitting tends to replace these modes by stable semi-vortices, with vorticities 0 in the larger component and 2 in the smaller one.