research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775

A semi-analytical approach for the characterization of ordered 3D nanostructures using grazing-incidence X-ray fluorescence

aMESA+ Institute for Nanotechnology, University of Twente, The Netherlands, bPhysikalisch-Technische Bundesanstalt, Berlin, Germany, and cNRC Kurchatov Institute, Moscow, Russia
*Correspondence e-mail: k.nikolaev@protonmail.com

Edited by S. M. Heald, Argonne National Laboratory, USA (Received 2 October 2019; accepted 3 December 2019; online 11 February 2020)

Following the recent demonstration of grazing-incidence X-ray fluorescence (GIXRF)-based characterization of the 3D atomic distribution of different elements and dimensional parameters of periodic nanoscale structures, this work presents a new computational scheme for the simulation of the angular-dependent fluorescence intensities from such periodic 2D and 3D nanoscale structures. The computational scheme is based on the dynamical diffraction theory in many-beam approximation, which allows a semi-analytical solution to the Sherman equation to be derived in a linear-algebraic form. The computational scheme has been used to analyze recently published GIXRF data measured on 2D Si3N4 lamellar gratings, as well as on periodically structured 3D Cr nanopillars. Both the dimensional and structural parameters of these nanostructures have been reconstructed by fitting numerical simulations to the experimental GIXRF data. Obtained results show good agreement with nominal parameters used in the manufacturing of the structures, as well as with reconstructed parameters based on the previously published finite-element-method simulations, in the case of the Si3N4 grating.

1. Introduction

Achievements in the field of science and technology related to the manufacturing of nanoscale devices are usually associated with the systematic decrease of the characteristic sizes of the structures within such devices. Such a decrease in characteristic sizes can lead to a strong performance dependency on minor variations in the device structure, its geometry and the elemental composition of different elements inside the structure. Prominent examples of such nanoscale device structures can be found in the microelectronics industry (Markov, 2014[Markov, I. L. (2014). Nature, 512, 147-154.]; Buitrago et al., 2016[Buitrago, E., Fallica, R., Fan, D., Kulmala, T. S., Vockenhuber, M. & Ekinci, Y. (2016). Microelectron. Eng. 155, 44-49.]). Understanding and improving the performance of such devices therefore requires the use of nanometrology techniques which, at best, are capable of reconstructing the geometry of the structure and the three-dimensional (3D) atomic concentration distributions of different elements. Such element-selective analysis can be performed using grazing-incidence X-ray fluorescence (GIXRF) (Soltwisch et al., 2018[Soltwisch, V., Hönicke, P., Kayser, Y., Eilbracht, J., Probst, J., Scholze, F. & Beckhoff, B. (2018). Nanoscale, 10, 6177-6185.]; Andrle et al., 2019[Andrle, A., Hönicke, P., Schneider, P., Kayser, Y., Hammerschmidt, M., Burger, S., Scholze, F., Beckhoff, B. & Soltwisch, V. (2019). Proc. SPIE, 11057, 110570M.]). GIXRF is based on the X-ray standing wave (XSW) which is excited due to the interference between incident and reflected radiation. Its position- and angle-dependent amplitude can substantially modulate the GIXRF intensities of an element depending on its location within the nanostructure. By varying the angle of incidence and/or incident photon energy, the location of the XSW field nodes and anti-nodes can be varied inside the nanostructure. Consequently, the emission of fluorescence radiation depends on the incident angle and the incident photon energy, as well as on the spatial distribution of the fluorescent atoms.

Measurement procedures and data analysis for one-dimensional (1D) depth distributions of fluorescent atoms have been well developed (Bedzyk et al., 1989[Bedzyk, M. J., Bommarito, G. M. & Schildkraut, J. S. (1989). Phys. Rev. Lett. 62, 1376-1379.]; Jörg & Kazimirov, 2013[Jörg, Z. & Kazimirov, A. (2013). The X-ray Standing Wave Technique: Principles and Applications, Vol. 7 of Series on Synchrotron Radiation Techniques and Applications. Singapore: World Scientific.]) and implemented for the study of epitaxial layers (Kröger et al., 2011[Kröger, I., Stadtmüller, B., Kleimann, C., Rajput, P. & Kumpf, C. (2011). Phys. Rev. B, 83, 195414.]), multilayers (Yakunin et al., 2014[Yakunin, S. N., Makhotkin, I. A., van de Kruijs, R. W. E., Chuev, M. A., Pashaev, E. M., Zoethout, E., Louis, E., Seregin, S. Y., Subbotin, I. A., Novikov, D. V., Bijkerk, F. & Kovalchuk, M. V. (2014). J. Appl. Phys. 115, 134303.]), Langmuir–Blodgett films (Bedzyk et al., 1989[Bedzyk, M. J., Bommarito, G. M. & Schildkraut, J. S. (1989). Phys. Rev. Lett. 62, 1376-1379.]; Novikova et al., 2003[Novikova, N. N., Zheludeva, S. I., Konovalov, O. V., Kovalchuk, M. V., Stepina, N. D., Myagkov, I. V., Godovsky, Y. K., Makarova, N. N., Tereschenko, E. Yu. & Yanusova, L. G. (2003). J. Appl. Cryst. 36, 727-731.]) and shallow ion implant profiles (Hönicke et al., 2010[Hönicke, P., Beckhoff, B., Kolbe, M., Giubertoni, D., van den Berg, J. & Pepponi, G. (2010). Anal. Bioanal. Chem. 396, 2825-2832.]), among others. However, if nanoscale devices, e.g. light-trapping structures in solar cells (Kröger et al., 2011[Kröger, I., Stadtmüller, B., Kleimann, C., Rajput, P. & Kumpf, C. (2011). Phys. Rev. B, 83, 195414.]), field emitter arrays (Fletcher et al., 2013[Fletcher, P. C., Mangalam, V. K. R., Martin, L. W. & King, W. P. (2013). J. Vac. Sci. Technol. B, 31, 021805.]) and nanorods (Malerba et al., 2015[Malerba, M., Alabastri, A., Miele, E., Zilio, P., Patrini, M., Bajoni, D., Messina, G. C., Dipalo, M., Toma, A., Proietti Zaccaria, R. & De Angelis, F. (2015). Sci. Rep. 5, 16436.]), are to be characterized, the calculation of the XSW is more complex. The depth atomic distribution profiles of such structures can still be analyzed in the framework of the conventional 1D XSW method with use of the effective layer approximation (Kennedy et al., 1999[Kennedy, G. P., Buiu, O. & Taylor, S. (1999). J. Appl. Phys. 85, 3319-3326.]). In this approximation, the atomic concentration distribution is averaged along the lateral directions. But this approach does not take into account the diffraction on the lateral structure of the sample and is therefore only applicable in the case of randomly distributed objects. Inherently, the information about lateral distribution is lost within the effective layer approach, and the effective atomic concentration profile can never fully explain the properties of such 2D or 3D devices.

In recent works (Soltwisch et al., 2018[Soltwisch, V., Hönicke, P., Kayser, Y., Eilbracht, J., Probst, J., Scholze, F. & Beckhoff, B. (2018). Nanoscale, 10, 6177-6185.]; Dialameh et al., 2018[Dialameh, M., Ferrarese Lupi, F., Hönicke, P., Kayser, Y., Beckhoff, B., Weimann, T., Fleischmann, C., Vandervorst, W., Dubček, P., Pivac, B., Perego, M., Seguini, G., De Leo, N. & Boarino, L. (2018). Phys. Status Solidi A, 215, 1700866.]), the sensitivity of GIXRF to the lateral distribution of atomic concentration in 2D and 3D structures of periodically arranged gratings and nanocolumns has been experimentally demonstrated. To achieve such sensitivity, a new experimental scheme has been employed, where measurements are carried out under different grazing-incidence and azimuthal-orientation angles. The incidence angle is typically varied within the range from zero to about 3αc where the critical angle αc is calculated for the bulk material. This procures formation of the standing wave in the structured layer. Variation of the azimuthal angle changes the distribution of XSW also in the lateral direction. The optical matrix method (Gibaud & Hazra, 2000[Gibaud, A. & Hazra, S. (2000). Curr. Sci. 78, 1467-1477.]) used for the analysis in conventional XSW (Yakunin et al., 2014[Yakunin, S. N., Makhotkin, I. A., van de Kruijs, R. W. E., Chuev, M. A., Pashaev, E. M., Zoethout, E., Louis, E., Seregin, S. Y., Subbotin, I. A., Novikov, D. V., Bijkerk, F. & Kovalchuk, M. V. (2014). J. Appl. Phys. 115, 134303.]) does not allow analysis of the lateral distribution of the XSW.

