# Mode switching dynamics in organic polariton lasing

Antti J. Moilanen<EMAIL_ADDRESS>Department of Applied Physics, Aalto
University School of Science, FI-00076 Aalto, Finland  Photonics Laboratory,
ETH Zürich, CH-8093 Zürich, Switzerland  Kristín B. Arnardóttir  DTU Fotonik,
Technical University of Denmark, 2800 Kongens Lyngby, Denmark  NanoPhoton-
Center for Nanophotonics, Technical University of Denmark, 2800 Kongens
Lyngby, Denmark  SUPA, School of Physics and Astronomy, University of St
Andrews, KY16 9SS St Andrews, United Kingdom  Jonathan Keeling  SUPA, School
of Physics and Astronomy, University of St Andrews, KY16 9SS St Andrews,
United Kingdom  Päivi Törmä<EMAIL_ADDRESS>Department of Applied
Physics, Aalto University School of Science, FI-00076 Aalto, Finland

(September 3, 2024)

###### Abstract

We study the dynamics of multimode polariton lasing in organic microcavities
by using a second-order cumulant equation approach. By inspecting the time
evolution of the photon mode occupations, we show that if multiple lasing
peaks are observed in time-integrated mode occupations, the reason can be
either bi-modal lasing or temporal switching between several modes. The former
takes place within a narrow range of parameters while the latter occurs more
widely. We find that the origin of the temporal switching is different in the
weak- and strong-coupling regimes. At weak coupling slope efficiency is the
determining factor, while for strong coupling it is changes in the eigenmodes
and gain spectrum upon pumping. This difference is revealed by investigating
the photoluminescence and momentum-resolved gain spectra. Our results
underscore the importance of understanding the time evolution of the
populations when characterizing the lasing behavior of a multimode polariton
system, and show how these features differ between weak and strong coupling.

##  I Introduction

In a system where lasing occurs, switching between lasing modes can occur, as
a function of parameters (e.g. excitation power) or as a function of time.
Mode switching is known to occur in semiconductor lasers that support several
densely spaced longitudinal modes. Multiple lasing peaks in the time-
integrated luminescence spectrum, as well as temporal switching between the
modes, have been observed in many different contexts Furfaro  _et al._ (2004);
Pedaci  _et al._ (2005); Shore and Kane (2005); Gordon  _et al._ (2008);
Wright  _et al._ (2017). The origin and nature of such mode switching differs
between systems. For example, the mode switching in edge-emitting
semiconductor lasers and vertical cavity surface emitting lasers (VCSELs) has
been attributed to spatial hole burning Herre and Barabas (1989); Lenstra and
Yousefi (2014); Valle  _et al._ (1995); Law and Agrawal (1997) and to four-
wave mixing Yacomotti  _et al._ (2004). Understanding the nature of mode
switching in a given system, and how to control it, is an important step
toward enabling applications.

Exciton-polaritons—polaritons for brevity—are the hybrid of a photon and an
exciton under strong light-matter coupling. Strong coupling can be engineered
by creating optical microcavities filled with an active layer of emitter
material. The active layer can take several forms, including inorganic
semiconductors, two-dimensional materials such as transition metal
dichalcogenides, or—the focus of this paper—organic molecules. Polariton
lasing (or polariton condensation) has been studied extensively during the
past two decades Deng  _et al._ (2010); Keeling and Berloff (2011); Keeling
and Kéna-Cohen (2020).

While the majority of the experimental and theoretical works have focused on
single-mode polariton lasing, multimode operation has been reported as well.
Switching of macroscopic polariton population between different spatial modes
or momentum states has been observed as a function of pump power Sun  _et al._
(2018); Kusudo  _et al._ (2013); Kim  _et al._ (2016); Maragkou  _et al._
(2010); Sawicki  _et al._ (2021). The switching between higher- and lower-
energy polariton modes has been attributed to polariton-polariton scattering,
which is stimulated by Bose enhancement at higher pump powers, driving the
relaxation of polaritons to the lowest-energy mode Doan  _et al._ (2005,
2006); Maragkou  _et al._ (2010); Sawicki  _et al._ (2021). Typically mode
switching has been observed only for specific detuning conditions. In coupled
wells or lattices of coupled polariton condensates, competition between
different modes has been studied experimentally Lai  _et al._ (2007); Kim  _et
al._ (2011); Jacqmin  _et al._ (2014); Amo and Bloch (2016); Scafirimuto  _et
al._ (2021); Boulier  _et al._ (2020) and numerically by using the generalized
Gross-Pitaevskii equations Eastham (2008); Wouters (2008). Recently, it was
noted that simultaneous polariton lasing can occur in a vertical lasing mode
and a horizontal guided mode in planar microcavities, and that such phenomena
could have taken place in several previous experiments where polariton lasing
in the vertical lasing mode has been reported with ZnO-based microcavities
Jamadi  _et al._ (2019). Mode competition dynamics has also recently been
explored in InGaAs based cavities Töpfer  _et al._ (2020). In closely related
experiments on photon Bose-Einstein condensates (BECs)—with weak light-matter
coupling in dye-filled microcavities—multimode condensation has been observed
experimentally and studied by second-order cumulant equations Vlaho and
Eckardt (2021); Vlaho  _et al._ (2019); Marelic  _et al._ (2016); Kurtscheid
_et al._ (2019). Similarly, for strongly-coupled plasmonic BECs, multiple
peaks near the ground state have been observed Väkeväinen  _et al._ (2020).

In our previous work Arnardottir  _et al._ (2020), we introduced a new model
for multimode polariton lasing in organic microcavities and computed steady-
state lasing phase diagrams as a function of pump strength and detuning of the
photon modes from the exciton. We showed that at certain pump strengths
multiple lasing peaks may be observed in the steady-state photon mode
occupations. Our second-order cumulant model goes beyond mean-field
descriptions, capturing the effects of non-lasing modes and fluctuations on
the system’s behavior. Retaining such fluctuations is important when
considering switching: the (initially small) populations of non-lasing modes
act as the seed for these modes to become macroscopically populated at later
times. In contrast, for mean-field approaches, explicit seed or noise terms
would be required.

In this paper, we investigate the dynamics of mode switching in multimode
polariton lasing, and analyze the gain spectrum in more depth than in
Arnardottir  _et al._ (2020), including momentum-resolved gain spectra. While
the model we consider is specific to organic molecules, the concepts we
discuss have relevance more widely. We find that in the case where several
lasing peaks are observed in the time-integrated mode occupations, the modes
are not necessarily lasing simultaneously, but the lasing mode may switch
(deterministically) as a function of time. Furthermore, even when only a
single lasing peak is observed, switching between different modes can occur
before the system reaches a steady state. We study mode switching in both weak
and strong coupling regimes.

Mode switching occurs in general because the mode which lases first does not
fully cross-saturate the gain for other modes. This raises a question of what
determines which other modes become competitive, and ultimately take over. At
weak coupling, the choice of which mode is selected can be explained by a
standard argument Siegman (1986) considering the different thresholds and
slope efficiencies (gradient of photon population vs pump) of each mode: When
several modes are above the threshold, the mode with highest gain will
eventually suppress the others. At strong coupling, the situation becomes more
complicated because of energy shifts. The polariton energies show a blue shift
as a function of pump strength, due to effective interactions that arise from
saturation of the gain medium. Importantly, we also show that the peak gain
shifts towards higher energies as the pump strength increases. This shift can
be attributed to dynamics of the vibrational degrees of freedom of the
molecules. The mode switching at strong coupling occurs as a result of
competition between both shifts in the gain and the polariton energies. Our
model considers spatially homogeneous systems, thus our results suggest that
mode switching does not necessarily require redistribution of real space
intensity.

The rest of the paper is organized as follows. In Section II, we present the
system and the model. In Section III, we study the temporal characteristics of
multimode lasing. In Section IV, we present the gain and photoluminescence
spectra as a function of pump strength. In Section V we consider pulsed
excitation, which is relevant for typical experiments on organic exciton-
polariton systems. Finally, in Section VI we discuss the interpretation and
implications of the results. Appendices provide technical details on the
cumulant equations and results for time-integrated populations.

##  II System and model

###  II.1 Density matrix equation of motion

We consider a planar microcavity that hosts multiple photon modes and an
active layer of emitters. We model the emitters as $N_{\mathrm{m}}$
vibrationally dressed two-level systems placed homogeneously in the cavity
plane. The system is thus described by the Tavis–Cummings–Holstein Hamiltonian
in the rotating wave approximation (RWA) Ćwik  _et al._ (2014); Strashko  _et
al._ (2018); Arnardottir  _et al._ (2020):

