Near- and far-field thermal radiation in plasmonic metamaterial emitters studied via fluctuation–dissipation FDTD

Abstract

Renewed interest in near-field thermal radiation stems from its potential in energy harvesting, nanoscale thermal management, and thermophotovoltaics, driving efforts to understand and engineer its fundamental mechanisms. Plasmonic materials are well recognized for their ability to enhance fields proximal to their surfaces; however, detailed quantitative assessment of the field intensity and energy density of thermal radiation in the near field (relative to the far field) and potential practical impact is lacking. This work presents a detailed study of a prototypical metamaterial design supporting plasmon resonance in the mid-infrared spectrum. Significant intensification of thermal radiation in near field is investigated using a numerical finite-difference time-domain method that incorporates randomly fluctuating dipoles in accordance with the fluctuation-dissipation theorem that enables excitation of high-order spatial harmonics. Electromagnetic response of the metamaterial thermal emitter to both thermal and optical excitations is evaluated and compared. Findings indicate that the near field intensity at a mere 500 nm from the metamaterial surface exhibits an enhancement of up to 38-fold in near-field electromagnetic energy density relative to the corresponding far-field value at 4.8 μm lying near the edge of Brillouin zone under thermal excitation. In contrast, 6-fold enhancement is observed at 4.35 μm lying near the center of the Brillouin zone under optical excitation. Further, thermal excitation yields a 3.5-fold increase in peak intensity compared to conventional blackbody radiation at 2 μm away from the outermost layer of the metamaterial. Moreover, photothermal simulations further reveal rapid temperature rises exceeding 600 K within nanoseconds under low-power pulsed excitation. These results demonstrate that thermal excitation enables access to dark electromagnetic modes that dominate near-field radiation and highlight the importance of excitation mechanism in designing metamaterial-based thermal photonic devices.

Graphical abstract

Keywords

Near-field thermal radiation / Plasmonic metamaterials / Fluctuation–dissipation theorem / Brillouin zone resonance / Dark electromagnetic modes / Photothermal response

1 Introduction

Thermal radiation is a fundamental electromagnetic process governing radiative heat transfer across length scales ranging from macroscopic systems to the nanoscale. Over the last hundred years, investigations into fundamental science and application of thermal radiation have been remarkably successful [1]. These studies, motivated by the need to improve effectiveness and productivity of large-scale thermal systems, have recognized the need for an in-depth understanding of the underlying microscopic mechanisms [2]. A grasp of how thermal radiation works at close distances is vital in many areas of inquiry, which include energy harvesting [3], thermal rectification [4], tailoring coherent thermal emission [5], radiative cooling [6], and surpassing the blackbody radiation limit in the near-field [7,8]. The use of advanced numerical modeling techniques, along with application of nanotechnology, have been essential for the accurate study of thermal radiation profiles of nanostructured devices. Most computational studies of thermal emission rely on optical excitation combined with Kirchhoff’s law to infer emissivity and far-field radiation. While effective for propagating modes [9−11], this approach inherently excludes non-radiative electromagnetic states that lack momentum matching to free space. As a result, near-field thermal radiation—where evanescent and dark modes can dominate—cannot be accurately captured using optical excitation alone. This limitation motivates the use of fluctuation–dissipation-based methods that explicitly account for thermally induced electromagnetic fluctuations. Thermal excitation can stimulate modes that are not triggered by optical means, which in turn can lead to significant discrepancies in near-field thermal radiation in contrast to that due to optical excitation [12]. Herein, we investigate nanoscale structures from both optical and thermal design perspectives with the underlying objective of modulating near and far field intensities in thermally activated systems.

At thermal equilibrium, a blackbody emits electromagnetic radiation spontaneously and continuously with a spectral dependence described by Planck’s and Kirchhoff’s laws such that with increasing temperature the radiant energy flux increases while the peak wavelength of the distribution decreases [13]. The prospect of achieving near-field–enhanced emitter has been investigated extensively over the last several decades with an eye to engineering surfaces where spectral radiance can exceed the blackbody limit over a certain range of frequencies albeit without violating the laws of thermodynamics [14]. Under equilibrium conditions, double-negative metamaterials with subwavelength structuring held at a fixed temperature have been shown theoretically to exhibit spectral emissivity exceeding unity over a certain range of wavelengths. Further, the study indicates that a metamaterial held at a sufficiently low temperature can perform as a “thermal black hole”, displaying a spectral absorptivity that exceeds unity [15]. These near-field–enhanced emittance and absorption outcomes are—associated with high-order highly-reactive spatial harmonics, predominantly dark modes, of the emitter’s fluctuating field—in conjugate-matched metamaterials which would otherwise not be evident in a simple planar absorbing material [16]. Thermal emission properties of such nanophotonic structures have been investigated using the fluctuation dissipation theorem in both frequency and time domains [17,18], however, there is a scarcity of comprehensive investigations of thermal radiation of metamaterials in both the near-field and far-field.

Recent advances in fluctuational electrodynamics have substantially expanded the numerical tools available for predicting thermal radiation from nanostructured and metamaterial emitters. Modern formulations now include fluctuating-volume-current and volume-integral approaches for inhomogeneous bodies [19,20] direct near-field radiative-transfer FDTD algorithms [21], Wiener-chaos/FDTD approaches for inhomogeneous thermal emitters [22] discrete-system Green’s-function methods for irregular three-dimensional geometries [21], and many-body/scattering-based frameworks for complex radiative heat-transfer systems [23]. These developments are important because near-field thermal radiation is governed not only by propagating modes, but also by evanescent, surface-polaritonic, and high-spatial-frequency electromagnetic states whose contribution cannot be fully inferred from far-field absorptivity alone [24−26].

Within this context, the FDT-FDTD approach adopted here is particularly suitable for the present multilayer plasmonic metamaterial because thermally fluctuating dipoles can be introduced directly into Maxwell’s equations throughout the dispersive metal-dielectric volume which builds on the work of Chan et al. [27] it models 3D periodic structures and further Maxwell’s equations are modified so as to consider all spatial harmonics to directly compute near- and far-field thermal radiation from a plasmonic metamaterial emitter. By comparing thermally and optically excited responses within the same structure, we isolate the role of dark electromagnetic modes in governing near-field enhancement. In addition, we couple electromagnetic absorption to transient heat transfer simulations to examine the photothermal response under pulsed excitation, providing a unified framework for thermal photonic device design.

2 Theory

2.1 Fluctuation-dissipation-based FDTD formulation

