Statistical information enhanced robust design method of optical thin film

: A novel merit function was constructed using the spectral coefficient average error and standard deviation, which can simultaneously optimize the expectation of spectral coefficient error and the envelope of standard deviation. Thus, a multi-objective optimization strategy based on Non-Dominated Sorting Genetic Algorithm and Sequential quadratic programming was proposed. By comparing result of wideband anti-reflection film, cut-off filter and Infrared dual-band filter designed by the conventional algorithm and the new algorithm, the control effect of the new algorithm on sensitivity of film parameters error was verified. The results show that the novel design method has the characteristics of time-efficient calculations and is capable of effectively improving the production yield of the film system, which has practical significance.


Introduction
Since the error between the theoretical and actual values of the optical film structure para-meters is inevitable, except for the spectrum properties, the sensitivity of spectral has been considered as well in the design of optical thin films. However, the sensitivity of spectral coefficients to errors is not fully considered in the conventional optimization algorithms [1][2][3][4][5], which a small film thickness error could also result in great variations in production yield. A robust film system with high production yield should be designed for the practical fabrication.
At present, the single object optimization method [6][7][8][9] has been commonly applied in the robust design of optical thin films, which has only one merit function F in single object (even if this function is the sum of several sections). In the early robust design, the corresponding calculation of the first-order partial derivative of the film [9][10][11] acted as an additional item of the merit function. However, the accuracy of the first order approximation is significantly limited as the true parameter error cannot be considered significantly small [12], and the merit function is commonly located near the local extreme value, where the gradient is close to or equal to zero. Among the robust design methods proposed recently [13][14][15][16][17][18], Monte Carlo method [12,13] adopts the statistical experimental method of sampling simulation and selects the sample average of merit function error as the mathematical expectation. It is capable of achieving effective results for various error distributions, unfortunately, the time taken to get an accurate result required for simulation is often intolerable. The analytical approximation method [15,16] is based on the second-order Taylor expansion of tolerance type merit function. The mathematical expectation of the merit function error is calculated with the fast calculation model of the first and second-order partial derivatives of the film system. It exhibits the advantages of fast calculation speed and high precision. It is applicable to the system which the film parameter errors produced in the fabrication process complies with the zero mean normal distribution. (e.g., quartz crystal oscillator monitoring system). The limitation is that it does not consider the standard deviation of the merit function, therefore it cannot fully represent the degree of error dispersion. In 2018, Shang Qi Kuang et al. proposed a multi-objective optimization robust design method. He constructed two merit functions for broad-band extreme ultraviolet (EUV) thin films [19] to calculate the ideal performance and robustness respectively. With this method, the effectiveness of the method was verified in elevating the fabrication yield of EUV mirrors through two types of design examples of broadband Mo / Si multilayers. The above-mentioned robust design methods are mainly based on the simulation or calculation of the mathematical expectation. However, the effect of the error standard deviation on the practical production is not well considered in these methods.
In this paper we propose new merit functions for robust design from the perspective of fabrication yield, calculating the mathematical expectation and standard deviation of errors, as an attempt to improve the yield about 5%∼15% for different layers. A multi-objective optimization strategy was found to avoid the pre-allocation of weights and the conflict between the original merit function and the additional items due to single objective optimization.

