# Instability-induced formation and non-equilibrium dynamics of phase defects in polariton condensates.

###### Abstract

We study, theoretically and numerically, the onset and development of modulational instability in an incoherently pumped spatially homogeneous polariton condensate. Within the framework of mean-field theory, we identify regimes of modulational instability in two cases: 1) Strong feedback between the condensate and reservoir, which may occur in scalar condensates, and 2) Parametric scattering in the presence of polarization splitting in spinor condensates. In both cases we investigate the instability induced textures in space and time including non-equilibrium dynamics of phase dislocations and vortices. In particular we discuss the mechanism of vortex destabilization and formation of spiraling waves. We also identify the presence of topological defects, which take the form of half-vortex pairs in the spinor case, giving an “eyelet” structure in intensity and dipole type structure in the spin polarization. In the modulationally stable parameter domains, we observe formation of the phase defects in the process of condensate formation from an initially spatially incoherent low-density state. In analogy to the Kibble-Zurek type scaling for nonequilibrium phase transitions, we find that the defect density scales with the pumping rate.

###### pacs:

71.36.+c,42.65.Sf,78.67.-n,03.75.Mn## I Introduction

The formation of complex spatio-temporal patterns and textures is a particularly intriguing phenomenon occurring in a diverse range of physical systems.Cross1993 (); StaliunasBook () The patterns commonly arise in nonlinear dissipative systems driven far from equilibrium. Amongst examples of such systems, exciton-polaritons in semiconductor microcavities have emerged as a hybrid light-matter system with strongly nonlinear properties. Effects such as Bose-Einstein condensation,Kasprzak2006 (); Balili2007 (); Lai2007 () superfluidity Carusotto2004 (); Amo2009a (); Amo2009b () and the formation of solitons Egorov2009 (); Egorov2010 (); Egorov2011 (); Sich2011 (); Ostrovskaya2012 () have been well-documented in the literature.Deng2010 (); Carusotto2013 ()

A particularly important nonlinear effect in the context of pattern formation in microcavities is the parametric scattering of resonantly excited polaritons in planar semiconductor microcavities.Savvidis2000 (); Tartakovskii2000 (); Stevenson2000 (); Saba2001 () The pair scattering of pairs of polaritons to different states in reciprocal space allows a homogeneous polariton field to spontaneously break translational symmetry. This effect enables the formation of ordered hexagonal/triangular lattices, predicted theoretically Saito2013 (); Egorov2014 (); Werner2014 () and observed experimentally Ardizzone2013 () under different excitation conditions, as well as lattices of breathing solitons.Egorov2013 ()

Under non-resonant/incoherent excitation, vortex lattices were predicted to occur in harmonic traps Keeling2008 (); Borgh2010 () and later observed in experiments involving multiple excitation spots.Cristofolini2013 () The formation of multi-lobed Manni2011 () and vortex-antivortex patterns Manni2013 () under ring shaped excitation, as well as sunflower ripples Christmann2012 () excited by a narrow pump spot have also been reported. In these examples, the translational symmetry of the system is already broken by the chosen shape of the pump-spot and/or the presence of a gradient in the potential of polaritons (as in, for example, the case of harmonic traps Balili2007 ()).

In this work we consider the possibility of spontaneous breaking of translational symmetry and pattern formation in planar microcavities excited by a spatially homogeneous incoherent pump. The pumping creates a reservoir of “hot” exciton-like polaritons, which form a polariton condensate through a stimulated scattering process. The translational symmetry breaking is triggered by linear instability of the homogeneous condensate to spatial modulations, and the nonlinear evolution of the unstable state leads to formation of spatial patterns. We consider two different mechanisms of such modulational instability (MI) in this system. The first one arises when the polariton condensate has a strong feedback effect on the reservoir, in the form of reservoir depletion due to stimulated scattering of reservoir excitons into polaritons. In this case, while polariton-polariton interactions are repulsive, the essentially saturable nature of exciton-polariton interactions may lead to effectively attractive nonlinearity for sufficiently low pump powers.Smirnov2014 () This effective focusing nonlinearity in the system naturally leads to MI of a spatially homogeneous state, which was established in several previous studies for both quasi-1D and 2D geometry.Wouters2007 (); Borgh2010 (); Carusotto2013 (); Littlewood06 (); Malpuech14 (); Smirnov2014 (); Li14 ().

Modulational instability is also known in spin-1 Bose-Einstein condensates of ultracold atoms due to parametric coupling Robins2001 (); Zhang2005 (); Li2005 (); Doktorov2007 () and nonlinear interactions between spin components,Matuszewski2010 () as well as in 1D exciton-polariton condensates with a spin (polarization) degree of freedom.Kamchatnov2013 () Although polariton systems are non-conservative and non-equilibrium, the two-component spin degree of freedom of polaritons does allow a second mechanism of MI, which works also in the defocusing regime where a strong condensate-reservoir feedback is unnecessary. A circularly polarized excitation splits the energy of the and states due to anisotropic interactions Vladimirova2010 (); Takemura2014 () occurring between the spin polarized reservoir and condensed polaritons. This splitting sets the foundation for a parametric scattering process as polaritons in an initially homogeneous state with wavevector on the upper spin-split branch can now scatter to degenerate non-zero wavevector states on the lower branch, reminiscent of experiments under resonant excitation in triple microcavities Diederichs2006 () or experiments in one-dimensional polariton systems.Xie2012 () While such a scattering process is not strictly allowed in isotropic cavities as it would violate spin conservation, the presence of sample anisotropy, which typically causes an additional linear polarization splitting and hybridization of the and branches, relaxes this limitation.

By considering the stability of the steady states of the system to weak perturbations, we find the zones of MI in the two different regimes. In the scalar case, where MI is derived from the condensate-reservoir feedback, we find that the homogeneous state breaks its translational symmetry and forms a turbulent state of phase dislocations, i.e., vortices. Unlike the previously studied cases, vortices do not form as the result of thermal fluctuations Giorgetti2007 () or scattering on disorder Liew2008 (). Rather, the spatial fragmentation of the initially homogeneous condensate due to the development of MI creates multiple interference between polariton flows generated by the randomly distributed sources, which leads to the development of multiple phase dislocations, similar to the scenario previously considered for multiple pump spots Berloff10 () and highly inhomogeneous trapped polariton condensates Borgh2010 (); Borgh2012 ().

