Efficient molecular emitter

We demonstrate extraordinarily spectrally selective narrowband mid-infrared radiation absorbance and thermal emittance with resonant peak FWHM ≤ 124 nm at λ = 5.73 μm, corresponding to a Q-factor of ∼ 92.3. This was achieved by harnessing mode coupling between a plasmonic metal-insulator-metal (MIM) metasurface and molecular vibrational mode resonances, with coupling constants ranging from η ∼ 3.9% to 6.6%. In addition, thermal radiation emissivity is in close accordance to the metamaterial absorbance spectrum, as described by Kirchhoff’s law of thermal radiation, and furthermore, emission was not angle dependent, unlike that exhibited by grating-based emitters. The experimentally investigated MIM structures remained stable up to a 250◦C heating temperature. MIM metamatrials with strong and spectrally tailored vibrational coupling behaviors represent a new paradigm in photo-thermal energy conversion. The experimentally observed pronounced resonant coupling behaviour was well described by finite-difference time-domain simulations of the plasmonic structure, where molecular vibration contributions were modeled using the Lorenz oscillator approximation.


Introduction
Harnessing and controlling radiative thermal energy continues to be a pressing scientific and technological challenge 1-3 , touching upon applications as diverse as molecular spectroscopy and bio-medical sensing, as well as enhancing energy conversion efficacy of solar cells. A blackbody remains the quintessential example of a thermal-to-radiative energy converter, however, even in its idealized state its utility is limited by a broad spectral energy distribution, dictated solely by its temperature. Yet, exerting greater control over radiated energy holds considerable unfulfilled promise 4,5 . For example, spectrally confining thermal radiation to a specific atmospheric transmission window has potential for facilitating the Earth's cooling, similar to how thermo-photovoltaic devices can benefit from a greater conversion efficiency when selective emitters are employed 6 . In contrast to such relatively broadband cooling-directed applications, chemical sensors require narrowband radiation tuned to a particular molecular vibrational mode to enable high selectivity and sensitivity detection. Sensors at the telecommunication spectral window at near-IR 7 are especially appealing due to capabilities to integrate them into CMOS (complementary metaloxide-semiconductor) platform. Hence, development of a thermal radiation emitter with spectral behaviour controlled more rigorously than lamp sources, and more affordable and ambient temperature compatible than quantum cascade lasers, is in considerable demand 8 .
It is expected that thermal radiation from a thick polymer film could be used for narrowband thermal emission purposes. It is straightforward to experimentally show that a 3 µm thick polyimide film can absorb almost 100% of incident thermal radiation at its primary molecular absorption band. However, the emission of radiation was not efficient at the the same wavelength, primarily due to strong re-absorption. The search for wavelength-selective thermal emitters is informed largely by Kirchhoff's law, according to which thermal radiation energy output is directly proportional to absorptivity in a particular spectral range and at a given temperature. One of the most direct ways of exploiting this effect relies on selection of natural materials or synthetic composites with appropriate electronic and phononic band structures, however, the best attainable performance and the degree of control that can be exerted in this case is fundamentally limited by physical chemistry considerations. A complimentary method of tailoring absorbance and, by extension, emittance behavior relies on wavelength scale structure patterning to manipulate electromagnetic mode propagation and confinement. Examples of such structures harnessed for thermal emitters include photonic crystals, and, in particular, different types of metamaterials -artificial wavelength-scale periodic arrangements of metal structures.
Of special note are various electromagnetic energy absober metasurfaces, which represent an intriguing approach to the creation of narrowband thermal radiation sources 9 . For instance, grating structures that allow for ultra-narrow band thermal emission out-coupled at a specific angle have been demonstrated, however, are subject to limitations due to that inherent angular dependence of resonance. On the other hand, plasmonic absorber metasurfaces likewise support spectrally narrow resonance modes, with widths only around 200-600 nm in mid-infrared (MIR) wavelength range, with the added benefit of their angular pattern of radiation being almost angle-independent 8 . Recent trends in narrow-band thermal radiation emitters show Q-factors ranging between 60 and 200 10-14 with specific materials/structures. A less than 100 cm −1 bandwidth is required for molecular detection, which is readily satisfied by molecular absorption enhanced using metasurface configurations [10][11][12][13][14] . The unique merit of using a molecular vibrational mode is the possibility for harnessing its spectral signature for targeted identification of the same molecular functional group, e.g. -OH or -NH and the like.
In a thermal emitter context a metasurface geometry fulfils advanced functions of free-space impedance matching as well as spatial and temporal electromagnetic energy concentration and confinement, however, thermal energy is both sourced from or absorbed by the vibrational degrees of freedom inherent to molecular structures of the constituent materials themselves. At a basic level molecular vibrations in a polymer film can be harnessed as an emitter, provided its thickness is sufficiently high 15 . For instance, a 300-µm-thick polyimide film can exhibit near 100% absorption and radiation emission. However, the thermal conductivity of a thick isotropic polymer film is low, which hampers rapid energy dissipation required for acceptable absorber/emitter performance. An alternative strategy utilizes plasmonic mode enhancement for attaining high photo-thermal absorbance and emission in lower thickness molecular films. Hence, in pursuit of exceedingly narrowband thermal energy radiation emitters, it can be expected that best results are attainable when narrowband metamaterial electromagnetic transducers are coupled to similarly narrowband material absorbers. In particular, the thermal radiation from molecular vibrational modes would have a suitably narrow spectral footprint to enable the targeted detection of analogous molecular vibrational signatures in a sensing device. Furthermore, thermal radiation extracted from molecular vibrations was explored as an approach for large-scale heat transfer between buildings and the surrounding atmosphere. For these types of applications thermal radiation out-coupling efficiency becomes the primary driver towards practical implementation. However, narrow-band MIR absorbers and emitters with suitably high spectral finesse and radiation out-coupling efficiency are yet to be experimentally demonstrated.
We report large emissivity and extraordinarily narrowband thermal radiation output from polyimide molecular vibrational modes, coupled to localized plasmon resonances supported by metal-insulator-metal (MIM) metasurfaces. The efficient coupling is facilitated by the strong localization of resonant radiation field in the thin-film polyimide insulator layer, as illustrated in Fig. 1. Coupling strength dependent interactions between MIM resonant modes and the modal structure of polyimide vibrational spectrum were modeled theoretically and systematically studied experimentally by varying top layer metal disk diameters [16][17][18] . Furthermore, metasurfaces were fabricated using highly scaleable and cost effective optical lithography methods, and the polyimide absorber film was selected for its compatibility with both semiconductor CMOS process and thermal treatments up to ∼ 300 • C. The investigated structures exhibit a remarkable agreement between theoretical and experimental findings, and provide a way to devise highly tailored planar MIR radiation sources which for their spectral signatures can draw upon a rich variety of molecular compounds.

