# Collective Effects of Organic Molecules based on Holstein-Tavis-Cummings
Model

Quansheng Zhang  Beijing Computational Science Research Center, Beijing
100193, China<EMAIL_ADDRESS>Ke Zhang  Helmholtz Institute Mainz, GSI
Helmholtzzentrum f$\ddot{u}$r Schwerionenforschung, 55099 Mainz, Germany
Johannes Gutenberg-University Mainz, 55099 Mainz, Germany

(August 27, 2024)

###### Abstract

We study the collective effects of an ensemble of organic molecules confined
in an optical cavity based on Holstein-Tavis-Cummings model. By using the
quantum Langevin approach and adiabatically eliminating the degree of freedom
of the vibrational motion, we analytically obtain the expression of the cavity
transmission spectrum to analyze the features of polaritonic states. As an
application, we show that the dependence for the frequency shift of the lower
polaritonic state on the number of molecules can be used in the detection of
the ultra-cold molecules. We also numerically analyze the fluorescence
spectrum. The variation of the spectral profile with various numbers of
molecules gives signatures for the modification of molecular conformation.

##  I Introduction

The interaction between molecule and light has played an important role in
photochemical reactivity [1], molecular spectroscopy [2, 3], chemical
fingerprinting [4], and the generation of non-classical states of light [5].
When molecules located in an optical cavity, the coupling strength can be
enhanced [6, 7, 8, 9]. In the strong coupling regime, the energy may oscillate
between molecules and cavity field at a faster rate than their relaxation
rates, inducing that the molecular and photonic states are hybridized into
polaritonic states [10]. In this situation, the potential energy landscape for
the molecular conformation will be modified, giving rise to a range of further
effects on the material properties of the molecules. As the modifications for
molecular conformation are tunable, it has opened up new routes for the
control of chemical reactivity [11, 12], energy and charge transport [13, 14],
and F$\ddot{\text{{o}}}$rster resonance energy transfer [15, 16, 17]. Besides,
the modifications for molecular conformation will bring influence on the
profile of molecular absorption and fluorescence spectra [18]. For the
electronic transitions, the coupling between electronic state and the
molecular conformation will lead to the structure of the emission (or
absorption) spectra observed in experiments [19, 9] with multiple peaks. The
observable molecular spectra will give signatures for the variation of
molecular conformation.

To provide an understanding of such phenomena, the Holstein-Tavis-Cummings
(HTC) model has been widely adopted to describe the light-electronic-vibration
problem [18, 20]. Furthermore, an approach based on the quantum Langevin
equation to solve the HTC model has been developed recently [17]. It provides
an alternative path to understand the Stokes and anti-Stokes processes [17],
Purcell effect under the influence of the phononic environments [21], and the
Floquent engineering [22] at the level of operators rather than states.
However, the previous works [17, 21, 22] mainly focus on a single molecule.
When multiple molecules confined in a cavity, the system will experience more
complex dynamics. In the case that many molecules are introduced to couple a
quantized cavity field, an efficient interaction between molecules will be
induced, leading to collective effects [23]. On the other hand, the
description of the electronic transition dressed by the vibrational motion
suggests that the vibrators for one molecule will be coupled to that for other
molecules. Considering above facts, a profound influence will be brought on
the polaritonic states, as well as the fluorescence (or absorption) spectrum.

Inspired of this, we provide here a further study of HTC model when multiple
molecules are presented. Going from a set of coupled standard Langevin
equations for the whole system, and adiabatically eliminating the degree of
freedom of the vibrational motion, a set of effective equations can be
derived, involving the average values of Pauli operators, photonic operators
and their correlations. From these equations, we study the collective effects
in the optical cavity via the cavity transmission spectrum, the steady
population and the fluorescence spectrum. To test the reasonableness of
approximation, we perform a numerical simulation for the case of two molecules
in the limit of weak coupling between the electronic states and molecular
conformation with help of the QuTiP package [24], and also take perturbative
treatment in the first order to compare against the perturbative method in
second order for general situations. We show that the profiles of cavity
transmission spectrum and the fluorescence spectrum can be modified via the
number of molecules, which will give signatures of the modification of
molecular conformation. Such phenomena will provide potential applications in
the detection of untracold molecules [25], as well as the control of chemical
reaction [12, 26], energy and charge transport [14, 20], and so on. In
contrast, we also give the analysis of a special case in the limit of without
the coupling between electronic states and molecular conformation, i.e. Tavis-
Cummings (TC) model.

The structure of this paper is organized as following. In Sec. II, we describe
the model and give the formal solution about the average value of the Pauli
lowering operator with the adiabatic elimination for the degree of freedom
vibrational motion. In Sec. III, we investigate the influence of the number of
molecules on the polaritonic features via the cavity transmitted field for HTC
model and TC model, respectively. In Sec. IV, we analyze the influence of the
number of molecules on the fluorescence spectrum for HTC model and TC model,
respectively. Finally, a summary is given in Sec. V.

##  II Model

Figure 1: (Color online) Model schematics. (a) An ensemble of $N$ organic
molecules are placed inside a Fairy-Perot microcavity with a single mode. The
dissipation for the cavity field occurs at rates $\kappa_{1}$ through the
mirror M1 and $\kappa_{2}$ through the mirror M2. (b) Diagram of molecular
level structure and relaxation processes. Each molecule is represented by two
harmonic oscillator potential surfaces of two electronic states with the shift
$\lambda$ in reaction coordinate. The cavity field resonantly couples to the
electronic transition $\ket{\textsl{g}}\leftrightarrow\ket{e}$ with the
coupling strength $g$. The relaxation processes of each molecule include the
electronic decay at rate $\Gamma$, and the pure dephasing at rate
$\Gamma_{\phi}$ and internal vibrational relaxation at rate $\gamma$.

As illustrated in Fig. 1, the proposed system of interest consists of an
ensemble of $N$ organic molecules, which confined in a Fabry-Perot
microcavity. Considering the organic molecules comprised of $z$ atoms, there
will be $3z-6$ (or $3z-5$ for linear molecules) normal modes of vibration in
the configuration space of the nuclear coordinates [20, 11]. When describing a
chemical reaction, often a few of these modes, which are associated with the
reaction coordinates, are investigated [27, 20, 11]. For a single reaction
coordinate, one can assume that the vibrational frequencies for the ground and
excited electronic states are the same. Thus, the system can be represented as
a two-level system (corresponding to its ground $\ket{\textsl{g}_{m}}$ and
excited $\ket{e_{m}}$ electric states) and a phonon mode (denoting the
harmonic vibrational degree of freedom with the annihilation operator
$b_{m}$).  Here, each molecule is labeled by the index $m$ (running from $1$
to $N$). The dynamic for a single free molecule is then governed by the
Holstein Hamiltonian [20]

$H_{m}=\nu
b_{m}^{\dagger}b_{m}+[\omega_{e}+\lambda\nu(b_{m}+b_{m}^{\dagger})]\sigma_{m}^{\dagger}\sigma_{m},$
(1)

where $\omega_{e}=\omega_{00}+\lambda^{2}\nu$ is the effective transition
frequency, with $\omega_{00}$ being the bare splitting between the ground and
excited electric states, and $\nu$ is the vibrational frequency for the
harmonic oscillator. The parameter $\lambda$ is related with the Huang-Rhys
factor $S$, i.e. $S=\lambda^{2}$. Here $S$ ranges from $0$ to $\sim 2$ [28,
29, 30], resulting from the displacement for the equilibrium position of the
vibrational mode between the excited and ground electronic states.
Furthermore, $\lambda$ also characters the coupling between the electronic
states and molecular conformation. In the limit of $\lambda=0$, one can
simplify the system for the molecule described by Eq. (1) will to be a two-
level system.

When the cavity field is single-mode, the molecule-cavity interaction can be
depicted by the TC Hamiltonian (we set $\hbar=1$)

$H_{\mathrm{int}}=\omega_{c}a^{\dagger}a+g\sum_{m=1}^{N}(a^{\dagger}\sigma_{m}+\mathrm{H.c.}),$
(2)

where the operator $a$ annihilates a cavity photon at frequency $\omega_{c}$,
$\sigma_{m}=\ket{\textsl{g}_{m}}\bra{e_{m}}$ annihilates the molecular
excitation and $g$ is the coupling strength for a single molecule. Notably,
this Hamiltonian implies the assumption that all molecules are equally coupled
to the cavity field via neglecting the influences resulting from the disorder
of transition dipole moment of molecules and molecular spatial distribution.
Such approximation will conveniently provide a qualitative description of the
collective behavior for organic polaritons [31].

Taking account of the free part of the molecules and the situation that the
cavity is probed by a classical field at frequency $\omega_{l}$ with the
coupling strength $\eta$, the total Hamiltonian for such system in the
rotating wave approximation can be given

$\displaystyle H=$
$\displaystyle\sum_{m}^{N}H_{m}+H_{\mathrm{int}}+i\eta(ae^{i\omega_{l}t}-a^{\dagger}e^{-i\omega_{l}t}).$
(3)

Applying a polaron transformation $H_{p}=U^{\dagger}HU$ with

$U=\exp[\sum_{m}^{N}\lambda\sigma_{m}^{\dagger}\sigma_{m}(b_{m}-b_{m}^{\dagger})],$
(4)

one can recast the Hamiltonian (3) into

$\displaystyle H_{p}=$ $\displaystyle\sum_{m}^{N}(\nu
b_{m}^{\dagger}b_{m}+\omega_{00}\sigma_{m}^{\dagger}\sigma_{m})+\omega_{c}a^{\dagger}a$
(5)
$\displaystyle+(\sum_{m}^{N}g\mathcal{D}_{m}^{\dagger}\sigma_{m}a^{\dagger}+i\eta
e^{i\omega_{l}t}a+\mathrm{H.c.}),$

where $\mathcal{D}_{m}=\exp[\lambda(b_{m}-b_{m}^{\dagger})]$ is a displaced
operator for $m$-th molecule. In this representation, the coupling between the
two-level system and cavity field is dressed by vibrations via the Frank-
Condon factors ${}_{\mathrm{vib}}\bra{n}\mathcal{D}\ket{m}_{\mathrm{vib}}$ for
the transition
$\ket{\textsl{g}}\ket{n}_{\mathrm{vib}}\rightarrow\ket{e}\ket{m}_{\mathrm{vib}}$,
where $\ket{n}_{\mathrm{vib}}$ (or $\ket{m}_{\mathrm{vib}}$) are Fock states
for vibrations.

To model the open system, we should consider the relaxation processes as
depicted in Fig. 1(b), including the loss of cavity field at rate
$\kappa=\kappa_{1}+\kappa_{2}$ (encompassing losses via both mirrors), the
spontaneous emission of the electronic transition at rate $\Gamma$, and the
pure dephasing of the electronic transition at rate $\Gamma_{\phi}$. Their
effects are described as Lindblad terms, and represented by
$\mathcal{L}_{\gamma_{\mathcal{O}}}[\mathcal{O}]=\gamma_{\mathcal{O}}(2\mathcal{O}\rho\mathcal{O}^{\dagger}-\mathcal{O}^{\dagger}\mathcal{O}\rho-\rho\mathcal{O}^{\dagger}\mathcal{O})$
with a collapse operator $\mathcal{O}$ and a corresponding decay rate
$\gamma_{\mathcal{O}}$. Finally the dynamics of the system can be described by
a master equation of the density matrix $\rho$

$\displaystyle\dot{\rho}=$
$\displaystyle-i[H,\rho]+\mathcal{L}_{\kappa}[a]+\sum_{m}^{N}\mathcal{L}_{\Gamma_{\phi}}[\sigma_{m}^{\dagger}\sigma_{m}-\sigma_{m}\sigma_{m}^{\dagger}]$
(6)
$\displaystyle+\sum_{m}^{N}\mathcal{L}_{\Gamma}[\sigma_{m}]+\sum_{m}^{N}\mathcal{L}_{\gamma}[b_{m}].$

The last term in Eq. (6) represents the vibrational relaxation process at rate
$\gamma$. Physically, when considering the coupling between the vibrator and a
reservoir of phonons (also described as Brownian noise dissipation model [17,
21]) in a solvent or surrounding medium, the vibrational relaxation of
molecules cannot be simply expressed in Lindblad form [22, 32, 21]. In this
work, for the sake of convenience to discuss, we still adopt the Lindblad
decay model to describe the vibrational relaxation process. This approach has
been used in several theoretical studies to explore the molecular spectroscopy
[22, 17], Purcell effect [21], energy transfer [17] and organic polariton
lasing [33]. In the limit $\lambda^{2}\gamma\ll\Gamma$, the Lindblad decay
model for the vibrational relaxation will be indistinguishable to the Brownian
noise model [22, 17]. Additionally, the spontaneous emission is also involved
in the vibrational motion, which can be ignorable due to its decay rate much
smaller than the vibrational relaxation rate $\gamma$.

