Near-field thermal radiative transfer between two coated spheres

TL;DR

This paper derives an expression for near-field thermal radiative transfer between coated spheres, revealing novel spectral features and conditions where layered spheres outperform homogeneous ones in heat transfer.

Contribution

It introduces a new analytical model for heat transfer between coated spheres and uncovers unique spectral phenomena and enhanced conductance due to layering effects.

Findings

Spectral peaks at non-resonant frequencies due to coatings

Near-field selective emission observed in metallic spheres with polar coatings

Layered spheres can surpass homogeneous spheres in total conductance

Abstract

In this work, we present an expression for the near-field thermal radiative transfer between two spheres with an arbitrary numbers of coatings. We numerically demonstrate that the spectrum of heat transfer between layered spheres exhibits novel features due to the newly introduced interfaces between coatings and cores. These features include broad super-Planckian peaks at non-resonant frequencies and near-field selective emission between metallic spheres with polar material coatings. Spheres with cores and coatings of two different polar materials are also shown to exceed the total conductance of homogeneous spheres in some cases.

Equations

$$Q_{A\to B}$$

$$\mathbf{P}(\tilde{\mathbf{r}})$$

$$\mathbf{E}(\tilde{\mathbf{r}})$$

$$\mathbf{H}(\tilde{\mathbf{r}})$$

$$⟨{\mathbf{J}}_{\mathbf{p}}^{\mathbf{e}}(\mathbf{r}){\mathbf{J}}_{\mathbf{q}}^{\mathbf{e}\ast }(\mathbf{r})⟩$$

$$⟨{\mathbf{J}}_{\mathbf{p}}^{\mathbf{m}}(\mathbf{r}){\mathbf{J}}_{\mathbf{q}}^{\mathbf{m}\ast }(\mathbf{r})⟩$$

$$⟨{\mathbf{J}}_{\mathbf{p}}^{\mathbf{e}}(\mathbf{r}){\mathbf{J}}_{\mathbf{q}}^{\mathbf{m}\ast }(\mathbf{r})⟩$$

$$G(\omega ,T)$$

$$G_t(T)={\int }_0^{\infty }\frac{d\omega }{2\pi }G(\omega ,T)={\int }_0^{\infty }d\lambda G(\lambda ,T),$$

$$\displaystyle\displaystyle T_{A\rightarrow B}^{e}(\omega)=2\Re\text{Tr}\oint\limits_{S_{B}}d\boldsymbol{r}\oint\limits_{S_{A}}d\boldsymbol{\widetilde{r}}\Bigg{[}$$

$$\frac{{\omega }^2}{c^2}\left[\hat{\mathbf{n}}(\mathbf{r})\times {\mu }_C{\overline{\overline{\mathbf{G}}}}_e(\mathbf{r},\tilde{\mathbf{r}})\right]\cdot {\left[\hat{\mathbf{n}}(\tilde{\mathbf{r}})\times {\epsilon }_C{\overline{\overline{\mathbf{G}}}}_m^T(\mathbf{r},\tilde{\mathbf{r}})\right]}^{\ast }$$

$$\displaystyle\displaystyle+\left[\boldsymbol{\widehat{n}}(\boldsymbol{r})\times\overline{\overline{\boldsymbol{G}}}_{M}(\boldsymbol{r},\boldsymbol{\widetilde{r}})\right]\cdot\left[\boldsymbol{\widehat{n}}(\boldsymbol{\widetilde{r}})\times\overline{\overline{\boldsymbol{G}}}_{E}^{T}(\boldsymbol{r},\boldsymbol{\widetilde{r}})\right]^{*}\Bigg{]},$$

$${\overline{\overline{\mathbf{G}}}}_e(\mathbf{r};\tilde{\mathbf{r}})=i{k}_C\sum _{m=-\infty }^{\infty }\sum _{l=\tilde{m}}^{\infty }\sum _{\nu =\tilde{m}}^{\infty }(-1{)}^m$$

$$\displaystyle\displaystyle\times\quad\left\{\begin{array}[]{r}\left[\begin{array}[]{r}D_{\nu,m}^{lM}\left[\boldsymbol{M}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(M)}(b)\boldsymbol{M}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\\[8.0pt] +D_{\nu,m}^{lN}\left[\boldsymbol{N}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(N)}(b)\boldsymbol{N}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\end{array}\right]\left[\begin{array}[]{l}\boldsymbol{M}_{l,-m}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})+R_{l}^{(M)}(a)\boldsymbol{M}_{l,-m}^{(3)}(k_{C}\boldsymbol{\widetilde{r}_{A}})\end{array}\right]\\[20.0pt] +\left[\begin{array}[]{r}J_{\nu,m}^{lM}\left[\boldsymbol{M}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(M)}(b)\boldsymbol{M}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\\[8.0pt] +J_{\nu,m}^{lN}\left[\boldsymbol{N}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(N)}(b)\boldsymbol{N}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\end{array}\right]\left[\begin{array}[]{l}\boldsymbol{N}_{l,-m}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})+R_{l}^{(N)}(a)\boldsymbol{N}_{l,-m}^{(3)}(k_{C}\boldsymbol{\widetilde{r}_{A}})\end{array}\right]\end{array}\right\},$$

$${\overline{\overline{\mathbf{G}}}}_E(\mathbf{r};\tilde{\mathbf{r}})=i{k}_C^2\sum _{m=-\infty }^{\infty }\sum _{l=\tilde{m}}^{\infty }\sum _{\nu =\tilde{m}}^{\infty }(-1{)}^m$$

$$\displaystyle\displaystyle\times\quad\left\{\begin{array}[]{r}\left[\begin{array}[]{r}D_{\nu,m}^{lN}\left[\boldsymbol{M}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(N)}(b)\boldsymbol{M}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\\[8.0pt] +D_{\nu,m}^{lM}\left[\boldsymbol{N}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(M)}(b)\boldsymbol{N}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\end{array}\right]\left[\begin{array}[]{l}\boldsymbol{M}_{l,-m}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})+R_{l}^{(M)}(a)\boldsymbol{M}_{l,-m}^{(3)}(k_{C}\boldsymbol{\widetilde{r}_{A}})\end{array}\right]\\[20.0pt] +\left[\begin{array}[]{r}J_{\nu,m}^{lN}\left[\boldsymbol{M}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(N)}(b)\boldsymbol{M}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\\[8.0pt] +J_{\nu,m}^{lM}\left[\boldsymbol{N}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(M)}(b)\boldsymbol{N}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\end{array}\right]\left[\begin{array}[]{l}\boldsymbol{N}_{l,-m}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})+R_{l}^{(N)}(a)\boldsymbol{N}_{l,-m}^{(3)}(k_{C}\boldsymbol{\widetilde{r}_{A}})\end{array}\right]\end{array}\right\},$$

$${\overline{\overline{\mathbf{G}}}}_m(\mathbf{r};\tilde{\mathbf{r}})=i{k}_C\sum _{m=-\infty }^{\infty }\sum _{l=\tilde{m}}^{\infty }\sum _{\nu =\tilde{m}}^{\infty }(-1{)}^m$$

$$\displaystyle\displaystyle\times\quad\left\{\begin{array}[]{r}\left[\begin{array}[]{r}J_{\nu,m}^{lN}\left[\boldsymbol{M}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(N)}(b)\boldsymbol{M}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\\[8.0pt] +J_{\nu,m}^{lM}\left[\boldsymbol{N}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(M)}(b)\boldsymbol{N}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\end{array}\right]\left[\begin{array}[]{l}\boldsymbol{M}_{l,-m}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})+R_{l}^{(N)}(a)\boldsymbol{M}_{l,-m}^{(3)}(k_{C}\boldsymbol{\widetilde{r}_{A}})\end{array}\right]\\[20.0pt] +\left[\begin{array}[]{r}D_{\nu,m}^{lN}\left[\boldsymbol{M}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(N)}(b)\boldsymbol{M}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\\[8.0pt] +D_{\nu,m}^{lM}\left[\boldsymbol{N}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(M)}(b)\boldsymbol{N}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\end{array}\right]\left[\begin{array}[]{l}\boldsymbol{N}_{l,-m}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})+R_{l}^{(M)}(a)\boldsymbol{N}_{l,-m}^{(3)}(k_{C}\boldsymbol{\widetilde{r}_{A}})\end{array}\right]\end{array}\right\},$$