Materials exhibit spontaneous electrical and magnetic moments due to quantum and thermal fluctuations which generate fluctuating electromagnetic fields both internally and externally that manifest as thermal radiation [28]. These fields can be modeled using Maxwell’s stochastic equations—essentially Maxwell’s equations augmented with random electric and magnetic current sources, akin to the approach for analyzing Brownian motion. Concurrently, the fluctuation-dissipation theorem (FDT) helps in calculating the mean spatial correlation of these random sources. Thus, the phenomenon of thermal radiation can be effectively understood through the integration of Maxwell’s stochastic equations with FDT.

This simulation of thermal electromagnetic emission is a finite-difference time-domain (FDTD) method that integrates the established Langevin method for modeling Brownian motion. Maxwell’s equations, in their current form, represent classical deterministic field equations. By introducing an element of randomness into these equations it is possible to generate the natural variability found in thermal fluctuations. Typically, there are three strategies by which randomness is introduced: direct integration of randomness in Newton’s equation of motion, incorporating a stochastic element in the displacement field, or embedding a random component in the free current density. Here we introduce randomness in the displacement field, D. The polarization response of a dispersive and absorbing medium is modeled using a Lorentz oscillator formalism, which forms the basis for incorporating thermal fluctuations within the FDTD framework [29].

(1)\frac{{\text{d}}^{2}P}{\text{d}{t}^2}+\gamma \frac{\text{d}P}{\text{d}t}+{\omega }_0^{2}P=\sigma E,

where γ, ω₀, and σ represent the damping coefficient, resonance frequency, and conductivity, respectively.

Considering thermal fluctuation, a random factor K(r,t) is introduced on the right-hand side of the above relation using the Langevin technique.

(2)\frac{{\text{d}}^{2}P}{\text{d}{t}^2}+\gamma \frac{\text{d}P}{\text{d}t}+{\omega }_0^{2}P=\sigma E+K(r,t).

Using the general harmonic solution, we obtain the solution for P:

(3)P(r,\omega )=\frac{\sigma E(r,\omega )}{{\omega }_0^2-{\omega }^2-\text{i}\gamma \omega }+\frac{K(r,\omega )}{{\omega }_0^2-{\omega }^2-\text{i}\gamma \omega }.

The displacement field, D, is related to polarization through the relationship D = E + 4πP, which comprises the external field E, a conventional non-stochastic component due to polarization, and an additional random element Q(r,\omega ) that we represent as follows:

(4)Q(r,\omega )=\frac{\text{4π}K(r,\omega )}{{\omega }_0^2-{\omega }^2-\text{i}\gamma \omega }.

2.2 Statistical attributes of thermal fluctuation

In order to analyze the fluctuating displacement field, the correlation function for Q needs to comply with the fluctuation-dissipation theorem as formulated by Levin et al. [30]:

(5)\begin{array}{rl}⟨Q_i(r,\omega )Q_j^{\ast }(r^{\prime },{\omega }^{\prime })⟩=& \frac{16{\text{π}}^3{c}^2\text{Im}\left[\epsilon \left(\omega \right)\right]}{{\omega }^3}\\ & \times I_0(\omega ,T){\delta }_{ij}{\delta }_{\omega {\omega }^{\prime }}\delta (r-r^{\prime }).\end{array}

Equation (5) describes the spatial and spectral correlation of thermally fluctuating displacement fields, indicating that thermal radiation originates from stochastic polarization currents distributed throughout the material volume. Importantly, these fluctuations excite all allowed electromagnetic states of the system, including non-radiative and highly confined modes that cannot be accessed through plane-wave optical excitation. where i=1,2,3 correspond to the components of Q, 〈...〉 denotes ensemble averaging, c is the speed of light, Im[ε(ω)] is the imaginary part of the permittivity including the polarization response in the absence of fluctuations, and I_0(\omega ,T)=(c\text{/(4π))}D(\omega )E(\omega ,T). In this expression, D\left(\omega \right)={\omega }^2/({\text{π}}^2{c}^3) is the free-space density of photon states and E(\omega ,T)=\hslash \omega /[\exp (\hslash \omega /(kT))-1] is the Bose−Einstein energy distribution function at absolute temperature T. The appearance of I_0(\omega ,T) in Eq. (5) exhibits the universality of the thermal fluctuations consistent with Planck’s formula and Kirchoff’s law. In the discrete limit, the fluctuations in K thus become [29]

(6)⟨K_i\left(r,\omega \right)K_j^{\ast }\left(r^{\prime },{\omega }^{\prime }\right)⟩=\frac{4{\text{π}}^2\sigma \omega \text{γ}}{{(\omega /c)}^2\text{Δ}V}{I}_0(\omega ,T){\delta }_{ij}{\delta }_{\omega {\omega }^{\prime }},

where \text{Δ}V represents the volume element employed in the simulation. To determine emissivity, the thermal emission intensity given in Eq. (6) is normalized against the free-space Planck radiation across the emission solid angle. The Bloch periodic boundary conditions are specified by the directions by k_x and k_y instead of the polar an azimuthal angle. In calculating the emissivity all frequency-dependent factors in the equation are effectively eliminated [29]. By exploiting the linearity of electromagnetism in linear material systems, the time-correlated function K(t)K^{\prime }(t^{\prime }) can be substituted with a stochastic excitation term to improve computational efficiency. In this approach, K′(t) is treated as a random variable with a uniform distribution over the interval ±\sqrt{\text{12π}\sigma \gamma /\left(N\text{Δ}V\right)}, where N represents the total number of time steps used in the Fourier transform and \mathrm{\Delta }V is the mesh element volume.

3 Design methodology and analytical treatment

With an eye to designing a representative metamaterial, we consider silica, the most commonly occurring dielectric material, and nickel is selected as the plasmonic constituent due to its thermal stability, mid-infrared plasmonic response, and compatibility with high-temperature operation, in contrast to noble metals which exhibit lower melting points and reduced thermal robustness, and commonly use in catalytic systems with economic viability. The cubic meta-atom geometry is chosen to support multiple resonant pathways and sharp field confinement while maintaining symmetry compatible with periodic boundary conditions. One of the vibrational modes of carbon dioxide [31]. The proposed design of the metamaterial thermal emitter is illustrated in Fig. 1a. The meta-atoms comprise of three nickel and silica bilayers of 25 nm and 50 nm thicknesses, respectively, with a cross-sectional profile of 2.7 μm square. The meta-atoms are distributed in a square lattice with a pitch of 3.2 μm. Further the meta-atoms lie on a silica sheet of 0.2 μm which in turn rests on a nickel ground plane of 0.3 μm thickness; the former is a spacer layer that supports resonant modes while the latter is a back reflector. The entire structure is supported on silicon. The optical properties of SiO₂ and Ni are taken from Malitson [32] and Palik [33], respectively.

