#  Effects of disorder on polaritonic and dark states in a cavity using the
disordered Tavis-Cummings model

Tarun Gera and K.L. Sebastian
Department of Inorganic and Physical Chemistry
Indian Institute of Science
Bangalore 560012
India

###### Abstract

We consider molecules confined to a microcavity whose dimensions are such that
an excitation of the molecule is nearly resonant with a cavity mode. We
investigate the situation where the excitation energies of the molecules are
randomly distributed with a mean value of $\epsilon_{a}$ and variance
$\sigma$. For this case, we find a solution that approaches the exact result
for large values of the number density $\mathscr{N}$ of the molecules. We find
the conditions for the existence of the polaritonic states, as well as
expressions for their energies. The polaritonic states are quite stable
against disorder. Analytical results are verified by comparison with
simulations. When $\epsilon_{a}$ is equal to that of the cavity state
$\epsilon_{c}$ (on resonance) the gap between the two polaritonic states is
found to increase with disorder, the increase being equal to
$2\frac{\sigma^{2}}{\sqrt{\mathscr{N}}|\tilde{V}|}$ where $\tilde{V}$ is the
coupling of a molecular excitation to the cavity state. An analytic expression
is found for the disorder induced width of the polaritonic peak. We results
for various densities of states, and the absorption spectrum. The dark states
that exist in the case $\sigma=0$ turn “grey” in presence of disorder with
their contribution to the absorption increasing with $\sigma$. We analyze the
effect of including lifetimes of the cavity and molecular states and find that
in the strong coupling regime, the width of the polaritonic peaks is dominated
by the lifetime effect and that disorder plays almost no role, if the Rabi
splitting is sufficiently large. We also consider the case where there is (a)
orientational disorder as well as (b) spatial variation of the cavity field
and find that they effectively amount to a renormlisation of the coupling.

##  I Introduction

In the recent past, states of molecules confined in a micro-cavity have
attracted much attention of chemists. The reason for this is that such
confinement can lead to very interesing modifications of transport properties
(energy and electron transfer reactions) and chemical reactivity. These happen
when an excited state of the molecule is almost resonant with a cavity mode.
In the strong coupling limit, modified chemical reactivity Thomas2016 ;
Thomas2019 ; Hirai2020 ; Peters2019 ; Munkhbat2018 , site selectivity
Thomas2019 , altered enzyme activities Vergauwe2019  and enhanced rate of
energy transfer Zhong2016 ; Zhong2017 ; Coles2014 ; Akulov2018  are some of
the processes that have been reported. These applications of light-matter
interactions in chemistry that have been reported confirms the presence of
hybrid light-matter states, for organic dye molecules at room temperature in
different states of matter confined in a cavity. The Ebbesen group has been
working on such systems and has reported Thomas2016 ; Thomas2019 ; Ebbesen2016
; Zhong2016 ; Zhong2017 ; Vergauwe2019 ; Hutchison2012 ; Orgiu2015  many
studies with various applications. Other research groups also have recently
published results Hirai2020 ; Peters2019 ; Pietron2019 ; Munkhbat2018 ;
Coles2014 ; Akulov2018 ; Casey2016  showing modified chemistry in such
cavities. One of the phenomena which is of great interest is the modification
in chemical reactivity, which can be achieved by making the cavity mode
resonate with different vibrational modes of the molecules. Many theoretical
models have been developed to understand this modified chemistry inside the
cavity. Transition state theory has been used to explain this phenomenon but
has not yet been able to provide a satisfactory explanation Galego2019 ;
Zhdanov2020 ; Jorge2020 . This shows the need of continuous exploration and
development of theories which could explain the experimental results. The
review articles Baranov2018 ; Flick2018 ; Herrera2020 ; Herrera2018 ;
Hertzog2019 ; Wang2021 ; Ribeiro2018 ; Feist2018  compile the past
developments, current status and unanswered questions in this rapidly
developing area of research.

When molecules are confined inside a cavity, there will be energy transfer
between excited states of the molecules and the cavity mode of the confined
electromagnetic field. If the rate of this energy exchange is fast enough
polaritonic modes are formed. For $N$ molecules confined to the cavity that
couple to the cavity mode through their transition dipoles, the coupling leads
to two polaritonic modes and $N-1$ other states. In the case where all the
molecules have identical excitation energies, these $(N-1)$ modes are
completely decoupled from the electromagnetic field and hence cannot be
excited by radiation. Consequently, they are referred to as the “dark modes”.
It is the formation of these mixed light+matter states that leads to the new
phenomena. The polaritonic states have been predicted to lead to enhanced
electron transfer rates for molecules confined in the cavity Herrera2016 . A
very interesting experimental demonstration is the energy transfer mediated by
polaritons Zhong2016 ; Coles2014 .

For molecules put inside a cavity there would be disorder in the energy levels
due to various reasons Scholes2020 . The causes may be inhomogeneities in the
material, different levels of aggregation, dispersal by solvent or solvent
fluctuations. It is well known that in general, transport properties are
strongly affected by the presence of disorder. The reason for this is that
disorder would in general cause the states to be localized (Anderson
localization) rather than delocalized. For models in which one has nearest
neighbour hopping, even a small disorder leads to localization in $1D$
(dimension) and $2D$. In $3D$ when the disorder exceeds a critical value the
system would undergo metal-insulator transition. Many aspects of disordered
systems confined to a cavity have been investigated, mostly by numerical
methods by Schachenmayer et. al. Botzung2020 . They consider the effect of the
cavity on the transport on disordered systems, with nearest neighbor hopping.
When the disorder is large enough, the states are localized and they
investigate the effect of the cavity on these localized states. They find that
the dark states are converted into “grey states”. Further, increasing the
coupling to the cavity increases delocalization. The states can be spread over
multiple sites, which are not neighbors in space, but are in energy. They find
that these “semi-localized” states obey a semi-Poisson statistics in the
energy level spacings. Also, the coupling through cavity modes is non-local
and hence the effects are not dependent on the dimensionality of the system.
The semi-localized states are found to be responsible for diffusive-like
dynamics. It has been suggested by Yuen-Zhou et al. Du2021 that the semi-
localized states contribute significantly to vibrational cooling, which leads
to significant enhancement of rate of electron transfer. In a new interesting
work, Celardo et al. Celardo2021  explore the long range hopping where
transport efficiency at first decreases, later increases and finally becomes
invariant as disorder increases.

Most of the investigations do not consider the position or orientational
dependence of the coupling of the molecule with the cavity mode, as they
assume the coupling to be a constant. In addition, the inhomogeneities of the
surrounding medium will affect the excitation energies, so that their values
would be distributed randomly. Such disorder has been investigated previously
in an interesting paper by Houdr$\acute{e}$ et. al. Houdre1996 . Our analysis
goes beyond theirs, in ways that are made clear in the Section IV of this
paper. In the present work we study the effect of these disorders on the
polaritonic and the dark states. Using the Tavis-Cummings Hamiltonian
Tavis1968 ; E1969  with added disorder, we analytically solve for the Green’s
function for the system in the limit $N>>1$. The density of states and
absorption cross-section also are analytically expressed in terms of the
Green’s function. In section III we consider the model where the coupling of
each molecule to the cavity state is considered to be the same, but their
excitation energies are randomly distributed. Having illustrated our approach
with this simple model, in section V we discuss the more general case in which
the coupling depends on the position as well as the orientation of the
transition dipole moment of the molecule with respect to the cavity field. To
see the effects of disorder clearly we have taken the lifetime of the excited
molecular state to be very long. We also discuss the inclusion of lifetime
effects. It is found that the polaritons have lifetimes which are weighted
averages of the lifetimes of the cavity and the molecular states. Further, on
resonance the disorder effects on the line-broadening are negligible if
$\Omega>>\frac{\gamma_{a}+\gamma_{c}}{2}>>\sigma$, where $\Omega$ is the rabi
splitting and $\gamma_{a}$ and $\gamma_{c}$ are the lifetimes of the molecular
and cavity states.

##  II The Hamiltonian and the Green’s operator

In the Tavis-Cummings Hamiltonian Tavis1968 ; E1969  with disorder in the
excitation energies, we take the excitation energy $\epsilon_{i}$ of the
$i^{th}$ molecule to be

$\epsilon_{i}=\epsilon_{a}+\xi_{i},$ (1)

where $\xi_{i}$ is a random variable.

The total Hamiltonian may then be written as

$\hat{H}=\epsilon_{c}|c\rangle\langle
c|+\sum_{i}^{N}(\epsilon_{a}+\xi_{i})|a_{i}\rangle\langle
a_{i}|+\sum_{i}^{N}(V_{i}|c\rangle\langle a_{i}|+V_{i}^{*}|a_{i}\rangle\langle
c|).$ (2)

