A high-accuracy complex-phase method of simulating X-ray propagation through a multi-lens system

Kshevetskii, S.; Wojda, P.; Maximov, V.

doi:10.1107/S1600577516013333

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 23| Part 6| November 2016| Pages 1305-1314

https://doi.org/10.1107/S1600577516013333

A high-accuracy complex-phase method of simulating X-ray propagation through a multi-lens system

S. Kshevetskii,^a ^* P. Wojda ^b and V. Maximov ^c

^aBaltic I. Kant Federal University, Kaliningrad, Russia, ^bGdansk University of Technology, Gdansk, Poland, and ^cSaint-Petersburg State University, St Petersburg, Russia
^*Correspondence e-mail: spkshev@gmail.com

Edited by G. Grübel, HASYLAB at DESY, Germany (Received 27 July 2015; accepted 19 August 2016; online 10 October 2016)

The propagation of X-ray waves through an optical system consisting of many X-ray refractive lenses is considered. For solving the problem for an electromagnetic wave, a finite-difference method is applied. The error of simulation is analytically estimated and investigated. It was found that a very detailed difference grid is required for reliable and accurate calculations of the propagation of X-ray waves through a multi-lens system. The reasons for using a very detailed difference grid are investigated. It was shown that the wave phase becomes a function, very quickly increasing with increasing distance from the optical axis, after the wave has passed through the multi-lens system. If the phase is a quickly increasing function of the coordinates perpendicular to the optical axis, then the electric field of the wave is a quickly oscillating function of these coordinates, and thus a very detailed difference grid becomes necessary to describe such a wavefield. To avoid this difficulty, an equation for the phase function is proposed as an alternative to the equation of the electric field. This allows reliable and accurate simulations to be carried out when using the multi-lens system. An equation for the phase function is derived and used for accurate simulations. The numerical error of the suggested method is estimated. It is shown that the equation for the phase function allows efficient simulations to be fulfilled for the multi-lens system.

Keywords: X-ray wave; X-ray optic; lens; non-uniform medium; focusing; numerical method; simulation; finite-difference; stability; numerical error.

1. Introduction

X-ray microscopy, a method of examining the internal structures of micro-objects, has been actively developed during the last 20 years. The application of biconcave lenses made of light metals is one of the most promising ways of focusing X-ray beams (Snigirev et al., 1996 ; Kohn et al., 2003 ; Kohn, 2002 ). Such lenses focus X-ray beams because the phase velocity of X-ray wave propagation in light metals slightly exceeds the speed of light. This is a consequence of the form of the complex refraction index of a material, n = $[1-\delta+i\beta]$ , where $[\delta]$ > 0. The imaginary part of n defines an absorption factor. When an X-ray wave propagates through the lenses the wave phase changes. The wave phase change is much larger for waves that propagate near the lens's aperture edges than for waves that propagate near the general optic axis. This leads to the effect of focusing of X-rays waves by lenses.

Recently, beryllium has been considered as one of the most promising materials for manufacturing lenses because this metal possesses low absorbing properties and a comparatively large refraction factor. However, lenses can also be manufactured from other materials.

Modern X-ray lenses have concave surfaces of elliptical, spherical or parabolic shape. In this paper the parabolic shape of lenses will be used. Such a form allows us to remove most of the aberrations inherent to lenses with spherical concave surfaces.

Since the refraction factor is small, the effect of focusing of one lens is weak, and often a system of many lenses is applied to achieve the essential focusing effect. If the lenses are arranged in a row, one after another, the lens system acts similarly as a single lens; such a system is known as a compound refractive lens (CRL). The focal distance of one lens is given by the formula F₁ = $[{{R}/({2\delta})}]$ , where R is the radius of curvature of the parabolic lens surface (Kohn, 2002). For a CRL consisting of N identical lenses, the focal distance F is approximately equal to F₁ /N.

Obtaining a high-quality image and the achievement of a high enlargement factor is one of the most important goals of the development of X-ray optics. As a consequence, the question of the influence of various lens defects on the image quality is now up to date. When we speak of defects, we mean any imperfections in the lens form and errors in the adjustment of lenses along the general optical axis of a system of lenses. We also consider any possible internal microscopic defects, i.e. cavities or inclusions of oxides in the lens's material.

This paper is devoted to the theory of numerical modeling of X-ray wave propagation in an optical system. The demands on numerical simulation methods essentially depend on our intentions. If one wants to explain the effects of focusing and image formation, then coincidence of the results of computation with experiment is obviously a necessary criterion of our simulations. Many publications of such a kind appear in the literature (Lengeler et al., 1999 ; Kohn, 2009 , 2012 ; Kohn et al., 2003).

If numerical simulations are used for scientific and engineering research (for example to investigate the influence of defects on focusing and imaging), the requirements for the theoretical validity of numerical methods is extremely high. Numerical simulations are often applied in cases where real experiments are difficult or impossible to carry out. Especially in these cases, it is good to possess some tools of control of simulation reliability and accuracy. The development of tools that monitor the accuracy of numerical modeling and investigation of the actual accuracy of the numerical simulation are important objectives of this study.

The computational complexity dependence on the optical system and on the experimental scheme is also an important issue. Normally the paraxial equation (2) (see §2) describes the X-ray wave propagation. We shall show that when the X-ray wave passes through a system of many lenses the solution of this equation is a rapidly oscillating function of distance from the main optical axis. Moreover, the oscillation frequency of the wavefield increases exponentially with the number of lenses. The higher the oscillation frequency then the more detailed the mesh required to describe such a function, and more computations are necessary. Although simulations for a small number of lenses may seem quite simple, the computational complexity increases rapidly with the number of lenses. Therefore, the simulation of X-rays propagating through many lenses may be an extremely complicated computational problem. Computer power is limited, and there is a specified limit for the quantity of lenses that can be used in simulations to obtain solutions with reasonable accuracy on the basis of the paraxial equation (Levy, 2000 ; Ishimaru, 1991 ; Babich & Buldyrev, 1991 $[Babich, V. M. & Buldyrev, V. S. (1991b). Short-Wavelength Diffraction Theory: Asymptotic Methods. Berlin: Springer Verlag.]$ ). Different approaches to how to increase the speed of calculation of synchrotron radiation are an interesting problem and investigation of this subject can be found in the literature (Bahrdt, 1997 ; Chubar & Elleaume, 1998 ).

