Unusual Dynamical Properties of Disordered Polaritons in Microcavities

We consider a microcavity containing $N$ quantum emitters labeled by $j$ as sketched in Fig. 1. To enable analytical calculations, we describe the electromagnetic field in the cavity with a single mode labeled by ($C$). The Hamiltonian describing the system reads as

$$\hat{H}=\hat{H}_{L}+\hat{H}_{M}+\hat{H}_{LM},$$

(1)

where the light, matter, and light-matter interaction Hamiltonians are given as

	$\displaystyle\hat{H}_{L}$	$\displaystyle=$	$\displaystyle E_{C}\hat{a}^{\dagger}\hat{a},$
	$\displaystyle\hat{H}_{M}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{N}E_{j}\hat{B}_{j}^{\dagger}\hat{B}_{j},$
	$\displaystyle\hat{H}_{LM}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{N}g_{j}\hat{a}\hat{B}_{j}^{\dagger}+\text{H.c.},$		(2)

respectively. The cavity mode is quantized by the photonic operator $\hat{a}$ and has energy $E_{C}$. The two-level systems, which are described by the operators $\hat{B}_{j}$, refer to general quantum emitters, such as atoms, charges, excitons, spins, electronic or vibrational levels of molecules. In the following, we focus on molecular excitations because of their experimental relevance, but our findings are generally valid for the other before mentioned systems. The excitation energies $E_{j}$ are distributed according to the Lorentz function

$\displaystyle P(E_{j})=\frac{1}{\pi}\frac{\sigma}{\left(E_{j}-E_{M}\right)^{2}+\sigma^{2}},$

(3)

where $E_{M}$ is the center of the probability distribution and $\sigma$ is its width, i.e., the disorder parameter. We emphasize that many of our results hold for arbitrary disorder distributions. The light-matter couplings can be expressed in terms of physical quantities as $g_{j}=\left(\frac{\hbar E_{C}}{2\epsilon_{0}V}\right)^{(1/2)}\mathbf{d}_{j}\cdot\mathbf{\hat{E}}$, where $V$ is the volume of the cavity, $\mathbf{d}_{j}$ is the dipole moment of the $j$-th molecule, and $\mathbf{\hat{E}}$ is the polarization of the cavity mode. For simplicity, we assume here a homogeneous coupling $g_{j}=g$, but generalizations to the inhomogeneous case are straightforward. Many calculations are performed in a thermodynamic limit which is defined as $g\rightarrow 0$ and $N\rightarrow\infty$ such that $g\sqrt{N}$ is constant. As $g\propto V^{-1/2}$, the thermodynamic limit implies a large cavity volume, but a constant molecule density.

The Hamiltonian in Eq. (1) is the celebrated Fano-Anderson model, which has been originally developed to understand the impact of continua on discrete levels and asymmetries in absorption spectra [41]. Aside from other applications, this and generalized Fano-Anderson models are also deployed to investigate transport through nanoscopic and mesoscopic systems [42, 43, 44, 45]. A multi-mode version of this model has been numerically studied in Ref. [46]. A recent investigation of the spectral and transport properties for a disorder distribution with compact support can be found in Refs. [47, 48], which make use of the exact expression of the eigenstates instead of the unifying Green’s function approach considered here. We note that the transport setup investigated in Refs. [47] is different from the one considered in our work in Sec. V.

The homogeneous system is defined for vanishing disorder $\sigma=0$. For simplicity, we focus here on the resonant system $E_{M}=E_{C}$. Using the excited states of the molecules $\left|e_{j}\right>_{M}=\hat{B}_{j}^{\dagger}\left|g\right>$, where $\left|g\right>$ is the collective ground state of the molecule ensemble, we define the collective excitations $\left|e_{k}\right>_{M}=\frac{1}{\sqrt{N}}\sum_{j}e^{ikj}\left|e_{j}\right>_{M}$ with $k=2\pi l/N$ and $l=0,...,N-1$. The Fock states of the cavity mode are denoted by $\left|n\right>_{L}$. In terms of these states, the two special eigenstates of the Hamiltonian

$$\left|\psi_{up/down}\right>=\frac{1}{\sqrt{2}}\left(\left|g\right>_{M}\left|1\right>_{L}\pm\left|e_{k=0}\right>_{M}\left|0\right>_{L}\right),$$

(4)

with energies $\epsilon_{0,1}=E_{M}\pm g\sqrt{N}$ are called the upper and lower polaritons, respectively. Their energy difference $\Omega_{R}=2g\sqrt{N}$, which increases with the square root of the molecule number, can be measured as a collective Rabi splitting $\Omega_{R}$ in spectroscopic experiments. The homogeneous state $\left|BS\right>_{M}=\left|e_{k=0}\right>_{M}$, which mixes with the cavity light field, is commonly denoted as the bright state. In contrast, the states $\left|DS_{k}\right>_{M}=\left|e_{k\neq 0}\right>_{M}$, which have energy $\epsilon_{k}=E_{M}$, completely decouple from the cavity light field and are denoted as dark states. We note that the bright and dark states can be also defined for an inhomogeneous coupling $g_{j}$. In this case, the bright state reads as $\left|BS\right>_{M}\propto\sum_{j}g_{j}\left|e_{j}\right>_{M}$, and the dark states are orthogonal to this state.

For the disordered system, the dark states are no longer eigenstates of the system, but all eigenstates have contributions of the cavity mode, the bright state and the dark states, as the matrix elements $\left<BS\right|\hat{H}_{M}\left|DS_{k}\right>\neq 0$ of the matter Hamiltonian in Eq. (2) become finite. The respective contributions of the cavity mode, the bright state and the dark states to the eigenstate is given respectively by the bright-state LDOS and the dark-state LDOS, which will be introduced in Sec. III. The mixing of bright and dark states influences the spectroscopic spectra in Sec. III and leads to the turnover in the relaxation and transport rates in Secs. IV and V.

As shown in Fig. 1, the Fano-Anderson model has a tree-shaped coupling structure, which can be solved analytically in Laplace space as shown in Appendix A. In doing so, we find an explicit expression of the single-particle retarded Green’s functions