$|c\rangle$ denotes the state where the cavity mode is excited to its first
excited state and it has an energy $\epsilon_{c}$. $|a_{i}\rangle$ denotes the
state in which the $i^{th}$ molecule is excited with all the other molecules
in their ground states. The coupling constant of the cavity mode to this
excitation is denoted by $V_{i}$ which is given by

$V_{i}=-\boldsymbol{E}(z_{i}).\boldsymbol{{\mu}}_{i}=-\sqrt{\frac{\epsilon_{c}}{2\epsilon_{o}{\cal
V}}}\mu_{eg}\sin(k_{c}z_{i})\cos(\theta_{i}),$ (3)

where ${\cal V}$ is the volume of the cavity and $\boldsymbol{{\mu}}_{i}$ is
transition dipole moment of the molecular excitation, having the magnitude
$\mu_{eg}$. The cavity field wave-vector and the permittivity of free space
are denoted as $k_{c}$ and $\epsilon_{o}$ respectively. The term
$\sin(k_{c}z_{i})$ accounts for the spatial variation of the cavity field and
$\theta_{i}$ is the angle between the transition dipole moment vector and the
electric field vector. We use the Green’s operator to analyze the problem
following the techniques of the classic paper of Anderson on the impurity
problem Anderson1961 . See also the book by Mattuck Mattuck2012 . The Green’s
operator for the Hamiltonian $\hat{H}$ is defined by
$\hat{G}(\omega)=(\omega-\hat{H})^{-1}$, where $\omega$ is a complex variable.
We shall calculate the matrix element $G_{ij}(\omega)$ of $\hat{G}(\omega)$
defined by

$G_{ij}(\omega)=\langle i|\frac{1}{\omega-\hat{H}}|j\rangle,$ (4)

where $|i\rangle$ and $|j\rangle$ are any two arbitrary states.
$G_{ij}(\omega)$ has the spectral representation

$G_{ij}(\omega)=\sum_{m}\frac{\langle
i|\varepsilon_{m}\rangle\langle\varepsilon_{m}|j\rangle}{\omega-\varepsilon_{m}},$
(5)

where $|\varepsilon_{m}\rangle$ is an eigenfunction of the Hamiltonian
operator having the eigenvalue $\varepsilon_{m}$, satisfying
$\hat{H}|\varepsilon_{m}\rangle=\varepsilon_{m}|\varepsilon_{m}\rangle.$ Note
that $m$ runs from $1$ to $(N+1)$. From this it is clear that any discrete
eigenvalue of $\hat{H}$ appears as a simple pole of $G_{ij}(\omega)$ and that
its residue at the pole at $\varepsilon_{m}$ is $\langle
i|\varepsilon_{m}\rangle\langle\varepsilon_{m}|j\rangle$, giving us an idea of
the amount of contribution that the $m^{th}$ eigenstate makes to $|i\rangle$
and $|j\rangle$. We now proceed to determine the matrix elements
$G_{ij}(\omega)$.

The Hamiltonian may be written in the matrix form as

$\underline{\underline{H}}=\left[\begin{array}[]{c|ccccccc}\epsilon_{c}&\;V_{1}&\;V_{2}&\;.&\;.&\;.&\;.&\;V_{N}\\\
\hline\cr V_{1}^{*}&\epsilon_{1}&0&.&.&.&.&0\\\
V_{2}^{*}&0&\epsilon_{2}&.&.&.&.&0\\\ .&.&.&.&.&.&.&.\\\
V_{N}^{*}&0&0&.&.&.&.&\epsilon_{N}\end{array}\right],$ (6)

which has the partitioned form

$\underline{\underline{H}}=\left[\begin{array}[]{c|c}\epsilon_{c}&\underline{V}\\\
\hline\cr\underline{V}^{\dagger}&\underline{\underline{\epsilon}}_{M}\end{array}\right].$
(7)

$M$ stands for molecules and $\underline{V}=[V_{1},V_{2},...,V_{N}]$ is a row
matrix having dimensions $1\times N$. $\underline{\underline{\epsilon}}_{M}$
is an $N\times N$ diagonal matrix of the excitation energies:

$\left(\underline{\underline{\epsilon}}_{M}\right)_{ij}=\epsilon_{i}\delta_{ij}.$
(8)

The matrix corresponding to $\hat{G}(\omega)$ operator is

$\underline{\underline{G}}=\left(\omega{\underline{\underline{I}}_{N+1}}-\underline{\underline{H}}\right)^{-1}.$
(9)

$\underline{\underline{I}}_{N+1}$ is the $(N+1)$ dimensional identity matrix.
$\underline{\underline{G}}$ too can be written in the same partitioned form as
in Eq. (7):

$\underline{\underline{G}}=\left[\begin{array}[]{c|c}G_{cc}&\underline{G}_{cM}\\\
\hline\cr\underline{G}_{cM}^{\dagger}&\underline{\underline{G}}_{M}\end{array}\right].$
(10)

Eq. (9) can be rearranged to get

$(\omega\underline{\underline{I}}_{N+1}-\underline{\underline{H}})\underline{\underline{G}}=\underline{\underline{I}}_{N+1}.$
(11)

Using Eq. (7), Eq. (9) and Eq. (10) we get

$\begin{pmatrix}\omega-\epsilon_{c}&-\underline{V}\\\
-\underline{V}^{\dagger}&\omega\underline{\underline{I}}-\underline{\underline{\epsilon}}_{M}\end{pmatrix}\begin{pmatrix}G_{cc}&\underline{G}_{cM}\\\
\underline{G}_{cM}^{\dagger}&\underline{\underline{G}}_{M}\end{pmatrix}=\begin{pmatrix}1&0\\\
0&\underline{\underline{I}}\end{pmatrix}.$ (12)

In the above $\underline{\underline{I}}$ is the $N\times N$ identity matrix.
Eq. (12) is equivalent to the following four equations:

$\displaystyle(\omega-\epsilon_{c})G_{cc}-\underline{V}\,\underline{G}_{cM}^{\dagger}$
$\displaystyle=$ $\displaystyle 1$ (13)
$\displaystyle(\omega-\epsilon_{c})\underline{G}_{cM}-\underline{V}\,\underline{\underline{G}}_{M}$
$\displaystyle=$ $\displaystyle 0$ (14)
$\displaystyle-\underline{V}^{\dagger}G_{cc}+(\omega\underline{\underline{I}}-\underline{\underline{\epsilon}}_{M})\underline{G}_{cM}^{\dagger}$
$\displaystyle=$ $\displaystyle 0$ (15)
$\displaystyle-\underline{V}^{\dagger}\underline{G}_{cM}+(\omega\underline{\underline{I}}-\underline{\underline{\epsilon}}_{M})\underline{\underline{G}}_{M}$
$\displaystyle=$ $\displaystyle\underline{\underline{I}}.$ (16)

Solving Eq. (15) for $\underline{G}_{cM}^{\dagger}$ gives

$\underline{G}_{cM}^{\dagger}=(\omega\underline{\underline{I}}-\underline{\underline{\epsilon}}_{M})^{-1}\underline{V}^{\dagger}G_{cc}.$
(17)

Using Eq. (17) in Eq. (13) and solving for $G_{cc}$ leads to

$G_{cc}=\left[\omega-\epsilon_{c}-\Sigma(\omega)\right]^{-1},$ (18)

with the self energy $\Sigma(\omega)$ defined by

$\Sigma(\omega)=\sum_{i}\frac{V_{i}^{*}V_{i}}{\omega-\epsilon_{i}}.$ (19)

We define the imaginary and real parts of the self energy by:

$\Sigma_{I}(\omega)=-\pi\sum_{i}|V_{i}|^{2}\delta(\omega-\epsilon_{i})$

and

$\Sigma_{R}(\omega)=-\frac{1}{\pi}{\cal
P}\int_{-\infty}^{\infty}d\epsilon\frac{\Sigma_{I}(\epsilon)}{\omega-\epsilon},$

where ${\cal P}$ stands for the principal value of the integral. It is useful
to introduce the density of states for the disordered molecular states by

$\rho_{d}(\xi)=\frac{1}{N}\sum_{i=1}^{N}\delta(\xi-\xi_{i}).$ (20)

In the following we first consider the simplest possible model where
$V_{i}=V,$ a constant. For this model, all the calculations are reported in
detail. After that we also discuss more general cases in section V.

##  III Model I: $V_{i}=V$

We note that the coupling $V\propto\frac{1}{\sqrt{{\cal{V}}}}$, where
${\cal{V}}$ is the volume of the cavity. We define $\mathscr{N}$ the number
density of molecules by $\mathscr{N}=N/{\cal{V}}$ is the number density of
molecules and volume independent coupling $\tilde{V}=V\sqrt{{\cal{V}}}$.

We use the notation $\omega^{+}=\omega+i\eta$ when $\omega$ is real and $\eta$
is taken to be an infinitestimal positive number. We then have

