Direct modeling of near field thermal radiation in a metamaterial

: The study of near field thermal radiation is gaining renewed interest thanks in part to their great potential in energy harvesting applications. It is well known that plasmonic or polaritonic materials exhibit strongly enhanced fields near the surface, but it is not trivial to quantitatively predict their impact on thermal radiation intensity in the near field. In this paper, we present a case study for a metamaterial that supports a surface plasmon mode in the terahertz region and consequently exhibits strongly enhanced near field thermal radiation at the plasmon resonance frequency. We implemented a finite-difference time-domain method that thermally excites the metamaterial with randomly fluctuating dipoles according to the fluctuation-dissipation theorem. The calculated thermal radiation from the metamaterial was then compared with the case of optical excitation by the plane wave incident on the metamaterial surface. The optical excitation couples only to the mode that satisfies the momentum matching condition while thermal excitation is not bound by it. As a result, the near field thermal radiation exhibits substantial differences compared to the optically excited surface plasmon modes. Under thermal excitation, the near field intensity at 1 µm away from metal surface of the metamaterial reaches a maximum enhancement of 43 fold over the far field at the frequency of the Brillouin zone boundary mode while the near field intensity under optical excitation reaches a maximum enhancement of 24 fold at the frequency of the Brillouin zone center mode. In addition, the peak near field intensity under thermal excitation shows a 4-fold enhancement over blackbody radiation with linear polarization radiation in the far field. The ability to precisely predict the local field intensity under thermal excitation is critical to the development of advanced


Introduction
Thermal radiation has been an intriguing topic over the past century.The study of blackbody radiation stimulated the development of quantum theory, which opened a new era of human history.The advance of nanotechnology in the past few decades has enabled researchers to study thermal radiation in the near field which exhibits many new phenomena, including coherent thermal radiation [1], surface wave excitation by a thermal source [2], and near-field radiation exceeding the blackbody radiation limit [3,4].This new research frontier has inspired a variety of applications, including thermal energy conversion [5], radiative cooling [6], nano-fabrication [7,8], and near-field imaging [9].
The thermal radiation research in nanostructures calls for a rigorous and efficient way of numerical modeling.In some studies [6], optical excitation by the far field plane wave is used as the excitation source.For studying far field thermal radiation, it is sufficient and computationally efficient to calculate absorption under far field illumination and invoke the Kirchhoff's law of thermal radiation.But the far field optical excitation method has the fundamental limitation that it can only excite the modes that satisfy the momentum matching condition.Modes that are not excited by optical excitation but that can be excited thermally could lead to significant deviations of the near field thermal radiation from the results under optical excitation.Even for the modes that are excited optically, the coupling efficiency may vary from geometry to geometry and is generally different from the case of thermal excitation.To accurately simulate the thermal radiation, therefore, one must simulate the radiation resulting from thermal fluctuations using the fluctuation-dissipation theorem [10].This approach has been implemented in both time domain [11,12] and frequency domain [13] to directly calculate thermal radiation from complex structures but these techniques have not yet been fully utilized to investigate the near field thermal radiation in artificially structured materials.
In this paper, we present the results of numerical study of the near field thermal radiation from a simple metamaterial thermal source using both far field plane wave optical excitation and thermal excitation.The metamaterial is designed to support surface plasmon polariton (SPP) modes at 5 THz at the Brillouin zone center and 4.6 THz at the Brillouin zone boundary, respectively.In the optical excitation model, only the Brillouin zone center mode can be excited, and exhibits a near field intensity enhancement peak at 5 THz with an enhancement factor of 24.However, the actual thermal radiation extracted from the thermal excitation model reveals a peak near field intensity enhancement at 4.6 THz with an enhancement factor of 43.From the thermal excitation simulations, we also calculate the energy density radiating from this metamaterial thermal source at a temperature of 600 K.The near field intensity was found to exceed the same temperature blackbody radiation intensity by approximately a factor of 4 in the frequency band of the surface plasmon mode.Our study confirms that a metamaterial structure supporting surface plasmon mode can exhibit near field thermal radiation intensity exceeding the blackbody radiation.It also shows the differences between optical and thermal excitation modeling and thereby the need for implementing direct thermal modeling for quantitative study of near field thermal radiation phenomena and their device applications.

