# Quantum dynamics of molecular ensembles coupled with quantum light: Counter-
rotating interactions as an essential component

Yi-Ting Chuang  Department of Chemistry, National Taiwan University, Taipei
10617, Taiwan  Institute of Atomic and Molecular Sciences, Academia Sinica,
Taipei 10617, Taiwan  Liang-Yan Hsu<EMAIL_ADDRESS>Institute of
Atomic and Molecular Sciences, Academia Sinica, Taipei 10617, Taiwan
Department of Chemistry, National Taiwan University, Taipei 10617, Taiwan
National Center for Theoretical Sciences, Taipei 10617, Taiwan

###### Abstract

The rotating-wave approximation to light-matter interactions is widely used in
the quantum electrodynamics Hamiltonian; however, its validity has long been a
matter of debate. In this article, we explore the impact of the rotating-wave
approximation on the quantum dynamics of multiple molecules in complex
dielectric environments within the framework of macroscopic quantum
electrodynamics. In general, we find that the energy shifts of the molecules
and the inter-molecule dipole-dipole interaction obtained in the weak coupling
regime are correct only when the counter-rotating interactions are considered.
Moreover, under the rotating-wave approximation, the energy shifts of the
ground-state molecules and a portion of the inter-molecule interaction are
discarded. Notably, in the near-field zone (short inter-molecular distance),
the reduction of inter-molecule interaction can reach up to 50 percent. We
also conduct a case study on the population dynamics of a pair of identical
molecules above a plasmonic surface. Through analytical and numerical
analysis, it is revealed that the rotating-wave approximation can profoundly
affect the dynamics of the molecules in both strong and weak coupling regimes,
emphasizing the need for careful consideration when making the rotating-wave
approximation in a multiple-molecule system coupled with quantum light.

††preprint: APS/123-QED

##  I Introduction

Over the past few decades, extensive experimental and theoretical research has
demonstrated that coupling molecules (or quantum emitters) with confined
electromagnetic fields in a well-designed photonic environment is a promising
avenue for modifying physical processes such as spontaneous emission [1, 2, 3,
4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16] and resonance energy transfer
[17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. Recently, there has been
a surge of interest in utilizing the concept of coupling molecules with
confined electromagnetic fields to manipulate chemical reactions [30, 31, 32,
33, 34, 35, 36]. The successful realization of such experiments has spurred
scientists to develop theoretical frameworks aimed at elucidating the
underlying physics behind observed phenomena and predicting novel effects.
Among these theoretical frameworks, cavity quantum electrodynamics (QED)
Hamiltonian are extensively used, especially the Tavis-Cummings model and its
extensions [37, 38, 39, 40, 41, 42, 43].

Qunatum electrodynamics serves as one of the most favorable foundations to
study the interactions between molecules and electromagnetic fields. Within
the framework of QED, one typically needs to make use of further
approximations to simplify the light-matter interactions, and the rotating-
wave approximation (RWA) [44, 45] is viewed as one of the most commonly
adopted approximations. To simply demonstrate the key concept of the RWA, we
consider a simple model consisting of a single molecule
($\omega_{\mathrm{m}}$) interacting with a single photonic mode
($\omega_{\mathrm{p}}$), i.e.,

$\displaystyle\hbar\omega_{\mathrm{m}}\hat{\sigma}^{(+)}\hat{\sigma}^{(-)}+\hbar\omega_{\mathrm{p}}\hat{a}^{\dagger}\hat{a}+\hbar
g\left[\hat{\sigma}^{(+)}+\hat{\sigma}^{(-)}\right]\left[\hat{a}^{\dagger}+\hat{a}\right].$

The spirit of the RWA lies in retaining the co-rotating interaction $\hbar
g\left[\hat{\sigma}^{(+)}\hat{a}+\hat{\sigma}^{(-)}\hat{a}^{\dagger}\right]$
while disregarding the counter-rotating interaction $\hbar
g\left[\hat{\sigma}^{(+)}\hat{a}^{\dagger}+\hat{\sigma}^{(-)}\hat{a}\right]$.
This choice is motivated by the fact that the co-rotating interaction
oscillates at a relatively low frequency
$\omega_{\mathrm{m}}-\omega_{\mathrm{p}}$ while the counter-rotating
interaction oscillates at a higher frequency
$\omega_{\mathrm{m}}+\omega_{\mathrm{p}}$. By neglecting the counter-rotating
interaction, the RWA simplifies the mathematical treatment of the system and
allows for an easier analysis. Typically, the RWA is considered a good
approximation when (i) the light-matter coupling strength $g$ is weak and (ii)
the detuning $\omega_{\mathrm{m}}-\omega_{\mathrm{p}}$ is small. However, when
molecules are coupled to multiple photonic modes (including off-resonant
modes), such as in free space (infinite photonic modes), condition (ii) seems
to be violated anyhow.

The validity of the RWA has long been a matter of debate. While there has been
extensive research on the impact of the RWA, it is worth noting that most of
these studies have primarily focused on scenarios involving only a single
photonic mode [46, 47, 48, 49, 50, 51], in free space [52, 53, 54, 55], or in
non-dispersive and non-absorbing (lossless) media [56]. Nevertheless, in
experimental setups, such as the Fabry-Pérot cavity, there exist infinite
numbers of photonic modes, including many off-resonant modes [57], that are
dressed by the surrounding photonic environments. The interactions between
molecules and these dressed photons (polaritons) are supposed to play a
crucial role in determining the chemical and physical properties of the
system. Therefore, restricting our analysis to a single photonic mode or the
free space scenario may result in misinterpretations of experimental
observations and impede our ability to accurately predict and quantify
polariton-coupled processes. In addition, it has been shown that plasmonic
materials can strongly affect physical processes, including spontaneous
emission [58, 16] and resonance energy transfer [59, 29], and it is natural to
consider dispersive and absorbing media in studying plasmonic effects. To
overcome the aforementioned issues, one needs advanced theoretical treatments
that incorporate the coupling of molecules with infinite photonic modes and
take into account the effects of (dispersive and absorbing) photonic
environments.

In this study, in order to demonstrate counter-rotating interactions as an
essential component, we investigate the quantum dynamics of multiple molecules
in dielectric environments using macroscopic quantum electrodynamics (MQED)
[60, 61, 62], which is an effective field theory for describing quantized
electromagnetic fields in any arbitrary inhomogeneous, dispersive, and
absorbing dielectric environment. We elucidate the importance of counter-
rotating interactions in the quantum dynamics of multiple molecules in both
strong and weak coupling regimes. In addition, we clarify that a reduction of
the resonant dipole-dipole interaction due to the RWA has been wrongly
interpreted as a non-negligible (quantum) effect in previous works [63, 64].
Note that our theoretical approach is not only restricted to the study of the
RWA, but can also be applied to investigate the combined effect of molecular
fluorescence and excitation energy transfer in any arbitrary photonic
environments, holding great promise for further advancements in areas such as
quantum optics, nanophotonics, and molecular engineering.

This article is organized as follows. In Sec. II, starting from the MQED
Hamiltonian (including without and with the use of the RWA), we derive the
dynamical equations that can describe the quantum dynamics of multiple
emitters in complex dielectric environments. In addition, the underlying
physics behind the dynamical equations is discussed. In Sec. III, we apply the
Markov approximation to the dynamical equations to study their behavior in the
weak coupling regime. In this regime, we obtain the energy shift, decay rate,
and inter-molecule dipole-dipole interaction, and compare our results to the
previous works. In Sec. IV, we apply our approach to investigate the role of
counter-rotating interactions in the quantum dynamics of a pair of identical
molecules above a plasmonic surface in both strong and weak coupling regimes.
In the last section, we provide a concise summary of this study.

##  II Theory

###  II.1 Hamiltonian

Considering a collection of two-level molecules coupled to polaritons (dressed
photons) in an arbitrary inhomogeneous, dispersive, and absorbing medium, the
total Hamiltonian $\hat{H}$ (without the RWA) under the electric-dipole
approximation in the multipolar coupling MQED [65, 45] can be expressed in
terms of the molecular Hamiltonian $\hat{H}_{\mathrm{M}}$, polaritonic
Hamiltonian $\hat{H}_{\mathrm{P}}$ and interaction Hamiltonian
$\hat{H}_{\mathrm{I}}$ as

$\displaystyle\hat{H}=\hat{H}_{\mathrm{M}}+\hat{H}_{\mathrm{P}}+\hat{H}_{\mathrm{I}},$
(1)

with

$\displaystyle\hat{H}_{\mathrm{M}}$
$\displaystyle=\sum_{\alpha}\hbar\omega_{\alpha}\hat{\sigma}^{(+)}_{\alpha}\hat{\sigma}^{(-)}_{\alpha},$
(2) $\displaystyle\hat{H}_{\mathrm{P}}$
$\displaystyle=\int\mathrm{d}\mathbf{r}\int_{0}^{\infty}\mathrm{d}\omega\,\hbar\omega\,\mathbf{\hat{f}}^{\dagger}(\bf{r},\omega)\cdot\mathbf{\hat{f}}(\bf{r},\omega),$
(3) $\displaystyle\hat{H}_{\mathrm{I}}$
$\displaystyle=-\sum_{\alpha}\hat{\bm{\mu}}_{\alpha}\cdot\hat{\mathbf{F}}(\mathbf{r}_{\alpha}).$
(4)

The molecular Hamiltonian $\hat{H}_{\mathrm{M}}$ in Eq. (2) describes the
total energy of the molecules, where $\omega_{\mathrm{\alpha}}$,
$\hat{\sigma}^{(+)}_{\alpha}$ and $\hat{\sigma}^{(-)}_{\alpha}$ denote the
electronic transition frequency, raising operator and lowering operator of the
$\alpha$-th molecule, respectively. The raising and lowering operators of the
$\alpha$-th molecule can be defined using the electronically ground state
$\ket{\mathrm{g}_{\alpha}}$ and excited state $\ket{\mathrm{e}_{\alpha}}$ of
$\alpha$ as follows:
$\hat{\sigma}^{(+)}_{\alpha}=\ket{\mathrm{e}_{\alpha}}\bra{\mathrm{g}_{\alpha}}$
and $\hat{\sigma}^{(-)}_{\alpha}$ =
$\ket{\mathrm{g}_{\alpha}}\bra{\mathrm{e}_{\alpha}}$. Note that we have
neglected the dipole self-interaction since it only contributes to a small
energy shift (free space Lamb shift) after renormalization [66]. The
polaritonic Hamiltonian $\hat{H}_{\mathrm{P}}$ in Eq. (3) describes the energy
of polaritons in the dielectric environment, where
$\mathbf{\hat{f}}^{\dagger}(\mathbf{r},\omega)$ and
$\mathbf{\hat{f}}(\mathbf{r},\omega)$ are the creation and annihilation
operators of the bosonic vector fields that satisfy the commutation relations:

$\displaystyle\left[\hat{f}_{k}\left(\mathbf{r},\omega\right),\hat{f}^{\dagger}_{k^{\prime}}\left(\mathbf{r}^{\prime},\omega^{\prime}\right)\right]$
$\displaystyle=\delta_{kk^{\prime}}\delta\left(\mathbf{r}-\mathbf{r}^{\prime}\right)\delta\left(\omega-\omega^{\prime}\right),$
$\displaystyle\left[\hat{f}_{k}\left(\mathbf{r},\omega\right),\hat{f}_{k^{\prime}}\left(\mathbf{r}^{\prime},\omega^{\prime}\right)\right]$
$\displaystyle=0.$

The interaction Hamiltonian $\hat{H}_{\mathrm{I}}$ in Eq. (4) describes the
couplings between the molecules and polaritons including counter-rotating
interactions, i.e., no RWA , where $\hat{\bm{\mu}}_{\alpha}$ and
$\hat{\mathbf{F}}(\mathbf{r}_{\alpha})$ are the transition dipole operator of
$\alpha$ and the field operator, respectively. The transition dipole operator
$\hat{\bm{\mu}}_{\alpha}$ can be written in terms of
$\hat{\sigma}^{(+)}_{\alpha}$ and $\hat{\sigma}^{(-)}_{\alpha}$ as follows,

$\displaystyle\hat{\bm{\mu}}_{\alpha}=\bm{\mu}^{\mathrm{eg}}_{\alpha}\hat{\sigma}^{(+)}_{\alpha}+\bm{\mu}^{\mathrm{ge}}_{\alpha}\hat{\sigma}^{(-)}_{\alpha},$
(5)

where $\bm{\mu}^{\mathrm{eg}}_{\alpha}=(\bm{\mu}^{\mathrm{ge}}_{\alpha})^{*}$
is the electronic transition dipole moment of $\alpha$. The field operator is
defined as:

$\displaystyle\hat{\mathbf{F}}(\mathbf{r}_{\alpha})=\hat{\mathbf{F}}^{(+)}(\mathbf{r}_{\alpha})+\mathrm{H.c.},$
(6)

with

$\displaystyle\hat{\mathbf{F}}^{(+)}(\mathbf{r}_{\alpha})=\int_{0}^{\infty}\differential{\omega}\int\differential{\mathbf{r}}\overline{\overline{\mathbf{g}}}(\mathbf{r}_{\alpha},\mathbf{r},\omega)\cdot\hat{\mathbf{f}}(\mathbf{r},\omega),$
(7)
$\displaystyle\overline{\overline{\mathbf{g}}}(\mathbf{r}_{\alpha},\mathbf{r},\omega)=i\sqrt{\frac{\hbar}{\pi\varepsilon_{0}}}\frac{\omega^{2}}{c^{2}}\sqrt{\mathrm{Im}\left[\varepsilon_{\mathrm{r}}(\mathbf{r},\omega)\right]}\,\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r},\omega),$
(8)

