Resonant dipole–dipole interaction in confined and strong-coupling dielectric geometries

Using the electromagnetic response function of an electric dipole located within a dielectric geometry, we derive the mathematical equivalence between the classical response and quantum mechanical resonant dipole–dipole interaction between two quantum objects (atoms, quantum dots, etc). Cooperative spontaneous emission likewise emerges from this equivalence. We introduce a practical numerical technique using finite difference time domain for calculating both dipole–dipole interaction and collective spontaneous emission in confined dielectric structures, where strong light–matter coupling might arise. This method is capable of obtaining resonant dipole–dipole interaction over a wide range of frequencies in a single run. Our method recaptures the results of quantum mechanical second order perturbation theory for weak light–matter coupling. In strong coupling situations such as near a photonic band edge, second order Rayleigh–Schrödinger perturbation theory leads to divergences, and instead Brillouin–Wigner perturbation theory is required. This is equivalent to the use of a variational wavefunction to describe the exciton transfer between initial and final states. We introduce a system of coupled classical oscillators, that describes resonant dipole–dipole interaction and vacuum Rabi splitting in the strong-coupling regime, and that provides an effective numerical scheme based on the finite difference time domain method. This includes the effects of quantum entanglement and the correlation of quantum fluctuations. We discuss the crossover to Forster energy transfer when quantum correlations between the dipoles are damped by strong environmental interactions.


Introduction
Forster (or Fluorescence) resonant energy transfer (FRET) has played a vital role in biomedical science in the past decade [1][2][3] and it has become an essential tool for investigating biomolecular structures. For example, FRET is nowadays used to quantify protein-protein, protein-DNA interactions [1] and it has been also exploited as a means to study protein conformational dynamics [1,2] (see figures 1(a) and (b)). FRET is a process where two atoms, molecules, quantum dots or chromophores exchange energy via resonant dipole-dipole interaction (RDDI). RDDI is believed to play an important role in photosynthesis [4]. In this case, there has been interest in the question whether RDDI transfer of excitation occurs coherently or incoherently [5]. In the incoherent Forster energy transfer, the dipole oscillators interact strongly with surrounding vibrational modes of the host structure on time scales short compared to RDDI time scale. In this case, phase information of the oscillating dipoles is lost during the transfer. On the other hand, very strong light matter interaction can lead to coherent RDDI transfer. Coherent resonant energy transfer (CRET) is important for the entanglement of qubits in quantum information processing systems [6]. CRET can also have important consequences, especially in low temperature quantum many-body physics. For example, it has been suggested [7] that excitons confined to photonic band gap (PBG)quantum well heterostructures may acquire photon-like effective mass, leading to effects such as Bose-Einstein condensates at elevated temperatures [7].
Depending on the separation between the two quantum dipoles, RDDI can be classified as far field or near field. In the former case, the interaction is mediated by real photons, while in the latter scenario, virtual photons dominate the exciton exchange. Most RDDI based bio-experiments operate in the near-field regime. More specifically, one chromophore in its excited state (called donor) emits a real/virtual photon which is in turn absorbed by another chromophore in its ground state (acceptor). Even though energy conservation is satisfied between the start and the end of the process, the intermediate state does not have to conserve energy. An important consequence of the nature of FRET and the uncertainty relationship is that the interaction time is very small and hence the virtual photon has to be absorbed after a very short propagation distance. In free space, the interaction is isotropic. In the near field regime, the rate of transfer (norm squared of RDDI matrix element) falls as R −6 , where R is the distance between the two molecules. A typical dynamic range for FRET experiments in chromophores is in the range of a few to tens of nanometers as depicted in figure 1(c).
It is possible to alter the strength and the directivity of RDDI by careful engineering of the electromagnetic (EM) modes of the environment in which the interaction takes place. By removing a band of EM modes as in the form of a PBG, one can vary the strength of RDDI on scales longer than the resonant wavelength [8,9]. However, on shorter length scales (tens of nanometers) RDDI is dominated by virtual photons far outside of the PBG and the interaction strength reverts to its value in ordinary vacuum [8]. Also directional stop-gaps in the EM spectrum or changes in the degeneracy properties of certain modes, can affect the local photonic density of states, leading to engineering of RDDI [10,11]. Experiments have verified that optical confinement structures directly affect RDDI [12][13][14].
The efficient numerical evaluation of RDDI in complex confined geometries is a necessity for engineering both CRET and FRET responses at will. Current computation schemes rely on direct calculation using second order Rayleigh-Schrödinger (RS) perturbation theory. This involves the calculation of a principal value integral, and care must be taken to obtain the correct numerical result near the singularity. When this integral oscillates, convergence problems often arise. If the integral is analytically tractable, this obstacle can be overcome using convergence factors [15]. However, when the dielectric geometry under consideration is complicated, this approach is inapplicable and more advanced techniques are required to solve the convergence problem. One possibility is the Longman algorithm [16] for treating oscillating integrals. Prior to RDDI calculations, the EM modes and photonic dispersion relation of the system may be required. For example in three-dimensional (3D) PBG structures, calculating the dispersion curve along the Brillion zone edges does not suffice and one needs complete information about the 3D band structure at all wave-vectors before computing RDDI. At each point in the band structure, the EM field has to be calculated and stored at the locations of the two atoms. Then, the RDDI is obtained as a 3D integral over all modes. These factors make computing RDDI inside complicated geometries very time consuming.
A numerically adequate method to solve some of the above problems was introduced in [17]. It is based on classical oscillators and uses the finite difference time domain (FDTD) method to compute RDDI in free space and some dielectric structures. However, it is mathematically equivalent to using RS second order perturbation theory and encounters divergences in the strong-coupling regime.
In this work, we demonstrate a computational method that simplifies the evaluation of RDDI, using the numerical equivalence between RDDI and the classical dipole response function [18]. Unlike previous works [17], we introduce a new approach capable of treating geometries with strong resonant optical cavities and photonic band edges where second order perturbation theory fails. Using a variational wavefunction expansion, we derive the time evolution of the expansion coefficients. This non-perturbative quantum mechanical analysis is then mapped to the problem of two interacting classical oscillators coupled to classical radiation. This recaptures phenomena such as vacuum Rabi splitting [19], quantum entanglement between the emitters and leads to an effective numerical scheme for investigating coherent RDDI dynamics under extreme strong-coupling conditions. Our paper is organized as follows. We begin with a brief review of the quantum mechanical processes of RDDI and collective spontaneous emission (CSE). In section 3 we re-derive, using the Green function approach, the equivalence between classical dipole emission and coherent RDDI. This enables the use of conventional EM numerical methods such as the FDTD method to compute RDDI. This method can be used to compute RDDI over a wide range of atomic transition frequencies in one single run. In section 4 we compare FDTD results with those obtained using RS perturbation theory (analytically or numerically) for free space idealized geometries, confined ideal metallic cavities and photonic crystals. In section 5 we investigate strong-coupling RDDI dynamics using the variational formulation and map this to a system of two interacting oscillators embedded in colored vacuum. We demonstrate the accuracy of this technique by applying it to an analytically solvable example and we demonstrate its power for resonant strong-coupling systems where conventional methods fail. Finally in section 6, we discuss the difference in RDDI dynamics when the two dipoles constitute a quantum mechanical pure state (CRET) and when each dipole experiences strong interactions with the environment (FRET). In the former case, the correlations of the quantum fluctuations between the two dipoles lead to more rapid oscillatory dynamics, whereas the latter case is described using mean-field approximation to the complete quantum dynamics.

