research papers
A highaccuracy complexphase method of simulating Xray propagation through a multilens system
^{a}Baltic I. Kant Federal University, Kaliningrad, Russia, ^{b}Gdansk University of Technology, Gdansk, Poland, and ^{c}SaintPetersburg State University, St Petersburg, Russia
^{*}Correspondence email: spkshev@gmail.com
The propagation of Xray waves through an optical system consisting of many Xray refractive lenses is considered. For solving the problem for an electromagnetic wave, a finitedifference 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 Xray waves through a multilens 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 multilens 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 multilens 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 multilens system.
Keywords: Xray wave; Xray optic; lens; nonuniform medium; focusing; numerical method; simulation; finitedifference; stability; numerical error.
1. Introduction
Xray microscopy, a method of examining the internal structures of microobjects, 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 Xray beams (Snigirev et al., 1996; Kohn et al., 2003; Kohn, 2002). Such lenses focus Xray beams because the phase velocity of Xray 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 = , where > 0. The imaginary part of n defines an absorption factor. When an Xray 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 Xrays 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 Xray 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_{1} = , 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_{1} /N.
Obtaining a highquality image and the achievement of a high enlargement factor is one of the most important goals of the development of Xray 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 Xray 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 Xray wave propagation. We shall show that when the Xray 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 Xrays 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). 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 highaccuracy computations for problems with many lenses, we propose another approach that is based on the solution of the equation for the complex phase instead of the solution of the paraxial equation for the electric field A(x,y,z) = . The phase function 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 . 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 Xray 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 Xray 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 Xray waves, partially coherent and nonstrictly monochromatic. We analyze both known and newly developed methods for solving the equations describing the Xray 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 Xray wave propagation
The propagation of a monochromatic electromagnetic wave with frequency in a medium with complex n = is described by the Helmholtz equation (Landau & Lifshitz, 1987; Goodman, 1996; Levy, 2000; Ishimaru, 1991),
Here E is an electric field; x, y, z are space coordinates; and k = , where c is the speed of light. In particular, this equation describes the propagation of Xray waves in a medium with a complex n, and an independence of n on the coordinates is assumed.
Let us consider the case where the wave propagates along the xaxis, and the characteristic scales l_{y} and l_{z} of the wave along the axes y and z are much larger than the l_{x} of the wave along the xaxis: , . 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 xvariable, we arrive at the paraxial equation (Goodman, 1996),
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 Xrays optics but also in standard optics, in the theory of propagation of radio waves and also in acoustics (Babich & Buldyrev, 1991a,b). 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 Xray waves in a material is under consideration, then the n = depends on frequency. Moreover, and . For example, for beryllium and for a 20 keV beam, we have = 0.852 × 10^{−6} and = 0.195 × 10^{−9}.
Xray 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 = ) into one general equation (Kshevetskii & Wojda, 2014),
Here b(x,y,z) = {k[n(x,y,z)^{2}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}^{2}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
R is the radius of curvature of the lens, S_{1} is the minimum thickness of the lens and S_{2} is the maximum thickness of the lens. The function x_{m,L}(y,z) describes the lefthand surface of the mth lens. The function x_{m,R}(y,z) describes the righthand surface of the mth lens. The function b(x,y,z) for a system of N beryllium lenses is as follows,
where 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 xaxis can be taken into account with the help of a shift operation if the lens's optical axis is parallel to the xaxis. If the lens's optical axis is not parallel to the xaxis, 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 to b(x,y,z) in equation (3). This function 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 Xray 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 Xray wave propagation through a sample. Therefore, equation (3) is a general equation which is worth insightful discussion and examination.
Substituting
into equation (3), where is a complex function {i.e the wave amplitude is equal to }, we obtain an equation for a complex phase of the wave,
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 Xray 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 Xray 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 finitedifference 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 xaxis 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 xaxis.
3. Analysis of the equation obtained from the paraxial approximation for Xray 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),
Here l^{ 2} = = ; F = 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) = , = 297 µm. For simplicity we use the value = 0. Even though taking into account the process of attenuation is quite straightforward, it is of no importance for our estimates.
When the multilens system is used, then the denominator l^{ 2} 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 Xray waves being equal to 12.4 keV, then l ≃ 3 × 10^{−6} 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 in (6) should be varied insignificantly at one step h of the difference mesh. This means that the conditions 1, 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 ≃ 2 × 10^{−8} m. Within the aperture, M 25000 mesh nodes along each spatial axis must be realised. The obtained estimate of the number of nodes required for a highquality solution of the problem is not exact but it is roughly true.
The phase 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 onedimensional 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 multilens 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.
Even in the case of onedimensional lenses, computer simulation of the problem with the number of nodes takes considerable time. The time required to solve the problem with twodimensional 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 phasefunction is a usual, largevalued, but nonoscillating function. Therefore, if we simulate the behavior of the phase function , then the difference grid requirements are significantly lower, which allows us to easily reach a high accuracy of calculation of . 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 highaccuracy calculation of the Xray wave propagation through the system of lenses. We expect that this equation is especially efficient for the multilens 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 Xray wave propagation through the sample.
4. Numerical simulation of the Xray wave propagation through the system of lenses and focusing of an Xray beam
4.1. A numerical method for solving equation (3)
On the basis of equation (3), we can solve the problem of Xray wave propagation through a system of lenses and focusing of an Xray beam. There are several methods of solving the problem of Xray wave propagation through a multilens 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 finitedifference methods and considering equation (3) in the parallelepiped of sufficiently large size. This parallelepiped includes all lenses and a space before the focus, by supposition. The scheme of and the grid is shown in Fig. 2. Finitedifference 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.
Let us denote by 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,
The boundary condition = 0 on the border is imposed. Here l is a normal to the boundary .
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
where
This relation is derived by multiplying equation (7) by and by adding up all obtained relationships and the complex conjugate relations over all nodes (n,m) within . 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 nonlinear 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 nonstandard approach. Let us use the substitution
in equation (7), where Ψ_{n,m} is a complex number {i.e the wave amplitude is equal to }, and use the expansions
After simplifying, we obtain a system of ordinary differential equations,
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 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 nonlinear 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/Xrayswavepropagation.
4.3. Estimation of accuracy of numerical simulations: the modified Runge rule for evaluation of errors of numerical simulations
The accuracy of a finitedifference 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 halfmaximum 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 highprecision quality of simulations than the position of the best focus. Our goal is the development of a simulation method with tools for controlling highprecision 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_{h1} and Z_{h2} are the results of approximate calculations of Z performed with the steps h_{1} and h_{2}, respectively (h_{2} < h_{1}). 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_{h1} + Ch_{1}^{k} + O(h_{1}^{k+1}) and Z = Z_{h2} + Ch_{2}^{k} + O(h_{2}^{k+1}). We neglect the small O(h_{1}^{k+1}), O(h_{2}^{k+1}) and equate these two different expressions for Z. Then, solving the resulting equation, we obtain
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 Xray propagation and Xray 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 ( = 6π × 10^{18} s^{−1}).
(iv) δ = 2.2156 × 10^{−6}, β = 3.1801 × 10^{−10}.
(v) Form of A before lenses: A(0,y) = (1/2πσ^{2})exp(−y^{2}/2σ^{2}), σ = 297 µm.
For the case of equation (7), the results of the simulations of Xray propagation through the onedimensional 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.

