Persistent polarization oscillations in ring-shape polariton condensates

We predict the limit cycle solution for a ring-shape bosonic condensate of exciton-polaritons confined in an optically induced rotating trap. The limit cycle manifests itself with polarization oscillations on a characteristic timescale of tens of picoseconds. The effect arises due to the interplay between orbital motion and the polarization degree of freedom. It is specific to spinor bosonic condensates and would be absent in a scalar case, where a bi-stability of stationary solutions would be observed instead. This work offers a tool of initialisation and control of qubits based on superpositions of polariton condensates characterised by different topologic charges.


I. INTRODUCTION
Strongly coupled exciton-photon systems are known for their pronounced optical nonlinearities, enabling the dynamic control of light.A typical representative of this class of systems is an optical microcavity with one or several semiconductor quantum wells characterised with resonant excitonic transitions [1].In such a structure, if the energies of confined photons and excitons are tuned close to resonance, their interaction leads to the formation of half-light, half-matter hybrid modes referred to as exciton polaritons.The field of polaritonics keeps attracting the enhanced attention of the research community, as it offers a convenient testing platform for various quantum coherent and nonlinear effects [2].A remarkable phenomenon that came to the focus of attention in the new century is the formation of dynamical Bose-Einstein condensates (BECs) of polaritons at exceptionally high temperatures [3,4].Similar to regular BECs formed at the thermodynamic equilibrium, polariton condensates clearly demonstrate such properties as macroscopic coherence and superfluidity, arising from the repulsive polariton-polariton interactions [5][6][7][8][9].The effects of polariton condensation can be employed for the creation of a new class of nanoscale coherent light sources known as polariton lasers [10][11][12].
Recently, it has been suggested that a superposition of polariton vortices with opposite topological charges can serve as a robust and scalable qubit [25][26][27][28].This makes especially important creation of a theoretical framework for the description of the properties of polariton vortex states, including their nonlinear dynamics [20,[29][30][31][32].
Central to this problem is the development of protocols allowing for the control of vortex topological charges [33][34][35][36][37][38][39] It should be noticed that while quantized vortices can exist in spatially uniform planar microcavities, the introduction of a confining potential can stabilize vortices and offer an additional tool to control their properties [40][41][42][43][44][45][46][47].In particular, introduction of chirality into the confining potential facilitates the excitation of polariton vortex states with a precise control over their topological charge [16,17,48,49].The problem of spinning of a polariton condensate by optically induced rotating potentials has been recently studied theoretically and experimentally [50][51][52][53].It has been demonstrated that in the case of a single component (scalar) condensate, in the vicinity of the condensation threshold, the formation of a condensate is governed exclusively by the properties of the fastest growing linear mode."The winner takes all" scenario [53,54], where the fastest growing mode is the first to reach the nonlinear stage at which it depletes the pump and suppresses the growth of the other modes, is realized in this regime.
In the present work, we analyze the mechanisms of mode selection in spinor condensates, where the polarization degree of freedom is accounted for.Polaritons of different polarizations exhibit distinct behaviors under the influence of spin-orbit interaction effects, which link the polarization of polaritons to their ballistic propagation within the microcavity plane.Among such effects are the optical spin Hall effect [55,56], zitterbewegung of polaritons [57,58], topological spin Meissner effect [59] and etc.Of particular interest are the effects of spin-orbit interaction in structures with annular geometry [59][60][61].Spin-orbit interaction occurs under conditions of splitting of polarization states of polaritons, which can be induced by both the inherent optical anisotropy of the structure [62,63] and external influences [64][65][66].In layered structures of optical microcavities, the most signif-icant source of spin-orbit interaction is the splitting of transverse electric (TE) and transverse magnetic (TM) polariton modes [67].
In this paper we show that in certain regimes, "the winner takes all" scenario can be realized in the spinor condensate as well.However, we also show that in the specific range of the rotation speeds, the superposition of two modes characterized by different frequencies and polarizations becomes stable.This gives rise to the limit cycle regime characterized by periodic polarization oscillations.Quite remarkably, these oscillations are not accompanied by the oscillations of the total polariton occupancy.This is demonstrated by direct numerical simulations of the two-dimensional (2D) generalized Gross-Pitaevskii equation describing the dynamics of a drivendissipative polariton condensate.Additionally, we develop a simple perturbation theory based on the coupled modes approach [17,18], which provides a qualitative description of the predicted phenomena.

II. THE FORMALISM
We consider the system consisting of an axially symmetric semiconductor microcavity that is incoherently pumped by two linearly polarized Laguerre-Gaussian optical beams.These beams possess different angular momenta, l 1 and l 2 , ∆l = l 1 − l 2 ̸ = 1, and are slightly detuned in frequency, so that Ω = ω 1 ̸ = ω 2 .The resulting spatial profile of the pump can be represented as P (t, r, θ) = P s (r) + P r (r) cos ∆l(θ + Ωt), (1) where r and θ denote the radial and angular coordinates, and P s (r) and P r (r) describe the radial distributions of the non-rotating and rotating components of the pump.
In further consideration, we threat the rotating pump excited by a superposition of two optical Laguerre-Gaussian beams having the angular momenta l 1 = 1 and l 2 = −1, thus ∆l = 2.The optical pump excites the reservoir of incoherent excitons.Due to the repulsive polariton-polariton interactions, the excitonic reservoir acts as an effective rotating complex trapping potential for the polariton condensate, with its imaginary part describing the stimulated relaxation of the excitons into the polariton state.We also account for an additional real rotationally symmetric trapping potential stemming from the patterning of a microcavity and stabilizing the polariton confinement.A schematic plot of the system is given in Fig. 1(a).
As typically the internal reservoir dynamics is much slower then the relaxation from it into the condensate, the reservoir can be adiabatically eliminated [68] and the system is described by the two coupled generalized Gross-Pitaevskii equations for order parameters Ψ ↑,↓ corresponding to two circular polarizations: The first term on the right-hand side describes the dispersion of polaritons characterized by effective mass m eff .The second term accounts for the presence of an external conservative potential V, created through microstruturing, together with a finite lifetime of polaritons characterized by an energy broadening γ.The third term corresponds to the TE-TM splitting, with normalized splitting magnitude β being proportional to the difference of the effective masses of TE and TM polarized polaritons [69].The fourth term describes the blueshift of a polariton condensate due to repulsive polariton-polariton interactions, with coefficient H describing the interaction between polaritons with same circular polarizations, H the interaction between polaritons with opposite circular polarizations.
The last term describes polariton-reservoir coupling, with its real part corresponding to the reservoir-induced polariton blueshift, and its imaginary part corresponding to the condensate gain stemming from the stimulated relaxation from the reservoir to the condensate.The exciton reservoir is driven by the linear polarized optical pump of intensity P , which excites equally both circular polarization components.The parameter G 1 is the scattering rate from the reservoir to the condensate, and Γ is the relaxation rate of the reservoir excitons.The coefficient G 2 defines the reservoir-induced blue shift of the condensate.The small red shift caused by the interaction of polaritons with excitons of opposite polarization is characterized by the coefficient G2 .
For further analysis, we introduce dimensionless variables, scaling time in units of the inverse polariton dissipation rate, t 0 = γ −1 , t → t/t 0 , the spatial coordinates in units of l 0 = ℏ/ 2γm eff , r → r/l 0 , and the order parameter Ψ ↑,↓ → 2G 1 / Γψ ↑,↓ .Then Eq. ( 2) can be written as following: where V = V/ℏγ is the normalized stationary potential, ϵ = G 1 /G 2 is the ration of the effective gain to the frequency shift caused by the reservoir, g = G2 /G 2 is the ration of the frequency shifts due to interactions of polaritons with the reservoir excitons of the opposite and the same polarizations, p = P G 2 /γΓ is the normalized incoherent pump, h = HΓ/2γG 1 and h = HΓ/2γG 1 are the relative strengths of the self-and cross-polarization polariton-polariton interactions.
In our simulations, we use the following values of the parameters, typical for polariton systems: β = 0.05, ϵ = 0.33, and g = −0.1.We focus on the scenario in which the nonlinearity stemming from the reservoir is dominating and take the coefficients h = 0.018 and h = −0.001.The corresponding values of dimensional parameters in Eq. ( 2) are given in [70].We also noted that in the vicinity of the condensation threshold, variations of the values of these parameters did not qualitatively change the picture.External confining potential was taken in the form: 25 and W V = 1.This corresponds to a ring-shape confinement with the height of the potential, radius, and width of the ring being about 0.28 meV, 9.89 µm, and 4.39 µm, respectively.
As discussed above, the pump is the combination of the static symmetric term p s (r) and the term p r (r) cos[2(θ − Ωt)] rotating with angular velocity Ω.The radial distribution of the static pump is taken in the following form: For our simulations, we take p s0 = 3.45, R s = 0.9 and W s = 0.25.This corresponds to a ring-shape pump of intensity 0.17 ps −1 µm −2 with radius and width of the ring being of about 3.96 µm and 1.1 µm, respectively.In this case, the confinement of polaritons results from the combination of the external conservative potential V and static potential induced by the pump.Selecting a smaller radius of the pump annulus compared to the radius of the external potential enhances the efficiency of the incoherent excitation of polaritons.Such a combination ensures that only condensates with angular indices of either +1 or −1 are formed in the system.
The radial distribution of the rotating pump component p r (r) is taken to be equal to that of the stationary one p s (r) = p r (r).In our simulations, we take the following values of the parameters of the rotating pump: p r0 = 0.3, R r = 0.9, and W r = 0.25.This corresponds to a ring-shape pump of intensity 0.015 ps −1 µm −2 with radius and width of the ring being of about 3.96 µm, and 1.1 µm, respectively.

III. RESULTS OF 2D MODELLING A. Stationary regime
We performed numerical studies of the dynamics of polariton states forming in the considered geometry in the presence of a complex rotating potential.We covered a broad range of the rotation velocities, wherein we observed the formation of the stationary states stably developing from a weak noise.Our results reveal that the density distributions of the polaritons in clockwise and counterclockwise circular polarizations (↑ and ↓) depend on the angular velocity of the potential rotation.For instance, at a relatively small rotation velocity Ω = 0.05, we observed the distinct lobes in the density distribution of ↑ polarization (as depicted in Fig. 2(a)), indicating that in this polarization, the condensate is formed by two counterpropagating waves of comparable amplitudes.The phase distribution of the order parameter, illustrated in Fig. 2(b), clearly demonstrates that the winding number for this polarization of the condensate is equal to 1. Let us note here that in the case of two counter-propagating waves, the topological charge is determined by the wave with the larger amplitude.
In the opposite polarization, ↓, the density distribution similarly exhibits two lobes, albeit notably less pronounced.This discrepancy indicates that in this polarization, the waves with angular momenta m = ±1 forming the state possess rather different amplitudes.It's important to highlight that the phase gradient is directed oppositely in the opposite circular polarization, meaning that polaritons of the ↑ and ↓ polarization flow in opposite directions.Furthermore, it's worth noting that the lobes in both polarizations are oriented along the symmetry axis, passing through the maxima of the rotating potential.
For higher rotation velocities, the distributions of the polarization intensities are different, as can be seen in Fig. 3, which shows the snapshots of the stationary polariton densities and phases for Ω = 0.15.Here the lobes in the ↑ polarization are less pronounced than in the ↓ polarization.Notably, the orientation of the lobes differs as well.Namely, the lobes in the ↑ polarization now align The angular velocity of the potential is Ω = 0.05.The dashed lines show the symmetry axes, with the white line passing through the minima and the gray line through the maxima of the rotating potential.
along the symmetry axis that traverses the minima of the potential.
For discussing the polarization properties of polaritons, it is convenient to introduce the normalized threecomponent Stokes vector ⃗ S = (S 1 , S 2 , S 3 ), which components are determined as follows: It's essential to note that we calculate the Stokes components in the rotating frame x ′ oy ′ , see Fig. 4(e).In this frame, the S 1 component corresponds to the x ′ / y ′ linear polarizations, while in the laboratory frame xoy, it characterizes the tangential and radial polarization components of the condensate.The other components, S 2 and S 3 , are responsible for diagonal/antidiagonal and clockwise/counterclockwise circular polarizations, respectively.
The azimuthal distribution of the densities, as well as the Stokes vector components of the polariton condensates in the stationary state, are depicted in Fig. 4 for the discussed rotation velocities of the potential Ω.It is evident that in the regime corresponding to the lower rotation velocity Ω = 0.05 (panels (a) and (c)), the polaiton field is predominantly radially polarized, with S 1 close to −1.Conversely, in the regime corresponding to the higher rotation velocity Ω = 0.15 (panels (b) and (d)), the polaiton field is primarily tangentially polarized, with S 1 close to +1.

B. Breathing regime
An important observation emerges: within a specific range of the angular velocities of the trapping potential Ω, the system supports only breathing steady state solutions.The typical dynamics of polaritons is illustrated in Fig. 6, where the evolution of the number of polaritons and the angular momentum of the polariton states in different polarizations is presented.As shown in panel (a), for relatively slow rotation, the number of particles grows in both polarizations until the steady state is reached.A similar behavior is observed for angular momenta, see panel (b).However, within the angular velocity range indicated by the pink rectangle in Figs.5(a,b), the dynamics takes on a markedly different character.This becomes evident in Fig. 6(c), where the steady state manifests oscillations in the number of polaritons in both polarizations.A spatial counterpart to this periodic inter-polarization transfer of polaritons has been recently observed in polariton waveguides [71].
The angular momenta in the ↑ and ↓ polarizations oscillate in phase, as depicted in panel (d), thus yielding pronounced oscillations in the total angular momentum of polaritons.Intriguingly, the numbers of polaritons in ↑ and ↓ polarization oscillate in anti-phase, leading to only minor oscillations in the total number of polaritons.
This phenomenon arises when the polarization oscillation period T osc is much shorter compared to the characteristic timescale of evolution of the polariton density T dns , T osc < T dns .The latter is determined by how much the pump intensity exceeds the condensation threshold intensity.In particular, for intensities that are sufficiently close to the threshold, this condition is al- ways satisfied.The total number of the polaritons thus cannot follow the rapid periodic transfer of the polaritons between the polarizations, and, therefore, the total number of the polaritons does exhibit pronounced oscillations.Some estimates are given in Appendix D to explain the discussed effect.
In the oscillatory regime, the temporal dynamics of the Stokes vector becomes quasi-periodic, see Fig. 7.The oscillations of the polarization measured at a fixed point in real space (specifically, at x = 0.8, y = 0) are accompanied by oscillations of intensity.The trajectory of the Stokes vector on the Poincaré sphere demonstrates intricate behavior and covers the entire sphere's surface, indicating that during the evolution, the polariton condensate undergoes oscillations between linear and circular polarization states.
With an increase in the rotation velocity, suppression of the oscillations becomes evident.As seen in Fig. 6(e,f), after some transitional processes, a stationary state (shown in Fig. 3) eventually emerges.Nevertheless, it's worth noting that the transitional stage is notably more intricate compared to that observed at lower angular velocities (sf.panels (a,b) and (e,f)).
In the next section, we develop a perturbation theory that allows us to analyze this phenomenon.

IV. COUPLED MODES APPROXIMATION
In order to comprehend the effects reported in the previous section, we developed a straightforward perturbation theory based on the observation that only the modes possessing angular indices of ±1 can be effectively excited in the system.We also restrict our consideration to the case where the spatial structure of the modes is predominantly defined by the real parts of stationary potentials, supposing that dissipative and nonlinear effects, the rotating potential and the TE-TM splitting can be treated as small perturbations.This allows us to represent the field as a sum of the four modes, each characterized by its polarization and angular momentum.The total field can be represented as follows: where ω 0 is the eigenfrequency of the unperturbed mode, and ρ(r) describes their radial structure of the mode.We normalize the latter as 4π rρ(r)dr = 1 ensuring that the total number of polaritons is equal to To eliminate explicit time dependencies of coefficients, we introduce new variables The dynamic equations can be then written as: The details of derivation of Eqs. ( 6) are given in Appendix A.
The effective losses of the modes γ 0 , the coupling strength between the modes of the same polarization but opposite angular momenta η, the effective TE-TM splitting σ, the dissipative nonlinearity α d , the conservative cubic nonlinearity α as well as the conservative crosspolarization nonlinearity α x coefficients can be derived by fitting the special series of 2D simulations.This procedure gives the following estimates for the coefficient values: γ 0 = 0.006, σ = 0.21, α d = 0.081, α = 0.281, α x = −0.027.The coefficient characterizing the rotating potential is expressed as η = η r + iη i , with η r = 0.052, η i = −0.017.
One should pay special attention to the last terms in the right-hand side of Eqs.(6a) and (6d).These terms elucidate the role of the TE-TM splitting in the polariton evolution in the ring geometry, which consists in the coupling of counter-propagating polariton flows via the polarization degree of freedom [18,52].Namely, it couples counterclockwise propagating (+) flows of ↑ polarized polaritons with clockwise propagating (−) flows of ↓ polarized polaritons (see Fig. 1(b)).
We can now define the angular momenta and the occupancies of the polariton state in ↑ and ↓ polarizations via the amplitudes of the polariton modes as follows: With the notation (7), the quantities M tot and N tot defined in the previous section keep their meaning.The results of the numerical solution of Eqs.It is evident that the results of the coupled mode approximation align well with the results of the full-scale 2D simulations.We have also observed that the fitting accuracy tends to improve at lower pump intensities, particularly when the system is just above the condensation threshold.It's worth noting that as the system approaches the condensation threshold, the time required for the formation of the stationary state becomes considerably elongated, rendering the simulations more timeconsuming.
It is worth mentioning that the transition to the oscillatory regime is also qualitatively reproduced by the coupled mode approximation.However, the range of the angular velocities where the oscillatory regime occurs (gray rectangles in Figs.5(a,b)), does not coincide with the one obtained from 2D simulations.This discrepancy could be attributed to the sensitivity of the dynamics of the system to the values of the parameters as well as to the neglect of certain effects in the coupled mode approach, such as the contributions of non-resonant modes and the frequency dependence of effective losses.
Another noteworthy aspect is that at low angular velocities of the rotating potential, the stationary states obtained from both 2D simulations and the coupled mode approach exhibit the same structural characteristics.However, at higher angular velocities, this congruence no longer holds, and the coupled mode approach predicts a distinct structure of the polariton condensates.A possible explanation is that the stationary state anticipated by the coupled mode approach exists in the 2D model, but it proves to be unstable with quite a relatively small increment.This instability might account for the scenario depicted in Fig. 6(e,f), where an oscillatory state initially forms with the averages quite close to that predicted by the coupled mode approach for the chosen parameters.Howeverm in the coupled mode approach, we observe the decay of the oscillations, whereas in full-scale 2D simulations, transition to a different state takes place.It is acknowledged that further research is required to comprehensively understand this intricate dynamic behavior, which is left for future exploration.Based on our numerical simulations, we propose that in the vicinity of the condensation threshold, the formation of the stationary state follows the "winner takes all" scenario.When the initial conditions involve a randomly distributed low-intensity field, the mode with the fastest growth rate is the first to reach the nonlinear stage.However, both the nonlinear depletion of the pump and the linear increment are small, and consequently, the nonlinearity does affect the structure of the modes.It is worth noting that the rotating pump couples the modes of the unperturbed problem, resulting in the eigenmodes becoming compositions of the eigenmodes of the axially symmetric problem.When we refer to the "structure of the eigenmode", we are describing the relative amplitudes of the waves that compose the eigenmode.
In the scenario under consideration, the fastest growing mode tends to dominate and suppress the other modes before its amplitude reaches the stationary value.Consequently, the structure of the nonlinear stationary state gets inherited from the fastest growing linear mode.This scenario holds well in the scalar case.For certain values of the rotation velocities, this concept also applies to the vector case, as detailed in Appendix B, where we analyze the linear eigenmodes.Further insights into the perturbative analysis of the nonlinear stage of the dynamics can be found in Appendix C, where we derive an approximate expression for the amplitude of the stationary state.
In the vector case, at specific rotation velocities Ω, the fastest growing mode exhibits a structure where a majority of polaritons are in the ↑ (or, for a different sign of Ω, in ↓) polarization.Consequently, the mode with dominant ↑ (↓) polarization cannot effectively suppress the mode where a significant portion of polaritons are in the ↓ (↑) polarization.This phenomenon is illustrated in Fig. 5(b), where the polariton portion in one polarization decreases before transitioning to the oscillation regime.To substantiate this observation, an analysis of the stability of the non-oscillating states is required.This analysis is presented in Appendix D. Additionally, in the same Appendix D, it is demonstrated that the oscillating states can be seen as the coexistence of two sub-states of different frequencies.The structures of the sub-state are similar to the structures of the fastest and second fastest growing modes.

V. CONCLUSION
We have studied the polarization dynamics of polariton condensates in the presence of a rotating potential.Our main attention was focused on low density condensates that are formed when the pump intensity is close to the condensation threshold.We have shown that during the initial linear stage of the condensation, the rotating potential leads to the formation of super-modes, resulting from the hybridization of modes affected by TE-TM splitting with those unaffected by it.The most pronounced hybridization occurs at the resonant angular velocity of the rotating potential, which matches the strength of TE-TM splitting.Due to the dissipative component of the rotating potential, the imaginary parts of these supermodes differ.Interestingly, the mode with the highest frequency also exhibits the highest growth rate.
We have demonstrated that in a non-conservative system, the "winner takes all" scenario can be realized.In this case, during the nonlinear stage of condensation, the fastest growing mode suppresses the other modes.Consequently, the stationary state is predominantly determined by the structure of the fastest growing linear mode, provided that the pump intensity is close to the condensation threshold and the density of the polariton state does not influence the field distributions.
The main result of this study is the finding of an original dynamics of the system characterised by persistent oscillations of the polarization.This limt cycle regime can be achieved even for the pump intensity slightly above the threshold.It occurs for certain angular velocities where the system is in a state where the fastest growing mode predominantly consists of polaritons of a specific circular polarization.Consequently, this mode cannot effectively deplete the gain created by the pump in the opposite polarization.As a result, the remaining pump is sufficient to support the growth of another super-mode.Thus, the final state comprises two coexisting super-modes with different polarizations.
The frequency difference between these modes leads to oscillations in the number of polaritons in clockwise and counterclockwise circular polarizations.The total angular momentum of polaritons undergoes periodic oscillations as well, but the total number of polaritons remains close to constant over time.
The results of direct numerical solutions of the twodimensional coupled generalized Gross-Pitaevskii equations are supported by a simple semi-analytical theory based on coupled mode approximation.
The eigenvalues of L correspond to the eigenfrequencies of the polariton modes, with the imaginary part indicating the decay rate of the mode.While exact solutions of Eq. (B1) can be obtained analytically, they are cumbersome and are not included in the text of this paper.
For further analysis, we rely on the numerically calculated dependencies of the eigenfrequencies on the angular velocity of the rotating potential.These dependencies are shown in Fig. 8.It is evident that there exist four modes characterized by distinct effective decay rates γ eff .These differing decay rates emerge due to the dissipative contribution introduced by the rotating potential, characterized by the coupling constant η.The effective decay rate γ eff encompasses a range of non-conservative processes that influence the mode.These processes include not only losses of polaritons but also the filling of the mode due to the presence of an external pump.Due to the pump, the dissipation rate γ eff can take negative values, which correspond to the temporal growth of the mode.From the physical point of view, this growth signifies the condensation of the incoherent reservoir excitons into the coherent polariton mode.It is worth mentioning that in our model, the mode with the highest frequency experiences the fastest growth.Given our focus on this mode in further analysis, we extract its associated eigenvector ⃗ X = (X ↑+ , X ↑− , X ↓+ , X ↓− ) T to emphasize its characteristics The structure of the stationary polariton state is expected to resemble that of the fastest growing mode when the latter effectively suppresses other modes, and the stationary condensate density remains low enough to alter the field distributions in the polariton state significantly.To validate the similarity between the structures of the fastest growing mode and the stationary polariton state, we conducted a comparison of their normalized angular momenta M ↑,↓,tot and population difference in the orthogonal circularly polarized components ∆N/N tot , as depicted in Figs.5(c,d).The indicated quantities for the fastest growing mode were calculated from Eq. ( 7) after substituting C with X.The dependencies of these quantities on the angular velocity Ω are shown in Fig. 5(c,d).It is evident that both the angular momenta and the population differences of the stationary nonlinear states are quite close to those of the fastest growing linear mode.This robust similarity strongly suggests that the former inherits the structure of the latter.We now extend our investigation by considering the evolution of polaritons within the coupled modes approximation while introducing nonlinear effects through the framework of perturbation theory.We amend Eq. (B1) by incorporating the nonlinear term ⃗ N ( ⃗ C) from the right hand side: The nonlinear term can be derived by the calculation of the nonlinear term and then by the projection on the modes of the axially symmetric conservative problem.
After a simple algebra one can obtain As discussed in the main part of the paper, the rotating potential couples the modes.To get analytical results on the nonlinear dynamics, it is convenient to write the equation analogous to C1 but for the amplitudes of the super-modes -the eigenmodes of the linear problem accounting for all real potentials.The dissipative and nonlinear effects will be treated perturbatively in terms of the amplitudes of the super-modes.Thus, in our treatment, we assume that the linear conservative terms included in the real part of the operator L, namely the TE-TM splitting (σ) and the conservative mode coupling originating from the rotating potential (η r ), significantly outweigh the increments (decrements) of the modes as well as the contributions from all nonlinear terms.
Then we split the right-hand side of Eq. (C1) into two components: the conservative linear part, which provides the fundamental solutions to a linear problem, and the supplementary part, encompassing dissipative and nonlinear effects that we treat as perturbations.Then, in the leading approximation order, the linear part is characterized by the operator including the TE-TM splitting (σ) and conservative coupling due to the rotating potential (η r ).The nonconservative effects are of the next perturbation order and characterized by the operator including the losses of the modes (γ 0 ) and gain from the rotation-induced mode mixing (η i ).In (C3) and (C4), η r = Re(η) and η i = Im(η).
We look for a solution in the following form: where ω k are the eigenvalues of L0 .One should underline that ω and ⃗ Y are purely real.For the sake of convenience, we normalize ⃗ Y such that | ⃗ Y | 2 = 1.We assume that the mode ⃗ Y 1 is the fastest growing and thus dominating mode.
Then we project Eq. (C1) onto the basis of eigenvectors ⃗ Y of the operator L0 = Re( L) and consider the terms L1 ⃗ C = Im( L) ⃗ C and ⃗ N ( ⃗ C) as perturbations.For our purposes, it is sufficient to derive the equation for the amplitude of the fastest growing super-mode, assuming that the other super-modes are negligibly small.The equation reads: where ⃗ N y1 is ⃗ N calculated for ⃗ Y 1 .From this equation, it is easy to find the stationary amplitude of the mode: where The correction to the frequency of the state is found as follows: The total number of polaritons is expressed through the amplitudes a k as N = k |a k | 2 .The perturbatively found dependency N on the angular velocity of the potential Ω is illustrated in Fig. 9(a) with a green dashed line.Clearly, perturbation theory provides a satisfactory approximation for the parameters employed in direct numerical simulations.We have verified that this agreement becomes better for lower values of the effective gain γ eff .
Appendix D: The stability of a single-mode solution and oscillating states The previously found single-mode solution may be unstable, with other modes potentially growing and significantly affecting the solution.To assess its stability, we linearize (C1) and derive an equation for small excitation ⃗ ξ on the background of the stationary solution a 1st ⃗ Y 1 .Excitations with eigenfrequencies substantially detuned from the frequency of the stationary solution, when the detuning is much greater than the characteristic timescales of the dissipative and nonlinear terms, can be described by the following equation: It is important to note that this equation does not account for parametric effects.
Considering M as a perturbation, we can search for a solution in the form ⃗ ξ k = ⃗ Y k e iω k t e iδ k t and derive expressions for δ k : One should note that (D2) is not applicable to excitaions having the same structure as the fastest growing mode ⃗ Y 1 , so this formula is valid only for k ̸ = 1.
The perturbatively obtained dependencies of the increments of the modes γ lp k = Im(δ k ) are shown in Fig. 9(b).The negative value of γ lp k indicates that the corresponding mode k grows in time.One can see that for the mode labeled as "3" in the figure (red curve), there is a range of angular velocities Ω, where it exhibits a positive increment.This implies that the stationary state cannot be defined by the fastest growing linear mode alone, and one can anticipate the oscillatory dynamics.
Note that the range of Ω for the oscillating behavior is shifted compared to what was observed in direct numerical simulations.This discrepancy can be attributed to the relatively large magnitude of the perturbation for the parameters used in direct numerical simulations.We have verified that most of this discrepancy arises from the nonlinear shift of the modes' eigenfrequencies, and that the agreement improves as pump intensity decreases.Now, let's consider whether the oscillatory state can be considered as coexistence of modes with the structures of the fastest growing linear mode and the mode that is not suppressed by the former.According to the perturbation theory, the second mode constituting the oscillating state should be similar to the linear mode labeled as "3" in Fig. 8.
Through direct numerical simulations, we have obtained the temporal spectra of the stationary states for various angular velocities of the potential Ω, as depicted in Fig. 10.In this plot, an interval of Ω exists where the system exhibits two-frequency dynamics, when some polaritons oscillate at a higher frequency δ 1 , while others oscillate at a lower frequency δ 2 .The white curves on this spectrum correspond to the eigenfrequencies of the modes "1" (fastest growing) and "3" shifted down by 0.04.This shift is approximately equal to the shift observed in these modes due to the nonlinearity effect.Notably, these spectral lines closely follow the dependencies of the linear eigenfrequencies on the potential rotation velocity.
In our numerical experiments, we have the flexibility to isolate the first or second spectral lines, enabling us to analyze the polarization and angular momentum of polariton states associated with each of the frequencies.In essence, we represent the stationary solution as ⃗ C = c  time, the dependencies plotted in panel (b) confirm that the polaritons with frequency δ 2 exhibit an order parameter similar to that of the mode that is not suppressed by the fastest growing mode, which is labeled as "3."Finally, let us provide some estimations to clarify the previously discussed effect, which is the limited influence of polariton transfer between polarizations on the total number of polaritons.The polarization oscillations arise due to the interaction between two coexisting modes, and therefore the period of the polarization oscillations T osc is inversely proportional to the frequency difference between these coexisting modes.For the modes involved, this frequency difference is primarily defined by the TE-TM splitting and can be estimated as ∆δ ≈ 2σ.On the other hand, the characteristic timescale for the evolution of the polariton number T dns is inversely proportional to |γ 0 +η i |.This allows us to estimate the ratio of these characteristic times as T dns / T osc ≈ 2|σ|/ |γ 0 + η i |.For the used parameters, this ratio is approximately 20, which explains why our numerical simulations do not show significant oscillations in the total polariton number despite the inter-polarization polariton transfer.

Figure 1 .
Figure 1.(a) Schematic of the excitation of a polariton condensate in a rotating optical potential.The potential is created by two optical beams with angular momenta l1 = 1 and l2 = −1, characterized by different frequencies, Ω = ω1 − ω2.The polariton condensate exists under the balance of the optical pump and decay.(b) Schematic of the coupling of counterpropagating polariton flows via polarization due to the TE-TM splitting.↑ polarized polaritons (blue arrows) propagating in the counterclockwise direction couple to ↓ polarized polaritons (red arrows) propagating in the clockwise direction.

Figure 2 .
Figure 2. The snapshots showing instantaneous profiles of the intensity and the phase distributions of the polariton field in opposite circular polarizations: ↑ for (a),(b) and ↓ for (c),(d).The angular velocity of the potential is Ω = 0.05.The dashed lines show the symmetry axes, with the white line passing through the minima and the gray line through the maxima of the rotating potential.

Figure 3 .
Figure 3.The snapshots showing instantaneous profiles of the intensity and the phase distributions of the polariton field in opposite circular polarizations: ↑ for (a),(b) and ↓ for (c),(d).The angular velocity of the potential is Ω = 0.15.The dashed lines show the symmetry axes, with the white line passing through the minima and the gray line through the maxima of the rotating potential.

Figure 4 .
Figure 4. (a) The angular distribution of the total density of the stationary state of the polariton condensate |ψ| 2 (black curves) as well as the condensate polarization components |ψ ↑ | 2 (red curves) and |ψ ↓ | 2 (blue curves) for the rotation angular velocities Ω = 0.05 (a) and Ω = 0.15 (b).(c), (d) The angular dependencies of the normalized Stokes parameters for the corresponding angular velocities.The dependencies are plotted for r = 0.8.The Stokes parameters are calculated in the rotating reference frame x ′ oy ′ , see panel (e).(e) The spatial distribution of the rotating potential.xoy is the laboratory reference frame, wile x ′ oy ′ is the rotating reference frame.

Figure 5 .
Figure 5.The dependencies of the normalized angular momenta M ↑,↓,tot (a),(c) and the normalized polarization occupation difference ∆N/Ntot = (N ↑ − N ↓ )/(N ↑ + N ↓ ) (b),(d) on the potential angular velocity Ω.The open circles in (a) and (b)show the quantities obtained from 2D numerical simulations.The solid lines in all panels correspond to numerical simulations based on nonlinear equations from the coupled mode theory.The dashed lines in (c) and (d) are for the fastest growing linear mode calculated from the coupled mode theory.Pink and gray rectangles indicate the ranges of the rotation velocity Ω, predicted from 2D numerical simulations (pink) and from the coupled mode approach (gray), where the steady state exhibits oscillatory behavior.
Figs. 5(a) and 5(a), respectively.The quantity ∆N/N tot actually represents the degree of circular polarization of the polariton condensate.It is evident that the results of the coupled mode approximation align well with the results of the full-scale 2D simulations.We have also observed that the fitting accuracy tends to improve at lower pump intensities, particularly when the system is just above the condensation threshold.It's worth noting that as the system approaches the condensation threshold, the time required for the formation of the stationary state becomes considerably elongated, rendering the simulations more timeconsuming.It is worth mentioning that the transition to the oscillatory regime is also qualitatively reproduced by the coupled mode approximation.However, the range of the angular velocities where the oscillatory regime occurs (gray rectangles in Figs.5(a,b)), does not coincide with the one obtained from 2D simulations.This discrepancy could be attributed to the sensitivity of the dynamics of the system to the values of the parameters as well as to the neglect of certain effects in the coupled mode approach, such as the contributions of non-resonant modes and the frequency dependence of effective losses.Another noteworthy aspect is that at low angular velocities of the rotating potential, the stationary states obtained from both 2D simulations and the coupled mode approach exhibit the same structural characteristics.However, at higher angular velocities, this congruence no longer holds, and the coupled mode approach predicts a distinct structure of the polariton condensates.A possible explanation is that the stationary state anticipated by the coupled mode approach exists in the 2D model, but it proves to be unstable with quite a relatively small increment.This instability might account for the scenario depicted in Fig.6(e,f), where an oscillatory state initially forms with the averages quite close to that predicted by the coupled mode approach for the

Figure 7 .
Figure 7. Evolution in time of the polariton densities (a) and the Stokes parameters (b) in the oscillatory regime for the polariton field at x = 0.8, y = 0. (c) The trajectory of the normalized Stokes vector on the Poincaré sphere.

Figure 8 .
Figure 8.The real (a) and imaginary (b) parts of the frequencies of the polariton eigenmodes as functions of the angular velocity of the potential Ω.
Appendix C: Perturbative treatment of the nonlinear problem

Figure 9 .
Figure 9. (a) The dependencies of the polariton number in the stationary states, calculated within the coupled modes approximation, on the rotation velocity Ω of the potential.The black solid line corresponds to the direct numerical simulations of Eq. (6), while the dashed green line is for the results of the perturbation theory treatment.The shaded rectangle marks the range of Ω where single-frequency solutions of Eqs.(6) do not exist.(b) The effective increments of the linear excitations on the backgroung of the stationary states calculated through the perturbation theory.The colors of the lines correspond to those in Fig.8.The shaded rectangle shows the range of Ω where one of the modes (red) has a positive increment.

1 ⃗ Z 1 eFigure 10 .
Figure 10.The temporal spectrum of the stationary state as a function of the angular velocity of the potential Ω.The simulations were conducted within the coupled mode approximation.The white curves show the red-shifted dependencies of the real parts of the eigenfrequencies of the two linear modes, including the fastest growing mode shown by the black curve and the state indicated by the red curve in Fig. 8(a).The eigenfrequencies are shifted downward by 0.04.

Figure 11 .
Figure 11.The structure of the components corresponding to the upper (a) and the lower (b) spectral lines shown in Fig. 10.The open circles correspond to the data extracted from the numerical simulations within the coupled modes model, while the solid lines correspond to the quantities calculated for the eigenmodes shown by black (for panel (a) ) and red (for panel (b) ) colors in Fig. 8.The field structure is characterized by the rations Ni n Nn = |Ci| 2 | ⃗ C| 2 shown as functions of Ω.