X-Ray Diffraction Tomography Recovery of the 3D Displacement-Field Function of the Coulomb-Type Point Defect in a Crystal.

A successive approach to the solution of the inverse problem of the X-ray diffraction tomography (XRDT) is proposed. It is based on the semi-kinematical solution of the dynamical Takagi-Taupin equations for the σ-polarized diffracted wave amplitude. Theoretically, the case of the Coulomb-type point defect in a single crystal Si(111) under the exact conditions of the symmetric Laue diffraction for a set of the tilted X-ray topography 2D-images (2D projections) is considered provided that the plane-parallel sample is rotated around the diffraction vector [[Formula: see text]20]. The iterative simulated annealing (SA) and quasi-Newton gradient descent (qNGD) algorithm codes are used for a recovery of the 3D displacement-field function of the Coulomb-type point defect. The computer recovery data of the 3D displacement-field function related to the one XRDT 2D projection are presented. It is proved that the semi-kinematical approach to the solution of the dynamical Takagi-Taupin equations is effective for recovering the 3D displacement-field function even for the one XRDT 2D projection.

Introduction

As is well-known[1-5], the X-ray diffraction topography method is a highly sensitive nondestructive diagnostics method for observing various crystal lattice defects like striations, grain boundaries, stacking faults, single clusters, and dislocations. All these defects change the positions of individual atoms of crystal cell relative to their regular position. The real crystal structure can be studied by the Laue diffraction topography method that provides the 2D image (2D projection) of crystal with defects. Based on the analytical approximate solutions of the X-ray optics problem by the stationary method (cf. the geometrical X-ray optics method in[6]), the theoretical fundamentals for an analysis of the X-ray diffraction scattering from macroscopical crystalline objects are elaborated and widely used (see[3,6,7], and references therein for details). In particular, in[7], in the case of the X-ray diffraction from the crystals with defects, the basic principles and theoretical analysis of the X-ray coherent diffraction intensity and incoherent (diffuse) scattering intensity, the latter is caused by the X-ray scattering from the statistical ensemble of crystal lattice defects, are comprehensively presented. In practice, to interpret and analyze the experimental 2D images of crystals with lattice defects, the defect images are compared with the corresponding ones that can be evaluated using numerical methods for solving the dynamical Takagi–Taupin equations[1-3,8-10].

In the last 20 years, the X-ray diffraction tomography (XRDT) method has been widely applied to investigate the real crystalline materials. In this method, the plate-parallel crystal sample is rotated around an axis perpendicular to a set of reflection crystal planes, namely: the rotation axis is chosen to coincide with the diffraction vector . Correspondingly, a set of the XRDT projections can be obtained at different rotation angle values. Each of these XRDT 2D projections corresponds to a certain orientation of the X-ray diffraction planes relative to the intrinsic Cartesian system of coordinates of the crystal sample as it is shown in Fig. 1.

Figure 1

Original coordinate system (X, S) for the inclined X-ray diffraction geometry; 0Y-axis is perpendicular to the (X, S)-plane. Φ is the crystal rotation angle around the diffraction vector parallel to the 0X-axis. The detector plane is perpendicular to the unit vector h = kh|/k along the diffracted wave propagation.

The first XRDT experiments have been successfully performed on the synchrotron radiation sources ESRF[11] and SPring-8[12]. Then, similar experiments were carried out on the laboratory setups using the X-ray characteristic line radiation[13,14]. Using a complete set of experimental XRDT 2D projections, the corresponding computer recovery of the 3D dislocation images has been carried out based on the modified conventional iterative tomography algorithm codes (cf.[13-15]).

It is obvious, that the 3D computer image recovery based on any experimental XRDT data is connected equally, if not to a greater extent, with the same difficulties as the interpretation of the XRDT 2D projections on the X-ray diffraction topograms. They are due to the complicated contrast mechanisms of the image formation related to various regions around defects in crystals[1-3]. In this respect, it is important to find out approximate analytical solutions of the dynamical Takagi–Taupin equations, which would allow describing some features of the X-ray diffraction topograms related to particular regions in the neighborhood of single defects in crystals. As is expected, this is a key factor for solving the inverse XRDT problem, in particular, for the 3D recovery of the displacement-field function near around single defects.

