Coherence of a non-equilibrium polariton condensate across the interaction-mediated phase transition

The physics of collective phenomena in quantum many-body systems is strongly connected to dimensionality, interactions and the intrinsic equilibrium or nonequilibrium nature of the system. A gas of identical bosons is a perfect example of the crucial role that these factors play in the onset of order¹. A three-dimensional (3D) gas of bosons can undergo Bose-Einstein condensation (BEC) transition, which results in a state of matter exhibiting long-range order. While thermal fluctuation preclude the formation of the long-range order in a two-dimensional (2D) gas, it can instead undergo the Berezinskii-Kosterlitz-Thouless (BKT) transition and exhibit superfluidity below a finite critical temperature. In the superfluid phase, the decay of the spatial correlations of the system is predicted to be a function of temperature and the superfluid density^2,3. In order to achieve a BEC phase transition in 2D, the bosons have to be non-interacting and trapped⁴. Both BEC and BKT transitions have been explored in the system of cold atomic gases by tuning the system’s dimensionality or the interaction strength^{5,6,7,8,9,10,11,12}.

2D bosonic systems, such as cavity photons¹³ or exciton polaritons^14,15, represent a class of systems with driving and dissipation, where the full thermalisation (energy relaxation) may not be achieved¹⁶ and the physics is inherently nonequilibrium. While cavity photons closely resemble an ideal (noninteracting) gas, exciton polaritons (polaritons herein) arise from the strong coupling of cavity photons and excitons in a semiconductor and interact strongly due to their excitonic component. Theoretical^17,18,19,20, and experimental^21,22 investigations demonstrated the relevance of the interaction-dominated BKT physics to polariton systems under nonresonant pumping, when the coherence in the system develops spontaneously. However, so far it remains unclear how the coherence in the intrinsically driven-dissipative polariton system depends on the interactions, how it is affected by confinement, and whether there exists, in analogy with conservative systems, any relation between the rate of decay of coherence and the superfluid fraction. These questions are also interesting because both the thermalisation (energy relaxation) and the phase transition in a polariton system are driven by interactions¹⁶.

In this work, we investigate the phase coherence, i.e. the first-order spatial correlation function, of a confined polariton condensate in an optical microcavity. The confinement allows us to explore the coherence features at high densities, when the interaction energy dominates over the kinetic energy. In the experiment, we vary the polariton interaction energy by varying the polariton density and the strength of interparticle interactions. The variation of the interaction strength is achieved by changing the energy detuning between the cavity photon and the exciton, which changes the excitonic/photonic fraction of the polariton. The larger excitonic fraction corresponds to stronger interparticle interaction. In turn, the condensate density is controlled by the power of the optical pump.

Our study reveals that the trapped, condensed polaritons exhibit different regimes of spatial coherence depending on the interaction strength. Full coherence across the spatial extent of the system is reached for stronger interacting (more excitonic) polaritons. The study of coherence dependence on the density for a fixed interaction strength (fixed excitonic fraction) shows similar behaviour, as the interaction energy grows with increasing density. However, at very high excitonic fractions and pump powers the coherence decreases due to an increasing role of various decoherence mechanisms. Nevertheless, at high densities and/or excitonic fraction the spatial decay of the first-order correlation function displays the power-law dependence consistent with a BKT phase. In the limit of weakly interacting (more photonic) polaritons, full thermalisation does not occur, and multi-mode condensation takes place instead of a ground state BEC predicted for a non-interacting trapped 2D bosonic system⁴. This confirms that the thermalisation in this system is driven by interactions and is inhibited in the photonic regime, making it difficult to reach the weakly or non-interacting BEC state. Since both condensation and thermalisation processes in polariton condensates are strongly interaction-driven, observation of a clean crossover from the interacting BKT to the non-interacting 2D BEC regimes by tuning the interaction strength (as observed in cold atoms using Feshbach resonances²³) is challenging.

To further characterise the transition to the BKT phase in this driven-disspative system, we introduce the concept of “coherent fraction” of the polariton condensate, and study its relation to the power-law exponent that describes the decay of the spatial first-order correlation function. Our analysis demonstrates growth and decay of the coherent fraction with the increasing interaction energy, consistent with the non-monotonic behaviour of the first-order spatial correlation function. Our findings also suggest a clear relationship between the power-law exponent, the coherent fraction and the effective temperature of a polariton system, in analogy with the relationship between the power-law exponent and superfluid fraction for quantum gases in thermal equilibrium²⁴.

