A master equation for strongly interacting dipoles

Here we identify sufficient conditions in order that the same master equation can be obtained from different microscopic Hamiltonians. This will be important when it comes to deriving the two-dipole master equation in the following sections. Let us consider a single dipole within the electromagnetic bath, and assume that there are only two relevant states ($| g\rangle ,| e\rangle$) of the dipole separated by energy ω₀ = ω_e − ω_g. Associated raising and lowering operators are defined by ${\sigma }^{+}=| e\rangle \langle g|$ and ${\sigma }^{-}=| g\rangle \langle e|$. The electromagnetic bath is described by creation and annihilation operators ${a}_{{\bf{k}}\lambda }^{\dagger },\,{a}_{{\bf{k}}\lambda }$ for a single photon with momentum k and polarisation λ. The photon frequency is denoted ${\omega }_{k}=| {\bf{k}}|$.

The energy of the dipole-field system is given by a Hamiltonian of the form H = H₀ + V, where

$$\begin{eqnarray}&&{H}_{0}={\omega }_{0}{\sigma }^{+}{\sigma }^{-}+\displaystyle \sum _{{\bf{k}}\lambda }{\omega }_{k}{a}_{{\bf{k}}\lambda }^{\dagger }{a}_{{\bf{k}}\lambda },\end{eqnarray} \tag{ 1 }$$

defines the free (unperturbed) Hamiltonian and V denotes the interaction Hamiltonian. Gauge-freedom within the microscopic description results in the freedom to choose a number of possible interaction Hamiltonians. We define the generalised-gauge transformation [26]

$$\begin{eqnarray}&&{R}_{\{{\alpha }_{k}\}}:= \exp [\hat{{\bf{d}}}\cdot {{\bf{A}}}_{\{{\alpha }_{k}\}}(0)],\qquad {{\bf{A}}}_{\{{\alpha }_{k}\}}({\bf{x}})=\displaystyle \sum _{{\bf{k}}\lambda }{\left(\displaystyle \frac{1}{2{\omega }_{k}{L}^{3}}\right)}^{\displaystyle \frac{1}{2}}{\alpha }_{k}{{\bf{e}}}_{{\bf{k}}\lambda }{a}_{{\bf{k}}\lambda }{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{x}}}+{\rm{H}}.{\rm{c}}.\end{eqnarray} \tag{ 2 }$$

In this expression the ${\alpha }_{k}$ are real and dimensionless, the e_kλ, λ = 1, 2 are mutually orthogonal polarisation unit vectors, which are both orthogonal to k, L³ is the volume of the assumed fictitious quantisation cavity, and $\hat{{\bf{d}}}$ denotes the dipole moment operator. By making the two-level approximation after having transformed the Coulomb gauge Hamiltonian using the unitary operator ${R}_{\{{\alpha }_{k}\}}$ we obtain the Hamiltonian H = H₀ + V where [26]

$$\begin{eqnarray}&&V=\left[\displaystyle \sum _{{\bf{k}}\lambda }{g}_{{\bf{k}}\lambda }\,{\sigma }^{+}({u}_{k}^{+}\,{a}_{{\bf{k}}\lambda }^{\dagger }+{u}_{k}^{-}\,{a}_{{\bf{k}}\lambda })+{\rm{h}}.{\rm{c}}.\right]+{V}^{(2)}\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&{V}^{(2)}=\displaystyle \sum _{{\bf{k}}\lambda }\displaystyle \frac{1}{2{L}^{3}}\,{\alpha }_{k}^{2}\,| {{\bf{e}}}_{{\bf{k}}\lambda }\cdot {\bf{d}}{| }^{2}+\displaystyle \frac{{e}^{2}}{2m}\,\tilde{{\bf{A}}}{(0)}^{2}.\end{eqnarray} \tag{ 4 }$$

The term ${V}^{(2)}$ is a self-energy term, which does not act within the two-level dipole Hilbert space and which depends on the field

$$\begin{eqnarray}&&\tilde{{\bf{A}}}({\bf{x}})=\displaystyle \sum _{{\bf{k}}\lambda }{\left(\displaystyle \frac{1}{2\omega {L}^{3}}\right)}^{\displaystyle \frac{1}{2}}(1-{\alpha }_{k}){{\bf{e}}}_{{\bf{k}}\lambda }\,{a}_{{\bf{k}}\lambda }{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{x}}}+{\rm{H}}.{\rm{c}}.\end{eqnarray} \tag{ 5 }$$

The coupling constant g_{k λ} and the (real) coefficients ${u}_{k}^{\pm }$ are defined as

$$\begin{eqnarray}&&{g}_{{\bf{k}}\lambda }=-{\rm{i}}{\left(\displaystyle \frac{{\omega }_{0}}{2{L}^{3}}\right)}^{\displaystyle \frac{1}{2}}\,{{\bf{e}}}_{{\bf{k}}\lambda }\cdot {\bf{d}},\qquad {u}_{k}^{\pm }=(1-{\alpha }_{k}){\left(\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}\right)}^{1/2}\mp {\alpha }_{k}{\left(\displaystyle \frac{{\omega }_{k}}{{\omega }_{0}}\right)}^{1/2}\end{eqnarray} \tag{ 6 }$$

where ${\bf{d}}$ and ${\omega }_{0}$ denote the two-level transition dipole moment and transition frequency, respectively. The real numbers α_k can be chosen arbitrarily. Choosing α_k = 0 yields the Coulomb-gauge Hamiltonian while choosing α_k = 1 yields the multipolar-gauge Hamiltonian. Letting α_k = 1 in equation (2) yields the well-known Power–Zienau–Woolley (PZW) transformation that relates the Coulomb and multipolar gauges. While the relation between the Coulomb and multipolar gauge has been discussed extensively [27–31], the PZW transformation is in fact a special case of a broader class of unitary gauge-fixing transformations [31–33]. More generally still, the freedom to choose the α_k within the canonical transformation (2) implies redundancy within our mathematical description and is henceforth referred to as generalised gauge-freedom. A third special case of equation (3) is afforded by making the choice α_k = ω₀/(ω₀ + ω_k), which specifies a symmetric mixture of Coulomb and multipolar couplings. This representation has proved useful in both photo-detection theory [31, 34] and open quantum systems theory [26], because within this representation ${u}_{k}^{+}\equiv 0$. The counter-rotating terms in the linear dipole-field interaction term V − V⁽²⁾ are thereby eliminated without use of the rotating-wave approximation.

Given the above arbitrary generalised-gauge description, it is clear that arbitrary matrix elements ${M}_{{fi}}(t)=\langle f| M(t)| i\rangle$ between eigenstates $| f\rangle$ and $| i\rangle$ of H₀ will not be the same when the evolution of the operator M is determined by different total Hamiltonians H and H', that have been obtained by making different choices of α_k in equation (3). In contrast on-energy-shell QED S-matrix elements are necessarily the same for two interaction Hamiltonians V and V', constrained such that the corresponding total Hamiltonians are related by a unitary transformation e^iT as [25, 28, 35]

$$\begin{eqnarray}&&H={H}_{0}+V,\qquad H^{\prime} ={{\rm{e}}}^{{\rm{i}}{T}}{{H}{\rm{e}}}^{-{\rm{i}}{T}}={H}_{0}+V^{\prime} .\end{eqnarray} \tag{ 7 }$$

For gauge-invariance to hold the unperturbed Hamiltonian H₀ must be identified as the same operator before and after the transformation by e^iT. Note however, that the unperturbed Hamiltonian H₀ given in equation (1) does not commute with the unitary transformation ${R}_{\{{\alpha }_{k}\}}$ given in equation (2) meaning that this H₀ represents a different physical observable depending on the choice of interaction. Despite this, S-matrix elements based on the partition H = H₀ + V are invariant, because H₀ in equation (1) does not explicitly depend on the α_k and is therefore the same for each different choice of generalised-gauge.

Having determined the conditions under which QED matrix elements are gauge-invariant we now turn our attention to deriving a master equation describing the two-level dipole within the radiation field. The conventional derivation of the second order quantum optical master equation, as found in [36] for example, does not at any point involve self-energy contributions due to the V⁽²⁾ term within the interaction V. In general however, this term does contribute to dipole level-shifts, as is shown in appendix A.1. In the general case that the temperature of the radiation field is arbitrary, the self-energy contributions from V⁽²⁾ can be incorporated into the master equation by defining the Hamiltonian

$$\begin{eqnarray}&&{\tilde{H}}_{d}=({\omega }_{0}+{\delta }_{e}^{(2)}-{\delta }_{g}^{(2)}){\sigma }^{+}{\sigma }^{-},\end{eqnarray} \tag{ 8 }$$

where the excited and ground state self-energy shifts are defined as

$$\begin{eqnarray}&&{\delta }_{n}^{(2)}=\mathrm{tr}(V| n\rangle \langle n| \otimes {\rho }_{F}^{\mathrm{eq}})=\mathrm{tr}({V}^{(2)}| n\rangle \langle n| \otimes {\rho }_{F}^{\mathrm{eq}}),\end{eqnarray} \tag{ 9 }$$

in which n = e, g and ${\rho }_{F}^{\mathrm{eq}}(\beta )={{\rm{e}}}^{-\beta {\sum }_{{\bf{k}}\lambda }{\omega }_{k}{a}_{{\bf{k}}\lambda }^{\dagger }{a}_{{\bf{k}}\lambda }}/\mathrm{tr}({{\rm{e}}}^{-\beta {\sum }_{{\bf{k}}\lambda }{\omega }_{k}{a}_{{\bf{k}}\lambda }^{\dagger }{a}_{{\bf{k}}\lambda }})$ with $\beta$ the inverse temperature of the radiation field. Since V⁽²⁾ is gauge-dependent we cannot include the self-energy shifts within the unperturbed Hamiltonian H₀ without ruining the gauge-invariance of any S-matrix elements obtained using the unperturbed states. Instead we replace the free system Hamiltonian in the usual Born–Markov master equation with the shifted Hamiltonian ${\tilde{H}}_{d}$ directly to obtain the second order master equation

$$\begin{eqnarray}&&\dot{\rho }=-{\rm{i}}[{\tilde{H}}_{d},\rho ]-{{\rm{e}}}^{-{{\rm{i}}{H}}_{0}t}{\int }_{0}^{\infty }{\rm{d}}{s}\,{\mathrm{tr}}_{F}[{V}_{I}(t),[{V}_{I}(t-s),{\rho }_{I}(t)\otimes {\rho }_{F}^{\mathrm{eq}}]]{{\rm{e}}}^{{{\rm{i}}{H}}_{0}t}.\end{eqnarray} \tag{ 10 }$$

This master equation automatically includes the level shifts due to V⁽²⁾ within the unitary evolution part, but the rest of the master equation is expressed in terms of the original partition H = H₀ + V. Using this partition where H₀ and V are given in equations (1) and (3) respectively, equation (10) yields the α_k-independent result

$$\begin{eqnarray}&&\dot{\rho }=-{\rm{i}}[{\tilde{\omega }}_{0}{\sigma }^{+}{\sigma }^{-},\rho ]+\gamma (N+1)\left({\sigma }^{-}\rho {\sigma }^{+}-\displaystyle \frac{1}{2}\{{\sigma }^{+}{\sigma }^{-},\rho \}\right)+\gamma N\left({\sigma }^{+}\rho {\sigma }^{-}-\displaystyle \frac{1}{2}\{{\sigma }^{-}{\sigma }^{+},\rho \}\right).\end{eqnarray} \tag{ 11 }$$

Here $N=1/({{\rm{e}}}^{\beta {\omega }_{0}}-1)$, ${\tilde{\omega }}_{0}={\omega }_{0}+{\rm{\Delta }}$ and

$$\begin{eqnarray}\gamma & = & 2\pi \displaystyle \sum _{{\bf{k}}\lambda }| \langle {\bf{k}}\lambda ,g| V| 0,e\rangle {| }^{2}\delta ({\omega }_{k}-{\omega }_{0})=\displaystyle \frac{{\omega }_{0}^{3}| {\bf{d}}{| }^{2}}{3\pi },\\ {\rm{\Delta }} & = & \displaystyle \int \displaystyle \frac{{{\rm{d}}}^{3}k}{{(2\pi )}^{3}}\displaystyle \sum _{\lambda }| {{\bf{e}}}_{\lambda }({\bf{k}})\cdot {\bf{d}}{| }^{2}(1+2{N}_{k})\displaystyle \frac{{\omega }_{0}^{3}}{{\omega }_{k}({\omega }_{0}^{2}-{\omega }_{k}^{2})},\end{eqnarray} \tag{ 12 }$$

where the continuum limit for wavevectors k has been applied and ${N}_{k}=1/({{\rm{e}}}^{\beta {\omega }_{k}}-1)$. Further details of the calculations leading to the final result for Δ in equation (12) are given in appendix A.1. We note that for N_k = 0 we have

