Critical fluctuations in a confined driven-dissipative quantum condensate

Phase fluctuations determine the low-energy properties of quantum condensates. However, at the condensation threshold, both density and phase fluctuations are relevant. While strong emphasis has been given to the investigation of phase fluctuations, which dominate the physics of the quantum system away from the critical point, number fluctuations have been much less explored even in thermal equilibrium. In this work, we report experimental observation and theoretical description of fluctuations in a circularly confined nonequilibrium Bose-Einstein condensate of polaritons near the condensation threshold. We observe critical fluctuations, which combine the number fluctuations of a single-mode condensate state and competition between different states. The latter is analogous to mode hopping in photon lasers. Our theoretical analysis indicates that this phenomenon is of a quantum character, while classical noise of the pump is not sufficient to explain the experiments. The manifestation of a critical quantum state competition unlocks possibilities for the study of condensate formation while linking to practical realizations in photonic lasers.


I. INTRODUCTION
The experimental realization of Bose-Einstein Condensates (BECs) enabled investigations of macroscopic quantum systems.One fundamental question relates to the role of fluctuations close to the critical point of a BEC transition.Critical quantum fluctuations arise from the uncertainty principle and play an essential role in a wide range of physical phenomena such as the structure of the universe, the Casimir effect and transitions between two competing phases of matter [1][2][3][4][5][6][7][8].Understanding such fluctuations is important for statistical physics of phase transitions and quantum critical phenomena.
Up to now, the focus has largely been on phase fluctuations, which are fundamental for understanding the physics away from the critical point [9].Yet, at the threshold of condensation, density fluctuations play a crucial role, while they have been much less explored.In the context of equilibrium systems, number fluctuations have been characterised to some extent in dilute atomic gases [10][11][12] and photon in a dye [13,14] condensates.For the out-of-equilibrium systems, however, the physics of critical number fluctuations has not been explored yet.
Here the scenario is more complicated; in addition to the usual density fluctuations in a single lowest energy state, we can witness a competition between different modes in the condensate formation process, which is particularly pronounced at vicinity of the critical point, where the condensate occupation is low.
Due to their strong nonlinear properties, microcavity exciton-polaritons (called here simply "polaritons") are the perfect candidates for the investigation of critical density fluctuations in non-equilibrium quantum systems.They are quasiparticles arising from the strong coupling of cavity photons and excitonic transitions, typically in semiconductor quantum wells (QWs) placed at the antinodes of an optical cavity.Polaritons can be viewed as photons with nonlinear interactions, or alternatively as excitons with very light effective mass (typically 10 −4 times the mass of an electron in vacuum).As interacting bosons, polaritons demonstrate quantum phenomena such as Bose-Einstein condensation [15,16], superfluidity [17] and quantized vortices [18][19][20][21][22][23] at temperatures from tens of Kelvin up to room temperatures.
In the present work, we undertake a detailed experimental and theoretical investigation of fluctuations near the condensation transition in a non-equilibrium system of non-resonantly pumped polaritons placed in an annular trap geometry.We explore theoretically the effect of quantum and classical noise on the system behaviour, and in particular on the mechanisms of the density fluctuations and mode competition.Apart from the fundamental interest, our results are also relevant for practical design of photonic lasers, where multimode switchingoriginating from their weakly-interacting nature-leads to "telegraphic noise" [24][25][26][27] instabilities.

II. RESULTS AND DISCUSSION
In order to investigate the fluctuations close to the critical point for condensation, we confined the polariton fluid in a circular trap, using a spatial light modulator to focus The pump laser acted as an incoherent source of excitons and polaritons, since it had photon energy well above the exciton energy, so that many phonons had to be emitted before the hot carriers cooled and turned into excitons and polaritons.This incoherent pumping therefore had two effects: it created both the potential energy profile of the trap, due to the repulsion of polaritons from the slow-moving excitons in the pump region, and it also acted as the source of the polaritons, as excitons converted into polaritons.As the pump power increased, the system went from non-condensed to a single-mode condensed regime, as reported in previous works.The large number of DBR periods gives the cavity a high Q-factor, resulting in a cavity lifetime of ∼135 ps and a polariton life time of ∼270 ps at resonance [28][29][30] (see Materials and Methods for a detailed description).
The sample was cooled in a continuous-flow coldfinger cryostat, which was held at a temperature of ∼5 K. To produce polaritons, the sample was excited nonresonantly with a continuous-wave (cw) laser, which was modulated by an optical chopper at 404 Hz with a duty cycle of 1.7% to prevent sample heating.The ring-shaped trap was generated by modulating the phase front of the laser using a spatial light modulator (SLM).
The diameter of the ring-shaped optical trap used in this manuscript was 45 µm, which is shown in Figure.S7(B).The cavity detuning from the exciton resonance was δ = −10.1 meV, corresponding to an exciton fraction |X| 2 = 0.3 for the lower polariton.The photoluminescence (PL) was collected using a microscope objective with a numerical aperture (NA) of 0.75, and was imaged onto the entrance slit of a spectrometer.The im-age was then sent through the spectrometer to a charged coupled device (CCD) for time integrated imaging.
Polariton condensates allow the study of a range of behavior as the detuning, and consequently the photon fraction of the polaritons, is changed.In the limit of high exciton fraction, the modes are closely spaced together in a near-continuum, and the interactions between the particles are strong, so that the dynamics are those of a close to equilibrium Bose-Einstein condensate [31,32].As the exciton fraction is decreased, the interactions become weaker and their mass lighter, so that they approach behavior that is more like a standard laser.The experiments we have performed here are at a middle ground; the polaritons are in the strong coupling limit, and travel tens of microns away from the pump region, but have resolvable transverse modes similar to those of a VCSEL.
We recorded several snapshots of the photoluminescence (PL) emission patterns from the trap for different pump powers.Each snapshot was time integrated for 100 ms.Since the pump laser was optically chopped, each snapshot is a sum of 42 images from separate laser pulses, each of which was about 42µs long.Crucially, the pump conditions are identical for all snapshots at each pump power, since we used a stabilized M Squared laser, with very low noise.Fig. S8(A) shows the density across the non-condensed-to-condensed transition as the pump power was increased; in Fig. S8(B) we show different snapshots of the PL in real space below, at and above threshold for condensation.Below threshold, polaritons remain in the vicinity of the pump spot and the PL looks spatially similar from one snapshot to the next as shown in Fig. S8(B), column I.For pump power above the con- densation threshold, the polariton population builds up inside the trap.When the density becomes larger than the critical density, stimulated scattering becomes significantly enhanced, giving rise to a condensate mode.Such modes have been seen in previous experiments with photonic detuning [33][34][35][36][37][38], and resemble pure lasing modes in vertical cavity, surface-emitting lasers (VCSELs) with imposed circular symmetry [24][25][26][27].
As expected, near the onset of condensation, we observe strong density fluctuations.Some of these fluctuations can be attributed to hopping of the condensate between different discrete eigenstates.Even though the pump conditions are identical from one snapshot to the next, the polariton gas near the threshold for condensation does not always condense to the same mode, but rather switches between several distinguishable modes as illustrated in Fig. S8(B), column II.Interestingly, although each snapshot is the result of summing 42 images from separate laser pulses, statistical clumping of certain modes -in small samples of 42-is able to reproduce the same mode switching effect (see Supplementary Information for more details).Notably, unlike optical multistability or bistability mechanisms, which can take place under resonant excitations at higher pump powers [39,40], the critical fluctuations we see are only present at the condensation threshold.Surprisingly, we find that the time scale of the fluctuations is very long compared to the intrinsic scattering time scale of the polaritons; features are seen to differ in images recorded with time integrations of 100 ms.This indicates that the system has metastable telegraph-type switching between modes, as have been seen for lasers [41][42][43].Here, this corresponds to the mode competition in the condensate formation process.In equilibrium BEC, it is the ground state of the system, which becomes macroscopically occupied.Since our system is driven-dissipative, and out-of-equilibrium, the process of the condensate formation is more complex.Which mode wins out to condense depends on the subtle interplay of gain and interactions effects.Eventually, as the pump power is increased away from the condensation threshold, the density fluctuations become less relevant and the polariton fluid occupies one mode for all different snapshots (see Fig. S8(B), column III).

