Extracting Thin-lm Optical Parameters from Spectrophotometric Data by Evolutionary Optimization

Extracting optical parameters from spectrophotometric measurements is a challenging task. In a photometric setup, an unknown thin-ﬁlm is subjected to an incident light beam for a range of admissible wavelengths, which outputs reﬂectance and transmittance spectra. The current work attempts to solve an inverse problem of extracting thin-ﬁlm thickness and complex refractive index from reﬂectance and transmittance spectra for an incident angle of light. The ﬁlm thickness is a scalar quantity, and the complex refractive index is composed of real and imaginary parts as functions of wavelengths. We leverage evolutionary optimization techniques to solve the underlying inverse problem, which determines the desired parameters associated with two optical dispersion models: ensemble of Tauc-Lorentz (TL) and ensemble of Gaussian oscillators, such that the generated spectra accurately ﬁt the input data. The optimal parameters involved in the adopted models are determined using efﬁcient evolutionary algorithms (EAs). Numerical results validate the effectiveness of the proposed approach in estimating the optical parameters of interest.


Introduction
The determination of a thin-film's thickness and complex refractive index over a broad spectral range is carried out using either ellipsometric or spectrophotometric analysis 1,2 . Ellipsometry measures how the polarization of a light beam changes on reflection from a surface. Spectrophotometry measures reflectance and/or transmittance of light through thin-films and substrates as function(s) of wavelength. Reflectance and/or transmittance spectral information are further fed to a numerical solver for extracting thickness and refractive index 3,4 . Thickness influences micro to macro-scale physics, including from electron mobility in thin-film transistors 5 to operating characteristics of thin-film gas sensors 6 . Refractive index based detectors play a significant role in characterizing the changes of several chemical and biological measurands, like gas concentration and water pH values 7 . Complex refractive index profiles are crucial in designing appropriate polymers to produce effective lenses and ultraviolet (UV)-absorbing coatings 8 . Thickness variations by adjusting the catalyst concentration or changing the heat-treatment process can be leveraged to shift the reflection band of a film from narrow to broad wavelength region, which is beneficial in fabricating dielectric reflectors for solar cells and band pass filters for optical instruments 9 . For compounded materials, an accurate estimation of a film's optical parameters is required to optimize the deposition process.
Computationally, it is quite nontrivial to achieve accurate estimates of all the optical parameters together 10 . A method that gives the most accurate solution for thickness, does not usually give the most accurate solution for complex refractive index 11,12 . Moreover, different methods perform well in different spectral ranges 12,13 . Two inherent difficulties in an inverse photometric problem are: (a) missing information induced by experimental uncertainty or film imperfections, and (b) multiple solutions due to the fact that the measured intensities might be compatible with other plausible combinations of optical constants 14 .
To dissolve the ambiguity in parameter extraction, previous researchers 12 exploited reflectance and transmittance spectral measurements at two distinct incident angles of light. Thus, there emerges a need of developing efficient methods for the precise determination of not only thickness but also complex refractive index. Several commercial packages of thin-film optical software, such as TFCalc, Filmwizard, Optilayer, and Essential McLeod, provide advanced modelling for optical dispersion laws and use fitting methods to determine the thickness and complex refractive index 11 . However, these packages cannot always find satisfactory solutions, especially when the initial candidate is far away from the desired solution in case of a strongly absorbing material like silicon in the visible range 11 . A recent development, RefDex 4, 15 provided an interactive fitting procedure to extract the complex refractive index of a film with known thickness from reflectance and transmittance data. This inverse problem is not unique and there exist multiple solutions that minimize the same loss function 4 , which mandates imposing physical constraints to determine meaningful solutions. The major concern is that the above-mentioned software-packages do not offer much flexibility on the choice of parameters, such as the number of oscillators and the bounds on model coefficients, involved in their respective optimization procedures.
Refractive index and extinction coefficient profiles can be substantially different for various materials 16 . For instance, consider some of the metal oxides; the complex refractive indices of the iron oxides, Hematite and Magnetite, exhibit more peaks and valleys than that of the copper oxides 17 . Due to this fact, a supervised learning approach confronts challenges in predicting refractive indices for a broad range of materials 18 . A supervised method that is trained with inorganic material data, might not give satisfactory prediction performance in the case of organic materials 19 . Thus, there arises a need for transfer learning in combination with supervised learning. On the contrary, an evolutionary optimization approach has the potential to bypass the requirement of intense training with labelled data. That is why, in the present work, we pose the underlying inverse problem as an optimization problem and solve it using evolutionary algorithms.
Earlier studies explored optimization-based fitting procedures to extract optical parameters from spectrophotometric data. Woollam et al. 20,21 proposed a tool that allowed sequential addition of various optical models and simultaneously minimized fitting errors to determine refractive index and thickness. A local optimizer, Levenberg-Marquardt (LM) algorithm 22 , was utilized in their approach, which gradually improved the solution accuracy. The optimizer starts from an initial guess and moves it to a feasible local minimum until another model addition becomes necessary for further reduction in the fitting error. Swarm-based evolutionary algorithms (EAs) are robust and good at finding the global optimum even for high-dimensional problems 23,24 , whereas single point based gradient descent variant algorithms face challenges due to local entrapment and require good initial guesses to reach the global optimum 22,25 . EAs have also been exploited to determine the optical parameters of interest 26 . Genetic algorithm (GA) and simulated annealing (SA) have been exploited in ellipsometric evaluations, where the traditional gradient-based LM method faces difficulty in tackling the related hilly error surfaces 27 . Gao et al. 11 applied GA and SA on an ensemble of Tauc-Lorentz oscillators to extract the refractive index real and imaginary parts of a film with known thickness from reflectance and transmittance spectra. However, their approach required a large population to determine the desired inverse solutions. In EAs, usually a small population size saves the cost of computing functions whereas a large population size enhances the search space exploration and helps in avoiding local optima at the cost of computational overhead. In recent years, covariance matrix adaptation evolution strategy (CMAES) has proved to be reliable in solving deterministic or stochastic global optimization problems even with a small population size due to the attributes like step-size adaptation, noise effect reduction, and invariance under coordinate systems 28 . The present work aims at determining a thin-film's complex refractive index and thickness with the help of EAs, such as GA and CMAES.
In this paper, we adopt EAs to solve the inverse problem of extracting a thin-film's optical parameters from spectrophotometric data of reflectance and transmittance spectra at just one incident angle of light. Further, the current solution approach can be applied to a variety of materials. The following first describes the related optical models along with the proposed optimization methodology to solve the underlying inverse problem, and then presents the numerical results obtained by the same. Finally, we discuss the major traits of the achieved inverse solutions and conclude the current work.