In the present study, an approximate analytical solution of the dynamical Takagi–Taupin equations is found out that seems to represent by itself a principal point to the inverse XRDT problem. Such the semi-kinematical approximation for the diffracted wave E(r) is used for recovering the 3D displacement-field function in the case of the Coulomb-type point defect (r0 is the radius-vector of a point defect in a crystal). By applying iterative simulated annealing (SA)[16] and quasi-Newton gradient descent (qNGD) algorithm codes[17,18], we will show that in the case of the symmetric Laue diffraction of the X-ray characteristic MoKα1-radiation with the diffraction vector  = [20] from a single crystal Si(111) the 3D displacement-field point defect function can be recovered.

Note that the first attempt to recover the 3D displacement-field point defect function was made in[19] using the so-called simultaneous algebraic reconstruction technique (SART) algorithm (cf.[13]). Unfortunately, the SART solutions turn out to converge only for a limited number of voxels (no more than 5 × 5 × 5 = 125 voxels) near around a point defect in a crystal.

Results and Discussion

The 2D diffraction imaging of a point defect. Semi-kinematical approximation

In this section, one derives the solution of the dynamical Takagi–Taupin equations in the semi-kinematical approximation, which seems to be effective for recovering the 3D displacement-field function using the inclined XRDT 2D projections data.

As is well-known[1-3], a direct image of defect in a crystal is due to the interbranch scattering of the X-ray Bloch waves in the strongly distorted region in the immediate vicinity of a defect, which can be interpreted as the kinematical scattering of Bloch waves propagating in a perfect crystal far from the defect central region. This assertion was confirmed by numerous numerical calculations (see, e.g.[2,3,10]).

For simplicity, in order to avoid cumbersome formulas and calculations, we restrict ourselves to the case of propagation of a σ-polarized X-ray wave-field in a non-absorbing (thin) crystal “on average” oriented in the exact Bragg position, , where 0 is the wave vector of the monochromatic X-ray plane wave incident onto a crystalline sample. In this case, the propagation of the total X-ray wave-field in a distorted crystal under conditions of the symmetric two-wave Laue diffraction is set in the form of the dynamical Takagi–Taupin equations (cf.[8,9])with the boundary conditions on the entrance crystal surface ()

Hereabove, in Eq. (1), , is the 3D displacement-field function describing the distorted crystal lattice of the crystal with defects and the oblique coordinates s0, s are directed along the wave-vectors of transmitted and diffracted waves, respectively.

Note that the boundary conditions (2) are automatically followed from the modified Eq. (1) with the added term in the right-hand side of the first line equation. Using the well-known substitutions and retaining the previous designations for the amplitudes of transmitted and diffracted waves, , one can easily show that the diffracted amplitude satisfies the inhomogeneous differential hyperbolic-type second-order equation in the partial derivatives over the dimensionless variables S0 and Sh (S0 = 0/Λ, Sh = h/Λ)

Hereafter, are the amplitudes of transmitted and diffracted waves in the crystalline medium; correspondingly, S0 and Sh are the oblique coordinates along the unit vectors of , . In the Eq. (3) the following notations for dimensionless variables S0 and Sh, complex parameters Γ, and are introducedwhere C is the X-ray polarization factor equal to 1 for σ-polarization and cos 2θ for π-polarization.

Further, we will choose a plane-parallel single crystal Si(111) sample with diffraction vector  = [20]; the wavelength λ of incident characteristic X-ray radiation equal to 0.071 nm of the characteristic X-ray Mo Kα1-radiation, photon energy of 17.48 keV. Correspondingly, the susceptibility coefficients Re(χ) = Re(χ−) = −1.921 × 10−6 and Im(χ) = Im(χ-) = 1.55 × 10−8; the Bragg angle θB = 10,65°; the X-ray extinction length , ℓ = 36.287 μm; C = 1 for the σ-polarized X-ray plane wave. The vector elastic displacement field of a Coulomb-type point defect, has a formwhere r0 is the radius vector of the point defect in a crystal.

