Creation and protection of entanglement in systems out of thermal equilibrium

Quantum systems may present correlations of both quantum and classical nature. Entanglement captures quantum correlations due to the non-separability of the system state [1–3]. The presence of these correlations is connected to the rise of non-local effects in quantum theory [4] and has been recognized as a key resource in several fields of quantum technology, including quantum computing [5], quantum cryptography [6], quantum teleportation [7] and quantum metrology [8]. A main obstacle to the concrete exploitation of quantum features in the above applications is the detrimental effect of environmental noise [9]. The unavoidable coupling with degrees of freedom of the surrounding environment generally leads to a decay of quantum coherence properties [10], preventing the possible exploitation of quantum correlations present in the system.

Considerable effort has been made to understand the effects of environmental noise on the dynamics of correlations present in an open quantum system [11–14], and to contrast the natural fragility of quantum coherence properties [15–17]. Reservoir engineering methods have pointed out the possibility of changing the perspective from reducing the coupling with the environment to modifying the environmental properties in order to manipulate the system of interest, thanks to its proper dissipative dynamics [18–20]. Other approaches exploit the effect of measurements and feedback to drive the systems toward a target state [21, 22].

A possible way to create quantum correlations between two systems is to make them interact with a common environment [23], which can also cause a revival of entanglement [24]. In the case of two emitters in a common vacuum or thermal electromagnetic field, in the absence of matter close to them, the mediated interaction plays a role over distances of the order of the common transition wavelength [25, 26]. It has been found that the presence of plasmonic waveguides near the emitters can allow a mediated interaction over larger distances [27] whose effect on the entanglement dynamics has been discussed [28]. However, at thermal equilibrium, the dynamic creation of entanglement eventually ceases at some time and the system thermalizes toward a thermal state which is a classical mixture. Steady entanglement can be instead generated by adding the action of an external driving laser [29].

The influence of several independent reservoirs at different temperatures, whose emission does not depend on their internal structure (material or geometry), has been considered in several contexts, including generation of entanglement in non-equilibrium steady states, both in the case of a few spins [30–32] and of a chain of spins [33–35] and in the context of quantum thermal machines [36, 37].

However, in a realistic configuration, the actual reflection and transmission properties of the bodies surrounding the quantum emitters should be taken into account, and may become particularly relevant if the emitters are placed close to the bodies (near-field effects). New possibilities emerging in such realistic systems out of thermal equilibrium have been recently pointed out in different contexts ranging from heat transfer [38, 39], to Casimir–Lifshitz forces [40–45]. There, radiation fields out of thermal equilibrium in configurations of a quite general nature have been characterized in terms of the correlators of the total field depending on the scattering matrices of the bodies composing the total system [46, 47]. In the case of single emitters in such environments, new tools exploiting the absence of thermal equilibrium to manipulate the atomic dynamics realizing inversion of population and cooling of internal atomic temperature have been pointed out [48, 49]. Recently, the case of two quantum emitters has also been analyzed, pointing out a new remarkable mechanism to generate and protect entanglement in a steady way in systems out of thermal equilibrium [50].

In this paper, we report a detailed investigation of this phenomenon by studying the internal dynamics of a system composed of two quantum emitters (real atoms or artificial ones such as quantum dots or superconducting qudits) placed in front of an arbitrary body embedded in a thermal radiation whose temperature is different from that of the body. The paper is organized as follows. In section 2, we describe the physical model under investigation and we derive a master equation for the general case of two N-level emitters. In section 3, we derive closed-form expressions for the functions governing the dynamics, in terms of the scattering matrices of the body and valid for arbitrary geometrical and material properties. In section 4, we develop these expressions in the case when the body is a slab of finite thickness. From section 5 onwards we specialize our analysis to the case of a two-qubit system, comparing cases in and out of thermal equilibrium. We point out the occurrence of peculiar phenomena emerging out of thermal equilibrium such as the generation of steady entanglement and a simple interpretation for this phenomenon is presented for a particularly interesting case. The general case of arbitrary values of the parameters is then discussed in section 6. In section 7, we draw our conclusions.

We consider a system made of two quantum emitters q = 1,2 interacting with an environment consisting of an electromagnetic field which is stationary and out of thermal equilibrium. This is generated by the field emitted by a body (M) at temperature T_M of arbitrary geometry and dielectric response and by the field emitted by the far surrounding walls (W) at temperature T_W, which is eventually transmitted and reflected by the body itself (see figure 1).

Zoom In Zoom Out Reset image size

T_M and T_W are kept fixed in time, realizing a stationary configuration for the electromagnetic field. The surrounding walls have an irregular shape and are distant enough from the body and the emitters so that their field can be treated at the emitters' locations, in the absence of the body, as blackbody radiation independent of their composition. This is not true for the field emitted by the body M which cannot be treated as a blackbody since its radiation depends on its actual properties such as its geometry and its dielectric function. The total Hamiltonian has the form

$$\begin{equation} H=H_{\rm S}+H_{\rm E}+H_{\rm I}, \end{equation} \tag{ 1 }$$

where H_S and H_E are the free Hamiltonians of the two emitters and of the environment. The interaction between the emitters and the field, in the multipolar coupling and in dipole approximation, is [51]

$$\begin{equation} H_{\rm I}=-\sum_{q} \mathbf{D}_q\cdot\mathbf{E}(\mathbf{R}_q) , \end{equation} \tag{ 2 }$$

where D_q is the electric-dipole operator of emitter q and E(R_q) is the electric field at its position R_q.

We first consider the general case in which each emitter has N_q internal levels |n〉_q where n∈{1,...,N_q} of frequency ω^q_n (ordered by increasing energy). Given two arbitrary levels n and m, their frequency difference is indicated by ω^q_nm = ω^q_n − ω^q_m and the transition matrix element of the dipole operator by d^q_mn = _q〈m|D_q|n〉_q. The free Hamiltonian of the two emitters is

$$\begin{equation} H_{\rm S}=\sum_{q=1}^2 H_q=\sum_{q, \epsilon_n^q} \epsilon_n^q \Pi(\epsilon_n^q), \end{equation} \tag{ 3 }$$

where $\Pi(\epsilon_n^q)=|n\rangle_{qq}\langle n|$ are the projectors associated with each eigenvalue ^q_n = ℏω^q_n (possibly degenerate) of H_q. The dipole operator of emitter q in the interaction picture, $\mathbf {D}_q(t)=\exp (\frac {\mathrm {i}}{\hbar }H_{\mathrm { S}}t)\mathbf {D}_q\exp (-\frac {\mathrm {i}}{\hbar }H_{\mathrm { S}}t)$, results in

$$\begin{equation} \mathbf{D}_q(t)=\sum_{^{\,m,n}_{n>m}}(\mathbf{d}_{m n}^q\ \sigma_{mn}^q\,{\rm e}^{-{\rm i}\omega_{nm}^qt}+\mathrm{h.c.}), \end{equation} \tag{ 4 }$$

where σ^q_mn = |m〉_qq〈n| and ω^q_nm ⩾ 0. By moving to the interaction picture, we obtain for H_I

$$\begin{equation} H_{\rm I}(t)=-\sum_{q} \mathbf{D}_q(t)\cdot\mathbf{E}(\mathbf{R}_q,t), \end{equation} \tag{ 5 }$$

where the time-dependent electric field is given by $\mathbf {E}(\mathbf {R}_q,t)=\exp (\frac {\mathrm {i}}{\hbar }H_{\mathrm { E}}t)E(\mathbf {R}_q)\exp (-\frac {\mathrm {i}}{\hbar }H_{\mathrm { E}}t)$. In the following each mode of the field is identified by the frequency ω, the transverse wave vector k = (k_x,k_y), the polarization index p (taking the values p = 1,2 corresponding to transverse-electric (TE) and transverse-magnetic (TM) polarizations respectively) and the direction of propagation ϕ = ± 1 (shorthand notation ϕ = ±) along the z-axis (see figure 1). In this approach, the total wavevector takes the form K^ϕ = (k,ϕk_z), where the z component of the wavevector k_z is a dependent variable given by $k_z=\sqrt {\frac {\omega ^2}{c^2}-k^2}$, where k = |k|. The explicit expression of the field at an arbitrary point R is

$$\begin{equation} \mathbf{E}(\mathbf{R},t)=2\mathrm{Re}\Bigg[\int_0^{+\infty}\frac{\mathrm{d}\omega}{2\pi}\mathrm{e}^{-\mathrm{i}\omega t}\mathbf{E}(\mathbf{R},\omega)\Bigg], \end{equation} \tag{ 6 }$$

where a single-frequency component reads

$$\begin{equation} \mathbf{E}(\mathbf{R},\omega)=\sum_{\phi,p}\int\frac{\mathrm{d}^2\mathbf{k}}{(2\pi)^2}{\rm e}^{{\rm i}\mathbf{K}^\phi\cdot\mathbf{R}}\hat{\boldsymbol{\epsilon}}^\phi_{p}(\mathbf{k},\omega)E^\phi_p(\mathbf{k},\omega), \end{equation} \tag{ 7 }$$

where E^ϕ_p(k,ω) is the field amplitude operator associated with the mode (ω,k,p,ϕ). For the TE and TM polarization vectors appearing in (7), we adopt the following standard definitions:

$$\begin{equation} \begin{array}{ll} \displaystyle \hat{\boldsymbol{\epsilon}}^\phi_{\mathrm{TE}}(\mathbf{k},\omega)=\hat{\mathbf{z}}\times\hat{\mathbf{k}}=\frac{1}{k}(-k_y\hat{\mathbf{x}}+k_x\hat{\mathbf{y}}),\\ \displaystyle \hat{\boldsymbol{\epsilon}}^\phi_{\mathrm{TM}}(\mathbf{k},\omega)=\frac{c}{\omega}\hat{\boldsymbol{\epsilon}}^\phi_{\mathrm{TE}}(\mathbf{k},\omega)\times\mathbf{K}^{\phi}=\frac{c}{\omega}(-k\hat{\mathbf{z}}+\phi k_z\hat{\mathbf{k}}), \end{array} \end{equation} \tag{ 8 }$$

where $\hat {\mathbf {x}}$, $\hat {\mathbf {y}}$ and $\hat {\mathbf {z}}$ are the unit vectors along the three axes and $\hat {\mathbf {k}}=\mathbf {k}/k$.

2.1. Master equation

The starting point to study the dynamics of the two emitters is, in the interaction picture, the von Neumann equation for the total density matrix ρ_tot(t):

$$\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\rho_{\mathrm{tot}}(t)=-\frac{\mathrm{i}}{\hbar}[H_{\rm I}(t),\rho_{\mathrm{tot}}(t)]. \end{equation} \tag{ 9 }$$

The reduced density matrix of the two emitters is given by ρ = Tr_E[ρ_tot], where Tr_E denotes the trace over the degrees of freedom of the environment. To derive a master equation for ρ we follow the procedure described in [9] for the case of one emitter by extending it to our system made of two emitters. We name ω_q an arbitrary transition frequency of emitter q (positive and negative). In general, several transitions can be characterized by the same frequency ω_q both because of degeneracy and/or the occurrence of equidistant levels. We rewrite each Cartesian component of the dipole operator, [D_q]_i (i = {x,y,z}), as

$$\begin{equation} [\mathbf{D}_q]_i= \sum_{\omega_q} \sum_{^{\:\:\:\:\{\epsilon^{q}_{n},\epsilon_m^q\}}_{\epsilon^{q}_{n}-\epsilon_m^q=\hbar \omega_q}} \Pi(\epsilon_m^q)[\mathbf{D}_q]_i\Pi(\epsilon^{q}_{n})=\sum_{\omega_q} A^q_i(\omega_q), \end{equation} \tag{ 10 }$$

where A^q_i(ω_q) and A^q †_i(ω_q) turn out to be the eigenoperators of H_q with frequencies −ω_q and +ω_q, respectively, i.e. [H_q,A^q_i(ω_q)] = −ω_qA^q_i(ω_q) and [H_q,A^q †_i(ω_q)] = + ω_qA^q †_i(ω_q). It also holds A^q †_i(ω_q) = A^q_i(− ω_q) and $\exp (\frac {\mathrm {i}}{\hbar }H_{\mathrm { S}}t) A^q_i(\omega _q)\exp (-\frac {\mathrm {i}}{\hbar }H_{\mathrm { S}}t)={\mathrm { e}}^{-{\mathrm { i}}\omega _q t}A^q_{i}(\omega _q)$. In the central term of (10), the first sum is over all of the frequencies ω_q while the second is over all of the couples of energy eigenvalues ^q_n and ^q_m of H_q such that ^q_n − ^q_m = ℏω_q. Following [9], it is useful to rewrite H_I(t) of (5) in terms of the eigenoperators A^q_i(ω_q) as

$$\begin{equation} H_{\rm I}(t)=-\sum_q\sum_{i,\omega_q}{\rm e}^{-{\rm i}\omega_q t}A^q_{i}(\omega_q) E_i(\mathbf{R}_q,t). \end{equation} \tag{ 11 }$$

From (10) it follows that the vector A^q(ω_q) = {A^q_x(ω_q),A^q_y(ω_q),A^q_z(ω_q)} is given by

$$\begin{equation} \mathbf{A}^q(\omega_q)= \sum_{^{\:\:\:\{m,n\}}_{\omega_{nm}^q=\omega_q}} \mathbf{d}_{mn}^q \sigma_{mn}^q=\mathbf{A}^{q\,\dagger}(-\omega_q), \end{equation} \tag{ 12 }$$

where the sum is over all of the couples n and m such that ω^q_nm = ω_q. By applying to the case of two emitters the standard procedure for the microscopic derivation of a master equation reported in [9], under Born, Markovian and rotating-wave approximations⁴, one can obtain (by using also the condition 〈E_i(R,t)〉 = 0) in the Schrödinger representation

$$\begin{eqnarray} &&\fl \frac{\mathrm{d}}{\mathrm{d} t}\rho =-\frac{\mathrm{i}}{\hbar} [H_{\rm S},\rho ]-\mathrm{i} \sum_{q, q', \omega}\sum_{i,i'} \{s_{ii'}^{qq'}(\omega) [A^{q\,\dagger}_{i}(\omega)A_{i'}^{q'}(\omega),\rho ] \nonumber \\ &&+\gamma_{ii'}^{qq'}(\omega) (A^{q'}_{i'}(\omega)\rho A^{q\,\dagger}_{i}(\omega) -\frac{1}{2}\{A^{q\,\dagger}_{i}(\omega)A^{q'}_{i'}(\omega),\rho \})\}, \end{eqnarray} \tag{ 13 }$$

where $\omega \gtreqless 0$, being terms with positive or negative ω associated, respectively, with downward and upward transitions. In the above equation, for q ≠ q' the sum $\sum _{q, q', \omega }$ is over all of the common frequencies ω_q = ω_q' = ω (this condition derives from the rotating wave approximation) while for q = q' it is over all of the transition frequencies of each emitter, and $\gamma _{ii'}^{qq'}(\omega )$, and $s_{ii'}^{qq'}(\omega )$ are defined by

