Reconstructing Detailed Line Profiles of Lamellar Gratings from GISAXS Patterns with a Maxwell Solver

TL;DR

This paper demonstrates that GISAXS combined with a Maxwell solver can accurately reconstruct detailed nanostructure line profiles, including critical dimensions, with sub-nanometer precision, enabling advanced dimensional metrology.

Contribution

It introduces a novel approach using Maxwell's equations and GISAXS to reconstruct detailed nanostructure profiles with high accuracy, validated by statistical methods.

Findings

Reconstructed line widths down to 55 nm.

Achieved sub-nanometer uncertainty in critical parameters.

Validated the method with MCMC sampling.

Abstract

Laterally periodic nanostructures were investigated with grazing incidence small angle X-ray scattering (GISAXS) by using the diffraction patterns to reconstruct the surface shape. To model visible light scattering, rigorous calculations of the near and far field by numerically solving Maxwell's equations with a finite-element method are well established. The application of this technique to X-rays is still challenging, due to the discrepancy between incident wavelength and finite-element size. This drawback vanishes for GISAXS due to the small angles of incidence, the conical scattering geometry and the periodicity of the surface structures, which allows a rigorous computation of the diffraction efficiencies with sufficient numerical precision. To develop dimensional metrology tools based on GISAXS, lamellar gratings with line widths down to 55 nm were produced by state-of-the-art e-beam lithography and then etched into silicon. The high surface sensitivity of GISAXS in conjunction with a Maxwell solver allows a detailed reconstruction of the grating line shape also for thick, non-homogeneous substrates. The reconstructed geometrical line shape models are statistically validated by applying a Markov chain Monte Carlo (MCMC) sampling technique which reveals that GISAXS is able to reconstruct critical parameters like the widths of the lines with sub-nm uncertainty.

Tables

Table 1. Table 1: The optimization results of the geometrical model parameters obtained from gratings G1 and G2 with the MCMC sampling, with corresponding confidence intervals ( ± 1 ​ σ plus-or-minus 1 𝜎 \pm 1\sigma ).

Parameter	Grating G1	Limits	Grating G2	Limits
pitch / nm	100	fixed	150	fixed
line height / nm	102.71 $\pm 0.12$	[95 105]	119.50 $\pm 0.11$	[115 125]
line width / nm	54.04 $\pm 0.56$	[50 60]	67.30 $\pm 0.31$	[60 70]
top r / nm	16.27 $\pm 0.39$	[0 20]	9.16 $\pm 0.39$	[0 20]
bot. r / nm	15.79 $\pm 0.24$	[0 30]	13.02 $\pm 0.58$	[0 30]
sidewall angle / ∘	90.91 $\pm 0.38$	[80 95]	84.73 $\pm 0.33$	[80 90]
${\sigma }_r$ (rms) / nm	1.35 $\pm 0.07$	[0 5]	1.87 $\pm 0.14$	[0 5]

Equations

$$\left(\begin{array}[]{c}q_{x}\\ q_{y}\\ q_{z}\end{array}\right)=\left(\begin{array}[]{c}k_{0}(\cos(\theta_{f})\cos(\alpha_{f})-\cos(\alpha_{i}))\\ k_{0}(\sin(\theta_{f})\cos(\alpha_{f}))\\ k_{0}(\sin(\alpha_{f})+\sin(\alpha_{i}))\end{array}\right)\\$$

$$I=I^{\text{model}}\cdot \exp \left(-{\sigma }_r^2{q}_y^2\right)\text{.}$$

$$\text{MANE}=\frac{1}{M}\bigg{(}\sum\limits_{N}\frac{|I^{\textrm{quasi}}_{N}-I^{\textrm{model}}_{N}|}{I^{\textrm{quasi}}_{N}}\bigg{)}\text{,}$$

$${\chi }^2=\sum _E{\tilde{\chi }}^2(\text{E})\text{,}$$

$$\tilde{\chi}^{2}(E)=\sum\limits_{m}\frac{\big{(}I_{m}^{\text{model}}(E)-I_{m}^{\text{meas}}(E)\big{)}^{2}}{{\sigma}^{2}(E)}\text{.}$$