Merit function of robust design
According to the automatic design method of optical thin films, the merit function refers to a quantitative index to determine the matching degree between the theoretical spectral coefficient and the target value, acting as the objective function in the process of optimization. The merit functions applied previously for conventional design are presented below [9]: Under the differences in the range of values of different spectral coefficients (e.g., 0 <T <1, 0 <φ<2π) and their different dimensions, the mentioned merit function can only be exploited to measure one spectral coefficient in many cases. These problems are solved with the following function where ∆ denotes the maximum permissible departure between Q(λ k ) and Q 0 (λ k ). The presented merit function is transformed into the ratio function, unifying the order of magnitude of the merit function and achieving dimensionless. When used, the square tolerance merit function with n = 2 is generally selected. In previous studies, the merit functions employed in robust design can fall to three types [9,15,18].
The merit functions are listed in Table 1, where m denotes the total number of layers; NN represents the total number of Monte Carlo simulations. The first type of merit function calculates the weighted sum of the first partial derivatives of different coatings by conventional merit function as an additional term of the merit function; the second type of merit function obtains the mathematical expectation of the merit function in the presence of error disturbance through Monte Carlo simulation. The third type of merit function is based on the second-order Taylor expansion of the conventional merit function. Under the assumption that the film parameter errors are independent and satisfy the zero mean normal distribution, the analytical solution of the mathematical expectation is taken as the merit function.
The accuracy of approximation of the first-order expansion employed in the first type of merit function is significantly limited, which cannot be applied to most robust designs. Though the accuracy of approximation of the second and third types of merit functions has been significantly improved, and each has its own advantages in the scope of application of the model and the calculation time, the key elements of the two types of merit functions are the mathematical expectation to determine the error of the merit function, but it can't reveal the dispersion degree of the error well. For the mentioned reason, building a merit function with complete statistical information to calculate the mathematical expectation and standard deviation can enhance the robustness of the design film system more effectively.
Given the spectral coefficient error attributed to the structural parameter error of the film, where Q λ denotes the practical spectral coefficient, Q ′ λ represents the theoretical value. Based on the second order Taylor expansion [20], it yields: where δd denotes the thickness error vector; δn represents the refractive index error vector. Due to the Hessian matrix being a symmetric matrix, thus For both sides of the Eq. (4), the mathematical expectation and standard deviation are obtained where E denotes the mathematical expectation, σ represents the standard deviation. Now, F, M and S functions are obtained.
For the practical industrial production, in addition to the expected value of spectral coefficient that determine the actual degree of the error dispersion, there is also the standard deviation of spectral coefficient, which owns great significance for fabrication. When the errors of the film parameters are independent of each other and comply with the normal distribution of the zero mathematical expectation, most of the expectation in Eq. (6) reaches to 0, the expectation of squared error is equated with the standard deviation. the higher order moments σ 2 (δd i 2 ), σ 2 (δn i 2 ) appear in the Eq. (7) and the statistic ( δd i −µ σ ) 2 follows the chi-square distribution χ 2 (1), hence the higher order moments can be expressed as σ 2 (δd i 2 ) = 2σ 4 (δd i ), σ 2 (δn i 2 ) = 2σ 4 (δn i ). Therefore, Eq. (6) and Eq. (7) can be given by Here the spectral coefficient error conforms to the normal distribution. In accordance with the statistical principle, the actual spectrum will be located in the standard deviation envelope with 68.26% probability and in the triple standard deviation envelope with 99.74% probability.
For other distribution of film parameter error, if some specific statistical parameters are known, the mathematical expectation and standard deviation of spectral coefficient increment can be determined by the mentioned formula as well by using Eq. (7), and the statistical distribution parameters under the corresponding distribution form can be subsequently obtained. For independent error distribution, it requires mathematical expectation and standard deviation; for correlated error distribution, it requires other high-order statistics.
According to the expectation and standard deviation of spectral coefficients, we propose two novel merit functions In Eq. (10), the coefficients a can be selected according to different statistical probabilities, which is corresponding to the error envelope of different probabilities. For normally distributed error, the coefficient is supposed to be set to 3, corresponding to the envelope with 99.74% probability. For other distributions, the coefficient a corresponding to a specific probability P is supposed to be determined according to the function P = ∫ a −a ϕ(x)dx, where ϕ(x) denotes the probability density function. The signs can be selected given the type of the designed film system. For instance, if the design target is the transmission spectrum T of antireflection film, we hope the lower limit of T could be higher, therefore the sigh is supposed to be selected as "-", and that of high reflection film is "+". For the case of a beamsplitter with a requirement of T = 50%, both equations using "+" and "-" signs respectively are supposed to be used. Since this design method is based on M.S functions, we name it M-S Robust method. With the fast analytical calculation method of the first and second order partial derivatives of the film parameters [21], the first-order partial derivative matrix and Hesse matrix of the presented formula can be calculated efficiently and accurately. Furthermore, Hesse matrix mostly calculates the continuous product of matrix where M i denotes the matrix of the i-layer film; D i represents the first-order partial derivative matrix of the i-layer film.
Matrix A can be rewritten as where Nested double cycle [21] is adopted to implement the calculation program: construct a four-dimensional matrix B, its elements B(:, :, i, j) = M i M i−1 . . . M j are adopted to record the continuous products of the film matrix, which can be reused in the calculation to implement the fast calculation of Hesse matrix.
In the calculation of spectral coefficients of optical thin films, according to the matrix calculating method, the main part of the algorithm is a multiplication of two-dimensional matrices. An optical film with m layers requires m times multiplication in the calculation of theoretical spectrum. For matrix B, calculation of each element requires only one-time multiplication, however, for the entire B matrix, calculation requires m(m+1)/2 times multiplications, this is due to the symmetry of Hesse matrix. When calculating the first-order and second-order partial derivatives using the B matrix, matrix multiplication calculation of 2 m and 2 m(m+1) times is required respectively. Therefore, the calculation time of f 1 and f 2 is about 3 m times that of theoretical spectrum, which is suitable for most calculation. Since the elements with different indexes i are independent of each other when calculating matrix B and Hessian matrix, parallel computing can be used to accelerate the algorithm.