Weak coupling resonant dipole-dipole interaction (RDDI) and collective spontaneous emission (CSE)
In this section we briefly review the physics of RDDI in the weak coupling regime. To do so, we consider two identical two level atoms (or quantum dots) embedded in a certain dielectric structure. We neglect non-radiative losses and dephasing effects associated with lattice vibrations that destroy coherence in radiative transfer. Furthermore, we assume that one atom is initially in its excited state while the other is in the ground state. In other words, the two atoms taken together have only one exciton. This scenario is sufficient for obtaining the values of both RDDIs and CSE inside any photonic structure.

Resonant dipole-dipole interaction
In the absence of strong light-matter coupling, RDDI is described by second order RS quantum mechanical perturbation theory [15]. Resonant and non-resonant contributions to RDDI can be described using two Feynman diagrams [15]. The interaction between matter and radiation is described using the multipolar Hamiltonian in the dipole approximation Here, the index λ runs over all possible parameters that characterize the photonic modes. For instance, in free space λ ≡ {k, σ } where k is the mode wavevector and σ denotes polarization. For a photonic crystal environment, λ ≡ {k, n} runs over all Bloch wavevectors and band numbers n. For localized cavity modes, λ is a discrete index that includes all different photonic states characterized by their mode indices. In equation (1) µ a,b are the atomic dipole moment operators of atoms located at positions r a,b and the electric field operator associated with eigenfrequencies ω λ is given by In addition, ε 0 is the free space permittivity,h is the Planck's constant and V is the volume (or the mode volume in case of localized EM modes). Finally, a λ and a + λ are the annihilation and creation operators and E λ (r ) is the normalized electric field vector associated with the photonic mode under consideration [15], obtained by solving Maxwell's equations under the appropriate boundary conditions. Assuming one atom is initially excited while the other is in its ground state, a straightforward treatment of Hamiltonian (1) using the RS perturbation theory yields the RDDI expression (transition matrix element for exciton exchange between the two atoms) [15] M (ω 10 In equation (2) µ a and µ b are the vectorial dipole matrix elements of the first and the second atoms, respectively, and they are chosen to be real. The atomic transition energy of both atoms is given by E 10 =hω 10 (the subscripts 1 and 0 refer to the excited and ground states of a bare atom, respectively).

Collective spontaneous emission
Spontaneous emission is the process in which an excited atom interacts with the vacuum fluctuation to emit a photon and return to its ground state. In free space, this is described by first order perturbation theory leading to the famous Fermi's golden rule. The rate of spontaneous emission is not an intrinsic atomic quantity, but can be tailored by controlling the local EM density of state (LDOS) surrounding the atom [20]. Scientific interest in both single and multiatom spontaneous emission has been rekindled by the possibility of strong light-matter coupling in optical cavities and photonic crystals [21,22] with sharp and abrupt features in the LDOS.
Here light emission is no longer described as a simple Purcell effect [20] and novel features may arise [23].
In this section we consider the case of two atoms sharing one exciton. Unlike the single atom case, the two atoms may form an entangled state composed of a superposition of the first atom excited and the second in ground state and vice versa. RDDI is the mechanism responsible for entanglement and the RDDI coupling energy determines the energy shift due to the formation of even or odd atomic superposition states [15]. More specifically, a system composed of a pair of two-level atoms and one photon/exciton can exist in one of the three states ψ 0,λ = |g 1 , g 2 , 1 λ or |ψ ± = 1 √ 2 {|g 1 , e 2 , 0 ± |e 1 , g 2 , 0 }. Here |ψ 0,λ is the state with both atoms in ground states with one photon in the EM mode characterized by a set of parameters λ, while |ψ ± consists of two atoms sharing one exciton and zero photons. In CSE the entangled atoms in state |ψ ± emit a photon and relax to the state |ψ 0,λ . The dynamics of this process is described using a general state |ψ = b(t) |ψ ± + λ a λ (t)|ψ 0,λ in the Schrödinger equation ih d|ψ dt = H |ψ . Here, the Hamiltonian is composed of atomic, radiation and interaction parts: , H A |ψ 0,λ = 2hω g |ψ 0,λ , H R |ψ ± = 0, H R |ψ 0,λ =hω λ |ψ 0,λ and H Iλ is defined by (1) and the expansion of the electric field operator. In the first two relations the atomic operator H A acts only on the atomic part of the state, while in the latter two H R affects only the radiation part. Herehω e is the energy of a bare excited atom,hω g is the atomic ground state energy and ψ 0,λ |H Iλ |ψ ± ∝ δ λ λ . Projecting the Schrödinger equation on the states |ψ ± and |ψ 0,λ , yields (in the absence of any permanent atomic dipole moment): Eliminating the rapid time-dependence using the substitution b(t) = w(t)exp {i (ω e + ω g ) t} and a λ (t) = c λ (t)exp {i(2ω g + ω λ )t} and employing the Wigner-Weisskopf approximation [24], we obtain w(t) = w(0)exp − 1 2 t exp (i t). In this expression, = 2π is the spontaneous emission decay rate is the Lamb shift, with P V denoting the principal value. Using the interaction Hamiltonian (1), we find that V ±0 A similar expression is obtained for V 0± λ and we finally obtain the two-atom spontaneous emission rate: represent the spontaneous emission rates from the individual atoms a and b. The third term describes the modification in emission rate due to dipole-dipole interaction and is called the collective (or cooperative) spontaneous emission rate where ω 10 = ω e − ω g . We show below that coll is directly related to the imaginary part of the classical EM response function calculated from Maxwell's equations. In free space, coll is sizable only when the distance between the two atoms is much smaller than the wavelength associated with the atomic transition energy. In such situations, and in the case of µ a = µ b , direct substitution of E λ (r a ) ≈ E λ (r b ) in the formula of emission rate gives = 0 when the system is in the anti-symmetric state |ψ − . Under these conditions emission is possible only when the system is initially prepared in the symmetric state |ψ + .

Computational method
In this section we relate RDDI and CSE to the classical response function of the electric dipole. While this connection was demonstrated previously [18] using correlation functions and the fluctuation-dissipation theorem, our derivation based on EM emission from an electric dipole is more transparent.
Consider the classical EM wave equation in a periodic dielectric medium with polarization current source [25] Here c is the speed of light and ε(r ) is the dielectric function. The operator is the electric field vector. In equation (4). p ex (r, t) is the source electric polarization. For a single radiating atom located at r , this becomes p ex (r, t) = p ex (t)δ(r − r ). Introducing the temporal Fourier transform of F{ p ex (t)} = S(ω), it follows that The solution of equation (5) can be written in terms of the transverse and longitudinal eigenmodes as follows: Since our analysis here focuses on periodic media, the sum will be written explicitly over all allowed Bloch wavevectors k and photonic bands n. In other words, here λ ≡ {k, n}. In equation (6), the transverse and longitudinal auxiliary field vectors are the eigenfunctions of the equations H Q T kn (r ) = (ω 2 kn /c 2 ) Q T kn (r ) and H Q L kn (r ) = 0, respectively (superscript T means transverse and L stands for longitudinal) satisfying periodic boundary conditions over some large volume V . In general, they are both 3D vector fields defined over coordinate space and they satisfy the conditions [25]: (6) denotes a tensor product and, for any two vectors u = (6) can also be written as Using the completeness relations and Q(r ) = √ ε(r ) E(r ), we finally obtain Here δ is infinitesimally small number we have added to ensure that we retain only the retarded part of the Green function [25]. The real and imaginary parts of (8) can be separated using Here P V is the principal value and δ is the Dirac delta function.
Next we focus our attention on the case when r = r . For a polarization current of the electric dipole pointing in the direction of the unit vector e a , i.e. S(r a , ω) = S a (ω) e a = S a e a , we define the response function , where we have used r = r b and r = r a . Here e b is another unit vector. The response function χ ba , representing the component of electric field generated at point r b in the direction of e b due to a unit polarization field located at r a and oscillating at frequency ω along the e a vector, is then given by In the above equation, Here I 1 and I 3 correspond to the ω + ω kn terms and I 2 and I 4 to the ω − ω kn terms. Since ω kn = ω −kn and E kn = E * −kn and the summation is over −∞ < k x , k y , k z < ∞ it follows that I 1,2,3,4 are also real quantities By comparing equation (10) with equation (2) we finally arrive at where ω 10 is the atomic transition frequency. Equation (11) establishes the connection between quantum mechanical RDDI denoted by M and the response function of a classical electric dipole. This implies that any of the widely available EM numerical techniques can be used to compute RDDI in a complex dielectric geometry.
A similar analysis leads to the relation between CSE of the two-atom system and χ ba . Given the fact that I 3,4 are both real quantities, it follows that Here we have used the fact that δ (ω + ω kn ) = 0 since ω kn is always positive.
Comparing equations (3) and (12) we obtain the desired relation between the classical response and CSE: Relations (11) and (13) facilitate the use of conventional EM techniques such as the FDTD method to evaluate RDDI and CSE. In order to calculate RDDI and CSE rates, we replace the first atom with a classical dipole current and then use FDTD to calculate the resulting radiation at the position of the second atom. The Fourier transform of temporal electric fields are used in equation (9) to obtain the response function. RDDI and CSE follow from equations (11) and (13). In our calculations, we assume a time-localized Gaussian dipole current to cover a wide spectral range. More specifically, p ex (t) = exp − (t − t shift ) 2 /T 2 0 e ex where t shift controls the start of the polarization current with respect to the simulation window and T 0 is proportional to the pulse width. We use units where the speed of light in vacuum, free space electric permittivity and magnetic permeability are all set to unity. Material dispersion is neglected in order to simplify our FDTD algorithm. If necessary, more involved FDTD codes [26] can be used in order to take such effects into account.
We note finally that the above analysis is easily reduced to free space (periodic medium with ε(r ) = 1 throughout the unit cell) and localized modes in cavity structures. The latter case follows by noting that an optical cavity with localized modes can be treated as a unit cell in a periodic structure made of identical copies of the cavity with a period much larger than the localization length of the cavity mode.

Numerical results for weak light-matter coupling
In this section we compare RDDI and CSE values obtained using the FDTD, with semianalytical results obtained using quantum mechanical perturbation theory. In order to validate our numerical scheme; we consider the simple cases of free space, one-dimensional (1D) and two-dimensional (2D) cavities, and 1D PBG structures.
We note that all simulations were run on a Intel quad processor (3.07 GHz), 64 bit machine, which has 24 GB of RAM memory.

Free space
The 1D idealization corresponds to the interaction of two parallel planar quantum wells or sheets of quantum dots. For simplicity we assume that the dipole transition moments of the excitons are parallel and in the plane of the sheets. Then the radiation modes mediating RDDI are 1D plane waves propagating in the direction e z perpendicular to the sheets (only photonic modes with electric field vectors parallel to the dipole moments of excitons contribute to these calculations): E k (r ) = exp (i k z z). A direct substitution in the expression for the RDDI and converting the discrete summation over k z to an integral yields: In this idealized 1D geometry, the transverse area A of the EM quantization box is also the area of the quantum wells having dipole moments of µ a and µ b while R is the separation between the quantum well layers. In free space, ω k = c |k z | and we write the atomic transition frequency ω 10 = ck 10 . Using contour integration in the complex plane, the above integral reduces to M = µ a µ b 2Aε 0 k 10 sin(k 10 R). Similarly, for an where A (for the 1D case) and L (for the 2D case) are the area/length of quantum wells/wires having a dipole moments of µ a,b .
idealized 2D parallel quantum wires of length L, separated by a distance R, RDDI is given where Y 0 is the Neumann function of order zero. In 3D the RDDI between a pair of quantum dots separated by distance R is given by [15] HereR i ≡R ·x i , wherê R is the unit vector between the dots andx i is a unit vector along the axis i.
Note the singular behavior of RDDI transition matrix element M for R = 0 in both 2D and 3D geometries. Figures 2(a)-(c) provide a comparison of the RDDI obtained using FDTD and the corresponding exact analytical expression. Nearly, perfect matching is found in all cases.
Next we compare CSE calculated using FDTD with the corresponding analytical expressions in 1, 2 and 3 dimensions. In 1D, starting from equation (3), we have , with photonic modes given by propagating plane waves E k (r ) = exp (i k z z). Substituting in the above equation and using R = Aε 0h k 10 cos(k 10 R). In 2D, a similar analysis in cylindrical coordinates yields coll = µ a µ b Lε 0h k 2 10 2 J 0 (k 10 R) where J 0 is the zeroth order Bessel function. Finally, in 3D coll = µ a µ b k 3 10 2π ε 0h sin(k 10 R) where A (for 1D geometry) and L (for 2D geometry) are the area/length of quantum wells/wires having a dipole moments of µ a,b .
where the two dipoles are separated by a distance R along the x-axis and their dipole moments are aligned along the z-axis. In both figures 2 and 3, the accuracy of FDTD is frequency dependent. Agreement with analytical results is perfect at lower frequency and it deteriorates slowly as frequency increases. In EM simulations, higher resolution is necessary for higher frequency. Another observation is that loss of accuracy with frequency occurs more rapidly in higher dimensions. This can be attributed to numerical dispersion of FDTD (deviation of computed wave velocity from its real value due to discretization effects) which increases with dimensionality. In order to avoid this effect, a smaller Courant factor, S = (c/n)( t/ x 2 + y 2 + z 2 ), is required where c/n is the speed of light in the medium, t is the time step and x, y and z are the discretizations along the main axes. However, a smaller Courant factor is directly translated into more time steps for the same spatial resolution and consequently more simulation time.
In all simulations, the pulse parameters were chosen to be t shift = 33 × (R/c) and T 0 = 0.33 × (R/c). The resolution in both 1D and 2D simulations were R/60 and 40 000 time steps used. On the other hand, for 3D simulations, the resolution was R/30. The total simulation time for the 1D, 2D and 3D cases were 6 s, 4 min and 1.86 h, respectively.

RDDI with resonant cavity modes
In order to test the utility of our numerical method, we consider RDDI between identical parallel sheets of quantum dots placed within a perfect 1D mirror cavity with their dipole moment vector parallel to the cavity walls, figure 4(a). The cavity EM modes are given by E m (x) = √ 2/L c sin(mπ x/L c ) where L c is the cavity length and m is an integer index that denotes the mode order.
ω+ω m for sheets located at x and x . For x = x , the RDDI expression takes the form where A and L are the transverse (in direction perpendicular to the schematic cross section shown in figure 4) area/length of quantum wells/wires having a dipole moments of µ a,b and L c is the cavity side length.
The above quantity can be easily numerically evaluated. Figure 5(a) shows a comparison of M obtained from the above calculations with FDTD simulations. Good agreement is observed. In the FDTD computations, the parameters of the exciting current pulse were taken to be t shift = 10 × (L c /c) and T 0 = 0.2 × (L c /c). In addition, the spatial resolution was taken to be L c /200 and 120 000 time steps were used with an overall running time of ≈ 10 s.
In the 2D scenario of figure 4(b), we consider two columns of quantum dots placed in two different positions inside an idealized two dimensional mirror cavity. The EM modes inside the cavity take the form E m 1 ,m 2 (x, y) = 2/ L x L y sin(m 1 π x/L x ) sin(m 2 π y/L y ) and the RDDI coupling is given by Again we calculate the above expression directly (for the same parameters used in 1D case and L x = L y ≡ L c = 1) and compare the results with those obtained from FDTD. Figure 5(b) reveals good agreement at most frequencies. The FDTD simulations parameters are identical to those used for the 1D case and the simulation time is ≈ 35 s.
Closer inspection of figure 5(b), however shows that FDTD simulations predict a false RDDI resonance at a normalized frequency f ≡ ωL c /2π c ≈ 5.7 (as indicated by the green arrow in the figure). Increasing the simulation time or spatial resolution does not eliminate this effect. The numerical error occurs since the 2D square cavity has four degenerate modes at f ≈ 5.7. The correct RDDI results can be recovered by modifying the cavity dimensions to L x = L c and L y = L c + δ and evaluating RDDI in the limit of δ → 0 according to L'Hopital's rule.

RDDI in photonic crystals
In this section we consider RDDI within PBG geometries [8,9]. We first study RDDI in a 1D multilayer structure. This model considers only coupling due to modes propagating normal to the structure and neglects off-normal wave vectors. This simplification enables comparison between FDTD numerical results and the RDDI computed semi-analytically using the Floquet-Bloch theorem. We then demonstrate the versatility of our numerical scheme by computing RDDI between two quantum dipoles embedded in a 2D photonic crystal structure.
The 1D multilayer structure, with lattice constant a (see figure 6(a)) consists of two alternating layers: air with thickness 0.8a and a hypothetical non-absorbing material with 0.2a thickness and ε(r ) = 12. In this Kronig-Penny model [27], the EM modes are solutions of the Helmholtz equation in each homogeneous region with Bloch boundary conditions. In other words, the solution in each layer is written as E (l) kn (z) = u (l) kn (z)exp (i kz) where l = 1, 2 denote the solution in the air and dielectric layers, respectively. Using the periodic boundary condition u (1) kn (−0.8 a) = u (2) kn (0.2 a), we obtain the band structure and the field profile corresponding to each Bloch momentum k and band n at any point inside the crystal. The resulting band structure is depicted in figure 6(b). Using the dispersion relation and the field profiles in equation (2), we calculate the idealized 1D RDDI between dipoles located at the centers of two different silicon layers separated by a distance of 5a (yellow arrows in figure 6(a)). Here, only one polarization is considered and the summation over the Bloch momenta is replaced by k = L 2π ∞ −∞ dk, where L is the overall length of the 1D system. In principle, Bloch modes from all photonic bands contribute to RDDI. In practice, for atoms separated by five lattice constants, only the five bands surrounding the transition frequency are required for accurate results. The field amplitudes are normalized [25] according to D Q * k m (z)Q kn (z) dz = Lδ mn δ k k where Q kn (z) = √ ε (z)E kn (z) and both m and n denote band numbers. Using the periodicity of u kn (z), the domain of integration can be reduced to a single unit cell cell ε r (z)u * kn (z)u kn (z) dz = a. A comparison between the fully numerical FDTD scheme and the semi-analytical Kronig-Penny analysis is shown in figure 6(c). Similar to the 1D free space case, here A is the overall transverse area of the quantum wells having dipole moments of µ a and µ b . Aaε 0 where A is the area of quantum wells (in a plane perpendicular to the 1D periodic structure) having dipole moments of µ a,b and a is the periodicity of the one dimensional PBG structure.
The FDTD results are obtained using a multilayer structure composed of 4001 unit cells excited with a Gaussian polarization current of T 0 = 0.4 × (a/c) pulse width and the simulation time was taken to be 5000 in units of a/c. In our computation, we used a large number of layers and the simulations were terminated before the radiated electric field can bounce back from the edge of the photonic crystal, ensuring that all Fabry-Perot resonances are eliminated from the calculations. The comparison in figure 6(c) shows a good agreement between FDTD calculations and those obtained using the Kronig-Penny model, except close to the band edge where RDDI becomes singular. A proper treatment of this situation requires the nonperturbative, strong-coupling analysis of section 5. We note that the RDDI value at the lower edge of the second band does not show a singular behavior in spite of the zero group velocity at that point. This is made possible since the local density of states at this band edge vanishes at the center of the dielectric layer.
Finally we study dipole-dipole interactions of two quantum radiators embedded inside the idealized 2D square lattice photonic crystal shown in figure 7(a). Here L represents the length of EM quantization box along the dielectric rods which is also the length of quantum wires (in a direction along the rods of the 2D photonic crystal) having a dipole moments of µ a,b .
As before, the cylinders are made of a high dielectric, non-absorbing material of ε(r ) = 12 in an air background, and are arranged on a square lattice. Each cylinder has a radius r = 0.2a where a is the lattice constant. In our idealized system, only wavevectors perpendicular to the cylinders are considered and the electric field is assumed parallel to the cylinders. The two quantum dipoles are located at the centers of the two cylinders marked by different colors/grayness. The band structure and RDDI values are depicted in figures 7(b) and (c), respectively. The resonant behavior at the upper edge of the first band is a direct result of vanishing group velocity and once again requires the more general treatment of section 5.
Very recently, another FDTD method for calculating RDDI was introduced [28]. While our Green's function method requires only one FDTD run in order to obtain RDDI and CSE values over wide range of frequencies, the authors of [28] require three different FDTD runs for every frequency point. Moreover computations in [28] yield CSE and a final Hilbert transform is needed to obtain RDDI. Both techniques compute RDDI corresponding to second order perturbation theory [15]. However, in situations where the eigenmodes of the system are discrete or when the principal value limit fails to exist, this perturbation approximation completely fails and gives infinities. A non-perturbative approach is required in these strong-coupling situations.

Strong light-matter coupling
RDDI in optical cavities with discrete modes or near certain band edges of photonic crystal structures is associated with strong light-matter coupling. In these situations, simple perturbation theory leads to divergent behavior. These singularities are also apparent in the FDTD method of section 4. Previous calculations were based on treating RDDI as a perturbation over the bare atoms. Moreover, these interactions were calculated using the RS perturbation theory [15]. Near an isolated discrete mode or close to dispersion band edges, the relevant atomic levels may be strongly dressed by the EM modes. Higher order approximations are required to describe these strong interactions. For example, a variational wavefunction approach was previously shown to recapture the phenomenon of vacuum Rabi splitting near a PBG [29]. This method is equivalent to the use of Brillouin-Wigner perturbation theory [30]. Recently, it was shown that the phenomenon of atomic Mollow splitting in a strong external field is recaptured using a self-consistent solution of the coupled Maxwell-Bloch (MB) equations [31].
In this section we extend our FDTD method to treat such strong coupling situations. Our generalized numerical method captures effects such as dressed atomic states and vacuum Rabi splitting, thereby removing the unphysical divergences of standard perturbation theory.

Non-perturbative analysis
Our analysis is based on the variational wavefunction method [29,32]. Consider the Schrödinger equation ih ∂|ψ ∂t = H |ψ , where the Hamiltonian H = H 0 + V consists of atomic and radiation terms H 0 (when the interaction is switched off) and the atom-photon interaction term V . We denote the states of interest as |a = |e, g, 0 , |b = |g, e, 0 and |k eλ = |e, e, 1 λ and |k gλ = |g, g, 1 λ . The first two represent states where one atom is excited with no photon while the last two describe states with one photon in the EM mode defined by the set of parameters λ while both atoms are excited or in their ground state. We consider a variational wavefunction consisting of a superposition of these states with corresponding time-dependent variational amplitudes Substituting (14) back in the Schrödinger equation and projecting the result on each of the above basis states yields four coupled differential equations for the variational coefficients: In the above equations, V α,β = α|V |k β = (V β,α ) * with α = {a, b} and β = {gλ, eλ}. Here E 10 is the energy of the excited atomic state with respect to atomic ground state energy which is taken to be zero. Note that our analysis does not invoke the rotating wave approximation that can lead to incorrect answers in certain situations [32]. Taking the Laplace transformation and denoting the Laplace transform of the coefficients by capital letters, we obtain the set of algebraic equations: ihs K eλ (s) = (2E 10 +hω λ ) K eλ (s) Here A(s), B(s), K gλ (s) and K eλ (s) are the Laplace transforms of a(t), b(t), k gλ (t) and k eλ (t), respectively. In what follows, we assume the system is initially prepared in state |a so that a(0) = 1.
In the above equations, each photonic mode is coupled to both atoms whereas there is no direct coupling between the atoms. Since the atom-atom interactions are mediated only by photons, we can eliminate the two middle equations and obtain a system of two coupled equations between the atoms: In (17), we define the complex interaction energy parameters where subscripts m and n denote either atom a or b and m = n describes self-interaction. The expression (18) is analogous to the perturbative RDDI expression (2) provided that we set ihs = E 10 . However, near optical resonances and photonic band edges, this so called pole approximation leads to unphysical singularities in RDDI. Our variational procedure circumvents these divergences by enabling effects such as vacuum Rabi splitting [19] and provides the correct physical behavior in such strong-coupling situations. In our more precise derivation of (18) mn can have both real and imaginary parts. Physically, the imaginary parts of mn (s) (m = n) describe the modification of the spontaneous emission rate of a single atom due to its interaction with the second atom. If the atoms are widely separated from each other (such that the rate of radiative decay becomes comparable or large compared to the rate of the resonant energy transfer), the coherent RDDI can be strongly suppressed relative to the value predicted by simple perturbation theory. This suppression of RDDI energy transfer may occur even for small separations if the bare atomic transition occurs in a spectral range where the EM density of states is very large. Similar damping effects are expected if the atoms are strongly coupled to phonons of the host material in which they are embedded. However, as we show below, in some idealized 1D systems the imaginary part of m,n (s) can lead to energy transfer even when the RDDI coupling obtained using RS perturbation theory of section 4 is zero (see figure 8(b) and related discussion).
Equations (17) can be written in the compact matrix form Similar to the derivation of equation (10), it is easy to verify that ab = ba , which we denote by a single complex number m . Likewise, we introduce the simplified notation for the selfenergies a = aa and b = bb . It follows that Since our system was initially prepared in state |a , the evolution of the system toward state |b is fully characterized by the coefficient b(t) obtained from the inverse Laplace transform of B. Thus knowledge of B suffices for obtaining RDDI dynamics. From (20), we have This non-perturbative result is valid in the absence of higher order processes such as multiphoton generation. We note that equation (21) can be also derived using an equivalent resolvant operator [9,33]. Our trial wavefunction (equation (14)) spans the one excitation sector and the three excitation sector. The two excitation sector (with one atom excited in the presence of one photon) is not included since states in this sector do not couple directly to the initial states of interest through the dipole interaction Hamiltonian considered here.
In the (perturbative) pole approximation (ihs = E 10 ), the parameters a , b and m become divergent near a photonic band edge (or in the presence of discrete modes where the summations cannot be transformed into integrations). However in our non-perturbative variational calculations (s remains a general complex number), these singularities are regulated by the dressing of the atomic levels by photons and the resulting repulsion of the atomic level from the cavity mode (vacuum Rabi splitting). In this case, one must calculate the exact inverse Laplace transform without any approximations or solve the time domain equations directly.

Equivalent classical oscillator model
We now introduce a classical EM model, mathematically equivalent to formula (21), that enables computation of dipole-dipole response using standard numerical techniques for solving Maxwell's equations. Consider a system of two coupled oscillators: where P a,b are the respective complex dipole moments of the oscillators, ω 10 is their natural oscillation frequency (chosen to match the atomic bare transition frequency ω 10 = E 10 /h) and ρ is a constant to be determined later. Also r a,b are the positions of the first and second oscillator, respectively. The total radiated electric field is given by E rad (t, r i ) = E ii + E i j , where E i j is the field produced at point i due to a dipole located at point j with i, j = {a, b}. In other words, the radiation field at each oscillator is the sum of the field radiated by the oscillator itself plus that produced by the other oscillator. The divergent static part of the field is not included here. The field amplitude E rad (t, r ) is required to satisfy Maxwell's equations in the confined strongcoupling geometry of interest. Since we use complex representations of the polarizations and electric fields, only first order time derivatives appear in equation (22). The imaginary parts of the fields can be eliminated to obtain second order differential equations for the physical field values. However equations (22) are more simple and computationally convenient. Assuming only one dipole is initially excited and equal dipole transition matrix elements, i.e. µ a = µ b ≡ µ, and taking the Laplace transformation of (22) yield where we have introduced the interaction energy parameter m (s) = µ i µ j χ i j (s) for i = j and i (s) = µ 2 i χ ii (s) in terms of the susceptibilities χ i j (s) = E i j (s)/P j (s). By choosing ρ = µ 2 , we obtain P b (s) P a (0) = ih m (s) (ihs −hω 10 − a (s)) (ihs −hω 10 The above formula is mathematically equivalent to equation (21) when the two interacting atoms experience the same photonic local density of states. It follows that, under this condition, a numerical solution of the coupled oscillator model (22) (together with the field equations for E rad (t, r )) in time domain provides a full description of the RDDI dynamics. Despite the classical nature of this equivalent oscillator model, it recaptures the effects of quantum entanglement of the two atoms as described in the state vector (14). In contrast, the semiclassical MB description [31] neglects the correlations of quantum fluctuations between the dipoles in favor of more dominant interactions with the environment. In order to numerically simulate RDDI with equivalent classical oscillators, the initial value of the dipole moment of the first oscillator is set to unity (corresponding to the excited atom) while that of the second oscillator is set to zero (to represent the ground state atom). As the system evolves with time according to equations (22) in conjunction with Maxwell's equations within the nanostructured dielectric medium, the first dipole will oscillate and emit EM fields. The field radiation is calculated using an FDTD algorithm with that particular dipole as a source term. The radiated field may scatter from the dielectric environment and then influence the oscillation of the initial dipole. This captures the self-interaction (atomic dressed state) effects described by aa (s) ≡ a (s) in equations (18)- (21). The radiated field (including self-dressing effects), upon arriving at the position of the second dipole acts as a driving force that induces oscillations. The resulting radiation from the second dipole is added to the total EM field through Maxwell's equation. The process continues and the fields due to both dipoles interact with the dielectric environment, cause self-dressing effects, and contribute to the excitation/deexcitation of each dipole. From the equivalence between equations (24) and (21), it is evident that the quantity |P a (t)/P a (0)| 2 corresponds to the quantum mechanical probability that an initially excited atom remains excited. Likewise, |P b (t)/P a (0)| 2 gives the time domain quantum dynamics of the energy transfer between the two atoms.
The numerical integration of equation (22) is performed using the Crank-Nicolson method [34]. This scheme has been shown to yield stable results in the context of MB equations [31] and the same applies here. A similar concept of coupled oscillators has been employed [35] to calculate the non-retarded dispersion force between atoms [36]. Given recent advances in fabricating 3D photonic crystals enabling strong light-matter coupling at visible frequencies [37], non-perturbative techniques may be essential in calculating RDDI inside such structures.

Illustrations and numerical results
We verify the efficiency of our two-oscillator model in computing dynamical evolution of atomic systems by first applying it in environments where the system's behavior is already known. We consider RDDI between two quantum dot sheets located in free space. The geometry is one dimensional and analytical solutions for this problem can be readily obtained. As in [31], we consider a planar layer of identical quantum dots, each of which is cubic with side length of L D = 4.652 nm and having an electric dipole moment of µ = 1.342 × 10 −28 Cm. The transition wavelength of each quantum dot is λ 10 = 1.5 µm and the dots are tightly assembled with a packing density of N A = 1/L 2 D . In section 4.1.1, we showed that the RDDI between two such quantum dot layers is given by m R ≡ M = µ 2 2A d ε 0 k 10 sin(k 10 R) where k 10 = ω 10 /c and R is the distance between the two layers. Also, A d = L 2 d is the average area of each quantum dot. The quantity M represents only the real part of the interaction coupling. In order to fully characterize the time domain dynamics of the RDDI process, the full complex coupling constant [8,9] must be considered. In the present geometry, the imaginary part is simply For m I = 0, the term inside the brackets reduces to sin( m R t/h) whereas for vanishing m R , the energy transfer between the two objects is given by |b(t)| 2 = exp (−2| s I | t/h)sinh( m I t/h) 2 . Since | s I | | m I |, the energy transfer probability given by |b(t)| 2 never exceeds unity.
Similarly, the probability amplitude of layer one to remain excited after initial excitation is given by We consider two specific situations, namely k 10 R = (n + 1/2)π and nπ with n being an integer. In the first, we have m I = 0 and in the latter m R = 0. For atomic transition wavelength λ 10 = 1.5 µm, the first scenario is achieved when R = λ 10 /4 = 0.375 µm. Under these conditions, m R /h = − s I /h ≈ 3.739 × 10 12 s −1 and m I = 0. Also note that in these 1D geometries, the Lamb shift s R = 0. The excitation probabilities then reduce to |b(t)| 2 = exp(−2| s I |t/h)sin 2 ( m R t/h) and |a(t)| 2 = exp(−2| s I |t/h)cos 2 ( m R t/h). Figure 8(a) depicts the close agreement between numerical and analytical results obtained for |a(t)| 2 (blue line starting at 1) and |b(t)| 2 (green line starting at zero) when R ≈ 0.375 µm. Due to spontaneous emission probability from each dot layer, the system ends up in the ground state with one photon emitted.
In the second scenario, when R = λ 10 /2 ≈ 0.75 µm, the excitation probabilities evolve as |b(t)| 2 = exp(−2| s I |t/h)sinh 2 ( m I t/h) and |a(t)| 2 = exp(−2| s I | t/h)cosh 2 ( m I t/h). Since s I = − m I neither expression exceeds unity. The dynamical evolution of the system in this case is illustrated in figure 8(b) where again good agreement between analytical and numerical results is evident.
An amusing characteristic of this choice of parameters is that the system initially decays and asymptotically approaches a probability of 1/4 that each atom remains excited. This feature is a peculiarity of 1D geometries where it is possible to cancel spontaneous emission from both atoms by destructive interference.
To further illustrate our technique, we revisit the problem of RDDI between two layers of quantum dots located inside a 1D metallic cavity (see figure 9(a)). We use the same parameters as before and we consider two different cavity lengths: L c = 1.5 µm and L c = 3 µm. In both cases, the layers are assumed to be positioned at L c /4 and 3L c /4. As before, the transition wavelength is λ 10 = 1.5 µm. For L c = 1.5 µm, the atomic transition frequency matches the second order cavity mode and the two layers are located at its maximum amplitude. Since both atoms couple strongly to the same mode, there is a significant dipole-dipole interaction. This is clearly observed in figure 9(b) where large and rapid amplitude oscillations occur.
On the other hand, when the cavity length is L c = 3 µm, the radiated fields match the fourth order cavity mode. Since both quantum dot layers are located at the nodes of this resonant mode, very weak interaction is expected. Figure 9(c) confirms this expectation. The first excited quantum dot layer remains mostly in its initial state while the second remains largely at its We note that for the resonant case of L c = 1.5 µm, high spatial resolution of L c /500 and a Courant factor of S = c dt/dx = 1/2 were required for the FDTD calculations to yield accurate non-divergent results.
We next test our oscillator model for single atom emission in the 1D photonic crystal environment studied in [31]. This structure is similar to that depicted in figure 6(a) except that here the high dielectric layers are made of SiO 2 having ε(r ) = 2.25 and it exhibits a PBG between 0.4080 < ω a/2π c < 0.4959. We consider a single atom located at the center of a high dielectric layer with ε(r ) = 2.25. The limiting case of single atom is numerically modeled by placing a classical oscillator at the location of that atom while removing the second oscillators to infinity. Vacuum Rabi splitting was predicted [31] for an atomic transition frequency of ω 10 a/2π c = 0.403. Figure 10(a) compares the transition probability computed for this system using our method (oscillatory decaying blue curve) and the same quantity in free space (purely decaying red curve). Figure 10(b) shows the Fourier transform of the radiated electric field at the output (calculated near the absorbing boundary layer). Clearly, vacuum Rabi splitting, similar to that predicted in [31], is observed. Since the classical oscillator model used here and the optical Bloch equations used in [31] are not identical, the exact ratio of the peak heights is not the same.
Finally we use our equivalent classical oscillator technique to compute the RDDI dynamics within the 1D photonic crystal shown in figure 6(a). The two atoms are located at the centers of the high dielectric layers five unit cells apart. As observed before, direct second order perturbation theory (as well as the previously introduced FDTD algorithm using the one oscillator model with no back action) predicts singular behavior for the coupling strength.
Here, the high dielectric layers have ε(r ) = 12 and the crystal exhibits a PBG in the range 0.2125 < ω a/2πc < 0.461. By choosing the lattice constant a = 318.75 nm, the lower band gap edge ω a/2π c = 0.2125 corresponds to a free space wavelength of 1.5 µm. Figure 11(a) Figure 11. RDDI between the two quantum layers indicated by yellow arrows in figure 6(a) when the atomic transition frequency (corresponding to λ 10 ≈ 3 µm) lies in the middle of the first band (a) and exactly (λ 10 = 1.5 µm) at the first band edge (b). In (a), the near linear dispersion leads to a behavior resembling free space. On the other hand, in (b) strong light-matter interaction is described without the unphysical singularities of the second order perturbation theory. Quantum dots parameters are presented in the text. Horizontal axes in both figures represent time t in seconds.
shows the system evolution when the atomic transition frequency is in the middle of the first band. Clearly, the evolution resembles that of free space, since in this range the dispersion is almost linear. On the other hand, figure 11(b) depicts the interaction dynamics when the atomic transition frequency coincides with the lower edge of the first band gap.
While second order perturbation theory predicts singular behavior, our time domain analysis predicts a well-defined strong interaction manifested by the highly oscillatory probabilities. Clearly our strong-coupling analysis, using an equivalent pair of classical oscillators, circumvents unphysical singularities of the previous perturbative analysis.

Phonon dephasing and non-radiative damping
Our analysis in the previous sections describes the interaction of two atoms sharing one exciton in the absence of any influence from the environment except from radiative decay. It does not account for depletion of population inversion due to dephasing effects or non-radiative decay arising from phonons in the host material. In general, the dynamics of energy transfer will depend on the interplay between RDDI and environmental effects. If the RDDIs between the two atoms are much stronger than the dephasing and damping effects, then the energy transfer process is coherent and is referred to as CRET. This is accurately modeled by the equivalent twooscillator model of section 5.2. On the other hand, if interaction with the environment dominates over RDDI, decoherence occurs and entanglement of the quantum state of the two atoms is lost. This is described using MB equations [31]. In this case the exact quantum correlation between the two atomic dipoles is replaced by mean-field approximation in which one dipole exerts a semi-classical, effective, electric-field on the other dipole.
We present below a brief description for two interacting atoms sharing one exciton in the Heisenberg picture. We derive the operator form of the MB equations and discuss a mean-field factorization in which the system is reduced to the semiclassical MB description. This modifies the CRET process to the more common process of FRET (where dephasing and damping effects are significant).
The full quantum mechanical Hamiltonian describing the interaction of multiple two-level atoms with a quantized EM field can be written as [38] H = λh ω λ a + λ a λ +h Here, σ z j is the atomic inversion operator while σ + j and σ j represent the atomic excitation and de-excitation operators (for the jth atom located at r j ), respectively. They satisfy the commutation relations: (27), a + λ and a λ are the creation and annihilation operators for photons in mode λ; ω 10 and ω λ are the atomic transition frequency and the photonic mode frequency, respectively and µ j is the vector transition dipole matrix element of atom j. Finally { E λ (r )} are a complete set of EM eigenmodes of the photonic structure (satisfying periodic boundary conditions over volume V ) and the coupling constant is given by The full Hamiltonian (27) includes (counter-rotating) terms usually ignored in the rotating wave approximation. Introducing the Pauli operators σ x j = σ j + σ + j and σ y j = i(σ j − σ + j ), equation (27) can be written as Here the electric field operatorÊ(r j ) = i λ h g a λ E * λ (r j ) − a + λ E λ (r j ) . In deriving (28) we assume that Heisenberg operators of different atoms commute and they in turn commute with the field operators [36]. Using the Heisenberg equation of motion (ih dÂ dt = [Â, H ]), we obtain The summation over j in (29d) provides atom-atom interaction as mediated by EM modes of the structure. Provided that operator products of the formÊ(r j ) σ y,z j are treated exactly, the quantum expectation values of equation (29) include the effects of entanglement between the atoms. However, if the system is interacting with a thermal bath of phonons, decoherence effects lead to a de-correlation of fluctuations between atomic and field operators. This motivates the mean field factorization Ê (r j ) σ y,z j = Ê (r j ) σ y,z j and reduces equations (29a)-(29c) to d dt where we identify the quantum expectation value of the field operators with classical EM fields satisfying Maxwell's equations In equation (31) we use the notation Ê (r, t) = E(r, t), and we identify H (r, t) and J (r, t) as the classical magnetic field and classical polarization current respectively. In this approximation, quantum entanglement between atoms is no longer accounted for. This is a departure from the equivalent classical oscillator model presented in section 5.2, where such entanglement effects play an important role in RDDI dynamics.
Dephasing effects due to interaction with phonons and non-radiative damping arising from anharmonic phonon decay to lower frequency phonons [39] can be phenomenologically incorporated in the optical Bloch equations [31] d dt Here T 2 and T 1 are the dephasing and non-radiative damping lifetimes, respectively. Numerical investigation of equations (32) together with (31) provides detailed dynamics of the system of interacting atoms under these conditions (decoherence, dephasing and non-radiative damping), and thus of great importance for FRET imaging experiments.
In order to gain insight into the energy transfer dynamics as described by equations (31) and (32), we consider the case of two 1D sheets of quantum dots inside a 1D perfect mirror cavity. This example was investigated in section 5.3 using the equivalent classical oscillator model (see figure 9) in the case of quantum mechanical pure state evolution. Here we analyze the same problem under decoherence conditions. Figure 12(a) shows the population inversions of the two layers as a function of time under the conditions T 1 = T 2 = ∞. We assume that the first layer of quantum dot is partially excited with σ y 1 = 0, σ z 1 = 0.9 and σ x 1 = √ 0.19, while the second atom is in its ground state with σ z 2 = −1, σ x 1 = 0 and σ y 1 = 0. Note that the energy transfer rate obtained in the semi-classical model is less than its pure state counterpart. More rapid oscillations occur when atomic fluctuations are quantum-mechanically correlated.
Finally, we study the dynamics when T 1 = T 2 = 10 −12 s. Under these conditions, energy transfer is accompanied by loss of population inversion as expected ( figure 12(b)). The character of the evolution depends on the competition between the dephasing/damping rate and the rate of energy transfer. If the energy transfer rate is faster than both dephasing and non-radiative damping rates, the two atoms will exchange energy several times before decaying to their ground states. On the other hand if decoherence effects occur on a time scale much less than energy transfer rate, the exchanges are less frequent or the first excited atom will decay to its ground state without any appreciable resonant energy transfer.
It would be of interest for future studies to develop models and numerical schemes capable of treating the gradual transition between CRET and FRET. A possible route towards this goal is to incorporate dephasing and/or damping effects in the quantum system from the beginning and to re-derive a new oscillator model that can capture their effects on the system dynamics. For instance one might try to model the effects of the environment using Brownian oscillators [40]. Such a study will not only be important for FRET imaging techniques, but will also shed light on the dynamics of decoherence.

Conclusion
We have presented a simple derivation of the mathematical equivalence between weak-coupling RDDI (and CSE) and the Green functions of Maxwell's equation. We have further demonstrated a powerful computational algorithm that provides an accurate description of coherent RDDI in complicated strong-coupling photonic structures using FDTD calculation. Our method supersedes simple perturbation theory in which singularities may appear at the edge of band gaps of photonic crystals or in the case of discrete cavity modes. Our analysis is based on the variational wavefunction approach (equivalent to resolvant operator method) applied to the onephoton, two-atom sector of the quantum electrodynamics Hilbert space.
Our numerical calculation utilizes an equivalent model of 'coupled classical oscillators' in which two oscillators are coupled to each other through their radiation fields and includes back reaction from the EM environment. We demonstrated that solving this model (in conjunction with Maxwell's equation) is equivalent to the results obtained from the variational wavefunction analysis for a quantum-mechanical pure state with entanglement between atoms. However, the two oscillators model has the advantage of being amenable to time domain analysis, providing an easy route to RDDI and CSE dynamics. Our numerical algorithm recaptures strong-coupling effects such as vacuum Rabi splitting inside photonic crystals. This enables accurate calculation of RDDI between two quantum objects under resonant, strong-coupling conditions, where the atomic transition frequency matches a cavity resonant frequency or a photonic band edge.
Finally we have compared our pure state quantum dynamics to what happens when decoherence is enabled by environmental interactions. This regime is described by imposing mean field approximations on the quantum system. This eliminates correlations in quantum fluctuations between the two atoms and fields. Under these conditions, The Heisenberg equations of motions for the atom-field system reduce to MB equations. We presented numerical results illustrating the slower energy exchange dynamics in FRET compared to CRET. The gradual transition between the fully coherent interactions and the semi-classical analysis is of great interest on both practical and theoretical fronts and merits further investigation.