$$\begin{eqnarray}{\rm{\Delta }} & = & {\tilde{\omega }}_{0}-{\omega }_{0}=\langle 0,e| V| 0,e\rangle +\displaystyle \sum _{{\bf{k}}\lambda }\displaystyle \frac{| \langle 0,e| V| {\bf{k}}\lambda ,g\rangle {| }^{2}}{{\omega }_{0}-{\omega }_{k}}-\langle 0,g| V| 0,g\rangle +\displaystyle \sum _{{\bf{k}}\lambda }\displaystyle \frac{| \langle 0,g| V| {\bf{k}}\lambda ,e\rangle {| }^{2}}{{\omega }_{0}+{\omega }_{k}}\\ & = & \displaystyle \int \displaystyle \frac{{{\rm{d}}}^{3}k}{{(2\pi )}^{3}}\displaystyle \sum _{\lambda }| {{\bf{e}}}_{{\bf{k}}\lambda }\cdot {\bf{d}}{| }^{2}\displaystyle \frac{{\omega }_{0}^{3}}{{\omega }_{k}({\omega }_{0}^{2}-{\omega }_{k}^{2})}.\end{eqnarray} \tag{ 13 }$$

The α_k-independence (generalised gauge-invariance) of the master equation (11) can be understood by noting that the spontaneous emission rate γ and level-shift Δ in equation (13) are gauge-invariant QED matrix elements that can be obtained directly using second order perturbation theory.

In summary, we have shown that the master equation obtained from different, unitarily equivalent microscopic Hamiltonians is the same provided it depends only on S-matrix elements. S-matrix elements are invariant if the bare Hamiltonian H₀ is kept the same for each choice of total Hamiltonian. In what follows this will be seen to be significant for the derivation of the master equation describing two strongly coupled dipoles.

Let us now turn our attention to obtaining the analogous result to equation (11) for the case of two identical interacting dipoles at positions R₁ and R₂ within a common radiation reservoir. The transition dipole moments and transition frequencies of the dipoles are independent of the dipole label and are denoted d and ω₀ respectively. The wavelength corresponding to ${\omega }_{0}$ is denoted with ${\lambda }_{0}$. An important quantity in the two-dipole system dynamics is the inter-dipole separation $R=| {{\bf{R}}}_{2}-{{\bf{R}}}_{1}|$. In terms of R we can identify in the usual way three distinct parameter regimes: $R\ll {\lambda }_{0}$ is the near zone in which R⁻³-dependent terms dominate, $R\sim {\lambda }_{0}$ is the intermediate zone in which R⁻²-dependent terms may become significant, and $R\gg {\lambda }_{0}$ is the far-zone (radiation zone) in which R⁻¹-dependent terms dominate. We give a general derivation of the standard two dipole master equation, in which the gauge freedom within the microscopic Hamiltonian is left open throughout. This reveals limitations within the conventional derivation using the multipolar gauge. To begin we define the two-dipole generalised PZW gauge transformation by [26, 34]

$$\begin{eqnarray}&&{R}_{\{{\alpha }_{k}\}}:= \exp \left({\rm{i}}\displaystyle \sum _{\mu =1}^{2}{{\bf{d}}}_{\mu }\cdot {{\bf{A}}}_{\{{\alpha }_{k}\}}({{\bf{R}}}_{\mu })\right),\end{eqnarray} \tag{ 14 }$$

where d_μ denotes the dipole moment operator of the μ'th dipole. We now transform the dipole approximated Coulomb gauge Hamiltonian using ${R}_{\{{\alpha }_{k}\}}$ and afterwards make the two-level approximation for each dipole. This implies that ${{\bf{d}}}_{\mu }={\bf{d}}({\sigma }_{\mu }^{+}+{\sigma }_{\mu }^{-})$, and that the Hamiltonian can be written H = H₀ + V with

$$\begin{eqnarray}&&{H}_{0}=\displaystyle \sum _{\mu =1}^{2}{\omega }_{0}{\sigma }_{\mu }^{+}{\sigma }_{\mu }^{-}+\displaystyle \sum _{{\bf{k}}\lambda }{\omega }_{k}{a}_{{\bf{k}}\lambda }^{\dagger }{a}_{{\bf{k}}\lambda }\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&V=\left[\displaystyle \sum _{{\bf{k}}\lambda }\displaystyle \sum _{\mu =1}^{2}{g}_{{\bf{k}}\lambda }{\sigma }_{\mu }^{+}({u}_{k}^{+}\,{a}_{{\bf{k}}\lambda }^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{R}}}_{\mu }}+{u}_{k}^{-}\,{a}_{{\bf{k}}\lambda }{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {{\bf{R}}}_{\mu }})+{\rm{h}}.{\rm{c}}.\right]+{V}_{\{{\alpha }_{k}\}}^{(2)}+{C}_{\{{\alpha }_{k}\}}.\end{eqnarray} \tag{ 16 }$$

Analogously to the single-dipole case the first term in equation (16) defines a linear dipole-field interaction component while the term ${V}_{\{{\alpha }_{k}\}}^{(2)}$ consists of self-energy contributions for each dipole and the radiation field;

$$\begin{eqnarray}&&{V}_{\{{\alpha }_{k}\}}^{(2)}=\displaystyle \sum _{{\bf{k}}\lambda }\displaystyle \frac{1}{{L}^{3}}\,{\alpha }_{k}^{2}\,| {{\bf{e}}}_{{\bf{k}}\lambda }\cdot {\bf{d}}{| }^{2}+\displaystyle \sum _{\mu =1}^{2}\displaystyle \frac{{e}^{2}}{2m}\,\tilde{{\bf{A}}}{({{\bf{R}}}_{\mu })}^{2},\end{eqnarray} \tag{ 17 }$$

where m is the dipole mass. Due to the two-level approximation the first term in (17) is proportional to the identity, while the second term is a radiation self-energy term. The field $\tilde{A}$ is defined as in the singe-dipole case by equation (5). The final term in equation (16) ${C}_{\{{\alpha }_{k}\}}$, has no analogue in the single-dipole Hamiltonian. This term gives a static Coulomb-like interaction between the dipoles, which is independent of the field;

$$\begin{eqnarray}&&{C}_{\{{\alpha }_{k}\}}=\displaystyle \sum _{{\bf{k}}\lambda }\displaystyle \frac{1}{{L}^{3}}\,({\alpha }_{k}^{2}-1)\,| {{\bf{e}}}_{{\bf{k}}\lambda }\cdot {\bf{d}}{| }^{2}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{R}}}{\sigma }_{1}^{x}{\sigma }_{2}^{x},\end{eqnarray} \tag{ 18 }$$

where ${\sigma }_{\mu }^{x}={\sigma }_{\mu }^{+}+{\sigma }_{\mu }^{-}$. In the Coulomb gauge (α_k = 0) the term ${C}_{\{{\alpha }_{k}\}}$ reduces to the usual dipole–dipole Coulomb interaction. In the multipolar gauge ${C}_{\{{\alpha }_{k}\}}$ vanishes, and the interaction in equation (16) therefore reduces to a sum of interaction terms for each dipole.

It is important to note that as in the single-dipole case ${R}_{\{{\alpha }_{k}\}}$ does not commute with H₀ given in equation (15) implying that H₀ represents a different physical observable for each choice of α_k. More generally, since ${R}_{\{{\alpha }_{k}\}}$ is a non-local transformation, which mixes material and transverse field degrees of freedom, the canonical material and field operators are different for each choice of α_k. This implies that the master equation for the dipoles will generally be different for each choice of α_k. We can, however, obtain a gauge-invariant result by ensuring that the master equation depends only on gauge-invariant S-matrix elements. These matrix elements are gauge-invariant despite the implicit difference in the material and field degrees of freedom within each generalised gauge. Usually the two-dipole master equation is derived using the specific choice α_k = 1 (multipolar gauge) for which the direct Coulomb-like coupling ${C}_{\{{\alpha }_{k}\}}$ vanishes identically. To obtain the same master equation for any other choice of ${\alpha }_{k}\ne 1$, we must include ${C}_{\{{\alpha }_{k}\}}$ within the interaction Hamiltonian V. The reason is that H₀ must be identified as the same operator for each choice of α_k in order that the gauge-invariance of the associated S-matrix holds.

We now proceed with a direct demonstration that the standard two-dipole master equation can indeed be obtained for any other choice of α_k, provided ${C}_{\{{\alpha }_{k}\}}$ is kept within the interaction Hamiltonian V. To do this we substitute the interaction Hamiltonian in equation (16) into the second order Born–Markov master equation in the interaction picture with respect to H₀, which is given by

$$\begin{eqnarray}&&{\dot{\rho }}_{I}(t)=-{\rm{i}}{{\rm{tr}}}_{F}[{V}_{I}(t),{\rho }_{I}(t)\otimes {\rho }_{F}^{{\rm{eq}}}]-{\int }_{0}^{\infty }{\rm{d}}{s}\,{{\rm{tr}}}_{F}[{V}_{I}(t),[{V}_{I}(t-s),{\rho }_{I}(t)\otimes {\rho }_{F}^{{\rm{eq}}}]],\end{eqnarray} \tag{ 19 }$$

where V_I(t) denotes the interaction Hamiltonian in the interaction picture and ρ_I(t) denotes the interaction picture state of the two dipoles. We retain contributions up to order e² and perform a further secular approximation, which neglects terms oscillating with twice the transition frequency ω₀. Transforming back to the Schrödinger picture and including the single-dipole self-energy contributions as in equation (10), we arrive after lengthy but straightforward manipulations at the final α_k-independent result

$$\begin{eqnarray}\dot{\rho } & = & -{\rm{i}}{\tilde{\omega }}_{0}\displaystyle \sum _{\mu =1}^{2}[{\sigma }_{\mu }^{+}{\sigma }_{\mu }^{-},\rho ]-{\rm{i}}{{\rm{\Delta }}}_{12}\displaystyle \sum _{\mu \ne \nu }^{2}[{\sigma }_{\mu }^{+}{\sigma }_{\nu }^{-},\rho ]+\displaystyle \sum _{\mu ,\nu =1}^{2}{\gamma }_{\mu \nu }\left[(N+1)\left({\sigma }_{\mu }^{-}\rho {\sigma }_{\nu }^{+}-\displaystyle \frac{1}{2}\{{\sigma }_{\mu }^{+}{\sigma }_{\nu }^{-},\rho \}\right)\right.\\ & & +\left.N\left({\sigma }_{\mu }^{+}\rho {\sigma }_{\nu }^{-}-\displaystyle \frac{1}{2}\{{\sigma }_{\mu }^{-}{\sigma }_{\nu }^{+},\rho \}\right)\right].\end{eqnarray} \tag{ 20 }$$

This equation is identical in form to the standard two-dipole master equation, which can be found in [4] for example. The coefficients within the master equation are as follows. The decay rates γ_μν are given by

$$\begin{eqnarray}{\gamma }_{\mu \mu } & = & \gamma ,\\ {\gamma }_{12} & = & {\gamma }_{21}=2\pi \displaystyle \sum _{{\bf{k}}\lambda }\langle 0,e,g| V| {\bf{k}}\lambda ,g,g\rangle \langle {\bf{k}}\lambda ,g,g| V| 0,g,e\rangle \delta ({\omega }_{0}-{\omega }_{k})={d}_{i}{d}_{j}{\tau }_{{ij}}({\omega }_{0},R),\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&{\tau }_{{ij}}({\omega }_{0},R)=\displaystyle \frac{{\omega }_{0}^{3}}{2\pi }\left(({\delta }_{{ij}}-{\hat{R}}_{i}{\hat{R}}_{j})\displaystyle \frac{\sin {\omega }_{0}R}{{\omega }_{0}R}+({\delta }_{{ij}}-3{\hat{R}}_{i}{\hat{R}}_{j})\left[\displaystyle \frac{\cos {\omega }_{0}R}{{({\omega }_{0}R)}^{2}}-\displaystyle \frac{\sin {\omega }_{0}R}{{({\omega }_{0}R)}^{3}}\right]\right).\end{eqnarray} \tag{ 22 }$$

In equation (21) and throughout we denote spatial components with Latin indices and adopt the convention that repeated Latin indices are summed. The quantity γ₁₂ denotes an R-dependent collective decay rate. The third equality in equation (21) wherein γ₁₂ has been expressed as a matrix element involving V makes the reason for the α_k-independence of this rate clear. We now turn our attention to the master equation shifts ${\rm{\Delta }}={\tilde{\omega }}_{0}-{\omega }_{0}$ and Δ₁₂. The single-dipole shift Δ includes all self-energy contributions, which have been dealt with in the same way as for the single-dipole master equation (see equation (10)). The shift Δ is therefore as in equation (12). Details of the calculation of the joint shift Δ₁₂ are given in appendix A.2 with the final result being

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{12}=\displaystyle \sum _{n}\displaystyle \frac{\langle f| V| n\rangle \langle n| V| {\rm{i}}\rangle }{{\omega }_{i}-{\omega }_{n}}=\int \displaystyle \frac{{{\rm{d}}}^{3}k}{{(2\pi )}^{3}}\displaystyle \sum _{\lambda }| {{\bf{e}}}_{\lambda }({\bf{k}})\cdot {\bf{d}}{| }^{2}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{R}}}\displaystyle \frac{{\omega }_{k}^{2}}{{\omega }_{0}^{2}-{\omega }_{k}^{2}}={d}_{i}{d}_{j}{V}_{{ij}}({\omega }_{0},R),\end{eqnarray} \tag{ 23 }$$

where $| f\rangle =| g,e,0\rangle$, $| i\rangle =| e,g,0\rangle$ and $| n\rangle$ are eigenstates of H₀ and