$${\sigma }^2(E)={\sigma }_{\text{num}}^2(E)+{\sigma }_{\text{N}}^2+{\sigma }_{\text{hom}}^2+{\sigma }_{\text{div}}^2\text{.}$$

$$P(\vec{x})\propto\exp\big{(}-{\chi}^{2}/2\big{)}\text{,}$$

Full text

I Introduction

Measurements on the length scale of several nm are challenged by the atomic granularity of matter and by structures which cannot easily be described with simple models. The continuously shrinking patterns in the semiconductor industry are at the forefront of technological development regarding the requirements for size reproducibility and regularity.

Scanning probe microscopy techniques (e.g. atomic force microscopy (AFM), scanning tunnelling microscopy (STM), scanning near-field optical microscopy (SNOM), scanning electron microscopy (SEM)) are powerful tools for the investigation of nanostructured surfaces and occupy the key positions for metrology tools in industry. However, in particular X-ray scattering is also an established technique in nanoscience. Grazing incidence small angle X-ray scattering (GISAXS) Guinier and Fournet (1955); Levine et al. (1989); Renaud et al. (2009) also offers comparably fast measurements, in addition to being destruction-free. With incidence angles close to the critical angle $\alpha_{c}$ of total external reflection, GISAXS is a technique with high surface sensitivity and is perfectly suited for in-situ applications. Due to the large photon beam footprint as compared to the nm pattern size, it directly yields statistical information on fluctuations such as structure roughness for a large structured area. However, the long beam footprint prevents probing specific small sample volumes in the illuminated area. A challenge for GISAXS is the characterization of structured surfaces with complex periodic and non-periodic sample layouts (e.g. photomasks). Changes in the design of the sample layout Pflüger et al. (2017) allow to overcome this problem, but this results also in an extreme loss in scattering intensity which eliminates one of the major advantages of GISAXS for in-situ applications. In contrast to GISAXS, transmission SAXS Wang et al. (2007); Jones et al. (2006); Sunday et al. (2015) (CDSAXS) is able to probe much smaller sample volumes, only limited by the overall beam dimension. The drawback of the transmission geometry results from the fact that a significant portion of the incident beam may be absorbed in the substrate. As the incidence angle moves away from normal incidence, the beam path length increases and for typical silicon wafers X-ray photon energies above 13 keV are needed for sufficient transmissionSunday et al. (2016). A combination of grazing incidence reflection and transmission measurements Lu et al. (2013) (GTSAXS) at the sample edges is often not practical because typical samples are not structured to the very edge. Therefore, GISAXS is ideally suited for the characterization of nano structured surfaces on non-homogeneous substrates.

Several groups have already performed GISAXS measurements on gratings for co-planar Tolan et al. (1995); Metzger et al. (1997); Jergel et al. (1999); Mikulík and Baumbach (1999) and conical scattering geometry Mikulík et al. (2001); Yan and Gibaud (2007). The diffraction patterns are well understood. The intensity modulation of the scattered orders is related to the line shape (form factor) of the grating. To reconstruct the line shape parameters, including the line width, the sidewall angle or the line height, the so-called inverse problem has to be solved. The corresponding calculations for GISAXS were mostly done using the distorted wave Born approximation (DWBA) including an analytic expression of the form factor Babonneau (2010); Lazzari (2002); Renaud et al. (2009); Rueda et al. (2012); Rauscher et al. (1995); Jiang et al. (2011); Hofmann et al. (2009); Suh et al. (2016); Meier et al. (2012). Arbitrarily shaped structures must be discretized and require numerical solvers Chourou et al. (2013). In the optical domain of the electromagnetic spectrum, modelling the light scattering by numerically solving the time-harmonic Maxwell’s equations with a higher-order finite-element method Pomplun et al. (2007); Kato et al. (2012) is well established. If the periodic structures are invariant in one dimension (i.e. along the grating lines), the computational domain can be reduced to two dimensions which significantly decreases the computational effort and also allows the calculation of rather large domains as compared to the incident wavelength. This enables, together with the shallow incident angle in GISAXS, the application of this approach to short wavelengths as in X-ray scattering. Solving an inverse problem means minimizing the difference between the measured scatter intensities and the simulated intensities by adapting the geometrical parameters. A statistical validation of the optimized models is possible with the Markov chain Monte Carlo (MCMC) approach Foreman-Mackey et al. (2013), which gives the opportunity to obtain parameter sensitivities and delivers confidence intervals.