The differential Eq. (3) describing the diffracted wave propagation through an imperfect crystal can be cast in the integral formwhere the dynamical diffracted wave amplitude within a perfect crystal and the kernel function take the form (see, e.g.[1,3] for details)respectively.

In the Eq. (7), and are the corresponding diffracted wave amplitude and the Green (point source) function in a perfect crystal; is the zero-order Bessel function of the first kind.

Further, we will use the image peculiarity of the XRDT 2D projections that are directly linked with strongly distorted regions near around a single defect core. As is shown in[1-3], in such the regions the X-ray diffraction scattering is, in general, kinematical owing to the interbranch scattering of the Bloch X-ray waves. Physically, it does mean that the so-called direct defect image contrast on the XRDT 2D projections is due to the diffracted wave propagation through the strongly distorted crystal region along the wave vector k. It immediately follows that the direct defect image on the XRDT 2D projections is formed due to the kinematical scattering.

The above assertion is equivalent to building the first-order perturbation theory solution of the dynamical Takagi-Taupin Eq. (1) if the corresponding zero-order approximation for the X-ray transmitted amplitude in a perfect crystal is used.

Thus, after some obvious straightforward manipulations with the integral Eq. (6), one obtains the following expressionfor the diffracted wave amplitude in the scope of the semi-kinematical approach, where the dynamical transmitted wave amplitude within a perfect crystal is defined byand the integral in the right-hand side of Eq. (8) is taken over the variable P along the direction of the diffracted wave propagation.

Formula (8) for the amplitude describes the kinematical scattering of the X-ray dynamical transmitted wave with the amplitude ain vicinity of the defect core, which represents by itself the pure phase object. That is why the formula (8) can be treated as the semi-kinematical approximate solution of the dynamical Takagi-Taupin Eq. (1).

Further, the theoretical formula (8) rewritten asis utilized for recovering the 3D displacement-field function together with the corresponding expression for function , namely:where Ф is the rotation angle around the diffraction vector (see Fig. 1). Note that the dimensionless X, Y, Z coordinates are linked with the intrinsic Cartesian coordinates in the crystal.

The formulas (10), (11) for the 3D displacement-field function are the basic theoretical expressions for solving the inverse problem XRDT under consideration. Finally, one has to determine the error-functional (the target function) and generate its optimization procedure. The latter has to apply to the target function in the standard form as followswhere the rectangular prism with dimensions 0 ≤ Z ≤ T, −X/2 ≤ X(T) ≤ X/2, −Y/2 ≤ Y(T) ≤ Y/2, summation in the right-hand side of Eq. (12) is carried out over an array of the XRDT 2D-projections, n is the number of the inclined XRDT 2D-projections; the dimensionless thickness T is chosen to be equal to unity.

Recovery of the 3D displacement-field function. The SA and qNGD algorithm codes

In this Section, we apply the iterative SA[16] and qNGD[17,18] algorithm codes, adapted to solving the inverse XRDT problem.

SA algorithm: setting the displacement-field function in analytical form

As is known, the SA algorithm[16] is applied to minimize nonlinear target χ2-functions. Essentially, it is one of the efficient methods for solving inverse problems with a large number of variables and combinatorial nature of iterative calculations. Starting with a model specified as the initial model and varying pseudo-randomly its parameters, the SA algorithm works until the current model fits best the data set ‘observed’ that means the minimum of the target function is achieved. The advantage of the SA algorithm seems to overcome local minima of the target function, which are the main obstacles for other nonlinear optimization methods, e.g., for the qNGD methods[17,18].

In our case, the position of a point defect is set by a radius vector r0 = nT/2, where n is the internal normal to the input crystal surface, Z = 0. Note, the crystal thickness T of a single crystal Si(111) is chosen such that the X-ray absorption in a sample can be neglected. In the first study stage, one needs to checkup that the SA algorithm code can be effective in the case of one XRDT 2D projection when only the term with Φ = 0 is taken into account in formula (12).