In the case of modulationally stable background, we show that multiple phase singularities can appear as a result of mean-field evolution of an initial white noise state, which mimics a pre-condensate state lacking spatial and phase coherence. Remarkably, formation of multiple vortices in this scenario seems to be analogous to, but not the same as, the Kibble-Zurek mechanism, which acts during the quench through a phase transition to the Bose-Einstein condensation (BEC) Kibble1976 (); Zurek1985 (). Indeed, the later describes the formation of boundary defects between different domains of condensate which develop an independent phase rather than inheriting it from the neighbouring spatial domains Zurek2009 (); Lamporesi2013 ().

We stress that the process of defect formation during nonequilibrium condensation of exciton-polaritons does not follow the scenario of the Kibble-Zurek mechanism Zurek1985 (); Dziarmaga2010 (); Matuszewski2014 (). The main difference is that in the latter, it is assumed that the system is initially in thermal equilibrium, and is driven out of equilibrium only in the vicinity of the phase transition Zurek1985 (); Dziarmaga2010 (). The process is divided into three phases, corresponding to adiabatic-impulse-adiabatic evolution. In nonequilibrium condensation, the system is far from equilibrium at the outset, and the transition to the quasi-equilibrium (condensed) state occurs only after crossing the critical point.

Nevertheless, the Kibble-Zurek mechanism and the defect formation in nonequilibrium systems have much in common. In both cases, defects are created due to symmetry breaking in separate parts of the system which cannot communicate in a finite time. In both cases, there is a competition of two timescales existing in the system, which results in the same algebraic forms of power-law scalings for the number of defects and their characteristic creation time Matuszewski2014 (). In the polariton condensation case, the quench time is replaced by the timescale of the formation of the condensate, which is controlled by the external pumping rate. We refer the reader to Sec. IV.4 and Ref. Matuszewski2014, for the detailed description of the process.

Regardless of the mechanism of the vortex formation, either as the result of the MI development or as a result of transition to BEC, we show that the presence of the incoherent reservoir affects substantially both stability of vortices and their collective dynamics even for the case of a stable homogeneous background. As a consequence, the vortices can lose their radial symmetry and develop either into spatially localized rotating phase dislocations or into non-localized spiraling waves.

A similar situation occurs in the spinor case, although multiple branches of modulationally unstable and stable solutions are present. Defects in the spin polarization of the condensate may appear even in the modulationaly stable regime. Such structures move randomly in the microcavity plane and are composed of half-vortex Rubo2007 () half-antivortex pairs,Manni2012 () exhibiting an associated dipole type spin texture. We predict that the density of vortices grows with increasing pump power similarly to the Kibble-Zurek scaling behaviour.

The paper is organised as follows. In Sec. II, we describe the mathematical model of a semiconductor microcavity operating in the strong coupling regime under incoherent homogeneous optical pump of a circular polarization. Then, in Sec. III, we study the stability and collective dynamics of phase dislocations in a single-component polariton condensate. Here the dynamics is mostly affected by the modulational instability originating from the strong feedback between the condensate and reservoir. In Sec. IV, we report a numerical analysis of the condensate dynamics in the presence of polarization splitting in spinor condensates. In Sec. IV.4, we study the defect formation and the scaling laws for their density in analogy to the Kibble-Zurek mechanism.

## Ii Theoretical Model

Let us begin by considering the incoherent excitation of a spinor polariton system, i.e., the system where the polarization degree of freedom is significant. The scalar case, which is valid when only one spin component is populated, is then easily obtained by removing one of the spin components. It is worth recalling that experiments with a circularly polarized optical pump have resulted in the excitation of a circularly polarized polariton condensate at the pump position, in both 2D Kammann2012 () and 1D Anton2014 () samples. A theoretical model can be based on the generalized Gross-Pitaevskii approach,Wouters2007 () where a condensate of exciton-polaritons can be described by the wavefunctions, and , of the and circularly polarized states, respectively:

(1) |

(2) |

Here we assume that the circularly polarized pumping creates a circularly polarized reservoir, , with dynamics described by the rate equation:

(3) |

where represents the pumping rate. Non-linear interactions between polaritons are characterized by and , representing the interaction strengths between parallel and antiparallel spins Vladimirova2010 (); Takemura2014 (), respectively. Similar parameters, and , characterize the blueshift caused by the circularly polarized reservoir. is the polariton decay rate. represents a linear polarization splitting, which has been reported in several experimental studies Klopotowski2006 (); Krizhanovskii2006 (); Amo2009 () and can take values of 0-0.2 meV.Klopotowski2006 (); Comitti2014 () Even larger values can be expected by the application of magnetic fields (in Voight configuration).

We note that in the limit , it is possible to proceed by adiabatic elimination of the reservoir dynamics.Borgh2010 (); Larre2013 () However we do not make such an approximation here since the reservoir dynamics is important for our further analysis.

## Iii Non-equilibrium Dynamics in the Scalar Case

### iii.1 Modulational instability of the homogeneous steady state.

If we consider a circularly polarized pump in a microcavity with negligible polarization splitting () then the model can be reduced to a scalar one () without the second polarization component (). First, we study a scalar steady-state homogeneous solution (HS) of the system of Eqs. (1)-(3) and discuss its stability. For the sake of generality we allow also the HSs with nonzero transversal momenta , which have the form of travelling waves

(4) |

where the condensate energy is given by . The HS becomes nontrivial provided that the external pump compensates for all losses and overcomes the threshold value:Wouters2007 () . The coherent exciton-polariton density and incoherent reservoir density are given by and , respectively.

The linear stability analysis of the homogeneous steady state of our scalar system and its modifications has been previously performed by many authors.Wouters2007 (); Carusotto2013 (); Littlewood06 (); Malpuech14 (); Smirnov2014 (); Li14 () For our choice of the system parameters, the linear stability analysis shows that the HS becomes modulationally (dynamically) unstable within a pump interval just above the threshold value of the pump [Fig. 1(a)] (details of the analysis are given in Appendix A). This MI is associated with the parametrical generation of field components with nonzero momenta . Figure 1(b) presents the linear growth rate of the unstable perturbations as a function of their momenta .