$$G_{X,Y}(t)\equiv i\Theta(t)\left<\hat{A}_{X}(t),\hat{A}_{Y}^{\dagger}\right>,$$

(5)

where $X\in\left\{C,j,BS,DS_{k}\right\}$, and $\Theta(t)$ denotes the Heaviside step function. Thereby, $\hat{A}_{C}^{\dagger}=\hat{a}^{\dagger}$, $\hat{A}_{j}^{\dagger}=\hat{B}_{j}^{\dagger}$ and $\hat{A}_{BS}^{\dagger}$ ($\hat{A}_{DS_{k}}^{\dagger}$ ) creates a bright state (dark state) excitation. In Laplace space, the Green’s functions read as [36, 37]

	$\displaystyle G_{C,C}(z)$	$\displaystyle=$	$\displaystyle\frac{1}{\left(z+iE_{C}(z)\right)},$
	$\displaystyle G_{C,j}(z)$	$\displaystyle=$	$\displaystyle-i\frac{g}{\left(z+iE_{C}(z)\right)\left(z+iE_{j}\right)}=G_{j,C}(z),$
	$\displaystyle G_{i,j}(z)$	$\displaystyle=$	$\displaystyle\frac{\delta_{i,j}}{z+iE_{j}}-\frac{g^{2}}{\left(z+iE_{i}\right)\left(z+iE_{C}(z)\right)\left(z+iE_{j}\right)},$

with the auxiliary function

$\displaystyle E_{C}(z)$

$\displaystyle=$

$\displaystyle E_{C}-i\sum_{j=1}^{N}\frac{g^{2}}{z+iE_{j}}.$

(7)

Up to factors $(z+iE_{j})$, the nominator of the third term in Eq. (LABEL:eq:greensFktDD) defines a polynomial $\mathcal{P}(z)$ of order $N+1$,

$$\mathcal{P}(z)=\left(z+iE_{C}(z)\right)\prod_{j}^{N}(z+iE_{j}),$$

(8)

which is equivalent to the characteristic polynomial of the Hamiltonian (1) when replacing $z\rightarrow-iE$. The inverse Laplace transformation can be expressed in terms of the roots of the characteristic polynomial, i.e., the poles of the Green’s function, as

$$G_{X,Y}(t)=\sum_{\alpha=1}^{N+1}A_{\alpha}e^{z_{\alpha}t},$$

(9)

where $A_{\alpha}=2\pi i\lim\limits_{z\rightarrow z_{\alpha}}(z-z_{\alpha})G_{X,Y}(z)$. Note that the pole of the first term in $G_{j,j}(z)$ of Eq. (LABEL:eq:greensFktDD) is not a pole of $G_{j,j}(z)$, as it cancels with a pole of the second term.

Here, we consider the system in the thermodynamic limit $N\rightarrow\infty$ and $g\rightarrow 0$ such that $g^{2}N=\text{const}$, and derive an effective Hamiltonian which exactly reproduces the dynamics of the single-particle Green’s function. The disorder-averaged Green’s function in the thermodynamic limit is defined by

$$\overline{G}_{X,Y}(z)\equiv\int dE_{1}...\int dE_{N}G_{X,Y}(z)P(E_{1})...P(E_{N}),$$

(10)

where $P(E_{j})$ is the disorder distribution function of $E_{j}$. When evaluating $\overline{G}_{X,Y}(z)$, $\sum_{j=1}^{N}\frac{g^{2}}{z+iE_{j}}\rightarrow\Pi(z)$ in Eq. (7) becomes a smooth function which depends on the statistics of $E_{j}$. Moreover, $\overline{G}_{i_{1},j_{1}}(z)=\overline{G}_{i_{2},j_{2}}(z)$ and $\overline{G}_{C,j_{1}}(z)=\overline{G}_{C,j_{2}}(z)$ for all $i_{1},i_{2},j_{1},j_{2}$, i.e., the disorder-averaged Green’s function is homogeneous. As a consequence, the Green’s function is block-diagonal in the basis of bright and dark states introduced in Sec. II.2, i.e., $\overline{G}_{e_{k_{1}},e_{k_{2}}}(z)\propto\delta_{k_{1},k_{2}}$ and $\overline{G}_{C,e_{k}}(z)\propto\delta_{k_{,}0}$. We use the disorder-averaged Green’s function to define an effective Hamiltonian by

$$\overline{G}(z)\equiv\frac{1}{z+iH_{eff}(z)}.$$

(11)

Note that the effective Hamiltonian depends on $z$ and the distribution of $E_{j}$. The block structure of $\overline{G}(z)$ in the basis of bright and dark states translates into a block structure of $H_{eff}(z)$, i.e.,

$$H_{eff}(z)=\left[\begin{array}[]{cccc}H_{eff}^{(C,BS)}(z)&&&\makebox(0.0,0.0)[]{\text{\Large 0}}\\ &H_{eff}^{(DS)}(z)&&\\ &&\ddots&\\ \makebox(0.0,0.0)[]{\text{\Large 0}}&&&H_{eff}^{(DS)}(z)\end{array}\right],$$

(12)

where $H_{eff}^{(C,BS)}(z)\in\mathbb{C}^{2\times 2}$ describes the interaction of the cavity mode and the bright state and $H_{eff}^{(DS)}(z)\in\mathbb{C}$ describes the dynamics of the dark states. As the effective dark state Hamiltonians are equal for all $k\neq 0$, we have suppressed the index $k$. The effective Hamiltonian Eq. (12) appears to suggest that the dynamics of the bright and dark states is decoupled. Yet, the influence of the dark states is incorporated in the non-Hermitian nature of $H_{eff}(z)$, which effectively leads to a dissipative dynamics.