Besides the regular diffraction pattern from a line grating, the diffuse scattering background in the GISAXS pattern reveals additional information about the surface. The appearance of long-range ordered superstructures in the scattering pattern yields information about the e-beam fabrication process. Roughness and imperfections of the surface structures produce a complex diffuse scattering background, which can be described by kinematic and dynamic scattering effects. These effects are correlated with the grating line shape and also allow the determination of geometrical parameters Meier et al. (2012); Soltwisch et al. (2016).

II Experimental Details

The basic geometry of GISAXS is shown in Fig. 1. A monochromatic X-ray beam idealized as a plane wave with the wave vector $\vec{k}_{i}$ impinges on the sample surface at a grazing incidence angle $\alpha_{i}$. The elastically scattered wave with the wave vector $\vec{k}_{f}$ propagates along the exit angle $\alpha_{f}$ and the azimuthal angle $\theta_{f}$. Here we use the common notation for the scattering vector $\vec{q}=\vec{k}_{f}-\vec{k}_{i}$

[TABLE]

using angular coordinates (see Fig. 1) and $k_{0}=2\pi/\lambda$ with $\lambda$ as the wavelength. Scattering from a periodically structured surface, e.g. a grating, leads to a characteristic diffraction pattern Tolan et al. (1995); Hofmann et al. (2009); Wernecke et al. (2012); Rueda et al. (2012); Soccio et al. (2014). The GISAXS pattern from a line grating observed with a 2D detector can be easily understood by a reciprocal space construction Mikulík et al. (2001); Yan and Gibaud (2007). The reciprocal space of a line grating, with pitch $p$, consists of truncation rods which are perpendicular to the surface and aligned with a spacing $\Delta q_{y}$ in the direction perpendicular to the scattering plane. The scattering pattern on the detector arises as a semicircle of discrete diffraction orders with $\Delta q_{y}=2\pi/p$ which follows from the intersection of the Ewald sphere with a radius of $k_{0}$ and the reciprocal space representation of the lattice at $q_{x}=0$ (see Fig. 1).

The GISAXS experiments were conducted at the four crystal monochromator beamline, operated by the Physikalisch-Technische Bundesanstalt (PTB), at the electron storage ring BESSY II in Berlin Krumrey and Ulm (2001). The beamline offers a photon energy from $E=$ 1.75 keV to $E=$ 10 keV. By using a beam-defining 0.52 mm pinhole about 150 cm before the sample position and a scatter guard 1 mm pinhole about 10 cm before the sample, we reach a beam spot size of about 0.5 x 0.5 mm at the sample with minimal parasitic scattering. A 6-axis goniometer installed in a UHV chamber allows the alignment of the scattering angle $\alpha_{i}$ with an accuracy of $\pm 0.001^{\circ}$. The grating lines were aligned in parallel orientation ($\varphi=0^{\circ}$) with respect to the incident beam with an accuracy of $\pm 0.002^{\circ}$. An in-vacuum PILATUS 1M hybrid pixel detector is installed on a movable sledge which allows the sample-detector distance to be varied from 1.7 m up to 4.5 m Wernecke et al. (2014). The detector features 20 bit counters for every pixel and consists of 10 separate modules. To extend the accessible photon energy range to the soft X-ray region (EUV), small angle X-ray scattering was also performed at the soft X-ray beamline Scholze et al. (2001). An Andor CCD camera was mounted on a UHV chamber with a fixed sample detector distance of $0.75$ m, which allows the performance of extreme ultraviolet small angle scattering (EUV-SAS) at steeper incident angles $\alpha_{i}$ of $7^{\circ}$.