Results and discussion
Plasmon coupling with vibrational modes of polyimide The MIM metasurface with a polyimide insulator layer investigated here is schematically illustrated in Figure 1, along with a scanning electron microscopy (SEM) image of a fabricated structure. A key prerequisite for effective coupling between molecular vibrations and plasmonic modes is ensuring both a spectral and spatial overlap between them. A layered metamaterial arrangement used here, in which polyimide is wedged between a 50 nm thickness top nanodisk structure and an underlying 200 nm thickness metal reflector, is exceptionally well-suited in this respect, as it ensures that the resonant optical field is primarily confined in the insulator, as illustrated by the finite-difference time-domain (FDTD) simulation derived plots in Fig. 1(a). Conversely, the spectral overlap can be tuned by altering the scale of the metamaterial pattern, where nanodisk diameter D where changed in the 500 nm to 1600 nm range and period P was fixed at 1.5 × D.
Experimental reflectance spectra of MIM samples with increasing nanodisk diameters, starting from D ≥ 900 nm, are given in Fig. 2(a), and reveal the absorbance of said structures. Here all the probed metamaterial patterns were fabricated on the same Au coated Si substrate, which allowed to ensure polyimide layer thickness and composition consistency. For D = 900 nm nanodisks, when the metamaterial resonance is mostly decoupled from polyimide vibrational modes, the absorption peak is situated in the λ = 5 µm wavelength region and possesses a relatively narrow 600 nm FWHM width. Likely the main factor limiting the metasurface quality factor is the scattering due to nanoscale surface roughness, derived from a conformal transfer of the the polyimide layer roughness to the overlaying gold disks during their subsequent deposition. On the other hand, the PMDA-ODA polyimide insulator molecular vibration spectrum is dominated by four sharp and strong absorption peaks at 1736, 1509, 1388, 1253 cm −1 -respectively attributable to carbonyl stretching, aromatic amine ring-breathing, C-N stretch, and imide deformation modes 19 -highlighted in Fig. 2(b). Polyimide has several bands whose orientation can be used for determining molecular alignment: 1780, 1745, 1520 and 1400 cm −1 , which are respectively associated with C=O symmetric (chain axis), C=O asymmetric (ring tilt), C 6 H 5 aromatic rings (phenyl) tangential vibrations, and C-N-C axial stretch [20][21][22] . The parallel (to film surface) bands are situated at 1491 cm −1 (ν C-C; skeletal vibration of 1,2,4-trisubstituted benzene ring in diamine moieties) and 1501 cm −1 (ν C-C; skeletal vibration of 1,4-disubstituted benzene in diamine moieties), while the perpendicular band is at 1724 cm −1 (ν C=O out of phase asymmetric stretching of carbonyl bonds in imide rings) [20][21][22] .
As the nanodisk diameter is increased, the MIM metasurface resonance is redshifted towards the spectral region of the different molecular vibrational modes and gradually begins to overlap with them. This is accompanied by an absorbance enhancement from vibrational modes by up to a factor of 50. When coupling comes into effect the polyimide molecular vibrational modes undergo hybridization with the plasmon resonance of MIM metasurface, leading to clear spectral shifts and splittings. Initially, when a polyimide peak is situated in close vicinity to the MIM mode, it gets slightly redshifted from its respective intrinsic molecular vibrational mode position. However, when plasmonic and vibrational modes overlap particularly strongly, peaks instead undergo a blueshift. These two stages are respectively termed the weak coupling and the strong coupling regime, as categorised by A. F. Kockum,et. al. 2 . The coupling strength is classified as strong or weak based on whether the light-matter coupling strength g is sufficiently large to exceed the losses in the system or not. The weak coupling regime is associated with the Purcell effect, responsible for a variety of photonic phenomena such as fluorescence enhancement, surface-enhanced Raman scattering (SERS) and surface-enhanced infrared absorption (SEIRA). On the other hand, the signature characteristic of strong coupling are vacuum Rabi oscillations that can manifest as an energy splitting in the dispersion relation [23][24][25][26][27] .
To further elucidate the spectral peculiarities of the vibrational mode coupling to the MIM metasurface structure, FDTD simulations were performed assuming a four-level Lorenz model for insulator layer permittivity, summarized in Fig. 2(c); details of the simulation method and parameters are provided in the Methods section. As plotted in Fig. 2(d,e), polyimide bands numerically synthesised for simulations exhibited linewidth and relative intensity values comparable to the experimental ones. Similarly, by adjusting fitting parameters the FDTD simulations qualitatively well reproduce the experimental absorbance measurement results. However, it must be noted that the experimental absorption peaks exhibit slightly narrower linewidths than numerically calculated ones, likely attributable to quantization of the molecular vibrational modes not accounted for by the classical Lorenz model.