It has been shown recently Smirnov2014 () that this MI is associated with the effective attractive nonlinearity induced by the saturation of the incoherent reservoir. Indeed, based on the intuition gained from paradigm nonlinear models, such as the Schrödinger equation with a Kerr-nonlinearity, one expects that the existence of the MI requires a focusing nonlinearity.Saito2001 (); AlKhawaja2002 (); Salasnich2003 (); Carr2004 () However, owing to repulsive interactions between excitons, the nonlinear behaviour of an exciton-polariton condensate is akin to that of optical waves in a defocusing media. This seeming contradiction clearly elucidates the influence of the open-dissipative nature of the system on the nonlinear behaviour,Smirnov2014 () which requires inclusion of an incoherent reservoir of “hot” excitons. To illustrate this influence we consider a nonlinear energy shift induced by both the coherent exciton-polaritons and the incoherent reservoir

(5) |

We note that the polariton density corresponds to a steady-state solution which, in general, is not necessarily given by the homogeneous value . The reservoir intensity follows this steady-state solution .

The first term of Eq. (5) describes the blue shift originating from the repulsive exciton-exciton interaction, whereas the second term describes the saturation of the reservoir. We define the effective nonlinearity coefficient as . Then the nonlinear response is effectively focusing provided that the coefficient is negative . Otherwise, if the nonlinear response is defocusing or repulsive [see thick solid lines in the Fig. 2(a)]. In the vicinity of the steady-state HS with intensity this coefficient takes the form (cf. Eq. (36) in Ref.[Smirnov2014, ]):

(6) |

The effective nonlinear coefficient [Eq. (6)] changes sign for a pump value . It means that in the vicinity of the HS the nonlinear response changes character from effectively focusing to effectively defocusing. As a result the modulationally unstable HS becomes stable exactly at this point . The condition gives a general criterion for the appearance of MI derived in Ref. [Smirnov2014, ]:

(7) |

Within the MI interval () indicated in Fig. 1(a), the HS experiences spontaneous translational symmetry breaking, resulting in the formation of non-uniform turbulent states of the condensate.

The result of a direct numerical calculation is shown in Fig. 1(c). In these calculations and those presented throughout the manuscript we use a square grid with at least points covering the plot area. We make use of an adaptive step Adams-Bashforth-Moulton procedure, which was previously found to be consistent with fixed step methods with an integration step of ps. To double check the results of our numerical simulations we repeated them with higher precision using up to grid points. Unless stated otherwise, periodic boundary conditions are applied.

The strongly non-equilibrium state in Fig. 1(c) includes one- and two-dimensional phase dislocations which move chaotically and overlap. Therefore this dynamics can be characterized as a “strong” turbulence regime with overlapping defects (see Ref. [Berloff10, ] and references therein). Two-dimensional Fourier transformation shows that the most part of the spatial momenta are bounded within the ring with the radius given by the cutoff momenta of the unstable modes [Fig. 1(d)]. Note that direct numerical simulations of the model (1,3) for different initial conditions did not reveal the formation of stationary periodical patterns, known for the coherently pumped polaritonic systems Saito2013 (); Ardizzone2013 (); Egorov2014 (); Werner2014 (). This is due to destabilization caused by fluctuations in the exciton reservoir. The modulational instability discussed later in the spinor case (section IV) does not depend on having a dynamic reservoir and can largely be reproduced assuming a static reservoir Keeling2008 (). This allows an explicit testing of the effect of the dynamic reservoir in the spinor case, where we find that periodic patterns are possible with a static reservoir but are prevented as soon as the reservoir density is allowed to evolve spatiotemporally. We assume that it is the freedom for density fluctuations to appear in the reservoir that lead to the disruption of regular patterns in the condensate also in the scalar case.

### iii.2 Single vortices in the dynamically stable regime.

The reservoir contribution becomes negligible in the limit of very strong pump . We expect that in this case the nonlinear dynamics is very similar to that known for the conservative systems, including the formation of the stable phase dislocations and vortices. Indeed the scalar version of the equations possesses vortex solutions within the stability interval of the HS for [see Figs. 2(b) and (c)]. Similar to the conservative case, one can define a characteristic length or an effective healing length in the vicinity of the HS :

(8) |

The healing length is a typical length scale over which can change significantly. It also gives the typical size of the vortices. More precisely the full width at half maximum (FWHM) of the vortices is approximately given by at least for a strong enough pumping [ in Fig. 2(b)]. It is worth to mention that the spatial oscillations of the real and imaginary parts of the amplitude profiles of the vortices go substantially beyond the healing length [Fig. 2(d)], especially for a moderate pumping in the vicinity of the MI. This indicates permanent energy and polariton exchange within the vortex profile. The presence of these intrinsic fluxes essentially changes the interaction dynamics between vortices FraserNJP2009 () providing a purely dissipative mechanism for their mutual repulsion.

The radial symmetry of the vortex phase is broken also for an inhomogeneous pump, for instance, a Gaussian pump Borgh2010 (); Ostrovskaya2012 (). As a result, the vortices have spiraling phases indicating again the permanent exchange of particles between different points within the resonator plane.

The stability analysis, discussed in the context of the HSs, can also be applied to the vortices, at least in the limit of very broad states. Indeed, following the healing length , the vortex width increases for smaller values of the pump and diverges in the vicinity of the MI [ in Fig. 2(b)]. The formal condition for the effective focusing nonlinearity () is always satisfied within the vortex profile close enough to its core. Therefore the steady-state condensate on a circle with a fixed radius and density can become unstable against spatially modulated perturbations. The periodical boundary conditions on the circle restrict the number of available momenta to the values , where is an integer. In the limit of large radius the values of the unstable momenta can be approximated by those values calculated for the HSs [see Fig. 1(b)]. Therefore the polariton condensate becomes unstable provided that at least one of the available momenta () is smaller than the cut-off value for MI. This allows estimation of the vortex diameter where the instability just sets in, i.e., . Therefore the vortices are stable provided that their diameters do not exceed the maximal value given by

(9) |

It is shown in Fig. 2(b) that the vortex diameter exceeds the maximal value within the pump interval ps in the vicinity of the MI domain. As we will see below (in the next subsection III.3) they become unstable and lose their radial symmetry.