This problem of GIXRF data analysis for well ordered structures has been addressed by Soltwisch et al. (2018[Soltwisch, V., Hönicke, P., Kayser, Y., Eilbracht, J., Probst, J., Scholze, F. & Beckhoff, B. (2018). Nanoscale, 10, 6177-6185.]), where the 2D structure of a lamellar Si3N4 grating has been analyzed. The experimentally measured GIXRF curves were analyzed by solving Maxwell's equations by means of a finite-element method (FEM) (Pomplun et al., 2007[Pomplun, J., Burger, S., Zschiedrich, L. & Schmidt, F. (2007). Phys. Status Solidi B, 244, 3419-3434.]). However, applicability of FEM is limited due to its high demand in computational effort. It quickly increases with the increase of the incident photon energy, the size and the dimensionality of the structure. The FEM simulations for the experiments on the 3D Cr nanocolumns published by Dialameh et al. (2018[Dialameh, M., Ferrarese Lupi, F., Hönicke, P., Kayser, Y., Beckhoff, B., Weimann, T., Fleischmann, C., Vandervorst, W., Dubček, P., Pivac, B., Perego, M., Seguini, G., De Leo, N. & Boarino, L. (2018). Phys. Status Solidi A, 215, 1700866.]) for instance are practically unrealisable.

Thus, in this study we provide an alternative approach for the calculation of the XSW field intensities within regular nanostructures by deriving semi-analytic equations based on the dynamical diffraction theory. We derive the solution of the Sherman equation (Sherman, 1955[Sherman, J. (1955). Spectrochim. Acta, 7, 283-306.]; Hönicke et al., 2010[Hönicke, P., Beckhoff, B., Kolbe, M., Giubertoni, D., van den Berg, J. & Pepponi, G. (2010). Anal. Bioanal. Chem. 396, 2825-2832.]) for the GIXRF intensity induced by XSW in the 3D periodic structure in linear-algebraic form. In order to test the new computational scheme, we perform numerical simulations for the same 2D lamellar grating as published by Soltwisch et al. (2018[Soltwisch, V., Hönicke, P., Kayser, Y., Eilbracht, J., Probst, J., Scholze, F. & Beckhoff, B. (2018). Nanoscale, 10, 6177-6185.]) and compare them with the results of FEM simulations and measurements. The semi-analytical nature of the derived equations allowed us to strongly reduce the computational effort, and to perform analysis of GIXRF also from a 3D nanostructured surface for the first time using the experimental data previously published by Dialameh et al. (2018[Dialameh, M., Ferrarese Lupi, F., Hönicke, P., Kayser, Y., Beckhoff, B., Weimann, T., Fleischmann, C., Vandervorst, W., Dubček, P., Pivac, B., Perego, M., Seguini, G., De Leo, N. & Boarino, L. (2018). Phys. Status Solidi A, 215, 1700866.]).

2. Theory

In Sections 2.1[link]–2.3[link] we consider the theoretical background of the dynamical diffraction theory in the many-beam approximation (Mikulík & Baumbach, 1999[Mikulík, P. & Baumbach, T. (1999). Phys. Rev. B, 59, 7632-7643.]) [in the literature also referred to as the rigorous coupled-wave analysis (Chateau & Hugonin, 1994[Chateau, N. & Hugonin, J.-P. (1994). J. Opt. Soc. Am. A, 11, 1321-1331.])]. In Section 2.5[link] we derive the solution of the Shermann equation in linear-algebraic form, which will further allow us to calculate GIXRF intensities of 2D and 3D structures.

2.1. Many-beam dynamical diffraction theory

The experimental geometry used by Dialameh et al. (2018[Dialameh, M., Ferrarese Lupi, F., Hönicke, P., Kayser, Y., Beckhoff, B., Weimann, T., Fleischmann, C., Vandervorst, W., Dubček, P., Pivac, B., Perego, M., Seguini, G., De Leo, N. & Boarino, L. (2018). Phys. Status Solidi A, 215, 1700866.]) and Soltwisch et al. (2018[Soltwisch, V., Hönicke, P., Kayser, Y., Eilbracht, J., Probst, J., Scholze, F. & Beckhoff, B. (2018). Nanoscale, 10, 6177-6185.]) for the GIXRF measurements is shown in Fig. 1(a)[link]. An X-ray beam impinges onto a sample surface under the grazing-incidence angle α and azimuthal angle ϕ. The excited fluorescence emission is measured using an energy-dispersive silicon drift detector D. To simulate the fluorescence intensity from the sample, the near-field (NF) distribution within the nanostructure must be calculated. The problem of the NF calculation is formulated by the Helmholtz equation,

[\left(\Delta+k_0^{\,2}\right)E({\bf{r}}) = -k_0^{\,2}\chi({\bf r})\,E({\bf r}). \eqno(1)]

Here, for simplicity we consider the Helmholtz equation in a scalar approximation, as the effect of polarization is negligible in grazing-incidence geometry in the X-ray spectral range; [E(\bf r)] is the electric field, the sample structure is represented by the dielectric susceptibility function χ(r) and k0 = 2π/λ is the wavenumber of the incident beam with wavelength λ. The Helmholtz equation can be solved using the FEM, kinematical diffraction theory or dynamical diffraction theory. With FEM being computationally challenging, and kinematical theory not sufficiently precise under grazing-incidence conditions (Mikulík & Baumbach, 1999[Mikulík, P. & Baumbach, T. (1999). Phys. Rev. B, 59, 7632-7643.]), we further consider the dynamical diffraction theory. Furthermore, to take into account the lateral structure of the sample one needs to consider the dynamical diffraction theory in the many-beam approximation (MBDDT).

[Figure 1]
Figure 1
(a) Sketch of the experimental geometry of GIXRF. D: energy-dispersive silicon drift detector; α: angle of incidence; k0: wave vector of the incident beam. (b) Sketch of a 2D structure (grating: periodic in x direction, finite depth in z direction), with the azimuthal rotation angle ϕ. (c) Sketch of a 3D periodic structure (nanocolumns: periodic in x and y directions, finite depth in z direction).

In the dynamical diffraction theory, equation (1)[link] is solved assuming that NF is represented as a Bloch wave,

[E({\bf r}) = \sum_{\bf h} E_{\bf h}(z) \exp\left(i{\bf k}_{{\bf h}\parallel} \cdot {\bf r}\right), \eqno(2)]

and the structure is represented as the Fourier series,

[\chi({\bf r}) = \sum_{\bf h} \chi_{\bf h} \exp(i{\bf h} \cdot {\bf r}), \eqno(3)]

where χh is the Fourier component,

[\chi_{{\bf h}} = {{1} \over {\Omega}} \int\!\!\!\int\chi(x,y) \, \exp\left(-i\,{\bf h}\cdot{\bf r}\right)\,{\rm{d}}S. \eqno(4)]

Here, integration is taken over the unit-cell area Ω for the corresponding reciprocal space vector,