$\Sigma(\omega^{+})=\mathscr{N}\tilde{V}^{2}\int_{-\infty}^{\infty}d\xi\frac{\rho_{d}(\xi)}{\omega+i\eta-\epsilon_{a}-\xi}.$
(21)

Separating out the real and imaginary parts of $\Sigma(\omega^{+})$ gives

$\Sigma(\omega^{+})=\mathscr{N}\tilde{V}^{2}{\cal{P}}\int_{-\infty}^{\infty}d\xi\frac{\rho_{d}(\xi)}{\omega-\epsilon_{a}-\xi}-i\pi\mathscr{N}\tilde{V}^{2}\int_{-\infty}^{\infty}d\xi\rho_{d}(\xi)\delta(\omega-\epsilon_{a}-\xi),$
(22)

where ${\cal{P}}$ denotes the principal value of the integral.

In this simplest model with $V_{i}=V$, the coupling with the cavity field is
the same for all molecules. Further we shall take $\xi_{i}$ to be identically
distributed Gaussian random variables with the probability distribution

$P(\xi)=\frac{1}{\sqrt{2\pi}\sigma}e^{\frac{-\xi^{2}}{2\sigma^{2}}}.$ (23)

In principle, one has to average over different realisations of the random
quantities $\xi_{i}$. However when ${\mathscr{N}}$ is large, we expect
$\rho_{d}(\omega)\rightarrow P(\xi).$ Hence in this limit one expects _the
ensemble averaged value_ of self energies to be given by:

$\Sigma_{I}(\omega)=\frac{\mathscr{N}\tilde{V}^{2}\pi}{\sqrt{2\pi}\sigma}e^{-\frac{(\omega-\epsilon_{a})^{2}}{2\sigma^{2}}}$
(24)

and

$\Sigma_{R}(\omega)=\frac{\sqrt{2}\mathscr{N}\tilde{V}^{2}}{\sigma}DawsonF\left(\frac{\omega-\epsilon_{a}}{\sqrt{2}\sigma}\right).$
(25)

$DawsonF(x)$ is the Dawson Integral (also referred to as Dawson function)
DawsonMathworld , defined by

$DawsonF(x)=e^{-x^{2}}\int_{0}^{x}dy\;e^{y^{2}}.$ (26)

Using these results in Eq. (18) we get

$G_{cc}(\omega)=\left[\omega-\epsilon_{c}-\frac{\sqrt{2}\mathscr{N}\tilde{V}^{2}}{\sigma}DawsonF\left(\frac{\omega-\epsilon_{a}}{\sqrt{2}\sigma}\right)+\frac{i\mathscr{N}\tilde{V}^{2}\pi}{\sqrt{2\pi}\sigma}e^{\frac{-(\omega-\epsilon_{a})^{2}}{2\sigma^{2}}}\right]^{-1}.$
(27)

###  III.1 The Density of States

We define the density of states for the cavity state by

$\displaystyle\rho_{c}(\omega)$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{N+1}\left|\langle
c|m\rangle\right|^{2}\delta(\omega-\varepsilon_{m})$ (28) $\displaystyle=$
$\displaystyle-\frac{1}{\pi}Im\\{G_{cc}(\omega^{+})\\}.$

$\rho_{c}(\omega)$ gives us an idea of the amount of participation of the
cavity state in the $m^{th}$ eigenstate of the system.

Now we calculate the matrix elements of Green’s operator for the molecular
states. Eq. (15) can be solved for $\underline{G}_{cM}$ to get

$\underline{G}_{cM}=(\omega-\epsilon_{c})^{-1}\underline{V}\,\underline{\underline{G}}_{M}.$
(29)

Using this result in Eq. (16) gives

$\underline{\underline{G}}_{M}(\omega^{+})=\left[\underline{\underline{I}}\omega^{+}-\underline{\underline{\epsilon}}_{M}-\frac{\underline{V}^{\dagger}\underline{V}}{\omega^{+}-\epsilon_{c}}\right]^{-1}.$
(30)

The operator corresponding to the matrix
$\underline{V}^{\dagger}\underline{V}$ can be written in the form
$\mathscr{N}\tilde{V}^{2}\left|mol\right\rangle\left\langle mol\right|$ with

$\left|mol\right\rangle=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\left|i\right\rangle.$
(31)

Hence we can write the operator corresponding to the matrix
$\underline{\underline{G}}_{M}(\omega^{+})$ as

$\hat{G}_{M}(\omega^{+})=\left[\omega^{+}-\hat{\epsilon}_{M}-\frac{\mathscr{N}\tilde{V}^{2}\left|mol\right\rangle\left\langle
mol\right|}{\omega^{+}-\epsilon_{c}}\right]^{-1}.$ (32)

The operator $\hat{\epsilon}_{M}$ is defined by its matrix elements $\langle
i|\hat{\epsilon}_{M}|i^{\prime}\rangle=\epsilon_{i}\;\delta_{ii^{\prime}}$.
Using the operator identity
$(\hat{A}-\hat{B})^{-1}=\hat{A}^{-1}+\hat{A}^{-1}\hat{B}(\hat{A}-\hat{B})^{-1}$
valid for any two operators $\hat{A}$ and $\hat{B}$, we can write

$\displaystyle\hat{G}_{M}(\omega^{+})$
$\displaystyle=\left[\omega^{+}-\hat{\epsilon}_{M}\right]^{-1}+\left[\omega^{+}-\hat{\epsilon}_{M}\right]^{-1}\frac{\mathscr{N}\tilde{V}^{2}\left|mol\right\rangle\left\langle
mol\right|}{\omega^{+}-\epsilon_{c}}\hat{G}_{M}(\omega^{+}).$ (33)

From this we find $\langle mol|\hat{G}_{M}(\omega^{+})$ to obey the equation

$\displaystyle\langle mol|\hat{G}_{M}(\omega^{+})$ $\displaystyle=\langle
mol|\left[\omega^{+}-\hat{\epsilon}_{M}\right]^{-1}+\frac{\Sigma(\omega^{+})}{(\omega^{+}-\epsilon_{c})}\langle
mol|\hat{G}_{M}(\omega^{+}).$ (34)

Solving for $\langle mol|\hat{G}_{M}(\omega^{+})$ gives

$\displaystyle\langle mol|\hat{G}_{M}(\omega^{+})=$
$\displaystyle(\omega^{+}-\epsilon_{c})G_{cc}(\omega^{+})\langle
mol|\left[\omega^{+}-\hat{\epsilon}_{M}\right]^{-1}.$

Using this back in Eq. (33) we get

$\hat{G}_{M}(\omega^{+})=\left[\omega^{+}-\hat{\epsilon}_{M}\right]^{-1}+\left[\omega^{+}-\hat{\epsilon}_{M}\right]^{-1}\mathscr{N}\tilde{V}^{2}\left|mol\right\rangle
G_{cc}(\omega^{+})\langle
mol|\left[\omega^{+}-\hat{\epsilon}_{M}\right]^{-1}.$ (35)

Using the above, the matrix element $G_{mol,mol}(\omega^{+})=\left\langle
mol\right|\hat{G}_{M}(\omega^{+})\left|mol\right\rangle$ can be calculated to
be:

$\displaystyle G_{mol,mol}(\omega^{+})$ $\displaystyle=$
$\displaystyle\left\\{1-\frac{\mathscr{N}\tilde{V}^{2}}{\omega^{+}-\epsilon_{c}}\left\langle
mol\right|(\omega^{+}-\hat{\epsilon}_{M})^{-1}\left|mol\right\rangle\right\\}^{-1}\left\langle
mol\right|(\omega^{+}-\hat{\epsilon}_{M})^{-1}\left|mol\right\rangle$ (36)
$\displaystyle=$
$\displaystyle(\omega^{+}-\epsilon_{c})G_{cc}(\omega^{+})\frac{1}{\mathscr{N}\tilde{V}^{2}}\Sigma(\omega^{+}).$
(37)

From $G_{mol,mol}(\omega^{+})$ we can calculate the density of states of
$|mol\rangle:$

$\rho_{mol}(\omega)=-\frac{1}{\pi}Im\\{G_{mol,mol}(\omega^{+})\\}.$ (38)

The total density of states may be defined by
$\rho_{T}(\omega)=\sum_{m=1}^{N+1}\delta(\omega-\varepsilon_{m})$ which may be
written as
$\rho_{T}(\omega)=-\frac{1}{\pi}Im\left(Tr{\hat{G}(\omega^{+})}\right)$. Using
the complete set $|c\rangle,|1\rangle,|2\rangle,....|N\rangle$ to calculate
the trace leads to

$\displaystyle\rho_{T}(\omega)$ $\displaystyle=$
$\displaystyle\rho_{c}(\omega)-\sum_{i=1}^{N}\frac{1}{\pi}Im\\{G_{ii}(\omega^{+})\\},$
(39)

