The field-induced interaction between non-resonant magnetic dipoles

TL;DR

This paper derives a general magnetic dipole-dipole interaction model mediated by the quantized electromagnetic field, applicable to non-resonant dipoles and including counter-rotating and mixed terms, without relying on the rotating-wave approximation.

Contribution

It provides a novel derivation of dipole-dipole interactions that accounts for non-resonant conditions and includes previously neglected interaction terms.

Findings

Interaction strength depends on dipole frequency difference.

Includes counter-rotating and mixed interaction terms.

No rotating-wave approximation needed.

Abstract

We make a general derivation for the magnetic dipole-dipole interaction based on the mediation of the quantized electro-magnetic field. Due to the interaction with the dipoles, the dynamics of the field is added by a dipole field, which finally gives rise to the dipole-dipole interaction. Different from previous studies, the rotating-wave-approximation is no longer needed throughout this derivation, and our result naturally gives the interaction for non-resonant dipoles. Moreover, our derivation also gives the counter-rotating interaction terms, and even the mixed interaction terms between the permanent and transition dipoles. We notice that this field-induced interaction is associated with the interference of the virtual/real photons emitted from the two dipoles, thus the interaction strength could be influenced by the frequency difference of the two dipoles.

Figures

Equations

$$V=\frac{\mu_{0}}{4\pi r^{3}}\big{[}\vec{m}_{1}\cdot\vec{m}_{2}-3(\vec{m}_{1}\cdot\hat{\mathrm{e}}_{\mathbf{r}})(\vec{m}_{2}\cdot\hat{\mathrm{e}}_{\mathbf{r}})\big{]}.$$

$${\hat{H}}_{\mathrm{int}}=-\sum _{i=1,2}{\hat{\mathfrak{m}}}_i(t)\cdot \hat{\mathbf{B}}({\mathbf{x}}_i,t).$$

$$\hat{\mathbf{B}}(\mathbf{x},t)=\sum_{\mathbf{k}\sigma}\mathrm{i}\hat{\mathrm{e}}_{\mathbf{k}\check{\sigma}}Z_{\mathbf{k}}\Big{[}\hat{a}_{\mathbf{k}\sigma}(t)e^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}}-{\rm H.c.}\Big{]},$$

$${\hat{H}}_{\mathrm{int}}^{(1)}$$

$$g_{i,\mathbf{k}\sigma }^{xy}$$

$$\hat{a}_{\mathbf{k}\sigma}(t)=\hat{a}_{\mathbf{k}\sigma}(0)e^{\text{--}\mathrm{i}\omega_{\mathbf{k}}t}-\sum_{xy}\frac{\mathrm{i}}{\hbar}\big{(}g_{1,\mathbf{k}\sigma}^{xy}\big{)}^{*}\int_{0}^{t}ds\,e^{\text{--}\mathrm{i}\omega_{\mathbf{k}}s}\hat{\tau}_{1}^{yx}(t-s).$$

$${\hat{\mathbf{B}}}_0=$$

$${\hat{\mathbf{B}}}_{\text{D1}}=$$

$$+H.c.$$

$${\hat{H}}_{1\to 2}=$$

$$=$$

$$\mathcal{D}_{1\rightarrow 2}^{yx,uv}(s):=\sum_{\mathbf{k}\sigma}-\frac{\mathrm{i}}{\hbar}\big{(}g_{1,\mathbf{k}\sigma}^{xy}\big{)}^{*}g_{2,\mathbf{k}\sigma}^{uv}\,e^{-\mathrm{i}\omega_{\mathbf{k}}s}:=-\mathrm{i}\int_{0}^{\infty}\frac{d\omega}{2\pi}\,J_{1\rightarrow 2}^{yx,uv}(\omega)e^{-\mathrm{i}\omega s}.$$

$$J_{1\to 2}^{yx,uv}(\omega )=$$

$$\displaystyle-3(\vec{m}_{1}^{yx}\cdot\hat{\mathrm{e}}_{\mathbf{r}})(\vec{m}_{2}^{uv}\cdot\hat{\mathrm{e}}_{\mathbf{r}})\,\big{[}\frac{1}{3}\eta^{2}\sin\eta+\eta\cos\eta-\sin\eta\big{]}\Big{\}},$$

