Noise-induced transition from superfluid to vortex state in two-dimensional nonequilibrium polariton condensates
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p Noise-induced transition from superfluid to vortex state in two-dimensionalnonequilibrium polariton condensates
Vladimir N. Gladilin and Michiel Wouters
TQC, Universiteit Antwerpen, Universiteitsplein 1, B-2610 Antwerpen, Belgium (Dated: September 30, 2019)We study the Berezinskii-Kosterlitz-Thouless mechanism for vortex-antivortex pair formationin two-dimensional superfluids for nonequilibrium condensates. Our numerical study is based on aclassical field model for driven-dissipative quantum fluids that is applicable to polariton condensates.We investigate the critical noise needed to create vortex-antivortex pairs in the systems, startingfrom a state with uniform phase. The dependence of the critical noise on the nonequilibrium andenergy relaxation parameters is analyzed in detail.
PACS numbers: 03.75.Lm, 71.36.+c
I. INTRODUCTION
With the advent of synthetic quantum systems, theinterest in driven-dissipative many-body systems hasgrown substantially in the last decade. Where particlesin ultracold atomic gases can to a very good approxi-mation be conserved, losses can be engineered [1] or areunavoidable in strongly interacting Bose gases [2]. In sys-tems based on electromagnetic degrees of freedom, bothin the microwave [3] and optical domain [4, 5], cavitylosses are often not negligible. In order to reach a steadystate, some driving of the system is then necessary tocompensate for the losses. This raises the question onthe modifications of the steady state with respect to thethermal equilibrium state in conservative systems.Here, we will consider the case of two-dimensionalweakly interacting bosons that are subject to single par-ticle losses, which are compensated by a nonresonantdrive. These systems are realized by microcavity polari-ton condensates [6], but there may be also the possibilityto construct them with ultracold atoms [1]. Microcav-ity polaritons are hybrid light-matter quasiparticles re-sulting from the strong coupling between an excitonictransition in a quantum well and a photonic mode inan enveloping microcavity. From their photonic compo-nent, they inherit a light effective mass, enhancing thespatial coherence, whereas from their excitonic compo-nent, they inherit interactions (see Refs. [7, 8] for theexperimental determination of the interaction constant).Under nonresonant excitation, a large density of excitonsis created, which subsequently relax to the lower polari-ton region. At equilibrium, the two-dimensional bose gasfeatures a Berezinskii-Kosterlitz-Thouless (BKT) phasetransition: for increasing temperature, thermally excitedvortex-antivortex pairs become unbound, resulting in theloss of superfluidity. A natural question is then how thistransition will be affected by losses and driving.From the experimental side, this physics was addressedin Ref. [9]. In their system with a long polariton life time,they did find a phase transition that was interpreted asthe binding to unbinding transition of vortex-antivortexpairs, reminiscent of the equilibrium situation. Already at the level of small phase fluctuations,nonequilibrium systems behave differently from theirequilibrium counterparts. Where for the latter, the phasedynamics is for small fluctuations to a good approxima-tion described by a linear equation, in the nonequilibriumcase, a nonlinear term appears, which brings them in theKardar-Parisi-Zhang universality class [10–16].Some first theoretical insight in the modification ofthe BKT transition due to the nonequilibrium conditioncan be gained by considering the modification of vorticeswhen going away from equilibrium. It turns out that thegain/losses introduce an additional current with the vor-tex core as its source. The phase profile is consequentlydeformed, resulting in a spiral wave [17]. This modifica-tion in the vortex flow field subsequently affects the inter-actions between vortices and antivortices: when the out-ward flows are more important than the usual azimuthalflows, the interaction between vortex and antivortex be-comes repulsive. These repulsive interactions hamper thevortex-antivortex recombination, enabling the formationof vortex-antivortex clusters with a very long life time[18]. A renormalization group based approach has shownthat these repulsions are fatal for the superfluid phase.The renormalization flow always goes toward the normalphase, even though this physics may manifest itself onlyat very large distances [10].In order to shed further light on the phase diagram of2D nonequilibrium polaritons, we resort here to numeri-cal simulations of the noisy generalized Gross-Pitaevskiiequation (ngGPE). This equation can be derived withinthe truncated Wigner approximation [19] or with theKeldysh field theory formalism [20, 21] and has beenwidely used to model nonresonantly pumped polaritoncondensates [22–27]. At equilibrium, when the Bose gasis fully characterized by its density, temperature and in-teraction constant, the critical temperature is found toequal T c ≈ π ~ n/ [ m ln(380 ~ /mg )], where n is the den-sity, m the mass and g the interaction strength [28].Away from equilibrium however, there are more micro-scopic parameters that enter the theoretical description.We investigate here how they affect the critical noisestrength (the analog of the temperature out of equilib-rium).In a previous theoretical study [18], based on a noisefree generalized GPE, we found that starting from aninitial state with a large number of vortices, several cansurvive in the steady state, because of the repulsive inter-actions between vortices and antivortices. We even foundthat these can form quite regular structures. With thisphysics understood, we will start here from the oppositeinitial condition with a homogeneous phase. In polaritoncondensates, such an initial condition can be achieved bysending a resonant pulse with a flat phase profile. As ex-pected, we find that only when sufficiently strong noiseis present, vortex-antivortex pairs can be formed in thesubsequent time evolution. We will show that the pairproduction shows a well defined noise threshold, allowingus to draw a phase diagram for the system.For systems that are far from equilibrium, we haveshown [18, 29] that the gGPE predicts a self-accelerationof vortices and production of new pairs in vortex-antivortex collisions, leading to chaotic dynamics. In thisparameter regime, we find that a moderate noise sup-presses this mechanism, leading to a stabilization of thesystem.The structure of the paper is as follows. In Sec. II,our model for nonequilibrium condensation is recapitu-lated. The phase diagram is discussed in Sec. III andconclusions are drawn in Sec. IV. II. MODEL
We consider nonresonantly excited two-dimensionalpolariton condensates. In the case of sufficiently fast re-laxation in the exciton reservoir, this reservoir can beadiabatically eliminated and the condensate is describedby the noisy generalized Gross-Pitaevskii equation [19–21, 27, 30, 31](i − κ ) ~ ∂ψ∂t = (cid:20) − ~ ∇ m + g | ψ | + i2 (cid:18) P | ψ | /n s − γ (cid:19)(cid:21) ψ + √ Dξ. (1)Here m is the effective mass and the contact interac-tion between polaritons is characterized by the strength g . The imaginary term in the square brackets on theright hand side describes the saturable pumping (withstrength P and saturation density n s ) that compen-sates for the losses ( γ ). We take into account the en-ergy relaxation κ in the condensate [32, 33]. The com-plex stochastic increments have the correlation function h ξ ∗ ( x, t ) ξ ( x ′ , t ′ ) i = 2 δ ( r − r ′ ) δ ( t − t ′ ).For polariton condensates, the validity of Eq. (1) isnot always straightfoward to justify. The repulsive inter-actions between the condensate and the exciton reservoirmay lead to an effective attraction between the polari-tons, leading to instability for a positive polariton mass[24, 34, 35]. This unstable state can be stabilized by anegative effective mass, that can be obtained in a po-lariton microcavity lattice [34]. We will not consider the instability physics in this work and assume positive massand positive interactions (but the physics remains unal-tered when the signs of the mass, interaction strengthand energy relaxation parameter ( κ ) are simultaneouslychanged).It is then convenient to rewrite Eq. (1) in a dimension-less form, by expressing the particle density | ψ | in unitsof n ≡ n s ( P/γ − ~ / ( gn ), lengthin units of ~ / √ mgn , and noise intensity in units of ~ n / (2 m ):(i − κ ) ∂ψ∂t = (cid:20) −∇ + | ψ | +i c − | ψ | ν | ψ | (cid:21) ψ + √ Dξ. (2)Besides the noise intensity D , equation (2) contains threeother dimensionless scalar parameters: κ characterizes,as described above, damping in the system dynamics, c = γ/ (2 gn s ) is a measure of the deviation from equilibrium,and ν = n /n s is proportional to the relative excess ofthe pumping intensity P over the threshold intensity. Inthe absence of noise, Eq. (2) has the steady state solution ψ = √ ne − int , where the density n satisfies κn = c (1 − n )1 + νn , (3)so that n is a decreasing function of κ .For weak noise, the fluctuations on top of this homo-geneous solution can in first approximation be analyzedwithin the Bogoliubov approach [36], but a more refinedanalysis has revealed that the dynamics of the phase fluc-tuations is in the KPZ universality class [12–15, 37, 38],where the phase nonlinearity cannot be neglected.Here, we will study the critical noise needed to cre-ate vortices. In order to address this problem, Eq.(2) is solved numerically using the same finite-differencescheme as in Ref. [29]. Specifically, we use periodicboundary conditions for a square of size L x = L y = 40with grid step equal to 0.2 corresponding to a cutoff inmomentum space equal to K c = 5 π . The model (1) is aclassical field model, that suffers from the usual ultravio-let catastrophe [28]. This implies that our results will becutoff dependent. Below, we will discuss this dependenceand how it will affect the comparison of our theoreticalresults with experiments.When analyzing the noise-induced BKT transition,each run starts from a uniform condensate distributionat D = 0. Then we apply noise with a fixed nonzerointensity D during a time interval t D . The detection ofvortex pairs in the presence of a strong noise is somewhattricky, but fortunately in the absence of noise their an-nihilation time is known to be rather long even at veryweak nonequilibrium [39, 40] and this time strongly in-creases with increasing c [18]. For these reasons, it ismore convenient to check the presence of vortex pairssometime later after switching off the noise. The usedtime delay (typically few our units of time) is sufficient FIG. 1. Distributions of the particle density n for c = 1 . ν = 1, κ = 0 . D = 0 .
015 and t D = 100 at the time moment t = 0 when the noise is switched off (a) and at t = 50 (b).Distribution of the phase θ of the order parameter at t = 5 isshown in panel (c). Panels (d), (e) and (f): same as in panels(a), (b) and (c), respectively, but for D = 0 . t D = 600and a shorter time delay after switching off the noise ( t = 5). for significant relaxation of the noisy component in thedensity and phase distributions of the condensate and, atthe same time, is too short for vortex pair annihilation.To determine the critical noise for the BKT transition, D BKT , we use the following criterion. If for a noise in-tensity D vortex pairs are present after a noise exposuretime t D (and hence D > D
BKT ), while for a certain noiseintensity D ′ < D no vortex pairs appear even at noiseexposures few times longer then t D , then D ′ lies eitherbelow D BKT or above D BKT and closer to D BKT thento D . Therefore, the critical noise intensity can be es-timated as D BKT = D ′ ± ( D − D ′ ). An illustration isdisplayed in Fig. 1. For D = 0 .
015 and noise exposuretime t D = 100, a noisy density distribution shown inFig. 1(a) evolves after switching off the noise into a pat-tern with clearly seen vortices and antivortices, whichpersist during a relatively long time [see Figs. 1(b) and(c)]. For a slightly lower noise intensity, D = 0 . t D = 600 ,after switching off the noise the corresponding noisy dis-tribution [Fig. 1(d)] rapidly relaxes towards a uniformvortex-free state [Figs. 1(e) and (f)]. According to ourapproach, the estimate for the critical noise in this caseis D BKT = 0 . ± . k c = 4 c = 0.3 D BK T (a)(b) c = 4 c = 0.3 n BK T d BK T k FIG. 2. Noise intensity D BKT corresponding to the BKT tran-sition [panel (a)] and average density of the condensate n BKT at D = D BKT [panel (b)] as a function of damping at ν = 1and two different values of the nonequilibrium parameter c .The error bars are smaller than the symbol size. Inset of panel(b): ratio d BKT = D BKT /n BKT for c = 4 (stars) and c = 0 . III. RESULTS AND DISCUSSION
In Fig. 2(a), the critical noise intensity correspondingto the BKT transition, D BKT , is shown as a functionof the damping parameter κ for the cases of moderate( c = 0 .
3) and strong ( c = 4) deviations from equilib-rium. In both cases the critical noise D BKT is seen toconsiderably increase with κ , despite the fact that theaverage density of the condensate at D = D BKT is amonotonously decreasing function of κ [see also Eq. (3)],this decrease being especially pronounced at smaller c [see Fig. 2(b)].The inset of Fig. 2(b) shows the ratio d BKT ≡ D BKT /n BKT as a function of κ . For the two very differ-ent values of the nonequilibrium parameter c , the curves d BKT ( κ ) appear to lie relatively close to each other. Onecan also notice that at c = 0 . d BKT on κ is nearly linear (except for the smallest values of κ ).At κ ≫ κ directly follows from an additional rescaling of time[ t → t/κ ; see Eq. (2)]. Our numerical results show thatthis linear dependence of d BKT on κ remains to goodapproximation unaltered even when κ < c = 0), the increase of the critical noisewith increasing κ can be understood by making the con-nection with the thermal equilibrium case. For vanish-ing nonequilibrium and κ ≫
1, our equation reduces tomodel A dynamics [41], which has a Boltzmann-Gibbssteady state distribution at temperature T = D/ κ . Inthis limit, the transition occurs at the equilibrium BKTtemperature [28], where T BKT = ηn , with η a numeri-cal cutoff dependent constant. For the critical noise, onethen obtains D BKT = 2 ηκn .Numerically, we have observed that for c →
0, the crit-ical noise obeys this relation, even for κ <
1. The factthat the dissipative part of the dynamics does not alterthe steady state can be understood from the following ar-gument. It is well established that pure Gross-Pitaevskiidynamics ( D = κ = c = 0) samples the thermal equilib-rium state in the microcanonical ensemble, at an energydetermined by the initial state. Similarly, the Langevindynamics [when omitting the i in the l.h.s. of Eq. (2)and taking c = 0] samples the phase space accordingto the canonical ensemble at a temperature determinedby the balance between noise and dissipation. Since thethermal state is the steady state of both the GP andLangevin dynamics, it is natural that the steady state ofthe combined dynamics is also at thermal equilibrium,with the temperature determined by the Langevin part.As implied by the results displayed in Fig. 2(a), thecritical noise intensity D BKT increases when movingaway from equilibrium. The dependence of D BKT onthe nonequlibrium parameter c is further illustrated inFig. 3(a) for the cases of weak ( κ = 0 .
01 and 0.1) andzero damping. At nonzero damping the increase of D BKT with c is partly related to the simultaneous increase of theaverage condensate density n BKT ( c ) shown in Fig. 3(b).The latter originates from the growing contribution of thepumping-loss term in Eq. (2) [the last term in the squarebrackets in Eq. (2)]. Indeed, this term tends to keep thecondensate density as close to 1 as possible. As a result,the degrading effect of damping or noise on the conden-sate density weakens with increasing c . This can be seen,e.g., from Fig. 2(b) by comparing to each other the valuesof n BKT at different c : for each fixed, not too small κ , thevalue of n BKT | c =4 is significantly larger than n BKT | c =0 . ,even though the density n BKT | c =4 corresponds to a con-siderably higher noise intensity than that for n BKT | c =0 . [see Fig. 2(a)].At the same time, as follows from the behavior of theratio d BKT ≡ D BKT /n BKT [see the inset of Fig. 3(a)],the increase of D BKT with c is considerably faster ascompared to that of n BKT . This means that the effectof pumping-loss processes on D BKT cannot be explainedsolely by the corresponding stabilization of the averagecondensate density against noise and damping. This be-comes even more evident when looking at the results for κ = 0. Indeed, while the critical noise D BKT demon-strates a clear rise with increasing c [see Fig. 3(a)], theaverage density n BKT ( c ) remains almost constant [seeFig. 3(b)]. At small c , the values of n BKT even slightly k = 0.1 k = 0.01 k = 0 D BK T (a) n BK T k = 0.1 k = 0.01 k = 0 c (b ) d BK T c d BK T c FIG. 3. Noise intensity D BKT corresponding to the BKT tran-sition [panel (a)] and average density of the condensate n BKT at D = D BKT [panel (b)] as a function of the nonequilibriumparameter c at ν = 1 and different values of the damping pa-rameter κ . The error bars are smaller than the symbol size.Inset of panel (a): ratio d BKT = D BKT /n BKT for κ = 0 . κ = 0 .
01 (circles). Inset of panel (b): ratio d BKT = D BKT /n BKT for κ = 0 (squares). The black linecorresponds to the linear fitting d BKT = ac , with a = 0 . decrease with c suggesting that the stabilizing effects ofthe pumping-loss processes do not completely compen-sate the average-density reduction caused by the increaseof the noise intensity at the BKT transition. However,an important aspect of the the pumping-loss processes,which is not directly reflected in the average density, isthat they impede formation of deep local suppressions ofthe condensate density. Those deep density suppressions,together with the appropriate phase gradient, are neces-sary prerequisites for vortex pair formation and hencefor the BKT transition. The discussed numerical resultsimply that just the stabilizing effect of the pumping-lossterm on the local density of the condensate governs theincrease of D BKT with c at κ →
0. At κ = 0 the depen-dence of D BKT on c becomes close to linear [see Fig. 3(a)],while for d BKT ( c ) a linear function d = ac with a = 0 . ψ ( x, t ) = p δn ( x, t ) e iθ ( x,t ) , one obtains inlinearized approximation for the Fourier components ∂∂t (cid:18) θ k n k (cid:19) = (cid:18) − κǫ k − ǫ k − ǫ k − c ν + 4 κ (cid:19) (cid:18) θ k n k (cid:19) + √ Dξ ( θ ) k √ Dξ ( n ) k ! , (4)where ǫ k = k . The white noise in Fourrier space hasthe correlation function h ξ ( θ ) ( k, t ) ξ ( θ ) ( k ′ , t ) i = 4 πδ ( k + k ′ ) δ ( t − t ′ ), analogous for ξ ( n ) , and h ξ ( θ ) ξ ( n ) i = 0.The steady state density fluctuations can be obtainedin closed form. In the limit of large k , they simplify to h n k n k ′ i = πδ ( k + k ′ ) 4 D (2 + κ ) κǫ k κ = 0 (5)2 πδ ( k + k ′ ) 8 D (1 + ν ) c κ = 0 (6)For the density fluctuations in real space, we then obtainfor a sufficiently large momentum cutoff K c h δn ( x ) i ∝ Dκ (2 + κ ) ln K c κ = 0 (7) D (1 + ν ) c K c κ = 0 (8)The local density fluctuations obtained here in the linearapproximation show behavior that is in line with whatwas observed numerically for the BKT transition. Den-sity fluctuations are first suppressed by the damping κ and for κ = 0 by the nonequilibrium parameter c .Note the very different dependence of both results onthe momentum cutoff. In the presence of a nonzerodamping, the cutoff dependence is a very weak logarithm,where in its absence, it becomes quadratic. This differ-ence is due to the fact that the momentum space densitytends to a constant in the κ = 0 case, where it decaysas k − for κ = 0. Even the latter is not fast enough toensure convergence in two dimensions, hence the loga-rithmic divergence of the density fluctuations.The values of D BKT , obtained from the numerical sim-ulations, correlate with the above analysis. At c = 1 . ν = 1, and κ = 0, an increase of the grid step from 0.2to 0.4 leads to a change of the numerically determined D BKT from 0 . ± . . ± . h n ( x ) i on K c in the absence of energy relaxation [see Eq. (8)]. At thesame time, for nonzero damping ( κ ≥ .
01) the relativedifference between the values of D BKT corresponding tothe grid step 0.4 and 0.2 is about 40% for κ = 0 .
01 anddoes not exceed 15% for κ ≥ .
1. The critical densityincreases somewhat when increasing the grid step, suchthat d BKT = D BKT /n BKT increases only by about 10%for κ ≥ .
1. Our analytical arguments and those of Ref.[28] suggest a logarithmic dependence on the cutoff, butwe were not able to test this scaling numerically due toa lack of accessible range in the grid spacing. Choos-ing the grid spacing much larger than 0.4, it becomes too coarse to give a good description of the vortex cores,while choosing a grid spacing smaller than 0.2 slows downthe calculations too much (the maximal kinetic energy in-creases quadratically with momentum cutoff). Moreover,decreasing the grid spacing below 0.2 would give a kineticenergy that is typically larger than what is required tojustify the lower polariton approximation.As seen from Eq. (2), the pumping-loss term increasesin magnitude when decreasing the parameter ν . One canexpect, therefore, that a decrease of ν leads to an increaseof the critical noise D BKT . Our simulations confirm thisexpectation but, at the same time, show that the influ-ence of ν on D BKT is relatively weak for nonvanishingdamping. At κ = 0 .
01 a decrease of ν by one order ofmagnitude (from 1 down to 0 .
1) results in an increase of D BKT approximately by 4% at c = 0 . c = 4. The corresponding increase of d BKT is about 20%for both c = 0 . ν on the critical noise is rather minor as compared to themuch stronger impact of κ and c .Let us now look in some more detail at the simultane-ous dependence of d BKT on κ and c . Taking into account,on the one hand, the nearly linear dependence d BKT on κ at large κ [see the inset in Fig. 2(b)] and, on the otherhand, the linear dependence on c , d BKT ≈ ac at κ → d BKT / ( κ + ac ). As seen from Fig. 4, all the resultsof our simulations lie in a rather narrow interval around1: 0 . < d BKT / ( κ + ac ) < .
3. This suggests that thesimple expression κ + 0 . c already can serve as a crudebut quite reasonable estimate for d BKT at any κ and c . Abetter approximation of the numerical results is obtained(see Fig. 4) by using the fitting function, that summarizesour numerical calculations rather accurately: d ∗ =0 . κ + 0 . c + κ . c . . κ . + 0 . c . + 1 . κ . c . (9)Since the numerical results for the dependence of d BKT on both κ and c are seen to be fitted quite well by d ∗ ( c, κ ),we may expect that this function can provide also mean-ingful results for inter- and extrapolation [see the insetin Fig. 4(b)].Finally, we address the question of whether the BKTtransition under consideration can be influenced by theeffect of vortex pair generation by moving vortices [18].This effect has been predicted to occur in stronglynonequilibrium condensates when self-accelerated vor-tices acquire sufficiently high velocities with respect tothe surrounding condensate. The results of [18] were ob-tained in the absence of noise. In our numerics contain-ing the noise term, we have observed that fluctuationsof the condensate density and currents tend to impedethe acceleration of vortices. Consequently, the genera-tion of vortex pairs in vortex collisions becomes impossi-ble at high noise intensities. This is illustrated in Fig. 5,which shows the upper boundary D m for the noise in-tensity range, where generation of new vortex pairs by d BK T / ( kk + a c ) , d * / ( + a c ) d BKT /( k + ac ) c = 4 c = 0.3 d */( k + ac ) c = 4 c = 0.3 k (a) d */( k + ac ) k c d BK T / ( kk + a c ) , d * / ( + a c ) (b) d BKT /( k + ac ) k = 0.1 k = 0.01 d */( k + ac ) k = 0.1 k = 0.01 c FIG. 4. Renormalized noise intensity d BKT / ( κ + ac ) (symbols)and its fitting by d ∗ / ( κ + ac ) (lines) as a function of the damp-ing parameter κ [panel (a)] and nonequilibrium parameter c [panel (b)].The inset shows the dependence of d ∗ / ( κ + ac ) onboth c and κ . Straight dashed lines correspond to the pa-rameter values covered by the curves in panels (a) and (b). self-accelerated vortices is possible, as a function of c at ν = 1 and κ = 0 .
01. The corresponding calculations areperformed, like in [18], starting with one pinned vortexpair in the region under consideration and simulating thedynamics during a time interval ∆ t ∼ D m on c is given in com-parison with the critical noise D BKT at the same κ and ν . As seen from Fig. 5, the region of c and D , wherepair generation by vortices is possible, lies at large c ( c > .
3) and has no overlap with the curve D BKT ( c ).Of course, the latter remains true for any κ ≥ .
01. In-deed, the critical noise D BKT increases with κ . At thesame time, damping slows down vortex motion and thusimpedes pair generation by moving vortices. As a result,the aforementioned region will shrink with increasing κ .Our simulations show that already at κ as small as 0.1the velocities of vortices are insufficient for pair gener-ation. So, we can conclude that at κ ≥ .
01 the BKTtransition is not affected by the processes of pair gener-ation by moving vortices. As concerns the limit of very D BKT D D c m FIG. 5. Noise intensity D BKT ( c ) corresponding to the BKTtransition in comparison with D m ( c ), the maximum noise in-tensity, which still allows for vortex pair generation by movingvortices, at ν = 1 and κ = 0 . weak damping κ < .
01, the nearly linear behavior of D BKT ( c ), obtained for κ = 0, manifests no visible pecu-liarities in the region of large c , where the processes ofpair generation by vortices become possible. This impliesthat, also in the limit κ →
0, the BKT transition is notconsiderably affected by these processes.
IV. CONCLUSIONS
We have investigated the critical noise for the spon-taneous formation of unbound vortex-antivortex pairs ina driven-dissipative bosonic system where particle lossesare compensated by nonresonant pumping. In this work,we have focused on the noise needed to form vortex-antivortex pairs starting from a uniform phase. In ourmodel, the critical noise strength depends – with a suit-able choice of units – on three dimensionless parameters:the energy relaxation rate ( κ ), the nonequilibrium pa-rameter ( c ) and the gain saturation parameter ( ν ).In the absence of energy relaxation, the nonequilibriumparameter determines the critical noise strength, but thiscritical value depends quadratically on the momentumspace cutoff, giving our numerical results limited predic-tivity for specific experiments. In the presence of suffi-ciently strong energy relaxation, the cutoff dependence ismuch weaker and experimentally relevant results can beextracted from the numerics in this case. In this regime,the critical noise strength increases both with energy re-laxation (as expected from equilibrium calculations) andwith nonequilibrium parameter c . The latter dependencecould seem counterintuitive, because with increasing c ,vortices and antivortices repel each other at large dis-tances, which could favor their unbinding. We interpretthe impeding of vortex-antivortex formation further awayfrom nonequilibrium as a consequence of the reduction ofthe density fluctuations. The effect of nonequilibrium onthe formation of vortices is therefore opposite to its effecton the annihilation : the life time of existing vortices isdramatically enhanced by nonequilibrium [18].Quantitatively, we have found that the effect of thenonequilibrium parameter is actually small for κ ≫ . c , which is satisfied when damping is not too smalland the system is not too far from equilibrium. Very farfrom equilibrium, self acceleration of vortices can lead tothe production of new vortex antivortex pairs. We haveshown here that this pair production ceases for increasingnoise. Our numerical simulations were performed for systemswith periodic boundary conditions whose size is compa-rable to the systems employed in current experiments.The study of the thermodynamical limit of infinite sys-tem size remains a challenge for both theoretical analysisand experimental investigation. ACKNOWLEDGEMENTS
We thank Sebastian Diehl for a stimulating discus-sion. This work was financially suppported by the FWOOdysseus program. [1] R. Labouvie, B. Santra, S. Heun, and H. Ott, Phys. Rev.Lett. , 235302 (2016).[2] P. Makotyn, C. E. Klauss, D. L. Goldberger, E. A. Cor-nell and D. S. Jin, Nat. Phys. , 116 (2014); C. Eigen,J. A. P. Glidden, R. Lopes, E. A. Cornell, R. P. Smithand Z. Hadzibabic Nature , 221 (2018).[3] M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch,and A. A. Houck, Phys. Rev. X , 011016 (2017).[4] I. Carusotto and C. Ciuti, Rev. Mod. Phys. , 299(2013).[5] F. Brennecke, R. Mottl, K. Baumann, R. Landig, T.Donner, and T. Esslinger, Proceedings of the NationalAcademy of Sciences , 11763 (2013).[6] J. Kasprzak, M. Richard, S. Kundermann, A. Baas,P. Jeambrun, J. Keeling, F. Marchetti, M. Szyma´nska,R. Andre, J. Staehli, V. Savona, P. B. Littlewood, B.Deveaud and Le Si Dang, Nature , 409 (2006).[7] A. Delteil, T. Fink, A. Schade, S. H¨ofling, C. Schneider,A. Imamo˘glu, Nat. Materials , 219 (2019).[8] S. R. K. Rodriguez, W. Casteels, F. Storme, N. Car-lon Zambon, I. Sagnes, L. Le Gratiet, E. Galopin, A.Lemaˆıtre, A. Amo, C. Ciuti, and J. Bloch, Phys. Rev.Lett. , 247402 (2017).[9] D. Caputo, D. Ballarini, G. Dagvadorj, C. S. Mu˜noz,M. De Giorgi,L. Dominici, K. West, L. N. Pfeiffer,G. Gigli, F. P. Laussy, M. H. Szyma´nska and D. San-vitto, Nat. Mat. ,145 (2018).[10] G. Wachtel, L. M. Sieberer, S. Diehl and E. Altman,Phys. Rev. B , 104520 (2016).[11] L. M. Sieberer, G. Wachtel, E. Altman and S. Diehl,Phys. Rev. B , 104521 (2016).[12] D. Squizzato, L. Canet and A. Minguzzi, Phys. Rev. B , 195453 (2018).[13] V. N. Gladilin, K. Ji, and M. Wouters, Phys. Rev. A ,023615 (2014).[14] L. He, L. M. Sieberer, E. Altman, and S. Diehl, Phys.Rev. B 92, 155307 (2015).[15] E. Altman, L. M. Sieberer, L. Chen, S. Diehl, and J.Toner, Phys. Rev. X 5, 011017 (2015).[16] L. M. Sieberer, M. Buchhold, and S. Diehl, Reports onProgress in Physics 79, 096001 (2016).[17] I. S. Aranson and L. Kramer, Rev. Mod. Phys. , 99(2002).[18] V. N. Gladilin and M. Wouters, arXiv:1810.06365. [19] M. Wouters and V. Savona, Phys. Rev. B , 165302(2009).[20] M. H. Szyman´ska, J. Keeling, and P. B. Littlewood, Phys.Rev. B , 195331 (2007).[21] L. M. Sieberer, S. D. Huber, E. Altman, S. Diehl, Phys.Rev. B , 134310 (2014).[22] K. Kalinin, P. G. Lagoudakis, N. G. Berloff, Phys. Rev.B , 094512 (2018).[23] P. Comaron, G. Dagvadorj, A. Zamora, I. Carusotto, N.P. Proukakis, M. H. Szyma´nska, Phys. Rev. Lett. ,095302 (2018).[24] N. Bobrovska, M. Matuszewski, K. S. Daskalakis, S. A.Maier, S. K´ena-Cohen, ACS Photonics , 111 (2017).[25] 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, Nat.Comm. , 2944 (2018).[26] H. Ohadi, Y. del Valle-Inclan Redondo, A. J. Ramsay, Z.Hatzopoulos, T. C. H. Liew, P. R. Eastham, P. G. Sav-vidis, J. J. Baumberg, Phys. Rev. B , 195109 (2018).[27] J. Keeling and N. G. Berloff, Phys. Rev. Lett. ,250401 (2008).[28] N. Prokof’ev, O. Ruebenacker, and B. Svistunov Phys.Rev. Lett. 87, 270402 (2001).[29] V. N. Gladilin and M. Wouters, New J. Phys. , 105005(2017).[30] M. Wouters and I. Carusotto, Phys. Rev. Lett. ,140402 (2007).[31] N. Bobrovska and M. Matuszewski, Phys. Rev. B ,035311 (2015).[32] L.P. Pitaevskii, Zh. Eksp. Teor. Fiz. , 408 (1958) [Sov.Phys. JETP , 282 (1959)].[33] M. Wouters, New J. Phys. , 075020 (2012).[34] F. Baboux, D. De Bernardis, V. Goblot, V.N. Gladilin, C.Gomez, E. Galopin, L. Le Gratiet, A. Lematre, I. Sagnes,I. Carusotto, M. Wouters, A. Amo, J. Bloch, Optica ,1163 (2018).[35] N. Bobrovska, E. A. Ostrovskaya and Michal Ma-tuszewski, Phys. Rev. B , 205304 (2014).[36] A. Chiocchetta and I. Carusotto, EPL, , 67007(2013).[37] K. Ji, V. N. Gladilin, and M. Wouters, Phys. Rev. B ,045301 (2015).[38] L. He and S. Diehl, New J. Phys. 19, 115012 (2017). [39] M. Kulczykowski and M. Matuszewski, Phys. Rev. B ,075306 (2017).[40] P. Comaron, G. Dagvadorj, A. Zamora, I. Carusotto, N.P. Proukakis, M. H. Szymanska, Phys. Rev. Lett. , 095302 (2018).[41] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys.49