where $\varepsilon_{0}$, $\varepsilon_{\mathrm{r}}(\mathbf{r},\omega)$, and
$c$ are the permittivity of free space, relative permittivity, and speed of
light in vacuum, respectively.
$\overline{\overline{\mathbf{g}}}(\mathbf{r}_{\alpha},\mathbf{r},\omega)$ is
an auxiliary tensor defined in terms of the dyadic Green’s function
$\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r},\omega)$ that
satisfies macroscopic Maxwell’s equations, i.e.,

$\left[\frac{\omega^{2}}{c^{2}}\varepsilon_{\mathrm{r}}(\mathbf{r}_{\alpha},\omega)-\nabla\times\nabla\times\right]\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r},\omega)=-\mathbf{\overline{\overline{I}}}_{3}\delta(\mathbf{r}_{\alpha}-\mathbf{r}),$
(9)

where $\mathbf{\overline{\overline{I}}}_{3}$ and
$\delta(\mathbf{r}_{\alpha}-\mathbf{r})$ are the $3\times 3$ identity matrix
and three-dimensional delta function, respectively. Note that the dyadic
Green’s function can be further decomposed as
$\overline{\overline{\mathbf{G}}}(\mathbf{r},\mathbf{r}^{\prime},\omega)=\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)+\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r},\mathbf{r}^{\prime},\omega)$,
where
$\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)$
represents the free-space dyadic Green’s function in the absence of the
dielectric bodies and
$\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r},\mathbf{r}^{\prime},\omega)$
represents the scattering dyadic Green’s function originating from the
presence of the dielectric bodies.

###  II.2 State Vector

To adequately include the effect of the counter-rotating interactions, we
extend the Wigner-Weisskopf wave function ansatz [67] and use the following
state vector to describe the total system [54, 68],

$\displaystyle\ket{\Psi(t)}=$
$\displaystyle\sum_{\alpha}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)e^{-iW^{\mathrm{E_{\alpha}},\left\\{0\right\\}}t}\ket{\mathrm{E_{\alpha}}}\ket{\left\\{0\right\\}}+\sum_{k=1}^{3}\int\differential{\mathrm{\bf{r}}}\int_{0}^{\infty}\differential{\omega}C^{\mathrm{G},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t)e^{-iW^{\mathrm{G},\left\\{1\right\\}}(\omega)t}\ket{\mathrm{G}}\ket{\left\\{1_{k}(\mathrm{\bf{r}},\omega)\right\\}}$
$\displaystyle\quad+\sum_{\alpha}\sum_{\beta>\alpha}\sum_{k=1}^{3}\int\differential{\mathbf{r}}\int_{0}^{\infty}\differential{\omega}C^{\mathrm{E_{\alpha\beta}},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t)e^{-iW^{\mathrm{E_{\alpha\beta}},\left\\{1\right\\}}(\omega)t}\ket{\mathrm{E_{\alpha\beta}}}\ket{\left\\{1_{k}(\mathrm{\bf{r}},\omega)\right\\}},$
(10)

with

$\displaystyle W^{\mathrm{E_{\alpha}},\\{0\\}}=\omega_{\alpha},$ (11a)
$\displaystyle W^{\mathrm{G},\\{1\\}}(\omega)=\omega,$ (11b) $\displaystyle
W^{\mathrm{E_{\alpha\beta}},\\{1\\}}(\omega)=\omega+\omega_{\alpha}+\omega_{\beta}.$
(11c)

The state vector includes the molecular (electronic) and photonic degrees of
freedom in the entire system. For the molecular part, the ket state
$\ket{\mathrm{G}}$ denotes that all the molecules are in their electronically
ground states, i.e.,
$\ket{\mathrm{G}_{\alpha}}=\ket{\mathrm{g}_{1}}\ket{\mathrm{g}_{2}}\dots\ket{\mathrm{g}_{\alpha}}\dots\ket{\mathrm{g}_{N}}$;
the ket state $\ket{\mathrm{E}_{\alpha}}$ denotes that $\alpha$ is in its
electronically excited state while the other molecules are in their
electronically ground states, i.e.,
$\ket{\mathrm{E}_{\alpha}}=\hat{\sigma}^{(+)}_{\alpha}\ket{\mathrm{G}}$; the
ket state $\ket{\mathrm{E}_{\alpha\beta}}$ denotes that both $\alpha$ and
$\beta$ are in their electronically excited states while the other molecules
are in their electronically ground states, i.e.,
$\ket{\mathrm{E}_{\alpha\beta}}=\hat{\sigma}^{(+)}_{\alpha}\hat{\sigma}^{(+)}_{\beta}\ket{\mathrm{G}}$.
For the photonic part, the ket state $\lvert\\{0\\}\rangle$ denotes the zero-
polariton state, and the ket state
$\lvert\\{1_{k}(\mathbf{r},\omega)\\}\rangle$ denotes the single-polariton
state, i.e.,
$\lvert\\{1_{k}(\mathbf{r},\omega)\\}\rangle=\hat{f}^{\dagger}_{k}\left(\mathbf{r},\omega\right)\lvert\\{0\\}\rangle$.
$C^{\mathrm{E_{\alpha},\left\\{0\right\\}}}(t)$,
$C^{\mathrm{G},\\{1_{k}\\}}(\mathbf{r},\omega,t)$ and
$C^{\mathrm{E_{\alpha\beta},\\{1_{k}\\}}}(\mathbf{r},\omega,t)$ are the
probability amplitudes of $\lvert\Psi(t)\rangle$ for the states
$\ket{\mathrm{E}_{\alpha}}\lvert\\{0\\}\rangle$,
$\ket{\mathrm{G}}\lvert\\{1_{k}(\mathbf{r},\omega)\\}\rangle$ and
$\ket{\mathrm{E}_{\alpha\beta}}\lvert\\{1_{k}(\mathbf{r},\omega)\\}\rangle$,
respectively; $\hbar W^{\mathrm{E}_{\alpha},\\{0\\}}$, $\hbar
W^{\mathrm{G},\\{1\\}}(\omega)$ and $\hbar
W^{\mathrm{E}_{\alpha\beta},\\{1\\}}(\omega)$ are the total energies of the
states $\ket{\mathrm{E}_{\alpha}}\lvert\\{0\\}\rangle$,
$\ket{\mathrm{G}}\lvert\\{1_{k}(\mathbf{r},\omega)\\}\rangle$ and
$\ket{\mathrm{E}_{\alpha\beta}}\lvert\\{1_{k}(\mathbf{r},\omega)\\}\rangle$,
respectively.

###  II.3 Quantum Dynamics: Equation of Motion

Figure 1: Schematic illustration of the quantum dynamics in Eqs. (12).

The quantum dynamics of the entire system without adopting the application of
the RWA can be obtained by solving the time-dependent Schrödinger equation
$i\hbar\partial\ket{\Psi(t)}/\partial t=\hat{H}\ket{\Psi(t)}$. After some
algebra, we obtain the equation of motion as follows (see Appendix A for the
details),

$\displaystyle\derivative{t}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=$
$\displaystyle-\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}\right]e^{-i\left(\omega-\omega_{\alpha}\right)t}e^{-i\left(\omega_{\alpha}-\omega\right)t^{\prime}}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}\right]e^{-i\left(\omega-\omega_{\alpha}\right)t}e^{-i\left(\omega_{\beta}-\omega\right)t^{\prime}}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\beta},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{eg}}_{\beta}\right]e^{-i\left(\omega+\omega_{\beta}\right)t}e^{-i\left(-\omega_{\beta}-\omega\right)t^{\prime}}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\beta},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{eg}}_{\alpha}\right]e^{-i\left(\omega+\omega_{\beta}\right)t}e^{-i\left(-\omega_{\alpha}-\omega\right)t^{\prime}}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t^{\prime}).$
(12)

The four terms in Eq. (12) corresponds to distinct physical processes:

1\. Spontaneous emission and reabsorption assisted by co-rotating
interactions: The first term on the right-hand side of Eq. (12) corresponds to
the physical process depicted by the blue arrow in FIG. 1. This process
involves an initially excited molecule $\alpha$ emitting a photon and
returning to its ground state. The emitted photon is then absorbed by $\alpha$
again, causing it to transition back to its excited state. In fact, this
process exactly corresponds to spontaneous emission, which has been derived
from previous studies [69, 70].

2\. Excitation energy transfer assisted by co-rotating interactions: The
second term on the right-hand side of Eq. (12) corresponds to the physical
process depicted by the red arrow in FIG. 1. In this process, an initially
excited molecule $\beta$ emits a photon, transitioning back to its ground
state. The emitted photon is subsequently absorbed by molecule $\alpha$,
resulting in the excitation of molecule $\alpha$. This process can be regarded
as excitation energy transfer assisted by co-rotating interactions [71].

3\. Virtual photon emission and reabsorption assisted by counter-rotating
interactions: The third term on the right-hand side of Eq. (12) corresponds to
the physical process depicted by the purple arrow in FIG. 1. It involves an
initially excited molecule $\alpha$ transitioning to a state where both
molecules $\alpha$ and $\beta$ are excited, accompanied by the emission of a
single polariton. Subsequently, the system transitions back to a state where
only $\alpha$ remains excited. This process is facilitated by the counter-
rotating interactions, which involve terms that describe the simultaneous
creation or annihilation of both photons and molecular excitations.