### iii.3 Collective dynamics of vortices

Beyond the MI instability interval for the non-trivial HS solution is stable. However, at the onset of condensation when it passes into the mean-field regime, the polariton state is spatially incoherent, with random phase. This leads to the generation of phase defects as the system goes through the condensate transition, in analogy to the Kibble-Zurek theory. To go beyond the mean-field approximation and describe first and second order spatial coherences when crossing the condensation threshold, one can make use of stochastic classical field approaches.Wouters2008 (); Wouters2009 () These describe the pre-condensate state as an ensemble of fluctuating white noise states (governed by a stochastic Gross-Pitaevskii equation). The projection onto the classical condensate state upon condensation selects a particular realization of the noise. After the condensate has formed we assume that fluctuations are weak in comparison to the condensate mean-field and can be neglected. This approach reproduces the typical establishment and coherent evolution of topological defects in polariton condensate experiments Lagoudakis2011 (); Manni2012 (). We note that the account of fluctuations throughout the evolution would be important for describing accurately spectral polariton properties or the polariton photoluminescence below threshold Wouters2008 (). Fluctuations can in principle shift the phase boundaries between stable and unstable regions Johne2009 (), however, these shifts are expected to be limited and have little further effect on the dynamics.

For a pumping substantially above the MI threshold the dynamics is dominated by the defocusing nonlinearity of a pure condensate (). The influence of dissipative dynamics of the reservoir becomes less important. Using the scalar version of the governing equations (1), (3) we calculated the condensate dynamics starting with a spatially incoherent state given by a small-amplitude white noise. Our initial condition mimics a particular realisation of the stochastic pre-condensate state. Similar to the conservative limit we observed the formation of a spatially coherent condensate accompanied by the spontaneous formation of the vortices which move chaotically and interact with each other. Two vortices with equal (opposite) topological charges repel (attract) each other. Two attracting vortices mutually annihilate if the distance between them becomes smaller than the healing length. Therefore, as discussed in Ref. [Berloff10, ], the whole number of vortices drops gradually with time and approaches zero, at least for a very strong pump ps.

In contrast, for a weaker pump, the number of phase singularities converges eventually to some constant value indicating the formation of a coherent state with a finite number of dislocations, i.e., a superfluid turbulence [Fig. 3(a)]. A snapshot of the intensity and phase profiles shows a state of well distinguishable vortices [Figs. 3(b),(c)]. The vortices move chaotically, interacting with their neighbours, and, in general, sustain a dynamical equilibrium. It is remarkable that the average separations between nearest vortices remains more or less constant for this particular interval of pump values. This means that there exists some equilibrium distance between vortices. Apparently the influence of the dissipative effects and the condensate flows (mentioned in the previous subsection III.2) are substantial.FraserNJP2009 () The out-going condensate flows from the vortex centers hinder attraction between vortices and their annihilation. For even weaker pump the dissipative effects become stronger and, as a consequence, the average distance between vortices under dynamical equilibrium becomes even smaller [Figs. 3(d),(e) and (f)]. We note that this dynamical equilibrium forms over a long time scale exceeding hundreds of nanoseconds. Therefore the Kibble-Zurek-like scaling law does not describe the number of vortices in this regime (see Sec. IV.4 below).

Even though the HSs are stable, the nonlinear dynamics of vortices is strongly affected by the reservoir saturation dynamics. The numerical modelling shows that the vortices themselves become unstable and develop into radially asymmetric rotating structures, as can be seen on the snapshot profiles in Figs. 3(e),(f) (a movie showing the time dynamics is available in the Supplemental Material SM_UnstableVortices ()). This is in agreement with the destabilization scenario for the vortices discussed in the previous subsection. Indeed the vortex size exceeds the maximal diameter given by Eq. (9) and therefore becomes unstable.

### iii.4 Formation of spiraling waves

In the vicinity of the MI the vortices are strongly unstable and the initial noise first develops into non-uniform dynamical states similar to those which appear for the modulationally unstable background [Figs. 4(a,b)]. In this “strong” turbulence regime Berloff10 () the characteristic distance between vortices is substantially smaller than their typical core size, so the vortices are not structured and the chaotic behavior is seen on the level of a single vortex. However, after some sufficiently long time of about several thousand polariton lifetimes, the system switches spontaneously into a more regular regime, characterized by the formation of a single spiraling topological dislocation [Fig. 4] (a movie showing the time dynamics is available in the Supplemental Material SM_UnstableVortices ()). This spiraling topological state drives away other phase dislocations in the system and covers eventually the whole computational window provided that ps. Similar spiraling waves are known for other non-equilibrium dissipative systems.Aranson1992 (); Huber1992 () In the simplest case they are solutions of the complex Ginzburg-Landau equation.

It is worth to articulate the differences between the spiraling topological states and the vortices. First, there is no rotational symmetry in the profiles [Figs. 4]. Second, these topological states are not stationary and experience uniform rotation of the density with a typical rotation period of about ps. Third, far from the center the profile converges to the homogeneous traveling wave solution given by Eq. (4) with the amplitude and a nontrivial momentum . Apparently this solution is characterized by a permanent radial flux of exciton-polaritons from the vortex center to the periphery and therefore resembles a point-like source of ring waves. We note that these spiraling waves can exist only in nonequilibrium dissipative systems.

An important question remains whether the spiraling state with a non-zero orbital angular momentum can emerge from the initially non-rotating turbulent state. Indeed, according to the conservation law of the total orbital angular momentum, the phase dislocations appear in pairs which is also valid for the turbulent state considered here. Apparently the local intensity fluctuations break the symmetry between the two dislocations within a pair in such a way that only one of them develops into the spiraling wave. The periodical boundary conditions (in and directions) were used in our numerical modelling. However, we have confirmed that the spiraling waves appear with the same probability independently on the computational window size and particular realisations of the initial seeding noise. Moreover, since some phase dislocations are always present [see Figs. 4(e) and (f)], the total orbital angular momentum of the condensate within the computational window remains zero. We performed additional numerical simulations of the condensate dynamics under a localized pump with a “flat-top” shape in the form of a super-Gaussian intensity distribution. It turned out that the spiraling waves appear also for the localized pump where the condensate density vanishes at the boundaries of the computational window. These calculations serve as a solid proof of the existence of the spiraling waves independently of a particular choice of the numerical boundary conditions.

