Noise-induced transition from superfluid to vortex state in two-dimensional nonequilibrium polariton condensates
Vladimir N. Gladilin Michiel Wouters
















TQC, Universiteit Antwerpen, Universiteitsplein 1,
B-2610 Antwerpen, Belgium




We study  the Berezinskii-Kosterlitz-Thouless mechanism for
vortex-antivortex pair formation in two-dimensional superfluids for
nonequilibrium condensates. Our numerical study is based on a
classical 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,
starting from a state with uniform phase. The dependence of the
critical noise on the nonequilibrium and energy relaxation
parameters is analyzed in detail.


03.75.Lm, 71.36.+c

Noise-induced transition from superfluid to vortex state in two-dimensional nonequilibrium polariton condensates
    Michiel Wouters
    December 30, 2023
================================================================================================================





§ INTRODUCTION


With the advent of synthetic quantum systems, the interest in
driven-dissipative many-body systems has grown substantially in the
last decade. Where particles in ultracold atomic gases can to a very
good approximation be conserved, losses can be engineered
<cit.> or are unavoidable in strongly interacting Bose gases
<cit.>. In systems based on electromagnetic degrees of
freedom, both in the microwave <cit.> and optical
domain <cit.>, cavity losses are often
not negligible. In order to reach a steady state, some driving of
the system is then necessary to compensate for the losses. This
raises the question on the modifications of the steady state with
respect to the thermal equilibrium state in conservative systems.

Here, we will consider the case of two-dimensional weakly
interacting bosons that are subject to single particle losses, which
are compensated by a nonresonant drive. These systems are realized
by microcavity polariton condensates <cit.>, but there
may be also the possibility to construct them with ultracold atoms
<cit.>. Microcavity polaritons are hybrid light-matter
quasiparticles resulting from the strong coupling between an
excitonic transition in a quantum well and a photonic mode in an
enveloping microcavity. From their photonic component, they inherit
a light effective mass, enhancing the spatial coherence, whereas
from their excitonic component, they inherit interactions (see Refs.
<cit.> for the experimental determination of
the interaction constant). Under nonresonant excitation, a large
density of excitons is created, which subsequently relax to the
lower polariton region. At equilibrium, the two-dimensional bose gas
features a Berezinskii-Kosterlitz-Thouless (BKT) phase transition:
for increasing temperature,  thermally excited vortex-antivortex
pairs become unbound, resulting in the loss of superfluidity. A
natural question is then how this transition will be affected by
losses and driving.

From the experimental side, this physics was addressed in Ref.
<cit.>. In their system with a long polariton life time,
they did find a phase transition that was interpreted as the binding
to unbinding transition of vortex-antivortex pairs, reminiscent of
the equilibrium situation.

Already at the level of small phase fluctuations, nonequilibrium
systems behave differently from their equilibrium counterparts.
Where for the latter, the phase dynamics is for small fluctuations
to a good approximation described by a linear equation, in the
nonequilibrium case, a nonlinear term appears, which brings them in
the Kardar-Parisi-Zhang universality class
<cit.>.

Some first theoretical insight in the modification of the BKT
transition due to the nonequilibrium condition can be gained by
considering the modification of vortices when going away from
equilibrium. It turns out that the gain/losses introduce an
additional current with the vortex core as its source. The phase
profile is consequently deformed, resulting in a spiral wave
<cit.>. This modification in the vortex flow field
subsequently affects the interactions between vortices and
antivortices: when the outward flows are more important than the
usual azimuthal flows, the interaction between vortex and antivortex
becomes repulsive. These repulsive interactions hamper the
vortex-antivortex recombination, enabling the formation of
vortex-antivortex clusters with a very long life time
<cit.>. A renormalization group based approach has shown
that these repulsions are fatal for the superfluid phase. The
renormalization flow always goes toward the normal phase, even
though this physics may manifest itself only at very large distances
<cit.>.