[{\bf h}_{x,y} = {{ 2\pi n_{x,y} } \over { D_{x,y} }} \, {\bf e}_{x,y}, \eqno(5)]

with nx,y the order of the diffraction index and Dx, y the periods along x or y directions. The parallel component of the wavevector of the hth diffraction order [{\bf k}_{{\bf h}\parallel}] = [{\bf k}_{0\parallel} + {\bf h}] is translationally invariant along the z direction; i.e. [{\bf k}_{{\bf h}\parallel}] is constant across all medias in a layered system for given h, while the vertical component is generally different in each medium and defined with the spherical dispersion equation,

[q_{{\bf h}z}^2 = k_0^{\,2} \left(1+\chi_0\right) - k_{{\bf h}\parallel}^{\,2}. \eqno(6)]

This equation is derived assuming that diffraction scattering is an elastic process: kh = (1 + χ0)k0 and assuming translational invariance of [{\bf k}_{{\bf h}\parallel}] mentioned above. Finally, substituting equations (2)[link], (3)[link] and (6)[link] into (1)[link], considering a property of the Fourier components, i.e. [\chi_{{\bf g}}\exp(i{\bf h}\cdot{\bf r})] = [\chi_{{\bf g} - {\bf h}}], results in a system of inhomogeneous linear ordinary differential equations (ODEs) of second-order,

[q_{{\bf h}z}^2 \, E_{\bf h}(z) + {{{\rm{d}}^2}\over{{\rm{d}}z^2}} E_{\bf h}(z) + k_0^{\,2} \sum_{{\bf g}\,\neq\,{\bf h}} \, E_{\bf g}(z) \chi_{{\bf g} - {\bf h}} = 0. \eqno(7)]

The general solution of such a system of ODEs is a linear combination of particular solutions of corresponding homogeneous ODEs, where the nth particular solution has the form of a standing wave with amplitudes Tn and Rn. Thus, the hth solution of equation (7)[link] has the form

[E_{\bf h}(z) = \sum_n \left[ T_n\exp\left(-ik_{nz}z\right)+R_n\exp\left(ik_{nz}z\right) \right] E_{hn}, \eqno(8)]

with linear combination coefficients Ehn. Therefore, the distribution of the NF is defined with equations (2)[link] and (8)[link]. Thus, the problem of NF calculation is reduced to finding kn,z, Ehn, Tn and Rn.

2.2. Characteristic equation

In this section we discuss the calculation of kn,z and Ehn. Variable kn,z has a physical meaning as the vertical component of the wavevector [see equation (8)[link]]. It defines the phase of the standing wave in the structured layer. One can assume that kn,z is defined with spherical dispersion kz = qz; however, under that assumption equation (7)[link] has no solutions. Therefore, values of kn,z deviate from spherical dispersion. To calculate the precise value of kn,z in the structured layer one can substitute equation (8)[link] into (7)[link]. The result is represented as the eigenvalues–eigenvectors problem,

[\left({\bf A} - k_{zn}^{\,2} {\bf I} \right) \,{\bf E}_n = 0, \eqno(9)]

where kzn 2 is an eigenvalue of matrix [{\bf A}] and En is an eigenvector composed of the coefficients Ehn from equation (8)[link]: [{\bf E}_n] = [(\ldots E_{-1,n},E_{0,n},E_{1,n}\ldots)^T]; [{\bf A}] is of the form

[{\bf A} = k_0^{\,2}{\bf C} - {\bf X}. \eqno(10)]

Matrix [{\bf C}] is the Toeplitz circulant matrix,

[{\bf C} = \left[ \matrix{ \ddots & & & & \cr & \chi_0 & \chi_{-1} & \chi_{-2} & \cr & \chi_1 & \chi_0 & \chi_{-1} & \cr & \chi_{2} & \chi_{1} & \chi_{0} & \cr & & & & \ddots \cr } \right], \eqno(11)]

and [{\bf X}] is the diagonal matrix with diagonal [(\ldots -k^2_{-1,\parallel},-k^2_{0,\parallel},-k^2_{1,\parallel} \ldots)]. Circulant matrices have a remarkable property: with increasing circulant matrix size, their eigenvalues asymptotically approach the exact values for an infinite matrix (Gray, 2006[Gray, R. M. (2006). Found. Trends Commun. Information Theory, 2, 155-239.]). Therefore, one can use a finite amount of Fourier components in equation (3)[link] to approximate the exact solution of equation (7)[link]. Consider a set of 2N+1 Fourier components [\{\chi_{-N}, \chi_{-N+1}, \ldots, \chi_0, \ldots \chi_{N-1}, \chi_{N} \}]. These Fourier components constitute a circulant matrix C of size [{\bf C} \in {\cal C}^{\,\,M \times M}], where M = N + 1. Solving the characteristic equation (9)[link] will give M eigenvalue–eigenvector pairs.

2.3. Boundary conditions

In this section we calculate the transmission Tn and reflection Rn amplitudes. Consider a sample as a stratified medium, consisting of layers. Tn and Rn are calculated in each layer using continuity conditions of the electric field and its first derivative. The continuity conditions (Born & Wolf, 2013[Born, M. & Wolf, E. (2013). Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. Elsevier.]) for the jth and (j + 1)th pair of layers can be written in a matrix form,

[{\bf P}^{\,(\,j)} \left[ \matrix{ {\bf T}^{\,(\,j)} \cr {\bf R}^{\,(\,j)} } \right] = {\bf P}^{\,(\,j+1)}{\bf Q}^{\,(\,j+1)} \left [\matrix{ {\bf T}^{\,(\,j+1)} \cr {\bf R}^{\,(\,j+1)} } \right]. \eqno(12)]

Here T and R are vectors composed of amplitudes Tn and Rn,

[T = \left(T_{-N/2}, T_{-N/2+1}, \ldots, T_0, \ldots T_{N/2-1}, T_{N/2} \right)^T. \eqno(13)]

Equation (12)[link] links amplitudes T(j), R(j) at the interface between the (j − 1)th and the jth layer, and amplitudes T(j+1), R(j+1) at the interface between the jth and the (j + 1)th layer. Matrix [{\bf P}] is the refraction matrix. For a structured layer it has the form

[{\bf P} = \left [\matrix{ {\bf E} & {\bf E} \cr -{\bf Ek}_z & {\bf E k}_z \cr } \right], \eqno(14)]

and for a homogeneous layer

[{\bf P} = \left [\matrix{ {\bf I} & {\bf I} \cr -{\bf k}_z & {\bf k}_z \cr } \right]. \eqno(15)]

Here, the matrix [{\bf E}] is composed of columns of eigenvectors and matrix kz is a diagonal matrix filled with kz,n. Refraction matrix [{\bf P}] is a 2 × 2 block matrix, thus [{\bf P} \in {\cal C}^{\,\,2M\times2M}]. Finally [{\bf Q}] is the propagation matrix,

[{\bf Q} = \left [\matrix{ {\bf Q}^+ & 0 \cr 0 & {\bf Q}^- \cr } \right], \eqno(16)]

where [{\bf Q}^\pm] are diagonal matrices with corresponding diagonals

[\left[ \ldots \exp\left(\mp ik_{-h,z}d_j\right), \exp\left(\mp ik_{0,z}d_j\right), \exp\left(\mp ik_{h,z}d_j\right) \ldots \right], \eqno(17)]

where dj is the thickness of jth layer. Although these equations can be used to calculate Tj and Ri, solving equation (12)[link] might be problematic due to the poorly conditioned transmission matrix in the case of a sufficiently large thickness of the sample and/or in the case of a sufficiently high number of Fourier components used in the calculation.

2.4. Numerical stability