Two lamellar gratings were manufactured by electron beam lithography. Grating G1 with a nominal pitch of 100 nm and a line width of 55 nm and grating G2 with a nominal pitch of 150 nm and a line width of 65 nm (see Table 1 and Fig. 2). The grating areas measure 1 mm by 15 mm (grating G1) and 0.51 mm by 4 mm (grating G2) with the lines oriented parallel to the long edge. To manufacture the gratings, a silicon substrate was spin coated with the positive resist ZEP520A (organic polymer). Pattern generation was done using a Vistec EBPG5000+ e-beam writer, operated with an electron acceleration voltage of 100 kV Senn et al. (2011). After resist development, the grating was etched into the silicon substrate via reactive ion etching, using the etching gases SF6 and C4F8. Finally, the remaining resist was removed with oxygen plasma treatment. Grating G1 was fabricated in 2013 and grating G2 in 2016 using a slightly different etching process. The resulting line shape is very sensitive to the conditions in the reactive ion etching chamber, which changed slightly due to other processes running in the same chamber between the manufacturing of the two samples, leading to the differences in geometry seen in Table 1 and Fig. 6.

III Superstructures

Besides the intense diffraction spots from the grating structure (see Fig. 3), the GISAXS pattern also shows a dominant superstructure in the ($q_{y},q_{z}$) pattern as additional diffraction semicircles around the main grating diffraction semicircle at $q_{x}=0$ (see Fig. 3 a)). The non-periodic shifting of the superstructures in the $q_{z}$ direction is clear evidence that the periodic modulation must be lateral along $q_{x}$. This behaviour is also illustrated in Fig. 1. A periodic modulation with the pitch $P$ of the grating in the lateral x direction leads to additional grating truncation planes with $\Delta q_{x}=2\pi/P$. This periodicity appears in the detector plane as additional GISAXS diffraction semicircles with the same $\Delta q_{y}$ periodicity for the diffraction orders but including the $\Delta q_{x}$ spacing of the lateral modulation. For the evaluation of the lateral periodicity, the corresponding ($q_{y},q_{x}$) projection must be calculated. With grazing incidence angles around 1*∘* and a photon energy of 6.5 keV, the observable scattering vector component $q_{x}$ is in the order of $\mu\text{m}^{-1}$, corresponding to $\mu$m sized structures on the sample surface (cf. Fig. 3 b)). This allows the extraction of the lateral pitch size (along the grating lines) as $\Delta x=$ 4.5 $\mu$m for the superstructure. This superstructure is directly related to the stitching field size of the e-beam writer (4.53$\times$ 4.53) $\mu\text{m}^{2}$. The lateral pitch size of the superstructure perpendicular to the scattering plane is not visible in the GISAXS scattering pattern, due to the different range for the $q_{y}$ scattering vector, but can be obtained with EUV-SAS at $\alpha_{i}$ = 7*∘* and photon energies around $E=80$ eV (see Fig.  3 c)). Thus, the measurement demonstrates that small angle scattering techniques are also able to extract long-range ordering of nm structured surfaces due to the elongated beam footprint. The finite-element approach in the X-ray spectral range, however, is limited to a two dimensional domain assuming a perfect infinite grating along the scattering direction. Such long-range perturbations can only be described in the model as a form of roughness.

IV Computational diffraction intensities

The reconstruction of the geometrical layout of nm sized structures by measuring the intensity distribution of the scattered photons means solving the inverse problem. In an iterative process, a theoretical model is used to compute the scattered intensity from a guess about the structure and compares the measured intensities with the simulation. In the last decade CDSAXS was developed as a powerful tool for the characterization of periodic nano structured surfaces. One major advantage of all transmission SAXS experiments, lies in the fact that a theoretical modelling with the first Born approximation (BA) is suitable for a surface shape reconstruction, as demonstrated by several groupsHu et al. (2004); Wang et al. (2009); Sunday et al. (2015, 2016); Lu et al. (2013). The incoming wave scatters only once at the target before forming the scattered wave. With an analytic description of the line shape (form factor) this modelling approach is relatively fast.

