## Abstract

Radiative heat transfer between uniform plates is bounded by the narrow range and limited contribution of surface waves. Using a combination of analytical calculations and numerical gradient-based optimization, we show that such a limitation can be overcome in complicated multilayer geometries, allowing the scattering and coupling rates of slab resonances to be altered over a broad range of evanescent wavevectors. We conclude that while the radiative flux between two inhomogeneous slabs can only be weakly enhanced, the flux between a dipolar particle and an inhomogeneous slab—proportional to the local density of states—can be orders of magnitude larger, albeit at the expense of increased frequency selectivity. A brief discussion of hyperbolic metamaterials shows that they provide far less enhancement than optimized inhomogeneous slabs.

© 2017 Optical Society of America

## 1. Introduction

Radiative heat transfer (RHT) between two bodies separated by a vaccuum gap and held at different temperatures is limited by Stefan–Boltzmann’s law in the far field, i.e. for gap distances *d* much larger than the thermal wavelength *ħc/k _{B}T* (several microns at ambient temperature). This limitation no longer applies in the near field, where the physics is dramatically altered by the presence of evanenscent tunneling of photons. For instance, the coupling of bound surface modes (e.g. phonon-polaritons in dielectrics) can result in orders-of-magnitude larger flux rates [1–5]. Frequency-selective bounds to RHT between arbitrarily shaped homogeneous bodies were recently established [6], demonstrating that typically considered configurations (e.g. two parallel slabs of ordinary materials) tend to be highly suboptimal. The possibility of increasing both far and near-field heat exchange through geometry has attracted a remarkable amount of attention over the last several deades due to its relevance in many technological applications, from thermophotovoltaic (TPV) energy conversion [1] to near-field thermal lithography [7].

In the far field, where the traditional, ray-optical form of Kirchoff’s law (equating emissivity and absorptivity) is applicable [8], computational advances have made it possible to exploit a variety of optimization strategies (exploiting a variety of methods, e.g. random-walk, genetic and particle swarm algorithms, and the Taguchi method) to realize, for instance, selective and/or wide-angle absorbers whose emissivity can come close to the blackbody limits [9–20], and whose objective is usually that of increasing the performance of a TPV device or solar cell. More recently, development of adjoint optimization techniques in combination with fast numerical EM solvers have allowed application of large-scale optimization [21] methods capable of efficiently tackling problems involving much higher number of degrees of freedom. For instance, these inverse-design techniques have been exploited to enhance the efficiency of solar-cell absorbers [12, 22], tailor the spectrum of incandescent sources [23], and to increase the functionality of photonic-crystal absorbers [24] and TPV systems [25]. Far less explored is the near field, where the possibility of tuning and amplifying RHT has been investigated mainly through parametric studies, i.e. shape optimization, involving only a few degrees of freedoms. For example, several authors have studied the role of thin films in amplifying the heat flux between planar objects [26–29]. Others have focused on enhancing desirable optical properties by examining variations with respect to Drude or Drude–Lorentz model parameters [28, 30–32], e.g. through doping [31, 33, 34]. Various geometries such as dielectric [35–37] and metallic [38–45] gratings, and even finite bodies [46–48] have been recently considered. Less restrictive inverse-design techniques have been recently employed to improve the performances of heat-assisted magnetic recording (HAMR) head [49]. Relatedly, some authors have also addressed the modulation and optimization of a closely related quantity, the near-field electromagnetic local density of states (LDOS). For instance, Ref. [50] performed a parametric study of the LDOS close to a multilayer arrangement of silicon carbide and silicon thin slabs as a function of distance and number of layers while Ref. [51] employed genetic algorithms to optimize the LDOS in proximity of a multilayer binary structure composed of alluminium and lossless dielectric layers.

In this paper, we apply adjoint, large-scale RHT optimization with the goal of enhancing RHT in two commonly studied scenarios, involving either a slab or a dipolar particle in proximity to a slab. In particular, we relax the typical assumption of homogeneous materials and consider instead arbitrary dielectric profiles (along the slab symmetry axis) in planar geometries. In order to fully address the richness of a non-uniform permittivity, we employ gradient-based optimization over a large number of degrees of freedom (number of layers ≳ 1000), a regime where previously explored techniques based on global, derivative-free optimization are bound to fail. We demonstrate that appropriately optimized, multilayer structures can lead to larger RHT compared to the best possible homogeneous thin films. Our results are motivated by and extend previously derived frequency-selective bounds for homogeneous, planar media [6], quantifying the maximum flux rates that can be achieved at any given frequency through the careful interference or “rate matching” of scattered and absorbed surface waves; such a condition can only be satisfied at a single wavevector in uniform slabs but can be much more broadband in inhomogeneous media. Note that the bound described in Ref. [6] is only an asymptotic bound, in the sense that it is derived under small separations. We find that with respect to uniform slabs, RHT between inhomogeneous slabs can only be weakly enhanced, with the optimized structures approaching the achievable bounds.

Much larger enhancements are possible in the dipole–plate geometry, where in principle RHT for a sufficiently small dipole—proportional to the local density of states (LDOS)—can be infinite [6, 52, 53]. We provide theoretical bounds for the flux contribution of each wavevector in semi-infinite, uniform media and discover structures that can approach these bounds over a broad range of wavevectors, limited mainly by the difficult task of achieving perfect and broadband absorption of waves close to the light line. Specifically, we find that the optimization procedure is able to produce finite but in principle arbitrarily large enhancements, limited only by the finite numerical discretization (number of layers) and sharp dielectric variations of the structures: these can reach two orders of magnitude at mid-range *ωd/c* ≳ 1 separations but at the expense of frequency selectivity. In the deep near field (*ωd/c* ≪ 1), on the other hand, the condition leading to maximal LDOS is satisfied by the ideal resonant plasmonic condition (Re[*ε*] = −1) in a uniform medium and is therefore bounded by previously derived bounds on homogeneous media [6]. We remark that while the problem of enhancing RHT between planar materials appears highly constrained and difficult to tackle, there is much more room for improvements and tunability when it comes to tailoring the LDOS in the vicinity of a multilayer body, which could also be of interest in other contexts, such as in near-field imaging [54,55].

Finally, although we only consider inhomogeneous slabs composed of dielectric media, we argue that similar enhancements are expected in the case of inhomogeneous, magnetically anisotropic materials. This is exemplified, for instance, in the case of hyperbolic metamaterials (HMM), which were recently considered in the study of radiative heat transfer in multilayer geometries [56]. In particular, focusing on the case of uniaxial media, we show that regardless of separation, HMMs suffer from the same limitations of uniform thin films and thus do not yield significant flux enhancements.

## 2. Formulation