4\. Excitation energy transfer assisted by counter-rotating interactions: The
fourth term on the right-hand side of Eq. (12) corresponds to the physical
process depicted by the green arrow in FIG. 1. In this process, an initially
excited molecule $\beta$ transitions to a state where both molecules $\alpha$
and $\beta$ are excited, with the presence of a single polariton.
Subsequently, the system transitions back to a state where only $\alpha$
remains excited. Similar to the third term, this process is facilitated by the
counter-rotating interactions and can be considered as another pathway of
excitation energy transfer.

###  II.4 Rotating-Wave Approximation

Under the RWA, we neglect the so-called counter-rotating terms in the
interaction Hamiltonian in Eq. (4). These counter-rotating terms involve
processes where photons and molecular excitation are simultaneously created or
annihilated. By discarding these terms, the interaction Hamiltonian is reduced
to

$\displaystyle\hat{H}^{\mathrm{RWA}}_{\mathrm{I}}$
$\displaystyle=-\sum_{\alpha=1}^{N}\left[\hat{\sigma}^{(+)}_{\alpha}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\hat{\mathbf{F}}^{(+)}(\mathbf{r}_{\alpha})+\mathrm{H.c.}\right],$
(13)

and the total Hamiltonian under the RWA is defined as

$\displaystyle\hat{H}^{\mathrm{RWA}}=\hat{H}_{\mathrm{M}}+\hat{H}_{\mathrm{P}}+\hat{H}^{\mathrm{RWA}}_{\mathrm{I}}.$
(14)

Note that for the total Hamiltonian under the RWA, we do not need to consider
the ket state
$\ket{\mathrm{E}_{\alpha\beta}}\lvert\\{1_{k}(\mathbf{r},\omega)\\}\rangle$
since it is not accessible without the presence of the counter-rotating
interactions (assume that there is no polariton at $t=0$); therefore, we can
simplify the state vector as:

$\displaystyle\ket{\Psi(t)}^{\mathrm{RWA}}=\sum_{\alpha}\tilde{C}^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)e^{-iW^{\mathrm{E_{\alpha}},\left\\{0\right\\}}t}\ket{\mathrm{E_{\alpha}}}\ket{\left\\{0\right\\}}+\sum_{k=1}^{3}\int\differential{\mathrm{\bf{r}}}\int_{0}^{\infty}\differential{\omega}\tilde{C}^{\mathrm{G},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t)e^{-iW^{\mathrm{G},\left\\{1\right\\}}(\omega)t}\ket{\mathrm{G}}\ket{\left\\{1_{k}(\mathrm{\bf{r}},\omega)\right\\}},$
(15)

where $\tilde{C}^{\mathrm{E_{\alpha},\left\\{0\right\\}}}(t)$ and
$\tilde{C}^{\mathrm{G},\\{1_{k}\\}}(\mathbf{r},\omega,t)$ are the probability
amplitudes of $\lvert\Psi(t)\rangle^{\mathrm{RWA}}$ for the states
$\ket{\mathrm{E}_{\alpha}}\lvert\\{0\\}\rangle$ and
$\ket{\mathrm{G}}\lvert\\{1_{k}(\mathbf{r},\omega)\\}\rangle$, respectively.

Solving the time-dependent Schrödinger equation
$i\hbar\partial\ket{\Psi(t)}^{\mathrm{RWA}}/\partial
t=\hat{H}^{\mathrm{RWA}}\ket{\Psi(t)}$, we obtain the equation of motion for
the quantum dynamics under the RWA as follows,

$\displaystyle\derivative{t}\tilde{C}^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=$
$\displaystyle-\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}\right]e^{-i\left(\omega-\omega_{\alpha}\right)t}e^{-i\left(\omega_{\alpha}-\omega\right)t^{\prime}}\tilde{C}^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}\right]e^{-i\left(\omega-\omega_{\alpha}\right)t}e^{-i\left(\omega_{\beta}-\omega\right)t^{\prime}}\tilde{C}^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t^{\prime}).$
(16)

Regarding Eq. (16), the first and second terms on the right-hand side
correspond to the physical processes represented by the blue and red arrows in
FIG. 1, respectively. The physical processes depicted by the purple and green
arrows are not included in Eq. (16) due to the absence of the counter-rotating
terms in the interaction Hamiltonian in Eq. (13).

##  III Weak Coupling Regime: Energy Shift, Decay Rate and Dipole-Dipole
Interaction

In this section, we will show counter-rotating interactions as an essential
component in describing dipole-dipole interactions, even in the weak light-
matter coupling regime. To clearly demonstrate the importance of counter-
rotating interactions, we further make the Markov approximation to Eq. (12)
and Eq. (16), which correspond to the equations of motion without and with
adopting the RWA, respectively, and then compare their physical meaning.

First of all, under weak light-matter interactions, we make the Markov
approximation to Eq. (12), and this equation can be simplified as (see
Appendix B):

$\displaystyle\derivative{t}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=$
$\displaystyle\quad-\frac{i}{\hbar}\left\\{\left[\Delta_{\mathrm{e}_{\alpha}}+\sum_{\beta\neq\alpha}\Delta_{\mathrm{g}_{\beta}}\right]-i\hbar\frac{\Gamma_{\alpha}}{2}\right\\}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)$
$\displaystyle\quad-\frac{i}{\hbar}\sum_{\beta\neq\alpha}\mathrm{V}_{\mathrm{DDI,\alpha\beta}}\,e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t),$
(17)

where $\Delta_{\mathrm{e(g)}}$ represents the energy shift of the excited
(ground) state of $\alpha$, $\Gamma_{\alpha}$ is the decay rate of the excited
state of $\alpha$, and $\mathrm{V}_{\mathrm{DDI,\alpha\beta}}$ denotes the
effective dipole-dipole interaction between $\alpha$ and $\beta$. The energy
shift $\Delta_{\mathrm{e(g)}}$ comprises two contributions: the free-space
Lamb shift $\Delta^{0}_{\mathrm{e(g)}}$ and Casmir-Polder potential
$\Delta^{\mathrm{Sc}}_{\mathrm{e(g)}}$, i.e.,

$\displaystyle\Delta_{\mathrm{e(g)}_{\alpha}}=\Delta^{0}_{\mathrm{e(g)}_{\alpha}}+\Delta^{\mathrm{Sc}}_{\mathrm{e(g)}_{\alpha}},$
(18)

with

$\displaystyle\Delta^{0}_{\mathrm{e(g)}_{\alpha}}=-\mathcal{P}\int_{0}^{\infty}\differential{\omega}\frac{\omega^{2}}{\pi\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}}{\omega-(+)\omega_{\alpha}},$
(19)
$\displaystyle\Delta^{\mathrm{Sc}}_{\mathrm{e(g)}_{\alpha}}=-\mathcal{P}\int_{0}^{\infty}\differential{\omega}\frac{\omega^{2}}{\pi\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}}{\omega-(+)\omega_{\alpha}},$
(20)

where $\mathcal{P}$ denotes the principal value. The energy shifts are
identical to those derived from perturbation theory [65, 72]. Note that the
free-space Lamb shift $\Delta^{0}_{\mathrm{e(g)}_{\alpha}}$ is divergent in
this context and requires proper treatment, such as renormalization. However,
the contribution of the free-space Lamb shift is typically small after
renormalization [66]. Therefore, for the sake of simplicity, this term can be
neglected in practical calculations (or considered as being included in the
transition energy of the molecule) [73, 69]. The decay rate is expressed as

$\displaystyle\Gamma_{\alpha}=\frac{2\omega_{\alpha}^{2}}{\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega_{\alpha})\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha},$
(21)

which is consistent with the spontaneous emission rate of a molecule in a
medium derived from Fermi’s golden rule [74]. The dipole-dipole interaction
$\mathrm{V}_{\mathrm{DDI,\alpha\beta}}$ can be divided into two components,
i.e.,

$\mathrm{V}_{\mathrm{DDI,\alpha\beta}}=\mathrm{V}_{\mathrm{RDDI,\alpha\beta}}+\mathrm{V}_{\mathrm{ORC,\alpha\beta}},$
(22)

with

$\mathrm{V}_{\mathrm{RDDI,\alpha\beta}}=\frac{-\omega_{\beta}^{2}}{\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega_{\beta})\cdot\bm{\mu}^{\mathrm{ge}}_{\beta},$
(23)

$\displaystyle\mathrm{V}_{\mathrm{ORC,\alpha\beta}}$
$\displaystyle=\int_{0}^{\infty}\differential{\omega}\frac{\omega^{2}}{\pi\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}}{\omega+\omega_{\beta}}$
$\displaystyle-\int_{0}^{\infty}\differential{\omega}\frac{\omega^{2}}{\pi\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}}{\omega+\omega_{\alpha}}.$
(24)

$\mathrm{V}_{\mathrm{RDDI,\alpha\beta}}$ represents the dipole-dipole
interaction between a pair of on-resonant ($\omega_{\alpha}=\omega_{\beta}$)
molecules, which is identical to the resonant dipole-dipole interaction in the
presence of dielectric bodies obtained using perturbation theory [75]. On the
other hand, $\mathrm{V}_{\mathrm{ORC,\alpha\beta}}$ can be regarded as a
correction to the dipole-dipole interaction when the molecules are off-
resonant ($\omega_{\alpha}\neq\omega_{\beta}$) since it is only non-zero when
$\omega_{\alpha}\neq\omega_{\beta}$. The calculation of
$\mathrm{V}_{\mathrm{ORC,\alpha\beta}}$ is complicated as it requires the
evaluation of integrals of the dyadic Green’s function over all the positive
frequencies. A useful technique is to transform the integral to the imaginary
axis, where the Dyadic Green’s function is much better behaved [63]. Using
this technique, one can derive the explicit form of the free-space off-
resonance correction $\mathrm{V}^{0}_{\mathrm{ORC,\alpha\beta}}$ [replace
$\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)$
with
$\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)$
in Eq. (24)], as shown in Appendix C. In this work, we will not conduct a
deeper discussion on the effect of the off-resonance correction; however, we
would like to emphasize that this correction may play an important role when
the frequency detuning between the molecular transitions is large.

Furthermore, we move on to the equation of motion with the RWA. Similarly, we
make the Markov approximation to Eq. (16), and this equation can be simplified
as

$\displaystyle\derivative{t}\tilde{C}^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=$
$\displaystyle\quad-\frac{i}{\hbar}\left[\Delta_{\mathrm{e}_{\alpha}}-i\hbar\frac{\Gamma_{\alpha}}{2}\right]\tilde{C}^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)$
$\displaystyle\quad-\frac{i}{\hbar}\sum_{\beta\neq\alpha}\tilde{\mathrm{V}}_{\mathrm{DDI,\alpha\beta}}\,e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}\tilde{C}^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t),$
(25)

where

$\tilde{\mathrm{V}}_{\mathrm{DDI,\alpha\beta}}=\mathrm{V}_{\mathrm{RDDI,\alpha\beta}}+\mathrm{V}_{\mathrm{QC,\alpha\beta}},$
(26)
$\mathrm{V}_{\mathrm{QC,\alpha\beta}}=\int_{0}^{\infty}\differential{\omega}\frac{\omega^{2}}{\pi\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}}{\omega+\omega_{\beta}}.$
(27)

Comparing Eq. (17) and Eq. (25), it is obvious that the quantum dynamics
without and with making the RWA under the Markov approximation exhibits two
major differences. The first major difference: the energy shifts of the
ground-state molecules ($\sum_{\beta\neq\alpha}\Delta_{\mathrm{g}_{\beta}}$)
are absent in the dynamical equation under the RWA. The second major
difference: from Eq. (22) and Eq. (26), one can find that the dipole-dipole
interactions between a pair of molecules exhibit disparities in the two
dynamical equations.

