11institutetext:  Joint Quantum Centre (JQC) Durham-Newcastle, School of
Mathematics, Statistics and Physics, Newcastle University, Newcastle upon
Tyne, NE1 7RU, United Kingdom
Institute of Physics, Polish Academy of Sciences, Al. Lotników 32/46, 02-668
Warsaw, Poland
INO-CNR BEC Center and Università di Trento, via Sommarive 14, I-38123 Povo,
Italy
Department of Physics and Astronomy, University College London, Gower Street,
London, WC1E 6BT, United Kingdom

Polaritons Static properties of condensates Vortices and turbulence

# Non-equilibrium Berezinskii-Kosterlitz-Thouless transition in driven-
dissipative condensates

P. Comaron 1122 I. Carusotto 33 M. H. Szymańska 44 N. P. Proukakis 1111223344

###### Abstract

We study the 2d phase transition of a driven-dissipative system of exciton-
polaritons under non-resonant pumping. Stochastic calculations are used to
investigate the Berezinskii-Kosterlitz-Thouless-like phase diagram for
experimentally realistic parameters, with a special attention to the non-
equilibrium features.

###### pacs:

71.36.+c

###### pacs:

03.75.H

###### pacs:

67.40.Vs

Phase transitions are ubiquitous in nature, both within the classical and
quantum realms. Dimensionality and symmetry are crucial ingredients for the
determination of the types of phase transition (PT) that a given system may
undergo. In a $3d$ system at thermal equilibrium, Bose particles can exhibit
off-diagonal long range order (ODLRO) when driven by a control parameter below
a specific critical temperature. This phenomenon is associated with the
appearance of a Bose-Einstein Condensate (BEC), predicted to occur in both
uniform and confined systems [1]. In $2d$ systems, instead, the presence of
thermal fluctuations destroys ODLRO, compromising the existence of a possible
PT to an ordered state at any finite temperature [2]. Nevertheless, it has
been shown that a different kind of PT to a quasi-condensate state may still
occur, with the decay of correlation functions going from an exponential to a
much slower algebraic law [3, 4]. This Berezinskii-Kosterlitz-Thouless (BKT)
transition can be pictorially understood in terms of the thermally activated
vortices, which change their spatial distribution when crossing the critical
temperature: at high-temperature they proliferate and are freely moving, at
low temperatures they are much less numerous and are bound in pairs, so their
detrimental impact on the coherence gets dramatically suppressed.

The physics becomes even more intriguing when one moves away from isolated
systems to driven-dissipative ones [2, 6, 7], whose stationary state is no
longer determined by thermal equilibrium, but by a non-equilibrium balance of
driving and dissipation. A most celebrated platform to study this physics is
based on exciton-polaritons in semiconductor microcavities, namely bosonic
quasiparticles that arise from the strong coupling between light and matter
excitations. These quasiparticles have a finite lifetime, which calls for some
external pumping to continuously compensate for losses [2]. As in standard
equilibrium BEC, for sufficiently high densities a macroscopic fraction of the
polariton gas condenses into a single momentum state and order develops across
the whole finite sample [8]. In spite of this apparent simplicity, the full
characterization of the PT and of its critical fluctuations in terms of
universality classes is still at the centre of an intense debate, in
particular given their intrinsically $2d$ nature. A strong attention has been
devoted, both experimentally [7, 1, 9, 10, 11] and theoretically [3, 12, 13,
14] to assess up to what point this PT can be described in terms of the
standard BKT theory of equilibrium systems. From the early days of this field,
dramatic consequences of non-equilibrium effects have been highlighted in
polariton systems, from the non-trivial shape of the condensate in real and
momentum spaces [15, 16] to the diffusive Goldstone mode in the collective
excitation spectrum of polariton condensates with small polariton lifetime
[17, 2]. Furthermore, the possibility of breaking BKT algebraic decay of
coherence in the quasi-ordered phase at very large distances has also been
pointed out in [18]. Except for specific cases [7], this occurs however on
length scales well beyond the experimental possibilities. Still, it has been
argued that for realistic system sizes the non-equilibrium character is
responsible for an algebraic decay of the spatial coherence with an exponent
exceeding the upper bound of 0.25 of equilibrium BKT theory [20, 7] and for
the ratio between spatial and temporal correlation exponents being equal to
$2$ [17, 21] instead of $1$ as in the case of an equilibrium-like system [9].

These theoretical predictions suggest that measurements of temporal coherence
are a key ingredient to characterize the nature of the PT: while early works
measured exponential or Gaussian decays of temporal coherence, not compatible
with a BKT transition [20, 22, 23, 24, 9, 25], and possibly related to single-
mode physics [26], power-law decay of temporal correlations have been reported
in recent works with improved samples [9]. Since the long polariton lifetime
in Ref. [9] exceeds other characteristic time scales, one can reasonably
assume the system to be in an equilibrium-like regime [9, 11]. On the other
hand, to date there is no direct numerical or experimental measurement of a
really non-equilibrium regime where the temporal and spatial algebraic
exponents are different.

Motivated by these open questions, in this work we undertake a detailed
numerical study of the PT exhibited by an incoherently pumped (IP) $2d$
polariton fluid under realistic experimental parameters. We numerically
investigate the non-equilibrium steady-state (NESS) phase diagram as a
function of the pump power and we characterize it in terms of the spatial and
temporal correlations, the spectrum of the collective excitation modes and the
spatial distribution of topological defects. Our predictions shine new light
on fundamental properties of the PT and on its non-equilibrium nature.

##  1 Theoretical modelling

We describe the collective dynamics of the polariton fluid through a
generalized stochastic Gross-Pitaevskii equation for the 2d polariton field as
a function of the position $\bm{\mathrm{r}}=(x,y)$ and time $t$, restricting
our investigation here to the simplest case of a spatially homogeneous system
with periodic boundary conditions. The equation describes the effective
dynamics of the incoherently-pumped lower polariton field
$\psi=\psi(\bm{\mathrm{r}},t)$ [27, 2] and includes the complex relaxation
processes by means of a frequency-selective pumping source [28, 29, 3]. The
model, which can be derived from both truncated Wigner (TW) and Keldysh field
theory [2, 17] reads ($\hbar=1$):

$id\psi=\bigg{[}-\frac{\nabla^{2}}{2m}+g|{\psi}|^{2}_{-}+\frac{i}{2}\bigg{(}\frac{P}{1+\frac{|{\psi}|^{2}_{-}}{n_{\mathrm{s}}}}-\gamma\bigg{)}\\\
+\frac{1}{2}\frac{P}{\Omega}\frac{\partial}{\partial t}\bigg{]}\psi{\,dt}+dW$
(1)