The master equation method is usually equivalent to the Langevin approach. For
the given master equation in the Lindblad form, described as
$\mathcal{L}_{\gamma_{\mathcal{O}}}[\mathcal{O}]$ with a collapse operator
$\mathcal{O}$ and decay rate $\gamma_{\mathcal{O}}$, we can map it onto a
Langevin form [34] as

$\displaystyle\frac{d}{dt}\mathcal{A}=$ $\displaystyle
i[H,\mathcal{A}]-\sum_{j}\left[\mathcal{A},\mathcal{O}_{j}^{\dagger}\right]\left\\{\gamma_{\mathcal{O}_{j}}\mathcal{O}_{j}-\sqrt{2\gamma_{\mathcal{O}_{j}}}\mathcal{O}_{j}^{\mathrm{in}}\right\\}$
(7)
$\displaystyle+\sum_{j}\left\\{\gamma_{\mathcal{O}_{j}}\mathcal{O}_{j}^{\dagger}-\sqrt{2\gamma_{j}}\mathcal{O}_{j}^{\mathrm{in}\dagger}\right\\}[\mathcal{A},\mathcal{O}_{j}],$

with an arbitrary system operator $\mathcal{A}$ and any collapse operator
$\mathcal{O}_{j}$ in the set {$b_{m}$, $\sigma_{m}$,
$\sigma_{m}^{\dagger}\sigma_{m}-\sigma_{m}\sigma_{m}^{\dagger}$, $a$}.

Under the condition of weak driving $\eta\ll\kappa$, the cavity field in
steady situation will have much less one photon and the molecules will pretty
much stay in the electronic ground state, i.e.
$\sigma_{m}^{\dagger}\sigma_{m}\ll 1$. In the rotating frame at the probe
frequency $\omega_{l}$, the quantum Langevin equations for the cavity field
$a$, and polarized operator
$\tilde{\sigma}_{m}=\mathcal{D}_{m}^{\dagger}\sigma_{m}$ are given as [17]

$\displaystyle\frac{d}{dt}a=$
$\displaystyle-(i\Delta_{c}+\kappa)a-ig\sum_{m}^{N}\sigma_{m}+\sqrt{2\kappa_{1}}A_{\mathrm{in}}$
(8a) $\displaystyle+\sqrt{2\kappa_{2}}a_{2}^{\mathrm{in}}(t),$
$\displaystyle\frac{d}{dt}\tilde{\sigma}_{m}\approx$
$\displaystyle-(i\Delta+\Gamma_{\bot})\tilde{\sigma}_{m}-ig\mathcal{D}_{m}^{\dagger}a+\sqrt{2\Gamma_{\bot}}\mathcal{D}_{m}^{\dagger}\sigma_{m}^{\mathrm{in}},$

with the detuning $\Delta_{c}=\omega_{c}-\omega_{l}$ and
$\Delta=\omega_{00}-\omega_{l}$. Here, $\Gamma_{\bot}=\Gamma+2\Gamma_{\phi}$
is the effective transverse relaxation rate. Notably, under the description of
the vibrational relaxation in the Lindblad form, an additional dephasing will
be caused to the polarized operator. However, we can neglect such dephasing as
this will lead to the disagreement with experimental observations of lifetime
for the electronic transition [17, 22, 21]. Additionally,
$A_{\mathrm{in}}(t)=\eta/\sqrt{2\kappa_{1}}+a_{1}^{\mathrm{in}}(t)$ denotes
the effective cavity input with zero-average input noise
$a_{1}^{\mathrm{in}}(t)$. Here we assume that the temperature for cavity field
is zero, i.e. $T_{\mathrm{cav}}=0$, then the non-vanishing correlation
$\braket{a_{1}^{\mathrm{in}}(t)a_{1}^{\mathrm{in}\dagger}(t^{\prime})}=\delta(t-t^{\prime})$
and
$\braket{a_{2}^{\mathrm{in}}(t)a_{2}^{\mathrm{in}\dagger}(t^{\prime})}=\delta(t-t^{\prime})$
will be obtained. The electronic transition is also affected by a white nose
input $\sigma_{m}^{\mathrm{in}}$ with nonzero correlation
$\braket{\sigma_{m}^{\mathrm{in}}(t)\sigma_{n}^{\mathrm{in}\dagger}(t^{\prime})}=\delta_{mn}\delta(t-t^{\prime})$.

Then we can formally integrate the Eq. (8) to get the solution for
$\tilde{\sigma}_{m}$ and subsequently for Pauli lowering operator

$\displaystyle\sigma_{m}(t)=$
$\displaystyle-\int_{0}^{t}dt_{1}e^{-(i\Delta+\Gamma_{\bot})(t-t_{1})}\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})[iga(t_{1})-\sqrt{2\Gamma_{\bot}}\sigma_{m}^{in}(t_{1})]+e^{-(i\Delta+\Gamma_{\bot})t}\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(0)\sigma_{m}(0).$
(9)

In a solvent or surrounding medium, the coupling between the vibrator and a
reservoir of phonons will lead to large relaxation rate for vibrational motion
(i.e., $\gamma\gg\Gamma_{\bot}$) [21, 35, 32]. Additionally, we also assume
that the vibrational relaxation rate much larger than cavity field (i.e.,
$\gamma\gg\kappa$). Under these conditions, the timescale of vibrational
relaxation will be much shorter than that for the relaxation of electronic
state, and also the decay of cavity field. Therefore, the vibrational state
may be considered to be in the steady state [32]. Via treating the vibrations
as a Markovian phonon bath [36, 32], the degrees of freedom for cavity field
and vibrational motion can be separated as

$\braket{a(t_{1})\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}\approx\braket{a(t_{1})}\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})},$
(10)

as well as that for Pauli lowering operators and vibrational motion

$\braket{\sigma_{n}(t_{1})\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}\approx\braket{\sigma_{n}(t_{1})}\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}.$
(11)

If these molecules are initially prepared in the ground electronic state, and
zero photon populated in cavity, the last term in Eq. (9) will be canceled
when we take average over the cavity field, the vibrational motion, and the
electronic state. Finally, the average dipole moment $\braket{\sigma_{m}}$ can
be derived

$\displaystyle\braket{\sigma_{m}}=-ig\int_{0}^{\infty}dt_{1}\mathcal{\mathcal{F}}_{m}(t-t_{1})\braket{a(t_{1})},$
(12)

with the definition
$\mathcal{\mathcal{F}}_{m}(t-t_{1})=\Theta(t-t_{1})\exp[-(i\Delta+\Gamma_{\bot})(t-t_{1})]\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}$,
where $\Theta(t)$ is the Heaviside step function. The two-time correlation for
the displacement operators
$\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}$ can be expressed
as (at $t>t_{1}$)

$\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}=e^{-\lambda^{2}}e^{\lambda^{2}e^{-(\gamma+i\nu)(t-t_{1})}},$
(13)

when the effective vibrational temperature satisfies
$k_{B}T_{\mathrm{vib}}/\hbar\nu\ll 1$ (see Appendix A). The important quantity
to be emphasized is the function
$\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}$, which implies
the influence of the vibrational mode on the excitons in organic molecules.
When $\lambda=0$, one can achieve the value of the function
$\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}=1$. Then the
interaction between electronic states and molecular conformation will be
removed. For finite $\lambda$, the coupling between the vibrational mode and
the excitons will result in some new fascinating physics.

Notably, the Pauli operator can feed back into itself via the cavity field,
resulting in higher order correlation between the displaced operators. The
adiabatic approximation depicted by Eq. (13) has neglected the contributions
from these higher order correlations. To explore the influence of these higher
order correlations, we have taken further discussion in Appendix B and C.

##  III Cavity Transmission

Figure 2: (Color online) Cavity-molecule spectroscopy. (a-f) Influence of the
number $N$ (increases as $2$, $20$, and $200$) of identical organic molecules
on the intensity of transmitted field at resonance $\omega_{c}=\omega_{00}$ at
$\lambda=0$ (blue dashed line), $\lambda=0.2$ (red solid line) and
$\lambda=0.6$ (black dotted line), respectively. (a) The green circle marks
shows the numerical simulation results for the case $\lambda=0.2$, $N=2$ with
the help of QuTiP. Intensity of transmitted field at resonance
$\omega_{c}=\omega_{00}$ as a function of the detuning $\Delta_{c}/\nu$ and
the number $N$ of organic molecules at $\lambda=0.2$ (g) and $\lambda=0.6$
(h), respectively. The spatial intensity profile of the transmitted field is
shown in arbitrary units but to scale. The green and blue dashed lines show
the frequency shift of the lower and upper polaritonic states various the
number $N$ of organic molecules, respectively, given by Eq. (19b). The
frequency shift of the lower polaritonic state is about $\omega_{-}\approx
g\sqrt{N}$ (yellow dotted line). The other parameters are $\nu/\Gamma=250$,
$\Gamma_{\phi}/\Gamma=2\kappa_{1}/\Gamma=2\kappa_{2}/\Gamma=1$,
$\gamma/\Gamma=50$, $g/\Gamma=5$, and $\eta/\Gamma=0.01$.

For the molecule-cavity system discussed in the previous section, the
electronic transition will be dressed by vibrational motion with the coupling
between the electronic states and molecular conformation. In the limit of
large vibrational relaxation (i.e. $\gamma\gg\kappa$ and
$\gamma\gg\Gamma_{\bot}$), the vibrational mode can be treated as a locale
phonon reservoir under the Markov approximation. Going from the Langevin
equation of polarized operator and adiabatically eliminating the degree of
freedom of vibrational motion, we formally derive the expression of the
average Pauli lowering operator $\braket{\sigma_{m}}$. It is noteworthy that
the first term in Eq. (12) represents a convolution. This character will
inspire us to utilize the Laplace transformation (defined as
$\braket{\bar{f}}(s)=\int_{0}^{\infty}dt\>f(t)e^{-st}$ for time dependent
function $f(t)$ at $t\geq 0$) to simplify the derivation. Then, Eq. (12) can
be reformed as

$\braket{\bar{\sigma}_{m}}=-ig\bar{\mathcal{F}}_{m}(s)\braket{\bar{a}},$ (14)

where $\bar{\mathcal{F}}_{m}$ is the Laplace transform of
$\mathcal{\mathcal{F}}_{m}(t)$, expressed as

$\bar{\mathcal{F}}_{m}=\sum_{k}\frac{\lambda^{2k}}{k!}\frac{e^{-\lambda^{2}}}{s+i(\Delta+k\nu)+(\Gamma_{\bot}+k\gamma)},$
(15)

with the Laplace transform variable $s$. Similarly, taking average over the
cavity field as well as the electronic degrees of freedom on both sides of Eq.
(8a), and applying the Laplace transformation, we finally obtain

$s\braket{\bar{a}}=-(i\Delta_{c}+\kappa)\braket{\bar{a}}-ig\sum_{m}^{N}\braket{\bar{\sigma}_{m}}+\frac{\eta}{s}.$
(16)

Considering the system with $N$ identical molecules, the Laplace form for the
average Pauli lowering operator $\braket{\bar{\sigma}_{m}}$ and the function
$\bar{\mathcal{F}}_{m}(s)$ would be the same for any molecule ($m$). Combining
the Laplace forms of the Pauli operator described by Eq. (14) and the cavity
field illustrated by Eq. (16), we can get

$\displaystyle\braket{\bar{\sigma}_{m}}$
$\displaystyle=-\frac{ig\eta\bar{\mathcal{F}}_{m}}{s[Ng^{2}\bar{\mathcal{F}}_{m}+i\Delta_{c}+\kappa+s]},$
(17a) $\displaystyle\braket{\bar{a}}$
$\displaystyle=\frac{\eta}{s[i\Delta_{c}+s+\kappa+Ng^{2}\bar{\mathcal{F}}_{m}]}.$
(17b)

From the final value theorem [37], we get the steady values

$\displaystyle\braket{\sigma_{m}}_{\mathrm{ss}}$
$\displaystyle=\lim_{s\rightarrow
0}s\braket{\bar{\sigma}_{m}}=-\frac{ig\eta\chi}{i\Delta_{c}+\kappa+Ng^{2}\chi},$
(18a) $\displaystyle\braket{a}_{\mathrm{ss}}$
$\displaystyle=\lim_{s\rightarrow
0}s\braket{\bar{a}}=\frac{\eta}{(i\Delta_{c}+\kappa)+Ng^{2}\chi},$ (18b)

where “ss” stands for the steady state situation, and
$\chi=\underset{s\rightarrow 0}{\lim}\bar{\mathcal{F}}_{m}$. The item
$Ng^{2}\chi$ in Eq. (18b) reflects the hybridization between photonic states
and molecule electronic states. One can approximately obtain the hybridized
decay rates and frequencies (under resonance conditions, i.e., fixed
$\Delta=\Delta_{c}=0$) [38, 10]

