## research papers

## Incorporation of interfacial roughness into recursion matrix formalism of dynamical X-ray diffraction in multilayers and superlattices^{1}

^{a}Atomicus OOO, Minsk, Belarus, and ^{b}Atomicus GmbH, Karsruhe, Germany^{*}Correspondence e-mail: igor.lobach@atomicus.by

Diffraction in multilayers in the presence of interfacial roughness is studied theoretically, the roughness being considered as a

Exact (within the framework of the two-beam theory) differential equations for field amplitudes in a crystalline structure with varying properties along its surface normal are obtained. An iterative scheme for approximate solution of the equations is developed. The presented approach to interfacial roughness is incorporated into the recursion matrix formalism in a way that obviates possible numerical problems. Fitting of the experimental rocking curve is performed in order to test the possibility of reconstructing the roughness value from a diffraction scan. The developed algorithm works substantially faster than the traditional approach to dealing with a (dividing it into a finite number of thin lamellae). Calculations by the proposed approach are only two to three times longer than calculations for corresponding structures with ideally sharp interfaces.Keywords: interfacial roughness; transition layers; multilayers; dynamical diffraction; recursion matrix formalism.

### 1. Introduction

The phenomenon of interfacial roughness in multilayers has been shown to have an influence on the functioning of electronic devices (*e.g.* see Zolotoyabko, 1998; Bolognesi *et al.*, 1992; Ming *et al.*, 1993; People & People, 1986; People *et al.*, 1984; Glaser *et al.*, 1990). Therefore it is important to have a means to estimate this nonideality of interfaces. In the case of grazing incidence in X-ray reflection (XRR), when most of the reflected power is produced by surface layers, interfacial roughness considerably affects measured XRR scans and has been extensively studied by this technique (*e.g.* see Feranchuk, Feranchuk *et al.*, 2003; de Boer, 1994, 1995, 1996; Feranchuk *et al.*, 2007; Sinha *et al.*, 1988; Feranchuk, Minkevich & Ulyanenkov, 2003; Benediktovitch *et al.*, 2014; Pietsch *et al.*, 2013).

In X-ray diffraction (XRD), just like in XRR, grazing-incidence geometry is sensitive to interfacial roughness (see Stepanov & Köhler, 1994). One may expect that in conventional geometries, for example – scans, the roughness effect will be subtle. This is true in many cases; however, sometimes scans are sensitive to interfacial roughness even in conventional geometries, for example for superlattices consisting of thin layers (Fullerton *et al.*, 1992). This is why we consider interfacial roughness in conventional XRD geometry in this paper.

It is possible to treat interfacial roughness as a ) by dividing it into many thin lamellae. However, this approach turns out to be very time consuming. Therefore, we decided to speed up the process by performing analytical calculations for transition layers and incorporation of interfacial roughness into the recursion matrix formalism (RMF) (see Stepanov *et al.*, 1998). We will show that it is possible to introduce only one additional matrix for each layer, describing the roughness of this layer. This boosts the speed of calculation, which is important for fitting of experimental curves (Ulyanenkov & Sobolewski, 2005).

This paper is organized in the following way. In §2, basic notations for the two-wave theory (DDT) in multilayered structures are introduced, which allow the to be treated as many thin lamellae (traditional numerical approach; see Stepanov & Köhler, 1994). In §3, several types of transition layers are analyzed and one of them is selected for further consideration. In §4, we obtain differential equations for field amplitudes in layers by assuming that the thickness of each lamella tends to zero. A natural ansatz for solution of the equations is proposed and a useful approximate solution is found. In the framework of this approach, it is shown that the effect of interfacial roughness can be described by an additional matrix in the RMF approach. §5 provides numerical examples and illustrates the fitting of one experimental rocking curve.

### 2. Summary of two-beam DDT in multilayered samples