In general the out-going radiation from the center of the topological solution repels the local inhomogeneities of the profile and other topological solutions. This gives an additional purely dissipative mechanism which enforces a long range spatial coherence in non-equilibrium systems operating in the regime of “strong” turbulence.

## Iv Non-equilibrium dynamics in the Spinor Case

### iv.1 Stability of Homogeneous States

In the presence of non-zero polarization splitting (), the population of component appears even for a fully polarized pumping. To emphasize the difference to the nonlinear dynamics discussed above (see Sec. III) we consider system parameters which do not satisfy the MI criteria (7) and, therefore, guarantee the stability of homogeneous solutions in the scalar limit.

Homogeneous stationary solutions can be found by substituting trial solutions in the form into Eqs. (1) and (2). This gives four stationary equations for real and imaginary parts of the amplitudes. These are supplemented by the requirement that the time derivative in Eq. (3) vanishes for a steady state, that is, . Noting that the phase reference of the system can be freely chosen, we can for simplicity set , which allows one to find the relation from one of the four stationary equations. The remaining three equations can be solved for the remaining unknown quantities: , and .

The dependence of the stationary HSs on the pumping power is shown in Figs. 5(a) and (b). Linear stability of the steady states can be determined by the standard extension of the Bogoliubov-de Gennes analysis Wouters2007 (); Carusotto2013 (); Littlewood06 (); Malpuech14 (); Smirnov2014 (); Li14 (); Kamchatnov2013 () by considering perturbations in the polariton and reservoir fields of the form and , respectively. The details of the derivation can be found in Appendix A.

The stationary states labelled and in Fig. 5(a,b) are unstable to fluctuations with . This “single-mode instability” indicates that even in a confined system (e.g., micropillar) the homogeneous steady state would be linearly unstable in the Lyapunov sense and any spatially homogeneous perturbation would grow. The stationary states labelled and are unstable to spatial modulations with non-zero wavevectors . These states would be stable in a confined system, with no spatial degrees of freedom, however, in a 2D planar system the states and undergo parametric scattering. The state labelled is fully stable, as the imaginary part of remains negative for all wavevectors (see Appendix A).

The “S”-shape of the curves , and is characteristic of multistability, which is a common feature of resonantly excited microcavities Gippius2004 (); Baas2004 (); Whittaker2005 (); Bajoni2008 () but less studied under the non-resonant or incoherent excitation Kyriienko2014 () that we consider here. While multistability is strictly only present in the confined system, since the and states are unstable in the presence of spatial degrees of freedom, they can still give rise to different (non-stationary) configurations under the same excitation conditions.

Fig. 5(c) shows the regions of MI in the system when the or branches are excited. The branch begins once the threshold for polariton condensation is passed (indicated by the grey rectangle), such that only weak pump intensities are needed to see the effects of MI.

### iv.2 Spin Textures due to Modulational Instability

Due to the presence of MI, we can expect the fragmentation of the homogeneous density of the condensate and spontaneous formation of spin textures, even in the presence of homogeneous pumping. Solving Eqs. (1)-(3) numerically when the system is excited just above threshold on the branch, we obtain the spin texture shown in Fig. 6. In analogy to the scalar case, the texture comprises a simultaneous modulation of the intensity and phase in the system, which oscillate with multiple frequencies. In addition to the scalar case, there is also an appearance of a non-uniform polarization, despite the fact that the pumping of the system is homogeneous in both intensity and polarization.

### iv.3 Spin Defects in the Dynamically Stable Regime

When the system is excited with a larger pump power, the system tends to follow the stable branch (S) in Fig. 5. In this case one can expect a spatially HS due to the stability of the (S) branch, however, defects present in the initial state after transition to condensation (in simulations taken as a low intensity white noise as in sections III.3 and III.4) are trapped in the system and stabilize with the structure shown in Fig. 7.

The defects have a non-trivial density structure with localized maxima inside of an otherwise circular shaped drop in density. These “eyelets” are not fixed in their locations, but move randomly in the plane with a typical speed on the order of m/ps. The eyelets also possess a characteristic spin polarization with a dipole type shape, as shown in Fig. 7(b). The structure of the eyelets can also be understood when looking at the phase distribution, which is shown in Figs. 7(c) and (d) for the and polarized components, respectively. Here it is clear that each eyelet is formed from a pair of half vortices Rubo2007 () – one appearing in each spin component.

In addition to the slow drift, the eyelets preserve their shape while undergoing a faster periodic rotation with a period close to (a movie of the motion is available in the Supplemental Material). Apparently this periodical motion of two bound vortices with opposite spins is induced by the polarization splitting in spinor condensates. It has been shown recently Egorov2014b () that similar spinor effects can evoke a uniform motion of bound polariton solitons in coherently driven microresonators.

### iv.4 Scaling Laws for the Defect Density

The number of vortices generated in the system outside of the MI region, during the mean-field evolution of an initially noisy, low-density state, is found to grow with the pump intensity, as shown in Fig. 8. The number of vortices is counted at a time after spatial coherence is established in the system, where the average polariton and reservoir density achieves a steady-state. Note that at very long-times there may be further recombination of vortex-antivortex pairs.

By taking advantage of the universality of dynamics of the system in the vicinity of the phase transition, we find approximate scaling laws governing the number of defects created during the transition, see solid line in Fig. 8(b). The scaling laws in the case of exciton-polaritons have similar forms as the ones obtained from the argument of Kibble and Zurek.Kibble1976 (); Zurek1985 () However, the dynamics of defect formation is different due to the fact that the initial state of the system is strongly out of equilibrium.Matuszewski2014 () While in the Kibble-Zurek mechanism the phase transition is assumed to begin in the initial state that is close to thermal equilibrium, here the white-noise initial state is dominated by fluctuations. In both cases, spontaneous symmetry breaking occurs differently in separate regions of space which cannot communicate on how the symmetry is broken due to the finite timescale of the process. On the borders between these separate regions defects can appear in the form of domain walls, vortices, or more sophisticated structures, depending on the dimensionality of the system and the form of the order parameter.Kibble1976 (); Zurek1985 (); Other_KZ ()

