---
abstract: 'We report on the influence of disorder on an exciton-polariton condensate in a ZnO based bulk planar microcavity and compare experimental results with a theoretical model for a non-equilibrium condensate. Experimentally, we detect intensity fluctuations within the far-field emission pattern even at high condensate densities which indicates a significant impact of disorder. We show that these effects rely on the driven dissipative nature of the condensate and argue that they can be accounted for by spatial phase inhomogeneities induced by disorder, which occur even for increasing condensate densities realized in the regime of high excitation power. Thus, non-equilibrium effects strongly suppress the stabilization of the condensate against disorder, contrarily to what is expected for equilibrium condensates in the high density limit. Numerical simulations based on our theoretical model reproduce the experimental data.'
author:
- Martin Thunert
- Alexander Janot
- Helena Franke
- Chris Sturm
- Tom Michalsky
- María Dolores Martín
- Luis Viña
- Bernd Rosenow
- Marius Grundmann
- 'Rüdiger Schmidt-Grund'
title: Cavity Polariton Condensate in a Disordered Environment
---

Introduction
============

The observation of a macroscopically coherent quantum state of exciton-polaritons, a so-called polariton Bose-Einstein condensate (BEC), [@Kasprzak.2006; @Balili.2007] has opened an active and challenging research field. Exciton-polaritons (for brevity polaritons) are mixed light-matter excitations in a microcavity (MC). [@Carusotto.2013; @Deng.2010] At finite quasi-particle density, several fascinating phenomena like superfluidity [@Amo.2009; @Amo.2009b; @Sanvitto.2010] and the formation of quantum vortices [@Lagoudakis.2008], were discovered. This allows for numerous novel applications like optical parametric oscillators [@Baumberg.2000], polariton lasers [@Schneider.2013; @Bhattacharya.2014] and logical elements [@Bajoni.2008; @Steger.2012; @Ballarini.2013; @Anton.2013b; @Sturm.2014], which are usually restricted to low temperatures. However, polariton BECs even at room-temperature were observed in MCs based on wide band gap materials like GaN [@Christopoulos.2007; @Christmann.2008; @Daskalakis.2013] and ZnO [@Lu.2012; @Li.2013; @Lai.2012] or organic materials [@Plumhof.2014], paving the way for technological applications. At the moment, experiments in these materials are significantly affected by disorder [@Christopoulos.2007; @Franke.2012; @Trichet.2013], and a thorough understanding of the impact disorder has on experimental observables in a polariton BEC is called for. In contrast to conventional BECs, occurring for example in cold atom systems, polaritons have a finite lifetime, which gives rise to unique properties of the condensate. Nonetheless, there remain similarities, for instance, in the absence of disorder quasi-long range order of a two-dimensional polariton condensate [@Roumpos.2012; @Chiocchetta.2013; @Spano.2012; @Spano.2013] and superfluidity is theoretically expected [@Wouters.2010; @Keeling.2011] and experimentally observed. [@Amo.2009; @Amo.2009b; @Sanvitto.2010] However, recent theoretical studies have revealed exciting differences between equilibrium and non-equilibrium condensates [@Sieberer.2013; @Sieberer.2014; @Tauber.2014; @Altman.2013; @Janot.2013]. For example, it is predicted that correlation functions for the condensate wave function decay exponentially [@Altman.2013] and that superfluidity vanishes in the presence of disorder. [@Janot.2013]

A polariton BEC is a steady state out of equilibrium where losses are compensated by external excitation. In the presence of disorder, spatial inhomogeneities of the condensate phase are induced. [@Janot.2013] If the phase fluctuates on length scales comparable to the condensate size, spatial correlations and phase rigidity are strongly reduced. In our work we will show that this leads to significant traces of disorder in the experimentally observed $k$-space intensity distribution, and theoretically demonstrate that the ratio of the condensate correlation length to the condensate size is independent of the condensate density. Consequently, in polariton condensates the stabilization against disorder fluctuations with increasing condensate density is strongly suppressed as compared to condensates in equilibrium. This prediction is supported by experimental investigations of the impact of disorder on a two-dimensional polariton BEC in a ZnO based MC. We measure the $k$-space intensity distribution as a function of excitation power, or rather condensate density, and observe significant disorder effects even at high densities. Numerical simulations allow to compare with experimental data confirming our theoretical predictions.

For an equilibrium BEC our observations would be unexpected, since an increasing density screens the disorder potential and leads to an ordered superfluid state [@Nattermann.2008; @Falco.2009; @Malpuech.2007]. Analogously, for a polariton BEC, interactions also can lead to superfluidity, as observed in clean samples [@Amo.2009; @Amo.2009b; @Sanvitto.2010]. However, as mentioned above, in the presence of disorder the polariton BEC is strictly speaking not a superfluid and long-range order is destroyed. [@Janot.2013] Thus, we expect and observe that disorder affects a dissipative polariton BEC much more than an equilibrium one. Several further observations found in literature seem to support this. For example, in one-dimensional CdTe MCs [@Manni.2011; @Stepnicki.2013] and ZnO MCs [@Trichet.2013] the spatial first-order correlation function of polariton BEC emission was analyzed in the presence of disorder and significant changes due to disorder were found. In the CdTe MCs the disorder effects remain present even with increasing excitation power, similarly to our findings in two-dimensional ZnO MCs. We note that the correlation length of the assumed disorder potential discussed in Ref.  is of the order of microns, which enables the trapping of the entire condensate. This is explicitly excluded in our model, since the disorder correlation length is assumed to be much smaller than the condensate size leading to spatial density and phase fluctuations of the condensate instead. Moreover, in various works on two-dimensional polariton BECs in CdTe based MCs disorder effects were also observed, leading to fluctuations within the far-field photoluminescence (PL) distribution [@Richard.2005] or the spatial first-order correlation function [@Kasprzak.2006]. Even frequency desynchronization between spatially separated condensate fragments can be induced, if the ratio between the disorder potential and the polariton interaction potential strength exceeds a critical value. [@Baas.2008; @Krizhanovskii.2009; @Wouters.2008b; @Eastham.2008] However, the dependence of the condensate density on the disorder effects was not analyzed within these works.

The paper is organized as follows: In Sec. \[sec:theory\_disorder\] we introduce our theoretical model. We discuss the disorder impact on a homogeneously and inhomogeneously excited condensate for a quasi-equilibrium (weak gain and loss) and driven dissipative (strong gain and loss) condensate, respectively. Furthermore, we provide a general argument that explains our experimental findings. These are presented in Sec. \[sec:experimental\_results\]. In Sec. \[sec:fit\_of\_exp\_data\] the theoretical predictions are confirmed by comparing experimental data to theoretical simulations. The summary and conclusion can be found in Sec. \[sec:conclusion\].

Theoretical Predictions {#sec:theory_disorder}
=======================

Model {#sec:theory_model}
-----

A phenomenological description of the dynamics of the polariton condensate wave function $\Psi(\vec{x},t)$ is given by an extended Gross Pitaevskii equation (eGPE) [@Wouters.2007; @Keeling.2008] $$\begin{aligned}
    \label{eq:eGPE}
    {\mathrm{i}}\hbar {\partial_{t}} \Psi = &\left( -\frac{\hbar^2}{2 m} \vec{\nabla}^2 + V(\vec{x}) + U \left|\Psi\right|^2 \right)\Psi \\
    \nonumber
        &+ {\mathrm{i}}\left(R(\vec{x}) - \Gamma \left|\Psi\right|^2\right) \Psi \ ,\end{aligned}$$ where $m$ is the effective mass of the lower polariton branch, $V$ an external potential and $U>0$ an onsite interaction constant. The function $R(\vec{x})$ describes the linear part of gain and loss due to inscattering from a reservoir of non-condensed polaritons and the finite lifetime of the condensate. The non-linearity $\Gamma |\Psi|^2$ implements a density dependent gain saturation with $\Gamma$ as gain depletion constant. Since the propagation of the reservoir polaritons can be neglected, the spatial shape of $R(\vec{x})$ can be related to the Gaussian profile of the excitation laser, namely $$\begin{aligned}
    \label{eq:reservoirfct}
    R(\vec{x}) = \hbar \gamma_{{\mathrm{c}}} \left(\frac{P}{P_{{\mathrm{th}}}} e^{-\vec{x}^2/\xi_P^2} - 1 \right) \ .\end{aligned}$$ The parameter $\gamma_{{\mathrm{c}}}$ is the condensate decay rate (inverse lifetime $\gamma_{{\mathrm{c}}} = 1/\tau$). The ratio $P/P_{{\mathrm{th}}}$ is the excitation power versus its value at threshold ${P_{{\mathrm{th}}}}$ at which condensation is observed first, and $\xi_P$ is the waist size of the Gaussian pump spot. We note that for the case of a spatially homogeneous excitation the eGPE  was successfully used to analyze a driven dissipative condensate. [@Sieberer.2013; @Altman.2013]

Because of interactions, the condensate energy is blueshifted by $n_0U$ where $n_0$ is the mean condensate density determined by the balance of gain and loss (for a definition of $n_0$ see [Eq. ]{}). The healing length $\xi \equiv \hbar / \sqrt{2 m n_0 U}$ is obtained by comparing kinetic and interaction energy of [Eq. ]{}.

The disordered environment is described by a random potential $V(\vec{x})$. We choose Gaussian-distributed delta-correlated disorder with zero mean and variance $\xi_V^2 V_0^2$, see Appendix \[sec:theory\_detailsmodel\] for details. We introduce an effective dimensionless disorder parameter, $$\begin{aligned}
  \label{eq:kappa}
  \kappa \equiv \frac{\xi_V \: V_0}{\xi\: n_0U} \ .\end{aligned}$$ An analysis of the gain and loss terms in [Eq. ]{} allows us to define a ’non-equilibrium parameter’ $$\begin{aligned}
  \label{eq:alpha}
  \alpha \equiv \frac{\Gamma}{U} \ .\end{aligned}$$ Its magnitude parametrizes the influence of gain and loss on the polariton BEC. For example, in the limit $\alpha \to 0$ (keeping $n_0$ finite) the equilibrium mean field description of a BEC is obtained, and, on the other hand, in the limit $\alpha \to \infty$ the condensate is totally dominated by gain and loss. In this work, we will focus on single-mode steady-state solutions and therefor make the ansatz $\Psi(\vec{x},t) = \Psi(\vec{x}) \exp( - {\mathrm{i}}\omega t)$, where $\hbar \omega$ is the condensate energy. However, in experimental realizations more than one condensate mode can exist. For any further details we refer to Appendix \[sec:theory\_detailsmodel\].

Disorder Effects
----------------

### Infinite condensate size

Before we discuss a finite size polariton BEC we would like to consider a homogeneously excited condensate ($\xi_P \to \infty$), such that the reservoir function [Eq. ]{} is a constant in space. We will i) review disorder effects on an equilibrium condensate [@Nattermann.2008; @Falco.2009], and ii) describe differences to a polariton BEC (driven dissipative condensate) [@Janot.2013].