The metasurface is analyzed for its excitation modes through optical and thermal excitation techniques described above. The optical parameters of the absorber denoted T, R, and A represent the spectral transmission, reflection, and absorption, respectively. We use impedance matching to maximize absorption and minimize reflection by aligning the impedance of the metasurface, Z, with that of free space (Z₀ = 375 Ω). We examine the dependence of the impedance Z on the S parameters (S₁₁ and S₂₁) of the scattering matrix (where R\left(\omega \right)={\left|S_{11}\right|}^2, T\left(\omega \right)={\left|S_{21}\right|}^2) and A\left(\omega \right)=1-R\left(\omega \right)-T\left(\omega \right)), and in particular to minimize the parameter S₁₁ [34].

The proposed structure, upon excitation by a plane wave optical source, exhibits a dip in reflection with a minimum in the mid-infrared region at λ = 4.35 μm, which corresponds to the asymmetric stretching C−O mode (see Fig. 1b). The reflectance at the resonant wavelength is reduced to a low value of 10%, indicating that this impedance aligned structure effectively supports excitation of plasmonic modes. To corroborate this observation, we calculate the metamaterial’s photonic band structure. The metasurface is stimulated with every conceivable mode of the system through the deployment of numerous randomly positioned broadband dipoles. The wave vector k is defined by the Bloch boundary conditions, necessitating one simulation for each k-vector; where a mode (or band) is present, the fields will continue to propagate indefinitely. At frequencies not associated with any mode, the fields rapidly dissipate due to destructive interference. Thus, by observing the resonant frequencies where the field persists in the simulation, we determine the band which is shown in Fig. 1c. The wave number k_y in a meta-unit cell (Fig. 1a) extends from the centre (value 0) to the edge of the first Brillouin zone where it has the maximum value of π/Λ, where Λ is the periodic spacing of 3.5 μm. An essentially constant horizontal line intersecting the y-axis at 4.35 μm represents the plasmonic mode of the structure, which is in agreement with the resonant dip observed in the reflectance spectrum (viewed under normal incidence). Also, there exists a pronounced mode along the diagonal, increasing in wavelength until it reaches the edge of Brillion zone. These polariton modes spanning across the band structure are able to produce intensely localized optical modes in the said frequency range, which are occasionally termed as dark modes [15]. However, these modes cannot be excited by direct coupling with free-space light due to momentum mismatch. This is confirmed by the lack of a resonant dip around the plasmon mode at 4.8 μm (Fig. 1b). In addition, Brillouin-zone momentum-matching analysis shows that the resonance near λ = 4.8 μm can couple to an off-normal far-field channel near θ = 43.3º, while the remaining high-k components contribute primarily to localized near-field enhancement rather than direct far-field emission (see Fig. S2).

To gain a thorough understanding of the structure’s performance, we analyzed the electric field behavior under both optical and thermal excitations. As can be seen in Fig. 1e, when the metamaterial structure is heated, the radiation emitted from all parts of the structure, including the surface and sharp edges of the nickel (Ni) layers, is due to the thermal properties and behavior of the materials. That is, at temperature all materials, Ni and SiO₂, emit thermal radiation based on their temperature and emissivity characteristics, and as such results in the presence of electromagnetic fields within the materials. As observed, a strong electric field appears at the outermost dielectric surface of the structure while away from the structure the highest field amplitude is observed at ~1 μm from the top of the meta-atom—exceeding the far-field intensity. Additionally, constructive and destructive interference, in accordance with Huygens’ principle, can be observed. This suggests that the thermal excitation leads to a more uniform and widespread distribution of the electric field across the various components of the structure, in contrast to the localized field enhancements observed under optical excitation.

However, under optical excitation (Fig. 1d), where the structure is excited by an incident electromagnetic wave, the response is primarily determined by the interaction of the refractive index of the materials with the specific wavelength of the optical source. In this case, the Ni layers – being metallic – for the most part will not exhibit significant electric field accumulation due to its higher conductivity and lower permittivity, whereas the dielectric properties of SiO₂ allow for stronger interaction with the incident electromagnetic wave at specific wavelengths giving rise to standing waves and hence concentration of the electric field [35]. Hence, the electric field accumulates principally at the sharp edges and on the surface of the Ni layers with minimal buildup of charges inside the Ni layers. The resonant modes supported by the dielectric material lead to field localization and enhancement in the SiO₂ regions. We emphasize that the behavior of the electric field within the metamaterial is highly dependent on the nature of the excitation. Optical excitation results in resonant field enhancements primarily within the dielectric SiO₂ layers owing to its interaction with the incident electromagnetic waves. Conversely, thermal excitation leads to a relatively uniform distribution of the electric field considering that the thermal radiation is defined by the thermal properties of the constituent materials. Understanding these differences is crucial for optimizing the performance of metamaterial structures for specific applications and operating conditions – optical or thermal.

We now investigate the response of the metasurface to thermal excitation, and in particular examine the near-field and far-field optical energy density represented by the term ε|E|². The contrast in the field densities is used to obtain the enhancement in the near-field calculated on the basis of the actual intensity in the far-field. Enhancement in the spectral intensity at several distances perpendicular to the plane tangential to the top-surface of the meta-atoms is shown in Figs. 1g−i. Given substantial variations in the local field intensity along the z-direction, the spectral field distribution shown at z is the average of the field extending from 0 to z for a meta-unit cell. The spectral field distributions at 500 nm, 2500 nm, and 80 μm away from the surface of the unit cell correspond to average intensities in the near, intermediate and far fields, respectively. The substantial enhancement in the average intensity of the near field compared to that in the far field for wavelengths ranging from 3.5 μm to 5.5 μm suggests the presence of confined modes due to the complex geometry of the metamaterial. Although these localized dark modes are electromagnetic eigenmodes of the metamaterial, they are weakly coupled to the normally incident far-field optical excitation used in the reflectance calculation because of momentum mismatch. Under thermal excitation, however, spatially distributed fluctuating dipoles can couple to these localized high-k modes, allowing them to influence the near-field thermal radiation spectrum and modify the apparent emissive response at specific wavelengths.

Dark modes are essential in understanding the complete resonant landscape of a metamaterial, including its thermal radiative properties, and are characterized by their subwavelength confinement. Considering the strong reactivity of the electromagnetic field linked to the dark modes and their high density of states, these modes are naturally narrowband. Proximal to the structure the enhancement intensity surpasses that of the bright modes, making the dark modes more prominent in the near field. For instance, at a wavelength of 4.35 μm, the intensity is lower than that at 4.8 μm at a distance of 500 nm from the metasurface. Further, the intensity enhancement under thermal excitation is approximately 38-fold in the near-field compared to the far-field; also, this dark mode is approximately 6 times greater in intensity than that in the bright mode which is situated at the center of Brillouin zone (Fig. 1c). As the distance from the metasurfae increases, the intensity of the dark mode diminishes while the bright mode intensity becomes dominant.