Physics to Mathematics
In our spectrophotometric analysis, the optical parameters to be extracted are: thickness (d), refractive index (n), and extinction coefficient (k), and the measured data are: reflectance (R) and transmittance (T ) spectra. A thin-film is situated in between air and substrate mediums in a photometric set-up, as shown by Figure 1. The film material has real and imaginary constituents in its complex refractive index profile. The real part (or refractive index) describes the propagation velocity of the incident light within the material, and the imaginary part (or extinction coefficient) concerns about how much of it gets absorbed in the medium. This study aims to solve the inverse problem of determining {d, n, k} from {R, T }. The following describes the use of optical oscillator models to emulate complex refractive index (n + ik) profiles required in the solution procedure.

Forward Calculations via Oscillator Models
The complex dielectric function, ε(ω) = ε 1 (ω) + i ε 2 (ω), is analytic in the upper half of the complex ω plane 16,29 , where ω = c λ denotes the frequency of the incident light with λ being its wavelength and c being the speed of light in air. The associated photon energy is represented by E = hω; h is Planck's constant. The analytic behavior of ε(ω) stems from the  principle of causality 16,30 . As a consequence, the imaginary and real parts of it are interconnected by the Kramers-Kronig relation 30 . Using this relation 29 , one can determine ε 1 from ε 2 as follows.
where P denotes Cauchy's principle value integral 16 . The complex dielectric function and complex refractive index are related by: ε 1 + iε 2 = (n + ik) 2 , which leads to By substituting Equation (2) into Equation (1), the dispersion relation 16,31 between the real (n) and imaginary (k) parts of a complex refractive index can be unravelled as To accurately evaluate n(ω) from k(ω), the above integration (3) requires to be solved for all frequencies, ranging from zero to infinity 16,30 . However, from an experimental perspective, it is only feasible to address a finite range of frequencies. A numerical integration with such a limited spectral range gives erroneous results 31 . Instead, if the functional form of k(ω) is known for all frequencies, then the functional form of n(ω) can be determined elegantly 16 . Therefore, oscillator models are utilized to generate efficient continuous function approximations. In this research, we take Tauc-Lorentz and Gaussian oscillator models into consideration, as described briefly in the following. Ensemble of Tauc-Lorentz Oscillators: Consider a material with the complex dielectric function: ε = ε 1 + iε 2 . According to an ensemble of N Tauc-Lorentz (TL) oscillators 11 , the imaginary part of the complex dielectric function can be expressed as where ω 0 i , ω g ,C i , A i represent the peak transition frequency, the band gap frequency, the broadening parameter, and the factor involving optical transition matrix elements for the i'th TL oscillator, respectively. To calculate ε 1 from ε 2 , let us now recall the Kramers-Kronig relation (1). In practice, the lower limit of the integral in Equation (1) is chosen as ω g instead of 0 because the Tauc-Lorentz model requires ε 2 to be zero for photon energies below the band gap 29 . Note that ε 1 (∞) > 1 is a high frequency dielectric constant to prevent ε 1 → 0 when ω < ω g .

