Extracting film thickness and optical constants from spectrophotometric data by evolutionary optimization

In this paper, we propose a simple and elegant method to extract the thickness and the optical constants of various films from the reflectance and transmittance spectra in the wavelength range of 350 − 1000 nm. The underlying inverse problem is posed here as an optimization problem. To find unique solutions to this problem, we adopt an evolutionary optimization approach that drives a population of candidate solutions towards the global optimum. An ensemble of Tauc-Lorentz Oscillators (TLOs) and an ensemble of Gaussian Oscillators (GOs), are leveraged to compute the reflectance and transmittance spectra for different candidate thickness values and refractive index profiles. This model-based optimization is solved using two efficient evolutionary algorithms (EAs), namely genetic algorithm (GA) and covariance matrix adaptation evolution strategy (CMAES), such that the resulting spectra simultaneously fit all the given data points in the admissible wavelength range. 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]. 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. For instance, Opti-Char allows an user to set lower and upper limits on the thickness, the refractive index, and the extinction coefficient, and it offers a variety of models like normal, anomal or arbitrary dispersion, and Sellmeier [15]; however, the number of variables and the optimization algorithm used in the estimation process remain unknown to an user [15]. The refractive index and extinction coefficient profiles can be substantially different for various materials [22,23]. For instance, consider some of the metal-oxides; the complex refractive indices of the iron-oxides, Hematite and Magnetite, exhibit more peaks and valleys (optical transitions) than that of the copper oxides [24]. Due to this fact, it is difficult for a supervised learning based estimation approach to adapt to a variety of materials [25]. Recently, a machine learning model has been developed to predict the correlation between spectral data and thickness; however, this approach was suitable for thickness characterization of dielectric materials but not materials with high extinction coefficients (like titanium nitride: TiN) [26]. To attain generality, there arises a need for supervised learning along with knowledge transfer techniques. On the contrary, an evolutionary optimization approach has potential to offer flexible solutions through swarm-based efficient search space exploration [27], while circumventing intense data curation and training requirements involved in traditional supervised learning techniques [26]. Hence, in the present work, we pose the underlying inverse problem as an optimization problem and solve it with evolutionary algorithms [27].
Earlier studies explored optimization-based fitting procedures to extract optical constants using spectrophotometric data. Woollam et al. [28,29] developed a solver that allowed sequential addition of various optical models to minimize fitting errors for extracting refractive index and thickness. A local optimizer, Levenberg-Marquardt (LM) algorithm [30], 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 [31][32][33], whereas single point based gradient descent variant algorithms face challenges due to local entrapment and require good initial guesses to reach the global optimum [30,34]. Moreover, multiple local search procedures from different starting points had been utilized in the Clustering Global Optimization (CGO) algorithm to estimate film thickness and optical constants from spectrophotometric data [15,35]. The estimation accuracy of CGO was good for SiO 2 and Ta 2 O 5 films, although it's efficiency depends on a balanced choice of the starting intervals [35]. EAs have also been exploited to determine the optical parameters of interest [36]. Genetic algorithm (GA) and simulated annealing (SA) have been exploited in the ellipsometric evaluations, where the traditional gradient-based LM method faces difficulty in tackling the related hilly error surfaces [37]. Gao et al. [15] applied GA and SA on an ensemble of Tauc-Lorentz oscillators to extract the real and imaginary parts of complex refractive index 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, a large population size enhances the search space exploration and helps in avoiding local optima at the cost of computational overhead.
The present work aims to determine various thin-films' complex refractive indices and thickness values with the help of EAs, such as GA and CMAES. Recently, covariance matrix adaptation evolution strategy (CMAES) has proved to be reliable in solving deterministic and 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 [27]. Moreover, CMAES exhibited efficient performance on high-dimensional and ill-conditioned optimization problems by utilizing an isotropic (rotation-invariant) evolution path [38]. Further, the performance of the proposed method is validated on metal-oxide and perovskite films, and its computational effectiveness is justified in comparison with the existing methods. The related optical dispersion models and the proposed optimization procedure to solve the addressed inverse problem are described in Physics to Mathematics section. The different types of data used in the present study are illustrated in Data Curation section. The achieved inverse solutions, thickness values and optical constants, of different films are presented in the Results section along with analysis and discussion. Finally, the benefits of the proposed approach are summarized in the Conclusion.

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. The film thickness (d) is a scalar quantity, and the complex refractive index (n + ik) is composed of real and imaginary constituents as functions of wavelengths. The real part of a complex refractive index (n) describes the propagation velocity of the incident light within the film material, and the imaginary part of it (k) concerns about how much of the light gets absorbed in the medium. This study aims to solve the inverse problem of determining {d, n, k} from {R, T}. In the following, we first explain the use of optical oscillator models to emulate complex refractive index and then dive into the problem formulation.