The contribution of dark/high-k modes can be estimated by comparing the full thermally excited near-field energy density with the far-field reference. At the thermal peak, the near-field enhancement reaches approximately 38 relatives to the far-field monitor. Thus, the excess near-field fraction is (38−1)/38×100% =97.4%, indicating that the local energy density is dominated by non-propagating or weakly radiative near-field channels. A second comparison can be made with the optically excited response, which reaches only approximately 6-fold enhancement. The peak-response difference (38−6)/38×100% = 84.2% indicates that most of the thermally excited near-field enhancement is absent from the bright-mode response accessible under normal-incidence optical excitation.

To ensure the numerical stability and convergence of the reported 38-fold enhancement, we conducted a rigorous grid convergence study by evaluating the system at three levels of spatial resolution: a baseline (Δ=1 nm), a coarser grid (1.5Δ), and a finer grid (0.5Δ). The convergence test for the near-field enhancement has been conducted, and the corresponding graph is included in the Supporting Information as Fig. S1. The peak enhancement factor exhibited a relative variation of less than 2% between the baseline and the refined mesh, confirming that the high-spatial-frequency components and evanescent modes are fully resolved.

Our meshing strategy specifically prioritizes the critical metal-dielectric (Ni-SiO₂) interfaces through localized mesh refinement. We utilized a spatial grid size of 1 nm, which is significantly smaller than the skin depth (δ) of nickel [36].

The foregoing discussion primarily examined differences in thermal radiation intensities between near-field and far-field. It is equally important to assess the radiative energy output from the heat source for a metasurface and a conventional blackbody at identical temperatures. Planck’s law defines the theoretical maximum thermal emission for a blackbody in the far field constrained by the density of free space states that facilitate the transmission of thermal energy. However, in the near field of a metasurface, additional modes can elevate the density of states above those in free space thus enabling radiative energy that can exceed blackbody radiation levels over spectral regions.

We examine the radiation intensities of the metasurface and conventional blackbody at 550 K considering that this temperature is commonly used for CO₂ mediated catalysis [37]. Applying Kirchhoff’s law, the emissivity is determined via modelling of far-field optical excitation which facilitates computation of the spectral radiation intensity in the near field. Matching of the far field intensity from the optical and thermal excitation provides o the scaling factors as a function of frequency at a given temperature. Accordingly, the near field and the far field radiation intensities due to thermal excitation are obtained. The average intensity of the near field at 2 μm from the outermost metasurface, the average intensity of the far field at 80 μm, and the far-field radiation level for a conventional blackbody are shown in Fig. 1j. We note that the electromagnetic field decays rapidly with distance from the metal surface with significant fluctuations in intensity, as observed in Fig. 1e. Noting that the 'far field' is defined as the region which is at least one wavelength in distance away from the source, beyond which the power decays proportionally with the square of the distance from the source [38], we observe relatively more uniform and slowly fluctuating fields (Fig. 1g). In the near field, the coherent nature of the modes leads to constructive and destructive phase-interference which gives to spectral enhancement of the field density. In this near-field scenario, the modes are tightly confined and effectively the fields do not propagate to the far field, thus dramatically influencing the electromagnetic environment at the nanoscale [39]. In Fig. 1k, the diagram illustrates the intensity of blackbody radiation at 550 K (indicated by the red curve) alongside the intensity of far-field radiation emanating from the heat source (represented by the green curve), which arises from the interaction between emissivity and blackbody radiation intensity. Considering that emissivity is capped at one for any frequency, the intensity of the far-field radiation consistently falls below that of the blackbody radiation across the full spectrum of frequencies. The observed near-field enhancement does not imply a violation of thermodynamic limits. The increased energy density arises from evanescent and non-propagating electromagnetic modes that decay exponentially with distance and do not contribute to far-field radiative power beyond the blackbody limit [40].

Equations (1)−(6) are written in the Gaussian electromagnetic-unit convention commonly used in the FDT-FDTD thermal-emission formulation. Then, all numerical simulation parameters, geometric dimensions, optical powers, thermal properties, heat-transfer quantities, and plotted results are reported in SI units. The transition to the heat-transfer model is therefore made explicitly by using SI quantities for absorbed energy, volumetric heat generation, thermal conductivity, heat capacity, density, and interfacial thermal conductance.

4 Heat transfer model

For a plane wave at normal incidence illuminating a metasurface – comprised of a 2D array of meta-atoms, the thermal generation associated with the absorbed light is given by the power density Q(r,t) [41]:

(7)Q\left(r,t\right)=\frac{1}{2}\text{Re}\left[J^{\ast }\left(r,t\right)\cdot E\left(r,t\right)\right],

where J^*(r,t) is the local induced current density due to the electric field E(r,t) = Re[E_ω(r). e^−iωt] inside the meta-atoms, which is converted into heat by Joule heating. This light-induced heat generation in the meta-atoms gives rise to a local temperature distribution T(r,t) [42]. The surface temperature distribution is analysed using the heat transfer equation [43]:

(8)C\rho \frac{\mathrm{\partial }}{\mathrm{\partial }t}T(r,t)-k\mathrm{\nabla }\cdot \left(\mathrm{\nabla }T(r,t)\right)=Q(r,t),

where T represents the temperature of the meta-atoms, a function of both space and time, C, ρ, and k are the spatially dependent specific heat capacity, density, and thermal conductivity of the material, respectively. The two components on the left side of the equation describe the temporal and spatial variation of thermal energy per unit volume, respectively. In solving the heat transfer equations, the boundary condition at the junction between the structure and air is defined as follows [44]:

(9)\overrightarrow{n}\cdot \left(-k\mathrm{\nabla }T\right)=h\left(T-T_{\text{r}}\right)+\epsilon \sigma \left(T^4-{T_{\text{r}}}^4\right),

where the term on the left indicates the outward conductive heat flux, while the first term on the right corresponds to the convective heat flux into the surrounding air and the second to the radiative heat flux. The normal to the surface is denoted by the vector \overrightarrow{n}, h is the coefficient for convective heat transfer of air, ε is the emissivity of the surface, σ is the Stefan−Boltzmann constant and T_r is the ambient room temperature. It is important to emphasize that as the size of the material approaches the nanoscale, the size becomes comparable to the grain size and thus surface effects become prominent. The movement of lattice phonons is notably restricted compared to that in larger specimen sizes, causing them to reflect more frequently. This reflection hinders heat transfer, significantly influencing thermal conductivity considering that the mean free path (MFP) of lattice phonons is comparable to the dimensions of the structure. The dependence of thermal conductivity on phonon mean free path and thus size at the nanoscale for different materials is presented in references [45−47]. Accordingly, the thermal properties of materials of interest – for the present study – as a function of geometric dimensions are given in Table 1 [48,49]. A key aspect of thermal transport in our metamaterial structure is the presence of a high density of interfaces and grain boundaries. These interfaces cause additional phonon scattering which in general can consider that they produce a temperature discontinuity ΔT at the interface, which is characterized by an interfacial (Kapitza) conductance, G_k [50]:

(10)\mathit{J}={\mathit{G}}_{\mathit{k}}\mathrm{\Delta }\mathit{T},

where J is the heat current. This is the analog of Fourier’s law. The inverse of the conductance, R = 1/G_k, is the Kapitza resistance. Originally, the term “Kapitza resistance” was applied only to solid/liquid interfaces at cryogenic temperatures, but it is now widely used to apply to solid/solid interfaces as well. For the present study, the experimentally determined Kapitza conductance for the nickel and SiO₂ interface is 28.2×10⁶ W/(K⋅m²) [35].

The metamaterial structure is stimulated with normally incident TM polarized light (where the electric field is along the x-axis) spanning the wavelength range of 3.5 μm to 5 μm. The absorption profile of the structure shown in Fig. 2a exhibits near-perfect absorption at 4.35 μm, highlighting the structure’s effectiveness in energy capture. Figure 2b illustrates the volumetric heat power density at this specific wavelength. Heat generation occurs principally at the metal-dielectric interface due to dissipative electron scattering interactions. Further, at the resonant wavelength of 4.35 μm the standing wave indicates a concentration of the magnetic field (H_y) within the dielectric (Fig. 2c), arising as a result of the circling electric displacement field in the adjacent nickel layers. This illustrates that the energy of the incident light is predominantly concentrated in the sub-wavelength layers of the meta-atoms [51].

The distributions of the light-induced heat dissipation, that is the heat sources, under normal incidence at 4.35 μm (resonant) and 4.8 μm (first resonant) are presented in Table 2. The results reveal that at the resonant wavelength the predominant heat source is concentrated in the lower most nickel layer of the meta-atom. For completeness, we show the distribution of the heat sources within the four distinct nickel resonators comprising the meta-atom. To examine the thermal response of the structure, the sample is subjected to a 20 μm diameter beam of incident light with a Gaussian distribution. The light fluence from a single pulse targeting the specimen is described as follows:

(11)F_{\left(r\right)}=\frac{2P}{\text{π}{w}^2{f}_r}\exp \left(-\frac{2r^2}{w^2}\right).

The full optical power of the incident light impinging the sample is P = 1.7 mW. The pulsed light source has a repetition rate of 25 kHz and a pulse width of 6 ns. At the focal point of the Gaussian beam, the fluence in each pulse is 0.043 J/cm². The optical energy contained within one unit cell is defined by the expression E_O(r) = P²\ast F_{\left(r\right)} , thus the optical energy within the meta-atom is E_O = 4.14 mJ at beam center. The thermal energy absorbed by an individual unit cell is calculated as E_t(r) = A_avgE_O(r), where A_avg is the average absorption coefficient of the structure – which is 0.40 for the range of wavelengths considered. This coefficient is calculated from the integral overlap of the power density spectrum of the light source and the metamaterial’s absorption spectrum from 3.5 μm to 5 μm, as illustrated in Fig. 2a. At the beam center, the thermal energy E_t(0) is therefore 1.65 mJ. This meta-atom light-induced heat source is modeled as a Gaussian pulsed heat source as follows:

(12)Q(r,t)=\frac{E_{\text{t}}\left(r\right)}{\text{Δ}V\sqrt{\text{π}}\tau }\exp \left(-\frac{{\left(t-t_0\right)}^2}{{\tau }^2}\right),

where ΔV represents the volume of the heat source, t₀ = 3 ns denotes the time delay of the pulse peak, and τ = 1.5 ns is the time constant of the light pulse. Figures 2d and 2e show the transient temperature profile of the proposed metamaterial structure at two distinct moments. Given the initial temperature of 293 K, the surface temperature of the meta-atom reaches 420 K within 8 ns of laser irradiation. At this point, the temperature distribution is observed to be predominantly within the nickel layers, a consequence of significant Joule heating. The temperature profiles of all the nickel layers as a function of time are shown in Fig. 2f. The lower thermal conductivity of the SiO₂ layers, situated between the Ni layers, clearly serve as thermal barriers, impeding the vertical heat transfer from the upper Ni layers to the bottom Ni layer. The first Ni layer, considering its direct exposure to light with the thermally insulating SiO₂ layer immediately below it, leads to rapid heating with a peak temperature of 1039 K in just 3.9 ns (noting that the structure is at room temperature initially). Although the bottom layer absorbs more energy, its proximity to a thick SiO₂ substrate which serves as a thermal barrier, moderates its temperature increase. Direct laser exposure of the top Ni layers, despite their lower energy absorption under optical excitation, has a significant bearing on the observed higher initial temperature spikes considering fewer immediate heat dissipation pathways compared to that of the bottom Ni layer. Over time, such as at 60 ns, the top layers of the square resonators, exhibit a slower temperature decrease compared to the bottom layers. This sustained higher temperature in the top Ni/SiO₂ layers relative to the bottom Ni layer during laser heating can be attributed to a combination of direct laser exposure, lower thermal mass, reduced heat dissipation due to insulation by SiO₂, and potentially more effective radiative cooling from the exposed surfaces [52]. Figure 2f shows the naturally relaxing temperature profile upon application of a single 6 ns laser pulse at 1.7 mW to the metamaterial unit cell. During this period, the third (from the top) nickel layer exhibits a higher peak temperature than the second nickel layer, consistent with the data in Table 2. These observations align with the thermal properties of the materials and the metamaterial’s structural configuration. Understanding these dynamics is crucial for accurately predicting and controlling the thermal behavior of such complex structures under various external conditions.

Also, it should be noted that the absorption profile in Table 2 is highly influenced by the Ni ground plane, which functions as a back reflector and forms a micro-cavity together with the SiO₂ spacer and the patterned multilayer meta-atom. This metal-insulator-metal-like configuration suppresses transmission and produces multiple reflections between the resonant meta-atom and the metallic ground plane. At the resonant wavelength of 4.35 μm, the cavity phase condition and impedance matching enhance the field residence within the spacer/meta-atom region, leading to stronger dissipative absorption in the Ni layers. The larger heat-generation fraction in the lower Ni layer is therefore consistent with cavity-assisted field localization near the spacer/back-reflector interface. In contrast, the response at 4.8 μm is associated more strongly with the Brillouin-zone-edge near-field mode, for which a larger fraction of the electromagnetic energy remains confined in non-radiative or weakly radiative channels rather than being efficiently coupled to the bright cavity mode [53].