With thick and non-homogeneous substrates the GISAXS measurement method is an interesting option for surface characterization. With decreasing incidence angle $\alpha_{i}$ the scattered amplitude increases while the transmission decreases. Incidence angles close to the external reflection allow to probe targets with a high surface sensitivity due to the comparably small penetration depths. This requires an extension of the BA, due to the fact that the BA neglects any multiple scattering effects. For GISAXS the distorted wave Born approximation is well developed in several Software packages (IsGisaxsLazzari (2002), BornAgainDurniak et al. (2014), HipGISAXSChourou et al. (2013), FitGisaxsBabonneau (2010)) which are able to deal with second order multiple scattering effects. They extend the BA with additional scattering contributions from different multiple scattering events. The DWBA is able to explain a lot of GISAXS measurements, but higher order scattering effects (e.g. higher order Yoneda scatteringSoltwisch et al. (2016)) are not covered in the typical implementation of the DWBA. The biggest advantages of both perturbation theories originate from the simplicity of the analytic form factor description and the resulting computational speed. This strong advantage reduces for arbitrary form factors, which must be discretized numericallyChourou et al. (2013). Beside those perturbation theories, solving the Maxwell’s equation is also possible with a finite-element approach. Arbitrarily shaped objects can be implemented easily as vector functions and parametrized. In contrast to the approximation methods, a Maxwell solver allows the computation of real far field scattering intensities, including the solution of the near field in the computational domain. Higher order scattering effects can be easily studied with the computed electromagnetic field distributions inside the scatter objects. Additionally, the computation of local field distributions gives rise to the simulation of depth depended absorption measurements, e.g. grazing incidence X-ray fluorescence (GIXRF). Coupled with the strong increase in the available computational power in the last decade, solving the Maxwell’s equations based on finite-elements is an interesting option in the X-ray spectral range. It is still limited to specific experimental settings (small angles, infinite gratings), i.e. an expansion to CDSAXS is not applicable at the moment, but possible future applications are only limited by the computational power.

In the following we will give only a rather compact introduction on the topic of finite-element discretization. A detailed summary can be found elsewhereMonk (2003); Pomplun et al. (2007).

X-rays are treated as electromagnetic plane waves of wavelength $\lambda=hc_{0}/E$ (with the Planck constant $h$, speed of light $c_{0}$ and photon energy $E$) which scatter on nano-structures. In this case, the set of Maxwell’s equations can be rewritten to a single, second order curl-curl equation for the electric fieldPomplun et al. (2007). The general idea of the finite element discretization is that the computational domain is subdivided into small patches (e.g. triangles). On these patches a vectorial ansatz function is defined usually with polynomials with a fixed order $p$. The approximate electric field solution is the superposition of these local ansatz functions. The numerical accuracy of the approximate electric field is a functional of the size of the finite-elements and the degree $p$ of the polynomials.

To compute diffraction efficiencies, we assume periodic structures which are invariant in one dimension (along the grating lines). Well-converged solutions are obtained with the package JCMsuite Pomplun et al. (2007) using a higher-order finite-element method.

We focus on a reconstruction of the form factor (the line shape) and not the structure factor (the grating pitch) which is well known from lithography Buhr et al. (2007). Therefore, we fixed the pitch of the structures to the nominal values in the reconstruction process. In Fig. 4, a computational domain and the corresponding near field calculation are shown for a typical GISAXS measurement geometry ($\alpha_{i}=0.86^{\circ}$, $E=5.5$ keV). The sketched line profile of the grating emphasizes the important structural parameters which were optimized in the reconstruction: the line width, the line height, the sidewall angle, the top corner rounding and the groove rounding. The line width of the structure is defined at the half-height. The parameter of the top corner rounding describes a circle with the respective radius. The bottom groove rounding is constructed with an elliptical shape. The major radius depends on the line width and the sidewall angle, the minor radius is parameterized for the reconstruction. Due to the reactive ion etching process, one could expect a strong groove rounding and only a minor rounding of the top corners because the top surface is still covered by the resist during etching (also visible in the SEM cross-section images, see Fig. 6). This model allows an adequate description of the expected line shape, with a low number of parameters.

Roughness or imperfections of the lamellar gratings cannot be modelled directly with the finite-element approach due to the small discretization length required in the X-ray region. The computational domain size would go beyond any available computer memory. To account for the line edge or line width roughness damping effects on the diffraction intensities $I$ of higher diffraction orders, an analytic correction must be applied. The line edge and line width roughness (LER/LWR) of the gratings were taken into account, in the optimization process, with the well-known analytic approach Mikulík and Baumbach (1999); Kato and Scholze (2010) based on Debye-Waller damping

[TABLE]

This allows the correction of the undisturbed computational diffraction intensities $I^{\textrm{model}}$ from the Maxwell solver in a post-process. The damping factor $\sigma_{r}$ of the rms line roughness was also included in the optimization process.

