# Organic polariton lasing and the weak- to strong-coupling crossover

Artem Strashko  SUPA, School of Physics and Astronomy, University of St
Andrews, St Andrews, KY16 9SS, United Kingdom  Peter Kirton  Vienna Center for
Quantum Science and Technology, Atominstitut, TU Wien, 1040 Vienna, Austria
Jonathan Keeling  SUPA, School of Physics and Astronomy, University of St
Andrews, St Andrews, KY16 9SS, United Kingdom

(August 27, 2024)

###### Abstract

Following experimental realizations of room temperature polariton lasing with
organic molecules, we present a microscopic model that allows us to explore
the crossover from weak to strong matter-light coupling. We consider a non-
equilibrium Dicke-Holstein model, including both strong coupling to
vibrational modes and strong matter-light coupling, providing the phase
diagram of this model in the thermodynamic limit. We discuss the mechanism of
polariton lasing, uncovering a process of self-tuning, and identify the
relation and distinction between regular dye lasers and organic polariton
lasers.

Bose-Einstein statistics underpin both the thermal equilibrium phenomenon of
Bose-Einstein condensation, and the non-equilibrium phenomenon of lasing.
Lying between these two extremes, there now exist several experimental
platforms; in particular exciton-polaritons (quasiparticles resulting from
strong coupling between photons and excitons) in semiconductor microcavities
at cryogenic temperatures Deng  _et al._ (2002); Kasprzak  _et al._ (2006);
Balili  _et al._ (2007), and photons in dye-filled microcavities at room
temperature Klaers  _et al._ (2010). Since these microcavities are imperfect,
they are sources of coherent light (as is a laser), but differ in mechanism
from photon lasing Deng  _et al._ (2010); Carusotto and Ciuti (2013). Indeed,
polariton lasing does not need electronic inversion, and so it is often stated
that it can provide coherent light sources with ultra-low thresholds.
Polariton lasing can also occur at room temperature in appropriate materials:
inorganic materials such as wide bandgap semiconductors Christopoulos  _et
al._ (2007); Guillet  _et al._ (2011); Li  _et al._ (2013) and two-dimensional
materials Waldherr  _et al._ (2018), and the focus of this Letter, organic
materials Kéna-Cohen and Forrest (2010); Plumhof  _et al._ (2014); Daskalakis
_et al._ (2014).

Polariton condensation in organic materials prompts interesting questions
regarding the mechanisms of polariton relaxation and lasing. Excitons in
organic materials are Frenkel excitons — electronic excitations of a molecule
or chromaphore delocalized by hopping Agranovich (2009). Excitons in organic
materials typically show complex absorption and emission spectra, due to
strong coupling between the electronic state and the nuclear configuration,
leading to rovibrational dressing. This causes a Stokes shift, so that
emission is at longer wavelengths than absorption. Spectral separation of
emission and absorption underpins the operation of dye lasers Schäfer (1990),
allowing gain without electronic inversion. Since both strong matter-light
coupling and large Stokes shifts are expected to reduce the lasing threshold,
how they act in concert is of both fundamental interest and practical
relevance.

Theoretical modeling of polariton condensates can follow a number of
approaches. To describe the macroscopic pattern formation and superfluid
hydrodynamics it mostly suffices to use the phenomenological complex Gross-
Pitaevskii equation Aranson and Kramer (2002); Carusotto and Ciuti (2013).
However, such order parameter equations are ubiquitous in non-equilibrium
systems breaking $U(1)$ symmetry, so similar equations also apply for a photon
laser Staliunas and Sanchez-Morcillo (2003); Berloff and Keeling (2013). Such
approaches are thus not well suited to understanding the relation between
polariton and photon lasing, or the particular properties of organic
polaritons. To answer such questions, a more relevant approach is to use
kinetic equations for the population of polaritons and excitonic reservoir,
with decay rates accounting for vibronically assisted processes Michetti and
La Rocca (2009); Fontanesi and La Rocca (2009); Mazza and La Rocca (2009);
Michetti  _et al._ (2015) relevant in organic materials. This approach has
been used to understand the onset of lasing in anthracene microcavities Mazza
_et al._ (2013) and J-aggregated dyes Coles  _et al._ (2011). It however
assumes polariton modes are well defined, so cannot access the weak- to
strong-coupling crossover. In this Letter we will work from a microscopic
Hamiltonian which can, in the appropriate limit, recover both the physics
described in such a kinetic model, as well as that of a dye laser in weak
coupling.

We consider the non-equilibrium Dicke-Holstein model, presented below, which
describes many molecules (with vibrationally dressed electronic transitions)
coupled to a common photon mode. The equilibrium phase diagram of this model
has been presented elsewhere Ćwik  _et al._ (2014). Similar models have been
used to understand the absorption and photoluminescence spectra of organic
polaritons Ćwik  _et al._ (2016); Herrera and Spano (2017a, b); Zeb  _et al._
(2017), and to explore whether strong matter-light coupling affects chemical
reactions Herrera and Spano (2016); Feist  _et al._ (2017); Ribeiro  _et al._
(2018). Despite this, the non-equilibrium dynamics of this model, when both
photons and vibrations couple strongly to the electronic state, has not been
explored. If coupling to vibrations is weak, one can eliminate these using
Bloch-Redfield-Wangsness theory Wangsness and Bloch (1953); Redfield (1955) to
capture the relaxation dynamics. For weak coupling to light, rate equations
are available Kirton and Keeling (2013, 2015); Hesten  _et al._ (2018);
Radonjić  _et al._ (2018), describing photon condensation. For a few
molecules, one can use exact numerical methods, treating the vibrational modes
as a non-Markovian dissipation process del Pino  _et al._ (2018a, b). However,
the thermodynamic limit with many molecules and strong coupling remains a
challenge, which we address here.

In this Letter we present the phase diagram of the non-equilibrium Dicke-
Holstein model, exploring the crossover from weak to strong matter-light
coupling. We uncover how the mechanism underlying polariton lasing evolves
with coupling strength. In the weak coupling limit, we recover the results of
models of photon condensation Kirton and Keeling (2013, 2015), while for
strong coupling our results are consistent with kinetic models Mazza and La
Rocca (2009). Finally, we address the practically significant question of
whether strong coupling reduces lasing threshold, and the optimal parameters
to realize polariton lasing.

Figure 1: (a) Cartoon illustrating our model: many molecules ($N$-level
systems) are coupled to a cavity mode. (b-d) Weak coupling behavior. (b)
Emission and absorption spectra of the molecules. (c) Weak coupling phase
diagram. The cyan dashed line marks the phase boundary without vibrational
dressing ($S=0$). (d) Dominant molecular transitions coupled into the lasing
mode at threshold. All energies are measured in units of $\varepsilon\equiv
1$, other parameters are $S=0.1$, $\omega_{v}=0.2$,
$\Gamma_{\downarrow}=\kappa=10^{-4}$, $\Gamma_{z}=0.03$, $\gamma_{v}=0.02$,
$k_{B}T_{v}=0.025$, $N_{v}=4$.