We study polariton condensates confined in an optically-induced circular trap^{25,26,27,28,29,30} in the steady-state regime. An off-resonant quasi-continuous wave (cw) excitation laser, shaped into a ring profile, pumps the microcavity sample with embedded GaAs quantum wells (see Methods), creating a similarly shaped excitonic reservoir that feeds and repels polaritons. As schematically shown in Fig. 1a, above a critical pump power (reservoir density), this configuration produces a trapped condensate inside the ring. The constant replenishment of the reservoir by the excitation laser ensures that the condensate density remains relatively constant despite the losses due to the finite polariton lifetime (~100 ps). The polaritons decay by emitting photons, which we detect in the experiment, and the intensity of the photoluminescence is directly proportional to the polariton density (Fig. 1b). The real space image (Fig. 1c), momentum space distribution, spectra, and first-order correlation functions of the photoluminescence are used to infer the information about the polariton condensate.

a Sketch of the excitation scheme where a ring-shaped laser excites the sample and creates a ring-shaped reservoir of excitons with the density n_r that traps the top-hat condensate with the density n_c. b Comparison of the cross-sections of the measured intensity of the condensate photoluminescence I (proportional to the polariton density) and the first-order correlation function ∣g⁽¹⁾∣ showing that the spatial extent of ∣g⁽¹⁾∣ is limited by the size of the condensate. Scale bar: 10 μm. c Example experimental image of the real-space distribution of the condensate. Scale bar: 10 μm. d Measured 2D g⁽¹⁾(r) of the condensate in (c), acquired using a modified Michelson interferometer (see ‘Methods’). Scale bar: 20 μm. e Example density dependence of the blueshift and healing length. The pseudo-color plot shows the spectra, with each spectrum normalised to 1 at its corresponding density. The excitonic fraction here is ∣X∣² = 0.34.

Full size image

At high densities of the condensate, when its interaction energy dominates over the kinetic energy, the density profile shows a relatively flat-top distribution (Fig. 1b). By virtue of the Thomas-Fermi approximation, this reflects a flat, box-like profile of the underlying trapping potential^25,27. The uniform profile of both the condensate and the reservoir-induced potential ensures that the densities of the condensate and the reservoir are spatially homogeneous near the centre of the trap. Note that this regime only works in our experiments for sufficiently large traps with a diameter around 30 μm. Smaller traps lead to quasi-harmonic potentials while very large traps need high excitation powers, beyond our experimental capability, to reach the high-density regime. In this work, we probe interaction energies up to 700 μeV (see Supplementary Note 2 for details).

A wedge in the microcavity sample arises during the fabrication process and leads to a spatial variation of the exciton-photon detuning^25,29, which translates to the variation of the excitonic fraction quantified by the Hopfield coefficient 0 ≤ ∣X∣²≤ 1 (see ’Methods’). The range of the excitonic Hopfield coefficient accessible in our sample is 0.21 ≤ ∣X∣²≤ 0.6. The excitonic fraction in the polariton quasiparticle is varied within this range by moving the laser excitation spot across the sample.

We use a modified Michelson interferometer³¹ to measure the two-dimensional first-order correlation function ∣g⁽¹⁾(Δx, Δy)∣ (see ‘Methods’), where (Δx, Δy) is the displacement from the autocorrelation point. An example of ∣g⁽¹⁾∣ is shown in Fig. 1b for an excitonic fraction of ∣X∣² = 0.37 with the peak density of n = 4.13 × 10³ μm⁻². It features a peak at the autocorrelation point Δr  = 0, which characterises short-range correlations. The slowly decaying shoulder at larger distances characterises the long-range order in this system which is the main interest of this work. By comparing the ∣g⁽¹⁾∣ to the condensate profile, one can clearly see in Fig. 1c, d that the long-range order is limited only by the size of the condensate. Note that despite the small inhomogeneities of the measured condensate profile, which can arise from both imperfections of the imaging setup and sample disorder, the measured ∣g⁽¹⁾∣ remains smooth, which strongly suggests that our measurement is not affected by the small perturbations of the steady-state density.

Density dependence of coherence at a fixed excitonic fraction

To understand how the condensate interaction energy affects the first-order correlation function, we performed a series of measurements at different excitation powers above the condensation threshold, resulting in different condensate densities n. Here, n is the average density in a small area around the centre of the condensate.

It is important to highlight that the size of the condensate (diameter (D approx 30 , mu m)) is an order of magnitude larger than the healing length (xi =hslash /sqrt{2gnm}), where g is the polariton-polariton interaction strength²⁵, m is the polariton mass, and ℏ is the Planck’s constant. The healing length ranges from 6 μm to 1 μm at the probed densities, while the measured condensate blueshift can reach 800 μeV as shown in Fig. 1e. The condensate interaction energy gn is less than this due to nonvanishing reservoir density inside the trap^27,32,33 (see Supplementary Note 2 for details). The small healing lengths relative to the system size guarantee that we are indeed probing the effect of interactions on the long-range coherence of the condensate. If ξ ~ D, the long-range universal behaviour in g⁽¹⁾ is masked by other short-range effects.

