Grazing-incidence small-angle X-ray scattering study of correlated lateral density fluctuations in W/Si multilayers

An inhomogeneity of material in W/Si multilayer structures was studied with grazing-incidence small-angle X-ray scattering. The experimental study revealed lateral density fluctuations in the Si spacer layers.


Introduction
Thin film periodic multilayer reflective coatings are used in spectroscopy and optical instruments in the photon energy range of a few tens of eV up to a few keV, i.e. from the extreme UV (XUV) to the soft X-ray range (Fewster, 1996;Louis et al., 2011). Depending on the application, various aspects of the structural quality of reflective coatings should be prioritized during the coating development process. As an example, often the reflectivity should be maximized over a certain bandwidth around a particular wavelength. However, for some applications like X-ray focusing optics the minimization of the offspecular diffuse scattering can be more important than maximization of the reflectivity. Primarily, diffuse scattering in multilayers is caused by interface roughness. In specific cases, diffuse scattering in multilayers may also be caused by a 3D distribution of defects in the volume of the structure (Pietsch et al., 2013). Analysis of the diffuse scattering provides information about the growth process including defects (Siffalovic et al., 2011). Such information can indicate directions for further optimization of the reflective properties of the multilayers.
In this study we focus on the characterization of the structural imperfections of a W/Si periodic multilayer with a 4.5 nm period. Such a multilayer is typically used for fluorescent spectroscopy analysis in the XUV range. A first multilayer structure characterization using scanning transmission electron microscopy (STEM) suggested that the Si spacer layers have lateral density fluctuations, different to what was expected from the growth model (Pelliccione & Lu, 2008). This observation was confirmed using grazing-incidence smallangle X-ray scattering (GISAXS). Numerical simulations of the GISAXS experimental data were performed to analyse the statistics of the distribution of the density fluctuations and the morphology of the interface roughness.
2. Sample preparation and characterization 2.1. Sample preparation W/Si multilayers were deposited on Si super-polished substrates using a magnetron sputtering system (Louis et al., 2011). A 6.35 mm-thick substrate was used in order to avoid a deformation (curvature) of the substrate caused by the tensile stress in the multilayer. A W target of 99.95% and Si target of 99.99% purity were used. Ar was used as a sputtering gas at a base pressure of 0.01 Pa. Magnetrons were set at a constant power mode equal to 86.7 and 213 W for the W and Si targets, respectively. The spatial uniformity of the coating was achieved by rotating the sample holder at 60 r min À1 during the deposition. The multilayer coating prepared contained 50 W/Si bilayers with a period thickness D = 4.5 nm and W layer thickness to period thickness ratio of À ¼ 0:23.

Preliminary sample characterization with XRR and HAADF-STEM
As a post-growth characterization of the manufactured sample, X-ray reflectivity (XRR) was used. The XRR curve was measured on an Empyrean laboratory diffractometer from Malvern Panalytical using characteristic Cu K 1 radiation from a long fine-focus line source. Monochromatization and primary collimation of the incident beam were achieved using a four-bounce asymmetrically cut germanium monochromator, which yields a beam divergence of 0.012 . Collimation at the detector side was achieved by an anti-scatter slit in combination with a tunable electronic receiving slit of a PIXcel3D detector. XRR data were analysed using the freeform approach (Zameshin et al., 2016). The resulting electrondensity depth profile was used to verify initial estimations of D and À, resulting in D = 4.46 nm and À = 0.2. The reconstructed parameters differ slightly from the design values because of the inter-diffusion between W and Si occurring during the deposition. These parameters are used further for the GISAXS numerical simulations presented in Section 4 for the qualitative investigation of the structure of the Si layer. XRR experimental data and analysis are presented in Appendix B2.
STEM measurements were performed in a Titan 80-300 equipped with a spherical aberration corrector (probe corrector), using an acceleration voltage of 300 kV. The microscope was equipped with an energy-dispersive X-ray Si(Li) spectrometer (EDAX), high-angle annular dark-field electron detector (Fischione) and Gatan image filter.
High-resolution bright-field transmission electron microscopy (TEM) measurements and an electron diffraction pattern (not shown here) indicated the presence of W nano-crystals in the W layers, with a crystallite size comparable with the W layer thickness, typical for nanometre-thickness metal layers grown by sputter deposition.
To study the microstructure of the Si layers, high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) measurements were performed. The HAADF-STEM images (see Fig. 1) revealed the presence of inhomogeneity inside the Si layers. Such inhomogeneity is unexpected, since to our knowledge amorphous Si is not known to develop any such density fluctuations during growth.
In Fig. 1 the Si layers appear darker than the W layers due to both atomic number contrast in HAADF-STEM and thickness effects as considered by Van den Broek et al. (2012). The inhomogeneity inside the Si layers appears to be quasiperiodic, and its distribution is statistically analysed in Appendix B1. The HAADF-STEM image, however, does not provide more detailed information on the density fluctuations including their positional correlations throughout the multilayer. To extract such information, a GISAXS study was performed, since GISAXS is highly sensitive to correlated structural imperfections [see the comprehensive review by Renaud et al. (2009)].