$${\hat{\tau }}_1^{yx}(t-s)\simeq {\hat{\tau }}_1^{yx}(t)\cdot e^{\mathrm{i}{\Omega }_1^{yx}s},$$

$${\hat{H}}_{1\to 2}=\sum _{xy,uv}{\xi }_{1\to 2}^{yx,uv}{\hat{\tau }}_1^{yx}{\hat{\tau }}_2^{uv}+H.c.,$$

$${\xi }_{1\to 2}^{yx,uv}=-\mathrm{i}{\int }_0^{t\to \infty }ds{\int }_0^{\infty }\frac{d\omega }{2\pi }{J}_{1\to 2}^{yx,uv}(\omega )e^{\mathrm{i}({\Omega }_1^{yx}-\omega )s}$$

$$=$$

$${\hat{H}}_{1\to 2}$$

$${\Lambda }_{1\to 2}^{yx,uv}$$

$$\displaystyle=K_{1\rightarrow 2}^{yx,uv}\big{(}\Omega_{1}^{yx}\big{)}-\frac{\mathrm{i}}{2}J_{1\rightarrow 2}^{yx,uv}\big{(}\Omega_{1}^{yx}\big{)},$$

$$K_{1\to 2}^{yx,uv}(\Omega ):=\mathbf{P}{\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }\frac{J_{1\to 2}^{yx,uv}(\omega )}{\Omega -\omega }$$

$$=$$

$$-$$

$${\hat{H}}_{1\leftrightarrow 2}$$

$$G_{yx,uv}^{(\mathtt{P})}$$

$$G_{yx,uv}^{(\mathtt{D})}$$

$${\hat{H}}_{xx,yy}:=(G_{xx,yy}^{(\mathtt{P})}+G_{xx,yy}^{(\mathtt{D})}){\hat{\tau }}_1^{xx}{\hat{\tau }}_2^{yy},$$

$${\hat{H}}_{xx,yy}=\frac{{\mu }_0}{4\pi r^3}[{\hat{\mathfrak{m}}}_1^{xx}\cdot {\hat{\mathfrak{m}}}_2^{yy}-3({\hat{\mathfrak{m}}}_1^{xx}\cdot {\hat{\mathrm{e}}}_{\mathbf{r}})({\hat{\mathfrak{m}}}_1^{yy}\cdot {\hat{\mathrm{e}}}_{\mathbf{r}})],$$

$${\hat{H}}_{\text{eg,ge}}=(G_{\text{eg,ge}}^{(\mathtt{P})}+G_{\text{eg,ge}}^{(\mathtt{D})}){\hat{\tau }}_1^{+}{\hat{\tau }}_2^{-}.$$

$$G_{\text{eg,ge}}^{(\mathtt{P})}$$

$$G_{\text{eg,ge}}^{(\mathtt{D})}$$

$${\hat{H}}_{\text{eg,eg}}$$

$$G_{\text{eg,eg}}^{(\mathtt{P})}$$

$$G_{\text{eg,eg}}^{(\mathtt{D})}$$

$${\hat{H}}_{\text{eg,gg}}$$

$$G_{\text{eg,gg}}^{(\mathtt{P})}$$

$$G_{\text{eg,gg}}^{(\mathtt{D})}$$

$$\hat{H}$$

Full text

I Introduction

The electric/magnetic dipole-dipole (DD) interaction exists in many different physical systems (Oh et al., 2008; Hu et al., 2018; Balasubramanian et al., 2008; Pervushin et al., 1997; Yang et al., 2010; El-Ganainy and John, 2013), such as the electric interaction between the dielectric defects in solid-state systems (Martinis et al., 2005a; Paik et al., 2011; Rigetti et al., 2012; Lisenfeld et al., 2016; Martinis et al., 2005b), and the magnetic interaction between the nitrogen-vacancy and the nuclear spins around (Zhao et al., 2012; Žnidarič et al., 2008; Doherty et al., 2013; Wrachtrup and Jelezko, 2006).

In classical electrodynamics, it is well-known that the magnetic dipole-dipole interaction is (Jackson, 1998)

[TABLE]