The second major difference naturally raises a question: which one is a
correct form of the dipole-dipole interaction, Eq. (22) or Eq. (26)? To answer
this question, we conduct the following discussion. First, when the molecules
are on-resonant ($\omega_{\alpha}=\omega_{\beta}$), the dipole-dipole
interaction within the RWA, i.e.,
$\tilde{\mathrm{V}}_{\mathrm{DDI,\alpha\beta}}$ in Eq. (25), is expressed as a
sum of the resonant dipole-dipole interaction
$\mathrm{V}_{\mathrm{RDDI,\alpha\beta}}$ [Eq. (23)] and an additional
correction term $\mathrm{V}_{\mathrm{QC,\alpha\beta}}$, which has been
reported in several previous studies [63, 28, 64]. It is worth noting that
this correction term is non-zero even when the molecules are on-resonant;
therefore, $\tilde{\mathrm{V}}_{\mathrm{DDI,\alpha\beta}}$ always deviates
from the resonant dipole-dipole interaction
$\mathrm{V}_{\mathrm{RDDI,\alpha\beta}}$. In some previous studies, this
deviation is considered an important effect that cannot be neglected in
describing inter-molecule interaction [63] and is regarded as a quantum
correction that cannot be obtained from classical electrodynamics [64].
However, we speculate that this so-called quantum correction term may only be
a product (an artifact) of the rotating-wave approximation since under the on-
resonance condition, this correction no longer exists when the counter-
rotating interactions are included. Second, when the molecules are off-
resonant ($\omega_{\alpha}\neq\omega_{\beta}$), the dipole-dipole interaction
without the RWA, i.e., $\mathrm{V}_{\mathrm{DDI,\alpha\beta}}$ in Eq. (22),
gives a reasonable correction $\mathrm{V}_{\mathrm{ORC,\alpha\beta}}$ because
the absolute value of $\mathrm{V}_{\mathrm{ORC,\alpha\beta}}$ is much smaller
than that of $\mathrm{V}_{\mathrm{QC,\alpha\beta}}$. Third, in the short-
distance (non-retarded) limit $\omega_{\beta(\alpha)}R/c\ll 1$, where
$\mathbf{r}_{\alpha}-\mathbf{r}_{\beta}\equiv R\mathbf{n}_{R}$, the dyadic
Green’s function becomes purely longitudinal [65], and the dipole-dipole
interactions can be approximated by their free-space and electrostatic
($\omega_{\alpha(\beta)}/c\rightarrow 0$) limits; therefore,

$\displaystyle\mathrm{V}_{\mathrm{DDI,\alpha\beta}}\approx
2\tilde{\mathrm{V}}_{\mathrm{DDI,\alpha\beta}}\approx\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}-3\left(\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathbf{n}_{R}\right)\left(\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathbf{n}_{R}\right)}{4\pi\varepsilon_{0}R^{3}}.$
(28)

Equation (28) clearly shows that $\mathrm{V}_{\mathrm{DDI,\alpha\beta}}$
converges to the conventional Coulomb dipole-dipole interaction, while
$\tilde{\mathrm{V}}_{\mathrm{DDI,\alpha\beta}}$ converges to only half of the
Coulomb dipole-dipole interaction. The omission of half of the Coulomb dipole-
dipole interaction within the framework of the RWA reinforces the significance
of incorporating the counter-rotating terms in the interaction Hamiltonian. In
addition to our study, in fact, the reduction of half of the dipole-dipole
interaction in free space due to the use of the RWA has also been reported
recently by Wubs et al [56]. Based on the above discussion, we can conclude
that counter-rotating interactions play a key role even in the weak light-
matter coupling regime.

In short, by examining the dynamical equation in the weak coupling regime [Eq.
(17)], we can obtain the energy shift, decay rate, and resonant dipole-dipole
interaction of molecules in a dielectric environment, and these results are
consistent with those derived from perturbation theory in previous works. This
agreement not only provides robust support for the validity of our theoretical
approach but also demonstrates counter-rotating interaction as an essential
component for dipole-dipole interactions.

##  IV Quantum Dynamics of A pair of Identical Molecules above a Plasmonic
Surface

Figure 2: Schematic illustration of a donor (A) and an acceptor (A) above a
plasmonic surface, where the donor-acceptor distance and the molecule-
dielectric distance are given by $d$ and $h$, respectively.

In this section, we numerically investigate the effect of counter-rotating
interactions on the quantum dynamics of multiple molecules in a complex
dielectric environment. For simplicity, we adopt a minimal model consisting of
a pair of identical molecules, i.e., a donor (D) and an acceptor (A), above a
plasmonic surface, with their transition dipole moments parallel to the normal
direction of the plasmonic surface, as depicted in FIG. (2). We denote the
distance between the donor and acceptor as $d$ and the distance between the
molecules and plasmonic surface as $h$. The plasmonic surface is modeled by
the following dielectric functions:

$\displaystyle\varepsilon_{\mathrm{r}}(\mathbf{r},\omega)=\begin{cases}1,&z>0,\\\
\varepsilon_{\mathrm{D}}(\omega),&z<0,\end{cases}$ (29)

where $\varepsilon_{\mathrm{D}}(\omega)=1-5/(\omega^{2}+0.1i\omega)$ is an
artificial Drude model. For the molecules, the transition frequency is given
by $\hbar\omega_{\mathrm{D}}=\hbar\omega_{\mathrm{A}}=3.525$ eV, which is in
resonance with the frequency of the surface plasmon polariton mode of the
chosen plasmonic surface, and the magnitude of the transition dipole moment is
given by
$\absolutevalue{\bm{\mu}_{\mathrm{D}}}=\absolutevalue{\bm{\mu}_{\mathrm{A}}}=10$
Debye. The dyadic Green’s functions of the system can be obtained through the
Fresnel method [76], and detailed information regarding the methodology can be
found in our previous works [77].

Figure 3: Excited-state population dynamics of a donor and an acceptor above a
plasmonic surface. (a) Population dynamics of the donor, (b) population
dynamics of the acceptor, and (c) total excited-state population dynamics when
$h=10$ nm and $d=4$ nm. (d) Population dynamics of the donor, (e) population
dynamics of the acceptor, and (f) total excited-state population dynamics when
$h=1$ nm and $d=1$ nm. FQD and MAQD denote full quantum dynamics and Markov-
approximated quantum dynamics, respectively.

In order to efficiently perform numerical calculations based on Eqs. (12),
(16), (17), and (25), we implement additional simplifications. For Eqs. (12)
and (16), we decompose the dyadic Green’s function as a sum of the free-space
contribution and scattering contribution, i.e.,
$\overline{\overline{\mathbf{G}}}(\mathbf{r},\mathbf{r}^{\prime},\omega)=\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)+\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r},\mathbf{r}^{\prime},\omega)$.
Since the coupling between the molecules and the free-space field is typically
weak, we can apply the Markov approximation only to the parts involving
$\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)$
and neglect the small free-space Lamb shift $\Delta^{0}_{\mathrm{e(g)}}$. This
simplification leads to our working equations, as shown in Eqs. (57) and (58).
These two equations will be referred to as the full quantum dynamics (FQD)
without and with the RWA in the subsequent discussion. For Eqs. (17) and (25),
we discard the small free-space Lamb shift $\Delta^{0}_{\mathrm{e(g)}}$,
resulting in our working equations as shown in Eqs. (59) and (60). These two
equations will be referred to as the Markov-approximated quantum dynamics
(MAQD) without and with the RWA in the subsequent discussion. The numerical
study will focus on how the counter-rotating interactions affects the quantum
dynamics in the weak and strong coupling conditions, where the (weak/strong)
coupling condition can be identified by the consistency between the population
dynamics obtained from the FQD and MAQD, e.g., the consistency of the
population dynamics of FQD and MAQD indicates the weak coupling condition [78,
79].

In FIG. 3(a)-(c), when $h=10$ nm and $d=4$ nm, the excited-state population of
the donor, i.e.,
$P^{\mathrm{E_{D},\\{0\\}}}(t)=\absolutevalue{C^{\mathrm{E_{D},\\{0\\}}}(t)}^{2}$,
the excited-state population of the acceptor, i.e.,
$P^{\mathrm{E_{A},\\{0\\}}}(t)=\absolutevalue{C^{\mathrm{E_{A},\\{0\\}}}(t)}^{2}$,
and the total excited-state population, i.e.,
$P^{\mathrm{E_{D},\\{0\\}}}(t)+P^{\mathrm{E_{A},\\{0\\}}}(t)$, obtained via
FQD and MAQD match each other; therefore, this circumstance can be identified
as a weak coupling condition. In this case, the individual population dynamics
of the donor $P^{\mathrm{E_{D},\\{0\\}}}(t)$ and acceptor
$P^{\mathrm{E_{A},\\{0\\}}}(t)$ calculated without and with the application of
the RWA are quite distinct from each other, as shown in FIG. 3(a) and (b).
More specifically, the oscillation frequency of the population dynamics
without the RWA is roughly two times the oscillation frequency of the
population dynamics with the RWA. Although the individual population dynamics
are different, it is surprising that the total population dynamics are still
identical in this case, as shown in Fig 3(c). This phenomenon can be further
confirmed by the analytical solution of MAQD (with the initial condition
$C^{\mathrm{E_{\alpha},\\{0\\}}}(t)=\delta_{\mathrm{\alpha D}}$), which can be
expressed as

$\displaystyle\begin{cases}P^{\mathrm{E_{D},\\{0\\}}}(t)=e^{-\Gamma
t}\left\\{\sinh^{2}\left[\mathrm{Im}(\mathrm{V}/\hbar)t\right]+\cos^{2}\left[\mathrm{Re}(\mathrm{V}/\hbar)t\right]\right\\},\\\
P^{\mathrm{E_{A},\\{0\\}}}(t)=e^{-\Gamma
t}\left\\{\sinh^{2}\left[\mathrm{Im}(\mathrm{V}/\hbar)t\right]+\sin^{2}\left[\mathrm{Re}(\mathrm{V}/\hbar)t\right]\right\\},\\\
P^{\mathrm{E_{D},\\{0\\}}}(t)+P^{\mathrm{E_{A},\\{0\\}}}(t)=e^{-\Gamma
t}\cosh\left[2\mathrm{Im}(\mathrm{V}/\hbar)t\right].\end{cases}$ (30)

In Eq. (30), $\Gamma$ and $\mathrm{V}$ are defined as: (i)
$\Gamma=\Gamma_{\mathrm{D}}=\Gamma_{\mathrm{A}}$ [Eq. 21] for both the
population dynamics without and with the RWA. (ii)
$\mathrm{V}=\mathrm{V}_{\mathrm{DDI,DA}}=\mathrm{V}_{\mathrm{DDI,AD}}$ [Eq.
(22)] for the population dynamics without the RWA, and
$\mathrm{V}=\tilde{\mathrm{V}}_{\mathrm{DDI,DA}}=\tilde{\mathrm{V}}_{\mathrm{DDI,AD}}$
[Eq. (26)] for the population dynamics with the RWA.