with $G_{ii}(\omega^{+})=\langle i|\hat{G}(\omega^{+})|i\rangle$.
$G_{ii}(\omega^{+})$ for $i\neq c$ may be easily evaluated using Eq. (35), to
get

$G_{ii}(\omega^{+})=(\omega^{+}-\epsilon_{i})^{-1}+V^{2}(\omega^{+}-\epsilon_{i})^{-2}G_{cc}(\omega^{+}),$

so that

$\sum_{i=1}^{N}G_{ii}(\omega^{+})=\sum_{i=1}^{N}(\omega^{+}-\epsilon_{i})^{-1}+V^{2}G_{cc}(\omega^{+})\sum_{i=1}^{N}(\omega^{+}-\epsilon_{i})^{-2}.$
(40)

A change in the total density of states for the molecules alone may be defined
by

$\Delta\rho_{M}(\omega)=-\frac{1}{\pi}Im\left(\sum_{i=1}^{N}G_{ii}(\omega^{+})\right)-\rho_{M}^{0}(\omega),$

where $\rho_{M}^{0}(\omega)$ is the density of molecular states in the case
where $\tilde{V}=0$. The above may be calculated using Eq. (40) to be

$\Delta\rho_{M}(\omega)=\frac{1}{\pi}Im\left(\frac{\partial\Sigma(\omega^{+})}{\partial\omega}G_{cc}(\omega^{+})\right)$
(41)

so that the change in the total density of states, due to the interaction of
the cavity mode with the molecules is

$\Delta\rho_{T}(\omega)=-\frac{1}{\pi}Im\\{G_{cc}(\omega^{+})\\}-\frac{1}{\pi}Im\\{\frac{\partial\Sigma(\omega^{+})}{\partial\omega}G_{cc}(\omega^{+})\\}.$
(42)

###  III.2 The Absorption Spectrum

The Hamiltonian for the interaction of the system with radiation is

$\hat{H}_{int}(t)=-\sum_{i}\text{$\boldsymbol{E}_{i}$}(t).\boldsymbol{\hat{\mu}}_{i}.$

The $\boldsymbol{\hat{\mu}}_{i}$ operator can be written in terms in terms of
transition dipole moment as:

$\boldsymbol{\hat{\mu}}_{i}=\boldsymbol{\mu_{i}}|g\rangle\langle
i|+\boldsymbol{\mu_{i}^{*}}|i\rangle\langle g|$

where $\boldsymbol{\mu_{i}}=\langle i|\boldsymbol{\hat{\mu}}_{i}|g\rangle$. We
can rewrite the Hamiltonian for interaction as:

$\hat{H}_{int}(t)=-\sum_{i}\boldsymbol{E}_{i}(t).\left(\boldsymbol{\mu_{i}}|g\rangle\langle
i|+\boldsymbol{\mu_{i}^{*}}|i\rangle\langle g|\right).$ (43)

$|g\rangle$ indicates the ground state in which all the molecules are in their
ground states. $|i\rangle$ is the site in which the $i^{th}$ molecule is
excited. $\boldsymbol{E}_{i}(t)$ is the harmonic electric field of the
frequency $\omega$ that the $i^{th}$ molecule experiences. In the spirit of
the model that assumes $V_{i}=V$ a constant, we first take the electric field
to be the same for all the molecules - i.e., neglect its space dependence so
that all the molecules experience the same electric field. Further, for
simplicity, we take all of them to be oriented in the direction of the
electric field. Under these assumptions,

$\hat{H}_{int}(t)=-\sum_{i}E(t)\left(\mu_{eg}|g\rangle\langle
i|+\mu_{ge}^{*}|i\rangle\langle g|\right)$

where $\mu_{eg}=|\boldsymbol{\mu_{i}}|$ and $E(t)$ is the electric field in
the direction of the transition dipole of $i^{th}$ molecule. The absorption
cross section $\alpha(\omega)$ for radiation of frequency $\omega$ is then

$\alpha(\omega)=\frac{\pi\omega}{\epsilon_{o}c\hbar}\;\;\sum_{f}\left|\left\langle
f\left|\hat{H}_{int}(0)\right|g\right\rangle\right|^{2}\delta(E_{f}-E_{g}-\hbar\omega).$

In the above, $E_{g}$ is the energy of the ground state. $|f\rangle$ denotes a
possible final state, having energy $E_{f}$. Following reference schatz2012
we write this as the Fourier transform of a correlation function:

$\alpha(\omega)=\frac{\omega}{2\epsilon_{o}c\hbar}\int_{-\infty}^{\infty}dte^{i\omega
t}\left\langle g\left|\hat{H}_{int}(t)\hat{H}_{int}\right|g\right\rangle,$

where
$\hat{H}_{int}(t)=e^{i\hat{H}t/\hbar}\hat{H}_{int}(0)e^{-i\hat{H}t/\hbar}$. On
using the expression for $H_{int}$ and simplifying we get

$\displaystyle\alpha(\omega)$
$\displaystyle=\frac{\omega}{2\epsilon_{o}c\hbar}\left|\boldsymbol{{\mu}}\right|^{2}\int_{-\infty}^{\infty}dte^{i\omega
t}\sum_{i,j}\left\langle i\left|e^{i\hat{H}t/\hbar}\right|j\right\rangle$
$\displaystyle=-\frac{\omega}{\epsilon_{o}c\hbar}N\left|\boldsymbol{{\mu}}\right|^{2}\;Im\left[G_{mol,mol}(\omega)\right].$
(44)

##  IV Results

###  IV.1 Existence of polaritonic states

We now analyze $G_{cc}(\omega)$ to see the effect of disorder on the
polaritonic states. At small disorder one expects the polaritonic states to
exist. However, as the disorder increases, one suspects that the polaritonic
states may disappear. For the sake of simplicity we consider only the case
where $\epsilon_{c}=\epsilon_{a}$, and it is easy to extend the analysis to
more general cases. In the following we make a quantitative analysis using the
figures 1 and 2, which differ in their values of ${\mathscr{N}}$. Both show
the plots of $\omega-\epsilon_{c}$, $\Sigma_{R}(\omega-\epsilon_{c})$ and
$\Sigma_{I}(\omega-\epsilon_{c})$ against $\omega-\epsilon_{c}$. A pole of
$G_{cc}(\omega)$ occurs at a point where the plot of $\omega-\epsilon_{c}$
intersects with that of $\Sigma_{R}(\omega-\epsilon_{c})$. Fig. 1 is for large
values of ${\mathscr{N}}$ and in this case there are three points where this
happens. The intersection at the origin has a large value for the imaginary
part $\Sigma_{I}(\omega-\epsilon_{c})$ implying that this a virtual state. The
other two intersections have only a small value of
$\Sigma_{I}(\omega-\epsilon_{c})$ indicating these are the long lived
polaritonic states. This figure is to be compared with Fig. 2 which is for
smaller values of ${\mathscr{N}}$ or larger value of $\sigma$. In this case
the polaritonic states do not exist. Only the virtual state is there. The plot
of $\Sigma_{R}(\omega-\epsilon_{c})$ has a slope of
$\frac{\mathscr{N}\tilde{V}^{2}}{\sigma^{2}}$ at $\omega=\epsilon_{c}$. Hence
the two polaritonic poles can exist only if
$\frac{\mathscr{N}\tilde{V}^{2}}{\sigma^{2}}>1$. For suitable values of the
parameters, it is also possible that only one polaritonic state exists.

Figure 1: This figure takes $\epsilon_{a}$ to be equal to $\epsilon_{c}$. The
plots of $\omega-\epsilon_{c},\;\Sigma_{R}(\omega-\epsilon_{c})$ and
$\Sigma_{I}(\omega-\epsilon_{c})$ against $(\omega-\epsilon_{c}).$ The
darkened small circles indicate points where $\omega-\epsilon_{c}$ intersects
with $\Sigma_{R}(\omega-\epsilon_{c})$ leading to the formation of two
discrete poles. These correspond to the polaritonic states. The intersection
at the origin indicates a virtual state as the imaginary part of the self
energy $\Sigma_{I}(0)$ is large at this point. Figure 2: This figure takes
$\epsilon_{a}$ to be equal to $\epsilon_{c}$. The plots of
$\omega-\epsilon_{c}$, $\Sigma_{R}(\omega-\epsilon_{c})$ and
$\Sigma_{I}(\omega-\epsilon_{c})$ against $(\omega-\epsilon_{c})$ for
increased disorder. There is only one point where $\omega-\epsilon_{c}$
intersects with $\Sigma_{R}(\omega-\epsilon_{c})$ and that is just the virtual
state.

