Novel theoretical approach to the GISAXS issue: the Green function formalism using the q-Eigenwaves propagating through a twofold rough-surfaced medium

To describe the 1D and 2D patterns of the grazing-incidence small-angle X-ray scattering (GISAXS) from a rough fractal surface, the novel integral equations for the amplitudes of reflected and transmitted waves are derived. To be specific, the analytical expression for the 2D total intensity distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dR_{tot} \left( {\theta ,\phi ;\theta_{0} } \right)}}{d\Omega }$$\end{document}dRtotθ,ϕ;θ0dΩ is obtained. The latter represents by itself a superposition of terms related to the GISAXS specular \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dR_{spec} \left( {\theta ;\theta_{0} } \right)}}{d\Omega }$$\end{document}dRspecθ;θ0dΩ and diffuse \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dR_{dif} \left( {\theta ,\phi ;\theta_{0} } \right)}}{d\Omega }$$\end{document}dRdifθ,ϕ;θ0dΩ patterns, respectively. Hereafter, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document}θ is the scattering meridian angle, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document}ϕ is the scattering azimuth angle; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta_{0}$$\end{document}θ0 is the angle of incidence. By using the above analytical expressions, the 1D and 2D GISAXS patterns are numerically calculated. Some new experimental measurements of the specular reflectivity curves Rspec(\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta_{0}$$\end{document}θ0) related to the fused quartz and crystal Si(111) samples have been carried out. Based on the theoretical approach developed, a direct least-squared procedure in a χ2-fit fashion has been used to determine the corresponding values of the root-mean-square roughness σ from the specular GISAXS reflectivity data.