For the Lorentz distribution in Eq. (3), the evaluation of the disorder-averaged Green’s function is straightforward and is equivalent to simply replacing $E_{j}\rightarrow E_{M}-i\sigma\equiv E_{M}^{(\sigma)}$ [49, 50]. Because of this simple replacement rule, the thermodynamic limit considered in Eq. (10) is actually equivalent to the disorder average, such that the following results are also valid for finite molecule numbers. The disorder-averaged Green’s functions are explicitly given as

	$\displaystyle\overline{G}_{C,C}(z)$	$\displaystyle=$	$\displaystyle\frac{1}{\left(z+iE_{C}(z)\right)},$		(13)
	$\displaystyle\overline{G}_{C,j}(z)$	$\displaystyle=$	$\displaystyle-i\frac{g}{\left(z+iE_{C}(z)\right)\left(z+iE_{M}^{(\sigma)}\right)}=\overline{G}_{j,C}(z),$
	$\displaystyle\overline{G}_{i,j}(z)$	$\displaystyle=$	$\displaystyle\frac{\delta_{i,j}}{z+iE_{M}^{(\sigma)}}$
		$\displaystyle-$	$\displaystyle\frac{g^{2}}{\left(z+iE_{M}^{(\sigma)}\right)\left(z+iE_{C}(z)\right)\left(z+iE_{M}^{(\sigma)}\right)},$

where now $E_{C}(z)=E_{C}+g^{2}N/(z+E_{M}^{(\sigma)})$. Transforming this into the basis of the bright and dark states, the blocks of the effective Hamiltonian in Eq. (12) are given as

$$H_{eff}^{(C,BS)}=\left[\begin{array}[]{cc}E_{C}&\frac{\Omega}{2}\\ \frac{\Omega}{2}&E_{M}^{(\sigma)}\end{array}\right]$$

(14)

with the Rabi frequency of the homogeneous system $\Omega=2g\sqrt{N}$ and $H_{eff}^{(DS)}=E_{M}^{(\sigma)}$. Note that the $z$-independence of the effective Hamiltonian is a consequence of the Lorentz distribution of $E_{j}$. The complex-valued eigenenergies of $H_{eff}^{(C,BS)}$ in Eq. (14)

	$\displaystyle\epsilon_{1,2}$	$\displaystyle=$	$\displaystyle\frac{E_{C}+E_{M}^{(\sigma)}}{2}$		(15)
		$\displaystyle\pm$	$\displaystyle\frac{1}{2}\sqrt{\left(E_{C}-E_{M}^{(\sigma)}\right)^{2}-\Omega^{2}}$

generalize the eigenenergies of the two polaritons for the disordered case. As it will become more clear when discussing spectral properties, the real part is related to the spectral position and the imaginary part to the spectral width of the absorption line shape. The Rabi splitting of the disordered system can be thus defined as $\Omega_{R}\equiv\text{Re}\left(\epsilon_{2}-\epsilon_{1}\right)$. The finite imaginary part appears due to the mixing of the bright and dark states, which eventually gives rise to a decaying occupation of the subsystem consisting of cavity mode and bright state.

Resonant system. For the resonant system $E_{C}=E_{M}$, the eigenenergies are given as

$\displaystyle\epsilon_{1,2}$

$\displaystyle=$

$\displaystyle E_{M}-i\frac{\sigma}{2}\pm\frac{1}{2}\left(\Omega^{2}-\sigma^{2}\right)^{\frac{1}{2}}$

(16)

and are depicted in Figs. 2(a) and 2(c). For very small or very large $\sigma$, they approximately read as

	$\displaystyle\sigma\ll\Omega:\qquad\epsilon_{\mu}$	$\displaystyle=$	$\displaystyle E_{M}-i\frac{\sigma}{2}\pm\left(\frac{\Omega}{2}-\frac{\sigma^{2}}{\Omega}\right),$		(17)
	$\displaystyle\sigma\gg\Omega:\qquad\epsilon_{\mu}$	$\displaystyle=$	$\displaystyle\begin{cases}E_{M}-i\sigma+i\frac{\Omega^{2}}{2\sigma}&\mu=1\\ E_{M}-i\frac{\Omega^{2}}{4\sigma}&\mu=2.\end{cases}$		(18)

These two limiting cases can be clearly seen in Figs. 2(a) and 2(c). For small $\sigma$, the two roots have different real parts, and equal imaginary parts. For large $\sigma$, the real parts are equal, but the imaginary parts differ. The agreement of either the real or the imaginary part is a consequence of $E_{C}=E_{M}$. Interestingly, for $\sigma=\Omega\equiv\sigma_{EP}\approx 0.09eV$, we observe an exceptional point, which has been frequently investigated recently [51]. In Secs. III.1 and III.2, we discuss signatures of the exceptional point in the absorption spectrum. In the following, we denote the region before and after the exceptional point as the underdamped (small $\sigma$) and overdamped (large $\sigma$) regime, respectively.

Off-resonant system. For the off-resonant system with $E_{M}<E_{C}$ in Figs. 2(b) and 2(d), the real and imaginary parts of both eigenvalues are well separated for all values of $\sigma$. Even for a small detuning of $E_{M}$ and $E_{C}$, the exceptional point observed in the resonant system does not exist. The eigenvalue $\epsilon_{1}$ is dominated by the molecular excitations, while $\epsilon_{2}$ is dominated by the cavity excitation. The imaginary part (describing the spectral width) of $\epsilon_{2}$ is much smaller than the imaginary part of $\epsilon_{1}$, as the light and matter become increasingly decoupled for larger disorder due to the decreasing density at $E_{j}\approx E_{C}$ in the off-resonant system.

First, we investigate the spectroscopic response of the cavity related to the perturbation operator $\hat{V}=\hat{a}+\hat{a}^{\dagger}$, which will be denoted as the cavity absorption spectrum in the following. In terms of the Green’s function, the cavity absorption can be expressed as

$\displaystyle\chi_{C}(\omega)$

$\displaystyle=$

$\displaystyle-\text{Im}\;G_{C,C}\left(-i\omega+0^{+}\right)=\pi\nu_{C}(\omega)$

(20)

and is thus directly proportional to the cavity LDOS. For the Lorentz distribution in the thermodynamic limit, we evaluate Eq. (20) using the disorder-averaged Green’s function $\overline{G}_{CC}(z)$ given in Eq. (13), which can be transformed into

	$\displaystyle\chi_{C}(\omega)$	$\displaystyle=$	$\displaystyle\sum_{\mu=1,2}\frac{-1}{\pi}\text{Im }\left[\frac{A_{\mu}}{\omega-\epsilon_{\mu}}\right],$
	$\displaystyle A_{\mu}$	$\displaystyle=$	$\displaystyle\frac{-i\epsilon_{\mu}+iE_{M}+\sigma}{-i\epsilon_{\mu}+i\epsilon_{\overline{\mu}}},$		(21)