The focal length is calculated as follows. Within the Runge–Kutta method, some mesh {x_{i}} of nodes on the xaxis, x_{i+1} = x_{i} + h_{x}, is used. The step h_{x} along the xaxis 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} kh^{2}. At each step x_{i} after the lens we find
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 Xrays propagating through onedimensional 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 highaccuracy simulations of the electromagnetic wave propagation through a system of 30 twodimensional beryllium lenses require long calculations with a powerful supercomputer.
Below we solve equation (5) using a finitedifference method. The results of the calculations are shown in Figs. 4 and 5. One can see that the real part of changes almost linearly, though it changes extremely fast after the set of lenses.
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.

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 highaccuracy 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 twodimensional 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 OXaxis without loss of accuracy. Therefore, we have a 10^{4}fold decrease in the amount of calculations in the case of twodimensional 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} kh^{2}). 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 Xray 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.
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 on the plane x = x_{i} completely determine the field 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 farwavefield after the focal spot, we recommend performing the calculation using the highaccuracy 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 highaccuracy 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 , known at some x_{n}, is input data to this method. This highaccuracy 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_{2}(M)] in onedimensional and twodimensional cases. For the finitedifference methods described in the paper, V = O(M_{1}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^{2}k); this means V = O[M_{1}^{ 3}L/(kd^{ 2})] in a onedimensional case and V = O[M_{1}^{ 2}L/(kd^{ 2})] in a twodimensional case, where d is the lens aperture. We use different notations M, M_{1} here because different fields are treated in computations and correspondingly different M, M_{1} 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^{5} in the onedimensional case and M = 10^{10} in the twodimensional case for highprecision 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^{6} in the onedimensional case and V ≃ 3 × 10^{11} in the twodimensional case.
For the finitedifference method developed in the paper, M_{1} = 4 × 10^{3} in the onedimensional case and M_{2} = 1. 6 × 10^{7} in the twodimensional case, for the same accuracy of computations. Correspondingly, if we take L = 0.3 m (distance to the best focus), then V ≃ 2 × 10^{6} in the onedimensional case and V ≃ 8 × 10^{9} in the twodimensional case.
We see that the complex phase method for problems with many lenses is somewhat more effective for highprecision computations than the standard FFT method. Nevertheless, for overcoming the difficulty of a great number of points for highprecision 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 complexphase equation method is new and does not have a history of improvement.
For simulation of Xray 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 Xray optics (Bahrdt, 1997; Chubar & Elleaume, 1998) to increase the calculation accuracy of a farwavefield after a focal spot could be another interesting aspect of investigation.
5. Conclusions
We have developed a highaccuracy method for simulating the propagation of Xray 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 finitedifference 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 highaccuracy 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 10^{−8} 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 twodimensional lenses, it reduces the amount of required computer memory by 625 times and reduces the computation time by more than 10^{4} times. It makes possible some simulations that are impossible without using this method.
The developed highaccuracy 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 Xray 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 ShortWavelength Difraction Theory. Oxford: Alpha Science International. Google Scholar
Babich, V. M. & Buldyrev, V. S. (1991b). ShortWavelength 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: McGrawHill. 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: ButterworthHeinemann. 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 InitialValue 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.