$${\overline{\overline{\mathbf{G}}}}_M(\mathbf{r};\tilde{\mathbf{r}})=i{k}_C^2\sum _{m=-\infty }^{\infty }\sum _{l=\tilde{m}}^{\infty }\sum _{\nu =\tilde{m}}^{\infty }(-1{)}^m$$

$$\displaystyle\displaystyle\times\left\{\begin{array}[]{r}\left[\begin{array}[]{r}J_{\nu,m}^{lM}\left[\boldsymbol{M}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(M)}(b)\boldsymbol{M}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\\[8.0pt] +J_{\nu,m}^{lN}\left[\boldsymbol{N}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(N)}(b)\boldsymbol{N}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\end{array}\right]\left[\begin{array}[]{l}\boldsymbol{M}_{l,-m}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})+R_{l}^{(N)}(a)\boldsymbol{M}_{l,-m}^{(3)}(k_{C}\boldsymbol{\widetilde{r}_{A}})\end{array}\right]\\[20.0pt] +\left[\begin{array}[]{r}D_{\nu,m}^{lM}\left[\boldsymbol{M}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(M)}(b)\boldsymbol{M}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\\[8.0pt] +D_{\nu,m}^{lN}\left[\boldsymbol{N}_{\nu m}^{(1)}(k_{C}\boldsymbol{r_{B}})+R_{\nu}^{(N)}(b)\boldsymbol{N}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})\right]\end{array}\right]\left[\begin{array}[]{l}\boldsymbol{N}_{l,-m}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})+R_{l}^{(M)}(a)\boldsymbol{N}_{l,-m}^{(3)}(k_{C}\boldsymbol{\widetilde{r}_{A}})\end{array}\right]\end{array}\right\},$$

$$C_{n,m}^{lM}$$

$$C_{n,m}^{lN}$$

$$\displaystyle\displaystyle\begin{split}D_{n,m}^{lM}&-\sum\limits_{\nu=\widetilde{m}}^{\infty}\left[\begin{array}[]{r}C_{\nu,m}^{lM}R_{\nu}^{(M)}(a)A_{nm}^{\nu m}(k_{C}d_{AB})\\ +C_{\nu,m}^{lN}R_{\nu}^{(N)}(a)B_{nm}^{\nu m}(k_{C}d_{AB})\end{array}\right]\\ &=A_{nm}^{lm}(k_{C}d_{AB}),\end{split}$$

$$\displaystyle\displaystyle\begin{split}D_{n,m}^{lN}&-\sum\limits_{\nu=\widetilde{m}}^{\infty}\left[\begin{array}[]{r}C_{\nu,m}^{lN}R_{\nu}^{(N)}(a)A_{nm}^{\nu m}(k_{C}d_{AB})\\ +C_{\nu,m}^{lM}R_{\nu}^{(M)}(a)B_{nm}^{\nu m}(k_{C}d_{AB})\end{array}\right]\\ &=B_{nm}^{lm}(k_{C}d_{AB}),\end{split}$$

$$F_{n,m}^{lM}$$

$$F_{n,m}^{lN}$$

$$\displaystyle\displaystyle\begin{split}J_{n,m}^{lM}&-\sum\limits_{\nu=\widetilde{m}}^{\infty}\left[\begin{array}[]{r}F_{\nu,m}^{lM}R_{\nu}^{(M)}(a)A_{nm}^{\nu m}(k_{C}d_{AB})\\ +F_{\nu,m}^{lN}R_{\nu}^{(N)}(a)B_{nm}^{\nu m}(k_{C}d_{AB})\end{array}\right]\\ &=B_{nm}^{lm}(k_{C}d_{AB}),\end{split}$$

$$\displaystyle\displaystyle\begin{split}J_{n,m}^{lN}&-\sum\limits_{\nu=\widetilde{m}}^{\infty}\left[\begin{array}[]{r}F_{\nu,m}^{lN}R_{\nu}^{(N)}(a)A_{nm}^{\nu m}(k_{C}d_{AB})\\ +F_{\nu,m}^{lM}R_{\nu}^{(M)}(a)B_{nm}^{\nu m}(k_{C}d_{AB})\end{array}\right]\\ &=A_{nm}^{lm}(k_{C}d_{AB}),\end{split}$$

$$T_{A\to B}^e\left(\omega \right)={\left(k_{C}a\right)}^2{\left(k_{C}b\right)}^2\sum _{m=-\infty }^{\infty }\sum _{l=\tilde{m}}^{\infty }\sum _{\nu =\tilde{m}}^{\infty }$$

$$\displaystyle\displaystyle\times\left[\begin{array}[]{r}\left[\epsilon_{\nu}^{(M)}(b)\left|D_{\nu,m}^{lM}\right|^{2}+\epsilon_{\nu}^{(N)}(b)\left|D_{\nu,m}^{lN}\right|^{2}\right]\epsilon_{l}^{(M)}(a)\\ \\ +\left[\epsilon_{\nu}^{(M)}(b)\left|J_{\nu,m}^{lM}\right|^{2}+\epsilon_{\nu}^{(N)}(b)\left|J_{\nu,m}^{lN}\right|^{2}\right]\epsilon_{l}^{(N)}(a)\end{array}\right],$$

$${ϵ}_{\nu }^{(P)}(w)$$

$$ϵ\left(\omega \right)=\sum _{m=-\infty }^{\infty }\sum _{l=\tilde{m}}^{\infty }\left[{ϵ}_l^{(M)}(w)+{ϵ}_l^{(N)}(w)\right].$$

$$G_{BB}(\lambda ,T)$$

$$G_{t,BB}(T)$$

$$CSC(\lambda )$$

$${\mathbf{M}}_{lm}^{(p)}(k\mathbf{r})=$$

$${\mathbf{N}}_{lm}^{(p)}(k\mathbf{r})=$$

$$+\frac{\sqrt{l(l+1)}}{kr}{z}_l^{(p)}(kr){\mathbf{V}}_{lm}^{(1)}(\theta ,\varphi ),$$

$${\mathbf{V}}_{lm}^{(1)}(\theta ,\varphi )$$

$${\mathbf{V}}_{lm}^{(2)}(\theta ,\varphi )$$

$${\mathbf{V}}_{lm}^{(3)}(\theta ,\varphi )$$

$$Y_{lm}(\theta ,\varphi )$$

Full text

I Introduction

Optical metamaterials are a class of artificial materials which exhibit electromagnetic behavior not otherwise observed in nature, such as negative refractive index,Smith et al. (2004); Zhang et al. (2005); Soukoulis et al. (2007); Valentine et al. (2008) cloaking, Alù and Engheta (2005); Schurig et al. (2006); Cai et al. (2007); Valentine et al. (2009) and superlensingPendry (2000); Fang et al. (2005); Smolyaninov et al. (2007); Zhang and Liu (2008) to name a few. Of particular interest are hyperbolic metamaterials (HMMs), those whose permissible wave-vector components form a hyperbolic isofrequency surface instead of the spherical surface found in typical isotropic materials. The simplest means of achieving HMM behavior is through layering different isotropic materials. When the layer thicknesses are much smaller than the free-space wavelength of light propagating through them, the entire nanocomposite can be viewed as a homogeneous material with hyperbolic effective optical properties.Halevi et al. (1999); Smith and Schurig (2003)

Hyperbolic metamaterials have found many uses in the field of near-field thermal radiative transfer. Appropriately designed HMMs have demonstrated the ability to tailor the spectrum of radiative transfer and to achieve heat transfer beyond that of Planck’s blackbody limit.Francoeur et al. (2011); Liu et al. (2011); Mason et al. (2011); Biehs et al. (2012); Guo et al. (2012); Guo and Jacob (2013); Liu and Shen (2013) To date, this has been achieved mostly by using layered planar surfaces.Biehs (2007); Fu and Tan (2009); Svetovoy et al. (2012) That configuration is attractive because the analytic solution to near-field thermal radiative transfer between two semi-infinite half spaces is well known, relatively straightforward to compute, and easily generalizes to include layered media.Polder and Van Hove (1971); Francoeur and Pinar Mengüç (2008); Francoeur et al. (2009); Song et al. (2016)