Correspondingly, formula (12) reduces towhere is the true X-ray topography 2D projection simulated according formula (10) with the true displacement-field function (11). The intensity is calculated according to formula (10) with trial displacement-field functions according to an iterative procedure of the target function optimization.

For simplicity, without the loss of generality, the scaling coefficient G can be set to unity in formula (11). Besides, to estimate the convergence of the iterative optimization procedure, one uses the error parameter CP, which is defined as

and it yields the quantitative estimate of the relative deviation of a current solution from the true solution .

The true XRDT 2D projection is depicted in Fig. 2. As is seen from Fig. 2, the calculated true XRDT 2D projection represents by itself the ‘double form’ image of the Coulomb-type point defect in accordance with the dipole behavior of the function (see Eq. (11), Fig. 3, and[19] for details).

Figure 2

True XRDT 2D imaging projection. The Coulomb-type point defect in a single crystal Si(111). The true displacement-field function with a descending index p = 1.5 is determined on the spatial grid nodes {15, 15, 15}, −0.5 ≤ {X, Y} ≤ 0.5. The rotation angle Φ = 0.

Figure 3

Cross-sections of the 3D displacement-field function in the spatial grid nodes {15, 15, 15}, respectively, −0.5 ≤ {X, Y} ≤ 0.5, at Z-levels: (a) Z = 0.375; (b) Z = 0.5. The upper and lower rows: the cross-sections of the true and recovered 3D displacement-field functions by the qNGD algorithm.

The initial model of the 3D displacement-filed function has been taken as the averaged superposition of functions with the descending indices {p}, i = {1 ÷ 4}. For each initial descending indices combination of p, i = {1 ÷ 4}, the operation range of descending indices in search has been chosen in the 3∙Random [Real, 0, 1] interval. For particular computer calculations, the spatial grid {41 × 41 × 41} for the Cartesian system coordinates (X, Y, Z) is used.

The recovery data of the SA algorithm code application for the 3D displacement-field function are listed in Table 1. As it follows from Table 1, in the case of the one descending index p for i = 1, the solution for the 3D displacement-field function has been obtained with the accuracy, when the error parameter CP reaches the value of 10−6. In the calculated cases of the two with combinations of descending indices p for i = {1 ÷ 3} and i = {1 ÷ 4}, the solution has been obtained with the error parameter CP of the order of 10−2.

Table 1

Initial and final recovery data for the 3D displacement-field function by the SA algorithm.

Spatial grid {i, j, k}	Initial			Final
	pi, {i = 1 ÷ 4}	Target Function	CP	pi, {i = 1 ÷ 4}	Target Function	CP
{41, 41, 41}	{0.9}	0.79	1	{1.5}	8 · 10−7	4 · 10−7
{41, 41,41}	{0.5, 1.0, 1.8}	0.48	0.96	{1.49, 1.49, 1.52}	3 · 10−4	1 · 10−2
{41, 41, 41}	{0.9,1.2,1.8,2.1}	1.44	1	{1.48,1.54,1.49, 1.49}	3 · 10−3	1.4 · 10−2

The initial 3D displacement-field function parameters {p}, i = {1 ÷ 4}. The 3D displacement-field function is determined in the analytical form.

SA, qNGD algorithm codes: numerical displacement-field function assignment

From the practical viewpoint, it is of great interest to perform the 3D displacement-field function recovering with its numerical assignment. In contrast to the case above described, we consider the 3D function values in each voxel of the spatial grid as parameters in search. Further, one uses the constraints for these functions that are due to the symmetry properties of the function over the variables X, Y, and Z − Z0. Besides, a requirement for a monotonic decrease of the function with an increase of variables |X|, |Y|, and |Z − Z0| is imposed as well.

Based on formulas (10), (11), the true XRDT 2D projection for a single crystal Si(111) sample calculated on the spatial grid {15, 15, 15} is displayed in Fig. 2. As is easily seen from Fig. 2, the XRDT 2D projection is symmetric towards the coordinate Y and shifted from the center along the X coordinate by T/2 · tan θB. The recovery data of the 3D function linked with the true XRDT 2D projection are listed in Table 2. The grid node numbers along the 0Z axis, in which the values of the function are altered in the optimization procedure, are shown in bold. Correspondingly, all the other values are chosen as the ones of the true function . As it follows from Table 2, in the case of the spatial grid {15, 15, 15} the SA algorithm code application provides the target function value to be reduced by more than six orders of magnitude, whereas the error parameter CP increases from 0.0658 to 0.118. This may occur because of the solution ambiguity of the inverse XRDT problem under consideration.

