Non-equilibrium Berezinskii-Kosterlitz-Thouless transition in driven-dissipative condensates^(a)

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 three-dimensional 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 two-dimensional (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 [5–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 [5]. 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 [1,7,9–11] and theoretically [3,12–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 has 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 [5,17]. 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 [19], 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 [19,20] 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 [9,20,22–25], and possibly related to single-mode physics [26], power-law decay of temporal correlations has 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.

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)$ [5,27] and includes the complex relaxation processes by means of a frequency-selective pumping source [28–30]. The model, which can be derived from both truncated Wigner (TW) and Keldysh field theory [5,17] reads $(\hbar=1)$

$$\begin{eqnarray} i\textrm{d} \psi &= \bigg[ - \frac{ \nabla^2}{2 m } + g|{\psi}|^2_{-} + \frac{i }{2} \bigg( \frac{P}{1+\frac{|{\psi}|^2_{-}}{n_\mathrm{s}}} - \gamma \bigg) \nonumber\\ & +\frac{1}{2}\frac{P}{\Omega}\frac{\partial}{\partial t} \bigg] \psi \textrm{d}t+ \textrm{d}W, \end{eqnarray} \tag{ 1 }$$

where m is the polariton mass, g is the polariton-polariton interaction strength, γ 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 Ω sets the characteristic scale of the frequency-dependence of gain. The renormalized density $|{\psi}|^2_{-} \equiv (|{\psi} |^2 - {1}/{{(2\textrm{d}V)}} )$ includes the subtraction of the Wigner commutator contribution (where $\textrm{d}V=a^2$ is the volume element of our 2D grid of spacing a). The zero-mean white Wiener noise $\textrm{d}W$ fulfils $\langle\textrm{d}W(\textbf{r},t)\textrm{d}W(\textbf{r}^\prime,t)\rangle = 0$, $\langle\textrm{d} W^*(\textbf{r},t)\textrm{d}W(\textbf{r}^\prime,t)\rangle = [(P+\gamma)/2] \delta_{\textbf{r},\textbf{r}^\prime}\textrm{d}t,$ 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} ({P}/{ \gamma} -1 )$ 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} [ {P}/ ( P g|\psi^{\mathrm{SS}}|^2/ \Omega + \gamma )-1]$ [31].

The effect of small excitations around the bare condensate steady-state solution can be described by means of the linearized Bogoliubov approximation [1,5]. 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 $\textbf{k}$-modes are decoupled, so we can move to Fourier space and define a $\textbf{k}$-dependent Bogoliubov matrix $\mathcal{L}_\textbf{k}$ [29,31],

$$\begin{equation} \mathcal{L}_\textbf{k} = \left( \begin{array}{c@{\quad}c} \Lambda ( {\epsilon_\textbf{k}} + \mu -i \Gamma ) & \Lambda (\mu -i \Gamma ) \\[4pt] \Lambda^{*} (-\mu -i \Gamma ) & \Lambda^{*} ( -{\epsilon_\textbf{k}}- \mu -i \Gamma) \end{array}\right), \end{equation} \tag{ 2 }$$

with $\Gamma = {\gamma} ( P-\gamma )/{2P}$, the free-particle dispersion $\epsilon_\textbf{k} = k^2/2m$, the interaction energy $\mu = g |\psi^{\mathrm{SS}}|^2$ and $\Lambda = ( \gamma_a + i \gamma_b )$ with $\gamma_a={1}/{[1+({P}/{(2 \Omega)} )^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{eqnarray} & \omega_\textbf{k}^{\pm} = -i [ {\gamma_a \Gamma }-{\gamma_b} ({\epsilon_\textbf{k}} + \mu)] \nonumber\\ & \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) }.\nonumber\\ \end{eqnarray} \tag{ 3 }$$

At high momenta k this spectrum recovers a single-particle behaviour with parabolic dispersion, while the frequency dependence of pumping results in an increasing linewidth for growing 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 [5], 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

$$\begin{equation} k_\mathrm{c} = \sqrt{ {2m} \left[ - \frac{ \Gamma {\gamma_b} }{\gamma_a}- \mu +\frac{\sqrt{ (\Gamma^2+\mu^2) ({\gamma_a}^2+{\gamma_b}^2)}}{\gamma_a} \right] }\qquad \end{equation} \tag{ 4 }$$

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