In what follows, we consider RHT in geometries involving slabs of arbitrarily varying dielectric profiles *ε*(*z*) along the direction *z* perpendicular to the slab–vacuum interfaces, depicted in Fig. 1(a). RHT in such a setup can be described via the fluctuational electrodynamics framework developed by Rytov, Polder, and van Hove (see Refs. [1,3] and references therein). Specifically, we extend a recently developed formulation of this problem [57, 58] that expresses the flux in scenarios involving two and three uniform bodies as a function of their reflection and transmission matrices. This approach, together with a semi-analytic expression for the reflection matrix of a slab of arbitrary *ε*(*z*), enables gradient-based optimizations of RHT. Since the system is time and translationally invariant in *x*–*y*, the reflection matrix
$\mathcal{R}$ is diagonal in the frequency *ω*, parallel wavevector **k*** _{β}*= (

*k*,

_{x}*k*), and polarization

_{y}*p*, and can be cast as the solution of a differential equation, derived as follows. Consider for each slab a local coordinate system such that

*z*= 0 lies at the interface between each slab and the vacuum gap (of size

*d*) and points away from the interface. Given a slab occupying the region [

*z, t*] (where 0

*< z < t*and

*t*is the possibly infinite thickness of the slab), let

*R*(

*z*) be the coefficient describing the reflection on the left side, i.e. at the interface

*z*. Adding a film of infinitesimal thickness Δ

*z*at

*z*, the reflection coefficient of the combined system (at the

*z*− Δ

*z*interface) is given by

*R*(

*z*−Δ

*z*) =

*ρ*(

*z*)+

*τ*

^{2}(

*z*)

*R*(

*z*)/(1−

*R*(

*z*)

*ρ*(

*z*)), where

*ρ*and

*τ*are the reflection and transmission coefficients of the film, respectively. Taking the limit Δ

*z*→ 0, one obtains the following nonlinear differential equation:

*R*(

*t*) = 0, describing the absence of the slab (thus a vanishing reflection coefficient) for

*z*=

*t*, completely specifies the reflectivity of the system. Here

*r*is the ordinary Fresnel reflection coefficient, ${k}_{zm}\left(z\right)=\sqrt{\epsilon \left(z\right){\left(\omega /c\right)}^{2}-{k}_{\beta}^{2}}$ the perpendicular wavevector inside the slab. Note that in the limiting case of uniform

*ε*, Eq. (1) yields the well-known solution $R\left(z\right)=r\left[1-{\mathrm{e}}^{2i{k}_{zm}\left(t-z\right)}\right]/\left[1-{r}^{2}{\mathrm{e}}^{2i{k}_{zm}\left(t-z\right)}\right]$, going to

*r*in the case of a semi-infinite slab (

*t*→ +∞). Equation 1 can be directly solved to obtain the reflection coefficient of a slab of arbitrarily varying

*ε*(

*z*): in the case of the two slabs of Fig. 1(a), one would also need to specify the boundary conditions

*R*(+∞) = 0 for both slabs A and B. Once the function

*R*(

*z*) is known for each slab, $\mathcal{R}=R\left(0\right)$ represents the reflection coefficient needed to calculate RHT and analyze the possible enhancements arising from a given

*ε*(

*z*), investigated below through optimization techniques.

We seek dielectric profiles *ε*(*z*) that maximize the heat flux
$H\left[\mathcal{R}\right]$ at a given frequency. In practice, given a choice of slab thickness, numerical evaluations require that the slab be discretized into segments, forming a multilayer geometry. We thus replace the function *ε*(*z*) with a piecewise-constant function that assumes values *ε _{i}* =

*ε*(

*z*), with all {

_{i}*ε*}

_{i}

_{i}_{=1},…,

*taken as variable degrees of freedom. Note that the size of individual layers typically needs to be very small in order to resolve the exponential decay of evanescent fields, with typical*

_{N}*N*≳ 100. Furthermore, while gradient information $\frac{\partial H}{\delta {\epsilon}_{i}}=\frac{\partial H}{\partial \mathcal{R}}\frac{\partial \mathcal{R}}{\partial {\epsilon}_{i}}$ is typically needed for large