We consider a single cavity mode, with frequency $\omega_{c}$, coupled to
$\mathcal{N}_{m}$ organic molecules as illustrated in Fig. 1(a). We model each
molecule as two electronic states (HOMO and LUMO levels), dressed by a bosonic
mode describing an internal molecular vibrational mode. This yields the Dicke-
Holstein model:

$H=\sum_{n}\left\\{\varepsilon\sigma_{n}^{z}+\omega_{v}\left[b^{\dagger}_{n}b_{n}^{\mathstrut}+\sqrt{S}\sigma_{n}^{z}(b^{\dagger}_{n}+b_{n}^{\mathstrut})\right]\right\\}\\\
+\omega_{c}a^{\dagger}a+g\left(a^{\dagger}+a\right)\sum_{n}\sigma_{n}^{x}+\frac{g^{2}\mathcal{N}_{m}}{\varepsilon}(a^{\dagger}+a)^{2},$
(1)

where $\sigma_{n}^{\alpha=x,y,z}$ are Pauli matrices describing the electronic
state of molecule $n$, $b^{\dagger}_{n}$ creates a vibrational excitation on
molecule $n$, and $a^{\dagger}$ creates a photon. To explore the ultrastrong
matter-light coupling regime we do not make a rotating-wave approximation and
we include a diamagnetic $A^{2}$ term, which prevents a ground state
superradiant transition (see Hepp and Lieb (1973); Rzażewski  _et al._ (1975);
Grießer  _et al._ (2016); De Bernardis  _et al._ (2018); Sup ). We note the
diamagnetic term could be eliminated by a Bogoliubov transform, yielding a
model without the $A^{2}$ term, but with a modified photon frequency
${\omega_{c}^{\text{eff}}=\sqrt{\omega_{c}\left(\omega_{c}+4g^{2}\mathcal{N}_{m}/\varepsilon\right)}}$,
and matter-light coupling
$g^{\text{eff}}=g\sqrt{\omega_{c}/\omega^{\text{eff}}_{c}}$. We characterize
matter-light coupling by the “bare” polariton splitting
$g\sqrt{\mathcal{N}_{m}}$, which would be the splitting between the upper and
lower polariton in the limit $S\to 0$.

To include incoherent pumping and decay processes, we use a Lindblad master
equation Breuer and Petruccione (2002) of the form:

$\dot{\rho}(t)=-i\left[H,\rho\right]+\kappa\mathcal{L}[a]+\sum_{n}\Big{(}\Gamma_{\downarrow}\mathcal{L}[\sigma_{n}^{-}]+\Gamma_{\uparrow}\mathcal{L}[\sigma_{n}^{+}]\\\
+\Gamma_{z}\mathcal{L}[\sigma_{n}^{z}]+\gamma_{\uparrow}\mathcal{L}[b^{\dagger}_{n}-\sqrt{S}\sigma^{z}]+\gamma_{\downarrow}\mathcal{L}[b^{\mathstrut}_{n}-\sqrt{S}\sigma^{z}]\Big{)},$
(2)

where $\mathcal{L}[X]=X\rho X^{\dagger}-\\{X^{\dagger}X,\rho\\}/2$. We include
electronic excitation and decay via non-cavity modes with rates
$\Gamma_{\uparrow},\Gamma_{\downarrow}$ respectively, dephasing with rate
$\Gamma_{z}$, and photon loss with rate $\kappa$. The final two terms describe
relaxation of the vibrational mode to thermal equilibrium at temperature
$T_{v}$, accounting for the electronic-state-dependent vibrational
displacement. These rates are thus $\gamma_{\uparrow}=\gamma_{v}n_{B}$,
$\gamma_{\downarrow}=\gamma_{v}(n_{B}+1)$ where
$n_{B}=[\exp(\omega_{v}/k_{B}T_{v})-1]^{-1}$. Throughout this Letter we
measure energies in units of the electronic transition energy, so
$\varepsilon\equiv 1$ by definition — typical physical values are
$\varepsilon\simeq 1$–$2$eV — other parameters are given in the caption of
Fig. 1. We choose parameters that are physical and demonstrate sidebands in
the molecular spectra, Fig. 1(b). The resulting weak coupling phase diagram,
Fig. 1(c), is straightforward to understand. The threshold
$\Gamma_{\uparrow}^{th}$ is reduced (i.e.
$\Gamma_{\uparrow}^{th}<\Gamma_{\downarrow}$) when the cavity frequency
matches the emission peak, and $\Gamma_{\uparrow}^{th}$ increases when the
cavity matches the absorption peak.

Our aim is to consider the large $\mathcal{N}_{m}$ limit of this model
exactly, allowing both for strong matter-light coupling and strong coupling
between vibrational and electronic states. In the large $\mathcal{N}_{m}$
limit, a mean-field (i.e. Maxwell-Bloch) treatment of the Dicke model becomes
exact Haken (1970); Kirton and Keeling (2017). There is however no decoupling
between electronic and vibrational degrees of freedom. Instead we define
generalized molecular operators, describing an $N$-level system as a whole,
with $N=2N_{v}$ for $N_{v}$ vibrational levels, where the choice of $N_{v}$
depends on the value of $S$. We write the operators for this $N$-level system
using a basis of generalized Gell Mann matrices $\lambda_{i}$ Stone and
Goldbart (2009), satisfying $\text{Tr}(\lambda_{i}\lambda_{j})=2\delta_{ij}$.
This enables us to write any operator as
$O=(\lambda_{i}/2)\text{Tr}\left(O\lambda_{i}\right)$ 111In writing this, we
have defined $\lambda_{0}=(2/N)\mathbf{1}_{N}$.. We may then write the
Hamiltonian as:

$H=\omega_{c}a^{\dagger}a+\sum_{n}\left[A_{i}+B_{i}(a^{\dagger}+a)\right]\lambda_{i}^{(n)}+\frac{g^{2}\mathcal{N}_{m}}{\varepsilon}(a^{\dagger}+a)^{2},$
(3)

with summation convention over $i$, where $A_{i},B_{i}$ can be found by
constructing the molecular operators in Eq. (1) as $N\times N$ matrices, and
taking traces. The Lindblad master equation can be rewritten in the same way:

$\dot{\rho}(t)=-i\left[H,\rho\right]+\kappa\mathcal{L}[a]+\sum_{\mu
n}\mathcal{L}\left[c_{i}^{\mu}\lambda_{i}^{(n)}\right].$ (4)

The sum over $\mu$ is over the five molecular dissipative channels in Eq. (2).
If these terms are written as $\sum_{\mu}\Gamma_{\mu}\mathcal{L}[J_{\mu}]$
then we have $c_{i}^{\mu}=\sqrt{\Gamma_{\mu}}\text{Tr}(J_{\mu}\lambda_{i})/2$.

The mean-field decoupling is realized by deriving the equations of motion for
variables $\alpha=\langle a\rangle,\ell_{i}=\langle\lambda_{i}\rangle$, where
$\langle O\rangle=\operatorname{Tr}(\rho O)$, and then making the mean-field
decoupling: $\langle a\lambda_{i}\rangle=\langle
a\rangle\langle\lambda_{i}\rangle$. This leads to the set of nonlinear coupled
differential equations:

$\displaystyle\partial_{t}\alpha$
$\displaystyle=-\left(i\omega_{c}+\frac{\kappa}{2}\right)\alpha-4i\frac{g^{2}\mathcal{N}_{m}}{\varepsilon}\mathrm{Re}[\alpha]-i\mathcal{N}_{m}B_{i}\ell_{i},$
(5) $\displaystyle\partial_{t}\ell_{i}$
$\displaystyle=\Big{(}\xi_{ik}+4f_{ijk}B_{j}\text{Re}[\alpha]\Big{)}\ell_{k}+\frac{4i}{N}c_{j}^{\mu}c_{k}^{\mu\ast}f_{ijk},$
(6)

where
$\xi_{ik}=2f_{ijk}A_{j}+ic_{l}^{\mu}c_{m}^{\mu\ast}(f_{ilp}\zeta_{mpk}+f_{mip}\zeta_{plk})$
with $\zeta_{ijk}\equiv\text{Tr}(\lambda_{i}\lambda_{j}\lambda_{k})/2$, and
$f_{ijk}\equiv\text{Tr}([\lambda_{i},\lambda_{j}]\lambda_{k})/4i$.

Mean-field theory shows a phase transition between a normal state with
$\alpha=0$, and a state with $\alpha\neq 0$ which we denote as a laser but
which may be either photon or polariton lasing. The phase boundary can be
found by considering when fluctuations about the normal state are unstable. We
thus write linearized equations of motion for fluctuations
$\alpha=\delta\alpha$, $\ell_{i}=\ell_{i,\text{ns}}+\delta\ell_{i}$, where
$\ell_{i,\text{ns}}$ is the normal state solution. Defining
$\upsilon=\left(\delta\alpha,\delta\alpha^{*},\mathbf{\delta\ell}\right)^{\intercal}$
we can write ${\partial_{t}\upsilon=\mathcal{M}\upsilon}$, and find the
eigenmodes, $\mathcal{M}v^{k}=\xi^{k}v^{k}$. The real part of $\xi^{k}$ gives
the growth (positive) or decay (negative) rate of a mode, while the imaginary
part gives its oscillation frequency.

From the eigenvector $\upsilon^{k}$ we can also find the contributions of
different molecular transitions to the given unstable mode (see Sup  for
details.) As an example of this, Fig. 1(d) shows the composition of the
unstable mode precisely at the lasing threshold, i.e. along the phase boundary
shown in Fig. 1(c). We see that where the threshold is low, the (1–0)
transition contributes most, while where the threshold is high, the (0–1)
transition dominates. In this Letter we show results with up to three
vibrational excitations, i.e. $N=8$. Our results are converged for this choice
as shown in Sup .

Using the methods above, we now explore how strong coupling to light modifies
the phase diagram, and understand the physics responsible for this
modification. Figure 2 shows the evolution of the
$\omega_{c},\Gamma_{\uparrow}$ phase diagram as matter-light coupling
increases, focusing on the low pumping regime. At moderate coupling,
$g\sqrt{\mathcal{N}_{m}}=0.1$, we see the same form as for weak coupling, with
a minimum (maximum) threshold at $\omega_{c}\simeq\varepsilon\mp\omega_{v}$.
As $g\sqrt{\mathcal{N}_{m}}$ increases, the most striking feature is that the
lobe at $\omega_{c}\approx\varepsilon-\omega_{v}$ bends and, at weak pumping,
extends to significantly higher cavity frequencies. This eventually leads to a
re-entrant phase diagram.

Figure 2: Evolution of phase diagrams with increasing coupling,
$g\sqrt{\mathcal{N}_{m}}$ (values as shown). Dash-dotted (yellow) lines
indicate the cuts shown in Fig. 3. Parameters as in Fig. 1. Figure 3: Nature
of lasing instability. Top (a–c,g–i): Real (right, solid cyan) and imaginary
(left, grayscale) parts of linearized eigenvalues $\xi$. The grayscale
indicates the photon component of that mode. The yellow dashed line highlights
which imaginary part corresponds to the mode that is unstable (i.e has a
positive real part). The pink dotted line shows the effective photon
frequency, $\omega_{c}^{\text{eff}}$. Bottom (d–f, j–l): Vibrational
composition of unstable mode. Parameters as in Fig. 1.

To understand the origin of this extension of the lobe to high cavity
frequencies, we explore the nature of the unstable mode. Figure 3 shows the
composition of the unstable mode, and the frequency of all modes for three
cuts across the phase diagram at fixed pump power at weaker and stronger
coupling. For orientation, we first summarize the more straightforward
behavior seen at weaker coupling, $g\sqrt{\mathcal{N}_{m}}=0.1$. In Fig. 3(a),
for low pumping $\Gamma_{\uparrow}=0.4\Gamma_{\downarrow}$, the frequencies of
the normal modes show a small polaritonic splitting where the effective photon
frequency 222i.e. this line plots the quantity
$\sqrt{\omega_{c}\left(\omega_{c}+4g^{2}\mathcal{N}_{m}/\varepsilon\right)}$,
hence its nonlinear dependence on $\omega_{c}$ at strong coupling crosses
molecular transitions, most clearly the zero phonon line at
$\text{Im}[\xi]=\varepsilon\equiv 1$, i.e the (0–0) transition. Lasing occurs
when the bare photon frequency is close to the (1-0) or (2-0) transitions (but
with different strengths for these two transitions). Indeed, as seen in Fig.
3(d), the unstable mode predominantly involves the (1–0) transition, crossing
over to the (0–0) transition as $\omega_{c}$ approaches the zero phonon line.
As pumping increases to $\Gamma_{\uparrow}=\Gamma_{\downarrow}$, Fig. 3(b,e),
saturation suppresses the polariton splitting, and lasing is possible over a
wider range of frequencies. Finally, at
$\Gamma_{\uparrow}=1.8\Gamma_{\downarrow}$, Fig. 3(c,f), we have electronic
inversion, and can achieve lasing even when absorption exceeds emission. Here
we have two lasing regions at low and high frequencies. At yet higher pump
strengths, these regions join to form a single region.

At stronger coupling, $g\sqrt{\mathcal{N}_{m}}=1.0$, the picture changes
dramatically. For weak pumping, Fig. 3(g) now shows a much clearer
anticrossing, but also shows a new feature: locking between the frequencies of
the polariton and of the (1–0) vibrational sideband. This frequency locking
persists over the range of cavity frequencies for which lasing occurs. As
such, although the bare photon frequency is high, the lower polariton mode is
at a lower frequency, and is self-tuned to allow feeding by the (1–0)
molecular transition Mazza  _et al._ (2013). Thus, at strong coupling, the
system can self-tune to support such feeding.

As pumping increases to $\Gamma_{\uparrow}=\Gamma_{\downarrow}$, the
polaritonic self-tuning effect reduces and we no longer have lasing over such
a wide range of cavity frequencies. As noted above, stronger pumping
suppresses the polaritonic splitting, and this prevents the self-tuning
effect. We nonetheless still have notable polaritonic splitting at low cavity
frequencies. The molecular transitions involved in lasing in Fig. 3(k) is
similar to that in Fig. 3(e), but now (0–1) and (1–1) transitions play a
larger role, as strong coupling admixes more electronic transitions into the
polariton. As we further increase pumping to
$\Gamma_{\uparrow}=1.8\Gamma_{\downarrow}$ we now have two distinct lasing
regimes: at low cavity frequency, there is polariton lasing, but at high
cavity frequencies, lasing is at the bare cavity photon frequency (see Sup
for an extended figure).

