Understanding the Nature of Mean-Field Semiclassical Light-Matter Dynamics: An Investigation of Energy Transfer, Electron-Electron Correlations, External Driving and Long-Time Detailed Balance

TL;DR

This paper compares different semiclassical light-matter interaction models to determine their reliability in simulating energy transfer and driven dynamics in molecular systems, revealing that simpler mean-field approaches can outperform full quantum treatments in some cases.

Contribution

It systematically evaluates multiple semiclassical schemes for light-matter interactions, highlighting their strengths and limitations in modeling energy transfer and external driving effects.

Findings

Hamiltonian #I performs best for resonance energy transfer.

Hamiltonian #II is stable but fails at short distances.

Full quantum treatment can overestimate energy under external driving.

Abstract

Semiclassical electrodynamics is an appealing approach for studying light-matter interactions, especially for realistic molecular systems. However, there is no unique semiclassical scheme. On the one hand, intermolecular interactions can be described instantaneously by static two-body interactions connecting different molecules plus a classical transverse E-field; we will call this Hamiltonian #I. On the other hand, intermolecular interactions can also be described as effects that are mediated exclusively through a classical one-body E-field without any quantum effects at all (assuming we ignore electronic exchange); we will call this Hamiltonian #II. Moreover, one can also mix these two Hamiltonians into a third, hybrid Hamiltonian, which preserves quantum electron-electron correlations for lower excitations but describes higher excitations in a mean-field way. To investigate which semiclassical scheme is most reliable for practical use, here we study the real-time dynamics of a pair of identical two-level systems (TLSs) undergoing either resonance energy transfer (RET) or collectively driven dynamics. While all approaches perform reasonably well when there is no strong external excitation, we find that no single approach is perfect for all conditions. Each method has its own distinct problems: Hamiltonian #I performs best for RET but behaves in a complicated manner for driven dynamics. Hamiltonian #II is always stable, but obviously fails for RET at short distances. One key finding is that, under externally driving, a full configuration interaction description of Hamiltonian #I strongly overestimates the long-time electronic energy, highlighting the not obvious fact that, if one plans to merge quantum molecules with classical light, a full, exact treatment of electron-electron correlations can actually lead to worse results than a simple mean-field treatment.

Figures

Tables

Table 1. Table 1: Synopsis of the main features of the semiclassical Hamiltonians for modeling light-matter interactions

Approach	Definition	Quantum e-e correlations	Computational complexity	Causality
#I FCI	Eqs. (4)-(5)	Fully accounted	$O\left(2^N\right)$	Violated
#I CIS	Eqs. (4)-(6a)	Partially accounted	$O\left(N^2\right)$	Violated
#II	Eqs. (9)-(10)	None	$O\left(N\right)$	Preserved
hybrid	Eqs. (13)-(16)	Partially accounted	$O\left(N^2\right)$	Violated

Equations

$$\frac{d}{dt}|{\Psi }_N(t)⟩=-\frac{i}{\hslash }{\hat{H}}_{\text{sc}}|{\Psi }_N(t)⟩.$$

$$\frac{\partial }{\partial t}\mathbf{B}(\mathbf{r},t)$$

$$\frac{\partial }{\partial t}\mathbf{E}(\mathbf{r},t)$$

$$\mathbf{J}(\mathbf{r},t)=\sum _{n=1}^N\frac{\partial }{\partial t}\text{Tr}\left(\hat{\rho }(t){\hat{\mathcal{P}}}^{(n)}(\mathbf{r})\right)$$

$${\hat{H}}_{sc}^I=\sum _{n=1}^N{\hat{H}}_s^{(n)}-\int d\mathbf{r}{\mathbf{E}}_{\perp }(\mathbf{r},t)\cdot {\hat{\mathcal{P}}}^{(n)}(\mathbf{r})+\sum _{n<l}{\hat{V}}_{\text{Coul}}^{(nl)}$$

$${\hat{V}}_{\text{Coul}}^{(nl)}=\frac{1}{{ϵ}_0}\int d\mathbf{r}{\hat{\mathcal{P}}}_{\parallel }^{(n)}(\mathbf{r})\cdot {\hat{\mathcal{P}}}_{\parallel }^{(l)}(\mathbf{r})$$

$$|{\Psi }_N(t)⟩\approx |{\Psi }_{\text{CIS}}(t)⟩=\sum _I{C}_I(t)|{\Psi }_I⟩$$

$$|{\Psi }_I⟩=D_{0,I}|{\Psi }_0^{\text{HF}}⟩+\sum _{i=L}^{\frac{N_e}{2}}\sum _{a=\frac{N_e}{2}+1}^M{D}_{I,i}^a|{\Psi }_i^a⟩$$

$$\frac{d}{dt}C_{I}(t)=-\frac{i}{\hbar}\sum_{J}\left\langle\Psi_{I}\bigg{|}\hat{H}_{\text{sc}}^{I}\bigg{|}\Psi_{J}\right\rangle C_{J}(t)$$

$${\hat{V}}_{\text{Coul}}^{(nl)}$$

$$\approx \frac{1}{{ϵ}_0}\int d\mathbf{r}\left[{\mathcal{P}}_{\parallel }^{(n)}(\mathbf{r},t)\cdot {\hat{\mathcal{P}}}_{\parallel }^{(l)}(\mathbf{r})+{\mathcal{P}}_{\parallel }^{(l)}(\mathbf{r},t)\cdot {\hat{\mathcal{P}}}_{\parallel }^{(n)}(\mathbf{r})\right]\text{(up to a constant)}$$

$${\hat{H}}_{sc}^{II}=$$

$${\hat{H}}_{\text{MF}}^{(n)}=$$

$$+\frac{1}{{ϵ}_0}\int d\mathbf{r}{\mathcal{P}}_{\parallel }^{(n)}(\mathbf{r},t)\cdot {\hat{\mathcal{P}}}_{\parallel }^{(n)}(\mathbf{r})$$

$$|{\Psi }_N(t)⟩=|{\psi }_1(t)⟩|{\psi }_2(t)⟩\cdots |{\psi }_N(t)⟩$$

$$\frac{d}{dt}|{\psi }_n(t)⟩=-\frac{i}{\hslash }{\hat{H}}_{\text{MF}}^{(n)}|{\psi }_n(t)⟩$$

$${\hat{H}}_{sc}={\hat{H}}_{sc}^{II}+\hat{\mathcal{Q}}\delta {\hat{V}}_{\text{Coul}}\hat{\mathcal{Q}}$$

$$\delta {\hat{V}}_{\text{Coul}}\equiv$$

$$+{\mathcal{P}}_{\parallel }^{(l)}(\mathbf{r},t)\cdot {\hat{\mathcal{P}}}_{\parallel }^{(n)}(\mathbf{r})]$$

$$\hat{\mathcal{Q}}\equiv \sum _{i\in W}|i⟩⟨i|$$

$$W_{0+1}=\{\text{ground or single excitonic states}\}$$

$$|{\Psi }_N(t)⟩=|{\Psi }_{\text{CIS}}(t)⟩\otimes |{\psi }_{he}(t)⟩$$

$${\hat{H}}_s^{(n)}=\hslash {\omega }_0{\hat{\sigma }}_{+}^{(n)}{\hat{\sigma }}_{-}^{(n)}$$

$${\hat{\mathcal{P}}}^{(n)}(\mathbf{r})={\mu }_{ge}{\mathbf{e}}_d^{(n)}\delta (\mathbf{r}-{\mathbf{r}}_n){\hat{\sigma }}_x^{(n)}$$

$$k_{\text{FGR}}=\frac{{\omega }_0^3{\mu }_{ge}^2}{3\pi {ϵ}_0{c}^3\hslash }$$

$$D_{\mu i,\nu j}=\left\langle\Psi_{N}\Big{|}\hat{a}^{\dagger}_{i\mu}\hat{a}_{j\nu}\Big{|}\Psi_{N}\right\rangle$$

$$U_{\text{s}}=\sum _{n=1}^N\text{Tr}\left(\hat{\rho }(t){\hat{H}}_s^{(n)}\right)$$

$$k_{\text{Eh}}={\rho }_{gg}{k}_{\text{FGR}}$$

$${\mathbf{f}}_{\parallel }(\mathbf{r})=\int d{\mathbf{r}}^{\prime }{\overleftrightarrow{\mathbf{\delta }}}_{\parallel }(\mathbf{r}-{\mathbf{r}}^{\prime })\mathbf{f}({\mathbf{r}}^{\prime })$$