We now make an analysis that is valid for arbitrary values of $\epsilon_{a}$
and $\epsilon_{c}$. We consider the case where disorder is not large so that
the situation in Fig. 1 occurs and there are two poles for $G_{cc}(\omega)$.
These poles occur for large values of $|\omega-\epsilon_{c}|$ and
$|\omega-\epsilon_{a}|$ . The Dawson Function has the asymptotic form

$DawsonF(x)=\frac{1}{2x}+\frac{1}{4x^{3}}+\frac{3}{8x^{5}}+O\left(\frac{1}{x^{7}}\right)\;\;\;\text{for
$|x|$ large.}$

Using this to the lowest order (i.e. only the $1/x$ term is retained), we find
that for $\left|\frac{\omega-\epsilon_{a}}{\sqrt{2}\sigma}\right|>>1$,
$\frac{\sqrt{2}\mathscr{N}\tilde{V}^{2}}{\sigma}DawsonF\left(\frac{\omega-\epsilon_{a}}{\sqrt{2}\sigma}\right)\simeq\frac{\mathscr{N}\tilde{V}^{2}}{(\omega-\epsilon_{a})}$
and hence

$G_{cc}(\omega)=\left(\omega-\epsilon_{c}-\frac{\mathscr{N}\tilde{V}^{2}}{(\omega-\epsilon_{a})}+i\frac{\mathscr{N}\tilde{V}^{2}\pi}{\sqrt{2\pi}\sigma}e^{\frac{-(\omega-\epsilon_{a})^{2}}{2\sigma^{2}}}\right)^{-1}.$

As $(\omega-\epsilon_{a})/(\sqrt{2}\sigma)>>1,$ the imaginary part in the
above is very small and may be neglected for a first analysis. Then
$G_{cc}(\omega)$ has poles at $\epsilon_{\pm}$ given by

$\epsilon_{\pm}=\frac{\left(\epsilon_{a}+\epsilon_{c}\right)}{2}\pm\sqrt{\mathscr{N}\tilde{V}^{2}+\left(\frac{\epsilon_{a}-\epsilon_{c}}{2}\right)^{2}}.$
(45)

These are the two polaritonic states $|+\rangle$ and $|-\rangle$. The
similarity of the above to simple molecular orbital theory coulson1979  is to
be noted. The residues of these at the poles are easily found to be

$|\langle
c|\pm\rangle|^{2}=\left[1+\frac{\mathscr{N}\tilde{V}^{2}}{(\epsilon_{\pm}-\epsilon_{a})^{2}}\right]^{-1}.$

For $\sqrt{\mathscr{N}}|\tilde{V}|\gg|\epsilon_{c}-\epsilon_{a}|$ these reduce
to

$\epsilon_{\pm}\simeq\frac{\left(\epsilon_{a}+\epsilon_{c}\right)}{2}\pm|\tilde{V}|\sqrt{\mathscr{N}}$

and

$|\langle c|\pm\rangle|^{2}\simeq\frac{1}{2}.$

At these poles, the imaginary part $\left[G_{cc}(\omega)\right]^{-1}$ can be
easily estimated to be $\approx
i\frac{\mathscr{N}\tilde{V}^{2}\pi}{\sqrt{2\pi}\sigma}e^{-\frac{\mathscr{N}\tilde{V}^{2}}{2\sigma^{2}}}$.
This gives the width of the polaritonic peaks. Note that the width of the peak
is decreasing exponentialy with ${\mathscr{N}}$ and hence its height would
increase exponentially.

From Eq. (45) we see the independence of $\epsilon_{\pm}$ upto first order in
$\sigma$. To see the effect of disorder, it is necessary to do an analysis up
to second order in $\sigma$. Under the assumption that the imaginary part is
small and keeping the next two terms too in the asymptotic expansion of
$DawsonF(x)$ we can find solutions for $\epsilon_{\pm}$ up to order
$\sigma^{2}$ using MATHEMATICA to be

$\epsilon_{s}=\frac{\left(\epsilon_{a}+\epsilon_{c}\right)}{2}+s\sqrt{\mathscr{N}\tilde{V}^{2}+\left(\frac{\epsilon_{a}-\epsilon_{c}}{2}\right)^{2}}+s\frac{\sigma^{2}\left((\frac{\epsilon_{a}-\epsilon_{c}}{2})\left(\frac{\epsilon_{a}-\epsilon_{c}}{2}+s\sqrt{\mathscr{N}\tilde{V}^{2}+(\frac{\epsilon_{a}-\epsilon_{c}}{2})^{2}}\right)+\mathscr{N}\tilde{V}^{2}\right)}{\mathscr{N}\tilde{V}^{2}\sqrt{\mathscr{N}\tilde{V}^{2}+(\frac{\epsilon_{a}-\epsilon_{c}}{2})^{2}}},$
(46)

where $s=\pm$. For $\epsilon_{c}=\epsilon_{a}$ this reduces to

$\epsilon_{\pm}=\epsilon_{c}\pm\left(\sqrt{\mathscr{N}}|V|+\frac{\sigma^{2}}{\sqrt{\mathscr{N}}|V|}\right).$
(47)

Thus when $\epsilon_{a}=\epsilon_{c}$, the gap between the two polaritonic
states is increased by the disorder, a conclusion which is valid to order
$\sigma^{2}$. Further, as the number of molecules increase, the effects of
disorder on the polaritonic states decrease. These results indicate that the
polaritons are quite stable even in the presence of disorder, provided the
number density of molecules is large enough, i.e., when
$\sqrt{\mathscr{N}}|\tilde{V}|\gg\sigma$, a condition that is satisfied
easily. It may be noted that an analysis valid up to first order in $\xi$ was
performed in the paper by Houdr$\acute{e}$ Houdre1996 , who concluded that
disorder has no effect on the separation between the polaritonic states. In
our notation, their analysis corresponds to doing calculations correct up to
first order in $\sigma$. From the above it is clear that at this order in
$\sigma$, disorder has no effect on the separation between the polaritonic
peaks, in agreement with them.

####  IV.1.1 Testing analytical results using simulations

We perform numerical calculations to determine the density of states using the
equations given in previous sections. For this we generated ${\mathscr{N}}$
values of $\xi_{i}$ using a Gaussian distribution with zero mean and standard
deviation $\sigma$. In principle, we should solve for the roots of
$(\omega-\epsilon_{c}-\Sigma_{R}(\omega))$ and then construct the plot. We can
easily bypass this by assuming that $\omega$ has a small imaginary part and
then calculating $-\frac{1}{\pi}Im[G_{cc}(\omega)].$ In our calculations, we
have taken the imaginary part to be equal to $0.001eV$. The resulting plot is
not very sensitive to the value of this imaginary part. All the other
quantities of interest are evaluated in a similar fashion. For each quantity
we average over 3000 realizations. In Fig. 3, 4 and 5 one can observe the good
agreement between the simulations and the theoretical results even for
$\mathscr{N}$ as small as $\mathscr{N}=1500\,m^{-3}$. Fig. 3 shows the close
agreement between theoretically calculated $\rho_{c}(\omega)$ with the one
obtained from simulations. In Fig. 4 we show plots of $\rho_{T}(\omega)$
against $\omega$ for two values of $\sigma$. Because of the division by
$\mathscr{N}$ the polaritonic peaks have a small height and are not clearly
visible.

Figure 3: Plot of $\rho_{c}(\omega)$ against $\omega$ for
$\mathscr{N}=1500\,m^{-3}$. Following values were used:
$\epsilon_{c}=2.0\,eV$, $\epsilon_{a}=2.0\,eV$, $\tilde{V}=3.56\times
10^{-3}\,eVm^{3/2}$ and $\eta=0.001\,eV$. The points represents the
simulations and the line plot is for the analytical result. Figure 4: Plot of
$\rho_{T}(\omega)/\mathscr{N}$ against $\omega$ for
$\mathscr{N}=1500\,m^{-3}$. Following values were used $\epsilon_{c}=2.0\,eV$,
$\epsilon_{a}=2.0\,eV$, $\tilde{V}=3.56\times 10^{-3}\,eVm^{3/2}$ and
$\eta=0.001\,eV$. The points represent the simulations and the line plot is
for the theoretical result. The polaritonic peaks barely visible in the plot,
due to the division by $\mathscr{N}$.

Fig. 5 shows the comparison of the analytical and simulated results for
$-\frac{1}{\pi}Im[G_{mol,mol}(\omega)]$. The agreement between the two is
good.

Figure 5: Plot of $-\frac{1}{\pi}Im[G_{mol,mol}(\omega)]$ against $\omega$ for
$\mathscr{N}=1500\,m^{-3}$. Following values were used:
$\epsilon_{c}=2.0\,eV$, $\epsilon_{a}=2.0\,eV$, $\tilde{V}=3.56\times
10^{-3}\,eVm^{3/2}$ and $\eta=0.001\,eV$. The solid points represent the
simulations and the line plot is for the analytical result.