V Uncertainty evaluation

The grazing incidence conical diffraction and the invariance of the grating in the direction of the scattering plane result in a standing wave field with a much larger period than the wavelength (see Fig. 4). This allows a significant increase in the size of the required discretization length $d$. It breaks the conventional rule of half the incident wavelength for the discretization to ensure numerical accuracy. This allows the use of a Maxwell solver based on the finite-element method to efficiently treat GISAXS applications.

To ensure numerical accuracy with the finite angle of incidence, the relative numerical error of the simulated diffraction intensities (far field) was calculated for different numerical precision settings. The two main numerical degrees of freedom are the spatial discretization length $d$ and the polynomial degree $p$ which defines the ansatz function used to approximate the fields. The numerical errors are defined as the difference of the actual diffraction intensities $I^{\textrm{model}}$ to the quasi-exact results $I^{\textrm{quasi}}$. The quasi-exact calculation is defined as the converged computation with the highest achievable numerical precision settings, typically limited by the amount of available computational memory. The numerical errors for typical GISAXS settings are shown in Figs. 5 a) and b) for a photon energy of $E=8$ keV and varying grazing incidence angles $\alpha_{i}$. Similar results can be obtained for the azimuthal angle $\varphi$. This figure reveals the coupling of $d$ and $p$ to the incident angle $\alpha_{i}$ and allows the trade-off between the numerical precision and the computational effort to be estimated. For example, a decrease in the incident photon energy would shift the numerical precision in the figure towards larger incident angles. In this rough estimation, the size of the scatter objects represents a natural border of the discretization length $d$, but modern implementations of the finite-element method allow to make use of adaptive meshing algorithmsBurger et al. (2015) which are able to tune the local discretization length, e.g. in critical regions like the vacuum-silicon interface. However, the obvious gap between the incident wavelength ($\lambda\ll 1$ nm) and the sufficient discretization length $d$ gives the opportunity for fast (below 1 s) and accurate simulations with a Maxwell solver based on the finite-element method. It should be noted that in Figs. 5 a) and b), only the mean absolute numerical error (MANE) of all diffraction orders is visualized

[TABLE]

with $M$ being the number of diffraction orders and $N$ the specific diffraction order. This hides the fact that the numerical error for the different diffraction orders is coupled with the geometrical layout and the actual experimental settings. Diffraction orders with an exit angle $\alpha_{f}$, which is close to the critical angle of the substrate or the grating effective layer (Yoneda band) Soltwisch et al. (2016), are most affected. The number of usable diffraction orders is critical for the reconstruction process, as it directly limits the possible complexity of the reconstructed model. To achieve an upper estimation of the numerical errors, we simulated 1000 randomly distributed geometrical grating layouts within the parameter boundaries for the extraction of the quasi-exact diffraction intensities. The identical grating layouts were simulated with the numerical discretization which will later be used in the specific grating reconstruction process (see Table 1). We then compared every diffraction intensity with the corresponding quasi-exact computation. The results of this numerical error estimation for the reconstruction of grating G2 are exemplified visually in Figs. 5 c) and d) for discretization lengths of $5$ nm and $6$ nm. The dashed red lines mark the $\sigma$ interval. These $\sigma$ values in the range of several percent reveal that numerical precision has a significant impact in the uncertainty evaluation. To achieve a better reconstruction result, the highest diffraction orders should be neglected in the fitting process due to the numerical instabilities at the critical angle (c.f. Fig. 3 a) Yoneda band). We also excluded the diffraction orders which are not used in the subsequent reconstruction process (e.g. zero order) of grating G1 and G2.