Zoom In Zoom Out Reset image size

We simulate the system dynamics by numerically integrating in time the stochastic differential equations for the polariton field shown in (1); numerical details and methods are reported in the Supplementary Material Supplementarymaterial.pdf (SM). 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 [9,19,30].

Zoom In Zoom Out Reset image size

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\textrm{d}{r}/L_xL_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 (see footnote  ¹ ).

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 [19] —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 Ω, this figure shows the phase diagram for three different choices of parameters, listed in the caption of the figure. For each analysed case, 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 \mu \text{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-antivortex pairs into a plasma of free vortices [18,19]. In fig. 2 the NESS average number of topological defects $\langle N_\mathrm{v} \rangle$ 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 the SM 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 $\langle N_\mathrm{v}\rangle \propto \xi^{-1/2}$ starts to decay as expected for a continuum PT with diverging correlation length ξ, a detailed study of which is presented in ref. [32]. 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, $\langle N_\mathrm{v}(P) \rangle$ 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.

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$_{\Omega = 50}$ case; the complementary study for the IP$_{\Omega = \infty}^{n_\mathrm{s}=1500}$ case with a frequency-independent pump is illustrated in the SM. Within our semi-classical model, the spatial and temporal two-point first-order correlation functions are defined, respectively, as

$$\begin{eqnarray} &&g^{(1)} (\Delta r) =\frac{\langle \psi^* ({r_0} + \Delta r,t) \psi^{} ({r_0},t) \rangle }{\sqrt{ \langle | \psi({r_0} + \Delta {r},t)|^2 \rangle \langle | \psi({r_0},t) |^2 \rangle}}, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &&g^{(1)} (\Delta t) = \frac{\langle \psi^* ({r_c},t_0) \psi^{} ({r_c},t_0 + \Delta t) \rangle }{\sqrt{\langle | \psi({r_c},t_0)|^2 \rangle \langle | \psi({r_c},t_0+ \Delta t) |^2 \rangle}}, \end{eqnarray} \tag{ 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.

Zoom In Zoom Out Reset image size

Inspired by earlier works [19,33], we characterize the behaviour of the steady-state correlation functions as a 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 function ² . 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 (see the SM). 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 (see the SM). In fig. 3 we superimpose on top of each correlation function g⁽¹⁾ the most accurate fitting curve, represented by red or blue dashed lines in the exponential or power-law cases, respectively.

Figure 4(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

$$\begin{equation} \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}}. \end{equation} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

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. 4) for both spatial and temporal correlation functions ³ .

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 ξ, 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. 4(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. 4 are coloured in a way to match the ones of fig. 3. In the disordered phase (left part of fig. 4(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 down ⁴ . In the quasi-ordered phase (right part of fig. 4(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} \mu \text{m}^{-1}$ is on the order of the momentum discretization $\Delta k = \pi/L=1.06 \times 10^{-2} \mu \text{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} \mu \text{m}^{-1}$ is extracted by simulating a system with a $(100 \times)$ smaller $\Delta k$ (see the SM): 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.

Zoom In Zoom Out Reset image size

Previous experimental [20] and theoretical [19] 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$_{\Omega=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 the SM, by enlarging the box by 1.5 and 2 times, power-law exponents are found to lie within the interval $0.25 < \alpha_s < 0.35$ (reported in fig. 4(b), (c) as a large errorbar for the brown point $\alpha_\mathrm{s}(P=1.0338)$). This confirms the argument that for a 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. 4(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}$ are due to its non-equilibrium nature or finite-size effects, or an interplay between the two.

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 [19] 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 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.

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 (PC and NPP)), and InterPol (through the National Science Center, Poland, Grant No. 366 2017/25/Z/ST3/03032 (PC) and the EPSRC Grant No. 368 EP/R04399X/1 (MHS)). We acknowledge EPSRC Grant No. EP/K003623/2 (MHS) and support from the Provincia Autonoma di Trento and from the Quantum Flagship Grant PhoQuS (820392) of the European Union (IC).