$$\begin{eqnarray}&&{V}_{{ij}}({\omega }_{0},R)=-\displaystyle \frac{{\omega }_{0}^{3}}{4\pi }\left(({\delta }_{{ij}}-{\hat{R}}_{i}{\hat{R}}_{j})\displaystyle \frac{\cos {\omega }_{0}R}{{\omega }_{0}R}-({\delta }_{{ij}}-3{\hat{R}}_{i}{\hat{R}}_{j})\left[\displaystyle \frac{\sin {\omega }_{0}R}{{({\omega }_{0}R)}^{2}}+\displaystyle \frac{\cos {\omega }_{0}R}{{({\omega }_{0}R)}^{3}}\right]\right).\end{eqnarray} \tag{ 24 }$$

As indicated by the second equality in equation (23) Δ₁₂ is nothing but the well-known gauge-invariant QED matrix-element describing resonant energy-transfer.

We have therefore obtained the standard result, equation (20), without ever making a concrete choice for the α_k. In order that the standard result is obtained the direct Coulomb-like interaction ${C}_{\{{\alpha }_{k}\}}$ must be kept within the interaction Hamiltonian V. This ensures that for all α_k the unperturbed Hamiltonian H₀ is that used in conventional derivations wherein α_k = 1. The S-matrix elements involving ${C}_{\{{\alpha }_{k}\}}$ that appear as coefficients in the master equation are then α_k-independent and are the same as those obtained in the conventional derivation. Our derivation makes it clear that when ${C}_{\{{\alpha }_{k}\}}$ is sufficiently strong compared with the linear dipole-field coupling term, its inclusion within the interaction Hamiltonian rather than H₀ may be ill-justified. The standard master equation may therefore be inaccurate in such regimes. This fact is obscured within conventional derivations that use the multipolar gauge α_k = 1, because in this gauge ${C}_{\{{\alpha }_{k}\}}$ vanishes identically. However, one typically still assumes weak-coupling to the radiation field in the multipolar gauge, and this leads to the standard master equation (20). If instead we adopt the Coulomb gauge α_k = 0 we obtain the static dipole–dipole interaction ${C}_{\{0\}}=C{\sigma }_{1}^{x}{\sigma }_{2}^{x}$, where in the mode continuum limit

$$\begin{eqnarray}&&C=\displaystyle \frac{{d}_{i}{d}_{j}}{4\pi {R}^{3}}({\delta }_{{ij}}-3{\hat{R}}_{i}{\hat{R}}_{j}).\end{eqnarray} \tag{ 25 }$$

This quantity coincides with the near-field limit of the resonant energy transfer element Δ₁₂ given in equation (23). In the near-field regime $R/\lambda \ll 1$, ${C}_{\{0\}}$ may be too strong to be kept within the purportedly weak perturbation V and the standard master equation, which only results when one treats C_{0} as a weak perturbation, should then break down. This will be discussed in more detail in the following section.

In the near-field regime $R/{\lambda }_{0}\ll 1$ the rate of spontaneous emission into the transverse field is much smaller than the direct dipolar coupling; $C/\gamma \gg 1$. Moreover, for a system of closely spaced Rydberg atoms, the electrostatic Coulomb interaction may be such that $C\sim {\omega }_{0}$. For example, given a Rydberg state with principal quantum number n = 50 we can estimate the maximum associated dipole moment as $(3/2){n}^{2}{a}_{0}e\sim {10}^{-26}\,\mathrm{Cm}$ where a₀ is the Bohr radius and e the electronic charge. For a $1\mu {\rm{m}}$ separation, which is approximately equal to $10{n}^{2}{a}_{0}$, the electrostatic dipole interaction $C/{\omega }_{0}\sim 1$ for ω₀ corresponding to a microwave frequency. In such situations it is not clear that the Coulomb interaction can be included within the perturbation V with the coupling to the transverse field. In the multipolar gauge where no direct Coulomb interaction is explicit the same physical interaction is mediated by the low frequency transverse modes, which must be handled carefully. A procedure which separates out these modes should ultimately result in a separation of the Coulomb interaction, which is of course already explicit within the Coulomb gauge. We remark that when considering realistic Rydberg atomic systems within the strong dipole–dipole coupling regime the validity of the two-level model should also be considered. However, moving beyond the two-level approximation is beyond the scope of this paper. Our aim is to consider general atomic, molecular and condensed matter systems strongly-coupled by dipole–dipole interactions for which two-level models are typically used [4, 37, 38]. Retaining the two-level model for each dipole allows us to succinctly compare with existing literature and thereby determine the relative difference produced by our non-perturbative treatment of dipole–dipole interactions.

In the Coulomb gauge the interaction Hamiltonian V coupling to the transverse radiation field is

$$\begin{eqnarray}&&V=\displaystyle \sum _{\mu =1}^{2}{\omega }_{0}{\sigma }_{\mu }^{y}{\bf{d}}\cdot {\bf{A}}({{\bf{R}}}_{\mu })+\displaystyle \sum _{\mu =1}^{2}\displaystyle \frac{{e}^{2}}{2m}{\bf{A}}{({{\bf{R}}}_{\mu })}^{2}\end{eqnarray} \tag{ 26 }$$

with ${\sigma }_{\mu }^{y}=-{\rm{i}}({\sigma }_{\mu }^{+}-{\sigma }_{\mu }^{-})$ and

$$\begin{eqnarray}&&{\bf{A}}({\bf{x}})=\displaystyle \sum _{{\bf{k}}\lambda }\sqrt{\displaystyle \frac{1}{2{\omega }_{k}{L}^{3}}}{{\bf{e}}}_{{\bf{k}}\lambda }{a}_{{\bf{k}}\lambda }^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{x}}}+{\rm{H}}.{\rm{c}}.\end{eqnarray} \tag{ 27 }$$

The contribution of the transverse field to Δ₁₂ in equation (21) is found using equation (26) to be

$$\begin{eqnarray}{\tilde{{\rm{\Delta }}}}_{12} & = & \displaystyle \sum _{n}\displaystyle \frac{\langle 0,e,g| V| n\rangle \langle n| V| 0,g,e\rangle }{{\omega }_{0}-{\omega }_{n}}\\ & = & -\displaystyle \frac{{\omega }_{0}^{3}{d}_{i}{d}_{j}}{4\pi }\left(({\delta }_{{ij}}-{\hat{R}}_{i}{\hat{R}}_{j})\displaystyle \frac{\cos {\omega }_{0}R}{{\omega }_{0}R}-({\delta }_{{ij}}-3{\hat{R}}_{i}{\hat{R}}_{j})\left[\displaystyle \frac{\sin {\omega }_{0}R}{{({\omega }_{0}R)}^{2}}-\displaystyle \frac{1-\cos {\omega }_{0}R}{{({\omega }_{0}R)}^{3}}\right]\right)\\ & = & \,{{\rm{\Delta }}}_{12}-C.\end{eqnarray} \tag{ 28 }$$

When the contribution $C=\langle 0,e,g| {C}_{\{0\}}| 0,g,e\rangle$ resulting from the direct Coulomb interaction is added to ${\tilde{{\rm{\Delta }}}}_{12}$ the fully retarded result Δ₁₂ is obtained. The two matrix elements Δ₁₂ and ${\tilde{{\rm{\Delta }}}}_{12}$ therefore only differ in their near-field components, which vary as R⁻³ and which we denote by ${{\rm{\Delta }}}_{12}^{\mathrm{nf}}$ and ${\tilde{{\rm{\Delta }}}}_{12}^{\mathrm{nf}}$, respectively. According to equations (23) and (28) the components ${{\rm{\Delta }}}_{12}^{\mathrm{nf}}$ and ${\tilde{{\rm{\Delta }}}}_{12}^{\mathrm{nf}}$ dominate at low frequencies ω₀. Since Δ₁₂ is evaluated at resonance ω_k = ω₀, it follows that within the multipolar-gauge the low ω_k modes within the system-reservoir coupling give rise to a strong dipole–dipole interaction in the form of ${{\rm{\Delta }}}_{12}^{\mathrm{nf}}\approx C$. In such regimes the multipolar interaction Hamiltonian cannot be classed as a weak perturbation. On the other hand, the matrix element ${\tilde{{\rm{\Delta }}}}_{12}$ is obtained using the Coulomb gauge interaction in equation (26), and is such that ${\tilde{{\rm{\Delta }}}}_{12}^{\mathrm{nf}}={{\rm{\Delta }}}_{12}^{\mathrm{nf}}-C\approx 0$. Within the Coulomb gauge the interaction equivalent to the low frequency part of the multipolar gauge system-reservoir coupling is a direct dipole–dipole Coulomb interaction C_{0}. This appears explicitly in the Hamiltonian, but has not been included within equation (26), which therefore represents a genuinely weak perturbation.

The collective decay rate ${\gamma }_{12}$ as given in equation (21) does not involve the direct Coulomb interaction C_{0} in any way, and can be obtained from the transverse field interaction in equation (26) or from the multipolar interaction. Crucially, in the near-field regime $R/{\lambda }_{0}\ll 1$ the terms γ, γ₁₂ and ${\tilde{{\rm{\Delta }}}}_{12}$, which result from the interaction in equation (26), are several orders of magnitude smaller than the direct electrostatic coupling C. Motivated by the discussion above, we include the Coulomb interaction within the unperturbed Hamiltonian, but continue to treat the interaction with the transverse field as a weak perturbation. This gives rise to a master equation depending on different S-matrix elements.

The unperturbed Hamiltonian H₀ = H_d + H_F is defined by

$$\begin{eqnarray}&&{H}_{d}=\displaystyle \sum _{\mu =1}^{2}{\omega }_{0}{\sigma }_{\mu }^{+}{\sigma }_{\mu }^{-}+C{\sigma }_{1}^{x}{\sigma }_{2}^{x},\qquad {H}_{F}=\displaystyle \sum _{{\bf{k}}\lambda }{\omega }_{k}{a}_{{\bf{k}}\lambda }^{\dagger }{a}_{{\bf{k}}\lambda },\end{eqnarray} \tag{ 29 }$$

where $C\in {\mathbb{R}}$ is given by equation (25). The corresponding interaction Hamiltonian is then given in equation (26). We begin by diagonalising H_d as

$$\begin{eqnarray}&&{H}_{d}=\displaystyle \sum _{n=1}^{4}{\epsilon }_{n}| {\epsilon }_{n}\rangle \langle {\epsilon }_{n}| ,\end{eqnarray} \tag{ 30 }$$

where

$$\begin{eqnarray}&&{\epsilon }_{1}={\omega }_{0}-\eta ,\qquad {\epsilon }_{2}={\omega }_{0}-C\qquad {\epsilon }_{3}={\omega }_{0}+C,\qquad {\epsilon }_{4}={\omega }_{0}+\eta ,\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}| {\epsilon }_{1}\rangle & = & \displaystyle \frac{1}{\sqrt{{C}^{2}+{({\omega }_{0}+\eta )}^{2}}}([{\omega }_{0}+\eta ]| g,g\rangle -C| e,e\rangle ),\qquad | {\epsilon }_{2}\rangle =\displaystyle \frac{1}{\sqrt{2}}(| e,g\rangle -| g,e\rangle ),\\ | {\epsilon }_{3}\rangle & = & \displaystyle \frac{1}{\sqrt{2}}(| e,g\rangle +| g,e\rangle ),\qquad | {\epsilon }_{4}\rangle =\displaystyle \frac{1}{\sqrt{{C}^{2}+{({\omega }_{0}-\eta )}^{2}}}(C| e,e\rangle -[{\omega }_{0}-\eta ]| g,g\rangle ),\end{eqnarray} \tag{ 32 }$$

with $\eta =\sqrt{{\omega }_{0}^{2}+{C}^{2}}$. Next we move into the interaction picture with respect to H₀ and substitute the interaction picture interaction Hamiltonian into equation (19). Moving back into the Schrödinger picture we eventually obtain

$$\begin{eqnarray}&&\dot{\rho }=-{\rm{i}}[{H}_{d},\rho ]+\displaystyle \sum _{\zeta ,\zeta ^{\prime} =\pm {\omega }_{\mathrm{1,2}}}\displaystyle \sum _{\mu ,\nu =1}^{2}[{{\rm{\Gamma }}}_{\mu \nu }(\zeta )({A}_{\nu \zeta }\rho {A}_{\mu \zeta ^{\prime} }^{\dagger }-{A}_{\mu \zeta ^{\prime} }^{\dagger }{A}_{\nu \zeta }\rho )+{\rm{h}}.{\rm{c}}.].\end{eqnarray} \tag{ 33 }$$

Here, ω₁ = η − C and ω₂ = η + C, while ${A}_{\mu (-\zeta )}={A}_{\mu \zeta }^{\dagger }$ and ${A}_{\mu {\zeta }_{n}}\equiv {A}_{\mu n}\,(n=1,2)$ with

$$\begin{eqnarray}{A}_{11} & = & a| {\epsilon }_{1}\rangle \langle {\epsilon }_{2}| +b| {\epsilon }_{3}\rangle \langle {\epsilon }_{4}| ,\qquad \,\,\,{A}_{12}=c| {\epsilon }_{1}\rangle \langle {\epsilon }_{3}| -d| {\epsilon }_{2}\rangle \langle {\epsilon }_{4}| ,\\ {A}_{21} & = & -a| {\epsilon }_{1}\rangle \langle {\epsilon }_{2}| +b| {\epsilon }_{3}\rangle \langle {\epsilon }_{4}| ,\qquad {A}_{22}=c| {\epsilon }_{1}\rangle \langle {\epsilon }_{3}| +d| {\epsilon }_{2}\rangle \langle {\epsilon }_{4}| ,\end{eqnarray} \tag{ 34 }$$

where

