# Mapping polariton Bose–Einstein condensate onto vibrational degrees of
freedom

Vladislav Yu. Shishkov1,2 Evgeny S. Andrianov1,2 Sergei Tretiak3,4 K. Birgitta
Whaley5,6 Anton V. Zasedatelev7<EMAIL_ADDRESS>1Dukhov
Research Institute of Automatics (VNIIA), 22 Sushchevskaya, Moscow 127055,
Russia  2Moscow Institute of Physics and Technology, 9 Institutskiy pereulok,
Dolgoprudny 141700, Moscow region, Russia  3Center for Integrated
Nanotechnologies, Los Alamos National Laboratory, Los Alamos, NM, 87545 USA
4Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, 87545
USA  5Department of Chemistry, University of California, Berkeley, California
94720, USA  6Berkeley Center for Quantum Information and Computation,
Berkeley, California 94720, USA  7Vienna Center for Quantum Science and
Technology (VCQ), Faculty of Physics, University of Vienna, Boltzmanngasse 5,
1090 Vienna, Austria

###### Abstract

We demonstrate a macro-coherent regime in molecular exciton-polariton systems,
where nonequilibrium polariton Bose–Einstein condensation coexists with
macroscopically occupied vibrational states. Strong vibronic coupling in
molecules induces an effective optomechanical interaction between cavity
polaritons and vibrational degrees of freedom of molecules, leading to
vibrational amplification in a resonant blue-detuned configuration. This
interaction maps out properties of the condensate to vibrational states,
offering a novel approach to achieve vibrational condensation with potential
applications in cavity-controlled chemistry, nonlinear and quantum optics.

††preprint: APS/123-QED

The formation of new eigenstates - cavity polaritons - through strong light-
matter interaction, grants remarkable control over the physical and chemical
properties of molecular materials. This opens new ways for the modification of
chemical reactions [1, 2, 3, 4], long-range energy transfer [5, 6, 7],
steering of singlet/triplet dynamics [8, 9, 10, 11], and enhancement of
nonlinear optical response [12, 13]. Macroscopic quantum phenomena, such as
nonequilibrium polariton Bose–Einstein condensation (BEC) [14, 15] and
superfluidity [16] in various organic materials at room temperature, have
established light-matter BECs as a versatile platform suitable for diverse
applications [17, 18, 19, 20]. Over the past decade, progress in exciton-
polariton BEC has given rise to novel architectures in molecular
optoelectronics, which include energy-efficient tunable coherent light sources
not reliant on population inversion [21, 22, 23], ultra-fast all-optical
transistors and logic gates [12, 24, 25], extreme photon nonlinearities for
emerging quantum technologies [13], and room-temperature topological states
immune to disorder [26, 27].

While exciton-polariton BECs are relatively well explored, vibrational
polariton condensation remains experimentally elusive. If realized, it will
bring molecular vibrations to a macroscopic quantum state, where all
vibrations share the same quantum properties, become indistinguishable and
behave as of macroscopic coherent matter wave demonstrating long range order.
However, achieving this state is challenging. Factors such as fast vibrational
relaxation, relatively high thermal fluctuations at room temperatures, the
absence of appropriate experimental methods, and insufficient quality of
cavity structures stand as current obstacles on a way towards vibrational BEC.
Amidst these challenges, the significance of vibrational condensates has
grown, particularly owing to the recent progress of cavity-controlled
chemistry. They offer an elegant solution to the dark state problem [28],
preventing undesirable decay into the dark state manifold. Recently Pannir-
Sivajothi and colleagues have harnessed the quantum statistics of BECs to
exploit benefits of light-matter condensates in driving chemical reactions,
thereby reporting distinct energetic and entropic advantages in electron
transfer [29].

Typically Raman-active molecular vibrations do not directly interact with the
cavity mode in exciton-polariton systems. However, they significantly
influence polariton dynamics through vibronic coupling with the excitonic
component. The strong vibronic coupling in molecular systems leads to distinct
replicas in absorption and emission spectra, as well as large Stokes shifts
[30, 31]. In polariton systems, intense high-energy vibrational resonances
enable efficient single-step relaxation [32, 33], reducing the threshold for
polariton condensation [12, 13]. Recent experimental reports showcase
unprecedented high nonlinearity, driven by bosonic stimulation through the
coupling to intense molecular vibrations [13]. However, the quantum state of
the vibrational degrees of freedom in exciton-polariton systems largely remain
an uncharted territory. In this study, we investigate vibrational states in
exciton-polariton systems that have strong vibronic coupling. We identify a
new regime, characterized by macroscopically occupied vibrational states that
coexist with a well-controlled exciton-polariton BEC. Our research introduces
novel approaches to achieve vibrational condensation and could serve as a
guideline for future experiments, taking advantage of strong vibronic coupling
as an inherent mechanism and general feature of molecular systems.

Figure 1: Schematic of the setup to generate macroscopic vibrational states
within exciton-polariton Bose–Einstein condensation (BEC) in a strong vibronic
regime. A coherent laser (Pump) in the blue-detuned configuration resonantly
excites bright excitonic states (Exc), causing vibrational amplification
through optomechanical coupling between molecular vibrations (Vib) and cavity
exciton-polariton BEC satisfying the phase-matching condition: $\hbar\bf
k_{\rm Pump||}=\hbar\bf k_{\rm BEC}+\hbar\bf k_{\rm Vib}$.

Here we consider a nonequilibrium microscopic model that describes a large
ensemble of organic molecules, each hosting a single exciton (${\hat{H}_{\rm
Exc}=\sum_{j=1}^{N_{\rm
mol}}{\hbar{\omega_{{\rm{exc}}}}\hat{\sigma}_{{\rm{Exc}}{j}}^{\dagger}{{\hat{\sigma}}_{{\rm{Exc}}{j}}}}}$,
where $N_{\rm mol}$ is the total number of molecules, $\omega_{\rm exc}$ is
the eigenfrequency of the excitons, $\hat{\sigma}_{{\rm{Exc}}{j}}^{\dagger}$
($\hat{\sigma}_{{\rm{Exc}}{j}}$) is the creation (annihilation) operator for
an exciton of a molecule located at the point ${{\bf r}_{j}}$) coupled
strongly to an optical cavity (${\hat{H}_{\rm Cav}=\sum_{\bf
k}\hbar{\omega_{{\rm Cav}|{\bf k}}}\hat{a}_{{\rm Cav}|{\bf
k}}^{\dagger}{{\hat{a}}_{{\rm Cav}|{\bf k}}}}$, where $\hat{a}_{{\rm Cav}|{\bf
k}}^{\dagger}$ ($\hat{a}_{{\rm Cav}|{\bf k}}$) is the creation (annihilation)
operator for a cavity photon with in-plane wavevector ${\bf k}$ and frequency
${\omega_{{\rm Cav}|{\bf k}}}$) and interacting with molecular vibrations
(${\hat{H}_{\rm Vib}=\sum_{j=1}^{N_{\rm
mol}}{\hbar{\omega_{\rm{Vib}}}\hat{b}_{{\rm Vib}j}^{\dagger}{{\hat{b}}_{{\rm
Vib}j}}}}$, where $\hat{b}_{{\rm Vib}j}^{\dagger}$ ($\hat{b}_{{\rm Vib}j}$) is
the creation (annihilation) operator of a molecular vibration with the
eigenfrequency $\omega_{\rm Vib}$ for the corresponding molecule). The
electronic excitation, initially confined to a single molecule, becomes
delocalised through interaction with the cavity field: $\hat{H}_{\rm Exc-
Cav}=\sum_{j=1}^{N_{\rm mol}}\sum_{{\bf k}}{\hbar\Omega_{j{\bf
k}}\left(\hat{\sigma}_{{\rm Exc}{j}}^{\dagger}{{\hat{a}}_{{\rm Cav}|{\bf
k}}}{e^{i{\bf k}{\bf r}_{j}}}+h.c.\right)}$, where $\Omega_{j{\bf k}}$ is Rabi
frequency of the $j$-th molecule and cavity mode with wavevector $\bf k$ [34].
This interaction implicitly involves vibronic coupling via the matter
component, represented by the term in the Hamiltonian: $\hat{H}_{\rm Exc-
Vib}=\sum_{j=1}^{N_{\rm mol}}{\hbar\Lambda\omega_{\rm Vib}\hat{\sigma}_{{\rm
Exc}j}^{\dagger}{{\hat{\sigma}}_{{\rm Exc}{j}}}\left({{{\hat{b}}_{{\rm
Vib}j}}+\hat{b}_{{\rm Vib}j}^{\dagger}}\right)}$, where $\Lambda$ is the
interaction constant between excitons and molecular vibrations being a square
root of the Huang–Rhys (HR) factor [30]. It quantifies interaction between the
electronic structure of a molecule and the displacement of its nucleus, which
subsequently extend significant influence on the polaritonic degrees of
freedom. The full system Hamiltonian thus includes molecular excitons along
with vibrational eigenenergies, cavity modes, and the aforementioned coupling
terms between subsystems (see details in Supplementary Materials, Section I).

$\hat{H}=\hat{H}_{\rm Exc}+\hat{H}_{\rm Cav}+\hat{H}_{\rm Vib}+\hat{H}_{\rm
Exc-Cav}+\hat{H}_{\rm Exc-Vib}.$ (1)

Considering the strong interaction between excitons, the cavity, and molecular
vibrations, we necessitate a transformation to the momentum representation of
the dressed states (see Supplemental Materials, Section I, for details). The
transformation unfolds in several stages. Firstly, we introduce vibrationally
dressed excitons and dressed molecular vibrations [35]. Subsequently, we
account for the light-matter interaction of the dressed excitons with the
cavity, yielding lower and upper polariton states [36]. Excitons accessible by
dipole coupling from the ground state, termed “bright excitons”, contrast with
the remaining “dark excitons”. Both bright excitons and polaritons are phase-
coherent, many-body delocalized states, which have a well-defined in-plane
momentum $\hbar\bf k$, matching the corresponding eigenstates of the cavity
and in-plane component of the incident wavevector of a pump beam $\bf k_{\rm
Pump||}$. Dark excitons, lacking well-defined momentum, represent a manifold
of localized states. Analogously, we separate the dressed molecular vibrations
into the bright and dark types. The bright vibrations are coherent and
delocalized states. Similar to the bright excitons they have a well-defined
momentum $\hbar\bf k_{\rm Vib}$ [37], unlike the localized dark vibrational
states.

Moving forward from a static eigenstate problem to the dynamic description of
the system, we deepen our analysis using the Lindblad approach, which accounts
for interactions with the environment [38, 34, 39]. We establish a theoretical
framework to calculate average values of observables, which includes
polariton, exciton, and vibration occupation numbers. Given that interactions
with the environment inevitably result in decoherence and energy relaxation,
it is very important to analyse potential relaxation and thermalization
processes in the system. Thus, intermolecular interaction and non-radiative
decay are the main mechanisms for excitons [40, 41]. In contrast, for
polaritonic states, cavity decay often dominates the energy and phase loss in
practical devices [17]. The coherence of vibrational states is typically
lifetime-limited by a fast inter- and intramolecular vibrational energy
redistribution (IVR) [42]. Thermalization is another essential process for
polariton BEC [43]. In organic polariton systems, thermalization can occur due
to intermolecular energy relaxation and via nonlinear interaction with low-
frequency vibrations. We describe all these relaxation and thermalization
processes by Lindblad superoperators $\hat{L}$ acting on the general density
matrix operator $\hat{\rho}$ (see details in Supplementary Materials, Section
II).

We consider coherent laser excitation of the bright excitonic states with in-
plane momentum $\hbar\bf k_{\rm ex}=\hbar\bf k_{\rm Pump||}$, adhering to the
resonant blue-detuned configuration $\hbar\omega_{\rm Pump}=\hbar\omega_{\rm
Exc}=\hbar\omega_{\rm Pol|\bf k=\bf 0}+\hbar\omega_{\rm Vib}$ as illustrated
in Figure 1. Here, $\hbar\omega_{\rm Pump}$ and $\hbar\omega_{{\rm Pol}|\bf
k=\bf 0}$ denote photon energy of the pump and the ground polariton state
energy. Next, we derive a master equation and develop a mean-field theory for
the system. The time evolution of average occupation numbers within the system
is described by the equations below, and further elaborated in Supplementary
Materials, Section III. This includes polaritons $n_{{\rm Pol}|{\bf k}}$ and
bright vibrational states $n_{{\rm Vib}|{\bf k}}$, associated with momenta
$\hbar\bf k$ and $\hbar\bf k_{\rm Vib}$ respectively, along with the
occupation for the resonantly pumped bright exciton state $n_{{\rm Exc}|{\bf
k}_{\rm ex}}$, all dark excitons $n_{\rm Exc_{D}}$, and vibrations $n_{\rm
Vib_{D}}$.