$\displaystyle H=$
$\displaystyle\sum_{\mathbf{k}}\omega_{\mathbf{k}}a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}+\sum_{n,{\mathbf{k}}}\left(g_{n,{\mathbf{k}}}a_{\mathbf{k}}\sigma_{n}^{+}+g_{n,{\mathbf{k}}}^{\ast}a_{\mathbf{k}}^{\dagger}\sigma_{n}^{-}\right)$
$\displaystyle+\sum_{n}\left[\frac{\varepsilon}{2}\sigma^{z}_{n}+\omega_{v}\left(b_{n}^{\dagger}b_{n}+\sqrt{S}\left(b_{n}^{\dagger}+b_{n}\right)\sigma_{n}^{z}\right)\right].$
(1)

Here $\varepsilon$ is the transition energy of the molecule and Pauli matrices
$\sigma_{n}^{z,+,-}$ correspond to the electronic states of molecule $n$. The
operator $b_{n}^{\dagger}$ creates a vibrational excitation of energy
$\omega_{v}$ on molecule $n$ while $a_{\mathbf{k}}^{\dagger}$ creates a photon
with energy $\omega_{\mathbf{k}}$ in the cavity mode with wave vector
${\mathbf{k}}$.

The photon energies are given by
$\omega_{\mathbf{k}}=\omega_{0}+{E_{\rho}}{(K_{x}^{2}+K_{y}^{2})}/{N_{\textrm{m}}}$,
where $K_{x,y}$ are integers labelling a discrete set of momentum vectors in
the two-dimensional plane of the cavity. The energy scale
$E_{\rho}=\pi^{2}\rho_{\text{2D}}{\hbar^{2}}/{2m}$ is defined in terms of the
molecular density $\rho_{2D}$ and effective photon mass $m$. The total number
of photon modes is truncated to $N_{\mathrm{ph}}$. The light-matter coupling
is given by the coupling constant
$g_{n,{\mathbf{k}}}=e^{-i\mathbf{k}\cdot\mathbf{r}_{n}}\Omega_{R}/\sqrt{N_{\mathrm{m}}}$,
where $\mathbf{r}_{n}$ is the position of the $n$th molecule. These plane-wave
cavity modes correspond to considering periodic boundary conditions. Coupling
of the electronic levels to vibrational modes results in broadening of
emission and absorption spectra and in the occurrence of a Stokes shift, which
is parametrized by the Huang-Rhys factor $S$ Jong  _et al._ (2015). In later
discussions it will be useful to use the notation $(n-m)$ to denote
transitions between the $n$th vibrational ground state and the $m$th
vibrational excited state; see Fig. 1(c,d). The total number of vibrational
states on each molecule will be truncated to $N_{\mathrm{v}}$.

Figure 1: Schematics of the system and model. (a) A planar microcavity with
multiple photon modes and a layer of emitter material. The cavity modes have
in-plane momentum $k$ and they decay by rate $\kappa$. The photons in the
cavity modes can strongly couple to the excitons, hybridizing to form exciton-
polaritons. (b) Discretized dispersion relation of the uncoupled cavity modes.
(c) Two-level system with both ground and excited state dressed by
$N_{\mathrm{v}}$ vibrational states. The bare transition energy between the
two electronic states is denoted by $\varepsilon$, the vibrational frequency
by $\omega_{v}$, and the Huang-Rhys factor by $S$. The rates of the drive and
dissipation processes involved are the external pumping $\Gamma_{\uparrow}$,
spontaneous decay $\Gamma_{\downarrow}$, dephasing $\Gamma_{z}$, and thermal
excitation and dissipation of the vibrational modes
$\gamma_{\uparrow,\downarrow}$. (d) Absorption and emission spectra of the
uncoupled emitters. The yellow shaded area corresponds to the energy range
which the uncoupled photon modes cover, starting from $\omega_{0}$ (0.65 eV).
The zero-phonon-line (0-0) and the emission shoulders related to the two first
vibrational transitions (1-0), (2-0) are labeled.

The incoherent processes that occur can be described by the Lindblad
formalism. The equation of motion for the density matrix of the system reads

$\partial_{t}\rho=-i[H,\rho]+\sum_{\mathbf{k}}\kappa\mathcal{L}[a_{\mathbf{k}}]+\sum_{n}\big{(}\Gamma_{\uparrow}\mathcal{L}[\sigma^{+}_{n}]+\Gamma_{\downarrow}\mathcal{L}[\sigma^{-}_{n}]\\\
+\Gamma_{z}\mathcal{L}[\sigma^{z}_{n}]+\gamma_{\uparrow}\mathcal{L}[b^{\dagger}_{n}+\sqrt{S}\sigma^{z}_{n}]+\gamma_{\downarrow}\mathcal{L}[b_{n}+\sqrt{S}\sigma^{z}_{n}]\big{)},$
(2)

with $\mathcal{L}[X]=X\rho X^{\dagger}-\frac{1}{2}(X^{\dagger}X\rho+\rho
X^{\dagger}X)$. The equations account for cavity loss at rate $\kappa$,
dephasing of the excitons $\Gamma_{z}$, and incoherent pumping and decay of
excitons with rates $\Gamma_{\uparrow}$ and $\Gamma_{\downarrow}$,
respectively 111We note that compared to our previous work Ref. Arnardottir
_et al._ (2020) we have corrected the sign of the $\sqrt{S}\sigma^{z}_{n}$
term in Eq. (2). This does not significantly change the results.. The last two
terms describe thermal relaxation and excitation with rates
$\gamma_{\downarrow}=\gamma_{v}(n_{b}+1)$ and
$\gamma_{\uparrow}=\gamma_{v}n_{b}$ where
$n_{b}=\left[\text{exp}(\omega_{v}/k_{B}T_{v})-1\right]^{-1}$ is the Bose-
Einstein distribution at temperature $T_{v}$. The values of all parameters
used for the simulations below are given in Table 1.

Table 1: Simulation parameters. Parameter | Value | Notes
---|---|---
$N_{\mathrm{m}}$ | $10^{8}$ |
$N_{\text{v}}$ | 4 |
$N_{\mathrm{ph}}$ | 9 |
$\varepsilon$ | 1 eV |
$\omega_{0}$ | 0.65 eV | $-2\omega_{v}<\omega_{0}-\epsilon<-\omega_{v}$
$E_{\rho}$ |  $5\times 10^{5}$ eV |
$\Omega_{\text{R}}$ | 0.1 eV & 0.4 eV |
$\omega_{v}$ | 0.2 eV |
S | 0.10 |
$\kappa$ |  $10^{-4}$ eV |  $\sim$ (40ps)-1
$\Gamma_{\downarrow}$ |  $10^{-4}$ eV |  $\sim$ (40ps)-1
$\Gamma_{z}$ | 0.03 eV |  $\sim$ (1ps)-1
$\gamma_{v}$ | 0.02 eV |
$k_{B}T_{v}$ | 0.025 eV | Room temperature

###  II.2 Second-order cumulant equations

The size of Hilbert space for $N_{\mathrm{m}}$ molecules, each with
$2N_{\mathrm{v}}$ internal states, is too large to directly simulate except
for small $N_{\mathrm{m}}$. To make the simulation computationally tractable
for realistic system sizes, we proceed from Eq. (2) by using the second-order
cumulant expansion. This provides a (closed) set of equations of motion for
first- and second-order correlations between molecule and photon operators
Kirton and Keeling (2017); Sánchez-Barquilla  _et al._ (2020); Arnardottir
_et al._ (2020).

In order to make the cumulant expansion, it is helpful to first relabel the
molecular degrees of freedom. The two electronic states $(\sigma)$ and the
$N_{\mathrm{v}}$ vibrational states ($b$) are combined into
$2N_{\mathrm{v}}$-level operators, which can be represented by generalized
Gell-Mann matrices, $\lambda$ (see Appendix A for details). The Hamiltonian
then becomes:

$H=\sum_{\mathbf{k}}\omega_{\mathbf{k}}a_{\mathbf{k}}^{\dagger}a_{\mathbf{k}}+\sum_{n}\left[A_{i}+\sum_{\mathbf{k}}(B_{i}a^{\dagger}_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{r}_{n}}+\textrm{h.c.})\right]\lambda_{i}^{(n)},$
(3)

where the vectors $A_{i}$ and $B_{i}$ include the terms in Eq. (1). Tensor
sums that run over Gell–Mann operator indices $i$ are implicit. In the new
basis, the master equation in Eq. (2) takes the form