The problem of numerical stability in the matrix formalism of dynamical diffraction theory was considered by Stepanov et al. (1998[Stepanov, S. A., Kondrashkina, E. A., Köhler, R., Novikov, D. V., Materlik, G. & Durbin, S. M. (1998). Phys. Rev. B, 57, 4829-4841.]). There the problem of numerical stability has been solved for the dynamical diffraction theory in a two-beam approximation (only χ−1, χ0 and χ1 have been taken into account). It has been solved by dividing matrices into 2 × 2 block matrices and solving equation (12)[link] separately for each block matrix by using recurrent formula. Although the recurrent matrix equations of Stepanov et al. (1998[Stepanov, S. A., Kondrashkina, E. A., Köhler, R., Novikov, D. V., Materlik, G. & Durbin, S. M. (1998). Phys. Rev. B, 57, 4829-4841.]) were derived for the two-beam case, they are generally applicable to the many-beam case. For brevity, we present these equations explicitly written for a three-layer system (see Fig. 2[link]), which is relevant to the experimental data we will consider further.

[Figure 2]
Figure 2
Sketch of the three-layer model. Arrows schematically depict the direction of propagation of the plane waves. Amplitudes of the plane waves are assembled into T and R vectors. [\widetilde{{\bf T}}] and [\widetilde{{\bf R}}] are amplitudes defined at the upper interface of the layer. T and R are defined at the bottom interface of the layer.

The continuity conditions for such a three-layer structure (vacuum – structured layer – substrate) are represented by the system of linear equations

[\eqalign{ \left[ \matrix{ {\bf T}^{\,({\rm v})} \cr {\bf R}^{({\rm v})} } \right] & = \left({\bf{P}}^{({\rm v})}\right)^{-1} {\bf{P}}^{({\rm g})}{\bf{Q}}^{({\rm g})} \left[\matrix{ {\bf T}^{\,({\rm g})} \cr {\bf R}^{({\rm g})} } \right], \cr \left [\matrix{ {\bf T}^{\,({\rm g})} \cr {\bf R}^{({\rm g})} } \right] & = \left({\bf P}^{({\rm g})}\right)^{-1} {\bf P}^{({\rm s})} \left [\matrix{{\bf T}^{\,({\rm s})} \cr {\bf R}^{({\rm s})} } \right]. } \eqno(18)]

One can rewrite that system as follows,

[\eqalign{ \left[ \matrix{ {\bf T}^{\,({\rm g})} \cr {\bf R}^{({\rm v})} } \right] & = {\bf M}^{\rm vg} \left [\matrix{ {\bf T}^{\,({\rm v})} \cr {\bf R}^{({\rm g})} } \right], \cr \left [\matrix{ {\bf T}^{\,({\rm s})} \cr {\bf R}^{({\rm g})} } \right] & = {\bf M}^{\rm gs} \left [\matrix{ {\bf T}^{\,({\rm g})} \cr {\bf R}^{({\rm s})} } \right]. } \eqno(19)]

Here matrix [{\bf M}] has the form of a block matrix,

[{\bf M} = \left [\matrix{ {\bf Q}^{-} {\bf V}_{11}^{\,-1} & -{\bf Q}^{-} {\bf V}_{11}^{\,-1} {\bf V}_{12}{\bf Q}^- \cr {\bf V}_{21} {\bf V}_{11}^{\,-1} & {\bf V}_{22}{\bf Q}^{-} - {\bf V}_{21} {\bf V}_{11}^{\,-1} {\bf V}_{12}{\bf Q}^- } \right], \eqno(20)]

where Vij is a matrix element of 2 × 2 block matrix V(vg, gs) = (P(v, g))−1P(g, s). Note that this equation does not include [{\bf Q}^+] whose elements are growing exponentially with respect to the thickness of the structured layer. Hence this matrix is numerically stable. Amplitudes [{\bf T}^{\,({\rm v})}] represent the incident beam, therefore

[{\bf T}^{\,({\rm v})} = (\ldots 0,1,0 \ldots)^T. \eqno(21)]

Additionally, for a sufficiently thick substrate we can assume

[{\bf R}^{({\rm s})} = (\ldots 0,0,0 \ldots)^T. \eqno(22)]

Taking into account these considerations, we derive equations for amplitudes in the structured layer,

[{\bf R}^{({\rm g})} = \left({\bf I} - {\bf M}_{21}^{({\rm gs})}{\bf M}_{12}^{({\rm vg})} \right)^{-1} {\bf M}_{21}^{({\rm gs})} {\bf M}_{11}^{({\rm vg})} {\bf T}^{\,({\rm v})} \eqno(23)]

and

[{\bf T}^{\,({\rm g})} = \left({\bf I} - {\bf M}_{12}^{({\rm vg})}{\bf M}_{21}^{({\rm gs})} \right)^{-1} {\bf M}_{11}^{({\rm vg})} {\bf T}^{\,({\rm v})}. \eqno(24)]

These amplitudes are calculated at the interface between the structured layer and the substrate (see Fig. 2[link]). One can calculate amplitudes at the vacuum-structured layer interface using

[\left [\matrix{ \widetilde{{\bf T}}^{({\rm g})} \cr \widetilde{{\bf R}}^{({\rm g})} } \right] = \left({\bf P}^{({\rm g})}\right)^{-1} {\bf P}^{({\rm v})} \left [\matrix{ {\bf T}^{\,({\rm v})} \cr {\bf R}^{({\rm v})} } \right]. \eqno(25)]

Finally, we need to rewrite equation (8)[link],

[E_{\bf h}(z) = \sum_n \left[ \widetilde T_n\exp\left(-ik_{nz}z\right) + R_n\exp\left(ik_{nz}[z+d]\right) \right] E_{hn}. \eqno(26)]

Here, both exponents decrease with respect to the depth, providing numerical stability.

2.5. X-ray fluorescence intensity

The fluorescence intensity Y can be calculated using the Sherman equation (Sherman, 1955[Sherman, J. (1955). Spectrochim. Acta, 7, 283-306.]), adapted for GIXRF (Hönicke et al., 2010[Hönicke, P., Beckhoff, B., Kolbe, M., Giubertoni, D., van den Berg, J. & Pepponi, G. (2010). Anal. Bioanal. Chem. 396, 2825-2832.]),

[Y \,\propto\, G(\alpha) \int\!\!\int\!\!\int \big|E({\bf r})\big|^2 p({\bf r}) \exp(-\mu \rho z) \,{\rm{d}}{\bf r}, \eqno(27)]

where [p({\bf r})] describes the density distribution of fluorescent atoms in the structure, and G(α) is the geometrical factor (Beckhoff, 2008[Beckhoff, B. (2008). J. Anal. At. Spectrom. 23, 845-853.]; Li et al., 2012[Li, W., Zhu, J., Ma, X., Li, H., Wang, H., Sawhney, K. J. & Wang, Z. (2012). Rev. Sci. Instrum. 83, 053114.]; Lubeck et al., 2013[Lubeck, J., Beckhoff, B., Fliegauf, R., Holfelder, I., Hönicke, P., Müller, M., Pollakowski, B., Reinhardt, F. & Weser, J. (2013). Rev. Sci. Instrum. 84, 045106.]). The integral is taken over the area of the elementary cell. The exponential term [\exp(-\mu\rho{z})] in equation (27)[link] takes into account the self-absorption of emitted fluorescent photons. Here, μ is the absorption coefficient and ρ is the effective density of the absorbing media. The integral in equation (27)[link] can be separated as follows [further, for brevity we do not explicitly write the multiplicative term G(α)],