$\displaystyle\Gamma_{\pm}$
$\displaystyle=\frac{\Gamma_{\mathrm{eff}}+\kappa}{2}\pm\mathfrak{R}\sqrt{\frac{(\Gamma_{\mathrm{eff}}+i\Delta_{\mathrm{eff}}-\kappa)^{2}}{4}-Ng^{2}},$
(19a) $\displaystyle\omega_{\pm}$
$\displaystyle=-\frac{\Delta{}_{\mathrm{eff}}}{2}\pm\mathfrak{I}\sqrt{\frac{(\Gamma_{\mathrm{eff}}+i\Delta_{\mathrm{eff}}-\kappa)^{2}}{4}-Ng^{2}},$
(19b)

where $\Gamma_{\mathrm{eff}}=\mathfrak{R}{\lim}1/\chi$ and
$\Delta_{\mathrm{eff}}=\mathfrak{I}{\lim}1/\chi$ denote the effective decay
rate and additional frequency shift for excitons, respectively.

In experiments, the hybridized decay rates and frequencies can be detected via
the cavity transmitted field. According to the input and output relation
$a_{2,\mathrm{out}}(t)=\sqrt{2\kappa_{2}}a(t)$, the transmission $\mathcal{T}$
via
$\mathcal{T}=\braket{a_{2,\mathrm{out}}}_{\mathrm{ss}}/\braket{A_{\mathrm{in}}}_{\mathrm{ss}}$
in steady state reign will be expressed as

$\mathcal{T}=\frac{2\sqrt{k_{1}k_{2}}}{(i\Delta_{c}+\kappa)+Ng^{2}\chi}.$ (20)

Figures 2(a-f) plot the intensity of the cavity transmission at resonance
$\omega_{c}=\omega_{00}$ in the various number $N$ of the organic molecules
for TC (see blue dashed line) and HTC (see red solid line and black dotted
line) model, respectively. The parameters we chosen are given in the caption
of Fig. 2, which associate with the recent theoretical researches [17, 21, 33,
22]. When multiple molecules emerged in a cavity, the system will be in the
collective behavior with a bright mode that couples to the cavity field with
an enhanced strength $g\sqrt{N}$ (for identical molecules). In the limit of
$\lambda=0$, [see blue dashed line in Fig. 2(a-c)], the cavity transmission
will be split into two symmetric peaks with the central frequency separation
$2\sqrt{Ng^{2}-(\Gamma_{\bot}+\kappa)^{2}/4}$ [39, 40, 38].

For finite $\lambda$, the coupling between electronic states and molecular
conformation will produce additional dephasing as well as the frequency shift
for excitons, resulting in an asymmetric cavity transmission profile [see red
solid line in Fig. 2(c) and black dotted line in Figs. 2(d-f)]. In the limit
of $\lambda^{2}\ll 1$ and $g\sqrt{N}\ll\nu$, the evolution of the transmitted
field for HTC model with the detuning $\Delta_{c}$ will perform similar
behaviors to the TC model, due to the small Franck-Condon overlap and the weak
collective effects. To clarify the influence of the number $N$ of organic
molecules on the hybridized frequencies as well as the decay rates, we
numerically show the spatial intensity profile of the transmitted field. It is
clearly to see that the evolution of the frequency for lower polariton with
the number $N$ of molecules obeys the rule given by Eq. (19b) [see green
dashed line in Fig. 2(g, h)]. For blue-detuned illumination, the spatial
intensity profile of the transmitted field seems more complicated [see Fig.
2(c-h)]. Typically, when the collective coupling strength $g\sqrt{N}>\nu$ for
large Franck-Condon overlap (i.e. $\lambda^{2}\sim 1$), one more polaritonic
peak appears in the cavity transmission profile [see black dotted line Fig.
2(f)]. For the left branch in Fig. 2(h), the relation between the number $N$
of molecules and the central frequency still obeys the rule given in Eq. (19b)
[see blue dashed line in Fig. (2)(h)]. This branch corresponds to the upper
polaritons. The appearance of the another branch is caused by the perturbation
for the excitation state of the vibrational motion, corresponding to the dark
polaritonic state [10].

In addition, we would like to point out here the cavity transmitted field can
be used to detect the ultra-cold molecules [25]. Considering the presence of
rovibrational states, several theoretical studies have modeled one molecule as
a several-level system [41, 25]. One can estimate the number of molecules via
the frequency separation between upper and lower polaritonic states. In
contrast to the previous work [25], the coupling between electronic states and
molecular conformation described by HTC model will make the cavity transmitted
field become more complicated [see Fig. 2(g-h)]. Particularly for the large
Franck-Conda overlap, it will be a challenge task to detect molecules
according to the method introduced by Ref. [25] due to the appearance of the
peak for dark polaritonic state, and the suppressed intensity for the upper
polaritonic state. Notably, a single peak appeared in the spatial intensity
profile of the transmitted field when $\Delta_{c}>0$, which is associated with
the lower polaritons. Considering the relation between the number $N$ of
molecules and the central frequency of the lower polaritonic peaks described
by Eq. (19b), we can approximately get the number of molecules as
$N\approx\omega_{-}^{2}/g^{2}$ [see yellow dotted line in Fig. 2(g, h)].

##  IV Fluorescence Spectrum

In the previous section, we have analyzed the polaritonic features via the
transmission field profile. With the increasing of the number of molecules,
qualitatively different phenomena arise while considering the collective
behavior for the emitters. On the other hand, the molecules can interact with
each other via the coupling with the quantized cavity field. This will give
rise to the modification of the molecular configuration [42], as well as the
transition between the electronic ground and excited states. In this section,
we will illustrate how the feature is reflected via the fluorescence spectrum.

Before the introduction of fluorescence spectrum, let us first give the
definition of the fluctuation operator. Since the expect values for Pauli
lowering operator $\braket{\sigma_{m}}_{\mathrm{ss}}$ and the cavity field
$\braket{a}_{\mathrm{ss}}$ have been obtained in steady state given by Eq.
(18), one can represent the fluctuation about this steady state by the
following operators

$\begin{aligned}
\delta\sigma_{m}&=\sigma_{m}-\braket{\sigma_{m}}_{\mathrm{ss}}\\\ \delta
a&=a-\braket{a}_{\mathrm{ss}}.\end{aligned},$ (21)

Considering the formal solution for Pauli operator given by Eq. (9), one can
easily derive the expression for the autocorrelation function

$\displaystyle\braket{\delta\sigma_{m}^{\dagger}(t)\delta\sigma_{m}(t+\tau)}=$
$\displaystyle-ig\int_{t}^{t+\tau}dt_{1}e^{-(i\Delta+\Gamma_{\bot})(t+\tau-
t_{1})}\braket{\delta\sigma_{m}^{\dagger}(t)\mathcal{D}_{m}(t+\tau)\mathcal{D}_{m}^{\dagger}(t_{1})\braket{a}_{\mathrm{ss}}}-\braket{\sigma_{m}}_{\mathrm{ss}}\braket{\delta\sigma_{m}^{\dagger}(t)}$
$\displaystyle-ig\int_{t}^{t+\tau}dt_{1}e^{-(i\Delta+\Gamma_{\bot})(t+\tau-
t_{1})}\braket{\delta\sigma_{m}^{\dagger}(t)\mathcal{D}_{m}(t+\tau)\mathcal{D}_{m}^{\dagger}(t_{1})\delta
a(t_{1})}$ (22)
$\displaystyle+e^{-(i\Delta+\Gamma_{\bot})\tau}\braket{\delta\sigma_{m}^{\dagger}(t)\mathcal{D}_{m}(t+\tau)\mathcal{D}_{m}^{\dagger}(t)\sigma_{m}(t)}.$

Under the large vibrational relaxation condition $\gamma\gg\kappa$ and
$\gamma\gg\Gamma$, the correlation time for vibrational motion would be much
shorter than the timescale of the correlation between cavity field and
molecules, as well as the intra-molecule correlations. Therefore, the four-
operator correlation functions can be separated as

$\displaystyle\braket{\delta\sigma_{m}^{\dagger}(t)\mathcal{D}_{m}(t+\tau)\mathcal{D}_{m}^{\dagger}(t_{1})\mathcal{O}(t_{1})}$
$\displaystyle\approx\braket{\delta\sigma_{m}^{\dagger}(t)\mathcal{O}(t_{1})}\braket{\mathcal{D}_{m}(t+\tau)\mathcal{D}_{m}^{\dagger}(t_{1})}$
(23)

with $\mathcal{O}(t)$ in the set $\\{$$a(t)$, $\sigma_{m}(t)$,
$\sigma_{n}(t)$$\\}$. In the limit of a long time scale, the system would be
in the steady state. Then the first line in Eq. (22) can be neglected as the
expectation value for fluctuation operator
$\braket{\delta\sigma_{m}}_{\mathrm{ss}}$ will be zero. The formulation for
the desired auto-correlation function is deduced with convolution integral

$\braket{\delta\sigma_{m}^{\dagger}(0)\delta\sigma_{m}(\tau)}_{\mathrm{ss}}=-ig\int_{0}^{\infty}dt_{1}\mathcal{F}_{m}(\tau-
t_{1})\braket{\delta\sigma_{m}^{\dagger}(0)\delta
a(t_{1})}_{\mathrm{ss}}-\braket{\delta\sigma_{m}^{\dagger}\delta\sigma_{m}}_{\mathrm{ss}}\mathcal{F}_{m}(\tau),$
(24)

where
$\braket{\delta\sigma_{m}^{\dagger}(0)\delta\mathcal{O}(\tau)}_{\mathrm{ss}}\equiv\underset{t\rightarrow\infty}{\lim}\braket{\delta\sigma_{m}^{\dagger}(t)\delta\mathcal{O}(t+\tau)}_{\mathrm{ss}}$
and
$\braket{\delta\sigma_{m}^{\dagger}\delta\sigma_{m}}_{\mathrm{ss}}=\braket{\sigma_{m}^{\dagger}\delta\sigma_{m}}_{\mathrm{ss}}-\braket{\delta\sigma_{m}}_{\mathrm{ss}}\braket{\sigma_{m}^{\dagger}}_{\mathrm{ss}}$.
Taking the Fourier transformation with the definition
$S_{\mathcal{O}}^{m}(\omega)=\int_{-\infty}^{\infty}d\tau\,\braket{\delta\sigma_{m}^{\dagger}(0)\delta\mathcal{O}(\tau)}_{\mathrm{ss}}e^{-i\omega\tau}$,
the correlation function given by Eq. (24) can be reformed in the Fourier
domain, i.e.,

$S_{\sigma_{m}}^{m}=-\mathcal{\tilde{F}}_{m}[igS_{a}^{m}-\braket{\delta\sigma_{m}^{\dagger}\delta\sigma_{m}}_{\mathrm{ss}}],$
(25)

where $\mathcal{\tilde{F}}_{m}(\omega)$ is the Fourier form for the factor
$\mathcal{F}_{m}(t)$.

Taking the similar method, one can respectively get the Fourier forms for the
correlation functions
$\braket{\delta\sigma_{m}^{\dagger}(0)\delta\sigma_{n}(\tau)}_{\mathrm{ss}}$
between molecules, and also the correlation function
$\braket{\delta\sigma_{m}^{\dagger}(0)\delta a(\tau)}_{\mathrm{ss}}$ between
cavity and molecules

$\displaystyle\begin{aligned}
S_{\sigma_{n}}^{m}=&-\mathcal{\tilde{F}}_{n}(igS_{a}^{m}-\braket{\delta\sigma_{m}^{\dagger}\delta\sigma_{n}}_{\mathrm{ss}}),\\\
S_{a}^{m}=&-ig\tilde{\mathcal{F}}_{a}(\sum_{n\neq
m}^{N}S_{\sigma_{n}}^{m}+S_{\sigma_{m}}^{m})+\tilde{\mathcal{F}}_{a}\braket{\delta\sigma_{m}^{\dagger}\delta
a}_{\mathrm{ss}},\end{aligned}$ (26)

where
$\braket{\delta\sigma_{m}^{\dagger}\delta\mathcal{O}}_{\mathrm{ss}}=\braket{\sigma_{m}^{\dagger}\mathcal{O}}_{\mathrm{ss}}-\braket{\mathcal{O}}_{\mathrm{ss}}\braket{\sigma_{m}^{\dagger}}_{\mathrm{ss}}$,
and $\tilde{\mathcal{F}}_{a}(\omega)=1/[i(\Delta_{c}+\omega)+\kappa]$.

For identical molecules, the functions $S_{\sigma_{n}}^{m}$ ($n\neq m$) are
same for any molecule pair ($m$, $n$) and $S_{\sigma_{m}}^{m}$ are
interchangeable for any molecule ($m$). As a result, Eqs. (25-26) can be
simply reformed as

$\mathcal{M}_{s}\bm{\mathcal{S}}_{\mathrm{Vec}}+\bm{\mathcal{S}}_{\mathrm{in}}=0,$
(27)

with the coefficient matrix