For quantum systems, usually the DD interaction is obtained by simply replacing the classical dipole moment $\vec{m}_{i}$ with the quantum operator $\hat{\mathfrak{m}}_{i}$.

But this interaction also can be considered in another way: in principle, the DD interaction should be induced by the mediation of the EM field. To carry out this idea, we can derive a Markovian master equation for two dipoles which are separated by a certain distance and locally interacting with the EM field, then the field-induced DD interaction appears as the off-diagonal Lamb shift in the unitary correction term (Lehmberg, 1970; Agarwal, 1974; Ficek et al., 1987; Craig and Thirunamachandran, 1998; Ficek and Swain, 2005).

For two-level dipoles, this derivation gave a transition interaction with the form $\hat{V}=g(\hat{\sigma}_{1}^{+}\hat{\sigma}_{2}^{-}+\hat{\sigma}_{1}^{-}\hat{\sigma}_{2}^{+})$ when their transition frequencies were exactly resonant. In this lengthy derivation, the rotating-wave-approximation (RWA) should be applied in a proper way, and the counter-rotating interaction terms cannot be obtained, which is different from the above direct exchange method from Eq. (1).

Even for this resonant transition interaction, it turned out the interaction strength $g$ does not exactly return to the above classical form (1), unless the distance between the two dipoles is much smaller than the transition wavelength $r/\lambda\ll 1$ (the Dicke limit) (Lehmberg, 1970; Agarwal, 1974; Ficek et al., 1987; Craig and Thirunamachandran, 1998; Ficek and Swain, 2005; Dicke, 1954).

We noticed that this difference can be understood if the permanent and transition dipoles are treated separately (Wang et al., 2018). By checking the dynamics of the EM field, we made another simplified way to derive the field-induced magnetic DD interaction, which gives both the permanent and transition DD interactions. The diagonal terms of the dipole operator correspond to the permanent dipole moments, and their interaction returns to the classical form directly. The off-diagonal terms are the transition dipole moments, which indeed do not have classical counterparts, thus their interaction intrinsically contains quantum corrections and does not need to return to the classical form exactly.

To make a more precise comparison with the interaction directly exchanged from the classical one, we still need to study the DD interaction between non-resonant dipoles, as well as the counter-rotating terms. In this paper, we make a more general derivation for the field-induced interaction between two magnetic dipoles with multiple energy levels. We calculate the dynamics of the quantized EM field, and it becomes the summation of the original vacuum field and the dipole field due to the interaction with the dipoles. Then the DD interaction is derived from the local interaction between the magnetic dipole and the incoming dipole field. The RWA is no longer needed through out this derivation, thus we can obtain the field-induced interaction for non-resonant dipoles, and even the counter-rotating interaction terms.

Different from the electric interactions, both the static and radiative magnetic interactions are related to the transverse modes of the EM field, and do not involve the longitudinal modes (Agarwal, 1974; Cohen-Tannoudji et al., 1989; Craig and Thirunamachandran, 1998). Therefore, our derivations can be well done under the Coulomb gauge and naturally gives both the permanent and transition DD interactions between two non-resonant dipoles. Moreover, our result also contains counter-rotating interaction terms and the mixed interaction terms between the permanent and transition dipoles, which are usually negligible under the RWA.

We notice that such field-induced DD interaction is associated with the interference of the virtual/real photons emitted from the two dipoles. As a result, the interaction strength could be influenced by the the frequency difference of the two dipoles.

This paper is organized as follows. In Sec. II, we introduce the model setup and study the dynamics of the quantized EM field induced by the magnetic dipoles. In Sec. III, we derive the general result for the field-induced interaction between two magnetic dipoles with multiple states. In Sec. IV, in the example of two-level dipoles, we discuss the different interaction types according to our general result. Sec. V is a brief summary.

II Field dynamics influenced by dipoles

We consider there are two magnetic dipoles placed in the EM field ($\hat{H}_{\text{{em}}}=\sum_{\mathbf{k}\sigma}\hbar\omega_{\mathbf{k}}\hat{a}_{\mathbf{k}\sigma}^{\dagger}\hat{a}_{\mathbf{k}\sigma}$). Each dipole has multiple energy levels, and the self-Hamiltonian reads $\hat{H}_{i}=\sum_{x}\text{{e}}_{i}^{(x)}|x\rangle_{i}\langle x|$ ($i=1,2$ for dipole-1,2), where $\text{{e}}_{i}^{(x)}$ and $|x\rangle_{i}$ are the eigen energy and the corresponding eigenstate.

