Ultranarrow-Band Wavelength-Selective Thermal Emission with Aperiodic Multilayered Metamaterials Designed by Bayesian Optimization.

We computationally designed an ultranarrow-band wavelength-selective thermal radiator via a materials informatics method alternating between Bayesian optimization and thermal electromagnetic field calculation. For a given target infrared wavelength, the optimal structure was efficiently identified from over 8 billion candidates of multilayers consisting of multiple components (Si, Ge, and SiO2). The resulting optimized structure is an aperiodic multilayered metamaterial exhibiting high and sharp emissivity with a Q-factor of 273. The designed metamaterials were then fabricated, and reasonable experimental realization of the optimal performance was achieved with a Q-factor of 188, which is significantly higher than those of structures empirically designed and fabricated in the past. This is the first demonstration of the experimental realization of metamaterials designed by Bayesian optimization. The results facilitate the machine-learning-based design of metamaterials and advance our understanding of the narrow-band thermal emission mechanism of aperiodic multilayered metamaterials.

Introduction

All materials emit or absorb thermal radiation. Therefore, in the exploration to utilize various thermal energy resources, tailoring thermal radiation plays a fundamentally important role.[1−3] While conventional thermal radiators typically exhibit broad-band, polarization-independent, and omnidirectional emission, the technology to control thermal radiation is rapidly progressing with the development of the fields of nanophotonics and metamaterials. Electromagnetic metamaterials are artificially engineered materials with characteristics tailored over a broad range of wavelengths.[4,5] Wavelength-selective narrow-band thermal emission control is a key technology with applications in high-efficiency thermophotovoltaics,[6−8] incandescent light sources,[9] biosensing,[10−12] microbolometers,[13,14] imaging,[15] and infrared heaters.[16] Different types of artificial nanostructures have been proposed in the past few decades: multilayer,[17,18] photonic crystal,[19−21] and metal–insulator–metal (MIM) metamaterials.[22−27]

Development of metamaterial thermal radiators generally requires high-cost nanofabrication. The reported narrow-band thermal radiator with the highest Q-factor to date (∼200) consists of 2D-grating-coupled surface phonon polaritons.[28] However, there is still a problem because there are large unwanted peaks and background in the emissivity spectra in the target wavelength range. This can be quantified by the low value of the figure of merit (defined to evaluate the radiator performance as will be shown later) due to the low wavelength selectivity. In addition, including another experimental demonstration of multiple quantum wells and a photonic crystal slab with a Q-factor of 107,[29] the complicated fabrication process faces practical problems because many applications of radiators require large surface area. In this sense, among various classes of metamaterials, multilayers with relatively less complication in fabrication have merit in scalability. Control of thermal emission by multilayer structures has been successfully demonstrated with a Fabry–Perot resonator with a Q-factor of 87[30] and a distributed Bragg reflector with a Q-factor of 36.[31] However, these structures are usually realized by simple and periodic design, despite the fact that periodic structures are a tiny subset of the entire possible range of multilayer structures. Several studies have reported control of light by such “aperiodic” multilayer structures,[32−36] but these results were obtained by numerical simulation. Furthermore, the optimal design of aperiodic multilayered metamaterials with desired thermal emission characteristics has been difficult because the search space, i.e., the number of possible candidates, becomes enormous.

The key technology to overcome this challenge is “materials informatics” (MI), which has the capability to efficiently identify materials with preferred properties. MI aims to identify the “best” materials with optimal structure and/or composition using unrecognized complex correlations in the data. It has been applied to find novel materials such as cathode materials for the lithium-ion batteries,[37] nitride semiconductors composed of earth-abundant materials,[38] piezoelectric materials,[39] and thermoelectric materials.[40−44] While these works have aimed to realize high-throughput screening of the best materials from the pool of stoichiometric compounds, another course of MI aims to create nanostructures by identifying the optimal geometry that maximizes the objective properties. This includes nanoparticles embedded in a matrix to modulate heat conduction,[45] solid–solid interfaces to identify energetically stable structures,[46] and multicore structures of plasmonic nanowires to control optical scattering and cloaking effects.[47]