3/12
The formulation of these optical functions was first proposed by Forouhi and Bloomer for amorphous semiconductors and insulators 16 , and later, extended for crystalline semiconductors and metals 30 . According to the seminal work by Jellison and Modine 29 , ε 1 (ω) can be derived from Equation (1) by exploiting the continuous approximation (4) of ε 2 (ω). These ε 1 and ε 2 are then used to solve Equation (2) for all ω, so that n and k can be deduced as a function of (3N + 3) parameters, where 3N parameters come from (ω 0 i ,C i , A i ) for i ∈ [1, N] and the rest three parameters are d, ω g , and ε 1 (∞). The number of decision variables involved in optimizing a TL ensemble model are (3N + 3). Ensemble of Gaussian Oscillators: The imaginary r-index profile according to a Gaussian oscillator model, is given by where, N is the number of Gaussian functions to be mixed; the parameters associated with the i'th Gaussian function are A i , µ i , σ i denoting its amplitude, mean and variance, respectively 32 . Next, n(ω) is determined by applying the Kramers-Kronig integration (3) to the continuous approximation (5) of k(ω), which takes shape as where n(∞) = 1 refers to k(∞) = 0 and n(∞) > 1 refers to k(∞) = 0 30, 32 . Here, Gaussian distributions are utilized to retrieve n and k directly instead of deriving from ε 1 and ε 2 . An ensemble of Gaussian oscillators (GO) composed of N Gaussian distributions, deals with 3N parameters as and two more parameters asn and d. Thus, the total number of decision variables involved in optimizing a GO ensemble model are (3N + 2). Note that the above optical characteristics n(ω) and k(ω) are expressed as n(λ ) and k(λ ), respectively, during the numerical implementation. The order of (n, k) sequences gets reversed, ascending to descending, when the independent variable ω changes to λ since ω = c λ . The forward calculation of {R calc (λ ), T calc (λ )} for a tuple {d, k(λ ), n(λ )} is carried out by the transfer-matrix method 13 . The transfer-matrix method is used to calculate the forward and backward propagating electric fields in smooth homogeneous films, and it relies on the superposition of electric fields. The overall transfer matrix is obtained by multiplying a matrix that quantifies the change in field due to light waves propagating through an interface (air-to-film) with another matrix that quantifies the change in field due to the same waves propagating within a layer (film).