$$\begin{equation} \begin{array}{ll} \displaystyle \gamma_{ii'}^{qq'}(\omega) =\Xi_{ii'}^{qq'}(\omega) +\Xi_{i' i }^{q' q \,*}(\omega) , \quad s_{ii'}^{qq'}(\omega) = \frac{\Xi_{ii'}^{qq'}(\omega) -\Xi_{i' i }^{q' q \,*}(\omega)}{2 i}, \\ \displaystyle \Xi_{ii'}^{qq'}(\omega)= \frac{1}{\hbar^2}\int_0^\infty \mathrm{d}s\,{\rm e}^{{\rm i}\omega s}\langle E_i(\mathbf{R}_q,t)E_{i'}(\mathbf{R}_{q'},t-s)\rangle , \end{array} \end{equation} \tag{ 14 }$$

where the field correlation functions enter in the function $\Xi _{ii'}^{qq'}(\omega )$. It follows that $\Xi _{ii'}^{qq'}(\omega )=\frac {1}{2}\gamma _{ii'}^{qq'}(\omega )+\mathrm {i} s_{ii'}^{qq'}(\omega )$, $[\gamma _{ii'}^{qq'}(\omega )]^{*}=\gamma _{i'i}^{q'q}(\omega )$ and $[s_{ii'}^{qq'}(\omega )]^{*}=s_{i'i}^{q'q}(\omega )$.

The initial state of the total system in (13) is assumed to be factorized, ρ_tot(0) = ρ(0)ρ_E. In the case ρ_E is a stationary state of the environment ([H_E,ρ_E] = 0) the correlation functions are homogeneous in time, that is 〈E_i(R_q,t)E_i'(R_q',t − s)〉 = 〈E_i(R_q,s)E_i'(R_q',0)〉, so that

$$\begin{equation} \gamma_{ii'}^{qq'}(\omega)=\frac{1}{\hbar^2}\int_{-\infty}^\infty \mathrm{d}s\,{\rm e}^{{\rm i}\omega s}\langle E_i(\mathbf{R}_q,s)E_{i'}(\mathbf{R}_{q'},0)\rangle \end{equation} \tag{ 15 }$$

does not depend on time. The functions defined in (14) appearing in the master equation (13) depend thus only on the field correlation functions 〈E_i(R_q,s)E_i'(R_q',0)〉, whose computation out of thermal equilibrium will be the subject of sections 3 and 4.

We now explicitly write the master equation (13) in the case of absence of degenerate and equidistant levels in each emitter, when the definition of the eigenoperators A^q(ω_q) (12) reduces to A^q(ω_q) = d^q_mnσ^q_mn (each ω_q corresponds to only one couple of energy eigenvalues {^q_m,^q_n}). For this purpose, we develop the sum over ω in (13), which for each |ω| runs over ω (downward transitions) and −ω (upward transitions), as $\sum _\omega f(\omega )=\sum _{\omega >0} f(\omega ) + \sum _{\omega >0} f(-\omega )$. From now on ω always indicates a positive frequency and we drop ' > 0' in the sums over ω. By introducing this new convention and using the explicit form for A^q_i(ω) (12), we can recast  (13) as