*Equilibrium condensate i):* The disorder potential attempts to pin the condensate into its minima, whereby the energy costs for density deformations (kinetic term in [Eq. ]{}) have to be compensated. The balance of pinning and kinetic energy determines the density Larkin length ${\mathcal{L}_{\mathrm{n}}}\approx \sqrt{\pi} \: \hbar^2 / m \: \xi_V V_0$ [@ImryMa.1975; @Nattermann.2008; @Falco.2009]. On the other hand, for a sufficiently large interaction energy $n_0U$ the disorder gets screened. [@Falco.2009] The ratio of healing to Larkin length, $\xi/{\mathcal{L}_{\mathrm{n}}}\sim \kappa$, describes this competition of disorder and interaction. For $\xi \ll {\mathcal{L}_{\mathrm{n}}}$ ($\xi \gg {\mathcal{L}_{\mathrm{n}}}$) the interaction energy is large (small) as compared to the disorder potential. Due to the fact that the interaction energy increases with increasing density (and $\xi \propto 1/\sqrt{n_0}$), $\xi/{\mathcal{L}_{\mathrm{n}}}$ decreases with increasing density, and disorder effects will fade away in this limit. Thus, for sufficiently high densities an equilibrium condensate will be ordered and superfluid. [@Falco.2009] *Non-equilibrium condensate ii):* In a driven system the mean density $n_0$ of the condensate is determined by a balance of gain and loss. Disorder induces density fluctuations about this mean value. In a region with reduced density, as compared to $n_0$, the gain mechanism tries to compensate the depletion, and more particles are scattered into the condensate than decay. On the other hand, in a region with increased density more particles decay than are injected from the reservoir. By virtue of the continuity equation, these local particle sources and sinks are connected by condensate currents. Because the density fluctuates randomly in space, a random distribution of sources and sinks forms and, thus, a random pattern of current flow is generated. The condensate current is proportional to the product of the density and the gradient of the condensate phase. Since the current is not constant, the phase cannot vary uniformly in space, and thus a random current configuration gives rise to a spatially fluctuating phase. We note that in this work the term ’fluctuations’ will be used for random spatial inhomogeneities. The correlation length, over which the phase typically varies by $2\pi$, is given by ${\mathcal{L}_{\phi}}\approx \sqrt{2}\pi \, {\mathcal{L}_{\mathrm{n}}}/ \alpha$. This scale can be obtained by a generalized Imry-Ma argument [@Janot.2013]: a condensate current flowing out of (or into) a region of diameter ${\mathcal{L}_{\phi}}$ is generated by an effective source (or sink) determined through an area average of multiple random sources and sinks. In contrast to an equilibrium condensate ($\alpha \to 0$ with ${\mathcal{L}_{\phi}}\to \infty$), the phase fluctuations occurring in the case ${\mathcal{L}_{\phi}}< \infty$ destroy the quasi-long-range order of the condensate. As a consequence of these phase fluctuations, the superfluid stiffness vanishes in the thermodynamic limit even for weak disorder, and a superfluid behavior is only present below a finite length scale, namely the superfluid depletion length ${\mathcal{L}_{s}}\approx \sqrt{2} \pi \, {\mathcal{L}_{\mathrm{n}}}/ \alpha^2$. [@Janot.2013]

### Finite condensate size

From the analysis above we conclude that in a disordered environment a condensate of size ${L_{c}}\ll {\mathcal{L}_{\phi}}$ will behave completely different from one of size ${L_{c}}\sim {\mathcal{L}_{\phi}}$.

In the following, we discuss these two scenarios sketched schematically in Fig. \[fig:sketch\_scenarios\]B.

For scenario I with ${L_{c}}\ll {\mathcal{L}_{\phi}}$ (called quasi-equilibrium in the following) the phase is correlated over the entire condensate region, and disorder induces mainly density fluctuations. As discussed above, the impact of disorder will decrease with increasing density, which should be directly observable by increasing the excitation power. Such kind of percolation transition from a disordered to an ordered regime was predicted (for a polariton BEC in equilibrium) in Ref. .

In the presence of gain and loss disorder induces phase fluctuations as explained above. For scenario II we assume that the phase correlation length ${\mathcal{L}_{\phi}}$ is comparable to the condensate size ${L_{c}}$, i.e. ${\mathcal{L}_{\phi}}\sim {L_{c}}$, such that spatial correlations and superfluidity are destroyed. The ratio ${L_{c}}/{\mathcal{L}_{\phi}}\propto (V_0 \ \xi_{{\mathrm{P}}} \xi_{{\mathrm{V}}} m / \hbar^2) \: (\Gamma / U)$ does not depend on the condensate density and, thus, is independent of the excitation power. A similar conclusion holds for the ratio ${L_{c}}/ {\mathcal{L}_{s}}$. As a consequence, a condensate stabilization with increasing density, as present in an equilibrium system, is strongly suppressed.

![(color online) Expectation value $\mu_P$ (upper row) and variance $\sigma^2_P$ (lower row) of the (normalized) intensity distribution $I_P(k)$. In order to compare the fluctuations for increasing excitation powers we present the ratio $\sigma^2_P(k) / \sigma^2_{P_0}(k)$ with $P_0 = 4 P_{{\mathrm{th}}}$. The left and right column depict a quasi-equilibrium (scenario I) and a driven dissipative (scenario II) condensate, respectively. For wavevectors $|k \xi| \lesssim 1$, scenario I shows a linear reduction of fluctuations with inverse excitation power, $\sigma^2_P(k) \propto 1/P$, while scenario II exhibits a suppressed stabilization with increasing excitation power. We averaged $1560$ disorder realizations and used: ${\mathcal{L}_{\mathrm{n}}}/{L_{c}}= 1 \: (4.5)$, ${\mathcal{L}_{\phi}}/{L_{c}}= 10 \: (2.5)$, ${\mathcal{L}_{s}}/{L_{c}}= 20 \: (0.3)$ for scenario I (II). []{data-label="fig:IkMeanVar"}](./{fig_MeanVar_Ik_cP4_ckappa0.4_0.1_alpha0.5_8_beta04}.pdf){width="1\hsize"}

In order to make our analysis more quantitative, we have studied theoretically the excitation power dependence of the two-dimensional $k$-space intensity $I_P({\vec{k}}) \propto \gamma_c |\Psi_{\vec{k}}|^2$ which can be directly compared to experimental data. To this end, [Eq. ]{} was simulated for many disorder realizations (see Appendix \[sec:theory\_numsim\] for details). We have extracted the expectation value, denoted by $\mu_P(k)$, and the variance, denoted by $\sigma^2_P(k)$, of the normalized intensity $I_P(\vec{k})$ by averaging over disorder configurations. We note that for a sufficiently large number of realizations, the disorder average restores radial symmetry, such that the expectation values $\mu_P(k)$ and $\sigma^2_P(k)$ depend on the magnitude $k = |\vec{k}|$ of wavevector only.

In Fig. \[fig:IkMeanVar\], the results for $\mu_P$ and $\sigma_P^2$ are shown for scenario I (left panels) and II (right panels). We find that the intensity $I_P$ vanishes for all wavevectors outside of the lower-polariton dispersion ( $k \gtrsim \xi^{-1}$) and that its average value does not change qualitatively as compared to a disorder-free system (cf. Ref. ). However, for a single snap-shot (see Fig. \[fig:snapshot\]) disorder breaks the radial symmetry and induces intensity fluctuations proportional to $\sigma_P$. For scenario I and for wavevectors $|k| \lesssim \xi^{-1}$, these fluctuations decay linearly with inverse excitation power, in agreement with the expectation $\sigma^2_P \sim \kappa^2 \propto 1/P$ for $\kappa \ll 1$. We note that regions with $k \approx \xi^{-1}$ show a high ratio $\sigma_P / \sigma_{\! P_{0}}$ (peaks in Fig. \[fig:IkMeanVar\] lower left panel). In this $k$-region, the emission intensity is increasing very rapidly with excitation power (see Fig. \[fig:IkMeanVar\] upper left panel), because of the repulsive potential hill created by the finite excitation spot. [@Wouters.2008] Thus, the increase of fluctuation strengths with excitation power for $k \approx \xi^{-1}$ is really due to the increase of emission power and does not yield information about the screening of the disorder potential for high condensate densities.

For scenario II, the stabilization with increasing excitation power is suppressed (see lower right panel of Fig. \[fig:IkMeanVar\]). As compared to scenario I, the decrease of $\sigma^2_P$ with increasing condensate density is weaker than $\sigma^2_P \propto 1/P$. These findings agree well with our argument provided above.

The reservoir of non-condensed polaritons interacts with the condensate and thus leads to an increase of the blueshift. [@Wouters.2007; @Wertz.2010; @Anton.2013c] Usually, this is accounted for by adding a potential term proportional to the reservoir density in [Eq. ]{}. [@Wouters.2007] Such a term will modify the emission frequency of the condensate (real part of [Eq. ]{}), however, does not change the non-equilibrium continuity equation (imaginary part of [Eq. ]{}). Hence, the mechanism of generating random condensate currents is not altered qualitatively by reservoir-condensate interaction and, thus, we believe that they can be safely neglected for our analysis.