In the case of dynamics described by (1)-(3), starting from the initial low density white-noise state, the growth of polariton density initially occurs without establishing any coherence (a movie showing the dynamics corresponding to the establishment of the state in Fig. 8a can be found in the supplemental material). This corresponds to the first stage of the process, where there is practically no nonlinearity and no -dependence of the growth rate. The defects are created in the second stage, when fluctuations are suppressed due to nonlinear interactions. In this process, regions of ordered phase (or patches of a condensate) appear out of the initial disordered strongly fluctuating phase. Due to spatial and polarization symmetry breaking, defects are created on the borders between these condensed regions.

We assume that in the emerging ordered regions described above the relative phase between the two polarization components is approximately equal to and densities of and are similar. This assumption is in agreement with the numerical data presented in Figs. 7(c) and 7(d), where the phase is equal to away from the eyelets, and is related to the fact that such configuration minimizes the energy in Eqs. (1-2). We substitute to obtain

(10) |

Next, we assume that the reservoir quickly adjusts to the change of ,

(11) |

Let us now assume that there exists a patch of approximately constant polariton field (or a condensate seed) at . We can describe small fluctuations around by

(12) |

where and are the chemical potential and the growth rate of the patch, and , represent the fluctuations. Expanding (11) in Taylor series around , we can rewrite (10) as

(13) |

with the accuracy of the order of . Here , , and where is the effective exciton decay rate. The above equation has the same form as the one derived in Ref. [Matuszewski2014, ], apart from the negligible dependence of on . The Bogoliubov-de Gennes modes are , and the mode frequencies

(14) |

where , the saturation parameter , and . This spectrum has the property that modes with high momenta are strongly damped (the imaginary part of the frequency is negative), in contrast to the initial linear dynamics where no -dependence of the imaginary part of the spectrum was present. We can define a characteristic momentum cutoff for the modes that are strongly damped, which scales with the parameters as . Any fluctuations with momenta higher than will be suppressed, while fluctuations with lower momentum can form regions of ordered condensate phase. To estimate the scaling of the number of defects, two different limiting cases can be considered. In the first case, we assume that the vortices are formed when the condensate density is already near to its equilibrium value . Here we obtain

(15) |

Note that the above scaling does not have a power-law form, which is due to the fact that the transition is effectively nonlinear.

In the opposite limit,Matuszewski2014 () we assume that the vortices are formed when the condensate density is still very small . In this case we obtain

(16) |

The two estimates [Eqs. (15) and (16)] are compared to the numerical results in Fig. 8(b). The numerically obtained scaling appears to be intermediate between the two extreme cases.

## V Conclusions

In this paper we presented a comprehensive theoretical study of non-equilibrium dynamics of polariton condensates in incoherently pumped semiconductor microcavities. We have anticipated two different destabilization mechanisms that govern nonlinear dynamics of this system. The first arises when the polariton condensate has a strong feedback effect on the reservoir of incoherent “hot” polaritons. The second one is associated with the parametric scattering in the presence of polarization splitting in a spinor condensate. Both mechanisms result in the formation of phase defects, i.e. vortices, triggered by the modulational instability of the homogeneous condensate.

In the scalar case, we have shown that the presence of the incoherent reservoir can affect substantially both the vortex stability and their mutual collective dynamics. In particular this can lead to the formation of rotating dislocations or delocalized spiraling waves.

In the spinor two-component condensate we have identified the presence of topological defects, which take the form of half-vortex pairs, giving an “eyelet” structure in intensity and dipole type structure in the spin polarization.

In the case when the phase defects are formed in the dynamically (modulationally) stable regime, as a result of the condensate formation from an initial spatially and phase-incoherent state, we find that the defect density scales with the pumping rate in analogy to the Kibble-Zurek type scaling for non-equilibrium phase transitions.

## Acknowledgements

T.C.H.L. acknowledges support from the Lee Kuan Yew Fellowship. O.A.E. acknowledges financial support by the Deutsche Forschungsgemeinschaft (DFG project EG344/2-1) and the Thuringian Ministry for Education, Science and Culture (TMESC project B514-11027). M.M. acknowledges
support from the National Science Center grant DEC-2011/01/D/ST3/00482. E.A.O. acknowledges funding by the Australian Research Council (ARC).

## Appendix A Linear Stability Analysis

The stability of the solutions can be checked using the standard approach of applying perturbations, in the form , . Substitution into Eqs. (1)-(3) and collecting terms linear in the small amplitudes and we have:

(17) | ||||

(18) |

(19) |

where we have defined:

(20) | ||||

(21) |

We consider perturbations in the polariton and reservoir fields of the form and , respectively. Here and are complex amplitudes, while is a real amplitude. is a complex eigenvalue to be determined. Stability of the original solutions occurs if such that the perturbation decays.

Substitution of and into Eqs. (17)-(19) and collection of terms oscillating as and yields a set of five coupled equations, which represent an eigenvalue problem for . In matrix form:

(22) |

where and .

The imaginary and real parts of the perturbation spectra, , are shown in Fig. 9(a,b). The energies of the stationary states are illustrated by thin horizontal lines in Fig. 9(a). Where the imaginary parts are positive, thick curves denote modulational instability at the given wavevectors . The stationary state labelled is unstable to fluctuations with .

## References