*N*[59], $\frac{\partial \mathcal{R}}{\partial {\epsilon}_{i}}$ can be straightforwardly obtained from Eq. (1). Because RHT can diverge with vanishing loss rate [6], while most plasmonic material have nonzero loss rate (unless doped with gain media [60]), we consider finite but uniform Im[

*ε*] thoughout the slabs, focusing only on optimizing with respect to {Re[

*ε*, which leaves modifications in the scattering rather than loss rate as the primary source of enhancement. Since this objective function is far from convex [61], we exploit local optimization algorithms [62, 63], with the optimum depending on the initial conditions [61] as well as the optimization procedures, albeit with no guarantee of finding global optimum.

_{i}## 3. Plate–plate

We first consider the scenario of two inhomogeneous, parallel slabs, depicted in Fig. 1(a), with both slabs A and B assumed to be in local thermal equilibrium at temperatures *T*_{A} and *T*_{B}, respectively. In this case, the well-known formalism for two slabs of uniform permittivities can be employed. The total heat transfer per unit surface is given by [27]
$H={\displaystyle \int \frac{\mathrm{d}\omega}{2\pi}\left[\mathrm{\Theta}\left({T}_{\mathrm{A}}\right)-\mathrm{\Theta}\left({T}_{\mathrm{B}}\right)\right]\mathrm{\Phi}\left(\omega \right)}$, where the monochromatic spectral component
$\mathrm{\Phi}\left(\omega \right)={\displaystyle {\sum}_{s,p}{\displaystyle \int \frac{\mathrm{d}{k}_{\beta}}{2\pi}{k}_{\beta}{\mathcal{Z}}_{s\left(p\right)}\left({k}_{\beta},\omega \right)}}$, with Θ(*T*) = *ħω/*[exp(*ħω/k _{B}T*)−1] denoting the Planck function. Here, we focus on the

*p*polarization, which supports surface modes and hence dominates RHT at short separations

*ωd/c*≲ 1. $\mathcal{Z}$ is known as the heat transmission coefficient, whose evanescent component is given by:

The extension of Eq. (2) to the case of inhomogeneous slabs consists in replacing the reflection operators with those from Eq. (1).

It is well known that energy transfer is optimal when the scattering and absorption decay rates of the surface modes described by Eq. (2) are equal [64]. This corresponds to a maximum transmissivity
$\mathcal{Z}=1$, realized at the wavenumber for which
${\mathcal{R}}_{\mathrm{A}}{\mathcal{R}}_{\mathrm{B}}^{*}={\mathrm{e}}^{2\mathrm{Im}\left[{k}_{z}\right]d}$. For two uniform semi-infinite slabs, such a “rate-matching” condition can only occur at a single *k _{β}*, depending on the separation and loss rate [6], in which case
$\mathcal{Z}$ exhibits a typical Lorentzian lineshape as a function of

*k*, whose peak lies close to a typical cutoff wavevector

_{β}*k*

_{max}, above which $\mathcal{Z}$ is exponentially suppressed. For small

*d*and loss rates and assuming operation close to the surface plasmon resonance of a uniform slab, i.e. Re[−1

*/χ*

_{spp}] = 1/2, such a cutoff can be approximated by [6] ${k}_{\mathrm{max}}\approx \frac{1}{2d}\mathrm{ln}[\frac{{|{\chi}_{\text{spp}}|}^{4}}{4{(\mathrm{Im}{\chi}_{\text{spp}})}^{2}}]$, where

*χ*=

*ε*− 1 is the susceptibility of the material. This leads to an upper bound on the RHT between uniform semi-infinite slabs, given by ${\mathrm{\Phi}}_{0}=\frac{1}{2\pi {d}^{2}}\mathrm{ln}[\frac{{|{\chi}_{\text{spp}}|}^{4}}{4{(\mathrm{Im}{\chi}_{\text{spp}})}^{2}}]$ [6].

Relaxing the assumption of uniform *ε* allows modes of different *k _{β}* to experience different scattering and absorption rates, potentially allowing rate-matching to not only persist over all

*k*≲

_{β}*k*

_{max}but even for

*k*

_{β}> k_{max}. However, in practice it is unfeasible to extend the optimization to the regime

*k*

_{β}> k_{max}. As a matter of fact, already in the case of two uniform slabs of finite thickness, the coefficient $\mathcal{Z}$ can in principle approach 1 at arbitrarily large

*k*, but only at the expense of exponentially diverging $\mathrm{Re}\left[\epsilon \right]=-\mathrm{Im}\left[\epsilon \right]{\mathrm{e}}^{{k}_{\beta}d}$ and vanishing thickness $t=\frac{2}{\mathrm{Im}\left[\epsilon \right]{k}_{\beta}}{\mathrm{e}}^{-{k}_{\beta}d}$ and bandwidths $\mathrm{\Delta}{k}_{\beta}={k}_{\beta}{\mathrm{e}}^{-{k}_{\beta}d}$, making such an intereference effect highly impractical if at all feasible to sustain over a wide range of

_{β}*k*. We find, however, that there exist structures that can achieve rate matching over all

_{β}*k*≲

_{β}*k*

_{max}and therefore whose flux is described by a larger achievable bound ${\tilde{\mathrm{\Phi}}}_{0}$, obtained by integrating Eq. (2) with $\mathcal{Z}=1$ up to

*k*

_{max}, given by:

The ratio
$\frac{{\tilde{\mathrm{\Phi}}}_{0}}{{\mathrm{\Phi}}_{0}}$ depends only on material loss, increasing with decreasing loss, as shown in Fig. 1(b) (black dotted line). In practice, however, such an enhancement tends to be relatively small because of the logarithmic power law and the fact that inhomogeneity seems to barely increase *k*_{max}, which is instead primarily determined by the choice of loss rate and *d*.

Next, we exploit optimization to discover inhomogeneous structures that can achieve or approach the monochromatic, achievable bounds
${\tilde{\mathrm{\Phi}}}_{0}$ above [Eq. (3)]. Although we consider the permittivities of the two slabs to be independent degrees of freedom, we find that the optimization always leads to a symmetric dielectric profile, which might be the easiest approach to guarantee the rate matching condition over a wide range of *k _{β}*, while asymmetric profiles are not, in principle, prohibited. The inset of Fig. 1(c) shows

*ε*(

*z*) for one such optimized slab, corresponding to the particular choice of Im[

*ε*] = 10

^{−2}and

*d*= 10 nm at the vacuum wavelength

*λ*= 8

*μ*m (frequency

*ω*= 2

*πc/λ*≈ 2.35 · 10

^{14}rad/s), with each dot representing the permittivity of a 1 nm-thick layer. The function Re[

*ε*(

*z*)] shows a strong variation near the slab–vacuum interface, approaching −1 away from the interface. This somewhat un-intuitive dielectric profile leads to nearly perfect $\mathcal{Z}=1$ for all

*k*

_{β}< k_{max}≈ 5

*/d*, shown in Fig. 1(c) (red solid line). In contrast, the transmissivity of either a uniform, semi-infinite slab of

*ε*= −1 (black solid line) or a finite slab of optimal thickness

*t*

_{opt}and permittivity

*ε*

_{opt}(blue solid line) exhibit $\mathcal{Z}~1$ over a smaller range of

*k*. Note that the transmissivity of the finite slab exhibits two peaks, which stem from a splitting of the modes localized around the surfaces of both slabs. Moreover, we find that these enhancements are robust with respect to frequency and layer thicknesses on the order of 10 nm. Figure 1(b) shows the enhancement factor associated with two different structures, optimized to maximize RHT at either

_{β}*d*= 10 nm (red lines) or

*d*= 500 nm (blue lines), as a function of the loss rate. We find that at small

*d*= 10 nm, the achievable enhancements agree well with the predictions of Eq. (3) while at larger

*d*= 500 nm and smaller loss rates, larger enhancements are observed; such a discrepancy is expected since the non-retarded approximation (assuming only contributions from

*k*≫

_{β}*k*

_{0}) employed in deriving Eq. (3) underestimates

*k*

_{max}at mid-range separations

*ωd/c*≳ 1. Even then, the flux rates of inhomogeneous slabs (solid lines) tends to be larger than those of uniform slabs (dashed lines).

While the configuration of isotropic (possibly inhomogeneous) parallel slabs explored above yields only a small enhancement factor, stemming primarily from the logarithmic power law and difficulty of increasing *k*_{max}, one might ask whether it is possible to further increase *k*_{max} by exploting more exotic media, e.g. electric and magnetic anisotropy. For the sake of simplicity, we restrict our discussion to uniaxial media, described by diagonal permittivity and permeability tensors given by:

For such media, RHT is still described by Eq. (2), but with a modified expression of perpendicular component of the wavevector inside the medium, which reads
${k}_{zm}=\sqrt{{\epsilon}_{\parallel}{\mu}_{\parallel}{k}_{0}^{2}-\frac{{\epsilon}_{\parallel}}{{\epsilon}_{\perp}}{k}_{\beta}^{2}}$ and
$\sqrt{{\epsilon}_{\parallel}{\mu}_{\parallel}{k}_{0}^{2}-\frac{{\mu}_{\parallel}}{{\mu}_{\perp}}{k}_{\beta}^{2}}$ for the *p* and *s* polarizations, respectively [65]. Moreover, the corresponding Fresnel reflection coefficients have to be modified and become
${r}^{p}=\frac{{\epsilon}_{\parallel}{k}_{z}-{k}_{zm}}{{\epsilon}_{\parallel}{k}_{z}+{k}_{zm}}$ and
${r}^{s}=\frac{{\mu}_{\parallel}{k}_{z}-{k}_{zm}}{{\mu}_{\parallel}{k}_{z}+{k}_{zm}}$ [65]. Since the reflection coefficients of the two polarizations are symmetric with respect to exchange of *ε* and *μ*, implying the existence of both electric or magnetic phonon-polaritons [66], one can focus on only one polarization, e.g. the *p* polarization. In the extreme near-field regime, where the non-retarded approximation is valid, the corresponding reflection coefficient is well approximated by

*k*

_{max}is to reduce Im[

*ε*] at the surface-resonance frequency, defined by the resonance condition Re[

*ε*

_{iso}] = −1, for which we have, assuming the same loss rate, Im[

*ε*

_{║}] = Im[

*ε*

_{⊥}] ≡ Im[

*ε*] ≪ 1,

Thus, one concludes that the anisotropy does not allow one to decrease losses and hence increase the cutoff wavevector *k*_{max}.

## 4. Dipole–slab

In this section, we study RHT between an inhomogeneous slab and a dipole, depicted in Fig. 2(a). Beginning with a brief overview of the formulation, we establish an approximate bound for RHT in the case of semi-infinite, homogeneous bodies and exploit optimization to show that multilayer structures can come close to approaching these limits over a broad range of wavevectors, leading to orders-of-magnitude larger RHT. Our work extends recent studies of near-field RHT between dipoles and HMMs or thin films [67] to the mid-field regime.

Consider a small sphere of radius *R* ≪ *d*, approximated as a dipolar particle of polarizability *α*, *d* being its distance from a planar substrate [see Fig. 2(a)]. In what follows, we focus on the off-resonance regime
$\alpha {\mathcal{D}}_{ll}\ll 1$ [68], in which the spectral transfer rate reads
$\mathrm{\Phi}\left(\omega \right)=4{\displaystyle {\sum}_{l=x,y,z}\mathrm{Im}\left[\alpha \left(\omega \right)\right]\mathrm{Im}\left[{\mathcal{D}}_{ll}\left(\omega \right)\right]}$, where
${\mathcal{D}}_{ll}$ denotes the Green’s function of the slab at the position of the dipole [2]. At short separations *ωd/c* ≲ 1, the relevant tensor components of the Green’s functions are:

*k*where the latter overestimates RHT since it also captures power radiating into the vacuum region [69].

_{β}< ω/cFirst, as in the plate–plate scenario, we investigate the possible enhancements in Φ(*ω*) or
$\mathcal{L}\left(\omega \right)$ that can arise from a spatially varying dielectric profile. Unlike the previous scenario, where
$\mathrm{\Phi}\left[\mathcal{R}\right]$ was a highly non-monotonic function of
$\mathcal{R}$, here the RHT integrand is linearly proportional to
$\mathrm{Im}\left[\mathcal{R}\right]$. Assuming small losses
$\mathrm{Im}\left[\epsilon \right]\ll 1$ and defining *k* = *ck _{β}/ω*, a useful figure of merit is the maximum
$\mathrm{Im}\left[\mathcal{R}\left(k\right)\right]$ and optimal permittivity of a uniform, semi-infinite slab at any given

*k*:

Both quantities are strongly divergent at *k* = 1, suggesting that the monochromatic LDOS
$\mathcal{L}\left(\omega \right)$ of an inhomogeneous slab can in principle be unbounded, with the main contribution to the divergence coming from wavevectors near the light cone *k _{β}* =

*ω/c*. We first observe that such a divergence would theoretically persist for any distance

*d*, since the separation enters Eq. (7) only as a parameter through the exponential factors. However, Eq. (9) shows that at least for a uniform medium, maximizing Im $\mathcal{R}$ in the limit

*k*→ 1 requires a perfect metal (

*ε*→ −∞), which can be shown to screen the response at other

*k*, resulting in vanishing Δ

*k*= 2(

*k*−1)

^{2}Im

*ε*and $\mathcal{L}\left(\omega \right)$. Consequently, the integrated response $\mathcal{L}\left(\omega \right)$ of a uniform slab is finite and maximized by a finite thickness and permittivity, which can be found numerically. More significant improvements, however, can be gained from a spatially varying

*ε*(

*z*), which provides additional degrees of freedom with which to simultaneously tune the scattering rate at different

*k*, allowing the response to approach the bounds described by Eq. (8) over wider bandwidths. Realizing such an enhancement presents, however, both conceptual and numerical challenges: waves approaching the light line have increasingly longer wavelength in the

*z*direction and are thus increasingly sensitive to spatial variations, requiring longer slabs and sharper variations in

*ε*(

*z*). Any numerical optimization strategy will thus benefit only from limited enhancements coming from

*k*≳ 1 due to the finite number of layers needed to resolve

*ε*(

*z*) and finite size of the system. Indeed, in the numerical optimization we find that $\mathcal{L}\left(\omega \right)$ improves slowly with increasing number of layers and resolution.

One should also consider that, as shown above in the simple case of a uniform slab, the permittivity that maximizes the LDOS at *k* ≫ 1 equals −1, while Eq. (9) requires *ε* = −∞ as *k* → 1. In practice, the optimal profile results from a tradeoff between these two conditions, since very high values of *ε* act to screen the response from other regions of the slab. Such distinct and challenging requirements make the optimization procedure highly nontrivial, increasing the computational cost of RHT calculations and slowing the convergence rate of the optimization algorithm, which can get easily trapped in multiple local optima. To illustrate these features, we perform separate optimizations with different constraints on the maximum allowed permittivity *ε*_{max} = max{|Re *ε*|}, which limits potential enhancements coming from waves near the light line.

Figure 2 reports optimizations of the evanescent contribution to
$\mathcal{L}\left(\omega \right)$ at a typical vacuum wavelength *λ* = 8 *μ*m and for *d* = 1 *μ*m. Figure 2(b) shows representative profiles Re[*ε*(*z*)] obtained under different *ε*_{max} = {5, 40} and at fixed Im[*ε*] = 10^{−3}. Noticeably, the lower profile exhibits rapid, subwavelength variations over small (tens to hundreds of nanometers) regions: as discussed earlier, these are needed to maximize
$\mathrm{Im}\left[\mathcal{R}\right]$ near *k _{β}* =

*ω/c*while avoiding screening effects at larger

*k*, with higher |Re[

_{β}*ε*]| occurring away from the interface for the same reason. This explains why larger

*ε*

_{max}lead to greater enhancements, illustrated in Figs. 2(c) and (d), which show $\mathrm{Im}\left[\mathcal{R}\right]$ and $\mathcal{L}$ as a function of

*k*. The results, which are also compared against those of uniform slabs of optimal thickness and permittivity (blue line), reveal that inhomogeneous structures can approach the bounds described by Eq. (8) (dashed black line) over much broader range of

_{β}*k*. Although producing a significant increase with respect to uniform slabs, the optimization fails as

_{β}*k*→

_{β}*ω/c*due to the practical limitations discussed above. Moreover, these enhancements will necessarily come at the expense of increased frequency selectivity, since waves near the light line are most sensitive to deviations in the long-range spatial pattern of the structure, here optimized to realize a specific interference pattern at

*ω*. This feature is apparent from the inset of Fig. 2(d), which shows the spectra $\mathcal{L}({\omega}^{\prime})$ of the optimized uniform and inhomogeneous slabs from above (neglecting material dispersion): namely, the contribution of lower

*k*states becomes increasingly restricted to frequencies

*ω*′ ≈

*ω*as

*k*→

*ω/c*. (Note that the factor of Im

*ε*in the abscissa is there because just as in the case of a uniform medium, the bandwidth of the Lorentzian-like spectrum is proportional to the loss rate.)

We now explore the enhancement factor from optimized inhomogeneous slabs for a wide range of separations *d* ∈ [0.5, 9] *μ*m. To begin with, Fig. 3(a) shows the maximum LDOS of the optimized inhomogeneous (solid lines) and uniform (dashed lines) slabs, normalized by the vacuum LDOS
${\mathcal{L}}_{0}\left(\omega \right)=\frac{{\omega}^{2}}{{\pi}^{2}{c}^{3}}$ [70], as a function of *d* and for multiple values of Im[*ε*] = {10^{−3}, 10^{−2}, 10^{−1}} (red, blue, and black lines respectively), with their ratio, the enhancement factor, depicted in Fig. 3(b). As shown, in both situations the LDOS increases rapidly with decreasing *ωd/c* and decreasing material losses. In particular, the enhancement factor [in principle infinite for any *d*, as suggested by Eq. (8)] increases up to a maximum value (dictated by the smallest participating, enhanced wavevector) and then decreases, approaching 1 as *d* → ∞. Essentially, at large *d*, the evanescent LDOS becomes increasingly dominated by wavevectors close to the light line (for which the optimization procedure fails). Consequently, beyond some separation the propagating contributions to
$\mathcal{L}\left(\omega \right)$ dominate. Such finite enhancements also imply that at small *d*, where the LDOS becomes increasingly dominated by large *kc/ω* ≫ 1 waves, the optimal slab is one satisfying the typical resonant condition of Re[*ε*] = −1 and hence the enhancement factor approaches 1.

For comparison, the green dotted line in Fig. 3(a) denotes the largest achievable far-field LDOS in planar media,

*R*

_{p}_{(}

_{s}_{)}| ≤ 1, in which case $|\mathrm{Re}[{R}_{p\left(s\right)}{\mathrm{e}}^{2i{k}_{z}d}]|\le 1$. More precisely, the limit in Eq. (10) can be derived by observing that for the

*p*polarization,

*s*polarization, leading to Eq. (10). Since it is derived under the assumption of an integrand maximized for any

*k*, $\mathrm{max}\left\{{\mathcal{L}}^{\text{prop}}\right\}$ provides an upper bound that is challenging to realize. However, as shown in Fig. 3, it is still smaller than the evanescent part of the LDOS of optimized inhomogeneous slabs over a wide range of separations. In particular, we remark that at the separations where the enhancement factor peaks, the evanescent contribution to the LDOS of the optimized homogeneous slabs is more than an order of magnutude larger than $\mathrm{max}\left\{{\mathcal{L}}^{\text{prop}}\right\}$.

_{β}We remark that in the scenario of inhomogenous slabs, the LDOS decreases smoothly with increasing separation and material loss, while in the case of optimal uniform slabs, material losses become irrelevant at distances *d* ≳ 3*c/ω*. In order to explain this feature, we observe that for Im[*ε*] = {10^{−3}, 10^{−2}} (red and blue lines respectively), the first-order derivative of the peak LDOS for uniform slabs is a discontinuous function of separation at a given
$\tilde{d}$, depending on Im[*ε*]. This results from an abrupt transition between two mechanisms of enhancement. In particular, depending on the separation, the uniform slab either maximizes the LDOS at some intermediate *k* ≫ 1 through a resonant Re[*ε*] ≈ −1 (small *d*) or near *k* ≈ 1 with Re[*ε*] → −∞ (large *d*). In the latter case, the imaginary part of the permittivity becomes irrelevant.

Since LDOS enhancements from inhomogeneous slabs prove significant at mid-range separations, one may wonder whether similar enhancements can be achieved in previously studied planar geometries. One such geometry are HMMs, which consist of alternating metal and dielectric layers and are known to exhibit hyperbolic dispersion [71, 72]. While it has been demonstrated that RHT between a dipole and a HMM in the deep near field is no larger than that of an appropriately designed uniform thin film [67], we analyze below whether this remains true in the mid-field regime. Consider for instance, a HMM of period *a* ≪ *λ, d* described as an effective, uniform anisotropic medium, with permittivities *ε*_{║} (surface-parallel) and *ε*_{⊥} and respectively, and by assumption the same (surface-perpendicular), having real parts *ε*_{1} and *ε*_{2} small imaginary part Im[*ε*] ≪ 1. The hyperbolic regions of the spectrum are those in which *ε*_{1}*ε*_{2} < 1, i.e. when the real parts have opposite signs. More specifically, one typically defines two categories of HMMs: type I with *ε*_{1} > 0 and *ε*_{2} < 0, and type II with *ε*_{1} < 0 and *ε*_{2} > 0. While in the extreme near-field there is no distinction between the two types, given that only the product *ε*_{1}*ε*_{2} matters [71], the two exhibit very different behavior in the mid-field. Focusing on the dominant, *p* polarization
$\left(\mathcal{R}\approx {r}^{p}\right)$, the relevant quantity in type-I HMMs is

*ε*]), resulting in much smaller RHT in the limit of small Im[

*ε*] ≪ 1. The behavior of type-II HMMs, on the other hand, varies depending on two different