Having illustrated the results using smaller values of $\mathscr{N}$, we now
extend our results to experimentally viable values, which are too large to be
simulated. As our analysis becomes better and better with increase in
$\mathscr{N}$ one expects the analytical results to be accurate for such large
values of $\mathscr{N}$. In Fig. 6 we plot $\rho_{c}(\omega)$ against $\omega$
for $\sigma=0.01\epsilon_{a},0.03\epsilon_{a}$ and $0.05\epsilon_{a}$. The
peak heights for the plot corresponding to $\sigma=0.01\,\epsilon_{a}$ value
is of the order of $10^{9}$ and hence the top portion has been cutoff. The
plots show that the polaritonic states are very much affected by the disorder,
as seen from the analytical estimate of their widths which is of the order of
$\approx\frac{\mathscr{N}\tilde{V}^{2}\pi}{\sqrt{2\pi}\sigma}e^{-\frac{\mathscr{N}\tilde{V}^{2}}{2\sigma^{2}}}$.

Figure 6: Plot for $\rho_{c}(\omega)$ v/s $\omega$ for $\mathscr{N}=1.15\times
10^{25}\,m^{-3}$ and for different values of $\sigma$. Following values were
used: $\epsilon_{c}=2.1\,eV$, $\epsilon_{a}=\epsilon_{c}$,
$\tilde{V}=4.06\times 10^{-14}\,eVm^{3/2}$.

In Fig. 7 we plot $\rho_{T}(\omega)/\mathscr{N}$ vs $\omega$ for
experimentally viable values for different values of $\sigma$.

Figure 7: Plot for $\rho_{T}(\omega)/\mathscr{N}$ v/s $\omega$ for
$\mathscr{N}=1.15\times 10^{25}\,m^{-3}$ and for different values of $\sigma$.
Following values were used: $\epsilon_{c}=2.1\,eV$,
$\epsilon_{a}=\epsilon_{c}$, $\tilde{V}=4.06\times 10^{-14}\,eVm^{3/2}$.

One does not see the polaritonic peaks in this plot because the height and
width of those two peaks are very small for all values of $\sigma$.

####  IV.1.2 Change in the density of states

Using the expression for the change in density of states given in Eq. (41) we
plot $\Delta\rho_{M}(\omega)$ in Fig. 8. Integrating the area under the curve
from $\omega=2.0\;eV$ to $\omega=2.2\;eV$ results in the value of $-1$. This
is what one would expect, because effectively out of the $N$ molecular states
$1$ is used up for the formation of the two polaritonic states. The plot has a
minimum at the middle. This is because it is the molecular states that are
near resonance with the cavity state that contribute the most to the
polaritonic states and hence the depletion of density of $\rho_{M}(\omega)$ is
maximum in the middle. In this figure, for $\sigma=0.01\,\epsilon_{a}$ the
polaritonic peaks are very high (height $\approx 10^{9}$) peak with a very
narrow (width $\approx 10^{-10}$).

Figure 8: Plot of change in molecular density of states
$\Delta\rho_{M}(\omega)$ vs $\omega$ for $\mathscr{N}=1.16\times
10^{25}\,m^{-3}$ and varied values of $\sigma$. Following values were used:
$\epsilon_{c}=2.1\,eV$, $\epsilon_{a}=\epsilon_{c}$, $\tilde{V}=4.06\times
10^{-14}\,eVm^{3/2}$.

####  IV.1.3 Absorption Cross Section

In Fig. 9 we plot the absorption cross-section given by Eq. (59) for
increasing values of $\sigma$. The shift of peak position with increasing
values of $\sigma$ is clearly seen. In addition the peaks get broader as
$\sigma$ increases. This indicates that with increase in disorder, the dark
states are turning “turning grey”. The dependence of width on $\sigma$ is
highly non-linear. The approximate relation between $\sigma$ and the width,
determined from the imaginary part of $[G_{cc}(\omega)]^{-1}$ at the poles, is
$\approx\frac{\mathscr{N}\tilde{V}^{2}\pi}{\sqrt{2\pi}\sigma}e^{-\frac{\mathscr{N}\tilde{V}^{2}}{2\sigma^{2}}}$.
When $\sigma<<2\sqrt{\mathscr{N}}\tilde{V}$ with
$\sigma=0.01\epsilon_{a}=0.02\,eV$ and
$2\sqrt{\mathscr{N}}\tilde{V}=0.275\,eV$, contribution of the dark states to
the absorption spectra is miniscule. The absorption is mostly from the two
polaritonic states and amounts to $99.86\%$ of the total and the absorption
due to the grey states amounts to $0.14\%$ when $\sigma=0.01\epsilon_{a}$. For
$\sigma=0.01\,\epsilon_{a}$, the FWHM for the peaks corresponding to the two
polaritonic states is of the order of $10^{-10}\,eV$ and peak height is of the
order of $10^{9}$. As we increase the value of $\sigma$ from $\sigma=0.01\,eV$
to $\sigma=0.05\epsilon_{a}=0.10eV$ the contribution of the dark states to the
absorption spectra steadily increases.

Figure 9: Plot of $\alpha$ against $\omega$ for $\mathscr{N}=1.15\times
10^{25}\,m^{-3}$ with varied values of $\sigma$. Following values were used:
$\epsilon_{c}=2.1\,eV$, $\epsilon_{a}=\epsilon_{c}$, $\tilde{V}=4.06\times
10^{-14}\,eVm^{3/2}$ and $|\boldsymbol{\mu}|=10\,D$. For the sake of clarity
of the figure, the y-axis is cutoff at the value of $20$. The peak height for
$\sigma=0.01\epsilon_{a}$ is of the order of $10^{10}$.

####  IV.1.4 Lifetime effects

Untill now we considered only the disorder as the source of line broadening.
In this section we add the effect of homogeneous broadening for the molecular
states $(\gamma_{a})$ and the cavity mode linewidth $(\gamma_{c})$. For this,
we replace $\epsilon_{c}$ by $\epsilon_{c}-i\gamma_{c}$ and $\epsilon_{a}$ by
$\epsilon_{a}-i\gamma_{a}$. Then the poles of the $G_{cc}(\omega)$ are at

$\epsilon_{\pm}=\frac{1}{2}\left(\epsilon_{a}+\epsilon_{c}-i(\gamma_{a}+\gamma_{c})\pm\sqrt{(-\epsilon_{a}+\epsilon_{c}+i(\gamma_{a}-\gamma_{c}))^{2}+4\mathscr{N}\tilde{V}^{2}}\right).$
(48)

The above expression shows that the polaritons have lifetimes which are
appropriately weighted averages of the lifetimes of the cavity and the
molecular states, and that for large $\mathscr{N}$, disorder effects on the
line-broadening are negligible.

For typical value of $\gamma_{a}=10\,ps$ $(0.3\,meV)$ and $\gamma_{c}=40\,fs$
$(0.1\,eV)$, we study the effect on the lineshape with increasing Rabi-
splitting ($\Omega=2\sqrt{\mathscr{N}}\tilde{V}$), shown in Fig. 10. In the
previous section we stated that as $\mathscr{N}$ increases the width of the
polaritonic peaks exponentially decreases
$\left(\frac{\mathscr{N}\tilde{V}^{2}\pi}{\sqrt{2\pi}\sigma}e^{-\frac{\mathscr{N}\tilde{V}^{2}}{2\sigma^{2}}}\right)$
with it, therefore as the value of $\Omega$ increases, the effects of disorder
vanishes exponentially and the contribution of the disorder to the linewidth
goes down. In figures 11 and 12 we plot the absorption cross-section with
varied values of $\sigma$ keeping $\Omega$ constant. We can observe from Fig.
11 that as $\sigma$ increases the peak position shifts and the lineshape is
modified. But as shown in Fig. 12 for a higher value of $\Omega$ there is
harldy any effect of $\sigma$ on the lineshape. Therefore when
$\epsilon_{a}=\epsilon_{c}$, under the condition
$\Omega>>\frac{\gamma_{a}+\gamma_{c}}{2}>>\sigma$ the effect of disorder is
negligible. This result echoes the findings of Houdré Houdre1996 .