Besides these numerical errors, small angle X-ray scattering is challenging for every experimental setup, especially if the accuracy of the alignment angles $\alpha_{i}$ and $\varphi$ is important. Under grazing incidence, small variations in the angles of incidence could change the measured diffraction efficiencies rapidly. This angular sensitivity is demonstrated in Fig. 5 e). The red dots are the simulated diffraction intensities (quasi-exact) of a typical lamellar silicon grating (grating G1) at $\alpha_{i}=0.7^{\circ}$ and $E=8$ keV. The changes in the diffraction intensities for angular variations in the order of $10^{-3}$ degrees are visualized with dashed lines. These angular uncertainties are close to the achievable motor step resolution in our experimental chamber. This is in principle not a big issue for the reconstruction process, because a subsequent calibration of the angles is still possible. But this calibration is also limited by the uncertainty from the estimation of the pixel size for the PILATUS detector Wernecke et al. (2014) and leads to a very similar uncertainty of $0.002^{\circ}$. To ensure an accurate reconstruction of the line shapes, both incidence angles must be included in the optimization process as independent parameters. This angular sensitivity in grazing incidence also influences the computational effort which is needed for an accurate simulation of measured diffraction intensities due to the divergence of the incident beam. The horizontal divergence of $0.01^{\circ}$ leads to an elongation of the diffraction peaks along the $q_{z}$ axis. The impact is directly visible in Fig. 3 a) for diffraction orders close to the horizon. To account for this angular distribution in the theoretical evaluation, we calculated 5 azimuthal incidence angles $\varphi$ for every simulation. The extracted diffraction intensities were weighted with a Gaussian distribution and the standard deviation of the incident beam divergence. Convergence studies with increasing numbers of incidence angles $\varphi$ (up to 50) reveal that the relative uncertainty $\sigma_{\text{div}}$ of the sampling with only 5 angles is in the region of 8%. A further reduction of this uncertainty is possible but would result in a significantly increased computational effort.

By using a photon counting detector, which follows Poisson statistics, for the measurements of diffraction intensities, the statistical measurement uncertainty $\sigma_{\text{N}}$ is typically an order of magnitude below the numerical uncertainty of the Maxwell solver. However, the homogeneity of the PILATUS detector at incident photon energies around 5 keV Wernecke et al. (2014) is not negligible, which leads to an additional uncertainty $\sigma_{\text{hom}}$ of 2.5% in the measured diffraction intensities. The distribution of the incident photon energies were taken into account, with the implementation of a Gaussian prior, with the bandwidth of $10^{-4}$ for each photon energy. For hard X-rays, the influence of uncertainties from the optical constants has almost vanished, especially for silicon far away (several keV) from any absorption edge. In the soft X-ray or EUV region, the situation could change dramatically and must be evaluated separately.

VI Reconstruction of the line shape

The general problem in a reconstruction process is the question of how much information is necessary to obtain a univocal solution. The constraints for the numbers of diffraction orders for a specific incident angle and photon energy are dictated by the grating pitch. This strictly limits the accessible information and can only be compensated for by mapping different regions of the reciprocal space (e.g. tuning the incidence angles or photon energy). To avoid the experimental uncertainties during angular scans ($\alpha_{i},\varphi$), we utilize the high stability of the four crystal monochromator and tune only the incident photon energy. The intensity of the specular reflection is often disturbed by the mismatch between the beam spot size and the grating target size. Parts of the beam are reflected from the surrounding substrate and not from the structured surface. This distorts the zero order diffraction intensity and prevents the extraction of absolute intensities. Fitting the relative diffraction intensities of the non-zero orders is thus more accurate for sample sizes below the size of the elongated beam footprint.

To compare the measured diffraction intensities with theoretical values, the Maxwell solver computes the near field solution for a specific parameter set depending on the model and extracts the theoretical photon flux for every diffraction order by post-processing using a Fourier transformation. The minimization functional of the optimization problem $\chi^{2}$ is defined as the total sum of the least-squares functionals for every photon energy $E$,

[TABLE]

where each of the functionals is defined as

[TABLE]

For the standard deviation $\sigma(E)$, the individual experimental and numerical uncertainty is propagated according to the Gaussian error propagation law,

[TABLE]

The numerical uncertainties $\sigma_{\text{num}}(E)$ were precalculated for the different photon energies which were used in the measurements and for both gratings, as described in the previous section. The discretization length $d$ and the polynomial degree $p$ were chosen such that the numerical uncertainties were $\sigma_{\text{num}}(E)\approx 2.6\%$ (grating G2) and $\sigma_{\text{num}}(E)\approx 3.6\%$ (grating G1).