Forward calculations via oscillator models
The complex dielectric function, �(ω) = � 1 (ω) + i � 2 (ω), is analytic in the upper half of the complex ω plane [23,39], where o ¼ c l denotes the frequency of the incident light with λ being its wavelength and c being the speed of light in the 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 [23,40]. Consequently, the imaginary and real parts of it are interconnected by the Kramers-Kronig relation [40,41]. Using this relation [39], one can determine � 1 from � 2 as follows.
where P denotes Cauchy's principle value integral [23]. The complex dielectric function and complex refractive index are related by: � 1 + i� 2 = (n + ik) 2 , which leads to Substituting Eq (2) into Eq (1), unfolds the dispersion relation [23,42] between the real (n) and imaginary (k) parts of a complex refractive index as To accurately evaluate n(ω) from k(ω), the above integration (3) requires to be solved for all frequencies ranging from zero to infinity [23,40]. 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 [42]. Instead, if the functional form of k(ω) is known for all frequencies, then the functional form of n(ω) can be determined elegantly [23]. Therefore, oscillator models are utilized to generate tractable continuous function approximations. In this research, Tauc-Lorentz and Gaussian oscillator models are taken into account, as described below.
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 [15], the imaginary part of the complex dielectric function can be expressed as where E 0i , E g , C i , A i represent the peak transition energy, the band gap energy, the broadening parameter, and the factor involving optical transition matrix elements for the i'th TL oscillator, respectively; ω g and ω 0i are the respective frequencies corresponding to the energies E g and E 0i .
To calculate � 1 from � 2 , let us now recall the Kramers-Kronig relation (1). In practice, the lower limit of the integral in Eq (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 [39]. Note that � 1 (1)>1 is a high frequency dielectric constant to prevent � 1 ! 0 when E < E g . The formulation of these optical functions was first proposed by Forouhi and Bloomer for amorphous semiconductors and insulators [23], and later, extended for crystalline semiconductors and metals [40]. According to the seminal work by Jellison and Modine [39], � 1 (ω) can be derived from Eq (1) by exploiting the continuous approximation (4) of � 2 (ω). These � 1 and � 2 are then used to solve Eq (2) for all ω, so that n and k can be deduced as a function of (3N + 3) parameters, where 3N parameters come from (E 0i , C i , A i ) for i 2 [1, N] and the rest three parameters are: d, ω g , and � 1 (1). The number of decision variables involved in optimizing an TL ensemble model are: (3N + 3).
Ensemble of Gaussian Oscillators: The imaginary r-index profile according to an ensemble of Gaussian oscillators, is given by where  (3) to the continuous approximation (5) of k(ω), which takes shape as where ErfiðÀ W i Þ ¼ 2 ffi ffi p p R À W i 0 exp x 2 dx represents the i th imaginary error function [44]; n(1) = 1 refers to k(1) = 0 and n(1)>1 refers to k(1)6 ¼0 [40,43]. Here, the Gaussian distributions are utilized to retrieve n and k directly instead of deriving them from � 1 and � 2 . An ensemble of Gaussian oscillators (GO) composed of N Gaussian distributions, deals with 3N parameters as A i , μ i , σ i for i 2 [1, N] and two more parameters as n and d. Thus, the total number of decision variables involved in optimizing an GO ensemble model are: (3N + 2).
Note that the above optical constants, n(ω) and k(ω), are expressed as n(λ) and k(λ), respectively, during the numerical implementation. Since o ¼ c l , the order of (n, k) sequences gets reversed when the independent variable is changed from ω to λ. The forward calculation of {R calc (λ), T calc (λ)} for a tuple {d, k(λ), n(λ)} is carried out by the transfer-matrix method [17]. The transfer-matrix method is used to calculate the forward and backward propagating electric fields in smooth homogeneous films, which relies on the superposition of the induced electric fields. The overall transfer matrix is obtained by multiplying a matrix that quantifies the change in field due to the 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 Fig 1. Optimization problem. The optimization problem is defined as In Eq (8), L is the total loss containing differences in the experimentally measured reflectance and transmittance data {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 process computationally efficient, we generate candidate solutions within a search space restricted by an intrinsic correlation [23,42] 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. However, our proposed optimization approach generates the candidate solutions from a search space restricted by the Kramers-Kronig relation and fits all the given spectral data points simultaneously.

Data curation
The proposed methodology is implemented on a diverse data set for two kinds of film materials: (A) metal-oxide, and (B) perovskite (MAPbI 3 ). Different types of data used in the present study are summarized in Table 1, and the details are described below.
Fully-synthetic dataset: There involves two steps in preparing the fully-synthetic data: (i) emulate refractive indices n and k by using a single Tauc-Lorentz oscillator (commonly used for metal oxide materials), and (ii) obtain the reflectance R and transmittance T by using the transfer-matrix method with the (n, k) profiles obtained in step (i) and the user-specified d values, where the wavelength range selected in the UV-Vis-NIR is 350-1000 nm. For this dataset, 1, 116 (n, k) profiles were simulated with a python implementation of the one-oscillator Tauc-Lorentz model in step (i), and 10 random thickness values within a range 20 − 2000 nm were used in step (ii) for obtaining (R, T) spectra.
Semi-synthetic dataset: There involves two steps in preparing the semi-synthetic data: (i) obtain refractive indices n and k from the literature [1] (for perovskite materials); (ii) obtain the reflectance R and transmittance T by using the transfer-matrix method with the (n, k) profiles obtained in step (i) and the user-specified d values, where the wavelength range selected in the UV-Vis-NIR is 350-1000 nm. For this dataset, a set of 18 distinct (n, k) profiles were selected from the literature in step (i), and 15, 640 (R, T) spectra were simulated with a python implementation of the transfer-matrix method by assigning thickness values d in a range 10 − 2010 nm with an increment of 1 nm. Unlike the fully-synthetic dataset, the semi-synthetic dataset only requires one-step simulation, i.e., the simulation of (R, T) spectra using the transfer-matrix method. Experimental dataset: The metal-oxide films, ITO and NiO, are sputtered via physical vapor deposition (PVD) on a glass substrate using an FHR SV-540 in-line sputtering tool. The MAPbI 3 perovskite film is deposited on a glass substrate with two process variables affecting thickness, which are the concentration of the perovskite precursor solution (PbI 2 and MAI with molar ratio of 1:1) and the spin coating speed. We conduct two characterizations: (1) spectrophotometry (UV-Vis) with an Agilent Cary 7000 UV-Vis-NIR Spectrophotometer [45], to obtain the optical reflection and transmission, and (2) profilometry with an KLA Tencor P-16 + Plus Stylus Profiler, to obtain the thickness of the deposited films.
The extinction coefficient of a material is its characteristic/intrinsic property. In the present study, we have considered a variety of materials, including hypothetical and realistic metaloxides and MAPbI 3 perovskites. Each material's extinction coefficient has a specific maximum and minimum. Considering various material samples used in the current data sets (fully-synthetic, semi-synthetic and experimental), the maximum and the minimum of all the extinction coefficients are 2.0 and 0.00, respectively. Note that the extinction coefficients in the semi-synthetic and experimental data sets are more realistic than that of the fully-synthetic data set, which never go beyond the peak value of 1.5.

Results
Implementation Aspects: During the implementation for extracting different films' thickness and optical constants, wavelength (λ) is considered as the independent variable. The 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. Algorithm Selection: 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 fully-synthetic data set. For each sample, three runs are considered with a maximum of 200 iterations and the best result is stored. Table 2 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 2 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. Upon selecting CMAES, it is applied onto the fully-synthetic, the semi-synthetic and the experimental data.

Performance on synthetic data
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). A loss value of 0.05 in Eq (8), 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).
Results using fully-synthetic data: CMAES is applied onto 100 samples randomly drawn from the entire synthetic dataset, and a detailed performance evaluation is presented in Table 3. Overall, the performance metrics, EE d , mEE n and mEE k , suggest less variance in the estimated thickness values and more variance in the estimated refractive indices. The modelbased optimization prove to be effective with TLO-2 and GO-3. The thickness estimates of various samples for the successful occasions, are shown in Fig 2. For three different inverse solutions, the estimated refractive index (n) and extinction coefficient (k) are shown in Fig 3, which are in conformity with their references. In Table 3, 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 an inverse solution where the real or imaginary refractive index estimates have high variances with respect to their reference profiles, although the corresponding thickness estimate is satisfactory. Results using semi-synthetic data: CMAES is applied onto 100 samples randomly drawn from the semi-synthetic dataset, and a detailed performance evaluation is presented in Table 4. To verify the robustness of the proposed method against measurement uncertainty, we inject normally distributed random disturbances into the reflectance and transmittance data of a d = 1350 nm thick perovskite film and run the optimizer to determine the respective inverse solutions. A performance comparison with different levels of measurement noise, i.e. noise m , is shown in Fig 5. As per Fig 5(a), the number of iterations required by the optimizer increases as the noise level increases from 5% to 10% of the (R, T) data. Without any noise, the optimizer takes 238 iterations to reach the global minimum with an estimation error <0.05. With 5% noise m , the optimizer is able to minimize the estimation error below the threshold 0.05 in 367 iterations, however, with 10% noise m , the error is 0.201 after 1500 iterations. The retrieved thickness estimates are: 1352.63, 1350.33, 1352.26 with 0%, 5%, 10% noise m , respectively. Due to the presence of noise m , Fig 5(c) indicates a downward shift in the estimated refractive index and an upward shift in the estimated extinction coefficient for wavelengths below 520 nm. Thus, the proposed method is robust against measurement noise up to certain extent.