$\partial_{t}\rho=-i[H,\rho]+\sum_{\mathbf{k}}\kappa\mathcal{L}[a_{\mathbf{k}}]+\sum_{\mu,n}\mathcal{L}\left[\gamma_{i}^{\mu}\lambda^{(n)}_{i}\right],$
(4)

where $\mu$ labels the drive and dissipation processes related to molecules:
(1) pump, (2) decay, (3) dephasing, and vibrational (4) excitation and (5)
decay.

We next write down equations of motion for first- and second-order
correlations of the operators
$\lambda^{(n)}_{i},a_{\mathbf{k}},a^{\dagger}_{\mathbf{k}}$. To produce a
closed set of equations of motion, all third- or higher-order correlations are
split into products of first and second-order correlations, by setting third-
order cumulants (e.g. $\langle ABC\rangle_{c}$ ) to zero. This allows to
decompose the third-order terms as $\langle ABC\rangle=\langle A\rangle\langle
BC\rangle+\langle B\rangle\langle AC\rangle+\langle C\rangle\langle AB\rangle$
Gardiner (2009). To further simplify the derivation, we make use of the
conservation of number of excitations implied by the RWA and divide the
operators $\lambda_{i}$ into three groups that either increase ($+$), decrease
($-$), or conserve ($z$) the electronic excitations. With this classification,
expectation values of operators that do not conserve the number of excitations
can be taken zero, such as the first-order terms $\langle
a_{\mathbf{k}}\rangle,\langle
a_{\mathbf{k}}^{\dagger}\rangle,\langle\lambda_{i_{-}}\rangle,\langle\lambda_{i_{+}}\rangle$.

Assuming a spatially homogeneous distribution (or random distribution with a
large number) of emitters in the microcavity, we may approximate summation
over sites by using
$\sum_{n}e^{i(\mathbf{k}-\mathbf{k}^{\prime})\cdot\mathbf{r}_{n}}=N_{\mathrm{m}}\delta_{\mathbf{k},\mathbf{k}^{\prime}}$,
allowing the use of momentum conservation. For the number-conserving Gell-Mann
operators we then obtain a site-independent form
$\ell_{i}=\langle\lambda_{i_{z}}^{(n)}\rangle$. For the two other groups of
Gell-Mann operators, $\lambda_{i_{+}}^{(n)}$ and $\lambda_{i_{-}}^{(n)}$, the
first order expectations are zero as they do not conserve the number of
excitations. We can write the Fourier components of the second-order
correlations as $d_{ij}^{\mathbf{k}}=\sum_{n,m\neq
n}e^{i\mathbf{k}\cdot(\mathbf{r}_{n}-\mathbf{r}_{m})}\langle\lambda_{i_{+}}^{(n)}\lambda_{i_{-}}^{(m)}\rangle/N_{\mathrm{m}}^{2}$
and $c_{i}^{\mathbf{k}}=\sum_{n}e^{i\mathbf{k}\cdot\mathbf{r}_{n}}\langle
a_{\mathbf{k}}\lambda_{i_{+}}^{(n)}\rangle/N_{\mathrm{m}}$. Photon mode
occupations are denoted by $n^{\mathbf{k}}=\langle
a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}\rangle$. The equations of motion for
the non-vanishing correlations are thus:

$\displaystyle\partial_{t}n^{\mathbf{k}}=-\kappa
n^{\mathbf{k}}-2N_{\mathrm{m}}\mathrm{Im}{B_{i}c_{i}^{\mathbf{k}}},$
$\displaystyle\partial_{t}\ell_{i}=\xi_{ij}\ell_{j}+\phi_{i}-8\mathrm{Re}{\beta_{ij}\sum_{k}c_{j}^{\mathbf{k}}},$
$\displaystyle\partial_{t}c_{i}^{\mathbf{k}}=\left[X_{ij}-\left(i\omega_{\mathbf{k}}+\frac{\kappa}{2}\right)\delta_{ij}\right]c_{j}^{\mathbf{k}}+2\beta^{\ast}_{ij}\ell_{j}n_{\mathbf{k}}-iB_{j}d_{ij}^{\mathbf{k}}$
$\displaystyle-\frac{i}{N_{\mathrm{m}}}\left(\zeta_{ij}\ell_{j}+\frac{B_{i}}{2N_{\mathrm{v}}}\right),$
$\displaystyle\partial_{t}d_{ij}^{\mathbf{k}}=X^{\ast}_{ip}d^{\mathbf{k}}_{pj}+X_{jp}d^{\mathbf{k}}_{ip}+2\ell_{p}\left(\beta_{ip}\tilde{c_{j}}^{\mathbf{k\ast}}+\beta^{\ast}_{jp}\tilde{c_{i}}^{\mathbf{k}}\right),$
(5)

where
$\tilde{c_{i}}^{\mathbf{k}}=c_{i}^{\mathbf{k}}-\sum_{\mathbf{q}}c_{i}^{\mathbf{q}}/N_{\mathrm{m}}$,
and the coefficients are given in Appendix A. Equations (5) describe the time
evolution of the strongly coupled system with multiple photon modes, capturing
the effect of fluctuations. Solving these equations yields the time evolution
of the photon mode populations, discussed in Section III.

Figure 2: Mode competition for weak coupling, $\Omega_{R}=0.1$eV. (a) Pump
dependence curves and (b) steady-state occupations at different $k$ modes. The
mode occupations in (b) are shown for selected pumps. The time-integrated
threshold curves are shown in Appendix B. On the right: (Top and middle row)
time evolution of photon population at different modes and (bottom row)
population inversion of the two-level system as a function of pump strength.
The pump strengths are (c-e) $\Gamma_{\uparrow}/\Gamma_{\downarrow}=0.30$,
(f-h) $0.60$, and (i-k) $0.70$. All other parameters are provided in Table 1.

###  II.3 Photoluminescence and gain spectra

When coupling between the cavity modes and the emitters is strong, polariton
lasing replaces photon lasing. In this limit, understanding the nature of mode
competition becomes more challenging, since both the gain profile and the
normal mode energies must be found self consistently, and change with pump
strength. In this section we present the numerical approach we will use to
extract both quantities, as used in Section IV.

The energies of the polariton modes can be found by computing the
photoluminescence spectrum, which is given by

$S_{\mathbf{k}}(\nu)=\int_{-\infty}^{\infty}dt\langle
a_{\mathbf{k}}^{\dagger}(t)a_{\mathbf{k}}(0)\rangle e^{i\nu t}.$ (6)

Here, we utilize the quantum regression theorem to calculate the two-time
correlators Breuer and Petruccione (2002). Solving $S_{\mathbf{k}}(\nu)$
requires the steady-state density matrix $\rho_{\textrm{ss}}$ which is used to
construct $\tilde{\rho}_{\mathbf{k}}(0)=a_{\mathbf{k}}\rho_{\textrm{ss}}$.
Time-evolving the effective density matrix
$\tilde{\rho}_{\mathbf{k}}(t)=e^{t\mathcal{L}}\tilde{\rho}_{\mathbf{k}}(0)$
and evaluating
$\text{Tr}\big{[}a_{\mathbf{k}}^{\dagger}\tilde{\rho}_{\mathbf{k}}(t)\big{]}$
gives the coupled differential equations for the two-time correlators Scully
and Zubairy (1997):

$\displaystyle\partial_{t}\langle
a^{\dagger}_{\mathbf{k}}(t)a_{\mathbf{k}}(0)\rangle$
$\displaystyle=\left(i\omega_{k}-\frac{\kappa}{2}\right)\langle
a^{\dagger}_{\mathbf{k}}(t)a_{\mathbf{k}}(0)\rangle+iN_{\mathrm{m}}B_{i}^{\ast}c_{i}^{\mathbf{k}}(t),$
$\displaystyle\partial_{t}c_{i}^{\mathbf{k}}(t)$
$\displaystyle=\xi_{ij}c_{j}^{\mathbf{k}}(t)+2f_{ijp}B_{j}\ell_{p}\langle
a^{\dagger}_{\mathbf{k}}(t)a_{\mathbf{k}}(0)\rangle,$ (7)

where
$c_{i}^{\mathbf{k}}(t)=\frac{1}{N_{\mathrm{m}}}\sum_{n}e^{i\mathbf{k}\cdot\mathbf{r}_{n}}\text{Tr}\left[\lambda_{i_{+}}^{(n)}\tilde{\rho}_{\mathbf{k}}(t)\right]$
and

