In this section, we discuss how the different parameters in the master equation can be obtained from the nature of surface polaritons for a system of two dielectric half-spaces separated by a vacuum gap and subsequently calculate the heat transfer for this system. Consider first a system of two coupled oscillators with natural frequency *ω*_{0}, half linewidth *γ*, coupling constant $\xi {\omega}_{0}^{2},$ and displacements *x*_{1}(*ω*), *x*_{2}(*ω*), where ${x}_{1/2}\mathrm{(}\omega \mathrm{)}={\displaystyle \int {e}^{i\omega t}{x}_{1/2}\mathrm{(}t\mathrm{)}dt},$ which are related as

$$\begin{array}{l}{\omega}^{2}{x}_{1}\mathrm{(}\omega \mathrm{)}={\omega}_{0}^{2}{x}_{1}\mathrm{(}\omega \mathrm{)}-2i\omega \gamma {x}_{1}\mathrm{(}\omega \mathrm{)}+\xi {\omega}_{0}^{2}{x}_{2}\mathrm{(}\omega \mathrm{)}\\ {\omega}^{2}{x}_{2}\mathrm{(}\omega \mathrm{)}={\omega}_{0}^{2}{x}_{2}\mathrm{(}\omega \mathrm{)}-2i\omega \gamma {x}_{2}\mathrm{(}\omega \mathrm{)}+\xi {\omega}_{0}^{2}{x}_{1}\mathrm{(}\omega \mathrm{)}\end{array}$$(7)The eigenfrequencies for such a system are given by

$${\omega}_{\pm}=\sqrt{{\omega}_{0}^{2}\mathrm{(}1\pm \xi \mathrm{)}-{\gamma}^{2}}-i\gamma \mathrm{;}$$(8)To find the equivalent coupling constant *ξ* between surface modes, we find the eigenmodes of the configuration of two flat surfaces separated by a vacuum gap using Maxwell’s equations and compare the expression to that of Eq. (8) for two coupled HOs. Consider two half-spaces separated by a gap of width *d* occupying the regions *z*<−*d*/2 and *z*>*d*/2 as shown in . From the dispersion relation for surface polaritons *ω*(**k**) and close to the surface polariton frequency, the surface mode is seen to acquire an electrostatic character with |**k**|→∞ [20]. In this electrostatic limit, satisfying the continuity of perpendicular component of displacement field gives the condition for surface modes as

$$\epsilon \mathrm{(}{\omega}_{\pm}\mathrm{)}=\mathrm{(}\begin{array}{c}\mathrm{(}1-{e}^{\mathrm{|}k\mathrm{|}d}\mathrm{)}\text{/}\mathrm{(}1+{e}^{\mathrm{|}k\mathrm{|}d}\mathrm{)}\\ \mathrm{(}1+{e}^{\mathrm{|}k\mathrm{|}d}\mathrm{)}\text{/}\mathrm{(}1-{e}^{\mathrm{|}k\mathrm{|}d}\mathrm{)}\end{array}$$(9)where *ε*(*ω*) is the dielectric function of the half-space. Taking the Lorentz model for the dielectric function:

$$\epsilon \mathrm{(}\omega \mathrm{)}={\epsilon}_{\infty}\mathrm{(}1+\frac{{\omega}_{L}^{2}-{\omega}_{T}^{2}}{{\omega}_{T}^{2}-{\omega}^{2}-i\omega \Gamma}\mathrm{)}$$(10)and solving for *ω* in Eq. (9) in the low-loss limit Γ → 0, and with some algebra, we obtain the eigenfrequencies of the coupled surface modes to be of the form:

$${\omega}_{\pm}=\sqrt{{{\omega}^{\prime}}_{0}{}^{2}\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)}\mathrm{(}1\pm {\xi}^{\prime}\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)}\mathrm{)}-\frac{{\Gamma}^{2}}{4}}-i\frac{\Gamma}{2}$$(11)where

$${{\omega}^{\prime}}_{0}\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)}=\sqrt{\frac{2{e}^{-\mathrm{|}k\mathrm{|}d}\mathrm{(}{\epsilon}_{\infty}\mathrm{(}{\omega}_{L}^{2}+{\omega}_{T}^{2}\mathrm{)}\text{\hspace{0.17em}}\mathrm{cosh}\mathrm{(}\mathrm{|}k\mathrm{|}d\mathrm{)}+\mathrm{(}{\epsilon}_{\infty}^{2}{\omega}_{L}^{2}+{\omega}_{T}^{2}\mathrm{)}\text{\hspace{0.17em}}\mathrm{sinh}\mathrm{(}\mathrm{|}k\mathrm{|}d\mathrm{)}\mathrm{)}}{{\mathrm{(}{\epsilon}_{\infty}+1\mathrm{)}\mathrm{}}^{2}-{\mathrm{(}{\epsilon}_{\infty}-1\mathrm{)}\mathrm{}}^{2}{e}^{-2|k\mathrm{|}d}}}$$(12)and