[\eqalignno{ Y \,\propto\, {}& \sum_{{\bf g},{\bf h}} \int\!\!\!\int p({\bf r}_{\parallel}) \exp\left[i\left({\bf k}_{{\bf g}\parallel} - {\bf k}_{{\bf h}\parallel}^*\right) \cdot {\bf r}\right]\,{\rm{d}}x\,{\rm{d}}y \cr& \times \sum_{m,n} \int \left[T_m\exp\left(-ik_{mz}z\right) + R_m\exp\left(ik_{mz}z\right) \right] \cr& \times \left[ T_n^{\,*}\exp\left(ik^{\,*}_{nz}z\right) + R_n^{\,*}\exp\left(-ik^{\,*}_{nz}d\right) \right] \cr& \times E_{gm}E_{hn}^{\,*} \exp(-\mu \rho z)\,{\rm{d}}z. &(28) }]

Such an integral separation imposes a restriction on the numeric density function: it must not be dependent on the z coordinate, p(r) ≡ p(rα); i.e. equation (28)[link] can only be used in cases when fluorescent atoms are distributed homogeneously along the z direction. Distribution in the xy-plane can be arbitrary. In the case of an inhomogeneous vertical distribution, one can discretize the structure along the z direction as a stack of sublayers and calculate equation (28)[link] for each sublayer.

Equation (27)[link] was rewritten in the form of equation (28)[link], so it can be conveniently represented in a linear algebraic language. The fluorescence intensity can be expressed as the sum of matrix elements [Y \,\propto\, \sum_{g,h}F_{gh}] of the matrix,

[{\bf F} = \boldPhi \circ \left({\bf E} \boldPsi {\bf E}^{\rm H} \right). \eqno(29)]

Here [\circ] represents an element-wise (Hadamard) multiplication and H represents a Hermitian transpose. Elements of matrix [\boldPhi] have the form

[\Phi_{hg} \,\equiv\, \int\limits_{-{D_x}/{2}}^{{D_x}/{2}} \!\!\! {\rm{d}}x \int\limits_{-{D_y}/{2}}^{{D_y}/{2}} \!\!\! {\rm{d}}y \,\,\, p({\bf r}_{\parallel}) \exp\left[i\left({\bf k}_{{\bf g}\parallel} - {\bf k}_{{\bf h\parallel}}^*\right) \cdot {\bf r}\right], \eqno(30)]

and elements of matrix [\boldPsi] have the form

[\eqalignno{ \Psi_{mn} \,\equiv\, {}& T_mT_n^{\,*} U(-k_{zm}+k_{zn}^{\,*}) +T_mR_n^{\,*} U(-k_{zm}-k_{zn}^*) \cr& +R_mT_n^{\,*} U(k_{zm}+k_{zn}^{\,*}) +R_mR_n^{\,*} U(k_{zm}-k_{zn}^{\,*}), &(31) }]

where

[U(q) \,\equiv \int\limits_{-d}^0 \exp(iqz) \exp(-\mu\rho{z}) \,{\rm{d}}z. \eqno(32)]

The [\boldPhi] matrix takes into account the distribution of fluorescent atoms and the electric field distribution in the lateral direction and the [\boldPhi] matrix takes into account photon absorption and the electric field distribution in the vertical direction. Equation (29)[link] allows the integral in equation (27)[link] to be calculated analytically, which is much more computationally efficient compared with the numerical integration.

3. Numerical simulations

3.1. 2D structure: Si3N4 lamellar grating

Here, we consider a 2D lamellar Si3N4 grating prepared using electron beam lithography. The original study with experimental data and numerical simulation of GIXRF intensity by means of FEM has been published by Soltwisch et al. (2018[Soltwisch, V., Hönicke, P., Kayser, Y., Eilbracht, J., Probst, J., Scholze, F. & Beckhoff, B. (2018). Nanoscale, 10, 6177-6185.]). The grating has a nominal period of Dx = 100 nm, the thickness of the structured layer is d = 90 nm and the line width is 40 nm.

The GIXRF measurements were carried out at the plane-grating monochromator (PGM) beamline (Senf et al., 1998[Senf, F., Flechsig, U., Eggenstein, F., Gudat, W., Klein, R., Rabus, H. & Ulm, G. (1998). J. Synchrotron Rad. 5, 780-782.]) for undulator radiation at the PTB laboratory (Beckhoff et al., 2009[Beckhoff, B., Gottwald, A., Klein, R., Krumrey, M., Müller, R., Richter, M., Scholze, F., Thornagel, R. & Ulm, G. (2009). Phys. Status Solidi B, 246, 1415-1434.]) of the BESSY II electron storage ring.

An incident photon energy of 520 eV was used. The PGM monochromator provides an energy resolution of δE/E < 5 × 10−4 in this energy range. The GIXRF intensities were obtained for the N-Kα fluorescence emission under various incidence angles α and azimuthal sample orientation angles ϕ [see Fig. 1(b)[link]]; ϕ = 0° corresponds to the conical orientation (Goray et al., 2018[Goray, L., Jark, W. & Eichert, D. (2018). J. Synchrotron Rad. 25, 1683-1693.]) of the sample grating. The recorded spectra from the silicon drift detector were deconvoluted using detector response functions in order to isolate the fluorescence signal from N-Kα from other spectral contributions. Further corrections, to take into account the detection efficiency and the geometrical factor (effective solid angle), were applied [see Soltwisch et al. (2018[Soltwisch, V., Hönicke, P., Kayser, Y., Eilbracht, J., Probst, J., Scholze, F. & Beckhoff, B. (2018). Nanoscale, 10, 6177-6185.]) for further details].

Best-fit simulations obtained by a sequential least-squares optimization algorithm and experimental GIXRF data are shown in Fig. 3[link] for various azimuthal orientation angles ϕ. For the simulation we use a simple box model in which the grating lines are treated as an array of boxes on top of the substrate [see Fig. 1(b)[link]]. Thus, the medium is divided into three areas: the vacuum, the structured layer in which the boxes are located, and the substrate. Within the box model, the sidewalls of the grating lines are considered to be parallel while the actual grating has a sidewall tilt angle. Based on the reconstruction of Soltwisch et al. (2018[Soltwisch, V., Hönicke, P., Kayser, Y., Eilbracht, J., Probst, J., Scholze, F. & Beckhoff, B. (2018). Nanoscale, 10, 6177-6185.]), this angle is not greater than β = 4°. In terms of the model it means that the Fourier transform in equation (3)[link] is changing along the z axis. To compensate for that in the simulations, within the one-layer model averaged Fourier components have been used, i.e. [\langle \chi_h \rangle = \chi_h\exp(-h^2\sigma^2/2)], with σ defined as half the projection of the sidewall on the x axis: [\sigma \equiv d\arctan{(\beta)}/2]. The best-fit line width (defined as the half-height width) is Dl = 39 nm and the best-fit sidewall tilt angle is β = 5°.

[Figure 3]
Figure 3
N-Kα GIXRF intensity, measured for various azimuthal orientation angles: (a) ϕ = 0° – conical, (b) ϕ = 0.2°, (c) ϕ = 1° and (d) ϕ = 3°. Red lines: numerical simulation; gray markers: experimental values.

Another feature of the actual sample that must be considered in the simulations is the effect of oxidation of surface and line edges. It affects the actual structure such that the concentration of fluorescent N atoms at the top part and at the line edges is strongly reduced. In the one-layer model, oxidation of the surface can be effectively incorporated by changing the integration limits in equation (32)[link], such that the integration in equation (32)[link] is taken only over a range where fluorescent atoms are present.

The best agreement with the experimental data was obtained with an effective surface layer thickness of dt = 3.3 nm at the top of the lines. The best fit suggests that N is not diluted at the line edges, since the reconstructed parameter of the effective edge thickness of the edges is ds = 0 nm. We note that this value ds is correlated with σ used in averaging of the Fourier components, and thus may not be representative. Also, note that these values only describes surface effects in terms of the absence of fluorescent N atoms, ignoring the gradual change in stoichiometry throughout the surface and the edges. It also neglects the change of optical properties of the structure due to oxidation. Best-fit parameter of the grating height, excluding the effective surface layer, is d = 88.7 nm. The average density of the line is [\rho_{{\rm Si}_3{\rm N}_4}] = 2.8 g cm−3 and the density of the substrate is ρSi = 2.22 g cm−3.