Figure 2 a presents the measured ∣g⁽¹⁾∣ at a fixed value of the excitonic fraction ∣X∣² = 0.37, which clearly shows that the long-range coherence increases with density. At low density, ∣g⁽¹⁾(Δx)∣ quickly decays with distance Δx. This spatial decay diminishes with increasing density until ∣g⁽¹⁾∣ becomes almost flat.

a Logarithmic plots of g⁽¹⁾ at ∣X∣² = 0.37 for different densities. Dots are experimental data and solid (dashed) lines are the exponential (power-law) fits (see Eq. (1)). b g⁽¹⁾(Δx) calculated using the stochastic model (2) with the experimental parameters. Extracted (c) power-law exponents α and (d) coherence lengths ℓ, respectively, as a function of density for different excitonic fractions ∣X∣². The black solid line and squares correspond to the pump strength in simulations. The horizontal dashed lines in (c, d) show the equilibrium limit for the algebraic exponent α =  0.25³ and the trap size, respectively. The error bars in (c, d) are estimated from the covariance matrix of the fit.

Full size image

To quantify the effect of interactions on the coherence length, we fit the experimental data using the functions⁴

According to the BKT theory of a spatially isotropic condensate, the fast, exponential (slow, power-law) decay of the first-order spatial coherence function should apply in the disordered (quasi-ordered) phase, below (above) the critical point P/P_th, where P_th is the pump power corresponding to the condensation threshold in our system³. For each data sets of Fig. 2a, we show the fitting to Eqs. (1) as solid blue and dashed red lines for the exponential and power-law fitting, respectively. The data are fitted considering the power-law exponents α and the coherence lengths ℓ as free parameters of the fit. The extracted quantities are then depicted in Fig. 2c, d. At larger densities, ℓ becomes much larger than the trap size (dashed line) while α decreases asymptotically to zero, signalling full coherence across the system.

We recall that in the canonical BKT transition³, the power-law decay of the coherence function is expected to exhibit a crossover in α from α_c = 0.25 at the critical point to 0 for large densities. In the driven-dissipative case, it has been previously shown that the actual α_c may assume larger values due to the non-equilibrium nature of the polariton system introduced by the pump fluctuations^34,35. Aiming to identify the transition point, we proceed by comparing the root mean square (RMS) of the fitting curves, a method employed in a previous work to ascertain the phase of the system¹⁷. The analysis of the RMS of the fits (reported in the Supplementary Note 3) reveals that even at very large densities, where power-law decay is expected, the exponential and power-law fits are still comparable. As a consequence, the ratio between the RMS of the fits cannot be used, in principle, to identify the critical point¹⁷. However, we note that in the regions where BKT behaviour is expected, the extracted coherence length is significantly larger—by about one order of magnitude or more—than the system size (see Fig. 2d). At the same time, the power-law exponent is nearly zero (see Fig. 2c). This observation suggests that the exponential behaviour is likely attributable to finite-size effects rather than the underlying physics of phase fluctuations.

We support our experimental finding with the numerical modelling of the polariton system using the stochastic driven-dissipative Gross-Pitaevskii equation³⁶. As described in “Methods”, we first evolve the model to reach a steady-state. Then, the first-order spatial correlation function (see Methods) is numerically computed and plotted in Fig. 2b for different densities. The resulting curves are then fitted to exponential and power-law functions. Parameters from the fits are extracted and compared to the experimental values in Fig. 2c, d. Overall, a good qualitative agreement between the theory and the experiment is found, especially for larger pump powers (condensate densities).

It is important to note that the slowest decay of correlations (i.e., the smallest α) is attained, both in experiment and in simulations, only up to n ~ 2 × 10³μm⁻². For larger densities, we observe a loss of coherence, followed by an increase of α and a reduction of ℓ. The plot of the numerically computed pump power required to achieve a particular value of the condensate density, shown by a black solid line in Fig. 2d, confirms that the system needs a dramatically stronger pump power to reach larger densities. Strong external drive leads to a strong decoherence of the condensate. Pump sources, in fact, are expected to introduce noise leading to dephasing at all powers^37,38,39. At stronger powers we expect these to be large. In addition, the dephasing may also arise from thermal fluctuations⁴⁰, the incoherent excitonic reservoir^41,42, and polariton-polariton interactions⁴³. In our numerical simulations, the pump noise is explicitly included in the noise term dW in Eq. (2) (see Methods), which accounts for the quantum fluctuations¹⁵.

Dependence of coherence on the excitonic fraction at a fixed polariton density

