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_BT (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 $ℛ$ is diagonal in the frequency ω, parallel wavevector k_β= (k_x, k_y), and polarization 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)+τ²(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:

 

Fig. 1 (a) Schematic of two inhomogeneous ε(z) slabs (A and B) separated by a vacuum gap of size d along the z direction. The radiative heat transfer (RHT) rate between them depends on their local temperatures T_A,B and reflection coefficients ${ℛ}_A{ }_{,B}$. Associated with each slab is a coordinate system centered at the slab–vacuum interface and pointing away from the gap. (b) Enhancement factor comparing RHT between optimized inhomogeneous (solid lines) or optimized uniform (dashed lines) slabs against that of semi-infinite uniform plates of Re[ε] = −1, at a fixed vacuum wavelength λ = 8 μm, as a function of material loss Im[ε], and for two representative separations d = 10 nm (red) and 500 nm (blue). The black dotted line shows the theoretical bound described by Eq. (3). (c) Transmission coefficient $\mathcal{Z}$ corresponding to uniform semi-infinite slabs of ε = 1 (black), finite slabs of optimal thickness t_opt and permittivity ε_opt (blue), and inhomogeneous slabs resulting from optimization (red), as a function of the dimensionless wavevector k_β d at fixed d = 10 nm and Im[ε] = 10⁻². The inset shows a typical dielectric profile Re[ε(z)] needed to achieve $\mathcal{Z}\approx 1$ over a broad range of k_β. To enhance the readability, two different x-scales are used in the ranges [0, 20] nm and [20, 100] nm, and the points resulting from numerical optimization are connected with segments.

Download Full Size | PDF

We seek dielectric profiles ε(z) that maximize the heat flux $H\left[ℛ\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_i), with all {ε_i}_i=1,…,_N 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 ≳ 100. Furthermore, while gradient information $\frac{\partial H}{\delta {\epsilon }_i}=\frac{\partial H}{\partial ℛ}\frac{\partial ℛ}{\partial {\epsilon }_i}$ is typically needed for large N [59], $\frac{\partial ℛ}{\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[ε_i, 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.

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_BT)−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 ${ℛ}_{\mathrm{A}}{ℛ}_{\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_{\max }\approx \frac{1}{2d}\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}\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⁻² and d = 10 nm at the vacuum wavelength λ = 8 μm (frequency ω = 2πc/λ ≈ 2.35 · 10¹⁴ 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₀) 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

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.

 

Fig. 2 (a) Schematic of an inhomogeneous ε(z) slab and a dipole separated by a vacuum gap of size d along the z direction. (b) Dielectric profile Re[ε(z)] corresponding to inhomogeneous slabs optimized to increase RHT from a dipole a distance d = 1 μm away from their z = 0 interface, at a fixed vacuum wavelength λ = 8 μm (frequency ω ≈ 0.785c/d). The profiles are obtained under different constraints on the maximum possible permittivity ε_max ≡ max|Re ε| = {5, 40} uniform (upper and lower figures) but correspond to the same Im[ε] = 10⁻³. (c) Imaginary part of the reflection coefficient $\mathrm{Im}[ℛ(k_{\beta })]$ at frequency ω as a function of k_β, for optimized inhomogeneous slabs [red and green solid lines, corresponding to the ε profiles in panel (b)], optimized uniform slabs (blue solid line), along with the theoretical bound (black dotted line) described by Eq. (8). (d) LDOS $ℒ(k_{\beta })$ (in SI units) at the location of the dipole and at frequency ω, with the same color convention of panel (c). The inset in (d) shows the k_β-integrated spectrum $ℒ({\omega }^{\prime })$ near ω as a function of the dimensionless frequency (ω′ − ω)/ω Im[ε], indicating that contributions from smaller k_β are increasingly sensitive to the wavelength.

Download Full Size | PDF

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:

First, as in the plate–plate scenario, we investigate the possible enhancements in Φ(ω) or $ℒ\left(\omega \right)$ that can arise from a spatially varying dielectric profile. Unlike the previous scenario, where $\mathrm{\Phi }\left[ℛ\right]$ was a highly non-monotonic function of $ℛ$, here the RHT integrand is linearly proportional to $\mathrm{Im}\left[ℛ\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[ℛ\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 $ℒ\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 $ℛ$ 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)² Im ε and $ℒ\left(\omega \right)$. Consequently, the integrated response $ℒ\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 $ℒ\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 $ℒ\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⁻³. 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[ℛ\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[ℛ\right]$ and $ℒ$ 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 $ℒ({\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 ${ℒ}_0\left(\omega \right)=\frac{{\omega }^2}{{\pi }^2{c}^3}$ [70], as a function of d and for multiple values of Im[ε] = {10⁻³, 10⁻², 10⁻¹} (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 $ℒ\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.

 

Fig. 3 (a) Peak LDOS of optimized inhomogeneous slabs (solid lines) and optimized uniform slabs (dashed lines), normalized by the vacuum LDOS ${ℒ}_0\left(\omega \right)=\frac{{\omega }^2}{{\pi }^2{c}^3}$, as a function of the dimensionless separation ωd/c, for multiple material loss rates Im[ε] = {10⁻¹, 10⁻², 10⁻³} (black, blue, and red, respectively). The green dotted line marks the largest possible LDOS in the far field, given by Eq. (10). (b) Enhancement factor comparing the peak LDOS of optimized inhomogeneous and uniform slabs.

Download Full Size | PDF

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

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 ≳ 3c/ω. In order to explain this feature, we observe that for Im[ε] = {10⁻³, 10⁻²} (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 ε₁ and ε₂ small imaginary part Im[ε] ≪ 1. The hyperbolic regions of the spectrum are those in which ε₁ε₂ < 1, i.e. when the real parts have opposite signs. More specifically, one typically defines two categories of HMMs: type I with ε₁ > 0 and ε₂ < 0, and type II with ε₁ < 0 and ε₂ > 0. While in the extreme near-field there is no distinction between the two types, given that only the product ε₁ε₂ matters [71], the two exhibit very different behavior in the mid-field. Focusing on the dominant, p polarization $\left(ℛ\approx r^p\right)$, the relevant quantity in type-I HMMs is

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].