The two magnetic dipoles interact with the EM field via the interaction (Wang et al., 2018; Weinberg, 2012)

[TABLE]

Here $\hat{\mathfrak{m}}_{i}$ is the magnetic dipole moment, $\mathbf{x}_{i}$ is the position of dipole-$i$. Under the Coulomb gauge, the quantized magnetic field $\hat{\mathbf{B}}\left(\mathbf{x},t\right)=\nabla\times\hat{\mathbf{A}}\left(\mathbf{x},t\right)$ reads as

[TABLE]

where $Z_{\mathbf{k}}:=\sqrt{\mu_{0}\hbar\omega_{\mathbf{k}}/2V}$ and $\hat{\mathrm{e}}_{\mathbf{k}\check{\sigma}}:=\hat{\mathrm{e}}_{\mathbf{k}}\times\hat{\mathrm{e}}_{\mathbf{k}\sigma}$. The index $\check{\sigma}$ refers to the polarization direction $\hat{\mathrm{e}}_{\mathbf{k}\check{\sigma}}$ orthogonal to $\hat{\mathrm{e}}_{\mathbf{k}\sigma}$.

Generally, the dipole moment operator $\hat{\mathfrak{m}}_{i}$ always can be expanded as $\hat{\mathfrak{m}}_{i}=\sum_{x,y}\vec{m}_{i}^{xy}\hat{\tau}_{i}^{xy}$, with $\vec{m}_{i}^{xy}:=\langle x|\hat{\mathfrak{m}}_{i}|y\rangle_{i}$ as the dipole moment amplitude, and $\hat{\tau}_{i}^{xy}:=|x\rangle_{i}\langle y|$ as the unitless transition operator.

The diagonal elements $\vec{m}_{i}^{xx}=\langle x|\hat{\mathfrak{m}}_{i}|x\rangle_{i}$ indicate the expectation value of the dipole moment on the level $|x\rangle_{i}$, thus they are permanent dipoles. The off-diagonal elements $\vec{m}_{i}^{xy}=\langle x|\hat{\mathfrak{m}}_{i}|y\rangle_{i}$ ($x\neq y$) are usually related to the photon radiation process, and they are transition dipoles (Wang et al., 2018).

When the magnetic dipoles are placed in the EM field, they change the original field dynamics. Then such influence propagates to each other, and that is how the interaction is generated between the two dipoles.

To make further calculation, we first consider that only dipole-1 is placed in the vacuum field, and study the field dynamics induced by dipole-1. With the above notations, such interaction $\hat{H}_{{\rm int}}^{(1)}=-\hat{\mathfrak{m}}_{1}(t)\cdot\hat{\mathbf{B}}(\mathbf{x}_{1},t)$ can be rewritten as

[TABLE]

Here the coefficients $g_{i,\mathbf{k}\sigma}^{xy}$ enclose the contributions from the EM field, the dipole moment $\vec{m}_{i}^{xy}$, and the position $\mathbf{x}_{i}$.

Under this interaction, the field dynamics is given by the Heisenberg equation $\frac{d}{dt}\hat{a}_{\mathbf{k}\sigma}=-\mathrm{i}\,\omega_{\mathbf{k}}\hat{a}_{\mathbf{k}\sigma}-\frac{\mathrm{i}}{\hbar}[\hat{a}_{\mathbf{k}\sigma},\hat{H}_{{\rm int}}^{(1)}]$, which gives

[TABLE]

Here the first term comes from the free field, and the second term is induced by the interaction with dipole-1.

Now we put this $\hat{a}_{\mathbf{k}\sigma}(t)$ back into the above magnetic field $\hat{\mathbf{B}}(\mathbf{x},t)$ [Eq. (3)], and it can be rewritten as $\hat{\mathbf{B}}(\mathbf{x},t)=\hat{\mathbf{B}}_{0}(\mathbf{x},t)+\hat{\mathbf{B}}_{\text{D1}}(\mathbf{x},t)$, where