On the basis of the above progress in geometry optimization, the methodology using Bayesian optimization has been extended to the design of nanostructures with optimal thermal conductance[48] and thermoelectric figure of merit.[49] There, to efficiently identify the optimal structures among the enormous number of candidates, phonon/electron transport calculations and machine learning/prediction are alternately conducted. The previous works have shown that such an approach can considerably accelerate nanostructure design for transport properties. As the method is not limited to phonons/electrons and is applicable to any other quasi-particles, this work aims to perform such optimization for polaritons and associated thermal radiation. It should be noted that for thermal radiation there have been reports on the optimal design of multilayer structures using a genetic algorithm,[50,51] but genetic algorithms do not involve machine learning/prediction. In addition, recent studies[52,53] reported numerical nanophotonics designs based on neural networks. The essential drawback of their approach is that it is “exploitation-only”. There is plenty of evidence that the exploitation-only approach cannot be more efficient than the approach balancing exploitation and exploration.[54] On the other hand, Bayesian optimization identifies an unknown function with respect to the descriptors with as few iterations as possible, where, at every iteration, learning and prediction based on a Gaussian process are performed. Our approach uses Bayesian inference to quantify uncertainties and takes the optimal balance between exploration and exploitation, and we have used it to solve an essentially more difficult problem than the ones solved using neural networks. Although the previous studies optimized the thickness of each layer only, we optimized how the three materials are arranged, i.e., our method optimizes the ordering of the materials as well. There are a huge number of possible orderings, which adds substantial difficulty to the optimization problem.

In this work, we computationally designed an ultranarrow-band wavelength-selective thermal radiator via Bayesian optimization methods[55] and experimentally demonstrated the optical characteristics of the designed multilayered metamaterials. Potential applications of this study include infrared sensors, infrared imaging, and infrared heaters since the target wavelength is in the mid-infrared range.

Results and Discussion

Figure a shows a schematic of the optimization method with MI combining electromagnetic simulation and Bayesian optimization. The designed metamaterial is divided into N unit layers with thickness dt. A unit layer can be either Ge, Si, or SiO2. The choice of compositions are commonly used semiconductor and dielectric materials for their high and low refractive indices, respectively. Since tungsten was chosen as the substrate, the substrate was considered opaque. Four basic elements are required when materials informatics is performed: the descriptor, calculator, evaluator, and optimization method. The descriptors are used to describe possible structure candidates during the optimization process. In this study, we used a text flag to indicate the state of each layer: “1”, “2”, and “3” represent the Ge, Si, and SiO2 layers, respectively. Such a simple descriptor has been shown to realize efficient optimization[48,49] in addition to being intuitive, general, and practical, which are important in the actual material development. As for the calculator, we employed the transfer matrix method (TMM) to calculate the emissivity spectra (see Methods).

Figure 1

(a) Schematic of the optimization method with material informatics combining electromagnetic simulation and Bayesian optimization. (b) Schematic of the ideal optical property of the narrow-band thermal radiator.

The desired optical property of the ultranarrow-band thermal radiator is shown in Figure b. The ideal radiator has a sharp and high thermal emission at a target wavelength λt with a bandwidth Δλ, and low thermal emission in the rest of the infrared wavelength region to reduce radiative heat loss. For the evaluator of designed multilayered metasurfaces, a figure of merit (FOM) is defined as follows:where ε is the spectral normal emissivity, Ibλ is the spectral blackbody intensity, and λmin and λmax are the minimum and maximum wavelengths considered for the optimization.

As we have N unit layers and three possible materials (Ge, Si, or SiO2), the total number of candidate structures is 3, which becomes enormous for a useful range of N. For these large-scale problems, efficiency of optimization becomes critical, and thus, we need a method that surpasses conventional optimization tools. For this, we employ Bayesian optimization using the open-source library COMBO (see section S1 in the Supporting Information).

As shown in Figure a, suppose that FOMs of n candidates are initially calculated, and we are to select the next ones to calculate. A Bayesian regression function is learned from n pairs of descriptors and FOMs (i.e., training examples). For all of the remaining candidates, a predictive distribution of FOMs is estimated. Finally, FOMs are calculated for the selected candidates, and they are added to the training examples. By repetition of this procedure, the calculation of FOMs is scheduled optimally, and the optimized structure can be found quickly. Here, one problem is that the Bayesian optimization requires large computational memory because it uses information on the text data for all of the candidates. Therefore, we employed a hierarchical method to reduce the required size of computational memory, as will be explained later (also see Figure S1).