where $m$ is the polariton mass, $g$ is the polariton-polariton interaction
strength, $\gamma$ is the polariton loss rate (inverse of the polariton
lifetime), $P$ the strength of the incoherent pumping providing the gain,
$n_{\mathrm{s}}$ is the saturation density, and $\Omega$ sets the
characteristic scale of the frequency-dependence of gain. The renormalized
density
$|{\psi}|^{2}_{-}\equiv\left(\left|{\psi}\right|^{2}-{1}/{{(2dV)}}\right)$
includes the subtraction of the Wigner commutator contribution (where
$dV=a^{2}$ is the volume element of our $2d$ grid of spacing $a$). The zero-
mean white Wiener noise $dW$ fulfils
$\left<dW(\textbf{r},t)dW(\textbf{r}^{\prime},t)\right>=0$,
$\left<dW^{*}(\textbf{r},t)dW(\textbf{r}^{\prime},t)\right>=[(P+\gamma)/2]\delta_{\textbf{r},\textbf{r}^{\prime}}dt,$
where the nonlinear density term is neglected since $|\psi|^{2}/n_{s}\ll 1$.
To describe the physics of the model we start by considering Eq. (1) at a
mean-field (MF) level, i.e. in the absence of the Wiener noise.As widely
discussed in the literature [27, 31], for the case of a frequency-independent
pump ($\Omega=\infty$), a condensate with density
$|\psi^{\mathrm{SS}}|^{2}=n_{\mathrm{s}}\left({P}/{\gamma}-1\right)$ is
expected to appear for pump strengths above threshold
$P>P_{\mathrm{MF}}=\gamma$ and to grow linearly in $P$ with a slope determined
by the saturation density $n_{\mathrm{s}}$. For a frequency-selective pump
($\Omega\neq\infty$), the NESS density loses its linear dependence on $P$, and
takes the slightly more complicated form
$|\psi^{\mathrm{SS}}|^{2}=n_{\mathrm{s}}\left[{P}/\left(Pg|\psi^{\mathrm{SS}}|^{2}/\Omega+\gamma\right)-1\right]$
[31].

The effect of small excitations around the bare condensate steady-state
solution can be described by means of the linearized Bogoliubov approximation
[1, 2]. By linearizing the deterministic part of Eq. (1) around the steady
state solution
$\psi(\textbf{r},t)=\psi^{\mathrm{SS}}+\delta\psi(\textbf{r},t)e^{-i\omega t}$
we obtain a pair of coupled Bogoliubov equations for the field
$\delta\psi(\textbf{r},t)$ and its complex conjugate
$\delta\psi^{*}(\textbf{r},t)$. Thanks to translational invariance, the
different k-modes are decoupled, so we can move to Fourier space and define a
k-dependent Bogoliubov matrix $\mathcal{L}_{\textbf{k}}$ [31, 29],

$\mathcal{L}_{\textbf{k}}=\left(\begin{array}[]{cc}\Lambda({\epsilon_{\textbf{k}}}+\mu-i\Gamma)&\Lambda(\mu-i\Gamma)\\\
\Lambda^{*}(-\mu-i\Gamma)&\Lambda^{*}(-{\epsilon_{\textbf{k}}}-\mu-i\Gamma)\\\
\end{array}\right)$ (2)

with $\Gamma={\gamma}\left(P-\gamma\right)/{2P}$, the free-particle dispersion
$\epsilon_{\textbf{k}}=k^{2}/2m$, the interaction energy
$\mu=g|\psi^{\mathrm{SS}}|^{2}$ and
$\Lambda=\left(\gamma_{a}+i\gamma_{b}\right)$ with
$\gamma_{a}={1}/{[1+\left({P}/{(2\Omega)}\right)^{2}]}$ and
$\gamma_{b}=-{P\gamma_{a}}/{2\Omega}$. The diagonalization of
$\mathcal{L}_{\textbf{k}}$ eventually leads to the double-branched excitation
spectrum

$\begin{split}\omega_{\textbf{{k}}}^{\pm}&=-i\left[{\gamma_{a}\Gamma}-{\gamma_{b}}({\epsilon_{\textbf{{k}}}}+\mu)\right]\pm\\\
&\sqrt{\Gamma^{2}{\gamma_{a}}^{2}+{\gamma_{b}}^{2}\mu^{2}-2\Gamma\gamma_{a}{\gamma_{b}}\left(\epsilon_{\textbf{{k}}}+\mu\right)-{\gamma_{a}}^{2}\epsilon_{\textbf{{k}}}\left(\epsilon_{\textbf{{k}}}+2\mu\right)}.\end{split}$
(3)

At high momenta $\mathbf{k}$ this spectrum recovers a single-particle
behaviour with parabolic dispersion, while the frequency-dependence of pumping
results in an increasing linewidth for growing $\mathbf{k}$. For small
$\textbf{k}\to 0$, the Goldstone mode describing long-wavelength twists of the
condensate phase and associated to the spontaneously broken U(1) symmetry
exhibits the diffusive behaviour typical of driven-dissipative systems [2],
rather than the sonic one characteristic of their equilibrium counterpart [1].

This physics is illustrated in Fig. 1, where the prediction of Eq. (3) is
plotted for increasing values of the pump strength $P$. The value of the
critical momentum

$k_{\mathrm{c}}=\sqrt{{2m}\left[-\frac{\Gamma{\gamma_{b}}}{{\gamma_{a}}}-\mu+\frac{\sqrt{\left(\Gamma^{2}+\mu^{2}\right)\left({\gamma_{a}}^{2}+{\gamma_{b}}^{2}\right)}}{{\gamma_{a}}}\right]}$
(4)

separating the diffusive behaviour from the sonic one at higher $\mathbf{k}$
increases as the system moves away from the threshold point $P_{\mathrm{MF}}$.

[width=]./fig1.pdf

Figure 1:  Real part (blue) and imaginary part (red) of the Bogoliubov
excitation spectrum (3) calculated for the parameters of the case IPΩ=50
(whose density is plotted in Fig. 2 as a blue curve), but where the relative
pump strength takes now the values $P/P_{\rm MF}~{}=~{}2,6,10,14,18,22$,
increased as indicated by the colour gradients. Corresponding real part for
$P/P_{\rm MF}~{}=~{}1.06$ is shown in Fig. 5(iv) as a blue curve.

##  2 The non-equilibrium Berezinskii-Kosterlitz-Thouless Phase diagram

We simulate the system dynamics by numerically integrating in time the
stochastic differential equations for the polariton field shown in (1);
numerical details are reported in [32]. In Fig. 2 we show the typical driven-
dissipative BKT PT-diagram of an incoherently-pumped polariton condensate, in
which the different observables are shown as a function of the pump strength
$P$. This is characterized by two distinct phases: a) a disordered phase
displaying a low density of polaritons, an exponential decay of spatial
correlations and a plasma of unbound free vortices; b) a superfluid phase
displaying a significant density of polaritons, an algebraic decay of spatial
correlations and a low density of vortices, mostly bound in vortex-antivortex
pairs [3, 7, 9].

Our first step in the investigation of the IP polariton PT was to clarify the
impact of fluctuations introduced by the stochastic noise on the average
density. The mean-field (stochastic) density
$|\psi|^{2}=\int|\psi({r})|^{2}d{r}/L_{x}L_{y}$
[$|\psi|^{2}=|{\psi}|^{2}_{-}$] is calculated by evolving Eq. (1) without
(with) the contribution of the Wiener noise. These two curves are plotted in
the inset of Fig. 2 as dashed black and solid blue curves, respectively.
Contrary to the mean-field case where $|\psi_{\mathrm{MF}}|=0$ in the
disordered phase $P<P_{\mathrm{MF}}$, within the stochastic framework the
density field is always non-zero, independently of the value of the pump
$P$111In the disordered phase, fluctuations are responsible for building up a
small but not negligible density of incoherent polaritons, the only zero-
density point coinciding with a vanishing pump strength $P=0$. In the quasi-
ordered phase the density grows considerably and asymptotically approaches the
mean-field prediction..