Narrowband MIR radiation emitters
Spectral plots of thermal radiation emitted by an externally heated polyimide MIM sample with D = 1500 nm diameter discs are given in Fig. 3. Emitted thermal radiation spectra, shown in Fig. 3(a), clearly mirror the reflectance, and by extension absorbance, exhibited by the metasurface. This shows that the plasmon-coupled molecular vibrational modes obey the Kirchhoff's law of thermal radiation. At the condition of strong coupling, and due to a relatively large bandwidth of the plasmon resonance, the absorption (emission) spectrum also became broader than that of the intrinsic vibrational mode absorption in polyimide. However, the 1736 cm −1 (∼ 5.76 µm) mode, subjected to a weak interaction regime, generates a spectrally narrow radiation line with an approximately 90 nm bandwidth. This is a marked improvement over the thermal radiation generated from MIM metasurfaces that typically exhibit a 200-600 nm bandwidth. A brief overview of reported MIM metasurface optical absorption and thermal radiation emittance performance is provided in Fig. 9 of Supporting Information.
Strong and weak coupling of a plasmonic resonance with molecular vibrational modes show considerable potential for enabling the creation of extraordinarily narrow band thermal radiation sources. The main prerequisite for that, however, is sufficient thermal stability on behalf of the metasurface structure, and particularly that of the insulator molecular compound. Polyimide is well known for its high thermal stability and heat resistance up to a 350 • C temperature. In our experiments, stable thermal radiation output was successfully measured for heating temperatures ranging from 50 to 250 • C. Some of the emitted radiation spectra are given in Fig. 3(b), where they are normalized to respective black-body spectral radiance. Above 300 • C, MIM structure became distorted due to softening and gradual decomposition of polyimide, and its plasmon resonance spectral position shifted, causing corresponding deleterious changes in thermal radiation output. However, when operated in the appropriate temperature range MIM metasurfaces exhibited a marked robustness in the angular dependence of emission upon heating, in contrast to grating-type surface radiators.
To realize even more spectrally narrow thermal radiation emitters it is important to investigate in detail the coupling between plasmon resonance and the molecular vibrational modes. Despite their disparate and complex physical origins, both the plasmon resonances supported by MIM metasurfaces as well as molecular absorption modes can be reasonably well approximated using the Lorenz oscillator model. If the two resonances have a good overlap both spatially and spectrally, the so called "strong coupling" conditions can be realised. Their characteristic feature is mode splitting into higher and lower energy branches that becomes apparent in dispersion spectra. Most of these strong coupling phenomena have been observed for J-aggregate dye molecules such as porphyrine and rhodamine, which have a high dipole moment in the aggregated form. In the same way, polymer materials are expected to have a relatively directional polymer chain which is favorable for realisation of larger dipole moments.
We have simulated a four-oscillator, as well as separate single-oscillator Lorenz models for each of the polyimide peaks interacting with the MIM plasmon mode as a function of the constituent metal disc diameters, as shown in Fig. 4. The plasmon resonance wavelength of the MIM metasurface has a linear dependence on disk diameter. However, when this plasmon resonance overlaps with molecular absorption, mode splitting becomes apparent. When the permittivity is defined by the single-oscillator Lorenz model (Fig. 4(c-f)), a conventional dispersion curve is obtained. However, when additional bands are included and neighboring Lorenz modes start to partially overlap each other, as in the experiment case of polyimide, the dispersion curve acquires a sigmoidal shape. This is a manifestation of the interaction observed as mode splitting. The normalized coupling strength η between the plasmon resonance and a molecular vibrational mode is described by 2 : where E − and E + are the energy values associated with the split bands, E 0 is the molecular vibrational mode energy, and therefore the splitting itself can be expressed as: where g 0 is the coupling strength and ω C is the angular frequency of molecular absorption. Based on experimental results, the η values for the weakly coupled modes I to IV (Fig. 5) are η I = 6.6%, η II = 4.8%, η III = 3.9%, and η IV = 4.1%, respectively. This result is in good correspondence with FDTD simulations. If the molecular absorption coefficient becomes larger, the η value is likewise increased. Such η values are typical in the MIR wavelength region [28][29][30] . However, since η is proportional to √ ω, the values obtained in the MIR spectral region tend to be 2 to 5 times smaller than those reported in the strong coupling case between organic molecules and plasmons at visible wavelengths 31,32 . Also, the absorption coefficient in polymers at IR is smaller than that of an electron system of dye molecules at visible wavelengths. Most of studies in visible wavelength region were carried out on J-aggregates of polphirine. A molar absorption coefficient of polphirine is ε = 2.6× 10 5 mol/dm −3 cm −1 , 33, 34 while its MIR absorption is in the range of ε = 10 3 -10 4 mol/dm −3 cm −1 . A quantum mechanical description of the molecular vibration mode is 35 : where ε is the molar absorption coefficient, Ψ i and Ψ j are wave functions of the molecular vibrational modes at the ground and excited states, respectively,μ is the dipole moment operator, µ 0 is the dipole moment of the equilibrium position, Q k is the small displacement, dτ is the elementary volume, ν is frequency of light, N A is Avogadro constant, h is Plank constant, c is speed of light, ε 0 is the permittivity of vacuum. The last term in Eqn. 3 is transition moment M, which is the dipole moment change of a molecular vibrational mode described as: Absorption in the molecular vibrational manifold occurs when ∂ µ/∂ Q k = 0. Therefore, if the change of molecular dipole moment becomes larger during oscillation, the absorption coefficient becomes enhanced. Polyimide main absorption band is related to C-N or C=O bond vibrations, as discussed above. A heterogeneous bonding generates a larger absorption coefficient than a CH functional group. The largest difference between UV-Vis and MIR absorption is that during UV-Vis absorption in an aggregated dye molecule (e.g., J-aggregate of polphyrine) a dipole extends over the span molecule, which makes the dielectric moment larger. Conversely, vibrational absorption at MIR wavelengths is weak due to harmonic oscillation of atoms which cause negligible interaction between functional groups, since they are independent of each other. Therefore, in a polymer chain, vibrational mode can be described as ε ∼ N × ε mon where ε mon is the monomer molar extinction, N is the number of polymerized monomers. When number N of oscillators are located in a cavity with transition moment M, the splitting strength η is proportional to √ N and M. We can conclude that one can observe splitting due to a high density of material even with a weak absorption.
To confirm the effect that coupling strength has on spectral signatures, a comparison between spectra was carried out with D = 1500 nm discs, outlined in Fig. 5. It should be noted that at least eight peaks in the absorption spectra of the vibrational modes of polyimide can be identified when Lorenzian curve fitting is applied. According to the rules associated with strong coupling, the interaction between a single plasmon resonance and one molecular absorption peak generates two spectral signatures. Furthermore, the MIM plasmon resonance has a sufficient spectral width to encompass all four vibrational peaks, hence, induces splitting of all four modes into a total of eight signatures (Fig. 5). The plasmon resonance overlaps most strongly with peaks II and III. In that case, split modes can be identified as having both upper and lower energy level branches. However, peaks I and IV, which represent a weaker overlap between the MIM plasmon resonance and molecular vibrational modes do not exhibit splitting behavior. If no interaction between the split modes is assumed, the modes in the four-oscillator Lorenz model are expected to appear at the same energy as described by the single-level Lorenz model. However, only the two peaks at the positions I and IV (Fig. 5), distinct from the crowded family of peaks mostly derived from II and III modes, match the positions predicted by their single-oscillator Lorenzian models. The other six closely spaced peaks formed from the original four-level Lorenz model have to be analysed in greater detail. Their complexity is consistent with the prominence of the interaction and mode mixing between the neighbouring split modes. In the case of strong coupling, all the peaks still have quality factors 7/16 Figure 6. Normalized spectral mapping of reflection and radiation obtained for the coupling constant η analysis from metasurfaces shown in Fig. 5.
Q ≥ 50, but initially more pronounced molecular vibrational modes become broadened due to the comparatively wide-band plasmon resonances as shown in Fig. s2.
The linewidth predicted by FDTD model was found broadened in the experiments. We used the four-polarization method to check for a presence of molecular orientation, which could explain narrowing of the bands 36 . Spectral maps at selected absorption bands showed that there was no molecular orientation in the MIM arrays with evaporated nano-thin polyimide (Fig.  s3). Interestingly, similar tendencies can be observed even using a spin coated amorphous polyimide film as shown in Fig. s4. Therefore, it could be concluded that a limitation of the transition probability of the quantized IR absorption levels is the most probable cause of the observed narrowing of the peaks when the molecular vibrations are strongly coupled to the collective plasmonic band. Also, it is expected that a higher coupling strength could be attained using oriented polyimide chain films. The resonance peaks also become shifted from their initial molecular absorption spectral positions. However, in the weak coupling regime, represented by mode I, the line width remains almost unchanged from an uncoupled molecular absorption peak. Therefore, the weak coupling condition with strong absorption enhancement has considerable potential for applications in molecular sensing where narrowband MIR emitter sources are in considerable demand. In this regard, how the strong coupling and mode splitting regime can be gainfully utilized as a narrowband source is less apparent.