$$\begin{eqnarray} &&\fl \frac{\mathrm{d}}{\mathrm{d} t}\rho =-\frac{\mathrm{i}}{\hbar}[H_{\rm S},\rho ]-\mathrm{i} \sum_{q, \omega} \{S^{qq}(\omega)[\sigma_{nn}^q,\rho ]+S^{qq}(-\omega)[\sigma_{mm}^q,\rho ]\}\nonumber\\ &&-\mathrm{i}\sum_{q\neq q', \omega}\Lambda^{qq'}(\omega)[\sigma_{mn}^{q\,\dagger}\sigma_{m'n'}^{q'},\rho ] +\sum_{q, q', \omega}\left\{\vphantom{\frac{1}{2}}\Gamma^{qq'}(\omega)(\sigma^{q'}_{m'n'}\rho \sigma_{mn}^{q\,\dagger}\right.\nonumber\\ &&\left.-\frac{1}{2}\{\sigma^{q\,\dagger}_{mn}\sigma^{q'}_{m'n'},\rho \})+\Gamma^{qq'}(-\omega)(\sigma^{q'\dagger}_{m'n'}\rho \sigma^{q}_{mn}-\frac{1}{2}\{\sigma^{q}_{mn} \sigma^{q'\,\dagger}_{m'n'},\rho \})\right\} , \end{eqnarray} \tag{ 16 }$$

where the sum $\sum _{q, q', \omega }$ in the second line is relative to all of the transition frequencies of each emitter for q = q' and only to the common transition frequencies for q ≠ q', (m,n) and (m',n') individuate, respectively, the transition of each emitter corresponding to the frequency ω, and we have defined the functions

$$\begin{eqnarray} &&\fl S^{qq}(\omega)=\sum_{i,i'} s_{ii'}^{qq}(\omega)[{\bf d}_{mn}^q]^{*}_{i} [{\bf d}_{mn}^{q}]_{i'},\quad S^{qq}(-\omega)=\sum_{i,i'}s_{ii'}^{qq}(-\omega)[{\bf d}_{mn}^q]_{i}[{\bf d}_{mn}^q]^{*}_{i'},\nonumber\\ &&\fl \Lambda^{qq'}(\omega)= \sum_{i,i'} [\mathbf{d}_{mn}^q]_i^{*} [\mathbf{d}_{m'n'}^{q'}]_{i'} [s_{ii'}^{qq'}(\omega)+s_{i'i}^{q' q}(-\omega)],\\ &&\fl \Gamma^{qq'}(\omega)= \sum_{i,i'}\gamma_{i i'}^{qq'}(\omega)[{\bf d}_{mn}^q]^{*}_{i} [{\bf d}_{m'n'}^{q'}]_{i'},\quad \Gamma^{qq'}(-\omega)=\sum_{i,i'}\gamma_{i i'}^{qq'}(-\omega)[{\bf d}_{mn}^q]_{i}[{\bf d}_{m'n'}^{q'}]^{*}_{i'}.\nonumber \end{eqnarray} \tag{ 17 }$$

We remark that it holds [Λ^qq'(ω)]* = Λ^q'q(ω) and [Γ^qq'(ω)]* = Γ^q'q(ω). In  (16), function Λ^qq'(ω) represents a coherent (dipole–dipole) interaction between the emitters mediated by the field while dissipative effects enter through the Γ functions. In particular, Γ^qq'(± ω) are individual (q = q') and common field-mediated collective (q ≠ q') emitter transition rates, related to both quantum and thermal fluctuations of the electromagnetic field at the emitters' position.

Here, we derive the field correlation functions needed to compute the functions in (17) for non-equilibrium configurations in the case of an arbitrary body and multilevel emitters. These functions will depend on the two temperatures T_M and T_W and on the material and geometrical properties of the body as well. We follow the derivation discussed in [47] in the more general case of two bodies and three temperatures and the derivation relative to a single quantum emitter in the presence of a single body and two temperatures [49]. Here, we extend the latter derivation to the case of two quantum emitters. Some of the computations involved are reported in appendix A.

The starting point is to decompose on the right side of the body where the emitters are located, the amplitude operators of total field modes propagating in the two directions z > 0 and z < 0 in terms of the fields emitted by the surrounding walls (W) and by the body (M). For a given set (ω,k,p), we have for the two directions

$$\begin{equation} E^+=E^{(\mathrm{M})+}+\mathcal{T}E^{\mathrm{(W)}+}+\mathcal{R}E^{\mathrm{\mathrm{(W)}}-}, \quad E^-=E^{\mathrm{(W)}-}, \end{equation} \tag{ 18 }$$

where we made the dependence on ω, k and p implicit. The total field E⁻ propagating toward the body (i.e. toward the left) is equal to the field emitted by the walls E^(W)− coming from the left, while the total field E⁺ propagating toward the right results from the field E^(M)+ directly produced by the body, the transmission through the body of the field E^(W)+ emitted by the walls coming from the left and the reflection by the body of the field E^(W)− coming from the right (see figure 2).

Zoom In Zoom Out Reset image size

The operators $\mathcal {R}$ and $\mathcal {T}$ are the reflection and transmission scattering operators associated with the right side of the body, whose explicit definition can be found for example in [47]. They connect any outgoing (reflected or transmitted) mode of the field to the entire set of incoming modes. By using (18) one can write the total field correlators in terms of the correlators of the fields emitted by each source.

The source fields have been characterized as in [47] by assuming that for the body M and the walls W a local temperature that remains constant in time can be defined and that the emission process of the body is essentially not influenced by the presence of the external radiation impinging on the body itself. This assumption leads to the hypothesis that the part of the total field emitted by the body is the same as it would be if the body were at thermal equilibrium with the environment at its own temperature so that the correlators of the field emitted by each body can still be deduced by using the fluctuation–dissipation theorem at its local temperature.

Under this assumption, the following symmetrized correlation functions [〈AB〉_sym = (〈AB〉 + 〈BA〉)/2] have been derived

$$\begin{eqnarray} &&\fl \langle E^{\mathrm{(M)}+}_p(\mathbf{k},\omega)E^{\mathrm{(M)}+\dagger}_{p'}(\mathbf{k}',\omega')\rangle_{\mathrm{sym}}=\frac{\omega}{2\epsilon_0c^2}N(\omega,T_{\mathrm{M}})2\pi\delta(\omega-\omega')\langle p,\mathbf{k}|\left(\mathcal{P}_{-1}^{\mathrm{(pw)}}\right.\nonumber\\ &&\left.\, -\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger} +\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(ew)}} -\mathcal{P}_{-1}^{\mathrm{(ew)}}\mathcal{R}^{\dagger}-\mathcal{T}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{T}^{\dagger}\right)|p',\mathbf{k}'\rangle, \nonumber\\ &&\fl \langle E^{\mathrm{(W)}\phi}_p(\mathbf{k},\omega)E^{\mathrm{(W)}\phi'\dagger}_{p'}(\mathbf{k}',\omega')\rangle_{\mathrm{sym}} =\frac{\omega}{2\epsilon_0c^2}N(\omega,T_{\mathrm{W}}) 2\pi\delta(\omega-\omega')\delta_{\phi,\phi'}\langle p,\mathbf{k}|\mathcal{P}_{-1}^{\mathrm{(pw)}}|p',\mathbf{k}'\rangle, \end{eqnarray} \tag{ 19 }$$

where we have introduced

$$\begin{equation} \begin{array}{ll} \hspace{-3pc}\displaystyle N(\omega,T)=\frac{\hbar\omega}{2}\coth\left(\frac{\hbar\omega}{2k_{\rm B}T}\right)=\hbar\omega\left[\frac{1}{2}+n(\omega,T)\right],\quad n(\omega,T)=(\mathrm{e}^{\frac{\hbar\omega}{k_{\rm B}T}}-1)^{-1},\\ \hspace{-3pc}\displaystyle \langle p,\mathbf{k}|\mathcal{P}_n^{\mathrm{(pw/ew)}}|p',\mathbf{k}'\rangle=k_z^n\langle p,\mathbf{k}|\Pi^{\mathrm{(pw/ew)}}|p',\mathbf{k}'\rangle. \end{array} \end{equation} \tag{ 20 }$$

In the above equation Π^(pw) and Π^(ew) are the projectors on the propagative (c k < ω, corresponding to a real k_z) and evanescent (c k > ω, corresponding to a purely imaginary k_z) sectors respectively. By combining (18) and (19), in appendix A a general expression for the total correlation functions in frequency space has been derived in (A.7). This expression can be used to compute the functions $\gamma _{ii'}^{qq'}(\omega), \gamma_{ii'}^{qq'}(-\omega)$ and $s_{ii'}^{qq'}(\omega )$ entering in (17), by exploiting their connection with the correlation functions between frequency components of the total field given in (A.8).

To move to the final expression of the functions in (17), we first rewrite the antinormally ordered correlation functions (A.7) as

$$\begin{eqnarray} &&\fl \langle E_i (\mathbf{R}_q,\omega) E_{i'}^\dagger(\mathbf{R}_{q'},\omega)\rangle=\frac{\hbar \omega^3}{3 \pi\epsilon_0 c^3}\{[1+n(\omega,T_{\mathrm{W}})] [\alpha_{\mathrm{W}}^{qq'}(\omega)]_{ii'}+ [1+n(\omega,T_{\mathrm{M}})][\alpha_{\mathrm{M}}^{qq'}(\omega)]_{ii'} \},\nonumber\\ \end{eqnarray} \tag{ 21 }$$

from which the normally ordered correlation functions are obtained by replacing [1 + n(ω,T_i)] with n(ω,T_i) and by taking the complex conjugate (this procedure derives from Kubo's prescription as explained in appendix A)

$$\begin{equation} \fl \langle E_i^\dagger(\mathbf{R}_q,\omega) E_{i'}(\mathbf{R}_{q'},\omega)\rangle=\frac{\hbar \omega^3}{3 \pi\epsilon_0 c^3}\{ n(\omega,T_{\mathrm{W}} [\alpha_{\mathrm{W}}^{qq'}(\omega)]_{ii'}^{*}+ n(\omega,T_{\mathrm{M}})][\alpha_{\mathrm{M}}^{qq'}(\omega)]_{ii'}^{*} \}, \end{equation} \tag{ 22 }$$

and where we have introduced two α functions which do not depend on the temperatures and on the dipoles, and do depend on the geometrical and material properties of the body through the operators $\mathcal {R}$ and $\mathcal {T}$:

$$\begin{eqnarray} &&\fl [\alpha_{\mathrm{W}}^{qq'}(\omega)]_{ii'}=\frac{3\pi c}{2 \omega}\sum_{p,p'}\int\frac{\mathrm{d}^2\mathbf{k}}{(2\pi)^2}\int\frac{\mathrm{d}^2\mathbf{k}'}{(2\pi)^2}{\rm e}^{{\rm i}(\mathbf{k}\cdot\mathbf{r}_q-\mathbf{k}'\cdot\mathbf{r}_{q'})} \langle p,\mathbf{k}|\{{\rm e}^{{\rm i}(k_z z_q-k_z^{'*}z_{q'})}\nonumber\\ &&\, \times [\hat{\boldsymbol{\epsilon}}_p^+(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}_{p'}^{+}(\mathbf{k}',\omega)]_{i'}^{*}(\mathcal{T}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{T}^{\dagger}+\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger})\nonumber\\ &&\,+{\rm e}^{{\rm i}(k_zz_q+k_z^{'*}z_{q'})}[\hat{\boldsymbol{\epsilon}}_p^+(\mathbf{k},\omega)]_i [\hat{\boldsymbol{\epsilon}}_{p'}^{-}(\mathbf{k}',\omega)]_{i'}^{*}\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}}\nonumber\\ &&\,+{\rm e}^{-{\rm i}(k_zz_q+k_z^{'*}z_{q'})}[\hat{\boldsymbol{\epsilon}}_p^-(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}_{p'}^{+}(\mathbf{k}',\omega)]_{i'}^{*}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger}\nonumber\\ &&+{\rm e}^{-{\rm i}(k_z z_{q}-k_z^{'*}z_{q'})}[\hat{\boldsymbol{\epsilon}}_p^-(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}_{p'}^{-}(\mathbf{k}',\omega)]_{i'}^{*}\mathcal{P}_{-1}^{\mathrm{(pw)}}\}|p',\mathbf{k}'\rangle, \nonumber\\ &&\fl [\alpha_{\mathrm{M}}^{qq'}(\omega)]_{ii'} =\frac{3\pi c}{2 \omega}\sum_{p,p'}\int\frac{\mathrm{d}^2\mathbf{k}}{(2\pi)^2}\int\frac{\mathrm{d}^2\mathbf{k}'}{(2\pi)^2}{\rm e}^{{\rm i}(\mathbf{k} \cdot \mathbf{r}_q-\mathbf{k}'\cdot\mathbf{r}_{q'})}\langle p,\mathbf{k}|\{ {\rm e}^{{\rm i}(k_z z_q-k_z^{'*}z_{q'})}\nonumber\\ &&\,\times [\hat{\boldsymbol{\epsilon}}_p^+(\mathbf{k},\omega)]_i [\hat{\boldsymbol{\epsilon}}_{p'}^{+}(\mathbf{k}',\omega)]_{i'}^{*} [(\mathcal{P}_{-1}^{\mathrm{(pw)}}-\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger}+\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(ew)}}\nonumber\\ &&\,-\mathcal{P}_{-1}^{\mathrm{(ew)}}\mathcal{R}^{\dagger}-\mathcal{T}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{T}^{\dagger})\,\}|p',\mathbf{k}'\rangle. \end{eqnarray} \tag{ 23 }$$

Functions [α^qq'_W(ω)]_ii' and [α^qq'_M(ω)]_ii' are real for q = q', while in general they are complex for q ≠ q' satisfying Im[α^qq'_W(ω)]_ii' = −Im[α^qq'_M(ω)]_ii'. The last property ensures that the function [α^qq'_W(ω)]_ii' + [α^qq'_M(ω)]_ii' is real as expected, being proportional to the imaginary part of the Green function (see (C.3)).

Now we can compute the transition rates in (17), by using (A.8), (21) and (22),

$$\begin{equation} \begin{array}{ll} \hspace{-3pc}\displaystyle \Gamma^{qq'}(\omega) = \sqrt{\Gamma_0^q(\omega)\Gamma_0^{q'}(\omega) } \{[1+n(\omega,T_{\mathrm{W}})]\alpha_{\mathrm{W}}^{qq'}(\omega) + [1+n(\omega,T_{\mathrm{M}})]\alpha_{\mathrm{M}}^{qq'}(\omega)\},\\ \hspace{-3pc}\displaystyle \Gamma^{qq'}(-\omega) = \sqrt{\Gamma_0^q(\omega)\Gamma_0^{q'}(\omega) } [n(\omega,T_{\mathrm{W}}) \alpha_{\mathrm{W}}^{qq'}(\omega)^{*} +n(\omega,T_{\mathrm{M}})\alpha_{\mathrm{M}}^{qq'}(\omega)^{*}], \end{array} \end{equation} \tag{ 24 }$$

where $\Gamma _0^q(\omega )=\frac {|{\bf d}_{mn}^q|^2\omega ^3}{3 \hbar \pi \epsilon _0 c^3}$ is the vacuum spontaneous-emission rate of transition |n〉_q → |m〉_q of emitter q and we have introduced the new functions

$$\begin{equation} \fl \alpha_{\mathrm{W}}^{qq'}(\omega) =\sum_{i,i'}[\tilde{{\bf d}}_{mn}^q]^{*}_{i} [\tilde{{\bf d}}_{m'n'}^{q'}]_{i'} [\alpha_{\mathrm{W}}^{qq'}(\omega)]_{ii'}, \quad \alpha_{\mathrm{M}}^{qq'}(\omega) =\sum_{i,i'}[\tilde{{\bf d}}_{mn}^q]^{*}_{i} [\tilde{{\bf d}}_{m'n'}^{q'}]_{i'} [\alpha_{\mathrm{M}}^{qq'}(\omega)]_{ii'}, \end{equation} \tag{ 25 }$$

being $[\tilde {{\bf d}}_{mn}^q]_{i}=[{\bf d}_{mn}^q]_{i}/|{\bf d}_{mn}^q|$. Differently from [α^qq'_W(M)(ω)]_ii', the functions α^qq'_W(M)(ω) depend on the choice of emitters' dipoles. In the case of two qubits, which will be treated in sections 5 and 6, there is only one transition for each emitter and the above equations (24) and (25) hold with the notation d^q_mn = d^q.

With regards to the function Λ^qq'(ω) we obtain, by using (A.8), (21) and (22),

$$\begin{equation} \fl \Lambda^{qq'}(\omega)= \frac{\sqrt{\Gamma_0^q(\omega)\Gamma_0^{q'}(\omega) }}{\omega^3} \sum_{i,i'} [\tilde{\mathbf{d}}_{mn}^q]_i^{*}[\tilde{\mathbf{d}}_{m'n'}^{q'}]_{i'} \mathcal{P}\int_{-\infty}^{+\infty}\frac{\omega'^3 \,\mathrm{d}\omega'}{2\pi}\frac{[\alpha_{\mathrm{W}}^{qq'}(\omega')]_{ii'} +[\alpha_{\mathrm{M}}^{qq'}(\omega')]_{ii'} }{\omega-\omega'} , \end{equation} \tag{ 26 }$$

where we used the properties [α^qq'_W(M)(ω')]_ii' = [α^q'q_W(M)(ω')]*_i'i and [α^qq'_W(M)(− ω')]_ii' = [α^qq'_W(M)(ω')]*_ii'. It follows that Λ^qq'(ω) does not depend on the presence or the absence of thermal equilibrium, being independent of the temperatures. By using the relation between α functions of (23) and the Green's function of the system in (C.3) derived in appendix C, the integration over frequencies in  (26) can be achieved by using the Kramers–Kronig relations connecting real and imaginary parts of Green's function:

$$\begin{eqnarray} \Lambda^{qq'}(\omega) &=& -\frac{1}{\hbar} \sum_{i,i'} [\mathbf{d}_{mn}^q]_i^{*} [\mathbf{d}_{m'n'}^{q'}]_{i'} \mathcal{P}\int_{-\infty}^{+\infty} \frac{\mathrm{d}\omega'}{\pi} \frac{\mathrm{Im} \,G_{ii'} (\mathbf{R}_q,\mathbf{R}_{q'},\omega')}{\omega'-\omega} \nonumber\\ & =& -\frac{1}{\hbar} \sum_{i,i'} [\mathbf{d}_{mn}^q]_i^{*} [\mathbf{d}_{m'n'}^{q'}]_{i'} \mathrm{Re} \, G_{ii'} (\mathbf{R}_q,\mathbf{R}_{q'},\omega). \end{eqnarray} \tag{ 27 }$$

We now specialize the derivation of the previous section to the case when the body is a slab of finite thickness δ, defined by the two interfaces z = 0 and −δ (see figure 3). In this simple case, explicit expressions for the transmission and reflection operators can be exploited [46, 47]. Owing to the translational invariance of a planar slab with respect to the xy plane, the slab reflection and transmission operators, $\mathcal {R}$ and $\mathcal {T}$, are diagonal and equal to

$$\begin{equation} \begin{array}{ll} \displaystyle \langle p,\mathbf{k}|\mathcal{R}|p',\mathbf{k}'\rangle=(2\pi)^2\delta(\mathbf{k}-\mathbf{k}')\delta_{pp'}\rho_{p}(\mathbf{k},\omega),\\ \displaystyle \langle p,\mathbf{k}|\mathcal{T}|p',\mathbf{k}'\rangle=(2\pi)^2\delta(\mathbf{k}-\mathbf{k}')\delta_{pp'}\tau_{p}(\mathbf{k},\omega), \end{array} \end{equation} \tag{ 28 }$$

where the Fresnel reflection and transmission coefficients modified by the finite thickness δ are given by (we recall that p = 1,2 corresponding to the TE and TM polarizations)

$$\begin{equation} \begin{array}{ll} \displaystyle \rho_{p}(\mathbf{k},\omega)=r_{p}(\mathbf{k},\omega)\frac{1-\mathrm{e}^{2\mathrm{i}k_{zm}\delta}}{1-r_{p}^2(\mathbf{k},\omega)\mathrm{e}^{2\mathrm{i}k_{zm}\delta}},\\ \displaystyle \tau_{p}(\mathbf{k},\omega)=\frac{t_{p}(\mathbf{k},\omega)\bar{t}_{p}(\mathbf{k},\omega){\rm e}^{{\rm i}(k_{zm}-k_z)\delta}}{1-r_{p}^2(\mathbf{k},\omega)\mathrm{e}^{2\mathrm{i}k_{zm}\delta}}. \end{array} \end{equation} \tag{ 29 }$$

In the previous equations we have introduced the z component of the K vector inside the medium

$$\begin{equation}k_{zm}=\sqrt{\varepsilon(\omega)\frac{\omega^2}{c^2}-\mathbf{k}^2}, \end{equation} \tag{ 30 }$$

ε(ω) being the dielectric permittivity of the slab, the ordinary vacuum–medium Fresnel reflection coefficients

$$\begin{equation}r_{\mathrm{TE}}=\frac{k_z-k_{zm}}{k_z+k_{zm}},\quad r_{\mathrm{TM}}=\frac{\varepsilon(\omega)k_z-k_{zm}}{\varepsilon(\omega)k_z+k_{zm}}, \end{equation} \tag{ 31 }$$

as well as both the vacuum–medium (denoted with t) and medium–vacuum (denoted with $\bar {t}$) transmission coefficients

$$\begin{equation} \begin{array}{ll} \displaystyle t_{\mathrm{TE}}=\frac{2k_z}{k_z+k_{zm}},\quad t_{\mathrm{TM}}=\frac{2\sqrt{\varepsilon(\omega)}k_{z}}{\varepsilon(\omega)k_z+k_{zm}},\\ \displaystyle \bar{t}_{\mathrm{TE}}=\frac{2k_{zm}}{k_z+k_{zm}},\quad\bar{t}_{\mathrm{TM}}=\frac{2\sqrt{\varepsilon(\omega)}k_{zm}}{\varepsilon(\omega)k_z+k_{zm}}. \end{array} \end{equation} \tag{ 32 }$$

Zoom In Zoom Out Reset image size

After replacing the matrix elements (28) in (23) we obtain for the α functions (we choose the x-axis along the vector r₁ − r₂ whose coordinates in the plane xy are then (r₁₂,0), being r₁₂ = |r₁ − r₂|),

$$\begin{eqnarray} &&\fl [\alpha_{\mathrm{W}}^{qq'}(\omega)]_{ii'}= \frac{3c}{8\pi\omega}\sum_{p}\Bigg\{\int_0^{\frac{\omega}{c}}\frac{\mathrm{d}k\,k}{k_z}[ {\rm e}^{{\rm i}k_z (z_q-z_{q'})} [N_p^{qq'}(k,\omega)]^{++}_{ii'} (|\rho_{p}(\mathbf{k},\omega)|^2+|\tau_{p}(\mathbf{k},\omega)|^2)\nonumber\\ &&+{\rm e}^{{\rm i}k_z(z_q+z_{q'})} [N_p^{qq'}(k,\omega)]^{+-}_{ii'}\rho_{p}(\mathbf{k},\omega)+{\rm e}^{-{\rm i}k_z(z_q+z_{q'})}[N_p^{qq'}(k,\omega)]^{-+}_{ii'}\rho_{p}(\mathbf{k},\omega)^{*}\nonumber\\ &&+{\rm e}^{-{\rm i}k_z (z_q-z_{q'})} [N_p^{qq'}(k,\omega)]^{--}_{ii'}]\Bigg\}, \nonumber\\ &&\fl [\alpha_{\mathrm{M}}^{qq'}(\omega)]_{ii'}= \frac{3c}{8\pi\omega}\sum_{p}\Bigg\{\int_0^{\frac{\omega}{c}}\frac{\mathrm{d}k\,k}{k_z} {\rm e}^{{\rm i}k_z (z_q-z_{q'})} [N_p^{qq'}(k,\omega)]^{++}_{ii'}(1-|\rho_{p}(\mathbf{k},\omega)|^2 -|\tau_{p}(\mathbf{k},\omega)|^2) \nonumber\\ &&-\mathrm{i}\int_{\frac{\omega}{c}}^{\infty}\frac{\mathrm{d}k\,k}{\mathrm{Im}(k_z)} \mathrm{e}^{-\mathrm{Im}(k_z)(z_q+z_{q'})} [N_p^{qq'}(k,\omega)]^{++}_{ii'} \bigl[\rho_{p}(\mathbf{k},\omega)-\rho_{p}(\mathbf{k},\omega)^{*}\bigr]\Bigg\} \end{eqnarray} \tag{ 33 }$$

where, using the fact that ρ_p(k,ω) and τ_p(k,ω) are independent of θ (the angle formed by k and the x-axis in the plane xy), we have introduced the angular integrals

$$\begin{equation} [N_p^{qq'}(k,\omega)]^{\phi\phi'}_{ii'}=\int_0^{2\pi} \frac{\mathrm{d}\theta}{\pi} [\hat{\boldsymbol{\epsilon}}_p^\phi(\mathbf{k},\omega)]_i [\hat{\boldsymbol{\epsilon}}_{p}^{\phi'}(\mathbf{k},\omega)]_{i'}^{*}\,{\rm e}^{{\rm i} k r_{qq'} \cos\theta}, \end{equation} \tag{ 34 }$$

where r₂₁ = −r₁₂. The matrix elements different from zero are, for p = 1, $[N_{1}^{qq'}]^{\phi \phi '}_{11}=\frac {2}{k r_{qq'}}J_1(k r_{qq'})$, $[N_{1}^{qq'}]^{\phi \phi '}_{22}=\frac {2}{k r_{qq'}}J_1(k r_{qq'})-2 J_2(k r_{qq'})$, while for p = 2 are

$$\begin{eqnarray} &&\fl [N_{2}^{qq'}]^{\phi\phi'}_{11} = \frac{2\phi \phi'c^2 |k_z^2|}{k r_{qq'}\omega^2}[J_1(k r_{qq'})- k r_{qq'} J_2(k r_{qq'})], \quad [N_{2}^{qq'}]^{\phi\phi'}_{13} = -\mathrm{i} \phi \frac{2c^2k k_z}{\omega^2}J_1(k r_{qq'} ),\nonumber\\ &&\fl [N_{2}^{qq'}]^{\phi\phi'}_{22} = \frac{2 \phi \phi' c^2 |k_z^2| }{k r_{qq'} \omega^2}J_1(k r_{qq'}), \quad [N_{2}^{qq'}]^{\phi\phi'}_{31} = -\mathrm{i} \phi' \frac{2c^2k k_z^{*}}{\omega^2}J_1(k r_{qq'} ), \\ &&\fl [N_{2}^{qq'}]^{\phi\phi'}_{33} = \frac{2c^2 k^2}{\omega^2}J_0(k r_{qq'} ),\nonumber \end{eqnarray} \tag{ 35 }$$

where J_n(x) is the nth-order Bessel function of the first kind. For x → 0, it is J₀(x) → 1, J₁(x) → 0, J₂(x) → 0 and J₁(x)/x → 1/2, so that [N^qq'₁₍₂₎(k,ω)]^ϕϕ'_ii' become diagonal and reduce to the vectors defined in (55) of [49] in the case of a single emitter.

To simplify the functions [α^qq'_W(ω)]_ii' and [α^qq'_M(ω)]_ii' in (33) we exploit the fact that the quantities [N^qq'₁(k,ω)]^ϕϕ'_ii' do not depend on ϕ and ϕ' and are real, and that in the propagative sector [N^qq'₂(k,ω)]⁺⁺_ii' = [N^qq'₂(k,ω)]⁻⁻ *_ii' and [N^qq'₂(k,ω)]⁺⁻_ii' = [N^qq'₂(k,ω)]⁻⁺ *_ii'. By using the angular integrals (34), equation (33) can thus be rewritten as

$$\begin{equation} \begin{array}{ll} \displaystyle [\alpha^{qq'}_{\mathrm{W}}(\omega)]_{ii'}=\frac{[A^{qq'}(\omega)]_{ii'}^{*}+[B^{qq'}(\omega)]_{ii'}+2[C^{qq'}(\omega)]_{ii'}}{2},\\ \displaystyle [\alpha^{qq'}_{\mathrm{M}}(\omega)]_{ii'}=\frac{[A^{qq'}(\omega)]_{ii'}-[B^{qq'}(\omega)]_{ii'}+2[D^{qq'}(\omega)]_{ii'}}{2}, \end{array} \end{equation} \tag{ 36 }$$

where we have introduced the integral matrices

$$\begin{eqnarray} &&\fl [A^{qq'}(\omega)]_{ii'} = \frac{3c}{4\omega} \sum_p\int_0^{\frac{\omega}{c}}\frac{k\,\mathrm{d}k}{k_z} {\rm e}^{{\rm i}k_z (z_q-z_{q'})} [N_p^{qq'}(k,\omega)]^{++}_{ii'},\nonumber\\ &&\fl [B^{qq'}(\omega)]_{ii'} = \frac{3c}{4\omega}\sum_p\int_0^{\frac{\omega}{c}}\frac{k\,\mathrm{d}k}{k_z}{\rm e}^{{\rm i}k_z (z_q-z_{q'})} [N_p^{qq'}(k,\omega)]^{++}_{ii'}(|\rho_{p}(k,\omega)|^2+|\tau_{p}(k,\omega)|^2),\\ &&\fl [C^{qq'}(\omega)]_{ii'} = \frac{3c}{4\omega}\sum_p\int_0^{\frac{\omega}{c}}\frac{k\,\mathrm{d}k}{k_z} \mathrm{Re}[ {\rm e}^{{\rm i}k_z(z_q+z_{q'})} [N_p^{qq'}(k,\omega)]^{+-}_{ii'} \rho_{p}(k,\omega)],\nonumber\\ &&\fl [D^{qq'}(\omega)]_{ii'} = \frac{3c}{4\omega}\sum_p\int_{\frac{\omega}{c}}^{+\infty}\!\!\!\!\frac{k\,\mathrm{d}k}{\mathrm{Im}(k_z)}\mathrm{e}^{-\mathrm{Im}(k_z)(z_q+z_{q'})} [N_p^{qq'}(k,\omega)]^{++}_{ii'}\mathrm{Im}[\rho_{p}(k,\omega)].\nonumber \end{eqnarray} \tag{ 37 }$$

For q = q', α^qq_W(ω) and α^qq_M(ω) coincide with the functions defined in (56) of [49] in the case of a single emitter. For q ≠ q', in the limit R_q' → R_q, [α^qq'_W(ω)]_ii' and [α^qq'_M(ω)]_ii' tend to their values in the case of a single emitter placed in R_q. In the limit, the distance between the two emitters goes to infinity, both α^qq'_W(ω) and α^qq'_M(ω) go to zero. In the limit of |r_qq'| → ∞, this is due to the fact that the functions [N^qq'_p(k,ω)]^ϕϕ' go to zero (for |x| → ∞, it is J₀(x) → 0, J₁(x) → 0 and J₂(x) → 0). In the limit |z_q − z_q'| → ∞, this is due to the presence of an oscillating function whose frequency goes to infinity in the integrals A(ω), B(ω) and C(ω) (this can be seen explicitly by integrating by parts) and to the presence of the exponential function going to zero in the integral D(ω).

Concerning the function Λ^qq'(ω), in order to develop its expression in (27) one has to compute the real part of the Green function in terms of the scattering operators $\mathcal {R}$ and $\mathcal {T}$. This is done in appendix C where a free term Λ^qq'₀(ω) (see (C.11)) remaining in the absence of matter has been isolated from a reflected part Λ^qq'_R(ω) (see (C.12)), Λ^qq'(ω) = Λ^qq'₀(ω) + Λ^qq'_R(ω). By using the expressions for $\mathcal {R}$ and $\mathcal {T}$ in the case of a slab and the angular integrals in (34), one can derive, starting from  (C.12):

$$\begin{eqnarray} &&\fl \mathrm{Re} \, G^{(\mathrm{R})}_{ii'}(\mathbf{R}_q,\mathbf{R}_{q'},\omega)=\frac{\mathrm{i} \omega^2}{4\epsilon_0c^2}\frac{1}{4 \pi}\sum_{p}\Biggl\{\int_0^{\frac{\omega}{c}}\frac{k\, \mathrm{d}k}{k_z} [{\rm e}^{{\rm i}k_z(z_q+z_{q'})} [N_p^{qq'}(k,\omega)]^{+-}_{ii'}\rho_{p}(\mathbf{k},\omega)\nonumber\\ &&-{\rm e}^{-{\rm i}k_z(z_q+z_{q'})}[N_p^{qq'}(k,\omega)]^{-+}_{ii'}\rho_{p}(\mathbf{k},\omega)^{*}-\mathrm{i}\int_{\frac{\omega}{c}}^{\infty}\frac{\mathrm{d}k\,k}{\mathrm{Im}(k_z)}\mathrm{e}^{-\mathrm{Im}(k_z)(z_q+z_{q'})}\nonumber\\ &&\times [N_p^{qq'}(k,\omega)]^{++}_{ii'} (\rho_{p}(\mathbf{k},\omega)+\rho_{p}(\mathbf{k},\omega)^{*})]\Biggr\}. \end{eqnarray} \tag{ 38 }$$

Equation (27) can be thus cast in the form

$$\begin{equation} \fl \Lambda^{qq'}(\omega)= \Lambda_0^{qq'}(\omega)+\sqrt{\Gamma_0^q(\omega)\Gamma_0^{q'}(\omega) } \sum_{i,i'}[\tilde{{\bf d}}_{mn}^q]^{*}_{i} [\tilde{{\bf d}}_{m'n'}^{q'}]_{i'} \bigr([C_2^{qq'}(\omega)]_{ii'}-[D_2^{qq'}(\omega)]_{ii'}), \end{equation} \tag{ 39 }$$

where we have introduced the integral matrices

$$\begin{equation} \begin{array}{lll} \hspace{-3pc}\displaystyle [C_2^{qq'}(\omega)]_{ii'}=\frac{3c}{8\omega}\sum_p\int_0^{\frac{\omega}{c}}\frac{k\,\mathrm{d}k}{k_z} \mathrm{Im}[ {\rm e}^{{\rm i}k_z(z_q+z_{q'})} [N_p^{qq'}(k,\omega)]^{+-}_{ii'} \rho_{p}(k,\omega)],\\ \hspace{-3pc}\displaystyle [D_2^{qq'}(\omega)]_{ii'}= \frac{3c}{8\omega}\sum_p\int_{\frac{\omega}{c}}^{+\infty}\!\!\!\!\frac{k\,\mathrm{d}k}{\mathrm{Im}(k_z)}\mathrm{e}^{-\mathrm{Im}(k_z)(z_q+z_{q'})} [N_p^{qq'}(k,\omega)]^{++}_{ii'}\mathrm{Re}[\rho_{p}(k,\omega)]. \end{array} \end{equation} \tag{ 40 }$$

We observe that the limit case when the body is absent is discussed in appendix B, where known expressions for Γ^qq'(ω) and Λ^qq'(ω) are retrieved.

From now on we specialize our investigation to the case of two emitters (qubits) characterized by two internal levels |1〉 ≡ |g〉 and |2〉 ≡ |e〉 with the same transition frequency ω = ω¹_e − ω¹_g = ω²_e − ω²_g. In this case, master equation (16) reduces to

$$\begin{eqnarray} &&\fl \frac{\mathrm{d}}{\mathrm{d} t}\rho = - \frac{\mathrm{i}}{\hbar} [H_{\rm S}+ \delta_{\rm S},\rho ]-\mathrm{i}\sum_{q\neq q'}\Lambda^{qq'}(\omega)[\sigma_{ge}^{q\,\dagger}\sigma_{ge}^{q'},\rho ]+\sum_{q, q'}\Gamma^{qq'}(\omega)\left(\vphantom{\frac{1}{2}}\sigma^{q'}_{ge}\rho \sigma_{ge}^{q\,\dagger}\right.\nonumber\\ &&\left. -\frac{1}{2}\{\sigma^{q\,\dagger}_{ge}\sigma^{q'}_{ge},\rho \}\right) +\sum_{q, q'}\Gamma^{qq'}(-\omega)\left(\sigma^{q'\dagger}_{ge}\rho \sigma^{q}_{ge}-\frac{1}{2}\{\sigma^{q}_{ge} \sigma^{q'\,\dagger}_{ge},\rho \}\right), \end{eqnarray} \tag{ 41 }$$

where we used [σ^q_gg,ρ] = −[σ^q_ee,ρ], so that

$$\begin{equation} \delta_{\rm S}=\sum_{q}\hbar [S^{qq}(\omega)-S^{qq}(-\omega)] \sigma_{ee}^q, \end{equation} \tag{ 42 }$$

and where the functions S^qq(± ω), Λ^qq'(ω) and Γ^qq'(± ω) are defined in (17) for the specific case {m,n} = {1,2} (in the following we use the notation d^q₁₂ = d^q).

We observe that master equation (41) can also describe the case in which the emitters' frequencies are close enough (but not identical) so that rotating wave approximation used in the derivation of (13) still holds. This typically occurs when the frequency difference is much smaller than the average frequency [26]. The operator δ_S (42) represents a shift of energy levels, the renormalized transition frequencies being equal to $\tilde {\omega }_e^1-\omega _g^1=\omega +S^{11}(\omega ) -S^{11}(-\omega )$ and $\tilde {\omega }_e^2-\omega _g^2= \omega +S^{22}(\omega ) -S^{22}(-\omega )$. In the following, these shifts will not play any role.

To discuss the properties of (41), we will use two different bases, the decoupled basis {|1〉 ≡ |gg〉,|2〉 ≡ |eg〉,|3〉 ≡ |ge〉,|4〉 ≡ |ee〉} and the coupled basis $\{|G\rangle \equiv |1\rangle ,|A\rangle \equiv (|2\rangle -|3\rangle )/\sqrt {2},|S\rangle \equiv (|2\rangle +|3\rangle )/\sqrt {2},|E\rangle \equiv |4\rangle \}$, where we have introduced the collective antisymmetric |A〉 and symmetric states |S〉. The coupled basis is the one diagonalizing the effective Hamiltonian, $H_{\mathrm { S}}+\delta _{\mathrm { S}}+\sum _{q\neq q'}\hbar \Lambda ^{qq'}(\omega )\sigma _{ge}^{q\,\dagger }\sigma _{ge}^{q'}$, appearing in the first line of (41). In particular, the sign of Λ¹²(ω) inverts the role of |A〉 and |S〉 in the eigenstates of the above effective Hamiltonian. The eigenvalues associated with |G〉, |A〉, |S〉, |E〉, are { − ℏω, −ℏ|Λ¹²(ω)|,ℏ|Λ¹²(ω)|,ℏω}, with respect to the energy E₀ = ℏ[ω¹_e + S¹¹(ω) − S¹¹(− ω) + ω²_e + S²²(ω) − S²²(− ω) − ω].

5.1. X states

In the decoupled basis we can distinguish elements along the two main diagonals of the two-qubit density matrix from the remaining ones because they are not connected through master equation (41). We thus focus our attention on the class of X states, having non-zero elements only along the main diagonal and anti-diagonal of the density matrix (we use the notation ρ_ij = 〈i|ρ|j〉),

$$\begin{equation} \rho_{\mathrm{X}} = \left(\begin{array}{@{}cccc@{}} \rho_{11} & 0 & 0 & \rho_{14} \\ 0 & \rho_{22} & \rho_{23} & 0 \\ 0 & \rho_{23}^{*} & \rho_{33} & 0 \\ \rho_{14} ^{*} & 0 & 0 & \rho_{44} \end{array}\right). \end{equation} \tag{ 43 }$$

Bell, Werner and Bell diagonal states belong to this class of states [55]. X-structure density matrices are found in a wide variety of physical situations and are also experimentally achievable [56]. For example, X states are encountered as eigenstates in all of the systems with odd–even symmetry such as in the Ising and the XY models [57]. Moreover, in many physical evolutions of open quantum systems an initial X structure is maintained in time [58], as it is in our case. Terms outside the two main diagonals initially populated would be eventually washed off asymptotically. In the following, the two-qubit state will always have an X structure.

5.2. Concurrence

We shall quantify the entanglement in the two-qubit dynamics by evaluating the concurrence, C(t) (C = 0 for separable states, C = 1 for maximally entangled states) [59]. For the X states it takes the form [58]

$$\begin{equation} \begin{array}{ll} C(t)=2\,\mathrm{max}\{0,K_1(t),K_2(t)\}, \\ K_1(t)=|\rho_{23}(t)|-\sqrt{\rho_{11}(t)\rho_{44}(t)}, \quad K_2(t)=|\rho_{14}(t)|-\sqrt{\rho_{22}(t)\rho_{33}(t)}. \end{array} \end{equation} \tag{ 44 }$$

The master equation (41) always induces an exponential decay for ρ₁₄(t), so that in the steady state only K₁(t) could be responsible for having C(∞) > 0.

To discuss new phenomena emerging out of thermal equilibrium, it will be instructive to rewrite K₁(t) in terms of the populations in the coupled basis (we use the notation ρ_IJ = 〈I|ρ|J〉 and ρ_I = 〈I|ρ|I〉):

$$\begin{equation} K_1(t)=\frac{1}{2}\sqrt{[\rho_{\mathrm{S}}(t)-\rho_{\mathrm{A}}(t)]^2+[\rho_{\mathrm{SA}}(t)-\rho_{\mathrm{AS}}(t)]^2} -\sqrt{\rho_{\mathrm{G}}(t)\rho_{\mathrm{E}}(t)}. \end{equation} \tag{ 45 }$$

We will see that out of thermal equilibrium, it is always ρ_AS(∞) = 0, but ρ_S(∞) and ρ_A(∞) can differ, so that K₁(∞) could be positive.

5.3. Thermal equilibrium

When T_W = T_M ≡ T, master equation (41) describes the thermalization toward the thermal equilibrium state, which is diagonal with the four steady populations given by

$$\begin{equation} \left(\begin{array}{@{}c@{}} \rho_{11}(\infty)\\ \rho_{22}(\infty)\\ \rho_{33}(\infty)\\ \rho_{44}(\infty)\\ \end{array}\right)_{\!\!\!\mathrm{eq}} =\frac{1}{Z_{\mathrm{eq}}}\left( \begin{array}{@{}c@{}} [1+ n(\omega, T)]^2 \\ n(\omega, T)[1+ n(\omega, T)] \\ n(\omega, T)[1+ n(\omega, T)]\\ n(\omega, T)^2 \end{array}\right), \end{equation} \tag{ 46 }$$

where Z_eq = [1 + 2n(ω,T)]². By moving to the coupled basis, the thermal state remains diagonal with ρ_S(∞) = ρ_A(∞) = ρ₂₂(∞) = ρ₃₃(∞).

As a mathematical remark, we note that the thermal state is always reached asymptotically except if the identities Γ¹¹(± ω) = Γ²²(± ω) = Γ¹²(± ω) = Γ²¹(± ω) ≡ Γ(± ω) are strictly verified. In this peculiar case, both in and out of thermal equilibrium, the steady state depends upon the initial state and may be entangled. In particular, it is diagonal in the coupled basis with populations equal to

$$\begin{equation} \left(\begin{array}{@c@} \rho_{\mathrm{G}}(\infty)\\ \rho_{\mathrm{A}}(\infty)\\ \rho_{\mathrm{S}}(\infty)\\ \rho_{\mathrm{E}}(\infty)\\ \end{array} \right) =\frac{1}{Z}\left(\begin{array}{@c@} \Gamma(\omega)^2 [1- \rho_{\mathrm{A}} (0)] \\ \rho_{\mathrm{A}} (0)\\ \Gamma(-\omega)\Gamma(\omega) [1- \rho_{\mathrm{A}} (0)]\\ \Gamma(-\omega)^2 [1- \rho_{\mathrm{A}} (0)] \end{array}\right), \end{equation} \tag{ 47 }$$

where Z = Γ(− ω)² + Γ(ω)Γ(− ω) + Γ(ω)². Apart from this case, at thermal equilibrium the steady state is always a thermal state, thus not entangled. We can see it by looking at the concurrence (44) which is zero being ρ₂₃(∞) = 0. This can also be seen in the coupled basis, where ρ_AS(∞) = 0 and ρ_S(∞) = ρ_A(∞), so that K₁(∞) (45) is negative.

5.4. Out of thermal equilibrium: an instructive case

When T_W ≠ T_M, qualitative differences emerge in the dynamics and in the steady states. To highlight these new features, we first consider a simple case where a clear physical interpretation in terms of |S〉 and |A〉 is available. This is the case when Γ¹¹(± ω) = Γ²²(± ω) ≡ Γ(± ω) and Γ¹²⁽²¹⁾(± ω) are real. These conditions are verified, for example, in the case of identical qubits, with d¹ = d² ≡ d, placed in equivalent positions with respect to the body (in the case of a slab, z₁ = z₂) and with d real and having components different from zero either only along the z-axis or only along the plane xy. In this case, master equation (41) gives in the coupled basis a set of rate equations for the populations, which are decoupled from the other density matrix elements:

$$\begin{equation} \begin{array}{lll} \displaystyle \dot{\rho}_{\mathrm{G}}=-(\Gamma_{\mathrm{A}} \,n_{\mathrm{A}}+\Gamma_{\mathrm{S}}\, n_{\mathrm{S}})\rho_{\mathrm{G}}+\Gamma_{\mathrm{A}}(1+n_{\mathrm{A}}) \rho_{\mathrm{A}}+\Gamma_{\mathrm{S}}(1+n_{\mathrm{S}}) \rho_{\mathrm{S}},\\ \displaystyle \dot{\rho}_{\mathrm{A}}=-\Gamma_{\mathrm{A}} (1+2 n_{\mathrm{A}}) \rho_{\mathrm{A}}+\Gamma_{\mathrm{A}}\, n_{\mathrm{A}} \rho_{\mathrm{G}}+\Gamma_{\mathrm{A}}(1+n_{\mathrm{A}}) \rho_{\mathrm{E}},\\ \displaystyle \dot{\rho}_{\mathrm{S}}=-\Gamma_{\mathrm{S}} (1+2 n_{\mathrm{S}}) \rho_{\mathrm{S}}+\Gamma_{\mathrm{S}}\, n_{\mathrm{S}} \rho_{\mathrm{G}}+\Gamma_{\mathrm{S}}(1+n_{\mathrm{S}}) \rho_{\mathrm{E}},\\ \displaystyle \dot{\rho}_{\mathrm{E}}=-[\Gamma_{\mathrm{A}} (1+2 n_{\mathrm{A}}) +\Gamma_{\mathrm{S}} (1+2 n_{\mathrm{S}}) ] \rho_{\mathrm{E}}+\Gamma_{\mathrm{A}}\, n_{\mathrm{A}} \rho_{\mathrm{A}} +\Gamma_{\mathrm{S}}\, n_{\mathrm{S}} \rho_{\mathrm{S}}. \end{array} \end{equation} \tag{ 48 }$$

Here, the coefficient Γ₀(ω) has been absorbed by the time variable in the derivative, which is now dimensionless, and we have used the relations

$$\begin{equation} \begin{array}{ll} \hspace{-3pc}\displaystyle \Gamma(\omega)-\Gamma^{12}(\omega) =\Gamma_0(\omega) \Gamma_{\mathrm{A}} (1+n_{\mathrm{A}}), \quad \Gamma(\omega)+\Gamma^{12}(\omega)=\Gamma_0(\omega) \Gamma_{\mathrm{S}} (1+n_{\mathrm{S}}),\\ \hspace{-3pc}\displaystyle \Gamma(-\omega)-\Gamma^{12}(-\omega) =\Gamma_0(\omega) \Gamma_{\mathrm{A}} \,n_{\mathrm{A}},\quad \Gamma(-\omega)+\Gamma^{12}(-\omega)=\Gamma_0(\omega) \Gamma_{\mathrm{S}} \,n_{\mathrm{S}} \end{array} \end{equation} \tag{ 49 }$$

with

$$\begin{equation} \begin{array}{ll} \displaystyle \Gamma_{\mathrm{A}}=\alpha_{\mathrm{W}}(\omega) -\alpha^{12}_{\mathrm{W}}(\omega)+\alpha_{\mathrm{M}}(\omega) -\alpha^{12}_{\mathrm{M}}(\omega),\\ \displaystyle \Gamma_{\mathrm{S}}=\alpha_{\mathrm{W}}(\omega) +\alpha^{12}_{\mathrm{W}}(\omega)+\alpha_{\mathrm{M}}(\omega) +\alpha^{12}_{\mathrm{M}}(\omega),\\ \displaystyle n_{\mathrm{A}} =\frac{1}{\Gamma_{\mathrm{A}}}\{ [\alpha_{\mathrm{W}}(\omega) -\alpha^{12}_{\mathrm{W}}(\omega)] n(\omega,T_{\mathrm{W}}) +[\alpha_{\mathrm{M}}(\omega) -\alpha^{12}_{\mathrm{M}}(\omega)]n(\omega,T_{\mathrm{M}}) \},\\ \displaystyle n_{\mathrm{S}} = \frac{1}{\Gamma_{\mathrm{S}}}\{ [\alpha_{\mathrm{W}}(\omega) +\alpha^{12}_{\mathrm{W}}(\omega)] n(\omega,T_{\mathrm{W}})+[\alpha_{\mathrm{M}}(\omega) +\alpha^{12}_{\mathrm{M}}(\omega)]n(\omega,T_{\mathrm{M}}) \}, \end{array} \end{equation} \tag{ 50 }$$

where α_W(M)(ω) ≡ α¹¹_W(M)(ω) = α²²_W(M)(ω). We remark that Λ^qq'(ω) does not enter in the rate equations (48), which are schematically represented in figure 4. We observe that with each decay channel from |E〉 to |G〉 we can associate distinct effective temperatures T_S and T_A confined between T_W and T_M in correspondence to the effective number of photons n_S and n_A, which have the property of being confined between n(ω,T_W) and n(ω,T_M) [49].

Zoom In Zoom Out Reset image size

Concerning the coherences in the second diagonal:

$$\begin{equation} \begin{array}{ll} \displaystyle \dot{\rho}_{\mathrm{A}\mathrm{S}}=-\frac{1}{2}\left[\Gamma_{\mathrm{A}} (1+2 n_{\mathrm{A}}) +\Gamma_{\mathrm{S}}(1+2 n_{\mathrm{S}}) -4 \mathrm{i}\Lambda^{12}(\omega)\right] \rho_{\mathrm{A}\mathrm{S}},\\ \displaystyle \dot{\rho}_{\mathrm{G}\mathrm{E}}=-\frac{1}{2}\left[\Gamma_{\mathrm{A}} (1+2 n_{\mathrm{A}}) +\Gamma_{\mathrm{S}}(1+2 n_{\mathrm{S}})-4\mathrm{ i} \Lambda^{12}(\omega) \right] \rho_{\mathrm{G}\mathrm{E}}, \end{array} \end{equation} \tag{ 51 }$$

which give for each coherence an exponential decay modulating oscillations due to Λ¹²(ω). The stationary solution of (48) is

$$\begin{equation} \left(\begin{array}{@{}c@{}} \rho_{\mathrm{G}}(\infty)\\ \rho_{\mathrm{A}}(\infty)\\ \rho_{\mathrm{S}}(\infty)\\ \rho_{\mathrm{E}}(\infty)\\ \end{array}\right)_{\!\!\!\mathrm{neq}} =\frac{1}{Z_{\mathrm{neq}}} \left(\begin{array}{@{}c@{}} (1+n_{\mathrm{A}})^2 (1+2n_{\mathrm{S}})\Gamma_{\mathrm{A}}+(1+2 n_{\mathrm{A}})(1+n_{\mathrm{S}})^2\Gamma_{\mathrm{S}}\\ n_{\mathrm{A}}(1+n_{\mathrm{A}})(1+2n_{\mathrm{S}}) \Gamma_{\mathrm{A}}+[ n_{\mathrm{A}}(1+2n_{\mathrm{S}}) \\ +n_{\mathrm{S}}^2(1+2n_{\mathrm{A}}) ] \Gamma_{\mathrm{S}} \\ n_{\mathrm{S}}(1+n_{\mathrm{S}})(1+2n_{\mathrm{A}}) \Gamma_{\mathrm{S}}+[ n_{\mathrm{S}}(1+2n_{\mathrm{A}}) \\ +n_{\mathrm{A}}^2(1+2n_{\mathrm{S}}) ] \Gamma_{\mathrm{A}} \\ n_{\mathrm{A}}^2 (1+2n_{\mathrm{S}})\Gamma_{\mathrm{A}}+(1+2 n_{\mathrm{A}})n_{\mathrm{S}}^2\Gamma_{\mathrm{S}} \end{array}\right), \end{equation} \tag{ 52 }$$

where Z_neq is the sum of the elements of the vector in the second line of the above equation. Out of equilibrium ρ₂₃(∞)_neq is different from zero and is given by

$$\begin{equation} \rho_{23}(\infty)_{\mathrm{neq}} =\frac{(n_{\mathrm{S}}-n_{\mathrm{A}})(\Gamma_{\mathrm{S}}+\Gamma_{\mathrm{A}})}{2 Z_{\mathrm{neq}}}, \end{equation} \tag{ 53 }$$

where we see easily how it tends to zero at thermal equilibrium when n_S = n_A. By using (52) in (44) and (45), we obtain for the steady concurrence:

$$\begin{eqnarray} &&\fl C(\infty)_{\mathrm{neq}}= \frac{2}{Z_{\mathrm{neq}}} [ |n_{\mathrm{S}}-n_{\mathrm{A}}|(\Gamma_{\mathrm{S}}+\Gamma_{\mathrm{A}}) /2\,-\sqrt{(1+n_{\mathrm{A}})^2 (1+2n_{\mathrm{S}})\Gamma_{\mathrm{A}}+(1+2 n_{\mathrm{A}})(1+n_{\mathrm{S}})^2\Gamma_{\mathrm{S}}}\nonumber\\ &&\times\sqrt{n_{\mathrm{A}}^2 (1+2n_{\mathrm{S}})\Gamma_{\mathrm{A}}+(1+2 n_{\mathrm{A}})n_{\mathrm{S}}^2\Gamma_{\mathrm{S}}} \; ] . \end{eqnarray} \tag{ 54 }$$

Simplifying Γ_S, C(∞)_neq becomes the function of the three dimensionless quantities Γ_A/Γ_S, n_S and n_A. This dependence is discussed in figure 5, where C = C(∞)_neq is depicted as a function of n_A for n_S = 0.001 and for several values of Γ_A/Γ_S, as indicated in the legend. We observe that by decreasing the value of Γ_A/Γ_S, higher values of C are reachable at higher values of n_A. The maximum value of C is 1/3, which can be obtained in the limits Γ_A/Γ_S → 0, n_S → 0 and n_A → ∞. In particular, the corresponding maximally entangled state, which is a statistical mixture of the ground and of the antisymmetric state with weights, respectively, equal to 2/3 and 1/3, has also been found in [35]. For smaller values of n_S the behavior remains almost identical, while by increasing its value, the values of C decrease progressively. We remark that identical behavior is found in the opposite case, i.e. when Γ_S/Γ_A → 0, the case in which the roles of states |S〉 and |A〉 are inverted. This can be achieved by looking for values of the various parameters such that α¹²_M(W)(ω) is negative and very close to α_M(W)(ω) in order to make the ratio Γ_S/Γ_A very small.

Zoom In Zoom Out Reset image size

Figure 5 describes the generation of steady entangled states emerging only in the absence of thermal equilibrium. Two main conditions must be fulfilled, the first being to have small values for Γ_A/Γ_S and the second to realize (quite) different effective temperatures for the two decay channels, which can be achieved only in the absence of thermal equilibrium.

Here, we report the numerical investigation concerning the case treated in section 4 when the body close to the emitters is a slab of finite thickness δ. According to (13), a relevant parameter involved in our investigation concerning the role of the body is the value of the dielectric permittivity at the common transition frequency of the two qubits. As a material we choose silicon carbide (SiC) whose dielectric permittivity ε(ω) is described by using a Drude–Lorentz model [53]

$$\begin{equation} \varepsilon (\omega) =\varepsilon_{\infty}\frac{\omega ^2-\omega_{\mathrm{l}}^2+\mathrm{i}\Gamma\omega}{\omega^2-\omega_{\rm r}^2+\mathrm{i}\Gamma\omega}, \end{equation} \tag{ 55 }$$

characterized by a resonance at ω_r = 1.495 × 10¹⁴ rad s⁻¹ and where ε_∞ = 6.7, ω_l = 1.827 × 10¹⁴ rad s⁻¹ and Γ = 0.009 × 10¹⁴ rad s⁻¹. This model implies a surface phonon–polariton resonance at ω_p = 1.787 × 10¹⁴ rad s⁻¹. A relevant length scale in this case is c/ω_r ≃ 2 μm while a reference temperature is ℏω_r/k_B ≃ 1140 K. We will assume that ε(ω) does not vary much in the interval of temperatures considered. In the following study, we explore a region of parameters much wider than that allowing the analytical description of section 5.4.

6.1. Steady configurations

We first focus on the properties of steady states. In particular, we are interested in the amount of entanglement present asymptotically, which is quantified by the concurrence (44). This analysis is supported by an analytical solution of the steady state of (41), which is not reported here, since it is particularly cumbersome.

In figure 6(a), we plot the maximum of steady concurrence obtained for an interval of transition frequencies ranging from 0.3 ω_r to 1.7 ω_r, in the case of δ = 0.01 μm. In our numerical sample, z₁ and z₂ may vary between 0.05 and 50 μm and r₁₂ between 0 and 15 μm. The two temperatures range in an interval such that the associated number of photons is between 0 and 3 (we checked that larger values are not needed). The red curve is relative to the case of dipoles oriented along the z-axis, while the green curve to the case of dipoles oriented along the x-axis (in figure 9 we will show that this is the best choice if we limit ourselves to directions lying in the xy plane). Higher values of concurrence are obtained immediately before/after the resonance frequency ω_r. The red configuration always gives better results except around the surface phonon–polariton frequency ω_p ≃ 1.2ω_r where choosing the dipole directions along the x-axis is the best choice. The values of the parameters corresponding to each maximum vary with frequency. The best configuration is always characterized by values of n(ω,T_W) close to zero and n(ω,T_M) between 1 and 3. Smaller values of n(ω,T_M) are needed in the green curve. The zone where to place the qubits is around 1 μm from the slab at 0.3ω_r, gradually decreasing (specially after ω_r) down to 0.25 μm at 1.7ω_r. For the red curve the best choice is always z₁ close to z₂ (in our numerical sample, we limit the minimal distance at the order of 0.1 μm) and r₁₂ = 0, while for the green curve it is z₁ = z₂ and r₁₂ small (of the order of 0.01 μm). This means that the best configuration is when the interatomic axis is aligned with the dipoles direction. For ω around ω_p we point out the occurrence of larger values of C for different values of δ, points indicated with a cross above the green curve. In the absence of large values of C in correspondence to the canonical choice of the parameters described above, small values of C become evident for a different set of parameters. This corresponds to larger values of δ (of the order of 1 μm or more), z₁ ≃ 2 μm, z₂ ≃ 4 μm and r₁₂ ≃ 0.5 μm. In part (b), we plot the dependence of C on δ for several values of ω as indicated in the figure. The maximum of C is always obtained close to δ = 0.01 μm, which is the value chosen in part (a), except around ω_r where much smaller values of δ are required. This explains why in part (a) concurrence decreases around ω_r for δ = 0.01 μm.

Zoom In Zoom Out Reset image size

In figure 7, we plot the steady concurrence as a function of the position of the second qubit z₂ and of the slab temperature T_M for four different values of r₁₂. From (a) to (d) the two-qubit distance [r²₁₂ + (z₁ − z₂)²]^1/2 increases leading to a progressive decrease of the values of concurrence generated. A maximum of C ≃ 0.24 is obtained in part (a) for z₂ ≃ 1.3 μm and T_M ≃ 1100 K. The white lines correspond to the case z₂ = z₁ for which equation (54) holds for concurrence. In part (a), the maximum along the white curve is C ∼ 0.222 in correspondence to Γ_A/Γ_S ∼ 4.6 × 10⁻⁷, n_S ∼ 0.02 and n_A ∼ 1.56 which correspond to effective temperatures for the two decay channels T_S ∼ 90 K and T_A ∼ 690 K. We observe that very high temperatures are considered in this plot only to highlight the entire region where steady entanglement is present. At unphysical temperatures (e.g. above the melting temperature), the plot is only indicative of what would occur if a different material was chosen such that similar values of ε(ω) (for SiC it is ε(0.3 ω_r) ∼ 10.3 + 0.00721i) were encountered at lower frequencies. In this case, similar behavior for steady concurrence at lower temperatures is expected. However, we remark that in our case values of C higher than 0.14 are already present at T_M ≃ 500 K in part (a).

Zoom In Zoom Out Reset image size

In figure 8, we discuss the behavior of n_S, n_A and Γ_A/Γ_S, appearing in (54), as a function of z. We plot in part (a) n_S and n_A as a function of z = z₁ = z₂ and compare them with the value of n(ω,T) computed at the minimal (here T_W = 30 K) and maximal (here T_M = 1200 K) temperature considered. The temperatures and the other parameters are equal to the ones giving the maximum of concurrence in figure 7(b) along the white line. In the inset (part (b)) we plot Γ_A/Γ_S as a function of z = z₁ = z₂. The plot proves that, near z ≃ 1 μm, both conditions to reach high values of C are satisfied: small values of n_S and Γ_A/Γ_S in correspondence with high enough values of n_A (see also figure 5 for a comparison).

Zoom In Zoom Out Reset image size

In figure 9, we analyze the dependence of steady concurrence on dipole orientations. In general, higher results are obtained when the two dipoles are parallel. We use again the set of parameters corresponding to the maximum in figure 7(b) along the white line, which is obtained for dipoles along the z-axis. We show how concurrence decreases by changing the dipole directions toward the x-axis (part (a)) always lying on the xz plane and then toward the y-axis (part (b)) always lying on the xy plane (see insets in figures 9(a) and (b)). From part (b) it emerges that aligning the dipoles direction with the interatomic axis (which here is the x-axis) is the optimal choice in the xy plane inducing a lack of symmetry in this plane between the x and y directions.

Zoom In Zoom Out Reset image size

6.2. Dynamics

Here, we discuss the dynamic behavior of two-qubit density matrix elements and concurrence out of thermal equilibrium, also making comparisons with the thermal equilibrium case. This analysis is performed by solving numerically the evolution governed by (41).

In figure 10, we plot several density matrix elements and concurrence as a function of dimensionless time Γ₀(ω)t. Parts (a) and (b) concern the case of maximally entangled initial states, respectively, the antisymmetric state |A〉 in (a) and the symmetric state |S〉 in (b) (see the values of the parameters in the caption of the figure). Quite different dynamic behavior is pointed out. While starting from |A〉 entanglement is just preserved at a high value, starting from |S〉 concurrence first decreases (going to zero) mainly because of the decrease of ρ_S and then revives because of the increase of ρ_A. Dynamic creation of entanglement is yet more evident in part (c) where the initial state is the factorized state |2〉. In this case, concurrence is initially zero and increases because of the mediated interaction between qubits. Oscillations of C and ρ₂₃ are linked to the behavior of ρ_AS(SA) which rapidly oscillate (see (51)) because of the large value of Λ¹²(ω) which here is equal to Λ¹²(ω)/Γ₀(ω) ≃ − 2.3 × 10⁴. The oscillations are present because ρ_AS(SA) are initially populated (ρ_AS(SA) = 1/2) while they are not in the cases plotted in parts (a) and (b). States |S〉 and |A〉 are initially equally populated and become different asymptotically.

Zoom In Zoom Out Reset image size

In figure 11, we compare the evolution of concurrence out of thermal equilibrium with the evolutions at equilibrium at the minimal temperature T_min = T_W = T_M = 30 K and at the maximal temperature T_max = T_W = T_M = 1300 K. Two initial maximally entangled configurations are compared, the antisymmetric state |A〉 in part (a) and the symmetric state |S〉 in part (b). At thermal equilibrium, concurrence vanishes on shorter times by increasing the temperature, while out of equilibrium steady entanglement is present. At equilibrium, a larger decay time is observed by starting from the antisymmetric state (see also figure 12 on this subject). In part (b), out of equilibrium, concurrence decays on the same equilibrium timescale, the two-qubit state becoming separable, but it reemerges successively. Both in (a) and (b) a large amount of the initial entanglement is thus asymptotically preserved. In part (c) the initial state is the factorized state |2〉. The main difference here is that concurrence presents strong oscillations (see comment on part (c) of figure 10). At thermal equilibrium, entanglement eventually vanishes on a timescale similar to the one of part (a), while out of equilibrium it is maintained after its creation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 12, we discuss the dependence of super- and sub-radiant effects from the presence/absence of thermal equilibrium. Here, super- and sub-radiance are connected to the occurrence of a decay rate larger or smaller than the one observed in the case of independent qubits, phenomenon due to the interaction of the qubits with a common environment and which depend on the nature of their initial state [26]. In particular, we compare the evolution of the ground state population ρ_G starting from |S〉 and |A〉 for two different values of r₁₂ (0.25 and 15 μm) at thermal equilibrium at T_W = T_M = 100 K (part (a)) and at T_W = T_M = 800 K (part (b)) and out of thermal equilibrium for T_W = 100 K, T_M = 800 K (part (c)). The figure evidences super-radiant behavior when the initial state is |S〉 and sub-radiant when it is |A〉. The faster or slower increase of ρ_G is due to the role of Γ₁₂(± ω) which in the two channel decay rates of (49) is summed to Γ(± ω) in the |S〉 case and subtracted in the |A〉 case. By increasing the value of r₁₂, Γ₁₂(± ω) decreases and the decay rates Γ^p(m)_A and Γ^p(m)_S, defined in the caption of figure 4, tend to the same value Γ(+ (−)ω)/Γ₀(ω), which is the decay rate in the case of single emitters. At thermal equilibrium, the asymptotic state is independent of the values of r₁₂ while this is not the case out of thermal equilibrium, as pointed out in part (c).

We finally remark that relevant differences are expected when the Markovian and the rotating wave approximation, here adopted, are not valid. In the non-Markovian regime another source of oscillations in the dynamics of concurrence typically emerges [13], while the effect of counter-rotating terms is known to modify the creation of entanglement between the two emitters [60].

In this paper, we have investigated a system made of two quantum emitters interacting with a common stationary electromagnetic field out of thermal equilibrium generated by an arbitrary body and by the surrounding walls held at fixed different temperatures. The environmental field is characterized by means of its correlation functions out of equilibrium which also depend on the scattering properties of the body. We have derived the expressions in the absence of thermal equilibrium of the various functions governing the dissipative dynamics of the two emitters and compared them with the ones holding at thermal equilibrium. This has been done in the case of emitters characterized by an arbitrary number of levels. We have then specialized our investigation to the case of two qubits discussing the new features emerging out of thermal equilibrium.

For a restricted parameter region we have analytically shown that the absence of equilibrium may lead to the generation of steady entangled states. This phenomenon has been interpreted in terms of different effective temperatures associated with two decay channels connecting the total excited and ground states via the symmetric and antisymmetric states, respectively. The two-qubit dynamics can be directed toward mixed states where the antisymmetric contribution is larger than the symmetric one (or vice versa), resulting in the presence of steady entanglement. In this specific case a value of 1/3 has been found as maximum for the concurrence, quantifying the steady entanglement.

We have then numerically investigated the general dependence of steady states and dynamics on the various parameters, without any restriction on the decay rates, in the case the body placed in proximity of the two qubits is a slab made of SiC. The dependence of steady entanglement on the two-qubit distance, their common transition frequency with respect to the slab resonances, the slab thickness, the dipoles orientations and the two involved temperatures have been discussed. Values of concurrence up to 0.24 have been found. Protection and/or generation of entanglement according to the nature of the two-qubit initial state, entangled or not, has been pointed out, also comparing entanglement dynamics in the presence or absence of thermal equilibrium. Higher values of steady concurrence are found for transition frequencies far from the slab resonances (ω/ω_r = 0.3) and small thickness (δ ≃ 0.01 μm). Remarkably, steady entanglement can be obtained by starting from configurations at thermal equilibrium and by increasing one of the two temperatures involved in the environment of the two qubits.

The possibility to observe the effects we discussed could be explored, for example, for emitters made by trapped atoms [42] or by artificial atoms such as quantum dots or superconducting qubits, placed in proximity to a substrate held at a temperature different from that of the cell surrounding the emitters and the substrate.

We thank R Messina for useful discussions and acknowledge financial support from the Julian Schwinger Foundation. MA is a member of the LabEx NUMEV.

Here, we connect the correlation functions to T_M and T_W and to the properties of the body as well. To this purpose, we first develop the connection between (14) and the correlation functions in frequency space. By using (6) and homogeneity in time, we have

$$\begin{eqnarray} &&\fl \Xi_{ii'}^{qq'}(\omega)=\frac{1}{\hbar^2}\int_{0}^\infty \mathrm{d}s\int_0^{+\infty}\frac{\mathrm{d}\omega'}{2\pi}\int_0^{+\infty}\frac{\mathrm{d}\omega''}{2\pi}[{\rm e}^{-{\rm i}(\omega''-\omega)s}\langle E_i(\mathbf{R}_q,\omega'')E_{i'}^\dagger(\mathbf{R}_{q'},\omega')\rangle\nonumber\\ &&+{\rm e}^{{\rm i}(\omega+\omega'')s}\langle E_i^\dagger(\mathbf{R}_q,\omega'')E_{i'}(\mathbf{R}_{q'},\omega')\rangle], \end{eqnarray} \tag{ A.1 }$$

where we have used 〈E_i(R_q,ω'')E_i'(R_q',ω')〉 = 〈E^†_i(R_q,ω'')E^†_i'(R_q',ω')〉 = 0. By using $\int _{0}^\infty \mathrm {d}s\exp (-\mathrm {i}\epsilon s)=\pi \delta (\epsilon ) - \mathrm {i} \mathcal {P}\frac {1}{\epsilon }$ (where $\mathcal {P}$ indicates the principal part of the integral), we obtain from the previous equation and (14) (we assume ω > 0)

$$\begin{eqnarray} \fl &&\gamma_{ii'}^{qq'}(\omega)= \frac{1}{\hbar^2} \int_0^{+\infty}\frac{ \mathrm{d}\omega'}{2\pi} \langle E_i(\mathbf{R}_q,\omega)E_{i'}^\dagger(\mathbf{R}_{q'},\omega')\rangle,\nonumber \\ \fl &&\gamma_{ii'}^{qq'}(-\omega)= \frac{1}{\hbar^2} \int_0^{+\infty}\frac{ \mathrm{d}\omega'}{2\pi}\langle E_i^\dagger(\mathbf{R}_q,\omega)E_{i'}(\mathbf{R}_{q'},\omega')\rangle,\nonumber \\ \fl &&s_{ii'}^{qq'}(\omega)= \frac{1}{\hbar^2}\mathcal{P}\int_0^{+\infty}\frac{\mathrm{d}\omega'}{2\pi}\int_0^{+\infty}\frac{\mathrm{d}\omega''}{2\pi} \Bigg[\frac{\langle E_i(\mathbf{R}_q,\omega'')E_{i'}^\dagger(\mathbf{R}_{q'},\omega')\rangle}{\omega-\omega''}\nonumber \\ \fl & &\hspace{6pc}+\frac{\langle E_i^\dagger(\mathbf{R}_q,\omega'')E_{i'}(\mathbf{R}_{q'},\omega')\rangle}{\omega+\omega''}\Bigg]. \end{eqnarray} \tag{ A.2 }$$

By using the decomposition in (7), we obtain

$$\begin{eqnarray} &&\fl \langle E_i(\mathbf{R}_q,\omega)E_{i'}^\dagger(\mathbf{R}_{q'},\omega')\rangle=\sum_{\phi,\phi',p,p'}\int\frac{\mathrm{d}^2\mathbf{k}}{(2\pi)^2}\int\frac{\mathrm{d}^2\mathbf{k}'}{(2\pi)^2} {\rm e}^{{\rm i}(\mathbf{K}^\phi\cdot\mathbf{R}_q-\mathbf{K'}^{\phi'*}\cdot\mathbf{R}_{q'})}\nonumber \\ &&\times [\hat{\boldsymbol{\epsilon}}_p^\phi(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}_{p'}^{\phi'}(\mathbf{k}',\omega')]_{i'}^{*}\langle E_p^\phi(\mathbf{k},\omega)E_{p'}^{\phi'\dagger}(\mathbf{k}',\omega')\rangle \end{eqnarray} \tag{ A.3 }$$

and

$$\begin{eqnarray} &&\fl \langle E_i^\dagger(\mathbf{R}_q,\omega)E_{i'}(\mathbf{R}_{q'},\omega')\rangle=\sum_{\phi,\phi',p,p'}\int\frac{\mathrm{d}^2\mathbf{k}}{(2\pi)^2}\int\frac{\mathrm{d}^2\mathbf{k}'}{(2\pi)^2}\,{\rm e}^{-{\rm i} (\mathbf{K}^{\phi*}\cdot\mathbf{R}_q-\mathbf{K'}^{\phi'}\cdot\mathbf{R}_{q'})}\nonumber\\ &&\times [\hat{\boldsymbol{\epsilon}}_p^\phi(\mathbf{k},\omega)]_i^{*}[\hat{\boldsymbol{\epsilon}}_{p'}^{\phi'}(\mathbf{k}',\omega')]_{i'}\langle E_p^{\phi\dagger}(\mathbf{k},\omega)E_{p'}^{\phi'}(\mathbf{k}',\omega')\rangle, \end{eqnarray} \tag{ A.4 }$$

where ω > 0. We observe that the last equation can be obtained by taking the complex conjugate of (A.3) after having interchanged the operators E^ϕ_p(k,ω) and E^ϕ'†_p'(k',ω').

We now combine equations (18) and (19) to obtain the symmetrized correlation functions of the amplitude operator of the total field in the region of interest

$$\begin{eqnarray} &&\fl\langle E^{+}_p(\mathbf{k},\omega)E^{+\dagger}_{p'}(\mathbf{k}',\omega')\rangle_{\mathrm{sym}}=2\pi\delta(\omega-\omega')\frac{\omega}{2\epsilon_0c^2} \langle p,\mathbf{k}|[N(\omega,T_{\mathrm{M}})(\mathcal{P}_{-1}^{\mathrm{(pw)}}-\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger}\nonumber\\ &&+\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(ew)}} -\mathcal{P}_{-1}^{\mathrm{(ew)}}\mathcal{R}^{\dagger}-\mathcal{T}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{T}^{\dagger}) +N(\omega,T_{\mathrm{W}})(\mathcal{T}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{T}^{\dagger}+\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger})]|p',\mathbf{k}'\rangle,\nonumber\\ &&\fl \langle E^{+}_p(\mathbf{k},\omega)E^{-\dagger}_{p'}(\mathbf{k}',\omega')\rangle_{\mathrm{sym}}=2\pi\delta(\omega-\omega')\frac{\omega}{2\epsilon_0c^2} N(\omega,T_{\mathrm{W}})\langle p,\mathbf{k}|\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}}|p',\mathbf{k}'\rangle,\nonumber\\ &&\fl \langle E^{-}_p(\mathbf{k},\omega)E^{+\dagger}_{p'}(\mathbf{k}',\omega')\rangle_{\mathrm{sym}}=2\pi\delta(\omega-\omega')\frac{\omega}{2\epsilon_0c^2}N(\omega,T_{\mathrm{W}})\langle p,\mathbf{k}|\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger}|p',\mathbf{k}'\rangle,\nonumber\\ &&\fl \langle E^{-}_p(\mathbf{k},\omega)E^{-\dagger}_{p'}(\mathbf{k}',\omega')\rangle_{\mathrm{sym}}=2\pi\delta(\omega-\omega')\frac{\omega}{2\epsilon_0c^2}N(\omega,T_{\mathrm{W}}) \langle p,\mathbf{k}|\mathcal{P}_{-1}^{\mathrm{(pw)}}|p',\mathbf{k}'\rangle. \end{eqnarray} \tag{ A.5 }$$

However, in order to develop equations (A.3) and (A.4) we need the non-symmetrized versions of these correlation functions. To compute them, we first remark that the source correlation functions reported in (19) have been derived by using thermal-equilibrium techniques at the temperature of each source individually (see [47] for a detailed discussion). It follows that we can use Kubo's prescription [52], according to which in order to obtain 〈AB〉 from 〈AB〉_sym the replacement N(ω,T_i) → ℏω[1 + n(ω,T_i)] must be performed, while 〈BA〉 results from the replacement N(ω,T_i) → ℏω n(ω,T_i).

By using (A.5) in (A.3) we obtain for the antinormally ordered correlation functions, the form

$$\begin{equation} \langle E_i (\mathbf{R}_q,\omega)E_{i'}^\dagger (\mathbf{R}_{q'},\omega')\rangle = 2\pi \delta (\omega-\omega') \langle E_i(\mathbf{R}_q,\omega)E_{i'}^\dagger(\mathbf{R}_{q'},\omega)\rangle, \end{equation} \tag{ A.6 }$$

being

$$\begin{eqnarray} &&\fl \langle E_i(\mathbf{R}_q,\omega)E_{i'}^\dagger(\mathbf{R}_{q'},\omega)\rangle=\frac{\hbar\omega^2}{2\epsilon_0 c^2}\sum_{p,p'}\int\frac{\mathrm{d}^2\mathbf{k}}{(2\pi)^2}\int\frac{\mathrm{d}^2\mathbf{k}'}{(2\pi)^2} {\rm e}^{{\rm i}(\mathbf{k}\cdot\mathbf{r}_q-\mathbf{k}'\cdot\mathbf{r}_{q'})}\nonumber\\ &&\,\times\langle p,\mathbf{k}|\{{\rm e}^{{\rm i}(k_z z_q-k_z^{'*}z_{q'})} [\hat{\boldsymbol{\epsilon}}_p^+(\mathbf{k},\omega)]_i [\hat{\boldsymbol{\epsilon}}_{p'}^{+}(\mathbf{k}',\omega)]_{i'}^{*} [[1+n(\omega,T_{\mathrm{M}})]\nonumber\\ &&\,\times(\mathcal{P}_{-1}^{\mathrm{(pw)}}-\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger} +\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(ew)}} -\mathcal{P}_{-1}^{\mathrm{(ew)}}\mathcal{R}^{\dagger} -\mathcal{T}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{T}^{\dagger})+[1+n(\omega,T_{\mathrm{W}})]\nonumber\\ &&\times (\mathcal{T}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{T}^{\dagger}+\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger})] +[1+n(\omega ,T_{\mathrm{W}})] [{\rm e}^{{\rm i}(k_zz_q+k_z^{'*}z_{q'})} \nonumber\\ &&\times [\hat{\boldsymbol{\epsilon}}_p^+(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}_{p'}^{-}(\mathbf{k}',\omega)]_{i'}^{*}\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}} +{\rm e}^{-{\rm i}(k_zz_q+k_z^{'*}z_{q'})}[\hat{\boldsymbol{\epsilon}}_p^-(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}_{p'}^{+}(\mathbf{k}',\omega)]_{i'}^{*}\nonumber\\ &&\times \mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger}+{\rm e}^{-{\rm i}(k_z z_{q}-k_z^{'*}z_{q'})}[\hat{\boldsymbol{\epsilon}}_p^-(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}_{p'}^{-}(\mathbf{k}',\omega)]_{i'}^{*}\mathcal{P}_{-1}^{\mathrm{(pw)}}]\}|p',\mathbf{k}'\rangle. \end{eqnarray} \tag{ A.7 }$$

We observe that the normally ordered correlation functions 〈E^†_i(R_q,ω)E_i'(R_q',ω')〉 of (A.4) are obtained from the two previous equations by replacing [1 + n(ω,T_i)] with n(ω,T_i) and by taking the complex conjugate. Equation (A.2) finally becomes

$$\begin{equation} \begin{array}{ll} \hspace{-5pc}\displaystyle \gamma_{ii'}^{qq'}(\omega)=\frac{1}{\hbar^2} \langle E_i(\mathbf{R}_q,\omega)E_{i'}^\dagger(\mathbf{R}_{q'},\omega)\rangle, \quad \gamma_{ii'}^{qq'}(-\omega)=\frac{1}{\hbar^2} \langle E_i^\dagger(\mathbf{R}_q,\omega)E_{i'}(\mathbf{R}_{q'},\omega)\rangle, \\ \hspace{-5pc}\displaystyle s_{ii'}^{qq'}(\omega)= \frac{1}{\hbar^2} \mathcal{P}\int_0^{+\infty}\frac{\mathrm{d}\omega'}{2\pi}\Bigg[\frac{\langle E_i(\mathbf{R}_q,\omega')E_{i'}^\dagger(\mathbf{R}_{q'},\omega')\rangle}{\omega-\omega'} +\frac{\langle E_i^\dagger(\mathbf{R}_q,\omega')E_{i'}(\mathbf{R}_{q'},\omega')\rangle}{\omega+\omega'}\Bigg]. \end{array} \end{equation} \tag{ A.8 }$$

Here, we treat explicitly the case when there is no body close to the two emitters. In this case, we have in (A.5) $\mathcal {T}=1$, $\mathcal {R} =0$, or equivalently in (28) ρ_p(k,ω) = 0 and τ_p(k,ω) = 1 [ε(ω) = 1]. It follows that the integral matrices in (37) reduce to [A^qq'(ω)]_ii' = [B^qq'(ω)]_ii' and [C^qq'(ω)]_ii' = [D^qq'(ω)]_ii' = 0, so that [α^qq'_M(ω)]_ii' = 0 and [α^qq'_W(ω)]_ii' (we choose the x-axis along the inter-atomic axis so that in (33) z_q = z_q' = 0 and R_q − R_q' = {r_qq',0,0}):

$$\begin{equation} \fl[\alpha^{qq'}_{\mathrm{W}}(\omega)]_{ii'}= \frac{3c}{4\omega} \sum_p\int_0^{\frac{\omega}{c}}\frac{k\,\mathrm{d}k}{k_z} \mathrm{Re}{ [N_p^{++}(k,\omega)]_{ii'}} = (\alpha^{qq'}_{\mathrm{W}})_\parallel \tilde{e}_\parallel + (\alpha_W)_\perp ( \tilde{e}_{\perp (y)}+\tilde{e}_{\perp (z)}) , \end{equation} \tag{ B.1 }$$

where $\tilde {e}_\parallel$, $\tilde {e}_{\perp (y)}$ and $\tilde {e}_{\perp (z)}$ are, respectively, unit vectors along the parallel and the perpendicular directions to the inter-atomic axis (x) and

$$\begin{equation} (\alpha^{qq'}_{\mathrm{W}})_\parallel = 3\left( \frac{ \sin \tilde{r} }{\tilde{r}^3 }-\frac{ \cos \tilde{r} }{\tilde{r}^2 }\right), \quad (\alpha^{qq'}_{\mathrm{W}})_\perp= \frac{3}{2}\left( \frac{ \sin \tilde{r} }{\tilde{r} }+\frac{\cos \tilde{r} }{\tilde{r}^2 }-\frac{ \sin \tilde{r} }{\tilde{r}^3 } \right), \end{equation} \tag{ B.2 }$$

where $\tilde {\mathbf {r}}=(\mathbf {R}_{q}-\mathbf {R}_{q'})\omega /c$ and $\tilde {r}= |\tilde {\mathbf {r}}|$. By using the two previous equations and (A.8) and (21), Γ^qq'(ω) of (17) can be cast in the form

$$\begin{eqnarray} &&\fl \Gamma^{qq'}(\omega) = \sqrt{\Gamma_0^q(\omega)\Gamma_0^{q'}(\omega) } [1+n(\omega,T_{\mathrm{W}})]\{ (\tilde{\mathbf{d}}_{mn}^{q\,*}\cdot \tilde{{\bf r}} ) (\tilde{\mathbf{d}}_{m'n'}^{q'}\cdot\tilde{{\bf r}}) (\alpha^{qq'}_{\mathrm{W}})_\parallel\nonumber\\ &&+ [ \tilde{\mathbf{d}}_{mn}^{q\,^{*}} \cdot \tilde{\mathbf{d}}_{m'n'}^{q'}- (\tilde{\mathbf{d}}_{mn}^{q\,*}\cdot\tilde{{\bf r}} ) \;(\tilde{\mathbf{d}}_{m'n'}^{q'}\cdot\tilde{{\bf r}})] (\alpha^{qq'}_{\mathrm{W}})_\perp \} \end{eqnarray} \tag{ B.3 }$$

which does not depend anymore on the reference system chosen to derive (B.2). In order to compare the previous result with the known expressions, let us consider the case of two qubits in vacuum with dipoles parallel between them (direction $\tilde {\mathbf {d}}$) with different modulus |d_q| ≠ |d_q'|. In this case, the previous equation reduces to the form (see for instance [25, 26])

$$\begin{eqnarray} &&\fl \Gamma^{qq'}(\omega) = \sqrt{\Gamma_0^q(\omega)\Gamma_0^{q'}(\omega) } [1+n(\omega,T_{\mathrm{W}})] \frac{3}{2} \Biggl\{ [1- (\tilde{\mathbf{d}}\cdot\tilde{\mathbf{r}})^2 ] \frac{\sin \tilde{r} }{\tilde{r} } [1- 3 (\tilde{\mathbf{d}}\cdot\tilde{\mathbf{r}})^2 ] \left( \frac{\cos \tilde{r} }{\tilde{r}^2 }-\frac{\sin \tilde{r} }{\tilde{r}^3 }\right)\Biggr\}.\nonumber\\ \end{eqnarray} \tag{ B.4 }$$

Here, we connect the approaches based on field correlations functions and on Green's function in order to develop the expression for Λ^qq'(ω) of (26). At thermal equilibrium the correlators of the total electromagnetic field outside the body follow from the fluctuation–dissipation theorem [54]

$$\begin{equation} \fl \langle E^{\mathrm{(tot)}}_i(\mathbf{R}_q,\omega) E_{i'}^{\mathrm{(tot)}\dagger}(\mathbf{R}_{q'},\omega')\rangle=2 \pi\delta(\omega-\omega') 2 \hbar [1+n(\omega,T)] \mathrm{Im} \,G_{ii'}(\mathbf{R}_q,\mathbf{R}_{q'},\omega). \end{equation} \tag{ C.1 }$$

In (C.1) G_ii'(R_q,R_q',ω) is the ii' component of the Green function of the system, solution of the differential equation (for two arbitrary points R and R')

$$\begin{equation} \left[\nabla_{\mathbf{R}}\times\nabla_{\mathbf{R}}-\frac{\omega^2}{c^2}\varepsilon(\omega,\mathbf{R})\right]\mathbb{G}(\mathbf{R},\mathbf{R}',\omega)=\frac{\omega^2}{\epsilon_0c^2}\,\mathbb{I}\,\delta(\mathbf{R}-\mathbf{R}') \end{equation} \tag{ C.2 }$$

$\mathbb {I}$ being the identity dyad and ε(ω,R) the dielectric function of the medium. The property (C.1) does not hold in the case of a non-equilibrium configuration. The comparison between (A.6) and (21) at equilibrium T_W = T_M = T, and (C.1) gives

$$\begin{equation} \mathrm{Im}\, G_{ii'} (\mathbf{R}_q,\mathbf{R}_{q'},\omega)=\frac{\omega^3}{3\pi \epsilon_0c^3}\frac{[\alpha_{\mathrm{W}}^{qq'}(\omega)]_{ii'} +[\alpha_{\mathrm{M}}^{qq'}(\omega)]_{ii'}}{2}. \end{equation} \tag{ C.3 }$$

Once we state this connection, which is used in (26), we need to compute the real part of the Green function to develop equation (27). Following appendix C of [47], the ii' component of the Green function for two arbitrary points R = {r,z} and R' = {r',z'} on the right side of the body reads like

$$\begin{eqnarray} &&\fl G_{ii'}(\mathbf{R},\mathbf{R}',\omega)=G^{(0)}_{ii'}(\mathbf{R},\mathbf{R}',\omega)+G^{\mathrm{(R)}}_{ii'}(\mathbf{R},\mathbf{R}',\omega),\nonumber\\ &&\fl G^{(0)}_{ii'}(\mathbf{R},\mathbf{R}',\omega)=\frac{\mathrm{i}\omega^2}{2\epsilon_0c^2}\sum_p\int\frac{\mathrm{d}^2\mathbf{k}}{(2\pi)^2}\exp[\mathrm{i}\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}')]\frac{1}{k_z}[\theta(z-z')[\hat{\boldsymbol{\epsilon}}^+_p(\mathbf{k},\omega)]_i\nonumber\\ &&\times[\hat{\boldsymbol{\epsilon}}^+_p(\mathbf{k},\omega)]_{i'} \exp[\mathrm{i}k_z(z\!-\!z')]+\theta(z'-z)[\hat{\boldsymbol{\epsilon}}^-_p(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}^-_p(\mathbf{k},\omega)]_{i'}\exp[\mathrm{i}k_z(z'\!-\!z)]],\nonumber\\ &&\fl G^{\mathrm{(R)}}_{ii'}(\mathbf{R},\mathbf{R}',\omega)=\frac{\mathrm{i}\omega^2}{2\epsilon_0c^2}\sum_{pp'}\int\frac{\mathrm{d}^2\mathbf{k}}{(2\pi)^2}\int\frac{\mathrm{d}^2\mathbf{k}'}{(2\pi)^2} \exp[\mathrm{i}(\mathbf{k}\cdot\mathbf{r}-\mathbf{k}'\cdot\mathbf{r}')] \nonumber\\ &&\times\frac{1}{k'_z}[\hat{\boldsymbol{\epsilon}}^+_p(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}^{-}_{p'}(\mathbf{k}',\omega)]_{i'} \exp[\mathrm{i}(k_zz+k'_zz')]\langle\mathbf{k},p|\mathcal{R}|\mathbf{k}',p'\rangle, \end{eqnarray} \tag{ C.4 }$$

where θ is the Heaviside step function (θ(x) = 1 for x > 0 and θ(x) = 0 elsewhere) and G_ii'(R,R',ω) has been divided into a free term, G⁽⁰⁾_ii'(R,R',ω), independent of the scattering operators, and a reflected one, G^(R)_ii'(R,R',ω), proportional to $\mathcal {R}$. With regards to the imaginary part of G_ii'(R_q,R_q',ω) it is possible to check starting from (C.4) that equation (C.3) is verified.

Concerning the real part of G_ii'(R_q,R_q',ω), to derive its expression we will make use of the properties of the polarization unit vectors

$$\begin{eqnarray} \fl \hat{\boldsymbol{\epsilon}}_{\mathrm{TE}}^\phi(-\mathbf{k},\omega)&=-\hat{\boldsymbol{\epsilon}}_{\mathrm{TE}}^\phi(\mathbf{k},\omega),\quad\hat{\boldsymbol{\epsilon}}_{\mathrm{TE}}^{-\phi}(\mathbf{k},\omega)=\hat{\boldsymbol{\epsilon}}_{\mathrm{TE}}^\phi(\mathbf{k},\omega), \quad (\hat{\boldsymbol{\epsilon}}_{\mathrm{TE}}^\phi(\mathbf{k},\omega))^{*}=\hat{\boldsymbol{\epsilon}}_{\mathrm{TE}}^\phi(\mathbf{k},\omega),\nonumber\\ \fl \hat{\boldsymbol{\epsilon}}_{\mathrm{TM}}^\phi(-\mathbf{k},\omega)&=\hat{\boldsymbol{\epsilon}}_{\mathrm{TM}}^{-\phi}(\mathbf{k},\omega),\quad (\hat{\boldsymbol{\epsilon}}_{\mathrm{TM}}^\phi(\mathbf{k},\omega))^{*}=\left\{\begin{array}{@{}ll@{}} \hat{\boldsymbol{\epsilon}}_{\mathrm{TM}}^\phi(\mathbf{k},\omega), & k_z\in\mathbb{R}, \\ \hat{\boldsymbol{\epsilon}}_{\mathrm{TM}}^{-\phi}(\mathbf{k},\omega), & k_z \notin \mathbb{R} \end{array}\right. \end{eqnarray} \tag{ C.5 }$$

and of the reciprocity relations of the scattering operators presented in appendix D of [47]

$$\begin{equation} \frac{1}{k^{'*}_z}(-1)^{p+p'}\langle-\mathbf{k}',p'|\mathcal{R}^{\phi\dagger}|-\mathbf{k},p\rangle=\frac{1}{k^{*}_z}\langle\mathbf{k},p|\mathcal{R}^{\phi\dagger}|\mathbf{k}',p'\rangle. \end{equation} \tag{ C.6 }$$

Starting from the free term G⁽⁰⁾_ii' in (C.4), its real part can be written as the sum of two terms coming, respectively, from the propagative and evanescent sectors (a change of variable from k to −k is done in the terms obtained by complex conjugation, we make use of (C.5) and we choose the interatomic axis along the z direction):

$$\begin{eqnarray} &&\fl \mathrm{Re}\, G^{(0)}_{ii'}(\mathbf{R},\mathbf{R}',\omega)_{\mathrm{PW}}=\frac{\mathrm{i }\omega^2}{4\epsilon_0c^2} \frac{1}{(2\pi)^2}\int_{0}^{\frac{\omega}{c}} \frac{k \,\mathrm{d} k}{k_z} [R(k,\omega)]_{ii'} \{ \theta(z-z') [\exp[\mathrm{i} k_z (z-z')] \nonumber\\ &&\, -\exp[-\mathrm{i} k_z (z-z')] ] +\theta(z'-z) [ \exp[\mathrm{i} k_z (z'-z)]-\exp[-\mathrm{i} k_z (z'-z)]]\} ,\nonumber\\ &&\fl \mathrm{Re}\, G^{(0)}_{ii'}(\mathbf{R},\mathbf{R}',\omega)_{\mathrm{EW}}=\frac{\mathrm{i} \omega^2}{4\epsilon_0c^2} \frac{2}{(2\pi)^2} \! \int_{\frac{\omega}{c}}^{+\infty}\! \frac{k\, \mathrm{d} k}{k_z} [R(k,\omega)]_{ii'} \nonumber\\ &&\,\times\{\theta(z-z')\exp[\mathrm{i} k_z (z-z')]+\theta(z'-z)\exp[\mathrm{i} k_z (z'-z)] \} , \end{eqnarray} \tag{ C.7 }$$

where we have used the angular integrals

$$\begin{equation} [R(k,\omega)]_{ii'}=\sum_p\int_0^{2\pi}\mathrm{d}\theta [\hat{\boldsymbol{\epsilon}}^{+(-)}_p(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}^{+(-)}_{p}(\mathbf{k},\omega)]_{i'} , \end{equation} \tag{ C.8 }$$

being the matrix R(k,ω) diagonal with [R(k,ω)]₁₁ = [R(k,ω)]₂₂ = π(2 − c²k²/ω²) and [R(k,ω)]₃₃ = 2πc²k²/ω². The integral in Re G⁽⁰⁾_ii'(R,R',ω)_PW gives two terms, one erasing exactly the integral in Re G⁽⁰⁾_ii'(R,R',ω)_EW, and the second being equal to (we distinguish diagonal elements perpendicular and parallel to the interatomic axis)

$$\begin{equation} \begin{array}{ll} \displaystyle \mathrm{Re}\, G^{(0)}_{\perp}(\mathbf{R},\mathbf{R}',\omega)=\frac{1}{4 \pi}\frac{ \omega^3}{\epsilon_0c^3} \left[\frac{(\tilde{r}^2-1) \cos \tilde{r} - \tilde{r} \sin \tilde{r}}{\tilde{r}^3} \right], \\ \displaystyle \mathrm{Re} \,G^{(0)}_{\parallel}(\mathbf{R},\mathbf{R}',\omega)=\frac{1}{2 \pi}\frac{ \omega^3}{\epsilon_0c^3} \left[\frac{ \cos \tilde{r} + \tilde{r} \sin \tilde{r}}{\tilde{r}^3}\right], \end{array} \end{equation} \tag{ C.9 }$$

where $\tilde {r}= |\mathbf {R}-\mathbf {R}'|\omega /c$ and we named Re G⁽⁰⁾_xx = Re G⁽⁰⁾_yy = Re G⁽⁰⁾_⊥ and Re G⁽⁰⁾_zz = Re G⁽⁰⁾_∥.

Function Λ^qq'(ω) of (27) can be thus decomposed into two parts, Λ^qq'(ω) = Λ^qq'₀(ω) + Λ^qq'_R(ω), connected to G⁽⁰⁾(R_q,R_q',ω) and G^(R)(R_q,R_q',ω), Λ^qq'₀(ω) being

$$\begin{eqnarray} &&\fl \Lambda_0^{qq'}(\omega)= -\sqrt{\Gamma_0^q(\omega)\Gamma_0^{q'}(\omega)} \frac{3}{4} \Biggl\{ [ \tilde{ \mathbf{d}}_{mn}^{q\,*} \cdot \tilde{\mathbf{d}}_{m'n'}^{q'} - (\tilde{\mathbf{d}}_{mn}^{q\,*}\cdot \tilde{{\bf r}} ) \;(\tilde{\mathbf{d}}_{m'n'}^{q'}\cdot \tilde{{\bf r}})] \nonumber\\ &&\times \left[\frac{(\tilde{r}^2-1) \cos \tilde{r} - \tilde{r} \sin \tilde{r}}{\tilde{r}^3} \right] +2 (\tilde{\mathbf{d}}_{mn}^{q\,^{*}}\cdot\tilde{{\bf r}} ) \;(\tilde{\mathbf{d}}_{m'n'}^{q'}\cdot\tilde{{\bf r}}) \left[\frac{ \cos \tilde{r} + \tilde{r} \sin \tilde{r}}{\tilde{r}^3}\right] \Biggr\}, \end{eqnarray} \tag{ C.10 }$$

which has been put in a form which does not depend anymore on the reference system chosen to derive equation (C.9). In the case of two qubits in a vacuum field in the absence of matter with electric dipoles parallel between them (direction $\tilde {\mathbf {d}}$) with |d_q| ≠ |d_q'|, Λ^qq'₀(ω) of (C.10) reduces to the known form [25, 26]

$$\begin{equation} \fl \Lambda_0^{qq'}(\omega)=\sqrt{\Gamma_0^q(\omega)\Gamma_0^{q'}(\omega) } \frac{3}{4} \left\{ [1- 3 (\tilde{\mathbf{d}}\cdot\tilde{\mathbf{r}})^2 ]\left( \frac{\sin \tilde{r} }{\tilde{r}^2 }+\frac{\cos \tilde{r} }{\tilde{r}^3 }\right) -[1- (\tilde{\mathbf{d}}\cdot\tilde{\mathbf{r}} )^2 ] \frac{\cos \tilde{r} }{\tilde{r} } \right\}, \end{equation} \tag{ C.11 }$$

where we used $\tilde {\mathbf {d}}\cdot \tilde{\mathbf{r}}=\tilde {d}_z$ and $\tilde {d}_x^2+\tilde {d}_y^2=1-(\tilde {d}_z)^2$.

We now consider the remaining part of the Green function in (C.4). By making a change of variable in the terms obtained by complex conjugation, (k,k') → (− k,−k'), and using the reciprocity relations of scattering operators in (C.6) and the properties of the polarization unit vectors (C.5), one can obtain

$$\begin{eqnarray} &&\fl \mathrm{Re}\, G^{(R)}_{ii'}(\mathbf{R},\mathbf{R}',\omega) =\frac{\mathrm{i} \omega^2}{4\epsilon_0c^2}\sum_{p,p'}\int\frac{\mathrm{d}^2\mathbf{k}}{(2\pi)^2}\int\frac{\mathrm{d}^2\mathbf{k}'}{(2\pi)^2}{\rm e}^{{\rm i}(\mathbf{k}\cdot\mathbf{r}-\mathbf{k}'\cdot\mathbf{r}')} \langle p,\mathbf{k}|\{{\rm e}^{{\rm i}(k_zz+k_z^{'*}z')}\nonumber\\ &&\times [\hat{\boldsymbol{\epsilon}}_p^+(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}_{p'}^{-}(\mathbf{k}',\omega)]_{i'}^{*} \mathcal{R}\mathcal{P}_{-1}^{\mathrm{(pw)}} -{\rm e}^{-{\rm i}(k_zz+k_z^{'*}z')}[\hat{\boldsymbol{\epsilon}}_p^-(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}_{p'}^{+}(\mathbf{k}',\omega)]_{i'}^{*}\mathcal{P}_{-1}^{\mathrm{(pw)}}\mathcal{R}^{\dagger}\nonumber\\ &&+ {\rm e}^{{\rm i}(k_z z-k_z^{'*}z')}[\hat{\boldsymbol{\epsilon}}_p^+(\mathbf{k},\omega)]_i[\hat{\boldsymbol{\epsilon}}_{p'}^{+}(\mathbf{k}',\omega)]_{i'}^{*} (\mathcal{R}\mathcal{P}_{-1}^{\mathrm{(ew)}}+\mathcal{P}_{-1}^{\mathrm{(ew)}}\mathcal{R}^{\dagger})\}|p',\mathbf{k}'\rangle. \end{eqnarray} \tag{ C.12 }$$