To further understand the dependence of the spatial coherence on interactions, we change the excitonic fraction ∣X∣² (i.e. the interaction strength, see Methods) and measure coherence for different polariton densities n. We then investigate the behaviour of spatial coherence as a function of ∣X∣² for a fixed polariton density n ~ 10³μm⁻². We analyse the data for the coherence, similarly to the analysis presented in Fig. 2, and display the results in Fig. 3.

a Logarithmic plots of |g⁽¹⁾| at constant density n ≈ 1 × 10³μm⁻² for different excitonic fractions ∣X∣². Dots are experimental data and solid (dashed) lines are the exponential (power-law) fits. Theoretical results are reported in (b). Extracted (c) power-law exponents α and (d) coherence lengths ℓ as a function of the excitonic fraction. The pump strength in simulations is shown by a black solid line and squares. Dashed lines in (c, d) correspond to the critical exponent α = 0.25 at the thermodynamic limit for the conservative case²⁴ (red), and to the system size (blue), respectively. The error bars shown (c, d) are estimated from the covariance matrix of the fit.

Full size image

Figure 3a shows that the system undergoes a crossover from low to high coherence as the excitonic fraction of the polariton increases, quantified by the fitting curves and parameters shown in Fig. 3a and c, d, respectively. This result clearly illustrates the interaction-driven nature of the coherence crossover for our system. At large excitonic fractions, the system shows a lack of coherence that we attribute to the dephasing mechanisms observed at large pump powers and discussed in the previous section. In fact, to reach similar densities, larger excitonic fractions require stronger pumps. This trend is clearly revealed by our simulations as shown in Fig. 3d, which depicts the growth of the required pump power as a function of the excitonic fraction when the system density is kept constant. Stronger dephasing may be induced at larger excitonic fractions by stronger polariton-polariton and polariton-reservoir interactions, which increase with the excitonic Hopfield coefficient. Furthermore, the trapping potential becomes tighter at larger excitonic fraction, resulting in a stronger overlap of the condensate with the reservoir, as previously shown in similar experiments²⁵. The largest ∣X∣² probed here already shows this effect.

The experimental results shown in Fig. 3a, c, d, are corroborated by the numerical simulations: the numerically extracted spatial coherence function, calculated at a fixed density (n ~ 10³μm⁻²) for different values of ∣X∣², are shown in Fig. 3(b). The extracted coherence lengths and power-law exponents are then compared to the experimental data in Fig. 3c, d. We find a very good agreement between the theory and experiment in the limits of both small and large excitonic fractions.

Both increasing the excitonic fraction of the polaritons at a fixed density and increasing the density at a fixed excitonic fraction, leads to the growth of the interaction energy. To better visualise the coherence behaviour in the space of these two parameters, in Fig. 4 we report the experimentally extracted α and ℓ as a function of density, Fig. 4a, and as a function of excitonic fraction, Fig. 4b. Panels (c) and (d) illustrate these parameters as a surface in the (n, ∣X∣²) space, showing the onset of coherence at lengths above the system size (shown as a dashed red line). In summary, in the confined system, coherence is enhanced when increasing the interaction energy by means of increasing either density or the excitonic fraction. The crossover to the full coherence takes place approximately at the system size, as predicted by the theory of trapped conservative quantum gases⁴. In the limits of both large densities and excitonic fractions, the polariton condensate experiences dephasing (due to the pump noise, polariton-polariton and polariton-exciton intercations), resulting in a reduction of coherence. As a result, the smallest α or largest ℓ occurs around ∣X∣² = 0.45 and n ~ 2 × 10³μm⁻².

Extracted power-law exponents α (top) and coherence lengths ℓ (bottom), as a function of (a) density and (b) excitonic fraction, for different experimental series. Simulation results are presented as green squares. c Power-law exponent α and (d) coherence length ℓ as carpet plots in the density-excitonic-fraction plane. Dashed lines in (c) and (d) correspond to the critical exponent α = 0.25 at the thermodynamic limit for the conservative case²⁴ (red), and to the system size (blue), respectively. The error bars in (a, b) are estimated from the covariance matrix of the fit.

Full size image

Coherent fraction and thermalisation