Determining analytic solutions to more complicated geometries opens up additional avenues of investigation for HMMs. Two geometries of interest are sphere-sphere and sphere-plane configurations. The formula for heat transfer between two homogeneous spheres has been determinedNarayanaswamy and Chen (2008); Mackowski and Mishchenko (2008); Krüger et al. (2012) and can be used to approximate the sphere-plane configurationOtey and Fan (2011) in the limit that one sphere is much greater than the other.Sasihithlu and Narayanaswamy (2014) The solution method in Ref. Narayanaswamy and Chen, 2008 cannot be easily extended to include coated spheres. Though the formalisms used in Refs. Mackowski and Mishchenko, 2008; Krüger et al., 2012, which are based on the T-matrix method, Waterman (1965); Peterson and Ström (1974) can in principle be used to calculate near-field radiation between coated spheres, the authors did not apply them for that purpose.

In this work, we present an expression for the near-field thermal radiative transfer between two spheres, either of which may have any number of coatings. The expression is derived using the framework of fluctuational electrodynamics and a Green’s function formalism. The key advance in this work is evaluating the radiative transfer between spheres using surface integrals instead of volume integrals, which permits the same formalism and numerical code to be used for coated as well as uncoated spheres, irrespective of the number of coatings. This approach is well suited for near-field radiative transfer analysis of coated spheres because, as in the case of two homogeneous spheres, the main geometric parameters of interest are still the center-to-center distance between the spheres (or alternatively the minimum gap between spheres) and the overall dimensions of each sphere, not the details of the coatings themselves. The details of the coatings on a sphere are encoded in the effective Mie reflection coefficients of the various vector spherical waves at the interface between the outer-most coating of that sphere and the intervening medium (usually vacuum).

We will show that coated spheres have some advantages over homogeneous spheres. Coated spheres exhibit broad, super-Planckian peaks (peaks exceeding that of blackbodies) which are due neither to surface phonon polariton nor surface plasmon polariton resonances and are not commonly observed in homogeneous spheres (see Fig. 2). We will also show that coated spheres with silver cores and silica coatings transition from silica-like behavior at gaps much smaller than the coating thickness to silver-like behavior for larger gaps (see Fig. 3). Last, we present coated spheres whose total conductance exceeds that of homogeneous spheres of either of its constitutive components (see Fig. 4). These advantages make coated spheres a promising platform for research and experimentation.

The structure of the paper is as follows: In Sec. II, the geometry of the two sphere problem is described. In Sec. III, the expression for thermal conductance is derived using the fluctuation-dissipation theorem, resulting in an expression which includes a transmissivity function for energy transfer. In Sec. IV, the transmissivity function is evaluated from the dyadic Green’s functions for the two coated spheres. Finally, in Sec. V, numerical results are presented for the conductance between two coated spheres.

II Geometry

The geometry of the problem is shown in Fig. 1. Two spheres, denoted $\displaystyle A$ and $\displaystyle B$, are composed of an arbitrary number of coatings, $\displaystyle N$ and $\displaystyle M$ coatings, respectively, atop a core. As shown in Fig 1(B), the outer radius of any coating $\displaystyle\rho$ ($\displaystyle 0<\rho\leq N$ for sphere $\displaystyle A$ and $\displaystyle 0<\rho\leq M$ for sphere $\displaystyle B$) in sphere $\displaystyle s$ ($\displaystyle s=A$ or $\displaystyle B$) is given by $\displaystyle r_{\rho,s}$ and the outer radius of the core is given by $\displaystyle r_{0,s}$. For simplicity, the outermost radii of spheres $\displaystyle A$ and $\displaystyle B$ are denoted $\displaystyle r_{N,A}=a$ and $\displaystyle r_{M,B}=b$, respectively. Though only the internal structure of sphere $\displaystyle A$ is shown in Fig. 1(B), sphere $\displaystyle B$ is similar, but with outer radius $\displaystyle b$ on layer $\displaystyle M$. The exterior region, denoted $\displaystyle C$, is vacuum.

As shown in Fig 1(A), a coordinate system is placed with its origin at the center of sphere $\displaystyle A$, henceforth referred to as the $\displaystyle A$-coordinate system. Similarly, a second coordinate system is placed with its origin at the center of sphere $\displaystyle B$ and referred to as the $\displaystyle B$-coordinate system. The two coordinate systems are oriented such that their $\displaystyle z$ axes are aligned along the vector connecting the origin of the $\displaystyle A$-coordinate system to that of the $\displaystyle B$-coordinate system. The spheres are separated by a center-to-center gap $\displaystyle d$. The minimum gap separating spheres $\displaystyle A$ and $\displaystyle B$ is $\displaystyle D=d-a-b$.

The same location may be represented by the position vector $\displaystyle\boldsymbol{r_{A}}$ in the $\displaystyle A$-coordinate system or $\displaystyle\boldsymbol{r_{B}}$ in the $\displaystyle B$-coordinate system. In order to differentiate between two different locations, we will use $\displaystyle\boldsymbol{r}$ and $\displaystyle\boldsymbol{\widetilde{r}}$, the need for which arises from the use of Green’s functions.

III Mathematical formulation

The radiative heat transfer from object A to object B, $\displaystyle Q_{A\rightarrow B}$, is given by

[TABLE]

where $\displaystyle\widehat{\boldsymbol{n}}(\boldsymbol{\widetilde{r}})$ and $\displaystyle\boldsymbol{P}(\boldsymbol{\widetilde{r}})$ are the unit outward normal of surface $\displaystyle B$ and the Poynting vector, respectively, at the location $\displaystyle\boldsymbol{\widetilde{r}}$. The Poynting vector is related to the cross-spectral density of the components of the electric and magnetic fields and is defined as

[TABLE]

where $\displaystyle\omega$ is the angular frequency, $\displaystyle\left<\cdot\right>$ denotes an ensemble average, $\displaystyle\left(\cdot\right)^{*}$ is the complex conjugate, and $\displaystyle\boldsymbol{E}(\boldsymbol{\widetilde{r}})$ and $\displaystyle\boldsymbol{H}(\boldsymbol{\widetilde{r}})$ are the Fourier-transformed electric and magnetic fields, respectively. The frequency dependence of $\displaystyle\boldsymbol{E}(\boldsymbol{\widetilde{r}})$, $\displaystyle\boldsymbol{H}(\boldsymbol{\widetilde{r}})$, and other fields and Green’s functions are suppressed for ease of notation. Inside each sphere, $\displaystyle\boldsymbol{E}(\boldsymbol{\widetilde{r}})$ and $\displaystyle\boldsymbol{H}(\boldsymbol{\widetilde{r}})$ are given by

[TABLE]

where $\displaystyle\boldsymbol{p}(\boldsymbol{r})=i\omega\mu_{0}\mu\boldsymbol{J^{e}}(\boldsymbol{r})$, $\displaystyle\boldsymbol{m}(\boldsymbol{r})=i\omega\varepsilon_{0}\varepsilon\boldsymbol{J^{m}}(\boldsymbol{r})$, $\displaystyle i$ is the imaginary unit, $\displaystyle\mu_{0}$ and $\displaystyle\varepsilon_{0}$ are the permeability and permittivity of free space, $\displaystyle\mu$ and $\displaystyle\varepsilon$ are the relative permeability and permittivity, and $\displaystyle\overline{\overline{\boldsymbol{G}}}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ is the dyadic Green’s function (DGF) relating sources and fields at $\displaystyle\boldsymbol{r}$ and $\displaystyle\boldsymbol{\widetilde{r}}$. The subscripts on the DGFs denote the electric ($\displaystyle e$) and magnetic ($\displaystyle m$) variants, and we define $\displaystyle\overline{\overline{\boldsymbol{G}}}_{E}(\boldsymbol{r};\boldsymbol{\widetilde{r}})=\nabla\times\overline{\overline{\boldsymbol{G}}}_{e}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ and $\displaystyle\overline{\overline{\boldsymbol{G}}}_{M}(\boldsymbol{r};\boldsymbol{\widetilde{r}})=\nabla\times\overline{\overline{\boldsymbol{G}}}_{m}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$, where $\displaystyle\nabla$ operates on functions involving $\displaystyle\boldsymbol{r}$ alone. More details about the DGF are given in Appendix B.1.

$\displaystyle\boldsymbol{J^{e}}(\boldsymbol{r})$ and $\displaystyle\boldsymbol{J^{m}}(\boldsymbol{r})$ are the Fourier transforms of the electric and magnetic current densities, respectively. The spectral densities of the components of $\displaystyle\boldsymbol{J^{e}}(\boldsymbol{r})$ and $\displaystyle\boldsymbol{J^{m}}(\boldsymbol{r})$ are related by the fluctuation-dissipation theorem of the second kind: Rytov (1967); Eckhardt (1982)