*k*regimes: $k>\sqrt{{\epsilon}_{2}}$, in which case $\mathrm{Im}\left[{\mathcal{R}}_{\text{II}}\right]\approx \mathrm{Im}\left[{\mathcal{R}}_{\mathrm{I}}\right]<1$ and $1\le k<\sqrt{{\epsilon}_{2}}$ (requiring

*ε*

_{2}> 1), in which case the peak

*ε*

_{1}→ −∞ and

*ε*

_{2}→ ∞. However, as before, such a large dielectric constant results in a significant screening effect and thus narrow bandwidth, whose full width at half maximum $\mathrm{\Delta}k\approx \frac{\left|{\epsilon}_{1}-1+3{\epsilon}_{2}\left({\epsilon}_{2}-1\right)\right|}{\sqrt{{\left(1-{\epsilon}_{1}{\epsilon}_{2}\right)}^{3}{\epsilon}_{2}\left(1-{\epsilon}_{1}\right)}}\mathrm{Im}\left[\epsilon \right]$. In the limit of a diverging permittivity, Δ

*k*→ 0 faster than $\mathrm{Im}\left[{\mathcal{R}}_{\parallel}\right]$ diverges, leading to vanishing RHT. Hence, one finds once again that HMMs are in principle not better than uniform, isotropic thin films at mid-field separations and are therefore never “discovered” by the optimization method.