$\mathcal{M}_{s}=\left(\begin{array}[]{ccc}-1&0&-ig\mathcal{\tilde{F}}_{m}\\\
0&-1&-ig\mathcal{\tilde{F}}_{n}\\\
-ig\mathcal{\tilde{F}}_{a}&-i(N-1)g\mathcal{\tilde{F}}_{a}&-1\end{array}\right),$

the vector for the Fourier form of correlation functions
$\bm{\mathcal{S}}_{\mathrm{Vec}}=($$S_{\sigma_{m}}^{m},$ $S_{\sigma_{n}}^{m}$,
$S_{a}^{m})^{T}$, and the input vector
$\bm{\mathcal{S}}_{\mathrm{in}}=($$\mathcal{\tilde{F}}_{m}\braket{\delta\sigma_{m}^{\dagger}\delta\sigma_{m}}_{\mathrm{ss}}$,
$\mathcal{\tilde{F}}_{n}\braket{\delta\sigma_{m}^{\dagger}\delta\sigma_{n}}_{\mathrm{ss}}$,
$\mathcal{\tilde{F}}_{a}\braket{\delta\sigma_{m}^{\dagger}\delta
a}_{\mathrm{ss}})^{T}$. To this end, the solution can be obtained via
$\bm{\mathcal{S}}_{\mathrm{Vec}}=-\mathcal{M}_{s}^{-1}\bm{\mathcal{S}}_{\mathrm{in}}.$
It should be noted that the above equation can also be derived via considering
the quantum regression theorem. Different with the traditional methods [43],
the convolution integral for $\braket{\sigma_{m}(t)}$ given by Eq. (12)
requires us to get the dynamic equation for the correlation function from the
perspective of the frequency domain. In addition, to obtain the value of the
vector $\bm{\mathcal{S}}_{\mathrm{Vec}}$, one need to calculate the
expectations $\braket{\sigma_{m}^{\dagger}\mathcal{O}}$ in steady state, which
have been discussed in Appendix C.

According to the definition of the collective operator
$J_{-}=\sum_{m=0}^{N}\sigma_{m}$, the fluorescence spectrum in steady
situation can be achieved by taking the Fourier transformation
$S=\underset{t\rightarrow\infty}{\mathrm{lim}}\mathfrak{R}\int_{-\infty}^{\infty}d\tau\,\braket{\delta
J_{-}^{\dagger}(t)\delta J_{-}(t+\tau)}e^{-i\omega\tau}$. Since
$\braket{\delta\sigma_{m}^{\dagger}(t)\delta\sigma_{m}(t+\tau)}$,
$\braket{\delta\sigma_{m}^{\dagger}(t)\delta\sigma_{n}(t+\tau)}$ ($m\neq n$)
are same for all molecules and all molecular pairs respectively, we can
further get the relation

$S=N\mathfrak{R}S_{\sigma_{m}}^{m}+N(N-1)\mathfrak{R}S_{\sigma_{n}}^{m}.$ (28)

Figure 3: (Colored online) The fluorescence spectrum at resonance
$\omega_{c}=\omega_{00}=\omega_{l}$ with various number $N$ of organic
molecules (increases as $2$, $20$, and $200$). We compare the HTC model for
different number $N$ of molecules, against the TC model (blue dashed line)
when (a-c) $\lambda=0.2$ and (d-f) $\lambda=0.6$. (a) The green circle markers
are numerical calculated with the help of QuTiP for the case $\lambda=0.2$,
and $N=2$. The other parameters are the same with Fig. 2.

Figure 3 illustrates the fluorescence spectrum at resonance
$\omega_{c}=\omega_{00}=\omega_{l}$. With the number $N$ of molecules
increasing, the spectral weight shows some interesting features. In the limit
of $\lambda=0$, the profile for the fluorescence spectrum will exhibit two
symmetric peaks with the central frequencies $\omega_{\pm}$ given by Eq. (19b)
(see blue dashed line). For finite $\lambda$, the coupling between electronic
and molecular conformation will bring us more amazing phenomena, e.g. the
Stokes scatting, the asymmetric profile for fluorescence spectrum [see red
solid line in Fig. 3(a-c) and black dotted line in Fig. 3(d-f)]. When $N=2$,
one more peak at frequency $\omega=-\nu$ appeared compared with TC (i.e.
$\lambda=0$) model [see red solid line in Fig. 3(a) and black dotted line in
Fig. 3(d)]. This peak corresponds to the Stokes process
$\ket{e}\ket{0}_{\mathrm{vib}}\rightarrow\ket{\textsl{g}}\ket{1}_{\mathrm{vib}}$,
with the transition rate characterized by the Franck-Condon factors
${}_{\mathrm{vib}}\bra{1}\mathcal{D}(\lambda)\ket{0}_{\mathrm{vib}}$ [44]. As
the number $N$ of organic molecules increases, the coverage of the Stokes line
in the frequency domain will be broadened [see black dotted line in Fig. 3(e,
f)]. Especially, at $N=20$, one can find that another peak appeared at
frequency $\omega=-2\nu$, corresponding to the Stoke process
$\ket{e}\ket{0}_{\mathrm{vib}}\rightarrow\ket{\textsl{g}}\ket{2}_{\mathrm{vib}}$.
The dependence of the profile for the Stokes line on the number of molecules
reflects that the molecular conformation has been modified via the collective
behavior. When the collective coupling strength $g\sqrt{N}>\nu$, the
contribution of dark polartonic state should be taken into account [see black
dotted line in Fig. 3(f)], which would bring new physics.

##  V Conclusion

In summary, we have studied the collective behavior for organic molecules
based on the HTC model. Via adiabatically eliminating the degree of freedom of
the vibrational motion, we have derived a set of linear equations in Laplace
(or Fourier) domain, involving the Pauli operators, the correlation between
cavity field and molecules, and the correlation between molecules. Based on
these equations, we have studied the cavity transmission spectrum for an
ensemble of identical organic molecules. Our results reveal that, different
with TC model, the coupling between electronic states and molecular
conformation will produce additional dephasing as well as the frequency shift
for excitons. This will result in some amazing phenomena, e.g. the suppressed
upper polaritonic peak, dark polaritonic peak with central frequency slightly
changed via increasing the number of molecules. In addition, the cavity
transmitted field is suggested to detect the ultra-cold molecules [25].
Considering the relationship between the hybridized frequencies and the number
of molecules, one can estimate the number of molecules via the frequency shift
for the lower polaritonic state. Moreover, we also study the steady population
on the excited state for total molecules and the fluorescence spectrum. With
the number of molecules increasing, the coverage of the Stokes line in the
frequency domain is broadened. This means that the probabilities for
transition from electronic excited state to the vibrational excited states in
the electronic ground state increased. In other words, the conformation of
molecules will be modified due to the collective effects, providing a
potential application in the control of the chemical reactivity [11, 12],
energy and charge transport [13, 14], and F$\ddot{\text{{o}}}$rster resonance
energy transfer [15, 16, 17].

###### Acknowledgements.

This work was supported by the Natural Science Foundation of China (under
Grants No. 12074030 and No. U1930402). Q. Z. acknowledges fruitful discussions
with Lei Du, Yong Li and Michael Reitz.

##  Appendix A Dynamic of the Vibrational Motion

In the limit of weak probe fields $\eta\ll\kappa$, the Heisenberg-Langevin
equation for the vibrational motion is given [17] by

$\frac{d}{dt}b=-i(\nu-i\gamma)b+\sqrt{2\gamma}b_{\mathrm{in}},$ (29)

where $b_{\mathrm{in}}$ is the annihilation operator of thermal noise with
zero-average value. The non-vanishing correlation functions are

$\displaystyle\begin{aligned}
\braket{b_{\mathrm{in}}^{\dagger}(t)b_{\mathrm{in}}(t_{1})}&=\bar{n}_{\mathrm{th}}\delta(t-t_{1}),\\\
\braket{b_{\mathrm{in}}(t)b_{\mathrm{in}}^{\dagger}(t_{1})}&=(1+\bar{n}_{\mathrm{th}})\delta(t-t_{1}),\end{aligned}$
(30)

where $\bar{n}_{\mathrm{th}}=1/[\exp(\hbar\nu/k_{B}T_{\mathrm{vib}})]-1]$ is
the average thermal photon numbers at the effective temperature
$T_{\mathrm{vib}}$. For brevity, let us assume that the temperature satisfies
$k_{B}T_{\mathrm{vib}}/\hbar\nu\ll 1$. According to the quantum regression
theorem [34, 43], the non-zero correlation for the boson operator $b$ at
$\tau>0$ is then obtained

$\braket{b(t)b^{\dagger}(t+\tau)}=e^{-(i\nu+\gamma)\tau}.$ (31)

Introducing the momentum operator $p=i(b^{\dagger}-b)/\sqrt{2}$, one can get
its correlation function

$\braket{p(t+\tau)p(t)}=\frac{1}{2}e^{-(i\nu+\gamma)\tau},$ (32)

and the variance $\braket{p^{2}}(t)=1/2$. Now, we can get the two-time
correlation function for displacement operator
$\mathcal{D}=\exp[\lambda(b-b^{\dagger})]=\exp(i\sqrt{2}\lambda p)$ as

$\displaystyle\braket{\mathcal{D}(t+\tau)\mathcal{D}^{\dagger}(t)}$ (33a)
$\displaystyle=$ $\displaystyle\braket{e^{i\sqrt{2}\lambda
p(t+\tau)}e^{-i\sqrt{2}\lambda p(t)}},$ $\displaystyle=$
$\displaystyle\braket{e^{i\sqrt{2}\lambda[p(t+\tau)-p(t)]}}e^{\lambda^{2}[p(t+\tau),p(t)]},$
$\displaystyle=$
$\displaystyle\braket{\sum_{k}\frac{(i\sqrt{2}\lambda)^{k}}{k!}[p(t+\tau)-p(t)]^{k}}e^{\lambda^{2}[p(t+\tau),p(t)]}.$
Taking use of the Isserlis’ theorem [22], one can simplify Eq. (33a) into
$\displaystyle\braket{\mathcal{D}(t+\tau)\mathcal{D}^{\dagger}(t)}$ (33b)
$\displaystyle=$ $\displaystyle
e^{-\lambda^{2}\braket{[(p(t+\tau)-p(t)]^{2}}}e^{\lambda^{2}[p(t+\tau),p(t)]},$
$\displaystyle=$ $\displaystyle
e^{-\lambda^{2}[2\braket{p^{2}}-2\braket{p(t+\tau)p(t)}]},$ $\displaystyle=$
$\displaystyle e^{-\lambda^{2}}e^{\lambda^{2}e^{-(i\nu+\gamma)(\tau)}}.$

Similarly, we can also get the four-time correlation function for displacement
operator at $t>t_{1}>t_{2}>t_{3}$,

$\displaystyle\braket{\mathcal{D}(t)\mathcal{D}^{\dagger}(t_{1})\mathcal{D}(t_{2})\mathcal{D}^{\dagger}(t_{3})}=$
$\displaystyle\braket{e^{i\sqrt{2}\lambda p(t)}e^{-i\sqrt{2}\lambda
p(t_{1})}e^{i\sqrt{2}\lambda p(t_{2})}e^{-i\sqrt{2}\lambda p(t_{3})}},$
$\displaystyle=$ $\displaystyle
e^{-\lambda^{2}\braket{[p(t)-p(t_{1})+p(t_{2})-p(t_{3})]^{2}}}e^{-\lambda^{2}[p(t)-p(t_{1}),p(t_{2})-p(t_{3})]}e^{\lambda^{2}[p(t),p(t_{1})]}e^{\lambda^{2}[p(t_{2}),p(t_{3})]},$
$\displaystyle=$ $\displaystyle
e^{-2\lambda^{2}[2\braket{p^{2}}-\braket{p(t)p(t_{1})}+\braket{p(t)p(t_{2})}-\braket{p(t)p(t_{3})}-\braket{p(t_{1})p(t_{2})}+\braket{p(t_{1})p(t_{3})}-\braket{p(t_{2})p(t_{3})}]}$
(34a) $\displaystyle=$ $\displaystyle
e^{-2\lambda^{2}}e^{2\lambda^{2}\braket{p(t)p(t_{1})}}e^{-2\lambda^{2}\braket{p(t)p(t_{2})}}e^{2\lambda^{2}\braket{p(t)p(t_{3})}}e^{2\lambda^{2}\braket{p(t_{1})p(t_{2})}}e^{-2\lambda^{2}\braket{p(t_{1})p(t_{3})}}e^{2\lambda^{2}\braket{p(t_{2})p(t_{3})}}.$

Submitting Eq. (32) into above equation, and taking the Taylor series
expansion, we finally get the expression for the four-time correlation in the
sum form

$\displaystyle\braket{\mathcal{D}(t)\mathcal{D}^{\dagger}(t_{1})\mathcal{D}(t_{2})\mathcal{D}^{\dagger}(t_{3})}$
(34b) $\displaystyle=$ $\displaystyle
e^{-2\lambda^{2}}\sum_{\\{k\\}}\frac{(-1)^{k_{2}+k_{5}}\lambda^{2(k_{1}+k_{2}+k_{3}+k_{4}+k_{5}+k_{6})}}{\underset{j}{\prod}k_{j}!}e^{-(k_{1}+k_{2}+k_{3})(i\nu+\gamma)(t-t_{1})}e^{-(k_{2}+k_{3}+k_{4}+k_{5})(i\nu+\gamma)(t_{1}-t_{2})}e^{-(k_{3}+k_{5}+k_{6})(i\nu+\gamma)(t_{2}-t_{3})}$

.

##  Appendix B Dynamic of Electronic States

From the Langevin equation for the cavity field $a$ described by Eq. (8a), we
can take the integration with respect of time to get the formal solution

$a=\int_{0}^{t}dt_{1}\,e^{-(i\delta_{c}+\kappa)(t-t_{1})}[-ig\sum_{m}\sigma_{m}(t_{1})+\eta+(\sqrt{2\kappa_{1}}+\sqrt{2\kappa_{2}})a_{\mathrm{in}}(t_{1})].$
(35)

By plugging Eq. (35) into the formal solution of the dipole operator
$\sigma_{m}(t)$ described by Eq. (9), and taking the averages, we find

$\displaystyle\braket{\sigma_{m}}=$ $\displaystyle-
ig\int_{0}^{t}dt_{1}e^{-(i\Delta+\Gamma_{\bot})(t-t_{1})}\int_{0}^{t_{1}}dt_{2}e^{-(i\Delta_{c}+\kappa)(t_{1}-t_{2})}\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})[-ig\sum_{n\neq
m}^{N}\sigma_{n}(t_{2})-ig\sigma_{m}(t_{2})+\eta]}.$ (36a) The item
$\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1}\sigma_{n}(t)}$
reflects the induced interaction between $m$-th and $n$-th molecule via the
quantized cavity field. The item
$\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})\sigma_{m}(t_{1})}$
relates to the modification of the frequency and decay rate. When we
considering the evolution for the Pauli operators, and submitting Eq. (9) into
Eq. (36a), the average value for Pauli operator is given by
$\displaystyle\braket{\sigma_{m}}=$ $\displaystyle ig^{3}\sum_{n\neq
m}^{N}\int_{0}^{t}dt_{1}e^{-(i\Delta+\Gamma_{\bot})(t-t_{1})}\int_{0}^{t_{1}}dt_{2}e^{-(i\Delta_{c}+\kappa)(t_{1}-t_{2})}\int_{0}^{t_{2}}dt_{3}\,e^{-(i\Delta+\Gamma_{\bot})(t_{2}-t_{3})}\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})\mathcal{D}_{n}(t_{2})\mathcal{D}_{n}^{\dagger}(t_{3})a(t_{3})}$
$\displaystyle+ig^{3}\int_{0}^{t}dt_{1}e^{-(i\Delta+\Gamma_{\bot})(t-t_{1})}\int_{0}^{t_{1}}dt_{2}e^{-(i\Delta_{c}+\kappa)(t_{1}-t_{2})}\int_{0}^{t_{2}}dt_{3}\,e^{-(i\Delta+\Gamma_{\bot})(t_{2}-t_{3})}\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})\mathcal{D}_{m}(t_{2})\mathcal{D}_{m}^{\dagger}(t_{3})a(t_{3})}$
$\displaystyle-\eta
g\int_{0}^{t}dt_{1}e^{-(i\Delta+\Gamma_{\bot})(t-t_{1})}\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})\int_{0}^{t_{1}}dt_{2}e^{-(i\Delta_{c}+\kappa)(t_{1}-t_{2})}.$
(36b)

