energies

: It is vital to further improve the design of TPV thermal emitters since the energy efﬁciency of thermophotovoltaic (TPV) systems is still not adequately high. In this paper, we propose a novel evaluator for the optimization of TPV thermal emitters, namely the percentage of effective ﬁgure (PEF) to replace the ﬁgure of merit (FOM). The associated algorithm, Bayesian optimization nesting simulated annealing (BOnSA), is developed to achieve better performance. By searching throughout the whole parameter space and then optimizing in a reduced space, BOnSA can lead to a satisfactory solution numerically for GaSb photovoltaic (PV) cells. When designing the emitter, the aperiodic material structure with an anti-reﬂection substructure and Fabry–Perot etalon is constructed from the material candidates. In particular, one of the optimal structures determined by BOnSA is {SiO 2 , ZnS, Ge, MgF 2 , W, Si, SiO 2 , W} with the value of PEF = 0.822, which is better than the previous work by comparison. Moreover, by applying BOnSA to various structures, we have obtained higher values of PEF with less time cost, which thus veriﬁes the efﬁciency and scalability of BOnSA. The results of our paper show that BOnSA provides an effective approach to the thickness optimization problem and that BOnSA is applicable in other relevant scenarios.


Introduction
The energy crisis has prompted scientists to seek new energy that can be applied on a large scale [1][2][3].In the past decades, solar energy is one of the most promising new energy sources [4][5][6][7] which can be utilized by thermophotovoltaic (TPV) systems [8][9][10].Theoretically, TPV systems can overcome the Shockley-Queisser limit [11,12] and obtain high energy conversion efficiency.However, nowadays, the energy conversion efficiency of TPV systems is still relatively low [13,14].One way to improve the efficiency is to put more effort into the design of TPV thermal emitters.Since photovoltaic (PV) cells in TPV systems can only convert photons of specific wavelengths into electric energy, it is crucial to improve the selectivity of TPV thermal emitters [15].Great efforts have been made to fabricate emitters with different materials including the luminescent bands of rare earth oxides [16,17] and optical crystals [18,19].However, these materials have their own shortcomings due to poor controllability and scarcity of raw materials [20].Thanks to the introduction of negative refractive index in the last century [21], thermal metamaterials have been deeply explored [22][23][24] and used extensively to engineer thermal emission [25][26][27][28].As for the design of TPV thermal emitters by means of metamaterials, there are in general two methods that have been implemented.One is metasurface design [29][30][31], and the other is multilayer structure design [32,33].The associated examples of these two design methods are listed in Table 1.Metasurface design varies significantly in different research works [34][35][36], whose corresponding structures are mostly designed by experience and lack reasonable explanation.By contrast, for multilayer structure design, the Fabry-Perot Energies 2023, 16, 416 2 of 16 etalon [37,38] and the periodic structure [39,40] are the most common.Although multilayer structures can be designed with good performance in most established work, the majority of researchers are inclined to use repeating segments for their materials, without subsequently providing physical explanation of the repeating segments [41][42][43] or the stacking order of chosen materials [44].It is undeniable that these two structures have contributed to the progress of thermal emitters, but the materials used in such structures are usually selected manually.Both lack an instructive guide for the design procedure.