Inverse Problem Formulation
We now explain how the underlying inverse problem is posed as an optimization problem. An overview of the associated forward and inverse processes is depicted in Figure 2. Optimization Problem: The optimization problem is defined as In Equation (7), L is the total loss containing differences in the experimentally measured reflectance and transmittance {R meas , T meas } from their theoretically calculated values {R calc , T calc }, for all wavelengths. The admissible range of the wavelength, i.e. {λ lb , λ ub }, depends on the experimental infrastructure. The decision variables associated with the optimization problem, are: thickness d, real r-index profile n(λ ), and imaginary r-index profile k(λ ). Thickness candidates are directly passed to the solver, whereas n and k candidates are passed in terms of the oscillator model parameters. The evaluation of {R calc (λ ), T calc (λ )} for a candidate solution {d, n(λ ), k(λ )} has been discussed in the preceding section. In order to make the optimization approach computationally efficient, we generate candidate solutions within a search space restricted by an intrinsic correlation 16,31 between the optical parameters of interest. The thickness candidates are chosen from a reasonable range of scalars. The imaginary and real refractive index candidate profiles are provided by the adopted oscillator models. A formulation over a range of wavelengths is computationally more tractable than a formulation at distinct wavelengths. In a discrete approach, the search space increases with the number of wavelengths and it is cumbersome to determine a physically meaningful refractive index (and/or extinction coefficient) profile out of every solution point at each wavelength, while satisfying the associated constraints. Whereas in our optimization approach, the parameter search space gets restricted by the Kramers-Kronig relation during candidate generation.

Data Curation
The proposed methodology is implemented on a diverse data set containing two kinds of thin-film materials: (A) metal oxide thin-films, and (B) perovskite thin-films.

Inverse
Refractive index (n) Extinction coefficient (k) . . . We obtain the synthetic data for two classes of the datasets, the Tauc-Lorentz (TL) dataset whose refractive index follows a single Tauc-Lorentz oscillator (commonly used for metal oxide materials), and the perovskite data set whose refractive index is from perovskite materials 33 . There are two steps involved in preparing the synthetic data: 1) obtain refractive indices n and k either from the simulation (for TL dataset) or from the literature 33 (for perovskite dataset); 2) obtain the reflectance R, and transmittance T by using the transfer-matrix method with the (n, k) profiles obtained in step (1) and user-specified d values, where the wavelength range selected in the UV-Vis-NIR is 350 -1000 nm. For TL dataset, Step (1): 1, 116 (n, k) spectra were simulated with a python implementation of a one-oscillator Tauc-Lorentz model.
Step (2): 15, 640 (R, T ) spectra were simulated with a python implementation of the transfer-matrix method by assigning 10 random thickness values d in the range of 20 − 2000 nm throughout the dataset. Unlike the TL dataset, the perovskite dataset only requires one-step simulation, i.e., the simulation of (R, T ) spectra using the transfer-matrix method, while a variety of 18 (n, k) profiles are picked from the literature 33 .
To generate an experimental data, the MAPbI 3 perovskite film is deposited on glass substrate with two process variables affecting thickness, which are the concentration (M = 1.25) of the perovskite precursor solution (PbI 2 and MAI with molar ratio of 1:1) and the spin coating speed (r.p.m.= 3000). We conduct two characterizations: (1) spectrophotometry (UV-Vis), with an Agilent Cary 7000 UV-Vis-NIR Spectrophotometer 34 , to obtain the optical reflection and transmission (2) profilometry, with a KLA Tencor P-16 + Plus Stylus Profiler, to obtain thickness of deposited films.