$\frac{dn_{\rm Exc|{\bf k}_{\rm ex}}}{dt}=-\gamma_{\rm Exc}(n_{\rm Exc|{\bf
k}_{\rm ex}}-\varkappa_{\rm Pump})+\gamma_{\rm Exc}^{\rm B-D}(n_{\rm
Exc_{D}}-\\\ n_{{\rm Exc}|{\bf k}_{\rm ex}})-\sum_{\bf k}\frac{G_{{\bf
k}}}{N_{\rm mol}}\left[n_{{\rm Exc}|{\bf k}_{\rm ex}}\left(n_{{\rm Pol}|{\bf
k}}+1\right)+\right.\\\ \left.n_{{\rm Vib}|{\bf k_{\rm ex}}-{\bf
k}}\left(n_{{\rm Exc}|{\bf k}_{\rm ex}}-n_{{\rm Exc}|{\bf k}_{\rm
ex}}\right)\right]$ (2) $\frac{dn_{\rm Exc_{D}}}{dt}=-\gamma_{\rm Exc}n_{\rm
Exc_{D}}+{\gamma_{\rm Exc}^{\rm B-D}\over N_{\rm mol}}(n_{{\rm Exc}|{\bf
k}_{\rm ex}}-n_{\rm Exc_{D}})-\\\ \sum_{\bf k}\frac{G_{\bf k}}{N_{\rm
mol}}\left[n_{\rm Exc_{D}}\left(n_{{\rm Pol}|{\bf k}}+1\right)+n_{\rm
Vib_{D}}\left(n_{\rm Exc_{D}}-n_{{\rm Pol}|{\bf k}}\right)\right]$ (3)
$\frac{dn_{{\rm Pol}|{\bf k}}}{dt}=-\gamma_{{\rm Pol}|{\bf k}}n_{{\rm
Pol}|{\bf k}}+G_{\bf k}\left[n_{\rm Exc_{D}}\left(n_{{\rm Pol}|{\bf
k}}+1\right)+\right.\\\ \left.n_{\rm Vib_{D}}\left(n_{\rm Exc_{D}}-n_{{\rm
Pol}|{\bf k}}\right)\right]+{G_{\bf k}\over N_{\rm mol}}\left[n_{\rm
Exc|\bf{k_{ex}}}\left(n_{{\rm Pol}|{\bf k}}+\right.\right.\\\
\left.\left.1\right)+n_{{\rm Vib}|{\bf k_{\rm ex}}-{\bf k}}\left(n_{{\rm
Exc}|{\bf k}_{\rm ex}}-n_{{\rm Pol}|{\bf k}}\right)\right]+\sum_{\bf
k^{\prime}}\left\\{\gamma_{\rm therm}^{\bf kk^{\prime}}\right.\\\
\left.\left(n_{{\rm Pol}|{\bf k}}+1\right){n_{{\rm Pol}|{\bf
k^{\prime}}}}-\gamma_{{\rm{therm}}}^{{\bf{k^{\prime}k}}}\left(n_{{\rm
Pol}|{\bf k^{\prime}}}+1\right){n_{{\rm Pol}|{\bf k}}}\right\\}$ (4)
$\frac{dn_{\rm Vib_{D}}}{dt}=-\gamma_{\rm Vib}\left(n_{\rm Vib_{D}}-n_{\rm
Vib}^{\rm th}\right)+{\gamma_{\rm Vib}^{\rm B-D}\over N_{\rm mol}}\sum_{\bf
k}\left(\right.\\\ \left.n_{{\rm Vib}|{\bf k}_{\rm ex}-{\bf k}}-n_{\rm
Vib_{D}}\right)+\sum_{\bf k}\frac{G_{\bf k}}{N_{\rm mol}}\left[n_{\rm
Exc_{D}}\left(n_{{\rm Pol}|{\bf k}}+1\right)+\right.\\\ \left.n_{\rm
Vib_{D}}\left(n_{\rm Exc_{D}}-n_{{\rm Pol}|{\bf k}}\right)\right]$ (5)
$\frac{dn_{{\rm Vib}|{\bf k}_{\rm ex}-{\bf k}}}{dt}=-\gamma_{\rm
Vib}\left(n_{{\rm Vib}|{\bf k}_{\rm ex}-{\bf k}}-n_{{\rm vib}}^{\rm
th}\right)+\gamma_{\rm Vib}^{\rm B-D}\\\ \left(n_{\rm Vib_{D}}-n_{{\rm
Vib}|{\bf k}_{\rm ex}-{\bf k}}\right)+\frac{G_{\bf k}}{N_{\rm
mol}}\left[n_{\rm Exc|{\bf k}_{\rm ex}}\left(n_{{\rm Pol}|{\bf
k}}+1\right)+\right.\\\ \left.n_{{\rm Vib}|{\bf k_{\rm ex}}-{\bf
k}}\left(n_{\rm Exc|{\bf k}_{\rm ex}}-n_{{\rm Pol}|{\bf k}}\right)\right]$ (6)

Here, $n_{{\rm Exc}|{\bf k}_{\rm ex}}$ refers to the number of dressed bright
excitons pumped directly by the resonant excitation, $n_{\bf k}$ is the
average number of lower polaritons in the state with $\hbar\bf k$ in-plane
momentum, $n_{\rm Exc_{D}}$ represents the number of dressed dark excitons,
$n_{\rm Vib}$ and $n_{\rm Vib_{D}}$ denote the number of bright and dark
molecular vibrations, respectively. The term $n_{\rm Vib}^{\rm
th}=(e^{\hbar\omega_{\rm Vib}/T}-1)^{-1}$ is the vibrational population in the
thermal equilibrium with an environment at temperature $T$. The parameters
$\gamma_{\rm Exc}$, $\gamma_{{\rm Pol}|{\bf k}}$ and $\gamma_{\rm Vib}$ are
the energy dissipation rates for excitons, lower polaritons with wavevector
$\bf k$ and molecular vibrations respectively. Thermalization rates
$\gamma_{{\rm{therm}}}^{\bf kk^{\prime}}$ and $\gamma_{{\rm{therm}}}^{\bf
k^{\prime}k}$ represent downward and upward energy relaxation in the polariton
subsystem, respectively, linked by the Kubo–Martin–Schwinger relation [44, 45,
39, 46]. Transition rates between bright and dark states $\gamma_{\rm
Exc}^{\rm B-D}$ and $\gamma_{\rm Vib}^{\rm B-D}$ correspond to the relaxation
of the bright excitonic and vibrational states to the manifold of localized
dark states due to IVR mechanisms leading to the dephasing $\Gamma_{\rm
Exc}=\gamma_{\rm Exc}^{\rm B-D}$ and $\Gamma_{\rm Vib}=\gamma_{\rm Vib}^{\rm
B-D}$, respectively. The resonant pumping rate $\varkappa_{\rm Pump}$, being
devided by total number of molecules $N_{\rm mol}$ has a clear physical
meaning of the stationary excitation number per molecule.

Figure 2: Phase diagrams for polariton and vibrational states demonstrating
the regime of coexisting polariton BEC and macroscopically occupied
vibrational states. (a),(b) - Polariton and vibrational occupation numbers as
the function of Huang-Rhys factor $\Lambda^{2}$ and the resonant pumping rate
$\varkappa_{\rm Pump}$ normalized to the number of molecules $N_{\rm mol}$.
Vibrational decoherence rate is assumed to be lifetime-limited $\Gamma_{\rm
Vib}=\gamma_{\rm Vib}/2=5\times 10^{-3}\omega_{\rm Vib}$. (c),(d) - Polariton
and vibrational occupation numbers as the function of vibrational decoherence
rate $\Gamma_{\rm Vib}$ and the resonant pumping rate $\varkappa_{\rm Pump}$.
Here we fix Huang-Rhys factor $\Lambda^{2}=1$. Black dashed line shows the
analytical value of the threshold for polariton and vibrational
macroscopically occupied states given by Eq. (8). The parameters here are
$N_{\rm mol}=10^{8}$, $\Omega_{R}=0.05$ $\rm eV$, $\gamma_{\rm Exc}=10^{-5}$
$\rm{eV}$, $\gamma_{\rm Exc}^{\rm B-D}=\Gamma_{\rm Exc}=10^{-2}$ $\rm{eV}$,
$\gamma_{\rm Vib}^{\rm B-D}=\Gamma_{\rm Vib}$, $\gamma_{{\rm{therm}}}^{\bf
kk^{\prime}}=10^{-5}$ $\rm eV$ for $|{\bf k}|<|{\bf k^{\prime}}|$ and $T=290$
$\rm{K}$, $\gamma_{{\rm Pol}|{\bf
k}}\approx\gamma_{\rm{Cav}|{\bf{k}}}=2.5\cdot 10^{-3}$ $\rm{eV}$,
$\omega_{{\rm Cav}|{\bf k}}=\omega_{{\rm Cav}|{\bf k}={\bf 0}}+\alpha_{\rm
Cav}{\bf k}^{2}$ with $\alpha_{\rm Cav}=2\cdot 10^{-3}$ ${\rm eV}\mu{\rm
m^{2}}$ and $S=500$ $\mu{\rm m^{2}}$. These parameters are consistent with the
recent experiments [12, 13].

The constant $G_{\bf k}$ plays a central role in polariton condensation and
generation of macroscopic vibrational states:

$G_{\bf k}=\frac{\Lambda^{2}\Omega_{R}^{2}\Gamma_{\rm Exc}}{(\omega_{\rm
Exc}-\omega_{{\rm Pol}|{\bf k}}-\omega_{{\rm Vib}})^{2}+(\Gamma_{\rm
Exc}/2)^{2}}$ (7)

where $\omega_{\rm Exc}=\omega_{\rm exc}-\Lambda^{2}\omega_{\rm Vib}$ and
$\Omega_{R}^{2}=\sum_{j=1}^{N_{\rm mol}}|\Omega_{j{\bf k}}|^{2}$ represents
the total Rabi frequency. Serving as the effective polariton-vibration
coupling constant, $G_{\bf k}$ characterizes the interaction strength between
exciton-polariton and vibrational states according to the term $\hbar G_{\bf
k}\hat{a}_{{\rm Pol}|{\bf k}}^{\dagger}{\hat{a}_{{\rm Pol}|{\bf
k}}\left({{{\hat{b}}_{{\rm Vib}|{\bf k}}}+b_{{{\rm Vib}|{\bf
k}}}^{\dagger}}\right)}$. Evidently, this coupling corresponds to a
conventional cavity optomechanical Hamiltonian [47], except that here, we deal
with a hybrid light-matter state (exciton-polariton) rather than a bare cavity
mode. Within this framework, molecular vibrational modes take on the function
of a mechanical oscillator. The optomechanical interaction steers population
dynamics across the cavity and mechanical degrees of freedom and depends on
the optomechanical back action rate, denoted as $\sim G_{\bf k}^{2}\times N$
[47]. In our case, the blue-detuned external laser drive $\hbar\omega_{\rm
Pump}=\hbar\omega_{{\rm Pol}|{\bf k}={\bf 0}}+\hbar\omega_{\rm Vib}$ imposes a
negative back action rate (anti-damping), resulting in strong vibrational
amplification under certain conditions, as illustrated in Figure 1.

Equations (2)–(6) provide a reduced microscopic model able to simulate
dynamics of electronic and vibrational states of molecules with strong
vibronic coupling in realistic cavities, given the right choice of
experimental parameters. We numerically solve Eqs. (2)–(6) and find a steady-
state population of polaritons at the ground state (${\bf k}={\bf 0}$) as well
as the population of the bright vibrational state mediating polariton
condensation that fulfils phase matching condition $\hbar\bf k_{\rm
Vib}=\hbar\bf k_{\rm Pump||}-\hbar\bf k_{\rm BEC}$ (refer to the Supplementary
Materials, Section IV for a discrete version of Eqs. (2)–(6)). Figure 2
presents density plots of the polariton and vibrational occupation numbers as
functions of the pumping rate $\varkappa_{\rm Pump}/N_{\rm mol}$, Huang-Rhys
factor $\sim\Lambda^{2}$, and vibrational decay rate $\gamma_{\rm Vib}$. An
analytic approximation to the steady-state solution reveals the threshold
pumping rate as a function of polariton decay rate and optomechanical coupling

$\varkappa_{\rm Pump}^{\rm thresh}=N_{\rm mol}\left.{\gamma_{{\rm Pol}|{\bf
k}}\over G_{{\bf k}}}\right|_{{\bf k}={\bf 0}},$ (8)

as detailed in Supplementary Materials, Section V. Notably, the threshold
value (8) $\varkappa_{\rm Pump}^{\rm thresh}\propto\Lambda^{-2}$ does not
depend on the vibrational decay rate $\gamma_{\rm Vib}$, which makes our
proposal robust against fast IVR mechanisms. Figures 2(a,c) show an average
occupation number in the polariton subsystem as a function of Huang-Rhys
factor and vibrational decay.

Figure 3: Average occupation number of the bright vibrational state (brown)
with the wavevector ${\bf k}_{\rm ex}$ and at the ground polariton state
(green) (a). (b) - Energy distributions for polariton states below (blue) and
above (red) the condensation threshold, where $\Delta\omega_{\bf
k}=\omega_{{\rm Pol}|{\bf k}}-\omega_{{\rm Pol}|{\bf k}={\bf 0}}$ . (c,d) -
Momentum distributions for vibrational and polariton states below (blue) and
above (red) the condensation threshold. The parameters are the same as for
Figure 2. Fitting curves in (b) and (d) represent the effective temperature
determined by the Boltzmann distribution in energy and momentum spaces,
respectively.