5 Conclusion

In summary, we demonstrated thermal excitation of plasmonic metamaterials enables access to dark electromagnetic modes that dominate near-field radiation but remain inaccessible under optical excitation. Using an FDT-FDTD framework, we quantitatively showed that these modes produce near-field energy density enhancements exceeding an order of magnitude relative to the far field. Coupled photothermal simulations further reveal rapid and localized temperature modulation under pulsed excitation. These findings establish excitation mechanism as a critical design parameter for thermal photonic devices and provide a computational framework for engineering near-field thermal emission in metasurfaces. Overall, the strong correlation between spectral position, excitation mechanism, and near-field intensity points to the importance of precise spatial and temporal control in the design of energy-harvesting or conversion platforms. The insights provided here not only contribute to the fundamental understanding of near-field thermal radiation but also pave the way for practical implementations in thermophotovoltaics, nanoscale heat engines, and active thermal modulation in metasurface-based devices.

6 Appendixes

6.1 A.1 Analytical local density of state validation of the FD-FDTD model

To demonstrate that the reported 38-fold enhancement is physically robust and free of numerical discretization errors, we have added a comprehensive convergence analysis. To validate our chosen grid resolution and time steps, we simulated the near-field thermal emission of a flat Al₂O₃ plate using identical discretization parameters. As shown in Fig. A1, our numerical results show excellent agreement with exact analytical calculations, confirming our discretization scheme is highly accurate. In the simulation of FDTD computational domain, we include a 50 nm thick (z direction) aluminum plate with size of 500 nm by 500 nm (x−y plane). Boundary condition to represent the semi-infinite aluminum plate applied in x−y direction and top and bottom are terminated with PML (perfectly matched layer) boundary. In the following, we test the simulation on a simple case with known theory. We employ the analytical theory to use the LDOS (local density of state) expression above a planar metal-vacuum interface. At the beginning, we look for the relation between electromagnetic energy density and LDOS at equilibrium temperature T, to convert LDOS to energy density as following [54]:

U\left(z,\omega \right)=\rho \left(z,\omega \right)\hspace{0.17em}\frac{\hslash \omega }{\exp \left(\hslash \omega /(k_{\text{B}}T)\right)-1},

where z denotes the distance from the surface, \omega is the angular frequency, \hslash is the reduced Planck constant, k_{\text{B}} is the Boltzmann constant, and T represents the absolute temperature. The quantity \rho (z,\omega ) corresponds to the local density of electromagnetic states (LDOS) at position z and frequency \omega, showing the local electromagnetic energy density as the product of the LDOS and the Bose–Einstein energy distribution factor. Then we write the vacuum LDOS, {\rho }_{\text{v}}\left(\omega \right)={\rho }_{\text{v}}\left(r,\omega \right)={\displaystyle \frac{{\omega }^2}{{\text{π}}^2{c}^3}}, where c is speed of light which is the far-field reference level used for normalization. The LDOS above a semi-infinite flat surface depends on parameters including distance from surface z, angular frequency \omega , in-plane wavevector integration variable, Fresnel reflection coefficients r_{12}^{\text{s}} and r_{12}^{\text{p}}, propagating and evanescent contributions as appears in the following expression above a plane interface:

\begin{array}{rl}\rho \left(z,\omega \right)=& \frac{{\rho }_{\text{v}}\left(\omega \right)}{2}\{{\int }_0^1\frac{\kappa \text{d}\kappa }{p}[2\\ & +{\kappa }^2\text{Re}(r_{12}^{\text{s}}{\text{e}}^{\text{2i}p\omega z/c}+r_{12}^{\text{p}}{\text{e}}^{\text{2i}p\omega z/c})]\\ & +{\int }_1^{\mathrm{\infty }}\frac{{\kappa }^3\text{d}\kappa }{\left|p\right|}\left[\text{Im}\left(r_{12}^{\text{s}}\right)\text{+Im}\left(r_{12}^{\text{p}}\right)\right]{\text{e}}^{-2\left|p\right|\omega z/c}\}.\end{array}

This expression is actually a summation over all possible plane waves with wave number k= ω/c(k,p), r_{12}^{\text{s}}, and r_{12}^{\text{p}} are the Fresnel reflection factors between media 1 and 2 in s and p polarizations, respectively, for a parallel wave vector. Also, p is being defined as:

p=\sqrt{1-{\kappa }^2},\mathrm{f}\mathrm{o}\mathrm{r}\kappa \le 1,

p\text{= i}\sqrt{1-{\kappa }^2},\mathrm{f}\mathrm{o}\mathrm{r}\kappa >1.

Here, κ ≤ 1 is propagating and for κ >1 evanescent waves. Now we can measure the enhancement factor (EF), normalize by far field so near-field enhancement factor is being defined as [55]:

EF\left(z,\omega \right)=\frac{U\left(z,\omega \right)}{U\left(z\to \mathrm{\infty },\omega \right)}=\frac{\rho \left(z,\omega \right)}{{\rho }_{\text{v}}\left(\omega \right)}.

To get numerical optical constants Drude material model parameters for aluminum are expressed as following:

\epsilon \left(\omega \right)=1-\frac{{\omega }_{\text{p}}^2}{\omega \left(\omega +\text{i}\gamma \right)},

and then uses for aluminum: {\omega }_{\text{p}}=15.3\text{eV} and \gamma =0.598\text{eV}.

6.2 A.2 Compare with FDTD

We simulate a finite aluminum slab in a domain designed to mimic a semi-infinite plate, then extract the energy density along z at \lambda =600 THz, treat the value far from the surface as the far-field reference, and compare the simulated EF with the analytical EF, as shown in Fig. A2, if you feed those random currents into Maxwell equations, you get electromagnetic fields everywhere. Finally Fig. A1 shows z distance from the metal surface and enhancement factor of electromagnetic energy density relative to far field. As we can see far away from the surface, EF\to 1, close to the surface EF > 1 and the enhancement decays as you move away from the interface. Different values for EF mean the interface supports additional near-field, especially evanescent, electromagnetic states that do not exist in free space. To put all together, high near the interface → decays with distance → reaches 1 far away.

The relative error was calculated at each distance as:

\text{Relative error}\left(z_i\right)=\frac{\left|E{F}_{\text{FDTD}}\left(z_i\right)-E{F}_{\text{Analytical}}\left(z_i\right)\right|}{E{F}_{\text{Analytical}}\left(z_i\right)}\times 100\mathrm{\%}.