Numerical Results
Implementation Aspects: During the implementation for extracting thin-film optical characteristics, wavelength (λ ) is considered as the independent variable. Reflectance and transmittance spectra are measured for a wavelength range of 350 to 1000 nm. The inverse problem of determining d, n(λ ), k(λ ) from R, T is not unique. So, there is a risk that an optimization algorithm unravels solutions that do not make sense physically. To mitigate this issue, we pose bounds on the decision variables involved in TLO and GO models such that negative n and k are always discouraged during the optimization process. The decision variables associated with the TLO ensemble model, i.e. A i , ω 0 i ,C i , ω g , ε ∞ , d, are selected from a bounded space: Choice of Algorithms: We apply two EAs, namely Genetic Algorithm (GA) and Covariance Matrix Adaptation Evolution Strategy (CMAES), for optimizing the parameters associated with the adopted optical dispersion models: TLO ensemble and GO ensemble. To validate the performance of CMAES and GA on TLO as well as GO ensembles, 20 samples are drawn randomly from the synthetic data set. For each sample, three runs are considered with a maximum of 200 iterations and the best result is stored. Table 1 presents a performance comparison between the employed local and global optimization algorithms, i.e. LM vs. GA vs. CMAES. The average loss is drastically higher in case of the local optimization algorithm than the global optimization algorithms because the former drives only a single point that may easily get trapped by local optima whereas the latter drives a population of candidates to efficiently explore the search space while avoiding local traps. Therefore, the performance of LM very much depends on the initial condition; a better starting point leads to a lower loss value. On the contrary, the performance of EAs like GA or CMAES is less sensitive to the initial population. Table 1 shows that the average loss achieved by CMAES is the lowest. This study bolsters the use of CMAES over GA in the underlying inverse problem. Performance Evaluation: CMAES is applied onto 100 samples of type A and 100 samples of type B thin-films, randomly drawn from the synthetic data set. For each sample, five runs are considered with a maximum of 1500 iterations and the best result is stored. We cut the run if CMAES reaches a loss value below 0.05 (stopping criteria). Note that a loss value of 0.05 in Equation (7), refers to a mean square error of 0.05/651 = 7.68E − 5 for (1000 − 350) + 1 = 651 wavelengths in the given range. A successful occasion is counted if the thickness estimation error is below 10% along with a minimum loss of ≤ 0.17 (saving criteria). Table 2 reports all the associated metrics, such as: EE d denotes the R 2 -score between the original thickness and estimated thickness values, mEE n denotes the median of the sequence of R 2 -scores between the original and estimated n(λ ) arrays, mEE k denotes the median of the sequence of R 2 -scores between the original and estimated k(λ ) arrays, mEE R denotes the mean of the sequence of R 2 -scores between the original and estimated R(λ ) arrays, mEE T denotes the mean of the sequence of R 2 -scores between the original and estimated T (λ ) arrays, SR denotes the success rate = number of successful occasions (n ss )/ total occasions (n s ), and mFE is the average number of function evaluations: (1/n ss ) * ∑ n ss i=1 FE i . In optical parameter estimation, optimizing the TLO ensemble model proves to be more efficient than optimizing the GO ensemble. For successful occasions, optimizing TLO ensemble requires less function evaluations than optimizing GM models.  Table 2. Estimation performance of the adopted oscillator ensemble models applied on 100 samples of type A and B films randomly picked from the total synthetic data set.
Moreover, we applied local and global optimization algorithms, i.e. LM, GA and CMAES, on TLO and GO models to extract a perovskite thin-film's (unknown n, k) optical characteristics from the experimentally measured spectrophotometric data, (R, T ). For each method, the respective optimizers are run ten times and the best results are recorded. The maximum number of iterations is set as 1000, however, we cut the run if the loss function value goes below 0.05. The corresponding thickness estimates and optimization loss outcomes are listed in Table 3. Figure 6 exhibits the fitting performance of the CMAES algorithm applied on two optical models, along with the extracted complex refractive index profiles. Interestingly, in   Table 2: (a) thickness estimation for type A thin-films, (b) thickness estimation for type B thin-films. this case, the optimized GO model gives better estimates of d, n(λ ), k(λ ) than the optimized TLO model.  Table 3. Estimation performance of various methods on the experimental data.
Inference from Results: Figure 3 shows all the thickness estimates of various samples for the successful occasions. According to Figure 3(a) and 3(b), the optimized TLO ensemble model finds more number of accurate thickness estimates (within 10 % of the original values) than the optimized GO ensemble model, which is also supported by Table 2. The performance of the adopted models in estimating refractive index (n) and extinction coefficient (k), can be visualized by the examples shown in Figures 4 and 5. It comes out that capturing n and k profiles simultaneously with thickness estimation is a challenging task. Therefore, mEE n and mEE k in Table 2 are slightly worse even though EE d is quite good. Among Tauc-Lorentz oscillator ensembles, TLO-2 and TLO-4 prove to be efficient and among Gaussian oscillator ensembles, GO-3 and GO-5 are effective. Overall, TLO ensemble optimization serves better than GO ensemble optimization. In Table 2, we mention the medians mEE n and mEE k as the associated R 2 -score sequences contain outliers that might skew the average of scores. An outlier refers to a solution where the real or imaginary refractive index estimates are inaccurate, although the corresponding thickness estimate is satisfactory. In Figure 4(b), n and k estimates exactly match with the original profiles throughout, however, in Figures 4(d), 5(b) and 5(d), there exist differences between the original and estimated n, k profiles especially in low-wavelength regime, partially including UV and visible spectral ranges. To further inspect this quantitatively, we split the entire wavelength range into two sectors and highlight the estimation performance of the well fitted models as per Table 2.
• TLO-2 for type A films: For a wavelength range of 350 − 500 nm, the estimation metrics are mEE n = 0.41617, mEE k = 0.98999; and for a wavelength range of 500 − 1000 nm, the estimation metrics are mEE n = 0.83078, mEE k = 0.99870.
• TLO-4 for type B films: For a wavelength range of 350 −500 nm, the estimation metrics are mEE n = −0.47954, mEE k = 0.72331; and for a wavelength range of 500 − 1000 nm, the estimation metrics are mEE n = 0.97758, mEE k = 0.95517.
The above reported results reveal that the fitting error in k(λ ) gets amplified into n(λ ) as it propagates through the Kramers-Kronig integration. In low-wavelength regime, mEE n and mEE k deteriorate seriously for perovskite materials. A possible reason is the presence of intra-band absorption giving rise to ε 2 for photon energies below the band gap, which can not be captured well by the Tauc-Lorentz model approximation due to an underlying assumption 29 : ε 2 (ω) = 0 for ω < ω g .
The numerical results justify that the TLO ensemble model performs well for the synthetic data of both type A and B thin-films. The estimation error metrics reported in Table 2   estimates simultaneously with d. The GO model can not be discarded based on its performance on the synthetic data in Table 2 because the same produces accurate estimates from the experimental data, as supported by Table 3 and Figure 6. It is worth noting that a local search algorithm, such as sequential least square programming (SLSQP from 'scipy.optimize'), can be applied to the estimates found by an EA to further improve the solution accuracy if needed.

