Measurement of layer thicknesses with an improved optimization method for depolarizing Mueller matrices

There are some commonly-used optimization techniques for the analysis of measured data in spectroscopic Mueller matrix ellipsometry (MME) used, for example, to calculate the layer thicknesses of samples under test. Concentrating on the metrological aspects of MME, we identified a non-optimal treatment of depolarization in all these techniques. We therefore recently developed an improved optimization method to adequately take depolarization in MME into account. In a further step, we also included statistical measurement noise and derived a likelihood function, which enabled us to apply both the maximum likelihood method and Bayesian statistics as well as the Bayesian information criterion for data evaluation. In this paper we concentrate on the application of this new method to measurements of SiO2-layer thicknesses on silicon. With a state-of-the-art SENTECH SENresearch 4.0 Mueller ellipsometer, we measured standard samples of different SiO2-layer thicknesses, whose calibrated thicknesses were between about 6 nm and 1000 nm. The MME results were compared to the calibration data. For all samples, an SiO2-SiO double-layer model turned out to be optimal. The measured total oxide layer thicknesses matched excellently with the calibration values, within the estimated range of uncertainties. All the results are presented here. This is the first comparison with traceable reference measurements demonstrating the validity of our novel MME analysis method.


Introduction
Spectroscopic ellipsometry (SE) is an indirect measuring technique widely used to measure the dimensional or optical Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. parameters of surface structures and surface layers. Our goal is a traceable and reliable SE for dimensional metrology of material parameters, layer systems and nanostructures on surfaces. To achieve this, a full uncertainty consideration is required, which must address uncertainty contributions arising from hardware, measurement and analysis.
A schematic for SE is shown in figure 1. To derive the parameters of the surface being measured, an inverse problem has to be solved, which is achieved by a fitting procedure. Within this procedure, all the spectral measurements for each angle of incidence (AOI) are combined into one merit function, which is optimized by comparing simulated and measured data. For the simulation of the optical response of unstructured surfaces, the Fresnel equations can be used. However, in the case of structured surfaces, Maxwell's equations have to be solved using a so-called Maxwell solver.
If there is no depolarization, i.e. for fully deterministic cases, the well-known Jones formalism can be used and simulated by the above-mentioned Fresnel or Maxwell solvers. However, if depolarization occurs in ellipsometric measurements, the Mueller-Stokes formalism is required to describe the depolarization, i.e. MME has to be used. The incoming and outgoing light is described by the Stokes vectors S in and S out and the polarizing action of the sample (including depolarization) is described by a 4 × 4 Mueller matrix M. This leads to the following equation: The full rigorous simulation of depolarizing cases would offer the best possible accuracy, but it incurs a very high computational cost. In practice, this usually requires too much effort or is even impossible. So, there is an important question: How to treat depolarization in the optimization process? To answer this question, we have recently developed a new analysis method for MME, which for the first time realistically takes into account the uncertainty contributions introduced by stochastic influences, such as (sample-induced) depolarization and measurement noise [1,2]. We present here the first application of our novel analysis method on MME layer thickness measurements. For this purpose, we measured a set of unstructured layer thickness standards with nominal layer thicknesses of 6 nm, 70 nm, 160 nm, 380 nm and 1000 nm, which have been calibrated by PTB using SE. However, the SE was itself calibrated using another identically-manufactured set of layer thickness standards with the same design, which were calibrated by PTB using traceable x-ray reflectometry (XRR) measurements [3][4][5]. Thus, the SE used for the calibration process is merely used as a comparator and the traceability of the calibration is essentially provided by XRR. The XRR calibration measurements were performed at the synchrotron radiation source BESSY II in Berlin.
In contrast, the analysis of our MME measurement data provides reference-free physical results and, at the least a reasonable estimation of the stochastic uncertainty contributions. Therefore, the comparison with the calibration values provides an excellent test of the validity of our new analysis method.
In the next chapter, we concisely discuss a major challenge of common merit functions and briefly describe our new merit function. The third chapter presents and discusses the application of our new method to the measured data of the layer thickness standards. Finally, we end with a short conclusion and outlook.