The formation of polariton BEC, in general, requires two conditions: 1 - the
rate of polariton thermalization overcomes the energy dissipation ($\sum_{{\bf
k}}\gamma_{\rm therm}^{{\bf k}{\bf 0}}n_{{\rm Pol}|{\bf k}}\gtrsim\gamma_{{\rm
Pol}|{\bf k}}$), and 2 - the total number of lower polariton surpasses the
critical number $\sum_{\bf k}n_{{\rm Pol}|{\bf k}}|_{\rm critical}\approx\nu
k_{B}T$ [48, 49]. In our system, the first condition is met well below the
threshold (8) while the second one is achieved at the threshold ${\sum_{\bf
k}n_{{\rm Pol}|{\bf k}}|_{\rm thresh}\sim\sum_{\bf k}n_{{\rm Pol}|{\bf
k}}|_{\rm critical}\sim 10^{2}-10^{3}}$ (see Supplementary Materials, Section
V). Followed by BEC formation of polaritons with in-plane wavevector $\bf
k=\bf 0$, the system exhibits macroscopic occupation of the vibrational
states. The momentum conservation law sets the in-plane wavevector of the
maroscopically occupated vabrational states to be ${\bf k}_{\rm Vib}={\bf
k}_{\rm ex}$. Figure 2b shows vibrational amplification towards the
macroscopically occupied state above a certain pumping rate, consistent with
the polariton BEC threshold in Eq. (8). The matching threshold values for
polariton and vibrational subsystems, unaffected by vibrational decay rate,
highlights the fact that macroscopic occupation of the vibrational state is a
result of polariton condensation. This regime of coexisting macroscopically
occupied exciton-polariton and vibrational states is a distinct feature of
exciton-polariton systems with strong vibronic interaction in the resonant
blue-detuned configuration. Although the decay of molecular vibrations does
not explicitly influence the threshold, it does have an impact on the average
vibration occupation number, as shown in Figure 2(d) and is also evidenced by
the relationship $n_{{\rm Vib}|{\bf k}_{\rm ex}}\propto\gamma_{\rm
Vib}^{\rm-1}$.

The polariton optomechanical interaction in our system allows the state of the
BEC to be mapped onto vibrational degrees of freedom. We consider two regimes
when being below (blue) and above (red) the condensation threshold. Figure 3a
shows the average occupation numbers of the bright vibrational state with the
wavevector ${\bf k}_{\rm ex}$ and the ground polariton state at the resonant
condition. Owing to the strong light-matter interaction, exciton-polaritons
show a parabolic dispersion relation near the ground state, unlike
dispersionless uncoupled vibrational states. Hence, further analysis of energy
distribution is restricted to the polariton subsystem. Figure 3b illustrates
energy distributions below and above the condensation threshold, with
polaritons exhibiting thermalization at an effective temperature $T_{\rm
eff}\approx 140K$. Above the threshold, the system undergoes a Bose–Einstein
distribution with the thermalized tail at $T_{\rm BEC}=0.6T_{\rm eff}$, a
distinctive cooling effect that comes from the nonequilibrium nature of
polariton BEC, recently demonstrated in Ref. [48] In momentum space, Figure
3c, vibrational states lock to the polariton distribution due to resonant
phase-matching condition. When the pumping rate exceeds the condensation
threshold, it collapses to the bright state with a well-defined in-plane
momentum $\hbar{\bf k}_{\rm ex}$ exhibiting thermalised distribution for
states at high in-plane momenta. This resembles the momentum distribution of
the polariton BEC depicted in Figure 3d, effectively mapping the vibrational
subsystem into a macroscopically occupied BEC-like state.

Next, we discuss evolution of the long-time vibronic coherences enabled by the
blue-detuned configuration in Figure 1. Our proposal involves two coherent
drives acting on the vibrational mode above condensation threshold: the laser
pump and polariton BEC. Similar to coherent Raman scattering [50] both waves
map out their coherence on the vibrational state. However, in contrast to
Raman-based methods, which imprint coherence only on vibrational degrees of
freedom, our approach enables a concurrent vibrational and exciton-polariton
macro-coherent regime. This necessitates resonant electronic excitation, for
which two-dimensional electronic spectroscopy serves as an accessible tool to
generate macroscopically coherent electron-vibrational (vibronic) states in
organic systems at room temperature [51]. Yet, ultra-broadband excitation
imposes certain constraints on coherence time, causing vibronic coherence to
dissipate within a hundred femtoseconds due to inherent dephasing of
electronic states [52, 53]. Our approach eliminates these effects by
leveraging exciton-polariton states, protected from natural decoherence by the
cavity. In fact, in our case the coherence time is limited by nonlinear
effects occurring at the exciton-polariton BEC [54, 55]. This enables an
improvement of at least three orders of magnitude in coherence time at room
temperature under narrow-band quasi-steady-state pumping conditions [54, 23].
Moreover, the coherence time can be further manipulated by engineering the
nonlinear potential landscape for the condensate [56, 57], potentially
extending it to the nanosecond level in future experiments, such as with novel
multi-component polariton systems [58, 59].

Figure 4: Probing macroscopic vibrational states by non-resonant Raman
spectroscopy of polariton BEC. The graph illustrates the angle-resolved and
pump-dependent distribution of the intensity ratio between anti-Stokes and
Stokes components of the Raman scattering, with parameters $\hbar\omega_{\rm
Probe}=2~{}{\rm eV}$, $\theta_{\rm Probe}=-30^{\rm o}$ and the rest parameters
are the same as for Figure 2. The inset shows proposed experimental
arrangement.

To validate the presence of macroscopically occupied vibrational states, we
propose using angle-resolved non-resonant Raman spectroscopy of polariton BEC.
The respective experimental arrangement is shown in Figure 4. We calculate the
intensity of the Stokes and anti-Stokes components of Raman scattering for a
non-resonant probe beam with an energy $\hbar\omega_{\rm Probe}$ and in-plane
momentum $\hbar\bf k_{\rm Probe}$. In the practical setups, the probe beam
must be synchronized with the onset of polariton condensation. Figure 4
demonstrates the ratio between the anti-Stokes and Stokes components as a
function of both pumping rate and detection angle. A coherent anti-Stokes
Raman scattering (CARS) signal emerges at a specific “magic angle” (MA),
denoted as $\theta_{\rm MA}$ (see Eq. (9)), due to resonant vibrational
amplification.

$\sin\theta_{\rm MA}={\omega_{\rm Pump}\sin\theta_{\rm Pump}+\omega_{\rm
Probe}\sin\theta_{\rm Probe}\over\omega_{\rm Probe}+\omega_{\rm Vib}}$ (9)

The magic angle provides a direct measurement of the momentum of the
vibrational mode $\hbar\bf k_{\rm ex}$ corresponding to the resonant condition
(see Supplementary Materials, Section VI for details). The formation of
macroscopically occupied vibrational states is manifested by the CARS signal,
which asymptotically approaches intensity of the Stokes component above
condensation threshold, as depicted in Figure 4.

In conclusion, we have introduced a new way of generating macroscopic
vibrational states through exciton-polariton condensation in the strong
vibronic regime. We have demonstrated that an effective strong optomechanical
coupling between vibrational and polariton degrees of freedom can be achieved
characterized by interaction strength $G_{\bf k}$. This coupling results in
parametric amplification of the resonant vibrational states above condensation
threshold in close analogy to the cavity optomechanics in the blue-detuned
configuration [47]. Our findings predict a unique macro-coherent regime in
strongly-coupled organic cavities. This has conceptual implications for
understanding coherence in molecular systems and its compatibility with
various chemical and physical processes [60]. We identified conditions where
vibronic light-matter BEC results in both macroscopic vibrational and exciton-
polariton states, aligning with cutting-edge experiments. Recent studies
further validate the vibronic mechanism in polariton condensation, revealing
significant energy advantage over non-resonant configurations [12, 13]. This
implies the likely existence of macroscopic vibrational states in these
experiments, awaiting verification. We further suggest a practical setup using
a non-resonant Raman probe to identify these states, if they are present.
Expanding upon the proposed scheme to include non-resonant vibrational control
methods is straightforward. Together with the recently developed concept of
the resonant seeded polariton condensation at room temperature, this will open
up new avenues for manipulating light-matter states across a broad spectral
range for nonlinear and quantum optics with molecular systems.

V.Yu.Sh. thanks Foundation for the Advancement of Theoretical Physics and
Mathematics ”Basis” for financial support. A.V.Z. acknowledges support from
the European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No 101030987 (LOREN). This work was
performed, in part, at the Center for Integrated Nanotechnologies, an Office
of Science User Facility at Los Alamos National Laboratory operated for the
U.S. Department of Energy (DOE) Office of Science.