$$\begin{eqnarray}a & = & \displaystyle \frac{{\omega }_{0}+\eta -C}{\sqrt{2({C}^{2}+{[{\omega }_{0}+\eta ]}^{2})}},\qquad b=\displaystyle \frac{{\omega }_{0}-\eta -C}{\sqrt{2({C}^{2}+{[{\omega }_{0}-\eta ]}^{2})}},\\ c & = & \displaystyle \frac{{\omega }_{0}+\eta +C}{\sqrt{2({C}^{2}+{[{\omega }_{0}+\eta ]}^{2})}},\qquad d=\displaystyle \frac{-{\omega }_{0}+\eta +C}{\sqrt{2({C}^{2}+{[{\omega }_{0}-\eta ]}^{2})}}.\end{eqnarray} \tag{ 35 }$$

The coefficients Γ_μν(ω) are defined by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mu \nu }(\omega )={\omega }_{0}^{2}{d}_{{\rm{i}}}{d}_{j}{\int }_{0}^{\infty }{\rm{d}}{s}\,{{\rm{e}}}^{{\rm{i}}\omega s}\langle {A}_{I,{\rm{i}}}({{\bf{R}}}_{\mu },s){A}_{I,j}({{\bf{R}}}_{\nu },0){\rangle }_{\beta },\end{eqnarray} \tag{ 36 }$$

where A_I(x, t) denotes the field A(x) in equation (27) once transformed into the interaction picture, and $\langle \cdot {\rangle }_{\beta }$ denotes the average with respect to the radiation thermal state at temperature β⁻¹. The Γ_μν(ω) are symmetric ${{\rm{\Gamma }}}_{\mu \nu }(\omega )={{\rm{\Gamma }}}_{\nu \mu }(\omega )$ and can be written

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mu \nu }(\omega )=\displaystyle \frac{1}{2}{\gamma }_{\mu \nu }(\omega )+{{\rm{i}}{S}}_{\mu \nu }(\omega ),\end{eqnarray} \tag{ 37 }$$

where

$$\begin{eqnarray}{\gamma }_{\mu \mu }(\omega ) & = & (1+N)\gamma \displaystyle \frac{\omega }{{\omega }_{0}},\qquad {\gamma }_{12}(\omega )=(1+N){d}_{i}{d}_{j}{\tau }_{{ij}}(\omega ,R)\displaystyle \frac{{\omega }_{0}^{2}}{{\omega }^{2}},\\ {S}_{\mu \mu }(\omega ) & = & \displaystyle \frac{\gamma }{2\pi }{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{k}\displaystyle \frac{{\omega }_{k}}{{\omega }_{0}}\left[\displaystyle \frac{1+{N}_{k}}{\omega -{\omega }_{k}}+\displaystyle \frac{{N}_{k}}{\omega +{\omega }_{k}}\right],\\ {S}_{12}(\omega ) & = & \displaystyle \frac{{d}_{i}{d}_{j}}{2\pi }{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{k}\,{\tau }_{{ij}}({\omega }_{k},R)\displaystyle \frac{{\omega }_{0}^{2}}{{\omega }_{k}^{2}}\left[\displaystyle \frac{1+{N}_{k}}{\omega -{\omega }_{k}}+\displaystyle \frac{{N}_{k}}{\omega +{\omega }_{k}}\right].\end{eqnarray} \tag{ 38 }$$

The frequency integrals in equation (38) are to be understood as principal values. The decay rates γ_μν(ω) in equation (38) coincide with those found in the standard master equation (20) when evaluated at ω₀, though are here evaluated at the frequencies ω_1,2. The quantities S_μν are related to the shifts Δ and Δ₁₂, defined in equations (12) and (23) respectively, by

$$\begin{eqnarray}&&{\rm{\Delta }}={S}_{\mu \mu }({\omega }_{0})-{S}_{\mu \mu }(-{\omega }_{0}),\qquad {{\rm{\Delta }}}_{12}-C={\tilde{{\rm{\Delta }}}}_{12}={S}_{12}({\omega }_{0})+{S}_{12}(-{\omega }_{0}).\end{eqnarray} \tag{ 39 }$$

In deriving equation (33) we have not yet performed a secular approximation, in contrast to the derivation of equation (20). However, naively applying a secular approximation that neglects off-diagonal terms for which $\zeta \ne \zeta ^{\prime}$ in the summand in equation (33) would not be appropriate, because this would eliminate terms that are resonant in the limit $C\to 0$. Instead we perform a partial secular approximation which eliminates off-diagonal terms for which ζ and ζ' have opposite sign. These terms remain far off-resonance for all values of C. The resulting master equation is given by

$$\begin{eqnarray}\dot{\rho } & = & -{\rm{i}}[{H}_{d},\rho ]\\ & & +\displaystyle \sum _{\zeta ,\zeta ^{\prime} ={\omega }_{\mathrm{1,2}}}\displaystyle \sum _{\mu ,\nu =1}^{2}[{{\rm{\Gamma }}}_{\mu \nu }(\zeta )({A}_{\nu \zeta }\rho {A}_{\mu \zeta ^{\prime} }^{\dagger }-{A}_{\mu \zeta ^{\prime} }^{\dagger }{A}_{\nu \zeta }\rho )+{{\rm{\Gamma }}}_{\mu \nu }(-\zeta )({A}_{\nu \zeta }^{\dagger }\rho {A}_{\mu \zeta ^{\prime} }-{A}_{\mu \zeta ^{\prime} }{A}_{\nu \zeta }^{\dagger }\rho )+{\rm{h}}.{\rm{c}}.].\end{eqnarray} \tag{ 40 }$$

We are now in a position to compare our master equation (40) with the usual result in (20). In the limit $C\to 0$ we have $\eta \to {\omega }_{0}$ so that ${\omega }_{\mathrm{1,2}}\to {\omega }_{0}$. The rates and shifts in equation (38) are then evaluated at ω₀ within equation (40). Also, the Hamiltonian H_d tends to the bare Hamiltonian ${\omega }_{0}({\sigma }_{1}^{+}{\sigma }_{1}^{-}+{\sigma }_{2}^{+}{\sigma }_{2}^{-})$, and furthermore we have that

$$\begin{eqnarray}&&\displaystyle \sum _{\zeta ={\omega }_{\mathrm{1,2}}}{A}_{\mu \zeta }\to {\sigma }_{\mu }^{-}.\end{eqnarray} \tag{ 41 }$$

Thus, taking the limit $C\to 0$ in equation (40) one recovers equation (20) with Δ₁₂ replaced by ${\tilde{{\rm{\Delta }}}}_{12}$ given in equation (28). However, since ${\tilde{{\rm{\Delta }}}}_{12}\to {{\rm{\Delta }}}_{12}$ when $C\to 0$, equations (40) and (20) coincide in this limit. For finite C, equation (40) offers separation-dependent corrections to the usual master equation and is the main result of this section.

It is important to note that while our master equation (40) is generally different to the usual gauge-invariant result (equation (20)) there is no cause for concern regarding the issue of gauge-invariance. As we have shown the standard master equation can be obtained when α_k = 0 provided one uses a partitioning of the Hamiltonian in the form $H={H}_{0}+{V}^{\mathrm{usual}}$ where H₀ is given by equation (15) and

$$\begin{eqnarray}&&{V}^{\mathrm{usual}}=\displaystyle \sum _{\mu =1}^{2}{\omega }_{0}{\sigma }_{\mu }^{y}{\bf{d}}\cdot {\bf{A}}({{\bf{R}}}_{\mu })+\displaystyle \sum _{\mu =1}^{2}\displaystyle \frac{{e}^{2}}{2m}{\bf{A}}{({{\bf{R}}}_{\mu })}^{2}+{C}_{\{0\}}.\end{eqnarray} \tag{ 42 }$$

Our master equation (40) has also been obtained by choosing α_k = 0, but our derivation makes use of the different partitioning $H={\tilde{H}}_{0}+V$ where ${\tilde{H}}_{0}$ is defined as in equation (29) and $V={V}^{\mathrm{usual}}-{C}_{\{0\}}$ is defined as in equation (26). The two different partitionings of the same Hamiltonian yield two different second order master equations.

As we have shown the standard Born–Markov-secular master equation (20) can be obtained for any other choice of α_k provided that the unperturbed Hamiltonian is always defined as in equation (15). Similarly a full secular approximation of our master equation (40) can also be obtained for any other choice of α_k provided the unperturbed Hamiltonian is always defined as in equation (29). We note further that the secular approximation is well justified within the near-field regime of interest $R\ll {\lambda }_{0}$. Let us consider for example the multipolar gauge obtained by choosing α_k = 1. In order to achieve the appropriate partitioning of the multipolar Hamiltonian for derivation of our master equation one must add C_{0} to the usual multipolar unperturbed Hamiltonian H₀ given in equation (15), and simultaneously subtract C_{0} from the usual multipolar interaction Hamiltonian. Using this repartitioning of the multipolar Hamiltonian the Born–Markov-secular master equation is found to coincide with our master equation (40) once a full secular approximation is performed within the latter. This derivation is however, more cumbersome than the Coulomb gauge derivation. Since the Coulomb energy is naturally explicit within the Coulomb gauge, the latter is the most natural gauge to choose for the purpose of including the relatively strong static interaction within the unperturbed Hamiltonian.

Any difference between the master equation (40) and the corresponding partially-secular result found using ${\alpha }_{k}\ne 0$ is contained entirely within non-secular contributions. These contributions are negligible within the regime of interest R ≪ λ₀ and have only been retained within equation (40) to facilitate comparison with the standard result equation (20). Moreover, in the far-field regime $R\gg {\lambda }_{0}$ the master equations (20) and (40) coincide, so the master equation (40) is also gauge-invariant within this regime.

For large inter-dipole separations R ≫ λ the master equations (20) and (40) coincide and they therefore yield identical physical predictions. However, in the near-zone R ≪ λ₀ the master equations generally exhibit significant differences. To compare the two sets of predictions we assume a vacuum field N = 0 and consider the experimental situation in which the system is prepared in the symmetric state $| {\epsilon }_{3}\rangle$. This state is a simultaneous eigenstate of the dipole Hamiltonian ${\omega }_{0}({\sigma }_{1}^{+}{\sigma }_{1}^{-}+{\sigma }_{2}^{+}{\sigma }_{2}^{-})$ appearing in the standard master equation (20), and of the Hamiltonian H_d appearing in our master equation (40). Experimentally, one expects to find that the system initially prepared in the state $| {\epsilon }_{3}\rangle$ decays into the stationary state. Theoretically, different stationary states are predicted by the two master equations (20) and (40), and the rates of decay into these respective stationary states are also different. Figures 1 and 2 compare the symmetric and stationary state populations found using master equations (20) and (40) when the system starts in the symmetric eigenstate $| {\epsilon }_{3}\rangle$. For small separations the ground and symmetric state populations obtained from our master equation (40) crossover earlier, which indicates more rapid symmetric state decay than is predicted by equation (20) (see figure 1). This gives rise to the different starting values at R = r_a of the curves depicted in figure 2. For larger separations the solutions converge and become indistinguishable for all times.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The different behaviour in figures 1 and 2 can be understood by looking at a few relevant quantities. The matrix element of the combined dipole moment between ground and symmetric eigenstates is found to be

$$\begin{eqnarray}&&{{\bf{d}}}_{31}=\langle {\epsilon }_{3}| {{\bf{d}}}_{1}+{{\bf{d}}}_{2}| {\epsilon }_{1}\rangle =2a{\bf{d}},\end{eqnarray} \tag{ 43 }$$

which is different to the usual transition dipole moment $\langle {\epsilon }_{3}| {{\bf{d}}}_{1}+{{\bf{d}}}_{2}| {gg}\rangle =\sqrt{2}{\bf{d}}$. Since $a\to 1/\sqrt{2}$ as $C\to 0$, the dipole moment d₃₁ reduces to $\sqrt{2}{\bf{d}}$ when $R\to \infty$. As R decreases, however, d₁₃ becomes increasingly large compared with $\sqrt{2}{\bf{d}}$. This is consistent with the more rapid decay observed in figure 1. A more complete explanation of this behaviour can be obtained by calculating the rate of decay of the symmetric state into the vacuum, which we denote γ_s. Using Fermi's golden rule, and the eigenstates given in equation (32), we obtain

$$\begin{eqnarray}&&{\gamma }_{s}({\omega }_{2})=2{c}^{2}{[{\gamma }_{\mu \mu }({\omega }_{2})+{\gamma }_{12}({\omega }_{2})]}_{N=0}.\end{eqnarray} \tag{ 44 }$$

Only when $C\to 0$, such that ${\omega }_{2}=\eta +C\to {\omega }_{0}$ and $c\to 1/\sqrt{2}$, does this decay rate reduce to that obtained when using the bare eigenstates $| i,j\rangle ,(i,j=e,g)$, which is

$$\begin{eqnarray}&&{\gamma }_{s,0}=\gamma +{\gamma }_{12}({\omega }_{0}),\end{eqnarray} \tag{ 45 }$$

where γ₁₂(ω₀) is given in equation (21). As shown in figure 3, for sufficiently small R the decay rate γ_s(ω₂) is significantly larger than γ_s,0.

Zoom In Zoom Out Reset image size

In contrast to the decay behaviour of the symmetric state, the predictions of the master equations (20) and (40) are the same if the system is assumed to be prepared in the anti-symmetric state $| {\epsilon }_{2}\rangle$, which like $| {\epsilon }_{3}\rangle$ is a simultaneous eigenstate of ${\omega }_{0}({\sigma }_{1}^{+}{\sigma }_{1}^{-}+{\sigma }_{2}^{+}{\sigma }_{2}^{-})$ and H_d. Both master equations predict that the population of the state $| {\epsilon }_{2}\rangle$ remains stationary, i.e., that it is a completely dark state. This can be understood by noting that the collective dipole moment associated with the anti-symmetric to stationary state transition vanishes when either stationary state, $| g,g\rangle$ or $| {\epsilon }_{1}\rangle$, is used. Finally, the predicted behaviour by our master equation (40) of the standard stationary state $| g,g\rangle$ is illustrated in figure 4. For an initial state $| {\epsilon }_{3}\rangle$ the population p_gg(t) of the state $| g,g\rangle$ at a given time t, is identical to that predicted by equation (20) only for sufficiently large R whereby $| {\epsilon }_{1}\rangle \approx | g,g\rangle$.

Zoom In Zoom Out Reset image size

In this section we apply our master equation (40) to calculate the emission spectrum of the two-dipole system initially prepared in the symmetric state $| {\epsilon }_{3}\rangle$. This provides a means by which to test experimentally whether our predictions are closer to measured values than the standard approach. The spectrum of radiation is defined according to the quantum theory of photodetection by [39, 40]

$$\begin{eqnarray}&&s(\omega )={\int }_{0}^{\infty }{\rm{d}}{t}{\int }_{0}^{\infty }{\rm{d}}{t}^{\prime} {{\rm{e}}}^{-{\rm{i}}\omega (t-t^{\prime} )}\langle {{\bf{E}}}_{s,{\rm{rad}}}^{(-)}(t,{\bf{x}})\cdot {{\bf{E}}}_{s,{\rm{rad}}}^{(+)}(t^{\prime} ,{\bf{x}}){\rangle }_{0},\end{eqnarray} \tag{ 46 }$$

where for simplicity we assume that the field is in the vacuum state. Since the master equations (20) and (40) yield different predictions for this experimentally measurable quantity, an experiment could be used to test which master equation is the most accurate.

The detector is located at position x with $x\gg R$, so that only the radiative component ${{\bf{E}}}_{s,\mathrm{rad}}$ of the electric source field, which varies as $| {\bf{x}}-{{\bf{R}}}_{\mu }{| }^{-1}$, need be used. This is the only part of the field responsible for irreversibly carrying energy away from the sources. The positive and negative frequency components of the radiation source field are given within both rotating-wave and Markov approximations by

$$\begin{eqnarray}&&{E}_{s,\mathrm{rad},i}^{(\pm )}(t,{\bf{x}})=\displaystyle \sum _{\mu =1}^{2}\displaystyle \frac{{\omega }_{0}^{2}}{4\pi {r}_{\mu }}({\delta }_{{ij}}-{\hat{r}}_{\mu ,i}{\hat{r}}_{\mu ,j}){d}_{j}{\sigma }^{\mp }({t}_{\mu }),\end{eqnarray} \tag{ 47 }$$

where ${{\bf{r}}}_{\mu }={\bf{x}}-{{\bf{R}}}_{\mu }$, and ${t}_{\mu }=t-{r}_{\mu }$ is the retarded time associated with the μ'th source. For $x\gg R$ we have to a very good approximation that ${{\bf{r}}}_{1}={{\bf{r}}}_{2}={\bf{r}}$, where r is the relative vector from x to the midpoint of R₁ and R₂. Substituting equation (47) into equation (46) within this approximation yields

$$\begin{eqnarray}&&s(\omega )=\mu {\omega }_{0}^{4}{\int }_{0}^{\infty }{\rm{d}}{t}{\int }_{0}^{\infty }{\rm{d}}{t}^{\prime} {{\rm{e}}}^{-{\rm{i}}\omega (t-t^{\prime} )}\displaystyle \sum _{\mu ,\nu =1}^{2}\langle {\sigma }_{\mu }^{+}(t){\sigma }_{\nu }^{-}(t^{\prime} )\rangle ,\end{eqnarray} \tag{ 48 }$$

where

$$\begin{eqnarray}&&\mu ={\left[\displaystyle \frac{1}{4\pi r}({\delta }_{{ij}}-{\hat{r}}_{i}{\hat{r}}_{j}){d}_{j}\right]}^{2}.\end{eqnarray} \tag{ 49 }$$

To begin with, let us use the standard master equation (20) to find the required two-time correlation function. Assuming that the system is initially prepared in the symmetric state $| {\epsilon }_{3}\rangle$, using the standard master equation (20) together with the method of calculation given in appendix A.3 we obtain the two-time correlation function

$$\begin{eqnarray}&&\displaystyle \sum _{\mu ,\nu =1}^{2}\langle {\sigma }_{\mu }^{+}(t){\sigma }_{\nu }^{-}(t^{\prime} )\rangle =2{{\rm{e}}}^{-\displaystyle \frac{{\gamma }_{s,0}}{2}(t+t^{\prime} )}{{\rm{e}}}^{{\rm{i}}({\tilde{\omega }}_{0}+{{\rm{\Delta }}}_{12})(t-t^{\prime} )},\end{eqnarray} \tag{ 50 }$$

where γ_s,0 and Δ₁₂ are given in equations (45) and (21), respectively. By direct integration of equation (50) one obtains the corresponding Lorentzian spectrum

$$\begin{eqnarray}&&{s}_{0}(\omega )=\displaystyle \frac{2{\omega }_{0}^{4}\mu }{{({\gamma }_{s,0}/2)}^{2}+{(\omega -[{\tilde{\omega }}_{0}+{{\rm{\Delta }}}_{12}])}^{2}}.\end{eqnarray} \tag{ 51 }$$

Let us now turn our attention to the spectrum obtained from our new master equation (40). We have seen that the solutions of equations (40) and (20) differ only in the near field regime $R\ll \lambda$. For sufficiently small R we have that ${\omega }_{0}\sim C$, and the frequency difference ${\omega }_{2}-{\omega }_{1}=2C\sim 2{\omega }_{0}$ is large. In this situation we can perform a full secular approximation within equation (40) to obtain the master equation

$$\begin{eqnarray}&&\dot{\rho }=-{\rm{i}}[{\tilde{H}}_{d},\rho ]+{ \mathcal D }(\rho ),\end{eqnarray} \tag{ 52 }$$

with

$$\begin{eqnarray}&&{\tilde{H}}_{d}={H}_{d}+\displaystyle \sum _{\omega =\pm {\omega }_{\mathrm{1,2}}}\displaystyle \sum _{\mu ,\nu =1}^{2}{S}_{\mu \nu }(\omega ){A}_{\mu \omega }^{\dagger }{A}_{\nu \omega },\end{eqnarray} \tag{ 53 }$$

and

$$\begin{eqnarray}&&{ \mathcal D }(\rho )=\displaystyle \sum _{\omega =\pm {\omega }_{\mathrm{1,2}}}\displaystyle \sum _{\mu ,\nu =1}^{2}{\gamma }_{\mu \nu }(\omega )\left[{A}_{\nu \omega }\rho {A}_{\mu \omega }^{\dagger }-\displaystyle \frac{1}{2}\{{A}_{\mu \omega }^{\dagger }{A}_{\nu \omega },\rho \}\right].\end{eqnarray} \tag{ 54 }$$

Solving this secular master equation allows us to obtain a simple expression for the emission spectrum.

The correlation function in equation (46) defines the radiation intensity when it is evaluated at t = t'. Naively calculating the quantity ${\sum }_{\mu ,\nu =1}^{2}\langle {\sigma }_{\mu }^{+}(t){\sigma }_{\nu }^{-}(t){\rangle }_{1}$ taken in the stationary state $| {\epsilon }_{1}\rangle$ yields a non-zero stationary intensity, because $| {\epsilon }_{1}\rangle$ is a superposition involving both $| g,g\rangle$ and the doubly excited bare state $| e,e\rangle$. A non-zero radiation intensity even in the stationary state is clearly non-physical. However, a more careful analysis recognises that when the radiation source fields are to be used in conjunction with equation (40) the optical approximations used in their derivation should be applied in the interaction picture defined in terms of the dressed Hamiltonian H_d given in equation (30). One then obtains the source field

$$\begin{eqnarray}&&{E}_{s,\mathrm{rad},i}^{(+)}(t,{\bf{x}})=\displaystyle \sum _{\mu =1}^{2}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{nm}}{n\lt m}}\displaystyle \frac{{\epsilon }_{{nm}}^{2}}{4\pi {r}_{\mu }}({\delta }_{{ij}}-{\hat{r}}_{\mu ,i}{\hat{r}}_{\mu ,j}){d}_{j}{\sigma }_{\mu ,{nm}}{\theta }_{{nm}}({t}_{\mu }),\end{eqnarray} \tag{ 55 }$$