The best-fit model and experimental data are qualitatively in good agreement. Qualitative agreement is also apparent on the GIXRF intensity (α, ϕ) maps shown in Fig. 4[link]. A full set of 48 experimental GIXRF curves taken along different azimuthal angles ϕ (from 0° to 2°) was interpolated on a (α, ϕ) grid [see Fig. 4(a)[link]]. The theoretical GIXRF map was calculated on the same (α, ϕ) grid using best-fit parameters from the data presented in Fig. 4[link].

[Figure 4]
Figure 4
Comparison between the experimental GIXRF N-Kα map (a) of the Si3N4 lamelar grating measured with the incidence photon energy E = 520 eV and the simulated GIXRF map (b) based on a best fit model. (c) Resonant lines in GIXRF map for Si3N4 grating structure, caused by interference between the reflected beam and the mth order of diffraction.

One can note a distinctive feature on the GIXRF map – resonant lines, which are visible both on the experimental data in Fig. 4(a)[link] and in the numerical simulations [see Fig. 4(b)[link]]. As a visual aid to notice these lines one can refer to the sketch in Fig. 4(c)[link]. In Fig. 4(c)[link] the position of the resonant lines is marked with black contour lines.

We assume that these lines are due to the interference between the reflected beam (zeroth order of diffraction) and a diffracted beam (mth order of diffraction). Therefore, the resonant lines must satisfy the Laue condition, which for this geometry can be formulated as kx 2+kz 2 = (kx+h)2. This formula geometrically corresponds to the Ewald sphere. For convenience we rewrite this equation in terms of the incidence and azimuthal angle,

[\sin\phi = {{\sin^2\alpha-\gamma^2}\over{2\gamma\cos\alpha}}, \eqno(33)]

where γ = λm/Dx. The contour lines in Fig. 4(c)[link] were calculated using this equation. Note that the resonant lines depend only on the lateral period of the structure Dx and the wavelength λ [see equation (33)[link]] – no other geometrical parameters are involved. Due to their explicit dependence on only the period of the structure, such lines might be used in the analysis of experimental data as a reference, to determine the lateral period of the structure, without needing a full structure reconstruction through model simulations. This equation can also be used to determine the angular range of interest for the measurement as the most interesting areas to probe are around these lines.

3.2. 3D structure: Cr nanocolumns

In this section we consider a periodic 3D nanocolumnar structure of Cr, manufactured using electron beam lithography (Altissimo, 2010[Altissimo, M. (2010). Biomicrofluidics, 4, 026503.]) on top of a SiO2 substrate. The structure of the sample is a regular square grid of box-shaped columns [see Fig. 1(c)[link]] on a substrate, with 300 nm × 300 nm lateral box dimensions and a Dx = Dy = 1 µm grid. The nominal height of the nanocolumns is d = 25 nm.

GIXRF measurements were carried out at the four crystal monochromator (FCM) beamline (Krumrey & Ulm, 2001[Krumrey, M. & Ulm, G. (2001). Nucl. Instrum. Methods Phys. Res. A, 467-468, 1175-1178.]) in the PTB laboratory (Senf et al., 1998[Senf, F., Flechsig, U., Eggenstein, F., Gudat, W., Klein, R., Rabus, H. & Ulm, G. (1998). J. Synchrotron Rad. 5, 780-782.]) of the BESSY II storage ring and reported by Dialameh et al. (2018[Dialameh, M., Ferrarese Lupi, F., Hönicke, P., Kayser, Y., Beckhoff, B., Weimann, T., Fleischmann, C., Vandervorst, W., Dubček, P., Pivac, B., Perego, M., Seguini, G., De Leo, N. & Boarino, L. (2018). Phys. Status Solidi A, 215, 1700866.]). The incident photon energy was set to E = 7 keV with an energy resolution of δE/E < 5 × 10−4. Numerical simulations are carried out similarly to those in Section 3.1[link]. The GIXRF experimental data and the best-fit obtained from dynamical diffraction theory simulations are shown in Fig. 5[link] for a selection of azimuthal angles.

[Figure 5]
Figure 5
Cr-Kα GIXRF intensity curves, measured for various azimuthal orientation angles (ϕ = 0° corresponds to the conical orientation). Red lines: numerical simulation; gray markers: experimental data.

Best-fit model parameters are: lateral period of the structure Dx = Dy = 1 µm, matching the same nominal values; lateral sizes of the nanocolumns are 300 nm × 300 nm, nanocolumn height d = 24 nm. The best-fit model suggests that there is no surface oxidation, dt = 0 nm; however, the effective thickness of the sidewalls is ds = 1.3 nm. The density of the nanocolumn material is equal to the nominal Cr density ρCr ≃ 7.2g cm−3, while the substrate density is [\rho_{{\rm SiO}_2}] = 2.4 g cm−3. Considering the large lateral period Dx = Dy = 1 µm (significantly larger than that of the Si3N4 lamellar grating structure), the sidewalls tilt is negligible, therefore σ = 0 nm, i.e. the best-fit model for the nanocolumn structure implies perfectly parallel sidewalls 〈χh〉 ≡ χh. The experimental GIXRF curves in Fig. 5[link] are in good agreement with numerical simulations.

It is important to note that in the case of grazing-incidence geometry the GIXRF curves calculated for the 3D structure could also be approximated with the use of an effective 2D model, albeit with reduced density. This is because in the grazing-incidence geometry the momentum transfer |ky| ≫ |kx|. In other words, measurements in grazing-incidence geometry are sensitive to the lower frequencies of the Fourier transform of the structure along the x direction and to the higher frequencies along the y direction, while the spacing between nodes in reciprocal space along the kx and ky directions are identical due to the symmetry Dx = Dy of the periodic structure.

Thus, GIXRF curves of Cr nanocolumns can be effectively represented in a first approach as a lamellar Cr grating with reduced density equal to the averaged density of the actual 3D structure. However, a direct comparison between 3D and 2D simulations (Fig. 6[link]) reveals some differences. For higher incident angles above the critical angle of total external reflection, the 2D model (dashed blue line in Fig. 6[link]) yields a monotonous angular dependence, while the experimental GIXRF curve clearly exhibits oscillatory behaviour in that angular range, with a maximum at [\alpha] ≃ 1.5°. In Fig. 6[link], curves are shown only for ϕ = 0°, but this oscillation in the range of higher incidence angles α is present in all experimental curves measured at different azimuthal orientations of the sample (see Fig. 5[link]). We attribute this oscillation to interference due to the periodicity of the structure along the y direction, which becomes more important at higher incident angles since the value of |ky| decreases with increasing incidence angle α and the measurement becomes more sensitive to the lower frequencies of the Fourier transform along the y direction. Such interference mode is not taken into account in the 2D simulations. Additionally, the 3D simulations show resonant peaks at α ≃ 1.15° and α ≃ 1.54° which are not resolved in experimental data. To observe these peaks, measurements with step sizes of δα = 0.01° should be resolved, which is experimentally feasible, as the resolution limit of modern synchrotron sample stage equipment is on the level of 0.001°.

[Figure 6]
Figure 6
Comparison of effective 2D and genuine 3D simulations of GIXRF Cr-Kα curve for 3D Cr nanocolumns structure in conical geometry (ϕ = 0°).

4. Discussion

In Table 1 we compare the structure parameters of the 2D lamellar Si3N4 grating, as reconstructed using the MBDDT simulations described in Section 2.5[link], with the nominal parameters used in fabrication of the grating. The results of of the MBDDT reconstruction are in good agreement with the nominal values.