Having analyzed the structure of the phase diagram, we next consider whether
strong matter-light coupling is a direct route to reduce the lasing threshold.
In Fig. 4(a–c) we plot phase diagrams vs $\Gamma_{\uparrow}$ and
$g\sqrt{\mathcal{N}_{m}}$ at three different bare photon frequencies, and we
compare these to the phase boundary predicted by weak coupling theory, see Sup
for details. In all cases, we see the threshold reduce with increasing
coupling. However, for cavity frequencies near the (1-0) transition, Fig.
4(a), the weak coupling prediction of the phase boundary matches the full
results well. For this optimal frequency, while threshold pumping does reduce
with increasing coupling, the threshold saturates as the system enters the
strong coupling regime (where $g\sqrt{\mathcal{N}_{m}}$ is larger than the
linewidth $\Gamma_{z}$). However, for cavity frequencies near the absorption
peak at the (0-1) transition, Fig. 4(c), one sees the full calculation
predicts a dramatically lower threshold than is predicted by weak coupling,
and the threshold continues to reduce even after entering the strong coupling
regime.

Figure 4: (a-c) Phase diagrams at various cavity frequencies. Red dashed line
shows the weak-coupling theory phase boundary, see Sup . (d) Minimal critical
pump strength (optimized over cavity frequency) as a function of matter-light
coupling. Red shading indicates weak- to strong-coupling crossover.

In Fig. 4(d) we summarize the data by plotting the minimal pumping strength
(optimized over cavity frequency Sup ) required for lasing at a given matter-
light coupling strength. While this optimized pump strength changes little
above $g\sqrt{\mathcal{N}_{m}}=0.1$, the results of Fig. 2 clearly show that
the range of cavity frequencies compatible with low threshold lasing does
increase significantly.

In conclusion, we have studied the phase diagram of the non-equilibrium Dicke-
Holstein model. We showed that at strong coupling, self-tuning of the
polariton to optimize feeding by a vibrational sideband leads to lasing over a
wide range of cavity frequencies. We also see that the minimum achievable
threshold is reduced as one approaches strong coupling, but does not change as
coupling is further increased. Our results and approach also open up a number
of possible future directions to explore. For example, one may go beyond mean-
field description with a cumulant expansion Meiser  _et al._ (2009); Meiser
and Holland (2010); Kirton and Keeling (2017, 2018) and explore the spectral
properties of emission and also luminescence below lasing threshold. Further,
by considering a multimode Dicke model, one can ask about the thermalization
of mode populations, analogous to that seen in the weak-coupling photon BEC
Klaers  _et al._ (2010); Kirton and Keeling (2013, 2015); Marelic and Nyman
(2015); Hesten  _et al._ (2018); Radonjić  _et al._ (2018). Finally, one may
consider more complex models, for example combining multiple vibrational modes
in the system, or using Redfield theory with the true system eigenstates to
better capture the physics of structured baths, leading to a more realistic
master equation Ciuti  _et al._ (2005); del Pino  _et al._ (2015). This may
allow one to capture physics that otherwise requires numerically intensive
non-Markovian simulations del Pino  _et al._ (2018a, b). Taken together, these
provide an efficient route to understand the full range of weak and strong
matter-light coupling for polariton condensation and lasing.

###### Acknowledgements.

AS acknowledges support from the EPSRC CM-CDT (EP/L015110/1). JK and PK
acknowledges financial support from EPSRC program ”Hybrid Polaritonics”
(EP/M025330/1). PK acknowledges support from the Austrian Academy of Sciences
(ÖAW).