## References

  * Hutchison  _et al._ [2012] J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, Angewandte Chemie International Edition 51, 1592 (2012).
  * Ahn  _et al._ [2023] W. Ahn, J. F. Triana, F. Recabal, F. Herrera, and B. S. Simpkins, Science 380, 1165 (2023).
  * Yuen-Zhou and Menon [2019] J. Yuen-Zhou and V. M. Menon, Proceedings of the National Academy of Sciences 116, 5214 (2019).
  * Garcia-Vidal  _et al._ [2021] F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Science 373, eabd0336 (2021).
  * Sáez-Blázquez  _et al._ [2018] R. Sáez-Blázquez, J. Feist, A. I. Fernández-Domínguez, and F. García-Vidal, Physical Review B 97, 241407 (2018).
  * Xiang  _et al._ [2020] B. Xiang, R. F. Ribeiro, M. Du, L. Chen, Z. Yang, J. Wang, J. Yuen-Zhou, and W. Xiong, Science 368, 665 (2020).
  * Georgiou  _et al._ [2021] K. Georgiou, R. Jayaprakash, A. Othonos, and D. G. Lidzey, Angewandte Chemie 133, 16797 (2021).
  * Martínez-Martínez  _et al._ [2018] L. A. Martínez-Martínez, M. Du, R. F. Ribeiro, S. Kéna-Cohen, and J. Yuen-Zhou, The Journal of Physical Chemistry Letters 9, 1951 (2018).
  * Martínez-Martínez  _et al._ [2019] L. A. Martínez-Martínez, E. Eizner, S. Kéna-Cohen, and J. Yuen-Zhou, The Journal of Chemical Physics 151, https://doi.org/10.1063/1.5100192 (2019).
  * Eizner  _et al._ [2019] E. Eizner, L. A. Martínez-Martínez, J. Yuen-Zhou, and S. Kéna-Cohen, Science Advances 5, eaax4482 (2019).
  * Polak  _et al._ [2020] D. Polak, R. Jayaprakash, T. P. Lyons, L. Á. Martínez-Martínez, A. Leventis, K. J. Fallon, H. Coulthard, D. G. Bossanyi, K. Georgiou, A. J. Petty, _et al._ , Chemical science 11, 343 (2020).
  * Zasedatelev  _et al._ [2019] A. V. Zasedatelev, A. V. Baranikov, D. Urbonas, F. Scafirimuto, U. Scherf, T. Stöferle, R. F. Mahrt, and P. G. Lagoudakis, Nature Photonics 13, 378 (2019).
  * Zasedatelev  _et al._ [2021] A. V. Zasedatelev, A. V. Baranikov, D. Sannikov, D. Urbonas, F. Scafirimuto, V. Y. Shishkov, E. S. Andrianov, Y. E. Lozovik, U. Scherf, T. Stöferle, _et al._ , Nature 597, 493 (2021).
  * Plumhof  _et al._ [2014] J. D. Plumhof, T. Stöferle, L. Mai, U. Scherf, and R. F. Mahrt, Nature materials 13, 247 (2014).
  * Cookson  _et al._ [2017] T. Cookson, K. Georgiou, A. Zasedatelev, R. T. Grant, T. Virgili, M. Cavazzini, F. Galeotti, C. Clark, N. G. Berloff, D. G. Lidzey, _et al._ , Advanced Optical Materials 5, 1700203 (2017).
  * Lerario  _et al._ [2017] G. Lerario, A. Fieramosca, F. Barachati, D. Ballarini, K. S. Daskalakis, L. Dominici, M. De Giorgi, S. A. Maier, G. Gigli, S. Kéna-Cohen, _et al._ , Nature Physics 13, 837 (2017).
  * Keeling and Kéna-Cohen [2020] J. Keeling and S. Kéna-Cohen, Annual Review of Physical Chemistry 71, 435 (2020).
  * Jiang  _et al._ [2022] Z. Jiang, A. Ren, Y. Yan, J. Yao, and Y. S. Zhao, Advanced Materials 34, 2106095 (2022).
  * Kavokin  _et al._ [2022] A. Kavokin, T. C. Liew, C. Schneider, P. G. Lagoudakis, S. Klembt, and S. Hoefling, Nature Reviews Physics 4, 435 (2022).
  * Luo  _et al._ [2023] S. Luo, H. Zhou, L. Zhang, and Z. Chen, Applied Physics Reviews 10, https://doi.org/10.1063/5.0121316 (2023).
  * Sannikov  _et al._ [2019] D. Sannikov, T. Yagafarov, K. Georgiou, A. Zasedatelev, A. Baranikov, L. Gai, Z. Shen, D. Lidzey, and P. Lagoudakis, Advanced Optical Materials 7, 1900163 (2019).
  * Wei  _et al._ [2021] M. Wei, A. Ruseckas, V. T. Mai, A. Shukla, I. Allison, S.-C. Lo, E. B. Namdas, G. A. Turnbull, and I. D. Samuel, Laser & Photonics Reviews 15, 2100028 (2021).
  * Putintsev  _et al._ [2020] A. Putintsev, A. Zasedatelev, K. E. McGhee, T. Cookson, K. Georgiou, D. Sannikov, D. G. Lidzey, and P. Lagoudakis, Applied Physics Letters 117, https://doi.org/10.1063/5.0019195 (2020).
  * Baranikov  _et al._ [2020] A. V. Baranikov, A. V. Zasedatelev, D. Urbonas, F. Scafirimuto, U. Scherf, T. Stöferle, R. F. Mahrt, and P. G. Lagoudakis, arXiv preprint arXiv:2005.04802 https://doi.org/10.48550/arXiv.2005.04802 (2020).
  * De  _et al._ [2023] J. De, X. Ma, F. Yin, J. Ren, J. Yao, S. Schumacher, Q. Liao, H. Fu, G. Malpuech, and D. Solnyshkov, Journal of the American Chemical Society 145, 1557 (2023).
  * Dusel  _et al._ [2021] M. Dusel, S. Betzold, T. H. Harder, M. Emmerling, J. Beierlein, J. Ohmer, U. Fischer, R. Thomale, C. Schneider, S. Hofling, _et al._ , Nano Letters 21, 6398 (2021).
  * Scafirimuto  _et al._ [2021] F. Scafirimuto, D. Urbonas, M. A. Becker, U. Scherf, R. F. Mahrt, and T. Stöferle, Communications Physics 4, 39 (2021).
  * Vurgaftman  _et al._ [2020] I. Vurgaftman, B. S. Simpkins, A. D. Dunkelberger, and J. C. Owrutsky, The journal of physical chemistry letters 11, 3557 (2020).
  * Pannir-Sivajothi  _et al._ [2022] S. Pannir-Sivajothi, J. A. Campos-Gonzalez-Angulo, L. A. Martínez-Martínez, S. Sinha, and J. Yuen-Zhou, Nature communications 13, 1645 (2022).
  * Atkins and Friedman [2011] P. W. Atkins and R. S. Friedman,  _Molecular quantum mechanics_ (Oxford university press, 2011).
  * Banwell and McCash [2017] C. N. Banwell and E. M. McCash,  _Fundamentals of molecular spectroscopy_ (Indian Edition, 2017).
  * Coles  _et al._ [2011] D. M. Coles, P. Michetti, C. Clark, W. C. Tsoi, A. M. Adawi, J.-S. Kim, and D. G. Lidzey, Advanced Functional Materials 21, 3691 (2011).
  * Somaschi  _et al._ [2011] N. Somaschi, L. Mouchliadis, D. Coles, I. Perakis, D. Lidzey, P. Lagoudakis, and P. Savvidis, Applied Physics Letters 99, https://doi.org/10.1063/1.3645633 (2011).
  * Scully and Zubairy [1997] M. O. Scully and S. Zubairy,  _Quantum optics_ (Cambridge University Press, 1997).
  * [35] _The operators for bare excitons and molecular vibrations, denoted $\hat{\sigma}_{{\rm Exc}j}$ and $\hat{b}_{{\rm Vib}j}$ , are related to the operators for dressed states ${\mathcal{\hat{S}}}_{{\rm Exc}j}$ and ${\mathcal{\hat{B}}}_{{\rm Vib}j}$ through the following relation: $\hat{\sigma}_{{\rm Exc}_{j}}={\mathcal{\hat{S}}}_{{\rm Exc}_{j}}e^{-\Lambda({\mathcal{\hat{B}}}^{\dagger}_{{\rm Vib}j}-{\mathcal{\hat{B}}}_{{\rm Vib}j})}$; $\hat{b}_{{\rm Vib}j}={\mathcal{\hat{B}}}_{{\rm Vib}j}-\Lambda{\mathcal{\hat{S}}}^{\dagger}_{{\rm Exc}_{j}}{\mathcal{\hat{S}}}_{{\rm Exc}_{j}}$ _.
  * [36] _The operators for lower and upper exciton-polariton states are defined as follows: $\hat{s}_{{\rm Low}|{\bf k}}=\hat{a}_{{\rm Cav}|{\bf k}}\cos\varphi_{\bf k}-\hat{c}_{{\rm Exc}|{\bf k}}\sin\varphi_{\bf k}$; $\hat{s}_{{\rm Up}|{\bf k}}=\hat{a}_{{\rm Cav}|{\bf k}}\sin\varphi_{\bf k}+\hat{c}_{{\rm Exc}|{\bf k}}\cos\varphi_{\bf k}$_.
  * [37] _Delocalized collective (bright) excitonic $\hat{c}_{{\rm Exc}|{\bf k}}$ and vibrational $\hat{c}_{{\rm Vib}|{\bf k}}$ states are defined as follows: $\hat{c}_{{\rm Exc}|{\bf k}}={1\over\Omega_{{\bf k}}}\sum_{j=1}^{N_{\rm mol}}\Omega_{j{\bf k}}{\mathcal{\hat{S}}}_{{\rm Exc}j}{e^{i{\bf k}{\bf r}_{j}}}$; $\hat{c}_{{\rm Vib}|{\bf k}}={1\over\sqrt{N_{\rm mol}}}\sum_{j=1}^{N_{\rm mol}}{\mathcal{\hat{B}}}_{{\rm Vib}j}e^{i{\bf k}{\bf r}_{j}}$_.
  * Carmichael [2009] H. Carmichael,  _An open systems approach to quantum optics: lectures presented at the Université Libre de Bruxelles, October 28 to November 4, 1991_ , Vol. 18 (Springer Science & Business Media, 2009).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Micha and Burghardt [2007] D. A. Micha and I. Burghardt,  _Quantum dynamics of complex molecular systems_ , Vol. 83 (Springer, 2007).
  * Ostroverkhova [2016] O. Ostroverkhova, Chemical reviews 116, 13279 (2016).
  * Martin and Hynes [2004] M. M. Martin and J. T. Hynes,  _Femtochemistry and femtobiology: ultrafast events in molecular science_ (Elsevier, 2004).
  * Bloch  _et al._ [2022] J. Bloch, I. Carusotto, and M. Wouters, Nature Reviews Physics 4, 470 (2022).
  * Kubo [1957] R. Kubo, Journal of the Physical Society of Japan 12, 570 (1957).
  * Martin and Schwinger [1959] P. C. Martin and J. Schwinger, Physical Review 115, 1342 (1959).
  * Kavokin  _et al._ [2017] A. V. Kavokin, J. J. Baumberg, G. Malpuech, and F. P. Laussy,  _Microcavities_ , Vol. 21 (Oxford university press, 2017).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Reviews of Modern Physics 86, 1391 (2014).
  * Shishkov  _et al._ [2022a] V. Y. Shishkov, E. S. Andrianov, A. V. Zasedatelev, P. G. Lagoudakis, and Y. E. Lozovik, Physical Review Letters 128, 065301 (2022a).
  * Shishkov  _et al._ [2022b] V. Y. Shishkov, E. Andrianov, and Y. E. Lozovik, Quantum 6, 719 (2022b).
  * Mukamel [1995] S. Mukamel,  _Principles of nonlinear optical spectroscopy_ (Oxford University Press, 1995).
  * Jonas [2003] D. M. Jonas, Annual review of physical chemistry 54, 425 (2003).
  * Mukamel [2000] S. Mukamel, Annual review of physical chemistry 51, 691 (2000).
  * Song  _et al._ [2015] Y. Song, C. Hellmann, N. Stingelin, and G. D. Scholes, The Journal of chemical physics 142, https://doi.org/10.1063/1.4916325 (2015).
  * Betzold  _et al._ [2019] S. Betzold, M. Dusel, O. Kyriienko, C. P. Dietrich, S. Klembt, J. Ohmer, U. Fischer, I. A. Shelykh, C. Schneider, and S. Höofling, ACS Photonics 7, 384 (2019).
  * Yagafarov  _et al._ [2020] T. Yagafarov, D. Sannikov, A. Zasedatelev, K. Georgiou, A. Baranikov, O. Kyriienko, I. Shelykh, L. Gai, Z. Shen, D. Lidzey, _et al._ , Communications Physics 3, 18 (2020).
  * Askitopoulos  _et al._ [2019] A. Askitopoulos, L. Pickup, S. Alyatkin, A. Zasedatelev, K. G. Lagoudakis, W. Langbein, and P. G. Lagoudakis, arXiv preprint arXiv:1911.08981 https://doi.org/10.48550/arXiv.1911.08981 (2019).
  * Baryshev  _et al._ [2022] S. Baryshev, A. Zasedatelev, H. Sigurdsson, I. Gnusov, J. D. Töpfer, A. Askitopoulos, and P. G. Lagoudakis, Physical Review Letters 128, 087402 (2022).
  * McGhee  _et al._ [2015] K. E. McGhee, M. Guizzardi, R. Jayaprakash, K. Georgiou, T. Jessewitsch, U. Scherf, G. Cerullo, A. Zasedatelev, T. Virgili, P. G. Lagoudakis, and D. G. Lidzey, Advanced Optical Materials , 2300262 (2015).
  * Putintsev  _et al._ [2023] A. D. Putintsev, K. E. McGhee, D. Sannikov, A. V. Zasedatelev, J. D. Töpfer, D. G. Lidzey, and P. G. Lagoudakis, arXiv preprint arXiv:2304.03617 https://doi.org/10.48550/arXiv.2304.03617 (2023).
  * Scholes  _et al._ [2017] G. D. Scholes, G. R. Fleming, L. X. Chen, A. Aspuru-Guzik, A. Buchleitner, D. F. Coker, G. S. Engel, R. Van Grondelle, A. Ishizaki, D. M. Jonas, _et al._ , Nature 543, 647 (2017).
  * Guha  _et al._ [2003] S. Guha, J. Rice, Y. Yau, C. M. Martin, M. Chandrasekhar, H. R. Chandrasekhar, R. Guentner, P. S. De Freitas, and U. Scherf, Physical Review B 67, 125204 (2003).

## Supplementary Materials

##  I  FUll Hamiltonian of the system

We consider a nonequilibrium microscopic model that describes an ensemble of
organic molecules. These molecules feature vibrationally dressed electronic
transitions and are coupled to cavity modes with differing in-plane momenta
$\hbar{\bf k}_{\parallel}$ (hereinafter $\hbar{\bf k}$). The electronic
excitations being localised at molecules are effectively treated as Frenkel-
type excitons, where each excitation can be described by the Pauli creation
and annihilation operators acting on a single molecule. We take into account
one vibrational state with eigenfrequency $\omega_{\rm Vib}$ in each molecule
The full Hamiltonian of the systems is Eq (1) in the main text

$\hat{H}=\sum\limits_{\bf k}\hbar{\omega_{{\rm Cav}|{\bf k}}}\hat{a}_{{\rm
Cav}|{\bf k}}^{\dagger}{{\hat{a}}_{{\rm Cav}|{\bf
k}}}+\sum\limits_{j=1}^{N_{\rm mol}}{\hbar{\omega_{{\rm
exc}}}\hat{\sigma}_{{\rm Exc}{j}}^{\dagger}{{\hat{\sigma}}_{{\rm
Exc}{j}}}}+\\\ \sum\limits_{j=1}^{N_{\rm mol}}{\hbar{\omega_{\rm
Vib}}\hat{b}_{{\rm Vib}j}^{\dagger}{\hat{b}_{{\rm
Vib}j}}}+\sum\limits_{j=1}^{N_{\rm mol}}{\hbar\Lambda\omega_{\rm
Vib}\hat{\sigma}_{{\rm Exc}j}^{\dagger}{{\hat{\sigma}}_{{\rm
Exc}{j}}}\left({{\hat{b}_{{\rm Vib}j}}+\hat{b}_{{\rm
Vib}j}^{\dagger}}\right)}+\\\ \sum\limits_{j=1}^{N_{\rm mol}}\sum\limits_{{\bf
k}}\hbar\Omega_{j{\bf k}}\left(\hat{\sigma}_{{\rm
Exc}{j}}^{\dagger}\hat{a}_{{\rm Cav}|{\bf k}}e^{i{\bf k}{\bf
r}_{j}}+\hat{\sigma}_{{\rm Exc}{j}}\hat{a}_{{\rm Cav}|{\bf
k}}^{\dagger}e^{-i{\bf k}{\bf r}_{j}}\right),$ (10)

From this Hamiltonian, we derive Eqs.Eqs. (2)–(6). To accomplish this, we
transition to the hybridized (or dressed) vibrational and exciton states and
subsequently to the polariton states through coupling to the cavity

###  I.1 Dressed vibrational and exciton states

Strong vibronic coupling is quite common feature of organic molecules.
Conjugated polymers is the one of particularly interesting class of molecular
systems that is characterized by a large Huang–Rhys factor. For example
ladder-type polymers such as poly(p-phenylene) also known as MELPPP typically
have Huang–Rhys factor $\Lambda$ close to unity [61]. The strong interplay
between electronic and vibrational degrees of freedom requires transformation
to the basis of dressed states. Therefore we substitute the initial operators
$\hat{\sigma}_{{\rm Exc}j}$ and $\hat{b}_{{\rm Vib}j}$ with the operators of
the dressed excitons ${\mathcal{\hat{S}}}_{{\rm Exc}j}$ and dressed vibrations
${\mathcal{\hat{B}}}_{{\rm Vib}j}$ as follows