[TABLE]

where $\displaystyle p,q=1,2,3$ are the Cartesian components of the current densities, $\displaystyle\Im{\left(\cdot\right)}$ denotes the imaginary component, $\displaystyle\hbar$ is the reduced Planck’s constant, $\displaystyle k_{b}$ is Boltzmann’s constant, $\displaystyle T$ is the thermodynamic temperature, and $\displaystyle\Theta=(\hbar\omega/2)\coth{\left(\hbar\omega/2k_{b}T\right)}$ is the average energy of a harmonic oscillator of frequency $\displaystyle\omega$ at temperature $\displaystyle T$.

Using Eqs. (1)-(4c), the linearized spectral conductance between two objects is given by Narayanaswamy and Zheng (2013a)

[TABLE]

where $\displaystyle X=\hbar\omega/2k_{b}T$, and $\displaystyle T_{A\rightarrow B}^{e}(\omega)$ is the transmissivity function for energy transfer. The total linearized conductance is obtained from

[TABLE]

where $\displaystyle c$ is the speed of light in vacuum and $\displaystyle\lambda=2\pi c/\omega$ is the free-space wavelength. The definition of total conductance in the second integral in Eq. (6) is the result of a change of integration variables from frequency to wavelength. In the discussion of the numerical results (Sec. V), the wavelength-dependent definition of spectral conductance will be employed.

IV Determination of transmissivity function

Since emission and absorption of electromagnetic waves are volumetric phenomena, it might seem intuitive to expect the expression for the transmissivity function for energy transfer to contain two volume integrals. Our previous workNarayanaswamy and Zheng (2013a) has shown that if the two volumes are isothermal, the properties of the vector Helmholtz equation [shown later as Eq. (30) in Appendix B] allow us to reduce the volume integrals into two surface integrals, so that the transmissivity function is given by

[TABLE]

where $\displaystyle\Re{\left(\cdot\right)}$ denotes the real part, $\displaystyle\text{Tr}\left(\cdot\right)$ denotes the trace, $\displaystyle\left(\cdot\right)^{T}$ is the transpose, and $\displaystyle\boldsymbol{r}$ and $\displaystyle\boldsymbol{\widetilde{r}}$ are locations on the surfaces $\displaystyle S_{B}$ and $\displaystyle S_{A}$ of spheres $\displaystyle B$ and $\displaystyle A$, respectively.

Although $\displaystyle\boldsymbol{r}$ and $\displaystyle\boldsymbol{\widetilde{r}}$ will be integrated over the surfaces of spheres $\displaystyle B$ and $\displaystyle A$, respectively, both position vectors must be defined within a medium. This gives us the option of defining the position vectors as either approaching the surfaces of the spheres from the inside or the outside of the spheres. Because the fluctuating charges responsible for emission and absorption are contained within the spheres, it would be natural to assume $\displaystyle\boldsymbol{r}$ and $\displaystyle\boldsymbol{\widetilde{r}}$ approach the spheres’ surfaces from the inside. This is the approach favored in Ref. Narayanaswamy and Chen, 2008 and is referred to as the “interior formula” by Narayanaswamy and Zheng.Narayanaswamy and Zheng (2013a) In the approach we adopt here, $\displaystyle\boldsymbol{r}$ and $\displaystyle\boldsymbol{\widetilde{r}}$ instead approach the surfaces of their respective spheres from region $\displaystyle C$, which is referred to as the “exterior method.”Narayanaswamy and Zheng (2013a) The advantage of this approach is that the same formalism can be used for uncoated spheres or spheres with any number of coatings.

Since $\displaystyle\boldsymbol{r}$ is integrated over the surface of sphere $\displaystyle B$ and $\displaystyle\boldsymbol{\widetilde{r}}$ is integrated over the surface of sphere $\displaystyle A$, it is most convenient to evaluate Eq. 7 using DGFs with $\displaystyle\boldsymbol{r}$ appearing as $\displaystyle\boldsymbol{r_{B}}$ and $\displaystyle\boldsymbol{\widetilde{r}}$ appearing as $\displaystyle\boldsymbol{\widetilde{r}_{A}}$. The appropriate DGFs when $\displaystyle\boldsymbol{\widetilde{r}},\boldsymbol{r}\in C$ are given by (see Appendix B.2 for a full discussion on their determination):

[TABLE]

where we define $\displaystyle\widetilde{m}=\max{\left\{\left|m\right|,1\right\}}$ for compactness, $\displaystyle k=(\omega/c)\sqrt{\varepsilon\mu}$ is the magnitude of the wavevector, $\displaystyle\boldsymbol{M}_{lm}^{(p)}(k\boldsymbol{r})$ and $\displaystyle\boldsymbol{N}_{lm}^{(p)}(k\boldsymbol{r})$ ($\displaystyle p=1$ or $\displaystyle 3$) are vector spherical waves (VSWs) of order $\displaystyle(l,m)$, and $\displaystyle R_{\nu}^{(M)}(r)$ and $\displaystyle R_{\nu}^{(N)}(r)$ are the Mie reflection coefficients at $\displaystyle r=\left|\boldsymbol{r}\right|$ for $\displaystyle\boldsymbol{M}_{\nu m}^{(1)}(k\boldsymbol{r})$ and $\displaystyle\boldsymbol{N}_{\nu m}^{(1)}(k\boldsymbol{r})$ waves, respectively (see Appendix C for definitions).

The value of $\displaystyle p$ determines the behavior of the VSWs in the radial direction. For $\displaystyle p=1$, $\displaystyle\boldsymbol{M}_{lm}^{(p)}(k\boldsymbol{r})$ is an incoming wave scaled radially by $\displaystyle z_{l}^{(1)}(kr)$, the spherical Bessel function of the first kind, where $\displaystyle r=|\boldsymbol{r}|$. For $\displaystyle p=3$, $\displaystyle\boldsymbol{M}_{lm}^{(p)}(k\boldsymbol{r})$ is an outgoing wave scaled radially by $\displaystyle z_{l}^{(3)}(kr)$, the spherical Hankel function of the first kind. The definition of $\displaystyle\boldsymbol{N}_{lm}^{(p)}(k\boldsymbol{r})$ and the relations between VSWs and vector spherical harmonics are given in Appendix A. $\displaystyle\boldsymbol{M}_{lm}^{(p)}(k\boldsymbol{r})$ and $\displaystyle\boldsymbol{N}_{lm}^{(p)}(k\boldsymbol{r})$ are related to one another by $\displaystyle k\boldsymbol{M}_{lm}^{(p)}(k\boldsymbol{r})=\boldsymbol{\nabla}\times\boldsymbol{N}_{lm}^{(p)}(k\boldsymbol{r})$ and $\displaystyle k\boldsymbol{N}_{lm}^{(p)}(k\boldsymbol{r})=\boldsymbol{\nabla}\times\boldsymbol{M}_{lm}^{(p)}(k\boldsymbol{r})$. For any two VSWs $\displaystyle\boldsymbol{P}$ and $\displaystyle\boldsymbol{Q}$, the term $\displaystyle\boldsymbol{P}\boldsymbol{Q}$, such as those appearing in Eqs. (8i)-(8aj), denotes the dyadic product of the two vectors. Lai et al.

$\displaystyle C^{M}$, $\displaystyle C^{N}$, $\displaystyle D^{M}$, $\displaystyle D^{N}$, $\displaystyle F^{M}$, $\displaystyle F^{N}$, $\displaystyle J^{M}$, and $\displaystyle J^{N}$ are unknown coefficients multiplying the VSWs. They are related to one another through a set of coupled linear equations, given by

[TABLE]

where $\displaystyle A_{nm}^{\nu m}$ and $\displaystyle B_{nm}^{\nu m}$ are well-known translation coefficients Chew (1992, 1995); Kim (2004); Dufva et al. (2008) and $\displaystyle d_{AB}=-d_{BA}=d$. The linear system is identical to that given by Mackowski (1991) (see Eqs. (11) and (12) of Ref. [Mackowski, 1991]) but for one small difference: the scattered field coefficients used in this work, such as $\displaystyle C_{n,m}^{lM}$ and $\displaystyle C_{n,m}^{lN}$, are analogous to $\displaystyle-a_{mn}^{i}/\alpha_{n}^{i}$ and $\displaystyle-b_{nm}^{i}/\beta_{n}^{i}$ in Mackowski’s work.