$\displaystyle\langle
a^{\dagger}_{\mathbf{k}}(t)a_{\mathbf{k}}(0)\rangle=\text{Tr}\left[a^{\dagger}_{\mathbf{k}}e^{t\mathcal{L}}a_{\mathbf{k}}\rho_{ss}\right].$

In the steady state, $\ell_{p}$ becomes constant. The linear equations can be
written in a matrix form as
$\partial_{t}\mathbf{C_{k}}=\mathcal{M}\mathbf{C_{k}}$, where the vector
$\mathbf{C_{k}}=[\langle
a^{\dagger}_{\mathbf{k}}(t)a_{\mathbf{k}}(0)\rangle,\\{c_{i}^{\mathbf{k}}(t)\\}_{i}]$
and the matrix $\mathcal{M}$ are obtained from Eqs. (7). The Fourier transform
can be written as

$\displaystyle S_{\mathbf{k}}(\nu)$
$\displaystyle=\int_{-\infty}^{\infty}e^{i\nu t+\mathcal{M}t}\mathbf{C}(0)dt$
$\displaystyle=(i\nu+\mathcal{M})^{-1}\mathbf{C}(0)=\sum_{i}\frac{\alpha_{i}}{\mu_{i}+i\nu}\left|r_{i}\right\rangle,$
(8)

where $\mu_{i}$ is the eigenvalue of $\mathcal{M}$ that corresponds to the
right $\left|r_{i}\right\rangle$ and left $\left\langle l_{i}\right|$
eigenvectors. The coefficient $\alpha_{i}$ is then given by
$\alpha_{i}=\left\langle l_{i}\right|\mathbf{C}(0)\rangle/\left\langle
l_{i}\right|r_{i}\rangle$.

The “gain” (i.e., emission minus absorption) can be studied by considering
two-time correlations of the molecules:

$\displaystyle G(\nu)$ $\displaystyle=\int_{-\infty}^{\infty}dte^{-i\nu
t}\langle[\sigma^{+(n)}(t),\sigma^{-(n)}(0)]\rangle$
$\displaystyle=V_{i}^{+}V_{j}^{-}\int_{-\infty}^{\infty}dte^{-i\nu
t}\langle[\lambda_{i}^{(n)}(t),\lambda_{j}^{(n)}(0)]\rangle,$ (9)

where $V_{i}^{\pm}=\frac{1}{2}\text{Tr}(\sigma^{\pm}\lambda_{i})$. In the
weak-coupling picture, this has a simple interpretation as the frequency-
resolved gain the photon modes see; at strong coupling the picture is more
complicated, but this nonetheless allows one to probe the molecular dynamics
separately from the photons.

By also considering the time-dependent correlations between different
molecules, one can further compute the momentum-resolved gain spectrum:

$\displaystyle G_{\mathbf{k}}(\nu)=$ $\displaystyle
V_{i}^{+}V_{j}^{-}\int_{-\infty}^{\infty}dte^{-i\nu t}$ (10)
$\displaystyle\sum_{n,m}e^{i\mathbf{k}\cdot(\mathbf{r}_{n}-\mathbf{r}_{m})}\langle[\lambda_{i}^{(n)}(t),\lambda_{j}^{(m)}(0)]\rangle.$

This quantity identifies the different gain seen by each cavity mode, and so
helps explain when mode switching occurs. Note that we can recover the form of
$G(\nu)$ in Eq. (II.3) by summing $G_{\mathbf{k}}(\nu)$ over all $\mathbf{k}$.

##  III Multimode lasing dynamics with continuous pump

###  III.1 Weak-coupling regime

Using the methods described in the previous section, we now explore the
dynamics of mode switching, starting with the weak-coupling regime, for which
we set $\Omega_{R}=0.1$eV. As shown in Fig. 2(a,b) (similar results also seen
in our previous work Arnardottir  _et al._ (2020)), for the parameters we use
the lasing mode switches between $k=5$ and $k=0$ as a function of pump
strength. Upon increasing the pump strength, the mode $k=5$, which is closest
to resonance with the $(1-0)$ transition, starts lasing first. When the pump
is further increased, the $k=0$ starts to lase, and the lasing at $k=5$ is
suppressed. As noted in Ref. Arnardottir  _et al._ (2020), the selection of
which mode wins mode competition in the weak coupling regime can be attributed
to the different threshold and slope efficiency—i.e., gradient of photon
population vs pump—of each mode. The modes near the vibrational shoulder
($k=4,5,6$) have the lowest threshold, but the mode near $k=0$ has the highest
slope efficiency (better visible in linear-scale threshold curves, see
Arnardottir  _et al._ (2020)). Thus at low pumping (low inversion), lasing
occurs at $k=5$, whereas at high pumping (high inversion), the mode $k=0$ with
the higher slope efficiency wins.

The time evolution of the system, as shown in Fig. 2(c-k), reveals that mode
switching occurs not only as a function of pump strength, but also in time.
Here, the term ‘switching’ corresponds to the time instant when the population
of one lasing mode surpasses the population of another mode. Fig. 2(c,f,i)
show for illustrative purpose the data shown in Fig. 2(d,g,j), which presents
the population of selected $k$ modes. Fig. 2(e,h,k) shows the population
inversion as a function of time.

At a low pump strength, in Fig. 2(c-d), lasing starts and remains at $k=5$
until steady state is reached. At higher pump strengths, as shown in Fig.
2(f-g) and Fig. 2(i-j), switching to $k=0$ mode occurs. At a high pump
strength, Fig. 2(i-j), lasing at $k=0$ completely suppresses the lasing at
$k=5$ after the switch. Only in a very narrow range of pump strengths, as
illustrated in Fig. 2(f-g), does the system show bi-modal lasing. In that
case, the $k=0$ mode slowly grows until it marginally prevails, and both modes
sustain a macroscopic population at steady state. The dashed lines in Fig.
2(d,g,j) show the photon population summed over all $k$ modes. It is
interesting to note that the sum over all $k$ modes shows only a single
threshold; whenever switching takes place the photon population is re-
distributed between different modes. One may note also that the population
inversion (as shown in Fig. 2(e,h,k)) appears to reach a stationary value
before switching occurs. After the switch, the inversion does however show a
small kink.

Figure 3: Metastable decay rate of the $k=5$ mode as a function of pump
strength, in the vicinity of the transition point. The decay rate is the
inverse of the time at which the lasing mode switch occurs, i.e., when the
population of mode $k=0$ exceeds that of mode $k=5$. Parameters as in Fig. 2.
Figure 4: Mode competition for strong coupling, $\Omega_{R}=0.4$eV. (a) Pump
dependence curves and (b) steady-state occupations at different $k$ modes. The
mode occupations in (b) are shown for selected pumps. The time-integrated
threshold curves are shown in Appendix B. (Top and middle row) time evolution
of photon population at different modes and (bottom row) population inversion
of the two-level system as a function of pump strength. The pump strengths are
(c-e) $\Gamma_{\uparrow}/\Gamma_{\downarrow}=0.60$, (f-h) $1.30$, and (i-k)
$1.50$. All other parameters are provided in Table 1.

In Figure 3, we show the rate at which the metastable $k=5$ state is replaced
by $k=0$. We see that the rate vanishes as one approaches the transition
point, i.e., the pump strength where bi-modal lasing is observed
($\sim\leavevmode\nobreak\ 0.6\Gamma_{\downarrow}$, Fig. 2(f-h)).

###  III.2 Strong coupling regime

At strong coupling ($\Omega_{R}=0.4$eV) we observe shifting of the lasing mode
towards higher $k$ on increasing the pump strength. As shown by the threshold
curves and spectra in Fig. 4(a-b), lasing is first triggered in the $k=5$ mode
but switches to $k=7$ at high pump strengths. In a narrow regime of pump
strengths, the $k=6$ mode has the highest occupation.

The time evolution of the mode occupations for strong coupling is presented in
Fig. 4(c-k). At a low pump strength, in Fig. 4(c-d), lasing starts at $k=5$
and switches to a lower $k=4$ mode. At higher pump strengths, in Fig. 4(f-g)
and Fig. 4(i-j), the lasing mode switches to higher $k$ modes. One notable
distinction from weak coupling is seen in the time evolution of the population
inversion. As seen in Fig. 4(h,k), the population inversion continues to grow
until the (final) mode switching occurs.

Both in the strong and weak coupling regimes, we have referred to the lasing
modes with the absolute value of the momentum vector, $k=|\mathbf{k}|$. Within
our model, as we have no source of disorder, all modes with a given
$|\mathbf{k}|$ are degenerate. In real experiments, the disorder that will be
present is likely to break the perfect degeneracy, leading to a particular
angular pattern being selected for the lasing mode.