Metasurface
High temperature selective emission An array of crosses above a backplane Pt, Al 2 O 3 , EBR [29] High efficiency thermal emitter optimization Topology optimized designs TiN, Si 3 N 4 [30] Rapidly generated globally optimized meta structure Topology optimized designs TiN, SiN [31] Light modulation and control A representative nanoantenna array Gold [34] Wavelength-selective near-infrared metasurface emitter A periodic metallic disk pattern on a dielectric film W, SiO 2 [35] Multilayer Tunable narrowband Silicon-based thermal emitter SiN/SiNO photonic crystal and a metallic W film on Si substrate SiN y O z , SiN x , Si, W [33] High efficiency rare-earth emitter for TPV system ErAG emitter with chirped-mirror (Q-matched Fabry-Perot modes) ErAG, other filters [37] Polariton-enhanced emittance for selective thermal emitter Metallic-dielectric multilayer structure W, SiO 2 [39] Achieve dual-band nonreciprocal thermal radiation A 1-D PC (AC) n , a 1-D PC (CA) m and a metal layer InAs, SiO 2 , Al [40] Selective TPV emitter with aperiodic multilayer structure Aperiodic multilayer structure W, Si, SiO 2 [41] A critical procedure in the design of thermal emitters is the optimal choice of materials with the best thickness of each layer.Therefore, the optimization algorithm and evaluator play a key role in this step.Regarding algorithms, besides traditional ones including simulated annealing (SA) [45] and genetic algorithm (GA) [46,47], a number of machine learning-related algorithms have come into our vision, such as Bayesian optimization (BO) [41,42].Different optimization algorithms have been applied to different work, but the evaluator the figure of merit (FOM), is commonly shared [41,42,48].However, the non-ideal performance of FOM can be easily noticed, especially when optimizing in the near-infrared region (NIR) [49].Moreover, the conventional BO costs much time and calculation particularly when using Gaussian process (GP) as its model [50], and has difficulties in jumping out of local optima [51].There is a drawback called dimension disaster in BO, which means the performance of BO in high dimension is almost the same as random search (RS) [52].The issues mentioned above drive us to pursue a more appropriate scheme, aimed for a better design of the thermal emitter and its optimization algorithm.
In our work, we propose a new evaluator to replace FOM and develop an optimization algorithm named Bayesian optimization nesting simulated annealing (BOnSA).As for determining the optimal structure for emitters, a novel method is put forward by referring to periodic materials and a Fabry-Perot etalon.In a word, a more general multilayer thermal emitter structure could be designed by inserting the periodic materials into the Fabry-Perot etalon and adding an anti-reflection substructure on the top of the Fabry-Perot etalon.Taking the time cost and final performance into consideration, we find that the best structure can be gained after one repeating segment.Therefore, the optimal structure we proposed is aperiodic.By performing RS to optimal structures with different material combinations, the algorithm can filter out several preliminary candidates and determine the material used in each layer of emitters for GaSb and InGaAsSb PV cells.Then BOnSA can be employed to them.To illustrate the scalability, BOnSA is also applied to optimize the thermal emitters in other structures.
The mathematical models and methods in our research are presented in Section 2 (more details are provided in Appendices A and B), where PEF is defined to serve as a different but efficient evaluator.By utilizing the transfer-matrix method, the mathematical model of n-layer multilayer TPV emitters can be obtained, characterizing the underlying dynamics behind the algorithm.By merging BO into SA, this combined optimization algorithm BOnSA is developed to take advantage of the strengths of BO and SA.The application and performance of BOnSA is discussed in detail in Section 3.After tuning preliminary parameters, we define the structure of TPV thermal emitters.The basic structure of emitters can then be divided into two parts, anti-reflection substructure and FP substructure.RS is employed to determine the initial values of PEF with regard to different combinations.The structure with the highest value of PEF is taken to the next step, which is the optimized by BOnSA.Due to the ingenious treatment of the local optima, the optimized structure can be quickly achieved.In the numerical simulations, it can be seen that the optimal structures perform well in different scenarios.Furthermore, we analyze and explain the reason why BOnSA can show a better performance.We also study the mechanism of optimal structures for thermal emitters, followed by conclusions given in Section 4.

PEF Evaluator
In this section, a new evaluator, the percentage of effective figure (PEF), will be specified.More precisely, PEF is defined to obtain a sharp peak in emission spectrum, which can be descried as follows in a mathematical sense: where λ 1 , λ 2 correspond to the lower and upper bounds in the given range of spectral emissivity, respectively.As displayed in Figure 1, the required emission spectrum encloses an ideal area in the target band (∆λ) and equals 0 out of the target band.However, the actual emission spectrum is composed of the effective area (S 2 ) in the target band, and the invalid areas (S 1 , S 3 ) out of the target band.Naturally, by taking the combination of these areas and then taking the ratio of the combination to the effective area, PEF can be written more succinctly as PEF = (S 2 − S 1 − S 3 )/ S 2 .
Energies 2023, 16, x FOR PEER REVIEW 4 of 17 emissivity spectrum into three parts.It can be easily observed from the definition that the larger the area of emissivity in the target band is the higher the value of PEF is (the better the performance of system is).Different from FOM, PEF measures the performance of TPV emitters at the level of geometric significance, which has nothing to do with the width and location of target band.

BOnSA Algorithm
Firstly, concerning the computation of spectral emissivity, we utilize the transfer-matrix method (TMM) [54], which has the following form: FOM has been widely used to evaluate the thermal performance of TPV emitters [41,42,48,53].However, it can be easily noticed that the ratio terms in FOM correspond to the same weight, which means even a high ratio in a rather short band could lead to a significant drop in FOM, particularly when optimizing in NIR.It is unreasonable that ratios in fairly short bands beyond the target range have the same effect as a high emissivity in the target band.One common practice to avoid the impact of ratios in a short band is to multiply a factor.However, this would raise a problem of tuning a new hyper-parameter.Therefore, compared to FOM, PEF could avoid the problem mentioned above due to its mathematical feature.Based on the waveform of target emissivity, PEF divides the emissivity spectrum into three parts.It can be easily observed from the definition that the larger the area of emissivity in the target band is the higher the value of PEF is (the better the performance of system is).Different from FOM, PEF measures the performance of TPV emitters at the level of geometric significance, which has nothing to do with the width and location of target band.