We further investigate a possible relationship between the coherence and the thermalisation properties of our system. We start by obtaining energy-resolved momentum-space images using a monochromator coupled with a CCD camera. The CCD intensity, I_CCD, of the central section of each image is collected as a function of momenta, and recorded for different densities (pump powers) and excitonic fractions (exciton-photon detunings). In Fig. 5a we report the normalised I_CCD for three cases: ∣X∣² = 0.21, 0.34, 0.60. For smaller ∣X∣² and at lower density, the system is in a multimode regime. At larger ∣X∣², the system displays a single-mode occupation at all densities. This strongly suggests that a larger excitonic fraction leads to a more efficient energy relaxation towards the ground state²⁵. In the context of conservative condensates, quasi-long-range coherence leads to the formation of different coherent components⁴⁴: i) the condensate fraction, defined as the normalised density of particles which populate the zero-momentum (k = 0) mode, ii) the superfluid density defined à la Nelson-and-Kosterlitz⁴⁵ and iii) the quasicondensate, defined as a condensate with a fluctuating phase⁴⁶. The superfluid density, in particular, can be directly linked to the power-law exponent α charactersising the BKT transition^3,24. For a non-equilibrium, open-dissipative condensate studied here, extraction of n_s is diffcult. Instead, we introduce the notion of ‘coherent density’ n_coh, defined as the fraction of particles below a certain momentum cut-off k_cut(∣X∣²), the latter obtained by Fourier-transforming the effective trapping potential that the system experiences at each value of the excitonic fraction.

a Momentum space distribution at different densities for three representative excitonic fractions ∣X∣² = 0.21, 0.34, 0.60. The averaged-particle-number momentum distribution at different densities: n ~ (b) 10³μm⁻², (c) 2 × 10³μm⁻² and (d) 3 × 10³μm⁻². e The coherent fraction extracted from the momentum distributions of polaritons. f–h At low densities, the energy-distributions of polariton number (dots and solid lines) are fitted with the Bose-Einstein distribution (dashed lines). RMS (root mean squares) of the fits are reported in (i). The low-density region of interest in panels (f–i) is highlighted in Fig. 4c. The error bars are estimated from the covariance matrix of the fit.

Full size image

In order to extract the coherent fraction of the condensate, we first obtain the particle distribution per momentum state, ({{mathcal{N}}}(k)), which is proportional to the CCD intensity, i.e., ({{mathcal{N}}}(k)={{mathcal{A}}}times {I}_{{{rm{CCD}}}}) taking into account the geometry of the system⁴⁷ (see Supplementary Note 2 for details). Figures 5b, c, d show I_CCD for different ∣X∣² at fixed densities n = 1, 2 and 3 × 10³μm⁻². From these distributions, we extract the coherent fraction n_coh/n_tot in the (high-density) single-mode regime, and plot its behaviour as a function of ∣X∣² at different values of the total density n_tot in Figure 5e. It is evident that the coherent fraction is increasing as a function of density and excitonic fraction (i.e., interaction strength). We note that a similar dependence on the interaction strength is expected for the superfluid fraction in a conservative system²⁴. For ∣X∣² larger than  ~0.35, the coherent fraction saturates at n_coh/n_tot ~ 0.9. At much larger excitonic fractions, the coherent fraction decreases, in agreement with the trend shown by the first-order correlation functions discussed in the previous sections.

In the second part of this analysis, we discuss thermalisation properties of the system. First, we move from momentum to energy space, employing the energy-momentum mapping provided by the dispersion relation of the lower-polariton branch ({E}_{{{rm{LP}}}}(k)=1/2[{E}_{{{rm{X}}}}(k)+{E}_{{{rm{C}}}}(k)-sqrt{{hslash }^{2}{Omega }^{2}+{delta }^{2}}]), where E_X(k) and E_C(k) are the energies of the exciton and cavity photon, respectively, δ(k) = E_C − E_X is the energy detuning, and ℏΩ is the Rabi splitting energy. We tentatively fit the data ({{mathcal{N}}}(E)) to the Bose-Einstein (BE) distribution of the average particle number ({{{mathcal{N}}}}_{{{rm{BE}}}}(E)=1/(exp ((E-mu )/{k}_{B}T)-1)), where T and μ are the temperature and chemical potential of the polaritons, respectively, and k_B is the Boltzmann constant. Here we treat the temperature T and chemical potential μ as free parameters. We are careful to apply this fitting only at low densities and high excitonic fraction, over a sufficiently populated tail of the modes which otherwise is obscured by a very bright (and peaked) main mode. This area is highlighted as a faint red shadowed region in Fig. 4c.

In Fig. 5f, g, h, we report the particle distributions in energy, with fitting curves as dashed lines. In panel (i) we show the RMS of the BE distributions; these give a quantitative tool for determining the goodness of the fit to a thermal distribution. Our results suggest that the fitting to a thermal tail improves as the densities and excitonic fractions increase. Note that at larger densities (n  ~ 10³μm⁻²), the fitting procedure deteriorates for larger ∣X∣²; this is in connection with the increase of the k = 0 peak which, at large excitonic fractions, dominates over other modes. Finally, we study the relation between the power-law exponent α and the coherent density n_coh, with the aim to draw a parallel with the relationship between α and the superfuid density n_s in conservative systems. In the canonical BKT theory, the two quantities are related by the expression α = 1/n_sλ² where (lambda =h/sqrt{(2pi m{k}_{B}T)}) is the thermal de-Broglie length⁴.