Results and discussion
To design a 1D periodic metamaterial thermal source, we first use the transverse magnetic (TM) polarized far field plane wave as the excitation source, and tune the reflectance spectrum so that the dip occurs at 5 THz.We choose this frequency band with an eye on eventual application in the rectenna technology [14].A rectenna directly rectifies the ac electric field of an incoming electromagnetic wave and produces dc electricity.One of the main challenges for the energy harvesting application is that the rectenna efficiency tends to be low at high operating frequencies while the thermal radiation intensity is low at low frequencies.A metamaterial thermal emitter offers a solution by significantly increasing the thermal radiation intensity in the near field at low frequencies.In this context, 5 THz is considered a good operating frequency.The simulation is conducted by the finite-difference time-domain (FDTD) method using the commercial software Lumerical.The metamaterial structure schematically shown in Fig. 1(a) consists of 1D periodic trenches in copper.The trenches are filled with SU-8 and there is also a thin layer of SU-8 on the metal surface.In the FDTD simulation of far field optical excitation, we set up a two-dimensional (2D) computational cell containing one unit cell along the direction of periodicity and terminate with periodic boundary condition (PBC).Above the metamaterial structure is an air region, which is terminated with the perfect matched layer (PML) boundary condition that absorbs all out-going waves exiting the computational domain.As for the material properties, we use the Drude dielectric function for copper [15], and a constant refractive index of 1.7 for SU-8 [16].In the final design, the period is 24 μm with a trench depth of 7.3 μm and width of 2 μm.The SU-8 overlayer is 1 μm thick.In the future integration with a photovoltaic (PV) converter, the PV device will be placed on the SU-8 overlayer and therefore we evaluate the field in the following discussion is at the SU-8 surface or 1 μm above the metal surface.
The reflectance spectrum of the metamaterial under TM polarized plane wave excitation is shown in Fig. 1(c).It has a strong resonance at 5 THz with reflectance suppressed down to 20%.This indicates good coupling between the plane wave excitation source and the surface plasmon mode.To further confirm the excitation of the surface plasmon mode, we calculate the photonic band structure of the metamaterial also using the FDTD method.In this calculation, we place multiple random electric dipoles near the metamaterial surface, and monitor the modes that persist long after the dipoles are turned off.The result is shown in Fig. 1(d).The wave number k y is for the direction along the periodicity (y direction as shown in Fig. 1(a)).The maximum value of k y is limited to the first Brillouin zone edge, / π Λ , where Λ = 24 μm is the periodicity.The straight line, which bounces back at the Brillouin zone boundary, is the light line in free space.The flat line that crosses the y-axis at 5 THz is the surface plasmon mode and the y-intercept corresponds to the resonance dip observed in the reflectance spectrum under normal incidence or k y = 0.At the Brillouin zone edge, the surface plasmon mode is located at 4.6 THz.The curve overlaps with the light line at low frequencies, and bends to the right side at higher frequencies until it reaches the resonance frequency.This surface plasmon polariton branch below the light line has the ability to produce strongly localized optical modes at desired frequencies and is sometimes called spoof surface plasmon [17].Near the resonance frequency, the slope of the dispersion curve decreases close to zero, indicative of a large local density of states and local electromagnetic field.However, due to the momentum mismatch, these modes cannot be excited by coupling light from free space.Thus, no resonance dip is observed around the spoof surface plasmon frequency of 4.6 THz in the reflectance spectrum in Fig. 1(c).To investigate the near field thermal radiation in a rigorous and quantitative manner, we implemented a direct thermal modeling method based on fluctuation-dissipation theorem following the approach proposed by Luo et al [11,18,19].Briefly, we introduced a random thermal fluctuation term ( )  (1) where γ is the damping coefficient, and 0 ω , σ are the resonance frequency and strength, respectively.Based on the fluctuation-dissipation theorem [20], the time correlation of ( ) t K in frequency domain is given as follows: where ( , ) I T ω represents the free space Planck blackbody radiation intensity at temperature T, K α β are the components of vector K, r is the spatial coordinate and ω is the angular frequency.Taking advantage of the linearity of electromagnetism in a linear material system, we replace the time correlated function should be implemented in 3-dimensional (3D) space, we used 3D FDTD modeling even though the structure is invariant in z direction.In the z direction, we extended the length to 2 μm, and set boundary condition to continuity.The cross-section on the x-y plane is shown in Fig. 1(b).To faithfully represent the random thermal fluctuations, we set up five simulations with different Bloch phase shifts along the y direction.The phase shift was varied at uniform interval from 0 to / π Λ , where Λ = 24 μm is the period of the metamaterial structure.For each Bloch phase shift, we averaged the results of seven simulations with different random number seeds to reduce the stochastic noise inevitable in this method [11].The final result is an average of all simulations with different Bloch phase shifts and random number seeds used to generate the ( ) t ′ K term.Taking into consideration the computational time and the FDTD simulation frequency resolution, we chose the FDTD time long enough to reach the temporal steady state and to have 0.2 THz resolution in frequency.The simulation is first tested for thermal emission from a flat aluminum surface for which an analytical solution exists.The simulation showed an excellent agreement with the analytical result.The details are presented in the Appendix.
For quantitatively comparisons between the optical and thermal excitation of surface plasmon mode, we first present the optical excitation results.In the following discussion, the effect of surface plasmon will be assessed by comparing the optical energy density, ε|E| 2 , and will be loosely termed as intensity enhancement, as it is proportional to the actual intensity in the far field, although one cannot rigorously define intensity for evanescent modes.Figure 2(a) shows the profile of intensity enhancement over the excitation source at 5 THz, the frequency of the surface plasmon at k y = 0.The intensity is largest near the trench and decays away from the metal surface and also from the trench center.At the largest intensity point, the intensity enhancement factor over the incident intensity reaches about 300.Inside the metal, the field is almost zero because copper has a very large negative permittivity and thus behaves like a perfect electric conductor in this frequency range.Since the local field intensity varies strongly along the y-direction, we averaged the field over one unit cell length for comparison's sake.The averaged intensity enhancement at 1 μm away from the copper surface is shown in Fig. 2(c) from frequency 3 THz to 7 THz.It has a clear peak at 5 THz with an enhancement factor of 24 over the far field intensity.The peak position aligns well with the frequency of the reflectance dip as well as the zone center mode (k y = 0) in the band structure.This is the well-known excitation of surface plasmon mode at k y = 0 via grating coupling.Similarly to the optical excitation case, we calculated the enhancement of near field optical intensity under thermal excitation.Once again, the near and far field intensities were calculated by averaging the intensity over the unit cell along y direction as shown in Fig. 1(b).The averaged near field intensity was evaluated at 1 μm above the metamaterial surface and the averaged far field intensity at 100 μm away.We confirmed that all near field modes decay sufficiently at 100 μm distance away from the metamaterial.Figure 2(b) shows the field profile of the thermally excited surface plasmon mode at 5.1 THz.Since no visual difference of intensity is observed in the x dimension from 30 μm and farther, the figure shows only the x range from −10 μm to 30 μm.Limited by our FDTD frequency resolution, 5.1 THz was the closest frequency to 5 THz, the zone center mode frequency at k y = 0.It is apparent that the mode shape is the same as that shown in Fig. 2(a) under far field TM polarized plane wave excitation.In both cases, most of the energy is confined near the trench and decays away from it.In direct thermal excitation, the plasmon mode is excited thermally, resulting in a different near field enhancement factor.Figure 2(b) shows strong field with randomness in intensity and spatial distribution in the copper region, which arise from the randomly fluctuating dipoles.In contrast, under far field plane wave excitation, no field resides inside the metal as shown in Fig. 2(a).Also, since the thermal source can couple into all available modes, the intensity enhancement profile in Fig. 2(b) must contain all surface plasmon modes with different k y values and polarizations.
Figure 2(d) shows the averaged near field intensity enhancement over the far field intensity.At the frequency of the zone center mode of 5 THz, the enhancement is around 6 fold.It is much weaker than that in the optical exaction case shown in Fig. 2(c).This means that the far field TM polarized plane wave excitation model overestimates the available thermal energy at the 5 THz zone center mode frequency.In contrast, the largest enhancement happens at the zone boundary mode frequency of 4.5 THz with an enhancement factor of 43.This peak enhancement at 4.6 THz does not show up in the optical excitation model because the far field plane wave excitation source cannot couple to the zone boundary mode due to momentum mismatch.In contrast, a thermal source can excite surface plasmon modes with various k y values and produce large near field thermal radiation over the whole surface plasmon frequency band.As shown in the band structure in Fig. 1(d), the first branch of the surface plasmon dispersion curve lies below the light line with wave number larger than light line.Therefore, the excited modes are all surface waves.The thermal energy coupled into these waves is trapped near the metamaterial surface and eventually gets absorbed back by copper.The mode profile of the zone boundary mode is shown in Fig. 3(b).The mode shape is very similar to the 5 THz zone center mode in Fig. 2(b), but the near field intensity enhancement over far field is much stronger.
All the previous discussions are focused on the relative thermal radiation intensity enhancement of the near field over the far field.It is also important to investigate the actual radiative energy that is produced by the metamaterial thermal source, and compare it with the same temperature blackbody.The blackbody radiation is the upper limit of thermal radiation to the far field from any object.It is ultimately limited by the free space density of state of propagating waves which convey thermal energy into the far field.In the near field, thanks to the additional surface modes, the density of states can exceed that of the free space and consequently the radiative energy exceeds that of blackbody radiation [4,21].It is stressed that this phenomenon happens only in the near field region and in the far field region the thermal radiation still obeys the Kirchhoff's law.To make a specific comparison with blackbody, we chose a temperature of 600 K.The thermal excitation we implement here does not explicitly include temperature information in it.We therefore use Kirchhoff's law of thermal radiation to specify the temperature.Specifically, using the emissivity calculated by the far field optical excitation simulation, we calculate the far field radiation intensity spectrum for the temperature 600 K.By equating this value with the far field intensity extracted from direct thermal excitation modeling, we obtain the frequency dependent scaling factors for this specified temperature over all frequencies of interest.We then use the same scaling factors to calibrate the near field.Once again, we use ε|E| 2 as relative intensity to directly compare the optical field in the near and far field regions.In Fig. 3(a), the black dashed line is the blackbody radiation intensity at 600 K.The red line is the far field radiation intensity from the metamaterial thermal source, which is the product of emissivity and the blackbody radiation intensity.Since the emissivity cannot exceed unity at any frequency, the far field radiation intensity is always weaker than blackbody radiation over the whole frequency band.It has a peak at around 5 THz due to the large emissivity arising from the out-coupling of the zone center mode.The blue line is the averaged relative intensity at 1 μm away from the metamaterial surface.It exceeds the blackbody radiation over the frequency band from around 4.2 THz to 5.2 THz due to the excitation of surface plasmon modes.The peak intensity is located near the zone boundary mode frequency of 4.6 THz with approximately 4 times higher intensity than the blackbody radiation.Judging from the intensity enhancement profile in Fig. 3(b), we expect even stronger radiation at positions closer to the copper surface.We also studied the polarization of the thermal emission from the metamaterial source.Here, we consider only the periodic boundary condition along y direction for the 0th diffraction order, which corresponds to the Brillouin zone center mode (k y = 0) in the band structure.As shown in Fig. 4(b), the far field radiation is strongly polarized in the y direction, especially near the surface plasmon resonance frequency of 5 THz.The intensity in the y polarization direction is over 1000 times stronger than the intensity in the x polarization direction.It indicates that the thermal radiation at this frequency is dominated by the excitation and subsequent out coupling of the surface plasmon mode.This is in good agreement with its reciprocal optical excitation, where a TM polarized source has a very strong coupling with the Brillouin zone center mode and consequently strong absorption at the surface plasmon frequency while a TE polarized source has very poor absorption.The polarization of the thermal radiation in the near field as shown in Fig. 4(a) is very different due to the presence of bound evanescent waves near the metamaterial surface.The ratio of intensity in y polarization over that in x polarization also reaches the peak at around 5 THz due to the Brillouin zone center mode.However, the value of the intensity ratio evaluated at 1 μm above the metal surface is much smaller, ranging between 1.5 and 2.5.This is in excellent agreement with the optical excitation simulations which show the polarization intensity ratio of 2.5 at 5 THz and confirms that the thermally excited plasmon modes retain the same polarization characteristics as the optically excited modes.