The structure considered in this section is presented in Fig. 1. In order to describe a by this structure one has to choose a change of Fourier components of susceptibility , , from layer to layer, such that the phase of oscillation of electron density does not encounter any jumps in the interface between neighboring layers. One can check by direct substitution that the following expression for the susceptibility in the *j*th layer is in agreement with the last requirement (Stepanov *et al.*, 1998):

with (layers are numbered by ). *L*^{( j )} is the thickness of the *j*th layer and *g*_{z}^{(m)} represents the *z* component of the vector in the *m*th layer . For see Fig. 1. In equation (2) , and are supposed to have the same phase as , and , respectively. Actually, these phases change slightly if one describes the transition between two real crystals owing to nonzero absorption, but this error is insignificant.

The problem of σ polarization only, because the results for π polarization can be readily obtained by multiplying and by the polarization factor , where is the The solution for the Fourier transform (over time) of the electric field in the *j*th layer is (see *e.g.* Stepanov & Köhler, 1994)

where *w* is the angular frequency in the Fourier transform and *i* = 1,2. represents the wavevector of the direct wave corresponding to the solution of the dispersion equation in the *j*th layer (in general, there are four solutions to the dispersion equation, but two of them correspond to negligible specularly reflected waves). *v*^{( j )}_{i} is the ratio of amplitudes of diffracted and direct waves in the *j*th layer. It is given by

where . represents the wavevector of the diffracted wave and is the wavevector of the incident wave. Actually, the component of the wavevector of the direct wave parallel to the crystal surface is fixed during solution of the dispersion equation. Therefore, we will call the normal component of the vector a solution of the dispersion equation. Hereinafter we denote the attenuating root by *k*_{1z} [] and the amplifying one by *k*_{2z} [].

The boundary conditions in the interface between the *j*th and ( *j*+1)th layers can be written as continuity of the electric field (see Benediktovitch *et al.*, 2014) separately for the waves with components parallel to the surface equal to and (:

Now let us switch to matrix notation, similar to that of Stepanov *et al.* (1998). We introduce a column vector of field amplitudes and two auxiliary matrices:

Making use of these matrices, equations (4) take the form

Even though now the equations look rather explicit, let us try to simplify them further by introducing new unknown amplitudes as a linear combination of current ones:

where *T* stands for `transmitted' and *D* stands for `diffracted'. *T*_{1} = 1 is the amplitude of the incident wave and is the amplitude of the wave diffracted from the whole sample. If we denote the vacuum below the sample as the (*n*+1)th layer, then *D*_{n+1} = 0, which means that no wave is incident on the lower surface of the sample in the direction opposite to the *z* axis. For the introduced amplitudes the following relations are true:

If we rewrite equation (9) in more detail (),

it becomes apparent that it is a linear system for unknown amplitudes of diffracted (*D*) and transmitted (*T*) waves. *D* can be expressed as

or, if we divide it by the increasing exponents exp( *i**k*^{(n)}_{2z}*L*^{(n)}),

In this way, the final expression for *D*,

does not contain any increasing exponents. Thus, consideration of thick layers will cause no numerical problems. Note that the approach presented here of precluding numerical errors is different from that reported by Stepanov *et al.* (1998). Clearly, one can consider any set of layers using the presented approach, not only transition layers. If one considers the more realistic problem where the (*n*+1)th vacuum layer is missing and the *n*th layer is the substrate (), then

### 3. Choice of model

In this section we consider several types of transition profiles and their influence on the diffracted intensity distribution. We investigate the following three types of transition profiles:

Type a

Type b

Type c

with

Additional indices 1 and 2 describe the properties of actual upper and lower layers, and , and are smooth functions providing the behavior shown in the lower part of Fig. 2:

where σ is the quantitative characteristic of roughness and *z*_{0} is the *z* coordinate of the interface between two actual layers (*z* of the upper surface of the sample is zero).

Type a (the Epstein profile; Epstein, 1930) causes an increase in diffracted intensity between the peaks (see Fig. 2*a*). The reason is that, if *g*_{z} changes smoothly, then at every angle of incidence (within the region between two peaks of different actual layers) there are sublayers in the for which the Bragg condition is fulfilled and they provide higher diffracted amplitude. In our opinion, the effect of interfacial roughness is likely to reduce the occurrence of diffraction in the lamellae close to the interface, rather than to shift it to other values of the incidence angle.

In cases of type b, the *z* component of changes stepwise and, as is seen in Fig. 2(*b*), the increase between the peaks disappears. However, near the left peak, one can see an enhancement, which is connected to the specific choice of dependence of and on *z*.

Type c (see Fig. 2*c*) can be understood as follows: The abrupt change of the *g*_{z} profile corresponds to an interface at which the average coherent lattice periodicity is changed. The drop in the susceptibilities and corresponds to a reduction of crystallographic order due to the static Debye–Waller factor conditioned by defects located at the interface. Such a situation can occur, for example, for the system described by Satapathy *et al.* (2005) where the nonuniform strain decays exponentially when moving away from the boundary. This type will be used further.

### 4. Differential equations for field amplitudes in transition layers: approximate solution

In the previous section we dealt with the problem of the *z*)' from now on], and the boundary conditions at all interfaces will be replaced by differential equations. Let us obtain them. Consider equation (8) in the case of infinitesimal sublayers. Making the substitutions , , in equation (8) results in