The analytical solution of the population dynamics in Eq. (30) provides plenty
of information. First, recall that $e^{-\Gamma t}$ is exactly the single-
molecule excited-state population, and $\cosh(x)$ is greater than or equal to
1. As a result, one can conclude that the inclusion of a second molecule can
slow down the decay of the total excited-state population, which is known as
the subradiance. Second, the total excited-state population only depends on
$\Gamma$ and $\mathrm{Im}(\mathrm{V})$, i.e., the imaginary part of the
dipole-dipole interaction. Thus, the total excited-state population is
unaffected by the RWA since
$\mathrm{Im}(\mathrm{V}_{\mathrm{DDI,DA}})=\mathrm{Im}(\tilde{\mathrm{V}}_{\mathrm{DDI,DA}})$.
Third, the individual excited-state population depends on both
$\mathrm{Im}(\mathrm{V})$ and $\mathrm{Re}(\mathrm{V})$, where
$\mathrm{Im}(\mathrm{V})$ modifies the decay behavior and
$\mathrm{Re}(\mathrm{V})$ generates the oscillatory pattern. Since
$\mathrm{Re}(\mathrm{V}_{\mathrm{DDI,DA}})\neq\mathrm{Re}(\tilde{\mathrm{V}}_{\mathrm{DDI,DA}})$,
the individual excited-state population is sensitive to the RWA. Moreover, the
two-times oscillation frequency can be explained by the fact that
$\mathrm{Re}(\mathrm{V}_{\mathrm{DDI,DA}})\approx
2\mathrm{Re}(\tilde{\mathrm{V}}_{\mathrm{DDI,DA}})$ in the short distance
limit in Eq. (28).

In FIG. 3(d)-(f), when $h=1$ nm and $d=1$ nm, the population dynamics obtained
via FQD and MAQD no longer match each other. In other words, we can identify
this circumstance as a strong coupling condition. In this case, not only the
individual population dynamics, i.e., $P^{\mathrm{E_{D},\\{0\\}}}(t)$ and
$P^{\mathrm{E_{A},\\{0\\}}}(t)$, but also the total excited population, i.e.,
$P^{\mathrm{E_{D},\\{0\\}}}(t)+P^{\mathrm{E_{A},\\{0\\}}}(t)$, are sensitive
to the RWA.

In summary, the numerical and analytical analysis of the population dynamics
in the specific system reveals the following findings. (i) Under weak coupling
conditions, the RWA can affect the excited-state population of individual
molecules while leaving the total excited-state population unchanged (in our
chosen system). However, we would like to emphasize that the total excited-
state population is unaffected by the RWA in the weak coupling regime is not a
general principle. In fact, in cases where the donor and acceptor are at
different heights, the donor and acceptor are non-identical, or when more than
two molecules are involved, the total excited-state population becomes
sensitive to the use of the RWA [56]. (ii) Under strong coupling conditions,
both the individual and total excited-state dynamics are highly sensitive to
the RWA. These results further emphasize the significance of counter-rotating
interactions in quantum dynamics in both strong and weak coupling regimes.

##  V Conclusion

In this study, we investigate the influence of counter-rotating interactions
on the quantum dynamics of multiple molecules in complex dielectric
environments within the framework of MQED. Our general theory of quantum
dynamics shows that the neglect of the counter-rotating interactions leads to
missing several crucial physical processes, e.g., virtual photon emission and
reabsorption. We summarize our main findings as follows. First, in the weak
coupling regime, the lack of these processes leads to the absence of energy
shifts in the ground-state molecules and incorrect dipole-dipole interactions
between molecule pairs. Second, our study clearly demonstrates that counter-
rotating interactions play an essential component in dipole-dipole
interactions. Our analysis points out that within the RWA, the dipole-dipole
interactions converge to only half of the conventional Coulomb interaction in
the short-distance (non-retarded) limit. Third, our numerical simulations
reveal that in the weak coupling regime, the absence of the counter-rotating
interactions can significantly influence the dynamics of individual molecules
while leaving the total excited-state population unchanged. Conversely, in the
strong coupling regime, both individual and total excited-state dynamics
exhibit sensitivity to the counter-rotating interactions. To sum up, through
the analysis of the dynamical equations and a specific case study of
population dynamics, we show that the counter-rotating interactions play
crucial roles in both strong and weak coupling regimes. Hence, it is
imperative to always exercise caution when making the rotating-wave
approximation. We believe that this work will provide important insights into
the study of light-matter interactions.

While we have demonstrated the significance of counter-rotating interactions,
there are still unresolved issues that warrant further investigations. First,
in this article, we have modeled molecules as two-level systems, neglecting
the influence of other molecular excited states. However, these excited states
play a crucial role in certain properties, such as energy shift calculations.
Second, our analysis includes up to two molecular excitations and one
polariton in the wavefunction ansatz. Nevertheless, in regimes of ultrastrong
and deep-strong coupling, higher excitation states, such as three-molecule
excitations and two-polariton states, may also impact quantum dynamics
significantly. These issues are important for exploring the quantum dynamics
of a collection of molecules coupled with quantum light. Finally, we hope that
this work could inspire further investigations into emerging quantum
electrodynamic phenomena in chemistry and molecular physics.

###### Acknowledgements.

Hsu thanks Academia Sinica (AS-CDA-111-M02) and the Ministry of Science and
Technology of Taiwan (110-2113-M-001-053 and 111-2113-M-001-027-MY4) for the
financial support.

##  Appendix A Derivation of Eq. (12)

To derive Eq. (12), we substitute the Hamiltonian $\hat{H}$ in Eq. (1) and the
state vector in Eq. (10) into the time-dependent Schrödinger equation
$i\hbar\partial\ket{\Psi(t)}/\partial t=\hat{H}\ket{\Psi(t)}$, and then we
obtain the following coupled differential equations,

$\displaystyle
i\hbar\partialderivative{t}C^{\mathrm{G},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t)e^{-iW^{\mathrm{G},\left\\{1\right\\}}(\omega)t}=-\sum_{\alpha}\left[\bm{\mu}^{\mathrm{ge}}_{\alpha}\cdot\overline{\overline{\mathcal{G}}}\vphantom{a}^{*}(\mathbf{r}_{\alpha},\mathbf{r},\omega)\right]_{k}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)e^{-iW^{\mathrm{E}_{\alpha},\left\\{0\right\\}}t},$
(31) $\displaystyle
i\hbar\partialderivative{t}C^{\mathrm{E_{\alpha\beta}},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t)e^{-iW^{\mathrm{E}_{\alpha\beta},\left\\{1\right\\}}(\omega)t}=$
$\displaystyle\quad-\left[\bm{\mu}^{\mathrm{eg}}_{\beta}\cdot\overline{\overline{\mathcal{G}}}\vphantom{a}^{*}(\mathbf{r}_{\beta},\mathbf{r},\omega)\right]_{k}e^{-iW^{\mathrm{E}_{\alpha},\left\\{0\right\\}}t}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)-\left[\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\overline{\overline{\mathcal{G}}}\vphantom{a}^{*}(\mathbf{r}_{\alpha},\mathbf{r},\omega)\right]_{k}e^{-iW^{\mathrm{E}_{\beta},\left\\{0\right\\}}t}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t),$
(32) $\displaystyle
i\hbar\derivative{t}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)e^{-iW^{\mathrm{E}_{\alpha},\left\\{0\right\\}}t}=$
$\displaystyle-\sum_{k=1}^{3}\int\differential{\mathbf{r}}\int_{0}^{\infty}\differential{\omega}\left[\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\overline{\overline{\mathcal{G}}}(\mathbf{r}_{\alpha},\mathbf{r},\omega)\right]_{k}e^{-iW^{\mathrm{G},\left\\{1\right\\}}(\omega)t}C^{\mathrm{G},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t)$
$\displaystyle-\sum_{\beta\neq\alpha}\sum_{k=1}^{3}\int\differential{\mathbf{r}}\int_{0}^{\infty}\differential{\omega}\left[\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\overline{\overline{\mathcal{G}}}(\mathbf{r}_{\beta},\mathbf{r^{\prime}},\omega)\right]_{k}e^{-iW^{\mathrm{E}_{\alpha\beta},\left\\{1\right\\}}(\omega)t}C^{\mathrm{E_{\alpha\beta}},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t),$
(33)

where $\left[\mathbf{v}\right]_{k}$ denotes the $k$-th component of the vector
$\mathbf{v}$. Note that we have ignored the two-polariton term and the three-
molecule-excitation term in obtaining Eqs. (31)-(33). Consider that there is
no polariton at $t=0$, i.e.,
$C^{\mathrm{G},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t=0)=C^{\mathrm{E_{\alpha\beta}},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t=0)=0$,
we can formally integrate Eqs. (31) and (32) and obtain

$\displaystyle
C^{\mathrm{G},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t)=\frac{i}{\hbar}\sum_{\alpha}\int_{0}^{t}\differential{t^{\prime}}\left[\bm{\mu}^{\mathrm{ge}}_{\alpha}\cdot\overline{\overline{\mathcal{G}}}\vphantom{a}^{*}(\mathbf{r}_{\alpha},\mathbf{r},\omega)\right]_{k}e^{-i\left(W^{\mathrm{E}_{\alpha},\left\\{0\right\\}}-W^{\mathrm{G},\left\\{1\right\\}}(\omega)\right)t^{\prime}}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime}),$
(34) $\displaystyle
C^{\mathrm{E_{\alpha\beta}},\left\\{1_{k}\right\\}}(\mathbf{r},\omega,t)=$
$\displaystyle\frac{i}{\hbar}\int_{0}^{t}\differential{t^{\prime}}\left[\bm{\mu}^{\mathrm{eg}}_{\beta}\cdot\overline{\overline{\mathcal{G}}}\vphantom{a}^{*}(\mathbf{r}_{\beta},\mathbf{r},\omega)\right]_{k}e^{-i\left(W^{\mathrm{E}_{\alpha},\left\\{0\right\\}}-W^{\mathrm{E}_{\alpha\beta},\left\\{1\right\\}}(\omega)\right)t^{\prime}}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle+\frac{i}{\hbar}\int_{0}^{t}\differential{t^{\prime}}\left[\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\overline{\overline{\mathcal{G}}}\vphantom{a}^{*}(\mathbf{r}_{\alpha},\mathbf{r},\omega)\right]_{k}e^{-i\left(W^{\mathrm{E}_{\beta},\left\\{0\right\\}}-W^{\mathrm{E}_{\alpha\beta},\left\\{1\right\\}}(\omega)\right)t^{\prime}}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t^{\prime}).$
(35)

In order to derive Eq. (12), we substitute Eqs. (34) and (35) into Eq. (33),
make use of the identity [65]

$\displaystyle\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r},\mathbf{r^{\prime}},\omega)=\int\differential{\mathbf{s}}\frac{\omega^{2}\mathrm{Im}\left[\varepsilon_{r}(\mathbf{s},\omega)\right]}{c^{2}}\overline{\overline{\mathbf{G}}}(\mathbf{r},\mathbf{s},\omega)\overline{\overline{\mathbf{G}}}\vphantom{G}^{\dagger}(\mathbf{r^{\prime}},\mathbf{s},\omega),$
(36)

apply the definitions of $W^{\mathrm{E}_{\alpha},\left\\{0\right\\}}$,
$W^{\mathrm{G},\left\\{1\right\\}}(\omega)$, and
$W^{\mathrm{E}_{\alpha\beta},\left\\{1\right\\}}(\omega)$ in Eqs (11a)-(11c),
and finally obtain

$\displaystyle\derivative{t}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=$
$\displaystyle\quad-\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}\right]e^{-i\left(\omega-\omega_{\alpha}\right)t}e^{-i\left(\omega_{\alpha}-\omega\right)t^{\prime}}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle\quad-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}\right]e^{-i\left(\omega-\omega_{\alpha}\right)t}e^{-i\left(\omega_{\beta}-\omega\right)t^{\prime}}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle\quad-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\beta},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{eg}}_{\beta}\right]e^{-i\left(\omega+\omega_{\beta}\right)t}e^{-i\left(-\omega_{\beta}-\omega\right)t^{\prime}}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle\quad-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\beta},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{eg}}_{\alpha}\right]e^{-i\left(\omega+\omega_{\beta}\right)t}e^{-i\left(-\omega_{\alpha}-\omega\right)t^{\prime}}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t^{\prime}).$
(37)