Fig. 1 .
Fig. 1.(a) 3D sketch of the metamaterial thermal emitter.Cu substrate and SU-8 coating are colored yellow and pink, respectively.(b) The unit cell of the metamaterial thermal emitter.The structural periodicity is 24 μm along y direction as denoted in the picture.Direction z is out of the plane. 1 μm thick SU-8 is coated on the copper substrate and fills the groove.Other parameters are stated in the text.Copper is in yellow, SU-8 is in red, and vacuum is in black.(c) The reflectance spectrum of the metamaterial structure under far field TM polarized plane wave excitation with normal incident angle.(d) The bandstructure of metamaterial in logarithmic color scale.The color bar represents the coupling strength of electric dipole to the modes.k y is the mode wave number along the periodicity in y direction.
N is the FDTD time steps, and V Δ is the mesh volume element used in the simulation.The final simulated thermal radiation is calibrated by requiring the far field frequency spectrum to be consistent with the Kirchhoff's law.To implement this scheme on the commercial FDTD software Lumerical, we worked with the Lumerical engineers to develop a new material model using their "flexible material plugin framework".The new material model simulates thermal radiation by adding random perturbation ( )t ′ Kto the polarization at each FDTD step.Because the random excitation term ( )t ′ K

Fig. 2 .
Fig. 2. (a) Electromagnetic intensity enhancement over far field TM polarized plane wave excitation at 5 THz.The white dash line depicts the copper surface.(b) Near field thermal radiation intensity enhancement over far field at 5.1 THz under thermal excitation.The enhancement strength is shown by the color bar.(c) Averaged near field (1 μm above copper surface) intensity enhancement under far field TM polarized plane wave optical excitation.(d) Averaged near field (1 μm above copper surface) intensity enhancement under thermal excitation.

Fig. 3 .
Fig. 3. (a) Averaged relative intensity at 100 μm away from the metamaterial thermal source surface (red line), at 1 μm above the copper surface of metamaterial (blue line), and far field free space radiation intensity of the blackbody thermal source (black dash line).All thermal source temperatures are set to 600K.(b) Thermal radiation intensity enhancement over far field at 4.7 THz frequency.

Fig. 4 .
Fig. 4. Polarization of the thermal radiation in the near field (a) and in the far field (b).Ix and Iy are the intensity in the x and y direction, respectively.