The power-law exponent and the coherent fraction of the condensate in our system are plotted in Fig. 6a, b, in the small-density and high excitonic component limit accessible from the previous analysis. Figure 6c clearly depicts the linear dependence α (propto) 1/n_cohλ² for the data sets obtained in our experiment. Note that this linear behaviour persists across almost two orders of magnitude for both α and 1/n_cohλ². It also shows a proportionality constant diverging from the expected value of 1 by two orders of magnitude, a deviation that may be due to the driven-dissipative and/or confined nature of the system.

a Algebraic exponents α and b coherent fraction n_coh in the low-density and high excitonic fraction regime (highlighted in Fig. 4c, d as a dashed area), for different densities (see legend in c). c Algebraic exponent shows an inverse proportionality to the coherent fraction, in analogy with the theory for conservative superfluids, where α = 1/n_sλ², where (lambda =h/sqrt{2pi m{k}_{B}T}) is the thermal de-Broglie wavelength, and n_s is the superfluid density. The dashed line corresponds to α = 10²/(n_sλ²). Linear scales plot is shown in the inset. The error bars are estimated from the covariance matrix of the fit.

Full size image

In this work, we investigated the influence of interactions and energy relaxation on the development of spatial coherence in a confined, high-density polariton system. This differs from past research, which identified the signatures of the BKT transition by examining the coherence of a low-density, untrapped polariton flow²². Our confinement technique also reduces the overlap with the incoherent excitonic reservoir. This is in contrast to previous studies^21,34 that used a top-hat excitation profile, with the condensate coherence suffering from the high degree of overlap with the reservoir. Using this confinement technique, we also avoid instabilities, like condensate fragmentation due to spatial hole-burning effect⁴⁸. A more detailed discussion on the BKT transition in nonequilibirum systems is contained in Supplementary Note 1.

Our experimental and numerical findings confirm that, in the driven-dissipative scenario of polariton condensation, the system exhibits behaviour consistent with a BKT-type phase transition, which is driven by interactions. As the interaction energy is tuned across a wide range (up to 700 μeV) by varying the density and interparticle interactions, we observe different regimes of coherence. In the high-density regime of stronger interacting, highly excitonic polaritons full coherence across the spatial size of the system can be reached. However, the coherence is degraded above certain densities and excitonic fractions. Our results suggest that this behaviour is driven by the dephasing due to the fluctuations of the laser pump, which are stronger at higher pump powers, as well as due to polariton-polariton interactions and interactions of polaritons with the reservoir excitons. This is exacerbated in the experiments due to the tightening of the reservoir-induced trap at higher excitonic fractions²⁵, resulting in stronger reservoir-condensate overlap.

In the low-density regime, we clearly see the short-range thermal component of the first-order correlation function⁴⁹. In the limit of highly photonic polaritons, we observe multi-mode condensation due to inefficient energy relaxation towards the ground state. We find that the fitting to BE statistics improves for stronger interactions and at larger densities, proving that thermalisation processes in non-equilibrium condensate take place via interactions. We note that this inhibited energy relaxation towards the ground state precludes the observation of a BEC phase in confined, 2D gas of weakly interacting, highly photonic polaritons.

Finally, we tentatively investigate the relationship between the power-law exponent of g⁽¹⁾(r), the condensate coherent fraction, and the effective temperature of the system. We find that the power-law exponent is directly proportional to the temperature of the system, which defines the thermal de Broglie wavelength of a polariton, and is inversely proportional to the coherent density, in analogy to the dependence of the power-law exponent on the superfluid density in the conservative case. However, we note that the system temperature in conservative systems is extracted from the fitting with a thermal distribution, while in a driven-dissipative system, full equilibration is not guaranteed and the BE fitting is useful only in a limited subset of the system’s parameters. To our knowledge, prior to this study, no investigation had explored this particular relation in the context of nonequilibrium, driven-dissipative quantum systems. Further experimental effort could address the issue of the extraction of the superfluid fraction (density) of the polariton condensate and its relation to the coherent fraction defined here.

Summarising, our work represents an important step towards the characterisation of interaction-mediated BKT phase transition and coherence in the driven-dissipative quantum system, across a range of interaction energy. Our findings are particularly relevant to studies of superfluid properties⁵⁰, universal behaviour⁵¹, and turbulent features⁵² in driven-dissipative quantum systems.

Sample and experimental setup