[TABLE]

Here $\hat{\mathbf{B}}_{0}(\mathbf{x},t)$ is the original free field when there are no dipoles. The presence of dipole-1 changes the EM field, and gives rise to $\hat{\mathbf{B}}_{\text{D1}}(\mathbf{x},t)$, which we call as the dipole field. When this dipole field propagates to dipole-2, the dipole-dipole interaction is induced.

III Interaction propagated by the dipole field

Now we consider dipole-2 is placed in the field ($\hat{\mathfrak{m}}_{2}=\sum_{uv}\vec{m}_{2}^{uv}\hat{\tau}_{2}^{uv}$), which interacts with EM field locally via $\hat{H}_{{\rm int}}^{(2)}=-\hat{\mathfrak{m}}_{2}(t)\cdot\hat{\mathbf{B}}(\mathbf{x}_{2},t)$. Clearly, the dipole-dipole interaction is the part induced by the above dipole field $\hat{\mathbf{B}}_{\text{D1}}(\mathbf{x},t)$ [Eq. (5)], which can be rewritten as

[TABLE]

In the last line, we collect the mode summation of $\mathbf{k}\sigma$ into a convolution kernel $\mathcal{D}_{1\rightarrow 2}^{yx,uv}(s)$, which is defined as

[TABLE]

Here we define J_{1\rightarrow 2}^{yx,uv}(\omega):=\frac{2\pi}{\hbar}\sum_{\mathbf{k}\sigma}\big{(}g_{1,\mathbf{k}\sigma}^{xy}\big{)}^{*}g_{2,\mathbf{k}\sigma}^{uv}\delta(\omega-\omega_{\mathbf{k}}) as the coupling spectral density, which can be used to turn the above summation of $\mathbf{k}\sigma$** **into a continuous integral (Breuer and Petruccione, 2002; Wang et al., 2018). For the free-space EM field it can be explicitly calculated out as (see Appendix A)

[TABLE]

where we denote $\eta:=kr=\omega r/c$, $c$ is the velocity of light, and $r=|\mathbf{x}_{1}-\mathbf{x}_{2}|$ is the distance between the two dipoles. We can see $J_{1\rightarrow 2}^{yx,uv}(\omega)=-J_{1\rightarrow 2}^{yx,uv}(-\omega)$ is an odd function, and $J_{2\rightarrow 1}^{uv,yx}(\omega)=J_{1\rightarrow 2}^{yx,uv}(\omega)$. When $\omega\rightarrow 0$, we have $J_{1\rightarrow 2}^{yx,uv}(\omega)\rightarrow 0$.

The interaction (6) naturally has a retarded form, and up to now it is still an exact result without any approximation. Here the kernel function $\mathcal{D}_{1\rightarrow 2}^{yx,uv}(s)$ characterizes the propagating process of the interaction, and it is a fasting-decaying function of $t$, because its (half-)Fourier image $J_{1\rightarrow 2}^{yx,uv}(\omega)$ has a quite wide spectrum.

Since $\mathcal{D}_{1\rightarrow 2}^{yx,uv}(t)$ is already of the second order of the field-dipole interaction $g_{i,\mathbf{k}\sigma}^{xy,uv}$, thus we make an approximation (Lehmberg, 1970; Wang et al., 2018)

[TABLE]

where $\Omega_{1}^{yx}:=(\text{{e}}_{1}^{(y)}-\text{{e}}_{1}^{(x)})/\hbar$ is the transition frequency between the two energy levels $|y\rangle_{1}$ and $|x\rangle_{1}$. That means, its dynamics is governed only by the self-Hamiltonian $\hat{H}_{1}$, without considering its interaction with the field of higher orders. Then the interaction (6) would become local in time, which is composed of summation terms $\sim\hat{\tau}_{1}^{yx}(t)\cdot\hat{\tau}_{2}^{uv}(t)$ multiplying by time-independent interaction strengths.

Further, since $\mathcal{D}_{1\rightarrow 2}^{yx,uv}(t)$ decays very fast to zero, we can extend the time integral in Eq. (6) to $t\rightarrow\infty$, and that reduces the interaction Hamiltonian as

[TABLE]