BOnSA Algorithm
Firstly, concerning the computation of spectral emissivity, we utilize the transfer-matrix method (TMM) [54], which has the following form: where E + k and E − k represent the forward and backward waves of the electric field in the k-th layer.More details about TMM can be found in Appendix A. The reflectivity and transmissivity are then given by R λ = (T 21 /T 11 ) 2 and T λ = (1/T 11 ) 2 , respectively.According to Kirchhoff's law of thermal radiation [55], we have that λ = A λ = 1 − R λ − T λ , where A λ is the absorptivity (the energy dissipation of scattering can be ignored).
Before proposing the algorithm, we show assumptions for our model.Firstly, the system is assumed to operate under suitable conditions (T = 1200 K) and the refractive index of relevant materials is available for the corresponding situation.Secondly, we assume that the direction of the incident light is perpendicular to the plane of the thermal emitter.This can be achieved by appropriate design of the optical structure.We also assume that influences of the factors outside the given range of spectral emissivity can be neglected.
It is widely shared that SA can escape from local optima [56,57], which helps make up for the deficiency of BO, but requires a lengthy time to achieve theoretical convergence.Therefore, the speed of SA could be significantly improved via incorporating BO during its operation process.Nevertheless, the GP modeling strategy adopted by typical BO becomes a stumbling block for efficiency.A convincing reason for abandoning the conventional BO based on GP is exactly its huge cost of time and space.
Therefore, we introduce BO based on tree-structured parzen estimator (TPE) [58,59].TPE divides the data into two categories with diverse density distributions, which are represented as follows: TPE updates the agent model by Gaussian mixture clustering model (GMM) [60] and adopts Expected Improvement (EI).EI has the following form for the minimum [60]: where γ = p(y < y * ) and M stands for the model-based algorithms approximation.Its solution for the minimum is x * = argmin x EI y * (x, M t−1 ).As for implementation, the algorithm has the setting of γ = 0.2 and y = −PEF.Due to its simplicity in contrast to GP, TPE not only runs faster and costs less, but also has excellent efficiency [60].
Energies 2023, 16, 416 As illustrated in Figure 2, the SA executes the external epoch.In each epoch, the Metropolis algorithm [61,62] selects n layers to generate a new solution, which follows a uniform distribution.If the new solution is rejected, the algorithm will go to the next epoch, otherwise BOnSA will narrow down the scope and undergo multiple BOs.It turns out that the greater the PEF of a new solution is, the higher expectation of PEF in its neighborhood will be.Through computation accomplished by BO, BOnSA will probably reach local optima.In fact, to jump out of the local region, it is reasonable to end BO by continuing the next epoch.After applying SA, the so-called greedy search will traverse the materials from top to bottom, optimizing the thickness of each layer.Completing the procedure mentioned above, a fulfilling solution can be obtained eventually.
where  = (  * ) and  stands for the model-based algorithms approximation.Its solution for the minimum is  * =  min EI * (,  ).As for implementation, the algorithm has the setting of  = 0.2 and  = −PEF.Due to its simplicity in contrast to GP, TPE not only runs faster and costs less, but also has excellent efficiency [60].
As illustrated in Figure 2, the SA executes the external epoch.In each epoch, the Metropolis algorithm [61,62] selects n layers to generate a new solution, which follows a uniform distribution.If the new solution is rejected, the algorithm will go to the next epoch, otherwise BOnSA will narrow down the scope and undergo multiple BOs.It turns out that the greater the PEF of a new solution is, the higher expectation of PEF in its neighborhood will be.Through computation accomplished by BO, BOnSA will probably reach local optima.In fact, to jump out of the local region, it is reasonable to end BO by continuing the next epoch.After applying SA, the so-called greedy search will traverse the materials from top to bottom, optimizing the thickness of each layer.Completing the procedure mentioned above, a fulfilling solution can be obtained eventually.SA [45] and BO [41,42] have been applied to optimize the design of thermal emitters.However, both of them have their own disadvantages.We thus combine the two algorithms by taking advantage of the strengths of each one to develop a new optimization algorithm BOnSA, which uses SA for global optimization and BO for local optimization.BOnSA effectively averts SA from redundant global search and prevents BO from being stuck at a local optimum.SA [45] and BO [41,42] have been applied to optimize the design of thermal emitters.However, both of them have their own disadvantages.We thus combine the two algorithms by taking advantage of the strengths of each one to develop a new optimization algorithm BOnSA, which uses SA for global optimization and BO for local optimization.BOnSA effectively averts SA from redundant global search and prevents BO from being stuck at a local optimum.