We use a high Q-factor GaAs/AlGaAs planar microcavity with 7-nm GaAs quantum wells embedded in a 3λ/2-cavity, similar to the samples used in refs. ^31,53. The thickness gradient of the sample enables control of the exciton-photon detuning by probing different positions on the sample. The 720-nm off-resonant continuous-wave (CW) excitation laser is spatially structured into a ring profile using an axicon, a pair of lenses, and a 50X 0.5-NA objective. The same objective collects the polariton emission which goes through a relay of lenses for real-space imaging with a sCMOS camera (Andor Zyla 4.2). The camera images the interference fringes used to measure the first-order correlation function. A monochromator (Princeton Instruments IsoPlane SCT 320) coupled with a EMCCD camera (Andor iXon Ultra 888) enables energy-resolved real and momentum-space imaging.

Tuning the excitonic fraction

The excitonic fraction ∣X∣² is tuned by probing areas of the samples with different photon-exciton detuning ({{rm{k} = 0}, , delta }). The sample has a gradient in detuning due to a wedge in the effective cavity length⁵³ resulting in a gradient in cavity photon energy. The detuning is calculated from the polariton dispersion, measured at very low excitation power (below condensation threshold) for each position on the sample, and using previously known parameters^25,27. We then extract the excitonic fraction from the detuning as (| X{| }^{2}=1/2(1+delta /scriptstylesqrt{{delta }^{2}+{Omega }^{2}})), where the Rabi splitting is about Ω = 15.9 meV²⁷. The exciton-photon detuning varies around 0.25 meV over the 30-μm diameter of the trap, which is equivalent to Δ∣X∣² = 0.004. Hence, ∣X∣² is approximately constant inside the trap.

Coherence measurement setup

We employ a modified interferometer³¹ consisting of a 50:50 beam splitter, a flat mirror, and a hollow retroreflector to measure the spatial first-order correlation function ∣g⁽¹⁾∣. The setup combines the polariton real-space emission I(r) = ∣ψ(r)∣² with its inverted image I( − r) = ∣ψ( − r)∣² onto a sCMOS camera where the interferogram image F(r) forms. The mirrors are aligned such that the autocorrelation point, (Δx, Δy) = 0, is at the centre of the real-space distribution. The first-order correlation function is then extracted using the equation

where V = (F_max − F_min)/(F_max + F_min) is the visibility of the interference fringes extracted by introducing sub-micron delays using a piezo actuator³¹. The mirrors are mounted on motorized stages to ensure near-zero temporal delay Δt ≈ 0. A piezo actuator provides a series of sub-wavelength delays near Δt ≈ 0 to accurately determine the spatial variation of the fringe visibility, V(r). To filter out low-frequency noise coming from the mechanical vibrations of the whole experimental setup, the CW laser is chopped down to  ~10–100-μs pulse lengths depending on the signal-to-noise ratio. This time scale is orders of magnitude longer than the polariton lifetime ( ~100-ps) and the pulse rise time ( ~50 ns), and hence is sufficient to capture the steady-state regime.

Modelling of the polariton field

We describe the effective dynamics of the polariton fluid through the equation of motions for the two-dimensional polariton field ψ = ψ(r, t) as a function of the position r = (x, y) and time t within the truncated Wigner approximation¹⁵, which reads (ℏ = 1)^36,54,55:

where m is the polariton mass, γ_c and γ_R are the polariton and excitonic reservoir decay rates, P(r) the circularly-shaped external drive, g_c is the polariton-polariton interaction strength. R is the scattering rate of reservoir particles into the condensate. The renormalised density (| psi {| }_{-}^{2}equiv left({leftvert psi rightvert }^{2}-1/2dVright)) includes the subtraction of the Wigner commutator contribution (where dV = a² is the element of volume of our 2D grid with the lattice spacing a). The zero-mean white Wigner noise dW fulfils the condition (langle d{W}^{* }({r},t)dW({{r}}^{{prime} },t)rangle =[((R/{gamma }_{R})P+{gamma }_{c})/2]{delta }_{{r},{{r}}^{{prime} }}dt). The equation above includes the phenomenological term proportional to the constant β, which quantifies the rate of energy relaxation in the system^{54,55,56,57,58}. The model Eq. (2) corresponds to the adiabatic approximation limit of the generalised polariton equations of motion coupled to an external reservoir³⁶, which is justified once the reservoir is able to adiabatically follow the evolution of the condensate and γ_c  ≫ γ_R⁵⁹.

Using the parameters adopted in ref. ⁴⁸ to model an experimental system similar to that described here, in this work we vary the excitonic fraction ∣X∣² in order to model the steady-state coherence across the range of interaction strengths. It is therefore crucial to note that most of the parameters in Eq. (2) depend on the excitonic Hopfield coefficient X. The latter also defines the value of the excitonic fraction ∣X∣², and depends on the exciton-photon detuning δ and the Rabi splitting ℏΩ as defined above. The polariton decay rate depends on the excitonic fraction as γ_c = (1 − ∣X∣²)γ_ph, where γ_ph = 1/τ_ph corresponds to the inverse of the photon lifetime. The constants g and g_R characterise the strengths of polariton-polariton and polariton-reservoir interactions, which become stronger for polaritons with a larger excitonic fraction. They can be estimated g_c = g_ex∣X∣⁴, g_R = g_ex∣X∣^2 60, with g_ex the exciton-exciton interaction. Following the discussion in ref. ⁴⁸, we assume that the stimulated scattering rate R from the reservoir into the polariton states is more efficient for more excitonic polaritons: from our results we find that (R={R}_{0}{({g}_{c}/{g}_{R})}^{2}), namely R scales quadratically with ∣X∣².