$\hat{\sigma}_{{\rm Exc}j}={\mathcal{\hat{S}}}_{{\rm
Exc}j}e^{-\Lambda({\mathcal{\hat{B}}}_{{\rm
Vib}j}^{\dagger}-{\mathcal{\hat{B}}}_{{\rm Vib}j})},$ (11) $\hat{b}_{{\rm
Vib}j}={\mathcal{\hat{B}}}_{{\rm
Vib}j}-\Lambda{\mathcal{\hat{S}}}^{\dagger}_{{\rm
Exc}j}{\mathcal{\hat{S}}}_{{\rm Exc}j}.$ (12)

This substitution preserve the commutation relations for the operators
${\mathcal{\hat{S}}}_{{\rm Exc}j}$ and ${\mathcal{\hat{B}}}_{{\rm Vib}j}$,
namely the dressed operators should have the same commutator relations as the
initial ones $\hat{\sigma}_{{\rm Exc}j}$ and $\hat{b}_{{\rm Vib}j}$.

The transition to the dressed state representation is equivalent to
diagonalization of the initial vibrational and excitonic parts of the
Hamiltonian (10) leading to

$\hat{H}=\sum\limits_{\bf{k}}\hbar\omega_{{\rm Cav}|{\bf k}}\hat{a}_{{\rm
Cav}|{\bf k}}^{\dagger}\hat{a}_{{\rm Cav}|{\bf k}}+\sum_{j=1}^{N_{\rm
mol}}\hbar\omega_{\rm Exc}{\mathcal{\hat{S}}}^{\dagger}_{{\rm
Exc}j}{\mathcal{\hat{S}}}_{{\rm Exc}j}+\sum_{j=1}^{N_{\rm
mol}}\hbar\omega_{\rm Vib}{\mathcal{\hat{B}}}_{{\rm
Vib}j}^{\dagger}{\mathcal{\hat{B}}}_{{\rm Vib}j}+\\\ \sum\limits_{j=1}^{N_{\rm
mol}}\sum\limits_{{\bf k}}\hbar\Omega_{j{\bf k}}\left(\mathcal{\hat{S}}_{{\rm
Exc}j}^{\dagger}\hat{a}_{{\rm Cav}|{\bf k}}e^{i{\bf k}{\bf
r}_{j}}+\mathcal{\hat{S}}_{{\rm Exc}j}\hat{a}_{{\rm Cav}|{\bf
k}}^{\dagger}e^{-i{\bf k}{\bf r}_{j}}\right)-\\\ \sum\limits_{j=1}^{N_{\rm
mol}}\sum\limits_{{\bf k}}\hbar\Lambda{\Omega_{j{\bf
k}}}\left(\mathcal{\hat{S}}_{{\rm Exc}j}^{\dagger}\mathcal{\hat{B}}_{{\rm
Vib}j}\hat{a}_{{\rm Cav}|{\bf k}}e^{i{\bf k}{\bf
r}_{j}}+\mathcal{\hat{S}}_{{\rm Exc}j}\mathcal{\hat{B}}_{{\rm
Vib}j}^{\dagger}\hat{a}_{{\rm Cav}|{\bf k}}^{\dagger}e^{-i{\bf k}{\bf
r}_{j}}\right),$ (13)

where we introduce the energy of the dressed exciton states $\omega_{\rm
Exc}=\omega_{\rm exc}-\Lambda^{2}\omega_{\rm Vib}$. In the Hamiltonian (13) we
decomposed the operator $e^{-\Lambda({\mathcal{\hat{B}}}_{{\rm
Vib}j}^{\dagger}-{\mathcal{\hat{B}}}_{{\rm Vib}j})}\approx
1-\Lambda({\mathcal{\hat{B}}}_{{\rm Vib}j}^{\dagger}-{\mathcal{\hat{B}}}_{{\rm
Vib}j})$ and omitted the non-resonant terms including $\hat{S}_{{\rm
Exc}j}\hat{B}_{{\rm Vib}j}\hat{a}^{\dagger}_{{\rm Cav}|{\bf k}}$. Thus, in the
Hamiltonian (13) we preserve only resonant terms with coefficients up to
$\Lambda^{2}n_{\rm Vib\\_per\\_mol}$, where $n_{\rm Vib\\_per\\_mol}$ is the
average occupation of the vibrational states per one molecule. After the
derivation of the main equations, we verify that $\Lambda^{2}n_{\rm
Vib\\_per\\_mol}\ll 1$ is always fulfilled for the system described in the
main text.

###  I.2 Bright and dark states

In the Hamiltonian (13) the wavevectors $\bf k$ belong to states of the
electromagnetic field in the cavity. In most experimental systems, the number
of states in an optical cavity is much smaller than the number of the
molecules. Thus, a Fabry–Pérot optical cavity with the fundamental mode
characterized by a wavelength $\lambda$ has the number of state equal to
$N_{\rm states}\approx\pi k_{\rm max}^{2}S/(2\pi)^{2}$, where $S$ is the area
of interest (typically active or illuminated area). In state-of-the-art BEC
experiments [13, 12] these parameters are the following: $\lambda\approx
500~{}{\mu m}$ and $S\approx 500~{}{\mu m}^{2}$, which result in $N_{\rm
states}\approx 10^{4}$ number of states. But, the number of molecules $N_{\rm
mol}$ in the same area is $\approx 10^{8}$, therefore $N_{\rm states}\ll
N_{\rm mol}$.

This analysis shows that the vast majority of exciton states do not actually
couple to the cavity and do not participate in the formation of the polariton
states. The exciton states, that coupled to the cavity we define as the bright
excitons, while the rest states classified as the dark ones. Bright excitons
are phase-coherent, many-body delocalized states with a well-defined in-plane
momentum $\hbar\bf k$, matching the corresponding eigenstates of the cavity.
In contrast the dark excitons, lacking well-defined momentum, represent a
manifold of localized states. Analogous to excitons — with the sole
distinction being their lack of dipole coupling to the cavity (we assume Raman
active transi) — we separate all molecular vibrations as either bright or dark
types. Bright vibrations represent coherent and delocalized states. Similar to
the bright excitons they have a well-defined momentum $\hbar\bf k_{\rm Vib}$,
unlike the localized dark vibrational states.

It is convenient to introduce the collective operator effectively describing
all bright excitons and bright vibrations

$\hat{c}_{{\rm Exc}|{\bf k}}={1\over\Omega_{\bf k}}\sum_{j=1}^{N_{\rm
mol}}\Omega_{j{\bf k}}{\mathcal{\hat{S}}}_{{\rm Exc}j}{e^{i{\bf k}{\bf
r}_{j}}}$ (14) $\hat{c}_{{\rm Vib}|{\bf k}}={1\over\sqrt{N_{\rm
mol}}}\sum_{j=1}^{N_{\rm mol}}{\mathcal{\hat{B}}}_{{\rm Vib}j}e^{i{\bf k}{\bf
r}_{j}}$ (15)

where $\Omega_{\bf k}=\sqrt{\sum_{j}|\Omega_{j{\bf k}}|^{2}}$ is the Rabi
frequency. These operators obey the commutation relations $\left[\hat{c}_{{\rm
Exc}|{\bf k}},\hat{c}_{{\rm Exc}|{\bf
k^{\prime}}}^{\dagger}\right]=\delta_{{\bf k},{\bf k^{\prime}}}$ and
$\left[\hat{c}_{{\rm Vib}|{\bf k}},\hat{c}_{{\rm Vib}|{\bf
k^{\prime}}}^{\dagger}\right]=\delta_{{\bf k},{\bf k^{\prime}}}$. Using
operators (14)–(15) and assuming $\Omega_{j{\bf k}}\approx\Omega_{j{\bf
k}^{\prime}}$ we transform the Hamiltonian (13)

$\hat{H}=\sum\limits_{\bf k}\hbar\omega_{{\rm Cav}|{\bf k}}\hat{a}_{{\rm
Cav}|{\bf k}}^{\dagger}\hat{a}_{{\rm Cav}|{\bf k}}+\sum\limits_{\bf
k}\hbar\omega_{\rm Exc}\hat{c}_{{\rm Exc}|{\bf k}}^{\dagger}\hat{c}_{{\rm
Exc}|{\bf k}}+\sum\limits_{\bf{k}}\hbar\omega_{\rm Vib}\hat{c}_{{\rm Vib}|{\bf
k}}^{\dagger}\hat{c}_{{\rm Vib}|{\bf k}}+\\\ \sum\limits_{\bf
k}\hbar\Omega_{\bf k}\left(\hat{c}_{{\rm Exc}|{\bf k}}^{\dagger}\hat{a}_{{\rm
Cav}|{\bf k}}+\hat{c}_{{\rm Exc}|{\bf k}}\hat{a}_{{\rm Cav}|{\bf
k}}^{\dagger}\right)-\\\ \sum\limits_{{\bf k},{\bf
k}^{\prime}}\hbar{\Lambda\Omega_{\bf k}\over\sqrt{N_{\rm
mol}}}\left(\hat{c}_{{\rm Exc}{\bf k^{\prime}}}^{\dagger}\hat{c}_{{\rm
Vib}|{\bf k^{\prime}}-{\bf k}}\hat{a}_{{\rm Cav}|{\bf k}}+\hat{c}_{{\rm
Exc}|{\bf k^{\prime}}}\hat{c}_{{\rm Vib}|{\bf k^{\prime}}-{\bf
k}}^{\dagger}\hat{a}_{{\rm Cav}|{\bf k}}^{\dagger}\right)+\\\ N_{\rm
mol}\hbar\omega_{\rm Exc}\hat{n}_{\rm Exc_{D}}+N_{\rm mol}\hbar\omega_{\rm
Vib}\hat{n}_{\rm Vib_{D}}+\hat{H}_{\rm dark\\_int},$ (16)

We denote the operator of the number of the dark excitons per one molecule as
$\hat{n}_{\rm Exc_{D}}$ and the similar operator for the dark vibrations as
$\hat{n}_{\rm Vib_{D}}$

$\hat{n}_{\rm Exc_{D}}={1\over N_{\rm mol}}\sum_{j=1}^{N_{\rm
mol}}{\mathcal{\hat{S}}}^{\dagger}_{{\rm Exc}j}{\mathcal{\hat{S}}}_{{\rm
Exc}j}-{1\over N_{\rm mol}}\sum\limits_{\bf{k}}\hat{c}_{{\rm Exc}|{\bf
k}}^{\dagger}\hat{c}_{{\rm Exc}|{\bf k}},$ (17) $\hat{n}_{\rm Vib_{D}}={1\over
N_{\rm mol}}\sum_{j=1}^{N_{\rm mol}}{\mathcal{\hat{B}}}_{{\rm
Vib}j}^{\dagger}{\mathcal{\hat{B}}}_{{\rm Vib}j}-{1\over N_{\rm
mol}}\sum\limits_{\bf k}\hat{c}_{{\rm Vib}|{\bf k}}^{\dagger}\hat{c}_{{\rm
Vib}|{\bf k}}.$ (18)

We also introduced the Hamiltonian of the indirect interactions between the
dark states and cavity modes $\hat{H}_{{\rm dark\\_int}}$

$\hat{H}_{\rm dark\\_int}=-\sum\limits_{j=1}^{N_{\rm mol}}\sum\limits_{\bf
k}\hbar\Lambda\Omega_{j{\bf k}}\left(\mathcal{\hat{S}}_{{\rm
Exc}j}^{\dagger}\mathcal{\hat{B}}_{{\rm Vib}j}\hat{a}_{{\rm Cav}|{\bf
k}}e^{i{\bf k}{\bf r}_{j}}+\mathcal{\hat{S}}_{{\rm
Exc}j}\mathcal{\hat{B}}_{{\rm Vib}j}^{\dagger}\hat{a}_{{\rm Cav}|{\bf
k}}^{\dagger}e^{-i{\bf k}{\bf r}_{j}}\right)+\\\ \sum\limits_{{\bf k},{\bf
k}^{\prime}}\hbar{\Lambda\Omega_{\bf k}\over\sqrt{N_{\rm
mol}}}\left(\hat{c}_{{\rm Exc}|{\bf k^{\prime}}}^{\dagger}\hat{c}_{{\rm
Vib}|{\bf k^{\prime}}-{\bf k}}\hat{a}_{{\rm Cav}|{\bf k}}+\hat{c}_{{\rm
Exc}|{\bf k^{\prime}}}\hat{c}_{{\rm Vib}|{\bf k^{\prime}}-{\bf
k}}^{\dagger}\hat{a}_{{\rm Cav}|{\bf k}}^{\dagger}\right).$ (19)

###  I.3 Polariton states

The strong interaction between excitons and electromagnetic field of the
cavity leads to the formation of the new collective states – exciton-
polaritons. It follows from the Hamiltonian (16), that only bright excitons
directly interacts with the cavity. Unlike bare photons in the cavity,
polaritons have the material part, that allows them to thermalize towards BEC
above condensation threshold. We introduce operators for the lower
$\hat{s}_{{\rm Low}|{\bf k}}$ and upper $\hat{s}_{{\rm Up}|{\bf k}}$
polaritons

$\hat{s}_{{\rm Low}|{\bf k}}=\hat{a}_{{\rm Cav}|{\bf k}}\cos\varphi_{{\bf
k}}-\hat{c}_{{\rm Exc}|{\bf k}}\sin\varphi_{{\bf k}}$ (20) $\hat{s}_{{\rm
Up}|{\bf k}}=\hat{a}_{{\rm Cav}|{\bf k}}\sin\varphi_{{\bf k}}+\hat{c}_{{\rm
Exc}|{\bf k}}\cos\varphi_{{\bf k}},$ (21)