First, we computationally designed narrow-band thermal radiators with three candidate materials (Ge, Si, and SiO2) for a target wavelength λt of 6.0 μm. The wavelengths Δλ, λmin, and λmax were set to 4 nm, 4 μm, and 8 μm, respectively. The number of layers N was fixed at 18. Variation of the total thicknesses of the multilayers, ttotal, was also considered within the range from 3.6 to 4.0 μm with an increment of 0.02 μm, giving 21 variations of ttotal. Therefore, the total number of possible candidates is 318 × 21 = 8 135 830 269. It should be noted that it was not possible to account for structures with translational and reversal symmetries prior to the calculation to reduce the number of candidates. In this case, the numbers of initial and predicted candidate structures were set to 200 and 400, respectively. The computational load for this calculation was so large that all of the candidates could not be evaluated. For the sake of saving the computational memory, the optimization was pursued in hierarchical steps; the overall candidates were randomly divided into 42 000 groups, and the optimization was first performed for each group, after which the global best structure was identified by ranking these 42 000 local best structures. The total computational time was about 24 days on our cluster machine with 24 parallel computation (UNI-i9X, TOWA Electric, Inc.). The computational memory size in this work was about 128 GB, which set the maximum total number of layers to be 18. This could be enlarged by using a computer with a larger memory, but as the FOM of the designed structure is already close to unity, there is in fact not much room left for noticeable improvement even if we further increased the number of layers. Therefore, one can see the current setup to be nearly optimal and free from hardware restrictions.

The resulting optimized structures are shown in Figure a. It is interesting to note that the optimized structure with the maximum FOM consists of only Ge and SiO2 layers despite the fact that the optimization was performed including Si also. The obtained structure is a counterintuitive aperiodic multilayer, which is explicitly different from conventional multilayered thermal radiators with periodic structures. The total thickness ttotal of the optimal multilayer in this case is 3.80 μm.

Figure 2

(a) Optimized structure of the narrow-band thermal emitter with three material candidates (Ge, Si, and SiO2). The optimal structure turned out to consist of only Ge and SiO2 layers. (b) Histories of the FOMs of 20 randomly selected groups. The global-maximum FOM was found in a certain group that is indicated by the thick red line.

Figure b shows the history of the maximum FOM with respect to the number of calculated structures. Here we randomly chose the cases of 20 groups with about 200 000 candidates each to show the optimization efficiency and its statistics. The maximum FOM could be realized within calculations of 168 000 000 structures on average, which means only 2.06% of the candidate structures needed to be calculated to identify the optimal structure.

We also designed two other types of narrow-band thermal radiators with different target wavelengths of 5.0 and 7.0 μm. For these cases, using the finding in the case of λt = 6.0 μm that the optimal structure consists only of two species (Ge and SiO2), the optimization was performed for these two species instead of the above three species, which reduced the number of candidates to 218 × 21 = 5 505 024. The bandwidth Δλ (=4 nm) and the evaluation range of wavelengths (λmin = 4.0 μm and λmax = 8.0 μm) were kept the same as in the three-species optimization, and the number of initial candidate structures and predicted candidate structures were reduced to 100 and 20, respectively. The resulting optimized structures for λt = 5.0 and 7.0 μm (Figure a,b) consist of aperiodic multilayers similar to that for λt = 6.0 μm. The total thicknesses of the multilayers for the corresponding samples are ttotal = 3.78 and 3.96 μm, respectively.

Figure 3

Optimized structures of the narrow-band thermal emitters for the target wavelengths of (a) 5.0 and (b) 7.0 μm.

The computational load for the two-species calculations was relatively small, so all of the candidates could be calculated to validate the optimal structure and efficiency. As a result, the optimal structures obtained by Bayesian optimization were confirmed to be exactly the same as the structures with maximum FOM among all of the candidates. We also confirmed from the probability distributions (see Figure S2) that the probability monotonically decreases as the FOM value approaches the maximum without noticeable local minima, indicating that the current problem is suited for Bayesian optimization.

Figure a shows the corresponding calculated spectral directional emissivities of the optimized structures. Extremely sharp and high emissivity can be realized with the optimized structures, and there are no extra peaks within the wavelength range of interest (from 4 to 8 μm). The corresponding emissivities of the peaks are unity, and their Q-factors are 217, 273, and 233 for λt = 5, 6, and 7 μm, respectively.

Figure 4