In order to achieve high-accuracy computations for problems with many lenses, we propose another approach that is based on the solution of the equation for the complex phase $[\Psi\left(x,y,z\right)]$ instead of the solution of the paraxial equation for the electric field A(x,y,z) = $[\exp[i\Psi\left(x,y,z\right)]]$ . The phase function $[\Psi\left(x,y,z\right)]$ is not a fast oscillating function, and therefore it does not require a detailed mesh for its digitization. The computational complexity of simulations for the phase function does not depend on the quantity of lenses, and one can easily achieve a high accuracy of numerical simulations for many lenses. An electric field may be elementarily restored when the complex phase function is known.

We derive our equation for the complex phase $[\Psi\left(x,y,z\right)]$ . Then we perform calculations with the help of the equation for the complex phase, investigate the accuracy of the method and demonstrate the advantages of the suggested method for the problem of focusing of X-ray waves by 30 beryllium lenses.

The simulation method is particularly advantageous when the case of many lenses is considered. As an extreme example, the simulation of focusing of X-ray waves by a system of 160 lenses has been performed.

In this paper, we set the incident wave in the form of a simple coherent wavefront, and we do not touch on the problem of real sources of X-ray waves, partially coherent and non-strictly monochromatic. We analyze both known and newly developed methods for solving the equations describing the X-ray wave propagation. The quality of approximate methods for solving equations is not dependent on wave sources, and therefore idealization of the wave source is used and attention is focused on the problems of solving equations rather than on the problem of computation for realistic sources. Since any wave can be represented as a superposition of simple monochromatic waves, one can use the obtained solutions as a base for building waves for any complex sources. However, such questions are not considered in this paper.

2. Basic equations of X-ray wave propagation

The propagation of a monochromatic electromagnetic wave with frequency $[\omega_{0}]$ in a medium with complex refractive index n = $[1-\delta+i\beta]$ is described by the Helmholtz equation (Landau & Lifshitz, 1987 ; Goodman, 1996 ; Levy, 2000; Ishimaru, 1991),

$[{{\partial^{2}E}\over{\partial x^{2}}}+{{\partial^{2}E}\over{\partial y^{2}}}+{{\partial^{2}E}\over{\partial z^{2}}}+k^{2}n^{2}E = 0 .\eqno(1)]$

Here E is an electric field; x, y, z are space coordinates; and k = $[{\omega_{0}}/{c}]$ , where c is the speed of light. In particular, this equation describes the propagation of X-ray waves in a medium with a complex refractive index n, and an independence of n on the coordinates is assumed.

Let us consider the case where the wave propagates along the x-axis, and the characteristic scales l_y and l_z of the wave along the axes y and z are much larger than the characteristic scale l_x of the wave along the x-axis: $[l_{x} \ll l_{y}]$ , $[l_{x} \ll l_{z}]$ . In this case, equation (1) can be greatly simplified. Substituting E(x,y,z) = exp(ikx)A(x,y,z) and neglecting the small term of the second derivative of the function A(x,y,z) with respect to the x-variable, we arrive at the paraxial equation (Goodman, 1996),

$[{{\partial A}\over{\partial x}}+{{k\left(n^{2}-1\right)}\over{2i}}\,A + {{1\,}\over{2\,ik}}\left({{\partial ^{2}}\over{\partial y^{2}}}+{{\partial^{2}}\over{\partial z^{2}}}\right)A = 0. \eqno(2)]$

Equation (2) was first proposed by M. A. Leontovich in 1944 to describe the propagation of a monochromatic electromagnetic wave within the paraxial approach (Leontovich, 1944 ). As a result, the approximate equation (2) is often called the paraxial equation. This equation is widely used not only in X-rays optics but also in standard optics, in the theory of propagation of radio waves and also in acoustics (Babich & Buldyrev, 1991a $[Babich, V. M. & Buldyrev, V. S. (1991a). Asymptotic Methods in Short-Wavelength Difraction Theory. Oxford: Alpha Science International.]$ ,b $[Babich, V. M. & Buldyrev, V. S. (1991b). Short-Wavelength Diffraction Theory: Asymptotic Methods. Berlin: Springer Verlag.]$ ). Equation (2) is also often called the parabolic equation due to the external similarity of this equation with the heat equation, despite the fact that the equation contains an imaginary unit in the coefficients and, strictly speaking, is not an equation of parabolic type.

In the case of electromagnetic wave propagation in vacuum or in air, n = 1, and the second term in equation (2) disappears.

If the propagation of X-ray waves in a material is under consideration, then the refractive index n = $[1-\delta+i\beta]$ depends on frequency. Moreover, $[\delta \ll 1]$ and $[\beta \ll \delta]$ . For example, for beryllium and for a 20 keV beam, we have $[\delta]$ = 0.852 × 10⁻⁶ and $[\beta]$ = 0.195 × 10⁻⁹.

X-ray waves penetrate through almost any material and propagate practically without reflection and absorption. Therefore, the boundary condition of continuity for function A on the border of different media is acceptable and widely used (Kohn, 2002, 2009; Lengeler et al., 1999; Kohn et al., 2003). It allows us to combine equation (2) for air (n = 1) and equation (2) for lens material (n = $[1-\delta+i\beta]$ ) into one general equation (Kshevetskii & Wojda, 2014 ),

$[{{\partial A}\over{\partial x}}+b\left(x,y,z\right) \,A+{{1}\over{2ik}}\left({{\partial^{2}}\over{\partial y^{2}}}+{{\partial ^{2}}\over{\partial z^{2}}}\right)A = 0. \eqno(3)]$

Here b(x,y,z) = {k[n(x,y,z)²-1]}/2i is a complex function of the coordinates depending on the material: b(x,y,z) = 0 for air and b(x,y,z) = [k(n_B²-1)]/2i inside the lens, where n_B is the refraction index for beryllium for a given frequency.