##  Appendix B Derivation of Eq. (17)

In the weak coupling regime, we apply the Markov approximation to Eq. (12),
i.e., we change
$C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime})\rightarrow
C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)$ and
$\int_{0}^{t}\differential{t^{\prime}}\rightarrow\int_{-\infty}^{t}\differential{t^{\prime}}$,
and make the substitution $\tau=t-t^{\prime}$
($\differential{\tau}=-\differential{t^{\prime}}$); then Eq. (12) becomes

$\displaystyle\derivative{t}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=$
$\displaystyle\quad-\left\\{\int_{0}^{\infty}\differential{\omega}\left[\int_{0}^{\infty}\differential{\tau}e^{-i\left(\omega-\omega_{\alpha}\right)\tau}\right]\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}\right]\right\\}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)$
$\displaystyle\quad-\sum_{\beta\neq\alpha}\left\\{\int_{0}^{\infty}\differential{\omega}\left[\int_{0}^{\infty}\differential{\tau}e^{-i\left(\omega-\omega_{\beta}\right)\tau}\right]\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}\right]\right\\}e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t)$
$\displaystyle\quad-\sum_{\beta\neq\alpha}\left\\{\int_{0}^{\infty}\differential{\omega}\left[\int_{0}^{\infty}\differential{\tau}e^{-i\left(\omega+\omega_{\beta}\right)\tau}\right]\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\beta},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{eg}}_{\beta}\right]\right\\}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)$
$\displaystyle\quad-\sum_{\beta\neq\alpha}\left\\{\int_{0}^{\infty}\differential{\omega}\left[\int_{0}^{\infty}\differential{\tau}e^{-i\left(\omega+\omega_{\alpha}\right)\tau}\right]\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\beta},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{eg}}_{\alpha}\right]\right\\}e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t).$
(38)

According to the Sokhotski–Plemelj theorem, we have

$\displaystyle\int_{0}^{\infty}\differential{\tau}e^{-i\left(\omega\pm\omega_{\alpha}\right)\tau}=\pi\delta\left(\omega\pm\omega_{\alpha}\right)-i\mathcal{P}\left(\frac{1}{\omega\pm\omega_{\alpha}}\right).$
(39)

Substituting Eq. (39) into Eq. (38), we obtain

$\displaystyle\derivative{t}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=$
$\displaystyle\quad-\left\\{\left[\frac{\omega_{\alpha}^{2}}{\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega_{\alpha})\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}\right]-i\mathcal{P}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}}{\omega-\omega_{\alpha}}\right]\right\\}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)$
$\displaystyle\quad-\sum_{\beta\neq\alpha}\left\\{\left[\frac{\omega_{\beta}^{2}}{\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega_{\beta})\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}\right]-i\mathcal{P}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}}{\omega-\omega_{\beta}}\right]\right\\}e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t)$
$\displaystyle\quad-\sum_{\beta\neq\alpha}\left\\{-i\mathcal{P}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\beta},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{eg}}_{\beta}}{\omega+\omega_{\beta}}\right]\right\\}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)$
$\displaystyle\quad-\sum_{\beta\neq\alpha}\left\\{-i\mathcal{P}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\beta},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{eg}}_{\alpha}}{\omega+\omega_{\alpha}}\right]\right\\}e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t).$
(40)

Using the following identities [65]

$\displaystyle\overline{\overline{\mathbf{G}}}(\mathbf{r},\mathbf{r}^{\prime},-\omega^{*})=\overline{\overline{\mathbf{G}}}\vphantom{G}^{*}(\mathbf{r},\mathbf{r}^{\prime},\omega),$
(41)
$\displaystyle\bm{\mu}_{1}\cdot\overline{\overline{\mathbf{G}}}(\mathbf{r}_{1},\mathbf{r}_{2},\omega)\cdot\bm{\mu}_{2}=\bm{\mu}_{2}\cdot\overline{\overline{\mathbf{G}}}(\mathbf{r}_{2},\mathbf{r}_{1},\omega)\cdot\bm{\mu}_{1},$
(42)
$\displaystyle\pi\omega^{2}\mathrm{Re}\overline{\overline{\mathbf{G}}}(\mathbf{r},\mathbf{r}^{\prime},\omega)$
$\displaystyle=\mathcal{P}\int_{-\infty}^{\infty}\differential{\omega}\frac{\omega^{2}\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r},\mathbf{r}^{\prime},\omega)}{\omega-\omega_{\alpha}}$
$\displaystyle=\mathcal{P}\int_{0}^{\infty}\differential{\omega}\frac{\omega^{2}\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r},\mathbf{r}^{\prime},\omega)}{\omega-\omega_{\alpha}}+\int_{0}^{\infty}\differential{\omega}\frac{\omega^{2}\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r},\mathbf{r}^{\prime},\omega)}{\omega+\omega_{\alpha}},$
(43)

and recalling the definitions of $\Delta_{\mathrm{e(g)}_{\alpha}}$ in Eq (18),
$\Gamma_{\alpha}$ in Eq (21), and $\mathrm{V}_{\mathrm{DDI,\alpha\beta}}$ in
Eq (22), we can transform Eq (40) into

$\displaystyle\derivative{t}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=-\frac{i}{\hbar}\left\\{\left[\Delta_{\mathrm{e}_{\alpha}}+\sum_{\beta\neq\alpha}\Delta_{\mathrm{g}_{\beta}}\right]-i\hbar\frac{\Gamma_{\alpha}}{2}\right\\}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)-\frac{i}{\hbar}\sum_{\beta\neq\alpha}\mathrm{V}_{\mathrm{DDI,\alpha\beta}}\,e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t).$
(44)

##  Appendix C Calculation of $\mathrm{V}_{\mathrm{ORC,\alpha\beta}}$ and
$\mathrm{V}_{\mathrm{QC,\alpha\beta}}$ on the Imaginary Axis

To obtain $\mathrm{V}_{\mathrm{ORC,\alpha\beta}}$ and
$\mathrm{V}_{\mathrm{QC,\alpha\beta}}$, we need to evaluate the integral

$\displaystyle
I=\int_{0}^{\infty}\differential{\omega}\frac{\omega^{2}}{\pi\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}}{\omega+\omega^{\prime}},$
(45)

where $\omega^{\prime}=\omega_{\alpha}$ or $\omega_{\beta}$. $I$ can be
alternatively expressed as:

$\displaystyle
I=\frac{1}{\pi\varepsilon_{0}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\left[\int_{0}^{\infty}\differential{\omega}\frac{\omega^{2}}{c^{2}}\frac{\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)}{\omega+\omega^{\prime}}\right]\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}.$
(46)

For simplicity, we define
$f(\omega)=\frac{\omega^{2}}{c^{2}}\frac{\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)}{\omega+\omega^{\prime}}$.
The integral in the square bracket can be evaluated using the contour integral
technique, i.e.,

$\displaystyle\int_{0}^{\infty}\differential{\omega}f(\omega)=\oint_{C}\differential{\omega}f(\omega)-\int_{C_{2}}\differential{\omega}f(\omega)-\int_{C_{3}}\differential{\omega}f(\omega),$
(47)

where the contour is shown in FIG. 4. The first term on the right-hand side of
Eq. (47) is zero since there is no singularity inside the contour $C$. The
second term on the right-hand side of Eq. (47) is also zero due to the
asymptotic behavior of the two-point dyadic Green’s function [65]

$\displaystyle\lim_{\absolutevalue{\omega}\rightarrow\infty}\left.\frac{\omega^{2}}{c^{2}}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\right|_{\mathbf{r}_{\alpha}\neq\mathbf{r}_{\beta}}=\overline{\overline{\mathbf{0}}}.$
(48)

Figure 4: The contour adopted in the integral. The total contour $C$ is equal
to $C_{1}+C_{2}+C_{3}$. Note that there is no singularity inside the contour.

Therefore, we have

$\displaystyle\int_{0}^{\infty}\differential{\omega}f(\omega)$
$\displaystyle=-\int_{i\infty}^{0}\differential{\omega}f(\omega)$
$\displaystyle=-\int_{i\infty}^{0}\differential{\omega}\frac{\omega^{2}}{c^{2}}\frac{\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)}{\omega+\omega^{\prime}}$
$\displaystyle=-\int_{i\infty}^{0}\differential{\omega}\frac{\omega^{2}\left(\omega-\omega^{\prime}\right)}{c^{2}}\frac{\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)}{\omega^{2}-\omega^{\prime
2}}$
$\displaystyle=-\int_{i\infty}^{0}\differential{\omega}\frac{\omega^{3}}{c^{2}}\frac{\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)}{\omega^{2}-\omega^{\prime
2}}$
$\displaystyle\quad+\int_{i\infty}^{0}\differential{\omega}\frac{\omega^{\prime}\omega^{2}}{c^{2}}\frac{\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)}{\omega^{2}-\omega^{\prime
2}}.$ (49)

Making the substitution $\kappa=-i\omega$,

$\displaystyle\int_{0}^{\infty}\differential{\omega}f(\omega)$
$\displaystyle=-\int_{0}^{\infty}\differential{\kappa}\frac{\kappa^{3}}{c^{2}}\frac{\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},i\kappa)}{\kappa^{2}+\omega^{\prime
2}}$
$\displaystyle\quad-i\int_{0}^{\infty}\differential{\kappa}\frac{\omega^{\prime}\kappa^{2}}{c^{2}}\frac{\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},i\kappa)}{\kappa^{2}+\omega^{\prime
2}}.$ (50)

Using Eq. (41), we can obtain that the dyadic Green’s function is purely real
on the imaginary axis, i.e.,
$\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},i\kappa)=\mathrm{Re}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},i\kappa)$,
and arrive at

$\displaystyle\mathrm{Im}\left[\int_{0}^{\infty}\differential{\omega}f(\omega)\right]=-\int_{0}^{\infty}\differential{\kappa}\frac{\omega^{\prime}\kappa^{2}}{c^{2}}\frac{\mathrm{Re}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},i\kappa)}{\kappa^{2}+\omega^{\prime
2}}.$ (51)

Now we can express the integral $I$ as [63, 64]:

$\displaystyle
I=-\int_{0}^{\infty}\differential{\kappa}\frac{\omega^{\prime}\kappa^{2}}{\pi\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Re}\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},i\kappa)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}}{\kappa^{2}+\omega^{\prime
2}}.$ (52)

Eq. (52) is a powerful tool in cases where the dyadic Green’s function on the
imaginary axis is available since
$\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},i\kappa)$
decays rapidly [63]. However, in complex dielectric environments, dyadic
Green’s function on the imaginary axis is difficult to obtain, and the
evaluation of $I$ through Eq. (45) is more convenient.

In free space, we have the explicit expression of the two-point dyadic Green’s
function, which reads:

$\displaystyle\left.\overline{\overline{\mathbf{G}}}_{{\mathrm{0}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\right|_{\mathbf{r}_{\alpha}\neq\mathbf{r}_{\beta}}=$
$\displaystyle\quad\frac{e^{ik_{0}R}}{4\pi
R}\Biggl{\\{}\vphantom{\frac{e^{R}}{R}}\left(\overline{\overline{\mathbf{I}}}_{3}-\mathbf{n}_{R}\otimes\mathbf{n}_{R}\right)$
$\displaystyle\quad\quad+\left(3\mathbf{n}_{R}\otimes\mathbf{n}_{R}-\overline{\overline{\mathbf{I}}}_{3}\right)\left[\frac{1}{(k_{0}R)^{2}}-\frac{i}{k_{0}R}\right]\Biggl{\\}},$
(53)

where $k_{0}=\omega/c$ and $\mathbf{r}_{\alpha}-\mathbf{r}_{\beta}\equiv
R\mathbf{n}_{R}$. Therefore, we can evaluate the integral using Eq. (52).
Inserting Eq. (53) into Eq. (52) and making the substitution $x=\kappa R/c$,
we obtain

$\displaystyle I^{0}$
$\displaystyle\equiv-\int_{0}^{\infty}\differential{\kappa}\frac{\omega^{\prime}\kappa^{2}}{\pi\varepsilon_{0}c^{2}}\frac{\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Re}\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},i\kappa)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}}{\kappa^{2}+\omega^{\prime
2}}$
$\displaystyle=-\frac{\omega^{\prime}}{4\pi^{2}\varepsilon_{0}cR^{2}}\Bigl{\\{}\left[\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}-\left(\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathbf{n}_{R}\right)\left(\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathbf{n}_{R}\right)\right]\mathcal{I}_{1}$
$\displaystyle\quad+\left[\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}-3\left(\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathbf{n}_{R}\right)\left(\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathbf{n}_{R}\right)\right]\left(\mathcal{I}_{2}+\mathcal{I}_{3}\right)\Bigl{\\}}.$
(54)

$\mathcal{I}_{1}$, $\mathcal{I}_{2}$, and $\mathcal{I}_{3}$ are the auxiliary
integrals defined as

$\displaystyle\mathcal{I}_{1}=\int_{0}^{\infty}\differential{x}\frac{x^{2}e^{-x}}{x^{2}+x^{\prime
2}},$ (55a)
$\displaystyle\mathcal{I}_{2}=\int_{0}^{\infty}\differential{x}\frac{xe^{-x}}{x^{2}+x^{\prime
2}},$ (55b)
$\displaystyle\mathcal{I}_{3}=\int_{0}^{\infty}\differential{x}\frac{e^{-x}}{x^{2}+x^{\prime
2}},$ (55c)

where $x^{\prime}=\omega^{\prime}R/c$. These auxiliary integrals can be
explicitly expressed in terms of the trigonometric functions and trigonometric
integrals as [80]

$\displaystyle\mathcal{I}_{1}$
$\displaystyle=-x^{\prime}\left[\mathrm{ci}(x^{\prime})\sin(x^{\prime})-\mathrm{si}(x^{\prime})\cos(x^{\prime})\right]+1,$
(56a) $\displaystyle\mathcal{I}_{2}$
$\displaystyle=-\mathrm{ci}(x^{\prime})\cos(x^{\prime})-\mathrm{si}(x^{\prime})\sin(x^{\prime}),$
(56b) $\displaystyle\mathcal{I}_{3}$
$\displaystyle=\frac{1}{x^{\prime}}\left[\mathrm{ci}(x^{\prime})\sin(x^{\prime})-\mathrm{si}(x^{\prime})\cos(x^{\prime})\right],$
(56c)

where the trigonometric integrals follow the definitions:

$\displaystyle\mathrm{ci}(x^{\prime})=-\int_{x^{\prime}}^{\infty}\differential{x}\frac{\cos(x)}{x},\,\,\mathrm{si}(x^{\prime})=-\int_{x^{\prime}}^{\infty}\differential{x}\frac{\sin(x)}{x}.$

The explicit form of $I^{0}$ is useful for the calculation of quantum
dynamics. In addition, the asymptotic behavior of $I^{0}$,
$V^{0}_{\mathrm{ORC},\alpha\beta}$, and $V^{0}_{\mathrm{QC},\alpha\beta}$ can
be analyzed through the expansion of the trigonometric integrals.

##  Appendix D Numerical Implementation of FQD and MAQD

To numerically implement FQD without the RWA, we start from Eq. (12), separate
$\overline{\overline{\mathbf{G}}}(\mathbf{r},\mathbf{r}^{\prime},\omega)$ as
$\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)+\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r},\mathbf{r}^{\prime},\omega)$,
apply the Markov approximation only to the part involve
$\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)$,
and discard the free-space Lamb shift $\Delta^{0}_{\mathrm{e(g)_{\alpha}}}$.
As a result, one can numerically calculate FQD without the RWA instead of
using Eq. (12) as follows,

$\displaystyle\derivative{t}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=$
$\displaystyle-\frac{\Gamma^{0}_{\alpha}}{2}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)-\frac{i}{\hbar}\sum_{\beta\neq\alpha}\mathrm{V}^{0}_{\mathrm{DDI,\alpha\beta}}\,e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t)$
$\displaystyle-\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}\right]e^{-i\left(\omega-\omega_{\alpha}\right)t}e^{-i\left(\omega_{\alpha}-\omega\right)t^{\prime}}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}\right]e^{-i\left(\omega-\omega_{\alpha}\right)t}e^{-i\left(\omega_{\beta}-\omega\right)t^{\prime}}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r}_{\beta},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{eg}}_{\beta}\right]e^{-i\left(\omega+\omega_{\beta}\right)t}e^{-i\left(-\omega_{\beta}-\omega\right)t^{\prime}}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{ge}}_{\beta}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r}_{\beta},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{eg}}_{\alpha}\right]e^{-i\left(\omega+\omega_{\beta}\right)t}e^{-i\left(-\omega_{\alpha}-\omega\right)t^{\prime}}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t^{\prime}),$
(57)

where
$\Gamma^{0}_{\alpha}=\frac{2\omega_{\alpha}^{2}}{\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega_{\alpha})\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}=\frac{\absolutevalue{\bm{\mu}^{\mathrm{ge}}_{\alpha}}^{2}\omega_{\alpha}^{3}}{3\pi\hbar\varepsilon_{0}c^{3}}$
is the spontaneous emission rate in free space, and
$\mathrm{V}^{0}_{\mathrm{DDI,\alpha\beta}}$ is obtained by substituting
$\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)$
with
$\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)$
in the expression of $\mathrm{V}_{\mathrm{DDI,\alpha\beta}}$.

To numerically implement FQD with the RWA, we start from Eq. (16) and adopt
the same procedure as in the above. As a result, we can numerically calculate
FQD with the RWA instead of Eq. (16) as follows,

$\displaystyle\derivative{t}\tilde{C}^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=$
$\displaystyle-\frac{\Gamma^{0}_{\alpha}}{2}\tilde{C}^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)-\frac{i}{\hbar}\sum_{\beta\neq\alpha}\tilde{\mathrm{V}}^{0}_{\mathrm{DDI,\alpha\beta}}\,e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}\tilde{C}^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t)$
$\displaystyle-\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r}_{\alpha},\mathbf{r}_{\alpha},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\alpha}\right]e^{-i\left(\omega-\omega_{\alpha}\right)t}e^{-i\left(\omega_{\alpha}-\omega\right)t^{\prime}}\tilde{C}^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t^{\prime})$
$\displaystyle-\sum_{\beta\neq\alpha}\int_{0}^{t}\differential{t^{\prime}}\int_{0}^{\infty}\differential{\omega}\left[\frac{\omega^{2}}{\pi\hbar\varepsilon_{0}c^{2}}\bm{\mu}^{\mathrm{eg}}_{\alpha}\cdot\mathrm{Im}\overline{\overline{\mathbf{G}}}_{\mathrm{Sc}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)\cdot\bm{\mu}^{\mathrm{ge}}_{\beta}\right]e^{-i\left(\omega-\omega_{\alpha}\right)t}e^{-i\left(\omega_{\beta}-\omega\right)t^{\prime}}\tilde{C}^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t^{\prime}),$
(58)

where $\tilde{\mathrm{V}}^{0}_{\mathrm{DDI,\alpha\beta}}$ is obtained by
substituting
$\overline{\overline{\mathbf{G}}}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)$
with
$\overline{\overline{\mathbf{G}}}_{0}(\mathbf{r}_{\alpha},\mathbf{r}_{\beta},\omega)$
in the expression of $\tilde{\mathrm{V}}_{\mathrm{DDI,\alpha\beta}}$.

To numerically implement MAQD without the RWA, we discard the free-space Lamb
shift $\Delta^{0}_{\mathrm{e(g)_{\alpha}}}$ in Eq. (17) and finally obtain

$\displaystyle\derivative{t}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=-\frac{i}{\hbar}\left\\{\left[\Delta^{\mathrm{Sc}}_{\mathrm{e}_{\alpha}}+\sum_{\beta\neq\alpha}\Delta^{\mathrm{Sc}}_{\mathrm{g}_{\beta}}\right]-i\hbar\frac{\Gamma_{\alpha}}{2}\right\\}C^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)-\frac{i}{\hbar}\sum_{\beta\neq\alpha}\mathrm{V}_{\mathrm{DDI,\alpha\beta}}\,e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}C^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t).$
(59)

To numerically implement MAQD with the RWA, we discard the free-space Lamb
shift $\Delta^{0}_{\mathrm{e_{\alpha}}}$ in Eq. (25) and finally arrive at

$\displaystyle\derivative{t}\tilde{C}^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)=-\frac{i}{\hbar}\left[\Delta^{\mathrm{Sc}}_{\mathrm{e}_{\alpha}}-i\hbar\frac{\Gamma_{\alpha}}{2}\right]\tilde{C}^{\mathrm{E_{\alpha}},\left\\{0\right\\}}(t)-\frac{i}{\hbar}\sum_{\beta\neq\alpha}\tilde{\mathrm{V}}_{\mathrm{DDI,\alpha\beta}}\,e^{-i\left(\omega_{\beta}-\omega_{\alpha}\right)t}\tilde{C}^{\mathrm{E_{\beta}},\left\\{0\right\\}}(t).$
(60)