The linear system is obtained by writing the DGFs as expansions of their VSW eigenfunctions in the coordinate systems of both spheres and then reexpanding some of the VSWs such that they may all be written into a consistent coordinate system. This allows for boundary conditions on the DGFs to be enforced. The details of this procedure may be found in Appendix B. It is important to note that the principles of this procedure are well established for dealing with pairs or even clusters of spheres. Previous works have used this technique in not only a wide range of electromagnetic scattering problems,Mackowski (1994); Xu and Gustafson (1996) but also in other wave scattering problems, such as those found in acoustics.Gumerov and Duraiswami (2002); Létourneau et al. (2017) Though the linear system of equations in Eqs. (9c)-(9p) is analogous to Waterman’s T-matrix representation, Mishchenko et al. (1996); Mishchenko and Martin (2013) additional difficulties are introduced due to our interest in near-field phenomena. The criterion for convergence when modeling near-field phenomena is more stringent than for far-field phenomena alone.Sasihithlu and Narayanaswamy (2011a)

For the DGFs given by Eqs. 8i-8aj, the transmissivity function can be evaluated as

[TABLE]

where

[TABLE]

and $\displaystyle P$ may be $\displaystyle M$ or $\displaystyle N$ and $\displaystyle w=a$ or $\displaystyle b$. Equation (13) obeys the reciprocity principle. See Appendix D for proof.

It is worth noting that our choices of symbol and definition of $\displaystyle\epsilon_{\nu}^{(P)}(w)$ were not arbitrary. The spectral emissivity of an isolated sphere of radius $\displaystyle w$ is given byKattawar and Eisner (1970)

[TABLE]

V Numerical results

In this section, all numerical results will be shown for $\displaystyle T=300$ K. In order to compute the transmissivity function efficiently, it must be rewritten into a more computationally efficient form. See Appendix E for details. Additionally, simulated values of spectral and total conductance will be normalized by the conductance between the same two spheres, assuming them to be blackbodies. The spectral conductance between two blackbodies is given by

[TABLE]

and the total conductance is given by

[TABLE]

where $\displaystyle A_{A}$ is the surface area of object A, $\displaystyle F_{A\rightarrow B}$ is the radiative view factor from object A to object B, and $\displaystyle\sigma=\pi^{2}k_{b}^{4}/60c^{2}\hbar^{3}$ is the Stefan-Boltzmann constant.

V.1 Dielectric coating atop metal core

Planar stratified HMMs have previously been investigated for heat transfer applications due to their broadband super-Planckian thermal emission properties.Guo et al. (2012); Guo and Jacob (2013) Though use of non-planar layered media is relatively rare in the study of near-field heat transfer, a thermal microelectromechanical systems (MEMS) device with a layer of polar material atop a curved chromium sensor has been used in extreme near-field experiments by Kim et al..Kim et al. (2015) It is important to note that a single layer of material does not make the device an HMM. Regardless, our work may still give insight into the behavior of the device. Despite their device having coatings, Kim et al. modeled their curved probe as homogeneous, composed of the polar material only. The authors provided only a post hoc justification of this assumption: the seeming agreement between modeled and measured results.

Numerical investigation of heat transfer can shed light on the validity of such an assumption. The simplest test case is to simulate the heat transfer between two identical single-coated spheres. For the simulations here we use a metallic core and a dielectric coating, composed of silverYang et al. (2015) and silica,Palik (1985) respectively. Varying the spheres’ dimensions, the position of the metal/dielectric interface and the separation gap allows for characterization of the impact of dielectric coatings atop metallic cores.

Figure 2 shows the effect of altering the position of the metal/dielectric interface on the spectrum of radiative heat transfer. The coated spheres have an outer radius, coating thickness, and core radius of $\displaystyle R$, $\displaystyle t$, and $\displaystyle R-t$, respectively. Their geometries are fixed such that $\displaystyle R=10$ $\displaystyle\mu$m and the minimum separation gap is $\displaystyle D=1$ $\displaystyle\mu$m.

The spectral conductance is shown in Fig. 2(A). For the case of a homogeneous silica sphere ($\displaystyle t/R=1$), the result is a relatively wideband distribution of frequencies contributing to the radiative transfer which exactly reproduces the results of previous work. Narayanaswamy and Chen (2008) The two well-known surface phonon polariton (SPhP) peaks are present at 8.75 $\displaystyle\mu$m and 20.3 $\displaystyle\mu$m [labeled in Fig. 2(A)]. As the silver core is allowed to grow, the spectrum incrementally changes into the case of two bare silver spheres. Spheres with $\displaystyle t/R=0.1$, $\displaystyle 0.05$, and $\displaystyle 0.01$ exhibit spectral conductances which appear to be roughly scaled versions of each other, the scaling proportional to the thickness of the coating.

The sequence in which the spectrum of radiative transfer for silica spheres transitions to that of silver spheres is not uniform across the spectrum. Although the magnitude of the spectral conductance of silver is always lower than that of silica, increasing the proportion of silver to silica may actually increase the spectral conductance for some wavelengths at some intermediate coating thicknesses. This is evident in the spectrum of spheres with $\displaystyle t/R=0.2$ at 12.5 $\displaystyle\mu$m and 23.5 $\displaystyle\mu$m and $\displaystyle t/R=0.1$ at the wavelength of 10 $\displaystyle\mu$m, where a broad super-Planckian peak, not associated with a SPhP, manifests. At these wavelengths, the conductance of the coated spheres exceeds that of pure silica spheres.

When a thin layer of polar material (actually, the class of materials is broader and any material which has narrow absorption bands may serve as such a thin-film material) is coated on a metallic substrate, the wavelength at which the magnitude of the dielectric function (or equivalently the complex refractive index) of the polar material reaches a minimum, $\displaystyle\lambda_{ENZ}$ (ENZ denoting epsilon near zero), takes special significance.Narayanaswamy et al. (2014) At this wavelength alone, the interface between the coating and vacuum behaves as a highly reflective mirror. The interface between the metallic substrate and the thin film is highly reflective at all wavelengths considered here because of the high dielectric function of metals for mid-infrared wavelengths.

Near $\displaystyle\lambda_{ENZ}$, electromagnetic waves experience reflective conditions at both interfaces, leading to a larger number of reflections than at other wavelengths, if the thin film is not too absorptive. The result of a greater number of reflections is the appearance of an optically thicker film. Because amorphous silica has a relatively high damping, these interesting effects manifest themselves in the near field only when the thickness becomes very small. Amorphous silica has a $\displaystyle\lambda_{ENZ}$ point at 7.95 $\displaystyle\mu$m [see Fig. 2(B)]. Hence, the stand-alone peak in Fig. 2(A) for $\displaystyle t/R=0.01$ at 8.06 $\displaystyle\mu$m is an epsilon near zero mode. As the thickness is increased, this peak can no longer be resolved because of its proximity to a SPhP peak.

Another class of peaks which appears in the spectrum of thermal radiative transfer of coated structures is the Fabry-Perot–like resonance. This type of resonance results from the interference of the multiple reflections of waves within a thin film. Because Fabry-Perot-like resonances require the constructive interference of waves, the location of the peak will drift as the thickness of the coating changes. This type of peak is evident in Fig. 2(A) for $\displaystyle t/R=0.1$ at 10.0 $\displaystyle\mu$m and $\displaystyle t/R=0.05$ at 9.55 $\displaystyle\mu$m.

The cumulative spectral contribution to conductance is shown in Fig. 2(C). The cumulative contribution at wavelength $\displaystyle\lambda$ is given by

[TABLE]

where $\displaystyle\lambda^{\prime}$ is a dummy integration variable. The slopes of the cumulative contribution curves indicate how relatively dominant a wavelength is in contributing to the total conductance. A greater slope indicates a greater relative contribution to the total conductance and vice versa.