where the coefficients $A_{\mu}$ can be evaluated in terms of the energies $\epsilon_{\mu}$ in Eq. (15). Details of the derivation can be found in Appendix B.1. We have introduced $\overline{\mu}$ as $\overline{\mu}\neq\mu$ for a compact notation. If $A_{\mu}$ is real-valued, the cavity absorption spectrum is given by two Lorentz functions centered at $\text{Re}\,\epsilon_{\mu}$ and having spectral width $\text{Im}\,\epsilon_{\mu}$. For complex-valued $A_{\mu}$, the shape deviates from the pure Lorentzian function. As the finite imaginary part of the eigenenergies describes the mixing of the bright and dark states, the dark states thus determine the functional form of the cavity LDOS.

In Figs. 3(a) and 3(d), we depict the cavity absorption spectrum for the resonant case $E_{M}=E_{C}$ in the underdamped and overdamped regimes, respectively. In the underdamped regime, we observe two Lorentzian peaks symmetrically located around $E_{M}$. Here, the eigenenergies fulfill $(-\epsilon_{1}+E_{M})=(\epsilon_{2}-iE_{M})^{*}$ as can be seen in Fig. 2, which explains the symmetry of the peaks when evaluating Eq. (21). In the overdamped regime, we observe a single peak with Lorentzian shape. Here, the eigenenergies fulfill $\epsilon_{1}+i\sigma/2=\epsilon_{2}-i\sigma/2$, such that there are actually two Lorentzian peaks centered at the same position $E_{M}$. The eigenenergy $\epsilon_{2}$ with a small imaginary part dominates the spectrum. Using Eq. (18) to evaluate $A_{1}$, we find that $A_{1}\rightarrow 0$ for $\sigma\rightarrow\infty$, such that the signatures of the eigenenergy $\epsilon_{1}$ are strongly suppressed. A physical interpretation of the eigenenergies in the overdamped regime is given in the next section.

Figures 3(g) and 3(j) depict the cavity absorption spectrum for the off-resonant system $E_{M}<E_{C}$ in the underdamped and overdamped regimes, respectively. Similar to the resonant system, we observe a peak close the position of the cavity frequency $E_{C}\approx\text{Re}\;\epsilon_{2}$. In Fig. 3(g), we mark a second peak related to the molecular eigenenergy $\epsilon_{1}$ located close to $E_{M}$, which is very small as the corresponding $A_{1}\rightarrow 0$ for large $\left|E_{C}-E_{M}\right|$. The observations in the overdampend regime in Fig. 3(j) are very similar, but the molecule peak is now completely suppressed.

Next, we consider the absorption spectrum which is related to the perturbation operator $\hat{V}=\sum_{j}D_{j}\left(\hat{B}_{j}+\hat{B}_{j}^{\dagger}\right)$, which we denote as matter absorption in the following. The coupling coefficients between the molecules and the probe field are $D_{j}=\mathbf{d}_{j}\cdot\mathbf{\hat{E}}_{p}$, where $\mathbf{\hat{E}}_{p}$ is the probe field. It measures directly the spectroscopic properties of the molecules instead of the cavity field as in Sec. III.1. Assuming a homogeneous $D_{j}=D$, we can express the matter absorption spectrum in terms of the single-particle Green’s functions as

	$\displaystyle\chi_{M}(\omega)$	$\displaystyle=$	$\displaystyle-D^{2}\sum_{i,j=1}^{N}\text{Im }G_{i,j}(-i\omega+0^{+})$		(22)
		$\displaystyle=$	$\displaystyle ND^{2}\pi\nu_{BS}(\omega).$

In contrast to the absorption spectrum of uncorrelated molecules, where only the diagonal elements of the Green’s function are taken into account, the absorption of the molecular polaritons includes all elements of the Greens’ function as the cavity induces strong coherences between the molecules.

The summations in Eq. (22) project the Green’s function onto the bright state introduced in Sec. II.2, such that the matter absorption spectrum is directly proportional to the bright-state LDOS. We note that the relation of the matter absorption and the bright-state density of states in Eq. (22) is also correct for inhomogeneous couplings as long as $g_{j}\propto D_{j}$, i.e., the polarizations of the cavity and the probe field are parallel.

For the Lorentz distribution in the thermodynamic limit, we can use the disorder-averaged Green’s function in Eq. (13) to evaluate the matter absorption, such that we find after some steps

	$\displaystyle\chi_{M}(\omega)$	$\displaystyle=$	$\displaystyle\sum_{\mu=1,2}\frac{-D^{2}}{\pi}\text{Im }\left[\frac{A_{\mu}}{\omega-\epsilon_{\mu}}\right],$
	$\displaystyle A_{\mu}$	$\displaystyle=$	$\displaystyle\frac{g^{2}}{\left(-\epsilon_{\mu}+iE_{M}+\sigma\right)\left(-i\epsilon_{\mu}+i\epsilon_{\overline{\mu}}\right)}.$		(23)

As explicitly demonstrated in Appendix B.2, the matter absorption spectrum can be expressed in terms of the cavity single-particle Green’s function as

$\displaystyle\chi_{M}(\omega)$

$\displaystyle=$

$\displaystyle ND^{2}\text{Im}\;G_{CC}\left(i\omega+iE_{M}+iE_{C}+\sigma\right),$

(24)

i.e., using a complex frequency $\omega\rightarrow-\omega+E_{M}+E_{C}+i\sigma$. Thus, the absorption spectra of the cavity mode and of the matter are directly related to each other.