## 5. Concluding remarks

We have studied optimization of frequency-selective near-field radiative heat transfer between either two slabs or between a dipolar particle and a slab, allowing inhomogenous dielectric profiles along the symmetry axis of the slabs. Our approach relied on application of large-scale, gradient-based optimization, which allows dealing with a large number of degrees of freedom, in combination with a simple-to-evaluate differential equation describing the reflectivity of a planar substrate of varying *ε*. In the plate–plate scenario, we extended previous theoretical bounds derived for homogeneous planar structures [6] and showed that inhomogeneous dielectric profiles enable rate matching (perfect absorption) of waves over a wide range of wavevectors, limited only by material loss rates. In spite of nearly perfect coupling, we find relatively low enhancement factors due to the logarithmic dependence of the amplification on material losses. Nevertheless, we observe that the plate–plate optimization is robust up to single-layer thicknesses of the order of tens of nanometers and with respect to frequency, making these predictions in principle experimentally feasible albeit challenging to test, e.g. by employing molecular beam epitaxy with vertically varying doping concentration and hence dielectric properties.

We found, however, that much larger enhancements can be achieved in the dipole–plate scenario, where RHT is proportional to the LDOS. While the LDOS in this configuration is in principle unbounded, we find that any RHT enhancement is in practice limited by a challenging and seemingly prohibitive compromise requiring perfect-metal behavior to increase the flux at small wavevectors (approaching the light line) and resonant, negative permittivities (Re[*ε*] ≈ −1) needed at larger wavevectors. Our optimization procedure was able to discover *ε* profiles able to partially satisfy these two criteria and hence achieve large near-field absorption over a wide range of wavevectors, leading to enhancement factors of up to two orders of magnitude at mid-range separations, limited only by the finite discretization and size of the multilayer structure as well as the existence of multiple local minima. As explained, because the source of the enhancement is increased absorption from waves near the light line, such increased RHT in the dipole–plate configuration necessarily comes at the expenses of increased frequency selectiviy and lack of robustness, depending sensitively on the precise dielectric arrangement and choice of wavelength, posing challenges to experimental realization, which constrains the achievable enhancement factor. Nevertheless, our numerical experiments provide a proof of principle that multilayer structures can overcome the wavevector–bandwidth limitations of uniform slabs in the near field, allowing the wavevector and spectral response of planar materials to be tailored. Moreover, our optimization technique on Green’s function through many degrees of freedom will also be relevant to Raman enhancement [73] and Casimir-Polder interactions [74]. Future work on optimization can be further extended to 2D-gratings through brute-force or with the help of analytic tools such as rigorous coupled wave analysis (RCWA) [45].

