Critical dynamics and tree-like spatiotemporal patterns in exciton-polaritoncondensates
Nataliya Bobrovska, Andrzej Opala, Paweł Miętki, Michał Kulczykowski, Piotr Szymczak, Michiel Wouters, Michał Matuszewski
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] A p r Critical dynamics and tree-like spatiotemporal patterns in exciton-polaritoncondensates
Nataliya Bobrovska, Andrzej Opala, Paweł Miętki, Michał Kulczykowski, Piotr Szymczak, Michiel Wouters, and Michał Matuszewski Institute of Physics Polish Academy of Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland Institute of Theoretical Physics, Faculty of Physics,University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland TQC, Universiteit Antwerpen, Universiteitsplein 1, B-2610 Antwerpen, Belgium
We study nonresonantly pumped exciton-polariton system in the vicinity of the dynamical in-stability threshold. We find that the system exhibits unique and rich dynamics, which leads tospatiotemporal pattern formation. The patterns have a tree-like structure, and are reminiscentof structures that appear in a variety of soft matter systems. Within the approximation of slowand fast time scales, we show that the polariton model exhibits self-replication point in analogy toreaction-diffusion systems.
I. INTRODUCTION
Semiconductor exciton-polaritons are quantum quasi-particles that exist in structures where strong light-matter coupling overcomes decoherence . Properties ofmicrocavity polaritons, which combine the extremelylow effective mass of confined photons with strong in-teractions of excitons, makes them an ideal candidatefor studying quantum fluids of light . Rapid progressin studies of these systems has led to observations ofremarkable phenomena, such as nonequilibrium Bose-Einstein condensation , quantum vortices , superflu-idity , and Berezinskii-Kosterlitz-Thouless phase transi-tion .Several recent experiments provided evidence of dy-namical instability in exciton-polariton condensates inthe case of nonresonant pumping . This instabil-ity is an inherent property of the open-dissipative GrossPitaevskii model, widely used for describing the dynam-ics of these systems . Signatures of instability were ob-served both in the case of organic microcavities , as wellas inorganic GaAs microcavities pumped continuously and with ultrashort pulses .Despite these experimental observations, most studiesof polariton fluids to date have focused on the stableregime of condensation. In particular, properties of thesystem close to instability threshold have not been a topicof a detailed study. This is of practical importance, sinceboth stable and unstable regimes of condensation havebeen demonstrated experimentally . It was pointedout that this regime can be characterized by interestingchaotic dynamics with unusual momentum distributionof fluctuations . Note that chaotic evolution has beenrecently predicted to occur also in a polariton model withresonant plane wave driving .Spatial pattern formation in polariton systems hasbeen studied in a number of different configurationsboth theoretically and experimentally . In this work,we investigate dynamical behavior close to the instabil-ity threshold in detail, and predict spatiotemporal pat-tern formation. We find that the dynamics results in tree-like structures in space-time coordinates, which ex-hibit branching, or self-replication. The behavior of thesystem becomes very similar to that occurring in cer-tain soft matter systems, including combustion , bacte-rial growth , chemical reactions , wetting films , orself-replicating pattern formation in general diffusion-reaction models . We describe the physical mechanismresponsible for branching, resulting from phase separa-tion into condensed and uncondensed regions. In analogyto reaction-diffusion systems, the existence of two timescales, corresponding to slow evolution and fast splittingdynamics, allows to understand the occurrence of self-replication and determine the threshold for its occurence.As a result, we find that polariton systems in the crit-ical regime display rich dynamics that is very differentfrom superfluid behavior observed in the stable regime.We also provide an analogy to extensively studied soft-matter systems.We discuss the experimental observation of splittingdynamics. We point out that while direct detectionof branching would be difficult in experiment due tothe chaotic nature of the process and the picosecondtime scales involved, it is possible to observe signaturesof branching in second-order spatiotemporal correlationfunction. This method allows to perform time-averagedexperiment in which many branching events occurring ina condensate over a long acquisition time contribute to anontrivial pattern of spatiotemporal correlations, whichcan be considered a smoking gun of branching dynamics. II. MODEL
We model evolution of an exciton-polariton conden-sate using the open-dissipative Gross-Pitaevskii equation(ODGPE) for the wavefunction ψ , coupled to the rateequation for the density of exciton reservoir, n R . In ourwork, we will focus mainly on the one dimensional case,when the condensate is trapped in a 1D microwire .Results in the two dimensional case are briefly discussedin Appendix B. The 1D evolution equations read i ~ ∂ψ∂t = − ~ D m ∗ ∂ ψ∂x + g C | ψ | ψ + g R n R ψ + i ~ Rn R − γ C ) ψ,∂n R ∂t = P ( x ) − ( γ R + R | ψ | ) n R , (1)where P ( x ) is the exciton creation rate determined bythe pumping profile, m ∗ is the effective mass of lower po-laritons, D = 1 − iA where A is a dimensionless constantaccounting for kinetic energy relaxation, γ C = τ − C and γ R = τ − R are the polariton and exciton loss rates relatedto their lifetimes τ C,R , and ( R, g i ) = ( R , g i ) / √ πd are the rates of stimulated scattering into the conden-sate and the interaction coefficients, rescaled in the one-dimensional case , where d is of the order of the mi-crowire width.In a model without noise, a nonzero homogeneous sta-tionary solution of the above model can be found in theform ψ ( x, t ) = ψ e − iµ t , n R ( x, t ) = n R . This solutionexists above threshold pumping P > P th = γ C γ R /R andis given by | ψ | = ( P/γ C ) − ( γ R /R ) , n R = γ C /R , and µ = g C | ψ | + g R n R . This homogeneous solution be-comes dynamically unstable (via Benjamin-Feir instabil-ity) in a certain parameter range, as predicted andrecently observed experimentally . The criterion forlinear stability in the case A = 0 was derived in PP th > g R g C γ C γ R . (2)In the case when A = 0 linear stability can be determinedby solving Bogoliubov eigenvalue problem numerically.We note that the dynamics predicted in this paper ap-pear to be quite general, and not limited to the modeldescribed above. As we demonstrate in Appendix A, thephysics described occurs as well in a model of polari-tons in a semimagnetic microcavity. This model doesnot include a reservoir explicitly, and the second degreeof freedom is provided by the magnetization of magneticions. III. RESULTS
Figures 1(a)-(c) present examples of numerical dynam-ics of the ODGPE model (1) (a) in the stable regime,(b) in the critical-unstable regime close to the stabilitythreshold of Eq. (2), and (c) in the deep unstable regime.We assume a small white noise in the polariton and reser-voir fields at t = 0 , and a constant homogeneous pump-ing P > P th for t > . In Figure 1(a) typical behaviorexpected for the stable regime is visible, where initialcondensate density fluctuation decays over time. In Fig-ure 1(c), an apparently random pattern of high densitypeaks is formed, as could be also expected in the unsta-ble case. On the other hand, in the intermediate case Fig. 1(b) the instability apparently leads to pattern for-mation and spatiotemporal chaos, which takes the formof tree-like branching of domains which are character-ized by low condensate density. We verified that suchpatterns appear in a relatively wide region of parameterspace in the vicinity of the critical threshold. However,the estimation of exact limits of this region is a non-trivial task which will be postponed for a future study.The tree-like patterns are reminiscent of those occurringin certain soft-matter systems . Below we draw ananalogy between diffusion-reaction systems and criticaldynamics of the polariton model. We note that some-what similar patterns were recently predicted to occur ina complex Ginzburg-Landau equation (CGLE) polaritonmodel, incorporating a carefully engineered complex pe-riodic potential . However this regime appears to resem-ble spatiotemporal intermittency regime of the CGLE rather than dynamics of diffusion-reaction systems.To investigate the dynamics of branching in more de-tail, we plot the evolution of condensate density andphase, together with the reservoir density for a singlebranching “event” in Fig. 2. The correspondence betweenthe regions of low condensate density and high reservoirdensity is a signature of phase separation, resulting fromthe repulsive polariton-reservoir interaction term g R inEqs. (1). Phase separation is the driving force of dynam-ical instability in a polariton system . Here, it leads tothe formation of well defined regions of high condensatedensity, separated from regions of high reservoir density,and the formation of separate branches visible in Fig. 2.At the same time, it does not lead to a complete decayof the condensate into small lumps, as in the deep unsta-ble regime of Fig. 1(c), since rather wide regions of al-most homogeneous condensate can still be distinguishedbetween the branches. The existence of such two qual-itatively different “phases” of low and high condensatedensity, corresponding to the branches and the regionsbetween them, can be justified by the existence of twostationary homogeneous solutions of Eqs. (1) ( a ) | ψ | = 0 , n R = Pγ R , (3) ( b ) | ψ | = Pγ C − γ R R , n R = γ C R , i. e. the zero solution and the nonzero stationary solu-tion. While both these (spatially infinite) solutions arenot stable in the unstable regime of condensation, the dy-namics of the system appears to locally follow the formof either (a) or (b). This is confirmed by the magnitudeof condensate and reservoir density in the branches andbetween them, which are close to values given by (a) and(b), respectively.The corresponding phase φ of the condensate wave-function ψ = | ψ | e iφ is shown in Fig. 2(c). Notice thatthe phase gradient in the time direction is different onthe left and right hand side of the branch, as followsfrom the different frequency of π rotations of the phasealong the time axis. This evidences the lack of phase co- FIG. 1: Spatiotemporal pattern formation. Upper panels show evolution of condensate density in (a) the stable regime, (b)the critical regime, close to instability threshold, and (c) the deep unstable regime, as indicated in the stability diagram below.Panel (b) reveals spontaneous formation of the tree-like spatiotemporal patterns. The grey star in the phase diagram belowcorresponds to the case where phase turbulence is observed . Purple stars correspond to cases with clear tree-like branchingevolution. Parameters are m ∗ = 3 . × − m e , τ R = 1000 ps, τ ( a ) C = 76 . ps, τ ( b ) C = 62 . ps, τ ( c ) C = 9 . ps, d = 4 µm , g DC = 1 . µ eV µm , g DR = 4 g DC , R D = 4 . × − µmps , A = 0 . . herence between the condensate regions on the two sides.In other words, the condensates which exist between thebranches form uncorrelated condensate islands with nomutual phase coherence, but with coherence within eachcondensate. The branches, on the other hand, are regionswhere there is almost no condensate density and no phasecoherence, which is visible as multiple phase discontinu-ities (spatiotemporal vortices) appearing in Fig. 2(c).The above observations, together with phase gradientsin x direction shown in Fig. 2, allow for the understand-ing of physical mechanism of branching. Density currentof polaritons can be calculated from the standard formula j = − i ~ / m ∗ ( ψ ∗ ∂ψ/∂x − c . c . ) , and is plotted in Fig. 2(d).A single branch before splitting is characterized by fluxof polaritons from inside the branch to the outside re-gions, as shown schematically in Fig. 3 (left). This resultsfrom the repulsive potential g R n R in Eq. (1), created bythe increased reservoir density in the (a) phase insidethe branch. Indeed, above threshold P > P th reservoirdensity is always higher in phase (a) than in phase (b).In the stable regime, this repulsive potential is screenedby the lower condensate density, which acts through the condensate self-interaction term g C | ψ | . However, as weenter the unstable regime, the reservoir-induced repul-sive potential begins to dominate, and leads to outflowof condensate density from the regions of increased reser-voir density, resulting in phase separation.The outside directed flow of polaritons from inside thebranch results in gradual increase of the spatial extentof the branch, as shown in the middle panels of Fig. 3,which is also visible as widening of the branch in timein Fig. 2(a). However, the spatial extent cannot increaseindefinitely, since the (a) phase inside the branch is nota stable state. When the branch becomes wide enough,dynamical instability sets in, leading to splitting of thebranch into two. The stability of the branch below acertain spatial extent of the branch and instability abovethis extent is a crucial property which makes the tree-likedynamics possible. When the branch splits, it developsa small high condensate density area in its center, whichgrows quickly thanks to the spontaneous scattering fromthe reservoir to the condensate. This is possible as thereservoir density is locally high, and the outflow of po-laritons is suppressed locally thanks to the flattening of FIG. 2: Example of (a) condensate density, (b) reservoir den-sity, (c) condensate wavefunction phase, and (d) condensatedensity current. The above figures correspond to a singlebranching event, selected from Fig. 1(b) (marked with a whitedashed box).FIG. 3: Schematic illustration of the physical mechanismleading to branch splitting. (a) The repulsive potential gen-erated by the reservoir-dominated branch leads to expulsionof polaritons from the branch, as well as growth of its dimen-sions. (b) When the branch becomes wide enough, the regionin its center with a flat section of the potential becomes a seedfor a new condensate island. (c) Condensate density quicklygrows, leading to separation of the two new branches. FIG. 4: Self-replication transition occurring when the the fi-nite system size is increased. Due to the periodic bound-ary conditions, this corresponds to increasing the distancebetween neighboring branches. Solid line in panel (a) showsthe position of a minimum (or two minima) of condensatedensity in a steady state, calculated in a box of size L . Thiscorresponds to a chain of equally spaced branches separatedby a distance L . Two representative states shown in panels(b) and (c) correspond to dashed lines in panel (a). The split-ting occurs as the distance to neighboring branches becomeslarger than L threshold ≈ . µ m. Parameters are g DC = 0 . µ eV µ m, g DR = 2 g DC , τ C = 4 ps, τ R = 3 . ps, R D = 0 . µmps d=2 µ m, P/P th = 1 . , A = 0 . . the effective potential as shown in Fig. 3 (middle). Thefast growth of condensate density leads to the formationof two separate branches as depicted in Fig. 3 (right).To describe the physics of splitting more quantita-tively, we employ the time scale separation method, in-troduced in the study of dynamics of self-replicating pat-terns in diffusion-reaction systems . This approach isbased on the assumption that the evolution occurs ontwo different time scales. The slow movement of branchesis occasionally interrupted by fast dynamics of splitting,or self-replication. Within this approach, the solutionsin the slow phase of motion can be found approximatelyby assuming a steady state which consists of a chain ofidentical branches or a single branch within a finite boxwith periodic boundary conditions . The threshold ofsplitting can be determined from stability properties ofthis periodic solution. Such an approximation, althoughclearly not adequate to exactly describe the dynamicsof non-periodic arrangement of branches as in Fig. 1(b),allows to gain insight into the main mechanism drivingthe branching dynamics and determine the approximatethreshold.As obtaining an exact analytical solution is not viablein our nonlinear system, we employ numerical methodbased on the evolution of Eq. (1) in a box of length FIG. 5: Second order correlation function g (2) ( d, τ ) is depicted for the stable (a), critical (b), and unstable (c) cases of Fig. 1.The critical case (b) with branching density patterns is characterized by nontrivial spatiotemporal correlations which cannotbe factorized into independent spatial and temporal parts. The characteristic “horn” features are signatures of branching indensity evolution from Fig. 1(b). In panel (d), cross-sections of g (2) ( d, τ ) in the case (b) are plotted for τ = 0 , an ps.Note that g (2) (0 , is approximately equal to unity in (a) and equal to two in (c), which correspond to a coherent state and aclassical random state (or thermal state), respectively. L with periodic boundary conditions. After sufficientlylong time of evolution, we obtain a stationary stable so-lution. To minimize transient effects and avoid possibleeffects of multistability, we perform the simulations adi-abatically, by feeding the result of one simulation as astarting point of another, with a slightly modified extentof the box L . This allows to follow one stable branch ofsolutions, and by changing L in both directions we candetect the possible effects of bistability. We show theresults of our investigation in Fig. 4. In panel (a) therange of investigated box sizes is shown on the verticalaxis, with the solid line showing the positions of eithera single minimum or two minima of the solution. Thesplitting of the minimum into two occurs when the boxsize is equal to about L threshold ≈ . µ m. The examplesof solutions with a single and two minima are shown inFigs. 4(b,c). These solutions resemble closely the densityprofiles obtained previously in a large system. The ob-tained threshold size of a branch before self-replication L threshold is also in good agreement with typical spatialscales on which branching occurs in full simulations. Atthe same time, we did not observe any region in whichthe two kinds of solutions shown in Figs. 4(b,c) would bestable for the same L . IV. DETECTION OF BRANCHING VIACORRELATIONS
Direct observation of branching shown in Fig. 1(b)would be a challenging task due to the short (picosec-ond) time scale of the dynamics. Although streak cam-eras can be used to observe polariton dynamics on suchtime scales, they usually require averaging over many re-peated realizations of the experiment or over a relativelylong acquisition time. Such methods would not provideevidence of branching due to the chaotic character of theprocess, in which patterns are expected to vary rapidlyand from shot to shot. We propose to circumvent thisproblem by measuring second-order spatiotemporal cor-relations instead of emission intensity. The reasoningbehind such approach is that even if branching occursat random positions and times, we can still recover itscharacteristic features in integrated correlation functions,since all branching events will contribute to it in a similarway. Second-order correlation function is defined as g (2) ( d, τ ) = R | ψ ( x, t ) | | ψ ( x + d, t + τ ) | dxdt ( R | ψ ( x, t ) | dxdt ) (4)where the spatial integral is taken over the size of thesystem. Time integration starts from the instant whenthe system achieves a quasi-stationary distribution, inwhich there are strong fluctuations, but observables havereached a steady state in a statistical sense. In practice,such state is established after several hundred picosec-onds of evolution, when average density saturates.In Figure 5 we visualize correlation functions corre-sponding to the three cases from Fig. 1. Clearly, sta-ble, critical and deep unstable cases are characterized byqualitatively different correlation functions. The charac-teristic horn-like shape of g (2) in Fig. 5(b) is an indica-tion of branching occurring in Fig. 1(b). Note that thehorns are directed both in positive and negative time di-rection, since g (2) ( d, τ ) as defined above is a time- andspace-symmetric function in the limit of infinite integra-tion time. It is important to note that only the criticalcase Fig. 5(b) is characterized by nontrivial spatiotem-poral correlations. In both Fig. 5(a) and (c) correlationscan be approximately factorized into spatial and tempo-ral functions, i. e. g (2) ( d, τ ) ≈ g (2) x ( d ) g (2) t ( τ ) , while suchfactorization is not possible in the case of Fig. 5(b). Thisis clearly shown in Fig. Fig. 5(d), where cross-sectionsof correlation function at three different values of τ areshown. On the other hand, we note that in the unstablecase of Fig. 5(c) temporal correlations g (2) t ( τ ) also havea nontrivial (non-Gaussian) character. V. CONCLUSIONS
In conclusion, we demonstrated that nonresonantlypumped exciton-polariton condensates at the thresholdof instability possess unique and rich dynamics, remi-niscent of self-replicating patterns encountered in manysoft-matter systems. We believe that these results pro-vide an interesting link between quantum coherent wavesystems and soft matter diffusion-reaction systems, whichmay stimulate further interaction between these areas ofphysics.
Acknowledgments
We thank Marzena Szymańska and Sebastian Diehlfor fruitful discussions. Support from National Sci-ence Centre, Poland Grants 2015/17/B/ST3/02273 and2016/22/E/ST3/00045 is acknowledged.
Appendix A: Branching in diluted magneticsemiconductor model
We discuss the generality of the observed effects.We find that branching appears not only in the open-dissipative Gross-Pitaevskii model with a reservoir, butalso in a model of semimagnetic exciton-polaritons inwhich the reservoir is absent. The role of the reservoir is played by the collective magnetization of manganese ionscoupled to the condensate.Recently, experimental investigations of semimag-netic microcavities (Cd − x Mn x Te) were performed in which quantum wells are doped with magnetic ions.In these cavities, phenomena such as giant Zeeman split-ting and polariton lasing were observed . Magnetiza-tion of a diluted magnetic semiconductor is given by theBrillouin function B J h M ( x, t ) i = n M g M µ B J B J (cid:18) g M µ B JB eff k B T (cid:19) , (A1)where n M is the magnetic ion concentration, g M is theirg-factor, J is the manganese total angular momentumequal to / , µ B is the Bohr magneton, k B is the Boltz-mann’s constant, T is the temperature of manganese ions. B eff = ( λ/ | ψ | is the effective magnetic field resultingfrom the presence of exciton-polariton condensate , withthe strength of the ion-polariton coupling is denoted with λ . FIG. 6: Example of branching evolution in the model ofsemimagnetic exciton-polaritons in which the second degreeof freedom is due to the collective magnetization of man-ganese ions rather than the reservoir. Parameters are τ M = . × − s , g = . × − meV m , n M = . × m − , B = , T = . , P − γ L = . × − meV , γ NL = . × − meV m , m ∗ = − m E , Rabi splitting Ω R = 5 meV. Within the description in terms of complex GinzburgLandau equation, which can be obtained from the fullopen-dissipative Gross-Pitaevskii model within the adi-abatic approximation , there is an additional effectivepotential caused by the ion-exciton interaction . i ~ ∂ψ∂t = − ~ m ∗ ∂ ψ∂x + g | ψ | ψ + iP ψ −− i γ L ψ − iγ NL | ψ | ψ − λM ψ, (A2)where the interaction between the polaritons with g , ex-ternal pumping with P and losses (linear and non-linear)with γ L and γ NL . We assume that circular pumping ishomogeneous and the condensate remains circularly po-larized, however the ion polarization is free to evolve.Moreover, we introduce the spin relaxation time ( τ M ) formagnetic ions. Then, the polariton evolution equationcouples to the equation for manganese magnetization re-laxation ∂M ( x, t ) ∂t = h M ( x, t ) i − M ( x, t ) τ M (A3)We found that within this model, there exist a large re-gion in parameter space in which tree-like branching oc-curs, and an example is shown in Fig. 6. Appendix B: Two-dimensional case
To investigate whether dimensionality is an importantfactor in the occurrence of branching, we perform a seriesof numerical simulations in the two-dimensional exten-sion of the model (1) i ~ ∂ψ∂t = − ~ D m ∗ ∇ ψ + g C | ψ | ψ + g R n R ψ ++ i ~ Rn R − γ C ) ψ,∂n R ∂t = P ( x ) − ( γ R + R | ψ | ) n R , (B1)where ∇ = ∂ /∂x + ∂ /∂y , and initial conditions forthe fields ψ and n R are the same as before. We foundthat for a similar range of parameters as in the 1D case,one can observe branching solutions as shown in Fig. 7.In panel (a), we show the density of the polariton con-densate at a given time t final . The dynamics of branchingis visible in panel (b), where a cross-section for y = 0 isshown, demonstrating the formation of spatiotemporalpatterns similar as in previous sections. The parametersof the simulation are given in the Figure caption. Thelower quality of figures is due to the increased numer-ical mesh spacing in the 2D case, which was necessarybecause of the limited computational resources. FIG. 7: Example of branching evolution in the two-dimensional model. (a) Spatial pattern of the density of thecondensate at t final = 10 ps. In panel (b), a cross-sectionof the density evolution at y = 0 is shown, demonstratingbranching patterns emerging from the initial state. Param-eters in physical units are m ∗ = 5 . × − m e , τ R = 10 ps, τ C = 9 ps, g C = 3 . µ eV µ m , g R = 6 . µ eV µ m , R = 5 . × − µm ps , L = 204 µ m, A = 0 . A. Kavokin, J. J. Baumberg, G. Malpuech, and F. P.Laussy,
Microcavities (Oxford University Press, 2007). I. Carusotto and C. Ciuti,Rev. Mod. Phys. , 299 (2013). H. Deng, G. Weihs, C. Santori, J. Bloch, and Y. Ya-mamoto, Science , 199 (2002). J. Kasprzak, M. Richard, S. Kundermann, A. Baas,P. Jeambrun, J. M. J. Keeling, F. M. Marchetti,M. H. Szymańska, R. André, J. L. Staehli, V. Savona,P. B. Littlewood, B. Deveaud, and L. S. Dang,Nature , 409 (2006). K. G. Lagoudakis, M. Wouters, M. Richard, A. Baas,I. Carusotto, R. André, L. S. Dang, and B. Deveaud-Plédran, Nat. Phys. , 706 (2008). D. Sanvitto, F. M. Marchetti, M. H. Szymanska, G. Tosi,M. Baudisch, F. P. Laussy, D. N. Krizhanovskii, M. S.Skolnick, L. Marrucci, A. Lemaître, J. Bloch, C. Tejedor,and L. Viña, Nature Physics , 527 EP (2010), article. L. Dominici, R. Carretero-González, A. Gian-frate, J. Cuevas-Maraver, A. S. Rodrigues, D. J.Frantzeskakis, G. Lerario, D. Ballarini, M. De Giorgi,G. Gigli, P. G. Kevrekidis, and D. Sanvitto,Nature Communications , 1467 (2018). G. Roumpos, M. D. Fraser, A. Löffler, S. Höfling,A. Forchel, and Y. Yamamoto, Nat. Phys. , 129 (2010). A. Amo, J. Lefrére, S. Pigeon, C. Adrados, C. Ciuti,I. Carusotto, R. Houdré, E. Giacobino, and A. Bramati,Nature Physics , 805 (2009). G. Dagvadorj, J. M. Fellows, S. Matyjaśkiewicz, F. M.Marchetti, I. Carusotto, and M. H. Szymańska,Phys. Rev. X , 041028 (2015). D. Caputo, D. Ballarini, G. Dagvadorj, C. Sánchez-Muñoz,M. De Giorgi, L. Dominici, K. West, L. N. Pfeiffer,G. Gigli, F. P. Laussy, M. H. Szymanska, and D. San-vitto, Nature Materials , 145 EP (2017), article. N. Bobrovska, M. Matuszewski, K. S. Daskalakis, S. A.Maier, and S. Kéna-Cohen, ACS Photonics , 111 (2018). E. Estrecho, T. Gao, N. Bobrovska, M. D. Fraser, M. Ste-ger, L. Pfeiffer, K. West, T. C. H. Liew, M. Matuszewski,D. W. Snoke, A. G. Truscott, and E. A. Ostrovskaya,Nature Communications , 2944 (2018). F. Baboux, D. D. Bernardis, V. Goblot, V. N. Glad-ilin, C. Gomez, E. Galopin, L. L. Gratiet, A. Lemaître,I. Sagnes, I. Carusotto, M. Wouters, A. Amo, and J. Bloch,Optica , 1163 (2018). M. Wouters and I. Carusotto,Phys. Rev. Lett. , 140402 (2007). N. Bobrovska and M. Matuszewski,Phys. Rev. B , 035311 (2015). S. S. Gavrilov, Phys. Rev. B , 195310 (2016). S. S. Gavrilov, Phys. Rev. Lett. , 033901 (2018). P. Cristofolini, A. Dreismann, G. Christmann,G. Franchetti, N. G. Berloff, P. Tsotsis, Z. Hat-zopoulos, P. G. Savvidis, and J. J. Baumberg,Phys. Rev. Lett. , 186403 (2013). G. Tosi, G. Christmann, N. G. Berloff, P. Tsotsis, T. Gao,Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg,Nature Communications , 1243 EP (2012), article. J. Keeling and N. G. Berloff, Phys. Rev. Lett. , 250401(2008). M. O. Borgh, J. Keeling, and N. G. Berloff, Phys. Rev. B , 235302 (2010). O. A. Egorov, A. Werner, T. C. H. Liew, E. A. Ostro-vskaya, and F. Lederer, Phys. Rev. B , 235302 (2014). F. Manni, K. G. Lagoudakis, T. C. H.Liew, R. André, and B. Deveaud-Plédran,Phys. Rev. Lett. , 106401 (2011). A. Gianfrate, L. Dominici, O. Voronych,M. Matuszewski, M. Stobinska, D. Ballarini,M. De Giorgi, G. Gigli, and D. Sanvitto,Light: Science & Applications , 17119 EP (2018),article. A. Werner, O. A. Egorov, and F. Lederer,Phys. Rev. B , 245307 (2014). N. H. Kwong, C. Y. Tsang, S. M. H. Luk,Y. C. Tse, C. K. P. Chan, P. Lewandowski,P. T. Leung, S. Schumacher, and R. Binder,Physica Scripta , 034006 (2017). G. Díaz-Camacho, C. Tejedor, and F. M. Marchetti,Phys. Rev. B , 245309 (2018). C. E. Whittaker, B. Dzurnak, O. A. Egorov, G. Buonaiuto,P. M. Walker, E. Cancellieri, D. M. Whittaker, E. Clarke,S. S. Gavrilov, M. S. Skolnick, and D. N. Krizhanovskii,Phys. Rev. X , 031033 (2017). T. Boulier, S. Pigeon, E. Cancellieri, P. Robin,E. Giacobino, Q. Glorieux, and A. Bramati,Phys. Rev. B , 024503 (2018). T. Boulier, H. Terças, D. D. Solnyshkov, Q. Glo-rieux, E. Giacobino, G. Malpuech, and A. Bramati,Scientific Reports , 9230 EP (2015), article. R. Hivet, E. Cancellieri, T. Boulier, D. Ballarini, D. San-vitto, F. M. Marchetti, M. H. Szymanska, C. Ciuti, E. Gia-cobino, and A. Bramati, Phys. Rev. B , 134501 (2014). T. Gao, E. Estrecho, G. Li, O. A. Egorov,X. Ma, K. Winkler, M. Kamp, C. Schneider,S. Höfling, A. G. Truscott, and E. A. Ostrovskaya,Phys. Rev. Lett. , 097403 (2016). O. Zik, Z. Olami, and E. Moses,Phys. Rev. Lett. , 3868 (1998). I. Golding, Y. Kozlovsky, I. Cohen, and E. Ben-Jacob,Physica A: Statistical Mechanics and its Applications , 510 (1998). K. J. Lee, W. D. McCormick, Q. Ouyang, and H. L. Swin-ney, Science , 192 (1993). S. M. Troian, X. L. Wu, and S. A. Safran,Phys. Rev. Lett. , 1496 (1989). W. N. Reynolds, J. E. Pearson, and S. Ponce-Dawson,Phys. Rev. Lett. , 2797 (1994). E. Wertz, L. Ferrier, D. D. Solnyshkov, R. Johne,D. Sanvitto, A. Lemaître, I. Sagnes, R. Grousson, A. V.Kavokin, P. Senellart, G. Malpuech, and J. Bloch,Nature Phys. , 860 (2010). N. Bobrovska, E. A. Ostrovskaya, and M. Matuszewski,Phys. Rev. B , 205304 (2014). L. A. Smirnov, D. A. Smirnova, E. A. Ostrovskaya, andY. S. Kivshar, Phys. Rev. B , 235310 (2014). T. C. H. Liew, O. A. Egorov, M. Matuszewski,O. Kyriienko, X. Ma, and E. A. Ostrovskaya,Phys. Rev. B , 085413 (2015). I. S. Aranson and L. Kramer,Rev. Mod. Phys. , 99 (2002). S. Yoon, M. Sun, Y. G. Rubo, and I. G. Savenko, ArXive-prints (2018), arXiv:1806.03070 [cond-mat.mes-hall] . M. van Hecke, Phys. Rev. Lett. , 1896 (1998). R. Mirek, M. Król, K. Lekenta, J.-G. Rous-set, M. Nawrocki, M. Kulczykowski, M. Ma-tuszewski, J. Szczytko, W. Pacuski, and B. Pietka,Phys. Rev. B , 085429 (2017). J.-G. Rousset, B. Pietka, M. Król, R. Mirek,K. Lekenta, J. Szczytko, W. Pacuski, and M. Nawrocki,Phys. Rev. B , 125403 (2017). J.-G. Rousset, B. Pietka, M. Król, R. Mirek, K. Lekenta,J. Szczytko, J. Borysiuk, J. Suffczyński, T. Kazimierczuk,M. Goryca, T. Smoleński, P. Kossacki, M. Nawrocki, andW. Pacuski, Appl. Phys. Lett. , 201109 (2015). B. Pietka, D. Zygmunt, M. Król, M. R. Molas,A. A. L. Nicolet, F. Morier-Genoud, J. Szczytko,J. Łusakowski, P. Zieba, I. Tralle, P. Stepnicki,M. Matuszewski, M. Potemski, and B. Deveaud,Phys. Rev. B , 075309 (2015). J. Gaj, R. Planel, and G. Fishman,Solid State Commun. , 435 (1979). I. A. Shelykh, T. C. H. Liew, and A. V. Kavokin,Phys. Rev. B , 201306 (2009). P. Miętki and M. Matuszewski,Phys. Rev. B96