The molecular absorption is depicted in Figs. 3(b) and 3(e) for the resonant system $E_{M}=E_{C}$. Because of Eq. (24), $\chi_{M}(\omega)$ has the same functional dependence as the cavity absorption Eq. (21) except for the complex-valued frequency. The discussions about the cavity absorption are thus also valid for the matter absorption. In the underdamped regime in Fig. 3(b), we find two peaks corresponding to the two polaritons which approximately have a Lorentzian shape, similar to the cavity absorption in Fig. 3(a). Interestingly, in Figs. 3(b), 3(e), 3(h), and 3(k) the matter absorption vanishes completely at $\omega=E_{M}$ due to level repulsion of the molecules which are in resonance with $E_{C}$. In contrast, the level repulsion for the cavity LDOS in Figs. 3(a), 3(d), 3(g), and 3(j) is not complete as the cavity mode is interacting with a continuous spectrum which leads to a coarse graining of the level repulsion.

The complete suppression of the absorption in Fig. 3 is reminiscent of the electromagnetically-induced transparency [52] (EIT) and the related vacuum-induced transparency (VIT) [53, 49] appearing in atomic and molecular three-level systems. In fact, the absorption suppression in Figs. 3(b), 3(e), 3(h), and 3(k) can be also understood as a destructive interference between two excited states. However, in contrast to the EIT and VIT, which considers individual three-level atoms or molecules coupled by a light field, in this paper, the light field itself is also a quantum state and couples the molecules collectively. Moreover, the suppressed absorption observed here is based on a non-local superposition of the quantum emitters (associated with the bright state) and is thus a collective effect, while the EIT and VIT are based on a destructive interference effect within individual atoms or molecules.

In the overdamped regime depicted in Figs. 3(e) and  3(k), the matter absorption strongly deviates from the cavity absorption despite their close relationship via Eq. (24). In contrast to the cavity absorption, both eigenenergies depicted in Fig. 2 are now relevant. The matter absorption is a superposition of two Lorentz functions which have different widths $\text{Im}\,\epsilon_{1}$ and $\text{Im}\,\epsilon_{2}$, but the same peak position of $\omega=E_{M}=E_{C}$. One peak has a positive amplitude, while the other has a negative amplitude, leading to a complete destructive interference at the cavity energy $E_{C}$. This complete destructive interference is a consequence of the level repulsion of the cavity mode and the molecular excitation energies with $E_{j}=E_{C}$. While the two peaks in the underdamped regime can be associated with the two polaritons, this picture breaks down in the overdamped regime. The distinct behavior can be understood in terms of the imaginary parts in Fig. 2. For an increasing disorder $\sigma$, the molecular excitations and the cavity mode gradually decouple. The spectral features of the molecular excitations thus approach the original Lorentz distribution with width $\sigma$, while the spectral width of the cavity continuously vanishes. In this regime, the polaritonic excitations are not well defined.

The off-resonant system $E_{M}<E_{C}$ in the underdamped regime depicted in Fig. 3 (h) is not symmetric. The molecule peak close to $\omega=E_{M}$ is significantly broader than the cavity peak close to $\omega=E_{C}$, which is directly related to the $\text{Im}\;\epsilon_{1,2}$ in Fig. 2. The cavity LDOS is thus only weakly influenced by the energetic disorder because of the energy splitting. The observations in the overdamped regime in Fig. 3(k) are similar to Fig. 3(e) but not symmetric.

Finally, we consider the total density of states to clarify the role of the bright and dark states introduced in Sec. II.2 in the presence of disorder. As the total density of states does not depend on the basis, it can be expressed either in the local basis of the molecules $j$ or in the basis of bright and darks states,

	$\displaystyle\nu(\omega)$	$\displaystyle\equiv$	$\displaystyle\nu_{C}(\omega)+\sum_{j}\nu_{j}(\omega)$		(25)
		$\displaystyle=$	$\displaystyle\nu_{C}(\omega)+\nu_{BS}(\omega)+\sum_{k=1}^{N-1}\nu_{DS_{k}}(\omega),$

i.e, it can be partitioned in the cavity LDOS, $\nu_{C}(\omega)$, the bright-state LDOS, $\nu_{BS}(\omega)$, and $N-1$ terms of the dark-state LDOS. In the thermodynamic limit, the cavity and bright-state LDOS converge to their respective limits given in Sec. III.1 and Sec. III.2, while the dark-state LDOS converge to $\nu_{DS_{k}}(\omega)\rightarrow(N-1)P(\omega)$. Thus, from the LDOS of the molecules in the non-interacting case $g=0$, which is $NP(\omega)$, one molecules is subtracted, which now forms the bright-state LDOS $\nu_{BS}(\omega)$.

In the homogeneous and resonant system discussed in Sec. II.2, we have $\nu_{C}(\omega)=\nu_{BS}(\omega)=\frac{1}{2}\left(\delta(\omega-E_{M}+g\sqrt{N})+\delta(\omega-E_{M}-g\sqrt{N})\right)$ and $\nu_{DS}(\omega)=(N-1)\delta(\omega-E_{M})$. Comparing these expressions with the LDOS considered in this section, we find that the disorder turns the delta functions into distribution functions with finite widths. In the homogeneous case the molecular excitations can be either classified as bright or dark states. Because of the disorder, bright and dark states are now mixed and both contribute to the formation of the eigenstates. On the average, the contributions are thereby proportional to $\nu_{BS}(\omega)$ or $\nu_{DS}(\omega)$, respectively. We recall that the dark states determine the functional shape of $\nu_{C}(\omega)$ and $\nu_{BS}(\omega)$, as the dark states are coupled to the bright state for a finite disorder.

In realistic spectroscopic experiments, one simultaneously measures the cavity and matter absorption, i.e., $\chi(\omega)=\alpha_{C}\cdot\chi_{C}(\omega)+\alpha_{M}\cdot\chi_{M}(\omega)$ with coefficients $\alpha_{C}$ and $\alpha_{M}$. For a fixed Rabi frequency $\Omega^{2}=4g^{2}N\propto N/V$, one can thus harness the scaling $\chi_{C}\propto\nu_{C}$ and $\chi_{M}\propto N\cdot\nu_{BS}\propto V\cdot\nu_{BS}$ to distinguish experimentally between both contributions via changing the cavity volume $V$ while keeping the molecule density constant.

To begin with, the occupation of the donor is given as

$$n_{1}(t)=\left|G_{1,1}(t)\right|^{2},$$

(26)