To further validate the computational scheme described in Section 2.5[link], FEM simulations have also been performed. The FEM simulations were performed using the JCMwave software (Pomplun et al., 2007[Pomplun, J., Burger, S., Zschiedrich, L. & Schmidt, F. (2007). Phys. Status Solidi B, 244, 3419-3434.]) for a box model based on the best-fit parameters in Table 1[link]. JCMwave is a rigorous Maxwell solver, which enables field simulations in structures of arbitrary shape. For the calculation of the finite-element solution, the computational domain is meshed into patches where a number of polynomial ansatz functions are defined. The finite-element side constraint of 4 nm and a polynomial degree of 4 have been used in the simulations and the GIXRF fluorescence intensities were calculated from electric fields as described by Soltwisch et al. (2018[Soltwisch, V., Hönicke, P., Kayser, Y., Eilbracht, J., Probst, J., Scholze, F. & Beckhoff, B. (2018). Nanoscale, 10, 6177-6185.]). A direct comparison of the MBDDT and the FEM simulations is shown in Fig. 7[link]. The GIXRF maps are symmetrical with respect to an axis at ϕ = 0°. In Fig. 7(a)[link] the map on the left is thus showing the MBDDT result, whereas that on the right shows the FEM result. Both simulation results are visually identical. In addition, the relative discrepancy is shown in Fig. 7(b)[link]. Here, the relative discrepancy is defined as

[\varepsilon_{ij} = {{ \left|Y_{ij}^{\,\rm (f)} - Y_{ij}^{\,\rm (m)}\right| }\over{ \max\left\{Y_{ij}^{\,\rm (f)}, Y_{ij}^{\,\rm (m)}\right\} }}, \eqno(34)]

where Yij (f,m) are the GIXRF intensities calculated in each (αi, ϕj) point using the FEM and MBDDT methods, respectively.

Table 1
Comparison of the 2D structure parameters of the Si3N4 lamellar grating as reconstructed by MBDDT with nominal parameters

  Nominal Simulation
Period Dx (nm) 100 100
Line height d (nm) 87 88.7
Line width (nm) 40 39
Effective surface thickness dt (nm) 3.3
Effective edge thickness ds (nm) 0
Sidewalls tilt (°) 4
Line density [\rho_{{\rm Si}_3{\rm N}_4}] (g cm−3) 3.2 2.8
Substrate density ρSi (g cm−3) 2.33 2.22
[Figure 7]
Figure 7
(a) Comparison of GIXRF maps as simulated by MBDDT (left-hand side) and FEM (right-hand side) approaches. (b) Relative discrepancy.

The absolute maximum of the relative discrepancy is 2.4% and the discrepancies are generally higher for the low incidence angles. It should be noted that the precision of the FEM calculation in this angular range may be limited due to the exponential decay of the evanescent waves. In general, the relative discrepancy is on a level of 1% for 80% of the points, proving the validity of the MBDDT approach for such GIXRF simulations.

An unambiguous comparison of the computational efficiency of FEM and MBDDT cannot be performed directly, since it depends on many factors such as photon energy, period of structure, numerical accuracy, and dimensionality of the problem. For the particular simulation shown in Fig. 7[link], the calculation time for one point of the GIXRF map on one CPU is on the level of 10 s for FEM and 0.01 s for MBDDT. Both MBDDT and FEM computations were performed on an NUMA computer with 160 CPUs [Intel(R) Xeon(R) CPU E7-4870 v2 @ 2.30 GHz]. Thus, a calculation of the GIXRF map of 200 × 100 size on this setup takes approximately 20 min and 2 s for FEM and MBDDT, respectively.

In Table 2[link] we compare the MBDDT-derived structure parameters of the Cr nanocolumns with their nominal parameters. A good agreement is obtained, especially since for the current case only a simple box model is used for the numerical simulations to describe the distribution of fluorescent atoms in the structure.

Table 2
Comparison of the 3D structure parameters of the Cr nanocolumns as reconstructed by MBDDT with nominal parameters

  Nominal Simulation
Period Dx, y (µm) 1 1
Column height d (nm) 25 24
Column width (nm) 300 300
Effective surface thickness dt (nm) 0
Effective edge thickness ds (nm) 1.3
Column density ρCr (g cm−3) 7.19 7.2
Substrate density [\rho_{{\rm SiO}_2}] (g cm−3) 2.65 2.4

With the MBDDT approach, it is also possible to take into account a structure with tilted sidewalls and surface oxidation. The model of the structure would needed to be discretized along the z direction, e.g. according to Pisarenco et al. (2016[Pisarenco, M., Quintanilha, R., van Kraaij, M. G. M. M. & Coene, W. M. J. (2016). J. Opt. Soc. Am. A, 33, 610-617.]). The sample can be approximated as a stack of homogeneous and/or structured layers, where each layer can have arbitrary structure parameters with the exception of the period, which must be maintained throughout the whole stack.

The main benefit in applications of the 3D XSW technique to the characterization of nanostructures is its sensitivity to the spatial distribution of the fluorescent atoms within the structure. Although the examples considered in this work exhibited a homogeneous lateral distribution of N and Cr within the grating line and nanocolumn, we can still demonstrate this sensitivity by performing simple calculations.

We use the same model for the Si3N4 lamellar grating as already shown earlier, but now we assume to have dopant atoms localized in a confined volume within the structure as shown in Figs. 8(a)–8(c)[link] instead of being homogeneously distributed, as assumed in Section 3.1[link]. The specific localization of the dopant atoms is depicted by green boxes and the resulting simulated GIXRF maps for the calculated fluorescence signal of the dopant atoms are shown. It can be observed that the corresponding GIXRF maps are highly sensitive to this variation. For an asymmetric distribution of fluorescent atoms, Fig. 8(c)[link], an asymmetry is also observed in the GIXRF maps. One may exploit such asymmetry, e.g. to distinguish chemical compositions of the left and right sidewall of the grating line. This may be useful in, for example, the characterization of gratings fabricated with multi-patterning (Weber et al., 2012[Weber, T., Käsebier, T., Szeghalmi, A., Knez, M., Kley, E. & Tünnermann, A. (2012). Microelectron. Eng. 98, 433-435.]) techniques.

[Figure 8]
Figure 8
Simulation of GIXRF maps for inhomogeneous distribution of fluorescent atoms within the lamellar grating structure. From (a) to (c): sketch of the structure; the green box depicts the localization of the dopant atoms. From (a′) to (c′): corresponding GIXRF maps.

5. Conclusions

A new computational scheme based on the dynamical diffraction theory has been developed and applied for the analysis of GIXRF experiments on 2D and 3D periodic nanostructures. It is capable of simulating GIXRF data from structures with specific element distributions both in-plane as well as in-depth. The computational scheme has been validated with a Maxwell solver based on the finite-element method and benchmarked on GIXRF experimental data obtained from Si3N4 2D lamellar gratings and Cr 3D nanocolumns. The reconstructed geometrical parameters of the lamellar grating derived from the elemental distribution are in good agreement with nominal values, as well as with parameters obtained from a previous study performed using a finite-element method. Furthermore, the parameters of the elemental distribution in the Cr 3D nanocolumns were reconstructed for the first time. A reconstruction of the geometrical parameters of this structure by means of FEM is practically impossible due to the required higher excitation photon energy, the larger period of the structures (and thus larger computational cell) and the 3D dimensionality of the sample. The obtained results of this reconstruction are in good agreement with the nominal parameters. Finally, we conclude that the MBDDT computational scheme can be used in conjunction with the GIXRF experimental technique as a powerful tool in element-selective nanometrology for 2D and 3D periodic structures.

Footnotes

Currently affiliated to Physikalisch-Technische Bundesanstalt, Berlin, Germany.

Acknowledgements