Performance on experimental data
Moreover, the model-based optimization with evolutionary algorithms, GA and CMAES, is used to extract the thickness and the optical constants of a 91 nm thick ITO film and a 99 nm thick MAPbI 3 film from the experimentally measured spectra. For each method, the respective optimizers are run ten times and the best results are recorded. The maximum number of allowed iterations is 1000, however, we cut the run if the estimation error goes below 0.05. For two different films, the resulting thickness estimates and the related optimization details are presented in Table 5. The evolution of estimation errors during the optimization process are shown in Fig 6. Even with a small population of 12 candidates, CMAES gives an estimation error of order 10 −2 in 500 iterations and it reduces faster than that of GA, as supported by Fig  6. For the second film in Table 5, the fitting performance and estimated complex refractive index are shown in Fig 7. In Fig 7(a), TLO fits the spectra better than GO throughout the wavelength range of 350 − 1000 nm. Fig 7(c) help visualizing that GO fits the spectra better than TLO in the wavelength range of 350 − 500 nm, with an appropriate broadening [46]. These results justify that the model-based optimization with CMAES can produce accurate d, n(λ), k(λ) estimates using real data with experimental uncertainty. Point-wise vs. Simultaneous Optimization: Earlier research [47] introduced a point-wise unconstrained minimization algorithm (PUMA) that considers finite number of points in a range of wavelengths and makes use of repeated calls to solve the underlying estimation problem. For instance, to extract the thickness and the optical constants of a thin-film from the experimental data, reported in Table 5, the following PUMA calls are made recursively. In this case, a quadratic error of 2.9356e − 2 with respect to 100 distinct wavelengths (points) in 350 − 1000 nm, is attained by PUMA. The estimated thickness (99 nm) is accurate with reference to the original thickness value (99 nm), though the estimated optical constants in Fig 7(d) do not capture fine transitions and they are not in good agreement with the index profiles found by our proposed method. Moreover, there is a lack of instruction on how to select a range for inflection points in PUMA calls and how many calls are sufficient to solve an inverse problem.
Comparison with Existing Package: Furthermore, the proposed method is applied on a literature data (R, T) to extract the optical constants and the thickness of a Si film [15]. In this case, a population of 50 candidates is leveraged and a minimum loss of 0.02 is selected as the stopping criteria for the CMAES algorithm. The optimal solution is presented in Table 6 and the estimated optical constants are shown in Fig 8. The benefits of the proposed method over the existing Optichar Software (OS) and Clustering Global Optimization (CGO) algorithm, are as follows: (i) In OS, a user has many options (normal, arbitrary dispersion, Sellmeier) for a model-based iterative optimization procedure, however, there are no information available on the decision variables and the optimization algorithm. The proposed method offers clarity on the optimization algorithm and the associated variables; (ii) The effectiveness of CGO depends on the choice of initial intervals of the associated decision variables [35]. It is difficult to explore all the minima with a small initial interval, and on the other hand, a broad initial interval might leave a sharp minimum undetected. Unlike CGO, the CMAES algorithm's performance is not so sensitive to its parameter bounds and the choice of bounds can be made without requiring much domain knowledge, as shown in Table 6; (iii) Also, CMAES performs faster than CGO and OS as it minimizes the error function with less evaluations (FE), as reported in Table 6. Table 4. Estimation performance of the adopted oscillator ensemble models applied on 100 samples of type B films, picked randomly from the semi-synthetic dataset. Note that the spectral data are in consistent with linear spectrophotometric measurements.