Nevertheless, even in the case of simple mesoscopic structure predominates, no one states that the GISAXS theory that is going out beyond both the distorted-wave Born approximation (DWBA) and the pure perturbation theory based on the parameter kσ θ 0 < 1 seems to be completed ( k = 2π/ , is the X-ray wavelength, θ 0 is the grazing-incidence angle).
There are a lot of the theoretical works concerning the GISAXS issue, amid which the works [10][11][12][13] are widely cited. In the work 11 , the static Debye-Waller factors of the reflected and transmitted GISAXS q 0 -waves coherently scattered have been derived. In the framework of the GISAXS perturbation theory for solving the integral stationary wave equation, authors of the work 12 have developed the Green function method proposed in 10 . In the work 13 , to solve the GISAXS issue, authors have developed the brilliant DWBA formalism widely used to treating the various GISAXS issues. In the paper 14 authors carefully have compared both the possibilities of the perturbation theory in the framework of the Green function formalism and the conventional DWBA. They have concluded that the conventional DWBA utilizes some simplified assumptions regarding an unperturbed wave field nearly a rough surface, which can be removed by using the Green function theory. Even without any calculations, one can see the similarity and difference of the above theoretical techniques. Their similarity and difference are the following. Both of them use the same set of conventional Fresnel solutions as the q-Eigenfunctions to find out the perturbation theory solutions of the GISAXS issue (cf. [14][15][16][17]. At the same time, the major difference each to other consists of the following. The conventional DWBA produces the perturbation theory solutions from the beginning, whereas the Green function formalism generates them at the stage of the rigorous asymptotic equations for reflected and transmitted wave amplitudes 14 . From that, it immediately follows that the range validity of the perturbation theory solutions built in the Green function framework cannot be less but only be more than the conventional DWBA solutions. In the works 15,16 the so-called self-consistent wave approach (SCWA) to solve the basic integral Green function wave equation has been suggested and analyzed. In particular, one has shown that the SCWA allows the GISAXS solutions to satisfy the optical theorem in the limit of large correlation lengths, kℓ → ∞ 16 . The next, to treat the GISAXS issue based on the Green function formalism, one has suggested the q-Eigenfunction wave field approach when the GISAXS issue solution in search represents by itself the direct 2D Fourier transform of the q-Eigenwaves propagating through the twofold homogeneous medium [15][16][17] .
Despite of a definite progress in the GISAXS theory developing, the above theoretical approaches [15][16][17] have been not successive by volume, especially, for large values of the r.m.s. roughness σ.
To date, it becomes clear that the Green function formalism needs to be advanced in a proper way aiming to make the GISAXS theory output to be more transparent from both the mathematical (math) and physical viewpoints. Due to the foregoing, the reformulation of the Green function formalism becomes to be of great interest.
A goal of the present study is to develop the GISAXS theory approach based on the Green function formalism to be modified. In other words, our aim is to highlight the modified Green function formalism differed from ones earlier proposed in [15][16][17] , and in fact, to reformulate the theory fundamentals in a proper way deriving the novel self-consistent (non-averaged) wave field equations to determine the reflected and transmitted GISAXS wave field.
As it will be shown underneath, reformulation of the Green function formalism allows to shorten the intermediate math calculations in comparison with similar ones in [15][16][17] , and makes the analytical formulae describing the specular (coherent) and non-specular (incoherent) X-ray scattering to be more transparent and understandable. Following to 17 , we will search the wave field E(x, z) propagating through the twofold rough-surfaced medium in the form where h(x) being the surface height at the 2D planar point coordinate x; the average <h(x)> is assumed to equal zero. δ 2 q − q 0 is the 2D delta function related to the incident X-ray plane wave Hereafter, the reflected and transmitted amplitudes, B(q) and C(q), are the inverse 2D q-Fourier transforms of Eigenwavefield (1), ( k z (q) = (k 2 − q 2 ) 1/2 ). Furthermore, to evaluate the specular (coherent) and non-specular (incoherent) GISAXS scattering, the Gaussian ensemble-statistics model of a randomly rough fractal surface [17][18][19][20] is assumed in terms of the r.m.s. roughness σ = h(x) 2 , the averaged roughness <h(x)> is assumed to be equal to zero and the two-point is the modified Bessel function of the second kind, x ≡ |x 1 − x 2 | ; ℓ is the correlation length, h is the fractal-surface model index (FSMI).
The present study is exposed as follows.
In "Theoretical fundamentals. The modified Green function formalism" section, formulary of the modified Green function formalism is described in details. In opposite to the conventional Green function formalism [12][13][14][15][16][17] , the artificial plane interface at coordinate z = T, |T|> > σ, is now introduced to define the conventional Fresnel solutions in terms of the q-Eigenfunction waves propagating in the direct and mirror-reversed scattering geometry. To derive the integral self-consistent (non-averaged) wave field equations for determining the reflected and transmitted wave amplitudes B(q) and C(q), the auxiliary parameter T is chosen as T < 0 and T > 0, respectively, and besides, it should be |T|> > σ in order to eliminate T-interfering effects due to high values of the random roughness field h(x).
In "Self-consistence Integral equations for probing the GISAXS issue" section, as a result of the modified Green function framework, the two separate self-consistent wave field equations have been derived for the wave field amplitudes B(q) and C(q) no linked each to other by using (z → ∓∞)-asymptotic of integral equation that governs the wave field (1) in search.
In "Zero-and first-order perturbation theory solutions of the GISAXS problem" section, in the scope of the averaged roughness approximation, analytical expressions for amplitude B(q) ≈ B 0 (q) + B 1 (q), which correspond both the specular, B 0 (q), and diffuse, B 1 (q), scattering channels from a randomly rough surface are found out. Concerning the specular amplitude B 0 (q) and non-specular amplitude B 1 (q), the corresponding analytical expressions for static specular and non-specular scattering factors are derived.
Based on the Gaussian ensemble-statistics, the analytical expressions for the 2D non-specular intensity distribution is the 2D intensity distribution (2D scan) in the angular range {θ, φ} concerning the fixed grazing-incidence angle θ 0 , and θ and φ are the polar and azimuth angles of the reflected GISAXS beam (see the GISAXS schematic in Fig. 1).
In "Numerical run-through for providing the GISAXS patterns: results and discussion" section, keeping in mind the experimental GISAXS data scans, the computer simulations of the scans R spe (θ 0 ) , dR dif d� (θ, φ) and the DSI dR dif dθ (θ) are presented for various parameters {kσ , kℓ, θ 0 /θ cr } and the FSMI h. Herein, in the case of the DSI dR dif dθ (θ; θ 0 ) the grazing incident angle θ 0 is assumed to be fixed and omitted for simplicity.
In "Concluding remarks" section, some results of the specular GISAXS measurements for the fused quartz and crystal Si(111) samples are presented. The corresponding values of the r.m.s. roughness σ are determined by applying standard least-square procedure in a χ 2 -fashion.

Theoretical fundamentals. The modified Green function formalism
In accordance with the basic idea of the Green function formalism, the integral wave equation can be cast in the form Figure 1. GISAXS layout. k 0 = q 0 + k z n is the incident wavevector; k R = q 0 − k z n and κ T = q 0 + κ z n are the wavevectors of the specularly reflected and transmitted waves. k R = q − k z (q)n and κ T = q + κ z (q)n are the wavevectors of the diffuse reflected and transmitted waves. n is the unit vector along the z-direction perpendicular to the flat surface z = 0. Detector D is the 2D CCD camera. The Green (a point-source) function G(r, r 1 ) is defined as for the twofold medium with the step-like abrupt electric susceptibility: Functions y 1 (z, q) and y 2 (z, q) are the two linearly independent Eigenfunctions of the differential wave equation over the coordinate z Noteworthy is the fact that in the Green function method the value of z-coordinate T may be arbitrary to some extent. In our case, we choose the value of |T| >> σ and sign of z-coordinate T being either T < 0 or T > 0 depending on a goal-problem-setting.
Accordingly, in our case, in the direct and mirror-reversed GISAXS geometry, solutions of Eq. (5) represent by themselves the conventional Fresnel functions for the step-like abrupt electric susceptibility χΘ(z − T) as follows respectively.
The reflection R 1 (q), R 2 (q), and transmission T 1 (q), T 2 (q) coefficients take the conventional Fresnel form where the z-components k z q and κ z q in Eqs. (6)-(7) are equal to they relate to the same vector q parallel to the plane z = 0, κ 2 = k 2 (1 + χ). The free term E 0 (r) in the right-hand side of integral wave Eq. (2) takes the evident form Hereafter, the following notations are introduced As it is seen from Eqs.
determines the behavior of the reflected and transmitted q-Eigenwaves in search.

Self-consistence integral equations for probing the GiSAXS issue
From Eq. (2) it immediately follows that in a limit of the z-coordinate tended to ∓∞ , the asymptotic GISAXS equations can be cast in the form 17 whereas the scattering amplitudes A R q, q 0 and A T q, q 0 are nothing else the inverse q-Fourier transforms of the corresponding integrals of the X-ray wave field in search (1) over z 1 -coordinate, namely 17 : www.nature.com/scientificreports/ To calculate the scattering amplitude A R q, q 0 , taking into account the specific z-coordinate behavior features of the X-ray wave field in search (1) and the y 1 z 1 , q -Eigenfunction (see the first-line Eq. (6)), the low limit of integral over z 1 -coordinate T it is convenient to be negative and value of |T|>> the r.m.s. roughness σ. To be specific, the last condition is needed to exclude T-interfering effects due to high-valued roughness field h(x) (remind that the T is nothing else the auxiliary parameter being called up to overcome math complexities of working in the framework of the standard Green function formalism 17 , where parameter T seems to be equal to zero. Then, combining the first-line Eqs. (10)- (11) and using the similar math-calculation scheme 17 , one obtains the integral equation to determine the reflected wave amplitude B(q), namely: It is seen that the Eq. (12) does not contain the auxiliary parameter T and provides the reflected wave amplitude B(q) in the form of integral self-consistent integral relationship, not linking with the transmitted wave amplitude C(q) in opposite to the standard Green function theory 17 , in which they form a system of the two integral equations linked each to other.
Similar to the above case, combining the second-line Eqs. (10), (11) and following the standard math scheme developed in 17 , one obtains the integral relationship for the transmitted wave amplitude C(q) can be derived in the form provided that the auxiliary parameters T needs to be chosen positive, T >> σ (evidently, Eq. (13) does not contain parameter T ).
In general, self-consistent (non-averaged) wave field Eqs. (12), (13) allow to developing the relatively simple procedure for solving the GISAXS issue of the q-Eigenwaves propagating through a twofold rough-surfaced medium. Notice that in the case of the step-like abrupt electric susceptibility χ(r) = χΘ(z), when the surface roughness field h(x) = 0, Eqs. (12), (13) definitely provide the conventional Fresnel solutions for the wave amplitudes B(q) and C(q). It is important that both the Eqs. (12) and (13) are free out of auxiliary parameter T and form the theoretical background of the GISAXS theory in terms of the q-Eigenwaves propagating through the twofold rough-surfaced medium.
Underneath, we will apply the self-consistent Eq. (12) to build up analytical solutions for describing the specular (coherent) and diffuse (incoherent) GISAXS from a randomly rough-surfaced medium.

Zero-and first-order perturbation theory solutions of the GISAXS problem
To build up the theoretical foundation of the specular and non-specular GISAXS phenomenon from a randomly rough surface, one applies the averaged roughness field approach. Integral Eq. (12) can be identically rewritten as follows Here, the averaged exponential factors W − (q, q 1 ), W + (q, q 1 ) are defined as.
where the symbol . . . denotes the general ensemble-statistics average procedure.
By using the Gaussian ensemble-statistics, the averages (15) can be analytically calculated, which are nothing else the exponential GISAXS factors (11) (15) W − q, q 1 = e i(κ z (q)−kz(q1))h(x) , W + q, q 1 = e i(κ z (q)+kz(q1))h(x) , Scientific RepoRtS | (2020) 10:11547 | https://doi.org/10.1038/s41598-020-68326-2 www.nature.com/scientificreports/ Confining ourselves by terms of the zero-and first-order perturbation theory for deriving the reflected wave amplitude B(q), a standard procedure yields the following expressions for the reflected wave amplitude B(q) Notice that the corresponding equation for specular wave amplitude (16) obtained is identical to the expression first derived in the work 11 . Accordingly, one can utilize analytical Eqs. (16) and (17) for the zero-and firstorder perturbation theory amplitudes B 0 (q) and B 1 (q) to determine the 1 D specular and 2D diffuse GISAXS scans.
As a result, expression for the reflected 2D intensity distribution can be written in the form 13,17 and represents by itself the superposition of the two terms, namely: (d� = cos θdθdφ is the elementary solid angle for the reflected beam, S 2 is an area of a rough surface impinged on by the incident X-ray beam). Directly following the math procedure scheme 13,17 , standard calculations yield the analytical equations, which describe the 1D specular intensity distributions and the diffuse 2D intensity distribution remind that B 1 q = 0.
Then, specular reflectivity R spe (θ 0 ) and diffuse scattering intensity (DSI) dR dif (θ)/dθ are equal to 17 The functions PSD 2D q − q 0 and PSD 1D q − q 0 and expressions for specular factor f R q 0 and nonspecular scattering factor R q, q 0 , are derived as described in Sections A and B of Supplemental Material (cf. 13,[15][16][17][18][19][20][21]. Thus, the resulting analytical Eqs. (20)-(23) provide the GISAXS issue solution in a framework of the modified Green function formalism and in terms of the assumed Gaussian ensemble-statistics model for a randomly oriented fractal surface [17][18][19][20] . Noteworthy is the fact that according to our calculations, the specular scattering factor f R (q 0 ), Eq. (20), is identical to Nevot-Croce's one 11 . Staying in the scope of the adopted surface model and the first-order perturbation theory used to solving self-consistent (non-averaged) wave field Eq. (12), one can state that expressions (20)-(23) provide in the proper way analytical description of the GISAXS issue and allow to evaluate the 1D specular R spec (θ 0 ) and 2D non-specular dR dif (θ ,φ) d� scans, and the 1D DSI dR dif (θ ) dθ scan as well.
(15a) W − q, q 1 = exp −0.5 k z q − k z q 1 2 s 2 , W + q, q 1 = exp −0.5 k z q + k z q 1 2 s 2 . (16) Scientific RepoRtS | (2020) 10:11547 | https://doi.org/10.1038/s41598-020-68326-2 www.nature.com/scientificreports/ Before proceeding further, it remains to see how well the developed approach works, its validity range, and prospects of solution of the GIAXS issue taking into account the high-order scattering effects. Now, one thing is clear that consideration of the GISAXS issue in the frame of high-order perturbation theory based on Kubo's cumulant average diagram technique (see, e.g., 22 ). By using such the technique, the self-consistent (non-averaged) wave field Eq. (12) can be converted to the Dyson-type and Bethe-Salpeter-type integral equations to describing the specular and non-specular GISAXS scattering, respectively. In this respect, by using the Kubo technique 22 , endeavor has been undertaken in the work 23 , in which the Dyson-type equation has been derived and analyzed in application to the specular GISAXS scattering. As was shown in 23 , specular scattering factor f R (q 0 ) becomes dependent on both the r.m.s. roughness σ and roughness correlation length ℓ. Clearly, application of the Kubo technique could enable to improving the present GISAXS theory approach clarifying at least, its validity range.
In general, real surface does not satisfy to the Gaussian ensemble-statistics and K 2 -correlation fractal surface model adopted in the present study. As was pointed out in [13][14][15][16][17] , such the model has been applied to clarify the GISAXS peculiar pattern. Here, we can only state that the modified Green function formalism does work and could be applied in a general case of the surface model provided that the ad hoc model of a rough surface is in a matter. For instance, instead of the play-in model, one can utilize the surface model, in which χ(r) = 0 for z < 0 whereas for z > 0 and χ(r) = < χ(z) > + δχ(r) at each point r = ( x, z) of the statistically distorted surface layer, where δχ(r) is the fluctuating part of χ(r), < δχ(r) > = 0.

numerical run-through for providing the GiSAXS patterns: results and discussion
In this section, based on formulae (20)-(23), the numerically simulated results for the 1D specular and 1D-2D non-specular GISAXS patterns from a randomly rough fractal surface based are presented. In particular, one presents the results of proper calculations for different dimensionless values of fractal surface parameters as kσ and kℓ and the FSMI h as well.
The numerically simulated 1D scans R spec (θ 0 ) versus the grazing-incidence angle θ 0 in the angular range 0 ≤ θ 0 ≤ 10θ cr are shown in Fig. 2 in the cases of the fused quartz and Si(111) samples. The corresponding parameters kσ are equal to 40.7999 (σ = 1 nm for fused quartz) and 26.9279 (σ = 0.66 nm for Si(111)). Simultaneously, in Fig. 2 the curves of Nevot-Croce's exponential factors (NCEFs) are depicted. It is evidently seen that the curve NCEF for the fused quartz (blue curve) goes below to the curve for Si(111) (black curve) due to the about two-times difference in their r.m.s. σ values.
Some examples of the 2D non-specular intensity distribution numerically simulated according to analytical formulae (21) and (23) are presented in Figs. 3a, 4a for the fused quartz and in Figs. 3b, 4b for Si(111), respectively. As it follows from Eqs. (21), (23), they contain the non-specular scattering factor R q, q 0 , the explicit expression of which is done in Section B of Supplemental Material (cf. [15][16][17][18]. It represents by itself the superposition of the exponential factors, which contain the r.m.s. roughness σ to be squared. Dependence of the 1D d� -scans on the correlation length ℓ and the FSMI parameter h are due to the PSD 1D (|q-q 0 |)-and PSD 2D (|q-q 0 |)-functions, the corresponding expressions of which are given in Section A of Supplemental Material (cf. 17 ).
Hereafter, having aim to demonstrate how it works, retrieval of the r.m.s. parameter σ from the 1D experimental R spec (q 0 ) data for the fused quartz and Si(111) samples that are displayed in Fig. 5.
Both the sample surfaces have been finished by mechanical polishing with iron oxide powders, and then, have been cleaned with ethanol just before the experimental measurements. The extracted NCF-functions versus the variable 4(θ 0 / θ cr ) 2 are displayed for interval of the grazing-incidence angles 5 θ cr ≤ θ 0 ≤ 12.5 θ cr in the insert in Fig. 5.
Scientific RepoRtS | (2020) 10:11547 | https://doi.org/10.1038/s41598-020-68326-2 www.nature.com/scientificreports/ procedure of random roughness field h(x), the integral Eq. (12) for the reflected wave amplitude B(q) has been exploited to search the zero-and first-order perturbation theory solutions of the GISAXS issue. These solutions are obtained in the analytical form of expressions (20)-(23) for determining the 1D specular (coherent) and 2D non-specular (incoherent) intensity distributions in the case of the GISAXS by a twofold rough-surfaced medium.
To develop the current GISAXS approach, we have advanced the Green function formalism and, besides, some key assumptions have been assumed. Namely, the Gaussian ensemble-statistics and Mandelbrot's fractal surface model have been exploited. Appropriately, the specular f R (q 0 ) and non-specular R q, q 0 scattering have been derived. The specular scattering factor f R (q 0 ) identically coincides with the corresponding expression earlier obtained in 11 , whereas in the case of the non-specular GISAXS, new expressions (20)-(23) for describing the 1D and 2D GISAXS intensity distributions have been derived in terms of the PSD 2D (|q − q 0 |)-and PSD 1D (|q − q 0 |)functions, respectively. At last, based on formulae (24), (25), the inverse issue of retrieving parameters {σ, ℓ, h} of Mandelbrot's rough surface model from the experimental GISAXS data has been shortly discussed. Accordingly, applying a direct fit-procedure in a χ 2 -sense 25 and using the theoretical representation of function NCF (24) in the case of the experimental specular GISAXS data collected for the fused quartz and Si(111) samples, retrieving the r.m.s. roughness parameter σ has been successfully realized (see Fig. 5).