- (1) M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993).
- (2) K. Staliūnas and V. J. Sánchez-Morcillo, Transverse patterns in Nonlinear Optical Resonators (Springer-Verlag Berlin Hidelberg, 2003).
- (3) 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 Le Si Dang, Nature 443, 409 (2006).
- (4) R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West, Science 316, 1007 (2007).
- (5) C. W. Lai, N. Y. Kim, S. Utsunomiya, G. Roumpos, H. Deng, M. D. Fraser, T. Byrnes, P. Recher, N. Kumada, T. Fujisawa, and Y. Yamamoto, Nature 450, 529 (2007).
- (6) I. Carusotto and C. Ciuti, Phys. Rev. Lett. 93, 166401 (2004).
- (7) A. Amo, J. Lefrere, S. Pigeon, C. Adrados, C. Ciuti, I. Carusotto, R. Houdre, E. Giacobino, and A. Bramati, Nature Phys. 5, 805 (2009).
- (8) A. Amo, D. Sanvitto, F. P. Laussy, D. Ballarini, E. del Valle, M. D. Martin, A. Lemaitre, J. Bloch, D. N. Krizhanovskii, M. S. Skolnick, C. Tejedor, and L. Viña, Nature 457, 291 (2009).
- (9) O. A. Egorov, D. V. Skryabin, A. V. Yulin, and F. Lederer, Phys. Rev. Lett. 102, 153904 (2009).
- (10) O. A. Egorov, A. V. Gorbach, F. Lederer, and D. V. Skryabin, Phys. Rev. Lett. 105, 073903 (2010).
- (11) O. A. Egorov, D. V. Skryabin, and F. Lederer, Phys. Rev. B 84, 165305 (2011).
- (12) M. Sich, D. N. Krizhanovskii, M. S. Skolnick, A. V. Gorbach, R. Hartley, D. V. Skryabin, E. A. Certa-Méndez, K. Biermann, R. Hey, and P. V. Santos, Nature Photon. 6, 50 (2012).
- (13) E. A. Ostrovskaya, J. Abdullaev, A. S. Desyatnikov, M. D. Fraser, and Y. S. Kivshar, Phys. Rev. A 86, 013636 (2012).
- (14) H. Deng, H. Haug, and Y. Yamamoto, Rev. Mod. Phys. 82, 1489 (2010).
- (15) I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).
- (16) P. G. Savvidis, J. J. Baumberg, R. M. Stevenson, M. S. Skolnick, D. M. Whittaker, and J. S. Roberts, Phys. Rev. Lett. 84, 1547 (2000).
- (17) A. I. Tartakovskii, D. N. Krizhanovskii, and V. D. Kulakovskii, Phys. Rev. B 62, 13298R (2000).
- (18) R. M. Stevenson, V. N. Astratov, M. S. Skolnick, D. M. Whittaker, M. Emam-Ismail, A. I. Tartakovskii, P. G. Savvidis, J. J. Baumberg, and J. S. Roberts, Phys. Rev. Lett. 85, 3680 (2000).
- (19) M. Saba, C. Ciuti, J. Bloch, V. Thierry-Mieg, R. André, Le Si Dang, S. Kundermann, A. Mura, G. Bongiovanni, J. L. Staehli, and B. Deveaud, Nature 414, 731 (2001).
- (20) H. Saito, T. Aioi, and T. Kadokura, Phys. Rev. Lett. 110, 026401 (2013).
- (21) O. A. Egorov, A. Werner, T. C. H. Liew, E. A. Ostrovskaya, and F. Lederer,Phys. Rev. B 89, 235302 (2014).
- (22) A. Werner, O. A. Egorov, and F. Lederer, Phys. Rev. B 89, 245307 (2014).
- (23) V. Ardizzone, O. Lewandowski, M. H. Luk, Y. C. Tse, N. H. Kwong, A. Lucke, M. Abbarchi, E. Boudin, E. Galopin, J. Bloch, A. Lemaitre, P. T. Leung, P. Roussignol, R. Binder, J. Tignon, and S. Schumacher, Sci. Rep. 3, 3016 (2013).
- (24) O. A. Egorov and F. Lederer, Phys. Rev. B 87, 115315 (2013).
- (25) J. Keeling and N. G. Berloff, Phys. Rev. Lett. 100, 250401 (2008).
- (26) M. O. Borgh, J. Keeling, and N. G. Berloff, Phys. Rev. B 81, 235302 (2010).
- (27) P. Cristofolini, A. Dreismann, G. Christmann, G. Franchetti, N. G. Berloff, P. Tsotsis, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Phys. Rev. Lett. 110, 186403 (2013).
- (28) F. Manni, K. G. Lagoudakis, T. C. H. Liew, R. André, and B. Deveaud-Plédran, Phys. Rev. Lett. 107, 106401 (2011).
- (29) F. Manni, T. C. H. Liew, K. G. Lagoudakis, C. Ouellet-Plamondon, R. André, V. Savona, and B. Deveaud, Phys. Rev. B 88, 201303 (2013).
- (30) G. Christmann, G. Tosi, N. G. Berloff, P. Tsotsis, P. S. Eldridge, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Phys. Rev. B 85, 235303 (2012).
- (31) L. A. Smirnov, D. A. Smirnova, E. A. Ostrovskaya, and Y. S. Kivshar, Phys. Rev. B 89, 235310 (2014).
- (32) M. H. Szymańska, J. Keeling, and P. B. Littlewood, Phys. Rev. Lett. 96, 230602 (2006).
- (33) D. D. Solnyshkov, H. Tercas, K. Dini, and G. Malpuech, Phys. Rev. A 89, 033626 (2014).
- (34) G. Li, M. D. Fraser, A. Yakimenko, and E. A. Ostrovskaya, arXiv:1408.1568 (2014).
- (35) M. Wouters and I. Carusotto, Phys. Rev. Lett. 99, 140402 (2007).
- (36) N. P. Robins, W. Zhang, E. A. Ostrovskaya, and Y. S. Kivshar, Phys. Rev. A 64, 021601(R) (2001).
- (37) W. Zhang, D. L. Zhou, M. S. Chang, M. S. Chapman, and L. You, Phys. Rev. Lett. 95, 180403 (2005).
- (38) L. Li, Z. Li, B. A. Malomed, D. Mihalache, and W. M. Liu, Phys. Rev. A 72, 033611 (2005).
- (39) E. V. Doktorov, V. M. Rothos, and Y. S. Kivshar, Phys. Rev. A 76, 013626 (2007).
- (40) M. Matuszewski, Phys. Rev. Lett. 105, 020405 (2010).
- (41) P.-É. Larré, N. Pavloff, and A. M. Kamchatnov, Phys. Rev. B 88, 224503 (2013).
- (42) M. Vladimirova, S. Cronenberger, D. Scalbert, K. V. Kavokin, A. Miard, A. Lemaître, J. Bloch, D. Solnyshkov, G. Malpuech, and A. V. Kavokin, Phys. Rev. B, 82, 075301 (2010).
- (43) N. Takemura, S. Trebaol, M. Wouters, M. T. Portella-Oberli, and B. Deveaud, Phys. Rev. B, 90, 195307 (2014).
- (44) C. Diederichs, J. Tignon, G. Dasbach, C. Ciuti, A. Lemaitre, J. Bloch, Ph. Roussignol, and C. Delalande, Nature 440, 904 (2006).
- (45) W. Xie, H. Dong, S. Zhang, L. Sun, W. Zhou, Y. Ling, J. Lu, X. Shen, and Z. Chen, Phys. Rev. Lett. 108, 166401 (2012).
- (46) L. Giorgetti, I. Carusotto, and Y. Castin, Phys. Rev. A 76, 013613 (2007).
- (47) T. C. H. Liew, Yu. G Rubo, and A. V. Kavokin, Phys. Rev. Lett. 101, 187401 (2008).
- (48) N. G. Berloff, Turbulence in exciton-polariton condensates, arXiv:1010.5225 (2010).
- (49) M. O. Borgh, G. Franchetti, J. Keeling, and N. G. Berloff, Phys. Rev. B, 86, 035307 (2012).
- (50) T. W. B. Kibble, J. Phys. A 9, 1387 (1976).
- (51) W. H. Zurek, Nature (London) 317, 505 (1985).
- (52) W. H. Zurek, Phys. Rev. Lett. 102, 105702 (2009).
- (53) G. Lamporesi, S. Donadello, S. Serafini, F. Dalfovo, and G. Ferrari, Nature Phys. 9, 656 (2013).
- (54) M. Matuszewski and E. Witkowska, Phys. Rev. B 89, 155318 (2014).
- (55) J. Dziarmaga, Adv. Phys., 59, 1063 (2010).
- (56) Yu. G. Rubo, Phys. Rev. Lett., 99, 106401 (2007).
- (57) F. Manni, K. G. Lagoudakis, T. C. H. Liew, R. André, V. Savona, and B. Deveaud, Nat. Commun. 3, 1309 (2012).
- (58) E. Kammann, T. C. H. Liew, H. Ohadi, P. Cilibrizzi, P. Tsotsis, Z. Hatzopoulos, P. G. Savvidis, A. V. Kavokin, and P. G. Lagoudakis, Phys. Rev. Lett. 109, 036404 (2012).
- (59) C. Antón, S. Morina, T. Gao, P. S. Eldridge, T. C. H. Liew, M. D. Martín, Z. Hatzopoulos, P. G. Savvidis, I. A. Shelykh, and L. Viña, arXiv:1410.8417 (2014).
- (60) L. Klopotowski, M. D. Martín, A. Amo, L. Viña, I. A. Shelykh, M. M. Glazov, G. Malpuech, A. V. Kavokin, and R. André, Solid State Commun. 139, 511 (2006).
- (61) D. N. Krizhanovskii, D. Sanvitto, I. A. Shelykh, M. M. Glazov, G. Malpuech, D. D. Solnyshkov, A. Kavokin, S. Ceccarelli, M. S. Skolnick, and J. S. Roberts, Phys. Rev. B 73, 073303 (2006).
- (62) A. Amo, T. C. H. Liew, C. Adrados, E. Giacobino, A. V. Kavokin, and A. Bramati, Phys. Rev. B 80, 165325 (2009).
- (63) V. S. Comitti, M. B. E. da Silva, and F. M. Matinaga, Solid State Comm. 193, 37 (2014).
- (64) P. É. Larré, N. Pavloff, and A. M. Kamchatnov, Phys. Rev. B 88, 224503 (2013).
- (65) H. Saito and M. Ueda, Phys. Rev. Lett. 86, 1406 (2001).
- (66) U. Al Khawaja, H. T. C. Stoof, R. G. Hulet, K. E. Strecker, and G. B. Partridge, Phys. Rev. Lett. 89, 200404 (2002).
- (67) L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. Lett. 91, 080405 (2003).
- (68) L. D. Carr and J. Brand, Phys. Rev. Lett. 92, 040401 (2004).
- (69) M. D. Fraser, G. Roumpos, and Y. Yamamoto, New J. Phys. 11 113048 (2009).
- (70) M. Wouters, I. Carusotto, and C. Ciuti, Phys. Rev. B 77, 115340 (2008).
- (71) M. Wouters and V. Savona, Phys. Rev. B 79, 165302 (2009).
- (72) K. G. Lagoudakis, F. Manni, B. Pietka, M. Wouters, T. C. H. Liew, V. Savona, A. V. Kavokin, R. André, and B. Deveaud-Plédran, Phys. Rev. Lett. 106, 115301 (2011).
- (73) R. Johne, N. S. Maslova, and N. A. Gippius, Solid State Commun., 149, 496 (2009).
- (74) See Supplemental Material at http:…. to view the time dynamics associated with Figs. 3, 4 and 7.
- (75) I. S. Aranson, L. Aranson, L. Kramer, and A. Weber, Phys. Rev. A 46, 2992(R) (1992).
- (76) G. Huber, P. Alstrøm, and T. Bohr, Phys. Rev. Lett. 69, 2380 (1992).
- (77) N. A. Gippius, S. G. Tikhodeev, V. D. Kulakovskii, D. N. Krizhanovskii, and A. I. Tartakovskii, Europhys. Lett. 67, 997 (2004).
- (78) A. Baas, J. P. Karr, H. Eleuch, and E. Giacobino, Phys. Rev. A 69, 023809 (2004).
- (79) D. M. Whittaker, Phys. Rev. B 71, 115301 (2005).
- (80) D. Bajoni, E. Semenova, A. Lemaitre, S. Bouchoule, E. Wertz, P. Senellart, S. Barbay, R. Kuszelewicz, and J. Bloch, Phys. Rev. Lett. 101, 266402 (2008).
- (81) O. Kyriienko, E. A. Ostrovskaya, O. A. Egorov, I. A. Shelykh, and T. C. H. Liew, Phys. Rev. B, 90, 125407 (2014).
- (82) O. A. Egorov and F. Lederer, Optics Lett. 39, 4029 (2014).
- (83) I. Chuang et al., Science 251, 1336 (1991); K. Pyka, J. Keller, H. L. Partner, R. Nigmatullin, T. Burgermeister, D.-M. Meier, K. Kuhlmann, A. Retzker, M. B. Plenio, W. H. Zurek, A. del Campo, and T. E. Mehlstübler, Nat. Commun. 4, 2291 (2013).