where

$\varphi_{\bf k}={1\over 2}{\rm arctg}\left(\frac{2\Omega_{\bf k}}{\omega_{\rm
Exc}-\omega_{{\rm Cav}|{\bf k}}}\right).$ (22)

After omitting the non-resonant terms, we obtain

$\hat{H}=\sum\limits_{\bf k}\hbar\omega_{{\rm Low}|{\bf k}}\hat{s}_{{\rm
Low}|{\bf k}}^{\dagger}\hat{s}_{{\rm Low}|{\bf k}}+\sum\limits_{\bf
k}\hbar\omega_{{\rm Up}|{\bf k}}\hat{s}_{{\rm Up}|{\bf
k}}^{\dagger}\hat{s}_{{\rm Up}|{\bf k}}+\sum\limits_{\bf k}\hbar\omega_{\rm
Vib}\hat{c}_{{\rm Vib}|{\bf k}}^{\dagger}\hat{c}_{{\rm Vib}|{\bf k}}-\\\
\sum\limits_{{\bf k},{\bf k}^{\prime}}{\hbar\Lambda\Omega_{\bf
k}\over\sqrt{N_{\rm mol}}}\cos\varphi_{\bf k}\cos\varphi_{\bf
k^{\prime}}\left[\left(\cos\varphi_{\bf k^{\prime}}\hat{s}_{{\rm Up}|{\bf
k^{\prime}}}^{\dagger}-\sin\varphi_{\bf k^{\prime}}\hat{s}_{{\rm Low}|{\bf
k^{\prime}}}^{\dagger}\right)\hat{c}_{{\rm Vib}|{\bf k^{\prime}}-{\bf
k}}\hat{s}_{{\rm Low}|{\bf k}}+h.c.\right]+\\\ N_{\rm mol}\hbar\omega_{\rm
Exc}\hat{n}_{\rm Exc_{D}}+N_{\rm mol}\hbar\omega_{\rm Vib}\hat{n}_{\rm
Vib_{D}}+\hat{H}_{\rm dark\\_int},$ (23)

where the new eigenfrequencies are:

$\omega_{{\rm Low}|{\bf k}}={\omega_{\rm Exc}+\omega_{{\rm Cav}|{\bf k}}\over
2}-\sqrt{{\left(\omega_{\rm Exc}-\omega_{{\rm Cav}|{\bf k}}\right)^{2}\over
4}+\Omega_{\bf k}^{2}}$ (24) $\omega_{{\rm Up}|{\bf k}}={\omega_{\rm
Exc}+\omega_{{\rm Cav}|{\bf k}}\over 2}+\sqrt{{\left(\omega_{\rm
Exc}-\omega_{{\rm Cav}|{\bf k}}\right)^{2}\over 4}+\Omega_{\bf k}^{2}}$ (25)

###  I.4 Interaction between polaritons, excitons and molecular vibrations

Even though dark excitons do not directly interact with the cavity to form
polariton states, they still significantly contribute to the polariton
transitions as evident from their commutation relations:

$\left[\hat{H}_{l{\rm dark\\_int}},\hat{n}_{\rm Exc_{D}}\right]\neq 0,$ (26)
$\left[\hat{H}_{l{\rm dark\\_int}},\hat{n}_{\rm Vib_{D}}\right]\neq 0,$ (27)
$\left[\hat{H}_{\rm dark\\_int},\sum\limits_{\bf k}\hbar\omega_{{\rm Low}|{\bf
k}}\hat{s}_{{\rm Low}|{\bf k}}^{\dagger}\hat{s}_{{\rm Low}|{\bf k}}\right]\neq
0.$ (28)

The interaction between molecular vibrations and excitons leads to the
transition of the energy from the dark excitons to the dressed vibrations and
lower polariton states. This interaction has a resonant nature and corresponds
to the following microscopic process: the annihilation of a dark exciton
followed by the creation of a lower polariton along with a dressed vibration.
Since the total number of dark exciton states greatly exceeds the number of
polariton states, this process has a become increasingly important below
condensation threshold.

##  II Interaction with the environment: pumping, relaxation, thermalization

The full state of the system can be characterized by its density matrix
$\hat{\rho}$. The density matrix is governed by the master equation [34, 38]

${\frac{d\hat{\rho}}{dt}}={\frac{i}{\hbar}}{\left[{\hat{\rho}},{\hat{H}}\right]}+\sum_{n}L_{\hat{A}_{n}}(\hat{\rho}),$
(29)

where $L_{\hat{A}_{n}}(\hat{\rho})$ is the Lindblad superoperators and
$\hat{A}_{n}$ is the relaxation operator,

$L_{\hat{A}}(\hat{\rho})=\hat{A}\hat{\rho}\hat{A}^{\dagger}-{1\over
2}\hat{\rho}\hat{A}^{\dagger}\hat{A}-{1\over
2}\hat{A}^{\dagger}\hat{A}\hat{\rho}$ (30)

The Lindblad superoperators can describe both dissipation and pumping of the
system.

Below we consider the relaxation processes in each of the subsystems
separately. For the sake of briefness, below we explicitly define only
relaxation operators $\hat{A}_{n}$.

The finite photon lifetime in the cavity as well as absorption/scattering
losses in the media result in dissipation of electromagnetic field in a cavity
mode $\bf k$ with the rate $\gamma_{{\rm Cav}|{\bf k}}$ This relaxation
process can be described by the relaxation operator $\sqrt{\gamma_{{\rm
Cav}|{\bf k}}}\hat{a}_{{\rm Cav}|{\bf k}}$.

Energy relaxation rate of dressed exciton states in the absence of a cavity is
defined by $\gamma_{\rm Exc}$. The corresponding relaxation operator for
$j$-th molecule is $\sqrt{\gamma_{\rm Exc}}\hat{\mathcal{\hat{S}}}_{{\rm
Exc}j}$. The dressed exciton states also undergo dephasing processes at the
rate $\Gamma_{\rm Exc}$ that corresponds to the relaxation operator
$\sqrt{\Gamma_{\rm Exc}}\hat{\mathcal{\hat{S}}}_{{\rm
Exc}j}^{\dagger}\hat{\mathcal{\hat{S}}}_{{\rm Exc}j}$. Usually, in molecular
systems at room temperature the exciton energy dissipation rate $\gamma_{\rm
Exc}$ is much smaller than the dephasing $\Gamma_{\rm Exc}$.

The dressed vibrational states dissipate their energy with the rate
$\gamma_{\rm Vib}$. Thermal fluctuations are important for the system at room
temperature, therefore the corresponding relaxation operators of dressed
vibrational states for the $j$-th molecule are $\sqrt{(1+n_{\rm Vib}^{\rm
th})\gamma_{\rm Vib}}\mathcal{\hat{B}}_{{\rm Vib}j}$ and $\sqrt{n_{\rm
Vib}^{\rm th}\gamma_{\rm Vib}}\mathcal{\hat{B}}_{{\rm Vib}j}^{\dagger}$, where
$n_{\rm Vib}^{\rm th}=1/({\rm exp}(\hbar\omega_{\rm Vib}/T)-1)$ is the mean of
the thermal distribution of molecular vibrations at an effective temperature
$T$.

The dissipation rate remains the same $\gamma_{\rm Vib}$ for both: dark and
bright vibrational states. The same applied for dark excitons, their energy
dissipation rate is equal to $\gamma_{\rm Exc}$. However, due to strong
interaction between cavity and bright exciton states the dissipation rates of
the polariton states is different. The dissipation rates of the lower
polariton states $\gamma_{{\rm Low}|{\bf k}}$ and for upper polariton states
$\gamma_{{\rm Up}|{\bf k}}$ are determined by the Hopfield coefficients and
equal to $\gamma_{{\rm Low}|{\bf k}}=\gamma_{{\rm Cav}|{\bf
k}}\cos^{2}\varphi_{\bf k}+\gamma_{{\rm Exc}}\sin^{2}\varphi_{\bf k}$ and
$\gamma_{{\rm Up}|{\bf k}}=\gamma_{\rm Exc}\cos^{2}\varphi_{\bf
k}+\gamma_{{\rm Cav}|{\bf k}}\sin^{2}\varphi_{\bf k}$. The corresponding
relaxation operators are $\sqrt{\gamma_{{\rm Low}|{\bf k}}}\hat{s}_{{\rm
Low}|{\bf k}}$ and $\sqrt{\gamma_{{\rm Up}|{\bf k}}}\hat{s}_{{\rm Up}|{\bf
k}}$.

Importantly polariton states undergo thermalization due to their matter
component. Thus, thermalization processes between lower polariton states with
the wavevector ${\bf k}_{1}$ and ${\bf k}_{2}$ can be described by the
relaxation operators $\sqrt{\gamma_{\rm therm}^{{{\bf k}_{1}}{{\bf
k}_{2}}}}{{\hat{s}}_{{\rm Low}|{{\bf k}_{2}}}}\hat{s}_{{\rm Low}|{{\bf
k}_{1}}}^{\dagger}$ and $\sqrt{\gamma_{\rm therm}^{{{\bf k}_{2}}{{\bf
k}_{1}}}}{{\hat{s}}_{{\rm Low}|{{\bf k}_{1}}}}\hat{s}_{{\rm Low}|{{\bf
k}_{2}}}^{\dagger}$. The relation between upward and downward thermalization
processes $\gamma_{{\rm therm}}^{{{\bf k}_{1}}{{\bf k}_{2}}}$ and
$\gamma_{{\rm therm}}^{{{\bf k}_{2}}{{\bf k}_{1}}}$ is determined by the
Kubo–Martin–Schwinger relation

$\gamma_{\rm therm}^{{{\bf k}_{1}}{{\bf k}_{2}}}=\gamma_{\rm therm}^{{{\bf
k}_{2}}{{\bf k}_{1}}}\exp\left(\frac{\omega_{{\rm Low}|{{\bf
k}_{2}}}-\omega_{{\rm Low}|{{\bf k}_{1}}}}{T}\right)$

where $T$ is the temperature of the environment.

We consider the following pumping scheme: the external laser source excites
the lower and upper polaritons with the wavevector ${\bf k}={\bf k}_{\rm ex}$.
This process can be described by the relaxation operators $\sqrt{\gamma_{{\rm
Low}|{\bf k}_{\rm ex}}\varkappa_{\rm Low}}\hat{s}_{{\rm Low}|{\bf k}_{\rm
ex}}$, $\sqrt{\gamma_{{\rm Low}|{\bf k}_{\rm ex}}\varkappa_{\rm
Low}}\hat{s}_{{\rm Low}|{\bf k}_{\rm ex}}^{\dagger}$, $\sqrt{\gamma_{{\rm
Up}|{\bf k}_{\rm ex}}\varkappa_{\rm Up}}\hat{s}_{{\rm Up}|{\bf k}_{\rm ex}}$
and $\sqrt{\gamma_{{\rm Up}|{\bf k}_{\rm ex}}\varkappa_{\rm Up}}\hat{s}_{{\rm
Up}|{\bf k}_{\rm ex}}^{\dagger}$ [48, 49]. We assume that ${\bf k}={\bf
k}_{\rm ex}$ is large enough so that $\omega_{{\rm Low}|{\bf k}_{\rm
ex}}\approx\omega_{\rm Exc}$ and ${|\omega_{{\rm Low}|{\bf k}_{\rm
ex}}-\omega_{\rm Exc}|\gg\Gamma_{\rm Exc}}$, also we assume that the external
source is tuned to the frequency $\omega_{\rm Exc}$ in agreement with state-
of-the-art experiments [13]. Therefore, we can omit the dynamics of the upper
polaritons only with the wavevector ${\bf k}_{\rm ex}$. We denote
$\varkappa_{\rm Pump}=\varkappa_{\rm Low}$.

The lower polariton states with large wavevectors possess energies nearly
identical to those of the dark excitons. Despite this, the total number of
dark exciton states significantly outnumbers that of the polariton states. As
a result, the dark excitons effectively act as a reservoir for the polariton
states with large wavevectors, especially for the lower polariton state with a
given wavevector ${\bf k}_{\rm ex}$. In equations (2)–(6) (main text) we
phenomenologically account for this mechanism, assuming that the corresponding
relaxation rate is $\gamma_{\rm Exc}^{\rm B-D}=\Gamma_{\rm Exc}$. Similar
principle applied to the bright and dark vibrational states. We set the
corresponding relaxation rate $\gamma_{\rm Vib}^{\rm B-D}$ to $\gamma_{\rm
Vib}/2$.

##  III Evolution of average occupation numbers

From the Hamiltonian (23) and master equation (29) we obtain the equations for
the average lower polaritons with the wavevector $\bf k$

$n_{{\rm Pol}|{\bf k}}=\langle\hat{s}_{{\rm Low}|{\bf
k}}^{\dagger}\hat{s}_{{\rm Low}|{\bf k}}\rangle,$ (31)

average number of the dark excitons per molecule

$n_{\rm Exc_{D}}=\langle\hat{n}_{\rm
Exc_{D}}\rangle\approx\sum_{j}\langle\hat{S}_{{\rm
Exc}j}^{\dagger}\hat{S}_{{\rm Exc}j}\rangle/N_{\rm mol},$ (32)