In contrast to the clear threshold shown by the MF curve, the smooth increase
of the density with pump strength shown by the stochastic theory requires a
more involved determination of the critical point. As done in previous works
[7] – and discussed in subsequent sections – our procedure to precisely
determine the critical point involves the functional form of the decay of
correlation functions, the behaviour of vortices in the vicinity of the
criticality and the appearance of the diffusive Goldstone mode in the
spectrum. Interestingly, we note in Fig. 2 that fluctuations are responsible
for an upward shift of the critical point $P_{\mathrm{BKT}}$ (vertical blue
line) with respect to the MF value $P_{\mathrm{MF}}$ (vertical gray line). In
order to unravel the dependence of $P_{\mathrm{BKT}}$ on the physical
parameters $n_{\mathrm{s}}$ and $\Omega$, this figure shows the phase diagram
for three different choices of parameters, listed in the caption of the
figure. For each case analysed, the critical point is highlighted with a
vertical coloured thick line. As general trends, we find that stronger
fluctuations in higher modes ($\Omega\rightarrow\infty$) and smaller
saturation densities ($n_{\mathrm{s}}\rightarrow 0$) lead to a larger shift of
$P_{\mathrm{BKT}}$ with respect to the mean-field $P_{\mathrm{MF}}$.

This feature can be understood by fixing one of the two parameters and
focusing on the other. On the one hand, for a frequency-independent pump
($\Omega=\infty$, green and violet lines), we note that increasing
$n_{\mathrm{s}}$ makes the BKT threshold $P_{\mathrm{BKT}}$ to shift closer to
$P_{\mathrm{MF}}$: the slope of the total density increases with
$n_{\mathrm{s}}$, so the critical density is reached at lower values of the
pump strength. On the other hand, for a fixed value of
$n_{\mathrm{s}}=500\mathrm{\mu m^{-2}}$ (blue and violet curves), the presence
of a frequency-selective pump leads to an effective thermal population of less
field modes. As a consequence, a weaker pump is sufficient to concentrate a
macroscopic population in the lowest modes, which has the effect of shifting
the threshold point back towards the mean-field value $P_{\rm MF}$.

As expected in a BKT-like picture, the IP phase transition can be pictorially
understood as being mediated by the unbinding of vortex-anti-vortex pairs into
a plasma of free vortices [7, 18]. In Fig. 2 the NESS average number of
topological defects $\left<N_{\mathrm{v}}\right>$ is plotted as a thick red
line for the parameters of the $\mathrm{IP_{\Omega=50}}$ case. Details on the
procedure we adopt to extract $N_{\mathrm{v}}$ are reported in [32] as well as
the illustrations of three exemplary configurations of vortices across the BKT
phase diagram. The low-pump disordered phase is characterized by a large
number of free vortices which are free to proliferate. As the pump power is
increased and approaches the threshold point, the number of vortices
$\left<N_{\mathrm{v}}\right>\propto\xi^{-1/2}$ starts to decay as expected for
a continuum PT with diverging correlation length $\xi$, a detailed study of
which is presented in Ref. [6]. At the critical point, around which the
process of vortices pairing starts to take place, the average number of
vortices is still non-zero but these are mostly grouped in vortex-antivortex
pairs. Due to vortex binding and annihilation processes,
$\left<N_{\mathrm{v}}(P)\right>$ shows a severe drop right above the critical
point and rapidly decreases to zero. Deep in the quasi-condensed phase when
$P\gg P_{\mathrm{BKT}}$, as the stochastic density grows and onset of
coherence appears, the dynamical annihilation processes are severe and
eventually leave the system free of defects.

##  3 Spatial-temporal coherence and the critical region

In this section we investigate the long-distance, late-time decay of the
spatial and temporal first-order correlation functions, $g^{(1)}(\Delta r)$
and $g^{(1)}(\Delta t)$ respectively. We focus here on the results for the
IPΩ=50 case; the complementary study for the
IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$ case with a frequency-independent
pump is illustrated in [32]. Within our semi-classical model, the spatial and
temporal two-point first order correlation functions are defined,
respectively, as

$\displaystyle g^{(1)}(\Delta r)=\frac{\langle\psi^{*}({r_{0}}+\Delta
r,t)\psi({r_{0}},t)\rangle}{\sqrt{\langle\left|\psi({r_{0}}+\Delta{r},t)\right|^{2}\rangle\langle\left|\psi({r_{0}},t)\right|^{2}\rangle}}\;,$
(5) $\displaystyle g^{(1)}(\Delta
t)=\frac{\langle\psi^{*}({r_{c}},t_{0})\psi({r_{c}},t_{0}+\Delta
t)\rangle}{\sqrt{\langle\left|\psi({r_{c}},t_{0})\right|^{2}\rangle\langle\left|\psi({r_{c}},t_{0}+\Delta
t)\right|^{2}\rangle}}\;,$ (6)

and are calculated at a sufficiently late time $t=t_{\mathrm{SS}}$ at which
the system has reached its NESS, and with $r_{c}=(L_{x}/2,L_{y}/2)$ being the
central point of the spatial grid. The numerical results for the spatial and
temporal correlations are illustrated in panels a) and b) of Fig. 3,
respectively.

Inspired by earlier works [7, 34], we characterize the behaviour of the
steady-state correlation functions as function of the pump strength $P$. In
Fig. 3 we show the transition from an exponential decay $g^{(1)}\sim
e^{-r/\xi}$ in the disordered phase, to a power-law decay $g^{(1)}\sim
r^{-\alpha}$ in the quasi-ordered phase, as expected for the spatial
correlation function of an equilibrium BKT transition. The same behaviour is
found for the temporal correlation function222Note that in our simulations the
accessible time duration are not long enough to observe the finite-size-
induced Schawlow-Townes decay [35, 26].. In order to identify whether a given
correlation function is characterised by either exponential or algebraic
decay, we have fitted each curve with both functions, paying particular
attention to ensure that all computational results are correctly converged
within the spatial and temporal windows selected for the fitting procedure
[32]. We have then calculated the Root-mean-square deviation (RMSD) of the
residuals of the fits within the fitting window selected and we have selected
the fit that minimizes the RMSD [32]. In Fig. 3 we superimpose on top of each
correlation function $g^{(1)}$, the most accurate fitting curve, represented
by red or blue dashed lines in the exponential or power-law cases,
respectively.

[width=]./fig2.pdf