The circularly symmetric pump profile P(r) is implemented as the sum of two contributions: i) the profile (P({{bf{r}}})={P}_{0}{gamma }_{c}/R[{e}^{-{left({{mathcal{E}}}/2{sigma }^{2}right)}^{2}}]), where P₀ is the pump strength, and ({{mathcal{E}}}={x}^{2}+{y}^{2}-{R_p}^{2}) is the circle equation with radius R_P, and ii) a weak flat pump ({P}={alpha }_{p}{P}_{0}{gamma }_{c}/R), with α_p <  1 positive and small, at the centre of the ellipse for r < R_P: with the latter, we account for the presence of a non-vanishing reservoir density at the centre of the trap. We note that the inclusion of such a term is important for enhancing the onset of coherence at the centre of the trap. We also find that the experimental data are adequately modelled by relaxation mechanisms^56,57,61 (controlled by the paramater β) being more efficient for larger excitonic fractions as β = β₀ + β₁∣X∣². In order to match the experimental measurements, we use the following parameters: ({m_{pol}})  = 3.6 × 10⁻⁵m_e with m_e the electron mass, τ_ph = 135 ps, γ_R = 10⁻³ps⁻¹, g_ex = 1.12 μeVμm², R₀ = 1.4 × 10⁻³μm² ps⁻¹. The parameters of the pump ring read: the radius R_P = R_in + (R_out − R_in)/2, with R_out = 15 μm and R_in = 13 μm, and σ_ring = 0.63 μm, corresponding to a FWHM of the laser of 1.5μm. For the relaxation coefficients we use β₀ ≈ 0.1 and β₁ ≈ 0.5. The weak flat pump at the centre of the trap has an amplitude parameter α_p = 0.1.

We simulate the dynamics of the polariton system by numerically integrating the stochastic differential equations Eq. (2) for the polariton field (see Supplementary Note 3). The numerical integration is performed on a two-dimensional lattice with Dirichlet boundary conditions. The lattice is composed of 128² grid points, with the lattice spacing Δx = 0.63μm. The lattice spacing is chosen to lie within the validity regime of the truncated Wigner methods used for the description of the stochastic field equations^15,55, but large enough to capture the macroscopic physics of the system. We notice that it also introduces a cut-off k_cut ∝ Δx⁻¹ in the momentum representation of the field. The integration of Eq. (2) in time is performed by using the XMDS2 software framework⁶². Specifically, we use a fixed time step which ensures stochastic noise consistency, and a fourth-order Runge-Kutta algorithm. The complex term dW has been implemented using the Wiener noise. We use a double precision SIMD-oriented Fast Mersenne Twister (dSFMT) algorithm to generate the random variable. We carefully check that dW is delta-correlated according to the formula given in the text above, by running numerical tests on the distribution of the noise which is let to evolve with ψ = 0, distributions analysed both in space and time.

Once the steady-state is achieved, to obtain a smooth first-order correlation function, g⁽¹⁾(r), we average over both the correlators at the x = 0 and y = 0 positions, using the expression ({g}^{(1)}(r)equiv [{g}_{x}^{(1)}(r)+{g}_{y}^{(1)}(r)]/2) where ({g}_{x}^{(1)}(r)={langle {psi }_{i,0}^{* }{psi }_{i+r,0}rangle }_{{{mathcal{N}}}}/scriptstylesqrt{{langle | {psi }_{i,0}{| }^{2}rangle }_{{{{N}}}}{langle | {psi }_{i+r,0}{| }^{2}rangle }_{{{{N}}}}}) and ({g}_{y}^{(1)}(r)={langle {psi }_{0,j}^{* }{psi }_{0,j+r}rangle }_{{{mathcal{N}}}}/scriptstylesqrt{{langle | {psi }_{0,j}{| }^{2}rangle }_{{{{N}}}}{langle | {psi }_{0,j+r}{| }^{2}rangle }_{{{{N}}}}}). Here, ψ_i,j = ψ(x_i, y_j) and the average ({langle ldots rangle }_{{{{N}}}}) is performed over the number ({{{N}}}) of stochastic realisations. All the results presented in our study are converged with respect to the number of stochastic realisations ({{{N}}}=100).