In the main text, we have roughly take the Markov approximation (under the
large vibrational relaxation condition $\gamma\gg\kappa,\gamma$) to separate
the degrees of freedom for cavity field and vibrational motion as

$\displaystyle\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})a(t_{1})}$
$\displaystyle\approx\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}\braket{a(t_{1})},$
(37a) and also the degrees of freedom for Pauli operator and vibrational
motion as
$\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})\sigma_{n}(t_{1})}\approx\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}\braket{\sigma_{n}(t_{1})}.$
(37b)

This approach implies that we have roughly take the following approximation

$\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})\mathcal{D}_{m}(t_{2})\mathcal{D}_{m}^{\dagger}(t_{3})}\approx\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}\braket{\mathcal{D}_{m}(t_{2})\mathcal{D}_{m}^{\dagger}(t_{3})}.$
(38)

In this case, all the displacement operators refer to the same vibrational
mode. To test the feasibility of this approximation, we make a further study.
According to the expression for the four-time correlation of the displacement
operator in Eq. (34a), we can take the Laplace transformation for Eq. (36b) to
obtain

$\braket{\bar{\sigma}_{m}}=ig^{3}\sum_{n\neq
m}\frac{\mathcal{\mathcal{\bar{F}}}_{m}\mathcal{\mathcal{\bar{F}}}_{n}}{i\Delta_{c}+\kappa}\braket{\bar{a}}+ig^{3}\bar{\mathcal{F}}_{2,m}\braket{\bar{a}}-i\eta
g\frac{\mathcal{\mathcal{\bar{F}}}_{m}}{s[i\Delta_{c}+\kappa]},$ (39)

where $\mathcal{\bar{F}}_{n}=\mathcal{\bar{F}}_{n}$ for identical molecules,
and

$\displaystyle\bar{\mathcal{F}}_{2,m}=$
$\displaystyle\sum_{\\{k\\}}\frac{1}{\underset{j}{\prod}k_{j}!}e^{-2\lambda^{2}}(-1)^{k_{2}+k_{5}}\lambda^{2(k_{1}+k_{2}+k_{3}+k_{4}+k_{5}+k_{6})}$
(40)
$\displaystyle\times\frac{1}{[s+i\Delta+i(k_{1}+k_{2}+k_{3})\nu+(k_{1}+k_{2}+k_{3})\gamma+\Gamma_{\perp}][\cdots(1)\cdots][s+i\Delta+i(k_{1}+k_{2}+k_{3})\nu+(k_{1}+k_{2}+k_{3})\gamma+\Gamma_{\perp}]}$

with
$[\cdots(1)\cdots]=[s+i\Delta_{c}+i(k_{2}+k_{3}+k_{4}+k_{5})\nu+(k_{2}+k_{3}+k_{4}+k_{5})\gamma+\kappa]$.

Considering the Laplace form of the cavity field described by Eq. (16), and
using the final value theorem, we finally get

$\displaystyle\braket{a}_{\mathrm{ss}}^{\prime}$
$\displaystyle=\frac{\eta\mathcal{F}_{c}^{0}(-1+Ng^{2}\mathcal{F}_{c}^{0}\mathcal{F}_{m}^{0})}{-1+Ng^{4}\mathcal{F}_{c}^{0}[\mathcal{F}_{2,m}^{0}+(N-1)\mathcal{F}_{c}^{0}(\mathcal{F}_{m}^{0})^{2}]},$
(41) $\displaystyle\braket{\sigma_{m}}_{\mathrm{ss}}^{\prime}$
$\displaystyle=-i\eta
g\frac{g^{2}\mathcal{F}_{c}^{0}[\mathcal{F}_{2,m}^{0}+(N-1)\mathcal{F}_{c}^{0}(\mathcal{F}_{m}^{0})^{2}]-\mathcal{F}_{c}^{0}\mathcal{F}_{m}^{0}}{-1+Ng^{4}\mathcal{F}_{c}^{0}[\mathcal{F}_{2,m}^{0}+(N-1)\mathcal{F}_{c}^{0}(\mathcal{F}_{m}^{0})^{2}]}.$

Figure 4: (Colored online) Cavity transmission at resonance
$\omega_{c}=\omega_{00}$ as a function of the detuning $\Delta_{c}/\nu$ with
the number $N=200$ of organic molecules at $\lambda=1$. The data is calculated
according to Eq. (41) (see red solid line) and Eq. (18b) (see blue dashed
line), respectively, The other parameters are the same with Fig. (2).

To verify the validity of the Markov approximation to adiabatically eliminate
the degree of freedom of vibrational motion, we plot the cavity transmission
spectrum at resonance $\omega_{c}=\omega_{00}$ with the number $N=200$ of
organic molecules at $\lambda=1$. As illustrated by Fig. 4, we can observe
that same behaviors performed when we consider the four-time correlation (see
red solid line) or just the two-time correlation (see blue dashed line).

##  Appendix C Steady Population

###  C.1 Dynamics

From the Heisenberg-Langevin equation in Eq. (7) (or the master equation in
Eq. (6), we can derive the dynamical equations for average number of photons
in cavity $\braket{n_{c}}=\braket{a^{\dagger}a}$, and the excited state
population for $m-$th dye molecule
$\braket{P_{m}}=\braket{\sigma_{m}^{\dagger}\sigma_{m}}$

$\displaystyle\frac{d}{dt}\braket{P_{m}}=$
$\displaystyle-\Gamma_{\|}\braket{P_{m}}-ig(\braket{a\sigma_{m}^{\dagger}}-\braket{a^{\dagger}\sigma_{m}}),$
(42a) $\displaystyle\frac{d}{dt}\braket{n_{c}}=$
$\displaystyle-2\kappa\braket{n_{c}}-ig\sum_{m}^{N}(\braket{a^{\dagger}\sigma_{m}}-\braket{a\sigma_{m}^{\dagger}})+\eta(\braket{a}+\braket{a^{\dagger}}),$
(42b)

where $\Gamma_{\|}=2\Gamma$ is the longitudinal relaxation rate. By taking
Laplace transformations for Eq. (42), we will have the results

$\displaystyle s\braket{\bar{P}_{m}}=$
$\displaystyle-\Gamma_{\|}\braket{\bar{P}_{m}}-ig(\braket{\overline{a\sigma_{m}^{\dagger}}}-\braket{\overline{a^{\dagger}\sigma_{m}}}),$
(43a) $\displaystyle s\braket{\bar{n}_{c}}=$
$\displaystyle-2\kappa\braket{\bar{n}_{c}}-ig\sum_{m}^{N}(\braket{\overline{a^{\dagger}\sigma_{m}}}-\braket{\overline{a\sigma_{m}^{\dagger}}})+\eta(\braket{\overline{a}}+\braket{\overline{a^{\dagger}}}).$
(43b)

The presence of the operator $a^{\dagger}\sigma_{m}$ in Eq. (43), relates to
the nonlinear behavior of the HTC model. In the limit of weak driving
$\eta\ll\kappa$, the cavity field in steady situation is approximately in
vacuum. Under this situation, the noise correlation between cavity field and
the molecules as well as the noise correlation between molecules will play
important roles. Therefore, we have to give the dynamics equations for the
operator $a^{\dagger}\sigma_{m}$. While considering the coupling between the
electronic states and molecular conformation, we start form the Langevin
equation for the cavity filed operator $a$ and the polarized operator
$\tilde{\sigma}_{m}$ as given in Eq. (8), to derive the dynamic equation for
$a^{\dagger}\tilde{\sigma}_{m}$

$\displaystyle\frac{d}{dt}a^{\dagger}\tilde{\sigma}_{m}=$
$\displaystyle(\frac{d}{dt}a^{\dagger})\tilde{\sigma}_{m}+a^{\dagger}\frac{d}{dt}\tilde{\sigma}_{m}$
(44a) $\displaystyle=$
$\displaystyle[i(\Delta_{c}-\Delta)-(\kappa+\Gamma_{\bot})]a^{\dagger}\tilde{\sigma}_{m}+ig\mathcal{D}_{m}^{\dagger}(\sigma_{m}^{\dagger}\sigma_{m}-a^{\dagger}a)+ig\mathcal{D}_{m}^{\dagger}\sum_{n\neq
m}^{N}\sigma_{n}^{\dagger}\sigma_{m}+\sqrt{2\kappa}A_{in}^{\dagger}\tilde{\sigma}_{m}+\sqrt{2\Gamma_{\bot}}\mathcal{D}_{m}^{\dagger}a^{\dagger}\sigma_{m}^{\mathrm{in}}.$
Integrating with respect to the time $t$ on both sides of Eq. (44a), then the
formal solution for the operator $a^{\dagger}\sigma_{m}$ will be obtained
$a^{\dagger}\sigma_{m}=\int_{0}^{t}dt_{1}e^{[i(\Delta_{c}-\Delta)-(\kappa+\Gamma_{\bot})](t-t_{1})}\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})\\{ig[P_{m}(t_{1})-n_{c}(t_{1})]+\eta\sigma_{m}(t_{1})\\}+\mathcal{X}_{m,\mathrm{ind}}^{t}+\mathcal{X}_{m,\mathrm{in}}^{t},$
(44b) with the dipole-dipole term
$\mathcal{X}_{m,\mathrm{d}-\mathrm{d}}^{t}=ig\int_{0}^{t}dt_{1}e^{[i(\Delta_{c}-\Delta)-(\kappa+\Gamma_{\bot})](t-t_{1})}\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})\sum_{n\neq
m}^{N}\sigma_{n}^{\dagger}\sigma_{m}(t_{1}),$ (44c) and the noise term
$\mathcal{X}_{\mathrm{in}}^{t}=\int_{0}^{t}dt_{1}e^{[i(\Delta_{c}-\Delta)-(\kappa+\Gamma_{\bot})](t-t_{1})}\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})[\sqrt{2\kappa}a_{in}^{\dagger}\sigma_{m}(t_{1})+\sqrt{2\Gamma_{\bot}}\mathcal{D}_{m}^{\dagger}a^{\dagger}\sigma_{m}^{\mathrm{in}}(t_{1})].$
(44d)