A. Theoretical Modelling
Near the critical threshold for condensation, mean-field theories are known to become inaccurate since the fluctuations of the order parameter about its mean value become large, and the state of the many-body system is no longer accurately described by a single-particle wave function.In order to model the experimental results, we describe the strongly fluctuating polariton field within the truncated Wigner (TW) approximation [45,46] in terms of a stochastic equation coupled to a rate equation for the excitonic reservoir.The full model and further details regarding numerical methods can be found in the section Materials and Methods.
We solve the stochastic equations with a Wiener noise at each time step for many realizations of the initial noise.Both initial noise and the noise during the time evolution simulate quantum fluctuations.We emphasize that the amplitude of the Wiener noise in the system is not tuned externally, but rather, it is introduced self-consistently in the model.In order to simulate the experimental measurements, the polariton dynamics is left to evolve for a long time (up to 0.5 µs) after the Non-Equilibrium Steady-State (NESS) is reached.Then, each realization is independently integrated over different time-windows of the NESS.In Fig. S8 we report typical results from our model.While the behaviour of the polariton density is depicted in Fig. S8(C), typical steady-state realizations are plotted in Fig. S8(D), for three typical cases: below, near and above the critical point (Fig. S8(D), column IV, V and VI respectively).Similar to the experimental cases explored in Fig. S8(A-B), below and far above threshold of condensation, the quantum noise plays no role, producing an ensemble of similar images for each realization (Fig. S8(C), columns IV and VI).Close to threshold, however, the quantum noise becomes important, giving rise to realizations with different modes (Fig. S8(D), column V).

B. Measuring critical mode competition
To gain a quantitative description of the critical fluctuation features exhibited by the polariton system, we calculated the normalized average image difference between different consecutive images from the experiments for a given set of parameters such as pump power.This can be written as where I n is the nth intensity image in real space and RMS stands for the root mean squared.This average image difference is then normalized for each pump power by the average intensity for that set.We plot this quantity for both the experiment and the numerical simulations in Fig. 3(A) and (D), respectively.While relatively small below and above the critical point, the average image difference I d is found to show a maximum exactly at the threshold region indicating that fluctuations are large near the condensation threshold.In order to verify that this behavior is due to the onset of competition between different modes, we make use of angle-resolved imaging to access the information in momentum space [47].Such a procedure allows us to extract the number of modes, together with their energy, over the whole set of images.In Fig. 3(B) we plot energy hopping for powers at, and above, the critical threshold, revealing the single-mode to many-mode transition.The quantity I(k ∥ , E) is shown in Fig. 3(C) for two typical configurations: single-mode (panel Ia) and multi-mode condensate (panel Ib).
Our numerical simulations are in good agreement with the experiment; the numerical average image difference has a maximum near the threshold (point III in Fig. 3(D)), while it drastically decreases at lower and higher pump powers (point IV).We can extract the different mode energies by Fourier transforming the computed wavefunction ψ = ψ(r, t) to ψ = ψ(k, ω).Experimentally, we did not measure the energy fluctuations in real time.Instead, we extracted the mode energy for each snapshot, allowing us to plot the energy as a function of snapshot number.We can, however, access this information from our numerical simulations; transforming over different temporal windows, we extract the energy peaks for each spectra and plot the energy as a function of time (Fig. 3(E)).