New merit function in a nutshell
It has been shown using simulated data of an SiO 2 layer on an Si substrate, in the ideal case of no experimental errors and rawdata noise (and hence with an uncertainty resulting from thickness variations within the spot size alone), that even for small amounts of depolarization, common merit functions lead to incorrect measurement uncertainty estimations that are obviously much too small [1]. Therefore, from our point of view, a merit function was needed that is able to return the preset uncertainties and hence, for the first time, also reliable uncertainty estimations in the case of real measured data. For this, we started with the well-known definition of a multivariate likelihood function (see equation (15) in [1]), which reaches its maximum for the best-fitting parameters. The related χ 2 -function reaches its minimum for the same parameters and was our starting point for further developments towards application in ellipsometry. Firstly, we integrated measured depolarization as a systematic part of the measurement process. To do so, we used Cloude's covariance matrix H [6]. After some computation, and benefitting from previous work in [7], where Cloude's covariance matrix H was already pointed out to be interpretable as a statistical quantity describing depolarization, we derived a merit function, which implicitly uses depolarization as a weighting for each residual contribution and also takes measurement noise into account [1]: Here j sim is a simulated Jones-element vector. H CF is the filtered H matrix using the so-called Cloude filter [6]. This merit function is used in the next chapter to analyse real data from the layer thickness standards.
A first test using the above-mentioned simulated measurement data again confirmed the performance of that function. In a further simulation step the measurement noise was additionally considered, which results in increasing uncertainties in the case of increasing noise levels and therefore again suggests plausible results. It is clear that, in reality, many more factors influence the measured data. However, from our point of view it is an appropriate method to make a first plausibility check of the new merit function, to establish whether it recovers preset uncertainties. If the new method fails even in such an ideal case, it will also fail in more complicated cases. For further details, please refer to [1].
Beside the quantities already discussed, i.e. (sampleinduced) depolarization and measurement noise, for a full uncertainty budget the determination of further quantities such as device-related contributions would be required. Moreover, for an even more precise description, the role of the degree of coherence during measurements (e.g. [8][9][10][11]) must be addressed as well. The implementation of a formalism representing the theoretically modeled values by M or H in our rigorous calculations, is a topic of current research.

Data analysis
We performed measurements on layer thickness standards with a Sentech SENresearch 4.0 Mueller ellipsometer (see figure 2), measuring in the spectral range of 190 nm to 1000 nm at N = 940 discrete wavelengths at an AOI of 70 • using a beam diameter of 2 mm.
For the data analysis we tested both a single (model 1: height h SiO2 ) and a double layer model (model 2: heights h SiO2 and h SiO ) via a maximum likelihood analysis. The dielectric parameters for Si and SiO 2 are taken from Herzinger et al [12] and for SiO from Palik [13]. The main results are shown in table 1. We used both the normalized χ 2 /N as well as a reduced Bayesian information criterion BIC red [1] derived from the Bayesian information criterion BIC, to select the best-fit model. For nearly all of the measured samples, we obtained the result that the double-layer model SiO 2 -SiO on Si is the optimal choice (bold values). This was confirmed in this example from applying both χ 2 /N as well as the BIC red criterion. Only for the nominal 6 nm thick layer did both models fit equally well. For the double layer model (model 2), very high normalized anticorrelations of the two fit parameters h SiO2 and h SiO were observed (see the last column of table 1) and hence this confirms our results in [1] and from the existing literature [14,15]. This is exemplarily shown in the correlation diagrams for the measurements of the nominal 6 nm and 380 nm layers, respectively (figure 3). The correlation diagrams for the measurements of the other samples provide quite similar results, consequently the derived uncertainties for the individual layers become large.  The obtained total heights h SiO2 + h SiO compare quite well with the height values obtained for h SiO2 with the single layer model. Particularly for the thin SiO layers used in model 2, the statistical uncertainty of the layer thickness exceeds the expected value of the layer height, so we obviously cannot It is therefore necessary to use Bayesian statistics. With this, we calculated the posterior distributions: since no external a priori information is available, we use a uniform height distribution as the prior distribution. The posterior distribution was sampled numerically with the Markov chain Monte Carlo (MCMC) technique using the random walk Metropolis algorithm [16,17]. For more details, see [18] and [19].
For the analysis of our MME measurements by applying the two-layer model, we get posterior distributions for each best-fit SiO 2 and SiO layer, each of which gives the statistical probability density function for the height values. As expected, all of the obtained posterior distributions are asymmetric and not Gaussian-distributed. An example is illustrated in figure 4, for the measurement of the 6 nm nominal layer thickness sample. This asymmetry is particularly introduced by the natural lower limit at a height of 0 nm. As already discussed in [1], such distributions are relatively difficult to interpret and the whole of the posteriors should be regarded as the result, because the uncertainties are only barely derivable.
Further on we have calculated the posterior distributions of the complete oxide heights h SiO2 +SiO . The results are shown in figures 5-9. For each posterior we have calculated 30 000 sample points.