where ${\sigma }_{\mu ,{nm}}={\sigma }_{\mu ,{nm}}^{+}+{\sigma }_{\mu ,{nm}}^{-}$ in which ${\sigma }_{\mu ,{nm}}^{\pm }$ denotes the nm'th matrix element of ${\sigma }_{\mu }^{\pm }$ in the basis $| {\epsilon }_{n}\rangle$. The transition frequencies associated with this basis are denoted ${\epsilon }_{{nm}}={\epsilon }_{n}-{\epsilon }_{m}$, and the raising and lowering operators are denoted ${\theta }_{{nm}}=| {\epsilon }_{n}\rangle \langle {\epsilon }_{m}| ,\,n\ne m$. The derivation of equation (55) is given in appendix A.4. According to equation (55) the annihilation (creation) radiation source field is now associated with lowering (raising) operators in the dressed basis $| {\epsilon }_{i}\rangle$ rather than in the bare basis $| n,m\rangle ,(n,m=e,g)$. Substitution of equation (55) into equation (46) yields

$$\begin{eqnarray}&&s(\omega )=\mu {\int }_{0}^{\infty }{\rm{d}}{t}{\int }_{0}^{\infty }{\rm{d}}{t}^{\prime} \displaystyle \sum _{\mu ,\nu =1}^{2}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{nm}}{n\lt m}}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{pq}}{q\lt p}}{\epsilon }_{{pq}}^{2}{\epsilon }_{{nm}}^{2}{\sigma }_{\mu ,{pq}}{\sigma }_{\nu ,{nm}}\langle {\theta }_{{pq}}(t){\theta }_{{nm}}(t^{\prime} )\rangle ,\end{eqnarray} \tag{ 56 }$$

where we have again assumed that ${{\bf{r}}}_{1}={{\bf{r}}}_{2}={\bf{r}}$. Unlike the correlation function in equation (92), when t = t' the correlation functions $\langle {\theta }_{{pq}}(t){\theta }_{{nm}}(t^{\prime} )\rangle ,\,p\gt q,\,m\gt n$ vanish in the stationary (ground) state θ₁₁. The radiation intensity is therefore seen to vanish in the stationary limit as required physically.

Taken in the symmetric state θ₃₃ the only non-zero two-time correlation function that contributes to equation (56) is found to be

$$\begin{eqnarray}&&{C}_{33}(t,t^{\prime} )=\langle {\theta }_{31}(t){\theta }_{13}(t^{\prime} )\rangle ={{\rm{e}}}^{-{c}^{2}[{\gamma }_{\mu \mu }({\omega }_{2})+{\gamma }_{12}({\omega }_{2})](t+t^{\prime} )}{{\rm{e}}}^{{\rm{i}}{\tilde{\omega }}_{2}(t-t^{\prime} )},\end{eqnarray} \tag{ 57 }$$

where

$$\begin{eqnarray}{\tilde{\omega }}_{2} & = & {\omega }_{2}+2({S}_{\mu \mu }(-{\omega }_{1})[{b}^{2}-{a}^{2}]+{S}_{12}(-{\omega }_{1})[{b}^{2}+{a}^{2}]\\ & & +{c}^{2}[{S}_{\mu \mu }({\omega }_{2})-{S}_{\mu \mu }(-{\omega }_{2})+{S}_{12}({\omega }_{2})-{S}_{12}(-{\omega }_{2})])\end{eqnarray} \tag{ 58 }$$

is the shifted symmetric to ground transition frequency. Integration of equation (57) according to equation (56) then yields the Lorentzian spectrum

$$\begin{eqnarray}&&s(\omega )=\displaystyle \frac{{(2a{\omega }_{2}^{2})}^{2}\mu }{{({\gamma }_{s}({\omega }_{2})/2)}^{2}+{(\omega -{\tilde{\omega }}_{2})}^{2}}.\end{eqnarray} \tag{ 59 }$$

Full details of the calculation of the spectrum in equation (59) are given in appendix A.4. In the limit of large separation $C\to 0$, which implies that ${\omega }_{2}\to {\omega }_{0}$, ${\tilde{\omega }}_{2}\to {\tilde{\omega }}_{0}+{{\rm{\Delta }}}_{12}$, ${\gamma }_{s}({\omega }_{2})\to {\gamma }_{s,0}$, and $a\to 1/\sqrt{2}$. As a result $s(\omega )\to {s}_{0}(\omega )$ for large R and the predicted spectra coincide. On the other hand, for sufficiently small R the spectrum s(ω) again offers separation-dependent corrections to the standard result s₀(ω).