## Funding

National Science Foundation (NSF) (DMR-1454836, DMR-1420541, EFRI-1005093); Princeton Center for Complex Materials (MRSEC).

## References and links

**1. **S. Basu, Z. Zhang, and C. Fu, “Review of near-field thermal radiation and its application to energy conversion,” Int. J. Ener. Res. **33**, 1203 (2009). [CrossRef]

**2. **A. Volokitin and B. N. Persson, “Near-field radiative heat transfer and noncontact friction,” Rev. Mod. Phys. **79**, 1291 (2007). [CrossRef]

**3. **K. Joulain, J.-P. Mulet, F. Marquier, R. Carminati, and J.-J. Greffet, “Surface electromagnetic waves thermally excited: Radiative heat transfer, coherence properties and casimir forces revisited in the near field,” Surf. Sci. Rep. **57**, 59 (2005). [CrossRef]

**4. **C. R. Otey, L. Zhu, S. Sandhu, and S. Fan, “Fluctuational electrodynamics calculations of near-field heat transfer in non-planar geometries: A brief overview,” J. Quant. Spectrosc. Radiat. Transf. **132**, 3 (2014). [CrossRef]

**5. **X. Liu, L. Wang, and Z. M. Zhang, “Near-field thermal radiation: Recent progress and outlook,” Nanoscale Microscale Thermophys. Eng. **19**, 98 (2015). [CrossRef]

**6. **O. D. Miller, S. G. Johnson, and A. W. Rodriguez, “Shape-independent limits to near-field radiative heat transfer,” Phys. Rev. Lett. **115**, 204302 (2015). [CrossRef] [PubMed]

**7. **W. Srituravanich, N. Fang, C. Sun, Q. Luo, and X. Zhang, “Plasmonic nanolithography,” Nano Lett. **54**, 1085–1088 (2004). [CrossRef]

**8. **L. Wang, S. Basu, and Z. Zhang, “Direct and indirect methods for calculating thermal emission from layered structures with nonuniform temperatures,” J. Heat Transf. **133**, 072701 (2011). [CrossRef]

**9. **I. Celanovic, F. OSullivan, M. Ilak, J. Kassakian, and D. Perreault, “Design and optimization of one-dimensional photonic crystals for thermophotovoltaic applications,” Opt. Lett. **29**, 863 (2004). [CrossRef] [PubMed]

**10. **J. Drevillon and P. Ben-Abdallah, “Ab initio design of coherent thermal sources,” J. Appl. Phys. **102**, 114305 (2007). [CrossRef]

**11. **A. David, H. Benisty, and C. Weisbuch, “Optimization of light-diffracting photonic-crystals for high extraction efficiency leds,” J. Displ. Technol. **3**, 133 (2007). [CrossRef]

**12. **M. Ghebrebrhan, P. Bermel, Y. Avniel, J. D. Joannopoulos, and S. G. Johnson, “Global optimization of silicon photovoltaic cell front coatings,” Opt. Express **17**, 7505 (2009). [CrossRef] [PubMed]

**13. **N. P. Sergeant, O. Pincon, M. Agrawal, and P. Peumans, “Design of wide-angle solar-selective absorbers using aperiodic metal-dielectric stacks,” Opt. Express **17**, 22800 (2009). [CrossRef]

**14. **Y.-B. Chen and K.-H. Tan, “The profile optimization of periodic nano-structures for wavelength-selective thermophotovoltaic emitters,” Int. J. Heat Mass Transf. **53**, 5542–5551 (2010). [CrossRef]

**15. **S. Chen, H. Cheng, H. Yang, J. Li, X. Duan, C. Gu, and J. Tian, “Polarization insensitive and omnidirectional broadband near perfect planar metamaterial absorber in the near infrared regime,” Appl. Phys. Lett. **99**, 253104 (2011). [CrossRef]

**16. **E. Nefzaoui, J. Drevillon, and K. Joulain, “Selective emitters design and optimization for thermophotovoltaic applications,” J. Appl. Phys. **11**, 084316 (2012). [CrossRef]

**17. **K. X. Wang, Z. Yu, V. Liu, Y. Cui, and S. Fan, “Absorption enhancement in ultrathin crystalline silicon solar cells with antireflection and light-trapping nanocone gratings,” Nano Lett. **12**, 1616–1619 (2012). [CrossRef] [PubMed]