$${\overleftrightarrow{\mathbf{\delta }}}_{\parallel }(\mathbf{r})=\sum _{i,j=x,y,z}{\delta }_{\parallel ij}(\mathbf{r}){\mathbf{e}}_i{\mathbf{e}}_j$$

$${\delta }_{\parallel ij}(\mathbf{r})$$

$$=\frac{1}{3}{\delta }_{ij}\delta (\mathbf{r})-\frac{\eta (\mathbf{r})}{4\pi |\mathbf{r}{|}^3}\left(\frac{3r_i{r}_j}{|\mathbf{r}{|}^2}-{\delta }_{ij}\right)$$

$${\overleftrightarrow{\mathbf{\delta }}}_{\perp }(\mathbf{r})=\sum _{i,j=x,y,z}{\delta }_{\perp ij}(\mathbf{r}){\mathbf{e}}_i{\mathbf{e}}_j$$

Full text

I Introduction

Recent experiments demonstrating collective phenomena with nanoscale light-matter interactionsTörmä and Barnes (2015); Jia et al. (2018); Schäfer et al. (2019) have highlighted the need for computational simulations of realistic molecular systemsThakkar et al. (2015); Du et al. (2018); Flick et al. (2017). Unfortunately, full quantum electrodynamical (QED) calculations scale unfavorably with the number of quantized photonic modes. Moreover, full QED is compatible only with full configuration interactions (CI) for the description of the matter system, such that QED also scales unfavorably with the number of molecules. Thus, mixed quantum–classical electrodynamics are a promising approach with reduced computational cost: one treats electronic/molecular subsystems with approximate quantum mechanics and describes light fully classically. For decades, semiclassical electrodynamical simulations have captured many exciting phenomena in the field of quantum optics and spectroscopyGross and Haroche (1982); Puthumpally-Joseph et al. (2014, 2015); Sukharev and Nitzan (2011); Neuhauser and Lopata (2007); Lopata and Neuhauser (2009a, b); Lisinetskaya and Mitrić (2014).

Nevertheless, semiclassical electrodynamics suffers from many well-known issues. First, vacuum fluctuations are ignored due to the classical treatment of EM fields, which are usually calculated via a mean-field (Ehrenfest) approximation. Owing to the failure of the classical EM field description, semiclassical electrodynamics cannot fully recover any pure quantum effects for a single electronic system, such as spontaneous emissionCrisp and Jaynes (1969); Milonni (1976). To date, many techniques have been proposed to account for this issue and this work is still continuingFlick et al. (2017); Li et al. (2018a); Chen et al. (2019); Hoffmann et al. (2019, 2018).