Determination of the Material Structure
The fundamental structure of TPV thermal emitters is designed with two substructures, namely the anti-reflection substructure and the Fabry-Perot etalon.The main feature of the Fabry-Perot etalon is the two W layers on top and bottom, as well as materials in the middle, with low refractive index and almost zero extinction coefficient.A brief explanation is that the Fabry-Perot etalon has been proved to provide enhancement to the selectivity of a filter [63], which greatly increases the emissivity of the system in the target band.The additional anti-reflection substructure on its top plays a role in reducing the reflection of near-infrared light in the target band, which can be compared to a hat put on the Fabry-Perot etalon.For convenience, the anti-reflection substructure and Fabry-Perot etalon are abbreviatedly referred to as the hat substructure and FP substructure, respectively.
The overall process of determining the optimal material structure is displayed in Figure 3.In the fundamental structure, we take the combination of these candidates to form various hat substructures inspired by the periodic multilayer emitters and we choose the FP substructure which contains the periodic materials.Taking the ascending time cost into consideration, the basic unit of periodic materials is set to two, and materials are selected from the candidates as MgF 2 , ZnS, Si and SiO 2 .For all the possible fundamental structures, RS is performed to confirm the preliminary PEF of each one.Then sorted by the values of PEF, the substructures that appear with high frequencies in the front will be selected for further optimization.Eventually, through applying BOnSA, the optimal structure distinguishes itself in the last candidates.
explanation is that the Fabry-Perot etalon has been proved to provide enhancement to the selectivity of a filter [63], which greatly increases the emissivity of the system in the target band.The additional anti-reflection substructure on its top plays a role in reducing the reflection of near-infrared light in the target band, which can be compared to a hat put on the Fabry-Perot etalon.For convenience, the anti-reflection substructure and Fabry-Perot etalon are abbreviatedly referred to as the hat substructure and FP substructure, respectively.
The overall process of determining the optimal material structure is displayed in Figure 3.In the fundamental structure, we take the combination of these candidates to form various hat substructures inspired by the periodic multilayer emitters and we choose the FP substructure which contains the periodic materials.Taking the ascending time cost into consideration, the basic unit of periodic materials is set to two, and materials are selected from the candidates as MgF2, ZnS, Si and SiO2.For all the possible fundamental structures, RS is performed to confirm the preliminary PEF of each one.Then sorted by the values of PEF, the substructures that appear with high frequencies in the front will be selected for further optimization.Eventually, through applying BOnSA, the optimal structure distinguishes itself in the last candidates.As shown in Figure 4a, the values of PEF can reach a high level when the repetition of periodic materials in FP substructures equals 1.Additionally, one can obtain the number of layers in the anti-reflection substructure.It can also be observed that, in terms of the four layers on top of W, anti-reflective layers do perform better than others, due to its larger ranging space and higher maximal value of PEF.To sum up, the optimal fundamental structure can be expressed as {x1, x2, x3, x4, W, y1, y2, W}, where xj (for j = 1, 2, 3, 4) and yj (for j = 1, 2) denote different materials.As shown in Figure 4a, the values of PEF can reach a high level when the repetition of periodic materials in FP substructures equals 1.Additionally, one can obtain the number of layers in the anti-reflection substructure.It can also be observed that, in terms of the four layers on top of W, anti-reflective layers do perform better than others, due to its larger ranging space and higher maximal value of PEF.To sum up, the optimal fundamental structure can be expressed as {x 1 , x 2 , x 3 , x 4 , W, y 1 , y 2 , W}, where x j (for j = 1, 2, 3, 4) and y j (for j = 1, 2) denote different materials.
By analyzing the values of PEF with regard to different combinations, the material in each layer can now be confirmed.After applying BOnSA to optimize the structures, the values of PEF for all the combinations are shown in Figure 4b.Apparently, the structure corresponding to {SiO 2 , ZnS, Ge, MgF 2 } as the anti-reflection substructure and {W, Si, SiO 2 , W} as the Fabry-Perot etalon has the best performance with the maximal value of PEF = 0.822.As a result, one of the optimal structures is determined as {SiO 2 , ZnS, Ge, MgF 2 , W, Si, SiO 2 , W}. Please note that Ge has been added to the candidates for antireflection substructure.Many previous studies have shown that one of the good choices for fabrication of PV cells is Ge [64,65], which possesses the features of low bandwidth and high refractive index.According to Kirchhoff's law of thermal radiation [55], Ge can help increase the emissivity (or absorptivity) of the target band, which has been verified by relevant experimentation [66].Moreover, in analogy to Si, a material that can be used to build PV cells [67] and multilayer thermal emitters [68], we believe that Ge should also be a good candidate, the feasibility of which has indeed been supported by previous work [42,69].Ge has also been proved to significantly improve the performance of thermal emitters in our work.By analyzing the values of PEF with regard to different combinations, the material in each layer can now be confirmed.After applying BOnSA to optimize the structures, the values of PEF for all the combinations are shown in Figure 4b.Apparently, the structure corresponding to {SiO2, ZnS, Ge, MgF2} as the anti-reflection substructure and {W, Si, SiO2, W} as the Fabry-Perot etalon has the best performance with the maximal value of PEF = 0.822.As a result, one of the optimal structures is determined as {SiO2, ZnS, Ge, MgF2, W, Si, SiO2, W}.Please note that Ge has been added to the candidates for anti-reflection substructure.Many previous studies have shown that one of the good choices for fabrication of PV cells is Ge [64,65], which possesses the features of low bandwidth and high refractive index.According to Kirchhoff's law of thermal radiation [55], Ge can help increase the emissivity (or absorptivity) of the target band, which has been verified by relevant experimentation [66].Moreover, in analogy to Si, a material that can be used to build PV cells [67] and multilayer thermal emitters [68], we believe that Ge should also be a good candidate, the feasibility of which has indeed been supported by previous work [42,69].Ge has also been proved to significantly improve the performance of thermal emitters in our work.