Using the near-field enhancement factor calculated via FDT-FDTD and analytical LDOS, the average relative error for distances z ≥ 0 is approximately 1.83%, with a maximum deviation of approximately 5.2%. The largest discrepancy occurs only at the first few nanometers from the interface, where the LDOS varies extremely rapidly and is most sensitive. The agreement over the remaining near-field range confirms that the stochastic source implementation, PML treatment, time-stepping, and spatial discretization reproduce the analytical LDOS trend before the method is applied to the patterned metamaterial.

References

[1] Rytov, S.M.: Theory of electric fluctuations and thermal radiation. AFCRCTR (1959)

[2] Xuan , Y. : An overview of micro/nanoscaled thermal radiation and its applications. Photon. Nanostructures 12(2), 93–113(2014)

[3] Shinohara , N. : Trends in wireless power transfer: WPT technology for energy harvesting, mllimeter-wave/THz rectennas, MIMO-WPT, and advances in near-field WPT applications. IEEE Microw. Mag 1(22), 46–59(2020)

[4] Chen , F. , Liu , X. , Tian , Y. , Zheng , Y. : Dynamic tuning of near‐field radiative thermal rectification. Adv. Eng. Mater 23(2), 2000825(2021)

[5] Greffet , J.J. , Carminati , R. , Joulain , K. , Mulet , J.P. , Mainguy , S. , Chen , Y. : Coherent emission of light by thermal sources. Nature 416(6876), 61–64(2002)

[6] Raman , A.P. , Anoma , M.A. , Zhu , L. , Rephaeli , E. , Fan , S. : Passive radiative cooling below ambient air temperature under direct sunlight. Nature 515(7528), 540–544(2014)

[7] Song , J. , Han , J. , Choi , M. , Lee , B.J. : Modeling and experiments of near-field thermophotovoltaic conversion: a review. Sol. Energy Mater. Sol. Cells 238, 111556(2022)

[8] Zhao , B. , Guizal , B. , Zhang , Z.M. , Fan , S. , Antezza , M. : Near-field heat transfer between graphene/hBN multilayers. Phys. Rev. B 95(24), 245437(2017)

[9] Mehrabi , S. , Rezaei , M.H. , Zarifkar , A. : Ultra-broadband solar absorber based on multi-layer TiN/TiO₂ structure with near-unity absorption. J. Opt. Soc. Am. B 36(9), 2602–2609(2019)

[10] Noreen , S. , Rehman , A.U. , Zubair , M. , Abbasi , Q.H. , Mehmoo , M.Q.: Design, fabrication and analysis of all-metal 3D metamaterial based electromagnetic absorber. In: Proceedings of 2nd International Conference on Microwave, Antennas & Circuits, pp. 1–5 (2025)

[11] Mehrabi , S. , Rezaei , M.H. , Rastegari , M.R. : High-efficient plasmonic solar absorber and thermal emitter from ultraviolet to near-infrared region. Opt. Laser Technol 143, 107323(2021)

[12] Park , K. , Zhang , Z. : Fundamentals and applications of near-field radiative energy transfer. Front. Heat Mass Transf 4(1), 1–26(2013)

[13] Johnson , C.: Mathematical Physics of Blackbody Radiation. Icarus iDucation (2012)

[14] Reiser , A. , Schächter , L. : Geometric effects on blackbody radiation. Phys. Rev. A 87(3), 033801(2013)

[15] Maslovski , S.I. , Simovski , C.R. , Tretyakov , S.A. : Overcoming black body radiation limit in free space: metamaterial superemitter. New J. Phys 18(1), 013034(2016)

[16] Luo , X. : Subwavelength artificial structures: opening a new era for engineering optics. Adv. Mater 31(4), 1804680(2019)

[17] Rytov , S.M. , Kravtsov , Y.A. , Tatarskii , V.I.: Principles of Statistical Radiophysics: Elements of random fields. Berlin: Springer (1989)

[18] Prost , J. , Joanny , J.F. , Parrondo , J.M.R. : Generalized fluctuation-dissipation theorem for steady-state systems. Phys. Rev. Lett 103(9), 090601(2009)

[19] Polimeridis , A.G. , Reid , M.T.H. , Jin , W. , Johnson , S.G. , White , J.K. , Rodriguez , A.W. : Fluctuating volume-current formulation of electromagnetic fluctuations in inhomogeneous media: incandescence and luminescence in arbitrary geometries. Phys. Rev. B Condens. Matter Mater. Phys 92(13), 134202(2015)

[20] Walter , L.P. , Tervo , E.J. , Francoeur , M. : Near-field radiative heat transfer between irregularly shaped dielectric particles modeled with the discrete system Green’s function method. Phys. Rev. B 106(19), 195417(2022)

[21] Didari , A. , Pinar Mengüç, M. : A design tool for direct and non-stochastic calculations of near-field radiative transfer in complex structures: the NF-RT-FDTD algorithm. J. Quant. Spectrosc. Radiat. Transf 197, 95–105(2017)

[22] Li , Z. , Li , J. , Liu , X. , Salihoglu , H. , Shen , S. : Wiener chaos expansion method for thermal radiation from inhomogeneous structures. Phys. Rev. B 104(19), 195426(2021)

[23] Biehs , S.A. , Messina , R. , Venkataram , P.S. , Rodriguez , A.W. , Cuevas , J.C. , Ben-Abdallah , P. : Near-field radiative heat transfer in many-body systems. Rev. Mod. Phys 93(2), 025009(2021)

[24] Shi , K. , Chen , Z. , Xing , Y. , Yang , J. , Xu , X. , Evans , J.S. , He , S. : Near-field radiative heat transfer modulation with an ultrahigh dynamic range through mode mismatching. Nano Lett 22(19), 7753–7760(2022)

[25] Habibi , M. , Beardo , A. , Cui , L. : Near-field thermal radiation as a probe of nanoscale hot electron and phonon transport. ACS Nano 19(6), 6033–6043(2025)

[26] Odebowale , A.A. , Berhe , A. , Ogundare , R.T. , Abdo , S. , Abdulghani , A. , Hattori , H.T. , Miroshnichenko , A.E. : Advances in radiative heat transfer: bridging far‐field fundamentals and emerging near‐field innovations. Adv. Funct. Mater 35(27), 2421051(2025)

[27] Chan , D.L.C. , Soljačić , M. , Joannopoulos , J.D. : Direct calculation of thermal emission for three-dimensionally periodic photonic crystal slabs. Phys. Rev. E Stat. Nonlin. Soft Matter Phys 74(3), 036615(2006)