**18. **C. Simovski, S. Maslovski, I. Nefedov, and S. Tretyakov, “Optimization of radiative heat transfer in hyperbolic metamaterials for thermophotovoltaic applications,” Opt. Express **21**, 14988 (2013). [CrossRef] [PubMed]

**19. **P. Wang and R. Menon, “Optimization of generalized dielectric nanostructures for enhanced light trapping in thin-film photovoltaics via boosting the local density of optical states,” Opt. Express **22**, A99 (2013). [CrossRef]

**20. **P. D. Anderson and M. L. Povinelli, “Optimized emission in nanorod arrays through quasi-aperiodic inverse design,” Opt. Lett. **40**, 2672 (2015). [CrossRef] [PubMed]

**21. **J. S. Jensen and O. Sigmund, “Topology optimization for nano-photonics,” Laser Photon. Rev. **5**, 308–321 (2011). [CrossRef]

**22. **V. Ganapati, O. D. Miller, and E. Yablonovitch, “Light trapping textures designed by electromagnetic optimization for subwavelength thick solar cells,” IEEE J. Photovolt. **4**, 175–182 (2014). [CrossRef]

**23. **O. Ilic, P. Bermel, G. Chen, J. D. Joannopoulos, I. Celanovic, and M. Soljačić, “Tailoring high-temperature radiation and the resurrection of the incandescent source,” Nat. Nanotech. **11**, 320 (2016). [CrossRef]

**24. **P. I. Borel, A. Harpoth, L. H. Frandsen, M. Kristensen, P. Shi, J. S. Jensen, and O. Sigmund, “Topology optimization and fabrication of photonic crystal structures,” Opt. Express **12**, 1996–2001 (2004). [CrossRef] [PubMed]

**25. **P. Bermel, W. Ghebrebrhan, M. Chan, Y. X. Yeng, M. Araghchini, R. Hamam, C. H. Marton, K. F. Jensen, M. Soljačić, J. D. Joannopoulos, S. G. Johnson, and I. Celanovic, “Design and global optimization of high-efficiency thermophotovoltaic systems,” Opt. Express **18**, A314 (2010). [CrossRef] [PubMed]

**26. **S.-A. Biehs, “Thermal heat radiation, near-field energy density and near-field radiative heat transfer of coated materials,” Eur. Phys. J. B **58**, 423 (2007). [CrossRef]

**27. **P. Ben-Abdallah, K. Joulain, J. Drevillon, and G. Domingues, “Near-field heat transfer mediated by surface wave hybridization between two films,” J. Appl. Phys. **106**, 044306 (2009). [CrossRef]

**28. **S. Basu and M. Francoeur, “Maximum near-field radiative heat transfer between thin films,” Appl. Phys. Lett. **98**, 243120 (2011). [CrossRef]

**29. **S. V. Boriskina, J. K. Tong, Y. Huang, J. Zhou, V. Chiloyan, and G. Chen, “Enhancement and tunability of near-field radiative heat transfer mediated by surface plasmon polaritons in thin plasmonic films,” in “*Photonics*,”, vol. 2 (Multidisciplinary Digital Publishing Institute, 2015), vol. 2, p. 659. [CrossRef]

**30. **X. J. Wang, S. Basu, and Z. M. Zhang, “Parametric optimization of dielectric functions for maximizing nanoscale radiative transfer,” J. Phys. D: Appl. Phys. **42**, 245403 (2009). [CrossRef]

**31. **Y. Zhao, G. H. Tang, and Z. Y. Li, “Parametric investigation for suppressing near-field thermal radiation between two spherical nanoparticles,” Int. Commun. Heat Mass Transf. **39**, 918 (2012). [CrossRef]

**32. **E. Nefzaoui, Y. Ezzahri, J. Drevillon, and K. Joulain, “Maximal near-field radiative heat transfer between two plates,” Eur. Phys. J. Appl. Phys. **63**, 30902 (2013). [CrossRef]

**33. **E. Rousseau, M. Laroche, and J.-J. Greffet, “Radiative heat transfer at nanoscale mediated by surface plasmons for highly doped silicon,” Appl. Phys. Lett. **95**, 231913 (2009). [CrossRef]

**34. **S. Basu, B. J. Lee, and Z. M. Zhang, “Near-field radiation calculated with an improved dielectric function model for doped silicon,” J. Heat Transf. **132**, 023302 (2010). [CrossRef]

**35. **X. L. Liu and Z. M. Zhang, “Graphene-assisted near-field radiative heat transfer between corrugated polar materials,” Appl. Phys. Lett. **104**, 251911 (2014). [CrossRef]

**36. **C. Fu and Z. Zhang, “Nanoscale radiation heat transfer for silicon at different doping levels,” Int. J. Heat Mass Transf . **49**, 1703 (2006). [CrossRef]

**37. **X. Liu, B. Zhao, and Z. M. Zhang, “Enhanced near-field thermal radiation and reduced casimir stiction between doped-si gratings,” Phys. Rev. A **91**, 062510 (2015). [CrossRef]

**38. **R. Guérout, J. Lussange, F. S. S. Rosa, J.-P. Hugonin, D. A. R. Dalvit, J.-J. Greffet, A. Lambrecht, and S. Reynaud, “Enhanced near-field thermal radiation and reduced casimir stiction between doped-si gratings,” Phys. Rev. B **85**, 180301 (2012). [CrossRef]

**39. **J. Dai, S. A. Dyakov, and M. Yan, “Enhanced near-field radiative heat transfer between corrugated metal plates: Role of spoof surface plasmon polaritons,” Phys. Rev. B **92**, 035419 (2015). [CrossRef]

**40. **J. Dai, S. A. Dyakov, and M. Yan, “Radiative heat transfer between two dielectric-filled metal gratings,” Phys. Rev. B **93**, 155403 (2016). [CrossRef]

**41. **J. Dai, S. A. Dyakov, S. I. Bozhevolnyi, and M. Yan, “Near-field radiative heat transfer between metasurfaces: A full-wave study based on two-dimensional grooved metal plates,” Phys. Rev. B **94**, 125431 (2016). [CrossRef]

**42. **Y. Yang and L. Wang, “Spectrally enhancing near-field radiative transfer between metallic gratings by exciting magnetic polaritons in nanometric vacuum gaps,” Phys. Rev. Lett. **117**, 044301 (2016). [CrossRef] [PubMed]

**43. **R. Messina, A. Noto, B. Guizal, and M. Antezza, “Radiative heat transfer between metallic gratings using adaptive spatial resolution,” Phys. Rev. B **95**, 125404 (2017). [CrossRef]

**44. **H. Chalabi, A. Alù, and M. L. Brongersma, “Focused thermal emission from a nanostructured sic surface,” Phys. Rev. B **94**, 094307 (2016). [CrossRef]