where the Green’s function is defined in Eq. (LABEL:eq:greensFktDD). To facilitate an analytical treatment, we apply the disorder average in the thermodynamic limit defined in Eq. (10). Yet, to account correctly for the microscopic dynamics of the donor occupation, the average is not applied to $E_{1}$. Specifying for the Lorentzian disorder, the Green’s function can be written in Laplace space as

$\displaystyle\overline{G}_{1,1}(z)$

$\displaystyle=$

$\displaystyle\frac{1}{z+iE_{1}}-\frac{g^{2}\left(z+iE_{M}+\sigma\right)}{\left(z+iE_{1}\right)\mathcal{P}(z)}.$

(27)

Thereby, the roots of the polynomial $\mathcal{P}(z)=\mathcal{P}_{0}(z)+\mathcal{P}_{1}(z)$ with

	$\displaystyle\mathcal{P}_{0}(z)$	$\displaystyle=$	$\displaystyle\left(z+iE_{C}\right)\left(z+iE_{M}+\sigma\right)\left(z+iE_{1}\right),$
			$\displaystyle+g^{2}N\left(z+iE_{1}\right),$
	$\displaystyle\mathcal{P}_{1}(z)$	$\displaystyle=$	$\displaystyle g^{2}\left(z_{\alpha}+iE_{M}+\sigma\right)$

determine the inverse Laplace transformation in Eq. (9). Note that $z=-iE_{1}$ is not a root of $\overline{G}_{1,1}(z)$ as the corresponding terms cancel exactly. The third-order polynomial is partitioned into the two parts $\mathcal{P}_{0}(z)$ and $\mathcal{P}_{1}(z)$. The first part can be solved analytically and we obtain the unperturbed roots $z_{1}^{(0)}=-i\epsilon_{1}$ and $z_{2}^{(0)}=-i\epsilon_{2}$ and $z_{3}^{(0)}=-iE_{1}$, with $\epsilon_{1,2}$ given in Eq. (15). The perturbative part $\mathcal{P}_{1}(z)$ gives a correction to $z_{\mu}^{(0)}$ in leading order of $g^{2}$ and is found by applying the PPT in Eq. (70) introduced in Appendix C.1. In doing so, we find that the two roots $z_{1}=z_{1}^{(0)}=-i\epsilon_{1}$ and $z_{2}=z_{1}^{(0)}=-i\epsilon_{2}$ remain unchanged in leading order, while the third root obtains a finite real part such that $z_{3}=-iE_{1}-\gamma(E_{1})+i\mathcal{O}\left(g^{2}\right)+\mathcal{O}\left(g^{4}\right)$ with

$\displaystyle\gamma(\omega)$

$\displaystyle=$

$\displaystyle g^{2}\nu_{C}(\omega),$

(28)

where the cavity LDOS in Eq. (21) can be explicitly written as

	$\displaystyle\nu_{C}(\omega)$	$\displaystyle=$	$\displaystyle\frac{1}{\pi}\frac{\kappa\sigma}{\left(-\omega+E_{C}-\kappa\left(-\omega+E_{M}\right)\right)^{2}+\kappa^{2}\sigma^{2}},$
	$\displaystyle\kappa$	$\displaystyle=$	$\displaystyle\frac{g^{2}N}{(-\omega+E_{M})^{2}+\sigma^{2}}.$		(29)

We note that Eq. (28) agrees with the result in Ref. [47]. The imaginary part shift $\propto g^{2}$ (i.e., the energy shift) can be neglected as it is not the leading order term. Using the roots to perform the inverse Laplace transformation, the Green’s function becomes as a function of time

	$\displaystyle G_{1,1}(t)$	$\displaystyle=$	$\displaystyle A_{1}e^{z_{1}t}+A_{2}e^{z_{2}t}+A_{3}e^{z_{3}t},$
	$\displaystyle A_{1}$	$\displaystyle=$	$\displaystyle\frac{g^{2}(z_{1}+iE_{M}+\sigma)}{\left(z_{1}+iE_{1}\right)^{2}(z_{1}-z_{2})}\rightarrow 0,$
	$\displaystyle A_{2}$	$\displaystyle=$	$\displaystyle\frac{g^{2}(z_{2}+iE_{M}+\sigma)}{\left(z_{2}+iE_{1}\right)^{2}(z_{2}-z_{1})}\rightarrow 0,$
	$\displaystyle A_{3}$	$\displaystyle=$	$\displaystyle\frac{g^{2}(-iE_{1}+iE_{M}+\sigma)}{\gamma(E)(-iE_{1}-z_{1})(-iE_{1}-z_{2})}\rightarrow 1.$		(30)

In the thermodynamic limit $N\rightarrow\infty,g\rightarrow 0$, which we have assumed throughout the calculation, the amplitudes $A_{1},A_{2}\rightarrow 0$, such that the occupation of $n_{1}(t)$ is solely determined by the real part of the root $z_{3}$. For consistency, we finally average $\gamma(E_{1})$ in Eq. (IV.1) over the donor energy $E_{1}$, which is distributed according to the Lorentz distribution in Eq. (3), such that the disorder-averaged relaxation rate reads as

	$\displaystyle\overline{\gamma}$	$\displaystyle=$	$\displaystyle\int dE_{1}\gamma(E_{1})P(E_{1})$		(31)
		$\displaystyle=$	$\displaystyle\frac{g^{2}}{\pi}\frac{\sigma+\frac{g^{2}N}{2\sigma}}{\left(E_{M}-E_{C}\right)^{2}+\left(\sigma+\frac{g^{2}N}{2\sigma}\right)^{2}},$

whose physical behavior will be discussed in the next section.

We note that the relaxation rate in Eq. (28) is exact in the thermodynamic limit $g\rightarrow 0,N\rightarrow\infty$ and holds also for general disorder distributions. The derivation is still valid for a Green’s function with an arbitrary number of poles. Possible branch cuts of the Green’s function vanish in the thermodynamic limit as the second term in Eq. (27) is proportional to $g^{2}$.

The energy-resolved relaxation rate in Eq. (IV.1) and the disorder-averaged relaxation rate in Eq. (31) are depicted in Figs. 4 and 5. The solid lines depict the analytical results, while the symbols depict the numerical results of the finite-size simulations, which are described in details in Appendix D. Several points are worthwhile to be discussed:

Disorder dependence. The energy-resolved relaxation rate $\gamma(E_{1})$ is shown in Fig. 4(a,b) as a function of disorder $\sigma$ for different values of $E_{1}$ for the resonant $E_{M}=E_{C}$ and off-resonant $E_{M}<E_{C}$ systems. We observe a turnover as a function of $\sigma$ for all values of $E_{1}$. For $E_{1}\neq E_{M}$ and/or $E_{M}\neq E_{C}$, we find from Eq. (IV.1) in the limiting cases

$$\gamma(E_{1})\propto\begin{cases}\frac{\sigma}{N}&\sigma\ll g\sqrt{N}\\ \frac{N}{\sigma}&\sigma\gg g\sqrt{N}\end{cases}.$$

(32)

Only in the special case $E_{1}=E_{M}=E_{C}$, $\gamma(E_{1})$ is monotonically decreasing as a function of $\sigma$ (not shown). As $\gamma(E_{1})\propto\nu_{C}(E_{1})$, the relaxation rate depends on the distance from the peak $\gamma(E_{1})\propto\left(E_{1}-\text{Re}\;\epsilon_{\mu}\right)^{2}$ and the peak width $\gamma(E_{1})\propto\text{Im}\;\epsilon_{\mu}$ of the peaks in Figs. 3(a), 3(d), 3(g), and 3(j).

For example, in the resonant case $E_{M}=E_{C}$ for an energy $E_{1}>E_{M}$ and $E_{1}-E_{M}\ll\Omega_{R}$, the right peak in Figs. 3(a) and 3(d) related to the eigenenergy $\epsilon_{2}$ determines $\gamma(E_{1})\propto\nu_{C}(E_{1})$. For small $\sigma$, we find from Eq. (18) that $\text{Re}\;\epsilon_{\mu}-E_{M}\propto\sqrt{N}$ and $\text{Im}\;\epsilon_{\mu}\propto\sigma$, which explains the behavior for small $\sigma$ in Eq. (32). Likewise for large $\sigma$, we find from Eq. (18) that $\text{Re}\;\epsilon_{\mu}=E_{M}=\text{const}$ and $\text{Im}\;\epsilon_{\mu}\propto N/\sigma$, which explains the behavior for large $\sigma$ in Eq. (32).

Energy dependence. Overall, the relaxation rate is larger for energies $\left|E_{1}-E_{C}\right|\approx 0$. For the selected energies in Figs. 4(a) and 4(b), the energies $E_{1}=0.9375\;\text{eV}$ and $E_{1}=1.025\;\text{eV}$ exhibit the largest relaxation rates, while the energies $E_{1}=0.85\;\text{eV}$ and $E_{1}=1.2\;\text{eV}$ exhibit the smallest relaxation rates. This is due to the shape of $\nu_{C}(E_{1})$, which decays with $1/\left(E_{1}-E_{C}\right)^{2}$ away from the peaks according to Eq. (28). We further observe that the maximum position in Figs. 4(a) and 4(b) as a function of $\sigma$ depends on the value of $E_{1}$. For the resonant system, the maxima are located around the same value, while for the off-resonant system, the maximum for the $E_{1}=0.9375\;\text{eV}$ curve is reached earlier than the other curves.

Disorder-averaged relaxation rate. The disorder-averaged relaxation rate is depicted in Figs. 4(c) and 4(d) as a function of $\sigma$ for different molecule numbers. In Fig. 4(c) for the resonant case, we observe a turnover which is a consequence of the turnovers of the energy-resolved relaxation rates in Fig. 4(a). According to the exact expression in Eq. (31), the disorder-averaged relaxation rate scales as

$\displaystyle\overline{\gamma}\propto\begin{cases}\frac{\sigma}{N}&\sigma\ll g\sqrt{N}\\ \frac{1}{\sigma}&\sigma\gg g\sqrt{N}\end{cases},$

(33)

in the two limiting cases, which is similar to Eq. (32) except for a factor $N$ for large disorder and can be clearly recognized in Fig. 4(c). For small $\sigma$, the width of the Lorentzian distribution increases with $\sigma$, which results in an increasing overlap with the cavity LDOS, whose spectral width also increases according to the imaginary part of the eigenvalues in Fig. 2. For large $\sigma$, the cavity LDOS $\nu_{C}(E_{1})$ becomes very narrow as discussed in Sec. III.1 and shown in Fig. 3(d), while the Lorentz distribution broadens and scales as $P(E_{1})\propto 1/\sigma$ for $E_{1}\approx E_{M}$. As a defining property, the integral of $\nu_{C}(E_{1})$ over all energies equals one, such that the overall relaxation rate scales as in Eq. (33).

The off-resonant case depicted in Fig. 4(d) exhibits two maxima for $N=100$ and $N=1000$. This is a consequence of the energy-resolved relaxation rate, whose peak position depends on the donor energy as can be seen in Fig. 4(b). Overall, the exact shape of the disorder-averaged dissipation rate sensitively depends on the shape of $\nu_{C}(E_{1})$ and thus on the system parameters.

Molecule-number dependence. The energy-resolved relaxation rate is depicted in Figs. 5(a) and 5(b) as a function of molecule number in the underdamped and overdamped regimes for the resonant system $E_{M}=E_{C}$. In both panels we observe a similar behavior. Interestingly, the rates exhibit a turnover in both regimes. The maximum position sensitively depends on the donor energy $E_{1}$, where the maximum for small $\left|E_{1}-E_{C}\right|$ is reached earlier than for large $\left|E_{1}-E_{C}\right|$. The limiting cases of $\gamma(E_{j})\propto N$ for small $N$ and $\gamma(E_{j})\propto 1/N$ for large $N$ can be directly inferred from Eq. (IV.1). For large $N$, the relaxation rate $\gamma(E_{1})$ becomes independent off $E_{1}$ as this limit is equivalent to a rescaling of the system energies $E_{1}/g\sqrt{N}\rightarrow 0$.