[28] Loomis , J.J. , Maris , H.J. : Theory of heat transfer by evanescent electromagnetic waves. Phys. Rev. B 50(24), 18517–18524(1994)

[29] Luo , C. , Narayanaswamy , A. , Chen , G. , Joannopoulos , J.D. : Thermal radiation from photonic crystals: a direct calculation. Phys. Rev. Lett 93(21), 213905(2004)

[30] Levin , M.L. , Rytov , S.M.: Theory of Equilibrium Thermal Fluctuations in Electrodynamics. Moscow: Nauka (1967)

[31] Ding , X. , Zhang , Z. , Sun , J. , Loh , J. , Ji , D. , Lu , J. , Liu , C. , Zhao , L. , Liu , W. , Zhao , J. , Tang , S. , Safari , M. , Cai , H. , Tu , W. , Kherani , N.,Hu , Z. , Ozin , G. , Zou , Z. , Wang , L. : Thermal radiative catalysis: selective dehydrogenation of ethane to ethylene by vibrationally excited carbon dioxide. Joule 7(10), 2318–2334(2023)

[32] Malitson , I.H. : Interspecimen comparison of the refractive index of fused silica. J. Opt. Soc. Am 55(10), 1205–1209(1965)

[33] Palik , E.D.: Handbook of Optical Constants of Solids. Cambridge: Academic Press (1998)

[34] Fu , P. , Chen , Z. , Yi , Z. , Cheng , S. , Ahmad , S. , Li , B. : Graphene‐based terahertz perfect absorber with broadband, wide‐angle, and dynamically tunable performance. Physica Status Solidi (RRL) 20(4), e70163(2026)

[35] Zhang , H. , Chen , Z. , Yi , Z. , Cheng , S. , Ahmad , S. , Tang , C. , Deng , J. : Three band THz perfect absorption device based on graphene metasurface structure. Physica E 181, 116542(2026)

[36] Ordal , M.A. , Bell , R.J. , Alexander , R.W. Jr, Long , L.L. , Querry , M.R. : Optical properties of Au, Ni, and Pb at submillimeter wavelengths. Appl. Opt 26(4), 744–752(1987)

[37] Li , K. , Niu , J. , Guo , B. , Liu , H. , Jin , Y. , Wang , J. , You , C. , Ran , J. : Understanding the transformation mechanism of carbon species on different Nickle facets in CO₂ reforming reaction. Mol. Catal 564, 114332(2024)

[38] Pedrotti , F.L. , Pedrotti , L.M. , Pedrotti , L.S.: Introduction to Optics. Cambrige: Cambridge University Press (2017)

[39] Kim , K. , Song , B. , Fernández-Hurtado , V. , Lee , W. , Jeong , W. , Cui , L. , Thompson , D. , Feist , J. , Homer Reid, M. , García-Vidal , F. , Cuevas , J. , Meyhofer , E. , Reddy , P. : Radiative heat transfer in the extreme near field. Nature 528(7582), 387–391(2015)

[40] Zhu, H., Ren, Y., Pan, H., Tang, G., Zhang, L., Jian-Sheng Wang, J.: Enhancing far-field thermal radiation by floquet engineering. arXiv preprint arXiv:2507.16688 (2025)

[41] Bohren , C.F. , Huffman , D.R.: Absorption and Scattering of Light by Small Particles. John Wiley & Sons (2008)

[42] Baffou , G. , Quidant , R. : Thermo‐plasmonics: using metallic nanostructures as nano‐sources of heat. Laser Photonics Rev 7(2), 171–187(2013)

[43] Christensen , M.G. , Adler-Nissen , J. : Simplified equations for transient heat transfer problems at low Fourier numbers. Appl. Therm. Eng 76, 382–390(2015)

[44] Isachenko , V.P. , Osipova , V.A. , Sukomel , A.S.: Heat Transfer. Nirali Prakashan (1980)

[45] Negi , A. , Rodriguez , A. , Zhang , X. , Comstock , A.H. , Yang , C. , Sun , D. , Jiang , X. , Kumah , D. , Hu , M. , Liu , J. : Thickness‐dependent thermal conductivity and phonon mean free path distribution in single‐crystalline barium titanate. Adv. Sci. (Weinh.) 10(19), 2301273(2023)

[46] Hosseini , S. , Safaei , M.R. , Goodarzi , M. , Alrashed , A.A.A.A. , Nguyen , T.K. : New temperature, interfacial shell dependent dimensionless model for thermal conductivity of nanofluids. Int. J. Heat Mass Transf 114, 207–210(2017)

[47] Jaiswal , R.L. , Pandey , B.K. : Modelling for the variation of thermal conductivity of metallic nanoparticles. Physica B 627, 413594(2022)

[48] Chien , H.C. , Yao , D.J. , Hsu , C.T. : Measurement and evaluation of the interfacial thermal resistance between a metal and a dielectric. Appl. Phys. Lett 93(23), 231910(2008)

[49] Yuan , S. , Jiang , P. : Thermal conductivity of nanoscale thin nickel films. Prog. Nat. Sci 15(10), 922–929(2005)

[50] Chernatynskiy , A. , Clarke , D.R. , Phillpot , S.R.: In Handbook of Nanoscience, Engineering, and Technology. Boca Raton: CRC Press (2018)

[51] Khuyen , B.X. , Tan , P.D. , Tung , B.S. , Hai , N.P. , Tuan , P.D. , Phong , D.X , Tung , D.K. , Anh , N.H. , Giang , H.T. , Vinh , N.P. , Tung , N.T. , Lam , V.D. , Chen , L. , Lee , Y.P. : Numerical optimization of metamaterial-enhanced infrared emitters for ultra-low power consumption. Photonics 12(6), 583(2025)

[52] Ijaz , S. , Kang , D. , Rana , A.S. , Kim , J. , Chani , M.T.S. , Zubair , M. , Abbassi , Q.H. , Mehmood , M.Q. , Rho , J. : Metasurface absorber–emitter pair-integrated high-efficiency thermophotovoltaic system. ACS Photonics 12(7), 3829–3839(2025)

[53] Zeng , N. , Chen , Z. , Yi , Z. , Cheng , S. , Ahmad , S. , Tang , C. , Gao , F. , Li , B. : Terahertz multi-band tunable refractive index sensing graphene absorber based on surface plasmon resonance. Physica B 734, 418608(2026)

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

[55] Blaber , M.G. , Arnold , M.D. , Ford , M.J. : Search for the ideal plasmonic nanoshell: the effects of surface scattering and alternatives to gold and silver. J. Phys. Chem. C Nanomater. Interfaces 113(8), 3041–3045(2009)

RIGHTS & PERMISSIONS

The author(s)