Comparison with calibration data
All posteriors obtained for the complete oxide heights h SiO2 +SiO are a particularly good approximation to a Gaussian distribution. Thus, we can interpret the width of the  obtained distributions as the statistical uncertainty contribution U stat of the MME measurements. The resulting mean values h MME = h SiO2+SiO and expanded uncertainties U stat (k = 2, for a confidence level of 95% [20]) are given in table 2. Please note that all the deterministic measurement errors are still missing, so that, for example, the influence of dispersion dataset uncertainties are yet to be included. However, the specification of uncertainty estimations for (n, k)-datasets are still not established, therefore the corresponding databases contain no uncertainties. Since the layer thicknesses and optical constants are correlated, the inaccuracy of the constants results in systematic errors in the determined thicknesses. Therefore, U stat is not yet a completely expanded uncertainty.
In table 2 we have listed the corresponding MME measurement results and statistical uncertainty estimations obtained for the complete oxide layer thicknesses and compared them with the certified calibration data h calib and U calib [21]. The uncertainties of the calibration data, as specified in the calibration certificates, comprise both the uncertainty of the XRR calibration and the uncertainty of the transfer using SE, as described above. The E n -value provides a very good comparative indicator of the agreement of the measurement results: [18] Although the obtained U stat are not yet providing the full measurement uncertainties of our MME measurements, we have used them to calculate the E n -values for this comparison. This already gives us a reasonable idea of the degree of conformity between the calibration values and our MME measurements, since the complete uncertainty U MME will always be larger than or equal to the statistical contributions: U MME ⩾ U stat . In any case, the observed E n -values that are well below one evidence an excellent and highly significant agreement. Only for the nominal 1000 nm layer the E n -value is a bit higher (0.63), but it still indicates a good agreement within the stated expanded uncertainties.

Analysis of depolarization behaviour
We further started to investigate the influence of depolarization on Mueller matrices especially in the short-wavelength range using as examples the (4,3)-M-elements of 6 nm, 160 nm and 1000 nm nominal layer thicknesses. The spectral behaviour of the polarization index P D [1,[23][24][25] for M 1,1 -normalized Mueller matrices with λ k being the eigenvalues of the corresponding Cloude covariance matrix H is shown below: as well as the probably more meaningful entropy S resulting from information theory [1,26] (P k : normalized eigenvalues): Also the associated χ 2 -distributions are illustrated. The selected (4,3)-M-elements (see figure 10) of the measured data for the stated layer thicknesses are now compared to the cases with (M CF the Cloude-filtered M [1]) and without depolarization information (M CF,nd ). To get M CF,nd , the M CF matrix is first transformed to its corresponding Cloude covariance matrix H, which has four eigenvectors and four eigenvalues in the interval (0,2). To get the nearest H without depolarization, the greatest eigenvalue is set to 2 and all the others are set to 0. With these new eigenvalues and the original eigenvectors we compute H now containing no depolarization (see [1]), which is then further transformed to M CF,nd [1]. One can see a significant deviation from M CF to M CF,nd especially in the short-wavelength range, which increases with increasing layer thickness. Therefore, the depolarization also increases with thickness and with reciprocal wavelength.
As a more quantitative investigation, both P D [23][24][25] and S [26], as plotted in figure 11, confirm that the depolarization increases with thickness and reciprocal wavelength. Furthermore, some characteristic depolarization peaks (already mentioned in e.g. [11,27]) are observable for the nominal 1000 nm standard in slightly decreasing energetic distances (from 0.55 eV to 0.40 eV), which would be equidistant without dispersion. For thinner layers (e.g. 6 nm and 160 nm) such peaks also occur, but shift to shorter wavelengths. This characteristic behaviour occurs due to the finite spectral resolution of the instrument, divergence of the ellipsometric light beam and non-uniformity of the oxide film thickness within the spot size [8,28,29], which mainly originates from the non-planarity of the air-SiO 2 and the SiO 2 -Si boundaries. Since the light beam therefore crosses slightly different oxide heights many times, this leads to thickness-dependent multiple interference (e.g. [30]). From this incoherent superposition the resulting beam yields an averaged measured Mueller matrix M (or equivalently a Cloude covariance matrix H) of local Mueller matrices over the area of the light spot of the ellipsometer.
A qualitative estimation confirmed that the characteristic peaks are introduced by interference effects.
Looking at the χ 2 -distributions of our merit function (see figure 12), which should optimally have a value of one over the whole spectral latitude, for the selected 6 nm, 160 nm and 1000 nm layer thicknesses one can see values very close to one for the 6 nm layer except in the UV range. Despite lots of depolarization peaks, in the case of the 1000 nm layer (see figure 11), please note the greatly weakened peaks occuring in the χ 2 -distribution.
For comparison with the standard M-optimization method, figure 13 illustrates the corresponding χ 2 -distributions for the same selection of layer thicknesses, which should optimally have a value of zero for the whole spectrum. For all layer thicknesses a stronger spectral dependence is observable, which increases from a relatively slight dependence at 6 nm with increasing values mainly in the UV range, again, to a strong dependence at 1000 nm. The depolarization peaks are still much more present in the case of the last thickness and therefore have a more significant perturbing influence on the data analysis than in our case.
The new merit function therefore overcomes the variations of the spectral contributions much better. This indicates once more, that with our merit function, the calculated uncertainties are more reliable than those computed by the current merit functions.
Other possible questions arising are: Can the depolarization due to the above-mentioned effects be neglected for modeling theoretical values in the case of 6 and 160 nm thicknesses, since it is relatively weak there? And: Does the depolarization in the case of the nominal 1000 nm thickness not have to be considered, to calculate the theoretical values for resilient results, owing to its characteristics visible in figure 11?
Regarding these questions, separate contributions to the discussion of different film thicknesses are presented below:  figure 11), again, there are usually only small deviations from one and zero, except in the UV range. The depolarization due to, inter alia, non-uniformity of the layer thickness, is therefore only rather weak. This could lead to the assumption that depolarization can be neglected in this case. In [1], however, we have already seen due to uncertainty simulations, that for the specification of a reliable measurement uncertainty estimation it is crucial to also take weak depolarization contributions into account, even for an oxide layer thickness of 2 nm. Hence, from our point of view, depolarization must additonally be considered in the simulated quantity of our merit function for such thinner layers as well. As soon as this is implemented it will confirm whether our current result for χ 2 and the excellent E n -values can be improved further.