Discussion on Results
Inverse problems involve one-to-many mappings and often they are ill-posed 35 , therefore, finding exact solutions to such problems is challenging. Further, the difficulty level rises when inverse solutions are determined from noisy measurements. The present work demonstrates that the proposed optimization approach in combination with EAs does a reasonably good job in extracting accurate inverse solutions. The optical characteristics including the nature of complex refractive index profiles vary across diverse materials, such as inorganic, organic, glass and other miscellaneous materials. The number of oscillator components to be used in an ensemble model depends on the type of the concerned material. It is not trivial to come up with an universal rule of selecting the same; the current choice is subjective to two types of film materials: metal oxide and perovskite. In the present study, the number of oscillator components used to tackle type B thin-films are higher than that of type A thin-films because perovskite materials' refractive index profiles are more complex than metal oxides. In general, our solver application can be extended to a variety of materials by adjusting the number of composition elements in the optical models. The proposed evolutionary optimization approach does not require any prior learning or memory-based mapping. The successful runs take only about a few minutes to find a solution in an i7-4600M CPU@2.90 GHz processor. The present approach utilizes the spectral data at just one incident angle of light, however, in future, more spectral information (at multiple incident angles) can be added to alleviate the effect of uncertainty in experimental measurements. Also, there remains a scope of devising advanced techniques to tackle surface roughness and layer inhomogeneity.

Conclusion
The proposed optimization methodology using an evolutionary algorithm elegantly solves the inverse problem of determining optical parameters from spectrophotometric data. We achieve a quite satisfactory level of accuracy in estimating thin-films' thickness values. The employed evolutionary algorithm, CMAES, proves to be efficient in finding global solutions without spending much function evaluations. Interestingly, our proposed evolutionary optimization approach is able to extract the optical characteristics of the perovskite thin-films with both real and imaginary parts of the complex refractive index. In overall parameter estimation, the optimized Tauc-Lorentz ensemble performs well and the optimized Gaussian ensemble performs descent for metal oxide as well as perovskite thin-films. To this end, our proposed method can extract various thin-films' optical characteristics, i.e. thickness, real and imaginary refractive index, simultaneously from spectrophotometric data over a broad range of wavelengths, and it offers generality to be applicable for diverse materials.