Experiment {#sec:experimental_results}
==========

In this section we discuss the experimentally observed behavior of the far-field PL emission pattern of a polariton condensate in a ZnO-based MC with pronounced structural disorder as a function of excitation power. For this experiment, the sample was excited using a pulsed Nd:YAG laser with a pulse duration of $500~{\mathrm{ps}}$. This is three orders of magnitude larger than the polariton relaxation time (0.4 ${\mathrm{ps}}$) which is determined from the spectral linewidth of the condensate emission. Thus, we can assume a quasi–continuous-wave excitation, which justifies the comparison with numerical simulations based on a steady state theory as will be discussed in Sec. \[sec:fit\_of\_exp\_data\]. Further details about the experimental setup can be found in Appendix \[sec:sample\_setup\]. The MC consists of a half wavelength ZnO cavity, which simultaneously acts as active medium, showing a quality factor of about 1000 and a maximum coupling strength of about 45 meV ($\Omega_{{\mathrm{Rabi}}} \approx 90~{{\mathrm{m}} {\mathrm{e}} {\mathrm{V}}}$) at $T = 10~{\mathrm{K}}$. By using a wedge-shaped cavity, the detuning between the cavity mode energy and the excitonic transition energy strongly varies with the lateral sample position. Structural investigations (atomic force microscopy, X-ray diffraction, cross-sectional transmission electron microscopy) yield a smooth but polycrystalline cavity layer, exhibiting a low interface roughness of $R_{{\mathrm{rms}}}$ = 1.9 nm. Furthermore, the cavity layer is preferentially $c$-plane oriented and laterally textured, containing large grains aligned in the growth direction reaching from the bottom to the top (grain sizes ranging from 20 nm up to 120 nm). Further information about the sample properties can be found in Ref. . Due to the textured structure we suppose that an electronic disorder potential is primarily caused by depletion of carriers, e.g. aluminum donor bound excitons [^1], due to interface band bending at grain boundaries [@Orton.1980]. (see Appendix \[sec:experiment\_disorderorigin\] for details).

![(color online) (a)-(d): Excitation power series of the far-field PL emission in a linear false color scale for $T = 10~{\mathrm{K}}$ and a detuning of $\Delta = -30~{{\mathrm{m}} {\mathrm{e}} {\mathrm{V}}}$. The excitation power is normalized to the condensation threshold $P_{{\mathrm{th}}}$. (e) $I_P(k)$ profiles of the condensate are extracted. These are integrated over the energy range $\Delta E$ marked by the white lines in the far-field emission pattern. The PL intensity for each power is normalized to the mean value of each $I_P(k)$ profile.[]{data-label="fig:Disorder"}](./{fig_Disorder_Images_10K_30meV_7.pdf}){width="0.9\hsize"}

Figures \[fig:Disorder\](a)-(d) show the excitation power dependence of the PL $k$-space emission pattern for $T = 10 \, {\mathrm{K}}$ and detuning $\Delta = -30 \, {{\mathrm{m}} {\mathrm{e}} {\mathrm{V}}}$. We deduce a polariton effective mass of $m=4.4 \times 10^{-5}~m_{{\mathrm{e}}}$ ($m_{{\mathrm{e}}}$: free electron mass) from the dispersion of the lower polariton branch (LPB) (not shown here). The excitation power density at condensation threshold is $P_{{\mathrm{th}}} = 79~{\mathrm{W}}{\mathrm{cm}}^{-2}$. Note that the determination of the excitation power density at threshold is quite complex, e.g. due to the coexistence of intense emission from uncondensed polaritons for $P \gtrsim {P_{{\mathrm{th}}}}$, but significant for the comparison with theoretical calculations discussed in Sec. \[sec:fit\_of\_exp\_data\]. Details for the experimental determination of ${P_{{\mathrm{th}}}}$ can be found in the Supplemental Material, Sec. \[sec:SM\_Pth\].

In all cases investigated here, the condensate emission is distributed dispersion-less at horizontal lines in $k$-space with maximum intensity between the LPB dispersion, which is visible in the far-field PL images (cf. Fig. \[fig:Disorder\]) for low excitation power $P \gtrsim {P_{{\mathrm{th}}}}$. This indicates a weak expansion of the condensed polaritons due to the background potential induced by the excitation spot, whose size is similar or even larger than the polariton propagation length. [@Franke.2012; @Wouters.2008] For the lowest excitation power shown here, $P = P_{{\mathrm{th}}}$, the emission intensity from the uncondensed polaritons and the condensate are of same order which prevents a clear distinction. With increasing excitation power the BEC states undergo a blueshift due to the increasing interaction potential, and we observe several states with different energy. Previous studies in the literature on this multimode behavior show that the emission from coexisting individual modes originates from different regions of the same condensate. [@Krizhanovskii.2006; @Baas.2008; @Krizhanovskii.2009] However, other studies on polariton condensates in a disordered environment found that long-range spatial coherence is still present for their energy-averaged emission [@Kasprzak.2006; @Richard.2005] indicating persistent correlations between different, possibly spatially separated condensate states.

For a wide range of excitation powers, condensate emission out of two energy ranges is observed, which are stable and energetically well separated. For a further analysis we select only one of these energy channels, in order to compare with numeric simulations of a single-mode condensate, cf. Sec. \[sec:fit\_of\_exp\_data\]. In Fig. \[fig:Disorder\](b)-(d) we marked the selected energy channel by two white dashed lines. This delimitation is defined by an energy range $\Delta E$ which corresponds to the excitation power dependent full width at half maximum of the condensate emission. Fig. \[fig:Disorder\](e) shows the far field emission profiles $I_P(k)$ for the selected energy channel and increasing excitation power, integrated over $\Delta E$. The $I_P(k)$ profiles show several randomly distributed inhomogeneities and differ strongly from the smooth and ideally radial symmetric distribution expected for a disorder-free sample. [@Wouters.2008] Remarkably, the intensity fluctuations persist even for high excitation power, i.e. high condensate densities. We note that the constant sharp stripes in the $I_P(k)$ profile at a specified *k* for all excitation powers are caused by imperfections of the setup, probably due to the microscope objective.

A similar finding with increasing excitation power was also observed for other detunings $\Delta = -50~{\mathrm{meV}},\ldots,-10~{\mathrm{meV}}$, and we conclude that our observation does not depend significantly on the particular choice of detuning within the mentioned range.

To investigate the temporal coherence properties of the condensate we used a Michelson interferometer in the plane mirror (PM) - retroreflector (RR) configuration to superimpose the PL emission of polaritons with opposite emission angles or rather wavevectors. For this experiment, the sample was excited by a frequency-tripled Ti:sapphire laser at 266 nm with a pulse duration of about 2 ps. Further details of the setup are provided in Appendix \[sec:sample\_setup\].

![(color online) (a)-(b): Energy-resolved far-field PL interference pattern in a linear false color scale for $T = 10~K$ and $\Delta = -47$ meV: (a) large path difference $\Delta s$: uncorrelated far-field PL emission (basically the sum of the emission from the individual interferometer arms), (b) path difference close to zero: distinct interference fringes indicating mutual temporal coherence of the polariton emission. The red dashed line indicates the energy of the investigated condensate state. (c) Normalized visibility of the interference fringes as a function of the path difference.[]{data-label="fig:EkSpace_Interference"}](./fig_Coherence_time_neu_v04a.pdf){width="\linewidth"}

The RR is mounted on a motorized linear stage, that allows us to vary the path difference $\Delta s$ between the emission collected from both interferometer arms. Fig. \[fig:EkSpace\_Interference\](a)-(b) show two selected interferograms of the $I(E,k)$ emission pattern for large (Fig. \[fig:EkSpace\_Interference\](a)) and short (Fig. \[fig:EkSpace\_Interference\](b)) path differences $\Delta s$, respectively. To investigate the temporal coherence properties of the polariton condensate we analyzed the normalized visibility of the interference fringes, $$\begin{aligned}
    \label{eq:Inorm}
    V_{{\mathrm{norm}}} & = \frac{I_ {\text{interf}} - I_{RR} - I_{PM}}{2\sqrt{I_{RR}  I_{PM}}} = g^1(\Delta t) \cos(\phi_{12}),\end{aligned}$$ as a function of the temporal delay $\Delta t = \frac{\Delta s}{c}$, where $c$ is the speed of light (cf. Fig. \[fig:EkSpace\_Interference\](c)). Here, $I_{{\mathrm{interf}}}$ is the intensity of the interference pattern, and $I_{{\mathrm{RR}}}$, $I_{{\mathrm{PM}}}$ are the intensities of the RR and PM arm, respectively, $g^1$ is the first-order coherence function and $\phi_{12}$ is the phase difference between the emission from the individual interferometer arms. By assuming a Gaussian decay of $g^1(\Delta t) = g^1(0) \exp[
- (\pi / 2) (\Delta t ^2 / \tau_{{\mathrm{coh}}}^2)]$[@Saleh.2008] we determined a coherence time of about $\tau_{{\mathrm{coh}} } = 8.7 \, {\mathrm{ps}}$. This is more than 50 times larger than the lifetime of the uncondensed polaritons of about 160 fs, which is deduced from the spectral linewidth of the polariton emission for $P \leq P_{{\mathrm{th}}}$ and for an energy range similar to the condensate energy at $P = P_{{\mathrm{th}}}$. Consequently, the coherence of the investigated quantum system is conserved during the multiple reabsorption and reemission processes, which can thus be identified as a condensate. We note that the experimentally estimated coherence time is a lower limit for the real value. We identify two experimental artifacts that restrict the determination of the real condensate’s coherence time, namely a spectrally and path difference dependent phase shift $\phi_{12}(\lambda, \Delta s)$ (artifact A) as well as a fast decay of the condensate emission intensity due to the short excitation pulses of about 2 ps that are used for the coherence time measurement (artifact B). By analyzing the impact of these artifacts quantitatively (cf. Supplemental Material, Sec. \[sec:SM\_limits\_cohTime\]), we roughly estimated the expected real values for the coherence time of $\tau^A_{\mathrm{corr}} = 10.3$ ps and $\tau^B_{\mathrm{corr}} = 14$ ps. By applying both corrections simultaneously, a maximum coherence time of $\tau_{\mathrm{corr}}$ = 24 ps was estimated.

For an ideal (homogeneous, disorder-free) condensate a linewidth of $\Delta E \approx 0.66 \,  h / \tau_{\mathrm{coh}} = 0.66 \, h / 8.7~{\mathrm{ps}} = 0.31~{\mathrm{meV}}$ would be expected for the condensate emission according to the Wiener-Khinchin theorem [@Saleh.2008] (and even less assuming the corrected values for $\tau_{\mathrm{coh}}$), where $h$ is the Planck constant. This is about a factor of 6.5 smaller than the observed minimum linewidth of 2 meV for the condensate emission in this experiment. Since the investigated condensate is a complex quantum system including spatial density and phase fluctuations we assume that the Wiener-Khinchin theorem cannot be applied here. We rather suppose that the mechanism which causes a broadening of the emission linewidth (e.g. repulsive particle interaction [@Porras.2003]) does not affect the coherence time to the same extent. This is supported by the quantitative discrepancy between the emission linewidth and the coherence time, which is observed also in a CdTe [@Kasprzak.2006] as well as in a ZnO MC. [@Lai.2012] We note that despite of the fast decay of polaritons, condensate emission can be observed up to 90 ps after the arrival of the exciting laser pulse, which thus allows for the experimental observation of coherence in the mentioned time range.

Summarizing, the experimental observations indicate a strong impact of disorder on the polariton BEC even at high excitation power well above the condensation threshold. As discussed in Sec. \[sec:theory\_disorder\] the suppression of disorder effects with increasing condensate density is strongly hindered for a polariton BEC. We assume that the interplay of gain-loss and disorder prevents a stabilization at high excitation power also in the experiment.

Comparison between Theoretical Model and Experiment {#sec:fit_of_exp_data}
===================================================

In the following, we will compare our experimental observations with numerical simulations. At threshold $P={P_{{\mathrm{th}}}}$ a cross-over from a non-condensed state to a polariton BEC takes place, typically indicated by a super-linear increase of the emission intensity. Such a transition is not very well described by the used eGPE . For this reason, the data analysis is done well above threshold, where both experimentally observed and theoretically calculated blueshift (condensate density) increase linearly with pump power. We note that the evolution of the experimentally measured polariton blueshift $\Delta E$ as a function of the excitation power shows two kinks at $P = 2~P_{{\mathrm{th}}}$ and $P = 4~P_{{\mathrm{th}}}$ (cf. Fig. \[fig:Blueshift\] in Appendix \[sec:experiment\_disorderorigin\]). We believe that the slope of $\Delta E$ for $P < 2~P_{{\mathrm{th}}}$ is predominantly caused by an electronic disorder potential, which starts to saturate for $P = 2~P_{{\mathrm{th}}}$, and that for $P \geq 4~P_{{\mathrm{th}}}$ the blueshift is governed by condensate-condensate interactions. Further discussions are presented in Appendix \[sec:experiment\_disorderorigin\] and references therein.

For the comparison between the theoretical model and the experimental data, the parameters of the eGPE  are chosen according to the experiment, see Tab. \[tab:parameters\]. [^2] We note that a quantitative determination of the disorder parameter from experiment is very challenging, cf. discussion in Appendix \[sec:experiment\_disorderorigin\], and we chose $\xi_{{\mathrm{V}}} V_0 \approx 0.15 \: \mu{\mathrm{m}} \, {{\mathrm{m}} {\mathrm{e}} {\mathrm{V}}}$ for simulations.

                          $m$                                  $\tau$            $\xi_{{\mathrm{P}}}$                $P_{{\mathrm{th}}}$                $P/P_{{\mathrm{th}}}$   $d \Delta \! E / d(P \! / \! {P_{{\mathrm{th}}}})$        
  --------------------------------------------------- ------------------------ ------------------------ ---------------------------------------------- ----------------------- ---------------------------------------------------- -- -- --
   $4.4 \cdot 10^{-5}\, {\mathrm{m}}_{{\mathrm{e}}}$   $0.4 \, {\mathrm{ps}}$   $2 \, \mu{\mathrm{m}}$   $79 \, {\mathrm{W}} \! / \! {\mathrm{cm}}^2$       $2,\ldots,20$                    $0.7 \: {\mathrm{meV}}$                      

  : Parameters extracted from experiment and corresponding parameters used for simulations as well as relevant length scales. For definitions see Appendix \[sec:theory\_detailsmodel\].[]{data-label="tab:parameters"}

   $\alpha$   $P/P_{{\mathrm{th}}}$   ${\xi_{{\mathrm{V}}} V_0}/{l_c \, \hbar \gamma_c}$   ${\xi_{{\mathrm{P}}}}/{l_{{\mathrm{c}}}}$   ${\mathcal{L}_{\mathrm{n}}}/ {L_{c}}$   ${\mathcal{L}_{\phi}}/ {L_{c}}$   ${\mathcal{L}_{s}}/ {L_{c}}$
  ---------- ----------------------- ---------------------------------------------------- ------------------------------------------- --------------------------------------- --------------------------------- ------------------------------
     $7$          $5,11,15,20$                             $0.125$                                            $3$                                       $5$                                  $3$                            $0.4$

  : Parameters extracted from experiment and corresponding parameters used for simulations as well as relevant length scales. For definitions see Appendix \[sec:theory\_detailsmodel\].[]{data-label="tab:parameters"}

![(color online) Snap shots of the normalized two-dimensional intensity distribution $I_P(\vec{k})$ of a polariton BEC with ${\mathcal{L}_{\phi}}\sim {L_{c}}$ in a disordered environment for increasing excitation power $P$. The used parameters are presented in Tab. \[tab:parameters\]. []{data-label="fig:snapshot"}](./{fig_Snap_Ik2d_L120_cP3_ckappa0.125_alpha7_snap8}.pdf){width="0.95\hsize"}

For a typical disorder realization, a series of numerically obtained snap shots of the two-dimensional intensity distribution $I_P(\vec{k})$ for increasing excitation power is shown in Fig. \[fig:snapshot\]. These images correspond to a polariton BEC described by scenario II. We clearly observe a disorder-induced deviation from the ideally radial distribution, which does not converge to a symmetric intensity distribution while increasing the excitation power. Such an asymmetry as well as its persistence is also observed experimentally, see Fig. \[fig:Disorder\], and thus agrees qualitatively with our simulations. We note that the experimental data represent the intensity distribution of one disordered sample, and correspond to a one-dimensional cut along a given line crossing the origin of the two-dimensional $k$-space distribution, for example the $x$-axis.

For a quantitative analysis we compare directly the experimental measurements with the numerically computed expectation value $\mu_{{\mathrm{P}}}$ and variance $\sigma^2_{{\mathrm{P}}}$ of the intensity distribution. To this end we symmetrize the experimental data $I_P(k) \to (I_P(k) + I_P(-k))/2$ with $k \geq 0$, and superimpose them with the results of the numerical simulations. Since the condensate density and healing length are hard to determine experimentally, we fix the scaling of $x$- and $y$-axis by a least-square fit. Fig. \[fig:Fit\] shows the result. We have excluded experimental data with wavevectors $k \geq 3 \: \mu {\mathrm{m}}^{-1}$, because a systematic artifact is present for all $k = 3,\dots,4 \ \mu {\mathrm{m}}^{-1}$ and for all excitation powers. [^3]

![(color online) Experimental $I_P(k)$ distribution (blue diamonds) compared to numerical simulations of scenario II for increasing excitation power $P$. The black solid lines depict the expectation value whereas the gray band indicates the standard deviation. For high momenta we have excluded systematically biased data (gray crosses). For each plot the mean squared error (MSE) and the goodness-of-fit value $\mathcal{Q}$ (see text) are computed. For used parameters see Tab. \[tab:parameters\]. []{data-label="fig:Fit"}](./{fig_figfit_T10_detuning-30_cP3_ckappa0.125_alpha7}.pdf){width="1.00\hsize"}

In order to quantify the agreement between theory and experiment we introduce the mean squared error (MSE) and the goodness-of-fit value $\mathcal{Q}$ (see Appendix \[sec:theory\_numsim\] for definitions). $\mathcal{Q}$ is a probability: if $\mathcal{Q} \sim 1$, the simulations describe the experimental data. On the other hand, if $\mathcal{Q} \ll 1$ the theoretical model does not reproduce the experiment. The experimental data are well described by simulations (of scenario II), cf. Fig. \[fig:Fit\]: the MSE is close to zero, and the goodness-of-fit value $\mathcal{Q}$ remains comparable to one for all pump powers studied. In contrast, trying to reproduce the experimental observations by simulations of a polariton BEC described by scenario I (quasi-equilibrium) fails, cf. Fig. \[fig:FitScenarioI\] of the Appendix \[sec:theory\_numsim\]. Thereto, we had chosen a non-equilibrium parameter $\alpha = 0.5$ and slightly increased the disorder strength. Then, the $\mathcal{Q}$-value drops from $\mathcal{Q} \approx 10^{-8}$ at $P\approx 5~P_{{\mathrm{th}}}$ to $\mathcal{Q} \approx 10^{-21}$ at $P = 20~P_{{\mathrm{th}}}$.

Summary and Conclusion {#sec:conclusion}
======================

In this work we have characterized a polariton condensate in a disordered environment. Our theoretical analysis shows that spatial fluctuations of the condensate phase, which are induced by the interplay of disorder and gain-loss of particles, do not depend on the mean condensate density. This leads to a reduced stabilization against disorder fluctuations with increasing density in contrast to an equilibrium condensate. To verify our prediction we have analyzed experimentally the photoluminescence emission of a ZnO based microcavity. Indeed, we find a lack of stabilization with increasing density in terms of pronounced intensity fluctuations within the $k$-space emission pattern even at high excitation power. This experimental finding can be reproduced by numerical simulations. From this we conclude that the polariton condensate in the microcavity is exposed to significant structural disorder, and that the persistence of disorder effects even at high excitation power, well above the condensation threshold, relies on the intrinsic non-equilibrium nature of polaritons. We note that these findings may also explain the observation of similar phenomena for polariton condensates in microcavities based on other materials, e.g. CdTe or GaN. [@Christopoulos.2007; @Baas.2008; @Krizhanovskii.2009]

Acknowledgement
===============

We acknowledge experimental support and fruitful discussions from the group of N. Grandjean at EPFL. AJ is supported by the Leipzig School of Natural Sciences BuildMoNa. This work was supported by Deutsche Forschungsgemeinschaft through project GR 1011/20-2, by Deutscher Akademischer Austauschdienst within the project PPP Spain (ID 57050448) and by Spanish, MINECO Project Nos. MAT2011-22997 and MAT2014-53119-C2-1-R.\
MT and AJ contributed equally to this work. MT performed the experimental part and AJ carried out the theoretical analysis.

Experimental setup {#sec:sample_setup}
==================

In order to investigate the optical properties of the polariton condensate, we applied two different photo-luminescence configurations, which have in common a non-resonant and pulsed excitation as well as a detection of the far-field emission. The setup to investigate the disorder effects on the polariton distribution and their dynamics as a function of the excitation power is described in Ref. . Here, the excitation was carried out by a pulsed Nd:YAG laser with pulse duration of 500 ps, whose Gaussian excitation spot covers a sample area of about 10 $\mu$m$^2$.

For the coherence measurements, the sample was excited via a frequency-tripled Ti:sapphire laser at 266 nm (repetition rate: 76 MHz, pulse length $\approx$ 2 ps). The PL signal of the Fourier plane was sent to a Michelson interferometer in the mirror-retroreflector configuration. The retroreflector image is a centrosymmetric counterpart of the mirror arm image. In the resulting interferogram we superimposed the signal with wavevector $\vec{k_{||}}$ with that of $-\vec{k_{||}}$. Interference maxima occur when the path difference between the individual beams, $\Delta L$ = $c \Delta t$, is an integer multiple of the PL emission wavelength, being $\Delta t$ the delay between the beams and $c$ the speed of light. With the help of a streak camera the relative delay between the two arms was set to zero for $k_{||} = 0$.

Real space measurements with a sufficient spatial resolution could not be performed due a to a spherical aberration induced by the cryostat window. For the conditions used in our experiments, namely the UV spectral range, a window thickness of 1.5 mm and the large range of collected emission angles of $\pm 23^\circ$, the resulting spatial distortion of the image is larger than structural fluctuations that we would like to resolve. Consequently, the distortion of the measured real space image prevents a precise investigation of the spatial distribution of the luminescence as well as spatially resolved correlation measurements. We note that far-field images are not affected by the cryostat window, which causes a parallel beam shift but does not change the angle of the transmitted rays.

Origin of disorder potential {#sec:experiment_disorderorigin}
============================

![(color online) Polariton blueshift $\Delta E$ as a function of the excitation power normalized to the condensation threshold for the observed ZnO MC. The red line corresponds to the polariton-polariton interaction whereas the blue line represents the blueshift due to an additional electronic background potential, which starts to saturated at about 2 $P_{{\mathrm{th}}}$ and is presumably totally saturated for $P > 4 {P_{{\mathrm{th}}}}$ .[]{data-label="fig:Blueshift"}](./fig_Blueshift_vs_excPower4.pdf){width="30.00000%"}

Due to the dual light-matter nature of the polaritons, the effective disorder potential can be of photonic as well as electronic origin. Photonic disorder can be caused by surface and interface roughness as well as thickness fluctuations within the MC structure. This leads to a spatial fluctuating cavity length and therefore to a variation of the cavity photon energy. Due to results of other ZnO-based MCs, a minimum potential strength of $V_C \geq$ 2 meV can be expected. [@Li.2013; @Orosz.2012] The corresponding correlation length $\xi_V$ is of the order of the photonic wavelength, of about 370 nm.

In the literature, usually electronic disorder is neglected. [@Trichet.2013; @Baas.2008; @Savona.2007; @Manni.2011; @Nardin.2009] In contrast to this, we assume a strong influence of an electronic background potential caused by randomly distributed excitonic states which are accumulated within the bulk of grains [@Orton.1980] or bound to impurities. This is supported by two facts: firstly, cross-sectional TEM analysis of a MC that is fabricated under the same conditions, provides a granular structure of the investigated ZnO MC with grain sizes ranging from 20 nm up to 120 nm. [@Franke.2012] Secondly, the slope of the polariton blueshift $\Delta E(P)$ is by a factor of about 6.3 larger for $P$ &lt; 2 $P_{\textnormal{th}}$ than above and even by a factor of 12.6 larger compared to the blueshift for $P > 4~{P_{{\mathrm{th}}}}$ (cf. Fig. \[fig:Blueshift\]). This can be explained by assuming an additional electronic background potential $\Delta E_b$, which may include localized states within a disorder potential or bound to impurities, as shown in Ref. . Since the concentration of these electronic defects is finite, their contribution to the condensate blueshift saturates for a certain excitation power or rather condensate density. Thus, the further blueshift for $P$ &gt; 4 $P_{\textnormal{th}}$ is restricted to condensate-condensate interaction.

We assume that the condensate blueshift for small excitation power $\Delta E$($P$ &lt; 2 $P_{\textnormal{th}})$ is primarily caused by its interaction with aluminum donor bound excitons ($D^0,X$) and that $\Delta E$ scales linearly with its concentration. As mentioned in Sec. \[sec:experimental\_results\], we suppose a depletion of bound excitons at grain boundaries and thus an accumulation of them within the grain bulk. [@Orton.1980] According to the model described in Ref.  the grain boundaries act like two back-to-back Schottky barriers and the carrier flow between grains is driven by thermionic emission over the Schottky barrier. In general, the average height and width of these barriers can be determined from the temperature-dependent evolution of the hall mobility. Unfortunately, this was not possible for our MC due to low current values, below the resolution limit of 1 nA, for temperatures below 200 K, caused by the small cavity thickness of about 100 nm as well as due to strong inhomogeneities of the current density, which may be caused by the cavity thickness gradient.

Assuming the mechanism of carrier depletion at grain boundaries to be the dominant one for the effective electronic disorder potential, its correlation length $\xi_V$ is similar to the grain size with values between 20 nm and 120 nm. This is about two orders of magnitude below the condensate size $L_c$, limited by the size of the pump spot and thus even lower than the assumed correlation length for photonic disorder of about 370 nm. Consequently, a trapping of the entire condensate within a minimum of the disorder potential can be excluded. We rather suppose that the disorder potential causes condensate density fluctuations and thus phase fluctuations due to the interplay of disorder and the non-equilibrium nature of the polariton condensate.

Details of the Model {#sec:theory_detailsmodel}
====================

A phenomenological description of the dynamics of the macroscopic polariton condensate wave function $\Psi(\vec{x},t)$ is given by an extended Gross-Pitaevskii Equation (eGPE), [@Wouters.2007; @Keeling.2008] $$\begin{aligned}
    \nonumber
    {\mathrm{i}}\hbar {\partial_{t}} \Psi = &\left( -\frac{\hbar^2}{2 m} \vec{\nabla}^2 + V(\vec{x}) + U \left|\Psi\right|^2 \right)\Psi \\
    \label{eq:eGPEdim}
        &+ {\mathrm{i}}\left(R(\vec{x}) - \Gamma \left|\Psi\right|^2\right) \Psi \ .\end{aligned}$$ The first part of the right hand side is the ordinary equilibrium GPE with $m$ as effective polariton mass of the lower polariton branch, $V(\vec{x})$ as external potential, and $U$ as repulsive onsite interaction potential. The second part models phenomenologically the gain and loss of condensed polaritons. Here, $R(\vec{x})$ describes the linear part of gain and loss, and the non-linearity implements a density dependent gain saturation with $\Gamma$ as gain depletion parameter. This provides a simplified description of the gain process from a reservoir, for example relaxation of high-momentum polaritons generated by incoherent excitation with an external laser beam, and the condensate decay due to its finite lifetime. Since the non-condensed polaritons have a short lifetime as compared to the lifetime of the condensate, we can safely neglect diffusion processes of these and relate the spatial extension of the reservoir to the Gaussian excitation profile of the laser beam. Then, $$\begin{aligned}
    \label{eq:reservoir}
    R(\vec{x}) = \hbar \gamma_{{\mathrm{c}}} \left(\frac{P}{P_{{\mathrm{th}}}} e^{- \vec{x}^{\:2}/\xi_P^2} -1 \right) \ ,\end{aligned}$$ with decay rate $\gamma_{{\mathrm{c}}} = 1/\tau$, where $\tau$ is the condensate lifetime, and waist size $\xi_P$ of the laser beam. The parameter  $P/P_{{\mathrm{th}}}$ is the excitation power normalized by its value at the threshold at which condensation is observed first. The disorder landscape is incorporated by the external potential $V(\vec{x})$. We use a $\delta$-correlated Gaussian distributed quenched disorder with vanishing mean value and variance, $$\begin{aligned}
    \label{eq:DisCorr}
    {\langle\!\langle V(\vec{x}) \rangle\!\rangle} = 0\ , \quad {\langle\!\langle V(\vec{x})V(\vec{y}) \rangle\!\rangle} = V_0^2 \xi_V^2 \delta(\vec{x}-\vec{y}) \ ,\end{aligned}$$ respectively. The average disorder strength is given by $V_0$ and its characteristic length is denoted by $\xi_V$. In the case of a spatially homogeneous excitation, i.e. $\xi_{{\mathrm{P}}} \to \infty$, our model  was first suggested in Ref. . As compared to Ref.  we do not consider the dynamics of the reservoir polaritons explicitly. However, the latter can be eliminated [@Carusotto.2013] for the typical case that the characteristic relaxation rate of the reservoir is much faster than the condensate decay rate [@Wouters.2008; @Carusotto.2013]. Then, an expansion to leading order in condensate density over reservoir density results in the eGPE . We note that a different theoretical approach may be suitable in order to describe propagation of a polariton BEC in a disorder-free environment [@Wertz.2012; @Malpuech.2014], which is not the aim of our work.

In the following we will discuss the model . The mean condensate density $n_0\equiv \frac{1}{\Omega_{{\mathrm{c}}}} \int_{\Omega_{{\mathrm{c}}}} |\Psi(\vec{x})|^2$ is found by averaging the second term of the right hand side of [Eq. ]{} over the condensate area $\Omega_{{\mathrm{c}}} \approx \pi \xi_{{\mathrm{P}}}^2$, and then demanding a balance of gain and loss, $$\begin{aligned}
            \label{eq:defn0}
    n_0 \approx \frac{\hbar \gamma_{{\mathrm{c}}}}{\Gamma} \left( \frac{P}{P_{{\mathrm{th}}}} -1 \right) \ .\end{aligned}$$ Since the interaction term in [Eq. ]{} is proportional to the density, we find an energy blueshift $n_0 U$. The healing length is extracted by comparing the kinetic energy term versus the interaction term in [Eq. ]{}, $$\begin{aligned}
  \label{eq:defhealinglength}
    \xi \equiv \sqrt{\frac{\hbar^2/2m}{n_0 U}} \ .\end{aligned}$$ Let us understand its physical relevance: For example, we assume a region in which the condensate has to vanish, $\Psi = 0$, however remains unperturbed everywhere else. Then, the healing length is the distance over which the condensate density changes from zero to $n_0$.

A dimensionless eGPE  takes the form $$\begin{aligned}
    \label{eq:eGPEdimless}
    {\mathrm{i}}\partial_t \psi = (-\nabla^2 + \vartheta(\vec{x}) + |\psi|^2) \psi + i \alpha( g_{\mathrm{R}}(\vec{x}) - |\psi|^2)\psi \ ,\end{aligned}$$ where density, length, energy and time are measured in units of $n_0$, $\xi$, $n_0 U$ and $\hbar/n_0U$, respectively. The ’non-equilibrium’ parameter $\alpha$ and the dimensionless reservoir function $g_{\mathrm{R}}$ are defined in [Eq. ]{} and [Eq. ]{}, respectively. With $\psi(\vec{x},t) \equiv \Psi(\vec{r},t)/\sqrt{n_0}$ we denote the dimensionless wave function and $\vartheta(\vec{x}) = V(\vec{x})/ n_0U$ is the disorder potential relative to the blueshift with $$\begin{aligned}
  \label{eq:disorder_stat}
    {\langle\!\langle \vartheta(\vec{x}) \rangle\!\rangle} = 0 \ , \qquad {\langle\!\langle \vartheta(\vec{x})\vartheta(\vec{y}) \rangle\!\rangle} = \kappa^2 \delta(\vec{x}-\vec{y}) \ .\end{aligned}$$ We have introduced two important dimensionless parameters, namely an effective disorder strength and a ’non-equilibrium’ parameter $$\begin{aligned}
  \label{eq:defkappaalpha}
  \kappa \equiv \frac{\xi_V \: V_0}{\xi\: n_0U} \ , \quad \text{and} \quad \alpha \equiv \frac{\Gamma}{U} \ .\end{aligned}$$ The first parameter $\kappa$ is also obtained by coarse graining the random disorder potential up to the healing length (assuming $\xi_V \ll \xi$). This process renormalizes the disorder strength by a factor $1/\sqrt{(\xi/\xi_V)^2}$. Then, the value $\xi_V V_0 / \xi$ is compared to the blueshift $n_0 U$. The second parameter $\alpha$ implements the non-equilibrium nature of polaritons. In the limit $\alpha \to 0$ (keeping $n_0$ constant) [Eq. ]{} reduces to the equilibrium GPE, whereas, for $\alpha \to \infty$ the condensate is totally dominated by gain and loss. The rescaled reservoir function yields $$\begin{aligned}
    \label{eq:defgP}
    g_{\mathrm{R}}(\vec{x}) =  \frac{(P/P_{{\mathrm{th}}}) \: e^{-x^2/x_P^2}-1}{P/P_{{\mathrm{th}}}-1} \ ,\end{aligned}$$ with $x_{{\mathrm{P}}} \equiv \xi_P/\xi$. For a steady state solution (single-mode condensate) we make the ansatz $$\begin{aligned}
  \label{eq:ssansatz}
    \psi(\vec{x},t) = \psi(\vec{x}) e^{-{\mathrm{i}}\omega t} = \sqrt{n(\vec{x})} e^{{\mathrm{i}}\phi(\vec{x}) - {\mathrm{i}}\omega t} \ ,\end{aligned}$$ where $\hbar \omega$ is the condensate energy.

We emphasize that both blueshift and healing length depend on the excitation power $P$ via $n_0$. Thus, $\kappa$ and $x_P$ depend on $P$, too. For our analysis it is useful to identify energy and length scale which are excitation power independent, namely the line width energy $\hbar \gamma_{{\mathrm{c}}}$ and the quantum correlation length $l_\gamma \equiv \sqrt{\hbar / 2 m \gamma_{{\mathrm{c}}}}$ (a non-equilibrium analogon of the thermal de Broglie wavelength) [@Trichet.2013], so that $\kappa$ and $x_P$ become functions of $\alpha, P/{P_{{\mathrm{th}}}}$ and sample parameters (see Table \[tab:parameters\]).

Numerical Simulations and Comparison with the Experiment {#sec:theory_numsim}
========================================================

*Numerical simulations –* Computing the condensate wave function by solving the eGPE  allows us to extract the real and $k$-space intensity distribution. We define the $k$-space intensity distribution according to $$\begin{aligned}
       \label{eq:defIk}
        I_P(\vec{k}) &\equiv \gamma_{{\mathrm{c}}} n_0 \: |\psi({\vec{k}})|^2\ ,\end{aligned}$$ where the momentum space wave function is defined via a two-dimensional discrete Fourier transform $\psi(\vec{k}_j) = (1/N^2) \sum_{\vec{x}_i} \psi(\vec{x}_i) e^{-{\mathrm{i}}\vec{k}_j \vec{x}_i}$, with $\vec{x}_i,\vec{k}_j$ being elements of a discrete lattice with $N$ lattice points in each spatial direction, such that $i,j = 1,\ldots,N^2$. We choose an appropriate set of simulation parameters extracted from the experiment (see Table \[tab:parameters\]) and solve [Eq. ]{} numerically. To this end we look for a steady state solution, see [Eq. ]{}, by solving the time evolution of the discretized wave function $\psi(\vec{x}_i,t)$. The latter is defined on a real-space square lattice with spacing $a=\xi$. We employ a variable order Adams-Bashforth-Moulton PECE algorithm [@NR.2007] to obtain the time evolution. First, we compute the steady state solution of the disorder-free system, $\vartheta = 0$. Then, we choose independent Gaussian distributed variables of vanishing mean and variance $\kappa^2$ for each lattice site and calculate the steady state solution of the disordered system. The time evolution of the disordered system is started from the disorder-free solution as initial condition. For each disorder realization the discretized two-dimensional wave function $\psi(\vec{x}_i)$ of the steady state is extracted, and then Fourier transformed in order to compute the two-dimensional $k$-space intensity $I_P(\vec{k}_j)$. Finally, we average over all disorder realizations and compute the expectation value and variance, $$\begin{aligned}
    \label{eq:meanIk}
    \mu_P(k) &\equiv {\langle\!\langle I_P(\vec{k}) \rangle\!\rangle}/{\langle\!\langle I_P(0) \rangle\!\rangle} \ ,\\      \label{eq:varIk}
    \sigma^2_P(k)  &\equiv {\langle\!\langle \left( I_P(\vec{k})-{\langle\!\langle I_P(\vec{k}) \rangle\!\rangle} \right)^2 \rangle\!\rangle}/{\langle\!\langle I_P(0) \rangle\!\rangle}^2\ ,\end{aligned}$$ respectively. Above, the bracket ${\langle\!\langle \ldots \rangle\!\rangle}$ denotes an average with respect to disorder, and we have normalized the mean and variance by the expectation value of the intensity at $k=0$. Since the excitation profile [Eq. ]{} is radial symmetric, Eqs. (\[eq:meanIk\], \[eq:varIk\]) are radial symmetric, too, assuming a sufficiently large number of disorder realizations.

![(color online) Experimental $I_P(k)$ distribution (blue diamonds) compared to numerical simulations of scenario I for increasing excitation power $P$. The black solid lines depict the expectation value whereas the gray band indicates the standard deviation. For high momenta we have excluded systematically biased data (gray crosses). For each plot the mean squared error (MSE) and the goodness-of-fit value $\mathcal{Q}$ are computed. Used simulation parameters: $\alpha = 0.5$, $\xi_{{\mathrm{V}}} V_0 / l_{{\mathrm{c}}} \hbar \gamma_{{\mathrm{c}}} = 0.4$, $\xi_P/l_{{\mathrm{c}}}=3$. These correspond to ${\mathcal{L}_{\mathrm{n}}}/{L_{c}}\approx 1.5$, ${\mathcal{L}_{\phi}}/{L_{c}}\approx 13$, ${\mathcal{L}_{s}}/{L_{c}}\approx 26$. []{data-label="fig:FitScenarioI"}](./{fig_figfit_T10_detuning-30_cP3_ckappa0.4_alpha0.5}.pdf){width="1.00\hsize"}

*Comparison with the experiment –* The numerically obtained mean and variance of the $k$-space intensity can be compared with the experimental data denoted by $I_{{\mathrm{ex}}}(k_x)$ here, cf. Sec. \[sec:experimental\_results\]. We note that these measurements represent a line-cut along an axis (e.g. the $x$-axis) of the two-dimensional intensity distribution and are measured for one disorder configuration determined by the disorder of the sample. We perform a spatial averaging step by symmetrizing the experimentally obtained intensity: $I_{{\mathrm{ex}}}(k_x) \to (I_{{\mathrm{ex}}}(k_x) + I_{{\mathrm{ex}}}(-k_x))/2$ and $k_x \geq 0$. In order to quantify the agreement between experiment and theory we introduce the chi-square value [@NR.2007] $$\begin{aligned}
    \label{eq:DefChiSqr}
    \chi^2_P = \sum \limits_{k_j \geq 0} \left( \frac{I_{{\mathrm{ex}}}(k_j)/a - \mu_{P}(b k_j)}{\sigma_{P}(b k_j)} \right)^2 \ .\end{aligned}$$ Since the condensate density $n_0$ and the healing length $\xi$ are hard to extract experimentally, we use two scaling parameters ($a,b$) instead. Both are determined by a least-squares fitting procedure. [@NR.2007]

In order to estimate the goodness-of-fit [@NR.2007] we extract the complement of the $\chi^2$-probability distribution function $F_{\chi^2}$, denoted by $\mathcal{Q} = 1-F_{\chi^2}(\chi^2_P)$, which is the probability that the simulations agree with the experimental data. If $\mathcal{Q} \ll 1$, the apparent discrepancies of model and data are unlikely to be random fluctuations, and we conclude that the model is not specified correctly, or that the fluctuation strength $\sigma_P$ is under-estimated. On the other hand, if $\mathcal{Q} \sim 1$, we conclude that the model describes the data correctly. Finally, we define the mean squared error: ${\mathrm{MSE}} \equiv (1/N^2) \: \sum_{k_j} (I_{{\mathrm{ex}}}(k_j)/a - \mu_P(k_j))^2$, which is a measure of how well the data match the simulated intensity distribution. The comparison of the experimental data and the numeric simulations of scenario I is shown in Fig. \[fig:FitScenarioI\], and the comparison with simulations of scenario II was shown and discussed in Sec. \[sec:fit\_of\_exp\_data\], Fig. \[fig:Fit\].

Supplemental Material {#supplemental-material .unnumbered}
=====================

Experimental Determination of Condensation Threshold Density {#sec:SM_Pth}
============================================================

For the comparison between the theoretical model and the experiment as discussed in Sec. \[sec:fit\_of\_exp\_data\] of the main text the determination of the condensation threshold ${P_{{\mathrm{th}}}}$ is crucial, since experimental as well as numerically simulated spectra shall be compared for similar ratios of $P / {P_{{\mathrm{th}}}}$.

![(a)Far-field emission pattern $I(E,k)$ for $P =  79~{\mathrm{W}}/{\mathrm{cm}}^{2}$. The black dotted line represents the energy barycenter for the LPB emission $E^{\mathrm{BC}}_{\mathrm{LPB}}(k)$ for each $k$ value. The white dashed lines highlight all relevant spectral positions which are plotted for increasing excitation power density in (b). Here, for the LPB the minimum of the parabolic dispersion as well as the energy barycenter for $k \approx 0$ was analyzed. The empty red circles represent BEC emission that appears as shoulder within the PL intensity spectra leading to larger uncertainties for the determination of the peak energy. (cf. Fig S.\[fig:SM:peakArea\](a))[]{data-label="fig:SM_Intro_spectr_contr"}](./SM_Intro_spectr_contr.pdf){width="0.7\linewidth"}

In general, we identify three spectral contributions within our far-field PL emission pattern, as shown in Fig. S\[fig:SM\_Intro\_spectr\_contr\], namely the lower polariton branch (LPB) emission, and two Bose-Einstein condensate (BEC) emission channels. Whereas the minimum of the LPB dispersion is almost excitation power independent ($E_{\mathrm{LPB}}^{\mathrm{min}} = 3.3207$ eV), the LPB emission gets significantly broadened towards higher energies for increasing excitation power. This is due to the fact that polaritons are located at different regions and are subject to different blueshifts due to the spatially inhomogeneous, e.g. Gaussian pump spot as well as the pronounced disorder potential. To consider this effect quantitatively, we calculated the energy barycenter of the LPB emission $E^{\mathrm{BC}}_{\mathrm{LPB}}(k) = (\int_{E} I(E) \cdot E \; dE) / (\int_{E} I(E) dE)$ for $k \approx 0$. To investigate the excitation power dependent evolution of the BEC emission, we performed a lineshape analysis of the PL spectra which are integrated over all observed $k$ values by assuming a Lorentzian lineshape. Both BEC emission channels show a large energy shift of about $E_{\mathrm{BEC 1}} = E_{\mathrm{BEC 2}} = 16\pm2$ meV with respect to $E_{\mathrm{LPB}}^{\mathrm{min}}$ for increasing excitation power density up to $P =  155~{\mathrm{W}}/{\mathrm{cm}}^{2}$. We assume that the initial large energy shift of $E^{\mathrm{BC}}_{\mathrm{LPB}}$ as well as of $E_{\mathrm{BEC 1,2}}$ is mainly caused by an electronic background potential that is discussed in detail in Appendix \[sec:experiment\_disorderorigin\] in the main text. For $P >  155~{\mathrm{W}}/{\mathrm{cm}}^{2}$ the slope of $E_{\mathrm{BEC 1,2}}$ is reduced, following the expectations for common BEC. In this regime, the electronic background potential starts to saturate and the condensate-condensate interaction becomes dominant.

Unfortunately, the experimental determination of the condensation threshold is accompanied by large uncertainties due to the interaction between both condensate modes as well as the superposition of intense LPB emission for a large range of excitation powers. Therefore, we discuss here the impact of the disorder on the determination of the threshold power. By doing so, we analyze at first the evolution of the total PL intensity with increasing excitation power, as usually done for a disorder free condensate. As will be discussed below, this method gives only an upper limit and therefore we apply two further methods: firstly, we examine the excitation power dependent evolution of the PL intensity for each spectral contribution separately. Secondly, we study the PL spectra $I_k(E)$ for each $k$ value and excitation power separately and deduce the FWHM for the individual emission channels as a function of $k$ and excitation power density.\

#### Evolution of the total Photoluminescence (PL) Intensity {#par:totalInt}

![(color online) Excitation power dependence of the total PL intensity. For small excitation power densities $P \leq 235~{\mathrm{W}}/{\mathrm{cm}}^{2}$ the behavior can be described by two power functions with different exponents below (red solid line) and above (blue solid line) the condensation threshold. The threshold can be deduced from the intersection between both functions: ${P_{{\mathrm{th}}}}= 130{\mathrm{W}}/{\mathrm{cm}}^{-2}$.[]{data-label="fig:SM:totalInt"}](./SM_totalInt_2.pdf){width="0.7\linewidth"}

As a first guess, we analyzed the excitation power dependence of the total PL intensity, integrated over all observed $k$ values and energies, as shown in Fig. S\[fig:SM:totalInt\]. The slope of the PL intensity increases abruptly for $P \geq 130~{\mathrm{W}}/{\mathrm{cm}}^{2}$. By assuming a power law behavior,[@kavokin.2007] the exponent increases from 1.4 for $P < 130~{\mathrm{W}}/{\mathrm{cm}}^{2}$ to 4.3 for $P > 130~{\mathrm{W}}/{\mathrm{cm}}^{2}$. Note that the estimated value of ${P_{{\mathrm{th}}}}= 130~{\mathrm{W}}/{\mathrm{cm}}^{2}$ is only an upper limit for the condensation threshold. This can be explained by the inhomogeneous shape of the (e.g. Gaussian) excitation spot profile. For excitation powers $P \gtrsim {P_{{\mathrm{th}}}}$ the critical density for polariton condensation is achieved within a small area only. In contrast to this, emission from uncondensed polaritons occurs for a much larger area, which superimposes the BEC emission. Thus, the BEC emission becomes dominant leading to the observed kink in the evolution of the PL intensity with increasing powers only for powers significantly larger than the real condensation threshold.\

#### Evolution of PL Emission for each BEC State {#subsec:BEC_peakArea}

![(color online) (a) Normalized PL spectra $I(E)$ integrated over all $k$ values for different excitation powers. The blue (dark yellow) arrow points at the PL spectra for which firstly a pronounced peak relating to the emission of BEC 2 (BEC 1) is observed. (b) Selected spectra that are significant for the determination of ${P_{{\mathrm{th}}}}$. (c) Exemplary PL spectrum for $P =  155~{\mathrm{W}}/{\mathrm{cm}}^{2}$ showing the FWHM of both BEC emission channels as well as the peak area. (d) Excitation power dependence of the PL intensity, integrated over the FWHM of the corresponding peak.[]{data-label="fig:SM:peakArea"}](./SM_peakArea3.pdf){width="0.9\linewidth"}

To analyze the impact of the superposition of the LPB and BEC emission on the determination of ${P_{{\mathrm{th}}}}$, we analyzed the excitation power dependent evolution of the PL emission for both contributions separately. For this purpose, we investigated the intensity spectra $I(E)$, integrated over all observed $k$ values, for different excitation power densities as shown in Fig. S\[fig:SM:peakArea\](a). For the lowest density ($ P = 16~{\mathrm{W}}/{\mathrm{cm}}^{2}$) only the LPB emission can be observed, with maximum intensity at the minimum of the LPB dispersion of $E_{{\mathrm{LPB}}} = 3.3207~{\mathrm{eV}}$. For increasing excitation power, the high energy edge of the LPB emission dominates due to the previously mentioned significant broadening towards higher energies. Already for $P = 79~{\mathrm{W}}/{\mathrm{cm}}^{2}$ an additional emission channel appears at $E_{\mathrm{BEC 2}} = 3.335~{\mathrm{eV}}$. This peak becomes pronounced and shows a strong narrowing for $P \geq 109~{\mathrm{W}}/{\mathrm{cm}}^{2}$. For further increasing excitation power a second pronounced BEC emission channel appears within the $I(k)$ spectra at smaller energies indicating the multimode BEC behavior, as discussed in Sec. \[sec:experimental\_results\] in the main text. Note that this emission channel is already observable for $P = 79~{\mathrm{W}}/{\mathrm{cm}}^{2}$ in the energy resolved $k$-space images as shown in Fig. S\[fig:SM\_Intro\_spectr\_contr\](a). However, it appears only as a small shoulder within the $k$-integrated intensity spectra for this excitation density range (cf. Fig. S\[fig:SM:peakArea\](a,b)). The energy position of the spectral contributions considered here (cf. Fig. S\[fig:SM\_Intro\_spectr\_contr\](b)) is indicated by the dashed lines in Fig. S\[fig:SM:peakArea\](a).

An exemplary lineshape analysis of the BEC 1 and BEC 2 contributions to the PL spectra is shown in Fig. S\[fig:SM:peakArea\](c) for an excitation power of $P = 155~{\mathrm{W}}/{\mathrm{cm}}^{2}$. The dependence on excitation power density of the PL peak area, integrated over the FWHM of the corresponding Lorentzian peaks, for both BEC emission channels separately as well as for the sum of both is compiled in Fig. S\[fig:SM:peakArea\](d) in a double-logarithmic scale. The PL intensity of BEC 2 starts to saturate for $P \approx 300~{\mathrm{W}}/{\mathrm{cm}}^{2}$ and even decreases for $P >  550~{\mathrm{W}}/{\mathrm{cm}}^{2}$, whereas the PL intensity of BEC 1 further increases for the excitation power density range shown here. This indicates an effective relaxation of polaritons from the high-energy BEC state 2 into the low-energy BEC state 1. Thus, both BEC emission channels are not independent but represent a system of coupled condensate states for which we estimate a single condensation threshold density. Similar to the previous method, we expect a kink in the evolution of the PL intensity at ${P_{{\mathrm{th}}}}$, however, only a discontinuity is barely visible at $P \approx 84~{\mathrm{W}}/{\mathrm{cm}}^{2}$ for the data set presented here (cf. Fig. S\[fig:SM:peakArea\](d)). We note that the PL spectra for $P <  109~{\mathrm{W}}/{\mathrm{cm}}^{2}$ do not show a clear peak for the energy range of the expected BEC emission due to the spectral overlap with the intense and spectrally broad LPB emission. Thus, the extracted peak area is subjected to large uncertainties in the mentioned range of excitation densities, which is highlighted by the gray dashed lines in Fig. S\[fig:SM:peakArea\](d). Due to this fact, the observed discontinuity in the evolution of the PL intensity is not fully reliable and only a range of possible values can be determined for ${P_{{\mathrm{th}}}}$. On the one hand, a peak at 3.335 eV is slightly visible for an excitation density of $P =  79~{\mathrm{W}}/{\mathrm{cm}}^{2}$, which leads to the observed discontinuity in the PL intensity evolution at $P \approx 84~{\mathrm{W}}/{\mathrm{cm}}^{2}$ and may indicate the onset of BEC emission. However, this peak may also be caused by effective polariton scattering into an already blueshifted state in the uncondensed regime. On the other hand, the peak gets pronounced and spectrally narrowed for $P =  109~{\mathrm{W}}/{\mathrm{cm}}^{2}$. This is a clear signature for polariton condensation. Conclusively, the range of values for ${P_{{\mathrm{th}}}}$ can be restricted to ${P_{{\mathrm{th}}}}= (84 - 109)~{\mathrm{W}}/{\mathrm{cm}}^{2}$ by using this method. A way to additionally reduce the impact of spectral overlapping between LPB and BEC emission and further specify ${P_{{\mathrm{th}}}}$ is discussed in the following section.\

#### $k$-dependent Evolution of FWHM {#par:FWHM_vs_P}

![(color online) Far-field emission pattern $I(E,k)$ for an exemplary excitation power of $P =  141~{\mathrm{W}}/{\mathrm{cm}}^{2}$. The parabolic shaped LPB emission as well as the dispersionless BEC emission channels are highlighted by white dashed lines. []{data-label="fig:SM:PL_Image"}](./SM_PL_Image.pdf){width="0.7\linewidth"}

For a more sophisticated investigation we analyzed the PL intensity spectra $I_k(E)$ for each $k$ value and excitation power separately. Thereby, we deduced the FWHM for both BEC emission channels as well as for the LPB emission. Fig. S\[fig:SM:PL\_Image\] shows exemplarily a far-field PL emission pattern for $P =  141~{\mathrm{W}}/{\mathrm{cm}}^{2}$ with logarithmic intensity scale. Here, all three emission channels, marked by white dotted lines, are energetically well separated for a large range of $k$ values.

![(color online) $k$-dependent FWHM of LPB and BEC emission for a large range of excitation power densities. Empty circles represent excitation power densities $P < {P_{{\mathrm{th}}}}$, whereas excitation power densities $P \geq {P_{{\mathrm{th}}}}$ are highlighted by full ones.[]{data-label="fig:SM:FWHM"}](./SM_FWHM.pdf){width="\linewidth"}

The $k$-dependent evolution of the FWHM with increasing excitation power density is shown in Fig. S\[fig:SM:FWHM\] for each emission channel separately. Note, that the missing data points correspond to PL spectra which show a strong spectral overlap of the emission channels and thus prevent a proper lineshape analysis. The broadening of the LPB emission increases with increasing absolute $k$ values due to an increasing excitonic fraction. For $P \geqslant  79~{\mathrm{W}}/{\mathrm{cm}}^{2}$ the FWHM of the BEC peak 2 is lower than the minimum FWHM of the LPB emission indicating the onset of polariton condensation. For the BEC peak 1 this situation is present for $P \geqslant  109~{\mathrm{W}}/{\mathrm{cm}}^{2}$. For both BEC emission channels, the narrowing saturates for $P =  141~{\mathrm{W}}/{\mathrm{cm}}^{2}$, indicating a maximum temporal coherence. Following the arguments about the interaction between both investigated BEC emission channels as discussed in the previous section, this coupled BEC system is characterized by a single threshold power of ${P_{{\mathrm{th}}}}=  79~{\mathrm{W}}/{\mathrm{cm}}^{2}$.\

#### Summary

The results for the determination of ${P_{{\mathrm{th}}}}$ are summarized for all three methods in Table \[tab:Pth\]. Using the typically used method by analyzing the total PL intensity as a function of excitation power, an upper value of ${P_{{\mathrm{th}}}}<  130~{\mathrm{W}}/{\mathrm{cm}}^{2}$ was estimated. By studying the PL intensity for each spectral contribution separately it was possible to reduce the impact of the disorder on the determined threshold value and further restrict the range of ${P_{{\mathrm{th}}}}$ to ${P_{{\mathrm{th}}}}= (84- 109)~{\mathrm{W}}/{\mathrm{cm}}^{2}$. The best minimization of the disorder influence on the determination of the threshold power for condensation was achieved by investigating the PL spectra for each $k$ value separately and deducing the FWHM for each spectral contribution. This method also differs from the other ones regarding its physical principle. The analysis of the PL intensity evolution for the total emission as well as for the individual BEC emission are based on an increasing rate of the parametric scattering process into the condensate state for $P \geq {P_{{\mathrm{th}}}}$ due to its bosonic nature. Thereby, the coexistence of a small area of condensed polaritons for $P \gtrapprox {P_{{\mathrm{th}}}}$ and a large area of uncondensed ones can cause a rather soft transition of the PL intensity evolution leading to large uncertainties for the determination of ${P_{{\mathrm{th}}}}$. Investigating the FWHM of the BEC emission channels rather gives insight into another property. As the FWHM of the BEC emission channels are inversely proportional to the temporal coherence of the corresponding system of particles the spontaneous build-up of coherence is observed, which is a basic property of a polariton condensate. Nevertheless, the convolution of a certain BEC emission peak with other emission channels leads to uncertainties for the determination of ${P_{{\mathrm{th}}}}$ using this method, too. In summary, we estimate ${P_{{\mathrm{th}}}}=  79~{\mathrm{W}}/{\mathrm{cm}}^{2}$ as the threshold value for polariton condensation in our MC for the investigated parameter set of $T = 10~{\mathrm{K}}, \Delta = -30~{\mathrm{meV}}$.

  method                                       ${P_{{\mathrm{th}}}}$ in ${\mathrm{W}}/{\mathrm{cm}}^{2}$   comments
  -------------------------------------------- ----------------------------------------------------------- -----------------------------------------------------------------------------------------------------------
  total PL intensity                           130                                                         upper limit, BEC emission superimposed by intense LPB emission
  PL intensity for each BEC emission channel   84 - 109                                                    superposition with LPB emission prevents reliable lineshape analysis for small excitation power densities
  $k$-dependent FWHM                           79                                                          provides best separation between LPB and BEC emission

  : Used methods for experimental determination of the condensation threshold density ${P_{{\mathrm{th}}}}$.[]{data-label="tab:Pth"}

Experimental limitations for Determination of the Coherence Time {#sec:SM_limits_cohTime}
================================================================

In this section, we discuss two experimental artifacts that lead to limitations in the determination of the condensate’s coherence time.

Spectrally dependent Phase Shift {#subsec:Tcoh_phase_shift}
--------------------------------

The phase difference $\phi_{12}$ between the emission from both interferometer arms depends on the emission wavelength $\lambda$ in the following way: $$\begin{aligned}
    \label{eq:interf_phase}
\phi_{12}(\lambda,\Delta s,\vec{r}) & =  \frac{2 \pi \Delta s} {\lambda} + (\vec{k_1} - \vec{k_2}) \vec{r} \notag \\ & = \frac{2 \pi \Delta s} {\lambda} +  \frac{2 \pi} {\lambda} |\vec{r}| (\sin(\alpha_1) - \sin(\alpha_2)).\end{aligned}$$ Here, $\Delta s$ is the path length difference between both interferometer arms, $\vec{k_1},\vec{k_2}$ are the wavevectors of the waves propagating along the individual interferometer arms, $\alpha_1,\alpha_2$ are the angle between the optical axis and the corresponding wavevectors $\vec{k_1},\vec{k_2}$ and $\vec{r}$ is the distance vector from the intersection point between both wavevectors (in the detector plane) and the point of interest of the resulting interference pattern. The geometry is sketched in Fig. \[fig:scheme\_PhaseDiff\] for the special case of $\alpha_2 = 0$. The second term of [Eq. ]{} defines the appearance of the observed interference pattern (fringe distance and orientation) for a specified $\Delta s$ and $\lambda$, whereas the first term causes an additional, spectrally dependent phase offset that increases linearly with increasing path length difference.

![(color online) Scheme of two superimposed, monochromatic plane waves with wavevectors $\vec{k}_1 = (2 \pi / \lambda) (- \sin (\alpha_1) \vec{e}_x + \cos (\alpha_1) \vec{e}_z) ,\vec{k}_2 = (2 \pi / \lambda) \vec{e}_z$ and equal wavelength $\lambda$ to illustrate the spatial phase distribution $\phi_{12}( \vec{r})$ (cf. second term in [Eq. ]{}) for the special case of $\alpha_2 = 0$. The solid (dashed) lines represent electric field maxima (minima) of the corresponding waves. Interference maxima (minima) within the detection plane occur at position $\vec{r}$ for which the electric field maximum of wave 2 is superimposed with a field maximum (minimum) of wave 1. The spatial phase distribution in the detection plane $\phi_{12}(x,y)$ = $\phi_{12}(x) = (2 \pi |\vec{r}| \sin(\alpha_1)) / \lambda$ is in accordance with the second term in [Eq. ]{} for $\alpha_2 = 0$. $\phi_{12}(x)$ can be found by relating the projection $|\vec{r}| \sin(\alpha_1)$ with the emission wavelength $\lambda$. This leads to a fringe period of $\Delta |\vec{r}| = \Delta x = \lambda /  \sin(\alpha_1)$. []{data-label="fig:scheme_PhaseDiff"}](./SM_scheme_PhaseDiff.pdf){width="0.8\linewidth"}

The spectral resolution is $\Delta \lambda = \lambda_{\mathrm{max}} - \lambda_{\mathrm{min}} = 0.08$ nm for our experiment. Therefore, accumulation of the intensity of the interference pattern over the spectral range $\Delta \lambda$ is accompanied by an integration over a range of phase differences

$$\begin{aligned}
\label{eq:interf_phasediff}
\Delta \phi_{12}(\lambda,\Delta s) = & 2 \pi   \left(  \frac{1}{\lambda_{\mathrm{min}}} - \frac{1}{\lambda_{\mathrm{max}}} \right)  \notag \\ & (\Delta s + |\vec{r}| (\sin(\alpha_1) - \sin(\alpha_2)).\end{aligned}$$

The second term in the second bracket in [Eq. ]{} can be neglected due to $|\vec{r}| (\sin(\alpha_1) - \sin(\alpha_2)) \ll \Delta s$ for almost the total range of path differences used here. The angle between the propagation directions of the two superimposed waves is about $0.14^\circ$ deduced from the fringe distance of the observed interference pattern, which represents an upper limit for $\sin(\alpha_1) - \sin(\alpha_2) \approx \alpha_1 - \alpha_2$ (in case of opposite signs for both angles). The radius of the observed interference pattern is in the range of $|\vec{r}|_{\mathrm{max}} \approx 1$ mm leading to $|\vec{r}|_{\mathrm{max}} \Delta \alpha_{\mathrm{max}} \approx 2.4~\mu$m. This is about three orders of magnitude below the maximum path length difference of $\Delta s_{\mathrm{max}} = 2.7$ mm.

The first term in [Eq. ]{} can be re-expressed in terms of $\Delta \phi_{12}(\lambda,\Delta s) =  2 \pi   \left(  \frac{1}{\lambda_{\mathrm{min}}} - \frac{1}{\lambda_{\mathrm{max}}} \right) \Delta s = 2 \pi   \left(  \frac{\lambda_{\mathrm{max}} - \lambda_{\mathrm{min}}}{\lambda_{\mathrm{max}} \lambda_{\mathrm{min}}} \right) \Delta s \approx   2 \pi \frac{\Delta \lambda}{\lambda^2} \Delta s$. Therefore, a significant phase shift of the order of $2\pi$ is induced, if the path difference $\Delta s$ is of the order of $\lambda^2 / \Delta \lambda$. This is the case for the measurement presented here, since $\lambda^2 / \Delta \lambda \approx 1.7$ mm. Thus, the experimentally observed decrease of the normalized visibility with increasing path difference (or temporal delay) is stronger than the pure reduction due to the decreasing temporal coherence $g^1(\Delta t)$, which is therefore under-estimated.

In order to quantify the impact of the spectrally dependent phase shift on the determined visibility we integrate the interference pattern $I_{\mathrm{interf}}$ over the range of phase differences $\Delta \phi_{12} = \phi_{12,{\mathrm{max}}} - \phi_{12,{\mathrm{min}}}$ that corresponds to a single CCD row $$\begin{aligned}
    I_{\mathrm{interf}}(\phi_{12,{\mathrm{min}}},\phi_{12,{\mathrm{max}}},\Delta s,\vec{r})  \notag \\ =  \int_{\phi_{12,{\mathrm{min}}}}^{\phi_{12,{\mathrm{max}}}} \! \frac{ I_{\mathrm{interf}}(\phi_{12},\Delta s,\vec{r}) \, \mathrm{d} \phi_{12}}{\phi_{12,{\mathrm{max}}} - \phi_{12,{\mathrm{min}}}}    
    \end{aligned}$$ where $\phi_{12}$ is given by Eq. \[eq:interf\_phase\]. Thereby, we assume a constant intensity distribution of the emission from the individual arms within the width of the single CCD row as well as total coherence between both intensity signals. Afterwards, we determined the normalized visibility according to Eq. \[eq:Inorm\] in the main text.

![(color online) (a) Decrease of the normalized visibility with increasing temporal delay of a theoretically assumed, totally coherent signal due to the spectrally dependent phase shift. The small fluctuations of the numerically determined data are caused by the finite number of points for the spatial interference pattern $I_{\mathrm{interf}}(\vec{r})$. (b) Corrected normalized visibility $V_{\mathrm{norm}}'$ as a function of the temporal delay between both intensity signals.[]{data-label="fig:SM_phaseShift_correction"}](./SM_phaseShift_correction.pdf){width="\linewidth"}

Following this procedure, we found empirically a sinusoidal decrease of $V_{\mathrm{norm}}$ with increasing path length difference or rather temporal delay, as shown in Fig. S\[fig:SM\_phaseShift\_correction\](a). By taking into account this systematic experimental error, a corrected coherence time of 10.3 ps could be deduced (cf. Fig. S\[fig:SM\_phaseShift\_correction\](b)).

Limited Excitation Pulse Width
------------------------------

For the coherence measurements we used pulsed excitation with a pulse length of about 2 ps. Note that for the PL experiments which are compared to theoretical simulations based on a steady-state theory, a different excitation laser with pulse length of about 500 ps was used. By means of previous time-resolved measurements of the investigated microcavity (MC) under similar excitation conditions, a condensate lifetime of about 4-8 ps was observed (not shown in this work). However, the PL intensity decreases exponentially after the excitation pulse vanishes, which strongly limits the determination of longer coherence times.

![(color online) (a) Monoexponential decay of two intensity signals $i_1(t)$ and $i_2(t)$ as well as of the superimposed interference signal $i_{\mathrm{interf}}(t)$. For the latter one, total coherence $g^1(\Delta t)$ and vanishing phase shift $\Delta \phi_{12}$ between both individual signals are assumed. (b) Comparison between the experimentally obtained normalized visibility and the adapted model $V_{{\mathrm{meas}}}  = V_{\mathrm{norm}}^{\mathrm{BEC}} \; V_{{\mathrm{norm}}}'$ (cf. Eq. \[eq:tau\_coh\].) Here, we considered the decrease of the calculated normalized visibility $V_{{\mathrm{norm}}}'$ or rather normalized temporal coherence with increasing temporal delay between both signals due to the finite pulse length of the condensate emission. We assumed a condensate lifetime of $\tau_{\mathrm{LT}} = 6~{\mathrm{ps}}$. (c) Similar procedure as in (b), but an additional correction due to the spectrally dependent phase shift was applied.[]{data-label="fig:SM_pulsedExc"}](./SM_pulsedExc.pdf){width="\linewidth"}

To quantify the impact of the finite pulse duration on the calculated coherence time, we consider two pulses $i_1(t)$ and $i_2(t)$ originating from both interferometer arms with equal amplitude $i_0$ and with a temporal delay of $\Delta t$. For simplicity, we describe the temporal evolution of both pulses by a mono-exponential decay

$$\begin{aligned}
\label{eq:exp_decay}
i_1(t) &= i_0 \:  \exp(-t/\tau_{\mathrm{LT}}) \\
i_2(t) &= \begin{cases} 0 & \text{for} \: t < \Delta t \\                         
                        i_0 \:  \exp(- (t+ \Delta t)/\tau_{\mathrm{LT}})  & \text{for} \: t \geq \Delta t  \end{cases}  
$$

while neglecting the onset time as shown in Fig. \[fig:SM\_pulsedExc\](a). For $t \geq \Delta t$ the intensity of the delayed signal $i_2(t)$ is by a factor of $A \equiv \exp(\Delta t / \tau_{\mathrm{LT}})$ larger than $i_1(t)$. The detection occurs time-integrated over millions of laser pulses, whereas each pulse acts as an individual statistical event. Thus, we have to consider the PL intensity integrated over the time interval between two consecutive pulses $T$:

$$\begin{aligned}
\label{eq:exp_decay_sum}
I_1 &= \int_0^T i_1(t) \: d{\mathrm{t}} & \overset{T \gg \tau_{\mathrm{LT}}}= i_0 \tau_{\mathrm{LT}}\\
I_2 &= \int_{\Delta t}^T i_2(t) \: d{\mathrm{t}}   & \overset{T \gg \tau_{\mathrm{LT}}}= i_0 \tau_{\mathrm{LT}} .
$$

The condition $T \gg \tau_{\mathrm{LT}}$ is fulfilled in our experiment, since the time interval $T$ is 13 ns, which is about three orders of magnitude larger than the lifetime of the condensate of about $\tau_{\mathrm{LT}} = 6~ {\mathrm{ps}}.$ To determine the coherence time, we calculate the normalized visibility $V_{{\mathrm{norm}}}$ of the interference pattern (cf.Eq. \[eq:Inorm\] in the main text), whose amplitude represents the temporal first order correlation function $g^1(\Delta t)$. This condition is only valid, if we calculate $V_{{\mathrm{norm}}}(t)$ for each point in time during the temporal decay of the PL intensity separately. But in our experiment, we firstly measure the temporally integrated intensity of the interference pattern $ I_{{\mathrm{interf}}} = \int_0^T  i_{{\mathrm{interf}}}(t) \: d{\mathrm{t}} = \int_0^T  (i_1(t) + i_2(t) + 2\sqrt{i_1(t) i_2(t)} g^1(\Delta t) \cos(\Delta \phi_{12}) ) \: d{\mathrm{t}}$ and calculate the normalized visibility $V_{\mathrm{norm}}'$ afterwards. To quantify the impact of the pulsed excitation we calculate $V_{{\mathrm{norm}}}'$ for an assumed superposition of two totally coherent signals without any phase shift ($g^1(\Delta t)~=~1,~\cos(\Delta \phi_{12} = 0$)). This leads to the following equation: $$\begin{aligned}
\label{eq:interf_tenp_integration}
I_{{\mathrm{interf}}}  &= \int_0^{\Delta t} i_1(t) \: d{\mathrm{t}} + \int_{\Delta t}^T  (i_1(t) + i_2(t) + \sqrt{i_1(t) i_2(t)} \: d{\mathrm{t}}  \notag \\
&= \int_0^{\Delta t} i_1(t) \: d{\mathrm{t}} + \int_{\Delta t}^T  
i_2(t) [\frac{1}{A} + 1 + \frac{2}{\sqrt{A}}] \: d{\mathrm{t}}  \notag  \\
&= i_0 \tau_{\mathrm{LT}} \: (1 - \frac{1}{A}) + i_0 \tau_{\mathrm{LT}} [\frac{1}{A} + 1 + \frac{2}{\sqrt{A}}] \notag \\
&= i_0 \tau_{\mathrm{LT}} \: (2 + \frac{2}{\sqrt{A}})  
$$ with $A$ being the intensity factor between both signals as defined previously in this paragraph. This value of $I_{{\mathrm{interf}}}$ is smaller than the expected maximum intensity for the superposition of two totally coherent signals with equal intensity of $I_{{\mathrm{interf}}}^{{\mathrm{max}}} = 4 i_0 \tau_{\mathrm{LT}}$, except for the trivial case of A = 1 that is fulfilled for $\Delta t  = 0$ only. Consequently, this leads also to a reduction of the normalized visibility $$\begin{aligned}
V_{{\mathrm{norm}}}' & = \frac{I_{{\mathrm{interf}}} - I_1 - I_2}
    {2\sqrt{I_1 I_2}}   \notag \\ & = \frac{1}{\sqrt{A}} = \exp(-\frac{\Delta t}{2 \tau_{\mathrm{LT}}})  & \leq 1 .\end{aligned}$$

Fig. \[fig:SM\_pulsedExc\](b) shows the evolution of the calculated normalized visibility $V_{{\mathrm{norm}}}'$ as a function of the temporal delay for $\tau_{\mathrm{LT}} = 6$ ps. Since we assumed total coherence and a vanishing phase shift between both signals $i_1(t)$ and $i_2(t)$, the reduction of $V_{{\mathrm{norm}}}'$ is exclusively caused by the temporal decay of the condensate emission due to the pulsed excitation. Therefore, $V_{{}{\mathrm{norm}}}'$ represents a correction function for the real value $V_{{}{\mathrm{norm}}}$. For the experimentally determined, uncorrected coherence time of $\tau_{\mathrm{coh}}' = 8.7~$ps we find a reduction of the normalized intensity by a factor of $V_{{}{\mathrm{norm}}}' = 0.49$. For comparison, the maximum value extracted from the experiment is about $V_{{}{\mathrm{norm}}} \approx 0.6$ (cf. Fig. \[fig:EkSpace\_Interference\] in the main text), whereas for uncorrelated emission a residual normalized visibility of $V_{{}{\mathrm{norm}}} \approx 0.04$ could be estimated (not shown here). For our MC condensation can only be achieved with pulsed excitation, thus we cannot determine the real coherence time directly from the measurement. However, with the help of the simplified model presented here, a corrected value of the coherence time can be estimated.

The normalized visibility obtained from the experiment $V_{\mathrm{meas}}$, is a convolution of the normalized visibility of the investigated condensate $V_{\mathrm{norm}}^{\mathrm{BEC}}$, for which a Gaussian decay is assumed for increasing temporal delay $\Delta t$, and of the correction function $V_{{\mathrm{norm}}}'$, leading to the following equation: $$\begin{aligned}
V_{{\mathrm{meas}}} & = V_{\mathrm{norm}}^{\mathrm{BEC}} \; V_{{\mathrm{norm}}}'  \notag \\ & 
= g^1(\Delta t =0) \exp(- \frac{\pi}{2} \frac{\Delta t^2}{\tau_{{\mathrm{coh}}}^2}) \exp(-\frac{\Delta t}{2 \tau_{\mathrm{LT}}})  \notag \\
 \label{eq:tau_coh} & =g^1(\Delta t =0) \exp(- \frac{\pi}{2} \frac{\Delta t^2}{\tau_{{\mathrm{coh}}}^2} - \frac{\Delta t}{2 \tau_{\mathrm{LT}}}).\end{aligned}$$

By comparing the experimental data with the corrected model obtained in Eq. \[eq:tau\_coh\] we can estimate a corrected coherence time of $\tau_{\mathrm{coh}} = 14$ ps. If we further apply corrections from the spectrally dependent phase shift, as discussed in the previous section, a maximum coherence time of $\tau_{\mathrm{coh}} = 24$ ps was estimated.

[^1]: aluminum is a common donor in ZnO layers and as a component of the used sapphire substrate and Bragg mirror layers it can easily diffuse into the ZnO cavity during the annealing process at high temperatures of $T~\approx~900~^\circ$

[^2]: The non-equilibrium parameter $\alpha$ can be extracted by the linewidth divided by the derivative of the condensate blueshift w.r.t. excitation power, $\alpha = \hbar \gamma_{{\mathrm{c}}} / (d \Delta E / d(P/{P_{{\mathrm{th}}}}))$. Within the error-bounds we chose $\alpha = 7$ in order to reproduce the experimental data.

[^3]: The shoulder in the experimental data at $k = 3,\dots,4 \ \mu {\mathrm{m}}^{-1}$ is visible for all excitation powers and its shape and magnitude is independent on the power. Therefore we can conclude that this shoulder is caused by residual effect and does not belong to the condensate emission. Different options are possible, e.g. the emission of uncondensed polaritons occupying the LP branch or artifacts from the experimental setup, e.g. transmission fluctuations from the microscope objective.