(a) Calculated spectral directional emissivities of the optimized structures obtained with Bayesian optimization and (b) measured spectral directional emissivities of the fabricated structures aimed at λt = 5.0 μm (red), 6.0 μm (blue), and 7.0 μm (green). (c) Cross-sectional TEM images of the fabricated sample for λt = 6.0 μm.

Finally, we experimentally fabricated the optimized structures by sputtering to demonstrate the feasibility of the structural optimization. Figure b shows the measured spectral directional emissivities of the fabricated structures. The three sharp peaks that correspond to the ones seen in the numerical simulations can be clearly observed, although the locations of the peaks are shifted by about 0.5 μm relative to the designed structures. The obtained peak emissivity values of the λt = 5, 6, and 7 μm samples are 0.76, 0.83, and 0.61, and the Q-factors are 132, 188, and 109, respectively. The reason for the discrepancies in the peak positions and emissivities/Q-factors of the designed and fabricated structures could be that the thicknesses of the constituent layers in the fabricated samples somewhat deviate from the designed values. Table quantifies the moderate but non-negligible differences between the layer thicknesses of the designed and fabricated structures obtained from the cross-sectional transmission electron microscopy (TEM) image for λt = 6.0 μm (Figure c) and the cross-sectional scanning electron microscopy (SEM) images for λt = 5.0 μm and λt = 7.0 μm (Figure S3). When we calculated the spectral directional emissivity for the layer thicknesses in the fabricated sample (Table ), the position of the peak approached the experimentally measured value (Figure S4). The remaining discrepancy can be attributed to the minor differences in the optical properties of the sputtered material and those used as inputs to the numerical simulation, since the optical properties may differ depending on fabrication conditions such as the deposition rate.

Table 1

Layer Thicknesses of the Designed and Fabricated Structures (in μm)

	λt = 5.0 μm		λt = 6.0 μm		λt = 7.0 μm
layer no.	sim.	exp.	sim.	exp.	sim.	exp.
1	0.42	0.42	0.42	0.43	0.44	0.44
2	0.63	0.61	0.63	0.69	0.66	0.62
3	0.42	0.43	0.42	0.45	0.44	0.44
4	1.05	0.97	0.85	0.91	0.88	0.84
5	0.63	0.63	0.85	0.87	0.44	0.45
6	0.63	0.58	0.63	0.65	0.22	0.22
7	–	–	–	–	0.44	0.44
8	–	–	–	–	0.44	0.41

To determine the sharpness of the interface, the atomic concentrations at the Ge–SiO2 interface were observed by energy-dispersive X-ray spectroscopy (EDX) (Figures S5 and S6), and the interface was confirmed to be sharp with small interdiffusion. Although fabrication with a more accurately calibrated sputtering process would improve the reproduction of the designed performance, which remains to be our future task, the key features in the designed structure, namely, ultranarrow-band emission with controlled peak wavelength, were clearly realized in the experiments. The obtained Q-factors are about 217–273 in the computational design and about 109–188 in the experiment, which are significantly higher than the values reported in the previous studies. In addition, the FOM of this work is significantly higher than in previous experimental work:[28] the FOM of previous work, evaluated with the same wavelength range around the target wavelengths, was only 0.02 for 0° and further decreased to −0.16 for 1°, which are considerably smaller than the FOM of 0.77 for the current structure aimed at 6 μm. Although the locally extracted Q-factor in the previous work reached 200,[28] the emissivity spectra had much larger background and unwanted peaks, and thus, our experiment exhibits significantly higher wavelength selectivity. To our knowledge, this is the first demonstration that narrow-band thermal radiators designed by machine learning can be realized in experiments.

We now discuss the mechanism of the enhanced emission in terms of the magnetic field profiles shown in Figure . The intensities of the magnetic profiles were normalized by the intensity of the normal incident light. In Figure a,b for λt = 5 and 6 μm, there are strong confinements of electromagnetic energy in the Ge layer. On the other hand, in Figure c, for λt = 7 μm, strong confinement can be observed in the SiO2 layer. These emissivity enhancements originate from localized modes, similar to defect modes of photonic crystals.[56] Defect modes of photonic crystals exist inside a photonic band gap; therefore, this phenomenon is usually observed with periodic structures (see section S3). However, it is interesting to note that we observed a similar localized mode inside the aperiodic multilayered metamaterials. In other words, two or more optimized defect layers are introduced into the photonic crystals that effectively serve to constitute a sharp peak in the emissivity. In particular, in Figure c, the defect layer corresponds to three layers of a thin SiO2 layer and upper and lower Ge layers sandwiching the SiO2 layer. Therefore, the aperiodic structure, when optimized, successfully suppresses the unnecessary emissivity peaks due to higher-order harmonics, or in other words, shifts the peaks to a shorter wavelength range. To quantify how much power is absorbed by the proposed structure, the power dissipation density w was calculated as[57]where ε0 is the permittivity of vacuum, ε″ is the imaginary part of the complex dielectric function, and ω is the angular frequency. The power dissipation densities, which are shown in Figure d–f, indicate the strong absorption at the tungsten substrate, although there is weak absorption within the SiO2 layer. Therefore, thermal energy dissipation mainly occurs in the metallic substrate because of the large optical loss.