Figure 2:  Non-equilibrium steady-state phase diagram showing mean-field (MF)
and averaged stochastic density (sGPE) (dashed and solid coloured curves,
respectively) in logarithmic-linear scale. For each set of parameters, we
associate a colour. Blue (labelled as IPΩ=50)$:\Omega=50\gamma=11.09ps^{-1}$
and $n_{\mathrm{s}}=500\mathrm{\mu m^{-2}}$. Green (labelled as
IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$): $\Omega=\infty$ and
$n_{\mathrm{s}}=1500\mathrm{\mu m^{-2}}$. Violet (labelled as IPΩ=∞):
$\Omega=\infty$ and $n_{\mathrm{s}}=500\mathrm{\mu m^{-2}}$. For each set of
parameters, the BKT threshold is shown as a vertical coloured line. The
vertical grey line shows the mean-field threshold $P_{\mathrm{MF}}=1$. The
average number of vortices $\left<N_{\mathrm{v}}\right>$ for the IPΩ=50 case
is depicted as a red curve. The inset shows comparison between the mean-field
(dashed black line) and the averaged stochastic density (blue solid line) for
the IP(Ω=50) case, plotted in linear-linear scale.

Figure 3:  Crossover from exponential to algebraic decay of the spatial (a)
and temporal (b) correlation functions, defined as in Eqs. (5) and (6). Thick
dashed red (blue) curves correspond to exponential (power-law) fitting, from
which the values of the correlation length $\xi$ and of the power-law
exponents $\alpha_{\mathrm{s}}$ and $\alpha_{\mathrm{t}}$ plotted in Fig. 8
were extracted. For each curve, we superimpose only the best fitting option.
Both fits are only shown for the curves which lie in the critical region. The
fits are restricted to the chosen fitting window, indicated by the gray
shadow.

Figure 4:  a) For the IPΩ=50 case, plot of the quantities
$\sigma_{\mathrm{s}}$ (solid red curve) and $\sigma_{\mathrm{t}}$ (dashed blue
line) defined in Eqs. (7), which identify the critical region (shaded blue
region) and the critical point $P_{\mathrm{BKT}}$ (vertical red line).
Squares, diamonds and circles correspond to points which fall before, within
and above the critical region, respectively. b) Log-linear plot of the spatial
($\alpha_{\mathrm{s}}$, filled circles) and temporal ($\alpha_{\mathrm{t}}$,
empty circles) exponents, extracted from the power-laws fits of Fig. 3 with
error bars. c) The correlation length $\xi$ (green solid curve) diverges when
approaching $P_{\mathrm{BKT}}$ from the disordered phase. Away from the
critical point into the quasi-ordered phase, the decay of the spatial
(temporal) algebraic exponent $\alpha_{\mathrm{s}}$ ($\alpha_{\mathrm{t}}$) is
shown as empty (filled) circles and red (blue) solid line. The excitations
spectrum for three characteristic values of the pump strength indicated with
i), ii) and iii) is shown in Fig. 5. The inset shows a plot of $\xi$ and
$\alpha_{\mathrm{s}}$ for the IPΩ=50 and
IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$ cases. In all panels above, the
colour of the markers corresponds to the one of the different curves in Fig.
3.  Figure 5:  Panels i-iii) show colorplots of the spectra obtained
numerically from the Fourier transform of $|\psi{(\textbf{r},t)}|^{2}$, for
the three points i) $P=1.03$, ii) $P=1.0338$ and iii) $P=1.06$ highlighted in
Fig. 8. Panel iv) shows a comparison with the real part of the analytical
dispersion in (3) (blue curves) for case iii). The inset shows a different
numerical simulation with a ($100\times$) larger box and spatial
discretization, able to capture the low-$k$ region and the diffusive branch of
the spectrum.

Fig. 8(a) shows how we characterise the critical region by means of the RMSD
ratios of the fits of the spatial (red solid curve) and temporal (dashed blue
curve) correlators, namely

$\sigma_{\mathrm{s}}=\frac{\mathrm{RMSD}^{\mathrm{pow}}_{\mathrm{s}}}{\mathrm{RMSD}^{\mathrm{exp}}_{\mathrm{s}}},\qquad\sigma_{\mathrm{t}}=\frac{\mathrm{RMSD}^{\mathrm{pow}}_{\mathrm{t}}}{\mathrm{RMSD}^{\mathrm{exp}}_{\mathrm{t}}}\,.$
(7)

By visually comparing the residuals of the exponential and power-law fits on
this figure, one can infer the position of the critical point as the point
where the two curves go through $1$. This point indicates the exponential-to-
power-law transition, which takes place for the same
$P_{\mathrm{BKT}}~{}\sim~{}1.0325$ (vertical red solid line in Fig. 8) for
both spatial and temporal correlation functions333 In both Fig. 3 and Fig.
8(a), there exists an intermediate regime in an interval of
$P_{\mathrm{BKT}}$, where the curves are neither exactly fitted by a power-law
or an exponential form. Therefore, we refer to the critical region as the
portion of the phase diagram located in between the last correlator showing
“clear” exponential decay (lower bound) and the first exhibiting “clear”
power-law decay (upper bound). In our case, these lower and upper bounds are
located at $P=1.032$ and $P=1.0338$, respectively. In both Figs. 2 and 8, the
critical region is highlighted with a blue shading..

This analysis of the decay of the correlation functions allows us to extract
quantities which are strictly linked to the physical nature of the PT. Namely,
the correlation length $\xi$, extracted from the exponential fit in the
disordered phase, and the algebraic exponents $\alpha_{s}$, $\alpha_{t}$ which
quantify the algebraic decay of space and time correlators in the quasi-
ordered one. In equilibrium systems, the former is known to be related to the
superfluid density [4]. These quantities are plotted in Fig. 8(c), as a
function of the pump strength $P$ and represented as solid green [$\xi(P)$],
solid red [$\alpha_{\mathrm{s}}(P)$] and blue thick [$\alpha_{\mathrm{t}}(P)$]
curves. Markers in Fig. 8 are coloured in a way to match the ones of Figs. 3.
In the disordered phase [left part of Fig. 8(c)] the coherence length $\xi(P)$
diverges when approaching the critical region from the left, as expected for a
continuum PT (in finite systems) undergoing critical slowing down444While the
dataset extracted is suitable for a qualitative description of the PT, a
possible quantitative extraction of critical exponents would require a more
advanced scaling analysis with larger sample sizes.. In the quasi-ordered
phase [right part of Fig. 8(b-c)], the exponents $\alpha_{\mathrm{s}}$ and
$\alpha_{\mathrm{t}}$ show a decreasing behaviour as the control parameter $P$
is increased, which is connected to the expected onset of coherence.

In Fig. 5 we plot three exemplary cases of excitation spectrum calculated from
the spatio-temporal Fourier Transform (FT) of $|\psi{(\textbf{r},t)}|^{2}$
across the PT. As expected, the system moves from a free-particle quadratic
dispersion below the transition [Fig. 5(i)] to a non-equilibrium spectrum, as
in (3), above the transition [Fig. 5(iii-iv)]. We find the analytical
Bogoliubov dispersion (3) to correctly describe our numerics: the agreement
between the peak of the numerical spectrum (colour map) and the analytical
prediction (blue curves) is explicitly illustrated for the last case in Fig.
5(iv). For this case, we find that the critical momentum
$k_{c}(P=1.06)=1.26\times 10^{-2}\mathrm{\mu m^{-1}}$ is on the order of the
momentum discretization $\Delta k=\pi/L=1.06\times 10^{-2}\mathrm{\mu m^{-1}}$
of the numerical simulation; as a consequence, in the main panels the
diffusive branch is hidden by the sonic behaviour of the dispersion at
$k>k_{c}$. The numerical value $k_{c}^{\textrm{num}}=1.0(5)\times
10^{-2}\mathrm{\mu m^{-1}}$ is extracted by simulating a system with a
($100\times$) smaller $\Delta k$ [32]: the low-$k$ part of the spectrum,
plotted as an inset of Fig. 5(iv), is now visible and in good agreement with
the analytically predicted curve.