In order to shed further light on the phase diagram of 2D
nonequilibrium polaritons, we  resort here to numerical simulations
of the noisy generalized Gross-Pitaevskii equation (ngGPE). This
equation can be derived within the truncated Wigner approximation
<cit.> or with the Keldysh field theory formalism
<cit.> and has been widely used to model
nonresonantly pumped polariton condensates <cit.>. At equilibrium, when
the Bose gas is fully characterized by its density, temperature and
interaction constant, the critical temperature is found to equal
T_c ≈ 2πħ^2 n/[m ln(380 ħ^2/mg)], where n is
the density, m the mass and g the interaction strength
<cit.>. Away from equilibrium however, there are more
microscopic parameters that enter the theoretical description. We
investigate here how they affect the critical noise strength (the
analog of the temperature out of equilibrium).

In a previous theoretical study <cit.>, based on  a noise
free generalized GPE, we found that starting from an initial state
with a large number of vortices, several can survive in the steady
state, because of the repulsive interactions between vortices and
antivortices. We even found that these can form quite regular
structures. With this physics understood, we will start here from
the opposite initial condition with a homogeneous phase. In
polariton condensates, such an initial condition can be achieved by
sending a resonant pulse with a flat phase profile. As expected, we
find that only when sufficiently strong noise is present,
vortex-antivortex pairs can be formed in the subsequent time
evolution. We will show that the pair production shows a well
defined noise threshold, allowing us to draw a phase diagram for the
system.

For systems that are far from equilibrium, we have shown
<cit.> that the gGPE predicts a self-acceleration of
vortices and production of new pairs in vortex-antivortex
collisions, leading to chaotic dynamics. In this parameter regime,
we find that a moderate noise suppresses this mechanism, leading to
a stabilization of the system.

The structure of the paper is as follows. In Sec. <ref>,
our model for nonequilibrium condensation is recapitulated. The
phase diagram is discussed in Sec. <ref> and conclusions
are drawn in Sec. <ref>.



§ MODEL 


We consider nonresonantly excited two-dimensional polariton
condensates. In the case of sufficiently fast relaxation in the
exciton reservoir, this reservoir can be adiabatically eliminated
and the condensate is described by the noisy generalized
Gross-Pitaevskii equation <cit.>

    ( i-κ)ħ∂ψ/∂ t =      [-ħ^2∇^2/2m +g |ψ|^2 . 
          .+ i/2(P/1+|ψ|^2/n_s-γ) ] ψ+√(D)ξ .