The second issue for semiclassical electrodynamics is that there is no unique semiclassical Hamiltonian, and inevitably, some inconsistency must arise because of the semiclassical ansatz. After all, how should we treat electron-electron interactions? Are they instantaneous and static? Are they mediated exclusively by the EM field or not? If one chooses a static picture, one assumes the electronic Hamiltonian is a combination of quantum two-body terms plus an electric dipole coupling term (which defines Hamiltonian #I in Ref. 21); here one finds that one can predict accurate short-range RET rate but at the cost of violating the long-range causality due to a quantum-classical mismatch of intermolecular interactionsLi et al. (2018b). By contrast, if one chooses to couple matter exclusively through the field, one assumes the electronic Hamiltonian will have only an extended dipole coupling term (which defines Hamiltonian #II); here one finds that one fails to capture any short-range RET rate quantitatively due to the lack of quantum electron-electron correlations but strictly preserves causalityLi et al. (2018b).

The problem of correlation versus causality is usually ignored in the literature. Nowadays, almost all calculations use Hamiltonian #II at the cost of inaccurate short-range interactionsZiolkowski et al. (1995); Sukharev and Nitzan (2011); DiLoreto and Rangan (2018); Jestädt et al. (2018). Of course, one means of improving Hamiltonian #II is to use density functional theory (DFT) for electronic structure. In principle, DFT/TD-DFT can give exact electronic structure while maintaining the single-body nature of the electronic Hamiltonian. Beyond DFT, however, there is very little work using explicitly correlated electronic wavefunctions interacting with both external and internal EM fields. In general, if one wants to use explicitly correlated electronic wavefunctions (to account for electron-electron correlations) while studying light-matter interactions, to date the usual premise has been to first diagonalize a molecular electronic Hamiltonian (with no explicit electric field but rather only with instantaneous Coulomb terms) and then allow the resulting many-body electronic states to interact with an external electric fieldKlamroth (2003); Krause et al. (2005); Greenman et al. (2010); Tremblay et al. (2008); Rohringer et al. (2006). As such, the electronic dynamics as induced by internally generated, dynamic electric fields is not usually accounted for. As a result, almost all standard approaches fail to capture some key effects of collective phenomena, for example modification of the spontaneous decay rate, the effect of the dielectric constant, or even the presence of an RET rateScheel et al. (1999); Lo et al. (2001); Cortes and Jacob (2018). For our purposes, we will not invoke DFT in the present paper, and our goal is to establish a clean benchmark of Hamiltonian #I and #II dynamics, and distinguish between purely mean-field electronic dynamics and explicitly correlated electronic dynamics. We will attempt to answer the following equations:

1. (i)

By including quantum electron-electron correlations, is Hamiltonian #I always superior to Hamiltonian #II in practice? 2. (ii)

Can we always improve semiclassical results for Hamiltonian #I by treating quantum electron-electron correlations at a higher level of accuracy? For example, in the context of a Hartree-Fork (HF) ground state and configuration interaction singles (CIS) excited states, does the performance always improve if we increase the size of our configuration interaction (CI) Hamiltonian to include higher excited CIs (e.g., doubly excited CIs)?

In order to answer these questions, we will investigate resonant energy transfer (RET) and collectively driven electronic dynamics for a minimalistic two-site model within the framework of mean-field Ehrenfest dynamics. While we have previously applied the same model to study the short-time RET rateLi et al. (2018b), we will now study long-time RET dynamics as well as the crucial effects of including an external driving field (using a standard dyadic Green’s function technique; see Appendix). Understanding this minimal model should pave the way for improving the currently available semiclassical methods.

This paper is organized as follows. In Sec. II, we introduce the framework of mean-field Ehrenfest dynamics as well as different semiclassical Hamiltonians. In Sec. III, we introduce the model and parameters for simulations. In Sec. IV, we present results for RET and driven dynamics, showing that some unexpected, anomalous behavior can emerge. In Sec. V, we explain the reasons for this reported anomaly. We conclude in Sec. VI.

II Method: Semiclassical Electrodynamics

As a brief review, we will now review the conventional semiclassical method — mean-field (Ehrenfest) dynamics — for propagating light-matter electrodynamics. First, according to which the matter side obeys the time-dependent Schrödinger equation,

[TABLE]

Here, $\lvert\Psi_{N}\rangle$ denotes the electronic wave function for $N$ molecules, and $\hat{H}_{\text{sc}}$ denotes the semiclassical Hamiltonian, which will be introduced later. Second, for the EM side, the classical Maxwell’s equations are evolved:

[TABLE]

where $\epsilon_{0}$ denotes the vacuum permittivity. Here, the current density $\mathbf{J}(\mathbf{r},t)$ is calculated by a mean-field approximation:

[TABLE]

Here, $\hat{\boldsymbol{\mathcal{P}}}^{(n)}$ denotes the polarization density operator for molecule $n$. Eqs. (1)-(3) are called the coupled Maxwell-Schrödinger equations. In this framework, the only remaining question is how to define the form of $\hat{H}_{\text{sc}}$.

II.1 Hamiltonian #I

For neutral and non-overlapping molecules that interact with the E-field solely, the standard semiclassical Hamiltonian readsMukamel (1999)

[TABLE]

Here, $\hat{H}_{s}^{(n)}$ denotes the molecular Hamiltonian for molecule $n$; molecules interact with each other through a classical transverse E-field $\mathbf{E}_{\perp}$, and electron-electron correlations between molecules are characterized by the intermolecular Coulomb operator

[TABLE]

We note that $\hat{V}_{\text{Coul}}^{(nl)}$ scales as $1/R^{3}$ (where $R$ denotes intermolecular separations), $\int d\mathbf{r}\ \mathbf{E}_{\perp}(\mathbf{r},t)\cdot\hat{\boldsymbol{\mathcal{P}}}^{(n)}(\mathbf{r})$ scales as $1/R$. Thus, $\hat{V}_{\text{Coul}}^{(nl)}$ dominates short-range intermolecular interactions, while $\mathbf{E}_{\perp}(\mathbf{r},t)$ dominates long-range intermolecular interactions. For usual Förster resonance energy transfer (FRET)Förster (1948), we usually account only for $\hat{V}_{\text{Coul}}^{(nl)}$, leading to a $1/R^{6}$ dependence of the energy transfer rate (which follows from a Fermi’s golden rule calculation).

According to Eqs. (4) and (5), the exchange operator between molecules is neglected, which is adequate when the wave functions between molecules do not overlap. In this manuscript, we will call Eq. (4) Hamiltonian #I. In general, for $N$ TLSs, the quantum two-body term $\hat{V}_{\text{Coul}}^{(nl)}$ introduces a great deal of the computational complexity. Hamiltonian #I formally should require a Hilbert space of size $2^{N}$. Thus, in practice, when modeling electrodynamics, one is forced to construct approximations to Hamiltonian #I of which there are many. We will now investigate two such variants with different electronic structure theories to propagate the time-dependent Schrödinger equation.

II.1.1 Time-Dependent Full Configuration Interaction

To fully account for $\hat{V}_{\text{Coul}}^{(nl)}$, if one has the means, one can propagate the time-dependent Schrödinger equation in a complete basis using Hamiltonian #I. These exact, molecular quantum dynamics are known as time-dependent full configuration interaction (TD-FCI). Obviously, TD-FCI is possible only for simple models, i.e., a few two-level systems (TLSs), as $2^{N}$ grows fast for large $N$.

II.1.2 Time-Dependent Configuration Interaction Singles

For large systems with many molecules, in order to reduce the computational cost, the most common treatment is to truncate Hamiltonian #I at the level of single excitations, also called the time-dependent configuration interaction singles (TD-CIS) method. Here, the time-dependent electronic wave function is expanded as:

[TABLE]

Given the CIS wavefunction that is defined in Eq. (6), one can propagate the wavefunction as

[TABLE]

where $\hat{H}_{\text{sc}}^{I}$ is already defined in Eq. (4).

II.2 Hamiltonian #II

Even simpler than TD-CIS, a more radical solution is to invoke the mean-field approximation (or Hartree approximation) for $\hat{V}_{\text{Coul}}^{(nl)}$ [Eq. (5)]:

[TABLE]

where ${\boldsymbol{\mathcal{P}}}^{(n)}_{\parallel}(\mathbf{r},t)$ denotes the longitudinal component of the classical polarization density for molecule $n$. Keen readers might well be confused about the mean-field treatment in Eq. (8a): Why not take $\hat{V}_{\text{Coul}}^{(nl)}\approx\frac{1}{2\epsilon_{0}}\int d\mathbf{r}\left[{\boldsymbol{\mathcal{P}}}^{(n)}_{\parallel}(\mathbf{r},t)\cdot\hat{\boldsymbol{\mathcal{P}}}^{(l)}_{\parallel}(\mathbf{r})+{\boldsymbol{\mathcal{P}}}^{(l)}_{\parallel}(\mathbf{r},t)\cdot\hat{\boldsymbol{\mathcal{P}}}^{(n)}_{\parallel}(\mathbf{r})\right]$ instead? The motivation behind Eq. (8a) is twofold: (i) Eq. (8a) allows us to define a semiclassical Hamiltonian that strictly preserves causality, as is shown below; (ii) the mean-field expansion in Eq. (8a) is already standard in the area of many-body physics; see Ref. 35 for a brief introduction. Because the last term in Eq. (8a) is just a time-dependent constant and will not alter the equations of motion for the molecular part, this term can be further neglected, leading to Eq. (8b).

By substituting Eq. (8b) into Eq. (4), we arrive at Hamiltonian #II:

[TABLE]

where

[TABLE]

Within Hamiltonian #II, molecules interact with each other only through a classical E-field $\mathbf{E}(\mathbf{r},t)$, and the last term above [in Eq. (10)] denotes the semiclassical self-polarization, which effectively renormalizes the energy levels of molecules slightly and does not significantly alter the overall dynamics. Hence, we will neglect the last term in our numerical simulations. Keen readers might wonder whether energy conservation is still valid if the last term is neglected — indeed, energy conservation can be guaranteed if we simply redefine the conserved quantityLi et al. (2018b).

II.2.1 Time-dependent Hartree Method

Given the one-body nature of Hamiltonian #II, the time-dependent Schrödinger equation can be evolved exactly with simple time-dependent Hartree (TDH) dynamics, i.e., the electronic wavefunction can be written as a Hartree product:

[TABLE]

where $\lvert\psi_{n}(t)\rangle$ denotes an effective one-body wavefunction for molecule $n=1,2,\cdots,N$. Following the variational principleBeck (2000), the equation of motion for each orbital $\lvert\psi_{n}(t)\rangle$ can be obtained as:

[TABLE]

where $\hat{H}^{(n)}_{\text{MF}}$ is defined in Eq. (10), and ${\boldsymbol{\mathcal{P}}}^{(n)}_{\parallel}(\mathbf{r},t)=\langle\psi_{n}(t)\rvert\hat{\boldsymbol{\mathcal{P}}}_{\parallel}^{(n)}(\mathbf{r})\lvert\psi_{n}(t)\rangle$.

II.3 A Hybrid Hamiltonian

While Hamiltonian #I treats only the transverse E-field classically, Hamiltonian #II treats all intermolecular interactions classically. Interestingly, we can write both of these Hamiltonians in a uniform way:

[TABLE]

where $\delta\hat{V}_{\text{Coul}}$ is defined as

[TABLE]

and $\hat{\mathcal{Q}}$ denotes a projection operator into a subspace ($W$) of the electronic states:

[TABLE]

When $W=\{\emptyset\}$, $\hat{\mathcal{Q}}=\mathbb{0}$, and Eq. (13) reduces to Hamiltonian #II; when $W$ is the entire electronic manifold of states $\mathcal{S}$, $Q=\mathds{1}$, and Eq. (13) reduces to Hamiltonian #I. By choosing an arbitrary subspace in between $\{\emptyset\}$ and $\mathcal{S}$, we can find intermediate Hamiltonians in between Hamiltonians #I and #II. Hence, Eqs. (13)-(15) form a generalized definition of a semiclassical Hamiltonian. Clearly, the choice of the subspace will play an important role in the quality of the Hamiltonian. In this manuscript, we define one intermediate subspace as:

[TABLE]

We call Eqs. (13)-(16) a hybrid Hamiltonian ($\hat{H}_{\text{sc}}^{hyb}$), in which there are two-body couplings ($\hat{V}_{\text{Coul}}^{(nl)}$) for the ground and singly excited states, but the double and higher excited states are entirely decoupled and reduced to mean-field interactions.

II.3.1 Time-dependent Hybrid Method

For the hybrid Hamiltonian, the many-body wave function can be expanded as follows:

[TABLE]

Here, $\lvert\Psi_{\text{CIS}}(t)\rangle$ characterizes the wave function for the CIS states, which is defined in Eq. (6), and $\lvert\psi_{he}(t)\rangle$ characterizes the wave function for higher excitations. On the one hand, we evolve $\lvert\Psi_{\text{CIS}}(t)\rangle$ by TD-CIS as in Eq. (7); on the other hand, because each higher excited state interacts with other states (i.e., CIS and other higher excited states) solely through a classical E-field, these states can be propagated independently with TDH as in Eq. (12). For example, for a pair of TLSs, the explicit form of the hybrid Hamiltonian is presented in Eq. (A32).

II.4 Summary of Semiclassical Hamiltonians

Fig. 1 is a cartoon of the four different semiclassical approaches (Hamiltonian #I FCI, #I CIS, #II, and a hybrid Hamiltonian) that have been introduced above. In this cartoon, we highlight how intermolecular interactions are described differently in these approaches.

Table 1 also summarizes the important features of these Hamiltonians, e.g., defining equations, whether or not quantum electron-electron correlations are accounted for, computational complexity as a function of molecular number ($N$), and whether long-range causality is preserved or not.

III Model

From now, natural units will be used: $[\hbar]=[c]=[\epsilon_{0}]=1$. We will perform calculations with a minimalistic quantum model — a pair of identical TLSs (labeled as D and A). The molecular Hamiltonian for molecule $n=D,A$ reads:

[TABLE]

where $\hbar\omega_{0}$ denotes the energy gap between ground state $\lvert ng\rangle$ and excited state $\lvert ne\rangle$ for molecule $n$, $\hat{\sigma}_{+}^{(n)}\equiv\lvert ne\rangle\langle ng\rvert$, and $\hat{\sigma}_{-}^{(n)}\equiv\lvert ng\rangle\langle ne\rvert$. After the long-wavelength approximation, $\hat{\boldsymbol{\mathcal{P}}}^{(n)}(\mathbf{r})$ reads:

[TABLE]

where $\hat{\sigma}_{x}^{(n)}=\lvert ng\rangle\langle ne\rvert+\lvert ne\rangle\langle ng\rvert$, $\mu_{ge}$ denotes the magnitude of the transition dipole moment, and $\mathbf{e}_{d}^{(n)}$ and $\mathbf{r}_{n}$ denote the unit vector along the transition dipole and the position of molecule $n$.

For our simulation parameters, we suppose that the TLSs are positioned symmetrically at $\mathbf{r}_{n}=(0,\pm\frac{R}{2},0)$, and their transition dipole moments are both oriented along the $z$-axis ($\mathbf{e}_{d}^{(n)}=\mathbf{e}_{z}$). We set $\omega_{0}=1$ and $\mu_{ge}=0.1$. In vacuum, the spontaneous emission rate for a single TLS is defined as

[TABLE]

With these parameters above, $k_{\text{FGR}}=1.6\times 10^{-3}$. To characterize the separation between TLSs, a dimensionless quantity $k_{0}R=\frac{\omega_{0}R}{c}$ is used. Finally, we will choose the intermolecular separation to be $k_{0}R=0.4$ (by default), corresponding to the dipole-dipole interaction $v_{dd}=\frac{\mu_{ge}^{2}}{4\pi\epsilon_{0}R^{3}}=1.2\times 10^{-2}$.

Since we operate in vacuum with no dielectric, we can calculate the time-dependent E-field by the dyadic Green’s function techniqueNovotny and Hecht (2006) instead of numerically solving Eq. (2) in a three-dimensional gridTaflove and Hagness (2005); see Appendix A for details. We numerically solve the reduced equation of motion for the molecular subsystem by a Runge-Kutta fourth-order propagatorButcher (2008) with the time step $\Delta t=0.01$.

IV Results

After introducing the model Hamiltonian and relevant dynamical methods, we will now perform simulations to mimic two different phenomena: (i) resonance energy transfer (RET) with no external EM field, and (ii) driven dynamics under an external driving cw field. For each case, four semiclassical treatments are considered: (i) Hamiltonian #I FCI; (ii) Hamiltonian #I CIS; (iii) Hamiltonian #II, and (iv) the hybrid Hamiltonian. To examine the performance of semiclassical approaches, we will compare them against either the time-dependent perturbative QED resultSalam (2010, 2018) or the results of Lehmberg-Agarwal master equation (LAME)Lehmberg (1970); Agarwal (1974) — the standard quantum approach for describing the dynamics of two-level systems (TLSs) in quantum optics; see Appendix B for details. Note that all LAME results presented below are calculated with FCI.

IV.1 Resonance Energy Transfer (RET)

For RET, no external driven field is considered. The donor (D) is initialized in a superposition state ($c_{g}\lvert Dg\rangle+c_{e}\lvert De\rangle$, where $|c_{g}|,|c_{e}|>0$), the acceptor (A) is initialized in the ground state. Here, we choose a superposition state for the donor so that we can initialize a time-dependent current density (and therefore EM field) without invoking any external EM fields. It is well known that Ehrenfest dynamics can depend (unphysically) on the initial state for the donor; for example, if $c_{e}=1$, Ehrenfest dynamics do not predict any spontaneous emission and are completely wrong. We consider two regimes: short-time dynamics, from which a RET rate ($k_{\text{ET}}$) can be extracted (see Appendix C for details), and long-time dynamics, in which dissipation effects become important.

RET rate

Fig. 2 plots the RET rate as a function of intermolecular separation ($k_{0}R$, where $k_{0}\equiv\frac{\omega_{0}}{c}$). Here, the perturbative QED calculation (black line) suggests that the RET rate obeys two mechanisms in different separation limits: at short range ($k_{0}R\ll 1$), the RET rate scales as $\frac{1}{R^{6}}$ due to dipole-dipole interactions, known as Förster resonance energy transfer (FRET)Förster (1948); at long range ($k_{0}R\gg 1$), the RET rate scales as $\frac{1}{R^{2}}$ because the transverse E-field dominates energy transferAndrews (1989). In general, all semiclassical approaches qualitatively predict these scalings but not quantitatively. For example, at short range, Hamiltonian #I [FCI (red circles) and CIS (blue triangles)] and the hybrid Hamiltonian (yellow squares) quantitatively agree with QED while Hamiltonian #II (cyan stars) predicts only a fraction of the true RET rate [proportional to the ground state population of the donor $\rho_{gg}^{(D)}(0)$; see Appendix C for an analytic proof].

At long range, not surprisingly, because all semiclassical approaches use a classical E-field and ignore vacuum fluctuations, none of the methods can predict the RET rate correctly (when $|c_{g}|\ll 1$). After all, in this limit, the correct physical mechanism is akin to spontaneous emission from one TLS followed by absorption by the other TLS. Interestingly, in this limit, Hamiltonian #I CIS predicts an RET rate with a larger error than Hamiltonian #I FCI; the underlying reason for this deterioration of accuracy is not obvious because, according to QED, excluding the doubly excited state should not alter the RET rate if the double is not populated initially (as is true for RET).

Long-time RET dynamics

Fig. 3 plots (from top to bottom) the long-time RET population dynamics for the donor and acceptor, as well as the impurity of the one-electron reduced density matrix (1-RDM) when the TLSs are close ($k_{0}R=0.4$). Here, the impurity of the 1-RDM is a measure to characterize how much the electronic states of different molecules are mixed. For example, when Hamiltonian #II is used, because the total wavefunction for a pair of TLSs can always be separated as a product of the wavefunctions for each TLS (which is certainly not true if other approaches are used), the impurity of 1-RDM is always zero (provided it starts at zero). Formally, the impurity of 1-RDM is calculated by $\text{Tr}\left(D\right)-\text{Tr}\left(D^{2}\right)$, where the matrix elements of the 1-RDM ($D$) are defined to be:

[TABLE]

Here, $\left\{\mu,\nu\right\}=\{1,2,\cdots,N\}$, $\{i,j\}=\{e,g\}$, $\hat{a}^{\dagger}_{\mu i}$ and $\hat{a}_{\mu i}$ are the creation and annihilation operators for state $\lvert\mu i\rangle$.

When the donor is weakly excited initially [$\rho_{ee}^{(D)}(0)=0.1$; left panel], all semiclassical approaches predict the same population dynamics (Figs. 3a,c) as the Lehmberg-Agarwal master equation (LAME, black line). These predictions agree with the consensus that a mean-field approximation should be valid when the donor is weakly excited, i.e., in the perturbative regime, where a classical E-field is good enough. When the donor is strongly excited [$\rho_{ee}^{(D)}(0)=0.9$; right panel], the semiclassical approaches can still predict some key features in population dynamics like oscillations (due to the dipole-dipole interaction), the dissipation, and the long-time slow decay of the dark state, but the dissipation rate is underestimated compared to LAME. In general, due to a lack of quantum dipole-dipole interactions, Hamiltonian #II (cyan solid) predicts slightly less accurate oscillation periods than other semiclassical approaches. More interestingly, for the impurity of the matter 1-RDM (Figs. 3e,f), we find that the more one properly accounts for quantum dipole-dipole interactions, the larger is the impurity of the matter subsystem as predicted by semiclassical dynamics [i.e., as far as the impurity of the matter subsystem, LAME $>$ Hamiltonian #I FCI(CIS) $>$ the hybrid Hamiltonian $>$ Hamiltonian #II $=0$]. In Fig. 3f, LAME predicts an impurity around $\frac{1}{2}$ at the long times, which can be understood as follows: for a pair of TLSs in vacuum, if the donor is fully excited, the final state for the TLSs plus the photonic field should be $\frac{1}{\sqrt{2}}\lvert gg;1\rangle+\frac{1}{\sqrt{2}}\lvert d;0\rangle$, where $\lvert gg;1\rangle$ denotes the TLSs in the ground state plus an emitted photon and $\lvert b;0\rangle$ denotes the TLSs in the bright state associated with no photon; thus, the corresponding reduced density matrix for the electronic degrees of freedom is $\sigma=\begin{pmatrix}\frac{1}{2}&0\\ 0&\frac{1}{2}\end{pmatrix}$, so that the impurity is $\frac{1}{2}$. By contrast, the fact that Ehrenfest is too pure (with an impurity much smaller than LAME) is a statement that additional decoherence is needed.

Apart from the RET dynamics of the two-level molecules, it is also worthwhile to study the frequency of the scattered E-field. Fig. 4 plots the frequency of the scattered E-field as a function of $\rho_{ee}^{(D)}(0)$ for RET dynamics. As predicted by LAME, the frequency of the E-field should not depend on $\rho_{ee}^{(D)}(0)$. However, we find that Hamiltonian #II disagrees with LAME and shows a slightly nonphysical behavior when $\rho_{ee}^{(D)}(0)$ gradually increases; by contrast, all other semiclassical approaches agree with LAME relatively well.

From the above RET results, we gather that Hamiltonian #II is slightly less accurate than the other semiclassical approaches especially when the donor becomes more than weakly excited, in which case one should include quantum dipole-dipole interactions.

IV.2 Collectively Driven Dynamics

Now, let us move to the case of collectively driven dynamics for a pair of TLSs prepared initially in ground state. The incident cw field takes the following form: $\mathbf{E}_{\text{in}}(\mathbf{r},t)=E_{0}\sin(\omega_{0}t-k_{0}x)\mathbf{e}_{z}$. To characterize the strength of the cw field, the Rabi frequency ($\Omega\equiv\mu_{ge}E_{0}$) is a good indicator: $\Omega<k_{\text{FGR}}$ ($\Omega>k_{\text{FGR}}$) represents a weak (strong) driving field. In general, for closely aggregated TLSs ($k_{0}R\ll 1$), because the spontaneous emission rate is strongly modified by intermolecular interactions ($v_{dd}$) instead of the vacuum value in Eq. (20), one would expect that semiclassical approaches should be valid as long as the Rabi frequency is much smaller than the dipole-dipole coupling ($\Omega\ll v_{dd}$). With this in mind, we check the results of driven dynamics as below.

Weakly driven dynamics

Fig. 5 plots the electronic energy (Fig. 5a-b) and the impurity of 1-RDM (Fig. 5c-d) for a pair of TLSs driven by a weak cw field ($\Omega=0.3k_{\text{FGR}}$) at both short times (left panel) and long times (right panel). Here, the electronic energy of the molecular subsystem ($U_{\text{s}}$) is defined as

[TABLE]

As explained above, we expect that all approaches (Hamiltonian #I FCI, #I CIS, #II, the hybrid Hamiltonian, and LAME) should predict the same dynamics for electronic energy. The surprising finding, however, is that after very long times ($t>200k_{\text{FGR}}^{-1}$), Hamiltonian #I FCI (red solid) predicts an unphysically large electronic energy compared to other approaches; see Fig. 5b. This unphysical behavior indicates (ironically) that a full accounting for quantum electron-electron correlations can actually be problematic even in the weak coupling limit. The reason for this anomaly will be addressed in Sec. V. For the impurity of the 1-RDM, as shown in Fig. 5d, while Hamiltonian #I FCI overestimates the impurity as compared with LAME, the hybrid Hamiltonian predicts similar steady-state impurity as LAME, and other semiclassical approaches predict much less impurity than LAME (note that here Hamiltonian #II still always predicts zero impurity).

Strongly driven dynamics

Fig. 6 plots the dynamics of the electronic energy and the impurity of 1-RDM when the cw field becomes stronger ($\Omega=2.0k_{\text{FGR}}<v_{dd}$). In this limit, because the contribution of the double is not negligible, as is shown in Fig. 6b(d), LAME predicts a much higher steady-state electronic energy (and impurity) than does Hamiltonian #I CIS, for which the double is truncated. Just as in Fig. 5, by including the double, Hamiltonian #I FCI overestimates the electronic energy significantly compared with LAME, reinforcing the notion that fully accounting for electron-electron correlation can be problematic (in both the weak and strong field limits). As far as the impurity of the matter 1-RDM (Fig. 6c,d.), the behaviors of the different semiclassical approaches are similar to what was found in the case of electronic energy, except for the fact that Hamiltonian #II always predicts zero impurity. Apparently, electronic FCI coupled to a classical EM field can predict nonphysical features, which conflicts with our intuition that including more electron-electron correlations should give better results.

Overall, for a reasonably strong field, no semiclassical approach can predict the steady-state electronic energy or the impurity of the matter 1-RDM correctly, which would naively conflict with the general consensus that semiclassical electrodynamics should be valid as long as the Rabi frequency ($\Omega$) is much smaller than the strength of the dipole-dipole coupling ($v_{dd}$). The validity of semiclassical electrodynamics is obviously complicated, and must depend on which Hamiltonian one uses. With this in mind, let us now digest the results above and consider why FCI behaves so poorly in Figs. 5 and 6.

V Discussion

From the results above in Figs. 2-6, our general conclusion is that no semiclassical method is perfect, but Hamiltonian #I CIS and the hybrid Hamiltonian seem to perform optimally and they are reasonably computationally efficient. Hamiltonian #II performs slightly worse (failing for RET and the impurity of 1-RDM). The most stunning conclusion is the drastic failure of Hamiltonian #I FCI under driven dynamics.

To better understand the failure of FCI in driven dynamics, consider the steady-state data in Fig. 7a. When the Rabi frequency ($\Omega$, $x$-axis) is much smaller than the dipole-dipole coupling ($v_{dd}$, the vertical magenta line), all semiclassical approaches predict similar steady-state population for the singles ($y$-axis) as compared to LAME. However, when we investigate the population of the double excitation (see Fig. 7b), conventional semiclassical approaches fail even when $\Omega\ll v_{dd}$. On the one hand, Hamiltonian #II always greatly underestimates the population of the double. If we restrict ourselves to the weak coupling limit ($\Omega\ll k_{\text{FGR}}$), such an underestimation is not very problematic because the population is so small as to have minimal effect on any physical observables. On the other hand, Hamiltonian #I FCI always overestimates the population for the double, leading to a nonphysically large electronic energy; see Fig. 7c. One may therefore hypothesize that the inclusion of the doubly excited state represents an important but risky proposal for semiclassical electrodynamics; overestimation of the double population is strongly correlated to the overestimation of the total electronic energy. Interestingly, the hybrid Hamiltonian does interpolate between Hamiltonian #I and #II, but there is minimal gain in accuracy when $\Omega\gg k_{\text{FGR}}$.

We can now answer the question above: why does FCI fail and predict an exorbitant accumulation of energy for the TLSs under a driving force? The root of this problem is the classical EM field. Note that, for a single TLS, due to the use of a classical EM field, Ehrenfest dynamics predicts a decay rate proportional to the ground state populationCrisp and Jaynes (1969); Milonni (1976); Li et al. (2018a):

[TABLE]

For a pair of closely aggregated TLSs ($k_{0}R\ll 1$, as considered in this manuscript), if one neglects the effect of the dark state and focuses on a three-level system with ground state $\lvert 0\rangle$, bright state $\lvert b\rangle$, and doubly excited state $\lvert 2\rangle$, the allowed optical transitions are $\lvert 0\rangle\leftrightarrow\lvert b\rangle$ and $\lvert b\rangle\leftrightarrow\lvert 2\rangle$ [and the Ehrenfest decay rates between these optical transitions also obey Eq. (23)]. For driven dynamics, with system initially in state $\lvert 0\rangle$, the quantum dipole-dipole interaction $\hat{V}_{\text{Coul}}^{(nl)}$ directly couples state $\lvert 0\rangle$ and state $\lvert 2\rangle$. Now, suppose we apply Hamiltonian #I with FCI. On the one hand, with driven dynamics, $\hat{V}_{\text{Coul}}^{(nl)}$ leads to an increase of the population for state $\lvert 2\rangle$; on the other hand, because initially $\rho_{bb}(0)=0$, according to Eq. (23), the decay rate from $\lvert 2\rangle$ to $\lvert b\rangle$ is greatly suppressed. As a result, state $\lvert 2\rangle$ will continuously accumulate the population, leading to an unphysically large electronic energy even in the weak coupling limit. In short, the exaggerated electronic energy predicted by Hamiltonian #I FCI (see Figs. 5-7) appears to come directly from the mismatch of the quantum electron-electron correlations and the classical EM field. Interestingly, this mismatch also causes the violation of long-range causalityLi et al. (2018b). The above discussion should be very general, valid for a pair of TLSs or in the case of many molecules: for driven systems, the population dynamics for higher excitations (beyond singles) cannot be correctly described by Hamiltonian #I FCI even when the driving field is very weak.

VI Conclusion

To conclude, in this manuscript, we have applied different semiclassical approaches to a minimalistic many-site model for light-matter interactions — a pair of identical TLSs. We find: (i) For the impurity of the 1-RDM, generally no semiclassical approach agrees with LAME very well; (ii) For RET dynamics, Hamiltonian #II is not an optimal candidate due to a lack of quantum dipole-dipole couplings; (iii) For collectively driven dynamics, all semiclassical approaches in Table 1 can correctly describe the population of single states when the Rabi frequency is much smaller than the dipole-dipole coupling ($\Omega\ll v_{dd}$); (iv) For collectively driven dynamics, even when $\Omega\ll v_{dd}$, Hamiltonian #I FCI always predicts a nonphysically large double population (and thus an incorrect electronic energy) due to a mismatch between quantum electron-electron correlations and a classical E-field; (v) A hybrid Hamiltonian can eliminate the reported anomaly for #I FCI in the weak field as well as outperform Hamiltonian #II with regard to RET. Nevertheless, the accuracy of the hybrid Hamiltonian is still far from quantitative.

For the moment, when using semiclassical electrodynamics to describe light-matter interactions, our recommendation is to use Hamiltonian #I CIS or the hybrid Hamiltonian as a trade-off between accuracy and computational cost. We must emphasize that (i) our present benchmark work was restricted to only a pair of TLSs, and (ii) no semiclassical algorithm performs quantitatively at all. In the future, these limitations must be addressed. On the one hand, for a large collection of molecules, more exciting collective phenomena should emerge and the performances of the different semiclassical approaches must be tested. On the other hand, and even more importantly, it is also natural to ask whether or not further algorithmic improvements can be made to the semiclassical methods above. For example, can we include some crucial aspects of spontaneous emission that are missed in a mean-field treatment and improve Hamiltonian #II? Recent experience Chen et al. (2019) suggests such improvements are possible and this work in ongoing.

VII Acknowledgement

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0019397. The research of A.N. is supported by the Israel-U.S. Binational Science Foundation. This research also used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

Appendix A Analytical and EM-Free Form of Semiclassical Hamiltonians

A.1 Longitudinal and Transverse Components

For a vector function $\mathbf{f}(\mathbf{r})=f_{x}(\mathbf{r})\mathbf{e}_{x}+f_{y}(\mathbf{r})\mathbf{e}_{y}+f_{z}(\mathbf{r})\mathbf{e}_{z}$, the longitudinal component is defined by

[TABLE]

where the dyadic longitudinal $\delta$-function $\overleftrightarrow{\boldsymbol{\delta}}_{\parallel}(\mathbf{r})$ is

[TABLE]

Here, $\mathbf{e}_{i}$ denotes a unit vector along direction $i=x,y,z$, and

[TABLE]

While the first definition, Eq. (A3a), is a natural definition of the longitudinal $\delta$-function, this expansion diverges at $|\mathbf{r}|=0$. To avoid such divergence, regularization is introduced, leading to the second definition, Eq. (A3b), in which $\eta(\mathbf{r})\equiv 0$ at $\mathbf{r}=0$ to suppress the divergence and $\eta(\mathbf{r})\equiv 1$ elsewhereCohen-Tannoudji et al. (1997).

Similar to Eq. (A2), the dyadic transverse $\delta$-function reads

[TABLE]

Note that $\delta_{\perp ij}(\mathbf{r})\equiv\delta_{ij}(\mathbf{r})-\delta_{\parallel ij}(\mathbf{r})$, so that the transverse component $\mathbf{f}_{\perp}(\mathbf{r})$ can be calculated by

[TABLE]

According to the definitions of the longitudinal and transverse $\delta$-functions, it is easy to show that $\int d\mathbf{r}\ \mathbf{f}_{\perp}(\mathbf{r})\cdot\mathbf{f}_{\parallel}(\mathbf{r})=0$ for all vector fields $\mathbf{f}(\mathbf{r})$.

A.2 Time-Dependent Dyadic Green’s Functions

If we assume that the electronic subsystem couples only to the E-field (as is true in this manuscript), it is more convenient to rewrite Maxwell’s equations [Eq. (2)] as

[TABLE]

A formal solution of the E-field reads

[TABLE]

where $\mathbf{E}_{\text{in}}(\mathbf{r},t)$ denotes the incoming field, and $\mathbf{E}^{(n)}(\mathbf{r},t)$ denotes the E-field that is emitted by molecule $n$, which can be further evaluated through the time-dependent dyadic Green’s function techniqueNovotny and Hecht (2006), i.e.,

[TABLE]

where $V$ denotes the integral volume that includes ${\boldsymbol{\mathcal{P}}}^{(n)}$. The time-dependent dyadic Green’s function $\overleftrightarrow{\mathbf{G}}(\mathbf{r},\mathbf{r}^{\prime},t,t^{\prime})$ is defined as

[TABLE]

where $k=\frac{\omega}{c}$. For a point source in a homogeneous environment, the time-dependent scalar Green’s function $G_{0}$ reads

[TABLE]

where $n=1$ in vacuum. By substituting Eqs. (A9) and (A10) into Eq. (A8), we arrive at a retarded expression of $\mathbf{E}^{(n)}$:

[TABLE]

Now, very often, within the content of electrodynamics with retardation, it is helpful to work with the time-independent dyadic Green’s function $\overleftrightarrow{\mathbf{G}}(\mathbf{r},\mathbf{r}^{\prime})$:

[TABLE]

where $R\equiv|\mathbf{r}-\mathbf{r}^{\prime}|$. Eq. (A12) can be rewritten as

[TABLE]

where $\overleftrightarrow{\mathbf{\eta}}_{1}$ and $\overleftrightarrow{\mathbf{\eta}}_{3}$ are defined as

[TABLE]

and $\hat{\mathbf{R}}_{i}$ denotes the unit vector along the direction of $\mathbf{R}_{i}=\mathbf{r}_{i}-\mathbf{r}_{i}^{\prime}$.

Because we will have different molecules at different sites, let us also introduce the following short-hand writing

[TABLE]

where $R_{nl}\equiv|\mathbf{r}_{n}-\mathbf{r}_{l}|$, $\mathbf{e}_{d}^{(n)}$ denotes the unit vector along the orientation of dipole $n$, $\eta_{1}^{(nl)}=\mathbf{e}_{d}^{(n)}\cdot\overleftrightarrow{\mathbf{\eta}}_{1}\mathbf{e}_{d}^{(l)}$ and $\eta_{3}^{(nl)}=\mathbf{e}_{d}^{(n)}\cdot\overleftrightarrow{\mathbf{\eta}}_{3}\mathbf{e}_{d}^{(l)}$. $G_{nl}$ in Eq. (A15) characterizes the magnitude of the light-matter coupling between the two unit dipoles at sites $n$ and $l$. The real and imaginary parts of $G_{nl}$ read

[TABLE]

Interestingly, when two dipoles overlap, i.e., $\mathbf{e}_{d}^{(n)}=\mathbf{e}_{d}^{(l)}$, $\eta_{1}^{(nl)}=\eta_{3}^{(nl)}=1$, and $R_{nl}\rightarrow 0$, a Taylor expansion of Eq. (A16) to leading order in $kR$ reduces to

[TABLE]

A.3 Analytical and EM-Free Form of Hamiltonian #I

For Hamiltonian #I [defined in Eq. (4)], it is unnecessary to evaluate $\mathbf{E}_{\perp}$ at all times. Instead, for neutral molecules (with no free charge), since the displacement field ($\mathbf{D}$) is transverse, i.e., $\mathbf{D}_{\parallel}=\mathbf{E}_{\parallel}+\frac{1}{\epsilon_{0}}{\boldsymbol{\mathcal{P}}}_{\parallel}=\mathbf{0}$, $\mathbf{E}_{\perp}$ can be rewritten as

[TABLE]

By substituting Eq. (A18) into Eq. (4), one obtains another form for Hamiltonian #I:

[TABLE]

At this point, let us evaluate all of the terms in Eq. (A19). If we make the long wave approximation, i.e., $\hat{\boldsymbol{\mathcal{P}}}^{(n)}(\mathbf{r})=\hat{\boldsymbol{\mu}}^{(n)}\delta(\mathbf{r}-\mathbf{r}_{n})\mathbf{e}_{d}^{(n)}$ (where $\hat{\boldsymbol{\mu}}^{(n)}\equiv\mu_{ge}\mathbf{e}_{d}^{(n)}\hat{\sigma}_{x}^{(n)}$ denotes the transition dipole operator for TLS $n$), and apply Eqs. (A1)-(A3), $\hat{V}_{\text{Coul}}^{(nl)}$ [Eq. (5)] is reduced to the dipole-dipole interaction form:

[TABLE]

where $\otimes$ denotes the Kronecker tensor product. Similarly, for $n\neq l$, the last term in Eq. (A19) can be simplified as

[TABLE]

where $\boldsymbol{\mu}^{(l)}(t)\equiv\text{Tr}\left(\hat{\rho}(t)\hat{\boldsymbol{\mu}}^{(l)}\right)=2\real{\rho_{ge}^{(n)}(t)}\mu_{ge}\mathbf{e}_{d}^{(n)}$, and $\rho_{ge}^{(n)}(t)$ denotes the coherence between the ground state and excited state for TLS $n$. For Hamiltonian #I, we can calculate $\rho_{ge}^{(n)}$ by $\rho_{ge}^{(n)}(t)=\text{Tr}\left(\hat{\rho}(t)\hat{\sigma}_{+}^{(n)}\right)$.

At this point, having evaluated all electronic matrix elements in Eq. (A19), for the sake of simplicity and efficiency, we would like to completely reduce Hamiltonian #I (when possible) into a Hamiltonian operating only on the electronic degrees of freedom, from which the electric and magnetic fields can be extrapolated analytically; this is, after all, the framework of the famous optical Bloch equation (OBE). To do so, let us evaluate the E-field using a Green’s function technique. For a TLS under the long wavelength approximation, ${\boldsymbol{\mathcal{P}}}^{(n)}(\mathbf{r},t)$ reads

[TABLE]

By substituting Eq. (A22) into Eq. (A9), we arrive at an analytical form for $\mathbf{E}^{(n)}(\mathbf{r},t)$:

[TABLE]

Between Eq. (A23a) and Eq. (A23b), we have neglected all retardation and assumed

[TABLE]

the time-independent Green’s function $\overleftrightarrow{\mathbf{G}}(\mathbf{r},\mathbf{r}_{n})$ is defined in Eq. (A12). Given Eq. (A23), the coupling between molecule $n$ and the E-field generated by molecule $l$ ($n\neq l$) is expressed as

[TABLE]

where $G_{nl}=G_{ln}$ is defined in Eq. (A15), $\real{G_{nl}}$ and ${\text{Im}\left[G_{nl}\right]}$ are defined in Eq. (A16). $\theta\left(t-\frac{R_{nl}}{c}\right)$ denotes the Heaviside step function which is required to preserve causality.

Finally, using Eqs. (A18) and (A21), the transverse interaction between sites $n$ and $l$ reads

[TABLE]

For the case of $n=l$, we apply Eqs. (A17), (A25c), and (A21), noting that the terms involving $\real{\rho_{ge}(t)}$ cancel. Then Eq. (A26b) becomes

[TABLE]

where $k_{\text{FGR}}$ is defined in Eq. (20).

Thus, in the end, provided we can make the assumption in Eq. (A24), we have obtained an analytical and EM-free form of Hamiltonian #I:

[TABLE]

where the analytical expressions of $\hat{V}_{\text{Coul}}^{(nl)}$, $\hbar\hat{\Omega}_{\perp}^{(nl)}$ and $\hbar\hat{\Omega}_{\perp}^{(nn)}$ are defined in Eqs. (A20), (A26) and (A27); $\hbar\hat{\Omega}_{\text{in}}^{(n)}$ denotes the coupling between molecule $n$ with the incoming field.

A.4 Analytical and EM-Free Form of Hamiltonian #II

By following the procedure above, we can also obtain an EM-free form for Hamiltonian #II:

[TABLE]

where $\hbar\hat{\Omega}^{(nl)}$ and $\hbar\hat{\Omega}_{\perp}^{(nn)}$ are defined in Eqs. (A25) and (A27). The analytical and EM-free forms of Hamiltonians #I and #II allow us to perform simulations of coupled light-matter interactions with negligible computational cost; the propagation of the EM fields on a grid is no longer necessary.

A.5 Hamiltonians for a Pair of Two-Level Systems

In this manuscript, we have presented results for a minimalistic many-site model — a pair of identical TLSs (labeled as D and A). For convenience, we now report the analytical and EM-free form of Hamiltonian #I for the pair of TLSs.

Let us define $\hbar\Omega^{(n)}_{\text{in}}$ and $\hbar\Omega^{(nl)}_{\perp}$ as the norms of the corresponding operators that have already been defined, i.e., $\hbar\hat{\Omega}^{(n)}_{\text{in}}=\hbar\Omega^{(n)}_{\text{in}}\hat{\sigma}_{x}^{(n)}$, $\hbar\hat{\Omega}^{(nl)}_{\perp}=\hbar\Omega^{(nl)}_{\perp}\hat{\sigma}_{x}^{(n)}$ (for $n,l=D,A$), where $\hbar\hat{\Omega}^{(nl)}_{\perp}$ is defined in Eqs. (A26) and (A27). Then, for the pair of TLSs,

[TABLE]

Here, $v_{dd}=\frac{\eta_{3}^{(DA)}}{4\pi\epsilon_{0}R_{DA}^{3}}$, $V_{D}=\hbar\Omega^{(D)}_{\text{in}}+\hbar\Omega^{(DD)}_{\perp}+\hbar\Omega^{(DA)}_{\perp}$, $V_{A}=\hbar\Omega^{(A)}_{\text{in}}+\hbar\Omega^{(AA)}_{\perp}+\hbar\Omega^{(AD)}_{\perp}$.

Similarly, Hamiltonian #II for a pair of TLSs reads

[TABLE]

where $V_{D}^{\prime}=\hbar\Omega^{(D)}_{\text{in}}+\hbar\Omega^{(DD)}_{\perp}+\hbar\Omega^{(DA)}$, and $V_{A}^{\prime}=\hbar\Omega^{(A)}_{\text{in}}+\hbar\Omega^{(AA)}_{\perp}+\hbar\Omega^{(AD)}$. Here, as above, we have defined $\hbar\Omega^{(nl)}$ as the norm of $\hbar\hat{\Omega}^{(nl)}$ [defined in Eq. (A25)], i.e., $\hbar\hat{\Omega}^{(nl)}=\hbar\Omega^{(nl)}\hat{\sigma}_{x}^{(n)}$ (for $n\neq l$).

Finally, for the hybrid Hamiltonian for a pair of TLSs [see Eq. (13)], the Hamiltonian reads:

[TABLE]

Appendix B The Lehmberg-Agarwal master equation (LAME)

For $N$ identical TLSs, the Lehmberg-Agarwal master equationLehmberg (1970); Agarwal (1974) (LAME) is the standard theory to describe the reduced dynamics of the electronic degrees of freedom in an open quantum environment. Formally, one can derive the LAME by taking the Born-Markov approximation from QED and a rotating wave approximation (RWA), leading to

[TABLE]

where the dissipative term $\mathcal{L}_{\text{L}}[\hat{\rho}_{\text{N}}]$ is called the Lindbladian:

[TABLE]

Here, $\hat{\rho}_{N}$ denotes the $N$-body density operator, $\hbar\hat{\Omega}^{(n)}_{\text{in}}$ denotes the coupling between the incoming E-field and molecule $n$, and the $a_{nl}$ and $b_{nl}$ terms describe the collective damping and the collective level shifts, which are defined as

[TABLE]

where the dimensionless intermolecular separation $x_{nl}$ is defined as $x_{nl}\equiv\omega_{0}R_{nl}/c$, and the Kronecker delta function $\delta_{nl}$ equals to $1$ if $n=l$ and equals [math] otherwise. As might be guessed from the structures of $a_{nl}$ and $b_{nl}$ (that contain both $1/R^{6}$ and $1/R^{2}$ terms as well as $k_{\text{FGR}}$), LAME can accurately capture the time-resolved RET dynamics between a pair of TLSs at both short and long range.

Although not the focus of this paper, when modeling dynamics with LAME, one key problem is that for a system with $N$ TLSs, the method requires one to build an exponentially large many-body density matrix of size $2^{N}$ during the course of a simulation and update the Lindbladian at every time step. As a result, LAME is usually applied only to a few TLSs. Furthermore, applying LAME for inhomogeneous systems (not in vacuum) is not obvious; and more generally, like any master equation, LAME is accurate only in the limit of weak light-matter coupling.

B.1 Connecting Ehrenfest Dynamics with LAME

Above, in Fig. 3, we have observed that LAME agrees with Ehrenfest #II if the TLSs are weakly excited. Let us now analytically show that LAME can indeed be connected to Ehrenfest #II after some approximations are made. We will start from the EM-free form of Hamiltonian #II [see Eq. (A29)] and take the RWA form of Hamiltonian #II.

We assume that $\text{Im}\left[\rho_{ge}\right]\approx\widetilde{\rho}_{ge}\text{Im}\left[e^{i\omega_{0}t}\right]$, where $\widetilde{\rho}_{ge}\equiv\rho_{ge}e^{-i\omega_{0}t}$ is a slowly varying variable compared with the time scale of $\omega_{0}^{-1}$. With this assumption, we obtain from Eq. (A27) the RWA form of $\hbar\hat{\Omega}_{\perp}^{(nn)}$:

[TABLE]

Similarly, for $\hbar\hat{\Omega}^{(nl)}$ ($l\neq n$) in Eq. (A25c), the corresponding RWA form reads

[TABLE]

where c.c. denotes the complex conjugate. Let us make the following definitions

[TABLE]

so that all light-matter couplings can be rewritten in a uniform way (for both $l=n$ and $l\neq n$)

[TABLE]

Because $\gamma_{nl}$ in Eq. (B6) is a c-number, the real and imaginary parts of $\gamma_{nl}$ contribute differently to the electronic dynamics, i.e., the imaginary part of $\gamma_{nl}$ leads to dissipation while the real part should be a level shift. Thus, it is necessary to separate the real and imaginary parts of $\gamma_{nl}$:

[TABLE]

With these definitions, we can further rewrite Eq. (B7) as

[TABLE]

and the RWA form of Hamiltonian #II becomes

[TABLE]

The expressions for the $a_{nl}^{\prime}$ and $b_{nl}^{\prime}$ coefficients reported here [Eq. (B8)] are exactly the same as the coefficients of LAME [$a_{nl}$ and $b_{nl}$ in Eq. (B3)], provided that the step function $\theta\left(t-\frac{r_{nl}}{c}\right)$ in Eq. (B6) (or causality) is ignored. The exact agreement between these distinct coefficients clearly suggests a consistency between Hamiltonian #II and LAME. In fact, one can show that, if the step function (or causality) is ignored, $\hat{H}_{\text{II, RWA}}^{(n)}$ is exactly the effective mean-field Hamiltonian $\hat{H}_{\text{mf}}^{(n)}$ of MF-LAME defined by

[TABLE]

Eq. (B11) can be derived from Eq. (B1) by supposing $\hat{\rho}_{N}(t)=\hat{\rho}^{(1)}(t)\otimes\cdots\otimes\hat{\rho}^{(N)}(t)$ is valid at any time $t$ and then tracing out $N-1$ degrees of freedom to form the one-body reduced density operator; see Ref. 46 for a detailed procedure.

Appendix C The RET Rate

In Fig. 2, we have compared the RET rate calculated by different approaches. For the sake of completeness, we will now briefly review the RET rate theory. Furthermore, we will also analytically calculate the short-time result of $\rho_{ee}^{(A)}(t)$ as propagated with Ehrenfest dynamics by Hamiltonian #II, confirming the numerical calculations in Fig. 2 and the finite-difference time-domain (FDTD) simulation results in Ref. 21.

C.1 Perturbative QED Result

According to the standard perturbative QED calculationsSalam (2012, 2018), in the weak coupling limit, the RET rate between a pair of TLSs [donor (D) + acceptor (A)] can be calculated by Fermi’s golden rule:

[TABLE]

Here, $\rho_{f}$ denotes the density of states for the final state, $M(\omega,R_{DA})$ denotes the transition matrix element between the final state and initial state, and $R_{DA}$ denotes the separation between the donor and acceptor. In order to evaluate $M(\omega,R_{DA})$, let us use the notation $\lvert nmk\rangle$ to represent the donor in state $\lvert n\rangle$, the acceptor in state $\lvert m\rangle$, and the photon in state $\lvert k\rangle$. Whereas the initial state $\lvert eg0\rangle$ and the final state $\lvert ge0\rangle$ do not couple directly through the EM field, one can show that they are coupled at second order. To do so, one simply expands the initial state to first order, i.e., $\lvert eg0\rangle\rightarrow\lvert\psi_{eg0}\rangle=\lvert eg0\rangle+\sum\limits_{k}\frac{\lvert ggk\rangle\left\langle ggk|\hat{H}_{\text{int}}|eg0\right\rangle}{\omega_{k}-\omega_{0}}$, where $\hat{H}_{\text{int}}$ denotes the interaction Hamiltonian. Then, to second order in the interaction, the key contribution should come from the state $\lvert ggk\rangle$ where the photon has frequency $\omega_{k}=\omega_{0}$, i.e., the photon energy should equal the energy gap for both the donor and acceptor. To second-order in the perturbation, the coupling matrix $M(\omega,R_{DA})$ takes the form (after a few integrations in 3D)Salam (2010):

[TABLE]

At this point, consider the density of states for the acceptor ($\rho_{f}$ in Eq. (C1)). If there are no vibrations (or other electronic degrees of freedoms), then over the time scale $\omega_{0}^{-1}\ll t\ll k_{\text{FGR}}^{-1}$, [where $k_{\text{FGR}}$ denotes the spontaneous emission (SE) rate for a single TLS] there can be no true rate of energy transfer. Instead, one will find large oscillations back and forth. At very short times, the excited state population for the acceptor is simply

[TABLE]

where $\rho_{ee}^{D}(0)$ denotes the initial excited state population for the donor.

C.2 Analytical RET Rate by Hamiltonian #II

According to Ehrenfest dynamics with Hamiltonian #II, the equations of motion for the acceptor (A) read

[TABLE]

where $\boldsymbol{\mu}=\mu_{ge}\mathbf{e}_{d}$. At short times, because the excited state population for the acceptor is much smaller than the donor, the EM fields that are felt by the acceptor predominately come from the donor. Thus, at short times, we can neglect the donor’s population decay, [i.e., $\rho_{ge}^{(D)}(t)\approx\rho_{ge}^{(D)}(0)e^{i\omega_{0}t}$], so that the light-matter coupling term for the acceptor is just $-\boldsymbol{\mu}\cdot\mathbf{E}=\hbar\Omega\approx\hbar\Omega^{(AD)}$, where $\hbar\Omega^{(AD)}$ is defined by $\hbar\hat{\Omega}^{(AD)}=\hbar\Omega^{(AD)}\hat{\sigma}_{x}^{(A)}$ [see Eq. (B5)]. Furthermore, according to the RWA in Eq. (B5),

[TABLE]

Because there is no non-Hamiltonian dissipative term in Ehrenfest dynamics, purity is strictly conserved, i.e.,

[TABLE]

where $\sqrt{\rho_{ee}^{(A)}\rho_{gg}^{(A)}}$ is slowing varying compared with the time scale of $2\pi/\omega_{0}$, and $\varphi$ is the initial phase for the acceptor. By further substituting Eqs. (C5) and (C6) into Eq. (C4b), we obtain

[TABLE]

For short times, the acceptor is not strongly excited, i.e., $\sqrt{\rho_{ee}^{(A)}\rho_{gg}^{(A)}}\approx\sqrt{\rho_{ee}^{(A)}}$, so that Eq. (C7) is easily integrated [with $\int_{0}^{t}\theta(t^{\prime}-T)dt^{\prime}=(t-T)\theta(t-T)$]:

[TABLE]

In other words, Ehrenfest predicts that the excited state population on the acceptor will be just $\rho_{gg}^{(D)}(0)$ times the perturbative QED result [$\rho_{ee,\text{QED}}^{(A)}(t)$ in Eq. (C3)]. When the donor is near the ground state, i.e., $\rho_{gg}^{(D)}(0)\rightarrow 1$, Ehrenfest dynamics [with Hamiltonian #II] exactly recovers the perturbative QED result. Note that the analytical derivations here exactly agree with our previous FDTD simulationsLi et al. (2018b).