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

Following the recent demonstration of the sensitivity of grazing-incidence X-ray fluorescence to the lateral structure of periodic nano-patterned devices, a computational scheme for the simulation of experimental data is presented. This can be used for the element-selective analysis of 3D atomic distributions in nano-patterned structures.


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;Buitrago et al., 2016). 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 threedimensional (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;Andrle et al., 2019). 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 ISSN 1600-5775 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 onedimensional (1D) depth distributions of fluorescent atoms have been well developed (Bedzyk et al., 1989;Jö rg & Kazimirov, 2013) and implemented for the study of epitaxial layers (Krö ger et al., 2011), multilayers (Yakunin et al., 2014), Langmuir-Blodgett films (Bedzyk et al., 1989;Novikova et al., 2003) and shallow ion implant profiles (Hö nicke et al., 2010), among others. However, if nanoscale devices, e.g. light-trapping structures in solar cells (Krö ger et al., 2011), field emitter arrays (Fletcher et al., 2013) and nanorods (Malerba et al., 2015), 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). 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;Dialameh et al., 2018), 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) used for the analysis in conventional XSW (Yakunin et al., 2014) 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), where the 2D structure of a lamellar Si 3 N 4 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). 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) 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;Hö nicke et al., 2010) 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) 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).

Theory
In Sections 2.1-2.3 we consider the theoretical background of the dynamical diffraction theory in the many-beam approximation (Mikulík & Baumbach, 1999) [in the literature also referred to as the rigorous coupled-wave analysis (Chateau & Hugonin, 1994)]. In Section 2.5 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.

Many-beam dynamical diffraction theory
The experimental geometry used by Dialameh et al. (2018) and Soltwisch et al. (2018) for the GIXRF measurements is shown in Fig. 1(a). 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, 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ðrÞ is the electric field, the sample structure is represented by the dielectric susceptibility function (r) and k 0 = 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), 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). In the dynamical diffraction theory, equation (1) is solved assuming that NF is represented as a Bloch wave, and the structure is represented as the Fourier series, where h is the Fourier component, Here, integration is taken over the unit-cell area for the corresponding reciprocal space vector, with n x, y the order of the diffraction index and D x, y the periods along x or y directions. The parallel component of the wavevector of the hth diffraction order k hk = k 0k þ h is translationally invariant along the z direction; i.e. k hk 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, This equation is derived assuming that diffraction scattering is an elastic process: k h = (1 + 0 )k 0 and assuming translational invariance of k hk mentioned above. Finally, substituting equations (2), (3) and (6) into (1), considering a property of the Fourier components, i.e. g expðih Á rÞ = gÀh , results in a system of inhomogeneous linear ordinary differential equations (ODEs) of second-order, 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 T n and R n . Thus, the hth solution of equation (7) has the form with linear combination coefficients E hn . Therefore, the distribution of the NF is defined with equations (2) and (8).
Thus, the problem of NF calculation is reduced to finding k n, z , E hn , T n and R n .

Characteristic equation
In this section we discuss the calculation of k n, z and E hn . Variable k n, z has a physical meaning as the vertical component of the wavevector [see equation (8)]. It defines the phase of the standing wave in the structured layer. One can assume that k n, z is defined with spherical dispersion k z = q z ; however, under that assumption equation (7) has no solutions. Therefore, values of k n, z deviate from spherical dispersion. To calculate the precise value of k n, z in the structured layer one can substitute equation (8) into (7). The result is represented as the eigenvalues-eigenvectors problem, where k 2 zn is an eigenvalue of matrix A and E n is an eigenvector composed of the coefficients E hn from equation (8): Matrix C is the Toeplitz circulant matrix, and X is the diagonal matrix with diagonal ð. . . À k 2 À1;k ; Àk 2 0;k ; Àk 2 1;k . . .Þ. Circulant matrices have a remarkable property: with increasing circulant matrix size, their eigenvalues asymptotically approach the exact values for an infinite matrix (Gray, 2006). Therefore, one can use a finite amount of Fourier components in equation (3) to approximate the exact solution of equation (7). Consider a set of 2N +1 Fourier components f ÀN ; ÀNþ1 ; . . . ; 0 ; . . . NÀ1 ; N g. These   Fourier components constitute a circulant matrix C of size C 2 C MÂM , where M = N + 1. Solving the characteristic equation (9) will give M eigenvalue-eigenvector pairs.