## References

  * Deng  _et al._ (2002) H. Deng, G. Weihs, C. Santori, J. Bloch, and Y. Yamamoto,  _Condensation of Semiconductor Microcavity Exciton Polaritons_ , Science 298, 199 (2002).
  * Kasprzak  _et al._ (2006) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. Keeling, F. Marchetti, M. Szymańska, R. Andre, J. Staehli, _et al._ ,  _Bose–Einstein condensation of exciton polaritons_ , Nature 443, 409 (2006).
  * Balili  _et al._ (2007) R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West,  _Bose-Einstein Condensation of Microcavity Polaritons in a Trap_ , Science 316, 1007 (2007).
  * Klaers  _et al._ (2010) J. Klaers, J. Schmitt, F. Vewinger, and M. Weitz,  _Bose–Einstein condensation of photons in an optical microcavity_ , Nature 468, 545 (2010).
  * Deng  _et al._ (2010) H. Deng, H. Haug, and Y. Yamamoto,  _Exciton-polariton Bose-Einstein condensation_ , Rev. Mod. Phys. 82, 1489 (2010).
  * Carusotto and Ciuti (2013) I. Carusotto and C. Ciuti,  _Quantum fluids of light_ , Rev. Mod. Phys. 85, 299 (2013).
  * Christopoulos  _et al._ (2007) S. Christopoulos, G. B. H. von Högersthal, A. J. D. Grundy, P. G. Lagoudakis, A. V. Kavokin, J. J. Baumberg, G. Christmann, R. Butté, E. Feltin, J.-F. Carlin, and N. Grandjean,  _Room-Temperature Polariton Lasing in Semiconductor Microcavities_ , Phys. Rev. Lett. 98, 126405 (2007).
  * Guillet  _et al._ (2011) T. Guillet, M. Mexis, J. Levrat, G. Rossbach, C. Brimont, T. Bretagnon, B. Gil, R. Butté, N. Grandjean, L. Orosz, _et al._ ,  _Polariton lasing in a hybrid bulk ZnO microcavity_ , Appl. Phys. Lett. 99, 161104 (2011).
  * Li  _et al._ (2013) F. Li, L. Orosz, O. Kamoun, S. Bouchoule, C. Brimont, P. Disseix, T. Guillet, X. Lafosse, M. Leroux, J. Leymarie, M. Mexis, M. Mihailovic, G. Patriarche, F. Réveret, D. Solnyshkov, J. Zuniga-Perez, and G. Malpuech,  _From Excitonic to Photonic Polariton Condensate in a ZnO-Based Microcavity_ , Phys. Rev. Lett. 110, 196406 (2013).
  * Waldherr  _et al._ (2018) M. Waldherr, N. Lundt, M. Klaas, S. Betzold, M. Wurdack, V. Baumann, E. Estrecho, A. Nalitov, E. Cherotchenko, H. Cai, E. A. Ostrovskaya, A. V. Kavokin, S. Tongay, S. Klembt, S. Höfling, and C. Schneider,  _Observation of bosonic condensation in a hybrid monolayer MoSe2-GaAs microcavity_ , (2018), 1805.03630 .
  * Kéna-Cohen and Forrest (2010) S. Kéna-Cohen and S. Forrest,  _Room-temperature polariton lasing in an organic single-crystal microcavity_ , Nat. Photon. 4, 371 (2010).
  * Plumhof  _et al._ (2014) J. D. Plumhof, T. Stöferle, L. Mai, U. Scherf, and R. F. Mahrt,  _Room-temperature Bose–Einstein condensation of cavity exciton–polaritons in a polymer_ , Nature materials 13, 247 (2014).
  * Daskalakis  _et al._ (2014) K. Daskalakis, S. Maier, R. Murray, S. Kéna-Cohen, _et al._ ,  _Nonlinear interactions in an organic polariton condensate_ , Nat. Mater 13, 271 (2014).
  * Agranovich (2009) V. M. Agranovich,  _Excitations in organic solids_ (Oxford University Press, Oxford, 2009).
  * Schäfer (1990) F. P. Schäfer, ed.,  _Dye Lasers_ , 3rd ed. (Springer-Verlag, Berlin, 1990).
  * Aranson and Kramer (2002) I. S. Aranson and L. Kramer,  _The world of the complex Ginzburg-Landau equation_ , Rev. Mod. Phys. 74, 99 (2002).
  * Staliunas and Sanchez-Morcillo (2003) K. Staliunas and V. J. Sanchez-Morcillo,  _Transverse patterns in nonlinear optical resonators_ , Vol. 183 (Springer-Verlag, Berlin, 2003).
  * Berloff and Keeling (2013) N. G. Berloff and J. Keeling,  _Universality in Modelling Non-equilibrium Pattern Formation in Polariton Condensates_ , in  _Physics of Quantum Fluids: New Trends and Hot Topics in Atomic and Polariton Condensates_, edited by A. Bramati and M. Modugno (Springer Berlin Heidelberg, Berlin, Heidelberg, 2013) pp. 19–38.
  * Michetti and La Rocca (2009) P. Michetti and G. C. La Rocca,  _Exciton-phonon scattering and photoexcitation dynamics in J-aggregate microcavities_ , Phys. Rev. B 79, 035325 (2009).
  * Fontanesi and La Rocca (2009) L. Fontanesi and G. C. La Rocca,  _Organic-based microcavities with vibronic progressions: Linear spectroscopy_ , Phys. Rev. B 80, 235313 (2009).
  * Mazza and La Rocca (2009) L. Mazza and G. C. La Rocca,  _Organic-based microcavities with vibronic progressions: Photoluminescence_ , Phys. Rev. B 80, 235314 (2009).
  * Michetti  _et al._ (2015) P. Michetti, L. Mazza, and G. C. La Rocca,  _Organic Nanophotonics, Nano-Optics and Nanophotonics_, edited by Y. S. Zhao (Springer, Berlin, 2015).
  * Mazza  _et al._ (2013) L. Mazza, S. Kéna-Cohen, P. Michetti, and G. C. La Rocca,  _Microscopic theory of polariton lasing via vibronically assisted scattering_ , Phys. Rev. B 88, 075321 (2013).
  * Coles  _et al._ (2011) D. M. Coles, P. Michetti, C. Clark, W. C. Tsoi, A. M. Adawi, J.-S. Kim, and D. G. Lidzey,  _Vibrationally Assisted Polariton-Relaxation Processes in Strongly Coupled Organic-Semiconductor Microcavities_ , Adv. Funct. Mater. 21, 3691 (2011).
  * Ćwik  _et al._ (2014) J. A. Ćwik, S. Reja, P. B. Littlewood, and J. Keeling,  _Polariton condensation with saturable molecules dressed by vibrational modes_ , Europhys. Lett. 105, 47009 (2014).
  * Ćwik  _et al._ (2016) J. A. Ćwik, P. Kirton, S. De Liberato, and J. Keeling,  _Excitonic spectral features in strongly coupled organic polaritons_ , Phys. Rev. A 93, 033840 (2016).
  * Herrera and Spano (2017a) F. Herrera and F. C. Spano,  _Dark Vibronic Polaritons and the Spectroscopy of Organic Microcavities_ , Phys. Rev. Lett. 118, 223601 (2017a).
  * Herrera and Spano (2017b) F. Herrera and F. C. Spano,  _Theory of nanoscale organic cavities: The essential role of vibration-photon dressed states_ , ACS Photonics 5, 65 (2017b).
  * Zeb  _et al._ (2017) M. A. Zeb, P. G. Kirton, and J. Keeling,  _Exact states and spectra of vibrationally dressed polaritons_ , ACS Photonics 5, 249 (2017).
  * Herrera and Spano (2016) F. Herrera and F. C. Spano,  _Cavity-Controlled Chemistry in Molecular Ensembles_ , Phys. Rev. Lett. 116, 238301 (2016).
  * Feist  _et al._ (2017) J. Feist, J. Galego, and F. J. Garcia-Vidal,  _Polaritonic chemistry with organic molecules_ , ACS Photonics 5, 205 (2017).
  * Ribeiro  _et al._ (2018) R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Campos-Gonzalez-Angulo, and J. Y. Zhou,  _Polariton Chemistry: controlling molecular dynamics with optical cavities_ , Chem. Sci. 9, 6325 (2018).
  * Wangsness and Bloch (1953) R. K. Wangsness and F. Bloch,  _The Dynamical Theory of Nuclear Induction_ , Phys. Rev. 89, 728 (1953).
  * Redfield (1955) A. G. Redfield,  _Nuclear Magnetic Resonance Saturation and Rotary Saturation in Solids_ , Phys. Rev. 98, 1787 (1955).
  * Kirton and Keeling (2013) P. Kirton and J. Keeling,  _Nonequilibrium Model of Photon Condensation_ , Phys. Rev. Lett. 111, 100404 (2013).
  * Kirton and Keeling (2015) P. Kirton and J. Keeling,  _Thermalization and breakdown of thermalization in photon condensates_ , Phys. Rev. A 91, 033826 (2015).
  * Hesten  _et al._ (2018) H. J. Hesten, R. A. Nyman, and F. Mintert,  _Decondensation in Nonequilibrium Photonic Condensates: When Less Is More_ , Phys. Rev. Lett. 120, 040601 (2018).
  * Radonjić  _et al._ (2018) M. Radonjić, W. Kopylov, A. Balaž, and A. Pelster,  _Interplay of coherent and dissipative dynamics in condensates of light_ , New J. Phys. 20, 055014 (2018).
  * del Pino  _et al._ (2018a) J. del Pino, F. A. Schröder, A. W. Chin, J. Feist, and F. J. Garcia-Vidal,  _Tensor network simulation of non-Markovian dynamics in organic polaritons_ , (2018a), 1804.04511 .
  * del Pino  _et al._ (2018b) J. del Pino, F. A. Schröder, A. W. Chin, J. Feist, and F. J. Garcia-Vidal,  _Tensor Network simulation of polaron-polaritons in organic microcavities_ , (2018b), 1807.00586 .
  * Hepp and Lieb (1973) K. Hepp and E. H. Lieb,  _Equilibrium statistical mechanics of matter interacting with the quantized radiation field_ , Phys. Rev. A 8, 2517 (1973).
  * Rzażewski  _et al._ (1975) K. Rzażewski, K. Wódkiewicz, and W. Żakowicz,  _Phase transitions, two-level atoms, and the A 2 term_ , Phys. Rev. Lett. 35, 432 (1975).
  * Grießer  _et al._ (2016) T. Grießer, A. Vukics, and P. Domokos,  _Depolarization shift of the superradiant phase transition_ , Phys. Rev. A 94, 033815 (2016).
  * De Bernardis  _et al._ (2018) D. De Bernardis, T. Jaako, and P. Rabl,  _Cavity quantum electrodynamics in the nonperturbative regime_ , Phys. Rev. A 97, 043820 (2018).
  * (45) See Supplemental Material for discussion of the $A^{2}$ term, further details of the linear stability calculation, convergence with vibrational states, and further analysis of the phase diagrams presented.
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2002).
  * Haken (1970) H. Haken,  _The semiclassical and quantum theory of the laser_ , in  _Quantum Opt._ , edited by S. M. Kay and A. Maitland (1970) p. 201.
  * Kirton and Keeling (2017) P. Kirton and J. Keeling,  _Suppressing and Restoring the Dicke Superradiance Transition by Dephasing and Decay_ , Phys. Rev. Lett. 118, 123602 (2017).
  * Stone and Goldbart (2009) M. Stone and P. Goldbart,  _Mathematics for physics: a guided tour for graduate students_ (Cambridge University Press, Cambridge, 2009).
  * Note (1) In writing this, we have defined $\lambda_{0}=(2/N)\mathbf{1}_{N}$.
  * Note (2) i.e. this line plots the quantity $\sqrt{\omega_{c}\left(\omega_{c}+4g^{2}\mathcal{N}_{m}/\varepsilon\right)}$, hence its nonlinear dependence on $\omega_{c}$ at strong coupling.
  * Meiser  _et al._ (2009) D. Meiser, J. Ye, D. R. Carlson, and M. J. Holland,  _Prospects for a Millihertz-Linewidth Laser_ , Phys. Rev. Lett. 102, 163601 (2009).
  * Meiser and Holland (2010) D. Meiser and M. J. Holland,  _Intensity fluctuations in steady-state superradiance_ , Phys. Rev. A 81, 063827 (2010).
  * Kirton and Keeling (2018) P. Kirton and J. Keeling,  _Superradiant and lasing states in driven-dissipative Dicke models_ , New J. Phys. 20, 015009 (2018).
  * Marelic and Nyman (2015) J. Marelic and R. A. Nyman,  _Experimental evidence for inhomogeneous pumping and energy-dependent effects in photon Bose-Einstein condensation_ , Phys. Rev. A 91, 033813 (2015).
  * Ciuti  _et al._ (2005) C. Ciuti, G. Bastard, and I. Carusotto,  _Quantum vacuum properties of the intersubband cavity polariton field_ , Phys. Rev. B 72, 115303 (2005).
  * del Pino  _et al._ (2015) J. del Pino, J. Feist, and F. J. Garcia-Vidal,  _Quantum theory of collective strong coupling of molecular vibrations with a microcavity mode_ , New J. Phys. 17, 053040 (2015).