To 1000 nm.
The behaviour of P D and S visible in figure 11 due to the finite spectral resolution of the instrument, etc, usually has a systematic deviation from one and zero. Of course, in this case, depolarization must certainly be considered for the calculation of theoretical values using suitable models realized by a formalism working with M or H [8,9,11]. This is part of the current research in our group but not covered in this manuscript. However, comparing the respective χ 2 -distributions in figures 12 and 13 again, we could already show evidence that our result is more reliable than that calculated by the established M-method. Furthermore, although the associated E n -value admittedly deviates about one order of magnitude from the others, it is nevertheless clearly smaller than one and therefore already indicates the presence of a significant overlap between our value and the calibrated value within the expanded uncertainties. Hence, the presented result for the 1000 nm oxide film is already resilient and shows that we are on the right track.

Summary and outlook
On our route to reliable uncertainty estimations and traceability in Mueller elipsometry, we concentrated first on stochastical uncertainty contributions induced by the measurement and analysis processes, namely measurement noise and depolarization.
We presented the first application of our improved analysis technique on MME measurements for the determination of the layer thicknesses of layer thickness standards with nominal values in the range of 6 nm to 1000 nm. We showed by a maximum likelihood analysis, that for all thickness standards an SiO 2 -SiO on Si model is a more adequate choice than a single SiO 2 layer model, as indicated by appling a modified BIC. We performed a Bayesian analysis via MCMC. The obtained height posterior distributions of the individual SiO 2 and SiO layers turned out to be non-Gaussian and strong anticorrelations between SiO 2 and SiO were observed. We also calculated the distribution of the total oxide layer thickness h SiO2 +SiO and obtained Gaussian distributions for all samples. We determined the mean values, standard deviations and also the statistical contribution to the expanded uncertainties U stat , respectively. These results of our MME measurements compare very well with the calibration values, so that all E n -values are significantly smaller than one. The inclusion of systematic uncertainty contributions might further increase the total uncertainty to a slight degree, thus decreasing the E n -values even further.
While the data analysis significantly favoured the two layer model (model 2) and it was notably possible to derive the combined oxide heights h SiO2 +SiO , due to the observed strong anticorrelation we could not derive the thicknesses of the individual SiO 2 and SiO-layers very meaningfully from the MME measurements. Further steps will be systematic investigations of the influence of dispersion dataset uncertainties and of further device-related uncertainty contributions (e.g. related to the retarder or spectrometer calibrations) on the MME measurement results to obtain a full uncertainty budget, as an important step towards self-traceable MME metrology. In addition, a more precisely mathematical description of coherence and incoherence during the measuring process, addressing the previously-mentioned finite spectral resolution, light beam divergence and non-uniformity of the layer thickness under test, will be of interest as an advancement of our current merit function. Since an averaged Mueller matrix M (or equivalently an averaged Cloudecovariance matrix H) is measured, an averaged M (or H) must also be calculated for the simulatable quantity in the modeling process. The local M can be calculated using the assumption that the thickness is constant within sufficiently small areas of spot size. However, they depend on the position of the point within the spot size for which the local M is calculated. A general formalism for this is presented, for example, in [8,9,11], but this is still not addressed in our merit function and hence it is currently a work in progress in our group.