Adjustment of polyimide thickness
The thickness of an insulating layer is an important parameter in tailoring the plasmonic response of a metasurface for best mode overlap with molecular absorbers and efficient coupling 37 . Among the polymer deposition methods spin coating makes it particularly easy to control thickness of polyimide compared with vacuum deposition methods. In addition to relatively thick 335 nm films, dilution of polyimide solution enabled the preparation of 45-nm-thick layers, which were be thinner even than those made by vapor deposition. The reflectance and radiation emittance spectroscopy data for such MIM metasurface embedded films are shown in Fig. 5. Spin coated thin films exhibited the same absorption peaks as compared to polyimide deposited by vacuum deposition and had similar peak splittings. Over the entire tested nano-disk diameter D range, a mirrormatching reflection and radiation spectra were observed, in close accordance to results shown in Fig, 3. Also, the peak splitting induced by coupling between plasmon resonance and molecular vibrational mode was likewise observed. To calculate peak separation precisely, a mapping was carried out for the normalized reflection (1 − R) and radiation in the spectral range from 1000 to 3000 cm −1 . The lower 45 nm thickness film exhibited a larger band splitting than the thick 335 nm polyimide film counterpart. The η value of the band-I was 9.8% for 45 nm and 7.6% for 335 nm films, respectively. The absolute value of η is larger than that for vacuum deposited polyimide. This is likely caused by slight differences of polymerization degree and/or absorption coefficient.