Multi-objective optimization strategy
Multi-objective optimization strategy in the present study avoids the pre-allocation of weights and the conflict between the original merit function and the additional items [21] due to single objective optimization. The multi-objective optimization strategy proposed in the present study first uses Non-Dominated Sorting Genetic Algorithm (NSGA-II) for initial optimization, and stops when reaching the preset maximum generation number or satisfies the terminating condition. Then, we transform the problem into an equivalent multi-objective programming problem, and the minimum maximum method based on Sequential quadratic programming (SQP) is used for local refined optimization.
NSGA-II [22,23] is a multi-objective optimization algorithm based on the genetic algorithm. Figure 1 shows the flow diagram of NSGA-II algorithm. It introduces elite retention strategy and fast non-dominated sorting strategy to NSGA, which can obtain Pareto frontier with multiple solutions: (1) Initialization. The initial population A 0 with n individuals is created, and the offspring population B 0 is generated by the genetic algorithm (crossover and mutation), and the two populations and population are merged into P t ; (2) Fast non-dominated sorting of A. The crowding degree of the individuals in each nondominated set is calculated. According to the non-dominated relationship and the crowding degree of the individuals, the tournament is conducted among the individuals. The winner is selected to form a novel population, and the next cycle process acts as the parent population; (3) Repeat step (2) until having reached the preset maximum generation number or satisfied the terminating condition.
SQP [24,25] acts as a significantly effective method to process small and medium-sized nonlinear programming. It transforms the problem to be optimized into a series of relatively simple quadratic programming subproblems and solves them. The advantage of this method is that it has local superlinear convergence. To implement multi-objective optimization, the maximum value of the two objective merit functions is optimized with the maximum minimum method to minimize the maximum value.
g Fig. 1. Flow diagram of NSGA-II algorithm

Designs
In our study, we independently write the program of merit function for robust design and use gamultiobj and fminimax functions involved in MATLAB optimization toolbox to implement the above algorithm flow. The function gamultiobj is based on NSGA-II, it reserves the nondominated set partially for the next generation in the elite retention and the default reservation number reaches 0.35. The function fminimax is optimized by the minimum maximum method by complying with SQP.

Design of visible and infrared dual band anti-reflection coatings
This design was one of the topics of the 1988 optical film design competition [26]. Its requirements are elucidated as follows: when designing antireflection coating in 420 -680 nm and 10500 -10700 nm, refractive index for substrate is 2.6 for the 400-700 nm band and 2.4 for the 10000-11000 nm band. In the robust design, for rigorous comparison, the same film layer number (28) is employed as their design, the sampling points number of 420 -680 nm and 10500 -10700 nm reached 27 and 11, the spectral tolerance takes up 0.1%, the lower limit of film thickness is 11 nm, the upper limit is 1.5D, where D denotes the thickness of Wu's design; the population size is 200, the crossover probability is 0.9, and the mutation probability reaches 0.1. In this scale problem, the computation time of each generation is about 0.79s, and the total time is about 6.5 minutes.

Design of robust cut-off filter
We design a cut-off filter with 100% transmittance in 400-650 nm band and 0% transmittance in 700-900 nm band at normal incidence. The refractive index of incident medium and substrate is 1.0 and 1.52, respectively. The high and low refractive index materials are Ti 3 O 5 and SiO 2 , respectively. Considering material dispersion, the fitted Cauchy dispersion coefficient is listed in Table 3. It is assumed that the geometric thickness errors of all coatings comply with the normal distribution of 1% standard deviation and the error in refractive index is not considered.
In the design, the target points of 400 -650 nm and 700 -900 nm are 26 and 21, respectively and the spectral accuracy coefficient is 1%. With the conventional design method, the initial structure of the film system refers to a stack of 21 layers with quarter wavelength L and h, and the reference wavelength reaches 700 nm. Combined with needle method optimization, the design result of 27 (27) is employed, the lower limit of film thickness is 0.5D, the upper limit is 1.5D, where D denotes the thickness of conventional design. The population size is 200 and the crossover probability is 0.9, the mutation probability reaches 0.1.

Merit function
Key elements Advantages Second order Taylor expansion of merit function The analytical calculation has high precision and high speed  Table 4 shows the theoretical reflectivity and transmittance of the two structures. In order to verify the effectiveness of robust design, we utilized Monte Carlo simulation method to simulate the two structures for 20000 times and counted the fabrication yield of each band. The standard deviation of film thickness is set to 1%. Figure 3 shows the theoretical transmission spectrum. Figure 4 and Fig. 5 show the quartiles of the simulation results and the spectral envelope of 95% of the samples. Table 5 shows the fabrication yield of each band.