The curves for spheres with $\displaystyle t/R\geq 0.2$ have relatively small slopes across most wavelengths. This is consistent with the fairly wideband behavior demonstrated in Fig. 2(A). As $\displaystyle t/R$ decreases, however, wavelengths differentiate into two categories: those with nearly zero slope and those with a very steep slope. For $\displaystyle t/R=0.10$ and $\displaystyle t/R=0.05$, the curves become nearly vertical at the SPhP wavelengths. In the most extreme case, for $\displaystyle t/R=0.01$, the curve is nearly vertical at 7.95 $\displaystyle\mu$m and 19.79 $\displaystyle\mu$m and nearly horizontal elsewhere. These wavelengths correspond to ENZ points of the silica layer. If the minimum separation gap between the spheres were to be decreased, SPhP peaks would grow and eventually dominate over ENZ peaks. The dominance of SPhP modes in the extreme near field is made clear in the discussion of Fig. 3.

As we have shown, adding a very thin layer of a material supporting surface polaritons to a metallic substrate creates a selective near-field emitter. (This was already known to be true in the far fieldGranqvist and Hjortsberg (1980); Narayanaswamy et al. (2014).) Experimental measurement of spheres with very thin coatings would allow for probing of resonant heat transfer, some of which may be due to SPhPs, while suppressing heat transfer at other wavelengths. Because SPhPs are known to dominate heat transfer between polar materials in the extreme near field, isolating the contributions from SPhPs by using a coated sphere would serve as a superior experimental method compared to measuring the effects of SPhPs with homogeneous spheres as in past experiments.Narayanaswamy et al. (2008); Shen et al. (2009); Guha et al. (2012)

Figure 3 shows the effect of varying the separation gap between spheres with outer radii of 5 $\displaystyle\mu$m on total conductance. According to classical radiative transfer, the distance dependence in the far field is due to changes in view factor. Indeed, for gaps such that $\displaystyle D/R\geq 2$, all cases are well approximated as graybodies, as indicated by the curves’ near-zero slopes. In that regime, the total conductance of spheres of constant radius increases as the fraction of silica increases.

As the separation gap decreases, the conductance between spheres with a silica coating begins to be dominated by the surface phonon polaritonic contributions. At a separation gap such that $\displaystyle D/R=0.004$, a coating of just 50 nm of silica can achieve 70% of the conductance of a fully silica sphere. This allows for the creation of spheres with silicalike behavior in the near field but tunable radiative transfer behavior in the far-field. As a simple rule of thumb, the conductance between two silica coated silver spheres exceeds 70% of that between two homogeneous silica spheres for $\displaystyle D/t\lesssim 1/4$. For larger gaps, the conductance is more like that of silver.

This observation partially validates the assumption made by Kim et al..Kim et al. (2015) Their device had a 100-nm silica coating atop an optically opaque chromium thermocouple. Their measurements were performed in the extreme near-field, at gaps ranging from 1 nm to 50 nm. For the smallest gaps, we have shown that the SPhP contributions will be dominant and modeling the whole body as homogeneous silica is a reasonable approximation. However, should the gaps of interest be larger or the materials not be dominated by SPhPs, more care must be taken to properly approximate near-field thermal radiative transfer of coated bodies.

V.2 Dielectric coating atop dielectric core

We have observed that two spheres with dielectric coatings and metallic cores have their total conductance effectively capped at that of two homogeneous spheres with the same dielectric material. A natural question is whether or not two coated spheres can ever exceed the total conductance of two homogeneous spheres which are composed of any of the coated spheres’ constitutive materials. Because of the importance of SPhPs in near-field radiative heat transfer, we simulate the conductance between two coated spheres whose cores and coatings both support SPhPs. As shown in Fig. 4(A), we simulate identical spheres with berylliaPalik (1985) cores and aluminaPalik (1985) coatings, or vice versa, with outer radii of 5 $\displaystyle\mu$m and a minimum separation gap of 100 nm. Those spheres simulated with an alumina coating and a beryllia core such that $\displaystyle 0.01\leq t/R\leq 0.5$ all exceed the total conductance between two homogeneous alumina spheres (which themselves exceed that of two homogeneous beryllia spheres). The maximum occurs at $\displaystyle t/R\approx 0.05$. Although spheres with beryllia coatings and alumina cores never exceed the total conductance of homogeneous alumina spheres, they too exhibit a slight local maximum at the same value of $\displaystyle t/R$. The maximum total conductance of the coated spheres outperforms homogeneous alumina spheres by 8.5%.

When looking at the spectral conductance of the homogeneous spheres and the coated sphere with the maximum conductance, it becomes apparent how the coated spheres are able to outperform the homogeneous spheres. As shown in Fig. 4(B), the coated spheres exhibit spectral features similar to features found in the spectra of their components. Most importantly, the coated spheres strongly reproduce the SPhP peaks of homogeneous alumina at 12.2 $\displaystyle\mu$m and 20.8 $\displaystyle\mu$m while capturing a portion of the enhancement due to the SPhP peak of homogeneous beryllia at 10.0 $\displaystyle\mu$m (SPhP peaks labeled in Fig. 4). This suggests that it may be possible to “stack” the effect of SPhPs at multiple wavelengths by choosing coatings of materials with spectrally spread SPhP peaks.

VI Conclusions

In this paper, we have presented a formula for calculating the near-field thermal radiative transfer between two spheres, which allows for the inclusion of any number of coatings. Since we lack the formalism to analyze the radiative transfer between two homogeneous spheres with anisotropic properties, we are unable, at this stage, to replace the two coated spheres with equivalent hyperbolic metamaterial spheres. Instead, the effective properties of the coated spheres are characterized by the Mie reflection coefficients.

We have demonstrated that spheres with metallic cores and coatings of polar materials also exhibit super-Planckian peaks at wavelengths not corresponding to surface phonon polaritons. Such spheres behave like polar materials in the extreme near field but like metals in the far field. This creates the possibility of measuring near-field radiative transfer due to surface phonon polaritons while suppressing contributions from other modes of heat transfer. We have also shown that a coating of alumina atop a beryllia core can outperform homogeneous alumina or beryllia for an optimal coating thickness, demonstrating that the surface phonon polariton contribution from material inside the sphere may also contribute to the conductance between two spheres.

This work will also be useful in the development of methods to estimate the heat transfer between curved surfaces such as proximity approximationsRousseau et al. (2009); Sasihithlu and Narayanaswamy (2011b); Song et al. (2015) or the thermal discrete dipole approximation.Edalatpour and Francoeur (2014, 2014) Future work should be devoted to developing a formalism of near-field radiative transfer between spheres which permits anisotropic dielectric properties.

VII Acknowledgments

This work was funded partially by ONR Grant N00014-12-1-0996 and NSF IGERT DGE-1069240.

Appendix A Mathematical definitions and useful relations

The vector spherical waves of interest are given by

[TABLE]

where $\displaystyle z_{l}^{(p)}(kr)$ is the spherical Bessel ($\displaystyle p=1$) or spherical Hankel ($\displaystyle p=3$) function of the first kind and $\displaystyle\zeta_{l}^{(p)}(kr)=\frac{1}{kr}\frac{d}{d(kr)}(krz_{l}^{(p)}(kr))$. $\displaystyle\boldsymbol{V}_{lm}^{(1)}(\theta,\phi)$, $\displaystyle\boldsymbol{V}_{lm}^{(2)}(\theta,\phi)$, and $\displaystyle\boldsymbol{V}_{lm}^{(3)}(\theta,\phi)$ are the vector spherical harmonics of order $\displaystyle(l,m)$ for polar and azimuthal angles $\displaystyle\theta$ and $\displaystyle\phi$, respectively. They are defined as

[TABLE]

where $\displaystyle Y_{lm}(\theta,\phi)$ is the scalar spherical harmonic of order $\displaystyle(l,m)$ and $\displaystyle\boldsymbol{\widehat{r}}$, $\displaystyle\boldsymbol{\widehat{\theta}}$, and $\displaystyle\boldsymbol{\widehat{\phi}}$ are the unit vectors of the spherical coordinate system. Scalar spherical harmonics are given by

[TABLE]

where $\displaystyle P_{l}^{m}(x)$ is the associated Legendre polynomial. Olver et al.

The Wronskian of $\displaystyle z_{n}^{(1)}(x)$ and $\displaystyle z_{n}^{(3)}(x)$ is given by

[TABLE]

The following integral identities are useful in evaluating surface integrals on the sphere:

[TABLE]

Appendix B Determination of the dyadic Green’s functions

B.1 Form of the dyadic Green’s functions

A dyadic Green’s function gives the vectorial response at a location due to a vector source at another location, the two positions being the arguments of the DGF. A convenient method to compute the DGFs in Eq. (7) is to expand them in terms of the eigenfunction solutions to the vector Helmholtz equation, given by