##  4 Discussion on the nature of the phase transition

Previous experimental [20] and theoretical [7] works showed that a spatial
power-law exponent exceeding the $\alpha_{\mathrm{s}}=0.25$ upper bound of the
equilibrium theory is a signature of the non-equilibrium nature of the PT.
While this is evidently the case for the
IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$ simulations, the numerically
obtained value of the exponent in the IPΩ=50 case never exceeds the
equilibrium upper bound. However, by enlarging the system size we find that
$g^{(1)}({\Delta}r)$ is converged in space over all the quasi-ordered pump
range except for the extreme point $P=1.0338$. This is expected, as in the
very vicinity of criticality finite size effects can be most important. As
shown in [32], by enlarging the box by 1.5 and 2 times, power-low exponents
are found to lie within the interval $0.25<\alpha_{s}<0.35$ [reported in Fig.
8(b,c) as a large errorbar for the brown point
$\alpha_{\mathrm{s}}(P=1.0338)$]. This confirms the argument that for non-
equilibrium driven-dissipative system $\alpha_{\mathrm{s}}$ can exceed the
upper equilibrium limit of $\alpha=0.25$ in the critical region, for both
frequency-independent and frequency-dependent pumping.

A key difference between equilibrium and non-equilibrium PTs is encoded in the
relation between the $\alpha_{\mathrm{s}}$ and $\alpha_{\mathrm{t}}$
exponents. In the equilibrium case, the sonic nature of the dispersion leads
to $\alpha_{\mathrm{s}}=\alpha_{\mathrm{t}}$. For a non-equilibrium driven-
dissipative condensate, the diffusive nature of the Goldstone mode suggests
instead that $\alpha_{\mathrm{s}}\sim 2\alpha_{\mathrm{t}}$ [17]. At first
sight, the prominent sonic branch visible in the spectrum of Fig. 5(iii) could
suggest that we are in a similar equilibrium-like scenario as in Ref. [11],
where almost equal values were measured for $\alpha_{\mathrm{s}}$ and
$\alpha_{\mathrm{t}}$, in strong contrast to our numerics. Looking at the
excitation spectrum in Refs. [9, 11] reveals that the critical momentum
$k_{c}(P=1.06)$ is there $2.53\times 10^{2}$ times smaller than the one
considered here, giving a characteristic length $2\pi/k_{c}$ that largely
exceeds the system size. This is due to the much longer lifetime displayed by
polaritons in those experiments and is responsible for the absence of an
observable diffusive region in the Goldstone mode. Our numerical study shows
instead a power-law decay of both spatial and temporal correlation function,
with an exponent ratio $\alpha_{\mathrm{s}}\sim 2\alpha_{\mathrm{t}}$ [Fig.
8(b)], suggesting a non-equilibrium nature of the condensate. However, due to
the inability to numerically simulate a large enough box to clearly highlight
the diffusive Goldstone mode, we cannot determine whether the different values
measured for $\alpha_{\mathrm{s}}$ and $\alpha_{\mathrm{t}}$ is due to its
non-equilibrium nature or finite-size effects, or an interplay between the
two.

##  5 Conclusions

In this paper we have undertaken a detailed numerical analysis to investigate
the non-equilibrium phase transition displayed by a polariton system under
incoherent pumping. We have characterized the non-equilibrium phase diagram
within both mean-field and stochastic pictures, confirming for realistic
system sizes a BKT-like scenario for non-equilibrium condensates featuring a
crossover between binding/unbinding of vortices and between an
exponential/power-law decay of correlations. Particular attention was given to
the role of fluctuations in the shift of the critical point with respect to
the mean-field picture and to the long-distance and late-time decay of the
spatial and temporal correlation functions. Our findings show that the non-
equilibrium driven-dissipative phase transition exhibits an algebraic exponent
exceeding the upper-bound equilibrium limit of $1/4$ in agreement with
previous experimental [20] and theoretical [7] works. A non-equilibrium nature
of the condensate is also suggested by the ratio
$\alpha_{\mathrm{s}}/\alpha_{\mathrm{t}}\sim 2$ of the algebraic decay
exponents of space and time correlators, extracted by our numerical
simulations and suggested by analytical calculations within a Keldysh
framework [17]. We note however that such an effect could be also due to a
possible interplay with finite-size effects. It would be of a great interest
to explore the interplay between non-equilibrium and finite-size effects in
spatial correlations in future works. Our results suggest that a complete
characterization of the non-equilibrium Berezinskii-Kosterlitz-Thouless phase
transition is within current experimental reach using polariton fluids.

###### Acknowledgements.

We acknowledge financial support from the Quantera ERA-NET cofund projects
NAQUAS (through the Engineering and Physical Science Research Council (EPSRC),
Grant No. EP/R043434/1 (P.C. and N.P.P.)), and InterPol (through the National
Science Center, Poland, Grant No. 366 2017/25/Z/ST3/03032 (P.C.) and the EPSRC
Grant No. 368 EP/R04399X/1 (M.H.S.)). We acknowledge EPSRC Grant No.
EP/K003623/2 (M.H.S.) and support from the Provincia Autonoma di Trento and
from the Quantum Flagship Grant PhoQuS (820392) of the European Union (I.C.).