where $\xi_{1\rightarrow 2}^{yx,uv}$ is calculated by111Here we used the relation $\int_{0}^{\infty}ds\,e^{\mathrm{i}\omega s}=\pi\delta(\omega)+\mathrm{i}\mathbf{P}\frac{1}{\omega}$, and $\mathbf{P}$ means the principal integral.

[TABLE]

By rearranging the summation indices, the interaction (6) can be rewritten as

[TABLE]

where K_{1\rightarrow 2}^{yx,uv}\big{(}\Omega_{1}^{yx}\big{)} is given by the function

[TABLE]

with the ratio $\eta=\Omega r/c$. Here we have utilized the relations $\Omega_{1}^{yx}=-\Omega_{1}^{xy}$, $J_{1\rightarrow 2}^{yx,uv}(-\omega)=-J_{1\rightarrow 2}^{yx,uv}(\omega)$, and $[J_{1\rightarrow 2}^{xy,vu}(\omega)]^{*}=J_{1\rightarrow 2}^{yx,uv}(\omega)$. Notice that $K_{1\rightarrow 2}^{yx,uv}(\Omega)=K_{1\rightarrow 2}^{yx,uv}(-\Omega)$ is an even function.

Following the same method as above, we should also consider the symmetric procedure that dipole-2 generates a dipole field $\hat{\mathbf{B}}_{\text{D2}}(\mathbf{x},t)$ and then interacts with dipole-1. That results to an interaction Hamiltonian $\hat{H}_{2\rightarrow 1}=\sum\varLambda_{2\rightarrow 1}^{uv,yx}\,\hat{\tau}_{2}^{uv}\hat{\tau}_{1}^{yx}$, which can be simply obtained by making proper variable exchanges from the above results. Finally, the total interaction between the two dipoles should be the average of the these two symmetric procedures $\hat{H}_{1\leftrightarrow 2}=\frac{1}{2}(\hat{H}_{1\rightarrow 2}+\hat{H}_{2\rightarrow 1})$, and we rewrite it as

[TABLE]

Here the interaction strength is divided into two parts: the principal term $G_{yx,uv}^{(\mathtt{P})}$ comes from the principal integral which is associated with virtual photon exchange, while the dissipation correction $G_{yx,uv}^{(\mathtt{D})}$ is associated with the real photon emission in the atom decay process (Li et al., 2015, 2016). Notice here $J_{i\rightarrow j}^{yx,uv}(\cdot)$ and $K_{i\rightarrow j}^{yx,uv}(\cdot)$ are not necessarily real numbers, which are determined by the vectors $\vec{m}_{i}^{yx}=\langle y|\hat{\mathfrak{m}}_{i}|x\rangle_{i}$.

Now we derived the general result for the field-induced interaction between two magnetic dipoles. No RWA was adopted, and there was no special requirement for the transition frequencies $\Omega_{i}^{yx}$. In the following, we will clarify the physical meaning of the different interaction terms in this result with the help of an example on a pair of two-level dipoles.

IV A pair of two-level dipoles

The above result (14) applies for general multi-state systems, and contains the interaction between any two matrix elements of the two dipole operators $\hat{\mathfrak{m}}_{i}=\sum\vec{m}_{i}^{xy}\hat{\tau}_{i}^{xy}$. Now we focus on an example that the two dipoles both have only two energy levels (Fig. 1), and discuss the different interaction types in the above general result (14) separately.

**Permanent dipoles: **

We first look at the interaction between the diagonal terms

[TABLE]

where $\hat{\tau}_{i}^{xx}=|x\rangle_{i}\langle x|$ for $x=\mathrm{e},\mathrm{g}$. In the interaction strength $G_{xx,yy}^{(...)}$, the dipole moment vector $\vec{m}_{i}^{xx}=\langle x|\hat{\mathfrak{m}}_{i}|x\rangle_{i}$ is the expectation value of the dipole operator $\hat{\mathfrak{m}}_{i}$ on the state $|x\rangle_{i}$, thus we say $\hat{\mathfrak{m}}_{i}^{xx}:=\vec{m}_{i}^{xx}\hat{\tau}_{i}^{xx}$ is the permanent dipole operator.