[TABLE]

where $\displaystyle\boldsymbol{F}(\boldsymbol{r})$ is the electric or magnetic field at location $\displaystyle\boldsymbol{r}$ and $\displaystyle k=(\omega/c)\sqrt{\varepsilon\mu}$. In spherical coordinates, the eigenfunctions of the vector Helmholtz equation are the vector spherical waves $\displaystyle\boldsymbol{M}_{lm}^{(p)}(k\boldsymbol{r})$ and $\displaystyle\boldsymbol{N}_{lm}^{(p)}(k\boldsymbol{r})$.

Each region has a DGF composed of two parts: a homogeneous DGF, $\displaystyle\overline{\overline{\boldsymbol{G}}}_{0}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$, corresponding to waves which travel directly from $\displaystyle\boldsymbol{\widetilde{r}}$ to $\displaystyle\boldsymbol{r}$, and a scattered DGF corresponding to waves which have experienced scattering at inhomogeneities. When expanding a DGF into its VSW eigenfunctions, the choice of coordinate system for $\displaystyle\boldsymbol{r}$ and $\displaystyle\boldsymbol{\widetilde{r}}$ becomes important because they appear in the arguments of vector spherical waves. Assuming that $\displaystyle\boldsymbol{r},\boldsymbol{\widetilde{r}}\in C$, $\displaystyle\overline{\overline{\boldsymbol{G}}}_{0}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ can be written in the $\displaystyle A$-coordinate system as

[TABLE]

where we define $\displaystyle\widetilde{m}=\max{\left\{\left|m\right|,1\right\}}$. With respect to the center of the $\displaystyle A$-coordinate system, the homogeneous DGF for $\displaystyle r_{A}>\widetilde{r}_{A}$ ($\displaystyle r_{A}<\widetilde{r}_{A}$) corresponds to outgoing (incoming) VSWs at $\displaystyle\boldsymbol{r_{A}}$. The double surface integral in Eq. (7) for computing $\displaystyle T_{A\rightarrow B}^{e}$ requires $\displaystyle\boldsymbol{r_{A}}\in S_{B}$ and $\displaystyle\boldsymbol{\widetilde{r}}_{A}\in S_{A}$. Hence, we must choose the branch of $\displaystyle\overline{\overline{\boldsymbol{G}}}_{0}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ for which $\displaystyle r_{A}>\widetilde{r}_{A}$.

The scattered DGF captures the collective effect of all scattering events at interfaces. The scattered DGF splits naturally into two parts: a part representing waves scattered off of a single sphere only, $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$, and a part representing waves scattered off of both spheres, $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$. $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ obviously includes multiple reflections between the two spheres. Because we chose to write Eq. (33) in the $\displaystyle A$-coordinate system, we must also express $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$, representing waves scattered by sphere A only, in the $\displaystyle A$-coordinate system. That part of the scattered DGF is related to the branch of $\displaystyle\overline{\overline{\boldsymbol{G}}}_{0}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ for $\displaystyle r_{A}<\widetilde{r}_{A}$. In that case, $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ can be thought of as arising from VSWs emitted at $\displaystyle\boldsymbol{\widetilde{r}}$ which travel inward before reflecting off of sphere A and proceeding to $\displaystyle\boldsymbol{r}$. It is given by

[TABLE]

where $\displaystyle R_{l}^{(M)}(a)$ and $\displaystyle R_{l}^{(N)}(a)$ are the Mie reflection coefficients at the surface of sphere A for $\displaystyle\boldsymbol{M}_{lm}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})$ and $\displaystyle\boldsymbol{N}_{lm}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})$ waves, respectively.

Some waves may reflect off of both spheres multiple times on their journey from $\displaystyle\boldsymbol{\widetilde{r}}$ to $\displaystyle\boldsymbol{r}$. The DGF which takes into account those multiple scatterings is given by

[TABLE]

where the coefficients on the VSWs are unknowns to be determined from the boundary conditions, which is discussed shortly.

All DGFs of the form discussed in this work are composed of dyadic productsLai et al. of VSWs. For any dyadic product, the VSW to the right can be any of the VSWs to the right in $\displaystyle\overline{\overline{\boldsymbol{G}}}_{0}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$, i.e., $\displaystyle\boldsymbol{M}_{l,-m}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})$, $\displaystyle\boldsymbol{N}_{l,-m}^{(1)}(k_{C}\boldsymbol{\widetilde{r}_{A}})$, $\displaystyle\boldsymbol{M}_{l,-m}^{(3)}(k_{C}\boldsymbol{\widetilde{r}_{A}})$, or $\displaystyle\boldsymbol{N}_{l,-m}^{(3)}(k_{C}\boldsymbol{\widetilde{r}_{A}})$. The vector to the left has to be an outgoing VSW in either of the coordinate systems in order to satisfy the far-field boundary conditions. Hence, the vector to the left can be a linear combination of the following VSWs: $\displaystyle\boldsymbol{M}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{A}})$, $\displaystyle\boldsymbol{N}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{A}})$, $\displaystyle\boldsymbol{M}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})$, and $\displaystyle\boldsymbol{N}_{\nu m}^{(3)}(k_{C}\boldsymbol{r_{B}})$. The expression in Eq. (45) takes into account all these possibilities.

The magnetic DGF takes the same form as the electric DGF but with a corresponding set of unknown magnetic coefficients, denoted with a tilde. That is to say that each $\displaystyle C$, $\displaystyle D$, $\displaystyle F$, and $\displaystyle J$ has a corresponding magnetic counterpart: $\displaystyle\widetilde{C}$, $\displaystyle\widetilde{D}$, $\displaystyle\widetilde{F}$, and $\displaystyle\widetilde{J}$. The same holds true for the Mie reflection coefficients — every $\displaystyle R$ has a counterpart $\displaystyle\widetilde{R}$. As will become apparent, only the waves in the DGFs representing scattering off of both sphere $\displaystyle A$ and sphere $\displaystyle B$ will contribute to the heat transfer between the two spheres (see Appendix B.2 for more details). This is completely analogous to the result for heat transfer between two semi-infinite half spaces.Narayanaswamy and Zheng (2013a)

In order to simplify the DGFs and determine the unknown coefficients, two steps must be taken. First, the scattered fields must be converted into a single coordinate system. Second, boundary conditions on the DGFs must be enforced at the interfaces between different media. This will allow us to express the boundary conditions in terms of the coefficients multiplying the VSWs.

The vector addition translation theorem is used to convert the coordinate system of VSWs.Chew (1992, 1995); Kim (2004); Dufva et al. (2008) To convert an outgoing VSW in the $\displaystyle B$-coordinate system into VSWs in the $\displaystyle A$-coordinate system, the following expression can be used:

[TABLE]

when the translation $\displaystyle\boldsymbol{r_{B}}-\boldsymbol{r_{A}}=\hat{z}d_{BA}$ is restricted to the $\displaystyle z$ axis and $\displaystyle\left|\boldsymbol{r_{A}}\right|<|d_{BA}|=d$.

Similarly, for conversion of outgoing VSWs from the $\displaystyle A$-coordinate system to VSWs in the $\displaystyle B$-coordinate system, we use

[TABLE]

where $\displaystyle\boldsymbol{r_{A}}-\boldsymbol{r_{B}}=\hat{z}d_{AB}$ and $\displaystyle\left|\boldsymbol{r_{B}}\right|<|d_{AB}|=d$.

At any location $\displaystyle\boldsymbol{r}$ on a boundary and for any $\displaystyle\boldsymbol{\widetilde{r}}$, the DGFs must satisfy

[TABLE]

Utilizing Eqs. (33)-(61) leads to a set of linear equations between the coefficients multiplying the various VSWs that appear in Eq. (45). Some details have been omitted here. For example, the DGFs used to derive the linear equations not only include those when $\displaystyle\boldsymbol{r}$ is in region $\displaystyle C$ [Eqs. (33)-(45)], but also those when $\displaystyle\boldsymbol{r}$ is inside spheres $\displaystyle A$ and $\displaystyle B$ (not given in this paper). The resultant set of linear equations is given in Eqs. (9c)-(9p) where

[TABLE]