## References

  * [1] Pitaevskii L. Stringari S. Bose-Einstein Condensation (Clarendon Press) 2003.
  * [2] Mermin N. D. Wagner H. Phys. Rev. Lett.1719661133.
  * [3] Berezinskii V. L. Sov. Phys. JETP321971493.
  * [4] Kosterlitz J. M. Thouless D. J. Journal of Physics C: Solid State Physics619731181.
  * [5] Carusotto I. Ciuti C. Reviews of Modern Physics852013299.
  * [6] Proukakis N. P., Snoke D. W. Littlewood P. B. (Eds.) Universal Themes of Bose-Einstein Condensation (Cambridge University Press) 2017\.
  * [7] Hohenberg P. C. Halperin B. I. Reviews of Modern Physics491977435.
  * [8] Kasprzak J., Richard M., Kundermann S., Baas A., Jeambrun P., Keeling J., Marchetti F., Szymańska M., Andre R., Staehli J. et al. Nature4432006409.
  * [9] Caputo D., Ballarini D., Dagvadorj G., Sánchez Muñoz C., De Giorgi M., Dominici L., West K., Pfeiffer L. N., Gigli G., Laussy F. P., Szymańska M. H. Sanvitto D. Nature Materials172017145 EP .
  * [10] Stepanov P., Amelio I., Rousset J.-G., Bloch J., Lemaître A., Amo A., Minguzzi A., Carusotto I. Richard M. Nature communications1020191.
  * [11] Ballarini D., Caputo D., Dagvadorj G., Juggins R., De Giorgi M., Dominici L., West K., Pfeiffer L. N., Gigli G., Szymańska M. H. et al. Nature communications1120201.
  * [12] Proukakis N. P., Gardiner S. A., Davis M. J. Szymańska M. H. (Eds.) Quantum Gases: Finite Temperature and Non-Equilibrium Dynamics Vol. 1 (World Scientific) 2013.
  * [13] Gladilin V. N. Wouters M. Phys. Rev. B1002019214506.
  * [14] Mei Q., Ji K. Wouters M. arxiv:2002.018062019.
  * [15] Richard M., Kasprzak J., Romestain R., André R. Dang L. S. Phys. Rev. Lett.942005187401.
  * [16] Wouters M., Carusotto I. Ciuti C. Phys. Rev. B772008115340.
  * [17] Szymańska M., Keeling J. Littlewood P. Physical review letters962006230602.
  * [18] Altman E., Sieberer L. M., Chen L., Diehl S. Toner J. Phys. Rev. X52015011017.
  * [19] Dagvadorj G., Fellows J., Matyjaśkiewicz S., Marchetti F., Carusotto I. Szymańska M. Physical Review X52015041028.
  * [20] Roumpos G., Lohse M., Nitsche W. H., Keeling J., Szymańska M. H., Littlewood P. B., Löffler A., Höfling S., Worschech L., Forchel A. Yamamoto Y. PNAS10920126467.
  * [21] Szymańska M., Keeling J. Littlewood P. Physical Review B752007195331.
  * [22] Krizhanovskii D. N., Sanvitto D., Love A. P. D., Skolnick M. S., Whittaker D. M. Roberts J. S. Phys. Rev. Lett.972006097402.
  * [23] Love A. P. D., Krizhanovskii D. N., Whittaker D. M., Bouchekioua R., Sanvitto D., Rizeiqi S. A., Bradley R., Skolnick M. S., Eastham P. R., André R. Dang L. S. Phys. Rev. Lett.1012008067404.
  * [24] Kim S., Zhang B., Wang Z., Fischer J., Brodbeck S., Kamp M., Schneider C., Höfling S. Deng H. Phys. Rev. X62016011026.
  * [25] Askitopoulos A., Pickup L., Alyatkin S., Zasedatelev A., Lagoudakis K., Langbein W. Lagoudakis P. G. arXiv:1911.089812019.
  * [26] Amelio I. Carusotto I. Phys. Rev. X102020041060.
  * [27] Wouters M. Carusotto I. Phys. Rev. Lett.992007140402.
  * [28] Wouters M., Liew T. C. H. Savona V. Phys. Rev. B822010245315.
  * [29] Chiocchetta A. Carusotto I. EPL102201367007.
  * [30] Comaron P., Dagvadorj G., Zamora A., Carusotto I., Proukakis N. P. Szymańska M. H. Phys. Rev. Lett.1212018095302.
  * [31] Wouters M. Carusotto I. Phys. Rev. Lett.1052010020602.
  * [32] See [link] for numerical details and methods.
  * [33] Zamora A., Dagvadorj G., Comaron P., Carusotto I., Proukakis N. P. Szymańska M. H. Phys. Rev. Lett.1252020095301.
  * [34] Comaron P., Larcher F., Dalfovo F. Proukakis N. P. Phys. Rev. A1002019033618.
  * [35] Keeling J., Szymańska M. H. Littlewood P. B. Keldysh Green’s function approach to coherence in a non-equilibrium steady state (Springer) 2010 pp. 293–329.

Supplementary Material for: Non-equilibrium Berezinskii-Kosterlitz-Thouless
transition in driven-dissipative condensates

In this Supplementary Material we present technical details related to our
numerical procedure and analysis for the study of a non-equilibrium driven-
dissipative phase transition. We discuss important numerical details,
numerical methods and convergence of spatial and temporal correlations and we
explain numerical methods for topological defect counting. Finally, show a
detailed analysis on the long-range behaviour close to the critical point.

##  6 Numerical calculations

We simulate the system dynamics by numerically integrating in time the
stochastic differential equations for the polariton field shown in Eq. (1) of
the main text. The numerical integration is performed on a 2d-grid with
periodic boundary conditions. The lattice is composed of $301^{2}$ grid
points, with total lengths $L_{\mathrm{x}}=L_{\mathrm{y}}=295.11\mu m$ and
lattice spacing $a=0.98\mu m$. Notice that the lattice spacing $a$, both
introduces a cut-off $\propto a^{-1}$ in the momentum representation of the
field, and is chosen between a lower bound given by the macroscopic scale of
the system, such as the healing length, and the upper bound $a^{2}\gg
g/\gamma$ given by the validity of the truncated Wigner methods used for the
description of the stochastic field equations [1, 2, 3]. The time dynamics of
the polariton field is here performed by integrating Eq. (1) in time with the
XMDS2 software framework [4]. Here we have used a fourth-order Runge-Kutta
algorithm with fixed time-step of $4.5\cdot 10^{-2}\mathrm{ps}$, which ensures
stochastic noise consistency. All the results expressed in the present work
are converged with respect to the number of stochastic realisations
$\mathcal{N}=200$.

In the present study, we use typical experimental parameters [5]: lifetime
$\tau=1/\gamma=4.5\mathrm{ps}$, $m=6.2\ 10^{-5}\ m_{\mathrm{e}}$, $g=6.82\
10^{-3}\ \mathrm{meV\mu m^{2}}$. In Fig. 1 of the main text the non-
equilibrium steady-state properties of the incoherently pumped polariton
system are discussed in term of three different sets of the values ($\Omega$,
$n_{\mathrm{s}}$). Such sets are labelled as IPΩ=50
($\Omega=50\gamma=11.09ps^{-1}$ and $n_{\mathrm{s}}=500\mathrm{\mu m^{-2}}$);
IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$ ($\Omega=\infty$ and
$n_{\mathrm{s}}=1500\mathrm{\mu m^{-2}}$) and IPΩ=∞ ($\Omega=\infty$ and
$n_{\mathrm{s}}=500\mathrm{\mu m^{-2}}$). The analysis on the coherence
properties in the main text is mainly centred around the case IPΩ=50 while the
complementary analysis for the case
IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$ is presented in this supplemental
document. For the sake of clarity, it is also worth noting that these sets of
parameters have been employed in other works of the authors. Specifically, in
Ref. Comaron et. al [3] the set of parameters IPΩ=50 and
IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$ have been used for the
investigation of late-time quench dynamics properties of driven-dissipative
condensate, while in Zamora et. al [6], the validity of the Kibble-Zurek
hypothesis for microcavity fluids is discussed; in the latter only the case
IPΩ=50 is considered.

The steady state values of the observables discussed in the analysis of the
main paper have been reached after evolving the MF (sGPE) model from a very
small initial value of the density (from a Gaussian white noise) to a steady
state with a non-zero coherent density, thereby simulating all cases starting
from a state without a condensate.