The two spectra ${s}_{0}(\omega )$ and s(ω) are compared in figures 5 and 6. As their relative widths are proportional to the rates, they are given in equations (45) and (44), respectively. These quantities have been plotted already in figure 3. The relative heights of the spectral peaks are ${s}_{0}({\tilde{\omega }}_{0}+{{\rm{\Delta }}}_{12})/\mu$ and $s({\tilde{\omega }}_{2})/\mu$ respectively, which are plotted in figure 7. This figure shows that the peak heights in the spectra begin to diverge as R decreases. At a separation of 15r_a, where r_a = n²a₀, n = 50 is a characteristic Rydberg atomic radius, the peak value of s(ω) is around two times larger than the peak value of s₀(ω) for the parameters chosen here. The positions of the peaks are ${\tilde{\omega }}_{0}+{{\rm{\Delta }}}_{12}$ and ${\tilde{\omega }}_{2}$, respectively, and these are plotted in figure 8. The ultra-violet cut-off chosen for the calculation of the single-dipole shift components corresponds to the inverse dipole radius wavelength, namely 2π c /r_a. This value is chosen for consistency with the electric dipole approximation that we have used throughout. For small R the spectrum s(ω) is blue-shifted relative to s₀(ω). Figure 8 shows that the ratio of peak positions approaches a constant value around two for very small R. These differences could in principle be detected in an experiment. At a separation of 20r_a, which is roughly 2.5 μm, for instance, the difference in shifted frequencies ${\tilde{\omega }}_{2}-{\tilde{\omega }}_{0}-{{\rm{\Delta }}}_{12}$ is around 1 Ghz for the parameters chosen in figure 1. This is similar in magnitude to the Lamb-shift in atomic hydrogen.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here we determine the contribution of self-energy terms to dipole level-shifts and demonstrate that the single-dipole transition shift is gauge-invariant. The self-energy term V⁽²⁾ is given in equation (4). The shifts arising from this term are divergent in the mode continuum limit ${\omega }_{k}\to \infty$, but this divergence is not unexpected within the non-relativistic dipole approximated treatment. It is typically handled through the introduction of an ultra-violet cut-off. In the treatment of the Lamb-shift in atomic hydrogen the Coulomb gauge self-energy V⁽²⁾ with α_k = 0 is independent of the atomic electron levels and is therefore ignored within the calculation of the measurable shift [35]. In the multipolar gauge V⁽²⁾ represents a polarisation self-energy term and when its contribution is combined with the remaining contribution to the shift coming from the linear part of the multipolar interaction Hamiltonian one obtains the same result as the Coulomb gauge treatment. In all cases mass renormalisation must also be performed to obtain the correct shift.

In the Coulomb gauge V⁽²⁾ does not contribute to the master equation transition shift of the two-level dipole, which is the difference between excited and ground state shifts. This is independent of whether the two-level approximation has been made. However, even within the Coulomb gauge it is important to note that one must generally account for all self-energy contributions when explicitly verifying that quantities are gauge-invariant. In particular, to verify that the ground and excited level-shifts are separately gauge-invariant, the contributions ${\langle n| }^{(2)}V| n\rangle ,\,n\,=\,e,g$ must be taken into account.

Using the Hamiltonian in equation (3), the standard Born-Markov master equation has form given in equation (11), in which the decay rate γ is independent of α_k. In this equation, the transition shift can be expressed as the difference between excited and ground state shifts as ${\tilde{\omega }}_{0}-{\omega }_{0}={\rm{\Delta }}={\delta }_{e}^{(1)}-{\delta }_{g}^{(1)}$ where

$$\begin{eqnarray}{\delta }_{e}^{(1)} & = & \displaystyle \int {{\rm{d}}}^{3}k\displaystyle \sum _{\lambda }\displaystyle \frac{| {{\bf{e}}}_{\lambda }({\bf{k}})\cdot {\bf{d}}{| }^{2}}{2{(2\pi )}^{3}}{\omega }_{0}\left[\displaystyle \frac{{{u}_{k}^{+}}^{2}{N}_{k}}{{\omega }_{k}+{\omega }_{0}}-\displaystyle \frac{{{u}_{k}^{-}}^{2}[1+{N}_{k}]}{{\omega }_{k}-{\omega }_{0}}\right],\\ {\delta }_{g}^{(1)} & = & \displaystyle \int {{\rm{d}}}^{3}k\displaystyle \sum _{\lambda }\displaystyle \frac{| {{\bf{e}}}_{\lambda }({\bf{k}})\cdot {\bf{d}}{| }^{2}}{2{(2\pi )}^{3}}{\omega }_{0}\left[\displaystyle \frac{{{u}_{k}^{-}}^{2}{N}_{k}}{{\omega }_{k}-{\omega }_{0}}-\displaystyle \frac{{{u}_{k}^{+}}^{2}[1+{N}_{k}]}{{\omega }_{k}+{\omega }_{0}}\right],\end{eqnarray} \tag{ 60 }$$

are α_k-dependent. This α_k-dependence is due to the lack of any contribution from the self-energy term V⁽²⁾ in equation (60).

The α_k-dependence within the master equation is eliminated when one accounts for the self-energy contributions and the effect of the two-level approximation, recalling that the latter was made after the transformation ${R}_{\{{\alpha }_{k}\}}$ was performed. More specifically it is possible to demonstrate that the single dipole master equation (10) is α_k-independent, and that it coincides with equation (11). First we note that we can continue to express the second line in equation (10) in terms of the original partition H = H₀ + V. Thus, provided H₀ is kept the same for each choice of the α_k the dissipative part of the master equation is α_k-independent.

It remains to show that when one adds the shift contributions ${\delta }_{e,g}^{(1)}$ coming from the second line in equation (10) to the corresponding self-energy contribution in equation (9) one obtains gauge-invariant total shifts. To this end let us first consider the Coulomb gauge α_k = 0. The total excited and ground state shifts are

$$\begin{eqnarray}&&{\delta }_{e,\mathrm{CG}}={\delta }_{e,\mathrm{CG}}^{(1)}+{\delta }_{\mathrm{CG}}^{(2)},\qquad {\delta }_{g,\mathrm{CG}}={\delta }_{g,\mathrm{CG}}^{(1)}+{\delta }_{\mathrm{CG}}^{(2)}.\end{eqnarray} \tag{ 61 }$$

The components ${\delta }_{e,\mathrm{CG}}^{(1)}$ are obtained by setting α_k = 0 in equation (60), while the remaining component

$$\begin{eqnarray}&&{\delta }_{\mathrm{CG}}^{(2)}=\displaystyle \frac{{e}^{2}}{2m}\int {{\rm{d}}}^{3}k\displaystyle \sum _{\lambda }\displaystyle \frac{| {{\bf{e}}}_{\lambda }({\bf{k}}){| }^{2}}{2{(2\pi )}^{3}{\omega }_{k}}(1+2{N}_{k})\end{eqnarray} \tag{ 62 }$$

is the Coulomb gauge self-energy shift due to the ${{\bf{A}}}_{{\rm{T}}}^{2}$ part of the Coulomb gauge interaction Hamiltonian. Since this term is independent of the dipole, the shift ${\delta }_{\mathrm{CG}}^{(2)}$ is the same for the ground and excited levels. The single-dipole transition shift Δ given in equation (12) in the main text can be expressed in terms of Coulomb gauge shifts as

$$\begin{eqnarray}&&{\rm{\Delta }}={\delta }_{e,\mathrm{CG}}-{\delta }_{g,\mathrm{CG}}={\delta }_{e,\mathrm{CG}}^{(1)}-{\delta }_{g,\mathrm{CG}}^{(1)}.\end{eqnarray} \tag{ 63 }$$

More generally, for arbitrary α_k the total ground and excited state level shifts are denoted δ_{e, g}. In what follows we will show that

$$\begin{eqnarray}&&{\delta }_{e}-{\delta }_{\mathrm{CG}}^{(2)}={\delta }_{e,\mathrm{CG}}^{(1)}\end{eqnarray} \tag{ 64a }$$

and

$$\begin{eqnarray}&&{\delta }_{g}-{\delta }_{\mathrm{CG}}^{(2)}={\delta }_{g,\mathrm{CG}}^{(1)},\end{eqnarray} \tag{ 64b }$$

from which it follows using equation (63) that δ_e − δ_g = Δ for all choices of α_k.

In order to show that equations (64a) and (64b) hold we must carefully account for the two-level approximation, which was performed after the gauge transformation ${R}_{\{{\alpha }_{k}\}}$. Let us consider a general shift of the m'th level of the dipole with the form

$$\begin{eqnarray}&&{\tilde{\omega }}_{m}={\omega }_{m}+\displaystyle \sum _{n}{\omega }_{{nm}}| {\bf{v}}\cdot {{\bf{d}}}_{{nm}}{| }^{2},\end{eqnarray} \tag{ 65 }$$

where v is arbitrary. If we restrict ourselves to two levels e and g, and if m = e in the above, then the sum includes only one other level n = g, so we get for the shift

$$\begin{eqnarray}&&\displaystyle \sum _{n}{\omega }_{{ne}}| {\bf{v}}\cdot {{\bf{d}}}_{{ne}}{| }^{2}=-{\omega }_{0}| {\bf{v}}\cdot {\bf{d}}{| }^{2},\end{eqnarray} \tag{ 66 }$$

where ${\omega }_{0}:= {\omega }_{{eg}}=-{\omega }_{{ge}}$ and ${\bf{d}}:= {{\bf{d}}}_{{eg}}={{\bf{d}}}_{{ge}}^{* }$. If instead m = g then the shift is

$$\begin{eqnarray}&&\displaystyle \sum _{n}{\omega }_{{ng}}| {\bf{v}}\cdot {{\bf{d}}}_{{ng}}{| }^{2}=+{\omega }_{0}| {\bf{v}}\cdot {\bf{d}}{| }^{2}.\end{eqnarray} \tag{ 67 }$$

The shift is clearly different in the m = e and m = g cases when considering a two-level system. However, for an infinite-dimensional dipole the shift is independent of m being given by

$$\begin{eqnarray}&&\displaystyle \sum _{n}{\omega }_{{nm}}| {\bf{v}}\cdot {{\bf{d}}}_{{nm}}{| }^{2}=\displaystyle \frac{{e}^{2}}{2m}| {\bf{v}}{| }^{2},\end{eqnarray} \tag{ 68 }$$

where we have made use of the identity

$$\begin{eqnarray}&&\displaystyle \sum _{n}{\omega }_{{nm}}{d}_{{nm}}^{i}{d}_{{mn}}^{j}={\rm{i}}\displaystyle \frac{{e}^{2}}{2m}\langle m| [{p}_{i},{r}_{j}]| m\rangle ={\delta }_{{ij}}\displaystyle \frac{{e}^{2}}{2m}.\end{eqnarray} \tag{ 69 }$$

The difference between the finite and infinite-dimensional cases arises because the proof of equation (69) rests directly on the CCR algebra $[{r}_{i},{p}_{j}]={\rm{i}}{\delta }_{{ij}}$, which can only be supported in infinite-dimensions. When the algebra is truncated to su(2), the same shift comes out level-dependent. Since the gauge transformation ${R}_{\{{\alpha }_{k}\}}$ is made on the infinite-dimensional dipole it is necessary to employ equation (68) in order to exhibit gauge-invariance of the shifts. Thus, in order to get the correct level-shifts within the two-level approximation, when dealing with the excited level shift m = e we use equations (66) and (68), which imply

$$\begin{eqnarray}&&{\omega }_{0}| {\bf{v}}\cdot {\bf{d}}{| }^{2}=-\displaystyle \frac{{e}^{2}}{2m}| {\bf{v}}{| }^{2},\end{eqnarray} \tag{ 70 }$$

but when dealing with the ground level shift m = g we use equations (67) and (68), which imply

$$\begin{eqnarray}&&{\omega }_{0}| {\bf{v}}\cdot {\bf{d}}{| }^{2}=\displaystyle \frac{{e}^{2}}{2m}| {\bf{v}}{| }^{2}.\end{eqnarray} \tag{ 71 }$$

We now proceed to verify that equations (64a) and (64b) hold. The complete shifts δ_{e, g} are obtained by taking the shifts in equation (60) and adding their respective self-energy contributions. Subtracting ${\delta }_{\mathrm{CG}}^{(2)}$ in equation (62) from δ_e and subsequently using equation (70), which is appropriate for the excited state shift, we obtain

$$\begin{eqnarray}{\delta }_{e}-{\delta }_{\mathrm{CG}}^{(2)} & = & \displaystyle \int {{\rm{d}}}^{3}k\displaystyle \sum _{\lambda }\displaystyle \frac{| {{\bf{e}}}_{\lambda }({\bf{k}})\cdot {\bf{d}}{| }^{2}}{2{(2\pi )}^{3}}\\ & & \times \left({\alpha }_{k}^{2}-{\alpha }_{k}({\alpha }_{k}-2)[1+2{N}_{k}]\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}+{\omega }_{0}\left[\displaystyle \frac{{{u}_{k}^{+}}^{2}{N}_{k}}{{\omega }_{k}+{\omega }_{0}}-\displaystyle \frac{{{u}_{k}^{-}}^{2}[1+{N}_{k}]}{{\omega }_{k}-{\omega }_{0}}\right]\right).\end{eqnarray} \tag{ 72 }$$

Using equation (6) we express the bracket within the integrand in this expression in terms of α_k. The part independent of N_k is