average number of the dark vibrations per molecule

$n_{\rm Vib_{D}}=\langle\hat{n}_{\rm
Vib_{D}}\rangle\approx\sum_{j}\langle\mathcal{\hat{B}}_{{\rm
Vib}j}^{\dagger}\mathcal{\hat{B}}_{{\rm Vib}j}\rangle/N_{\rm mol}.$ (33)

We separately consider the pumped bright excitons and lower polariton states
with the wavevector ${\bf k}_{\rm ex}$ and denote the corresponding average

$n_{{\rm Exc}|{\bf k}_{\rm ex}}=\langle\hat{s}_{{\rm Low}|{\bf k}_{\rm
ex}}^{\dagger}\hat{s}_{{\rm Low}|{\bf k}_{\rm ex}}\rangle.$ (34)

We also consider the average number of bright vibrations, that assist the
transition between the lower polariton state with the wavevector ${\bf k}_{\rm
ex}$ to the lower polariton state with the wavevector $\bf k$,

$n_{{\rm Vib}|{\bf k_{\rm ex}}-{\bf k}}=\langle\hat{c}_{{\rm Vib}|{\bf k}_{\rm
ex}-{\bf k}}^{\dagger}\hat{c}_{{\rm Vib}|{\bf k}_{\rm ex}-{\bf k}}\rangle.$
(35)

Assuming that $\Gamma_{\rm Exc}\gg\gamma_{\rm Exc},\gamma_{{\rm Low}|{\bf
k}}$, $\gamma_{{\rm Vib}}$ and that $\Omega_{\bf k}$ do not depend on $\bf k$
we obtain the following equations

$\frac{dn_{{\rm Exc}|{\bf k}_{\rm ex}}}{dt}=-\gamma_{{\rm Low}|{\bf k}_{\rm
ex}}(n_{{\rm Exc}|{\bf k}_{\rm ex}}-\varkappa_{\rm Pump})+\gamma_{\rm
Exc}^{\rm B-D}(n_{\rm Exc_{D}}-n_{{\rm Exc}|{\bf k}_{\rm ex}})-\\\ \sum_{\bf
k}\frac{\tilde{G}_{{\bf k}}}{N_{\rm mol}}\left[n_{{\rm Exc}|{\bf k}_{\rm
ex}}\left(n_{{\rm Pol}|{\bf k}}+1\right)+n_{{\rm Vib}|{\bf k_{\rm ex}}-{\bf
k}}\left(n_{{\rm Exc}|{\bf k}_{\rm ex}}-n_{{\rm Pol}|{\bf k}}\right)\right]$
(36) $\frac{dn_{\rm Exc_{D}}}{dt}=-\gamma_{\rm Exc}n_{\rm
Exc_{D}}+{\gamma_{\rm Exc}^{\rm B-D}\over N_{\rm mol}}(n_{{\rm Exc}|{\bf
k}_{\rm ex}}-n_{\rm Exc_{D}})-\sum_{\bf k}\frac{G_{\bf k}}{N_{\rm
mol}}\left[n_{\rm Exc_{D}}\left(n_{{\rm Pol}|{\bf k}}+1\right)+n_{\rm
Vib_{D}}\left(n_{\rm Exc_{D}}-n_{{\rm Pol}|{\bf k}}\right)\right]$ (37)
$\frac{dn_{{\rm Pol}|{\bf k}}}{dt}=-\gamma_{{\rm Low}|{\bf k}}n_{\bf k}+G_{\bf
k}\left[n_{\rm Exc_{D}}\left(n_{{\rm Pol}|{\bf k}}+1\right)+n_{\rm
Vib_{D}}\left(n_{\rm Exc_{D}}-n_{{\rm Pol}|{\bf k}}\right)\right]+\\\
{\tilde{G}_{\bf k}\over N_{\rm mol}}\left[n_{{\rm Exc}|{\bf k}_{\rm
ex}}\left(n_{{\rm Pol}|{\bf k}}+1\right)+n_{{\rm Vib}|{\bf k_{\rm ex}}-{\bf
k}}\left(n_{{\rm Exc}|{\bf k}_{\rm ex}}-n_{{\rm Pol}|{\bf
k}}\right)\right]+\\\ \sum_{\bf k^{\prime}}\left\\{\gamma_{\rm therm}^{\bf
kk^{\prime}}\left(n_{{\rm Pol}|{\bf k}}+1\right){n_{{\rm Pol}|{\bf
k^{\prime}}}}-\gamma_{{\rm{therm}}}^{{\bf{k^{\prime}k}}}\left(n_{{\rm
Pol}|{\bf k^{\prime}}}+1\right){n_{{\rm Pol}|{\bf k}}}\right\\}$ (38)
$\frac{dn_{\rm Vib_{D}}}{dt}=-\gamma_{\rm Vib}\left(n_{\rm Vib_{D}}-n_{\rm
Vib}^{\rm th}\right)+{\gamma_{\rm Vib}^{\rm B-D}\over N_{\rm mol}}\sum_{\bf
k}\left(n_{{\rm Vib}|{\bf k_{\rm ex}}-{\bf k}}-n_{\rm Vib_{D}}\right)+\\\
\sum_{\bf k}\frac{G_{\bf k}}{N_{\rm mol}}\left[n_{\rm Exc_{D}}\left(n_{{\rm
Pol}|{\bf k}}+1\right)+n_{\rm Vib_{D}}\left(n_{\rm Exc_{D}}-n_{{\rm Pol}|{\bf
k}}\right)\right]$ (39) $\frac{dn_{{\rm Vib}|{\bf k_{\rm ex}}-{\bf
k}}}{dt}=-\gamma_{\rm Vib}\left(n_{{\rm Vib}|{\bf k_{\rm ex}}-{\bf k}}-n_{{\rm
Vib}}^{\rm th}\right)+\gamma_{\rm Vib}^{\rm B-D}\left(n_{\rm Vib_{D}}-n_{{\rm
Vib}|{\bf k_{\rm ex}}-{\bf k}}\right)+\\\ \frac{\tilde{G}_{\bf k}}{N_{\rm
mol}}\left[n_{{\rm Exc}|{\bf k}_{\rm ex}}\left(n_{{\rm Pol}|{\bf
k}}+1\right)+n_{{\rm Vib}|{\bf k_{\rm ex}}-{\bf k}}\left(n_{{\rm Exc}|{\bf
k}_{\rm ex}}-n_{{\rm Pol}|{\bf k}}\right)\right]$ (40)

where $\tilde{G}_{\bf k}=G_{\bf k}\sin^{2}\varphi_{{\bf k}_{\rm ex}}$ and
$\gamma_{{\rm Low}|{\bf k}_{\rm ex}}=\gamma_{{\rm Cav}|{\bf k}_{\rm
ex}}\cos^{2}\varphi_{{\bf k}_{\rm ex}}+\gamma_{{\rm Exc}}\sin^{2}\varphi_{{\bf
k}_{\rm ex}}$.

Under the following assumption: $\varphi_{{\bf k}_{\rm ex}}\approx\pi/2$ we
substitute $\tilde{G}_{\bf k}\approx G_{\bf k}$, $\gamma_{{\rm Low}|{\bf
k}_{\rm ex}}\approx\gamma_{\rm Exc}$ in Eqs. (36)–(40), that leads to Eqs.
(2)–(6) represented in the main text. We also denote $\gamma_{{\rm Pol}|{\bf
k}}=\gamma_{{\rm Low}|{\bf k}}$ and $\omega_{{\rm Pol}|{\bf k}}=\omega_{{\rm
Low}|{\bf k}}$.

Here, we consider the cavity with quadratic dispersion relation

$\omega_{{\rm Cav}|{\bf k}}=\omega_{{\rm Cav}|{\bf k}={\bf 0}}+\alpha_{\rm
Cav}{\bf k}^{2}$ (41)

where under the blue-detuned resonant configuration $\omega_{\rm
Exc}-\omega_{{\rm Pol}|{\bf k}={\bf 0}}\approx\omega_{\rm Vib}$, from Eq. (24)
it becomes

$\omega_{{\rm Pol}|{\bf k}}=\omega_{{\rm Pol}|{\bf k}={\bf 0}}+\alpha_{\rm
Pol}{\bf k}^{2}$ (42)

with

$\alpha_{\rm Pol}=\alpha_{\rm Cav}{\omega_{\rm Vib}^{2}\over\omega_{\rm
Vib}^{2}+\Omega_{R}^{2}}$ (43)

##  IV  Discretization of the equations (2)–(6) from the main text

To discretize the system we map its momenta states onto the energy space,
moving from the continuum of wavevectors to a discrete set of frequencies. We
consider a discrete finite set of frequencies (${\omega_{0}}$, ${\omega_{1}}$,
${\omega_{2}}$, …, ${\omega_{N}}$) with a sampling rate equal to
$\delta\omega$, such that ${\omega_{j+1}}-{\omega_{j}}=\delta\omega$. The
total number of frequencies that we account for is $N+1$. For the sake of
simplicity we introduce the ground state frequency as
${\omega_{0}}={\omega_{{\rm Pol}|{\bf k}={\bf 0}}}$. We denote the frequency
intervals $R_{j}=(\omega_{j}-\delta\omega,\omega_{j}]$ and introduce
occupation numbers and parameters of the discretized system

$n_{{\rm Pol}0}=n_{{\rm Pol}|{\bf k}={\bf 0}},$ (44) $n_{{\rm
Pol}j}=\sum_{{\bf k}:\;\omega_{{\rm Pol}|{\bf k}}\in R_{j}}n_{{\rm Pol}|{\bf
k}},\;\;1<j<N,$ (45) $n_{{\rm Vib}0}=n_{{\rm Vib}|{\bf k}={\bf 0}},$ (46)
$n_{{\rm Vib}j}=\sum_{{\bf k}:\;\omega_{{\rm Pol}|{\bf k}}\in R_{j}}n_{{\rm
Vib}|{\bf k_{\rm ex}}-{\bf k}},\;\;1<j<N,$ (47) $\gamma_{{\rm
Pol}j}=\,\left.\gamma_{{\rm Pol}|{\bf k}}\right|_{\omega_{{\rm Pol}|{\bf
k}}=\omega_{j}},$ (48) $\gamma_{\rm therm}^{jm}=\left.\gamma_{\rm therm}^{\bf
k^{\prime}k^{\prime\prime}}\right|_{\omega_{{\rm Pol}|{\bf
k^{\prime}}}=\omega_{j},\,\omega_{{\rm Pol}|{\bf
k^{\prime\prime}}}=\omega_{m}},$ (49)
$G_{j}=\frac{\Lambda^{2}\Omega_{R}^{2}\Gamma_{\rm Exc}}{(\omega_{\rm
Exc}-\omega_{j}-\omega_{\rm Vib})^{2}+(\Gamma_{\rm Exc}/2)^{2}},$ (50)
$\tilde{G}_{j}=G_{j}\sin^{2}\varphi_{{\bf k}_{\rm ex}},$ (51)

The total number of the lower polariton states within the frequency interval
$\delta\omega$ close to $\omega_{j}$ we denote as $D_{j}$ and define as
follows

$D_{0}=1$ (52) $D_{j}=\sum_{{\bf k}:\;\omega_{{\rm Pol}|{\bf k}}\in
R_{j}}1,\;\;1<j<N.$ (53)

One can show that $D_{j}=S\delta\omega/(4\pi\alpha_{\rm Cav})$ for $1<j<N$
where $S$ is the pumped area.

As a result we obtain the following discrete version of equations above

$\frac{dn_{{\rm Exc}|{\bf k}_{\rm ex}}}{dt}=-\gamma_{\rm Exc}(n_{{\rm
Exc}|{\bf k}_{\rm ex}}-\varkappa_{\rm pump})+\gamma_{\rm Exc}^{\rm B-D}(n_{\rm
Exc_{D}}-n_{{\rm Exc}|{\bf k}_{\rm ex}})-\\\
\sum_{m=0}^{N}\frac{\tilde{G}_{m}}{N_{\rm mol}}\left[n_{{\rm Exc}|{\bf k}_{\rm
ex}}\left(n_{{\rm Pol}m}+D_{m}\right)+n_{{\rm Vib}m}\left(n_{{\rm Exc}|{\bf
k}_{\rm ex}}-{n_{{\rm Pol}m}\over D_{m}}\right)\right]$ (54)

$\frac{dn_{\rm Exc_{D}}}{dt}=-\gamma_{\rm Exc}n_{\rm Exc_{D}}+{\gamma_{\rm
Exc}^{\rm B-D}\over N_{\rm mol}}(n_{{\rm Exc}|{\bf k}_{\rm ex}}-n_{\rm
Exc_{D}})-\sum_{m=0}^{N}\frac{G_{m}}{N_{\rm mol}}\left[n_{\rm
Exc_{D}}\left(n_{m}+D_{m}\right)+n_{\rm Vib_{D}}\left(D_{m}n_{\rm
Exc_{D}}-n_{m}\right)\right]$ (55)