##  IV Gain spectrum

To better understand why mode switching occurs in response to pumping, we next
discuss how gain and photoluminescence spectra change as a function of pump
power. To avoid the subtleties associated with time-dependent spectra Eberly
and Wódkiewicz (1977), we only consider the steady-state spectra.

###  IV.1 Weak coupling regime

Figure 5: Momentum-integrated gain spectra $G(\nu)$ near the lasing threshold
for weak coupling, $\Omega_{R}=0.1$eV. The top panel shows $G(\nu)$, while the
bottom panel shows $-G(\nu)$, both on a logarithmic scale. The pump strengths
are color coded as in the legend and the inset which shows the corresponding
part of the threshold curve. Solid lines in the panels (and filled circles in
the inset) are results above the lasing threshold. Threshold is around 0.21
$\Gamma_{\uparrow}/\Gamma_{\downarrow}$, marked by the vertical line in the
inset. All other parameters as in Table 1.  Figure 6: Momentum-integrated gain
spectra $G(\nu)$ as a function of pump strength for the weak coupling case
$\Omega_{R}=0.1$eV. This figure shows an expanded range of energies and pump
strengths compared to Fig. 5.

To introduce the properties of the gain spectrum, we start by considering the
total (momentum-integrated) spectrum $G(\nu)$ of the molecules defined in Eq.
(II.3). Figure. 5 shows $G(\nu)$ for weak coupling, at a narrow range of pump
strengths (see legend) near the initial lasing threshold ($\sim$ 0.21
$\Gamma_{\uparrow}/\Gamma_{\downarrow}$). Since $G(\nu)$ can be either
positive (net gain from the molecules) or negative (net absorption), we show
both $G(\nu)$ and $-G(\nu)$ on logarithmic scales. Below threshold (dashed
lines) $G(\nu)$ is negative at all frequencies. As the pump strength increases
a region of positive net gain develops. Positive gain from the molecules does
not immediately lead to lasing, as it is not sufficient to overcome the cavity
losses. When the gain at the frequency of the $k=5$ cavity mode (around 0.75
eV) becomes sufficiently high to exceed cavity loss, lasing occurs. Figure 6
then shows how the gain spectrum continues to develop, considering a wider
range of frequencies (covering all relevant cavity modes) and a wider range of
pump strengths. From this figure, one sees that between pump strengths 0.60
and 0.70 $\Gamma_{\uparrow}/\Gamma_{\downarrow}$ the peak of the gain switches
to a lower energy mode at 0.63-0.64 eV.

We can further study how the gain affects different cavity modes by
considering $G_{\mathbf{k}}(\nu)$. The momentum-resolved gain spectra are
shown in Fig. 7(a-c) for the three selected pumps that were studied in Section
III (Fig. 2). Here only the positive gain values are plotted. The results
confirm that the peak of the gain shifts from the mode $k=5$ to $k=0$ as the
pump is increased, and there is a narrow region of pump strengths where the
gain is similar for the two modes. We calculated also the photoluminescence
spectrum $S_{\mathbf{k}}(\nu)$ to see how the mode energies change upon
increasing the pump strength. As shown in Fig. 7(d-f), under weak coupling the
modes are only slightly shifted from the bare photon dispersion, and so in
this case the mode locations remain practically unchanged as the pump is
increased.

Figure 7: (a-c) Momentum-resolved gain spectra $G_{\mathbf{k}}(\nu)$ and (d-f)
photoluminescence spectra $S_{\mathbf{k}}(\nu)$ for selected pump strengths
for the weak coupling case, $\Omega_{R}=0.1$eV. In (a-c) the numbers next to
the peaks label the $k$ of the mode(s) that are lasing in the steady state. In
(d-f) the white crosses mark the uncoupled cavity modes. The pump strengths
are (a,d) $0.30\leavevmode\nobreak\ \Gamma_{\uparrow}/\Gamma_{\downarrow}$,
(b,e) $0.60\leavevmode\nobreak\ \Gamma_{\uparrow}/\Gamma_{\downarrow}$, and
(c,f) $0.70\leavevmode\nobreak\ \Gamma_{\uparrow}/\Gamma_{\downarrow}$; all
other parameters as in Table 1.

###  IV.2 Strong coupling regime

In the strongly coupled system one should now consider lasing arising from
scattering into the polariton modes. As discussed elsewhere Mazza  _et al._
(2013); Strashko  _et al._ (2018), one can still consider a gain spectrum for
feeding into the polariton mode. This gain spectrum corresponds to the
reservoir of dark excitonic states. At strong coupling, the gain spectrum
shows more complex changes as a function of pump strength. Figure 9 shows the
evolution of total gain $G(\nu)$ with pump strength. As before, below
threshold (dashed line) the gain is negative, while above threshold one sees
various regions with net gain. As pumping is further increased, the net gain
moves toward higher energies. The shift can be explained by considering two
effects that arise due to strong coupling and pumping.

First, gain is increased at higher energies as a result of the vibrational
sidebands becoming successively inverted with stronger pumping. To understand
this we must consider both the relative populations of the ground and excited
electronic manifolds, and the populations of the vibrational sublevels
labelled by $n$ as illustrated in Fig. 8. Considering the ground state
electronic manifold, we may state the following about the populations of the
sublevels. Unless pumping becomes very strong, the vibrational sublevels will
be populated thermally (according to a Boltzmann distribution), so the $n=0$
sublevel is most populated and the higher $n$ are less populated. Turning now
to the population of the electronic excited states, we may note that the
population of the excited state manifold is zero at weak pumping and grows
with pumping.

Therefore, with increasing pumping the population of the excited state
manifold will start to exceed the populations of the ground-state-manifold
sublevels, such that the first sideband to get inverted corresponds to
transitions from the 0th state of the excited manifold to the highest $n$
state of the ground state manifold (the ($n-0$) transition), see Fig. 8. That
first inverted transition is also the transition involving the smallest energy
difference, thus it leads to gain at lower energies. As the pumping increases
further, sidebands with successively smaller $n$ additionally become inverted,
thus one starts to get gain at higher energies. At the same time the excited
state sublevels with $m>0$ can also become inverted. Both these processes mean
($n-m$) transitions with larger energy difference are inverted, shifting the
gain toward higher energy Strashko  _et al._ (2018).

Figure 8: Schematic cartoon of successive inversion of sublevel transitions
when vibrational sublevels are thermally populated. Blue and red solid lines
represent the $m$ and $n$ levels for transitions $n-m$ from excited state
manifold ($m$) to the ground state manifold ($n$) of an organic emitter. As
pump strength is increased, the first transition to undergo population
inversion is the $4-0$ transition, which corresponds to transitions from the
0th state of the excited state manifold to the 4th state of the ground state
manifold.

The second effect to be considered is that, not only the gain but also the
energies of the strongly-coupled normal modes evolve with pump strength. As
shown by the momentum-resolved results in Fig. 10(a-c), for instance when the
pump strength is increased from 1.30 to $1.50\leavevmode\nobreak\
\Gamma_{\uparrow}/\Gamma_{\downarrow}$, the gain spectrum moves from 0.78–0.82
eV to 0.84–0.86 eV. The photoluminescence spectra in Fig. 10(d-f) show that
the dips associated with the polaritonic mode frequencies move to higher
energy as the whole dispersion blue shifts upon increasing pump. The blue
shift originates from saturation of the two-level systems Arnardottir  _et
al._ (2020), reducing the effective coupling strength, ultimately recovering
the uncoupled photon mode dispersion at very strong pumping.

Summarizing, the shifting of the lasing mode toward higher energies upon
increasing pumping in the strong coupling case is caused by combined effects
of the gain moving toward higher energies, as well as the blueshifting of the
polariton mode frequencies.

Figure 9: Total gain spectra $G(\nu)$ as a function of pump strength for
strong coupling, $\Omega_{R}=0.4$eV. As in Figs. 5 and 6, the two panels show
$G(\nu)$ and $-G(\nu)$ on a logarithmic scale. The pump strengths are color
coded according to the legend, and all other parameters are as in Table 1.
Figure 10: (a-c) Momentum-resolved gain and (d-f) photoluminescence spectra
for selected pump strengths for the strong coupling case $\Omega_{R}=0.4$eV.
(a-c) The numbers next to the gain peaks label the $k$ of the mode(s) that are
lasing in the steady state. (d-f) The white crosses mark the uncoupled cavity
modes. The color scale is logarithmic. The pump strengths are (a,d)
$0.60\leavevmode\nobreak\ \Gamma_{\uparrow}/\Gamma_{\downarrow}$, (b,e)
$1.30\leavevmode\nobreak\ \Gamma_{\uparrow}/\Gamma_{\downarrow}$, and (c,f)
$1.50\leavevmode\nobreak\ \Gamma_{\uparrow}/\Gamma_{\downarrow}$.