Design of robust infrared dual-band filter
We designed an infrared dual-band filter with 100% transmittance in 3160-3460nm and 7570-7770nm band and 0% transmittance in 3700-6000nm band at normal incidence The refractive index of incident medium and 1.0. The low refractive index material is ZnS, the high refractive index material and structure are both Ge. The fitted Cauchy dispersion coefficient is listed in Table 6. It is assumed that the geometric thickness errors of all coatings also comply with the normal distribution of 1% standard deviation and the error of refractive index is not considered.  In the design, each band has 50 target points and the spectral accuracy coefficient is 1%. With the conventional design method, the result of 26 layers is achieved according to Liu's [28] (26) was employed, he lower limit of film thickness is 0.5D, the upper limit is 1.5D, where D denotes the thickness of conventional design. The population size is 200 and the crossover probability is 0.85, the mutation probability reaches 0. 15 Table 7 shows the theoretical reflectivity and transmittance of the two structures. In effort to confirm the effectiveness of robust design, we also utilized Monte Carlo simulation method to simulate the two structures for 20000 times and counted the fabrication yield of each band. The standard deviation of film thickness is also set to 1%. Figure 6 and Fig. 7 show the quartiles of the simulation results and the spectral envelope of 95% of the samples. Table 8 shows the fabrication yield of each band. Note, in the above design, the variable range in NSGA-II algorithm was set according to the results of traditional design, i.e. using the traditional design as the reference point. This results in similarities between the two designs. However, the design obtained by the robust design method demonstrated improved performance as compared to traditional designs. And the calculation times of robust designs are listed in Table 9.

Analysis
In the design of dual band antireflection coatings, Fig. 2 illustrates the theoretical spectrum of the film system in Table 1, the mathematical expectation of the spectrum in the presence of error and the envelope of standard deviation (the actual spectrum falls in it with 68.26% probability) under the assumption of the relative error of the film geometric thickness of the normal distribution type with the standard deviation of 1%. As indicated from the comparison of the data in Table 1 and Fig. 2, the mathematical expectation and standard deviation envelope of spectrum in 420 -550 nm band are slightly worse than those of Wu, whereas the range of 550 -680 nm is better than that of Wu. The average value of the whole wavelength range is higher than that of Wu and the average value is elevated by 0.0011% and 0.0018%, respectively. In the range of 10500 -10700 nm, it is better than that of Wu and the average value increased by 0.0031% and 0.0034%, respectively. The spectral sensitivity of the two bands is lower than that of Wu, demonstrating that the results achieved with the proposed method are more robust.
In the design of cut-off filter, Fig. 4 and Fig. 5 show the spectrums of the two design results. It can be seen that the spectral area of interest: cut-off band, is very similar between the two design methods. In particular, the robust design method achieves calculation results in much reduced time with high fabrication yield. However, the ripples seen for robust design method is more pronounced compared to conventional design method, this is a comprise to achieve higher fabrication yield and more efficient calculation while maintaining maximum optimization near cut-off band. The data in Table 5 shows that in the 20000 Monte Carlo simulations, the fabrication yield of robust design results has increased by about 5% on average in the transmission band and about 2% in the stop band.
In the design of infrared dual-band filter, the spectral performance of the conventional design is similar to that of the robust design, however the robust design has better performance in the Monte Carlo simulations. The data in Table 8 shows that in 20000 Monte Carlo simulations, the fabrication yield of robust design results has increased by about 9.5% on average in the transmission band, especially in 3160-3460nm band, the fabrication yield that the transmission of each wavelength is better than 98% is improved by 12%. Figure 6 and Fig. 7 show the data of each wavelength in detail. In the range of 3400-3460nm, the quartile of transmittance increased by 2%. In the range of 7570-7770nm, the quartile of reflectivity decreased by 0.3%. Although the fabrication yield of cut-off band decreased by about 1%, from 99.99% to 98.68%, the fabrication yield of other bands increased significantly. This shows that robust design method has a good application in the design of infrared dual-band filter.
The above analysis shows that the robust design method proposed in this paper can effectively improve the fabrication yield of various thin films. However, the improvement effect is limited for films whose spectrum has a band that is significantly sensitive to errors, such as cut-off filters. For this kind of thin films, the weight setting at different wavelengths is necessary to be reconsidered.

Conclusion
In accordance with the principle of statistics, two types of novel assessment functions including the complete discrete information of spectral coefficients are built. A multi-objective optimization strategy based on NSGA-II and SQP is proposed to simultaneously optimize the mathematical expectation of spectral coefficient error and the standard deviation envelope. In addition, the robustness of the design results is verified by performing the design experiments of visible infrared dual band antireflection film, cut-off filter and infrared dual-band filter. As demonstrated by the results, the dual band antireflection coating developed with this robust design method exhibits stronger robustness and the spectral error envelope fluctuation of the cut-off filter is smaller. Furthermore, the method is easy to implement by software and capable of effectively improving the production yield of the film system with low calculation time and space cost.