GISAXS experiment
GISAXS measurements were done on the bending-magnet beamline Langmuir of the synchrotron radiation source Siberia-2 at the Kurchatov Institute (Korchuganov et al., 2012). Monochromatization at the beamline is carried out by a thermally stabilized two-bounce Si monochromator with (111) reflection. Higher harmonics of the monochromated beam are suppressed with quartz and tungsten X-ray mirrors. The synchrotron beam was collimated with three sets of slits. The resulting beam size is 50 Â 300 mm and the corresponding average direct-beam intensity is approximately 3 Â 10 7 counts s À1 . The vertical beam divergence is 4 arcsec and horizontal beam divergence 20 arcsec.
Experimental data were measured with a Pilatus 100k 2D detector. Measurements were taken at the wavelength = 0.1 nm in 12 exposures of 15 min each, in order to avoid detector saturation. The angle of incidence was set to 0 = 0.4 , in between total external reflection and the first Bragg peak. The sum of 12 GISAXS measurements is shown in Fig. 2 on a logarithmic scale, where the colour scale of Fig. 2(a) is chosen to emphasize high-intensity scattering, whereas the same data are given in Fig. 2(b) with the colour scale adjusted to emphasize lower-intensity details.
The most pronounced feature in the GISAXS scattering is the resonant diffuse scattering sheets (Kaganer et al., 1996). The series of these sheets is marked by arrows '1' in Fig. 2(a). These sheets are caused by the correlated interface roughness (Holý & Baumbach, 1994). In the literature, resonant diffuse scattering sheets are also referred to as Bragg sheets (Siffalovic et al., 2009). The scattering geometry for the Bragg sheets is shown in Fig. 2(a 0 ). Here, k in , k sc are the wavevectors for the incident and scattered beam, respectively, and q ¼ k sc À k in is the reciprocal-space vector. Bragg sheets are located at q z ' 2m=D.
A second feature marked by arrows '2' in Fig. 2(a) is known in the literature as the Bragg singularity lines (Kaganer et al., 1995). These Bragg singularity lines are minima caused by the destructive interference of the scattered waves with each other [k sc and k 0 sc in Fig. 2(a 00 )], while Bragg sheets are due to the Bragg interference of incident and scattering waves. Bragg singularity lines are located at positions for which the exit angle of the scattered beam sc satisfies the Bragg condition 2D sinð sc Þ ¼ m for the multilayer structure.
In addition to the features discussed above one can observe two other features in Fig. 2(b). The first feature is the 'halo' in between Bragg sheets [marked with arrows '3' in Fig. 2 The second feature is the 'shadow' in between the third and the fourth Bragg sheets [marked with arrows '4' in Fig. 2 In Section 4 we will show that these 'halo' and 'shadow' features can be reproduced by simulations taking into account scattering on correlated density fluctuations of the material in the Si layers. Within this approach the position of the 'shadow' is due to the size of the density fluctuations and the 'halo' is due to the correlation in the positions of fluctuations.

GISAXS theoretical background
Calculations of diffuse X-ray scattering are performed using a perturbation theory (Sinha et al., 1988). A rigorous formulation of the second-order perturbation theory is given in Kaganer et al. (1996). There, the theory is formulated in terms of the reciprocity theorem of electrodynamics (Landau et al., 1984). In this theorem, the scattering length f of the diffuse scattering wave E diff ¼ f =r expðÀikrÞ has the form where the function ðrÞ ¼ ðrÞ À 0 represents a deviation of the actual dielectric susceptibility ðrÞ of a structure from the value 0 of an ideal structure. Equation (1) is written in a scalar approximation. This 0 value is used in the calculation of the wavefields in the structure: the field E in is induced by the incident beam, and the field E sc by the diffusely scattered wave. Equation (1) is commonly referred to as the distortedwave Born approximation (DWBA). We have considered here s-polarized radiation, as this was used in the measurements. Therefore, the fields E in and E sc are represented as scalar functions.
The intensity of the diffuse scattering is described with the scattering differential cross section: Here, averaging is applied to account for the random nature of imperfections ðrÞ of the structure. ðr) is considered to be a stochastic variable, describing a spatial distribution and structure of the imperfections. In this work we will consider the form of equation (2) explicitly written for interface roughness (in Section 3.1) and the 3D paracrystal of the density fluctuations (in Section 3.2). Finally, the wavefields E in;sc are considered as plane waves with phase terms expðAEik in;sc Á rÞ. Therefore, considering the integration in equation (1), it is convenient to represent the imperfection of the structure in the form of a Fourier transform: ðqÞ, where q ¼ k sc À k in . The way of defining the probability density function of ðqÞ defines a recipe for the simulation of the diffuse scattering on various types of imperfections. We now briefly review the theoretical models.

Scattering on correlated interface roughness
For the calculation of the diffuse scattering on interface roughness of the multilayer, the wavefield within the jth layer is considered to be E ðjÞ ðzÞ ¼ E ðjÞ t expðik z zÞ þ E ðjÞ r expðÀik z zÞ. The second summation in equation (3), where the indices l; m; n; p denote the transmitted (t) and reflected (r) waves, has 16 terms. Each term has a correlation function C jk ðqÞ ¼ h i ðqÞ Ã j ðqÞi that describes the interface roughness morphology (Pelliccione & Lu, 2008) of each j; k pair of interfaces.
In our numerical simulations, we employ Ming's model (Ming et al., 1993) for the calculation. For the characterization of roughness morphology, Ming's model involves the following parameters: root-mean-square (r.m.s.) roughness amplitude , lateral correlation length , Hurst parameter H which characterizes the jaggedness of a sample and the vertical correlation length L vert . For a detailed description of these parameters, see Ming et al. (1993) and Siffalovic et al. (2011).

Scattering on the density fluctuations
Here we consider lateral density fluctuations as an array of spheroids included in a homogeneous matrix of the Si spacer layer. In the STEM image (see Fig. 1), one can notice that the fluctuations are confined within the Si layers. We note that the vertical sizes are close to the value of the Si layer thickness and are not correlated with the positions along the vertical direction. One can also note that the density fluctuations appear as spheroids of comparable sizes. Additionally, we assume that the density fluctuations are isotropically arranged in the lateral plane inside the Si layers and we neglect correlations between their sizes and positions. Based on these considerations, the density fluctuations are best described using a 3D paracrystal model employing the mono-dispersion and decoupling approximations.
A comprehensive model for the simulation of diffuse scattering on a 3D paracrystal (Eads & Millane, 2001) is given in Buljan et al. (2012), where scattering from quantum dots is investigated. There, an expression for the differential cross section is derived using a decoupling approximation. In addition to the decoupling approximation, we simplify the expression for the differential cross section given in Buljan et al. (2012) assuming the mono-dispersion approximation: That equation is derived assuming Here F R ðrÞ is the shape function of the density fluctuation. It is equal to unity inside the density fluctuation located at R and equal to zero elsewhere. Taking the Fourier transform of equation (5) in equation (1) results in two functions:F FðqÞ and GðqÞ. The form factor functionF FðqÞ is a Fourier transform of a single density fluctuation shape function F R ðrÞ located at the origin R ¼ 0. This function is deterministic due to the mono-dispersion approximation. We considered a spheroid shape of the density fluctuation which has a form factor (Lazzari, 2002) that can be approximated byF where ¼ q k d L þ q z d z ; d L and d z are the lateral diameter and height of the spheroid, respectively, and V sp is the spheroid volume.
The correlation function GðqÞ is a sum of phase displacements related to positions of the density fluctuations: Explicit mathematical expressions for this model are given in Appendix A1. The characteristic parameters of this model are: lateral mean distance a L between each two neighbouring density fluctuations, dispersion L of the lateral mean distance, vertical mean distance a z and dispersion z defined analogously to the lateral parameters, lateral d L and vertical d z sizes of the density fluctuations. Thus, parameters a L , a z , L , z describe positions of the density fluctuations and d L , d z describe the size of the density fluctuations.

Analysis of correlated interface roughness
The diffuse scattering due to interface roughness was analysed by fitting model simulations, as described in Section 3.1, to three line extractions of the data: the first line extraction is taken in the plane of incidence (along the q z direction in q y = 0 nm À1 ), while the second and third line extractions are in the Bragg sheet planes (along q y in q z = 1.5 nm À1 and q z = 2.9 nm À1 ). Experimental data and best fits are shown in Fig. 3.
The best-fit parameters are: Hurst parameter H = 1, roughness r.m.s. amplitude = 0.2 nm, vertical correlation length L vert = 48.3D and lateral correlation length = 8 nm. From the good agreement between simulations of the interface roughness and the experimental data, we conclude that the diffuse scattering along the plane of incidence [ Fig. 3(a)] and in the Bragg sheets [Figs. 3(b) and 3(c)] is primarily caused by interface roughness, which is well described by Ming's model.
The interface roughness at = 0.2 nm is considered small. The value of H = 1 suggests that interfaces are smooth with no jaggedness. The vertical correlation length (L vert = 48.3D) is close to the full stack thickness, indicating a high degree of roughness replication from interface to interface.
The obtained parameters allowed us to simulate the scattering caused by the interfaces for the full range of experimental values of q y , q z . The results are shown in Fig. 4(a). As a visual aid, Fig. 4(b) shows the experimental data, together with a contour line based on the simulation of interface roughness scattering that is shown in Fig. 4(a). We attribute the scattering intensity inside this contour line to the correlated interface,

Analysis of correlated density fluctuations
The intensity of the scattering from density fluctuations is strongly dependent on the correlation function and the form factor that describe the distribution of the fluctuations. To build the initial model for fitting one can already assume the characteristic parameters of the density fluctuation distribution.
In the approximation of an ideally ordered distribution of density fluctuations (described in Appendix A1), the scattering will have a peak at q y ' 2=a L . In Fig. 5 one can see maxima at q y ' AE0.8 nm À1 . Thus, as an initial guess we assume that the mean lateral distance a L = 8 nm. Tails of the curves in Fig. 5 are primarily determined by the form factor (see Appendix A2). Thus, simulating the slopes of the line-extraction curves in Fig. 5 (q y range from 0.9 to 1.5 nm À1 ) as an initial guess, we assume the sizes of the density fluctuations d L = 4 nm, d z = 2 nm. Statistical parameters of the density fluctuations were estimated using a best fit to the GISAXS line extractions of experimental data. The confidence intervals of these parameters were estimated using a Hessian matrix calculated at the local minima of the best fit. Statistical parameters of the density fluctuations were also estimated using STEM data (see Appendix B1). The results are shown in Table 1.
Visible in Table 1 is the excellent agreement between STEM and GISAXS in estimates of the positional parameters (a L , a z , L and z ). It is noted that the mean lateral distance of the density fluctuations matches very well with the interface roughness lateral correlation a L ' = 8 nm. It hints that formation of the density fluctuations affects interface roughness morphology. The mean vertical distance a z matches very well with the period of the multilayer D = a z = 4.5 AE 0.2 nm and the dispersion of the mean vertical distance z is lower than the r.m.s. roughness amplitude estimated by fitting an interface roughness model, z < (see Table 1). Thus, we conclude that the density fluctuations are confined to the Si layer and do not penetrate into the W layer. This observation is consistent with the STEM image shown in Fig. 1.
Considering the absolute value of the scattered intensity, the size of the density fluctuations and correlated parameters, using equation (4)   Experimental data (markers) and model simulations based on interface roughness (solid line), for line extractions of data at (a) q y = 0 nm À1 , (b) q z = 1.5 nm À1 , (c) q z = 2.9 nm À1 . the density contrast: Á = 0.26 AE 0.05 g cm À3 , which is approximately 11% of bulk Si density in normal conditions. Using the GISAXS best-fit parameters, the intensity of the diffuse scattering was simulated for the full area of the 2D detector shown in Fig. 4(c). In Fig. 4(c) one notices the interesting effect that scattering on the density fluctuations not only causes an enhancement of the scattered intensity between Bragg sheets ('halo'), but also a 'shadow' of the scattering by destructive interference. This 'shadow' is defined by the shapes of the density fluctuations in the statistical ensemble. The position of the shadow in Fig. 4(c) is consistent with experimental data in Fig. 4(b). Numerical examples of how 'shadow' and 'halo' features change with the variation of the parameters of the paracrystal model are discussed in detail in Appendix A2. The model used for the numerical simulations is one of the simplest to describe the arrangement of the density fluctuations in the periodical structure. However, the good agreement of simulations with measurement data justifies the approximation used.
Comparing the absolute scattered intensities from interface roughness and density fluctuations, we conclude that diffuse scattering in the studied W/Si multilayer system is primarily due to the interface roughness, with only approximately 8% of the total scattered intensity being caused by the density fluctuations. The formation of the density fluctuations in the Si layers of a W/Si multilayer system is unexpected. Additionally, comparison of the mean lateral distance of the density fluctuations and lateral correlation length of the interface roughness hints that density fluctuations affect the interface roughness morphology. Thus, the analysis of these density fluctuations is of interest for understanding the multilayer growth model.
Our hypothesis is that the density fluctuations in the Si layers are formed due to the interaction with the high-energy back-scattered ions present during the magnetron sputtering deposition. Interacting with the sample surface during the Si layer growth, these high-energy ions distribute energy to the Si layers allowing the formation of the lower-density phase, i.e. density fluctuations. Preliminary analysis of various W/Si multilayers deposited with various doses of ion assistance show that higher ion currents result in stronger 'shadow' and 'halo' effects. For a more detailed investigation of the density fluctuations, the formalism of physical kinetics (Lifshitz et al., 1981) can be used. That theory allows one to analyse the dynamics of the statistical parameters of density fluctuations. In this article, we estimated the final statistical parameters of the density fluctuations, which can be used in further analysis in the formalism of physical kinetics as a subject of further research.

Conclusions
HAADF-STEM and GISAXS were used to study density fluctuations inside Si layers in periodic W/Si multilayers of nanoscale-thickness films. The fluctuations are ordered vertically with a dispersion of z ' 0.11 nm and mean distance a z ' 4.5 nm which is equal to the period of the multilayer sample. In the lateral direction, i.e. within the Si layer, these density fluctuations have a mean mutual distance of a L ' 8 nm, while the dispersion in the lateral direction is L ' 3.2 nm. The density fluctuations are strongly confined within the Si layers and have reduced density (Á ' 0.26 g cm À3 ). This study exemplifies the level of detail on growth phenomena that can be found using a combination of STEM and GISAXS analysis.

APPENDIX A Statistical characterization of density fluctuations A1. Correlation function of the paracrystal model
For the simulations in Section 4.2, the distribution of the density fluctuations is taken into account using a paracrystal model. Considering the HAADF-STEM data (Fig. 1) we assume that fluctuations have short-range ordering: the positions are correlated more strongly for neighbouring density fluctuations than for distant density fluctuations. Short-range ordering is taken into account using cumulative position errors (Buljan et al., 2012). Within that approach the correlation function for a 1D paracrystal lattice has the form where = =ð Ã À Þ and = ½1 À ðÞ N =ð1 À Þ, = expðÀiq Á aÞ is the phase displacement related to the paracrystal lattice vector a, = hexpðq Á dÞi is attributed to the Experimental data (markers) and model simulations based on interface roughness (red line) and density fluctuations (blue line), for line extractions of data at (a) q z = 2.2 nm À1 , (b) q z = 3.6 nm À1 . position dispersion, d is the cumulative position error vector and is an effective parameter related to the number of irradiated particles N and X-ray absorption: Assuming that the positions of the density fluctuations are normally distributed, can be calculated using a general formula (Payne & Clemens, 1993): where r is the vector whose components are position dispersion parameters along the x, y and z axes. In the correlation function [equation (8)] absorption of X-rays is taken into account. Note here that equation (8) takes into account 3D displacement of an object from its original positions defined as a 1D paracrystal lattice.
Let us now consider limiting cases. In the case of no absorption, ! N, the correlation function converges to Both functions, equation (8) and equation (11), have limits: lim jqj!0 G ! N 2 and lim jqj!1 G ! N. In equation (4) one can notice that dependence of the differential scattering cross section on q y is defined with the correlation function and the form factor, equation (6). For jqj ! 1 the correlation function is constant. Thus, in this case the scattering intensity is primarily due to the form factor, equation (6). In other words, we can assume that the shape of the scattering 'tails' in line extractions of the experimental data shown in Fig. 5 is solely defined by the form factor. This assumption was used as an initial guess in the fitting procedure in Section 4.2. Finally, in the limiting case of no dispersion jdj ! 0, we obtain G ¼ sin 2 ðq Á aN=2Þ sin 2 ðq Á a=2Þ : In this case the correlation function converges to the interference function. In other words, within the paracrystal model, scattering on an ideal distribution of fluctuations is similar to diffraction from an ideal crystal. For the numerical simulations in Section 4.2 we approximate the correlation function of the 3D distribution as where G i is the correlation function of the 1D paracrystal model, equation (8), with 3D lattice basis vectors chosen as a 1 ¼ ða L ; 0; 0Þ; a 2 ¼ ð0; a L ; 0Þ; a 3 ¼ ð0; 0; a z Þ; ð14Þ and the dispersion is parameterized with vectors: Strictly speaking, equation (13) does not take into account position cross-correlations between different lattice vectors, i.e. displacements of an object along a i and a j , i 6 ¼ j, are considered as statistically independent. To compensate for this in the simulations we performed azimuthal averaging of the correlation function: equation (13) is numerically integrated with respect to rotation of basis vectors around the z axis. In the following sections we discuss how variation of the statistical parameters around the best-fit solution affects the correlation function and the form factor.

A2. Variation of statistical parameters
The correlation function, equation (13), calculated for the best-fit parameters (Table 1) with various values of the lateral mean distance parameter a L is shown in Figs. 6(a), 6(b) and 6(c) and the correlation function with various values of the vertical mean distance a z is shown in Figs. 6(d), 6(e) and 6(f).
It has been shown in Appendix A1 that for the limiting case q ! 1 the correlation function converges to a constant value. Indeed, in Fig. 6 one can see that the correlation function is constant near the edges of the frame. Increasing and decreasing the value of the lateral mean distance a L [Figs. 6(a), 6(b) and 6(c)] results in a change of the 'halo' along the q y axis. The width of the 'halo' is proportional to $ 2=a L . A change of the vertical mean distance a z [Figs. 6(d), 6(e) and 6(f)] results in a shift of Bragg sheets and the 'halo' along the q z direction. It is important to note here that Bragg sheets in Figs. 6(d), 6(e) and 6(f) are not caused by interface roughness. Here, Bragg sheets are due to the periodicity of the density fluctuations in the z direction. In the experimental data in Fig. 4(a) these peaks are not resolvable due to the interface roughness contribution to the diffuse scattering [ Fig. 4(b)]. However, the vertical position of fluctuations can be estimated by the shape of the 'halo' feature.
In the HAADF-STEM image (Fig. 1) one can see that density fluctuations are confined within the Si spacer layer. The same conclusion can be drawn from best-fit parameters: mean vertical distance a z is equal to the period of the multilayer D ¼ a z ¼ 4.5 nm and dispersion of the mean vertical distance z is lower than the r.m.s. roughness amplitude estimated by fitting the interface roughness model, z < . Let us consider the case in which fluctuations start to penetrate the W layer. One can simulate that situation by increasing parameter z . In that case the 'halo' features are blurred significantly already at the dispersion values of z ¼ 0.5 nm [see Fig. 7(f)].
In the opposite case of low values of distance dispersion z = 0.11 nm and L = 0.8 nm [see Fig. 7(a)] the correlation function is similar to the diffraction on the crystal powder. That observation is consistent with equation (12) where the limiting case jrj ! 0 is considered.
The form factor function, equation (6), calculated for various size parameters is shown in Fig. 8. One can observe in Fig. 8 that the position of the 'shadow' from the best-fit model [see Fig. 4(c)] is not as broad as it is in the experimental data [see Fig. 4(a)]. This might be due to the mono-dispersion approximation. Finally, one can observe in Fig. 8 that the value of the form factor decreases along the q z direction. Since the scattering differential cross section distribution along q y is defined only with the correlation function and the form factor [see equation (4)] and the correlation function is constant near the edge of the frame in Fig. 6, the tails in Fig. 5 are defined by the form factor.

APPENDIX B Details of preliminary sample characterization B1. Statistical analysis of HAADF-STEM image
In this section we describe the procedure that has been employed for estimation of the statistical parameters in Table 1 from the STEM image. The boundaries of fluctuations were manually marked on the STEM image, four points for each dark object in Fig. 9. The position of each fluctuation has been calculated as a centre of mass of these four points. Lateral and vertical sizes have been calculated for each fluctuation as maxðx i Þ À minðx i Þ and maxðy i Þ À minðy i Þ, respectively, where     HAADF-STEM image of a W/Si multilayer. Each density fluctuation is manually marked with four dots for the statistical analysis.
To estimate confidence intervals of these parameters we used a t distribution. Position parameters a L , a z , L and z were considered to be normally distributed, while size parameters d L and d z were considered to be gamma-distributed, since variation of size parameters is large and negative values are nonphysical. Negative a L and a z parameters are nonphysical as well; however, since variation of these parameters is smaller than that of the size parameters, we therefore consider a normal distribution to be a valid assumption for these parameters.

B2. XRR analysis
XRR curve fitting has been used for a preliminary characterization of the sample. XRR was measured on a laboratory diffractometer using the characteristic Cu K radiation (incident-beam wavelength = 0.154 nm). The free-form approach was used for XRR fitting [see Zameshin et al. (2016) for a detailed description]. The XRR curve was measured up until the incidence angle max = 8 . That corresponds to the maximal vertical momentum transfer q ðmaxÞ z = 11.4 nm À1 . According to the sampling theorem, the minimal discretization step of the optical constant profile ( profile) is d min = =q ðmaxÞ z . Fitting by a model with a discretization step lower than d min does not provide relevant information about the structure of the sample. Therefore, the model of the W/Si bilayer period contained 16 sublayers with individual thickness of approximately 0.28 nm each. In the model all 50 periods were assumed to be identical. The measured XRR curve was fitted by a simulated curve using the derivativefree bound-constrained optimization algorithm BOBYGA (Powell, 2009). The 2 functional (Hughes & Hase, 2010) was used as a goodness-of-fit criterion. The measured XRR curve and the best-fit simulated curve are shown in Fig. 10(a). Residuals normalized on measurement uncertainty are shown in Fig. 10(b). The best-fit model corresponds to the goodnessof-fit criterion of 2 = 2.3. The profile of the best fit is shown in Fig. 10(c). Model uncertainties [blue error corridors in Fig.  10(c)] were calculated using the formalism of an inverse Hessian.
In the model, the surface of the sample was considered separately from the period to take oxidation into account. That allowed us to improve the goodness of fit in between Bragg peaks [see the inset in Fig. 10(a)]. As a reference, the vertical size of the density fluctuations d z is shown in Fig.  10(a). One can notice that the vertical size of the density fluctuation fits perfectly in the Si layer, which is consistent with GISAXS and TEM data.