Boundary conditions
In this section we calculate the transmission T n and reflection R n amplitudes. Consider a sample as a stratified medium, consisting of layers. T n and R n are calculated in each layer using continuity conditions of the electric field and its first derivative. The continuity conditions (Born & Wolf, 2013) for the jth and ( j + 1)th pair of layers can be written in a matrix form, Here T and R are vectors composed of amplitudes T n and R n , Equation (12) 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 P is the refraction matrix. For a structured layer it has the form and for a homogeneous layer Here, the matrix E is composed of columns of eigenvectors and matrix k z is a diagonal matrix filled with k z, n . Refraction matrix P is a 2 Â 2 block matrix, thus P 2 C 2MÂ2M . Finally Q is the propagation matrix, where Q AE are diagonal matrices with corresponding diagonals where d j is the thickness of jth layer. Although these equations can be used to calculate T j and R i , solving equation (12) 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.

Numerical stability
The problem of numerical stability in the matrix formalism of dynamical diffraction theory was considered by Stepanov et al. (1998). 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) separately for each block matrix by using recurrent formula. Although the recurrent matrix equations of Stepanov et al. (1998) 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), which is relevant to the experimental data we will consider further.
The continuity conditions for such a three-layer structure (vacuum -structured layer -substrate) are represented by the system of linear equations One can rewrite that system as follows, Here matrix M has the form of a block matrix, where V ij is a matrix element of 2 Â 2 block matrix V (vg, gs) = (P (v, g) ) À1 P (g, s) . Note that this equation does not include Q þ whose elements are growing exponentially with respect to the thickness of the structured layer. Hence this matrix is numerically stable. Amplitudes T ðvÞ represent the incident beam, therefore Additionally, for a sufficiently thick substrate we can assume Taking into account these considerations, we derive equations for amplitudes in the structured layer, These amplitudes are calculated at the interface between the structured layer and the substrate (see Fig. 2). One can calculate amplitudes at the vacuum-structured layer interface using e T T ðgÞ e R R ðgÞ Finally, we need to rewrite equation (8), Here, both exponents decrease with respect to the depth, providing numerical stability.

X-ray fluorescence intensity
The fluorescence intensity Y can be calculated using the Sherman equation (Sherman, 1955), adapted for GIXRF (Hö nicke et al., 2010), where pðrÞ describes the density distribution of fluorescent atoms in the structure, and G() is the geometrical factor (Beckhoff, 2008;Li et al., 2012;Lubeck et al., 2013). The integral is taken over the area of the elementary cell. The exponential term expðÀzÞ in equation (27) 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) can be separated as follows [further, for brevity we do not explicitly write the multiplicative term G()], 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 k ); i.e. equation (28) 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) for each sublayer.
Equation (27) was rewritten in the form of equation (28), so it can be conveniently represented in a linear algebraic language. The fluorescence intensity can be expressed as the sum of matrix elements Y / P g;h F gh of the matrix, Here represents an element-wise (Hadamard) multiplication and H represents a Hermitian transpose. Elements of matrix U have the form and elements of matrix W have the form where UðqÞ Z 0 Àd expðiqzÞ expðÀzÞ dz: The U matrix takes into account the distribution of fluorescent atoms and the electric field distribution in the lateral direction and the U matrix takes into account photon absorption and the electric field distribution in the vertical direction. Equation (29) allows the integral in equation (27) to be calculated analytically, which is much more computationally efficient compared with the numerical integration.

