A master equation for strongly interacting dipoles

We consider a pair of dipoles for which direct electrostatic dipole-dipole interactions may be significantly larger than the coupling to transverse radiation. We derive a master equation using the Coulomb gauge, which naturally enables us to include the inter-dipole Coulomb energy within the system Hamiltonian rather than the interaction. In contrast, the standard master equation for a two- dipole system, which depends entirely on well-known gauge-invariant S-matrix elements, is usually derived using the multipolar gauge, wherein there is no explicit inter-dipole Coulomb interaction. We show using a generalised arbitrary-gauge light-matter Hamiltonian that this master equation is obtained in other gauges only if the inter-dipole Coulomb interaction is kept within the interaction Hamiltonian rather than the unperturbed part as in our derivation. Thus, our master equation, while still gauge-invariant, depends on different S-matrix elements, which give separation-dependent corrections to the standard matrix elements describing resonant energy transfer and collective decay. The two master equations coincide in the large separation limit where static couplings are negligible. We provide an application of our master equation by finding separation-dependent corrections to the natural emission spectrum of the two-dipole system.


Introduction
Dipole-dipole interactions are central to several important effects in atomic and molecular physics. Early studies by Eisenschitz, London and Förster [1,2] treated dipolar interactions as perturbative effects arising from direct electrostatic coupling.
Molecular quantum electrodynamics (QED) extends these treatments by incorperating retardation effects due to finite signal propagation. As was first shown by Casimir and Polder [3], a striking retardation effect occurs at large separations R/λ 1 where the R −6 dependence of the dispersion energy is increasingly replaced by an R −7 dependence.
In order to study the dynamics of systems of interacting dipoles open quantum systems theory has proven useful [4]. The master equation formalism can be used to obtain dynamical information about state populations and coherences, and to obtain fluorescence spectra [5,6,7,8]. As will be confirmed in this work, the standard second-order Born-Markov-secular master equation describing two dipoles within a common radiation reservoir depends entirely on well-known quantum electrodynamic (QED) matrix elements. These matrix elements describe dipole-dipole coupling and decay with retardation effects included. This master equation is obtained by treating the direct electrostatic coupling between the dipoles as a perturbation along with the coupling to transverse radiation. However, it is clear that if the former is sufficiently strong this approach may not be justified, in analogy with the case of externally imposed interactions [9]. Here we consider a system of free dipoles strongly coupled by dipoledipole interactions. Our focus is on discerning the full dependence of the physics on the interdipole separation. We also delineate how microscopic gauge-freedom effects the ensuing master equation derivation.
An important class of systems strongly coupled by dipole-dipole interactions are Rydberg atoms, which have been of interest for some time [10]. In recent years dipole-dipole interactions of Ryberg atoms have been the subject of numerous experimental and theoretical works [11,12,13,14,15,16,17,18,19,20,21,22,23]. Recently the first experimental confirmation of Förster resonant energy transfer was demonstrated using two Rydberg atoms separated by 15µm [14]. This type of resonant energy transfer is an important mechanism within photosynthesis, whose quantum nature is of continued interest within open quantum systems theory [24]. Dipoledipole interactions of Rydberg atoms also offer promising means of implementing quantum gates in which adjacent Rydberg states are treated as effective two-level systems and dipolar interactions are tuned with the use of lasers [21].
Such adjacent Rydberg states are typically separated by microwave transitions, which for small enough separations can be matched or even exceeded by the electrostatic dipole-dipole interaction strength divided by . Thus, a novel regime of strong electrostatic coupling occurs, in which the usual weak-coupling theory is expected to break down. A repartitioning of the Hamiltonian is necessary in order to identify a genuinely weak system-reservoir interaction, which can then constitute the starting point for perturbation theory. More specifically, we include the direct inter-dipole Coulomb energy within the unperturbed part of the Hamiltonian and only treat the coupling to transverse radiation as a weak perturbation. The master equation we derive exhibits a different dependence on the inter-dipole separation, and this has important consequences for the predicted physics. The rates of collective decay and resonant energy transfer are altered, as are the properties of the light emitted by the system.
There are five sections in this paper. We begin in Section 2 by reviewing the standard one and two dipole master equations in the Born-Markov and secular approximations. We show how the standard two-dipole master equation can be obtained for various choices of gauge for the microscopic Hamiltonian. Our purpose is to clearly identify limitations in the standard derivation, which is usually always performed using the multipolar Hamiltonian [4]. This concrete form of the Hamiltonian is the form obtained by choosing the multipolar gauge, also known as the Poincaré gauge [25]. In Section 3 we derive an alternative master equation describing the two-dipole system, which only reduces to the standard result in the limit of vanishing direct electrostatic coupling between the dipoles. This occurs in the limit of large separation. In Section 4 we solve the master equation derived in Section 3 and compare the solution with that of the standard master equation. We also obtain corrections to the emission spectrum of the two-dipole system. Finally in Section 5 we summarise our findings. We assume natural units = 0 = c = 1 throughout.