$\frac{dn_{{\rm Pol}j}}{dt}=-\gamma_{{\rm Pol}j}n_{{\rm
Pol}j}+G_{j}\left[n_{\rm Exc_{D}}\left(n_{{\rm Pol}j}+D_{j}\right)+n_{\rm
Vib_{D}}\left(D_{j}n_{\rm Exc_{D}}-n_{{\rm Pol}j}\right)\right]+\\\
{G_{j}\over N_{\rm mol}}\left[n_{{\rm Exc}|{\bf k}_{\rm ex}}\left(n_{{\rm
Pol}j}+D_{j}\right)+n_{{\rm Vib}j}\left(n_{{\rm Exc}|{\bf k}_{\rm
ex}}-{n_{{\rm Pol}j}\over D_{j}}\right)\right]+\\\
\sum_{m=0}^{N}\left\\{\gamma_{\rm therm}^{{jm}}\left(n_{{\rm
Pol}j}+D_{j}\right){n_{{\rm Pol}m}}-\gamma_{{\rm therm}}^{mj}\left(n_{{\rm
Pol}m}+D_{m}\right)n_{{\rm Pol}j}\right\\}$ (56) $\frac{dn_{\rm
Vib_{D}}}{dt}=-\gamma_{\rm Vib}\left(n_{\rm Vib_{D}}-n_{\rm Vib}^{\rm
th}\right)+{\gamma_{\rm Vib}^{\rm B-D}\over N_{\rm mol}}\sum_{m=0}^{N}(n_{{\rm
Vib}m}-D_{m}n_{\rm Vib_{D}})+\\\ \sum_{m=0}^{N}{G_{m}\over N_{\rm
mol}}\left[n_{\rm Exc_{D}}\left(n_{{\rm Pol}m}+D_{m}\right)+n_{\rm
Vib_{D}}\left(D_{m}n_{\rm Exc_{D}}-n_{{\rm Pol}m}\right)\right]$ (57)
$\frac{dn_{{\rm Vib}j}}{dt}=-\gamma_{\rm Vib}\left(n_{{\rm Vib}j}-D_{j}n_{\rm
Vib}^{\rm th}\right)+\gamma_{\rm Vib}^{\rm B-D}(D_{j}n_{\rm Vib_{D}}-n_{{\rm
Vib}j})+\\\ \frac{G_{j}}{N_{\rm mol}}\left[n_{{\rm Exc}|{\bf k}_{\rm
ex}}\left(n_{{\rm Pol}j}+D_{j}\right)+n_{{\rm Vib}j}\left(n_{{\rm Exc}|{\bf
k}_{\rm ex}}-{n_{{\rm Pol}j}\over D_{j}}\right)\right],$ (58)

where $j=0,1,...,N$. The total amount of the equations is $2N+7$.

To solve the equations of motion for polariton and vibrational states we
generate the system of differential equations by discretizing the wavevector
space. These differential equations include the pumped state at the lower
polariton branch ${{\bf{k}}_{{\rm ex}}}$, the set of equidistant $101$ modes
in frequency space that fill the bottom of lower polariton dispersion within
the momenta range of [$|k_{\rm min}|=0$ $\mu{\rm m}^{-1},|k_{\rm max}|=2$
$\mu{\rm m}^{-1}$] as well as $101$ modes for bright vibrational states, an
equation for dark excitons and an equation for dark vibrations.

##  V  Analytic approximation to the condensation threshold

To facilitate further derivation of the analytic approximation we introduce
the following notations

$I=n_{{\rm Exc}|{\bf k}_{\rm ex}}+N_{\rm mol}n_{\rm Exc_{D}},$ (59)
$P=\sum_{\bf k}n_{{\rm Pol}|{\bf k}}.$ (60) $V=N_{\rm mol}n_{\rm
Vib_{D}}+\sum_{\bf k}n_{{\rm Vib}|{\bf k}},$ (61)

From the Eqs. (2)–(6) in the main text we obtain

${dI\over dt}=-\gamma_{\rm Exc}(I-\varkappa_{\rm Pump})-{1\over N_{\rm
mol}}I(PG_{{\bf k}={\bf 0}}+G_{S})-{1\over N_{\rm mol}}J$ (62) ${dP\over
dt}=-\gamma_{{\rm Pol}|{\bf k}={\bf 0}}P+{1\over N_{\rm mol}}I(PG_{{\bf
k}={\bf 0}}+G_{S})+{1\over N_{\rm mol}}J$ (63) ${dV\over dt}=-\gamma_{\rm
Vib}(V-N_{\rm mol}n_{\rm Vib}^{\rm th})+{1\over N_{\rm mol}}I(PG_{{\bf k}={\bf
0}}+G_{S})+{1\over N_{\rm mol}}J$ (64)

were we assume $\tilde{G}_{\bf k}\approx G_{\bf k}$, $\gamma_{{\rm Pol}|{\bf
k}_{\rm ex}}\approx\gamma_{\rm Exc}$, $\sum_{\bf k}G_{\bf k}n_{{\rm Pol}|{\bf
k}}\approx G_{{\bf k}={\bf 0}}P$, $\sum_{\bf k}\gamma_{{\rm Pol}|{\bf
k}}n_{{\rm Pol}|{\bf k}}\approx\gamma_{{\rm Pol}|{\bf k}={\bf 0}}P$,
$\sum_{\bf k}G_{\bf k}n_{{\rm Vib}|{\bf k}_{\rm ex}-{\bf k}}\approx G_{{\bf
k}={\bf 0}}V$, $N_{\rm mol}+\sum_{\bf k}1=N_{\rm mol}+N_{\rm states}\approx
N_{\rm mol}$ and denote

$J=\sum_{\bf k}G_{\bf k}\left[n_{{\rm Vib}|{\bf k}_{\rm ex}-{\bf k}}(n_{{\rm
Exc}|{\bf k}_{\rm ex}}-n_{{\rm Pol}|{\bf k}})+N_{\rm mol}n_{\rm
Vib_{D}}(n_{\rm Exc_{D}}-n_{{\rm Pol}|{\bf k}})\right],$ (65) $G_{S}=\sum_{\bf
k}G_{\bf k}.$ (66)

Figure 5: Numerical simulations (a) and analytic results (b) obtained for the
total number of lower polaritons. The dashed black line shows the analytic
approximation to the threshold values given by Eq. (68).

Here we consider a situation when the number of vibrations per molecule is
much less than one.In this case, the term $J/N_{\rm mol}$ is much smaller than
the term $I(PG_{{\bf k}={\bf 0}}+G_{S})/N_{\rm mol}$ in the Eqs. (62)–(64).
Therefore we can neglect this term and obtain the stationary solution for the
total amount of polaritons

$\sum_{\bf k}n_{{\rm Pol}|{\bf k}}={\gamma_{\rm Exc}\over 2\gamma_{{\rm
Pol}|{\bf k}={\bf 0}}}(\varkappa_{\rm Pump}-\varkappa_{\rm Pump}^{\rm
thresh})+\sqrt{\left[{\gamma_{\rm Exc}\over 2\gamma_{{\rm Pol}|{\bf k}={\bf
0}}}(\varkappa_{\rm Pump}-\varkappa_{\rm Pump}^{\rm
thresh})\right]^{2}+\varkappa_{\rm Pump}{\gamma_{\rm Exc}\over\gamma_{{\rm
Pol}|{\bf k}={\bf 0}}}{G_{S}\over G_{{\bf k}={\bf 0}}}},$ (67)

where the threshold pumping rate is defined as

$\varkappa_{\rm Pump}^{\rm thresh}=N_{\rm mol}{\gamma_{{\rm Pol}|{\bf k}={\bf
0}}\over G_{{\bf k}={\bf 0}}}+{G_{S}\over G_{{\bf k}={\bf 0}}}{\gamma_{{\rm
Pol}|{\bf k}={\bf 0}}\over\gamma_{\rm Exc}}$ (68)

For the system considered in the main text the second term of the right-hand
side of the Eq. (68) is much smaller than the first one. Therefore in the main
text we preserve only the first term.

We present the analysis of the Eq. (68) in the main text.

On Figure 5 we compare the numerical simulations for $\sum_{\bf k}n_{{\rm
Pol}|{\bf k}}$ and the analytic result Eq. (67). One can see, that analytic
approximation is in good agreement with the numerical simulations. The values
$\sum_{\bf k}n_{{\rm Pol}|{\bf k}}$ at the threshold exhibit an abrupt change,
a feature that is more pronounced when comparing numerical simulations to
analytics. We believe this discrepancy comes from our analytic approach where
we neglect the term $J/N_{\rm mol}$ in the Eqs. (62)–(64). This term, in
particular, describes the influence of the microscopic occupation of the
bright vibrational states on the BEC dynamics.

At the threshold (68) the total number of the lower polaritons is given by

$\sum_{\bf k}n_{{\rm Pol}|{\bf k}}|_{\rm thresh}=\sqrt{\varkappa_{\rm
Pump}^{\rm thresh}{\gamma_{\rm Exc}\over\gamma_{{\rm Pol}|{\bf k}={\bf
0}}}{G_{S}\over G_{{\bf k}={\bf 0}}}}.$ (69)

The formation of a polariton BEC generally requires that the total number of
lower polaritons surpasses the critical number below

$\sum_{\bf k}n_{{\rm Pol}|{\bf k}}|_{\rm critical}\approx\nu k_{B}T$ (70)

where $\nu$ is the density of states obtained from Eq. (42)

$\nu\approx{S\over 4\pi\alpha_{\rm Cav}}.$ (71)

For the realistic parameters of the system, considered in the main text we get
the following total polariton number

$\sum_{\bf k}n_{{\rm Pol}|{\bf k}}|_{\rm thresh}\sim\sum_{\bf k}n_{{\rm
Pol}|{\bf k}}|_{\rm critical}\sim 10^{2}-10^{3}.$ (72)

Thus, crossing of the threshold (68) is immediately followed by the BEC
formation.

Interestingly, $\sum_{\bf k}n_{{\rm Pol}|{\bf k}}|_{\rm
thresh}\propto\Lambda^{-1}$ and $\sum_{\bf k}n_{{\rm Pol}|{\bf k}}|_{\rm
critical}$ does not depend on $\Lambda$, whereas $\sum_{\bf k}n_{{\rm
Pol}|{\bf k}}|_{\rm critical}\propto T$ and $\sum_{\bf k}n_{{\rm Pol}|{\bf
k}}|_{\rm thresh}$ does no depend on $T$, given that $\Gamma_{\rm Exc}$
remains constant. Therefore, by adjusting the Huang-Rhys factor and the
temperature of the sample, one can distinguish between two regimes: the
polariton laser and the BEC.

##  VI  Coherent anti-Stokes Raman Scattering on polariton BEC: magic angle
conditions

Figure 6: Phase matching condition for the formation of polariton BEC and
macroscopically occupied vibrational states (a), phase matching condition for
the Raman probe (b) and schematic layout of the CARS measurement (c).

Here we calculate the magic angle of the CARS signal induced by Raman
scattering from polariton BEC. We consider the probe beam with the frequency
$\omega_{\rm Probe}$ and the incident angle $\theta_{\rm Probe}$. We denote
the frequency of the pumping beam as $\omega_{\rm Pump}$ and its incident
angle as $\theta_{\rm Pump}$. We use the following assumption $\omega_{\rm
Pump}=\omega_{\rm ex}$. Given the in-plane wavevector of the bright vibrations
${\bf k}_{\rm Vib}$ and in-plane wavevector of BEC ${\bf k}_{\rm BEC}$ we can
express the energy conservation law in the following form

$\omega_{\rm Pump}=\omega_{\rm BEC}+\omega_{\rm Vib}$ (73) $\omega_{\rm
CARS}=\omega_{\rm Probe}+\omega_{\rm Vib}$ (74)

as well as in-plane phase matching condition

$({\bf k}_{\rm Pump})_{||}=({\bf k}_{\rm BEC})_{||}+{\bf k}_{\rm Vib}$ (75)
$({\bf k}_{\rm CARS})_{||}=({\bf k}_{\rm Probe})_{||}+{\bf k}_{\rm Vib}$ (76)

where $||$ denotes an in-plane component of the corresponding wavevector and
$({\bf k}_{\rm Pump})_{||}={\bf k}_{\rm ex}$. Under the assumption that all
vectors $({\bf k}_{\rm Probe})_{||}$, ${\bf k}_{\rm Vib}$ and $({\bf k}_{\rm
Pump})_{||}$ lie in the same plane, one can treat them as c-numbers.
Consequently, the following angles $\theta_{\rm Probe}$, $\theta_{\rm Pump}$
and $\theta_{\rm CARS}$ can be introduced (see Figure 6)

$({\bf k}_{\rm Pump})_{||}={\omega_{\rm Pump}\over c}\sin\theta_{\rm Pump}$
(77) $({\bf k}_{\rm Probe})_{||}={\omega_{\rm Probe}\over c}\sin\theta_{\rm
Probe}$ (78) $({\bf k}_{\rm CARS})_{||}={\omega_{\rm CARS}\over
c}\sin\theta_{\rm CARS}$ (79)

The bright vibrational state with maximal occupation corresponds to the
wavevector ${\bf k}_{\rm Vib}={\bf k}_{\rm ex}=({\bf k}_{\rm Pump})_{||}$, in
agreement with Eq. (75) where ${\bf k}_{\rm BEC}\approx{\bf 0}$. In the
specific configuration of the polariton BEC we denote the corresponding angle
$\theta_{\rm CARS}$ as $\theta_{\rm MA}$ and refer to it as the “magic angle”.
Based on Eqs. (73)–(79) and condition ${\bf k}_{\rm BEC}={\bf 0}$ we derive
the expression for $\theta_{\rm MA}$ as presented in Eq. (9) of the main text.