Conclusion
Control over extraordinarily narrow spectral bandwidth thermal radiation absorption and emission, with a quality factor of Q ≥ 90, is demonstrated by exploiting weak and strong coupling of molecular vibrational modes with a plasmon resonance supported by a MIM metamaterial. Reciprocity between absorption and emission strength and spectral line-shapes is in accordance with the Kirchhoff's law. The thermal radiation emission spectral features exhibit nearly the same widths as the molecular vibrational modes they were derived from. The uniformity of angle-resolved radiation pattern, its narrow spectral signatures and radiation intensity are all highly compatible with sensing applications. This metamaterial IR radiation source exhibits properties superior compared to alternative grating-based thermal radiation emitter surfaces. Polyimide is among the most thermally stable polymer materials, however, thermal radiance could be enhanced further if polymer materials that can maintain their properties at even higher temperatures can be utilized. Recent advances in production of high-temperature stable polymers culminated with a variety of polyimde stable up to 500 • C 38 . Most polymeric materials exhibit IR absorption in the 4 to 6 µm spectral region, located around the MIR transmittance window of air. Hence, in principle, their metal-polymer composites can be harnessed towards applications in radiative cooling. Furthermore, narrowband emission spectral fingerprint can be tuned by specifically selecting molecular compounds with particular vibrational mode sets. Alternatively, plasmon enhancement can be applied not only to organic materials but to inherently robust inorganic materials such as SiO 2 , SiC and other thermally stable materials 30,39 .