Gauge-invariant master equations 2.1. Single-dipole Hamiltonian and master equation
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 , |e ) of the dipole separated by energy ω 0 = ω e − ω g . Associated raising and lowering operators are defined by σ + = |e g| and σ − = |g e|. The electromagnetic bath is described by creation and annihilation operators a † kλ , a kλ for a single photon with momentum k and polarisation λ. The photon frequency is denoted ω k = |k|.
The energy of the dipole-field system is given by a Hamiltonian of the form defines the free (unperturbed) Hamiltonian and V denotes the interaction Hamiltonian. Gaugefreedom within the microscopic description results in the freedom to choose a number of possible interaction Hamiltonians. We define the generalised-gauge transformation [26] α k e kλ a kλ e ik·x + H.c. (2) In this expression the α k are real and dimensionless, the e kλ , λ = 1, 2 are mutually orthogonal polarisation unit vectors, which are both orthogonal to k, L 3 is the volume of the assumed fictitious quantisation cavity, andd denotes the dipole moment operator. By making the twolevel approximation after having transformed the Coulomb gauge Hamiltonian using the unitary operator R {α k } we obtain the Hamiltonian and 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 The coupling constant g kλ and the (real) coefficients u ± k are defined as where d and ω 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 Eq. (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,28,29,30,31], the PZW transformation is in fact a special case of a broader class of unitary gauge-fixing transformations [31,32,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 Eq. (3) is afforded by making the choice α k = ω 0 /(ω 0 +ω 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 ≡ 0. The counter-rotating terms in the linear dipole-field interaction term V −V (2) 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 f i (t) = f | M (t) |i between eigenstates |f and |i of H 0 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 Eq. (3). In contrast onenergy-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 [35,25,28] For gauge-invariance to hold the unperturbed Hamiltonian H 0 must be identified as the same operator before and after the transformation by e iT . Note however, that the unperturbed Hamiltonian H 0 given in Eq. (1)  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 Ref. [36] for example, does not at any point involve self-energy contributions due to the V (2) term within the interaction V . In general however, this term does contribute to dipole level-shifts, as is shown in appendix 6.1. In the general case that the temperature of the radiation field is arbitrary, the self-energy contributions from V (2) can be incorporated into the master equation by defining the Hamiltoniañ where the excited and ground state self-energy shifts are defined as in which n = e, g and ρ eq F (β) = e −β kλ ω k a † kλ a kλ /tr(e −β kλ ω k a † kλ a kλ ) with β the inverse temperature of the radiation field. Since V (2) is gauge-dependent we cannot include the selfenergy shifts within the unperturbed Hamiltonian H 0 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 Hamiltoniañ H d directly to obtain the second order master equatioṅ This master equation automatically includes the level shifts due to V (2) within the unitary evolution part, but the rest of the master equation is expressed in terms of the original partition H = H 0 + V . Using this partition where H 0 and V are given in Eqs. (1) and (3) respectively, Eq. (10) yields the α k -independent resulṫ Here N = 1/(e βω 0 − 1),ω 0 = ω 0 + ∆ and where the continuum limit for wavevectors k has been applied and N k = 1/(e βω k − 1). Further details of the calculations leading to the final result for ∆ in Eq. (12) are given in Appendix 6.1. We note that for N k = 0 we have 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 Eq. (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 0 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.

Arbitrary gauge derivation of the standard two-dipole master equation
Let us now turn our attention to obtaining the analogous result to Eq. (11) for the case of two identical interacting dipoles at positions R 1 and R 2 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 ω 0 respectively. The wavelength corresponding to ω 0 is denoted with λ 0 . An important quantity in the two-dipole system dynamics is the interdipole separation R = |R 2 − R 1 |. In terms of R we can identify in the usual way three distinct parameter regimes: R λ 0 is the near zone in which R −3 -dependent terms dominate, R ∼ λ 0 is the intermediate zone in which R −2 -dependent terms may become significant, and R λ 0 is the far-zone (radiation zone) in which R −1 -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 Power-Zienau-Woolley gauge transformation by [34,26] where d µ denotes the dipole moment operator of the µ'th dipole. We now transform the dipole approximated Coulomb gauge Hamiltonian using R {α k } and afterwards make the two-level approximation for each dipole. This implies that d µ = d(σ + µ + σ − µ ), and that the Hamiltonian and Analogously to the single-dipole case the first line in Eq. (16) defines a linear dipole-field interaction component while the term V (2) {α k } consists of self-energy contributions for each dipole and the radiation field; where m is the dipole mass. Due to the two-level approximation the first term is proportional to the identity, while the second term is a radiation self-energy term. The fieldÃ is defined as in the singe-dipole case by Eq. (5). The final term in Eq. (16) C {α k } , has no analog in the single-dipole Hamiltonian. This term gives a static Coulomb-like interaction between the dipoles, which is independent of the field; where σ x µ = σ + µ + σ − µ . In the Coulomb gauge (α k = 0) the term C {α k } reduces to the usual dipole-dipole Coulomb interaction. In the multipolar gauge C {α k } vanishes, and the interaction in Eq. (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 {α k } does not commute with H 0 given in Eq. (15) implying that H 0 represents a different physical observable for each choice of α k . More generally, since R {α 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 gaugeinvariant 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 {α k } vanishes identically. To obtain the same master equation for any other choice of α k = 1, we must include C {α k } within the interaction Hamiltonian V . The reason is that H 0 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 {α k } is kept within the interaction Hamiltonian V . To do this we substitute the interaction Hamiltonian in Eq. (16) into the second order Born-Markov master equation in the interaction picture with respect to H 0 , which is given byρ 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 2 and perform a further secular approximation, which neglects terms oscillating with twice the transition frequency ω 0 . Transforming back to the Schrödinger picture and including the single-dipole selfenergy contributions as in Eq. (10), we arrive after lengthy but straightforward manipulations at the final α k -independent resulṫ This equation is identical in form to the standard two-dipole master equation, which can be found in Ref. [4] for example. The coefficients within the master equation are as follows. The decay rates γ µν are given by where In Eq. (21) and throughout we denote spatial components with Latin indices and adopt the convention that repeated Latin indices are summed. The quantity γ 12 denotes an R-dependent collective decay rate. The third equality in Eq. (21) wherein γ 12 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 ∆ =ω 0 − ω 0 and ∆ 12 . 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 [cf Eq. (10)]. The shift ∆ is therefore as in Eq. (12). Details of the calculation of the joint shift ∆ 12 are given in appendix 6.2 with the final result being where |f = |g, e, 0 , |i = |e, g, 0 and |n are eigenstates of H 0 and As indicated by the second equality in Eq. (23) ∆ 12 is nothing but the well-known gauge-invariant QED matrix-element describing resonant energy-transfer.
We have therefore obtained the standard result, Eq. (20), without ever making a concrete choice for the α k . In order that the standard result is obtained the direct Coulomb-like interaction C {α k } must be kept within the interaction Hamiltonian V . This ensures that for all α k the unperturbed Hamiltonian H 0 is that used in conventional derivations wherein α k = 1. The S-matrix elements involving C {α 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 {α k } is sufficiently strong compared with the linear dipolefield coupling term, its inclusion within the interaction Hamiltonian rather than H 0 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 {α 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 This quantity coincides with the near-field limit of the resonant energy transfer element ∆ 12 given in Eq. (23). In the near-field regime R/λ 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.

Derivation of an alternative master equation
In the near-field regime R/λ 0 1 the rate of spontaneous emission into the transverse field is much smaller than the direct dipolar coupling; C/γ 1. Moreover, for a system of closely spaced Rydberg atoms, the electrostatic Coulomb interaction may be such that C ∼ ω 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 ∼ 10 −26 Cm where a 0 is the Bohr radius and e the electronic charge. For a 1µm separation, which is approximately equal to 10n 2 a 0 , the electrostatic dipole interaction C/ω 0 ∼ 1 for ω 0 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 twolevel 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 [37,38,4]. 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 and The contribution of the transverse field to ∆ 12 in Eq. (21) is found using Eq. (26) to bẽ When the contribution C = 0, e, g| C {0} |0, g, e resulting from the direct Coulomb interaction is added to∆ 12 the fully retarded result ∆ 12 is obtained. The two matrix elements ∆ 12 and ∆ 12 therefore only differ in their near-field components, which vary as R −3 and which we denote by ∆ nf 12 and∆ nf 12 , respectively. According to Eqs. (23) and (28) the components ∆ nf 12 and∆ nf 12 dominate at low frequencies ω 0 . Since ∆ 12 is evaluated at resonance ω k = ω 0 , 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 ∆ nf 12 ≈ C. In such regimes the multipolar interaction Hamiltonian cannot be classed as a weak perturbation. On the other hand, the matrix element∆ 12 is obtained using the Coulomb gauge interaction in Eq. (26), and is such that∆ nf 12 = ∆ nf 12 − C ≈ 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 Eq. (26), which therefore represents a genuinely weak perturbation.
The collective decay rate γ 12 as given in Eq. (21) does not involve the direct Coulomb interaction C {0} in any way, and can be obtained from the transverse field interaction in Eq. (26) or from the multipolar interaction. Crucially, in the near-field regime R/λ 0 1 the terms γ, γ 12 and∆ 12 , which result from the interaction in Eq. (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 where C ∈ R is given by Eq. (25). The corresponding interaction Hamiltonian is then given in Eq. (26). We begin by diagonalising H d as and with η = ω 2 0 + C 2 . Next we move into the interaction picture with respect to H 0 and substitute the interaction picture interaction Hamiltonian into Eq. (19). Moving back into the Schrödinger picture we eventually obtaiṅ Here, where The coefficients Γ µν (ω) are defined by where A I (x, t) denotes the field A(x) in Eq. (27) once transformed into the interaction picture, and · β denotes the average with respect to the radiation thermal state at temperature β −1 . The Γ µν (ω) are symmetric Γ µν (ω) = Γ νµ (ω) and can be written where with N k = 1/(e βω k − 1). The frequency integrals in Eq. (38) are to be understood as principal values. The decay rates γ µν (ω) in Eq. (38) coincide with those found in the standard master equation (20) when evaluated at ω 0 , though are here evaluated at the frequencies ω 1,2 . The quantities S µν are related to the shifts ∆ and ∆ 12 , defined in Eqs. (12) and (23) respectively, by In deriving Eq. (33) we have not yet performed a secular approximation, in contrast to the derivation of Eq. (20). However, naively applying a secular approximation that neglects offdiagonal terms for which ζ = ζ in the summand in Eq. (33) would not be appropriate, because this would eliminate terms that are resonant in the limit C → 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 bẏ We are now in a position to compare our master equation (40) with the usual result in (20). In the limit C → 0 we have η → ω 0 so that ω 1,2 → ω 0 . The rates and shifts in Eq. (38) are then evaluated at ω 0 within Eq. (40). Also, the Hamiltonian H d tends to the bare Hamiltonian , and furthermore we have that Thus, taking the limit C → 0 in Eq. (40) one recovers Eq. (20) with ∆ 12 replaced by∆ 12 given in Eq. (28). However, since∆ 12 → ∆ 12 when C → 0, Eqs. (40) and (20)

Discussion: gauge-invariance of the new master equation
It is important to note that while our master equation (40) is generally different to the usual gauge-invariant result [Eq. (20)] there is no cause for concern regarding the issue of gaugeinvariance. 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 usual where H 0 is given by Eq. (15) and Our master equation (40) has also been obtained by choosing α k = 0, but our derivation makes use of the different partitioning H =H 0 + V whereH 0 is defined as in Eq. (29) and V = V usual − C {0} is defined as in Eq. (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 Eq. (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 Eq. (29). We note further that the secular approximation is well justified within the nearfield regime of interest R λ 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 0 given in Eq. (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 α k = 0 is contained entirely within non-secular contributions. These contributions are negligible within the regime of interest R λ 0 and have only been retained within Eq. (40) to facilitate comparison with the standard result Eq. (20). Moreover, in the far-field regime R λ 0 the master equations (20) and (40) coincide, so the master equation (40) is also gauge-invariant within this regime.

Solutions and emission spectrum 4.1. Solutions
For large inter-dipole separations R λ 0 the master equations (20) and (40) coincide and they therefore yield identical physical predictions. However, in the near-zone R λ 0 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 | 3 . This state is a simultaneous eigenstate of the dipole Hamiltonian ω 0 (σ + 1 σ − 1 + σ + 2 σ − 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 | 3 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. Figs. 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 | 3 . 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 Eq. (20) (see Fig. 1). This gives rise to the different starting values at R = r a of the curves depicted in Fig. 2. For larger separations the solutions converge and become indistinguishable for all times.
The different behaviour in Figs. 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 which is different to the usual transition dipole moment 3 | d 1 +d 2 |gg = √ 2d. Since a → 1/ √ 2 as C → 0, the dipole moment d 31 reduces to √ 2d when R → ∞. As R decreases, however, d 13 becomes increasingly large compared with √ 2d. This is consistent with the more rapid decay observed in Fig. 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 Eq. (32), we obtain Only when C → 0, such that ω 2 = η + C → ω 0 and c → 1/ √ 2, does this decay rate reduce to that obtained when using the bare eigenstates |i, j , (i, j = e, g), which is where γ 12 (ω 0 ) is given in Eq. (21). As shown in Fig. 3, for sufficiently small R the decay rate γ s (ω 2 ) is significantly larger than γ s,0 . Figure 1. The populations of the stationary state (| 1 or |g, g ) and symmetric state | 3 are plotted as functions of t for fixed separation R = 10r a , where r a = n 2 a 0 is a characteristic Rydberg atomic radius, with n = 50 and a 0 the Bohr radius. We have assumed no thermal occupation of the field, N = 0 and that the transition dipole moment d is orthogonal to the separation vector R. The transition frequency is chosen in the microwave regime ω 0 = 10 10 . We use p g and p s to denote the stationary and symmetric state populations respectively, and we use p and p 0 to denote populations obtained from master equations (40) and (20), respectively. In the case of Eq. (40) the stationary state is | 1 whereas in the case of Eq. (20) the stationary state is simply |g, g . The initial condition chosen is p s = 1. The ground and symmetric state populations obtained from the master equation (40) crossover significantly earlier than those obtained from Eq. (20).
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 antisymmetric state | 2 , which like | 3 is a simultaneous eigenstate of ω 0 (σ + 1 σ − 1 + σ + 2 σ − 2 ) and H d . Both master equations predict that the population of the state | 2 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 or | 1 , is used. Finally, the predicted behaviour by our master equation (40) of the standard stationary state |g, g is illustrated in Fig. 4. For an initial state | 3 the population p gg (t) of the state |g, g at a given time t, is identical to that predicted by Eq. (20) only for sufficiently large R whereby | 1 ≈ |g, g .

Emission spectrum
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 | 3 . 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   Fig. 1. Figure 4. The population of the state |g, g found using Eq. (40) is plotted as a function of separation R for various times. All remaining parameters are as in Fig. 1. The dashed lines give the corresponding populations found using Eq. (20), which are insensitive to variations in R over the range considered. For large R the two sets of solutions agree. In particular, the steady state population p g,g (∞) = | g, g| 1 | 2 is equal to unity only for sufficiently large R.
photodetection by [39,40] 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 R, so that only the radiative component E s,rad of the electric source field, which varies as |x − R µ | −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 where r µ = x − R µ , and t µ = t − r µ is the retarded time associated with the µ'th source. For x R we have to a very good approximation that r 1 = r 2 = r, where r is the relative vector from x to the midpoint of R 1 and R 2 . Substituting Eq. (47) into Eq. (46) within this approximation yields where 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 | 3 , the standard master equation (20) together with the method of calculation given in appendix 6.3 we obtain the two-time correlation function where γ s,0 and ∆ 12 are given in Eqs. (45) and (21), respectively. By direct integration of Eq. (50) one obtains the corresponding Lorentzian spectrum Let us now turn our attention to the spectrum obtained from our new master equation (40). We have seen that the solutions of Eqs. (40) and (20) differ only in the near field regime R λ 0 . For sufficiently small R we have that ω 0 ∼ C, and the frequency difference ω 2 − ω 1 = 2C ∼ 2ω 0 is large. In this situation we can perform a full secular approximation within Eq. (40) to obtain the master equationρ and Solving this secular master equation allows us to obtain a simple expression for the emission spectrum. The correlation function in Eq. (46) defines the radiation intensity when it is evaluated at t = t . Naively calculating the quantity 2 µ,ν=1 σ + µ (t)σ − ν (t) 1 taken in the stationary state | 1 yields a non-zero stationary intensity, because | 1 is a superposition involving both |g, g and the doubly excited bare state |e, e . 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 Eq. (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 Eq. (30). One then obtains the source field where σ µ,nm = σ + µ,nm + σ − µ,nm in which σ ± µ,nm denotes the nm'th matrix element of σ ± µ in the basis | n . The transition frequencies associated with this basis are denoted nm = n − m , and the raising and lowering operators are denoted θ nm = | n m | , n = m. The derivation of Eq. (55) is given in Appendix 6.4. According to Eq. (55) the annihilation (creation) radiation source field is now associated with lowering (raising) operators in the dressed basis | i rather than in the bare basis |n, m , (n, m = e, g). Substitution of Eq. (55) into Eq. (46) yields where we have again assumed that r 1 = r 2 = r. Unlike the correlation function in Eq. (92), when t = t the correlation functions θ pq (t)θ nm (t ) , p > q, m > n vanish in the stationary (ground) state θ 11 . The radiation intensity is therefore seen to vanish in the stationary limit as required physically. Taken in the symmetric state θ 33 the only non-zero two-time correlation function that contributes to Eq. (56) is found to be is the shifted symmetric to ground transition frequency. Integration of Eq. (57) according to Eq. (56) then yields the Lorentzian spectrum Full details of the calculation of the spectrum in Eq. (59) are given in appendix 6.4. In the limit of large separation C → 0, which implies that ω 2 → ω 0 ,ω 2 →ω 0 + ∆ 12 , γ s (ω 2 ) → γ s,0 , and a → 1/ √ 2. As a result s(ω) → s 0 (ω) 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 0 (ω).
The two spectra s 0 (ω) and s(ω) are compared in Figs. 5 and 6. As their relative widths are proportional to the rates, they are given in Eqs. (45) and (44), respectively. These quantities have been plotted already in Fig. 3. The relative heights of the spectral peaks are s 0 (ω 0 +∆ 12 )/µ and s(ω 2 )/µ respectively, which are plotted in Fig. 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 2 a 0 , n = 50 is a characteristic Rydberg atomic radius, the peak value of s(ω) is around two times larger than the peak value of s 0 (ω) for the parameters chosen here. The positions of the peaks areω 0 + ∆ 12 and ω 2 , respectively, and these are plotted in Fig. 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 0 (ω). Fig. 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ω 2 −ω 0 − ∆ 12 is around 1 Ghz for the parameters chosen in Fig. 1. This is similar in magnitude to the Lamb-shift in atomic Hydrogen.
s(ω+ω2) s(ω2) Figure 5. The spectra s(ω) and s 0 (ω) are plotted with R = 10r a and with all remaining parameters as in Fig. 1. For this separation the peak heights and centres are quite different as shown in Figs. 7 and 8, respectively. Here, for illustrative purposes, the spectra have both been centred at zero and normalised by their respective peak heights. Figure 6. The spectra s(ω) and s 0 (ω) are plotted with R = 50r a and with all remaining parameters as in Fig. 1. For this separation the positions of the peak centres remain quite different on the frequency scale set by the width γ s (ω 2 ) ≈ γ s,0 . Here, for illustrative purposes, the spectra have both been centred at zero. However, for this value of R the peak heights are effectively the same. Therefore, the spectra have been normailsed by the same peak value s 0 (ω 0 + ∆ 12 ).

Conclusions
In this paper we have derived a partially secular master equation valid for arbitrarily separated dipoles within a common radiation field at arbitrary temperature. The equation is intended for the modelling of dipolar systems in which static dipole-dipole interactions are strong compared with the coupling to transverse radiation. This situation can arise in systems of Rydberg atoms and other molecular systems [10,11,12,13,14,15,16,17,18,19,20,21,22,23,41,42,43,44]. We have shown that the standard gauge-invariant two-dipole master equation can only be derived in gauges other than the multipolar gauge if the direct inter-dipole Coulomb energy is included within the interaction Hamiltonian rather than the unperturbed part. Our arbitrary gauge approach makes a particular limitation of this method clear. Specifically, the usual approach can only be justified when the direct Coulomb interaction is weak along with the coupling to transverse radiation. In situations in which this is not the case our master equation, which is based on a repartitioning of the Hamiltonian into unperturbed and interaction parts, yields significant corrections to previous results. In addition to corrections to the decay of the excited states of the system, we have found corrections to the natural emission spectrum of the initially excited system. In principle, spectroscopy could be used to determine which predictions are closer to the measured values. A possible extension of our result would be to include an external driving Hamiltonian that represents coherent irradiation. The techniques employed here could then be used to calculate the fluorescence spectrum of the driven system. Acknowledgment: This work was supported by the Engineering and Physical Sciences Research Council. We thank Jake Iles-Smith and Victor Jouffrey for useful discussions.
6. Appendix 6.1. Self-energy contributions and the Gauge-invariance of the single dipole-shift 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 (2) is given in Eq. (4). The shifts arising from this term are divergent in the mode continuum limit ω k → ∞, 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 (2) 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 (2) 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 (2) 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 n| (2) V |n , n = e, g must be taken into account.
Using the Hamiltonian in Eq. (3), the standard Born-Markov master equation has the form given in Eq. (11), in which the decay rate γ is independent of α k . The transition shift expressed as the difference between excited and ground state shifts asω 0 − ω 0 = ∆ = δ are α k -dependent. This α k -dependence is due to the lack of any contribution from the self-energy term V (2) in Eq. (60). The α k -dependence within the master equation is eliminated when one accounts for the selfenergy contributions and the effect of the two-level approximation, recalling that the latter was made after the transformation R {α 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 Eq. (11). First we note that we can continue to express the second line in Eq. (10) in terms of the original partition H = H 0 + V . Thus, provided H 0 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 δ (1) e,g coming from the second line in Eq. (10) to the corresponding self-energy contribution in Eq. (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 CG .
The components δ (1) e,CG are obtained by setting α k = 0 in Eq. (60), while the remaining component is the Coulomb gauge self-energy shift due to the A 2 T part of the Coulomb gauge interaction Hamiltonian. Since this term is independent of the dipole, the shift δ (2) CG is the same for the ground and excited levels. The single-dipole transition shift ∆ given in Eq. (12) in the main text can be expressed in terms of Coulomb gauge shifts as More generally, for arbitrary α k the total ground and excited state level shifts are denoted δ e,g .
In what follows we will show that δ e − δ and from which it follows using Eq. (63) that δ e − δ g = ∆ for all choices of α k . In order to show that Eqs. (64a) and (64b) hold we must carefully account for the two-level approximation, which was performed after the gauge transformation R {α k } . Let us consider a general shift of the m'th level of the dipole with the form 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 where ω 0 := ω eg = −ω ge and d := d eg = d * ge . If instead m = g then the shift is n ω ng |v · d ng | 2 = +ω 0 |v · d| 2 .
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 where we have made use of the identity The difference between the finite and infinite-dimensional cases arises because the proof of Eq. (69) rests directly on the CCR algebra [r i , p j ] = iδ ij , which can only be supported in infinite-dimensions. When the algebra is truncated to su(2), the same shift comes out leveldependent. Since the gauge transformation R {α k } is made on the infinite-dimensional dipole it is necessary to employ Eq. (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 Eqs. (66) and (68), which imply but when dealing with the ground level shift m = g we use Eqs. (67) and (68), which imply We now proceed to verify that Eqs. (64a) and (64b) hold. The complete shifts δ e,g are obtained by taking the shifts in Eq. (60) and adding their respective self-energy contributions. Subtracting δ (2) CG in Eq. (62) from δ e and subsequently using Eq. (70), which is appropriate for the excited state shift, we obtain Using Eq. (6) we express the bracket within the integrand in this expression in terms of α k . The part independent of N k is In this expression we identify the coefficient of α 2 k as and the coefficient of 2α k as Thus, Eq. (73) is α k -independent. The remaining part is The N k -dependent parts of δ e − δ CG can be dealt with in a similar manner. The coefficient of α 2 k in the N k -dependent part of the bracket within the integrand of the expression for δ e − δ Similarly, the coefficient of α k is The remaining N k -dependent part is Combining Eqs. (73) and (79) we obtain the α k -independent result which completes the proof of Eq. (64a). The shift appearing on the left-hand-side of Eq. (64b) is found using Eq. (71) to be Similar calculations to those above for the excited state yield the final result This completes the proof that the transition shift δ e − δ g is α k -independent and that it equals ∆ given in Eq. (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 γ, γ 12 and ∆ 12 in the standard two-dipole master equation (20). These coefficients are immediately seen to coincide with gauge-invariant matrix elements.

Calculation of the standard joint shift
The joint shift ∆ 12 resulting from the arbitrary gauge master equation derivation is given by Using Eq. (6) all α k -dependence can be shown to vanish in the same way as with the single-dipole shifts dealt with in Appendix 6.1. The final result is Evaluating the angular integral and polarisation summation yields The integral is regularised by introducing a convergence factor e − ω k under the integral, and finally taking the limit → 0 + . We substitute τ ij given in Eq. (22) into Eq. (85) and evaluate the resulting integrals term by term. The integral arising from the first part of τ ij is We now make the substitution z = ω k R, and make a suitable choice of contour C such that by the residue theorem we obtain Thus, the part of the shift ∆ 12 arising from the first part ( which we recognise as the R −1 component of ∆ 12 in Eq. (21). The remaining parts of Eq. (85) can be evaluated in a similar way, which yields the final result given in Eq. (21).

Method of calculation of the spectrum
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] We define the super-operator Λ byρ(t) = Λρ(t) using the master equation [Eq. (20) or Eq. (40)]. Since Λ is time-independent, from the initial condition F (0, 0) ≡ I we obtain the general solution 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, ..., 16}, which is closed under Hermitian conjugation. The trace defines an inner-product O, O = tr(O † O ) with respect to which the basis x i is assumed to be orthonormal. We identify two resolutions of unity as i tr(x † i ·)x i = I = i tr(x i ·)x † i , which imply that any operator O can be expressed as Expressing both sides of the equatioṅ F (t) = ΛF (t) in the basis x i yields the relatioṅ where Eq. (90) can be written in the matrix formḞ = ΛF(t) whose solution is expressible in the matrix exponential form F(t) = e Λt . A general two-time correlation function of system operators can then be expressed using Eq. (89) as where O i = tr(Ox i ), ρ i = tr(x † i ρ) and O ij = tr(x † i O x j ). Choosing the basis {x i } to be the operators obtained by taking the outer products of the bare states |n, m , (n, m = e, g), the above machinery can be used to obtain the correlation function (50).

Derivation of spectrum associated with the new master equation
The mode expansion for the transverse field canonical momentum Π T is This operator represents a different physical observable for each choice of α k , because it does not commute with the generalised gauge transformation R {α k } . Similarly the photonic operators a kλ 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; Π T (x) = −E(x), x = R µ . The positive frequency (annihilation) and negative frequency (creation) components of the electric field are therefore defined for x = R µ by 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 Since the dipole moment operators d µ commute with the transformation R {α k } they represent the same physical observable for each choice of α k . This implies that Eq. (95) can be expressed in terms of Coulomb gauge raising and lowering operators σ ± µ 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 σ ± µ in the dressed basis | n to obtain a kλ,s (t) = ω k 2L 3 e kλ · d where σ µ,nm = σ + µ,nm + σ − µ,nm , nm = n − m , and θ nm = | n m |. We now perform a rotatingwave approximation, which eliminates terms that are rapidly oscillating within the interaction picture defined by the dressed Hamiltonian H d given in Eq. (30). Substitution of the resulting expression into Eq. (94) yields in the mode continuum limit whereθ nm (t ) denotes the operator θ nm (t ) transformed into the interaction picture with respect to H d , and r µ = x − R µ . Performing the angular integration and polarisation summation, and retaining only the radiative component yields Finally, using the Markov approximation = 2πf ( mn )e i nmt [δ(t − (t − r µ )) − δ(t − (t + r µ ))], valid for a suitably behaved function f , we obtain the final result Eq. (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 6.3 as that obtained by taking the outer-products of the basis states | n given in Eq. (32). Using Eq. (92) we define the array of correlation functions where p is restricted to values such that x p is diagonal, and where the matrix X m has elements (X m ) jk = tr(x † j x m x k ). Taken in the symmetric state θ 33 the correlations appearing in Eq. (56) are all elements of the array C nm (t, t ), which is given by Eq. (100) with x p = θ 33 . We choose a labelling whereby the x i are given by In this case the only non-zero off-diagonal element of C nm (t, t ) is C 1,11 (t, t ) where x 1 = θ 11 and x 11 = θ 33 . Furthermore the diagonal elements C nn (t, t ) are zero unless n is odd. It follows that the only non-vanishing correlations in Eq. (56) are C 33 (t, t ) and C 77 (t, t ). Moreover, since 2 µ,ν=1 σ µ,32 σ ν,23 = 0 only the term involving C 33 (t, t ) contributes. This term describes correlations associated with the symmetric to ground state transition and is given by Eq. (57) in the main text. Integration of this correlation function according to Eq. (46) then yields the spectrum in Eq. (59).