##  V Multimode lasing dynamics with pulsed pump

Next we consider the case of pulsed excitation, which is relevant for typical
experiments on organic exciton-polariton systems. We excite the system with a
4 ps pulse that has a Gaussian time dependence, and study the dynamics of mode
switching as a function of pump strength (amplitude). All other parameters are
as in Section III.

###  V.1 Weak coupling regime

In the weak coupling regime, at low pump strength (Fig. 11(c-d)), lasing
occurs at $k=5$ until decay causes lasing to cease. At higher pump strengths,
as shown in Fig. 11(f-g) and Fig. 11(i-j), the lasing mode switches from $k=5$
to $k=0$ similar to the continuous pump case discussed in Section III.
However, here the lasing mode switches from $k=0$ back to $k=5$ as the
population decays. Throughout the studied pump range, the higher $k=6$ mode
gains large population at the onset of lasing. This is also visible in the
time-integrated threshold curve and populations in Fig. 11(a-b): occupation of
the $k=6$ mode follows quite closely that of the $k=5$ mode. As shown by Fig.
11(h,k), at higher pump strengths the two-level systems have very high
inversion before the onset of lasing, which is likely connected to excess
gain, enough to trigger lasing in several modes.

As noted above, considering time-dependent gain spectra present some issues
with interpretation Eberly and Wódkiewicz (1977). The steady-state spectra
presented in Section IV can however provide some intuition on how the gain and
photoluminescence might shift under time-dependent population inversion. For
example, at weak coupling, the peak of the gain switched from $k=5$ to $k=0$
mode when the (continuous) pump strength was increased from 0.30 to 0.70. Here
in the pulsed case, the lasing peak switches in a qualitatively similar way
(both as a function of pump strength as well as in time).

Figure 11: Mode competition with pulsed excitation for weak coupling,
$\Omega_{R}=0.1$eV. (a) Time-integrated pump dependence curves and (b)
occupations at different $k$ modes for selected pumps. On the right (c-k):
(Top and middle row) time evolution of photon population at different modes
and (bottom row) population inversion of the two-level system at various pump
strengths. The excitation pulse has a Gaussian time dependence centered at
around 16 ps with full width at half maximum of around 4 ps. The pulse profile
is shown in (e,h,k) by the red dashed line (in arbitrary units). The peak pump
strengths are (c-e) $\Gamma_{\uparrow}/\Gamma_{\downarrow}=5.00$, (f-h)
$20.0$, and (i-k) $30.0$. All other parameters are provided in Table 1. Note
that here the time axis is in linear scale.

###  V.2 Strong coupling regime

At strong coupling, for pulsed excitation we observe again rich switching
dynamics. At low pump strength, as shown in Fig. 12(c-d), lasing starts at
$k=4$ mode and switches to a higher $k=5$ mode before the population decays.
At higher pump strengths, lasing starts in the high $k=7$ mode and undergoes
switching via $k=5$ to $k=4$ and back to $k=5$. The duration of the $k=7$
lasing coincides with the time that $\langle\sigma_{z}\rangle$, shown in Fig.
12(h,k), remains clamped close to $0$ by polariton lasing. This duration grows
upon increasing the pump strength. The switching sequence to the lower $k$
modes, 5–4–5, occurs during the decay of $\langle\sigma_{z}\rangle$ in a time
scale that we observe to be roughly independent of the pump strength (we
confirmed this up to pump amplitude
$\Gamma_{\uparrow}/\Gamma_{\downarrow}=150$). This time scale is around 40 ps
which corresponds to the decay times of spontaneous emission and cavity modes
(Table 1).

Similar to the weak coupling case, for a short time at the onset of lasing,
multiple modes are highly occupied, which contributes to the time-integrated
populations shown in Fig. 12(a-b).

Figure 12: Mode competition with pulsed excitation for strong coupling,
$\Omega_{R}=0.4$eV. (a) Time-integrated pump dependence curves and (b)
occupations at different $k$ modes for selected pump powers. On the right
(c-k): (Top and middle row) time evolution of photon population at different
modes and (bottom row) population inversion of the two-level system at various
pump strengths. The excitation pulse has a Gaussian time dependence centered
at around 16 ps with full width at half maximum of around 4 ps. The pulse
profile is shown in (e,h,k) by the red dashed line (in arbitrary units). The
peak pump strengths are (c-e) $\Gamma_{\uparrow}/\Gamma_{\downarrow}=5.00$,
(f-h) $20.0$, and (i-k) $30.0$. All other parameters are provided in Table 1.
Note that here the time axis is in linear scale.

##  VI Discussion

By simulating the time evolution of the photon population, we have explored
the dynamics of mode switching, behavior which could not be seen if only the
steady-state result is studied. For example, in Fig. 2(c,f), the steady-state
result would show single mode ($k=0$) lasing whereas the time evolution
reveals that a higher-energy mode ($k=5$) is lasing at earlier times.
Furthermore, when comparing simulations to experiments dealing with fast
dynamics, the experimental measurement is often restricted to observing the
time-integrated spectrum, which may lead to incorrect interpretations of the
properties of the system (see Appendix B). Especially in experiments with fast
pulsed excitation, it is common to record time-integrated spectra. Here we
have shown that for pulsed excitation the lasing mode can temporally switch
such that the time-integrated spectra will exhibit several lasing peaks (or a
broadened peak with sub-structure, depending on the mode spacing and
resolution of the measurement apparatus). Resolving the time evolution is
therefore important in describing the behavior of such systems both in
simulations and in experiments.

From the gain calculations we can generally conclude that lasing occurs at a
mode (or modes) where there is a sharp positive peak in the gain spectrum. The
results suggest that if several modes have positive gain, the mode with the
highest positive gain peak will win out the lasing mode competition. In the
weak coupling case the highest peak in the gain spectrum simply switches from
a higher $k=5$ mode into the $k=0$ mode, given by the different thresholds and
slope efficiencies of each mode, In the strongly coupled system, the observed
shift of the lasing mode towards higher $k$ modes originates from the shifting
of both the gain and the polariton modes toward higher energies as pump is
increased.

We believe that observing mode competition effects such as those shown here is
realistic in current experiments. In this paper we used a loss rate of 0.1 meV
for both the cavity modes and the excitons. This corresponds to lifetimes of
$\sim 40$ps (see Table 1 for all simulation parameters). Under continuous
pumping, the lasing starts at around $t=100$ ps ($t=10$ ps) in the weak
(strong) coupling case. As shown in Fig. 3, the time that it takes for the
mode switching to occur depends strongly on the pump strength. Nevertheless,
we can take examples from our data. In the weak (strong) coupling case with
continuous pumping, at pump strength
$\Gamma_{\uparrow}/\Gamma_{\downarrow}=0.70$ ($1.30$), the switching occurs at
around $t=1000$ ps ($t=100$ ps), i.e., after 10 times the time required for
lasing to start, or 25 times (2.5 times) the cavity mode lifetime. Another
example can be drawn from the pulsed excitation scheme, where the lasing mode
is shown to undergo several switching events within 80–120 ps, or 2–3 times
the cavity lifetime, from the arrival of the excitation pulse. The temporal
switching of lasing modes could be resolved by, e.g., a streak camera.

We note that the cavity lifetime used in the simulations is rather long for an
organic polariton system. Cavity lifetimes of the order 10–100 ps are typical
for inorganic polariton systems, whereas in organic systems, the lifetimes are
typically in the picosecond scale Keeling and Kéna-Cohen (2020). The choice of
simulation parameters in the present study is based on our earlier works
Strashko  _et al._ (2018); Arnardottir  _et al._ (2020), and we stress here
that the results are similar for shorter lifetimes.

Multimode polariton lasing could potentially be of use for technological
applications. Mode switching phenomena could be used in optical information
processing in a similar fashion as polariton switches Gao  _et al._ (2012);
Amo  _et al._ (2010) and transistors Zasedatelev  _et al._ (2019); Ballarini
_et al._ (2013) (based on on/off-switching of the output intensity of a
single-mode condensate) have been proposed. Furthermore, when a system is at a
point where multiple modes are lasing in the steady state, even a small
disturbance in the operating conditions may push the system from one
configuration to another, which could be utilized in sensing applications.

###### Acknowledgements.