Table 2

Initial and final recovery data of the 3D displacement-field function by the qNGD and SA algorithms.

Spatial grid {i, j, k}	Initial			Final
	p ini	Target function	CP	Target function	CP
{15, 15, 1–6|7–9|10–15}	qNGD algorithm
	1.55	2.52 · 10−11	1.77 · 10−2	5.67 · 10−14	1.81 · 10−2
	1.45	1.21 · 10−11	1.71 · 10−2	6.88 · 10−21	4.90 · 10−6
	SA algorithm
	1.55	1.67 · 10−2	0.0161	1.31 · 10−5	0.0171
{15, 15, 1–5|6–10|11–15}	qNGD algorithm
	1.55	2.82 · 10−11	2.92 · 10−2	2.14 · 10−22	2.70 · 10−3
	1.45	1.47 · 10−11	2.79 · 10−2	9.56 · 10−28	6.00 · 10−4
	SA algorithm
	1.55	1.56 · 10−2	0.0273	2.94 · 10−6	0.029
{15, 15, 15}	qNGD algorithm
	1.55	2.42 · 10−11	0.0706	5.03 · 10−24	0.119
	1.45	1.24 · 10−11	0.0665·	7.99 · 10−17	0.107
	SA algorithm
	1.55	1.42 · 10−2	0.0658	2.63 · 10−8	0.118

The 3D displacement-field function is determined in the numerical form.

We compare the above3D displacement-field function recovery results with the corresponding ones by using the nonlinear qNGD algorithm based on the Levenberg–Marquardt scheme[20]. Further, one utilizes the open access NL2SNO program code as an implementation of the qNGD algorithm (see[18] for details). The corresponding optimization results of the target function (13) obtained by the qNGD algorithm are listed in Table 2. It is seen that for the spatial grids, i.e.: {15, 15, 1–6|7–9|10–15} and {15, 15, 1–5|6–10|11–15}, the final target function values of order of 10−22 and the error parameter CP ones of the order of 5 · 10−5 are obtained. It should be noted that in the case of the spatial grid {15, 15, 15}, the values of the error parameter CP are much more than the initial error CP ones due to the solution ambiguity of the inverse XRDT problem under consideration.

The recovered 3D displacement-field function is illustrated in Fig. 3, where the cross-sections of the 3D function are presented for the two values of Z equal to 0.375, 0.5, respectively.

The upper row (Fig. 3a,b) gives cross-sections of the true 3D function , the lower row gives cross-sections of the 3D function recovered by the qNGD algorithm code. The corresponding 2D cross-section images match each other better for the value of Z = 0.5. Generally, this indicates that the 3D function recovery data are much better in the nearer vicinity, Z = 0.5, of the Coulomb-type point defect core.

Conclusions

The theoretical semi-kinematical approach has been developed that describes the X-ray diffraction propagation in the vicinity of the point defect core in a crystal. To solve the inverse XRDT problem the iterative SA and qNGD algorithm codes have been applied.

In the case of the Coulomb-type point defect the recovery of the 3D displacement-field function has been obtained. With certain limitations on a class of the searched functions specified in both the analytical and numerical forms, the iterative SA and qNGD algorithm codes in use work well even for the one XRDT 2D projection.

At this point, it should be noted once more that the present semi-kinematical approach embraces the inverse XRDT problem even if alternative algorithm codes have to be applied. On the other hand, a convergence of the optimization procedure for the target function based on the SA and qNGD algorithm codes might be improved by exploiting a number of the true XRDT 2D projections and, besides, some modifications of the SA and qNGD algorithm codes in use. The question of whether the elaborated theoretical approach overall is effective or whether it should be given up in favor of other approaches, which should be applied instead, remains a good topic for future work.