Figure 10: Plot for absorption cross-section $\alpha$ v/s $\omega$ for varied
value of $\Omega$. Following values were used: $\epsilon_{c}=2.1\,eV$,
$\epsilon_{a}=\epsilon_{c}$, $\tilde{V}=4.06\times 10^{-14}\,eVm^{3/2}$,
$\sigma=0.01\epsilon_{a}$, $\gamma_{a}=0.3\,meV$ and $\gamma_{c}=0.1\,eV$.
Figure 11: Plot for absorption cross-section $\alpha$ v/s $\omega$ for varied
value of $\sigma$. Following values were used: $\epsilon_{c}=2.1\,eV$,
$\epsilon_{a}=\epsilon_{c}$, $\tilde{V}=4.06\times 10^{-14}\,eVm^{3/2}$,
$\Omega=256\,meV$, $\gamma_{a}=0.3\,meV$ and $\gamma_{c}=0.1\,eV$. Figure 12:
Plot for absorption cross-section $\alpha$ v/s $\omega$ for varied value of
$\sigma$. Following values were used: $\epsilon_{c}=2.1\,eV$,
$\epsilon_{a}=\epsilon_{c}$, $\tilde{V}=4.06\times 10^{-14}\,eVm^{3/2}$,
$\Omega=812\,meV$, $\gamma_{a}=0.3\,meV$ and $\gamma_{c}=0.1\,eV$.

##  V Model II: $V_{i}=V(\theta_{i},z_{i})$

The next situation we consider has orientational and position dependent
coupling i.e., $V_{i}=V\cos(\theta_{i})\sin(kz_{i}/L)$. $\theta_{i}$ is the
angle between the electric field and transition dipole moment of the $i^{th}$
molecule. The $\sin(kz_{i}/L)$ term accounts for the spatial variation of the
electric field in the cavity of length $L$. $\xi_{i}$ is assumed to be
identically distributed Gaussian random variables with the probability
distribution same as the one considered in Model I. As a result, we get

$\displaystyle\langle\Sigma(\omega)\rangle=\frac{V^{2}}{4\pi
L}\int_{0}^{\pi}\int_{0}^{L}\int_{0}^{2\pi}d\theta_{i}d\phi_{i}dz_{i}\sum_{i}\frac{\sin(\theta_{i})\cos(\theta_{i})^{2}\sin(kz_{i})^{2}}{\omega-\epsilon_{i}}.$
(49)

Averaging over the angles and length will give us:

$\displaystyle\langle\Sigma(\omega)\rangle=\frac{V^{2}}{6}\sum_{i}\frac{1}{\omega-\epsilon_{i}}.$
(50)

This means that effectively the only result of orientational and positional
averaging is to reduce the coupling per molecule by a factor of $\sqrt{6}$.

###  V.1 The Absorption Spectrum for Model II

Using the general interaction Hamiltonian given by Eq. (43), the absorption
cross-section can be written as:

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\frac{\omega}{2\epsilon_{o}c\hbar}\int_{-\infty}^{\infty}e^{i\omega
t}\left\langle g\left|\hat{H}_{int}(t)\hat{H}_{int}\right|g\right\rangle dt$
(51)

For this case the definition of $|mol\rangle$ needs to be modified. Now it is
given by

$|mol\rangle=F\sum_{i}\boldsymbol{\hat{\mu}}_{i}.\boldsymbol{E}_{i}|g\rangle,$
(52)

where $\boldsymbol{\hat{\mu}}_{i}=\boldsymbol{\mu_{i}}|g\rangle\langle
e_{i}|+\boldsymbol{\mu_{i}^{*}}|e_{i}\rangle\langle g|$ and $F$ is the
normalization constant given by:

$F^{2}\sum_{i}\left|\boldsymbol{{\mu}}_{i}.\boldsymbol{E}_{i}\right|^{2}=1$
(53)

Using this definition for $|mol\rangle$ we can rewrite Eq. (51) as

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\frac{\omega}{2\epsilon_{o}c\hbar
F^{2}}\int_{-\infty}^{\infty}e^{i\omega t}\langle
mol\left|e^{-iHt/\hbar}\right|mol\rangle e^{iE_{g}t/\hbar}dt$ (54)
$\displaystyle=$ $\displaystyle-\frac{\omega}{\epsilon_{o}c\hbar
F^{2}}Im\left\\{G_{mol,mol}(\omega)\right\\}.$ (55)

Since $F$ contains the information over orientation and spatial distribution
we substitute it back and average over length and angle,

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle-\frac{\omega}{\epsilon_{o}c\hbar}\sum_{i}\left|\boldsymbol{{\mu}}_{i}.\boldsymbol{E}_{i}\right|^{2}Im\left\\{G_{mol,mol}(\omega)\right\\}$
(56)

To determine $G_{mol,mol}(\omega)$ we re-write Eq. (30)

$\underline{\underline{G}}_{M}(\omega^{+})=\left[\omega^{+}-\underline{\underline{\epsilon}}_{M}-\frac{\underline{V}^{\dagger}\underline{V}}{\omega^{+}-\epsilon_{c}}\right]^{-1}.$
(57)

Using Eq. (52) we can write the operator corresonding to
$\underline{V}^{\dagger}\underline{V}$ in the form $F^{-2}|mol\rangle\langle
mol|$, this gives us,

$\hat{G}_{M}(\omega^{+})=\left[\omega^{+}-\hat{\epsilon}_{M}-\frac{F^{-2}\left|mol\right\rangle\left\langle
mol\right|}{\omega^{+}-\epsilon_{c}}\right]^{-1}.$ (58)

Similar to the procedure for model I we can follow the steps to get:

$\displaystyle{G}_{mol,mol}(\omega^{+})=$
$\displaystyle\left\\{1-\frac{F^{-2}\langle
mol|(\left[\omega^{+}-\hat{\epsilon}_{M}\right]^{-1})|mol\rangle}{\omega^{+}-\epsilon_{c}}\right\\}^{-1}\langle
mol|\left[\omega^{+}-\hat{\epsilon}_{M}\right]^{-1}|mol\rangle$
$\displaystyle=$
$\displaystyle\left\\{1-\frac{1}{\omega^{+}-\epsilon_{c}}\sum_{i}\frac{\left|\boldsymbol{{\mu}}_{i}.\boldsymbol{E}_{i}\right|^{2}}{\omega^{+}-\epsilon_{i}}\right\\}^{-1}F^{2}\sum_{i}\frac{\left|\boldsymbol{{\mu}}_{i}.\boldsymbol{E}_{i}\right|^{2}}{\omega^{+}-\epsilon_{i}}$

The average over orientation and length are performed with an approximation
that since the terms are all squares of
$\left|\boldsymbol{{\mu}}_{eg,i}.\boldsymbol{E}_{i}\right|$, we can
individually average over each term. On averaging,

$\displaystyle{G}_{mol,mol}(\omega^{+})=$
$\displaystyle\left\\{1-\frac{1}{\omega^{+}-\epsilon_{c}}\langle\Sigma(\omega)\rangle\right\\}^{-1}\frac{6}{NV^{2}}\langle\Sigma(\omega)\rangle$

Here, $\langle\Sigma(\omega)\rangle$ is used as given in Eq. (50). With this
result the absorption cross-section averaged over the length and orientation
is,

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle-\frac{\omega}{6\epsilon_{o}c\hbar}N\left|\boldsymbol{{\mu}}\right|^{2}Im\left\\{G_{mol,mol}(\omega)\right\\}.$
(59)

##  VI Conclusions

A general solution for Tavis-Cummings model for a system with energetic
disorder which follows a Gaussian distribution is provided in this work.
Analytical solutions are obtained in the limit of large $\mathscr{N}$ for a
simple system where coupling is same for all the molecules and for a general
system where coupling depends in the position and orientation of the
molecules. Conditions for the existence of the polaritonic states is derived.
A general expression for the polaritonic energies correct upto second order in
disorder is derived. For the case where $\epsilon_{a}=\epsilon_{c}$, it is
found that increase in disorder leads to an increase in the Rabi splitting. In
general, polaritonic states are found to be very stable against disorder. They
would always exist if the number density of molecules is sufficiently large.
Our calculation of the absorption spectrum shows that disorder causes the dark
states to turn grey. However, in the case where Rabi splitting is sufficiently
large, disorder plays almost no role in the line width and in this limit the
line width is dominated by lifetime effects. Orientational and position
dependence of the coupling is shown to effectively renormalize the coupling
strength.

##  VII Acknowledgements

The authors thank Prof. Srihari Keshavamurthy (IITK) for his wonderful
lectures titled “Chemistry with Quantum Light”. This work is a result of an
after lecture discussion with him and Prof. Madhav Ranganathan (IITK) and the
authors are grateful to both of them.