Sample fabrication
Polyimide film was formed using vapor deposition polymerization of pyromellitic dianhydride (PMDA) and 4,4'-diaminodiphenyl ether (ODA) on the 200 nm Au layer coated Si substrate 40,41 . Temperature of the source was set to 150 • C (PMDA) and 110 • C (ODA), while the substrate temperature was set to 30 • C. After 30 min of deposition, thermal imidization treatment has been carried out at 150 • C in air for 1 hour. Polyimide layer formation was confirmed using infrared reflection-absorption spectroscopy (IR-RAS). Thickness of polyimide was 60 nm.

9/16
Another method of polymer film formation involved spin coating the polyimide precursor varnish made of poly(pyromellitic dianhydride-co-4,4'-oxydianiline), amic acid solution 15.0-16.0 wt.% in N-methylpyrrolidone (NMP) (Pyre-M.L. RC-5019, Sigma-Ardrich). After spin coating on a Si substrate, the film was annealed at 200 • C for 4 hours in a vacuum chamber. A dehydration condensation between -NH and -OH was induced, resulting in PMDS-ODS polyimide formation. When an undiluted solution was spin coated on a Si substrate at 3000 rpm, a 3 µm thickness polyimide film was obtained. To optimize and lower the thickness to be suitable for MIM metasurface creation, solution was diluted with NMP. In this study, films of 45 nm and 335 nm thicknesses were made.
Metasurfaces were fabricated using conventional electron beam lithography (EBL), electron beam evaporation of metals and a lift-off process. First, a 200 nm thickness Au layer was thermally evaporated on the Si substrate with a 5 nm Cr adhesion layer, followed by 60 nm of vacuum polymerization of polyimide, carried out onto the Au film without any adhesion layer. Then, spin coating of positive tone EB resist, ZEP 520A 1:2 diluted by anisole, was spin coated at 1000 rpm for 10 seconds and 3000 rpm for 40 sec. Then it was prebaked at 130 • C on a hotplate for 10 min and used for EBL writing at 100 kV acceleration voltage, 5 nA of irradiation current with 350 µC/cm 2 area dose. After E-beam drawing the pattern was developed in ZED-N50 for 1 min, followed by thermal evaporation and subsequent lift-off for creation of 50 nm thickness Au disk MIM pattern. The same EBL pattern definition and lift-off process was used for creating metasurfaces with polyimide deposited by spin coating.

Characterization of reflection / radiation of metasurfaces
Optical characterization of reflection spectra of metasurface has performed using a microscope (IRT-1000) combined with a Fourier transform infrared spectrometer FTIR (FT-IR 6200). A 350 nm thickness Au film, which has 0 transmission and ∼98% reflectance, was used as a calibration reference. Thermal radiation emitted from metasurfaces was acquired using a similar FTIR system that was customized to use the MIM device as an external light source . Black body radiation (94% radiation blackbody paint) is used for the references of thermal radiation measurements. Four-level Lorenz model for electric permittivity can be described using the following equation:

Strong coupling modeling by FDTD
where the high (infinite) frequency permittivity ε ∞ = 3.68 and a reduction parameter f = 0.025 was used to attain parity with experimental results of the vibrational mode coupled plasmon resonance of MIM metasurface. The peak intensity ratio of g k , Lorenz parameter of central frequency ω 0,k , and full width at half maximum (FWHM) of the line width γ k was obtained from experimental IR-RAS spectra of polyimide, shown in Fig. 2(c) and Fig. 4 as well as summarized in Table 1.