##  I Supplementary Material for: ”Organic polariton lasing and the weak- to
strong-coupling crossover”

##  II Role of diamagnetic $A^{2}$ term

This section discusses the role of the diamagnetic $A^{2}$ term, allowing one
to disentangle the effects of this term vs other effects of strong coupling on
the form of the phase diagram.

As noted in the text, this $A^{2}$ term is included in order to avoid a
ground-state superradiant transition, following the results of Rzażewski  _et
al._ (1975). Without the diamagnetic term, the ground state of the Dicke model
at strong coupling is a state with a macroscopic occupation of the photon Hepp
and Lieb (1973), even in the absence of pumping. The inclusion of the $A^{2}$
term, with a coefficient determined by the oscillator strength sum rule
prevents this transition occurring Rzażewski  _et al._ (1975), so the ground
state remains stable at all coupling strengths. Recently the question of
whether a phase transition is in fact possible in the ground state has been
re-opened, with a consensus developing that a transition is possible, but
driven by dipole-dipole interactions between atoms when in the appropriate
geometry Grießer  _et al._ (2016); De Bernardis  _et al._ (2018). We neglect
the possibility of this ferromagnetic transition here, and include an $A^{2}$
term so that only the transition induced by incoherent pumping is present.

Figure S1: Role of $A^{2}$ term, fixed coupling, parameters as in Fig. 1.
White region denotes normal state, black denotes lasing, and red the
superradiant phase. Note the different $x$-axis scale in panel (f).

Figure S1 compares the phase diagrams (vs bare cavity frequency and pumping)
with (left) and without (right) the diamagnetic $A^{2}$ term. In the phase
diagrams without the $A^{2}$-term a superradiant phase is present at low
photon frequencies and low pumping. As discussed elsewhere Kirton and Keeling
(2018), this superradiant state can be distinguished from the lasing state by
the nature of the instability and the spectrum of the emission: ground state
superradiance is a stationary steady state, while the lasing state is at
finite frequency, and thus is time dependent in the lab frame. As a
consequence, the superradiant state instability involves a single unstable
frequency (regions marked red on Fig. S1), while the lasing state instability
involves a complex conjugate pair of unstable modes (regions marked black on
Fig. S1).

Regarding the effect of the $A^{2}$ term on the shape of the polariton lasing
region, we see from Fig. S1 that the lasing phase remains a similar shape,
although the characteristic frequency ranges are significantly rescaled at
strong coupling. As pointed out in the Letter, the $A^{2}$ term has a twofold
effect. Firstly, it leads to renormalization of the effective photon frequency
$\omega_{c}^{\text{eff}}$, which leads to a $g$-dependent effective detuning.
Secondly, it reduces the effective exciton-photon coupling. In the ultrastrong
coupling shown in Fig. S1(e,f), the combination of these effects results in
the effective decoupling of the photon from the molecular excitations, causing
the lasing region to shrink to a significantly lower frequency range. As the
shape is similar, we can be clear the bending of the lasing lobe is not
associated with the $A^{2}$ term.

##  III Effect of vibrational truncation

In Fig. S2 we show the dependence of the phase diagram on the number of
vibrational states $N$. We show results for $S=0.1$ and for two values of
matter-light coupling, $g\sqrt{\mathcal{N}_{m}}=0.1$ and
$g\sqrt{\mathcal{N}_{m}}=1.0$, corresponding to the top left and bottom right
panels of Fig. 2 in which we used $N_{v}=4$. For the weaker coupling,
$g\sqrt{\mathcal{N}_{m}}=0.1$, as $N_{v}$ increases the phase diagram retains
a similar structure, but one can see the appearance of new peaks at
$\omega_{c}\approx 0.6,0.4,0.2$ corresponding to (2-0), (3-0) and (4-0)
transitions respectively. For stronger coupling,
$g\sqrt{\mathcal{N}_{m}}=1.0$, we see a similar emergence of extra peaks, and
a slight evolution of the extended lasing lobe. In both cases we see that the
results do not differ significantly between $N_{v}=4$ and $N_{v}=5$,
confirming we can safely use $N_{v}=4$ in our Letter.

Figure S2: Dependence of phase diagrams on the number of vibrational levels
retained, $N_{v}$. Other parameters as in Fig. 1.