AJM and PT acknowledge support by the Academy of Finland under project numbers
303351, 307419, 327293, 318987 (QuantERA project RouTe), 318937 (PROFI), and
320167 (Flagship Programme, Photonics Research and Innovation (PREIN)), and by
Centre for Quantum Engineering (CQE) at Aalto University. AJM acknowledges
financial support by the Jenny and Antti Wihuri Foundation and ETH Zurich
Postdoctoral Fellowship. KBA and JK acknowledge financial support from EPSRC
program “Hybrid Polaritonics” (EP/M025330/1). KBA acknowledges support from
The RSE Saltire Research Award.

##  Appendix A Details of generalized Gell-Mann matrices

In this appendix we provide mathematical details of the generalized Gell-Mann
matrix representation of the molecules.

Gell-Mann matrices which have previously been used, e.g., in describing the
strong interaction in particle physics Gell-Mann (1962); Stone and Goldbart
(2009). Their commutation and product rules are given by

$\displaystyle\lambda_{i}=\lambda_{i}^{\dagger},\qquad\text{Tr}\lambda_{i}=0,\qquad\text{Tr}[\lambda_{i}\lambda_{j}]=\delta_{ij},$
$\displaystyle[\lambda_{i},\lambda_{j}]=2if_{ijk}\lambda_{k},\qquad\lambda_{i}\lambda_{j}=\zeta_{ijk}\lambda_{k}+\frac{2}{N_{\lambda}}\delta_{ij},$

where we may write $\zeta_{ijk}=if_{ijk}+t_{ijk}$ and $N_{\lambda}$ is the
dimension of the matrices. The tensors $f_{ijk}$ and $t_{ijk}$ are the
symmetric and antisymmetric structure constants, respectively. Gell-Mann
matrices allow for expanding any linear $N_{\lambda}\times N_{\lambda}$
operator in this basis as

$A=\frac{1}{N_{\lambda}}\text{Tr}[A]\mathbbm{1}_{N_{\lambda}}+\frac{1}{2}\text{Tr}[A\lambda_{i}]\lambda_{i}.$
(11)

Using the above results, we may identify the coefficients in the second-order
cumulant equations (Eq. (5)) as follows:

$\displaystyle\phi_{i}=\frac{2i}{N_{\mathrm{v}}}f_{ipr}\sum_{\mu}\gamma_{p}^{\mu}\gamma_{r}^{\mu\ast},$
$\displaystyle\beta_{ij}=\frac{1}{2}B_{p}(f_{i_{x}pj}-if_{i_{y}pj}),$
$\displaystyle X_{ij}=\xi_{i_{x}j_{x}}+i\xi_{i_{y}j_{x}},$
$\displaystyle\xi_{ij}=2f_{ipj}A_{p}+i\left(f_{ips}\zeta_{rsj}+f_{ris}\zeta_{spj}\right)\sum_{\mu}\gamma_{p}^{\mu}\gamma_{r}^{\mu\ast}.$

##  Appendix B Time-integrated threshold curves

In this appendix we present threshold curves and population distributions
obtained by time-integrating the photon mode populations until steady state
has been reached. The results shown in Fig. 13 represent a common case in
experimental measurements of a system with fast dynamics, where, instead of
having access to the temporal evolution of the photon populations, mode
populations are collected by integrating over a finite time. Our results show
that observing multiple peaks in the time-integrated mode populations can
refer to temporal switching between modes.

Figure 13: Time-integrated results with continuous excitation for weak and
strong coupling. (a,c) Pump dependence curves and (b,d) occupations at
different $k$ modes. The mode occupations in (b,d) are shown for selected pump
powers. The weak coupling case ($\Omega_{R}=0.1$eV) is shown on the top row
and the strong coupling case ($\Omega_{R}=0.4$eV) on the bottom row.