Here m is the effective mass and the contact interaction between
polaritons is characterized by the strength g. The imaginary term
in the square brackets on the right hand side describes the
saturable pumping (with strength P and saturation density n_s)
that compensates for the losses (γ). We take into account the
energy relaxation κ in the condensate <cit.>. The
complex stochastic increments have the correlation function ⟨ξ^*(x,t) ξ(x',t') ⟩=2 δ( r- r')
δ(t-t').

For polariton condensates, the validity of Eq. (<ref>) is not
always straightfoward to justify. The repulsive interactions between
the condensate and the exciton reservoir may lead to an effective
attraction between the polaritons, leading to instability for a
positive polariton mass <cit.>.
This unstable state can be stabilized by a negative effective mass,
that can be obtained in a polariton microcavity lattice
<cit.>. We will not consider the instability physics in
this work and assume positive mass and positive interactions (but
the physics remains unaltered when the signs of the mass,
interaction strength and energy relaxation parameter (κ) are
simultaneously changed).

It is then convenient to rewrite Eq. (<ref>) in a dimensionless
form, by expressing the particle density |ψ|^2 in units of
n_0≡ n_s(P/γ-1 ), time in units of ħ/(gn_0),
length in units of ħ/√(2mgn_0), and noise intensity in
units of ħ^3 n_0/(2m):

    ( i-κ)∂ψ/∂ t=      [-∇^2
    +|ψ|^2  ψ^2/ψ^2. 
          .+ ic1-|ψ|^2/1+ν |ψ|^2] ψ
    +√(D)ξ .

Besides the noise intensity D, equation (<ref>) contains
three other dimensionless scalar parameters: κ characterizes,
as described above, damping in the system dynamics,
c=γ/(2gn_s) is a measure of the deviation from equilibrium,
and ν=n_0/n_s is proportional to the relative excess of the
 pumping intensity P over the threshold intensity. In the
absence of noise, Eq. (<ref>) has the steady state solution
ψ = √(n) e^- i n t, where the density n satisfies

    κ n = c(1-n)/1+ν n,

so that n is a decreasing function of κ.

For weak noise, the fluctuations on top of this homogeneous solution
can in first approximation be analyzed within the Bogoliubov
approach <cit.>, but a more refined analysis has
revealed that the dynamics of the phase fluctuations is in the KPZ
universality class <cit.>, where the phase nonlinearity cannot be
neglected.

Here, we will study the critical noise needed to create vortices. In
order to address this problem, Eq. (<ref>) is solved
numerically using the same finite-difference scheme as in
Ref. <cit.>. Specifically, we use periodic boundary conditions
for a square of size L_x=L_y=40 with grid step equal to 0.2
corresponding to a cutoff in momentum space equal to K_c=5 π.
The model (<ref>) is a classical field model, that suffers from
the usual ultraviolet catastrophe <cit.>. This implies
that our results will be cutoff dependent. Below, we will discuss
this dependence and how it will affect the comparison of our
theoretical results with experiments.

When analyzing the noise-induced BKT transition, each run starts
from a uniform condensate distribution at D=0. Then we apply noise
with a fixed nonzero intensity D during a time interval t_D. The
detection of vortex pairs in the presence of a strong noise is
somewhat tricky, but fortunately in the absence of noise their
annihilation time is known to be rather long even at very weak
nonequilibrium <cit.> and this time strongly
increases with increasing c <cit.>. For these reasons,
it is more convenient to check the presence of vortex pairs sometime
later after switching off the noise. The used time delay (typically
few our units of time) is sufficient for significant relaxation of
the noisy component in the density and phase distributions of the
condensate and, at the 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 intensity D vortex pairs are present after a noise exposure
time t_D (and hence D>D_ BKT), while for a certain noise
intensity D^'<D no vortex pairs appear even at noise
exposures few times longer then t_D, then D^' lies either
below D_ BKT or above D_ BKT and closer to D_
BKT then to D. Therefore, the critical noise intensity can be
estimated as D_ BKT=D^'± (D-D^'). An
illustration is displayed in Fig. <ref>. For D=0.015 and
noise exposure time t_D=100, a noisy density distribution shown in
Fig. <ref>(a) evolves after switching off the noise into a
pattern with clearly seen vortices and antivortices, which persist
during a relatively long time [see Figs. <ref>(b) and (c)].
For a slightly lower noise intensity, D=0.014, and significantly
longer noise exposure time, t_D=600 , after switching off the
noise the corresponding noisy distribution [Fig. <ref>(d)]
rapidly relaxes towards a uniform vortex-free state
[Figs. <ref>(e) and (f)]. According to our approach, the
estimate for the critical noise in this case is D_
BKT=0.014± 0.001.








§ RESULTS AND DISCUSSION  


In Fig. <ref>(a), the critical noise intensity corresponding to
the BKT transition, D_ BKT, is shown as a function of the
damping parameter κ for the cases of moderate (c=0.3) and
strong (c=4) deviations from equilibrium. In both cases the
critical noise D_ BKT is seen to considerably increase with
κ, despite the fact that the average density of the
condensate at D=D_ BKT is a monotonously decreasing function
of κ [see also  Eq. (<ref>)], this decrease being
especially pronounced at smaller c [see Fig. <ref>(b)].



The inset of Fig. <ref>(b) shows the ratio d_ BKT≡
D_ BKT/n_ BKT as a function of κ. For the two very
different values of the nonequilibrium parameter c, the curves
d_ BKT(κ) appear to lie relatively close to each other.
One can also notice that at c=0.3 the dependence of d_ BKT
on κ is nearly linear (except for the smallest values of
κ).

At κ≫ 1 linear scaling of the noise intensity with
κ directly follows from an additional rescaling of time
[t→ t / κ; see Eq. (<ref>)]. Our numerical
results show that this linear dependence of d_ BKT on
κ remains to good approximation unaltered even when κ
<1.

At equilibrium (c=0), the increase of the critical noise with
increasing κ can be understood by making the connection with
the thermal equilibrium case.  For vanishing nonequilibrium and
κ≫ 1, our equation reduces to model A dynamics
<cit.>, which has a Boltzmann-Gibbs steady state
distribution at temperature  T = D/2κ. In this limit, the
transition occurs at the equilibrium BKT temperature
<cit.>, where T_ BKT = η n, with η a
numerical cutoff dependent constant. For the critical noise, one
then obtains D_ BKT =2 ηκ n.

Numerically, we have observed that for c→ 0, the
critical noise obeys this relation, even for κ < 1. The fact
that the dissipative part of the dynamics does not alter the steady
state can be understood from the following argument. It is well
established that pure Gross-Pitaevskii dynamics (D=κ=c=0)
samples the thermal equilibrium state in the microcanonical
ensemble, at an energy determined by the initial state. Similarly,
the Langevin dynamics [when omitting the i in the l.h.s. of Eq.
(<ref>) and taking c=0] samples the phase space according to
the canonical ensemble at a temperature determined by the balance
between noise and dissipation. Since the thermal state is the steady
state of both the GP and Langevin dynamics, it is natural that the
steady state of the combined dynamics is also at thermal
equilibrium, with the temperature determined by the Langevin part.

As implied by the results displayed in Fig. <ref>(a), the
critical noise intensity D_ BKT increases when moving away
from equilibrium. The dependence of D_ BKT on the
nonequlibrium parameter c is further illustrated in
Fig. <ref>(a) for the cases of  weak (κ=0.01 and 0.1) and
zero damping. At nonzero damping the increase of D_ BKT with
c is partly related to the simultaneous increase of the average
condensate density n_ BKT(c) shown in Fig. <ref>(b). The
latter originates from the growing contribution of the pumping-loss
term in Eq. (<ref>) [the last term in the square brackets in
Eq. (<ref>)]. Indeed, this term tends to keep the condensate
density as close to 1 as possible. As a result, the degrading effect
of damping or noise on the condensate density weakens with
increasing c. This can be seen, e.g., from Fig. <ref>(b) by
comparing to each other the values of n_ BKT at different
c: for each fixed, not too small κ, the value of n_
BKT|_c=4 is significantly larger than n_ BKT|_c=0.3,
even though the density n_ BKT|_c=4 corresponds to a
considerably higher noise intensity than that for n_
BKT|_c=0.3 [see Fig. <ref>(a)].


At the same time, as follows from the behavior of the ratio d_
BKT≡ D_ BKT/n_ BKT [see the inset of
Fig. <ref>(a)], the increase of D_ BKT with c is
considerably faster as compared to that of n_ BKT. This means
that the effect of pumping-loss processes on D_ BKT cannot be
explained solely by the corresponding stabilization of the average
condensate density against noise and damping. This becomes even more
evident when looking at the results for κ=0. Indeed, while
the critical noise D_ BKT demonstrates a clear rise with
increasing c [see Fig. <ref>(a)], the average density n_
BKT(c) remains almost constant [see Fig. <ref>(b)]. At small
c, the values of n_ BKT even slightly decrease with c
suggesting that the stabilizing effects of the pumping-loss
processes do not completely compensate the average-density reduction
caused by the increase of 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, is
that they impede formation of deep local suppressions of the
condensate density. Those deep density suppressions, together with
the appropriate phase gradient, are necessary prerequisites for
vortex pair formation and hence for the BKT transition. The
discussed numerical results imply that just the stabilizing effect
of the pumping-loss term on the local density of the condensate
governs the increase of D_ BKT with c at κ→ 0. At κ = 0 the dependence of D_ BKT on
c becomes close to linear [see Fig. <ref>(a)], while for
d_ BKT(c) a linear function d=ac with a=0.003 provides a
nearly perfect approximation [see Fig. <ref>(b)].

This interpretation of our numerical results can be elucidated by a
linear analysis of the fluctuations. Writing the field as
ψ(x,t)=√(1+δ n(x,t))e^i θ(x,t), one obtains in
linearized approximation for the Fourier components

    ∂/∂ t(
    [ θ_k; n_k ]) =      (
    [        -κϵ_k     -ϵ_k/2-1;        2 ϵ_k - 2c/1+ν+4 κ ]) (
    [ θ_k; n_k ]) 
          +(
    [   √(D)ξ^(θ)_k; √(4 D)ξ^(n)_k ]),

where ϵ_k=k^2. The white noise in Fourrier space has the
correlation function ⟨ξ^(θ)(k,t)
ξ^(θ)(k',t) ⟩ = 4πδ(k+k') δ(t-t'),
analogous for ξ^(n), and ⟨ξ^(θ)ξ^(n)⟩=0.

The steady state density fluctuations can be obtained in closed
form. In the limit of large k, they simplify to



⟨n_k n_k' ⟩= 
   2 πδ(k+k')
4 D (2+κ^2)/κϵ_k     κ≠ 0
   

   2 πδ(k+k') 8 D  (1+ν)/c      κ = 0

For the density fluctuations in real space, we then obtain for a
sufficiently large momentum cutoff K_c
⟨δn^2(x) ⟩∝
    D/κ(2+κ^2) lnK_c     κ≠ 0 
   

    D(1+ν)/c K_c^2      κ = 0 

The local density fluctuations obtained here in the linear
approximation show behavior that is in line with what was observed
numerically for the BKT transition. Density fluctuations are first
suppressed by the damping κ and for κ=0 by the
nonequilibrium parameter c.

Note the very different dependence of both results on the momentum
cutoff. In the presence of a nonzero damping, the cutoff dependence
is a very weak logarithm, where in its absence, it becomes
quadratic. This difference is due to the fact that the momentum
space density tends to a constant in the κ=0 case, where it
decays as k^-2 for κ≠ 0. Even the latter is not fast
enough to ensure convergence in two dimensions, hence the
logarithmic divergence of the density fluctuations.


The values of D_ BKT, obtained from the numerical
simulations, correlate with the above analysis. At c=1.2, ν=1,
and κ=0, an increase of the grid step from 0.2 to 0.4 leads
to a change of the numerically determined D_ BKT from
0.00355± 0.0001 to 0.0139 ± 0.0001 that is in line with the
quadratic dependence of ⟨ n^2(x) ⟩ on K_c in the
absence of energy relaxation [see Eq. (<ref>)]. At the same
time, for nonzero damping (κ≥ 0.01) the relative
difference between the values of D_ BKT corresponding to the
grid step 0.4 and 0.2 is about 40% for κ=0.01 and does not
exceed 15% for κ≥ 0.1. The critical density increases
somewhat when increasing the grid step, such that d_ BKT=
D_ BKT/n_ BKT increases only by about 10% for κ≥ 0.1. Our analytical arguments and those of Ref.
<cit.> suggest a logarithmic dependence on the cutoff,
but we were not able to test this scaling numerically due to a lack
of accessible range in the grid spacing. Choosing 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 down the calculations too much (the maximal
kinetic energy increases quadratically with momentum cutoff).
Moreover, decreasing the grid spacing below 0.2 would give a kinetic
energy that is typically larger than what is required to justify the
lower polariton approximation. 

As seen from Eq. (<ref>), the pumping-loss term increases in
magnitude when decreasing the parameter ν. One can expect,
therefore, that a decrease of ν leads to an increase of the
critical noise D_ BKT. Our simulations confirm this
expectation but, at the same time, show that the influence of ν
on D_ BKT is relatively weak for nonvanishing damping. At
κ=0.01 a decrease of ν by one order of magnitude (from 1
down to 0.1) results in an increase of D_ BKT approximately
by 4% at c=0.3 and by 12% at c=4. The corresponding increase
of d_ BKT is about 20% for both c=0.3 and 4. In other
words, the effect of ν on the critical noise is rather minor as
compared to the much stronger impact of κ and c.







Let us now look in some more detail at the simultaneous 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. <ref>(b)] and, on the
other hand, the linear dependence on c, d_ BKT≈ ac at
κ→ 0, it seems reasonable to consider the
renormalized quantity d_ BKT/(κ+ac). As seen from
Fig. <ref>, all the results of our simulations lie in a rather
narrow interval around 1: 0.6<d_ BKT/(κ+ac)<1.3. This
suggests that the simple expression κ + 0.003c already can
serve as a crude but quite reasonable estimate for d_ BKT at
any κ and c. A better approximation of the numerical
results is obtained (see Fig. <ref>) by using the fitting
function, that summarizes our numerical calculations rather
accurately:

    d^*=       0.609 κ+0.003c 
          +κ^1.411
    c^0.411/2.646κ^0.610  + 0.0706 c^0.889 + 1.622
    κ^1.514 c^1.308

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
meaningful results for inter- and extrapolation [see the inset in
Fig. <ref>(b)].


Finally, we address the question of whether the BKT transition under
consideration can be influenced by the effect of vortex pair
generation by moving vortices <cit.>. This effect has
been predicted to occur in strongly nonequilibrium condensates when
self-accelerated vortices acquire sufficiently high velocities with
respect to the surrounding condensate. The results
of <cit.> were obtained in the absence of noise. In our
numerics containing the noise term, we have observed  that
fluctuations of the condensate density and currents tend to impede
the acceleration of vortices. Consequently, the generation of vortex
pairs in vortex collisions becomes impossible at high noise
intensities. This is illustrated in Fig. <ref>, which shows the
upper boundary D_ m for the noise intensity range, where
generation of new vortex pairs by self-accelerated vortices is
possible, as a function of c at ν=1 and κ=0.01. The
corresponding calculations are performed, like in
 <cit.>, starting with one pinned vortex pair in the
region under consideration and simulating the dynamics during a time
interval Δ t ∼ 1000 after depinning of vortices, now in
the presence of noise.


In Fig. <ref>, the dependence of D_ m on c is given in
comparison with the critical noise D_ BKT at the same
κ and ν. As seen from Fig. <ref>, the region of c
and D, where pair generation by vortices is possible, lies at
 large c (c > 3.3) and has no overlap with the curve
D_ BKT(c). Of course, the latter remains true for any κ≥ 0.01. Indeed, the critical noise D_ BKT increases with
κ. At the same time, damping slows down vortex motion and
thus impedes 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.1 the
velocities of vortices are insufficient for pair generation. So, we
can conclude that at κ≥ 0.01 the BKT transition is not
affected by the processes of pair generation by moving vortices. As
concerns the limit of very weak damping κ< 0.01, the nearly
linear behavior of D_ BKT(c), obtained for κ=0,
manifests no visible peculiarities in the region of large c, where
the processes of pair generation by vortices become possible. This
implies that, also in the limit κ→ 0, the BKT transition
is not considerably affected by these processes.



§ CONCLUSIONS 


We have investigated the critical noise for the spontaneous
formation of unbound vortex-antivortex pairs in a driven-dissipative
bosonic system where particle losses are 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 our model,
the critical noise strength depends – with a suitable choice of
units – on three dimensionless parameters: the energy relaxation
rate (κ), the nonequilibrium parameter (c) and the gain
saturation parameter (ν).

In the absence of energy relaxation, the nonequilibrium parameter
determines the critical noise strength, but this critical value
depends quadratically on the momentum space cutoff, giving our
numerical results limited predictivity for specific experiments. In
the presence of sufficiently strong energy relaxation, the cutoff
dependence is much weaker and experimentally relevant results can be
extracted from the numerics in this case. In this regime, the
critical noise strength increases both with energy relaxation (as
expected from equilibrium calculations) and with nonequilibrium
parameter c. The latter dependence could seem counterintuitive,
because with increasing c, vortices and antivortices repel each
other at large distances, which could favor their unbinding. We
interpret the impeding of vortex-antivortex formation further away
from nonequilibrium as a consequence of the reduction of the density
fluctuations. The effect of nonequilibrium  on the
formation of vortices is therefore opposite to its effect
on the annihilation:  the life time of existing vortices is
dramatically enhanced by nonequilibrium <cit.>.

Quantitatively, we have found that the effect of the nonequilibrium
parameter is actually small for κ≫ 0.003 c, which is
satisfied when damping is not too small and the system is not too
far from equilibrium. Very far from equilibrium, self acceleration
of vortices can lead to the production of new vortex antivortex
pairs. We have shown here that this pair production ceases for
increasing noise.

Our numerical simulations were performed for systems with
periodic boundary conditions whose size is comparable to
the systems employed in current experiments.
The study of the thermodynamical limit of infinite
system size remains a challenge for both theoretical analysis and
experimental investigation.





§ ACKNOWLEDGEMENTS


We thank Sebastian Diehl for a stimulating discussion. This work was
financially suppported by the FWO Odysseus program.



ott16 R. Labouvie, B. Santra, S. Heun, and H. Ott, Phys. Rev. Lett. 116, 235302 (2016).

bose-unit  P. Makotyn, C. E. Klauss, D. L. Goldberger, E. A. Cornell and D. S.  Jin, Nat. Phys. 10, 116 (2014); C. Eigen, J. A. P. Glidden, R. Lopes, E. A. Cornell, R. P. Smith and Z. Hadzibabic
Nature 563, 221 (2018).

fitzpatrick17 M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck, Phys. Rev. X 7, 011016 (2017).

carusotto-ciuti I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).

brenecke13  F. Brennecke, R. Mottl, K. Baumann, R. Landig, T. Donner, and T. Esslinger, Proceedings of the National Academy of
Sciences 110, 11763 (2013).

Kasprzak06
J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun,
J. Keeling, F. Marchetti, M. Szymańska, R. Andre, J. Staehli, V.
Savona, P. B. Littlewood, B. Deveaud and Le Si Dang, Nature 443, 409 (2006).


delteil19
A. Delteil, T. Fink, A. Schade, S. Höfling, C. Schneider, A.
Imamoğlu, Nat. Materials 18, 219 (2019).

rodriguez17
S. R. K. Rodriguez, W. Casteels, F. Storme, N. Carlon Zambon, I.
Sagnes, L. Le Gratiet, E. Galopin, A. Lemaître, A. Amo, C. Ciuti,
and J. Bloch, Phys. Rev. Lett. 118, 247402 (2017).

caputo2016
D. Caputo, D. Ballarini, G. Dagvadorj, C. S. Muñoz,
M. De Giorgi,L. Dominici, K. West, L. N. Pfeiffer, G. Gigli, F. P.
Laussy, M. H. Szymańska and D. Sanvitto, Nat. Mat. 17,145
(2018).

wachtel  G. Wachtel, L. M. Sieberer, S. Diehl and E. Altman, Phys. Rev. B 94, 104520 (2016).

sieberer L. M. Sieberer, G. Wachtel, E. Altman and S. Diehl, Phys. Rev. B 94, 104521 (2016).



squizzato18 D. Squizzato, L. Canet and A. Minguzzi, Phys. Rev. B 97, 195453 (2018).

gladilin14 V. N. Gladilin, K. Ji, and M. Wouters, Phys. Rev. A 90, 023615 (2014).

he15 L. He, L. M. Sieberer, E. Altman, and S. Diehl, Phys.
Rev. B 92, 155307 (2015).

altman15 E. Altman, L. M. Sieberer, L. Chen, S. Diehl, and J.
Toner, Phys. Rev. X 5, 011017 (2015).

sieberer16 L. M. Sieberer, M. Buchhold, and S. Diehl, Reports on
Progress in Physics 79, 096001 (2016).

aranson2002 I. S. Aranson and L. Kramer, Rev. Mod. Phys. 74, 99 (2002).

gladilin18 V. N. Gladilin and M. Wouters, arXiv:1810.06365.

wouters09 M. Wouters and V. Savona, Phys. Rev. B 79, 165302 (2009).

szymanska07 M. H. Szymanśka, J. Keeling, and P. B. Littlewood, Phys. Rev. B 75, 195331 (2007).

sieberer14 L. M. Sieberer, S. D. Huber, E. Altman, S. Diehl, Phys. Rev. B 89, 134310 (2014).

kalinin18 K. Kalinin, P. G. Lagoudakis, N. G. Berloff, Phys. Rev. B  97, 094512 (2018).

comaron18 P. Comaron, G. Dagvadorj, A. Zamora, I. Carusotto, N. P. Proukakis, M. H. Szymańska,  Phys. Rev. Lett. 121, 095302 (2018).

bobrovska17 N. Bobrovska, M. Matuszewski, K. S. Daskalakis, S. A. Maier, S. Kéna-Cohen, ACS Photonics 5, 111 (2017).

estrecho18 E. Estrecho, T. Gao, N. Bobrovska, M. D. Fraser, M. Steger, L. Pfeiffer, K. West, T. C. H. Liew, M. Matuszewski, D. W. Snoke, A. G. Truscott and E. A. Ostrovskaya, Nat. Comm. 9, 2944 (2018).

ohadi18 H. Ohadi, Y. del Valle-Inclan Redondo, A. J. Ramsay, Z. Hatzopoulos, T. C. H. Liew, P. R. Eastham, P. G. Savvidis, J. J. Baumberg, Phys. Rev. B 97, 195109 (2018).

16 J. Keeling and N. G. Berloff, Phys. Rev. Lett. 100,
250401 (2008).

prokofev01 N. Prokof'ev, O. Ruebenacker, and B. Svistunov Phys. Rev. Lett. 87, 270402 (2001).

30 V. N. Gladilin and M. Wouters, New J. Phys. 19, 105005 (2017).

15 M. Wouters and I. Carusotto, Phys. Rev. Lett. 99, 140402
(2007).

a31 N. Bobrovska and M. Matuszewski, Phys. Rev. B 92, 035311
(2015).

38 L.P. Pitaevskii, Zh. Eksp. Teor. Fiz. 35, 408 (1958) [Sov. Phys.
JETP 35, 282 (1959)].

39 M. Wouters, New J. Phys. 14, 075020 (2012).

baboux18 F. Baboux, D. De Bernardis, V. Goblot, V.N. Gladilin, C. Gomez, E. Galopin, L. Le Gratiet, A. Lemaître, I. Sagnes, I. Carusotto, M. Wouters, A. Amo, J. Bloch, Optica 5, 1163 (2018).

bobrovska14 N. Bobrovska, E. A. Ostrovskaya and Michal Matuszewski, Phys. Rev. B 90, 205304 (2014).




chiocchetta13 A. Chiocchetta and I. Carusotto,  EPL, 102, 67007 (2013).


ji15 K. Ji, V. N. Gladilin, and M. Wouters, Phys. Rev. B 91,
045301 (2015).

he17  L. He and S. Diehl,  New J. Phys. 19, 115012 (2017).




kulczy M. Kulczykowski and M. Matuszewski, Phys. Rev. B 95, 075306 (2017).

Comaron17 P. Comaron, G. Dagvadorj, A. Zamora, I. Carusotto, N. P. Proukakis,
M. H. Szymanska,  Phys. Rev. Lett. 121, 095302 (2018).

hohenberg-halperin P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49,435 (1977).