**45. **V. Fernández-Hurtado, F. J. García-Vidal, S. Fan, and J. C. Cuevas, “Enhancing near-field radiative heat transfer with si-based metasurfaces,” Phys. Rev. Lett. **118**, 203901 (2017). [CrossRef] [PubMed]

**46. **O. Ramezan Choubdar and M. Nikbakht, “Radiative heat transfer between nanoparticles: Shape dependence and three-body effect,” J. Appl. Phys. **120**, 144303 (2016). [CrossRef]

**47. **H. Chalabi, E. Hasman, and M. L. Brongersma, “Effect of shape in near-field thermal transfer for periodic structures,” Phys. Rev. B **91**, 174304 (2015). [CrossRef]

**48. **H. Chalabi, E. Hasman, and M. L. Brongersma, “Near-field radiative thermal transfer between a nanostructured periodic material and a planar substrate,” Phys. Rev. B **91**, 014302 (2015). [CrossRef]

**49. **S. Bhargava and E. Yablonovitch, “Lowering hamr near-field transducer temperature via inverse electromagnetic design,” IEEE Trans. Magn. **51**, 3100407 (2015). [CrossRef]

**50. **Y. Li, C. J. Zhang, T.-B. Wang, J.-T. Liu, T.-B. Yu, Q.-H. Liao, and N.-H. Liu, “Modulation of electromagnetic local density of states by coupling of surface phonon-polariton,” Mod. Phys. Lett. B **31**, 1750050 (2017). [CrossRef]

**51. **P. Ben-Abdallah, K. Joulain, J. Drevillon, and G. Domingues, “Tailoring the local density of states of nonradiative field at the surface of nanolayered materials,” Appl. Phys. Lett. **94**, 153117 (2009). [CrossRef]

**52. **O. D. Miller, A. G. Polimeridis, M. H. Reid, C. W. Hsu, B. G. DeLacy, J. D. Joannopoulos, M. Soljačić, and S. G. Johnson, “Fundamental limits to optical response in absorptive systems,” Opt. Express **24**, 3329 (2016). [CrossRef] [PubMed]

**53. **S. I. Maslovski, C. R. Simovski, and S. A. Tretyakov, “Overcoming black body radiation limit in free space: metamaterial superemitter,” New J. Phys. **18**, 013034 (2016). [CrossRef]

**54. **S. A. Ramakrishna, J. Pendry, M. Wiltshire, and W. Stewart, “Imaging the near field,” J. Mod. Opt. **50**, 1419 (2003). [CrossRef]

**55. **X. Zhang and Z. Liu, “Superlenses to overcome the diffraction limit,” Nat. Mater. **7**, 435 (2008). [CrossRef] [PubMed]

**56. **S.-A. Biehs and P. Ben-Abdallah, “Near-field heat transfer between multilayer hyperbolic metamaterials,” Z. Naturforsch. A **72**, 115–127 (2017). [CrossRef]

**57. **R. Messina and M. Antezza, “Scattering-matrix approach to casimir-lifshitz force and heat transfer out of thermal equilibrium between arbitrary bodies,” Phys. Rev. A **84**, 042102 (2011). [CrossRef]

**58. **R. Messina and M. Antezza, “Three-body radiative heat transfer and casimir-lifshitz force out of thermal equilibrium for arbitrary bodies,” Phys. Rev. A **89**, 052104 (2014). [CrossRef]

**59. **R. H. Byrd, J. Nocedal, and R. A. Waltz, “Knitro: An integrated package for nonlinear optimization,” in “*Large-scale nonlinear optimization*,” (Springer, 2006), p. 35. [CrossRef]

**60. **C. Khandekar, W. Jin, O. D. Miller, A. Pick, and A. W. Rodriguez, “Giant frequency-selective near-field energy transfer in active–passive structures,” Phys. Rev. B **94**, 115402 (2016). [CrossRef]

**61. **S. Boyd and L. Vandenberghe, *Convex optimization* (Cambridge university, 2004). [CrossRef]

**62. **K. Svanberg, “A class of globally convergent optimization methods based on conservative convex separable approximations,” SIAM J. Ooptimiz. **12**, 555 (2002). [CrossRef]

**63. **J. Nocedal, “Updating quasi-newton matrices with limited storage,” Math. Comput. **35**, 773 (1980). [CrossRef]

**64. **B. Liu and S. Shen, “Theory of the thermal radiation from optical antennas in far-and near-fields,” https://arXiv:1509.00939 (2015).

**65. **L. Hu and S. Chui, “Characteristics of electromagnetic wave propagation in uniaxially anisotropic left-handed materials,” Phys. Rev. B **66**, 085108 (2002). [CrossRef]

**66. **M. Francoeur, S. Basu, and S. J. Petersen, “Electric and magnetic surface polariton mediated near-field radiative heat transfer between metamaterials made of silicon carbide particles,” Opt. Express **19**, 18774 (2011). [CrossRef] [PubMed]

**67. **O. D. Miller, S. G. Johnson, and A. W. Rodriguez, “Effectiveness of thin films in lieu of hyperbolic metamaterials in the near field,” Phys. Rev. Lett. **112**, 157402 (2014). [CrossRef] [PubMed]

**68. **A. Volokitin and B. Persson, “Dissipative van der waals interaction between a small particle and a metal surface,” Phys. Rev. B **65**, 115419 (2002). [CrossRef]

**69. **A. Oskooi and S. G. Johnson, “Electromagnetic wave source conditions,” https://arXiv:1301.5366 (2013).

**70. **R. Messina, J.-P. Hugonin, J.-J. Greffet, F. Marquier, Y. De Wilde, A. Belarouci, L. Frechette, Y. Cordier, and P. Ben-Abdallah, “Tuning the electromagnetic local density of states in graphene-covered systems via strong coupling with graphene plasmons,” Phys. Rev. B **87**, 085421 (2013). [CrossRef]

**71. **Y. Guo, C. L. Cortes, S. Molesky, and Z. Jacob, “Broadband super-planckian thermal emission from hyperbolic metamaterials,” Appl. Phys. Lett. **101**, 131106 (2012). [CrossRef]

**72. **S.-A. Biehs, M. Tschikin, R. Messina, and P. Ben-Abdallah, “Super-planckian near-field thermal emission with phonon-polaritonic hyperbolic metamaterials,” Appl. Phys. Lett. **102**, 131106 (2013). [CrossRef]

**73. **X. Ling, W. Fang, Y.-H. Lee, P. T. Araujo, X. Zhang, J. F. Rodriguez-Nieva, Y. Lin, J. Zhang, J. Kong, and M. S. Dresselhaus, “Raman enhancement effect on two-dimensional layered materials: graphene, h-bn and mos2,” Nano Lett. **14**, 3033–3040 (2014). [CrossRef] [PubMed]

**74. **A. Seyedzahedi, A. Moradian, and M. Setare, “Intensifying the casimir force between two silicon substrates within three different layers of materials,” Phys. Lett. A **380**, 1475–1480 (2016). [CrossRef]