For $S=0.1$, we can in fact expect from simple arguments that the dominant
physics comes from single vibrational excitation. For an isolated molecule,
neglecting strong coupling, the probability for a ($n$–0) transition is
$P_{n}=e^{-S}S^{n}/n!$, thus for $S=0.1$, we would expect results to be
dominated by $n=0,1$. The behavior shown in Fig. S2 shows that despite this,
higher vibrational states do play some role even for these small values of
$S$.

##  IV Nature of the lasing instability

In this section we provide further details on the calculation of the
vibrational composition of the lasing mode, and present figures illustrating
the evolution of the spectrum.

###  IV.1 Composition of linear stability eigenmodes

To determine the molecular transitions corresponding to a given unstable mode,
we note that from the eigenvector
$\upsilon=\left(\delta\alpha,\delta\alpha^{*},\mathbf{\delta\ell}\right)^{\intercal}$
of the stability matrix $\mathcal{M}$, we can extract the matter part
$\delta\ell_{i}$ and so construct the corresponding molecular density matrix,
$\rho(t)=\rho_{\text{ns}}+\mathcal{A}(\delta\ell_{i}e^{\xi
t}+\delta\ell_{i}^{\ast}e^{\xi^{\ast}t})\lambda_{i}/2$ where $\mathcal{A}$ is
an arbitrary amplitude and $\rho_{\text{ns}}$ is the normal state density
matrix. (For simplicity of notation, we neglect the superscripts on
$\delta\ell$ and $\xi$ labeling eigenmodes.) The complex conjugates appearing
here are required in order to guarantee Hermiticity (given the Gell Mann
matrices are defined to be Hermitian). This step is crucial since the
equations of motion mix $\ell_{i}$ and $\ell_{i}^{\ast}$. To find the
amplitude of the oscillatory component of a given element of the density
matrix $\rho_{ij}$, we first define the matrix $r=\delta\ell_{i}\lambda_{i}/2$
and then note the oscillatory component of $\rho_{ij}$ takes the form
$\delta\rho_{ij}=r_{ij}e^{\xi t}+r_{ji}^{\ast}e^{\xi^{\ast}t}$.

While the diagonal components of a density matrix give the molecular state
populations, nondiagonal ones correspond to coherences. In particular, the
block of the density matrix which is off-diagonal in terms of electronic
states gives the weights of molecular transitions involved in lasing. For the
diagonal components of the density matrix we can immediately write
$\delta\rho_{ii}=2|r_{ii}|\cos\big{(}\xi^{\prime\prime}t+\text{Arg}(r_{ii})\big{)}e^{\xi^{\prime}t}$
where $\xi^{\prime},\xi^{\prime\prime}$ denote real and imaginary parts of the
eigenvalue respectively. For the off diagonal components, the behavior is more
complex as in general $|r_{ij}|\neq|r_{ji}|$; we thus expect that the quantity
$\delta\rho_{ij}$ traces an elliptical spiral in the Argand plane,
$\delta\rho_{ij}=[r_{A}\cos(\xi^{\prime\prime}t+\phi)+ir_{B}\sin(\xi^{\prime\prime}t+\phi)]e^{i\theta+\xi^{\prime}t}$.
One may readily show that the semi-major and semi-minor axes of this ellipse
are given by $r_{A,B}=|r_{ij}|\pm|r_{ji}|$, while
$\phi,\theta=[\text{Arg}(r_{ij})\pm\text{Arg}(r_{ji})]/2$. Given this
behavior, we define the amplitude of a given component by the semi-major axis.
The contribution of a given molecular transition $(n-m)$ to the lasing mode
thus corresponds to the amplitude
$|r_{n\downarrow,m\uparrow}|+|r_{m\uparrow,n\downarrow}|$.

###  IV.2 Evolution of spectrum with coupling

In Fig. 3 we showed the evolution of the spectrum with increasing pumping.
Here, in Fig. S3, we complement this by showing how the spectrum evolves with
increasing matter-light coupling at weak pumping. We clearly see how the
polariton splitting increases and self-tuning develops with increasing
coupling. We may also note that the real part of the unstable eigenvalue
increases, so the lasing instability will develop on a shorter timescale.

Figure S3: Evolution of spectrum as coupling $g\sqrt{\mathcal{N}_{m}}$
increases, for weak pumping $\Gamma_{\uparrow}=0.4\Gamma_{\downarrow}$. Line
colors and styles as in Fig. 3, parameters as in Fig. 1.

###  IV.3 Transition from polariton to photon lasing

In Fig. S4 we show the spectrum for the same parameters as in Fig. 3(i)
(strong coupling and strong pumping), but over an extended axis range.This
figure illustrates the transition between the two different lasing states
which occurs in this regime. For low cavity frequencies (around
$\omega_{c}\lesssim 1.5$) the system supports self-tuning and we get polariton
lasing. However at larger cavity frequencies, the larger detuning in
combination with strong pumping suppresses strong coupling leading to “normal”
lasing at effective photon frequency $\omega_{c}^{\text{eff}}$. As this regime
leads to effectively weak coupling photon lasing, we find the growth rate for
this second lasing instability is much slower than in the polariton lasing
regime.

We note that in this figure, because we use $N_{v}=4$, we see only a limited
range of vibrational sidebands. A calculation with large $N_{v}$ would show a
sequence of sidebands up to higher frequency (but with a weight that decreases
as one goes to higher frequencies).

Figure S4: Spectrum for $g\sqrt{\mathcal{N}_{m}}=1.0$,
$\Gamma_{\uparrow}=1.8\Gamma_{\downarrow}$ on an enlarged axis. Line colors
and styles as in Fig. 3, parameters as in Fig. 1.

##  V Phase boundary vs coupling strength

In this section we discuss further the phase boundaries and minimum lasing
threshold shown in Fig. 4. We provide further cuts of the phase diagram at
different cavity frequencies, present the optimal cavity frequency at a given
coupling strength, and provide details of the weak coupling calculation that
is shown in Fig. 4(a-c).

###  V.1 Evolution of phase boundary with cavity frequency

We first provide further detail on the how the phase boundary evolves with
cavity frequency. Figure 4(a-c) presents the phase diagrams at three specific
frequencies; to show more clearly how the phase boundary evolves with
frequency, the top panel of Fig. S5 shows a sequence of curves, giving the
evolution of the phase boundary with bare photon frequency. Note that the
phase boundary has different topology for $\omega_{c}\leq 0.9$ and
$\omega_{c}\geq 1$. For low frequency, there is a single lasing region (but
with a hole appearing at large coupling around
$\Gamma_{\uparrow}\simeq\Gamma_{\downarrow}$). Around $\omega_{c}=1$, this
hole reaches the phase boundary, dividing the lasing region into two: a large
region at high pumping, and a smaller region, below inversion and at strong
coupling.

Figure S5: Top: Critical pump strength as a function of coupling strength for
various bare cavity frequencies $\omega_{c}$. Bottom: The cavity frequencies
$\omega_{c}$ at which the threshold has a local minimum as a function of
coupling $g\sqrt{\mathcal{N}_{m}}$. The thick gray line highlights which
frequency corresponds to the global minimum. Parameters as in Fig. 1.

###  V.2 Optimal cavity frequency