## References

  * (1) Anoop Thomas, Jino George, Atef Shalabney, Marian Dryzhakov, Sreejith J. Varma, Joseph Moran, Thibault Chervy, Xiaolan Zhong, Eloise Devaux, Cyriaque Genet, James A. Hutchison, and Thomas W. Ebbesen.  Ground-state chemical reactivity under vibrational coupling to the vacuum electromagnetic field.  Angewandte Chemie International Edition, 55(38):11462–11466, August 2016.
  * (2) A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet, J. Moran, and T. W. Ebbesen.  Tilting a ground-state reactivity landscape by vibrational strong coupling.  Science, 363(6427):615–619, February 2019.
  * (3) Kenji Hirai, Rie Takeda, James A. Hutchison, and Hiroshi Uji-i.  Modulation of Prins Cyclization by Vibrational Strong Coupling.  Angewandte Chemie - International Edition, 59(13):5332–5335, 2020\.
  * (4) Vanessa N. Peters, Md Omar Faruk, Joshua Asane, Rohan Alexander, D’angelo A. Peters, Srujana Prayakarao, Sangeeta Rout, and M. A. Noginov.  Effect of strong coupling on photodegradation of the semiconducting polymer p3ht.  Optica, 6(3):318, March 2019.
  * (5) Battulga Munkhbat, Martin Wersall, Denis G. Baranov, Tomasz J. Antosiewicz, and Timur Shegai.  Suppression of photo-oxidation of organic chromophores by strong coupling to plasmonic nanoantennas.  Science Advances, 4(7):9552, July 2018.
  * (6) Robrecht M. A. Vergauwe, Anoop Thomas, Kalaivanan Nagarajan, Atef Shalabney, Jino George, Thibault Chervy, Marcus Seidel, Eloise Devaux, Vladimir Torbeev, and Thomas W. Ebbesen.  Modification of enzyme activity by vibrational strong coupling of water.  Angewandte Chemie, 131:15468–15472, 10 2019.
  * (7) Xiaolan Zhong, Thibault Chervy, Shaojun Wang, Jino George, Anoop Thomas, James A. Hutchison, Eloise Devaux, Cyriaque Genet, and Thomas W. Ebbesen.  Non-radiative energy transfer mediated by hybrid light-matter states.  Angewandte Chemie - International Edition, 55(21):6202–6206, 2016\.
  * (8) Xiaolan Zhong, Thibault Chervy, Lei Zhang, Anoop Thomas, Jino George, Cyriaque Genet, James A. Hutchison, and Thomas W. Ebbesen.  Energy Transfer between Spatially Separated Entangled Molecules.  Angewandte Chemie - International Edition, 56(31):9034–9038, 2017\.
  * (9) David M. Coles, Niccolo Somaschi, Paolo Michetti, Caspar Clark, Pavlos G. Lagoudakis, Pavlos G. Savvidis, and David G. Lidzey.  Polariton-mediated energy transfer between organic dyes in a strongly coupled optical microcavity.  Nature Materials, 13(7):712–719, May 2014.
  * (10) Katherine Akulov, Dan Bochman, Adina Golombek, and Tal Schwartz.  Long-distance resonant energy transfer mediated by hybrid plasmonic–photonic modes.  The Journal of Physical Chemistry C, 122(28):15853–15860, June 2018\.
  * (11) Thomas W. Ebbesen.  Hybrid light–matter states in a molecular and material science perspective.  Accounts of Chemical Research, 49(11):2403–2412, October 2016.
  * (12) James A. Hutchison, Tal Schwartz, Cyriaque Genet, Eloise Devaux, and Thomas W. Ebbesen.  Modifying chemical landscapes by coupling to vacuum fields.  Angewandte Chemie International Edition, 51(7):1592–1596, January 2012.
  * (13) E. Orgiu, J. George, J. A. Hutchison, E. Devaux, J. F. Dayen, B. Doudin, F. Stellacci, C. Genet, J. Schachenmayer, C. Genes, G. Pupillo, P. Samorì, and T. W. Ebbesen.  Conductivity in organic semiconductors hybridized with the vacuum field.  Nature Materials, 14(11):1123–1129, September 2015.
  * (14) J. J. Pietron, K. P. Fears, J. C. Owrutsky, and B. S. Simpkins.  Electrochemical modulation of strong vibration–cavity coupling.  ACS Photonics, 7(1):165–173, December 2019.
  * (15) Shaelyn R. Casey and Justin R. Sparks.  Vibrational strong coupling of organometallic complexes.  Journal of Physical Chemistry C, 120:28138–28143, 12 2016.
  * (16) Javier Galego, Clàudia Climent, Francisco J. Garcia-Vidal, and Johannes Feist.  Cavity Casimir-Polder Forces and Their Effects in Ground-State Chemical Reactivity.  Physical Review X, 9(2):1–22, 2019.
  * (17) Vladimir P. Zhdanov.  Vacuum field in a cavity, light-mediated vibrational coupling, and chemical reactivity.  Chemical Physics, 535(March):110767, 2020.
  * (18) Jorge A. Campos-Gonzalez-Angulo and Joel Yuen-Zhou.  Polaritonic normal modes in transition state theory.  The Journal of Chemical Physics, 152:161101, 4 2020.
  * (19) Denis G. Baranov, Martin Wersall, Jorge Cuadra, Tomasz J. Antosiewicz, and Timur Shegai.  Novel nanostructures and materials for strong light-matter interactions.  ACS Photonics, 5:24–42, 1 2018.
  * (20) Johannes Flick, Nicholas Rivera, and Prineha Narang.  Strong light-matter coupling in quantum chemistry and quantum photonics.  Nanophotonics, 7(9):1479–1501, September 2018.
  * (21) Felipe Herrera and Jeffrey Owrutsky.  Molecular polaritons for controlling chemistry with quantum optics.  Journal of Chemical Physics, 152, 3 2020.
  * (22) Felipe Herrera and Frank C. Spano.  Theory of Nanoscale Organic Cavities: The Essential Role of Vibration-Photon Dressed States.  ACS Photonics, 5(1):65–79, 2018.
  * (23) Manuel Hertzog, Mao Wang, Jurgen Mony, and Karl Borjesson.  Strong light-matter interactions: A new direction within chemistry.  Chemical Society Reviews, 48:937–961, 2 2019.
  * (24) Derek S. Wang and Susanne F. Yelin.  A roadmap toward the theory of vibrational polariton chemistry.  ACS Photonics, 8:2818–2826, 10 2021.
  * (25) Raphael F. Ribeiro, Luis A. Martinez-Martinez, Matthew Du, Jorge Campos-Gonzalez-Angulo, and Joel Yuen-Zhou.  Polariton chemistry: controlling molecular dynamics with optical cavities.  Chemical Science, 9:6325–6339, 2018.
  * (26) Johannes Feist, Javier Galego, and Francisco J. Garcia-Vidal.  Polaritonic chemistry with organic molecules.  ACS Photonics, 5:205–216, 1 2018.
  * (27) Felipe Herrera and Frank C. Spano.  Cavity-Controlled Chemistry in Molecular Ensembles.  Physical Review Letters, 116(23):1–6, 2016.
  * (28) Gregory D. Scholes.  Polaritons and excitons: Hamiltonian design for enhanced coherence: Hamiltonian design for coherence.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476, 10 2020.
  * (29) T. Botzung, D. Hagenmuller, S. Schutz, J. Dubail, G. Pupillo, and J. Schachenmayer.  Dark state semilocalization of quantum emitters in a cavity.  Physical Review B, 102, 10 2020.
  * (30) Matthew Du and Joel Yuen-Zhou.  Can dark states explain vibropolaritonic chemistry?  Arxiv, 4 2021.
  * (31) Nahum C. Chavez, Francesco Mattiotti, J. A. Mendez-Bermudez, Fausto Borgonovi, and G. Luca Celardo.  Disorder-enhanced and disorder-independent transport with long-range hopping: Application to molecular chains in optical cavities.  Physical Review Letters, 126, 4 2021.
  * (32) R. Houdr$\acute{e}$, R. P. Stanley, and M. Ilegems.  Vacuum-field rabi splitting in the presence of inhomogeneous broadening: Resolution of a homogeneous linewidth in an inhomogeneously broadened system vacuum-field rabi splitting in the presence of inhomogeneous broadening: Resolution of a homogeneous linewidth in an inhomogeneously broadened system.  Physical Review A, 53:2711–2715, 1996.
  * (33) Michael Tavis and Frederick W. Cummings.  Exact solution for anN-molecule—radiation-field hamiltonian.  Physical Review, 170(2):379–384, June 1968.
  * (34) Michael Tavis and Frederick W. Cummings.  Approximate solutions for anN-molecule-radiation-field hamiltonian.  Physical Review, 188(2):692–695, December 1969.
  * (35) P. W. Anderson.  Localized magnetic states in metals.  Physical Review, 124(1):41–53, October 1961.
  * (36) R.D. Mattuck.  A Guide to Feynman Diagrams in the Many-Body Problem: Second Edition.  Dover Books on Physics. Dover Publications, 2012.
  * (37) https://mathworld.wolfram.com/dawsonsintegral.html.
  * (38) G.C. Schatz and M.A. Ratner.  Quantum Mechanics in Chemistry.  Dover Books on Chemistry. Dover Publications, 2012.
  * (39) C.A. Coulson and R. McWeeny.  Coulson’s Valence.  Oxford Chemistry Series. Oxford University Press, 1979.