##  7 First-order spatial and temporal correlations: methods and numerical
convergence

In Fig. 6 we repropose the data shown in Fig. 3 of the main text, plotted both
in log-log (left) and log-linear (right) scale.

Figure 6: Spatial (top row) and temporal (bottom row) correlations, as shown
in Fig. (2) of the main text, both in log-log (left) and linear-log scale
(right).  Figure 7: Binding/unbinding crossover of topological defects. The
three panels depict the phase ${\vartheta}=\arg{\\{\psi\\}}$ for (left)
$P=0.95$ with $\left<N_{\mathrm{v}}\right>\sim 1.3\times 10^{3}$, (centre)
$P=1.0338$ with $\left<N_{\mathrm{v}}\right>\sim 76$ and (right) $P=1.06$ with
$\left<N_{\mathrm{v}}\right>\sim 0$. Vortices and anti-vortices are reported
as red dots and green triangles, respectively.

Figure 8: Spatial correlations for larger system sizes. Top panel: spatial
correlations for $P=1.0340$, immediately outside the critical region in the
quasi-ordered phase, calculated for different lattice size
$L=295.11\mathrm{\mu m}$ (red points) and $L=444.42\mathrm{\mu m}$ (blue
points), for the case IPΩ=50. Bottom panel: as top panel but for $P=1.0338$
(middle), borderline case of the critical region, calculated for different
lattice size $L=295.11\mathrm{\mu m}$ (red points), $L=444.42\mathrm{\mu m}$
(blue points) and $L=590.22\mathrm{\mu m}$ (green points). Red points
correspond to the spatial correlation for the “standard” box size, employed
throughout the analysis of this chapter. In both panels power-law (solid
line), stretched-exponential (dashed line) and exponential (dotted line) fits
are reported. The exponent $\alpha$ as a function of the system size $L$ is
reported in the insets of each panel.

Spatial first order correlation function. In Fig. 6(a,c) we show the spatial
first order correlation function $g^{(1)}(\Delta r)$ as a function of the pump
strength $P$, calculated by means of Eq. (8) of the main text. To extract such
a quantity, we first compute the two-point correlation function for each
vertical and horizontal array of a single realisation of the wave function
$\psi(\mathrm{\textbf{r}})$; we then average and normalize before additionally
averaging over the $\mathcal{N}$ realizations, and computing corresponding
error bars. We have moreover tested our results for i) convergence in lattice
size, ii) convergence in stochastic numerical realisations iii) convergence in
computational method. More details are reported in the Supplemental Materials
of Ref. [3].

Temporal first order correlation function. Fig. 6(b,d) shows the temporal
first order correlation function $g^{(1)}(\Delta t)$, calculated by means of
Eq. (9) of the main text. First, we make sure we reach the steady state at a
given pump power $P$ and we then continue running the simulation for a long
enough time (in our case we push the system evolution up to $\Delta
t^{max}\sim 1\mathrm{ns}$). In Eq. (9) of the main text, $\Delta t$ correspond
to the relative time between the two central points of the two fields
considered. In the vicinity of the critical point, we expect strong density
and phase fluctuations. Therefore, a solid averaging procedure is crucial to
achieve smooth and clean results for the correlator curves; this is important
for extracting exponents as precise as possible. In order to reach this, we
average over realisations $\mathcal{N}$, but also over different times
$t_{i}>t_{0}$ with same temporal windows, checking that different $t_{i}$ are
sufficiently far apart so that any memory/correlations of the system are lost.
Therefore, a number of numerical results have been carried out. Specifically,
we test: i) convergence over averaging upon different initial temporal points
$t_{0}$ chosen to be sufficiently far apart in time (see aforementioned
discussion), ii) convergence in the choice of the initial $t_{0}$ (also in
order to confirm that steady-state is reached for the minimum $t_{0}$ chosen),
iii) convergence in the number of stochastic paths $\mathcal{N}$.

We now discuss the procedure for evaluating the exponents and their errorbars.
Fig. 4(a) of the main text illustrates how the critical region is inferred by
the decay of the two-point correlators; here the Root Mean Deviation (RMSD) of
the residuals of the fits are calculated as

$\mathrm{RMSD}=\sqrt{\frac{1}{n}\sum_{n}\left(d_{i}-f_{i}\right)^{2}},$ (8)

where $d_{i}$ and $f_{i}$ are the data points and the fit points at a position
$j$ in space within the fitting window, with $j=1,\dots,n$. For both spatial
and temporal correlators, the exponents plotted in Fig. 4(b,c) of the main
text are extracted by choosing a sensible choice of the spatial (temporal)
window of the fits (gray shadows in Fig. 3). These are chosen so to avoid
discretisation effects at small scales (times) and finite-size effects at
large scales (late-time fluctuations given by the stochastic nature of the
problem). We check the robustness of such a fitting-window choice, and extract
relative errorbars for the power-law exponents (e.g. in Fig. 4(b,c) of the
main text), by quantifying how the spatial (temporal) exponents change by
varying the fitting window limits by $20\%$ of their default values.

##  8 Numerical counting of topological defects

The number of topological defects is evaluated from the phase gradients around
closed paths of each grid point. First, a Gaussian (low-pass) filter is
applied to the wave function, in order to remove all the high frequency noise
components. This removes all noise with wavelength (in pixels) smaller than,
or of the order of, the standard deviation of the filter’s Gaussian kernel. We
proceed by extracting the phase, $\vartheta$, and gradient of the phase,
$\nabla\vartheta=v$, using finite difference methods. Finally, a vortex is
identified (together with its location and charge) in a specific grid point
when the circulation $\Gamma=\int_{C}vd{r}\gtrsim 2\pi$, around a close path
$C$ of double the size of the vortex healing length.

##  9 Binding/unbinding crossover of topological defects

In Fig. 2 of the main text, the average number of topological defects
$\left<N_{\mathrm{v}}\right>$ at steady-state are plotted as a thick red line.
The condensate phase ${\vartheta}=\arg{\\{\psi\\}}$ of three exemplary
configurations of vortices, below, in the vicinity of and above the BKT
threshold are illustrated in Fig. 7. The left panel corresponds to a low-pump
disordered phase (left), characterized by a high density of free vortices
which are free to proliferate. In spite of the high density of topological
defects $d_{\mathrm{v}}=N_{\mathrm{v}}/(L_{x}L_{y})$ at this stage, we have
ensured that the vortex density remains always smaller that the polariton
density $\left<|\psi(t)|^{2}\right>$, in order for our physical description of
vortices to be meaningful. Already for $P\sim 0.9$, we find that
$d_{\mathrm{v}}\sim 0.01\mathrm{\mu m^{-2}}$ and
$d_{\mathrm{v}}/\left<|\psi(t)|^{2}\right>\sim 0.01$. At the threshold point,
$P\sim P_{\mathrm{BKT}}$, (centre) the process of vortices pairing starts to
take place and is depicted. At high pump values $P\gg P_{\mathrm{BKT}}$, in
the ordered phase (right), the system is left free of defects.

##  10 Numerical extraction of the critical momentum