The Emission Spectrum of Different Optimal Structures
By employing BOnSA, the best optimal structure for a GaSb PV cell can be obtained in the above section, with the thickness of each layer being {140.0,60.0, 129.6, 423.6, 10.5, 50.6, 8.4, 205.0} nm.As a complement, BOnSA is applied to another structure with great optimization potential, whose structure is confirmed as {MgF2, ZnS, Ge, SiO2, W, ZnS, SiO2, W} with thickness being {148.4

The Emission Spectrum of Different Optimal Structures
By employing BOnSA, the best optimal structure for a GaSb PV cell can be obtained in the above section, with the thickness of each layer being {140.0,60.0, 129.6, 423.6, 10.5, 50.6, 8.4, 205.0} nm.As a complement, BOnSA is applied to another structure with great optimization potential, whose structure is confirmed as {MgF 2 , ZnS, Ge, SiO 2 , W, ZnS, SiO 2 , W} with thickness being {148.4,59.1, 129.6, 403.3, 17.3, 93.2, 23.4,201.6} nm.We also optimize the structures in the references [41,70] to see the performance of BOnSA proposed in this paper.The emission spectra of different structures are illustrated in Figure 5a.
In terms of GaSb PV cells, it is clear that the emissivity is lower in the band higher than λ pv , which can significantly reduce the thermalization loss of the TPV system.Moreover, in the required band ∆λ, the emissivity curve of optimal structure can also maintain a steeper and plumper peak.According to Kirchhoff's law of thermal radiation, low reflectivity can be obtained only in the target area, resulting in a higher emissivity.Generally, as shown by the solid line in Figure 5a, the shape of the curve in the target area is like a rectangular pulse, which is close to the desired one.Thus, it can be concluded that the structure proposed by us performs well in the target band.Similarly, BOnSA can also be employed to the structure designed for the InGaAsSb PV cell, which can be seen in Figure 5b.Compared to the structure obtained in the original paper [41], the emissivity of the optimal structure has relatively high values (nearly 1.0) via BOnSA in our paper with less fluctuation in the target band, which indicates the TPV thermal emitter can absorb as much energy as possible that can then be consumed by PV cells.
optimize the structures in the references [41,70] to see the performance of BOnSA proposed in this paper.The emission spectra of different structures are illustrated in Figure 5a.In terms of GaSb PV cells, it is clear that the emissivity is lower in the band higher than  , which can significantly reduce the thermalization loss of the TPV system.Moreover, in the required band Δ, the emissivity curve of optimal structure can also maintain a steeper and plumper peak.According to Kirchhoff's law of thermal radiation, low reflectivity can be obtained only in the target area, resulting in a higher emissivity.Generally, as shown by the solid line in Figure 5a, the shape of the curve in the target area is like a rectangular pulse, which is close to the desired one.Thus, it can be concluded that the structure proposed by us performs well in the target band.Similarly, BOnSA can also be employed to the structure designed for the InGaAsSb PV cell, which can be seen in Figure 5b.Compared to the structure obtained in the original paper [41], the emissivity of the optimal structure has relatively high values (nearly 1.0) via BOnSA in our paper with less fluctuation in the target band, which indicates the TPV thermal emitter can absorb as much energy as possible that can then be consumed by PV cells.

The Optimization Trajectories of BOnSA and Other Algorithms
The optimization process for one of the optimal structures ({SiO2, ZnS, Ge, MgF2, W, Si, SiO2, W}) is depicted in detail in Figure 6.Here the dots represent all the local PEF in trials and the lines stand for the historical optimal PEF.Based on the distribution histograms of PEF, the local optima generated by BOnSA is much better than GA and RS.Compared to BOnSA, although GA can give a global optimum, its efficiency is completely defeated by BOnSA since GA costs much more time.
In the early stage of optimization, BOnSA quickly reaches a temporary stable level, thanks to the boost via multiple BOs.The graph of historical optimal PEF has several ascending phases, resulting from the boost via BO in small regions.It can be understood in the sense that BO based on TPE can quickly approach the local optima, while SA will transfer from one area that may have been full excavated to another area to be exploited.In these simulations, it takes 3 to 5 min to obtain a good solution on our laptop (AN515-55-57C4, Acer, Inc., Suzhou, China), which also proves the efficiency of BOnSA.

The Optimization Trajectories of BOnSA and Other Algorithms
The optimization process for one of the optimal structures ({SiO 2 , ZnS, Ge, MgF 2 , W, Si, SiO 2 , W}) is depicted in detail in Figure 6.Here the dots represent all the local PEF in trials and the lines stand for the historical optimal PEF.Based on the distribution histograms of PEF, the local optima generated by BOnSA is much better than GA and RS.Compared to BOnSA, although GA can give a global optimum, its efficiency is completely defeated by BOnSA since GA costs much more time.

Sensitivity Analysis
In terms of the sensitivity analysis in comparison to other work, BOnSA shows strong reliability.
As demonstrated in Figure 7, we analyze the sensitivity in our optimal structure by changing the thickness or varying the incident angle slightly.The gray area indicates the variation of emissivity when the thickness of one or more layers varies 1 to 3 nm, as a simulation of deviation that may arise in fabrication.Apparently, the slight change under the shadow almost does not affect its stability a little.Similarly, as the incident angle varies from 0 to 20°, the PEF value of the optimized structure changes by only 1 × 10 − , which illustrates the stability of the structure as well.In the early stage of optimization, BOnSA quickly reaches a temporary stable level, thanks to the boost via multiple BOs.The graph of historical optimal PEF has several ascending phases, resulting from the boost via BO in small regions.It can be understood in the sense that BO based on TPE can quickly approach the local optima, while SA will transfer from one area that may have been full excavated to another area to be exploited.In these simulations, it takes 3 to 5 min to obtain a good solution on our laptop (AN515-55-57C4, Acer, Inc., Suzhou, China), which also proves the efficiency of BOnSA.As demonstrated in Figure 7, we analyze the sensitivity in our optimal structure by changing the thickness or varying the incident angle slightly.The gray area indicates the variation of emissivity when the thickness of one or more layers varies 1 to 3 nm, as a simulation of deviation that may arise in fabrication.Apparently, the slight change under the shadow almost does not affect its stability a little.Similarly, as the incident angle varies from 0 to 20 • , the PEF value of the optimized structure changes by only 1 × 10 −4 , which illustrates the stability of the structure as well.

Sensitivity Analysis
In terms of the sensitivity analysis in comparison to other work, BOnSA shows strong reliability.
As demonstrated in Figure 7, we analyze the sensitivity in our optimal structure by changing the thickness or varying the incident angle slightly.The gray area indicates the variation of emissivity when the thickness of one or more layers 1 to 3 nm, as a simulation of deviation that may arise in fabrication.Apparently, the slight change under the shadow almost does not affect its stability a little.Similarly, as the incident angle varies from 0 to 20°, the PEF value of the optimized structure changes by only 1 × 10 − , which illustrates the stability of the structure as well.By contrast, BOnSA has also been applied to other structures designed for the GaSb PV cell, such as {Si3N4, W, Si3N4, W} [70] and {SiO2, Si, W, Si, W} [41].As shown by the dotted lines in Figure 7, when it comes to the low band close to  − Δ/2, a sharp decline appears and heavily affects the efficiency of the GaSb PV cell in the other two works.However, our approach can maintain a high PEF in the whole target band.By contrast, BOnSA has also been applied to other structures designed for the GaSb PV cell, such as {Si 3 N 4 , W, Si 3 N 4 , W} [70] and {SiO 2 , Si, W, Si, W} [41].As shown by the dotted lines in Figure 7, when it comes to the low band close to λ c − ∆λ/2, a sharp decline appears and heavily affects the efficiency of the GaSb PV cell in the other two works.However, our approach can maintain a high PEF in the whole target band.

Main Novelties
One of the reasons for the difference between our work and others lies in the choice of evaluator used in BOnSA.Although FOM is derived from the emission rate on average in the physical sense, shorter wavelength contributes less to the performance index, which leads to a shift to longer wavelength.However, PEF is derived according to the underlying mathematics that pays more attention to the entire region, which can overcome the abovementioned shortcomings of FOM.This can be proved by applying BOnSA and PEF to the structures as shown in Figure 5a of our manuscript.By employing the optimization approach proposed in our work, the value of PEF for the {Si 3 N 4 , W, Si 3 N 4 , W} structure increases from 0.537 to 0.677, and PEF of {SiO 2 , Si, W, Si, W} structure increases from 0.374 to 0.642, where the corresponding thicknesses are {76.43, 12.14, 51.31, 180} nm and {108.57, 30, 15.71, 15.71, 200} nm, respectively.Therefore, the results when using PEF and BOnSA are superior in NIR.It is worth mentioning that it takes merely several optimization steps to gain the optimal outcome.
In addition, our results possess merits in terms of multilayer structure design and optimization.In some previous work on multilayer thermal emitters, the elaboration of substructures has been given little attention.Meanwhile, regarding the optimization of material thicknesses in each layer, the majority of the previous work tend to use a single algorithm, which is demonstrated in Table 2. None Enumeration Enumeration and simulation [72] In contrast, the BOnSA algorithm proposed by us is a combination of two effective optimization algorithms.In particular, the epoch of SA is set to be a small finite value, which allows SA to find more local optima in a huge parameter space with reduced time cost.On this basis, BO searches for the global optimum in the neighborhood of local optima.Therefore, the two algorithms can cooperate together building on their respective strengths which leads to the optimization with better performance.

Mechanism Explanation
To explore the mechanism for the optimal structures for thermal emitters, we plot the normalized electric and magnetic field intensities.As illustrated in Figure 8a, the electric field is mainly concentrated in the layer of SiO 2 and MgF 2 .The reason is that SiO 2 on the top layer acts as an anti-reflection substructure to reduce the reflection in the range of target wavelengths.observed inside the aperiodic multilayer structure.Together with the help of the Fabry-Perot etalon, the reflectivity and transmissivity can be sharply declined, and the emissivity can be improved in the visible light area and NIR.The optimal structure of TPV thermal emitters designed via BOnSA in turn helps us in understanding the selective thermal emission mechanism of aperiodic multilayered metamaterials better.Equipped with optimal TPV thermal emitters designed by BOnSA, the TPV systems can obtain high energy conversion efficiencies.Furthermore, the efficiencies in such systems based on near-field TPV can be further improved due to the near-field effects [73,74].

Calculation of Conversion Efficiency
We now calculate the conversion efficiency of TPV systems by the energy conversion relationship.The energy of the incident photon can be obtained by the emission spectrum What is more, it is obvious that a strong magnetic field enhancement exists in Ge as the third layer, as shown in the Figure 8b.The emissivity enhancement originates from localized modes, which are similar to the confined Tamm modes [37,38].Tamm modes of multilayer structure exist inside a photonic band gap, and these modes localized at the surface are usually excited in the periodic system.It is exciting that Tamm modes can be observed inside the aperiodic multilayer structure.Together with the help of the Fabry-Perot etalon, the reflectivity and transmissivity can be sharply declined, and the emissivity can be improved in the visible light area and NIR.The optimal structure of TPV thermal emitters designed via BOnSA in turn helps us in understanding the selective thermal emission mechanism of aperiodic multilayered metamaterials better.Equipped with optimal TPV thermal emitters designed by BOnSA, the TPV systems can obtain high energy conversion efficiencies.Furthermore, the efficiencies in such systems based on near-field TPV can be further improved due to the near-field effects [73,74].

Calculation of Conversion Efficiency
We now calculate the conversion efficiency of TPV systems by the energy conversion relationship.The energy of the incident photon can be obtained by the emission spectrum of the following form: The energy absorbed by the PV cell is given by where λ is the emissivity of the structure, η q,λ is the external quantum efficiency of GaSb PV cell and E bλ is the blackbody radiation intensity.Therefore, the conversion efficiency can be calculated as P abs / P in .With regard to the structures designed for GaSb cell, we have obtained the specific conversion efficiencies at 1200 K, as shown in Table 3.According to Table 3, it is obvious that by means of the optimization approach BOnSA, conversion efficiencies of the corresponding structures have been evidently improved in comparison to the original results in the references.Meanwhile, it can be seen that the selection of materials has a strong impact on the conversion efficiency.
Remark.It is undeniable that there are shortcomings of the algorithm BOnSA.Firstly, the results of the algorithm are partially dependent on the adjustment of the hyperparameters, which turns out to be tricky.Secondly, the emission spectra obtained by BOnSA are still depressed in some regions concerning the wavelength as shown in Figure 5. Though BOnSA can be further improved, it can be extended to solving related problems.For example, in general, in a problem where the physical model and the associated evaluator are known, BOnSA can be applied to seek for optimal parameters.It is promising that BOnSA may also be used to solve black-box problems.

Conclusions
In this paper, we propose a new evaluator PEF to replace FOM and employ BOnSA to carry on the design of TPV thermal emitters for GaSb PV cells as well as InGaAsSb PV cells.Through mathematical modeling and simulation, the optimal material structure for emitters with anti-reflection substructure and Fabry-Perot etalon presents excellent physical characteristics.
Following the rules given in this paper, the aperiodic multilayer structure with better performance can be designed by considering the optical properties of materials in simulations.Specifically, our optimal structure can be described as {SiO 2 , ZnS, Ge, MgF 2 , W, Si, SiO 2 , W}, with the corresponding thickness being {140.0,60.0, 129.6, 423.6, 10.5, 50.6, 8.4, 205.0} nm, whose PEF reaches the value 0.822.
In addition, to demonstrate the robustness and reliability of such an approach, we apply BOnSA to other recent work on the structure design of TPV thermal emitters for GaSb PV cells.The numerical simulation shows that BOnSA leads to a better solution as the value of PEF is increased by more than 20% with less time cost.One can thus conclude that BOnSA has strong adaptability and a low calculation requirement.
Furthermore, with the help of normalized electric and magnetic field intensities, the mechanism of optimal TPV thermal emitters reveals the rationality of the structure from the underlying physics point of view.
In a word, through our optimization, we obtain optimal emitters with excellent optical characteristics, which can provide enlightening guidance for the construction of TPV thermal emitters.
As for the calculation of the electromagnetic physical quantities in the model, one can adopt the transfer-matrix method (TMM) [54,75,76], finite difference time domain method (FDTD) [77] and the rigorous coupled-wave analysis method (RWCA) [78].In this work, we perform TMM, which has the following total transmission moment of electric field amplitude [54]: where E + k and E − k represent the forward and backward traveling waves of the electric field in the k-th layer, which can be vividly seen in Figure A1.As we are concerned with the interface in multilayers, the dynamical matrix of the j-th layer is denoted by D j .The product of dynamical matrices can be expressed as D −1 j D j+1 , and the propagation matrix of the j-th layer is denoted by P j .After applying matrix operations, we can obtain a transfer-matrix for the traveling waves.What is more, the reflectivity and transmissivity of the multilayer structure can be characterized by the parameters of the transfer-matrix.In addition, according to Kirchhoff's law, we can obtain the emission spectrum of thermal emitters.
For the sake of standardization and clarity, we give the notations in this paper.As for an n-layer multilayer structure, d j (1 ≤ j ≤ n) corresponds to the thickness of the j-th layer material.∆λ stands for the target bandwidth, and λ c is the center wavelength of ∆λ.D = [d 1 , d 2 , . . . ,d n ] represents the ensemble of parameters to be op- timized, which spans the parameter space P.

Appendix B
Python was chosen for realizing the algorithm in numerical simulations, which provides a wide variety of open-source and powerful packages such as NumPy [79] and SciPy [80].As for the implement of BO, among several alternatives [81,82], we chose HyperOpt [83], due to its highly customized API and TPE approach.
In terms of hyper-parameters tuning for GaSb cell, the range for center wavelength λ c is configured from 1.0 µm to 1.2 µm.Since when T = 300 K, the band gap wavelength of GaSb cell λ pv = 1.708 µm.Thickness of each layer ranges from 5 nm to 305 nm.λ 1 and λ 2 are set to be 0.4 µm and 5 µm.Bandwidth ∆λ is set to 1.3 µm.
The total epoch of SA is 400.Here, the convergence rate is removed from SA because it would take lengthy calculation time to reach the theoretical limit and it is replaced with η = 0.4, which means the algorithm will terminate if no better PEF appears after epoch × η times.Additionally, the constant κ used to adjust the probability of acceptance is set to 0.4 here.With this setting for BO, BO will search around the local field with a certain size.The ratio of searching field to total space is written as β, whose value is 0.1.For the sake of convenience, all the parameters considered in the optimization are listed in Table A1.

Figure 1 .
Figure 1.Mathematical interpretation of PEF. represents the effective area in the target band. and  represent the areas out of the target band, to the left and right of  respectively.

Figure 1 .
Figure 1.Mathematical interpretation of PEF. S 2 represents the effective area in the target band.S 1 and S 3 represent the areas out of the target band, to the left and right of S 2 respectively.

Figure 2 .
Figure 2. Detailed description of BOnSA flow.In "Generate and Select solutions", the circles represent solutions generated by means of a uniform distribution.The yellow circles represent the solutions that have been selected, while the blue circles are those that have not been selected during this process.The red arrows demonstrate a visual analogy of parameter selection in the optimization algorithm.

Figure 2 .
Figure 2. Detailed description of BOnSA flow.In "Generate and Select solutions", the circles represent solutions generated by means of a uniform distribution.The yellow circles represent the solutions that have been selected, while the blue circles are those that have not been selected during this process.The red arrows demonstrate a visual analogy of parameter selection in the optimization algorithm.

Figure 3 .
Figure 3.The process of determining optimal structure.

Figure 3 .
Figure 3.The process of determining optimal structure.

Figure 4 .
Figure 4. (a) The boxplot of values of PEF corresponding to structures with different hats and repetition of periodic materials.(b) The heatmap of different FP substructures and hat substructures.

Figure 4 .
Figure 4. (a) The boxplot of values of PEF corresponding to structures with different hats and repetition of periodic materials.(b) The heatmap of different FP substructures and hat substructures.

Figure 5 .
Figure 5.The emission spectrum of different optimal structures including the optimization of the structures in [41,70] using BOnSA for GaSb PV cell (a) and for InGaAsSb PV cell (b).

Figure 5 .
Figure 5.The emission spectrum of different optimal structures including the optimization of the structures in [41,70] using BOnSA for GaSb PV cell (a) and for InGaAsSb PV cell (b).

Figure 6 .
Figure 6.The optimization trajectories using three different algorithms with their historical PEF distribution histograms.

Figure 6 .
Figure 6.The optimization trajectories using three different algorithms with their historical PEF distribution histograms.
In terms of the sensitivity analysis in comparison to other work, BOnSA shows strong reliability.

Figure 6 .
Figure 6.The optimization trajectories using three different algorithms with their historical PEF distribution histograms.

Figure 7 .
Figure 7. Analysis of the emission spectrum regarding the optimal structure obtained via BOnSA in comparison with other works.

Figure 7 .
Figure 7. Analysis of the emission spectrum regarding the optimal structure obtained via BOnSA in comparison with other works.

Figure 8 .
Figure 8. (a,b) Contour of normalized electric and magnetic field intensities for our optimal structure.

Figure 8 .
Figure 8. (a,b) Contour of normalized electric and magnetic field intensities for our optimal structure.

Table 1 .
Two common design methods together with relevant research.

Table 2 .
Some previous studies of multilayer thermal emitters.

Table 3 .
Conversion efficiency of the system before and after optimization by BOnSA.

Table A1 .
Physical and algorithmic parameters considered in this paper.