The coefficients for the magnetic DGFs, $\displaystyle\widetilde{C}$, $\displaystyle\widetilde{D}$, $\displaystyle\widetilde{F}$, and $\displaystyle\widetilde{J}$ [see discussion following Eq. (45)], are related to the coefficients of the electric DGFs by interchanging $\displaystyle\varepsilon$ and $\displaystyle\mu$ and noticing $\displaystyle\widetilde{R}_{\nu}^{(M)}(r)=R_{\nu}^{(N)}(r)$ and $\displaystyle\widetilde{R}_{\nu}^{(N)}(r)=R_{\nu}^{(M)}(r)$. Accordingly, we get

[TABLE]

B.2 Simplification of dyadic Green’s function

The derivation of the simplified electric DGF, previously presented as Eq. (8i) in its final form, is given below. The electric DGF is given by

[TABLE]

where $\displaystyle\overline{\overline{\boldsymbol{G}}}_{0}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$, $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$, and $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ are given by Eqs. (33)-(45), and represent scattering off of neither sphere, only sphere $\displaystyle A$, and spheres $\displaystyle A$ and $\displaystyle B$, respectively. For reasons which will become apparent shortly, we will group $\displaystyle\overline{\overline{\boldsymbol{G}}}_{0}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ and $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ together. $\displaystyle\overline{\overline{\boldsymbol{G}}}_{0}(\boldsymbol{r};\boldsymbol{\widetilde{r}})+\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ is the DGF for an isolated sphere $\displaystyle A$ in the absence of sphere $\displaystyle B$.

Using Eqs. (33) and (36), we get

[TABLE]

Using Eqs. (62)-(69), Eq. (45) simplifies to

[TABLE]

In order to conveniently evaluate $\displaystyle T_{A\rightarrow B}^{e}$ [Eq. (7)], it will prove useful to convert the VSWs on the left in the dyadic products of the DGF to the $\displaystyle B$-coordinate system. To do so, we employ the vector addition translation theorem, given in Eqs. (54) and (57). Accordingly, we get

[TABLE]

Equations (9g)-(9h) and (9o)-(9p) may then be used to eliminate summation over the index $\displaystyle\nu$ of the $\displaystyle C$ and $\displaystyle F$ coefficients in Eq. 123. After further simplification, we get

[TABLE]

Adding Eq. (85) and Eq. (133), we obtain the Eq. (8i) for $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$. The rationale for splitting $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ into $\displaystyle\overline{\overline{\boldsymbol{G}}}_{0}(\boldsymbol{r};\boldsymbol{\widetilde{r}})+\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ and $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ should be clearer now: only the coefficients of VSWs appearing in $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}^{(sc)\prime\prime}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ manifest in the final expression for $\displaystyle T_{A\rightarrow B}^{e}$ [Eq. (7)].

The form of the magnetic DGF, $\displaystyle\overline{\overline{\boldsymbol{G}}}_{m}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$, is identical, but for $\displaystyle\widetilde{D},\widetilde{J},\widetilde{R}$ in place of $\displaystyle D,J,R$ in Eq. (8i). See Eqs. (70)-(77) for relations between coefficients of VSWs in $\displaystyle\overline{\overline{\boldsymbol{G}}}_{m}$ and $\displaystyle\overline{\overline{\boldsymbol{G}}}_{e}$. $\displaystyle\overline{\overline{\boldsymbol{G}}}_{E}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ and $\displaystyle\overline{\overline{\boldsymbol{G}}}_{M}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ are obtained from their definitions: $\displaystyle\overline{\overline{\boldsymbol{G}}}_{E}(\boldsymbol{r};\boldsymbol{\widetilde{r}})=\nabla\times\overline{\overline{\boldsymbol{G}}}_{e}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$ and $\displaystyle\overline{\overline{\boldsymbol{G}}}_{M}(\boldsymbol{r};\boldsymbol{\widetilde{r}})=\nabla\times\overline{\overline{\boldsymbol{G}}}_{m}(\boldsymbol{r};\boldsymbol{\widetilde{r}})$

Appendix C Computation of effective Mie reflection coefficients

The expressions for Mie reflection coefficients, $\displaystyle R_{\nu}^{(M)}(r_{\rho,s})$ and $\displaystyle R_{\nu}^{(N)}(r_{\rho,s})$ (see Sec. II for the geometry as well as definition of terms), are well known from the literature for both uncoated and single-coated spheres. Kaiser and Schweiger (1993); Bohren and Huffman (2004); Zhou and Hu (2006); Zheng and Ghanekar (2015) The full, multilayered Mie reflection coefficients can be determined recursively in a manner similar to that of Fresnel reflection coefficients for planar stratified media.Narayanaswamy and Zheng (2013b) The recurrence relation is given by

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In Eqs. (136)-(141), $\displaystyle Z=Z_{0}\sqrt{\mu/\varepsilon}$ is the electromagnetic impedance for a dielectric material and $\displaystyle Z_{0}=\sqrt{\mu_{0}/\varepsilon_{0}}$ is the impedance of free space. Additionally, we define $\displaystyle k_{N+1,A}=k_{M+1,B}=k_{C}$. The recursion relation is terminated by $\displaystyle R_{\nu}^{(M)}(r_{0,s})=-u_{\nu}(r_{0,s})$ and $\displaystyle R_{\nu}^{(N)}(r_{0,s})=-v_{\nu}(r_{0,s})$.

Appendix D Reciprocity of transmissivity function

The transmissivity function obeys the reciprocity relation $\displaystyle T_{A\rightarrow B}^{e}\left(\omega\right)=T_{B\rightarrow A}^{e}\left(\omega\right)$. Proof of this property for Eq. (13) is a useful check of the validity of our derived transmissivity function. Because of the properties of electric and magnetic fields, DGFs must obey the following reciprocity relations: Tai1994

[TABLE]

To use these relations, the locations of the $\displaystyle\boldsymbol{r}$ and $\displaystyle\boldsymbol{\widetilde{r}}$ must be interchanged while holding the locations of the spheres fixed. The DGFs from the right-hand sides of Eqs. (142)-(144) can be written with their own set of unknown coefficients, which we accent with a “ ˘ ” symbol. They arise when considering heat transfer from sphere $\displaystyle B$ to $\displaystyle A$. From Eqs. (142)-(144),

[TABLE]

Using Eqs. (145)-(148), showing $\displaystyle T_{A\rightarrow B}^{e}\left(\omega\right)=T_{B\rightarrow A}^{e}\left(\omega\right)$ follows naturally from Eq. (13). We get

[TABLE]

Here, the $\displaystyle l$ and $\displaystyle\nu$ indices and $\displaystyle m$ and $\displaystyle-m$ swap roles, but the overall sum is equal.

Appendix E Computational implementation

To be of any practical use, the transmissivity function must be numerically calculable. A number of numerical problems are introduced, however, due to the presence of spherical Bessel and Hankel functions of high order, which may experience issues with underflow and overflow on computers.Sasihithlu and Narayanaswamy (2014) To avoid such problems, we make two changes in the calculation of the transmissivity function.

First, instead of using the set of coefficients multiplying the VSWs in Eqs. (9c)-(9p), we use the related coefficients defined in Eqs. 62-69 denoted with a $\displaystyle(1)$ superscript. These coefficients are analogous to the unknown coefficients used in Ref. [Narayanaswamy and Chen, 2008] and the scattering coefficients used in Ref. Mackowski and Mishchenko, 2008. Second, we introduce prefactors to the translation and VSW coefficients which help stabilize the linear system used to solve for the VSW coefficients. The details of this procedure are discussed in greater detail by Sasihithlu and Narayanaswamy.Sasihithlu and Narayanaswamy (2014) The prefactors take the form of a ratio of spherical Bessel or Hankel functions. For example, the prefactor for $\displaystyle D_{\nu,m}^{lM^{(1)}}$ is $\displaystyle z_{l}^{(1)}(k_{C}a)/z_{\nu}^{(1)}(k_{C}b)$. The linear system is modified in such a way that we solve directly for the value of $\displaystyle\left[z_{l}^{(1)}(k_{C}a)/z_{\nu}^{(1)}(k_{C}b)\right]D_{\nu,m}^{lM^{(1)}}$, and other such coefficients. The expression for the transmissivity function may then be modified such that the VSW coefficients, with the appropriate prefactors, appear explicitly.

After some manipulation, we get

[TABLE]

where

[TABLE]

and $\displaystyle w=a$ or $\displaystyle b$.

The above results are obtained using the relation

[TABLE]

and simplified using Wronskian relations [see Eq. (25)].