The lower panel of Fig. S5 shows the cavity frequencies at which the threshold
has local minima. As could already be seen in Fig. 2, the threshold shows
multiple local minima at vibronic replicas
$\omega_{c}=\varepsilon-n\omega_{v}$. At the lowest coupling, the minimal
threshold is near $\omega_{c}=\varepsilon$ because this corresponds to the
peak emission — for very low coupling, the critical condition is for emission
to beat cavity loss. Note that the threshold lasing power diverges around
$g\sqrt{\mathcal{N}_{m}}\simeq 0.0024$, i.e. for the parameters we use, there
is no lasing at weaker coupling than this. When coupling increases to
$g\sqrt{\mathcal{N}_{m}}\approx 0.0035$, we get another local minimum emerging
at $\omega_{c}\approx\varepsilon-\omega_{v}$ corresponding to the (1-0)
transition. For $g\sqrt{\mathcal{N}_{m}}\gtrsim 0.004$ this new minimum
becomes the global minimum, i.e. it leads to the lowest threshold pumping as
shown in Fig. 4(d). As we enter strong coupling regime $g\gtrsim 0.03$ we get
one more local minimum driven by (2-0) transition and then another one due to
(3-0) transition. All these local minima of pumping can be seen in Fig. 2.

Increasing the coupling further (into the strong coupling self-tuning regime)
we see a continuous evolution of (2-0) and (3-0) branches to lower frequency,
while the (1-0) branch moves to higher frequencies. The evolution of the (1-0)
transition frequency can be associated with the self-tuning effect which
occurs at strong coupling, i.e. self-tuning allows the lasing threshold for
this transition to move out to high cavity frequencies. In contrast, the (2-0)
and (3-0) transitions do not show any strong-coupling self-tuning effect — the
lasing transitions at these frequencies appear similar to weak coupling
lasing. The reason these lasing frequencies shift to lower frequencies is the
renormalization of effective photon frequency by the $A^{2}$ term. i.e., the
resonance condition is in fact
$\omega_{c}^{\text{eff}}\simeq\varepsilon-n\omega_{v}$, and the dependence of
$\omega_{c}^{\text{eff}}$ on $g\sqrt{\mathcal{N}_{m}}$ means the resonance
moves to lower bare cavity frequencies as coupling increases.

In the ultrastrong coupling regime, the (1-0) and (2-0) branches give almost
the same threshold pumping around $\Gamma_{\uparrow}^{th}\approx
0.22\Gamma_{\uparrow}$, while the (3-0) transition has
$\Gamma_{\uparrow}^{th}\approx 0.38\Gamma_{\uparrow}$, as can also be seen in
Fig. 2.

###  V.3 Comparison to weak coupling theory

In Fig. 4, we compared the strong coupling phase diagram to a weak coupling
prediction of the phase boundary. Here we provide details of this weak
coupling prediction. The weak coupling theory of the Dicke-Holstein model has
been discussed in Ref. Kirton and Keeling (2013, 2015). In those works, it is
shown that treating matter-light coupling perturbatively, one can derive a
weak coupling master equation:

$\dot{\rho}(t)=-i\left[\tilde{H}_{0},\rho\right]+\kappa\mathcal{L}[a]+\sum_{n}\Big{(}\Gamma_{\downarrow}\mathcal{L}[\sigma_{n}^{-}]+\Gamma_{\uparrow}\mathcal{L}[\sigma_{n}^{+}]\\\
+\Gamma_{E}\mathcal{L}[a^{\dagger}\sigma^{-}]+\Gamma_{A}\mathcal{L}[a\sigma^{+}]\Big{)}\rho,$
(S1)

where the free Hamiltonian is
$\tilde{H}_{0}=(\omega_{c}-\varepsilon)a^{\dagger}a$ and the rates of emission
and absorption processes are given by

$\Gamma_{A,E}=g^{2}\int_{-\infty}^{\infty}dte^{\pm
i\omega_{c}t}\langle\sigma^{-}(t)\sigma^{+}(0)\rangle_{0}.$ (S2)

The two-time correlation function $\langle\ldots\rangle_{0}$ is calculated for
free molecules, neglecting the matter-light coupling. In Refs. Kirton and
Keeling (2013, 2015), a different form of dissipation was assumed for the
vibrational modes, allowing an explicit form for $\Gamma_{A,E}$ through a
Keldysh path integral. The form of dissipation considered here and written in
Eq. (2) does not allow such a form, however the integrals in $\Gamma_{A,E}$
can be numerically calculated from the master equation using the quantum
regression theorem Breuer and Petruccione (2002).

From this equation, we can use the same mean-field decoupling discussed
earlier, i.e. $\langle a\sigma^{\pm}\rangle=\langle
a\rangle\langle\sigma^{\pm}\rangle$ and $\langle a^{\dagger}a\rangle=|\langle
a\rangle|^{2}$ to yield a set of non-linear coupled equations:

$\displaystyle\partial_{t}\langle a\rangle=$
$\displaystyle\frac{\mathcal{N}_{m}}{4}\Big{[}\Gamma_{A}\left(1-\langle\sigma^{z}\rangle\right)-\Gamma_{E}\left(1+\langle\sigma^{z}\rangle\right)\Big{]}\langle
a\rangle$
$\displaystyle-\left[i(\omega_{c}-\varepsilon)+\frac{\kappa}{2}\right]\langle
a\rangle,$ (S3) $\displaystyle\partial_{t}\langle\sigma^{z}\rangle=$
$\displaystyle\Big{(}1-\langle\sigma^{z}\rangle\Big{)}\Big{(}\Gamma_{\uparrow}+\Gamma_{A}|\langle
a\rangle|^{2}\Big{)}$
$\displaystyle-\Big{(}1+\langle\sigma^{z}\rangle\Big{)}\Big{(}\Gamma_{\downarrow}+\Gamma_{E}|\langle
a\rangle|^{2}\Big{)}.$ (S4)

We can then again consider linear stability of the normal state, by
considering fluctuations $\langle a\rangle=\alpha$ and
$\langle\sigma^{z}\rangle=\langle\sigma^{z}\rangle_{ns}+z$ as in the Letter.
One may see that $\alpha$ and $z$ are uncoupled, so the corresponding
equations can be solved trivially, and the linearized eigenvalue is given by

$\xi=-i(\omega_{c}-\epsilon)-\frac{\kappa}{2}+\frac{\mathcal{N}_{m}}{4}(\Gamma_{E}-\Gamma_{A})+\frac{\mathcal{N}_{m}}{4}(\Gamma_{E}+\Gamma_{A})\langle\sigma^{z}\rangle_{\text{ns}},$
(S5)

where
$\langle\sigma^{z}\rangle_{\text{ns}}=(\Gamma_{\uparrow}-\Gamma_{\downarrow})/(\Gamma_{\uparrow}+\Gamma_{\downarrow})$.
One can extract threshold pumping as a function of coupling $g$, by solving
$\text{Re}[\xi]=0$, which gives:

$\Gamma_{\uparrow}^{th}=\frac{\mathcal{N}_{m}\Gamma_{A}+\kappa}{\mathcal{N}_{m}\Gamma_{E}-\kappa}\Gamma_{\downarrow}.$
(S6)