A heuristic optimization method is ideally suited to explore a large parameter space of the inverse problem, but requires massive parallelization to reduce the computational effort to a reasonable time. With a well-optimized Maxwell solver and a commercially available workstation, we were able to solve 106 structures in less than one day. We used the particle swarm optimization Kennedy (2010) method, which ideally delivers the global minimum of the total $\chi^{2}$ functional with an acceptable computational effort. However, no information about parameter sensitivity or the quantification of confidence intervals is delivered with a particle swarm optimization. To solve this issue, we applied an affine invariant Markov chain Monte Carlo sampling technique Foreman-Mackey et al. (2013). The likelihood is given by

[TABLE]

where $\vec{x}$ is the set of parameters of the model. As a starting point, a random set of parameters was generated with respect to predefined boundaries around the global minimum of the particle swarm optimization. The confidence intervals and mean values were calculated by evaluating the probability distribution for each parameter, after the MCMC procedure converged. The confidence intervals ($\pm 1\sigma$) given here represent percentiles of the number of samples found in the interval defined by the upper and lower bounds.

The results of the MCMC optimization are summarized in Fig. 6 and Table 1. In Fig. 6, the extracted diffraction intensities for both gratings are shown as a function of the scattering vector component $q_{y}$ and photon energy, and compared to the simulated diffraction intensities of the MCMC optimization. Diffraction orders which were influenced by the detector gaps (dark regions in the GISAXS images, see Fig. 3) were manually taken out of the reconstruction (missing points in Fig. 6). The good agreement between the measurement and simulation is evidence that the chosen model is able to describe the diffraction patterns from the grating structures. Also, the reconstructed geometry fits well to the expected nominal values. This is illustrated in Fig. 6 by a comparison of the reconstructed line shape (c.f. red line) with cross-section SEM images which we obtained from witness samples. A direct comparison depends on the homogeneity of the grating sample, because the SEM images show only a rather small part of the grating in contrast to the GISAXS measurements, which capture the whole structured surface. However, the reconstructed grooves show a significant rounding which is expected from the masked ion beam etching process. This leads to well-defined corners on the top and more rounded structures inside the grooves. This behaviour corresponds well with the reconstructed line shapes of both gratings.

A comparison of the reconstructed parameter values of both gratings and the estimated confidence intervals can be found in Table 1.

The confidence intervals of the height are rather small with $\pm 120$ pm as compared to expected process variations of the etching. One could argue that the estimated confidence intervals correspond to the well-defined mean value of the grating line height over the large illuminated beam footprint area. The measured diffraction intensities can be understood as superposition of the diffraction patterns from several different gratings with slightly different parameters due to stochastic roughness or other systematic deviations in the etching process. This effect was accounted for by the inclusion of the Debye-Waller damping in our model. It should, however, be mentioned that this Debye-Waller damping approach may not hold for any kind of process variation as it was derived for line edge roughness in a binary grating model. In particular, the line height roughness does not disturb the two-dimensional periodicity of the surface, as line edge roughness does, and might have a different impact. The estimated confidence intervals are therefore only valid within the framework of the model presented here. For a complete uncertainty evaluation in GISAXS, the impact of these modelling assumptions must be further investigated.

VII Summary

We investigated lamellar gratings etched in silicon with pitch sizes of 100 nm and 150 nm using GISAXS in conical orientation. The very distinct diffraction patterns from the rather similar grating line shapes highlight the high sensitivity of GISAXS to structure details. A dominant superstructure in the diffraction patterns is correlated to artefacts of the e-beam writing process and may in future be exploited to further tune the lithography process. The conical diffraction geometry combined with grazing incidence angles allows the use of the finite-element approach with a Maxwell solver in the X-ray spectral range. The Maxwell solver in conjunction with the finite-element method is a versatile tool to solve the inverse problem for highly periodic structures because it allows arbitrary geometric models of the line shape to be parameterized. A detailed reconstruction of the line shape is possible due to a significant reduction of the computational effort by adapting the discretization length to the grazing angle of incidence. It reveals geometric details, such as strong rounding inside the grooves, which are often not accessible with other direct non-destructive measurement methods like AFM. A first parameter sensitivity evaluation based on the Markov chain Monte Carlo method is presented.

VIII Acknowledgment

We thank the European Metrology Research Programme (EMRP) for financial support within the Joint Research Project IND17 ”Scatterometry”. The EMRP is jointly funded by the EMRP participating countries within EURAMET and the European Union.