## References

  * Furfaro  _et al._ (2004) L. Furfaro, F. Pedaci, M. Giudici, X. Hachair, J. Tredicce, and S. Balle,  _Mode-switching in semiconductor lasers_ , IEEE J. Quantum Electron. 40, 1365 (2004).
  * Pedaci  _et al._ (2005) F. Pedaci, M. Giudici, J. Tredicce, and G. Giacomelli,  _Experimental analysis of mode-hopping in bulk semiconductor lasers_ , Appl. Phys. B 81, 993 (2005).
  * Shore and Kane (2005) K. A. Shore and D. M. Kane,  _Unlocking dynamical diversity: optical feedback effects on semiconductor lasers_ (John Wiley & Sons, Chichester, 2005).
  * Gordon  _et al._ (2008) A. Gordon, C. Y. Wang, L. Diehl, F. X. Kärtner, A. Belyanin, D. Bour, S. Corzine, G. Höfler, H. C. Liu, H. Schneider, T. Maier, M. Troccoli, J. Faist, and F. Capasso,  _Multimode regimes in quantum cascade lasers: From coherent instabilities to spatial hole burning_ , Phys. Rev. A 77, 053804 (2008).
  * Wright  _et al._ (2017) L. G. Wright, D. N. Christodoulides, and F. W. Wise,  _Spatiotemporal mode-locking in multimode fiber lasers_ , Science 358, 94 (2017).
  * Herre and Barabas (1989) P. Herre and U. Barabas,  _Mode switching of Fabry-Perot laser diodes_ , IEEE J. Quantum Electron. 25, 1794 (1989).
  * Lenstra and Yousefi (2014) D. Lenstra and M. Yousefi,  _Rate-equation model for multi-mode semiconductor lasers with spatial hole burning_ , Opt. Express 22, 8143 (2014).
  * Valle  _et al._ (1995) A. Valle, J. Sarma, and K. A. Shore,  _Dynamics of transverse mode competition in vertical cavity surface emitting laser diodes_ , Opt. Commun. 115, 297 (1995).
  * Law and Agrawal (1997) S. Law and G. Agrawal,  _Effects of spatial hole burning on gain switching in vertical-cavity surface-emitting lasers_ , IEEE J. Quantum Electron. 33, 462 (1997).
  * Yacomotti  _et al._ (2004) A. M. Yacomotti, L. Furfaro, X. Hachair, F. Pedaci, M. Giudici, J. Tredicce, J. Javaloyes, S. Balle, E. A. Viktorov, and P. Mandel,  _Dynamics of multimode semiconductor lasers_ , Phys. Rev. A 69, 053816 (2004).
  * Deng  _et al._ (2010) H. Deng, H. Haug, and Y. Yamamoto,  _Exciton-polariton Bose-Einstein condensation_ , Rev. Mod. Phys. 82, 1489 (2010).
  * Keeling and Berloff (2011) J. Keeling and N. G. Berloff,  _Exciton–polariton condensation_ , Contemp. Phys. 52, 131 (2011).
  * Keeling and Kéna-Cohen (2020) J. Keeling and S. Kéna-Cohen,  _Bose–Einstein Condensation of Exciton-Polaritons in Organic Microcavities_ , Annu. Rev. Phys. Chem. 71, 435 (2020).
  * Sun  _et al._ (2018) Y. Sun, Y. Yoon, S. Khan, L. Ge, M. Steger, L. N. Pfeiffer, K. West, H. E. Türeci, D. W. Snoke, and K. A. Nelson,  _Stable switching among high-order modes in polariton condensates_ , Phys. Rev. B 97, 045303 (2018).
  * Kusudo  _et al._ (2013) K. Kusudo, N. Y. Kim, A. Löffler, S. Höfling, A. Forchel, and Y. Yamamoto,  _Stochastic formation of polariton condensates in two degenerate orbital states_ , Phys. Rev. B 87, 214503 (2013).
  * Kim  _et al._ (2016) S. Kim, B. Zhang, Z. Wang, J. Fischer, S. Brodbeck, M. Kamp, C. Schneider, S. Höfling, and H. Deng,  _Coherent Polariton Laser_ , Phys. Rev.X 6, 011026 (2016).
  * Maragkou  _et al._ (2010) M. Maragkou, A. J. D. Grundy, E. Wertz, A. Lemaître, I. Sagnes, P. Senellart, J. Bloch, and P. G. Lagoudakis,  _Spontaneous nonground state polariton condensation in pillar microcavities_ , Phys. Rev. B 81, 081307 (2010).
  * Sawicki  _et al._ (2021) K. Sawicki, T. J. Sturges, M. Ściesiek, T. Kazimierczuk, K. Sobczak, A. Golnik, W. Pacuski, and J. Suffczyński,  _Polariton lasing and energy-degenerate parametric scattering in non-resonantly driven coupled planar microcavities_ , Nanophotonics 10, 2421 (2021).
  * Doan  _et al._ (2005) T. D. Doan, H. T. Cao, D. B. Tran Thoai, H. Haug, D. Thoai, and H. Haug,  _Condensation kinetics of microcavity polaritons with scattering by phonons and polaritons_ , Phys. Rev. B 72, 085301 (2005).
  * Doan  _et al._ (2006) T. Doan, H. Cao, D. Thoai, and H. Haug,  _Microcavity polariton kinetics for bosonic condensation and lasing in II-VI compound materials_ , Phys. Rev. B 74, 115316 (2006).
  * Lai  _et al._ (2007) C. W. Lai, N. Y. Kim, S. Utsunomiya, G. Roumpos, H. Deng, M. D. Fraser, T. Byrnes, P. Recher, N. Kumada, T. Fujisawa, and Y. Yamamoto,  _Coherent zero-state and $\pi$-state in an exciton–polariton condensate array_, Nature 450, 529 (2007).
  * Kim  _et al._ (2011) N. Y. Kim, K. Kusudo, C. Wu, N. Masumoto, A. Löffler, S. Höfling, N. Kumada, L. Worschech, A. Forchel, and Y. Yamamoto,  _Dynamical d-wave condensation of exciton–polaritons in a two-dimensional square-lattice potential_ , Nat. Phys. 7, 681 (2011).
  * Jacqmin  _et al._ (2014) T. Jacqmin, I. Carusotto, I. Sagnes, M. Abbarchi, D. Solnyshkov, G. Malpuech, E. Galopin, A. Lemaître, J. Bloch, and A. Amo,  _Direct Observation of Dirac Cones and a Flatband in a Honeycomb Lattice for Polaritons_ , Phys. Rev. Lett. 112, 116402 (2014).
  * Amo and Bloch (2016) A. Amo and J. Bloch,  _Exciton-polaritons in lattices: A non-linear photonic simulator_ , C. R. Phys. 17, 934 (2016).
  * Scafirimuto  _et al._ (2021) F. Scafirimuto, D. Urbonas, M. A. Becker, U. Scherf, R. F. Mahrt, and T. Stöferle,  _Tunable exciton–polariton condensation in a two-dimensional Lieb lattice at room temperature_ , Commun. Phys. 4, 39 (2021).
  * Boulier  _et al._ (2020) T. Boulier, M. J. Jacquet, A. Maître, G. Lerario, F. Claude, S. Pigeon, Q. Glorieux, A. Amo, J. Bloch, A. Bramati, and E. Giacobino,  _Microcavity Polaritons for Quantum Simulation_ , Adv. Quantum Technol. 3, 2000052 (2020).
  * Eastham (2008) P. R. Eastham,  _Mode locking and mode competition in a nonequilibrium solid-state condensate_ , Phys. Rev. B 78, 035319 (2008).
  * Wouters (2008) M. Wouters,  _Synchronized and desynchronized phases of coupled nonequilibrium exciton-polariton condensates_ , Phys. Rev. B 77, 121302 (2008).
  * Jamadi  _et al._ (2019) O. Jamadi, F. Réveret, D. Solnyshkov, P. Disseix, J. Leymarie, L. Mallet-Dida, C. Brimont, T. Guillet, X. Lafosse, S. Bouchoule, F. Semond, M. Leroux, J. Zuniga-Perez, and G. Malpuech,  _Competition between horizontal and vertical polariton lasing in planar microcavities_ , Phys. Rev. B 99, 085304 (2019).
  * Töpfer  _et al._ (2020) J. D. Töpfer, H. Sigurdsson, S. Alyatkin, and P. G. Lagoudakis,  _Lotka-Volterra population dynamics in coherent and tunable oscillators of trapped polariton condensates_ , Phys. Rev. B 102, 195428 (2020).
  * Vlaho and Eckardt (2021) M. Vlaho and A. Eckardt,  _Nonequilibrium mode competition in a pumped dye-filled cavity_ , Phys. Rev. A 104, 063709 (2021).
  * Vlaho  _et al._ (2019) M. Vlaho, H. A. M. Leymann, D. Vorberg, and A. Eckardt,  _Controlled two-mode emission from the interplay of driving and thermalization in a dye-filled photonic cavity_ , Phys. Rev.Research 1, 033191 (2019).
  * Marelic  _et al._ (2016) J. Marelic, L. F. Zajiczek, H. J. Hesten, K. H. Leung, E. Y. X. Ong, F. Mintert, and R. A. Nyman,  _Spatiotemporal coherence of non-equilibrium multimode photon condensates_ , New Journal of Physics 18, 103012 (2016).
  * Kurtscheid  _et al._ (2019) C. Kurtscheid, D. Dung, E. Busley, F. Vewinger, A. Rosch, and M. Weitz,  _Thermally condensing photons into a coherently split state of light_ , Science 366, 894 (2019).
  * Väkeväinen  _et al._ (2020) A. I. Väkeväinen, A. J. Moilanen, M. Nečada, T. K. Hakala, K. S. Daskalakis, and P. Törmä,  _Sub-picosecond thermalization dynamics in condensation of strongly coupled lattice plasmons_ , Nat. Commun. 11, 3139 (2020).
  * Arnardottir  _et al._ (2020) K. B. Arnardottir, A. J. Moilanen, A. Strashko, P. Törmä, and J. Keeling,  _Multimode Organic Polariton Lasing_ , Phys. Rev. Lett. 125, 233603 (2020).
  * Siegman (1986) A. E. Siegman,  _Lasers_ (University science books, 1986).
  * Ćwik  _et al._ (2014) J. A. Ćwik, S. Reja, P. B. Littlewood, and J. Keeling,  _Polariton condensation with saturable molecules dressed by vibrational modes_ , Eur. Lett. 105, 47009 (2014).
  * Strashko  _et al._ (2018) A. Strashko, P. Kirton, and J. Keeling,  _Organic Polariton Lasing and the Weak to Strong Coupling Crossover_ , Phys. Rev. Lett. 121, 193601 (2018).
  * Jong  _et al._ (2015) M. d. Jong, L. Seijo, A. Meijerink, and F. T. Rabouw,  _Resolving the ambiguity in the relation between Stokes shift and Huang–Rhys parameter_ , Phys. Chem. Chem. Phys. 17, 16959 (2015).
  * Note (1) We note that compared to our previous work Ref. Arnardottir  _et al._ (2020) we have corrected the sign of the $\sqrt{S}\sigma^{z}_{n}$ term in Eq. (2). This does not significantly change the results.
  * 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).
  * Sánchez-Barquilla  _et al._ (2020) M. Sánchez-Barquilla, R. E. F. Silva, and J. Feist,  _Cumulant expansion for the treatment of light–matter interactions in arbitrary material structures_ , J. Chem. Phys. 152, 034108 (2020).
  * Gardiner (2009) C. Gardiner,  _Stochastic Methods: A Handbook for the Natural and Social Sciences_ , Springer Series in Synergetics (Springer, Berlin Heidelberg, 2009).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2002).
  * Scully and Zubairy (1997) M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge Univ. Press, Cambridge, 1997).
  * Eberly and Wódkiewicz (1977) J. H. Eberly and K. Wódkiewicz,  _The time-dependent physical spectrum of light $\ast$_, J. Opt. Soc. Am. 67, 1252 (1977).
  * 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).
  * Gao  _et al._ (2012) T. Gao, P. S. Eldridge, T. C. H. Liew, S. I. Tsintzos, G. Stavrinidis, G. Deligeorgis, Z. Hatzopoulos, and P. G. Savvidis,  _Polariton condensate transistor switch_ , Phys. Rev. B 85, 235102 (2012).
  * Amo  _et al._ (2010) A. Amo, T. C. H. Liew, C. Adrados, R. Houdré, E. Giacobino, A. V. Kavokin, and A. Bramati,  _Exciton–polariton spin switches_ , Nat. Photonics 4, 361 (2010).
  * 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,  _A room-temperature organic polariton transistor_ , Nat. Photonics 13, 378 (2019).
  * Ballarini  _et al._ (2013) D. Ballarini, M. D. Giorgi, E. Cancellieri, R. Houdré, E. Giacobino, R. Cingolani, A. Bramati, G. Gigli, and D. Sanvitto,  _All-optical polariton transistor_ , Nat. Commun. 4, 1 (2013).
  * Gell-Mann (1962) M. Gell-Mann,  _Symmetries of Baryons and Mesons_ , Phys. Rev. 125, 1067 (1962).
  * Stone and Goldbart (2009) M. Stone and P. Goldbart,  _Mathematics for physics: a guided tour for graduate students_ (Cambridge University Press, Cambridge, 2009).