DOI: 10.1103/PhysRevB. 99.041403 Radiative heat transfer (RHT) between objects separated by near-field distances (on the order of or shorter than the thermal wavelength) exhibits a number of remarkable features. Primarily, evanescent contributions, absent in the far field, cause the rate of heat exchange to scale inversely with separation down to nanometer scales, leading to flux rates many orders of magnitude larger than those predicted by the Stefan-Boltzmann (blackbody) law [1][2][3][4]. Further, this increased flux can be enhanced and controlled by resonant electromagnetic surface modes [5][6][7][8], allowing heat to be concentrated into narrow and designable spectral bandwidths. These properties, in principle, provide a means of significantly improving the degree to which heat can be manipulated compared to thermal conduction, leading to the consideration of applications and devices based on near-field RHT in various contexts, with proposals ranging from thermophotovoltaic energy capture [9][10][11][12][13][14] to high-density heat management [15][16][17], and heatassisted magnetic recording [18,19].
Yet, a concrete understanding of what can be accomplished with near-field RHT remains elusive. Simple high-symmetry structures where analytic solutions are possible provide valuable insight, but appear to be far from ideal. In particular, the most well-studied platform for implementing selective RHT enhancement [20,21], involving parallel metal plates [22][23][24][25] supporting surface resonances (plasmon or phonon polaritons) [26,27], has critical deficiencies. First, as dictated by the Planck distribution, there is a natural wavelength scale for observing significant thermal radiative effects near ambient * arod@princeton.edu temperatures that spans the near-to mid-infrared (1-10 μ m) spectrum [28][29][30]. Typical (low-loss) optical materials do not support polariton resonances at these wavelengths, and often lack sufficient thermal stability to withstand long-term operations [31][32][33]. Second, the tightest known limits of RHT between extended structures, recently derived using energy conservation and reciprocity arguments [34], reveal that both practical material (|χ | 2 / Im χ ) and geometric ( λ/d ) factors enable enhancements relative to blackbody emission as large as
orders of magnitude larger than what is achievable with ideal planar media, suggesting that dramatic improvements are possible through the use of nanostructured surfaces [34]. (In this expression, d stands for separation, and χ material susceptibility, assumed to be the same in both bodies.) Moreover, the |χ | 4 / Im[ χ ] 2 scaling in (1) indicates that materials exhibiting a strong metallic response, far from the typical planar surface-plasmon polariton (SPP) condition Re[ χ ] = -2, have a much greater potential for enhancement. To date, however, this behavior has yet to be observed, and tested geometries employing nonresonant materials [35][36][37][38][39] have yet to surpass the optimal rates associated with planar bodies. In this Rapid Communication, we provide direct evidence of this material scaling enhancement in periodic systems, at least for a wide range of values of χ . Building on our earlier examination of RHT between multilayer bodies [40], we now employ inverse design [41] to investigate RHT between generalized two-dimensional (2D) gratings (restricting the analysis to realistic materials and fabricable structures). Focusing on lossy metals far from the SPP condition at infrared wavelengths, we observe that while F does not approach (1), the predicted material scaling is nevertheless present. For the specific example of tungsten (W), we also find a grating geometry that selectively achieves 50% of the RHT of a pair of ideal, lossless metal plates satisfying the SPP condition, for separations as small as 2% of a tunable design wavelength. These predictions represent RHT enhancements of nearly two orders of magnitude compared to corresponding planar objects, confirming the potential of even relatively simple structuring for selectively enhancing RHT.
The application of inverse design to selective RHT enhancement between extended structures is complicated by several characteristics. First, near-field RHT is controlled by evanescent electromagnetic fluctuations. The large density of these states makes it challenging to apply traditional resonant nanophotonic strategies for enhancing far-field emission over narrow spectral windows [42,43]. Moreover, the characteristically large field amplitudes and subwavelength features of evanescent states make them sensitive to small variations in structural and material properties [44], and correspondingly accurate modeling of RHT requires fine numerical resolution [45][46][47][48]. Second, unlike the far field, RHT cannot be decomposed into approximately equivalent independent subproblems. Alterations in the structure of any one object affects the response of the entire system, meaning that the scattering properties of all bodies must be controlled simultaneously. Finally, Maxwell's equations depend nonlinearly on the dielectric properties and shapes of all bodies, making the optimization nonconvex [49] and any a priori guarantee of globally optimal solutions impossible.
Consequently, tractable general approaches for calculating near-field RHT have only recently been realized [45][46][47]50], and nearly all previously studied geometries have been designed via trial-and-error approaches exploiting brute-force searches over a handful of high-symmetry design parameters [51]. Beginning with bulk metamaterials [39,[52][53][54][55][56][57], thin films [58][59][60][61][62], plasmonic materials [50,55,[62][63][64], and, more recently, metallic metasurfaces [35,65,66] and gratings [38,64,[67][68][69], selective RHT enhancement has primarily been achieved by tuning the permittivity response, either real or effective, to create or mimic surface resonances [70]. (Notably, a recent silicon metasurface design [35] was predicted to have a larger integrated RHT than planar SiO 2 , which exhibits low-loss surface phonon polaritons, down to gap distances of 10 nm.) Other similarly high-symmetry approaches have sought to increase the photonic density of states by exploiting interference (hybridization [71]) among the localized plasmons of individual nanostrutures [72]. Building from simple shapes, tunable RHT rates have been demonstrated in nanobeam (triangular, ellipsoidal, and rectangular unit cells) and nanoantenna arrays exploiting both Mie [73] and Fano [74,75] resonances. Although conceptually promising, these approaches have been found to have diminishing returns at small separations (relative to the thermal wavelength).
In the case of periodic gratings, RHT can be expressed as
where the integration is carried out over Bloch vectors k in the first Brillouin zone (BZ), and the scattering properties of the structures are captured by the transfer function T ( ω, k) described by (6). To maximize the transfer function, both the density and coupling efficiency of the participating states [76,77] must be made as large as possible at all k. A similar problem arises in the design of far-field emitters [44,[78][79][80][81][82], where resonant structures are often intuitively designed to match the absorption and coupling rates of a wide range of externally excited states, leading to a complete suppression of scattered fields. But the task here is more complicated, as the electromagnetic field must be regulated over a much larger range of wavevectors. For intuition-based structures with a few tunable parameters, there does not seem to be enough design freedom to reach this level of control, with the range of rate-matched states occurring in current designs [36] falling short of those achievable in planar high-index dielectrics ( n > 3). Without a viable means of addressing these questions, we turn to inverse-design techniques [41].
As an initial step towards the broader development of this area, and for computational convenience, we focus exclusively on two-dimensional gratings (fields and dielectrics are translationally invariant along the third dimension). This choice has major consequences for the underlying physical processes. Particularly, in moving from three to two dimensions, the geometric and material scaling of the density of states decreases and, as a result, the maximum RHT rate between two ideal planar metals, Re[ χ ( ω )] = -2, becomes finite in the limit of vanishing loss [34],
with F ( ω ) = ( ω ) / 0 ( ω ), and 0 ( ω
the spectral emission rate per unit area of a two-(three-) dimensional planar blackbody. In contrast,
exhibits stronger geometric and material enhancement factors. Consequently, achieving a strong material response χ ( ω ) at the desired frequency window along with broadband rate matching through nanostructuring is expected to have a more significant impact in three dimensions. As (3) is a useful standard for comparing the efficacy of any given design, it will be used as a normalization throughout the remainder of this Rapid Communication.
Keeping these considerations in mind, we now describe a computational method that allows fast computations of RHT between arbitrarily shaped gratings of period , separated by a vacuum gap d . Within fluctuational electrodynamics, the calculation of RHT consists of determining the absorbed power within a body B , (
, resulting from thermally excited current sources originating within a different body A . Given a discretized computational grid and assuming local media, these sources obey the fluctuation-dissipation relation [83],
Here, ( ω, T ) = hω/ ( e hω/k B T -1) is the Planck function,
• • • a thermal ensemble average, {i, j } the index of a given location or pixel within the computational grid, and {γ , β} = {x, y, z} the vector polarizations. Equation ( 4) can be used in conjunction with knowledge of the electric Green's function G of the system [84] to obtain
Writing this in matrix form, with superscripts denoting projections onto the respective body, and G denoting the matrix form of the electric Green's function, it follows that RHT can be written as a Frobenius norm,
The main challenge in evaluating ( 6) lies in the need to repeatedly evaluate and multiply G AB , the inverse of a sparse matrix. Direct application of either sparse-direct [85] or iterative solvers [86] would demand extraordinary computational resources, especially in three dimensions. In particular, without additional simplifications, within a particular numerical discretization, the number of computations at each iteration of an optimization required is at least the rank of the matrix Im[ χ B ], or three times the number of pixels in B (polarizations). However, because G AB does not describe fields created by current sources within the same body (but only disjoint bodies), it admits a low-rank approximation [87]. Hence, ( 6) can be well approximated by a singular value decomposition of the matrix
requiring only a small set of singular values {σi }. Applying the fast randomized singular value decomposition (SVD) algorithm [88], detailed in the Supplemental Material [89], we find that typically no more than eight singular values are needed to reach an error estimate better than 1 × 10 -3 , reducing the number of required matrix inversion calculations to 16. (As derived in the Supplemental Material [89], this trace formulation also enables fast gradient computations via the adjoint method [41].)
The inverse problem is then to maximize i |σi | 2 with respect to variations in χ . Such an optimization can be carried out in the framework of topology optimization using the adjoint method [41], allowing a huge range of design parameters (each pixel within the optimization domain). We find, however, that local, gradient-based optimization leads to slow convergence to fabricable structures and comparatively suboptimal designs. To avoid these difficulties, we instead considered a range of shape optimizations [41]. While limiting the space of discoverable structures, this choice allows for application of statistical Bayesian algorithms [90,91] in combination with a gradient-based optimization. Specifically, the susceptibility profile over the periodic computational domain is described by the product χ i = χ α f α (xi ; {p α }, ζi ), where χ denotes the susceptibility of the metal, and each f α is a shape function characterized by geometric parameters {pα}. For improved convergence, f α also contains a smoothing kernel, allowing gradual variations between metal and nonmetal regions. To obtain a binary structure, the {ζi } smoothing parameters are The depicted optimization proceeds from an unstructured planar system (0) to structured gratings ( N ) by successively introducing additional ellipsoidal degrees of freedom N to the design space. Performance is qualified by the magnitude of RHT at a single design frequency ω opt = 2 πc/ 3 μ m ( λ opt = 3 μ m), where tungsten behaves as a highly lossy metal far from the planar surface plasmon resonance. The gap separation and the period of gratings are d ( ) = 0 . 02 (0 . 04) λ opt , respectively. Even in the case of a single ellipsoid, producing a relative simple grating, heat transfer is enhanced by a factor of 40, while the more complicated design (3) is observed to come within 50% of the ideal planar limit. reduced at each successive iteration of the optimization until the shape functions are piecewise constant [92]. In what follows, this procedure, implemented with a simple 2D finitedifference frequency-domain (FDFD) Maxwell solver [93], is applied to selectively enhance RHT at the thermal wavelength λ opt = 3 μ m corresponding to peak emission at T = 1000 K. The surface-surface separation between the two bodies is fixed to be d = 0 . 02 λ opt (deeply in the subwavelength regime).
We begin by considering tungsten gratings, which at λ opt exhibits a highly metallic response χ ≈ -170 + 37 i far from the planar SPP condition. Figure 1 depicts the spectral enhancement factor F/F 2D pl for both unstructured plates (0) and optimized gratings ( N ) obtained by successively increasing the number of (ellipsoidal) shapes allowed in each unit cell, N . The spectra of the two optimized gratings, illustrated as insets [94], both peak at λ opt (black dashed line), with magnitudes F/F 2D pl = { 0 . 40 , 0 . 53} increasing with the number of ellipsoids {1 , 3} [95]. On the one hand, enhancements of this magnitude for lossy metals, |χ |/ Im( χ ) 1, are considered challenging [96] at small separations, d λ opt . On the other hand, large RHT is known to be possible in ultrathin films through the interference of coupled SPPs [97]. However, to reach the magnitudes obtained here, unrealistically small thicknesses 1 nm are needed. In contrast, no feature in the gratings of Fig. 1 is smaller than 10 nm. Notably, although the optimization is carried out at a single frequency, the (There is no visible difference for gratings 6-9 and so only one of these is shown.) The lower inset depicts the RHT of tungsten grating 7 when the susceptibility is varied without altering the structure. The dashed line represents the susceptibility scaling predicted by recently derived limits on RHT [34].
discovered enhancement peaks are always broadband due to the high level of material absorption. Consequently, the frequency-integrated RHT at T = 1000 K exhibited by grating (3) is found to be roughly 10% larger than that of two planar silicon carbide (black line), a low-loss polaritonic material. To explain this enhancement, Fig. 1 (right) examines the transfer function T ( ω, k ) versus frequency and wave number k . In moving from grating 1 to 3, the color plot demonstrates (frequency axis) that additional modes are successively created and pushed towards ω opt , enhancing the density of states. Owing to the large size of the BZ (small periodicity = 0 . 04 λ opt ), we find that the range of rate matching achieved here is considerably larger than that observed in previously examined grating structures [68].
Another key finding is depicted in Fig. 2, which plots RHT enhancement for representative optimizations [98] across an array of material and geometry combinations, as a function of the material scaling factor |χ | 4 / Im[ χ ] 2 of (1). Three differ-ent classes of design are explored: collections of ellipsoids (circles), single polyline interfaces (squares), and Fourier curves (triangles). Uniformly, every one of these structures is observed to enhance RHT by at least an order of magnitude compared to the corresponding planar systems and within factors of unity of the ideal planar bound (3). Regardless of the particular parametrization, over the range of examined Re[ χ ] and Im[ χ ], a clear linear trend in F (red dotted line) as a function of |χ | 4 / Im[ χ ] 2 is observed. This linear scaling becomes increasingly difficult to observe at larger values of χ , where larger resolutions are needed to accurately capture resonant behavior and the optimization requires increasingly larger numbers of iterations to find structures along the fit line. As should be expected, material scaling consistent with (1) is seen only for gratings optimized specifically for each particular value of χ ( ω opt ). For a fixed geometry [Fig. 2 (insets)], varying either Re[ χ ] or Im[ χ ] shifts the resonance away from ω opt and diminishes RHT. The appearance of this linear trend indicates that aspects of the arguments leading to the energy bounds in Ref. [34] are coming into play. However, the fact that the slope of the fit line is too flat to pass through the origin indicates that a substantial increase in available design space has failed to significantly bridge the magnitude gap, and that approaching these bounds (at least in 2D) may prove challenging with practical structures. Conversely, it should be emphasized that our results do not preclude the existence of structures with much larger enhancements. First, although the design space we have investigated is substantially larger than previous work, it is still relatively limited. Second, the complexity of the structures and degree of enhancement are limited by the spatial resolution of the chosen discretization, 0 . 0005 λ opt (1 / 40th of the gap size).
To summarize, in investigating potential radiative heat transfer enhancements through inverse design, we have found evidence supporting the material scaling recently predicted by shape-independent bounds [34]. While the observed heat transfer rates are still far from matching the magnitudes predicted by general bounds, we have found that RHT rates between fabricable tungsten gratings (a highly lossy metal), for subwavelength gap separations as small as 2% of the design wavelength, can selectively approach 50% of the rate achieved by ideal planar materials (lossless metals satisfying the SPP condition) in the infrared. The results represent nearly two orders of magnitude greater RHT rates in structured compared to planar materials. It remains to be seen to what degree similar strategies might enhance RHT in three dimensions, where the photonic density of states is larger and its associated dependence on material losses significantly stronger.