## References

  * Bonifacio  _et al._ [1971a] R. Bonifacio, P. Schwendimann, and F. Haake, Phys. Rev. A 4, 302 (1971a).
  * Bonifacio  _et al._ [1971b] R. Bonifacio, P. Schwendimann, and F. Haake, Phys. Rev. A 4, 854 (1971b).
  * Morawitz [1979] H. Morawitz, IBM J. Res. Dev. 23, 517 (1979).
  * Seke [1986a] J. Seke, Phys. Rev. A 33, 739 (1986a).
  * Seke [1986b] J. Seke, Phys. Rev. A 33, 4409 (1986b).
  * Kulina  _et al._ [1986] P. Kulina, C. Leonard, and R. H. Rinkleff, Phys. Rev. A 34, 227 (1986).
  * Bashkirov  _et al._ [1990] E. Bashkirov, F. Le Kien, and A. Shumovsky, Phys. A: Stat. Mech. Appl. 167, 935 (1990).
  * Alber [1992] G. Alber, Phys. Rev. A 46, R5338 (1992).
  * Kozierowski and Chumakov [1995] M. Kozierowski and S. M. Chumakov, Phys. Rev. A 52, 4194 (1995).
  * Temnov and Woggon [2005] V. V. Temnov and U. Woggon, Phys. Rev. Lett. 95, 243602 (2005).
  * Martín-Cano  _et al._ [2010] D. Martín-Cano, L. Martín-Moreno, F. J. García-Vidal, and E. Moreno, Nano Lett. 10, 3129 (2010).
  * Goban  _et al._ [2015] A. Goban, C.-L. Hung, J. D. Hood, S.-P. Yu, J. A. Muniz, O. Painter, and H. J. Kimble, Phys. Rev. Lett. 115, 063601 (2015).
  * Zhang  _et al._ [2016] Y. Zhang, Y. Luo, Y. Zhang, Y.-J. Yu, Y.-M. Kuang, L. Zhang, Q.-S. Meng, Y. Luo, J.-L. Yang, Z.-C. Dong, and J. G. Hou, Nature 531, 623 (2016).
  * Kim  _et al._ [2018] J.-H. Kim, S. Aghaeimeibodi, C. J. K. Richardson, R. P. Leavitt, and E. Waks, Nano Lett. 18, 4734 (2018).
  * Luo  _et al._ [2019] Y. Luo, G. Chen, Y. Zhang, L. Zhang, Y. Yu, F. Kong, X. Tian, Y. Zhang, C. Shan, Y. Luo, J. Yang, V. Sandoghdar, Z. Dong, and J. G. Hou, Phys. Rev. Lett. 122, 233901 (2019).
  * Wang and Hsu [2023] S. Wang and L.-Y. Hsu, J. Phys. Chem. C 127, 12904 (2023).
  * Gersten and Nitzan [1984] J. I. Gersten and A. Nitzan, Chem. Phys. Lett. 104, 31 (1984).
  * Gersten and Nitzan [1985] J. I. Gersten and A. Nitzan, Surf. Sci. 158, 165 (1985).
  * Hua  _et al._ [1985] X. M. Hua, J. I. Gersten, and A. Nitzan, J. Chem. Phys. 83, 3650 (1985).
  * Kobayashi  _et al._ [1995] T. Kobayashi, Q. Zheng, and T. Sekiguchi, Phys. Rev. A 52, 2835 (1995).
  * Andrew and Barnes [2000] P. Andrew and W. L. Barnes, Science 290, 785 (2000).
  * Zhang  _et al._ [2007] J. Zhang, Y. Fu, and J. R. Lakowicz, J. Phys. Chem. C 111, 50 (2007).
  * Govorov  _et al._ [2007] A. O. Govorov, J. Lee, and N. A. Kotov, Phys. Rev. B 76, 125308 (2007).
  * Reil  _et al._ [2008] F. Reil, U. Hohenester, J. R. Krenn, and A. Leitner, Nano Lett. 8, 4128 (2008).
  * Lessard-Viger  _et al._ [2009] M. Lessard-Viger, M. Rioux, L. Rainville, and D. Boudreau, Nano Lett. 9, 3066 (2009).
  * Lunz  _et al._ [2011] M. Lunz, V. A. Gerard, Y. K. Gun’ko, V. Lesnyak, N. Gaponik, A. S. Susha, A. L. Rogach, and A. L. Bradley, Nano Lett. 11, 3341 (2011).
  * Zhong  _et al._ [2017] X. Zhong, T. Chervy, L. Zhang, A. Thomas, J. George, C. Genet, J. A. Hutchison, and T. W. Ebbesen, Angew. Chem. Int. Ed. 56, 9034 (2017).
  * Ren  _et al._ [2017] J. Ren, T. Chen, B. Wang, and X. Zhang, J. Chem. Phys. 146, 144101 (2017).
  * Wu  _et al._ [2018] J.-S. Wu, Y.-C. Lin, Y.-L. Sheu, and L.-Y. Hsu, J. Phys. Chem. Lett. 9, 7032 (2018).
  * Hutchison  _et al._ [2012] J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, Angew. Chem. Int. Ed. 51, 1592 (2012).
  * Thomas  _et al._ [2016] A. Thomas, J. George, A. Shalabney, M. Dryzhakov, S. J. Varma, J. Moran, T. Chervy, X. Zhong, E. Devaux, C. Genet, J. A. Hutchison, and T. W. Ebbesen, Angew. Chem. Int. Ed. 55, 11462 (2016).
  * Thomas  _et al._ [2019] 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, Science 363, 615 (2019).
  * Lather  _et al._ [2019] J. Lather, P. Bhatt, A. Thomas, T. W. Ebbesen, and J. George, Angew. Chem. Int. Ed. 58, 10635 (2019).
  * Vergauwe  _et al._ [2019] R. M. A. Vergauwe, A. Thomas, K. Nagarajan, A. Shalabney, J. George, T. Chervy, M. Seidel, E. Devaux, V. Torbeev, and T. W. Ebbesen, Angew. Chem. Int. Ed. 58, 15324 (2019).
  * Thomas  _et al._ [2020] A. Thomas, A. Jayachandran, L. Lethuillier-Karl, R. M. Vergauwe, K. Nagarajan, E. Devaux, C. Genet, J. Moran, and T. W. Ebbesen, Nanophotonics 9, 249 (2020).
  * Sau  _et al._ [2021] A. Sau, K. Nagarajan, B. Patrahau, L. Lethuillier‐Karl, R. M. A. Vergauwe, A. Thomas, J. Moran, C. Genet, and T. W. Ebbesen, Angew. Chem. Int. Ed. 60, 5712 (2021).
  * Galego  _et al._ [2015] J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. X 5, 041022 (2015).
  * Herrera and Spano [2016] F. Herrera and F. C. Spano, Phys. Rev. Lett. 116, 238301 (2016).
  * Wu  _et al._ [2016] N. Wu, J. Feist, and F. J. Garcia-Vidal, Phys. Rev. B 94, 195409 (2016).
  * Ribeiro  _et al._ [2018] R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Campos-Gonzalez-Angulo, and J. Yuen-Zhou, Chem. Sci. 9, 6325 (2018).
  * Semenov and Nitzan [2019] A. Semenov and A. Nitzan, J. Chem. Phys. 150, 174122 (2019).
  * Spano [2020] F. C. Spano, J. Chem. Phys. 152, 204113 (2020).
  * Du and Yuen-Zhou [2022] M. Du and J. Yuen-Zhou, Phys. Rev. Lett. 128, 096001 (2022).
  * Scully and Zubairy [1997] M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * Vogel and Welsch [2006] W. Vogel and D. Welsch,  _Quantum Optics_, 3rd ed. (John Wiley & Sons, 2006).
  * Zaheer and Zubairy [1988] K. Zaheer and M. S. Zubairy, Phys. Rev. A 37, 1628 (1988).
  * Klimov  _et al._ [2001] A. B. Klimov, J. L. Romero, and C. Saavedra, Phys. Rev. A 64, 063802 (2001).
  * Larson [2007] J. Larson, Phys. Scr. 76, 146 (2007).
  * Irish [2007] E. K. Irish, Phys. Rev. Lett. 99, 173601 (2007).
  * Werlang  _et al._ [2008] T. Werlang, A. V. Dodonov, E. I. Duzzioni, and C. J. Villas-Bôas, Phys. Rev. A 78, 053805 (2008).
  * Zueco  _et al._ [2009] D. Zueco, G. M. Reuther, S. Kohler, and P. Hänggi, Phys. Rev. A 80, 033846 (2009).
  * Agarwal [1971] G. S. Agarwal, Phys. Rev. A 4, 1778 (1971).
  * Knight and Allen [1973] P. L. Knight and L. Allen, Phys. Rev. A 7, 368 (1973).
  * Friedberg and Manassah [2008] R. Friedberg and J. T. Manassah, Phys. Lett. A 372, 2514 (2008).
  * Fleming  _et al._ [2010] C. Fleming, N. I. Cummings, C. Anastopoulos, and B. L. Hu, J. Phys. A Math. Theor. 43, 405304 (2010).
  * Jørgensen and Wubs [2022] M. A. Jørgensen and M. Wubs, J. Phys. B: At. Mol. Opt. Phys. 55, 195401 (2022).
  * Wei  _et al._ [2021] Y.-C. Wei, M.-W. Lee, P.-T. Chou, G. D. Scholes, G. C. Schatz, and L.-Y. Hsu, J. Phys. Chem. C 125, 18119 (2021).
  * Wang  _et al._ [2020a] S. Wang, G. D. Scholes, and L.-Y. Hsu, J. Chem. Phys. 11, 5948 (2020a).
  * Hsu  _et al._ [2017] L.-Y. Hsu, W. Ding, and G. C. Schatz, J. Chem. Phys. 8, 2357 (2017).
  * Gruner and Welsch [1996] T. Gruner and D.-G. Welsch, Phys. Rev. A 53, 1818 (1996).
  * Dung  _et al._ [1998] H. T. Dung, L. Knöll, and D.-G. Welsch, Phys. Rev. A 57, 3931 (1998).
  * Scheel  _et al._ [1998] S. Scheel, L. Knöll, and D.-G. Welsch, Phys. Rev. A 58, 700 (1998).
  * Dzsotjan  _et al._ [2011] D. Dzsotjan, J. Kästel, and M. Fleischhauer, Phys. Rev. B 84, 075419 (2011).
  * Varguet  _et al._ [2021] H. Varguet, A. A. Díaz-Valles, S. Guérin, H. R. Jauslin, and G. Colas des Francs, J. Chem. Phys. 154, 084303 (2021).
  * Buhmann [2012a] S. Y. Buhmann,  _Dispersion Forces I_, Springer Tracts in Modern Physics, Vol. 247 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012).
  * Power  _et al._ [1959] E. A. Power, S. Zienau, and H. S. W. Massey, Phil. Trans. R. Soc. A 251, 427 (1959).
  * Weisskopf and Wigner [1930] V. Weisskopf and E. Wigner, Z. Phys. 63, 54 (1930).
  * Svidzinsky  _et al._ [2010] A. A. Svidzinsky, J.-T. Chang, and M. O. Scully, Phys. Rev. A 81, 053821 (2010).
  * Dung  _et al._ [2000] H. T. Dung, L. Knöll, and D.-G. Welsch, Phys. Rev. A 62, 053804 (2000).
  * Wang  _et al._ [2019] S. Wang, G. D. Scholes, and L.-Y. Hsu, J. Chem. Phys. 151, 014105 (2019).
  * Wang  _et al._ [2022] S. Wang, Y.-T. Chuang, and L.-Y. Hsu, J. Chem. Phys. 157, 184107 (2022).
  * Buhmann [2012b] S. Y. Buhmann,  _Dispersion Forces II_, Springer Tracts in Modern Physics, Vol. 248 (Springer Berlin Heidelberg, 2012).
  * Scheel  _et al._ [1999] S. Scheel, L. Knöll, and D.-G. Welsch, Phys. Rev. A 60, 4094 (1999).
  * L. Novotny and Hecht [2012] L. Novotny and B. Hecht,  _Principle of Nano-Optics_ (Cambridge University Press, New York, 2012).
  * Dung  _et al._ [2002] H. T. Dung, L. Knöll, and D.-G. Welsch, Phys. Rev. A 65, 043813 (2002).
  * Chew [1995] W. C. Chew,  _Waves and fields in inhomogeneous media_ , Vol. 16 (Inst of Electrical &, 1995).
  * Lee  _et al._ [2021] M.-W. Lee, Y.-T. Chuang, and L.-Y. Hsu, J. Chem. Phys. 155, 074101 (2021).
  * Wang  _et al._ [2020b] S. Wang, G. D. Scholes, and L.-Y. Hsu, The Journal of Physical Chemistry Letters 11, 5948 (2020b).
  * Chuang  _et al._ [2022] Y.-T. Chuang, S. Wang, and L.-Y. Hsu, J. Chem. Phys. 157, 234109 (2022).
  * I.S. Gradshteyn [2015] I. R. I.S. Gradshteyn,  _Table of Integrals, Series, and Products_ (Elsevier, 2015).