Numerical simulations
3.1. 2D structure: Si 3 N 4 lamellar grating Here, we consider a 2D lamellar Si 3 N 4 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). The grating has a nominal period of D x = 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 planegrating monochromator (PGM) beamline (Senf et al., 1998) for undulator radiation at the PTB laboratory (Beckhoff et al., 2009) 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)]; = 0 corresponds to the conical orientation (Goray et al., 2018) 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 research papers efficiency and the geometrical factor (effective solid angle), were applied [see Soltwisch et al. (2018) for further details].
Best-fit simulations obtained by a sequential least-squares optimization algorithm and experimental GIXRF data are shown in Fig. 3 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)]. 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), this angle is not greater than = 4 . In terms of the model it means that the Fourier transform in equation (3) 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. h h i ¼ h expðÀh 2 2 =2Þ, with defined as half the projection of the sidewall on the x axis: d arctan ðÞ=2. The best-fit line width (defined as the half-height width) is D l = 39 nm and the best-fit sidewall tilt angle is = 5 .
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), such that the integration in equation (32) 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 d t = 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 d s = 0 nm. We note that this value d s 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 Si 3 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. A full set of 48 experimental GIXRF curves taken along different azimuthal angles (from 0 to 2 ) was interpolated on a (, ) grid  Comparison between the experimental GIXRF N-K map (a) of the Si 3 N 4 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 Si 3 N 4 grating structure, caused by interference between the reflected beam and the mth order of diffraction.
[see Fig. 4(a)]. The theoretical GIXRF map was calculated on the same (, ) grid using best-fit parameters from the data presented in Fig. 4.
One can note a distinctive feature on the GIXRF mapresonant lines, which are visible both on the experimental data in Fig. 4(a) and in the numerical simulations [see Fig. 4(b)]. As a visual aid to notice these lines one can refer to the sketch in Fig. 4(c). In Fig. 4(c) 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 k 2 x þ k 2 z = ðk x þ hÞ 2 . This formula geometrically corresponds to the Ewald sphere. For convenience we rewrite this equation in terms of the incidence and azimuthal angle, where = m/D x . The contour lines in Fig. 4(c) were calculated using this equation. Note that the resonant lines depend only on the lateral period of the structure D x and the wavelength [see equation (33)] -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.

3D structure: Cr nanocolumns
In this section we consider a periodic 3D nanocolumnar structure of Cr, manufactured using electron beam lithography (Altissimo, 2010) on top of a SiO 2 substrate. The structure of the sample is a regular square grid of box-shaped columns [see Fig. 1(c)] on a substrate, with 300 nm Â 300 nm lateral box dimensions and a D x = D y = 1 mm 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) in the PTB laboratory (Senf et al., 1998) of the BESSY II storage ring and reported by Dialameh et al. (2018). 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. The GIXRF experimental data and the best-fit obtained from dynamical diffraction theory simulations are shown in Fig. 5 for a selection of azimuthal angles. Best-fit model parameters are: lateral period of the structure D x = D y = 1 mm, 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, d t = 0 nm; however, the effective thickness of the sidewalls is d s = 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 SiO 2 = 2.4 g cm À3 . Considering the large lateral period D x = D y = 1 mm (significantly larger than that of the Si 3 N 4 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 i h . The experimental GIXRF curves in Fig. 5 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 |k y | ) |k x |. 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 k x and k y directions are identical due to the symmetry D x = D y of the periodic structure.
Thus, GIXRF curves of Cr nanocolumns can be effectively represented in a first approach as a lamellar Cr grating with 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. reduced density equal to the averaged density of the actual 3D structure. However, a direct comparison between 3D and 2D simulations (Fig. 6) reveals some differences. For higher incident angles above the critical angle of total external reflection, the 2D model (dashed blue line in Fig. 6) yields a monotonous angular dependence, while the experimental GIXRF curve clearly exhibits oscillatory behaviour in that angular range, with a maximum at ' 1.5 . In Fig. 6, 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). 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 |k y | 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 .

Discussion
In Table 1 we compare the structure parameters of the 2D lamellar Si 3 N 4 grating, as reconstructed using the MBDDT simulations described in Section 2.5, 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, FEM simulations have also been performed.
The FEM simulations were performed using the JCMwave software (Pomplun et al., 2007) for a box model based on the best-fit parameters in Table 1. 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 finiteelement 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). A direct comparison of the MBDDT and the FEM simulations is shown in Fig. 7. The GIXRF maps are symmetrical with respect to an axis at = 0 . In Fig. 7(a) 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). Here, the relative discrepancy is defined as where Y ðf;mÞ ij are the GIXRF intensities calculated in each ( i , j ) point using the FEM and MBDDT methods, respectively.
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 Comparison of effective 2D and genuine 3D simulations of GIXRF Cr-K curve for 3D Cr nanocolumns structure in conical geometry ( = 0 ).  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 In Table 2 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.
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). 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 Si 3 N 4 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) instead of being homogeneously distributed, as assumed in Section 3.1. The specific localization of the dopant atoms is depicted by green boxes and the resulting simulated GIXRF maps for the calculated fluores-cence 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), 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) techniques.

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 Si 3 N 4 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

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 0 ) to (c 0 ): corresponding GIXRF maps. 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.