$${\xi}^{\prime}\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)}=\frac{{\epsilon}_{\infty}\mathrm{(}{\omega}_{L}^{2}-{\omega}_{T}^{2}\mathrm{)}}{{\epsilon}_{\infty}\mathrm{(}{\omega}_{L}^{2}+{\omega}_{T}^{2}\mathrm{)}\text{\hspace{0.17em}}\mathrm{cosh}\mathrm{(}\mathrm{|}k\mathrm{|}d\mathrm{)}+\mathrm{(}{\epsilon}_{\infty}^{2}{\omega}_{L}^{2}+{\omega}_{T}^{2}\mathrm{)}\text{\hspace{0.17em}}\mathrm{sinh}\mathrm{(}\mathrm{|}k\mathrm{|}d\mathrm{)}}\mathrm{;}$$(13)In the limit of large gaps |**k**|*d*→∞, ${{\omega}^{\prime}}_{0}\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)}\approx {\omega}_{0},$ where *ω*_{0} is the surface polariton frequency for a single half-space given by ${\omega}_{0}=\sqrt{\mathrm{(}{\epsilon}_{\infty}{\omega}_{L}^{2}+{\omega}_{T}^{2}\mathrm{)}\text{/}\mathrm{(}{\epsilon}_{\infty}+1\mathrm{)}\mathrm{}}$ and the coupling parameter *ξ*′(**k**, *d*) ≈ 0. For the case when *ε*_{∞}=1, these expressions reduce to the simple forms: ${{\omega}^{\prime}}_{0}=\sqrt{\mathrm{(}{\omega}_{L}^{2}+{\omega}_{T}^{2}\mathrm{)}\text{/}2}$ and ${\xi}^{\prime}\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)}={e}^{-\mathrm{|}k\mathrm{|}d}\mathrm{(}{\omega}_{L}^{2}-{\omega}_{T}^{2}\mathrm{)}\text{/}\mathrm{(}{\omega}_{L}^{2}+{\omega}_{T}^{2}\mathrm{)}$[15]. This method of finding the eigenmodes of the surface modes in the electrostatic limit is well known and has been previously employed, for example, in the context of deriving the van der Waals force [12] and the steady-state van der Waals heat transfer [15]. By comparing Eqs. (11) and (8), we can relate the parameters of the HO model with those of the surface modes:

$${\omega}_{0}\to {{\omega}^{\prime}}_{0}\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)};\text{\hspace{0.17em}}\xi \to {\xi}^{\prime}\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)};\text{\hspace{0.17em}}\gamma \to \Gamma \text{/}2$$(14)We can thus use the results from Section 2 to find the heat transfer between the two flat surfaces due to the coupling between the surface modes. As noted earlier, we add the contributions from all the HO modes (labeled by **k**), observing that heat transfer does not result in change in momentum and that only surface modes of the same in-plane wave-vector component interact across the two half-spaces. The heat flux between two half-spaces *P*(*t*, *d*) (W m^{−2}) is then given by

$$P\mathrm{(}t\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)}={\displaystyle \int \frac{1}{4{\pi}^{2}}\u3008P\u3009\text{\hspace{0.17em}}{d}^{2}k}$$(15)where 〈*P*〉(*t*, **k**, *d*) is the rate of heat transfer (W) between two HOs obtained from Eq. (4). The relation between *g* in Eq. (4) and ξ in Eq. (8) can be derived by comparing the interaction Hamiltonians in the quantum mechanical picture and the classical model. In the classical picture, the potential energy of the Hamiltonian of the coupled HO system is given by

${H}_{\text{PE}}=\frac{1}{2}{\omega}_{0}^{2}{x}_{1}^{2}+\frac{1}{2}{\omega}_{0}^{2}{x}_{2}^{2}+{\omega}_{0}^{2}\xi {x}_{1}{x}_{2}$

where the interaction term is given from ${H}_{\text{I}}={\omega}_{0}^{2}\xi {x}_{1}{x}_{2}.$ The equivalent quantum mechanical model can be obtained by substituting the displacements *x*_{1} and *x*_{2} in terms of the commutator operators *a*, *a*^{†}, *b*, and *b*^{†} using the standard relations ${x}_{1}=\mathrm{(}a+{a}^{\u2020}\mathrm{)}\sqrt{\hslash \text{/}\mathrm{(}2{\omega}_{0}\mathrm{)}},$ and ${x}_{2}=\mathrm{(}b+{b}^{\u2020}\mathrm{)}\sqrt{\hslash \text{/}\mathrm{(}2{\omega}_{0}\mathrm{)}}.$ Including only terms resulting in photon exchange, we get

${\widehat{H}}_{\text{I}}={\omega}_{0}\xi \mathrm{(}{a}^{\u2020}b+a{b}^{\u2020}\mathrm{)}\frac{\hslash}{2}.$

Comparing to the interaction Hamiltonian in the quantum mechanical picture from Eq. (1) (for a particular **k** value) gives us *g*=*ω*_{0}*ξ*/2. Substituting for the values of *ω*_{0} and *ξ* from Eq. (14) gives

$$g\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)}=\frac{{{\omega}^{\prime}}_{0}\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)}{\xi}^{\prime}\mathrm{(}k\mathrm{,}\text{\hspace{0.17em}}d\mathrm{)}}{2}$$(16)