The lenses are taken into account in the model as follows. Let x_m be the position of the center of the mth lens. Then the surfaces of the parabolic lenses are described by functions

$[\eqalign{ x_{m,{\rm{L}}}(y,&z)= \cr&\left\{ \matrix{ \eqalign{x_{m}&-({R}/{2})(y^2+z^2)\cr&-({S_1}/{2}_{\vphantom{\Big|}})}\hfill & \,\,\,\eqalign{{\rm{for}}\,\, &({R}/{2})(y^2+z^2)+({S_1}/{2})\cr&\,\lt\,{S_2}/{2}_{\vphantom{\Big|}},}\hfill \cr \eqalign{x_m&-({S_2}/{2})\cr&\hphantom{+???}}\hfill & \,\,\,\eqalign{{\rm{for}}\,\, &({R}/{2})(y^2+z^2)+({S_1}/{2})\cr&\geq({S_2}/{2}),}\hfill } \right. \cr\cr x_{m,{\rm{R}}}(y,&z)= \cr&\left\{ \matrix{ \eqalign{x_{m}&+({R}/{2})(y^2+z^2)\cr&+({S_1}/{2}_{\vphantom{\Big|}})}\hfill & \,\,\,\eqalign{{\rm{for}}\,\, &({R}/{2})(y^2+z^2)+({S_1}/{2})\cr&\,\lt\,({S_2}/{2}_{\vphantom{\Big|}}),}\hfill \cr \eqalign{x_m&+({S_2}/{2})\cr&\hphantom{+???}}\hfill & \,\,\,\eqalign{{\rm{for}}\,\, &({R}/{2})(y^2+z^2)+({S_1}/{2})\cr&\geq({S_2}/{2}).}\hfill } \right. }]$

R is the radius of curvature of the lens, S₁ is the minimum thickness of the lens and S₂ is the maximum thickness of the lens. The function x_m,L(y,z) describes the left-hand surface of the mth lens. The function x_m,R(y,z) describes the right-hand surface of the mth lens. The function b(x,y,z) for a system of N beryllium lenses is as follows,

$[b\left(x,y,z\right)= {{k\left(n_{\rm{B}}^2-1\right)}\over{2i}} \sum\limits_{m\,=\,1}^N \eta\left[x-x_{m,{\rm{L}}}(y,z)\right]\eta\left[x_{m,{\rm{R}}}(y,z)-x\right],]$

where $[\eta(x)]$ is the Heaviside function. The deviations of the lens's surfaces from a parabolic shape can be taken into account with the help of the addition of small perturbations to the continuous functions x_m,L(y, z), x_m,R(y,z). The discrepancy between the lens's optical axis and the x-axis can be taken into account with the help of a shift operation if the lens's optical axis is parallel to the x-axis. If the lens's optical axis is not parallel to the x-axis, then we can take into account this defect with the help of a rotation operation. Finally, if the lenses include cavities or inclusions, we can take into account these defects by adding some function $[\alpha(x,y,z)]$ to b(x,y,z) in equation (3). This function $[\alpha(x,y,z)]$ will take into account heterogeneity of the lens material.

Equation (3) is very convenient because it allows us to carry out thorough calculations when we solve the problem of the X-ray wave propagation through a system of lenses and in a vacuum. This equation allows us to easily take into account any inclusions in the lens material or other lens defects. Equation (3) can also be used to describe the X-ray wave propagation through a sample. Therefore, equation (3) is a general equation which is worth insightful discussion and examination.

Substituting

$[A\left(x,y,z\right) = \exp\left[i\Phi\left(x,y,z\right)\right] \eqno(4)]$

into equation (3), where $[\Phi\left(x,y,z\right)]$ is a complex function {i.e the wave amplitude is equal to $[\exp[-\Im(\Phi)]]$ }, we obtain an equation for a complex phase of the wave,

$[\eqalignno{ & {{\partial\Phi\left(x,y,z\right)}\over{\partial x}}-ib\left(x,y,z\right) + {{1}\over{2\,ik}} \Bigg[\left({{\partial^{2}}\over{\partial y^{2}}} + {{\partial^{2}}\over{\partial z^{2}}}\right)\Phi\left(x,y,z\right) \cr& \qquad + i\left\{\left[{{\partial\Phi\left(x,y,z\right)}\over{\partial y}}\right]^{2} + \left[{{\partial\Phi\left(x,y,z\right)}\over{\partial z}}\right]^{2}\right\} \Bigg] = 0 . & (5)}]$

If the experiment scheme is such that some sample is placed before or after the system of lenses, then equation (5) may also be used for simulation of the X-ray wave propagation through the sample.

S. M. Rytov was the first to propose equation (5) in 1937 when he considered diffraction of light on ultrasonic waves in a medium. In the Rytov theory (Rytov et al., 1988 ), the suggestion that the medium parameters are slowly changing with the coordinates has been used for deriving an equation in the form of (5). When we derive equation (5), the assumption of slow variation of the medium parameters was not required, so equation (5) follows directly from equation (3). The reduction of the requirements apparently became possible due to the high penetration capability of X-ray waves, which leads to a slow variation of wave parameters even when the medium parameters are changed abruptly.

S. M. Rytov developed a perturbation theory for the solution of equation (5), the method of smooth perturbations. Perturbation theory starts from a linearized equation, the contribution of the nonlinear terms being taken into account in the following orders of the perturbation theory. In any case, the Rytov perturbation theory is hardly acceptable for solving our problem because the nonlinear terms play a leading role in our problem of focusing. We solve equation (5) with a finite-difference method, without any simplifications.

If we neglect the refractive term in (5) and if we set b(x,y,z) = 0, then we arrive at the geometrical optics equation written in some special approximation. In this approximation, only the wave along the positive direction of the x-axis is taken into account and the scales of the wave along the y and z axes are much larger than the scale of the wave along the x-axis.

3. Analysis of the equation obtained from the paraxial approximation for X-ray waves propagation

The last term of the paraxial equation (3) is comparatively small (Kohn, 2009, 2012). If we neglect it, we obtain an approximate formula (6) for the function A(x,y,z). This formula describes the behavior of the wave A(x,y,z) that has passed through a system of N lenses (Kohn, 2009, 2012),

$[\eqalign{ A(x,y,z)&\approx A(0,y,z)\exp\left[i\varphi\left(y,z\right)\right],\cr \varphi\left(y,z\right)&= {{\left(y^{2}+z^{2}\right)}/{l^{\,2}}}.} \eqno(6)]$

Here l² = $[{{cR}/{\omega_{0}\delta N}}]$ = $[({{2c}/{\omega_{0}}})F]$ ; F = $[{{R}/({2\delta N})}]$ is the focus distance for a system which consists of N lenses; R is the radius of curvature of the parabolic lenses. Here A(0,y,z) is an electric field before the lenses. For estimation purposes we used function A(0,y,z) = $[\exp[-y^{2}+z^{2}/(2\sigma^{2})]]$ , $[\sigma]$ = 297 µm. For simplicity we use the value $[\beta]$ = 0. Even though taking into account the process of attenuation is quite straightforward, it is of no importance for our estimates.

When the multi-lens system is used, then the denominator l² may be small because the frequency is high and the radius of curvature of the lenses is small.

For an optical system of 30 beryllium lenses with R = 50 µm and for the energy of the X-ray waves being equal to 12.4 keV, then l ≃ 3 × 10⁻⁶ m.

When equation (3) is solved numerically, then some mesh with nodes (y_n,z_m) is always used in order to write the field A(x,y,z) at these nodes. It is a standard step which is independent of the numerical method applied for solving the problem. We denote the mesh space step by h. It is obvious that in the numerical simulations the phase $[\varphi]$ in (6) should be varied insignificantly at one step h of the difference mesh. This means that the conditions $[|({\partial\varphi}/{\partial z})h|]$ $[\ll]$ 1, $[|({\partial\varphi}/{\partial y})h|]$ $[\ll]$ 1 should be satisfied within the limits of the lens aperture. If we take into account the lens aperture d ≃ 0.5 mm, we obtain the condition $[h \ll ({l^{\,2}}/{d})]$ ≃ 2 × 10⁻⁸ m. Within the aperture, M $[\gg]$ 25000 mesh nodes along each spatial axis must be realised. The obtained estimate of the number of nodes required for a high-quality solution of the problem is not exact but it is roughly true.

The phase $[\varphi]$ becomes a fast growing function when the wave passes through the system of 30 lenses. Correspondingly, the function A(x,y,z) becomes a fast oscillating function, and therefore we have to use a very detailed difference mesh for presentation of such a function. The field Re(A) for the wave that has passed through the system of one-dimensional lenses depends on the number of lenses in the system; it is shown in Fig. 1, which clearly clarifies why a very detailed mesh should be used for the multi-lens case. A very detailed grid is required to approximate a fast oscillating function, but it is not connected to any peculiarities of the applied mathematical methods.

Figure 1
The field Re(A) immediately after the wave has passed through one beryllium lens (graph on the left), the system of ten beryllium lenses (graph in the middle) and the system of 30 beryllium lenses (graph on the right). The calculations are made with the use of equation (6)

. R = 50 µm and the energy of the X-ray waves is equal to 12.4 keV.

Even in the case of one-dimensional lenses, computer simulation of the problem with the number of nodes $[M \gg 25000]$ takes considerable time. The time required to solve the problem with two-dimensional lenses is at least M times longer.

We suggest using equation (5) for simulations because, while the function A(x,y,z) is quickly oscillating, the phase-function $[\Phi\left(x,y,z\right)]$ is a usual, large-valued, but non-oscillating function. Therefore, if we simulate the behavior of the phase function $[\Phi\left(x,y,z\right)]$ , then the difference grid requirements are significantly lower, which allows us to easily reach a high accuracy of calculation of $[\Phi]$ . Then we can easily calculate the electric field A(x,y,z) by means of formula (4).

Owing to the reasons described above, equation (5) easily allows us to fulfill a high-accuracy calculation of the X-ray wave propagation through the system of lenses. We expect that this equation is especially efficient for the multi-lens case.

If the experiment scheme is such that some sample is placed before the system of lenses, then equation (5) should also be used for simulating the X-ray wave propagation through the sample.

4. Numerical simulation of the X-ray wave propagation through the system of lenses and focusing of an X-ray beam

4.1. A numerical method for solving equation (3)

On the basis of equation (3), we can solve the problem of X-ray wave propagation through a system of lenses and focusing of an X-ray beam. There are several methods of solving the problem of X-ray wave propagation through a multi-lens system. For example, Kohn (2002) developed and applied an efficient method based on replacement of the system of lenses with one long lens with an average refraction index. Also, some approximate solution can be obtained by the consequent combining of the solutions of equation (2) without the refraction term with the solutions of equation (2) for n = 1 (Kohn, 2009). All of these approaches are acceptable, but they are not exact, and tools for controlling their accuracy and reliability are necessary. So, we attempt to develop accurate methods and to obtain reliable estimates of the method accuracy.

We solve equation (3) numerically using finite-difference methods and considering equation (3) in the parallelepiped $[\Omega]$ of sufficiently large size. This parallelepiped $[\Omega]$ includes all lenses and a space before the focus, by supposition. The scheme of $[\Omega]$ and the grid is shown in Fig. 2. Finite-difference methods may not be the most efficient but they are quite flexible and universal. They are also good for obtaining general information about the solution behavior.

Figure 2
Schematic representation of the field $[\Omega]$ of numerical integration. The figure shows only two lenses for simplicity. In simulations all 30 lenses are included in the parallelepiped $[\Omega]$ . The wave before the system of lenses has the form A(0,y,z) = $[({1}/{2\pi\sigma^2})\exp(-y^{2}+z^{2}/{2\sigma^2})]$ . Then the wave propagates in a heterogeneous medium, and lenses create the heterogeneity of the medium. The mesh is schematically pictured with parallel lines in crimson. The real mesh is more dense and contains 100000 parallel lines.

Let us denote by $[\partial\Omega]$ the boundary of the parallelepiped that is intersected by the y and z axes. We denote by S the square that is the intersection of the parallelepiped with the plane Oyz. Within the square S we define a difference grid, with the nodes and a constant step h along the axes y and z.

We approximate equation (3) by the following system of ordinary differential equations,

$[\eqalignno{ {{\partial A_{n,m}}\over{\partial x}}&+{b\left(x,y_{n},z_{m}\right)}\,A_{n,m} \cr& + {{1}\over{2ik}}\Bigg({{A_{n+1,m}-2A_{n,m}+A_{n-1,m}}\over{h^{2}}} \cr& + {{A_{n,m+1}-2A_{n,m}+A_{n,m-1}}\over{h^{2}}}\Bigg) = 0. &(7)}]$

The boundary condition $[{\partial A}/{\partial l}\vert_{\partial\Omega}]$ = 0 on the border $[\partial\Omega]$ is imposed. Here l is a normal to the boundary $[\partial\Omega]$ .

The stability and convergence of the suggested method (7) has been proved by Kshevetskii & Wojda (2014). The proof of the stability is standard for equations of this type and is based on the relation

$[{{\partial}\over{\partial x}}\left\Vert A\left(x\right)\right\Vert \leq 0,]$

where

$[\left\Vert A\left(x\right)\right\Vert = \left[\textstyle\sum\limits_{\Omega} \left\vert A_{n,m}\left(x\right)\right\vert^{2}h^{2}\right]^{1/2}.]$

This relation is derived by multiplying equation (7) by $[A_{nm}^{\ast}h^{2}]$ and by adding up all obtained relationships and the complex conjugate relations over all nodes (n,m) within $[\Omega]$ . Equation (7) is linear and approximates the original equation (3). Then convergence of the numerical solution to an exact solution of equation (3) automatically follows from the proven stability, in agreement with the Lax theorem (Richtmayer & Morton, 1967 ).

The system of ordinary differential equations (7) can be solved using any standard numerical method appropriate for this purpose: for example, we can apply the Runge–Kutta method of the second order of accuracy.

4.2. A numerical method for solving equation (5)

Equation (5), despite its many advantages, has one major drawback: it is more complicated from a mathematical point of view than equation (3). First of all, it is a non-linear equation, and, as far as we know, the mathematical theory and numerical methods for solving such nonlinear equations have not been developed until now.

To construct a numerical method that is converging, we use the following non-standard approach. Let us use the substitution

$[A_{n,m}=\exp\left(i\Psi_{n,m}\right)\eqno(8)]$

in equation (7), where Ψ_n,m is a complex number {i.e the wave amplitude is equal to $[\exp[-{\rm{Im}}(\Psi_{n,m})]]$ }, and use the expansions

$[\eqalign{ & {\eqalign{ \exp\left(i\Psi_{n\mp1,m}\right) & = \exp\left(i\Psi_{n,m}\right) \exp\left[i\left(\Psi_{n\mp1,m}-\Psi_{n,m}\right)\right] \cr& \simeq \exp\left(i\Psi_{n,m}\right)\Big[1+i\left(\Psi_{n\mp1,m}-\Psi_{n,m}\right)\cr& \quad-{{1}\over{2}}\left(\Psi_{n\mp1,m}-\Psi_{n,m}\right)^{2}\Big]_{\vphantom{\Big|}},}} \cr& {\eqalign{ \exp\left(i\Psi_{n,m\mp1}\right) & = \exp\left(i\Psi_{n,m}\right) \exp\left[i\left(\Psi_{n,m\mp1}-\Psi_{n,m}\right)\right] \cr& \simeq \exp\left(i\Psi_{n,m}\right)\Big[1+i\left(\Psi_{n,m\mp1}-\Psi_{n,m}\right)\cr& \quad-{{1}\over{2}}\left(\Psi_{n,m\mp1}-\Psi_{n,m}\right)^{2}\Big].}}}\eqno(9)]$

After simplifying, we obtain a system of ordinary differential equations,

$[\eqalignno{ & {{\partial \Psi _{n,m}}\over{\partial x}}-ib\left(x,y_{n},z_{m}\right) \cr& + {{1}\over{2ik}} \Bigg\{\!\Bigg(\!{{\Psi_{n+1,m}-2\Psi_{n,m}+\Psi_{n-1,m}}\over{h^{2}}} + {{\Psi_{n,m+1}-2\Psi_{n,m}+\Psi_{n,m-1}}\over{h^{2}}}\Bigg) \cr& + i\Bigg[{{\left(\Psi_{n+1,m}-\Psi_{n,m}\right)^{2}+\left(\Psi_{n,m}-\Psi_{n-1,m}\right)^{2}}\over{2h^{2}}} \cr& + {{\left(\Psi_{n,m+1}-\Psi_{n,m}\right)^{2}+\left(\Psi_{n,m}-\Psi_{n,m-1}\right)^{2}}\over{2h^{2}}}\Bigg]\Bigg\}=0. &(10)}]$

Convergence and stability of the system (10), at least starting from a certain sufficiently small h, follows from stability and convergence of the system (7), provided that $[\Psi(x,y,z)]$ is continuous in the variables y and z.

The system of ordinary differential equations (10) can be solved by the usual Runge–Kutta method of the second order of accuracy, as is done in this paper. For non-linear equations, different versions of the Runge–Kutta method may differ considerably in their effectiveness. Here we use the simplest version of the Runge–Kutta method of the second order of accuracy that allows us to obtain a solution if h is sufficiently small. However, mathematical estimates show that it is possible to develop a much more efficient method for solving equation (10) than we use here. We hope to develop a more effective method for solving equations (10) or (5) in the near future.

The code of all the programs for the calculations has been posted at https://sourceforge.net/projects/X-rays-wave-propagation.

4.3. Estimation of accuracy of numerical simulations: the modified Runge rule for evaluation of errors of numerical simulations

The accuracy of a finite-difference simulation can be found by means of comparison of calculations made under various steps of difference meshes. The Runge rule (Kraus & Langer, 2007 ) of practical estimation of a numerical error is popular in numerical calculations. Initially the Runge rule was formulated for approximate calculations of integrals; however, this rule is also applicable for estimating the accuracy of other quantities. The Runge rule is applicable to estimating the accuracy of quantities independent of coordinates whereas the electric field we calculate depends on coordinates. Both the position of the best focus and the FWHM (full width at half-maximum of the intensity distribution perpendicular to the lenses axis) are independent of coordinates and are important characteristics of an optical system; therefore, they can be used for estimation of simulation error with the help of the Runge rule. The position of the best focus is obtained for such an x where the intensity distribution reaches maximum (or, equivalently, where the FWHM takes minimum).

There are other important characteristics of optical systems; for example, the intensity distribution of the wavefield. Although the distribution of the field intensity is obviously calculated in our work, it cannot be used as a tool for automatically controlling the computational accuracy because it depends on coordinates. Numerical experiments show that the FWHM of the intensity distribution is significantly more sensitive to the simulation accuracy than the position of the best focus. Therefore, the FWHM is a better tool for automatic control of high-precision quality of simulations than the position of the best focus. Our goal is the development of a simulation method with tools for controlling high-precision quality. Two utilized characteristics, the FWHM and the position of the best focus, are enough to achieve this goal. However, in the future the tools for controlling the simulation accuracy can be improved, and other wavefield characteristics independent of coordinates can be included in monitoring and adjustment of the simulation accuracy.

In a classical variant of the Runge rule, calculations carried out with a fixed and double steps are used. In our case, the steps in both cases are not very different. Therefore, we modify the classical Runge rule to a more general case when the steps are not necessarily equal, i.e. h and h/2.

Let us consider an approximate calculation of some Z value that is an exact value of any of the searched characteristics. We assume that Z_h₁ and Z_h₂ are the results of approximate calculations of Z performed with the steps h₁ and h₂, respectively (h₂ < h₁). Let us also assume that the main error term Q has a structure Q = Ch^k, where h is a step, C is some constant and k denotes an accuracy order of the method. Then we can write Z = Z_h₁ + Ch₁^k + O(h₁^k+1) and Z = Z_h₂ + Ch₂^k + O(h₂^k+1). We neglect the small O(h₁^k+1), O(h₂^k+1) and equate these two different expressions for Z. Then, solving the resulting equation, we obtain

$[Q(h_1)= \left\vert Ch_{1}^{k} \right\vert = \left\vert {{Z_{h_{1}}-Z_{h_{2}}} \over {\left({h_{2}}/{h_{1}}\right)^{k}\,-\,1}} \right\vert. \eqno(11)]$

The numerical methods (7) and (10) are methods of the second order of accuracy; therefore we have to take k = 2 in equation (11).

4.4. Estimation of accuracy of numerical simulations: comparison of computational errors and volumes of calculations for both methods

We have provided computation of X-ray propagation and X-ray focusing using the parameters given below:

(i) Number of beryllium lenses N = 30.

(ii) Radius of curvature of the lenses R = 50 µm.

(iii) Energy = 12.4 keV ( $[\omega_0]$ = 6π × 10¹⁸ s⁻¹).

(iv) δ = 2.2156 × 10⁻⁶, β = 3.1801 × 10⁻¹⁰.

(v) Form of A before lenses: A(0,y) = (1/2πσ²)exp(−y²/2σ²), σ = 297 µm.

For the case of equation (7), the results of the simulations of X-ray propagation through the one-dimensional lenses depending on the number of nodes and the grid step h are shown in Table 1. Results of computation obtained with the help of equation (7) with the described parameters are shown in Fig. 3.

Table 1
Results of computation of the focal distance and FWHM for different values of h

	Number of points	h	FWHM at 0.34 m after CRL	Focal distance
1	6 × 10⁴	0.06 µm	22.8864 µm	0.368319 m
2	8 × 10⁴	0.0125 µm	20.0654 µm	0.367792 m
3	10⁵	0.01 µm	18.5089 µm	0.367297 m

Figure 3
Dependence of the absolute value [|A|]

on the distance from the optical axis. Blue: incident wave; red: results 4 cm after the set of lenses; yellow: results 20 cm after the set of lenses; green: results 36 cm after the set of lenses. Calculation provided for 100000 points.

The focal length is calculated as follows. Within the Runge–Kutta method, some mesh {x_i} of nodes on the x-axis, x_i+1 = x_i + h_x, is used. The step h_x along the x-axis must be chosen in such a way that the computational error associated with the finiteness of h_x is insignificant and does not affect the result. (This computational error is controlled with the help of the Runge rule, but is formulated with respect to the step h_x.) The requirement that the error associated with the finiteness of h_x is small leads to the fact that h_x $[\ll]$ kh². At each step x_i after the lens we find

$[a_{i}= \max\limits_{(n,m)} \left\vert A_{nm} \left(x_{i}\right) \right\vert^{2}.]$

After finalization of the calculations, we find x_i at which the maximum of a_i is achieved. The position of the best focus is determined by the value of x_i.

The results of simulations of X-rays propagating through one-dimensional lenses which depend on the number of nodes and the grid step h are shown in Table 1.

We apply the Runge rule (11) to the results of computations in Table 1 and we find out that the error of computation of the FWHM is approximately 20% and the error of computation of the best focal distance is 0.4%. Therefore, even the number M = 100000 of mesh nodes is insufficient for obtaining the FWHM with fully satisfactory precision.

We see that the computational error in the value of FWHM is about 50 times greater than the computational error in the position of the best focus; that is, the FWHM is more demanding to simulation quality, while the position of the best focus can be found even by rough calculation. Consequently, the FWHM is more suitable for controlling the accuracy of computations than the position of the best focus.

The number M = 100000 is chosen here because the calculations were performed using a personal computer (PC). The simulations for greater M are very difficult using a PC because simulations for M = 100000 already take several hours.

Thus, we come to the conclusion that more than 100000 points of a difference mesh should be used along each axis to accurately calculate the parameters of a focal spot to two significant figures. It also means that high-accuracy simulations of the electromagnetic wave propagation through a system of 30 two-dimensional beryllium lenses require long calculations with a powerful supercomputer.

Below we solve equation (5) using a finite-difference method. The results of the calculations are shown in Figs. 4 and 5. One can see that the real part of $[\Phi]$ changes almost linearly, though it changes extremely fast after the set of lenses.

Figure 4
Imaginary part of $[\Phi]$ (it can also be described by $[-\ln|A|]$ on y). Blue: incident wave; red: results 4 cm after the set of lenses; yellow: results 20 cm after the set of lenses; green: results 36 cm after the set of lenses. Calculation provided for 4000 points.

Figure 5
Real part of $[\Phi]$ . Blue: incident wave; red: results 4 cm after the set of lenses; yellow: results 20 cm after the set of lenses; green: results 36 cm after the set of lenses. Calculation provided for 4000 points.

For investigation of the calculation error, we performed simulations with 4000 points and 8000 points and the calculation error was evaluated using the Runge rule (11). The results of the calculation of the FWHM and the position of the best focus are shown in Table 2.

Table 2
Results obtained with numerical solution of equation (5) for the FWHM at distance 0.34 m from the system of lenses and the focal distance

Number of points	h	FWHM at 0.34 m after CRL	Focal distance
4000	0.25 µm	17.3315 µm	0.36575 m
8000	0.125 µm	17.3056 µm	0.365938 m

Analysis of the calculation results in Table 2 shows that the calculation error of the best focus distance is less than 0.07%. The calculation error of FWHM at a distance of 34 cm is less than 0.2%.

Thus, using equation (5) we achieved high-accuracy results, despite the fact that the grid used in this simulation contains 25 times fewer points along each direction than in the previous case. For the case of two-dimensional lenses, equation (5) allows 625 times fewer grid points to be used. The benefits are extremely high and application of equation (5) lives up to expectations.

Moreover, thanks to the use of a 25 times larger step h with respect to the previous case, we can also use at least 25 times fewer steps along the OX-axis without loss of accuracy. Therefore, we have a 10⁴-fold decrease in the amount of calculations in the case of two-dimensional lenses.

For each lens we have at least 4000 points in the direction perpendicular to wave propagation and at least 160 space steps in the longitudinal direction (in accordance with the condition h_x $[\ll]$ kh²). This enables the forms of lenses to be modeled with high precision. To clearly demonstrate the benefits of the proposed new method, we have solved the problem of propagation of the X-ray waves through 160 beryllium lenses. In this problem, the wavefield is oscillating more rapidly than in the case of 30 lenses, and the problem is computationally extremely difficult if one solves it with the help of equation (3). However, the use of the equation for the complex phase (5) makes it possible to solve this problem using a PC. The results of simulating the problem with 160 lenses are shown in Fig. 6.

Figure 6
Imaginary part of $[\Phi]$ (it can be also described by $[-\ln|A|]$ on y). Focal spot for 160 lenses (focal distance is 7 mm after the last lens). Calculation provided for 4000 points.

4.4.1. Recommended range of applicability of the complex phase method

The focal spot is very small, and only a small number of mesh nodes belong to the focal spot. Accordingly, the resolution of the mesh in the focal spot is low, and the numerical approximation error increases significantly during passage through the focal spot. The Runge–Kutta methods are sequential; that is, the values of $[\Psi(x_{i},y,z)]$ on the plane x = x_i completely determine the field $[\Psi(x_{i+1},y,z)]$ at the plane x = x_i+1. If a sequential method loses accuracy at any x_i, then the calculation accuracy cannot be high for the following points. Therefore, despite the fact that the parameters of the focal spot are calculated with high accuracy, the numerical simulation accuracy significantly drops at the distances where the wave has passed through the focal spot. For this reason, we do not recommend the use of the numerical method (10) to perform calculations at distances greater than the focal length.

If one is interested in the far-wavefield after the focal spot, we recommend performing the calculation using the high-accuracy method for equation (10) only up to distances where the wave has completely overcome all of the lenses. Then we recommend using the obtained complex phase after passing all the lenses as input data for further calculations with use of another high-accuracy method. This recommended method is based on an analytical solution of equation (2) for n = 1 and describes the propagation of the wave in a vacuum. The method is based on the Filon method of calculating rapidly oscillating integrals. It is specially designed so that the complex phase $[\Psi(x_{n},y,z)]$ , known at some x_n, is input data to this method. This high-accuracy calculation method was published by Kshevetskii & Wojda (2015 ). The use of the complex phase allows numerical calculations to be performed with high accuracy for a large number of lenses for any distance.

The question of comparing the efficiency of the phase equation method with that of Fourier optics methods which are widely used now is interesting. The volume V of calculations in the fast Fourier transform (FFT) method is estimated as O[Mlog₂(M)] in one-dimensional and two-dimensional cases. For the finite-difference methods described in the paper, V = O(M₁L/h_x), where L is the distance for which calculations need to be performed. In our most trivial algorithm, h_x might be estimated as h_x = O(h²k); this means V = O[M₁³L/(kd²)] in a one-dimensional case and V = O[M₁²L/(kd²)] in a two-dimensional case, where d is the lens aperture. We use different notations M, M₁ here because different fields are treated in computations and correspondingly different M, M₁ are necessary for achievement of the same accuracy by different methods.

We have shown in this paper that for 30 beryllium lenses it is necessary to take M = 10⁵ in the one-dimensional case and M = 10¹⁰ in the two-dimensional case for high-precision computations with any method where the electric field is calculated; in particular, it might be true also for the standard FFT method. In this case, for the FFT method, V ≃ 2 × 10⁶ in the one-dimensional case and V ≃ 3 × 10¹¹ in the two-dimensional case.

For the finite-difference method developed in the paper, M₁ = 4 × 10³ in the one-dimensional case and M₂ = 1. 6 × 10⁷ in the two-dimensional case, for the same accuracy of computations. Correspondingly, if we take L = 0.3 m (distance to the best focus), then V ≃ 2 × 10⁶ in the one-dimensional case and V ≃ 8 × 10⁹ in the two-dimensional case.

We see that the complex phase method for problems with many lenses is somewhat more effective for high-precision computations than the standard FFT method. Nevertheless, for overcoming the difficulty of a great number of points for high-precision calculations with the FFT method, Chubar et al. (2007 , 2008 ) and Canestrari et al. (2014 ) suggested some powerful formal transformation allowing mitigating very effectively the effect of the fast phase increase and to reduce M considerably by means of analytical treatment of the quadratic phase term. This transformation, however, demands preliminary knowledge of approximate values of the wavefront radius and the center. We hope also that the effectiveness of the phase equation method will increase with time; the modern FFT is a result of a long evolution and optimization, but the complex-phase equation method is new and does not have a history of improvement.

For simulation of X-ray wavefront propagation through a CRL, the SRW code exists which utilizes the Fourier optics and compatible methods (Bahrdt, 1997; Chubar & Elleaume, 1998).

The combined use of the complex phase method with methods conventional for X-ray optics (Bahrdt, 1997; Chubar & Elleaume, 1998) to increase the calculation accuracy of a far-wavefield after a focal spot could be another interesting aspect of investigation.

5. Conclusions

We have developed a high-accuracy method for simulating the propagation of X-ray waves through many lenses and focusing the waves. The simulation method is equipped with tools for monitoring the accuracy of the simulations. The simulation method uses a finite-difference approximation of the proposed equation for the complex phase of the wave.

We have also estimated the number of mesh nodes necessary to obtain a high-accuracy solution by conventional methods based on the solution of the parabolic equation. We have shown that for the case of 30 beryllium lenses the necessary dimension of the mesh is more than 100000 nodes along each of the axes perpendicular to the optical axis of the lens system (the mesh step must satisfy the condition h $[\le]$ 10⁻⁸ m). The use of the equation for the complex phase allows the mesh dimension to be reduced by at least 25 times along each axis perpendicular to the main optical axis. In the case of two-dimensional lenses, it reduces the amount of required computer memory by 625 times and reduces the computation time by more than 10⁴ times. It makes possible some simulations that are impossible without using this method.

The developed high-accuracy method of simulation is designed to investigate the influences of various defects in the lenses (cavities, inclusions, violations of forms) on focusing and image. Also, the new method allows solving problems with the complex shapes of lenses, or calculating the propagation of X-ray waves through a sample with complex internal structure.

Acknowledgements

We thank Professor S. Leble for his advice and assistance in the research. We also thank A. Snigirev for fruitful discussions and offer of research topics. We are grateful to V. Kohn for discussions and useful advice and comments. The work is supported by Ministry of Education and Science of the Russian Federation (Contracts No. 14.Y26.31.0002 and 02.G25.31.0086).

References

Babich, V. M. & Buldyrev, V. S. (1991a). Asymptotic Methods in Short-Wavelength Difraction Theory. Oxford: Alpha Science International. Google Scholar
Babich, V. M. & Buldyrev, V. S. (1991b). Short-Wavelength Diffraction Theory: Asymptotic Methods. Berlin: Springer Verlag. Google Scholar
Bahrdt, J. (1997). Appl. Opt. 36, 4367–4381. CrossRef PubMed CAS Web of Science Google Scholar
Canestrari, N., Chubar, O. & Reininger, R. (2014). J. Synchrotron Rad. 21, 1110–1121. Web of Science CrossRef IUCr Journals Google Scholar
Chubar, O., Couprie, M., Labat, M., Lambert, G., Polack, F. & Tcherbakoff, O. (2008). Nucl. Instrum. Methods Phys. Res. A, 593, 30–34. Web of Science CrossRef CAS Google Scholar
Chubar, O. & Elleaume, P. (1998). Proceedings of the Sixth European Particle Accelerator Conference (EPAC98), pp. 1177–1179. Bristol: Institute of Physics. Google Scholar
Chubar, O., Polack, F., Couprie, M.-E., Labat, M., Lambert, G. & Tcherbakoff, O. (2007). Proceedings of FEL 2007, pp. 192–195, Novosibirsk, Russia. Google Scholar
Goodman, J. W. (1996). Introduction to Fourier Optics, 2nd ed. New York: McGraw-Hill. Google Scholar
Ishimaru, A. (1991). Electromagnetic Wave Propagation, Radiation and Scattering. Englewood Cliffs: Prentice Hall. Google Scholar
Kohn, V. G. (2002). J. Environ. Toxicol. Public Health Lett. 76, 701–704. Google Scholar
Kohn, V. G. (2009). J. Surf. Invest. 5, 32–39. Google Scholar
Kohn, V. G. (2012). J. Synchrotron Rad. 19, 84–92. Web of Science CrossRef CAS IUCr Journals Google Scholar
Kohn, V., Snigireva, I. & Snigirev, A. (2003). Opt. Commun. 216, 247–260. Web of Science CrossRef CAS Google Scholar
Kraus, J. & Langer, U. (2007). Editors. Lectures on Advanced Computational Methods in Mechanics. p. 164. Berlin: de Gruyter. Google Scholar
Kshevetskii, S. P. & Wojda, P. R. (2014). Proc. SPIE, 9209, 92090Q. Google Scholar
Kshevetskii, S. P. & Wojda, P. (2015). Math. Appl. 43, 193–206. Google Scholar
Landau, L. D. & Lifshitz, E. M. (1987). The Classical Theory of Field, 4th ed. Oxford: Butterworth-Heinemann. Google Scholar
Lengeler, B., Schroer, C., Tümmler, J., Benner, B., Richwin, M., Snigirev, A., Snigireva, I. & Drakopoulos, M. (1999). J. Synchrotron Rad. 6, 1153–1167. Web of Science CrossRef IUCr Journals Google Scholar
Leontovich, M. A. (1944). Izv. Akad. Nauk. SSSR Ser. Phys. 8, 16. Google Scholar
Levy, M. (2000). Parabolic Equation Methods for Electromagnetic Wave Propagation. The Institution of Electrical Engineers, London, UK. Google Scholar
Richtmayer, R. & Morton, K. (1967). Difference Methods for Initial-Value Problems. New York: Interscience. Google Scholar
Rytov, S. M., Kravtsov, Yu. A. & Tatarskiy, V. I. (1988). Introduction to Statistical Radiophysics. Berlin/Heidelberg: Springer. Google Scholar
Snigirev, A., Kohn, V., Snigireva, I. & Lengeler, B. (1996). Nature (London), 384, 49–51. CrossRef CAS Web of Science Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.