The disorder-averaged relaxation rate $\overline{\gamma}$ is depicted in Figs. 5(c) and 5(d). For small $N$, the rate $\overline{\gamma}$ is approximately constant in agreement with Eq. (33). For large $N$ (i.e., $g\sqrt{N}\gg\sigma$), the factor $1/N$ can be explained by the shape of $\nu_{C}(\omega)$ discussed in Sec. III.1: for the resonant system, the polariton peaks in Fig. 3(a) are located around $E_{peak}=E_{M}\pm g\sqrt{N}$ and decay as $\nu_{C}(E)\propto 1/(E-E_{peak})^{2}$ with increasing distance from the peak center. When the donor energy is distributed around $E_{M}$ for small $\sigma/g\sqrt{N}$, then the cavity mode overlap contributes the extra factor $1/N$.

Finite-size simulation. In Figs. 4 and 5, the analytical solutions are compared with finite-size simulations, which are represented as symbols. Overall, both calculations agree very well with each other, except for deviations for very small $\sigma$ in Figs. 4(a) and 4(b). The deviations occur for donor energies $E_{1}$ away from the center of the Lorenz distribution $E_{M}$. For these energies, the density of states is very low, such that the continuum limit, which has been applied in the derivation, is not justified. This problem is even more prominent when the donor energy $E_{1}$ is close to $E_{C}$ as can be seen for $E_{1}=1.125\;\text{eV}$ and $2.0\;\text{eV}$ in Figs. 4(b), as the level repulsion leads to a depletion of the cavity LDOS in this region. However, these deviations are not visible in the disorder-averaged rate as the energies in the sparse spectral regime hardly contribute to the disorder average. In Fig. 5 we clearly observe the improvement of the analytical calculation when approaching the thermodynamic limit $N\rightarrow\infty$.

To calculate the transport dynamics, we couple the Hamiltonian in Eq. (1) to an additional reservoir at the acceptor $j=N$ as depicted in Fig. 1(c). The amended Hamiltonian is then given by

	$\displaystyle H_{tr}$	$\displaystyle=$	$\displaystyle H+\hat{H}_{R}+\hat{H}_{AR},$
	$\displaystyle\hat{H}_{R}$	$\displaystyle=$	$\displaystyle\sum_{k}\epsilon_{k}\hat{c}_{k}^{\dagger}\hat{c}_{k},$
	$\displaystyle\hat{H}_{AR}$	$\displaystyle=$	$\displaystyle\sum_{k=1}^{N_{R}}g_{R}\hat{B}_{N}^{\dagger}\hat{c}_{k}+\text{H.c.},$		(34)

where the reservoir is described by the operators $\hat{c}_{k}$, which are coupled with strengths $g_{R}$ to the acceptor. The transport behavior of this and similar Hamiltonians has been numerically investigated in Refs. [32, 33, 34, 54, 35].

Initially, the system is excited at the donor $j=1$ such that the relevant Green’s functions are

	$\displaystyle G_{C,1}(z)$	$\displaystyle=$	$\displaystyle-i\frac{g}{\left(z+iE_{C}(z)\right)\left(z+iE_{1}\right)}$
		$\displaystyle+$	$\displaystyle i\frac{g^{3}}{\left(z+iE_{C}(z)\right)^{2}\left(z+iE_{N}(z)\right)\left(z+iE_{1}\right)},$
	$\displaystyle G_{j\neq N,1}(z)$	$\displaystyle=$	$\displaystyle\frac{\delta_{j,1}}{z+iE_{1}}$
		$\displaystyle-$	$\displaystyle\frac{g^{2}}{\left(z+iE_{j}\right)\left(z+iE_{C}(z)\right)\left(z+iE_{1}\right)}$
			$\displaystyle+\frac{g^{4}}{\left(z+iE_{j}\right)\left(z+iE_{C}(z)\right)^{2}\left(z+iE_{N}(z)\right)\left(z+iE_{1}\right)},$
	$\displaystyle G_{N,1}(z)$	$\displaystyle=$	$\displaystyle-\frac{g^{2}}{\left(z+iE_{N}(z)\right)\left(z+iE_{C}(z)\right)\left(z+iE_{1}\right)},$

where the two auxiliary functions are defined by

	$\displaystyle E_{C}(z)$	$\displaystyle=$	$\displaystyle E_{C}+\sum_{j=1}^{N-1}\frac{g^{2}}{z+iE_{j}},$		(36)
	$\displaystyle E_{N}(z)$	$\displaystyle=$	$\displaystyle E_{N}+\sum_{k}\frac{g_{R}^{2}}{z+iE_{k}}+\frac{g^{2}}{z+iE_{C}(z)}.$		(37)

For an appropriate multiplication with the factors $(z+iE_{j})$ and $(z+iE_{k})$, the denominators of the respective last terms in Eq. (V.1) define the polynomial

	$\displaystyle\mathcal{P}(z)$	$\displaystyle=$	$\displaystyle\left(z+iE_{N}(z)\right)\left(z+iE_{C}(z)\right)$		(38)
		$\displaystyle\times$	$\displaystyle\prod_{j}^{N}(z+iE_{j})\prod_{k}^{N_{R}}(z+iE_{k}),$

which is equivalent to the characteristic polynomial of the Hamiltonian in Eq. (34) (with a replacement of $z\rightarrow-iE$). It is not possible to find all roots of this polynomial, which has order $N+N_{R}+1$. As in Sec. II.4, we can take the thermodynamic limit in the second term of Eq. (37). Under the assumption that $E_{k}$ is distributed according to the Lorentz distribution, we can simply replace $E_{k}\rightarrow E_{R}+\Sigma$, where $E_{R}$ denotes the center of the reservoir distribution and $\Sigma$ is its width. Yet, as the transport rate is determined by the microscopic dynamics of the molecular excitations, the disorder average is not applied to the auxiliary function $E_{C}(z)$. After the disorder average of the reservoir states, the effective characteristic polynomial reads

	$\displaystyle\mathcal{P}(z)$	$\displaystyle=$	$\displaystyle\left(z+iE_{C}(z)\right)\left(z+iE_{N}(z)\right)(z+iE_{R}+\Sigma)$		(39)
		$\displaystyle\times$	$\displaystyle\prod_{j}^{N-1}(z+iE_{j}),$

whose $N+2$ complex roots $z_{\mu}$ determine the inverse Laplace transformation in Eq. (9).