Since these individuals molecules are introduced to couple a quantized cavity
field, an efficient interaction between molecules will be induced as described
by Eq. (44c). To find the expression of the dipole-dipole term
$\mathcal{X}_{m,\mathrm{d}-\mathrm{d}}^{t}$, one need to consider the
expression of $\sigma_{n}^{\dagger}\sigma_{m}$,

which can be derived through the following Langevin equation
$\frac{d}{dt}\tilde{\sigma}_{n}^{\dagger}\tilde{\sigma}_{m}=-2\Gamma_{\bot}\tilde{\sigma}_{n}^{\dagger}\tilde{\sigma}_{m}+ig\mathcal{D}_{m}^{\dagger}\mathcal{D}_{n}(a^{\dagger}\sigma_{m}-a\sigma_{n}^{\dagger})+\sqrt{2\Gamma_{\bot}}\mathcal{D}_{m}^{\dagger}\mathcal{D}_{n}(\sigma_{n}^{\dagger}\sigma_{m}^{\mathrm{in}}+\sigma_{n}^{\mathrm{in}\dagger}\sigma_{m}).$
(45a) The evolution of $\sigma_{n}^{\dagger}\sigma_{m}(t)$ is then similarly
obtained as
$\sigma_{n}^{\dagger}\sigma_{m}=ig\int_{0}^{t}dt_{1}e^{-2\Gamma_{\bot}(t-t_{1})}\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})[a^{\dagger}\sigma_{m}(t_{1})-a\sigma_{n}^{\dagger}(t_{1})]\mathcal{D}_{n}(t_{1})\mathcal{D}_{n}^{\dagger}(t)+\mathcal{Y}_{\mathrm{in}}^{t},$
(45b) with the noise term
$\mathcal{Y}_{\mathrm{in}}^{t}=\sqrt{2\Gamma_{\bot}}\int_{0}^{t}dt_{1}e^{-2\Gamma_{\bot}(t-t_{1})}\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})[\sigma_{n}^{\dagger}\sigma_{m}^{\mathrm{in}}(t_{1})+\sigma_{n}^{\mathrm{in}\dagger}\sigma_{m}(t_{1})]\mathcal{D}_{n}(t_{1})\mathcal{D}_{n}^{\dagger}(t).$
(45c)

From now on, the formal solution for the interactions between cavity and
molecules characterized by Eq. (44b) as well as the interactions between
molecules depicted by Eq. (45b) have been obtained. To get the excited state
population for total molecules, one can take the Markov approximation to
adiabatically eliminate the vibrational motion without (or with) submitting
Eq. (45b) into Eq. (44b), which relates to the perturbative treatment in the
first (or second) order.

###  C.2 Perturbative Treatment in the 1st Order

In this subsection, let us take the perturbative treatment in the first order.
Under the Markov approximation with the large vibrational relaxation, we can
adiabatically eliminate the vibrational motion as introduced in Sec. II to get
the dynamic equations for the item $\braket{a^{\dagger}\sigma_{m}}$ and
$\braket{\sigma_{n}^{\dagger}\sigma_{m}}$ in Laplace space

$\displaystyle\braket{\overline{a^{\dagger}\sigma_{m}}}=$ $\displaystyle
ig\bar{\mathcal{F}}_{m}^{\prime}(\braket{\bar{P}_{m}}-\braket{\bar{n}_{c}}+\sum_{n\neq
m}^{N}\braket{\overline{\sigma_{n}^{\dagger}\sigma_{m}}})+\eta\bar{\mathcal{F}}_{m}^{\prime}\braket{\bar{\sigma}_{m}}$
(46a) $\displaystyle\braket{\overline{\sigma_{n}^{\dagger}\sigma_{m}}}=$
$\displaystyle
ig\bar{\mathcal{F}}_{mn}(\braket{\overline{a^{\dagger}\sigma_{m}}}-\braket{\overline{a\sigma_{n}^{\dagger}}}),$
(46b)

where $\bar{\mathcal{F}}_{m}^{\prime}$ is the Laplace form for
$\mathcal{F}_{m}^{\prime}(t)=\exp[(i\Delta_{c}-i\Delta-\kappa-\Gamma_{\bot})t]\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(0)}$
with the expression

$\bar{\mathcal{F}}_{m}^{\prime}=\sum_{k}\frac{\lambda^{2k}}{k!}\frac{e^{-\lambda^{2}}}{s+i(\Delta-\Delta_{c}+k\nu)+(\kappa+\Gamma_{\bot}+k\gamma)},$
(47)

and $\bar{\mathcal{F}}_{mn}$ is the Laplace form for
$\mathcal{F}_{mn}(t)=\exp(-2\Gamma_{\bot}t)\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(0)}\braket{\mathcal{D}_{n}(0)\mathcal{D}_{n}^{\dagger}(t)}$
with the expression

$\mathcal{\bar{F}}_{mn}=\sum_{k_{m},k_{n}}\frac{\lambda^{2(k_{m}+k_{n})}}{k_{m}!k_{n}!}\frac{e^{-2\lambda^{2}}}{s+i(k_{m}-k_{n})\nu+[2\Gamma_{\bot}+(k_{m}+k_{n})\gamma]}.$
(48)

For identical molecules, the Laplace forms for the population of the excited
electric state $\braket{\overline{\sigma_{m}^{\dagger}\sigma_{m}}}$, the
expect value for the interaction between cavity and molecule
$\braket{\overline{a^{\dagger}\sigma_{m}}}$ as well as the function
$\bar{\mathcal{F}}_{m}^{\prime}$ are the same with that for any molecule
($m$), and the expect value for the intra-molecule interaction
$\braket{\overline{\sigma_{n}^{\dagger}\sigma_{m}}}$, as well as the function
$\mathcal{\bar{F}}_{mn}$ are the same for any molecule pair ($m$, $n$). As a
result, we can simplify Eqs. (43), and (46), and transform them into matrix
form

$\mathcal{M}_{1}\bm{V}_{1}+\eta\bm{V}_{1,\mathrm{in}}=0,$ (49)

with the drift matrix

$\mathcal{M}_{1}=\left(\begin{array}[]{ccccc}-2\kappa-s&iNg&-iNg&0&0\\\
ig\bar{\mathcal{F}}_{m}^{\prime*}&-1&0&-ig\bar{\mathcal{F}}_{m}^{\prime*}&-i(N-1)g\bar{\mathcal{F}}_{m}^{\prime*}\\\
-ig\bar{\mathcal{F}}_{m}^{\prime}&0&-1&ig\bar{\mathcal{F}}_{m}^{\prime}&i(N-1)g\bar{\mathcal{F}}_{m}^{\prime}\\\
0&-ig&ig&-\Gamma_{\|}-s&0\\\
0&-ig\bar{\mathcal{F}}_{mn}&ig\bar{\mathcal{F}}_{mn}&0&-1\end{array}\right),$
(50)

the vector $\bm{V}_{1}=(\braket{\bar{n}_{c}},$
$\braket{\overline{a\sigma_{m}^{\dagger}}},$
$\braket{\overline{a^{\dagger}\sigma_{m}}},$ $\braket{\bar{P}_{m}},$
$\braket{\overline{\sigma_{m}^{\dagger}\sigma_{n}}})^{T}$ and the input vector
$\bm{V}_{1,\mathrm{in}}=(\braket{\overline{a}}+\braket{\overline{a^{\dagger}}},$
$\bar{\mathcal{F}}_{m}^{\prime*}\braket{\bar{\sigma}_{m}^{\dagger}},$
$\bar{\mathcal{F}}_{m}^{\prime}\braket{\bar{\sigma}_{m}},$ $0,$ $0)^{T}$.

We then obtain the solutions
$\bm{V}_{1}=-\eta\mathcal{M}_{1}^{-1}\bm{V}_{1,\mathrm{in}}$ by direct matrix
inversion. According to the final value theorem, the excited state population
for all molecules is finally given by

$\mathcal{P}_{N}=\sum_{n=1}^{N}\braket{P_{m}}.$ (51)

If we set $N=1$, there will be only one molecule couples to the cavity and the
inter-molecule interaction induced by cavity [depicted by Eq. (44c)] can be
neglected. In such situation, the expect value for excited state population in
steady state will be simplified as

$\displaystyle\mathcal{P}_{1}$ $\displaystyle=\frac{\eta
g^{2}\mathfrak{R}\braket{a}_{\mathrm{ss}}\mathfrak{R}\mathcal{F}_{m}^{\prime
0}-\eta g\kappa\mathfrak{I}\\{\mathcal{F}_{m}^{\prime
0}\braket{\sigma_{m}}_{\mathrm{ss}}\\}}{\Gamma\kappa+g^{2}(\Gamma+\kappa)\mathfrak{R}\mathcal{F}_{m}^{\prime
0}},$ (52)

where $\mathcal{F}_{m}^{\prime 0}=\underset{s\rightarrow
0}{\lim}\bar{\mathcal{F}}_{m}^{\prime}$.

Figure 5 illustrates the excited state population for all molecules at
resonance $\omega_{c}=\omega_{00}$ with and without the incoherent coupling
between molecules. We compare the HTC model (red solid line) for varying
number of molecule $N$, against the TC model (blue dashed line) in the limit
$\alpha=0$ (a-d) and $\alpha=1$ (e-h). For finite $\lambda$, the transfer rate
of population from upper polaritonic states to lower polaritonic states and
dark polaritonic states increased with the number of organic molecules
increasing. Consequently the population for upper polaritonic states will be
reduced, and the population for dark polaritonic states is going to be
increased. In addition, the transfer rate for the anti-Stokes process
$\ket{\textsl{g}}\ket{0}_{\mathrm{vib}}\rightarrow\ket{e}\ket{1}_{\mathrm{vib}}$
is enhanced with the increasing of the number $N$ of organic molecules. This
reflects that the molecular conformation has been modified via the collective
behavior.

Figure 5: (Colored online) Steady population. The excited state population for
all molecules at $\omega_{c}=\omega_{00}$ in steady state as a function of the
detuning $\Delta_{c}/\nu$ with different number $N$ (increases as $2$, $20$
and $200$) of organic molecules at $\lambda=0.6$ (red solid line) and
$\lambda=0$ (blue dashed line) when (a-c) $\alpha=0$ and (d-f) $\alpha=1$. (a,
d) The black dotted line shows the numerical simulation results for the case
$\lambda=0.6$, $N=2$ with the help of QuTiP. $P_{1}^{\prime}$ is the unit
value given by Eq. (52) at $\Delta=0,$ $\lambda=0$. The other parameters are
the same with Fig. 2.

###  C.3 Perturbative Treatment in the 2nd Order

To test the reasonableness of the approximation as previously mentioned, we
will take perturbative treatment in higher order. Substituting Eq. (45b) into
Eq. (44c), and taking the averages, we find

$\displaystyle\braket{\mathcal{X}_{m,\mathrm{d}-\mathrm{d}}^{t}}\approx$
$\displaystyle-g^{2}\sum_{n\neq
m}^{N}\int_{0}^{t}dt_{1}\,e^{[i(\Delta_{c}-\Delta)-(\kappa+\Gamma_{\bot})](t-t_{1})}\int_{0}^{t_{1}}dt_{2}\,e^{-2\Gamma_{\bot}(t_{1}-t_{2})}\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{2})}\braket{\mathcal{D}_{n}(t_{2})\mathcal{D}_{n}^{\dagger}(t_{1})}$
$\displaystyle\times[\braket{a^{\dagger}\sigma_{m}}(t_{2})-\braket{a\sigma_{n}^{\dagger}}(t_{2})],$
(53)

Here we have adiabatically eliminated the vibrational motion under the large
vibrational relaxation condition (i.e. $\gamma\gg\kappa$ and
$\gamma\gg\Gamma)$. Noticing that the correlation function
$\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{2})}$ (and
$\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(t_{1})}$) depends on the
time difference $t-t_{2}$ (and $t-t_{1}$) as described in Eq. (13), we can
reformed Eq. (53) as