$$\begin{eqnarray}&&{\alpha }_{k}^{2}-{\alpha }_{k}({\alpha }_{k}-2)\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}-\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}-{\omega }_{0}}\left[{(1-{\alpha }_{k})}^{2}\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}+{\alpha }_{k}^{2}\displaystyle \frac{{\omega }_{k}}{{\omega }_{0}}+2{\alpha }_{k}(1-{\alpha }_{k})\right].\end{eqnarray} \tag{ 73 }$$

In this expression we identify the coefficient of ${\alpha }_{k}^{2}$ as

$$\begin{eqnarray}&&1-\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}-{\omega }_{0}}\left(\displaystyle \frac{{\omega }_{k}}{{\omega }_{0}}+\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}-2\right)-\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}=1-\displaystyle \frac{{\omega }_{k}-{\omega }_{0}}{{\omega }_{k}}-\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}=0,\end{eqnarray} \tag{ 74 }$$

and the coefficient of 2α_k as

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}-\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}-{\omega }_{0}}\left(1-\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}\right)=\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}-\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}-{\omega }_{0}}\displaystyle \frac{{\omega }_{k}-{\omega }_{0}}{{\omega }_{k}}=0.\end{eqnarray} \tag{ 75 }$$

Thus, equation (73) is α_k-independent. The remaining part is

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{0}^{2}}{{\omega }_{k}({\omega }_{0}-{\omega }_{k})}.\end{eqnarray} \tag{ 76 }$$

The N_k-dependent parts of ${\delta }_{e}-{\delta }_{\mathrm{CG}}^{(2)}$ can be dealt with in a similar manner. The coefficient of ${\alpha }_{k}^{2}$ in the N_k-dependent part of the bracket within the integrand of the expression for ${\delta }_{e}-{\delta }_{\mathrm{CG}}^{(2)}$ is

$$\begin{eqnarray}&&-2\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}+\left(\displaystyle \frac{{\omega }_{0}^{2}}{{\omega }_{k}}+{\omega }_{k}\right)\left(\displaystyle \frac{1}{{\omega }_{0}+{\omega }_{k}}-\displaystyle \frac{1}{{\omega }_{k}-{\omega }_{0}}\right)+{\omega }_{0}\left(\displaystyle \frac{1}{{\omega }_{0}+{\omega }_{k}}+\displaystyle \frac{1}{{\omega }_{k}-{\omega }_{0}}\right)\\ &&=\,2\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}\left[-1+\displaystyle \frac{1}{{\omega }_{k}^{2}-{\omega }_{0}^{2}}(-{\omega }_{0}^{2}-{\omega }_{k}^{2}+2{\omega }_{k}^{2})\right]=0.\end{eqnarray} \tag{ 77 }$$

Similarly, the coefficient of α_k is

$$\begin{eqnarray}&&4\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}+2\displaystyle \frac{{\omega }_{0}^{2}}{{\omega }_{k}}\left(\displaystyle \frac{1}{{\omega }_{k}-{\omega }_{0}}-\displaystyle \frac{1}{{\omega }_{k}+{\omega }_{0}}\right)-2\left(\displaystyle \frac{1}{{\omega }_{k}-{\omega }_{0}}-\displaystyle \frac{1}{{\omega }_{k}+{\omega }_{0}}\right)\\ &&\quad =\,4\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}\left[1-\displaystyle \frac{1}{({\omega }_{0}+{\omega }_{k})({\omega }_{k}-{\omega }_{0})}({\omega }_{k}^{2}-{\omega }_{0}^{2})\right]=0.\end{eqnarray} \tag{ 78 }$$

The remaining N_k-dependent part is

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{0}^{2}}{{\omega }_{k}}\left(\displaystyle \frac{1}{{\omega }_{k}+{\omega }_{0}}+\displaystyle \frac{1}{{\omega }_{0}-{\omega }_{k}}\right).\end{eqnarray} \tag{ 79 }$$

Combining equations (73) and (79) we obtain the α_k-independent result

$$\begin{eqnarray}&&{\delta }_{e}-{\delta }_{\mathrm{CG}}^{(2)}=\int {{\rm{d}}}^{3}k\displaystyle \sum _{\lambda }\displaystyle \frac{| {{\bf{e}}}_{\lambda }({\bf{k}})\cdot {\bf{d}}{| }^{2}}{2{(2\pi )}^{3}}\displaystyle \frac{{\omega }_{0}^{2}}{{\omega }_{k}}\left(\displaystyle \frac{[1+{N}_{k}]}{{\omega }_{0}-{\omega }_{k}}+\displaystyle \frac{{N}_{k}}{{\omega }_{0}+{\omega }_{k}}\right)={\delta }_{e,\mathrm{CG}}^{(1)},\end{eqnarray} \tag{ 80 }$$

which completes the proof of equation (64a).

The shift appearing on the left-hand side of equation (64b) is found using equation (71) to be

$$\begin{eqnarray}{\delta }_{g}-{\delta }_{\mathrm{CG}}^{(2)} & = & \displaystyle \int {{\rm{d}}}^{3}k\displaystyle \sum _{\lambda }\displaystyle \frac{| {{\bf{e}}}_{\lambda }({\bf{k}})\cdot {\bf{d}}{| }^{2}}{2{(2\pi )}^{3}}\\ & & \times \left({\alpha }_{k}^{2}+{\alpha }_{k}({\alpha }_{k}-2)[1+2{N}_{k}]\displaystyle \frac{{\omega }_{0}}{{\omega }_{k}}+{\omega }_{0}\left[\displaystyle \frac{{{u}_{k}^{-}}^{2}{N}_{k}}{{\omega }_{k}-{\omega }_{0}}-\displaystyle \frac{{{u}_{k}^{+}}^{2}[1+{N}_{k}]}{{\omega }_{k}+{\omega }_{0}}\right]\right).\end{eqnarray} \tag{ 81 }$$

Similar calculations to those above for the excited state yield the final result

$$\begin{eqnarray}&&{\delta }_{g}-{\delta }_{\mathrm{CG}}^{(2)}=-\int {{\rm{d}}}^{3}k\displaystyle \sum _{\lambda }\displaystyle \frac{| {{\bf{e}}}_{\lambda }({\bf{k}})\cdot {\bf{d}}{| }^{2}}{2{(2\pi )}^{3}}\displaystyle \frac{{\omega }_{0}^{2}}{{\omega }_{k}}\left(\displaystyle \frac{[1+{N}_{k}]}{{\omega }_{0}+{\omega }_{k}}+\displaystyle \frac{{N}_{k}}{{\omega }_{0}-{\omega }_{k}}\right)={\delta }_{g,\mathrm{CG}}^{(1)}.\end{eqnarray} \tag{ 82 }$$

This completes the proof that the transition shift δ_e − δ_g is α_k-independent and that it equals Δ given in equation (12).

We remark that the need to account for the self-energy contributions along with the effect of the two-level truncation is a peculiarity of the single-dipole shift term Δ. The same need does not arise in the case of the remaining coefficients γ, γ₁₂ and Δ₁₂ in the standard two-dipole master equation (20). These coefficients are immediately seen to coincide with gauge-invariant matrix elements.

The joint shift Δ₁₂ resulting from the arbitrary gauge master equation derivation is given by

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{12}=\int \displaystyle \frac{{{\rm{d}}}^{3}k}{{(2\pi )}^{3}}\displaystyle \sum _{\lambda }| {{\bf{e}}}_{{\bf{k}}\lambda }\cdot {\bf{d}}{| }^{2}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{R}}}\left({\alpha }_{k}^{2}-1-\displaystyle \frac{{\omega }_{0}}{2}\left[\displaystyle \frac{{{u}_{k}^{+}}^{2}}{{\omega }_{k}+{\omega }_{0}}+\displaystyle \frac{{{u}_{k}^{-}}^{2}}{{\omega }_{k}-{\omega }_{0}}\right]\right).\end{eqnarray} \tag{ 83 }$$

Using equation (6) all α_k-dependence can be shown to vanish in the same way as with the single-dipole shifts dealt with in appendix A.1. The final result is

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{12}=\int \displaystyle \frac{{{\rm{d}}}^{3}k}{{(2\pi )}^{3}}\displaystyle \sum _{\lambda }| {{\bf{e}}}_{\lambda }({\bf{k}})\cdot {\bf{d}}{| }^{2}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{R}}}\displaystyle \frac{{\omega }_{k}^{2}}{{\omega }_{0}^{2}-{\omega }_{k}^{2}}.\end{eqnarray} \tag{ 84 }$$

Evaluating the angular integral and polarisation summation yields

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{12}=\displaystyle \frac{1}{\pi }{\int }_{0}^{\infty }{\rm{d}}{\omega }_{k}\,{d}_{i}{d}_{j}{\tau }_{{ij}}({\omega }_{k},R)\displaystyle \frac{{\omega }_{k}}{{\omega }_{0}^{2}-{\omega }_{k}^{2}}.\end{eqnarray} \tag{ 85 }$$