Analysis and discussion
In Fig 2(a) and 2(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 3. The thickness estimates are slightly dispersed around their original values for films with d > 1.5 μm. At higher film thickness, the optical constants exhibit a nonlinear dependence commonly observed in earlier research findings [49]. In Fig 3(b), n and k estimates exactly match with the original profiles throughout, however, in Figs 3(d), 4(b) and 4(d), there exist slight differences between the original and estimated n, k profiles especially in the low-wavelength regime partially including ultraviolet-visible spectral range. To

PLOS ONE
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 Tables 3 and 4. • 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 the low-wavelength regime, mEE n and mEE k deteriorate for perovskite materials, which could be due to the interaction between multiple inter-band optical transitions [50] or film inhomogeneity [4].
The error metrics reported in Tables 3 and 4 reveal the challenge in achieving a good accuracy of n(λ), k(λ) estimates simultaneously with d. The achieved results justify that the interband transitions for metal-oxide and perovskite materials are well captured by Tauc-Lorentz oscillators. Moreover, the model-based optimization with CMAES can handle measurement uncertainty and extract optical parameters from experimental data, as supported by Table 5 and Fig 7. It is worth noting that a local search algorithm, such as sequential least square   Table 5.
https://doi.org/10.1371/journal.pone.0276555.g006 programming (SLSQP from 'scipy.optimize'), can be applied to the estimates found by an EA to further improve the solution accuracy if needed. Inverse problems involve one-to-many mappings and often they are ill-posed [38,51], therefore, finding exact solutions to such problems is challenging. Further, the difficulty level rises when the inverse solutions are extracted from noisy measurements. The present work demonstrates that the proposed model-based optimization with CMAES does a reasonably good job in extracting accurate inverse solutions. The nature of complex refractive index profiles varies across diverse materials, such as inorganic, organic, and other miscellaneous materials. Thus, the number of oscillator components to be used in an optical ensemble model depends on the type of the concerned material(s), and it is difficult 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. To tackle multiple optical transitions in case of MAPbI 3 perovskite films, the number of oscillator components is selected as higher than that of metal-oxides films. In general, our solver application can be extended to a variety of materials by adjusting the number of composition elements in the optical dispersion 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 an inverse solution with an i7-4600M CPU@2.90 GHz processor. The present approach utilizes the spectral data at just one  Table 6. Performance comparison between the proposed method and the existing methods based on OptiChar software and CGO algorithm. Here, the estimation error is defined according to the reference [15]. Note that the § error is calculated with reference to the spectral data [15] generated using a plot digitizer [48], while considering 651 points in the wavelength range of 400 − 900 nm. incident angle of light, however, in future, additional spectral information (at multiple incident angles) can be considered to alleviate the effect of uncertainty in the experimental measurements.

Conclusion
The proposed model-based optimization with covariance matrix adaptation evolution strategy (CMAES) succeeds in extracting the thickness and the optical constants of metal-oxide and perovskite films from spectrophotometric data using tangible parameters. The employed evolutionary algorithm, CMAES, proves to be efficient in finding globally optimal solutions without spending much function evaluations, and it is robust against measurement uncertainties as well. Overall, the proposed method finds accurate thickness estimates and it can estimate complex refractive indices with multiple transitions.