In Fig. 5 of the main text we discuss three exemplary cases of excitation
spectrum calculated from the spatio-temporal Fourier Transform (FT) of
$|\psi{(\textbf{r},t)}|^{2}$ across the PT. For the latter case $P=1.06$
[panel (iv)] the critical momentum $k_{c}(P=1.06)=1.26\times
10^{-2}\mathrm{\mu m^{-1}}$ (calculated by means of Eq. (5) of the main text)
is found on the order of the momentum discretization $\Delta
k=\pi/L=1.06\times 10^{-2}\mathrm{\mu m^{-1}}$ of the numerical simulation. In
order to extract numerically the value of $k_{c}$, a system $100\times$ larger
in size (but with same number of grid points) is simulated and plotted in the
inset of Fig. 5(iv). Finally, with such a choice of parameter, we can access a
system with ($100\times$) smaller $\Delta k$. We note that while with such a
spatial discretization the model is no more able to capture the microscopic
physics of the system, the momentum discretization $\Delta k$ is sufficiently
small to measure a critical momentum $k_{c}^{\textrm{num}}=1.0(5)\times
10^{-2}\mathrm{\mu m^{-1}}$, in good agreement with the analytically predicted
curve [Eq. (6) of the main text)]. In future works, it would be of a great
interest simulating smaller numerical spatial spacing in such a large systems,
so to be able to quantitatively study the effects of the microscale physics
(accessible by means of a smaller numerical spatial spacing) on the value of
the numerical critical momentum $k_{c}^{\mathrm{num}}$.

Figure 9:  Vortex-antivortex pairs for the borderline case. Filtered phase
with detected vortices (red dots) and anti-vortices (green dots) are displayed
for the case $P=1.0338$, for $L=295.11\mathrm{\mu m}$ (left) and
$L=444.42\mathrm{\mu m}$ (right), for the case IPΩ=50. We note that only
vortex-antivortex pairs are present in both configurations, suggesting the
system is in a quasi-ordered phase. Colormap as in Fig. 7.

##  11 Long-range behaviour at the critical point

Figure 10: Exponential and algebraic decay crossover of spatial correlations
for IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$. As Fig. 3(a) of the main
text, but for the case IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$, in log-
linear (top) and log-log (bottom) scales. The values of the pump powers $P$
are reported in the legend.

In this section we describe the behaviour of the correlation function in the
vicinity of a critical point in the (quasi)ordered phase; we discuss presence
of possible finite-size effects and the critical values of the algebraic
exponent around the BKT threshold. As in the main text, we focus on the case
IPΩ=50 for which we identify the upper bound of the critical region to be
located at $P=1.0338$, when the first (spatial and temporal) correlator shows
“clear” power-law decay. Differently from the
IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$ case, for the system size and
parameter set considered in the main text analysis, the maximum exponent
extracted does not exceed the equilibrium upper-bound limit, i.e.
$\alpha_{s}<0.25$.

We then proceed by investigating finite-size scaling properties around this
limit. Fig. 8 shows comparison of spatial correlation functions
$g^{(1)}(\Delta r)$ calculated for different box sizes, for $P=1.0338$ (top),
the borderline case of the ordered regime, and for a point just outside the
critical region $P=1.034$ (bottom). Red points correspond to the correlations
computed with a lattice size $L=295.11\mu m$ and $N=301$ points (as in main
text). Blue and green points (when shown) correspond to $L=444.42\mu m$ and
$L=590.22\mu m$, i.e. to a box size enlarged by $50\%$ and $100\%$,
respectively. The different curves are then fitted with exponential and power-
law functions, as in the analysis performed in the main text. We note that
while for the case $P=1.034$, the curves are converged within the fitting area
(grey faint region), confirming the power-law decay for larger boxes, the case
$P=1.0338$ exhibits a further decrease of the slope as the system size
increases, suggesting that the correlators are indeed not converged, with
finite size effects being still relevant and responsible for the changing of
the curve shapes. It is worth noting that such a behaviour is consistent with
observations within the OPO framework [7].

However, we note that for the limit case $P=1.0338$ the spatial correlators
for experimentally-relevant system sizes still exhibit an algebraic decay
(solid lines) in favour of an exponential behaviour (dashed curves). Thus, we
proceed evaluating the exponent $\alpha$ extracted from the power-law fit for
the case $P=1.0338$. In the insets of Fig. 8, we plot the exponent $\alpha$ as
a function of the box size for each case. We find that the default box size
shows $\alpha=0.23$ (just below the equilibrium threshold) while, for larger
boxes $\alpha$ exceeds the limit of $1/4$, with $\alpha=0.28$ and
$\alpha=0.33$ for $L=444.42\mu m$ and $L=590.22\mu m$, respectively. To
consolidate our conclusions on the BKT-like nature of the system in this limit
case, we investigate its topological defects configuration. A quasi-ordered
regime is indeed confirmed by the presence of only vortex pairs in the system,
as illustrated in Fig. 9 for one of the many realizations, for both smaller
and larger boxes.

Concluding, in this section we demonstrated that in the vicinity of the
critical point, within experimentally-relevant sizes, the system manifests an
exponent $\alpha>1/4$, therefore suggesting its nature to be of the “non-
equilibrium” type, in agreement with previous works [5, 7].

Figure 11: As Fig. 4(a) of the main text, but for the case
IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$

##  12 Spatial coherence for the frequency-independent pump case
IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$

In this section we present the fitting procedure and analysis for the
frequency-independent pump case IP${}_{\Omega=\infty}^{n_{\mathrm{s}}=1500}$.
Similarly to Fig. 3(a) of the main text for the frequency-selective case, Fig.
10 shows the exponential-to-algebraic long-range decay crossover of the
spatial and temporal correlation functions for increasing pump power $P$.
Thick dashed red (blue) curves correspond to exponential (power-law) fitting,
from which values of correlation length $\xi$ and power-law exponents
$\alpha_{\mathrm{s}}$ are extracted and plotted in the inset of Fig. 4(c). For
each curve, we superimpose only the appropriate fit, over the chosen fitting
window (plotted as gray shadow); we show both fits for curves which lie in the
critical region. From this analysis, we find that the critical region lower
and upper boundaries correspond to the values $P=1.0865$ and $P=1.0866$,
respectively. In Fig. 11 we plot the RMDS of the fits.

## References

  * [1] Sinatra A. Castin Y. Phys. Rev. A782008053615.
  * [2] Carusotto I. Ciuti C. Reviews of Modern Physics852013299.
  * [3] Comaron P., Dagvadorj G., Zamora A., Carusotto I., Proukakis N. P. Szymańska M. H. Phys. Rev. Lett.1212018095302.
  * [4] Graham R., Hope J. Johnsson M. Comp. Phys. Commun.1842013201.
  * [5] Nitsche W. et al. Phys. Rev. B902014205430.
  * [6] Zamora A., Dagvadorj G., Comaron P., Carusotto I., Proukakis N. P. Szymańska M. H. Phys. Rev. Lett.1252020095301.
  * [7] Dagvadorj G., Fellows J., Matyjaśkiewicz S., Marchetti F., Carusotto I. Szymańska M. Physical Review X52015041028.