Figure 5

(a–c) Contour plots of normalized magnetic field intensity and (d–f) power dissipation density for target wavelengths of (a, d) 5.0 μm, (b, e) 6.0 μm, and (c, f) 7.0 μm.

Because of the localized mode of the electromagnetic wave, the proposed emitter has an angular dependence of the optical properties (Figure S8). Isotropic thermal emission is preferred in certain applications such as infrared heaters. In this design, the angular dependences of the optical properties of transverse magnetic and transverse electric polarization within 20° are small, as the spectral shifts were only about 1%, which therefore is acceptable for practical applications. It should be noted that it is also possible to include the angular dependence in the FOM for preferred angular dependence, which will be explored in the future. The obtained results enhance our understanding of the narrow-band thermal emission mechanism of aperiodic multilayered metamaterials and facilitate the effective design of new metamaterials via Bayesian optimization.

Conclusion

We computationally designed ultranarrow-band wavelength-selective thermal radiators via Bayesian optimization methods and experimentally demonstrated the optical characteristics of the designed multilayered metamaterials. The optimized structures could be found within calculations of only a few percent of the total numbers of candidate structures. The optimized structure for each target wavelength consists of aperiodic multilayers that give rise to sharp and near-unity emissivity. The designed structures were experimentally realized with reasonable accuracy, and the obtained structures exhibit Q-factors significantly larger than in previous works based on empirical design. Post-analysis of the magnetic fields of the structures revealed that the aperiodic multilayers can result in highly effective localization. The current work demonstrates the effectiveness, feasibility, and accuracy of developing narrow-band thermal emission materials using Bayesian optimization. In addition, the follow-up analysis of the mechanism demonstrates that such a materials informatics approach is also useful to enhance our understanding of narrow-band thermal emission.

Methods

Safety Statement

No unexpected or unusually high safety hazards were encountered.

Electromagnetic Simulation

The TMM was used to solve Maxwell’s equations, allowing the calculation of the spectral radiative properties of multilayered metamaterials.[58] The spectral directional emissivity could be obtained by applying Kirchhoff’s law, i.e., ελ = 1 – Rλ, where Rλ is the reflectance obtained from the TMM simulation. The dielectric functions of SiO2, Si, Ge, and W were obtained from tabulated data.[59]

Bayesian Optimization

Bayesian optimization is a design algorithm based on machine learning[60] and a well-established technique for black-box optimization.[55] Bayesian prediction models are employed to predict the black-box function, where the uncertainty of the predicted function is also evaluated as predictive variance. The next candidate for the experiment is selected on the basis of predicted values and variances. Bayesian optimization has been recognized as an important technique in machine learning research because of successful hyperparameter tuning in deep learning algorithms. Bayesian optimization can be applied not only to materials sciences but also to various kinds of problems. However, the precondition is that each candidate point is represented as a numerical vector of identical dimensionality (i.e., descriptor).

Sample Fabrication and Reflectivity Measurement

The narrow-band thermal radiators designed on the basis of the Bayesian optimization method were experimentally fabricated and characterized. SiO2 and Ge layers were alternately deposited on a tungsten substrate by a magnetron sputtering machine. An FTIR spectrometer (iS50R, Thermo Scientific Nicolet) was used for reflectivity measurements, with an opaque gold film as a reference. In order to avoid atmospheric absorption, the measurements were conducted with flowing nitrogen gas. The incident angle was arranged within 1°, and therefore, the measured spectral reflectivity data could be considered as near normal reflectivity. Once the reflectivity was obtained, the spectral directional emissivity was obtained by applying Kirchhoff’s law.