$\displaystyle\braket{\mathcal{X}_{m,\mathrm{ind}}^{t}}=$
$\displaystyle-g^{2}e^{-\lambda^{2}}\sum_{n\neq
m}^{N}\int_{0}^{t}dt_{1}\,e^{[i(\Delta_{c}-\Delta)-(\kappa+\Gamma_{\bot})](t-t_{1})}\int_{0}^{t_{1}}dt_{2}\,e^{-2\Gamma_{\bot}(t_{1}-t_{2})}e^{\lambda^{2}e^{-(\gamma+i\nu)(t-t_{2})}}\braket{\mathcal{D}_{n}(t_{2})\mathcal{D}_{n}^{\dagger}(t_{1})}$
$\displaystyle\times[\braket{a^{\dagger}\sigma_{m}}(t_{2})-\braket{a\sigma_{n}^{\dagger}}(t_{2})].$
$\displaystyle=$ $\displaystyle-g^{2}\sum_{n\neq
m}^{N}\sum_{k}\int_{0}^{\infty}dt_{1}\int_{0}^{\infty}dt_{2}\,\mathcal{F}_{m}^{(k)}(t-t_{1})\mathcal{F}_{mn}^{(k)}(t_{1}-t_{2})[\braket{a^{\dagger}\sigma_{m}}(t_{2})-\braket{a\sigma_{n}^{\dagger}}(t_{2})],$
(54)

where $\tilde{\Delta}_{k}=\Delta_{c}-\Delta+k\nu$,
$\tilde{\Gamma}_{1,k}=\kappa+\Gamma_{\bot}+k\gamma$,
$\tilde{\Gamma}_{2,k}=2\Gamma_{\bot}+k\gamma$,
$\mathcal{F}_{m}^{(k)}(t)=\Theta(t)\lambda^{2k}\exp(-\lambda^{2})\exp[(i\tilde{\Delta}_{k}-\tilde{\Gamma}_{1,k})t]/{k!},$
and
$\mathcal{F}_{mn}^{(k)}(t_{1}-t_{2})=\Theta(t_{1}-t_{2})\exp[-(ik\nu+\tilde{\Gamma}_{2,k})(t_{1}-t_{2})]\braket{\mathcal{D}_{n}(t_{2})\mathcal{D}_{n}^{\dagger}(t_{1})}.$
Since Eq. (54) represents a convolution, one can take use the Laplace
transformation to rewritten Eq. (53) into the Laplace space:

$\braket{\bar{\mathcal{X}}_{m,\mathrm{ind}}^{t}}=-g^{2}\sum_{n\neq
m}^{N}\mathcal{\bar{F}}_{mn}[\braket{\overline{a^{\dagger}\sigma_{m}}}-\braket{\overline{a\sigma_{n}^{\dagger}}}],$
(55)

where

$\displaystyle\mathcal{\bar{F}}_{mn}=$
$\displaystyle\underset{k}{\sum}\mathcal{\bar{F}}_{m}^{(k)}\mathcal{\bar{F}}_{mn}^{(k)},$
(56) $\displaystyle=$
$\displaystyle\sum_{k_{m},k_{n}}\frac{\lambda^{2(k_{m}+k_{n})}}{k_{m}!k_{n}!}\frac{e^{-2\lambda^{2}}}{[s+i(\Delta-\Delta_{c}+k_{m}\nu)+(\kappa+\Gamma_{\bot}+k_{m}\gamma)][s+i(k_{m}-k_{n})\nu+2\Gamma_{\bot}+(k_{m}+k_{n})\gamma]}.$

Via substituting Eq. (9) into Eq. (44b), and adiabatically eliminating the
degree of freedom of vibrational motion, we similarly obtain the evolution of
$\braket{a^{\dagger}\sigma_{m}}$ in Laplace space

$\braket{\overline{a^{\dagger}\sigma_{m}}}=ig[\bar{\mathcal{F}}_{m}^{\prime}(\braket{\overline{P}_{m}}-\braket{\bar{n}_{c}})-\eta\mathcal{\bar{F}}_{cm}^{\prime}\braket{\bar{a}}]-g^{2}\sum_{n\neq
m}^{N}\mathcal{\bar{F}}_{mn}[\braket{\overline{a^{\dagger}\sigma_{m}}}-\braket{\overline{a\sigma_{n}^{\dagger}}}],$
(57)

where $\bar{\mathcal{F}}_{m}^{\prime}$ is the Laplace form for
$\mathcal{F}_{m}^{\prime}(t)=\exp[(i\Delta_{c}-i\Delta-\kappa-\Gamma_{\bot})t]\braket{\mathcal{D}_{m}(t)\mathcal{D}_{m}^{\dagger}(0)}$
with the expression

$\bar{\mathcal{F}}_{m}^{\prime}=\sum_{k}\frac{\lambda^{2k}}{k!}\frac{e^{-\lambda^{2}}}{s+i(\Delta-\Delta_{c}+k\nu)+(\kappa+\Gamma_{\bot}+k\gamma)},$
(58)

and

$\mathcal{\bar{F}}_{cm}^{\prime}=\sum_{k}\frac{\lambda^{2k}}{k!}\frac{e^{-\lambda^{2}}}{[s+i(\Delta-\Delta_{c}+k\nu)+(\kappa+\Gamma_{\bot}+k\gamma)][s+i(\Delta+k\nu)+(\Gamma_{\bot}+k\gamma)]}.$
(59)

For identical molecules, the Laplace forms for the population of the excited
electric state $\braket{\overline{\sigma_{m}^{\dagger}\sigma_{m}}}$, the
expect value for the interaction between cavity and molecule
$\braket{\overline{a^{\dagger}\sigma_{m}}}$ as well as the function
$\bar{\mathcal{F}}_{m}^{\prime}$ are the same with that for any molecule
($m$), and the expect value for the intra-molecule interaction
$\braket{\overline{\sigma_{n}^{\dagger}\sigma_{m}}}$, as well as the function
$\mathcal{\bar{F}}_{mn}$ are the same for any molecule pair ($m$, $n$). As a
result, we can simplify Eqs. (43), and (57) into matrix form

$\mathcal{M}_{2}\bm{V}_{2}+\eta\bm{V}_{2,\mathrm{in}}=0,$ (60)

with the drift matrix

$\mathcal{M}_{2}=\left(\begin{array}[]{cccc}-2\kappa-s&-iNg&iNg&0\\\
-ig\bar{\mathcal{F}}_{m}^{\prime}&-(N-1)g^{2}\mathcal{\bar{F}}_{mn}-1&(N-1)g^{2}\mathcal{\bar{F}}_{mn}&ig\bar{\mathcal{F}}_{m}^{\prime}\\\
ig\bar{\mathcal{F}}_{m}^{\prime*}&(N-1)g^{2}\mathcal{\bar{F}}_{mn}^{*}&-(N-1)g^{2}\bar{\mathcal{F}}_{mn}^{*}-1&-ig\bar{\mathcal{F}}_{m}^{\prime*}\\\
0&ig&-ig&-\Gamma_{\|}-s\end{array}\right),$ (61)

the vector $\bm{V}_{2}=(\braket{\bar{n}_{c}},$
$\braket{\overline{a^{\dagger}\sigma_{m}}},$
$\braket{\overline{a\sigma_{m}^{\dagger}}},$ $\braket{\bar{P}_{m}})^{T}$ and
the input vector
$\bm{V}_{2,\mathrm{in}}=(\braket{\bar{a}}+\braket{\overline{a^{\dagger}}},-ig\bar{\mathcal{F}}_{cm}^{\prime}\braket{\bar{a}},ig\bar{\mathcal{F}}_{m}^{\prime*}\braket{\overline{a^{\dagger}}},0)^{T}$.
We then obtain the solutions
$\bm{V}_{2}=-\eta\mathcal{M}_{2}^{-1}\bm{V}_{\mathrm{2,in}}$ by direct matrix
inversion. Finally, the excited state population for total molecules is given

$\mathcal{P}^{\prime\prime}_{N}=\sum_{n=1}^{N}\braket{P_{m}}.$ (62)

If we set $N=1$, the expect value for excited state population in steady state
will be simplified as

$\displaystyle\mathcal{P}^{\prime\prime}_{1}$ $\displaystyle=\eta
g^{2}\frac{\mathfrak{R}\braket{a}_{\mathrm{ss}}\mathfrak{R}\mathcal{F}_{m}^{\prime
0}+\kappa\mathfrak{R}\\{\mathcal{F}_{cm}^{\prime
0}\braket{a}_{\mathrm{ss}}\\}}{\Gamma\kappa+g^{2}(\Gamma+\kappa)\mathfrak{R}\mathcal{F}_{m}^{\prime
0}},$ (63)

where $\mathcal{F}_{m}^{\prime 0}=\underset{s\rightarrow
0}{\lim}\bar{\mathcal{F}}_{m}^{\prime}$, and $\mathcal{F}_{cm}^{\prime
0}=\underset{s\rightarrow 0}{\lim}\bar{\mathcal{F}}_{cm}^{\prime}$.

Figure 6: (Colored online) Steady population. The excited state population for
total molecules at $\omega_{c}=\omega_{00}$ in steady state as a function of
the detuning $\Delta_{c}/\nu$ with the number $N=200$ of organic molecules at
$\lambda=1$. We compare the perturbative method in 1st order (blue dashed
line) against the 2nd order (red solid line). $P_{1}^{\prime}$ is the unit
value given by Eq. (52) at $\Delta=0,$ $\lambda=0$. The other parameters are
the same with Fig. 2.

To test the reasonableness of the Markov approximation to adiabatically
eliminate the degree of freedom for the vibrational motion, we plot the
excited steady population for total molecules various the detuning
$\Delta_{c}/\nu$ at $\omega_{c}=\omega_{00}$ with the number $N=200$ of
organic molecules at $\lambda=1$. We can observe that similar behaviors
performed when we take the perturbative method in second order (see red solid
line) and first order (see blue dashed line).