C. Mode hopping during relaxation
Let us now focus on the mode hopping behaviour of the different cases shown in Fig. 3(E).At pump powers well above threshold (Fig. 3(E), panel IV), the system undergoes energy fluctuations and hopping during only the early dynamics of the system, and then settles down to a single mode.Such a behaviour is linked to the relaxation dynamics expected to take place in coupled condensate-reservoir systems.In previous works, the relaxation dynamics and the dynamical instabilities have been found to be a direct effect of the polaritonreservoir repulsive interaction [48][49][50][51][52], the latter being closely related to the polariton hole-burning effect [46].In Fig. 4(B), corresponding to the case above threshold (point IV in Fig. 3(D)), the strong density fluctuations take place up to times 20 < τ RD < 100 ns before eventually relaxing down.It is interesting to note that at higher pump powers (Fig. 4 relaxation dynamics with characteristic time of the order of a few nanoseconds.

D. Critical mode competition
At pump powers close to the condensation threshold the situation changes substantially.As depicted in Fig. 3(E), panel III the system spends time in three or more distinguishable states along its entire dynamics.Fig. 4(A) shows that the mode competition leads to to strong density fluctuations.In Fig. 3(F) we plot the I(k, E) spectra of two typical modes (further details on the spectral analysis of the mode competition are reported in the Supplementary Information).This numerical observation helps us to give a quantitative picture of the system's behavior at the critical point: different stochastic realizations (ini-tialized with different numerical seeds) lead to a very different and independent dynamics.Each stochastic path exhibits strong fluctuations in density and hopping between different modes and energies.It is important to note here that such fluctuations persist for timescales up to ∼ 20ns, which are orders of magnitude larger than the polariton lifetime, leading to persistent "intermittent" oscillations, which resemble telegraphic behavior observed in photonic systems [24][25][26][27].Importantly, we note that, as opposed to the configurations at larger pumps discussed in the previous section, the fluctuations persist over the whole NESS dynamics.
While one would expect that by time-averaging over long periods, such fluctuations would eventually be washed out, we find that imaging of mode competition is accessible even when integrating up to times comparable to the NESS duration (see Fig. S6 and discussion in the SI).The mode-competition effect can therefore be measured by making use of such "statistical clumping" of certain modes in small ensembles of measures; we discuss this in more detail in the Supplemental Information.We also investigated the effects of disorder on the mode hopping.Introducing a static disorder profile, with amplitude and spatial correlations matching the experiment, we observe an enhancement of the image difference I d .A more detailed discussion is reported in the Supplemental Information.
Finally, we note that in the experiment the measured time scales, of the order of a few microseconds, are longer than the ones observed in the theoretical analysis.We find that such a discrepancy can be explained by introducing low frequency pump modulations, mimicking external perturbations that are intrinsic to the experimental setup, but not accounted in the earlier theoretical modeling.We analysed how such time scales depend on external perturbations; our results are discussed in the Supplementary Information.We introduced small periodic modulations to the pump, of the order of the timescale of the numerical mode hopping (∼ 20ns).We observed that the switching time scales can be extended up to the total duration of the numerical dynamics, i.e. on the order of microseconds, therefore matching the switching timescales of the experimental measurement.

E. Quantum fluctuations vs classical noise
We now discuss the role of quantum noise in the onset of critical mode competition.We can see the importance of quantum fluctuations by comparing the numerical simulations in the beyond-mean-field TW picture with the same model treated at mean-field (MF) level.In the latter case, to account for classical variations in the driving strength, we include real-valued spatiotemporal fluctuations to the driving profile.Details of both models are reported in the Material and Methods section.
We can calculate the quantity I d from the outcomes of the two models, which are plotted in Fig. 3(D) as a blue solid curve (TW) and a pink dashed curve (MF).
Let us now restrict ourselves to the critical region, corresponding to the grey shadowed area in Fig. 3(D).While the TW (blue) curve exhibits a large peak in the vicinity of the critical point, which results from the mode competition physics discussed in the previous sections, the MF (pink) curve shows clearly an absence of such a peak.The peak near the threshold for condensation can be reproduced only by the self-consistent quantum fluctuations accounted within the TWA.Inspection of the density profiles in the two cases confirms the absence of mode hopping in the latter case.We demonstrate here that fluctuations in both density and phase are essential to explain the critical mode hopping; "classical" amplitude fluctuations in the pump strength alone are not sufficient.While we cannot account for all potential sources of classical fluctuations (e.g.phonons, high energy excitons), we believe that they are orders of magnitude weaker than the quantum noise coming from photon decay and pumping (see Supplementary Information for further discussion).

F. Dynamical instabilities at larger pumps
While we have focused so far on the critical region, it is interesting to also comment on the physics taking place at higher pump powers (i.e., higher polariton densities).In Fig. 3(D), we note the presence of a second peak at P /P TW th ∼ 1.92, i.e. point V. Our investigation of the polariton and reservoir density profiles suggests that this peak originates from the mechanism of dynamical instability [53,54].As discussed in more detail in the Supplementary Information, these can be reproduced by moving to a mean-field formulation of Eq. ( 2)-(3), while adding a classical-noise variations to the pump profile (pink dashed curve of Fig. 3(D)).Interestingly, these instabilities are found to be responsible for the large peak in the mean-field (MF) curve too.At very large pump powers (P /P TW th ∼ 3), the quantity I d eventually decreases.The density profiles shown in Fig. S8(D) reveal that at point VI in Fig. 3(D), the image difference is due to strong density fluctuations extended over the whole confined region; a different behaviour when compared to the critical mode competition observed at point III.We expect that such density fluctuations would eventually smooth out -and resemble more the experimental profiles-when considering multiple time-windows of the order of tens of microseconds as in the experiment.

III. CONCLUSIONS
In this work, we investigated experimentally and numerically the density fluctuations in a circularly-confined non- equilibrium Bose-Einstein condensate of polaritons.As in some experiments with mode hopping in lasers, at the vicinity of the critical region we observe "telegraphic" switching with long periods of metastability.While the phenomena of mode-hopping is a standard feature in photonic systems such lasers and VCSELs [24][25][26][27], in this work for the first time we are able to demonstrate that, in the photonic limit of an interacting polariton fluid, the mode competition arises from the interplay between interactions and quantum noise.Our results suggest that, analogously to other activation problems, quantum fluctuations are responsible for switching between different states.
In addition, we distinguish between the mode competition at the critical point, and the relaxation mechanisms expected to take place at larger drivings.We explore in detail the effects of external perturbations and the role of dynamical instabilities at larger pump powers.Most importantly, our numerical analysis demonstrates that beyond-mean-field quantum effects are crucial for the onset of mode competition and that it cannot be replicated by a model with just classical noise (i.e., fluctuations in the pump power).Such findings are determinant for understanding of critical density fluctuations in confined quantum systems.Our work paves the way to further investigation of the mechanism of condensate formation in both equilibrium and non-equilibrium settings.

A. Sample design
The microcavity used in this work consists of a total of 12 GaAs quantum wells (7 nm thick) with AlAs barriers embedded between two distributed Bragg reflector (DBR) mirrors, which together form an optical cavity (Fig. S7(A)).The DBRs are made of alternating layers of AlAs and Al 0.2 Ga 0.8 As.The top DBR is composed of 32 pairs and the bottom DBR is composed of 40 pairs.The quantum wells are in groups of 4, with each group placed at one of the three antinodes of the 3λ/2 cavity.The large number of DBR periods gives the cavity a high Q-factor, resulting in a cavity lifetime of ∼135 ps and a polariton life time of ∼270 ps at resonance [28][29][30].The long cavity life time allows polaritons to propagate over macroscopic distances of up to millimeters.The wavelength of the pump laser was tuned to the second reflective minimum (725 nm), about 113 meV above the lower polariton resonance.

B. Numerical model and parameters
Within the TW approximation [45,46], the strongly fluctuating polariton field can be described in terms of a stochastic equation of motion coupled to the rate equation for the excitonic reservoir, which read where ψ(r, t) is the polariton field, n R (r, t) the excitonic reservoir density and m the polariton mass.γ c and γ R define the decay rates of condensed polaritons and the excitonic reservoir respectively.Following Ref. [46], the former term can be estimated as γ c = (1−|X| 2 )γ ph , where γ ph = 1/τ ph corresponds to the inverse of the photon lifetime and X is the Hopfield coefficient [55].The excitonic fraction where Ω is the Rabi splitting and the detuning δ corresponds to the difference between the photonic and excitonic energies.In Eq. ( 2) the Wigner noise accounts for quantum fluctuations and is space-and time-correlated as ⟨dW c (r, t)dW c (r ′ , t) ⟩ = 0, ⟨dW c (r, t)dW * c (r ′ , t)⟩ = (γ c + Rn R (r))/2 δ r,r ′ dt.Moreover, the renormalized density |ψ| 2 W ≡ |ψ| 2 − 1/2a 2 includes the subtraction of the Wigner commutator contribution, where a is the numerical 2D grid lattice spacing.The relaxation parameter β sets the amount of energy relaxation and increases as the detuning δ = E c − E ex increases.The constants g c and g R are the strengths of polariton-polariton and polariton-reservoir interactions respectively; they can be estimated as g c = g ex |X| 4 , g R = g ex |X| 2 , with g ex the exciton-exciton interaction.The parameter R = R 0 g c /g R quantifies the stimulated scattering rate of the reservoir excitons into the polariton condensate [46].In order to match the experimental geometry, the pump power profile P r (r

C. Mean-field with classical fluctuations in the pump
In the main text we discuss the role of quantum fluctuations in the mode competition by comparing the results obtained from the TW equations, i.e.Eq. ( 2) and Eq. ( 3), with their mean-field approximation and only a classical noise (introduced as fluctuations in the pump power).
The polariton equations of motion in the mean-field approximation read: To explore the effects of classical fluctuations in the driving, we introduce spatial-temporal variation in the pump power: where dW(r,t) is a zero-mean, real-valued Gaussian noise both in space and time.This formulation excludes the non-classical correlations present in the Truncated Wigner approximation.Here, the term χ quantifies the amount of classical fluctuations on top of the driving pump P r (r).
In typical experimental setups, the strength of the classical noise is usually ≈ 1% of the pump intensity.In our numerical analysis, we push the strength of the noise up to 10% of the total intensity of the pump, corresponding to χ = 0.1.Even at these strong classical noise amplitudes, we do not see any peak around the threshold P TW th ≈ 1 on the I d purple curve plotted in Fig. 3(D), confirming that the I d peak measured at criticality is due to the presence of the Wigner quantum noise included in the TW approximation, as discussed in the main text.We refer the reader to the Supplementary Information for further numerical details.
As explained in the main text, experimentally, we did not measure the energy fluctuations in real time.Instead, we took many snapshots of the polariton gas using angleresolving imaging for each pump power and obtained the mode energy for each snapshot image I(k, E).This was extracted by first summing the image I(k, E) horizontally (i.e integrating over k) to obtain the intensity I(E), and then numerically estimating the energy for each I(E) peak as shown in Fig. S1 (B).The columns in Fig. S1 from left to right relate to below, near and above the threshold for condensation (which correspond to points I, II, and III in Fig. 2 in the main text).We then plot the extracted energy as a function of the snapshot number (Fig. S1 (C)).To help visualize the frequency of each mode, we use the points in Fig. S1 (C) to create a histogram (Fig. S1 (D)).We find that the polariton gas is multimode near the phase transition and single mode away from the phase transition (i.e above and below the threshold for condensation).
We follow a similar procedure for the numerical data, which allows us to extract the energy fluctuations (Fig. S2).Fig. S2, columns II, III, V, and VI correspond to points III, IV, V and VI of Fig. 3 in the main text, respectively.These points are discussed in the main text while here we include the full data for completeness.We first integrate the real space solution ψ(r, t) over a 20ns temporal window.The real space wavefunction ψ(r) is then Fourier transformed to ψ(k).|ψ(r)| 2 and its Fourier transform |ψ(k)| 2 are shown in Fig. S2(A) and Fig. S2(B) respectively for the last 20 ns temporal window of the dynamics.To find the energy resolved spectra, we Fourier transform the wavefunction ψ = ψ(r, t) to ψ = ψ(k, ω) over the same 20ns temporal window.We then follow the same method mentioned earlier to extract the energy fluctuations by detecting the energy peaks from the numerical I(k, E) images (see Fig. S2(D-F)).We find that below the condensation threshold there is a continuum of energy states (see Fig. S1 (E), column I).Near the critical regime, the condensate switches between three main states with almost equal distribution * Address correspondence to: haa108@pitt.edu(see Fig. S2 (F), column II).Importantly, these critical fluctuations persist over the total simulated time as shown in Fig. S1 (E), column II.Well above threshold, the condensate undergoes switching during early dynamics and then settles down to a single mode (Fig. S1(E), column III).Since the switching happens only during early dynamics, the histogram shows a single dominant mode as illustrated in Fig. S1(F), column III.Fig. S1(E), columns V and IV show the effect of the dynamical instability (a term which refers to a spatially fragmented condensate, caused by effective attractive interactions, mediated by the excitonic reservoir [56,57]) taking place at higher pump powers.In these cases, the condensate density starts to build up in the centre of the ring, over the whole confined region.What we see then are the spatial fluctuations in the condensate density (spatially non-uniform/fragmented condensate) in a single energy mode.This effect is different from the one observed at criticality; a detailed explanation of this phenomenon is given in the next sections.

II. STATISTICAL CLUMPING
The experimental emission patterns in Fig. 2 of the main text are a sum of 42 images from separate pulses, while in the simulations each realization corresponds to a single pulse.In this section, we show that the same mode hopping effect can be reproduced by statistical clumping of certain modes in small samples of 42.To model the statistical clumping, we first start with three modes.We then randomly choose one of these modes N pulse times and compute the sum of these randomly chosen N pulse modes.This sum of the N pulse images is therefore equivalent to a snapshot in the experiment (in the experiment, N pulse = 42).Experimentally, we recorded 300 snapshots for each pump power.Therefore, we repeat the above process in the simulations to generate 300 snapshots.We then use Eq. ( 1) from the main text to compute the average image difference of these snapshots.To illustrate the significance of statistical clumping, we vary N pulse (i.e the number of images that are being summed over) and we plot the image difference as a function of N (see Fig. S3).As seen in Fig. S3, statistical clumping of three modes in small sample of 42 (the red circle shown in the plot) shows a noticeable image difference.As expected, when the number N pulse becomes very large, the image difference tends to zero.However, the experimental regime is in small samples of N pulse = 42 where mode-hopping can still be observed.We find that the image difference I d gets reduced by roughly 0.2 when N pulse = 42 compared to the case of N pulse = 1.Therefore, the effect of summing over images from separate pulses reduces I d but the overall mode hopping effect can still be observed in a small sample of N pulse = 42.
Although the results we have presented in Fig. 2 of the main paper show the condensate in only a single mode, very often we also observe the polariton gas in a linear combination of these modes.However, since we take many snapshots of the condensate, occasionally we also observe it in a single mode even after averaging over 42 chopped pulses.Figure S4 shows an example of polariton condensate in a single model (left and middle panels) and a linear combination of these modes (right panel) for the same pump power P /P th = 1.

III. ENERGY SPACING
In the experiment, we observe that energy spacing between the switching modes of the system is about 200 µeV , while in the simulations, the observed energy spacing is 2 µeV .In this section, we show that in the simulations the switching happens between states with energy close to the ground state of the system, while experimentally the switching occurs between higher energy    modes.
To find which of the single particle states the condensates chooses, we solve the Schrodinger equation without nonlinearities in 2D with an annular potential.
We then compare which of the solutions |ψ n,m (r r r)| 2 look spatially similar (i.e the same number of pedals or ripples) to the experimental and theoretical real-space images.The eigenstates are calculated using potential of the form V (r r r) = V 0 exp [−(r − r 0 ) 2 /2σ], where r 0 is the radius and σ is the width of the potential ring.We solve Eq. (S1) numerically to find |ψ n,m (r r r)| 2 and E n,m .
The quantum number n determines the number of nodes in the radial direction, while m determines the nodes around the ring.The number of radial nodes is given by n − 1 and the number of nodes around the ring is given by 2m. .We find that in the experiment the condensate switches between three main modes ψ (n=1,m=16) , ψ (n=5,m=1) and ψ (n=3,m=16) , which are shown in Fig. 2, column II of the main text.In the simulations, the switching happens between modes ψ (n=1,m=8) , ψ (n=2,m=8) and ψ (n=1,m=1) .The energy spacing is roughly proportional to z (m final , n final ) 2 − z (m initial , n initial ) 2 , suggesting that switching between higher energy modes will give larger energy spacing, where z (m, n) is the n-th zero of the regular Bessel function J m (z).Of course, the real-space images of the condensate both in the experiment and the simulations include nonlinearities.However, the solution to the linear Schrodinger equation still allows us to find the quantum number of the single particle states, which participate in condensation.The states obtained from solving the Schrodinger equation are close to the real space distributions obtained from the experiment and the simulations with only small deviations due to relatively weak polariton interactions.
In our numerical simulations, we find that the modes, which participate in the mode switching process, depend on various parameters of the model such as the relaxation and interaction strengths.The values of these parameters are not exactly known for semiconductor microcavities [58].Exploring in principle a large range of possible parameter space is both time-consuming and computationally expensive.We, therefore, focus on reproducing the experiment qualitatively i.e. showing that the mode-hopping is pronounced at the phase boundary and is greatly suppressed away from it.

IV. DENSITY CALIBRATION
In Fig. 2A and Fig. 3A of the paper, we have used the photon counting method to calibrate the density of the polaritons.Photon counting allows us to relate the number of counts on the CCD camera to the number of photons detected.This was done by matching the laser wavelength to the polariton emission wavelength (776.5 nm).The laser was sent to a mirror at the sample plane, which reflects the laser through the same optical path that was used in the experiment.This reflected beam was then imaged with the CCD camera.The CCD count is proportional to the number of photons, which can be written as: where ∆t is the integration time of the camera, N ph is the number of photons detected and η is the efficiency factor.The number of photons that the camera detects during a time ∆t is then: where P is the measured power of the laser, h is Planck constant, c is the speed of light and λ is the wavelength of the laser.This then allows us to find a single efficiency factor to convert from CCD counts to number of photons being detected.The efficiency factor is given by: This efficiency factor was used to calibrate the total density of the polaritons, which is given by: where I CCD is the CCD count, τ is the average radiative lifetime of the polaritons and A obs is the observed area on the sample from which the light was collected.Since the excitation laser beam was chopped in the experiment with a duty cycle d = 1.7%, then the total number of polaritons is where Here |C k ∥ | 2 is the photon fraction and τ cav is the cavity lifetime.

V. NUMERICAL ANALYSIS
In the main text we introduce the numerical model used to simulate the experimental observations, and discuss the mode switching near the phase transition.In this section, we give more details about the numerical analysis; we first discuss the model and the results from the mean-field equations with noise added to the pump profile (Sec.V A), we then explain the effects of the system relaxation dynamics on the mode competition (Sec.V B), and finally we explore introducing classical external perturbation simulated by a periodic pumping (Sec.V C).
A. Numerical details of the mean-field modeling with a noisy pump profile.
In the main paper, we investigate the role of quantum critical fluctuations on the observation of mode-hopping features in the polariton system by comparing different models.The mean-field modeling with a fluctuating pump profile is described in the Methods Section of the main text.In this section, we give further numerical details regarding the simulations of this model.
In Fig. S5 we plot the results obtained by solving Eqs. ( 4)-( 5) of the Methods Section in the main text.The nonequilibrium steady state is reached after evolving the system from a random initial condition for 1µs.We average over N = 100 realizations.The time evolution of the total averaged density ⟨|ψ(t)| 2 ⟩ N is plotted in the left panel of Fig. S5 for different pump values.Once the system reaches a nonquilibrium steady state, we extract the average steady-state density shown in the top right panel of Fig. S5.We then calculate the image difference I d by means of Eq. (1) of the main text and plot it in the bottom right panel of Fig. S5, and in Fig. 3(D) of the main text.

B. Effects of dynamical relaxations and instabilities on the mode-hopping.
We proceed by discussing the effects of the relaxation mechanisms on the mode swithing.As anticipated in the main text, and shown in Fig. S2, in certain regimes above threshold (P TW th ≈ 1.5) the system is seen to undergo swithing between different energy modes during only the early dynamics of the system, before settling down to a single mode.Comparison of the time evolution of the polariton and reservoir densities at this particular pump power regime (shown in Fig. (4) of the main text) suggests that this early-time mode competition originates from the relaxation dynamics in the condensatereservoir coupled system [48][49][50][51]46].As shown in Fig. (4) of the main text, different stochastic realisations (initialized with a different numerical seeds) lead to different dynamics.Interestingly, both at criticality and at early times of the polariton dynamics above threshold, they show "intermittent" oscillations, which resemble a telegraphic behavior.Such intermittency may be an effect of introducing the relaxation term β in the equations of motions (Eqs.( 2)-(3) of the main text): preliminary results on the observation of intermittent behavior has been reported in uniform systems [59].We leave the numerical investigation of this behavior, for the case of confined systems, for future work.
To investigate further the effect of the relaxation dynamics on the mode switching over the whole phase diagram, we extract the quantity I d for different integrationtime windows, and investigate real-space density distributions of the polariton condensate.In Fig. S6, we plot three different I d , calculated within the TW approximation, by time-integrating over different temporal windows, from different starting times till the end of the evolution, i.e. t = 0.5µs.Specifically, we integrate be-tween: 0.02µs < T < 0.5µs (red curve); 0.05µs < T < 0.5µs (black curve); 0.2µs < T < 0.5µs (blue curve).The latter corresponds to the blue curve Fig. (3) (D) of the main text.We note that at P /P TW th ≈ 1.6, a large peak appears when considering longer times of integration.Comparison with the density evolution plots in Fig. (4) of the main text suggests that this peak is due to the inclusion of the switching between modes originating from the relaxation mechanisms.
Next, we compare the TW approximation results in Fig. S6 with the image difference I d , calculated using classical noise added to the pump profile in the MF approximation.(This is shown in Fig. S6 as a pink curve; the same curve is reported in Fig. 3 (D) of the main text.)Although the MF result does not exhibit a peak at criticality, as discussed in Sec.V A, it is interesting to note the presence of a maximum at larger pump powers.An investigation of the real-space density images along the time evolution suggests that this is due to the so called modulational instability [51], which is a mean-field effect.Attractive interactions between polaritons, caused by the presence of the reservoir, lead to the condensate fragmentation and so fluctuating density, which persist over the whole dynamics as shown in the left panel of Fig. S5.This indicates a possible relation between the observed MF peak and the peak at P /P TW th ≳ 1.8 observed in the TW formulation (blue curve).

C. Modelling classical external perturbation with periodic pumping
In the experimental setup, external perturbations e.g.temperature variations or laser modulations may have a large impact on the critical timescales of the system; such perturbations are not included in our model.In this section we discuss the role of the dynamics of the external pumping on the switching time-scales of our theoretical modelling.Specifically, our analysis shows that the timescale of the mode switching can be extended up to the whole dynamics (of the order of microseconds) by introducing small periodic modulations of the amplitude of the pump.We model this phenomenon by means of a periodic pump defined as For our analysis, we choose timescales comparable to the temporal mode-hopping scales of the numerical simulations, namely τ mod = 20 ns.To probe different regimes of modulations, we vary A in the range A ∈ [0, 10 −2 , 10 −1 ].
Investigation of the real-space density distributions shows that a choice of A = 0.01 extends the switching time scale to over the whole time-dynamics, namely 0.5µs.Moreover, comparison of the quantity I d between modulated and non-modulated cases (A = 0) shows no appreciable differences (< 10%) in the critical region.Indeed, the polariton density follows the modulation of the pump, yet the density profiles do not show large changes when compared to the unperturbed case.To test our results, we do numerics for values up to A = 0.1.An investigation of the polariton density distribution shows very large oscillations; moreover, the image difference increases substantially compared to the non-oscillating case.

D. Effect of static disorder on mode switching
We have studied the effect of static disorder on mode switching by introducing a time-independent term V (r) to the Hamiltonian of Eq.( 2) of the main text.We assumed that the random potential has a mean amplitude and root mean square fluctuation given by ⟨V (r)⟩ = 0 and √ ⟨V 2 (r)⟩ = 10 µeV respectively.The correlation length of this potential is 2 µm.We have run 50 different stochastic dynamics with the same random static noise profile but with different random noise in the initial wave function.We then averaged over these 50 different stochastic dynamics and repeated the procedure for 50 different random static noise profiles.
We have found that the static disorder enhances the image difference I d by a factor of 4.However, we note that the number of modes that undergo mode hopping is still the same when compared to the case with only quantum fluctuations and no static disorder (i.e the results in the main text shown in Figure 3).Therefore, we conclude that static disorder only introduces density variation of the polariton condensate rather than affecting the number of modes.

E. Convergence Analysis
In this section, we discuss the numerical convergence in the spatial grid spacing a.The TWA method, employed in this work, is computationally very efficient and is able to describe quantum fluctuations.However, TWA has a relatively narrow window of applicability: it can be shown [45] that TWA is valid only for lattice spacings greater than the square root of the ratio between interactions and losses: γ LP ≫ g LP /a 2 .This condition essentially sets a lower limit for a.This originates from the derivation of TWA equations, where one assumes the 3rd order derivative in Fokker-Plank equation to be small (at least smaller than the 2nd order term which we keep) so that can be ignored [45].
In our convergence and numerical validity tests, we kept the same set of parameters as in the main text and only varied a and pump power P .We have run 50 stochastic realizations for different spatial spacing, namely a = 0.8, 1.02, 1.1, 1.17, 1.55 µm.We have then investigated the time-averaged density profiles for three different pump power values, those corresponding to points III, IV and VI in Fig. 3(D) of the main text.
For the pump power corresponding to point IV (V) in Fig. 3, we find that the condensate is of a single mode (multi-mode) for all a values.At the critical point (i.e., the pump power corresponding to point III), we observed that mode switching only happens for a = 1.1, 1.17 µm, while for the small grid spacing (a = 0.8 µm and a = 1.02 µm), and for the larger grid spacing (a = 1.55 µm), we observed a single-mode condensate in each stochastic realization.
For the parameters considered in our work, the TWA condition reads a ≫ (g LP /γ LP ) 1/2 = 0.85 µm.a = 0.8 and a = 1.02 can therefore be excluded.At large grid spacing (i.e a ≥ 1.55 µm), the spatial discretization becomes too coarse for the typical scales of the system.Due to a too small high-momentum cut-off (k max = π/a), momenta contributions that are still significant in the system are suppressed.

VI. QUANTUM FLUCTUA-TIONS VS. CLASSICAL NOISE
Let us first define what we mean by quantum fluctuations.The definition of quantum noise which we use in this study is the same as in the seminal book by Crispin Gardiner and Peter Zoller [60].Fluctuations in quantum system arising from interactions with the external world leading to either drive or dissipation are defined as quantum noise.In contrast, fluctuations caused by finite temperature are referred to as classical noise.In the truncated Wigner approximation (TWA), quantum fluctuations of the open system manifest themselves as an additive white noise in space and time present at every time step of the dynamics.This arises from the 2nd order derivative term of the Fokker-Planck (FP) equations from where TWA is derived.It has to be mentioned that TWA has also been used, especially in the context of cold atoms, to describe finite temperatures [61].But, as it can be found in literature, in closed system the 2nd order term is missing in FP equations, and so there is no additive dynamical noise.Instead, the finite temperature fluctuations are encoded in the initial conditions.The dynamical equations remain mean-field and the thermal fluctuations are incorporated in the random/noisy initial conditions, where the observables are computed by averaging over different initial conditions.
Since our system is open (driven/dissipative) and weakly interacting, it is not in thermal equilibrium and the temperature is not defined.If we were to encode some temperature in the initial condition, as in cold atoms, this would have no effect on the steady state as due to continues drive and dissipation (i.e noise at every time step) the influence of any initial conditions will be washed out.Our steady state is initial condition independent.
Based on timescales, we have identified two realistically relevant processes that can induce the mode hopping.The quantum noise acting on photons due to their finite lifetime in the cavity, and the potential variations in the pump power (classical noise).Note that the real-valued noise corresponding to spatio-temporal variations in the pump gives fluctuations in the density only.There are no fluctuations in the phase of the condensate.What we can say for certain is that we need fluctuations in the phase to explain the experimental results.Fluctuations in density alone are not sufficient.That is why we believe that fluctuations in pump power are not the cause of the effects described here.There potentially could be another source of noise which gives phase fluctuations (external to polariton population) such as phonons or high energy excitons.However, we show in the next sections that interactions with phonons or high energy excitons are orders of magnitude weaker than the dissipation coming from the finite photon lifetime.

A. Interaction with phonons
Modeling thermal fluctuations within the G-P model is a difficult task.To show polariton-phonon interaction being weak, we take a different approach to estimate the polariton-phonon interaction, namely a quantum Boltzmann equation.To simulate the dynamics of the polaritons, we make use of the semiclassical Boltzmann equation, which reads: Here n ⃗ k is the occupation number, τ ⃗ k is the characteristic lifetime, P ⃗ k is the pumping term.The W (i) are the interaction terms for particle-particle collisions.The laser generation is modeled by using a time-independent pump term P that is longer than the total simulated time (i.e c.w pumping).Since the laser generation is nonresonant with energy much greater than the polariton energies, we assumed that the free electrons and holes created in the pump process lead each polariton state being pumped with equal probability.We solve Eq. (S8) numerically by including polariton-phonon interaction and polariton-polariton interaction to find n ⃗ k .The updated n k (t) is then used to find the new W ⃗ k→ ⃗ k ′ (t) until a steady-state distribution is reached.The interactions terms for the polariton-polariton and polaritonphonon interactions are given in Ref [62].The occupation number of the polartions obtained from Eq. (S8) as a function of their energy is shown in Fig. S7 Figure S7 shows several important aspects about the exciton-polaritons in general.First, the three regions of the energy distribution have very different properties.At high energy, which corresponds to the excitonic range of the spectrum, there is a thermal tail which fits a Maxwell-Boltzmann distribution since these high energy excitons interact very efficiently with phonons.The temperature of this distribution is very close to the lattice temperature.However, at low energy, which corresponds to the Figure S7: The occupation of the lower polaritons as a function of their energy.Red: when only the polaritonphonon interaction is considered.Blue: when both the lower polariton-phonon and polariton-polariton interactions are considered.
. polaritonic region, the polaritons have a nearly flat distribution when only the polariton-phonon interaction is included and therefore do not have a well defined temperature since they interact very weakly with the lattice.For this reason, we believe that the dissipation coming from the finite photon lifetime is orders of magnitude larger than anything connected with interactions with phonons.When the polariton-polariton interaction is included, the polaritons have a nearly thermal distribution with a temperature well above the lattice temperature.What we find is that polaritons thermalize with each other and not with the lattice.This is consistent with Quantum Boltzmann solutions previously reported in the literature [62].
To be more concrete about the strength of this polariton-phonon interaction, we calculate the polaritonphonon out-scattering rate per particle in a momentum state k (Ref.[63], Chapter 4): where ± correpsond to phonon emission (+) and phonon absorption (−).Here n ⃗ k is the occupation number of the lower polariton and n phon ⃗ q is the occupation number of the phonons, which we assume are in thermal equilibrium with n phon ⃗ q given by the Planck distribution n phon ⃗ q = 1/ (e ̵ hωq/k B T − 1).The polariton-phonon interaction is based on hydrostatic deformation potential, which takes the form [64], Figure S8: The scattering time of polaritons with phonons as a function of the polariton energy calculated from Eq. S9 .
where V , ρ, u are volume, density and longitudinal sound velocity respectively.X k is the Hopfield coefficient and a e and a h are the deformation coefficients of the conduction and valence band for GaAs respectively.I

⊥(∥)
e(h) are the overlap integrals between the exciton and phonon mode.The scattering time is plotted in Fig. S8.The scattering time in the polaritonic region is in the order of 10 −6 seconds, which is negligible compared to the photon decay which is in the order of a few hundred picoseconds.The scattering time with phonons is four orders of magnitude weaker than the photon decay.In the excitonic region, the scattering time is in the order of tens of picoseconds, which is why they can efficiently thermalize with phonons.

B. High energy excitons
In the experimental setup, the high energy excitons stay mostly near vicinity of the pump spot since they are 10 4 heavier than the polaritons.The diffusion length of these high energy excitons is typically ⩽ 1 µm while the diameter of the annular trap is 45 µm.We therefore expect that the density of these high energy excitons inside the trap where the condensate is forming is nearly zero.The bottleneck excitons have some polaritonic character and do not exhibit the nearly stationary behavior [65]; this poses an interesting question, however modeling such excitons is difficult and beyond the scope of this work.

Figure 1 :
Figure 1: Optical trapping using annular pump in a GaAs microcavity.(A) Microcavity structure used in the experiment.The dark and light blue layers represent the DBRs used to confine the light field.The yellow layer indicates the quantum wells, where bound excitons are formed and the red lines indicate the intensity of the confined optical field.(B) Normalized real-space image of the excitation laser generated by a spatial light modulator.(C) A slice of the intensity of the excitation ring at y=0 line through the center of the excitation ring pattern shown in (B).

Figure 2 :
Figure 2: Critical mode competition and density patterns in experiment (Left panels (A,B) and theoretical modeling (Right panels (C,D)).Panels (B,D) are selected snapshots of the polariton density distributions, where each column corresponds to a fixed pump power below (I column), near (II column) and well above (III column) the condensation threshold.The exact P /P th for each column is indicated by the black circles in (A) and (C).In panel (D), the single-realization density profiles are integrated in time after the system has reached the non-equilibrium steady state with the integration period of τ int = 20ns.

Figure 3 :
Figure 3: Image difference, energy hopping and spectra from experiment (left(A,B,C)) and theoretical modelling (right (D,E,F) panels).(A) Normalized average image difference I d given by Eq. 1 (blue circle and blue solid curve) and the polariton density (gray solid line) as a function of P /P th .(B) The extracted energy peaks for each snapshot, (I) near and (II) far from the critical point.In some cases the condensate is single-mode (purple dots), while in others is multi-mode (blue dots).(C) Representative I (k, E) snapshots from which the energy peaks were extracted.We show two examples for a single (Ia) and a multi-mode state (Ib) close to criticality.(D) The same as in panel (A), but from theoretical mean-field (MF) and Truncated-Wigner (TW) analysis.Polariton densities are plotted with dashed (solid) gray line for MF (TW).The average image difference I d is shown as magenta (blue) circles and magenta (blue) solid line for MF (TW).In panels (A) and (D) the critical region is marked as a gray shaded area; consistent with previous works [44].The threshold power was identified where the plot changes slope in a log-log scale.(E) The extracted energy peaks along the time evolution of a single realization of the model (2-3), within (III) and above (IV) the critical region.(F) The same as panel (C), but for the theoretical analysis showing two time-integrated snapshots with different number of modes.
modeled as a tilted elliptical geometry, with a rotated axis of angle θ.E is the tilted-ellipse equationE = Ax 2 + Cy 2 + Bxy − 1 where A = (cos 2 θ)/a 2 + (sin 2 θ)/b 2 , B = 2 sin θ cos θ and C = (sin 2 θ)/a 2 + (cos 2 θ)/b 2 ,with a and b the major and minor axis respectively.The model (2)-(3) is solved numerically adopting an explicit Runge-Kutta method of orders forth and fifth on a two-dimensional numerical grid with N = 128 2 points with grid-spacing a = 1.17µm.We solve the dynamical equations with the following experimental parameters: m = 4.2 × 10 −5 m e with m e the electron mass,

31 Figure S1 :
Figure S1: Experimental energy analysis (A) Representative I(k, E) snapshots from which the energy maxima were extracted.(B) I(E) obtained by horizontally summing I(k, E).The red circles illustrate the energy peaks extraction.(C) Switching dynamics for each snapshot.(D) A histogram of the occurrence of modes with given energy.The first column corresponds to below, the second to near, and the last column to above the condensation threshold.

Figure S2 :
Figure S2: Theoretical energy analysis.(A) Real space |ψ(r)| 2 and (B) momentum space |ψ(k)| 2 densities.(C) |ψ(k, E)| 2 from Fourier transforming ψ = ψ(r, t) to ψ = ψ(k, ω).In (A-C) we have integrated over the last 20ns of the dynamics.(D) I(E) from integrating I(k, E) over k.The red circles illustrate the energy peaks extraction.(E) The energy peaks as a function of time.We periodically sample different temporal windows of 20ns.Colors of the dots correspond to: single mode (purple); two modes (blue); three and more modes (green).(F) A histogram of the number of realisations within given energy interval.The red Gaussian curves in the last two histograms are guide to the eye.

Figure S3 :
Figure S3: Statistical clumping of modes.Normalized image difference by randomly choosing 3 modes N pulses times, and summing over them.This process is repeated 300 times, each corresponding to a single snapshot.We use Eq.(1) in the main text to compute I d .The red circle at N pulses = 42 shows the regime we work with experimentally.

Figure S4 :
FigureS4: Linear combination of modes.Multiple snpashots of the polarion density profile in real space taken at the threshold power.In the left and middle panels, the condensate is single mode while for the right panel the condensate is in a linear combination of these two modes.

Figure S5 :
Figure S5: Mean-field approximation.(left panel) Time evolution of the total density ⟨|ψ(t)| 2 ⟩ N from Eqs. (4)-(5) of the Methods and Material section in the main text evraged over N = 100 realizations, plotted in a linear-log scale.Different colors are for different pump powers as in the right panel.(top right panel) The extracted steadystate densities as a function of pump power normalised to the condensation threshold P /P TW th on a log-log scale.The same as the gray point-dashed curve in Fig. 3 of the main text.(bottom right panel) Normalized image difference I d in a log-linear scale obtained by time-integrating the dynamics shwon in the left panel over the interval 0.5µs < T < 1µs.

Figure S6 :
Figure S6: Image difference.The image difference I d , calculated within the Truncated-Wigner (Eqs.(2)-(3) of the main text) marked as solid blue and mean-field approximations (Eqs.(4)-(5) of the main text) marked as dashed pink line.The same data is shown in Fig. 3(D) of the main text.The red and black solid lines correspond to the I d for the TW, when integrating over different temporal windows T .