For these permanent dipole moments, the transition frequencies are $\Omega_{i}^{xx}=0$, thus the dissipation correction vanishes $G_{xx,yy}^{(\mathtt{D})}=0$, and the principal term $G_{xx,yy}^{(\mathtt{P})}$ directly returns the classical result [see $\Omega\rightarrow 0$ in Eqs. (13, 14)]. Finally, such interaction terms can be written as

[TABLE]

which has exactly the same form with the classical one.

**Resonant transition dipoles: **

Now we consider the interaction term between resonant transition dipoles. For the transition operators $\hat{\tau}_{1}^{+}:=|\mathrm{e}\rangle_{{}_{1}}\langle\mathrm{g}|$, $\hat{\tau}_{2}^{-}:=|\mathrm{g}\rangle_{{}_{2}}\langle\mathrm{e}|$, with resonant transition frequencies $\Omega_{1}^{\text{eg}}=-\Omega_{2}^{\text{ge}}:=\Omega$, the corresponding interaction term is

[TABLE]

From Eq. (14) we can see the two summation terms in $G_{\text{eg,ge}}^{(\mathtt{D})}$ exactly cancel each other [since $J_{2\rightarrow 1}^{uv,yx}(\omega)=J_{1\rightarrow 2}^{yx,uv}(\omega)$], thus the dissipation term always vanishes ($G_{\text{eg,ge}}^{(\mathtt{D})}=0$). The two terms in the principal term equal to each other, and that gives $G_{\text{eg,ge}}^{(\mathtt{P})}=K_{1\rightarrow 2}^{\text{eg,ge}}(\Omega)$. This result is the same as previous calculations based on master equations (Lehmberg, 1970; Agarwal, 1974; Ficek et al., 1987; Ficek and Swain, 2005; Wang et al., 2018). This principal interaction strength does not return to the classical form directly, but returns at the Dicke limit $r/\lambda\rightarrow 0$, where $\lambda$ is the wavelength of the transition frequency $\Omega$.

Following the same discussion, there is also a conjugate term for the transition dipole interaction $\hat{H}_{\text{ge,eg}}=G_{\text{ge,eg}}^{(\mathtt{P})}\,\hat{\tau}_{1}^{-}\hat{\tau}_{2}^{+}$, with $G_{\text{ge,eg}}^{(\mathtt{P})}=[G_{\text{eg,ge}}^{(\mathtt{P})}]^{*}$.

Non-resonant transition dipoles:

Further, we consider the case that frequencies of the two transition dipoles are different $\Omega_{1}\neq\Omega_{2}$. In this case, the above result (14) gives the principal term and the dissipation correction as

[TABLE]

Since $\Omega_{1}\neq\Omega_{2}$, the two terms in $G_{\text{eg,ge}}^{(\mathtt{D})}$ can no longer cancel each other completely. Notice that, when deriving the GKSL (Lindblad) master equation for an atomic system with coherent transitions, there is also such a similar unitary correction term like $G_{\text{eg,ge}}^{(\mathtt{D})}$ (see Eq. (12) in Ref. (Li et al., 2015) or Eq. (17) in Ref. (Li et al., 2016)), which comes from the dissipation terms of the master equation. Therefore, here we say $G_{xy,uv}^{(\mathtt{D})}$ is the dissipation correction. Similarly, the principal term becomes the average of $K_{1\rightarrow 2}^{\text{eg,ge}}(\Omega_{1})$ and $K_{1\rightarrow 2}^{\text{eg,ge}}(\Omega_{2})$.

From this result we can see the two dipoles emit both real and virtual photons, and they have interference with each other. But due to the frequency difference, they cannot achieve the complete interference as the above resonance case, thus the two terms in $G_{\text{eg,ge}}^{(\mathtt{D})}$ cannot cancel each other exactly.

**Counter-rotating terms and mixed terms: **

It is worth noticing that the above interaction (14) also contains counter-rotating terms, such as

[TABLE]

Here the principal interaction strength $G_{\text{eg,eg}}^{(\mathtt{P})}$ is quite similar with the above transition interaction (18) (except the ordering difference in the superscript). The dissipation term $G_{\text{eg,eg}}^{(\mathtt{D})}$ always exists for nonzero $\Omega_{i}$.

Moreover, indeed the interaction (14) also contains “mixed” interaction terms between the transition and permanent dipoles, e.g.,