This work is part of the research programme of the Industrial Focus Group XUV Optics, being part of the MESA+ Institute for Nanotechnology and the University of Twente (https://www.utwente.nl/xuv). It is supported by ASML, Carl Zeiss SMT AG and Malvern Panalytical, as well as the Province of Overijssel and the Netherlands Organization for Scientific Research (NWO).

Funding information

This project has received funding from the Electronic Component Systems for European Leadership Joint Undertaking under grant agreement No 826589 – MADEin4. This Joint Undertaking receives support from the European Union's Horizon 2020 research and innovation programme and Netherlands, France, Belgium, Germany, Czech Republic, Austria, Hungary, Israel.

References

First citationAltissimo, M. (2010). Biomicrofluidics, 4, 026503.  Web of Science CrossRef PubMed Google Scholar
First citationAndrle, A., Hönicke, P., Schneider, P., Kayser, Y., Hammerschmidt, M., Burger, S., Scholze, F., Beckhoff, B. & Soltwisch, V. (2019). Proc. SPIE, 11057, 110570M.  Google Scholar
First citationBeckhoff, B. (2008). J. Anal. At. Spectrom. 23, 845–853.  Web of Science CrossRef CAS Google Scholar
First citationBeckhoff, B., Gottwald, A., Klein, R., Krumrey, M., Müller, R., Richter, M., Scholze, F., Thornagel, R. & Ulm, G. (2009). Phys. Status Solidi B, 246, 1415–1434.  Web of Science CrossRef CAS Google Scholar
First citationBedzyk, M. J., Bommarito, G. M. & Schildkraut, J. S. (1989). Phys. Rev. Lett. 62, 1376–1379.  CrossRef PubMed CAS Web of Science Google Scholar
First citationBorn, M. & Wolf, E. (2013). Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. Elsevier.  Google Scholar
First citationBuitrago, E., Fallica, R., Fan, D., Kulmala, T. S., Vockenhuber, M. & Ekinci, Y. (2016). Microelectron. Eng. 155, 44–49.  Web of Science CrossRef CAS Google Scholar
First citationChateau, N. & Hugonin, J.-P. (1994). J. Opt. Soc. Am. A, 11, 1321–1331.  CrossRef Web of Science Google Scholar
First citationDialameh, M., Ferrarese Lupi, F., Hönicke, P., Kayser, Y., Beckhoff, B., Weimann, T., Fleischmann, C., Vandervorst, W., Dubček, P., Pivac, B., Perego, M., Seguini, G., De Leo, N. & Boarino, L. (2018). Phys. Status Solidi A, 215, 1700866.  Web of Science CrossRef Google Scholar
First citationFletcher, P. C., Mangalam, V. K. R., Martin, L. W. & King, W. P. (2013). J. Vac. Sci. Technol. B, 31, 021805.  Web of Science CrossRef Google Scholar
First citationGibaud, A. & Hazra, S. (2000). Curr. Sci. 78, 1467–1477.  CAS Google Scholar
First citationGoray, L., Jark, W. & Eichert, D. (2018). J. Synchrotron Rad. 25, 1683–1693.  Web of Science CrossRef IUCr Journals Google Scholar
First citationGray, R. M. (2006). Found. Trends Commun. Information Theory, 2, 155–239.  CrossRef Google Scholar
First citationHönicke, P., Beckhoff, B., Kolbe, M., Giubertoni, D., van den Berg, J. & Pepponi, G. (2010). Anal. Bioanal. Chem. 396, 2825–2832.  Web of Science PubMed Google Scholar
First citationJörg, Z. & Kazimirov, A. (2013). The X-ray Standing Wave Technique: Principles and Applications, Vol. 7 of Series on Synchrotron Radiation Techniques and Applications. Singapore: World Scientific.  Google Scholar
First citationKennedy, G. P., Buiu, O. & Taylor, S. (1999). J. Appl. Phys. 85, 3319–3326.  Web of Science CrossRef CAS Google Scholar
First citationKröger, I., Stadtmüller, B., Kleimann, C., Rajput, P. & Kumpf, C. (2011). Phys. Rev. B, 83, 195414.  Google Scholar
First citationKrumrey, M. & Ulm, G. (2001). Nucl. Instrum. Methods Phys. Res. A, 467–468, 1175–1178.  Web of Science CrossRef CAS Google Scholar
First citationLi, W., Zhu, J., Ma, X., Li, H., Wang, H., Sawhney, K. J. & Wang, Z. (2012). Rev. Sci. Instrum. 83, 053114.  Web of Science CrossRef PubMed Google Scholar
First citationLubeck, J., Beckhoff, B., Fliegauf, R., Holfelder, I., Hönicke, P., Müller, M., Pollakowski, B., Reinhardt, F. & Weser, J. (2013). Rev. Sci. Instrum. 84, 045106.  Web of Science CrossRef PubMed Google Scholar
First citationMalerba, M., Alabastri, A., Miele, E., Zilio, P., Patrini, M., Bajoni, D., Messina, G. C., Dipalo, M., Toma, A., Proietti Zaccaria, R. & De Angelis, F. (2015). Sci. Rep. 5, 16436.  Web of Science CrossRef PubMed Google Scholar
First citationMarkov, I. L. (2014). Nature, 512, 147–154.  Web of Science CrossRef CAS PubMed Google Scholar
First citationMikulík, P. & Baumbach, T. (1999). Phys. Rev. B, 59, 7632–7643.  Google Scholar
First citationNovikova, N. N., Zheludeva, S. I., Konovalov, O. V., Kovalchuk, M. V., Stepina, N. D., Myagkov, I. V., Godovsky, Y. K., Makarova, N. N., Tereschenko, E. Yu. & Yanusova, L. G. (2003). J. Appl. Cryst. 36, 727–731.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationPisarenco, M., Quintanilha, R., van Kraaij, M. G. M. M. & Coene, W. M. J. (2016). J. Opt. Soc. Am. A, 33, 610–617.  Web of Science CrossRef Google Scholar
First citationPomplun, J., Burger, S., Zschiedrich, L. & Schmidt, F. (2007). Phys. Status Solidi B, 244, 3419–3434.  Web of Science CrossRef CAS Google Scholar
First citationSenf, F., Flechsig, U., Eggenstein, F., Gudat, W., Klein, R., Rabus, H. & Ulm, G. (1998). J. Synchrotron Rad. 5, 780–782.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationSherman, J. (1955). Spectrochim. Acta, 7, 283–306.  CrossRef CAS Web of Science Google Scholar
First citationSoltwisch, V., Hönicke, P., Kayser, Y., Eilbracht, J., Probst, J., Scholze, F. & Beckhoff, B. (2018). Nanoscale, 10, 6177–6185.  Web of Science CrossRef CAS PubMed Google Scholar
First citationStepanov, S. A., Kondrashkina, E. A., Köhler, R., Novikov, D. V., Materlik, G. & Durbin, S. M. (1998). Phys. Rev. B, 57, 4829–4841.  Web of Science CrossRef CAS Google Scholar
First citationWeber, T., Käsebier, T., Szeghalmi, A., Knez, M., Kley, E. & Tünnermann, A. (2012). Microelectron. Eng. 98, 433–435.  Web of Science CrossRef CAS Google Scholar
First citationYakunin, S. N., Makhotkin, I. A., van de Kruijs, R. W. E., Chuev, M. A., Pashaev, E. M., Zoethout, E., Louis, E., Seregin, S. Y., Subbotin, I. A., Novikov, D. V., Bijkerk, F. & Kovalchuk, M. V. (2014). J. Appl. Phys. 115, 134303.  Web of Science CrossRef Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775
Follow J. Synchrotron Rad.
Sign up for e-alerts
Follow J. Synchrotron Rad. on Twitter
Follow us on facebook
Sign up for RSS feeds