Further, it is necessary to take into account that under the assumption of infinitesimal sublayer thickness

where indices *j* are omitted and the dependence on *z* is assumed. Also, . By equation (8) with we obtain

Equation (19) is the differential equation for field amplitudes in a crystalline sample whose properties vary along the *z* direction. Its domain of applicability coincides with that of two-beam DDT with two diffraction roots [diffraction in each infinitesimal sublayer of the sample should be properly described by DDT in order that equation (19) be correct].

Our present goal is to find how roughness influences the transition matrix of one layer in a structure. In the case of nonzero roughness it is unclear what an interface is. We introduce it formally as the point where *g*_{z} jumps from one value to the other (as in type c). As we consider one layer, we choose new coordinates, where *z* = 0 in the upper interface. The thickness of the layer is *L*, the susceptibilities not perturbed by roughness are , and , and the susceptibilities perturbed by roughness are , and . We will also introduce a function : , . The vector is and is constant. Note that by the last definitions we can consider roughness in both the upper and the lower surfaces of the layer (corresponding modifications to type c will be made below). Further, we make ordinary approximations of two-beam DDT (search for roots in the vicinity of *k*_{0z}) and obtain an expression for the matrix . Then, after some algebra, equation (19) takes the form

*k*_{0z} is the *z* component of , which is the wavevector of the incident wave. Similarly, *k*_{g0z} is the *z* component of . Equations (20) are similar to Takagi–Taupin equations (Takagi, 1962, 1969; Taupin, 1964). However, the validity of the obtained equations is not restricted to a smooth variation of susceptibilities. In the case of constant susceptibilities, equation (20) can be solved by Euler's method, searching for a solution in the form of exponents, and the solution is as follows:

which is precisely what one would expect in the case of an ideal crystal. *E*_{1}^{(0)} and *E*_{2}^{(0)} are constants of integration. The meanings of *v*_{1}, *v*_{2}, *k*_{1z} and *k*_{2z} are the same as in equations (3) and (4) and they are expressed through , , , and .

We will seek a solution of equation (20) in the form of equations (21) with varying . Substitution of the ansatz into equation (20) gives the following:

with

One can see that if then the solution is . Equation (22) has an exact analytical solution in the case that , and are linear functions of , where σ is the roughness parameter. There is a plausible model consistent with the last requirement, but the solution is cumbersome and has numerical problems. It is expressed through the Laguerre polynomial and hypergeometric function (not shown here). However, it appears that an approximate solution to equation (22) is very accurate at realistic values of roughness. The relative error with respect to the exact solution is comparable to the accuracy of two-beam DDT. Namely, one can find a first-order approximation by integration of equation (22), assuming is constant on the right-hand side:

In further consideration we will need the value of :

We call the roughness matrix. Let us now consider equation (6) as the boundary conditions at the interface between two actual layers. One has the field at the lower interface of the *j*th layer on the left-hand side. In the case of nonzero roughness this field can be expressed through the above-mentioned ansatz using equation (25). Then equation (6) takes the form

where is the roughness matrix of the *j*th layer. Further, it is convenient to denote and obtain the diffracted amplitude from a multilayer structure of actual layers by equation (13), having made the substitution

in equations (12). Thus, the formal solution to the problem with nonzero roughness is found.

Below we will find an expression for and in the case of a

of type c. Let us modify it to account for roughness in both interfaces of the layer:Type c*

where upper indices U and L imply the layer above or below the one under consideration (see Fig. 3).

In this case are continuous functions and is almost continuous. Actually, a small discontinuity has a subtle influence as only the integrals of these functions will be used. Performing the integration in equation (25) in the case of type c* one obtains

where `' means the Hadamard product, that is (no summation over *i* and *j*).

with . However, in the final equation for the amplitude diffracted from the sample one needs an expression for :

This is expressed *via* the following quantities ():

It appears that and diverge with . However, the final expression for is not divergent since it involves the convergent quantities and in equations (31). Thus, to avoid numerical problems one should use of equation (30) using expressions from equations (31). We can check the limiting case . It is seen from equations (31) that this leads to as well. for the substrate layer is readily obtained by considering the limit in equations (31) and (30). Lastly, some expressions in equations (31) are convergent only if . However, this condition is not fulfilled only at very large values of the roughness parameter (hundreds of nanometres). The derivation of the roughness matrix for amorphous layers is analogous and also much simpler than for crystalline layers. Therefore, we do not present it here.

### 5. Numerical examples and fitting of experimental rocking curve

Let us prove that the suggested approximate solution is accurate. A comparison of exact and approximate scans for a simple test sample is given in Fig. 4. Hereafter, by exact solution we mean the one obtained by division of the sample into many thin lamellae of finite thickness. The values of the roughness parameter are very large in this example (up to the thickness of the upper layer). They are chosen in this way in order to find the limits of applicability of the approximate solution. In Fig. 4 one can see that the dashed lines start to deviate slightly from the solid lines at values of the roughness parameter of about 10 nm. For realistic values (up to several ångströms) the exact and approximate lines are indistinguishable, which proves the reliability of the approximate solution.

Fig. 5 shows calculated scans of a based on the type c model. Note that some peaks are enlarged in the case of nonzero roughness, whereas others are reduced. Actually, this is in accordance with diffraction in each period of the That is, if the diffracted amplitude from one period is decreased owing to roughness at a certain incidence angle, then the amplitude of the wave diffracted from the is reduced too and *vice versa*.

We also performed a test fitting of an experimental curve by the proposed method of calculation of rocking curves taking into account roughness (see Fig. 6). In addition, the fitting was done by the conventional nominal model with ideally sharp interfaces. In the case of the model with roughness, the variable parameters were the vertical scale, the layer's thickness, and the roughness parameter at the interface between the layer and the substrate. In the case of the nominal model, the only variable parameter was the vertical scale; the layer's thickness was fixed and equal to 5 nm (this follows from the period of thickness fringes). The best fit was determined by finding the curve with the smallest squared deviation (in logarithmic scale) from the experimental rocking curve. The experimental curve was taken from the work of Ulyanenkova *et al.* (2013).

The area near the substrate peak where the largest discrepancy is observed is dominated by the shape of the substrate peak. Since it is much higher in intensity than the layer, the deviation of its shape from DDT predictions affects the curve significantly. The deviations can be due to instrumental effects or due to the influence of defects [see, for example, the paper by Shreeman & Matyi (2010), where the fit of the region near the substrate for a similar structure was obtained in the frame of statistical theory].

In our case we are within the DDT formalism and the focus was on the intensity from the layer. Hence, the left part of the profile where the thickness fringes are visible is of interest. Here we can see that the interfacial roughness model follows the curve closer than the sharp interface model.

### 6. Conclusions

In this paper we obtained a general result, the application of which is not limited to interfacial roughness. We presented the exact (in the framework of DDT) equations for wave amplitudes in a crystal whose properties vary along its surface normal [see equation (20)]. These equations could be potentially applied to many other problems apart from interfacial roughness. For instance, one could use these equations in crystals with strain along the surface normal. This might be a subject of further research.

Another independent result is that equation (20) has an exact solution for a plausible model of the describing interfacial roughness [when the Fourier components of susceptibility are linear functions of ]. However, this solution is impractical because it engenders copious numerical problems. This is why it was not exploited in this paper.

Nevertheless, we have developed an approximate iterative method of solution of equation (20) [see equations (24) and (25)]. In the numerical examples in §5, it is shown that for realistic values of the roughness parameter the approximate solution is indistinguishable from the exact one. However, in practice it is calculated much faster than the exact one. The speed is comparable (two to three times longer) to that of corresponding calculations in the model with ideally sharp interfaces. As to the lamellar approach, to achieve sufficient accuracy we had to use about 100 lamellae and this made the calculations several hundred times longer with regard to the ideal interface model. Generally, the roughness matrix approach improves the speed of calculations by about a factor of *N*, where *N* is the number of lamellae used to model the Clearly, this is a boon when performing fitting of experimental rocking curves. The method of taking into account roughness was incorporated into RMF in multilayers [see equation (26)]. In order to take advantage of the speed of the approximate solution one has to calculate analytically an explicit form of the roughness matrix for the used model of the The results of such a calculation for the model considered in this paper are given in equations (28), (30) and (31).

A trial fitting of the experimental curve showed that the roughness matrix approach with

model type c provides better agreement with experimental rocking curves in the region of thickness fringes far from the peak of the substrate. In the vicinity of the substrate both models (nominal and type c) deviate from the experimental curve. The disagreement may be caused by instrumental effects or by the influence of defects. However, this single fit does not provide enough evidence to make any final conclusions. In the future, additional experiments should be performed, possibly with superlattices, where interfacial roughness is thought to have a stronger influence on rocking curves.It would also be desirable to make some additional independent measurements of the roughness of studied samples, for example by XRR or electron microscopy.

### Footnotes

^{1}This article will form part of a virtual special issue of the journal, presenting some highlights of the 13th Biennial Conference on High-Resolution X-ray Diffraction and Imaging (XTOP2016).

### Acknowledgements

The authors would like to acknowledge many useful discussions with Professor Dr Ilya Feranchuk.

### References

Benediktovitch, A., Feranchuk, I. & Ulyanenkov, A. (2014). *Theoretical Concepts of X-ray Nanoscale Analysis.* Berlin: Springer. Google Scholar

Boer, D. K. G. de (1994). *Phys. Rev. B*, **49**, 5817–5820. CrossRef Web of Science Google Scholar

Boer, D. K. G. de (1995). *Phys. Rev. B*, **51**, 5297–5305. Google Scholar

Boer, D. K. G. de (1996). *Phys. Rev. B*, **53**, 6048–6064. CrossRef Web of Science Google Scholar

Bolognesi, C., Kroemer, H. & English, J. (1992). *Appl. Phys. Lett.* **61**, 213–215. CrossRef CAS Web of Science Google Scholar

Epstein, P. (1930). *Proc. Natl Acad. Sci. USA*, **16**, 627–637. CrossRef PubMed CAS Google Scholar

Feranchuk, I. D., Feranchuk, S., Komarov, L., Sytova, S. & Ulyanenkov, A. (2003). *Phys. Rev. B*, **67**, 235417. Web of Science CrossRef Google Scholar

Feranchuk, I. D., Feranchuk, S. & Ulyanenkov, A. (2007). *Phys. Rev. B*, **75**, 085414. Web of Science CrossRef Google Scholar

Feranchuk, I. D., Minkevich, A. & Ulyanenkov, A. (2003). *Eur. Phys. J. Appl. Phys.* **24**, 21–26. Web of Science CrossRef CAS Google Scholar

Fullerton, E. E., Schuller, I. K., Vanderstraeten, H. & Bruynseraede, Y. (1992). *Phys. Rev. B*, **45**, 9292–9310. CrossRef CAS Web of Science Google Scholar

Glaser, E., Trombetta, J., Kennedy, T., Prokes, S., Glembocki, O., Wang, K. & Chern, C. (1990). *Phys. Rev. Lett.* **65**, 1247–1250. CrossRef PubMed CAS Web of Science Google Scholar

Ming, Z., Krol, A., Soo, Y., Kao, Y., Park, J. & Wang, K. (1993). *Phys. Rev. B*, **47**, 16373–16381. CrossRef CAS Web of Science Google Scholar

People, R., Bean, J., Lang, D., Sergent, A., Störmer, H. L., Wecht, K., Lynch, R. & Baldwin, K. (1984). *Appl. Phys. Lett.* **45**, 1231–1233. CrossRef CAS Web of Science Google Scholar

People, R. & People, R. (1986). *IEEE J. Quantum Electron.* **22**, 1696–1710. CrossRef Web of Science Google Scholar

Pietsch, U., Holy, V. & Baumbach, T. (2013). *High-Resolution X-ray Scattering: From Thin Films to Lateral Nanostructures.* New York: Springer Science and Business Media. Google Scholar

Satapathy, D., Kaganer, V., Jenichen, B., Braun, W., Däweritz, L. & Ploog, K. (2005). *Phys. Rev. B*, **72**, 155303. Web of Science CrossRef Google Scholar

Shreeman, P. K. & Matyi, R. J. (2010). *J. Appl. Cryst.* **43**, 550–559. Web of Science CrossRef CAS IUCr Journals Google Scholar

Sinha, S., Sirota, E. B., Garoff, S. & Stanley, H. (1988). *Phys. Rev. B*, **38**, 2297–2311. CrossRef CAS Web of Science Google Scholar

Stepanov, S. & Köhler, R. (1994). *J. Appl. Phys.* **76**, 7809–7815. CrossRef CAS Web of Science Google Scholar

Stepanov, S., Kondrashkina, E., Köhler, R., Novikov, D., Materlik, G. & Durbin, S. (1998). *Phys. Rev. B*, **57**, 4829–4841. Web of Science CrossRef CAS Google Scholar

Takagi, S. (1962). *Acta Cryst.* **15**, 1311–1312. CrossRef CAS IUCr Journals Web of Science Google Scholar

Takagi, S. (1969). *J. Phys. Soc. Jpn*, **26**, 1239–1253. CrossRef CAS Web of Science Google Scholar

Taupin, D. (1964). *Bull. Soc. Fr. Mineral. Cristallogr.* **87**, 469–511. CAS Google Scholar

Ulyanenkov, A. & Sobolewski, S. (2005). *J. Phys. D Appl. Phys.* **38**(10A), A235–A238. Web of Science CrossRef Google Scholar

Ulyanenkova, T., Myronov, M., Benediktovitch, A., Mikhalychev, A., Halpin, J. & Ulyanenkov, A. (2013). *J. Appl. Cryst.* **46**, 898–902. Web of Science CrossRef CAS IUCr Journals Google Scholar

Zolotoyabko, E. (1998). *J. Appl. Cryst.* **31**, 241–251. Web of Science CrossRef CAS IUCr Journals 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.