The integral is regularised by introducing a convergence factor ${{\rm{e}}}^{-\epsilon {\omega }_{k}}$ under the integral, and finally taking the limit $\epsilon \to {0}^{+}$. We substitute τ_ij given in equation (22) into equation (85) and evaluate the resulting integrals term by term. The integral arising from the first part of τ_ij is

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{\epsilon \to {0}^{+}}{\int }_{0}^{\infty }{\rm{d}}{\omega }_{k}\,{\omega }_{k}^{3}{{\rm{e}}}^{-\epsilon {\omega }_{k}}\displaystyle \frac{\sin {\omega }_{k}R}{{\omega }_{0}^{2}-{\omega }_{k}^{2}}=\displaystyle \frac{1}{2{\rm{i}}}\mathop{\mathrm{lim}}\limits_{\epsilon \to {0}^{+}}{\int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{k}\,{\omega }_{k}^{3}{{\rm{e}}}^{-\epsilon {\omega }_{k}}\displaystyle \frac{{{\rm{e}}}^{{\rm{i}}{\omega }_{k}R}}{{\omega }_{0}^{2}-{\omega }_{k}^{2}}.\end{eqnarray} \tag{ 86 }$$

We now make the substitution z = ω_k R, and make a suitable choice of contour C such that by the residue theorem we obtain

$$\begin{eqnarray}&&\displaystyle \frac{1}{2{{\rm{i}}{R}}^{2}}\mathop{\mathrm{lim}}\limits_{\epsilon \to {0}^{+}}{\int }_{C}{\rm{d}}{z}\,\displaystyle \frac{{z}^{3}{{\rm{e}}}^{{\rm{i}}{z}-\epsilon z/R}}{{({\omega }_{0}R)}^{2}-{z}^{2}}=-\displaystyle \frac{\pi }{2}{\omega }_{0}^{2}\cos {\omega }_{0}R.\end{eqnarray} \tag{ 87 }$$

Thus, the part of the shift Δ₁₂ arising from the first part (R⁻¹ component) of τ_ij is

$$\begin{eqnarray}&&-\displaystyle \frac{{\omega }_{0}^{2}}{4\pi R}({\delta }_{{ij}}-{\hat{R}}_{i}{\hat{R}}_{j}){d}_{i}{d}_{j}\cos {\omega }_{0}R,\end{eqnarray} \tag{ 88 }$$

which we recognise as the R⁻¹ component of Δ₁₂ in equation (21). The remaining parts of equation (85) can be evaluated in a similar way, which yields the final result given in equation (21).

We denote the dynamical map governing evolution of the reduced density matrix by F(t, t'), which is such that F(t, t')ρ(t') = ρ(t). A general two-time correlation function for arbitrary system observables O and O' can be written [36]

$$\begin{eqnarray}&&\langle O(t)O^{\prime} (t^{\prime} )\rangle =\mathrm{tr}({OF}(t,t^{\prime} )O^{\prime} F(t^{\prime} )\rho ).\end{eqnarray} \tag{ 89 }$$

We define the super-operator Λ by $\dot{\rho }(t)={\rm{\Lambda }}\rho (t)$ using the master equation (equation (20) or equation (40)). Since Λ is time-independent, from the initial condition $F(0,0)\equiv I$ we obtain the general solution $F(t,t^{\prime} )={{\rm{e}}}^{{\rm{\Lambda }}(t-t^{\prime} )}$. For convenience we write F(t,0) = F(t), so that $F(t,t^{\prime} )=F(t-t^{\prime} )$.

In order to calculate the two-time correlation functions we first find a concrete representation of the maps Λ and F(t). For this purpose we introduce a basis of operators denoted $\{{x}_{i}:i=1,\,\ldots ,\,16\}$, which is closed under Hermitian conjugation. The trace defines an inner-product $\langle O,O^{\prime} \rangle =\mathrm{tr}({O}^{\dagger }O^{\prime} )$ with respect to which the basis x_i is assumed to be orthonormal. We identify two resolutions of unity as ${\sum }_{i}\mathrm{tr}({x}_{i}^{\dagger }\cdot ){x}_{i}=I={\sum }_{i}\mathrm{tr}({x}_{i}\cdot ){x}_{i}^{\dagger }$, which imply that any operator O can be expressed as $O={\sum }_{i}\mathrm{tr}({x}_{i}^{\dagger }O){x}_{i}={\sum }_{i}\mathrm{tr}({x}_{i}O){x}_{i}^{\dagger }$. Expressing both sides of the equation $\dot{F}(t)={\rm{\Lambda }}F(t)$ in the basis x_i yields the relation

$$\begin{eqnarray}&&{\dot{F}}_{{jk}}(t)=\displaystyle \sum _{l}{{\rm{\Lambda }}}_{{jl}}{F}_{{lk}}(t),\end{eqnarray} \tag{ 90 }$$

where

$$\begin{eqnarray}&&{F}_{{jk}}(t)=\mathrm{tr}[{x}_{j}^{\dagger }F(t){x}_{k}],\qquad {{\rm{\Lambda }}}_{{jl}}=\mathrm{tr}[{x}_{j}^{\dagger }{\rm{\Lambda }}{x}_{l}].\end{eqnarray} \tag{ 91 }$$

Equation (90) can be written in the matrix form $\dot{{\bf{F}}}={\boldsymbol{\Lambda }}{\bf{F}}(t)$ whose solution is expressible in the matrix exponential form ${\bf{F}}(t)={{\rm{e}}}^{{\boldsymbol{\Lambda }}t}$. A general two-time correlation function of system operators can then be expressed using equation (89) as

$$\begin{eqnarray}&&\langle O(t)O^{\prime} (t^{\prime} )\rangle =\displaystyle \sum _{{ijkl}}\mathrm{tr}({{Ox}}_{i}){F}_{{ij}}(t-t^{\prime} )\mathrm{tr}({x}_{j}^{\dagger }O^{\prime} {x}_{k}){F}_{{kl}}(t^{\prime} )\mathrm{tr}({x}_{l}^{\dagger }\rho )={{\bf{O}}}^{{\rm{T}}}{\bf{F}}(t-t^{\prime} ){\bf{O}}^{\prime} {\bf{F}}(t^{\prime} ){\boldsymbol{\rho }},\end{eqnarray} \tag{ 92 }$$

where ${O}_{i}=\mathrm{tr}({{Ox}}_{i})$, ${\rho }_{i}=\mathrm{tr}({x}_{i}^{\dagger }\rho )$ and ${O}_{{ij}}^{{\prime} }=\mathrm{tr}({x}_{i}^{\dagger }O^{\prime} {x}_{j})$. Choosing the basis $\{{x}_{i}\}$ to be the operators obtained by taking the outer products of the bare states $| n,m\rangle ,\,(n,m=e,g)$, the above machinery can be used to obtain the correlation function (50).

The mode expansion for the transverse field canonical momentum ${{\boldsymbol{\Pi }}}_{{\rm{T}}}$ is

$$\begin{eqnarray}&&{{\boldsymbol{\Pi }}}_{{\rm{T}}}(t,{\bf{x}})=-{\rm{i}}\displaystyle \sum _{{\bf{k}}\lambda }\sqrt{\displaystyle \frac{{\omega }_{k}}{2{L}^{3}}}{{\bf{e}}}_{{\bf{k}}\lambda }{a}_{{\bf{k}}\lambda }(t){{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{x}}}+{\rm{h}}.{\rm{c}}.\end{eqnarray} \tag{ 93 }$$

This operator represents a different physical observable for each choice of ${\alpha }_{k}$, because it does not commute with the generalised gauge transformation ${R}_{\{{\alpha }_{k}\}}$. Similarly the photonic operators ${a}_{{\bf{k}}\lambda }$ are implicitly different for each choice of α_k. In the multipolar gauge the field canonical momentum coincides with the total electric field away from the sources; ${{\boldsymbol{\Pi }}}_{{\rm{T}}}({\bf{x}})=-{\bf{E}}({\bf{x}}),\,{\bf{x}}\ne {{\bf{R}}}_{\mu }$. The positive frequency (annihilation) and negative frequency (creation) components of the electric field are therefore defined for ${\bf{x}}\ne {{\bf{R}}}_{\mu }$ by

$$\begin{eqnarray}&&{{\bf{E}}}^{(+)}(t,{\bf{x}})={\rm{i}}\displaystyle \sum _{{\bf{k}}\lambda }\sqrt{\displaystyle \frac{{\omega }_{k}}{2{L}^{3}}}{{\bf{e}}}_{{\bf{k}}\lambda }{a}_{{\bf{k}}\lambda }(t){{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{x}}},\qquad {{\bf{E}}}^{(-)}(t,{\bf{x}})={{\bf{E}}}^{(+)}{(t,{\bf{x}})}^{\dagger },\end{eqnarray} \tag{ 94 }$$

where a_kλ is the photon annihilation operator within the multipolar gauge. For a system of two dipoles the integrated Heisenberg equation for the multipolar photon annihilation operator yields the source component

$$\begin{eqnarray}&&{a}_{{\bf{k}}\lambda ,s}(t)=\sqrt{\displaystyle \frac{{\omega }_{k}}{2V}}\displaystyle \sum _{\mu =1}^{2}{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{R}}}_{\mu }}{\int }_{0}^{t}{\rm{d}}{t}^{\prime} {{\rm{e}}}^{-{\rm{i}}{\omega }_{k}(t-t^{\prime} )}{{\bf{e}}}_{{\bf{k}}\lambda }\cdot {{\bf{d}}}_{\mu }(t^{\prime} ).\end{eqnarray} \tag{ 95 }$$

Since the dipole moment operators d_μ commute with the transformation ${R}_{\{{\alpha }_{k}\}}$ they represent the same physical observable for each choice of α_k. This implies that equation (95) can be expressed in terms of Coulomb gauge raising and lowering operators ${\sigma }_{\mu }^{\pm }$ in the two-level approximation, despite the implicit difference between these operators and their counterparts defined within the multipolar gauge. We subsequently express the Coulomb gauge operators ${\sigma }_{\mu }^{\pm }$ in the dressed basis $| {\epsilon }_{n}\rangle$ to obtain

$$\begin{eqnarray}&&{a}_{{\bf{k}}\lambda ,s}(t)=\sqrt{\displaystyle \frac{{\omega }_{k}}{2V}}{{\bf{e}}}_{{\bf{k}}\lambda }\cdot {\bf{d}}\displaystyle \sum _{\mu =1}^{2}{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {{\bf{R}}}_{\mu }}{\int }_{0}^{t}{\rm{d}}{t}^{\prime} {{\rm{e}}}^{-{\rm{i}}{\omega }_{k}(t-t^{\prime} )}\displaystyle \sum _{{nm}}{\sigma }_{\mu ,{nm}}{\theta }_{{nm}}(t^{\prime} ),\end{eqnarray} \tag{ 96 }$$

where ${\sigma }_{\mu ,{nm}}={\sigma }_{\mu ,{nm}}^{+}+{\sigma }_{\mu ,{nm}}^{-}$, ${\epsilon }_{{nm}}={\epsilon }_{n}-{\epsilon }_{m}$, and ${\theta }_{{nm}}=| {\epsilon }_{n}\rangle \langle {\epsilon }_{m}|$. We now perform a rotating-wave approximation, which eliminates terms that are rapidly oscillating within the interaction picture defined by the dressed Hamiltonian H_d given in equation (30). Substitution of the resulting expression into equation (94) yields in the mode continuum limit

$$\begin{eqnarray}&&{{\bf{E}}}_{s}^{(+)}(t,{\bf{x}})\\ &&\quad ={\rm{i}}\displaystyle \int {{\rm{d}}}^{3}k\displaystyle \sum _{\lambda }\displaystyle \frac{{\omega }_{k}}{2{(2\pi )}^{3}}{{\bf{e}}}_{\lambda }({\bf{k}})[{{\bf{e}}}_{\lambda }({\bf{k}})\cdot {\bf{d}}]\displaystyle \sum _{\mu =1}^{2}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{nm}}{n\lt m}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {{\bf{r}}}_{\mu }}{\displaystyle \int }_{0}^{t}{\rm{d}}{t}^{\prime} {{\rm{e}}}^{-{\rm{i}}{\omega }_{k}(t-t^{\prime} )}{{\rm{e}}}^{{\rm{i}}{\epsilon }_{{nm}}t^{\prime} }\,{\sigma }_{\mu ,{nm}}\,{\tilde{\theta }}_{{nm}}(t^{\prime} ),\end{eqnarray} \tag{ 97 }$$

where ${\tilde{\theta }}_{{nm}}(t^{\prime} )$ denotes the operator ${\theta }_{{nm}}(t^{\prime} )$ transformed into the interaction picture with respect to H_d, and ${{\bf{r}}}_{\mu }={\bf{x}}-{{\bf{R}}}_{\mu }$. Performing the angular integration and polarisation summation, and retaining only the radiative component yields

$$\begin{eqnarray}&&{E}_{s,\mathrm{rad},i}^{(+)}(t,{\bf{x}})\\ &&\quad =\,\displaystyle \frac{{\rm{i}}}{4{\pi }^{2}}\displaystyle \sum _{\mu =1}^{2}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{nm}}{n\lt m}}({\delta }_{{ij}}-{\hat{r}}_{\mu ,i}{\hat{r}}_{\mu ,j}){{\rm{d}}}_{j}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{k}{\displaystyle \int }_{0}^{t}{\rm{d}}{t}^{\prime} {\omega }_{k}^{2}\displaystyle \frac{\sin ({\omega }_{k}{r}_{\mu })}{{r}_{\mu }}{{\rm{e}}}^{-{\rm{i}}{\omega }_{k}(t-t^{\prime} )}{{\rm{e}}}^{{\rm{i}}{\epsilon }_{{nm}}t^{\prime} }{\sigma }_{\mu ,{nm}}{\tilde{\theta }}_{{nm}}(t^{\prime} ).\end{eqnarray} \tag{ 98 }$$

Finally, using the Markov approximation

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{k}\,f({\omega }_{k}){{\rm{e}}}^{{\rm{i}}({\omega }_{k}+{\epsilon }_{{nm}})t^{\prime} }[{{\rm{e}}}^{-{\rm{i}}{\omega }_{k}(t-{r}_{\mu })}-{{\rm{e}}}^{-{\rm{i}}{\omega }_{k}(t+{r}_{\mu })}]\\ &&\quad \approx \,f({\epsilon }_{{mn}}){\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{k}\,{{\rm{e}}}^{{\rm{i}}({\omega }_{k}+{\epsilon }_{{nm}})t^{\prime} }[{{\rm{e}}}^{-{\rm{i}}{\omega }_{k}(t-{r}_{\mu })}-{{\rm{e}}}^{-{\rm{i}}{\omega }_{k}(t+{r}_{\mu })}]\\ &&\quad =\,2\pi f({\epsilon }_{{mn}}){{\rm{e}}}^{{\rm{i}}{\epsilon }_{{nm}}t^{\prime} }[\delta (t^{\prime} -(t-{r}_{\mu }))-\delta (t^{\prime} -(t+{r}_{\mu }))],\end{eqnarray} \tag{ 99 }$$

valid for a suitably behaved function f, we obtain the final result equation (55) given in the main text.

To calculate the spectrum according to our master equation (40) we choose the basis of operators {x_i} used within the general method laid out in appendix A.3 as that obtained by taking the outer-products of the basis states $| {\epsilon }_{n}\rangle$ given in equation (32). Using equation (92) we define the array of correlation functions

$$\begin{eqnarray}&&{C}_{{nmp}}(t,t^{\prime} )=\langle {x}_{n}^{\dagger }(t){x}_{m}(t^{\prime} ){\rangle }_{{x}_{p}}={({\bf{F}}(t-t^{\prime} ){{\bf{X}}}_{m}{\bf{F}}(t^{\prime} ))}_{{np}},\end{eqnarray} \tag{ 100 }$$

where p is restricted to values such that x_p is diagonal, and where the matrix ${{\bf{X}}}_{m}$ has elements ${({{\bf{X}}}_{m})}_{{jk}}=\mathrm{tr}({x}_{j}^{\dagger }{x}_{m}{x}_{k})$. Taken in the symmetric state θ₃₃ the correlations appearing in equation (56) are all elements of the array ${C}_{{nm}}(t,t^{\prime} )$, which is given by equation (100) with x_p = θ₃₃. We choose a labelling whereby the x_i are given by

$$\begin{eqnarray}{x}_{i} & = & | {\epsilon }_{1}\rangle \langle {\epsilon }_{i}| ,\,\,\,i=1,\,\ldots ,\,4,\\ {x}_{i} & = & | {\epsilon }_{2}\rangle \langle {\epsilon }_{i-4}| ,\,\,\,i=5,\,\ldots ,\,8,\\ {x}_{i} & = & | {\epsilon }_{3}\rangle \langle {\epsilon }_{i-8}| ,\,\,\,i=9,\,\ldots ,\,12,\\ {x}_{i} & = & | {\epsilon }_{4}\rangle \langle {\epsilon }_{i-12}| ,\,\,\,i=13,\,\ldots ,\,16.\end{eqnarray} \tag{ 101 }$$

In this case the only non-zero off-diagonal element of ${C}_{{nm}}(t,t^{\prime} )$ is ${C}_{\mathrm{1,11}}(t,t^{\prime} )$ where x₁ = θ₁₁ and x₁₁ = θ₃₃. Furthermore the diagonal elements ${C}_{{nn}}(t,t^{\prime} )$ are zero unless n is odd. It follows that the only non-vanishing correlations in equation (56) are C₃₃(t, t') and C₇₇(t, t'). Moreover, since ${\sum }_{\mu ,\nu =1}^{2}{\sigma }_{\mu ,32}{\sigma }_{\nu ,23}=0$ only the term involving C₃₃(t, t') contributes. This term describes correlations associated with the symmetric to ground state transition and is given by equation (57) in the main text. Integration of this correlation function according to equation (46) then yields the spectrum in equation (59).