## References

  * Klessinger and Michl [1995] M. Klessinger and J. Michl,  _Excited States and Photo-Chemistry of Organic Molecules,_ (Wiley, 1995).
  * Aroca [2013] R. F. Aroca, Phys. Chem. Chem. Phys. 15, 5355 (2013).
  * Doppagne  _et al._ [2017] B. Doppagne, M. C. Chong, E. Lorchat, S. Berciaud, M. Romeo, H. Bulou, A. Boeglin, F. Scheurer, and G. Schull, Phys. Rev. Lett. 118, 127401 (2017).
  * Wrigge  _et al._ [2008] G. Wrigge, I. Gerhardt, J. Hwang, G. Zumofen, and V. Sandoghdar, Nat. Phys. 4, 60 (2008).
  * Kühn  _et al._ [2006] S. Kühn, U. Håkanson, L. Rogobete, and V. Sandoghdar, Phys. Rev. Lett. 97, 017402 (2006).
  * Lidzey  _et al._ [1998] D. G. Lidzey, D. D. C. Bradley, M. S. Skolnick, T. Virgili, S. Walker, and D. M. Whittaker, Nature 395, 53 (1998).
  * Holmes and Forrest [2007] R. Holmes and S. Forrest, Org. Electron. 8, 77 (2007).
  * Tischler  _et al._ [2005] J. R. Tischler, M. S. Bradley, V. Bulović, J. H. Song, and A. Nurmikko, Phys. Rev. Lett. 95, 036401 (2005).
  * Kéna-Cohen  _et al._ [2008] S. Kéna-Cohen, M. Davanço, and S. R. Forrest, Phys. Rev. Lett. 101, 116401 (2008).
  * Herrera and Owrutsky [2020] F. Herrera and J. Owrutsky, J. Chem. Phys. 152, 100902 (2020).
  * Herrera and Spano [2016] F. Herrera and F. C. Spano, Phys. Rev. Lett. 116, 238301 (2016).
  * Du  _et al._ [2019] M. Du, R. F. Ribeiro, and J. Yuen-Zhou, Chem 5, 1167 (2019).
  * Schachenmayer  _et al._ [2015] J. Schachenmayer, C. Genes, E. Tignone, and G. Pupillo, Phys. Rev. Lett. 114, 196403 (2015).
  * Hagenmüller  _et al._ [2017] D. Hagenmüller, J. Schachenmayer, S. Schütz, C. Genes, and G. Pupillo, Phys. Rev. Lett. 119, 223601 (2017).
  * Zhong  _et al._ [2016] X. Zhong, T. Chervy, S. Wang, J. George, A. Thomas, J. A. Hutchison, E. Devaux, C. Genet, and T. W. Ebbesen, Angew. Chem. Int. Edit. 55, 6202 (2016).
  * Reitz  _et al._ [2018] M. Reitz, F. Mineo, and C. Genes, Scientific Reports 8, 9050 (2018).
  * Reitz  _et al._ [2019] M. Reitz, C. Sommer, and C. Genes, Phys. Rev. Lett. 122, 203602 (2019).
  * Zeb  _et al._ [2018] M. A. Zeb, P. G. Kirton, and J. Keeling, ACS Photonics 5, 249 (2018).
  * Klaers  _et al._ [2010] J. Klaers, J. Schmitt, F. Vewinger, and M. Weitz, Nature 468, 545 (2010).
  * May and Kuhn [2011] V. May and O. Kuhn,  _Charge and Energy Transfer Dynamics in Molecular Systems_ , 3rd ed. (John Wiley & Sons, Ltd, 2011).
  * Reitz  _et al._ [2020] M. Reitz, C. Sommer, B. Gurlek, V. Sandoghdar, D. Martin-Cano, and C. Genes, Phys. Rev. Research 2, 033270 (2020).
  * Reitz and Genes [2020] M. Reitz and C. Genes, J. Chem. Phys. 153, 234305 (2020).
  * Zhang  _et al._ [2020] Y. Zhang, J. Aizpurua, and R. Esteban, ACS Photonics 7, 1676 (2020).
  * Johansson  _et al._ [2013] J. Johansson, P. Nation, and F. Nori, Comput. Phys. Commun. 184, 1234 (2013).
  * Sawant  _et al._ [2018] R. Sawant, O. Dulieu, and S. A. Rangwala, Phys. Rev. A 97, 063405 (2018).
  * Schwartz  _et al._ [2011] T. Schwartz, J. A. Hutchison, C. Genet, and T. W. Ebbesen, Phys. Rev. Lett. 106, 196405 (2011).
  * Miller  _et al._ [1980] W. H. Miller, N. C. Handy, and J. E. Adams, J. Chem. Phys. 72, 99 (1980).
  * Calvo  _et al._ [2020] J. Calvo, D. Zueco, and L. Martin-Moreno, Nanophotonics 9, 277 (2020).
  * Sánchez-Carrera  _et al._ [2010] R. S. Sánchez-Carrera, M. C. R. Delgado, C. C. Ferrón, R. M. Osuna, V. Hernández, J. T. L. Navarrete, and A. Aspuru-Guzik, Organic Electronics 11, 1701 (2010).
  * Herrera and Spano [2017] F. Herrera and F. C. Spano, Phys. Rev. A 95, 053867 (2017).
  * Herrera and Spano [2018] F. Herrera and F. C. Spano, ACS Photonics 5, 65 (2018).
  * Ćwik  _et al._ [2016] J. A. Ćwik, P. Kirton, S. De Liberato, and J. Keeling, Phys. Rev. A 93, 033840 (2016).
  * Arnardottir  _et al._ [2020] K. B. Arnardottir, A. J. Moilanen, A. Strashko, P. Törmä, and J. Keeling, Phys. Rev. Lett. 125, 233603 (2020).
  * Crispin Gardiner [2004] P. Z. Crispin Gardiner,  _Quantum noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics_ , 2nd ed., Springer series in synergetics (Springer-Verlag, Berlin Heidelberg, 2004).
  * Nitzan  _et al._ [1975] A. Nitzan, S. Mukamel, and J. Jortner, J. Chem. Phys. 63, 200 (1975).
  * Radonjić  _et al._ [2018] M. Radonjić, W. Kopylov, A. Balaž, and A. Pelster, New J. Phys. 20, 055014 (2018).
  * Li  _et al._ [2006] J.-H. Li, J.-B. Liu, A.-X. Chen, and C.-C. Qi, Phys. Rev. A 74, 033816 (2006).
  * Plankensteiner  _et al._ [2019] D. Plankensteiner, C. Sommer, M. Reitz, H. Ritsch, and C. Genes, Phys. Rev. A 99, 043843 (2019).
  * Bernardot  _et al._ [1992] F. Bernardot, P. Nussenzveig, M. Brune, J. M. Raimond, and S. Haroche, EPL 17, 33 (1992).
  * Khitrova  _et al._ [2006] G. Khitrova, H. M. Gibbs, M. Kira, S. W. Koch, and A. Scherer, Nat. Phys. 2, 81 (2006).
  * Wellnitz  _et al._ [2020] D. Wellnitz, S. Schütz, S. Whitlock, J. Schachenmayer, and G. Pupillo, Phys. Rev. Lett. 125, 193201 (2020).
  * Galego  _et al._ [2019] J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. X 9, 021057 (2019).
  * Carmichael [1999] H. J. Carmichael,  _Statistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations_ , Theoretical and Mathematical Physics (Springer-Verlag, Berlin Heidelberg, 1999).
  * Sauer  _et al._ [2011] M. Sauer, J. Hofkens, and J. Enderlein,  _Handbook of Fluorescence Spectroscopy and Imaging: From Single Molecules to Ensembles_ (Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2011).


* (27) Gorban, A.N. Maxallent: Maximizers of all entropies and uncertainty of uncertainty. Comput. Math. Appl. 2013, 65, 1438–1456. Corrected postprint: arXiv:1212.5142 [physics.data-an].
  * (28) Golub, G.H.; Van Loan, C.F. Matrix Computations; Johns Hopkins University Press: Baltimore, MD, USA, 1996.
  * (29) Gorban, A.N. Equilibrium encircling. Equations of Chemical Kinetics and Their Thermodynamic Analysis; Nauka: Novosibirsk, USSR, 1984. (in Russian)
  * (30) Yablonskii, G.S.; Bykov, V.I.; Gorban, A.N.; Elokhin, V.I. Kinetic Models of Catalytic Reactions; Series “Comprehensive Chemical Kinetics,” Volume 32; Elsevier: Amsterdam, The Netherlands, 1991.
  * (31) Grmela, M.; Öttinger, H.C. Dynamics and thermodynamics of complex fluids. I. Development of a general formalism. Phys. Rev. E 1997, 56, 6620.
  * (32) Grmela, M. Fluctuations in extended mass-action-law dynamics. Physica D 2012, 241, 976–986.
  * (33) Giovangigli, V.; Matuszewski, L. Supercritical fluid thermodynamics from equations of state. Physica D 2012, 241 649–670.
  * (34) Gorban, A.N.; Shahzad, M. The Michaelis-Menten-Stueckelberg Theorem. Entropy 2011, 13, 966–1019. Corrected postprint: arXiv:1008.3296 [physics.chem-ph].
  * (35) Hangos, K.M. Engineering model reduction and entropy-based Lyapunov functions in chemical reaction kinetics. Entropy 2010, 12, 772–797.
  * (36) Orlov, N.N.; Rozonoer, L.I. The macrodynamics of open systems and the variational principle of the local potential. J. Franklin Inst-Eng. Appl. Math. 1984, 318, 283–341.
  * (37) Feinberg, M. On chemical kinetics of a certain class. Arch. Rat. Mechan. Anal. 1972, 46, 1–41.
  * (38) Gorban, A.N.; Mirkes, E.M.; Yablonsky, G.S. Thermodynamics in the limit of irreversible reactions. Physica A 2013, 392, 1318–1335. arXiv:1207.2507 [cond-mat.stat-mech].
  * (39) Boltzmann, L. Neuer Beweis zweier Sätze über das Wärmegleichgewicht unter mehratomigen Gasmolekülen. Sitzungsberichte der Kaiserlichen Akademie der Wissenschaften in Wien. 1887, 95, 153-164. (in German)
  * (40) Lorentz, H.-A. Über das Gleichgewicht der lebendigen Kraft unter Gasmolekülen. Sitzungsberichte der Kaiserlichen Akademie der Wissenschaften in Wien. 1887, 95, 115–152. (in German)
  * (41) Stueckelberg, E.C.G. Theoreme $H$ et unitarite de $S$. Helv. Phys. Acta 1952, 25, 577–580.
  * (42) Horn, F.; Jackson, R. General mass action kinetics. Arch. Ration. Mech. Anal. 1972, 47, 81–116.
  * (43) Szederkényi, G.; Hangos, K.M. Finding complex balanced and detailed balanced realizations of chemical reaction networks. J. Math. Chem. 2011, 49, 1163–1179. arXiv:1010.4477 [q-bio.MN].
  * (44) Kalpazidou, S.L. Cycle Representations of Markov Processes; Book series: Applications of Mathematics, Volume 28; Springer: New York, NY, USA, 2006.
  * (45) Rockafellar, R.T. Convex Analysis; Princeton University Press: Princeton, NJ, USA, 1997.
  * (46) Kolokoltsov, V.N. Nonlinear Markov Processes and Kinetic Equations; Book series Cambridge Tracts in Mathematics, Volume 182; Cambridge University Press: Cambridge, UK, 2010.
  * (47) Greenberg, H.J.; Pierskalla, W.P. A review of quasi-convex functions. Oper. Res. 1971, 19, 1553–1570.
  * (48) Ponstein, J. Seven kinds of convexity. SIAM Rev. 1967, 9, 115–119.
  * (49) Matoušek, J. On directional convexity. Discret. Comput. Geom. 2001, 25, 389–403.
  * (50) Hwang, F.K.; Rothblum, U.G. Directional-quasi-convexity, asymmetric Schur-convexity and optimality of consecutive partitions. Math. Oper. Res. 1996, 21, 540–554.
  * (51) Marshall, A.W.; Olkin, I.; Arnold, B.C. Inequalities: Theory of Majorization and its Applications; Springer: New York, NY, USA, 2011.
  * (52) Gorban, A.N.; Yablonskii, G.S. Extended detailed balance for systems with irreversible reactions. Chem. Eng. Sci. 2011, 66, 5388–5399. arXiv:1101.5280 [cond-mat.stat-mech].
  * (53) Gorban, A.N. Equilibrium Encircling. Equations of Chemical kinetics and their Thermodynamic Analysis; Nauka: Novosibirsk, Russia, 1984. (in Russian)
  * (54) Gorban, A.N. Invariant sets for kinetic equations. React. Kinet. Catal. Lett. 1979, 10, 187–190.
  * (55) Gorban, A.N. Thermodynamic Tree: The Space of Admissible Paths. SIAM J. Applied Dynamical Systems 2013, 12 (1), 246–278. arXiv:1201.6315 [cond-mat.stat-mech].
  * (56) Jaynes, E.T. Gibbs vs Boltzmann entropies. American Journal of Physics 1965, 33 (5), 391–398. Reprint online: http://lptms.u-psud.fr/membres/trizac/Ens/M2MQPL/Jaynes_entropy.pdf Accessed June 03 2014.
  * (57) Gorban, A.N.; Bykov, V.I.; Yablonskii, G.S., Essays on chemical relaxation, Novosibirsk: Nauka, 1986. (in Russian)
  * (58) Lebowitz, J.L. Boltzmann’s entropy and time’s arrow. Physics Today 1993, 46, 32–32. Reprint online: http://users.df.uba.ar/ariel/materias/FT3_2008_1C/papers_pdf/lebowitz_370.pdf Accessed June 03 2014
  * (59) Goldstein, S.; Lebowitz, J. L. On the (Boltzmann) entropy of non-equilibrium systems. Physica D: Nonlinear Phenomena 2004, 193 (1), 53–66. arXiv:cond-mat/0304251 [cond-mat.stat-mech]
  * (60) Csiszár, I.; Matúš, F. Information Projections Revisited. IEEE Transactions On Information Theory 2003, 49 (6), 1474–1490. Reprint online: http://yaroslavvb.com/papers/cziszar-information.pdf Accessed June 03 2014
  * (61) Gorban, A.N.; Kaganovich, B.M.; Filippov, S.P.; Keiko, A.V.; Shamansky, V.A.; Shirkalin, I.A. Thermodynamic Equilibria and Extrema: Analysis of Attainability Regions and Partial Equilibria. Springer, Berlin-Heidelberg-New York, 2006.
  * (62) Gorban, A.N.; Karlin, I.V. Uniqueness of thermodynamic projector and kinetic basis of molecular individualism. Physica A 2004, 336 (3-4), 391–432. arXiv:cond-mat/0309638
  * (63) Gorban, A.N.; Gorban, P.A.; Karlin, I.V. Legendre integrators, post-processing and quasiequilibrium. J. Non-Newtonian Fluid Mech. 2004, 120, 149–167. arXiv:cond-mat/0308488
  * (64) Gorban, A.N.; Karlin, I.V. Method of invariant manifold for chemical kinetics. Chem. Eng. Sci. 2003 58, 4751–4768. arXiv:cond-mat/0207231 [cond-mat.stat-mech]
  * (65) Gorban, A.N.; Karlin, I.V. Invariant Manifolds for Physical and Chemical Kinetics. Lect. Notes Phys. 660, Springer, Berlin–Heidelberg, 2005.
  * (66) Gorban, A.N. On the problem of boundary equilibria. React. Kinet. Catal. Lett. 1980, 15 (3), 315–319. (in Russian)
  * (67) Zylka, Ch. A note on the attainability of states by equalizing processes. Theor. Chim. Acta 1985, 68, 363–377.
  * (68) M. Gopalkrishnan, E. Miller, A. Shiu, A geometric approach to the Global Attractor Conjecture, 2013, arXiv:1305.5303 [math.DS]