[TABLE]

where $K_{1\rightarrow 2}^{\text{eg,gg}}(0)$ has the classical interaction form.

Usually, these counter-rotating and mixed interaction terms do not exhibit significant physical effects since they oscillate too fast ($\sim e^{\pm\mathrm{i}(\Omega_{1}+\Omega_{2})t}$, or $\sim e^{\pm\mathrm{i}\Omega_{i}t}$) comparing with the typical time scale determined by the interaction energy $G_{xy,uv}^{(...)}$, thus averagely they do not have significant contributions and can be omitted by the RWA. When the detuning between the two frequencies is quite large (e.g., $\Omega_{1}\gg\Omega_{2}$), we have $\Omega_{1}-\Omega_{2}\simeq\Omega_{1}+\Omega_{2}\simeq\Omega_{1}$, then these oscillating terms could be comparable with the non-resonant transition terms.

From this example of two-level dipoles, we can see the above general result (14) indeed contains different types of magnetic dipole-dipole interactions, and they are derived altogether based on the field propagation. Clearly, the generalization to multi-state systems is straightforward.

In the Dicke limit ($\Omega_{1,2}r/c\rightarrow 0$), we can see all the dissipation corrections approach zero, and all the principal terms approach $K_{1\rightarrow 2}^{yx,uv}(0)$. Thus, after summing up all the above terms, the interaction can be written as

[TABLE]

which just returns the classical form (1). In this sense, such field-induced interaction is consistent with that directly exchanged from the classical interaction. If the distance between the two dipoles is not short enough, the quantum corrections in the off-diagonal interaction terms (the transition terms, counter-rotating terms, and the mixed terms) should be considered.

V Summary

In this paper we made a general derivation for the magnetic dipole-dipole interaction based on the propagation of the quantized EM field. In this derivation, the two dipoles may have multiple energy levels, and we only adopted the second-order perturbation and the Markovian approximation, but did not require the resonance condition. Since the RWA is no longer needed, our result naturally gives the DD interaction for non-resonant dipoles, and even the counter-rotating interaction terms, as well as the mixed interaction terms between the permanent and transition dipoles. As long as the dipole operator $\hat{\mathfrak{m}}=\sum\vec{m}_{xy}|x\rangle\langle y|$ is written down for a specific physical system, all these interaction terms can be derived based on the field mediation.

From this result, we can notice that this field-induced DD interaction is associated with the inference of the virtual/real photons emitted from the two dipoles. As a result, the interaction strength could be influenced by the frequency difference of the two dipoles. If the two dipole moments has the same frequency, our result naturally returns to the result in previous studies. In principle, our derivation may also apply if the EM field has special properties (e.g., by applying certain periodic modulation or placing the system in a cavity with finite size (Wang et al., 2018; Baranov et al., 2016; Donaire et al., 2017; Tang et al., 2018; Davis et al., 2019)), or if the mediation field is not the EM field (e.g., the phonon field).

*Acknowledgements *- This study is supported by the National Natural Science Foundation of China (Grant No. 11675014), the Ministry of Science and Technology of China (2013YQ030595-3), and the BIT Research Fund Program for Young Scholars. S.-W. Li appreciated quite much for the helpful discussions with N. Wu and D. Xu in Beijing Institute of Technology.

Appendix A The coupling spectral density in free space

Here we show the derivation of the coupling spectral density $J_{1\rightarrow 2}^{yx,uv}(\omega)$ in free space, which is define by

[TABLE]

The summation over $\mathbf{k},\sigma$ is changed into integral by

[TABLE]

Thus the coupling spectral density $J_{1\rightarrow 2}^{yx,uv}(\omega)$ is given by (denoting $\mathbf{r}:=\mathbf{x}_{2}-\mathbf{x}_{1}$)

[TABLE]

Finally this integral gives (denoting $k:=\omega/c$, see Ref. (Wang et al., 2018))

[TABLE]

From this expression, we can see $J_{1\rightarrow 2}^{yx,uv}(-\omega)=-J_{1\rightarrow 2}^{yx,uv}(\omega)$ is an odd function of $\omega$. And we also have the following relation when exchanging the indices:

[TABLE]