High-accuracy computation of hard X-ray focusing and imaging for refractive optics

A mathematical apparatus for solving problems of X-ray wave propagation through complex optical systems, when the lens thickness can change with jumps, is developed and presented. The developed method is based on the use of the superposition of oriented Gaussian beams, which satisfy the Helmholtz equation with high accuracy. The wave propagation in air and through kinoform and ordinary lenses is considered. Focusing and imaging properties are compared for both types of X-ray optics. The diffraction effects arising due to thickness jumps in the kinoform lenses and the influence of these jumps on the X-ray focusing and imaging are investigated. The prospect of using the developed theory for X-ray optics applications is discussed.

1. Introduction

Research with X-ray radiation is rapidly expanding, with many new highly coherent and extreme-brightness radiation sources being created globally (Cho, 2020). Today the new generation of synchrotrons plays a leading role in a wide range of scientific and technological applications. However, high-quality coherent-relative optics are still vital for the full utilization of highly coherent and intense light from diffraction-limited X-ray sources. Additionally, the development of precise and high-efficiency X-ray optics requires extreme precision and almost nanometre-accuracy from both methods and instruments. Therefore, the critical issue remains the ability to predict the interaction of radiation with real materials and elements of X-ray optics (Roth et al., 2017; Lyatun et al., 2020) mathematically, with controlled calculation accuracy for highly coherent radiation.

Currently, coherent imaging based on X-ray optics is developing rapidly due to the need for quick visualization of micro- and nano-scale objects with high spatial resolution. To visualize and study the internal structure of small objects, it is necessary to have the appropriate optics to enable us to control, collimate, transport and focus hard X-rays with high accuracy and precision. Such optics for hard X-rays were proposed in 1996 (Snigirev et al., 1996), which now copes with the task of utilizing coherent radiation in the world's leading synchrotrons. However, almost all the materials used for X-ray optics are not homogeneous, which leads to a sharp change in the phase of the radiation as it passes through the optical material (Roth et al., 2017). In this paper, we propose an approach that provides a highly accurate solution to the problem of calculating the scattering of short-wave radiation by elements of refractive X-ray optics. Consideration is given to the electrodynamic approach.

Here, we will investigate the mathematical apparatus developed for the analysis of the optical properties of two types of X-ray optics with different shape profiles – concave parabolic and kinoform. The latter type is of particular interest because it has a sharp stepwise change in the lens profile. When an X-ray wave propagates through a kinoform lens, wavefield oscillations of significant amplitudes but very fine scales occur at the lens thickness jumps. With distance from the lens, these very fine oscillations of the wavefield disappear, but they result in changing the physical system properties as a whole. We expect similar physical effects when X-ray waves propagate through lenses with internal defects (inclusions, cavities) or through porous materials (Goikhman et al., 2015). This sparks interest in the consideration of kinoform lenses as simpler systems of a similar kind. We believe that the study of the propagation of X-ray waves through such inhomogeneous structures is important for the development of the theoretical apparatus for X-ray optics. At present, such systems are very difficult in their theoretical description and calculations. This situation necessitates interest in the theoretical consideration of such problems.

The technique of describing the propagation and refraction of X-ray waves, applied and developed in this paper, is based on the so-called Gaussian beams. The technique of using rays in optics is presented in the literature (Kogelnik, 1965; Keller & Streifer, 1971; Deschamps, 1972; Deschamps et al., 1983). Gaussian beams are also used in radiolocation (Chabory et al., 2005, 2010; Ghannoum et al., 2009). The theoretical application of Gaussian beams to problems of X-ray optics was presented by Wojda & Kshevetskii (2019).

It is also worth pointing out that the literature contains a description of the asymptotic solutions of the Helmholtz equation (Maslov, 1994), which can be interpreted as beams, and which are more general in comparison with the Gaussian beams used in this work. Gaussian beams have been used to solve various short-wave physical problems (Popov, 2002; Anikin et al., 2019) but not for X-rays.

Initially, Gaussian beams were understood in the literature as exact solutions of the paraxial equation; they are widely used in modelling electromagnetic wave propagation (Good­man, 1996; Kohn, 2012). In this approach, electromagnetic waves are described with a superposition of Gaussian beams. Some development of this technique, adequate to the peculiarities of X-ray waves, was suggested by Wojda & Kshev­etskii (2019). Instead of conventional Gaussian beams, it was suggested to use Gaussian beams multiplied by the plane wave of the carrier frequency. Moreover, the use of oriented Gaussian beams, in which the directions of wave propagation in beams coincides with the local X-ray wave propagation directions at the wavefront, was suggested. The condition is that the direction of propagation of the beams should be perpendicular to the front of the modelled wave and determines the tilt angles of the beams. Thus, the developed approach combines the properties of geometrical optics and wave optics. As soon as Gaussian beams oriented in space are used, it is already more appropriate to consider such wave beams not as solutions of the paraxial equation but as approximate solutions of the more general and more exact Helmholtz equation.

The authors consider the proposed method to be highly accurate not so much because the Helmholtz equation is being solved but rather to emphasize that some automatic procedure for finding the required number of Gaussian beams for high-quality calculations is suggested. Over time, this computation tuning procedure may be improved. However, a reasonable automatic choice of the method parameters ensuring the accuracy and reliability of simulations is proposed and implemented in this paper.

A famous book by Gerrard & Burch (1975) is devoted to the mathematical apparatus of conventional optics. In a significant part of this book, the method of beams, in reality, a certain version of geometric optics, is presented. Compared with the book, this paper uses a wave approach. The Gaussian wave beam is given in Gerrard & Burch (1975) and is discussed in detail. However, the methods based on the superposition of Gaussian beams being developed and used in this paper are not developed or discussed in Gerrard & Burch (1975). In the calculations of this paper we use the representation of the wavefield in terms of the sum of thin Gaussian beams. The amplitudes and directions of the introduced Gaussian beams are calculated from the electric wavefield digitized on the mesh on the plane, while in Gerrard & Burch (1975) the wave propagation paths from the source are calculated for the entire optical system at once based on the optical system configuration and its important optical characteristics. In this paper, in front of each lens and after each lens, in front of the sample and after the sample, the amplitudes and propagation directions of wave beams are calculated anew from the analysis of the wavefield, and the symmetry of the optical system is not necessary. The proposed method in this paper allows us to take into account irregular surfaces or internal defects in the lens or sample. Finally, in this work, X-ray waves not only pass through the lenses but also scan the sample and draw its image, and the sample is included in the problem.

Currently, the Fourier method is often used to calculate the propagation of X-ray waves in air. Programs that use it include, inter alia, Synchrotron Radiation Workshop (Chubar & Elleaume, 2013; Kohn, 2020). The Fourier method is especially popular due to the very fast standard computer FFT library allowing the Fourier transform to be performed very quickly. Despite the obvious differences between the proposed Gaussian beam method and the Fourier method, these two different approaches have much in common. One can interpret the Gaussian beam method as some kind of generalized Fourier method. Within this, some special non-trivial set of basis functions, adapted to the specific considered problems, is used instead of the complex exponentials of the usual method. The adaptation of the basis functions to specific considered problems determines the method's effectiveness.

Additionally, the focusing of X-rays by beryllium and diamond lenses will be considered in this paper since beryllium is a widespread optical material and diamond is the most promising material for X-ray optics and new ultra-bright sources with high thermal load because of its high heat conductivity (Isakovic et al., 2009).

The proposed method allows calculating the smallest effects on the focal plane resulting from X-ray diffraction at the edges of the lens aperture (Wojda & Kshevetskii, 2019; Medvedskaya et al., 2020). Hence the idea to investigate the diffraction effects on the thickness jumps of a kinoform lens arose. This work is devoted to the development of a more advanced mathematical apparatus for theoretical studies and the application of the developed methods to the calculation of X-ray focusing by a kinoform lens system. We hope that the developed approach can be easily adapted to real microstructured and porous materials of X-ray optics (Goikhman et al., 2015; Lyatun et al., 2020).

2. Basic equations and starting points

X-ray waves usually propagate at small angles to some specific main direction of X-ray propagation. We further consider the propagation of monochromatic X-ray waves, which are used in X-ray microscopy. The monochromaticity of waves allows an equation to be based on the Helmholtz equation, which directly follows from Maxwell's equations where the wave frequency ω₀ is fixed (Levy, 2000; Ishimaru, 1991; Babich & Buldyrev, 1991a,b),

In equation (1), E is the electric field; r = (x, y, z), k₀ is the wavenumber, c is the speed of light; n = 1 − δ + iβ is a complex refractive index, β is an absorption coefficient, and δ is a refractive coefficient (β, δ are non-negative).

In the case of electromagnetic wave propagation in vacuum or in air, n = 1. If the propagation of X-ray waves in a material is under consideration, then the refractive index n = 1 − δ + iβ depends on frequency. Moreover, δ << 1 and β << δ. For example, for beryllium and for a 20 keV beam, we have δ = 0.852 × 10⁻⁶ and β = 0.195 × 10⁻⁹.

The electric field E(r) is assumed to be continuous at the air–material interface.

In X-ray optics, an approximate so-called paraxial equation is often used, which follows from (1) and which was first derived by Leontovich (1944). The paraxial equation takes into account that waves propagate mainly along one specific direction. The paraxial equation is often solved using Fourier transform methods; finite-difference methods have also been proposed (Kshevetskii et al., 2016). However, finite-difference methods require a very dense computational grid, which makes the calculations time-consuming.

In this paper, the Helmholtz equation (1) is directly solved, without simplifications. The solution method is based on the use of oriented Gaussian beams and was proposed by Wojda & Kshevetskii (2019). The problem under consideration has some specifics and is more complicated than in Wojda & Kshevetskii (2019), so further development of the method will be required.

2.1. X-ray propagation in a vacuum: approximation of the X-ray wave electric field with a superposition of oriented Gaussian beams

Waves in X-ray optics propagate at small angles to the main direction of propagation. Let the X-rays propagate mainly along the OX axis, not exactly along this axis but at a small angle to it, and let the vector e₁ indicate the direction of wave propagation. Let the vector j be directed along the axis OZ. We denote e₂ = j × e₁, e₃ = e₁ × e₂. The approximate solution of the Helmholtz equation within which the wave propagates in the direction indicated by the vector e₁ can be written as follows (see, for example, Wojda & Kshevetskii, 2019),

Here r = (x, y, z) and r₀ = (x₀, y₀, z₀) is a point in space through which a Gaussian beam passes and this point lies on the symmetry axis of the Gaussian beam. σ is a parameter that determines the Gaussian beam width. The product (r − r₀) e_j is a dot product of the vectors r − r₀ and e_j, j = 1, 2, 3.

Under the condition (k₀σ)² ≪ 1, which is obviously satisfied by a large margin in our study, the formula (2) gives an almost exact particular solution to the Helmholtz equation (1). Therefore, it was proposed in Wojda & Kshevetskii (2019) to construct a general solution of equation (1) in the form of a superposition of particular solutions (2) with appropriately chosen parameters.

Let us consider the following general problem. Let a known electric field E₀(y, z) of an electromagnetic wave be given on the plane x = x₀. Using the Helmholtz equation, it is necessary to calculate the electric field in the entire half-space .

We introduce in space a system of lines given by the equations (y = y_m, z = z_n), where y_m and z_n are real numbers, and y_m+1 = y_m + h, z_n+1 = z_n + h, where h is the step of the introduced mesh κ. In Wojda & Kshevetskii (2019), it is proposed to search for a general solution to the problem of the X-ray wave propagation in the form

Here e_1,mn are the unit vectors that determine the directions of wave propagation in each beam. They are chosen so that they should be locally perpendicular to the wavefront of the wave and directed along the local direction of wave propagation. The formulas for e_1, mn are given in Section 3; they depend on the form in which the wavefield is represented. The coefficients E_0mn of the superposition (3) have been calculated in Wojda & Kshevetskii (2019); it is shown that E_0mn = E₀(y_m, z_n).

In Wojda & Kshevetskii (2019) it is shown that it is wise to choose σ = h.

Further in the text, we consider both cases: the planar one and the three-dimensional one. For the planar case, the Gaussian beam looks like

The wavevector e₁ indicates the direction of propagation of X-ray waves in the beam.

The general solution for a planar case of the problem of X-ray wave propagation has the form

For the planar case, the electric field from the X-ray source is given by the formula

whereas for the three-dimensional case

In the works of Popov (1982, 1987, 2002) and Grikurov & Popov (1983), the wavefield of short waves (not X-ray waves, but described by the Helmholtz equation) in an inhomogeneous medium was considered. The wavefield was generated by a point source. For the solution, some representation of the wavefield as a sum over Gaussian beams similar to (3) and (5) (or as an integral over angular variables) was used. The main difference in the case here under consideration is obvious in that we do not have an analytically specified point source of waves but instead we have an arbitrarily specified electric field of X-ray waves digitized at points r_mn on the plane x = x₀ in the three-dimensional case, or digitized on the grid on the line x = x₀ in the two-dimensional case. This electric wavefield is decomposed into a superposition of Gaussian beams based on the fact that these Gaussian beams, starting from the plane x = x₀, behave like Gaussian functions at the section x = x₀ and play the role of basis functions. Expansion coefficients and direction vectors of Gaussian beams were derived by Wojda & Kshevetskii (2019). In (3) and (5), E_0mn and E_0m are the values of the digitized electric wavefield at the points r_mn on the plane in the three-dimensional problem or at the points r_n on the straight line in the two-dimensional problem, respectively. The direction vectors of Gaussian beams will be given later.

3. X-ray propagation through optical elements

Let the X-rays propagate mainly along the axis OX. We will use some simplifications. We will not describe in detail the X-ray propagation in the thickness of each lens, but we will apply some formulas giving the wavefield after each lens. Also, we neglect very small corrections as the wave can only really propagate at a very small angle to the OX axis. Finally, we use only sufficiently thin Gaussian beams, which allow us to take into account the lens shapes in detail. In Wojda & Kshevetskii (2019), it is shown that when the X-ray wave passed through a lens it acquired a local phase shift,

In (8), E(x_C, y, z)⁻ is the electric field in front of an optical element and E(x_C, y, z)⁺ is the field after the optical element. That is, if an X-ray wave passes through an optical element with a local thickness F(y, z), then at the exit from the optical element the X-ray wave acquires an additional phase ϕ(y, z) = k₀(−δ + i β) F(y, z). It is also equivalent to the function describing the electric field gaining a factor .

Also, the wave changes its propagation direction. The changes in the components of the vector e₁ are given by the formulas (Wojda & Kshevetskii, 2019)

The formulas (8) and (9) are applied to Gaussian beams. Namely, on the surface in front of each lens, the electric field is represented as (3). After the lens, each Gaussian beam acquires a factor according to (8) and changes its propagation direction according to (9).

In specific calculations, the partial derivatives in the formulas (9) for propagation directions of Gaussian beams are approximately calculated using finite-difference methods. For example,

3.1. Influence of jumps in the thickness of a kinoform lens and lens edges on X-ray wave propagation: substantiation of calculation methods in the case of a non-smooth lens surface

Idealized lenses may not have a width, and then their form can be expressed as follows: F(y, z) = R(y² + z²) + W_sm. However, all real objects have limited sizes. For a lens, this is particularly evident in that the lens aperture is limited. The thickness of a lens can be written using the formula F_L(y,z) = min{R (y²+z²)+W_sm,W}, where R denotes the lens curvature, W is the maximum thickness of the lens and W_sm is the minimum thickness of the lens (see Fig. 1).

Figure 1 The shape of a parabolic biconcave lens and a kinoform lens. For the biconcave parabolic lens FL(y,z) = min{R (y^2+z^2)+Wsm,W}. For the kinoform lens , WK is the kinoform lens thickness, l = 0, 1,…, W/WK − 1.

Formulas (9) contain derivatives of lens thickness. Some of these derivatives may not exist if the point (y_m, z_n) lies on the lens thickness jump. In the case of a conventional lens, this is not of great importance, since in experiments the intensity of X-ray radiation incident on the lens edge or the lens aperture edge is small in comparison with the maximum radiation intensity. The publication by Wojda & Kshevetskii (2019) shows that the lens edges and the lens aperture edges generate diffraction effects.

Such diffraction effects are even more significant for kinoform lenses because the kinoform thickness at the jumps changes rapidly, and the radiation intensity falling onto the sharp thickness change is comparable with the maximum intensity.

According to the Huygens–Fresnel principle, spherical waves should form in places where there is a sharp change in the kinoform thickness. However, it is impractical to use this principle in calculations, since this requires a huge amount of calculations. The use of Gaussian beams, whose propagation directions coincide with local wave propagation directions, allows calculating the solution with high accuracy using a relatively small number of Gaussian beams. This greatly reduces the volume of required computations.

However, the question of how to correct the formulas (9) and (10) so that they can be applied to kinoform lenses arises. Additionally, it is necessary to mathematically substantiate the applicability of the proposed mathematical apparatus as a whole to kinoform lenses.

First, we correct the formulas (9) and (10). Kinoform lenses can be viewed as consisting of concentric rings whose thickness is described with differentiable functions. We use this circumstance and replace the finite-difference approximations of the derivatives in (10) with any other finite-difference approximations within which the used points lie within one such concentric ring. Then the formulas (9) and (10) will be correct even in the case of kinoform lenses.

Second, let us consider the applicability of the proposed mathematical apparatus in general to kinoform lenses. It is enough to consider only one jump in lens thickness. The essence of the estimates given below is that very narrow Gaussian beams are used when passing through the lens, and there are a lot of them. Only a small number of beams fall into the lens thickness jump, and, as can be seen, they practically do not affect the overall wave pattern. So, let the thickness jump of a three-dimensional kinoform lens take place on a circle of radius R_R. Consider the ring (R_R − ε, R_R + ε). The upper bound Q of the contribution of this ring to the total wave electric field is

where |E|_M is the electric field amplitude. In obtaining this estimate, we took into account that very narrow beams are used when crossing the lens, and they intersect weakly. Therefore, the contribution to the wavefield near the lens at any point is practically the contribution of a single Gaussian beam. With increasing distance from the lens, the beams expand significantly, and together they create a wave pattern. However, the beam amplitudes decrease when the beam widths increase, and therefore the contribution of the ring to the overall wave pattern remains the same as near the lens.

Let us consider the limit , . [Consideration of the limit is possible from a mathematical point of view. However, the Gaussian beam (2) is an approximate asymptotic solution to the Helmholtz equation (1) only if (k₀σ)² ≫ 1. In our calculations, (k₀σ)² > 10⁹. Therefore, we assume that σ is very small compared with the lens size and compared with the distance between thickness jumps, but in the considered limit.] It is also reasonable to take . One can see that in the considered limit. This also means that, although the thickness jump certainly affects the wave pattern, nevertheless fine details near the thickness jump do not affect the wavefield. For example, if you slightly smooth the thickness jump, the result will be very similar.

4. Calculation of local propagation directions of an electromagnetic wave

The use of oriented Gaussian beams allows one to calculate the propagation of waves with high accuracy. To apply the method, one must first represent the electromagnetic wave in the form of a superposition of oriented Gaussian beams. In particular, one needs to efficiently calculate the propagation directions of the Gaussian beams, which by definition coincide with the local wave propagation directions.

In solving the problem, two different situations are encountered in which the directions of local wave propagation need to be calculated. In the first case, the electromagnetic wave is represented with a relatively simple analytical expression (for example, this situation takes place when the wave from the source falls on the lens). We represent the X-ray wavefield in the vicinity of the point r_mn = (x₀, y_m, z_n) in the form

where ϕ(x, y, z) is a real local phase of the wave. In this representation of the wavefield, the local phase function ϕ(x, y, z) describes the entire dependence of the wavefield on the coordinates. From (12), the following formula follows,

Therefore, the phase function gradient can be found based on (13),

Here denotes the real part of a complex number. The phase function gradient (14) gives the local propagation directions of the electromagnetic wave, and therefore it is reasonable to define the propagation directions of the introduced Gaussian beams by the formula

Partial derivatives on the right-hand side of (14) can be approximately calculated based on the values of the field E on the grid, using finite-difference approximations of the derivatives.

In reality, unfortunately, the formula (15) is only useful in simple cases. For example, when an X-ray wave from a source falls on a lens, the wave propagates almost along one direction, and this case greatly simplifies the calculations. At the same time, usually, an X-ray wavefield is described with very rapidly changing functions. The wave propagation directions at different points can vary quite significantly. Altogether, this may bring some situations when local directions of wave propagation cannot be calculated with high accuracy based on the formula (15) without use of dense calculation mesh.

In this case, a description of such a rapidly changing wavefield with a superposition of Gaussian beams is adequate, and some acceptable formula is necessary.

The approximate relation follows from (2),

This approximate relation is valid because in (2) the derivatives of the first exponential factor are much larger than the derivatives of the remaining factors. So we have the formula

Equation (17) follows from equations (13), (14) and (16). Further, using the formulas (14) and (15), we come to the formulas for calculating the local wave propagation directions,

and

where e_1, mn = (e_1, mn, x, e_1, mn, y, e_1, mn, z).

The formulas (18) allow us to decompose the wavefield, represented as a superposition of Gaussian beams, into a new superposition of Gaussian beams. For example, if a wave presented in the form of a superposition of Gaussian beams falls on a lens, then one needs to apply such a decomposition because the Gaussian beam widths grow with distance, but one needs to pass through the lens only with very narrow Gaussian beams in order to take into account the lens shape in detail. The formulas (18) allow one to perform such a decomposition.

5. Modelling of propagation and focusing of X-ray waves, imaging with X-ray waves: experimental conditions and technique for solving the problems

We consider two cases: planar and three-dimensional propagation of X-rays through some optical systems. The propagation and focusing of X-ray waves are studied for the following conditions. Let us consider the situation when a source of monochromatic X-ray waves with a wavelength λ₀ = 10⁻¹⁰ m produces in the plane x = x₀ () a Gaussian distribution of electric wavefield intensity with the width = 100 µm. More specifically, on the plane x = x₀, the X-ray source creates an electric field

in the planar case and

in the three-dimensional case. Then the X-ray radiation propagates further to a distance of 40 m.

Then, for the planar case, the propagation of X-ray waves through a lens (lens system) made of beryllium and the wave propagation after the lens (lens system) is considered. The general scheme for calculating the X-ray wave propagation is shown in Fig. 2. For beryllium lenses and the given X-ray wavelength, β_Be = 3.1801 × 10⁻¹⁰ and δ_Be = 2.2156 × 10⁻⁶. The lens curvature radius is 50 µm, and the smallest width is 30 µm.

Figure 2 Optical system for the planar case (top) and the general scheme for calculating the X-ray wave propagation (bottom).

For the three-dimensional case, the following problem is considered: propagation of an X-ray wave from a given wavefield on the x = x₀ plane towards the grating at a distance of 40 m, propagation through a copper grating (β_Cu = 8.8546 × 10⁻⁷ and δ_Cu = 1.1015 × 10⁻⁵), then wave propagation to a distance of d_X meters, up to a diamond lens (β_D = 2.8366 × 10⁻⁹ and δ_D = 4.7423 × 10⁻⁶), propagation through the lens and, finally, propagation to the screen at a distance of d_Y from the lens. The general scheme for calculating the X-ray wave propagation is shown in Fig. 3.

Figure 3 Optical system for the three-dimensional case (top) and the general scheme for 3D calculations of the X-ray wave propagation (bottom).

The planar problem is solved as follows. The specified wave electric field obtained experimentally from the source of monochromatic X-ray waves is digitized on a mesh in the x = x₀ plane, at a large distance to the grating and lenses. Then, using this digitized wave electric field, the amplitudes and local propagation directions of the introduced system of Gaussian beams centred at the mesh points are calculated. Summing over all Gaussian beams, an expression for the electric field of X-rays in the region is constructed using the formulas (4), (5) and (18). The resulting wave electric field, which is the sum of the Gaussian beams, is digitized on a mesh in the x = x₁ plane. These values of the wave electric field on the mesh are used to calculate the amplitudes and local propagation directions of the new introduced Gaussian beams for the region x ≥ x₁. Then, the changes in the amplitudes and propagation directions of Gaussian beams due to their propagation through the optical element (lens/kinoform) or sample (copper grating) is calculated using the formulas (8) and (9). The resulting electric field of the wave, which is the sum over Gaussian beams, is obtained for x ≥ x₁. Next, the values of this electric field are calculated on the mesh in the x = x₂ plane. Therefore, when you put it all together, the wave propagation through the sample or lens is calculated using the formulas (4), (5), (18), (8) and (9), and the electric field of the wave in front of the next optical element on the mesh on the plane x = x₂ is calculated. The calculated values of the electric field on the mesh are used to calculate the amplitudes and local propagation directions of the system of new Gaussian beams for the region x ≥ x₂. Thus, an expression for the wave electric field is obtained, which is the sum over Gaussian beams for x ≥ x₂. The resulting electric wavefield for the region x ≥ x₂ is digitized on a mesh in the plane x = x₃, and so on. This is repeated until the wavefield after the last lens/kinoform is calculated. The resulting wavefield after the last lens/kinoform, which is the sum over the Gaussian beams, is used to calculate the results for |E| after all optical elements. In particular, the resulting wavefield after all lenses/kinoforms is used to determine the best focus position or to image the grating.

The three-dimensional problem is solved in a similar way. However, instead of formulas (4) and (5), the formulas (2) and (3) are used.

The propagation of X-ray waves through the sample (grating) is calculated using the same formulas as the propagation of waves through the lens, but only the coefficients are taken for copper.

5.1. Setting up a computer program

Setting up a computer program includes finding out and regulating the accuracy of the approximate solution obtained. This issue has some complications.

In general, the question of the actual calculation error is one of the most difficult in the theory of numerical calculations. Usually, for the mathematical justification of a numerical method, the asymptotic convergence of an approximate solution to the exact one is proved. However, asymptotic convergence is not a real error estimate, but it is the only possibility for achieving the required accuracy. In papers, usually, instead of evaluating the calculation error, the rate of convergence of an approximate solution to an exact one or the calculation results showing that for different values of the method parameter the calculation results are close to each other are given. Due to these objective mathematical difficulties, most papers do not contain numerical estimates of actual errors of the results obtained. Instead, for example, an approximate solution is compared with an experiment (although this is not the same thing, because the model can have drawbacks).

The calculation accuracy strongly depends on subtle details of the problem (for example, whether the equation coefficients are differentiable functions or not). When solving similar problems, but differing in such important details, the accuracy can differ significantly. Also, if one solves the same problem with different simulation methods, the accuracies of the obtained approximate solutions can differ.

Finally, the answer to the question of accuracy also depends on what we mean by accuracy. For example, Fourier methods converge on average over the area, and this convergence is often considered sufficient. Accordingly, the calculation accuracy can be estimated on average over the region. However, having high accuracy on the average over the area, there can be large deviations of the approximate solution from the exact one at and near some points. This depends on solution smoothness, and the solution smoothness depends on the equation coefficient smoothness, and on the smoothness of the field at the boundary.

Despite these theoretical difficulties, in practice it is often possible to suggest some simple asymptotic working formulas for the calculation error under conditions when the deviation of the numerical solution from the exact one is small. Thus, in the paper by Kshevetskii et al. (2016), some simple formula was proposed for the calculation error for used finite-difference methods. The formula was developed on the basis of the well known Runge's rule, and the calculation errors were evaluated. For the Gaussian beam method used in this work, simple formulas based on Runge's rule were also proposed for calculation errors (Wojda & Kshevetskii, 2019).

In Wojda & Kshevetskii (2019), the accuracy of the focusing simulation was investigated and it was shown that the maximum amplitude of the electric field in the focal spot and the half-width of the focal spot is very sensitive to the numerical simulation quality, while the focal length is easily calculated even with rough simulations.

The resulting image quality mainly depends on the focusing quality. Therefore, to ensure high simulation accuracy, we control the errors of calculating the maximum electric field amplitude in the focal spot and the focal spot width at half of the maximum amplitude as the most sensitive to the calculation quality.

As a whole, the algorithm for the accuracy control in the case of ideal lenses is as follows:

(i) First of all, one needs to find the number of Gaussian beams required for qualitative focusing. It is necessary to take so many Gaussian beams so that the error in calculating the maximum electric field amplitude at the focal spot and the width of the focal spot at half of the maximum amplitude does not exceed 1%. To estimate these errors, the method and the formula for error from Wojda & Kshevetskii (2019) is used, degree k = 2.

(ii) Further, if we calculate the grating image, then the number of Gaussian beams should not be less than in the focusing calculation. To find the required number of beams, it is appropriate to proceed as in the case of focusing calculations. Several bright points in the grating image are selected, and the number of beams is found in order that the error in the electric field amplitude at these points does not exceed 1%.

The number of Gaussian beams is found, and the calculations are tuned.

Now let us take a look at kinoform lenses. In principle, the algorithm for setting up calculations with kinoforms is the same. However, a kinoform lens has a rather smeared, wide focal spot, elongated along the main optical axis, and uneven. The maximum field amplitude in the focal spot is poorly expressed, and the general blurring of the focal spot also interferes with the adjustment (this can lead to poor image quality due to physical, not mathematical, reasons). However, an additional rule helps: the number of Gaussian beams should be no less than in a similar case for ideal lenses. In the formula for estimating the error (Wojda & Kshevetskii, 2019), it is necessary to take k = 1 because it follows from formula (11) that the error is O(h) on the lens thickness jumps.

In reality, the calculation parameters are as follows. For the planar case, the calculations were carried out for 2000 points with a space step h = 0.5 µm. For three-dimensional ones the calculation was carried out for 400 × 400 points with a space step h = 1.25 µm. The step values differ in simulations because only one lens is used in 3D modelling, while in the planar case up to 30 lenses are taken into account. The computational complexity (the number of operations) of the presented method using Gaussian beams is estimated as N_I × N_F. Here N_I is the number of points digitizing the electric field after the lenses, and N_F is the number of points at which we calculate the electric field after the lenses at the distance we are interested in. After the X-rays have propagated through a system of several tens of lenses, the wavefield becomes fast oscillating, and a rather dense numerical grid is required for its digitization, which determines the number N_I. The number N_F is determined only by the desirable resolution of the picture of interest.

5.2. Comparison of the calculation method with some other existing ones

At present, Fourier methods are often used for calculations in X-ray optics. It is known that the computational complexity of the fast Fourier transform (FFT) is estimated as Mlog₂(M), where M is the number of points used in the FFT method (Cooley et al., 1967). The FFT method is very fast. However, the FFT algorithm is somewhat lacking in robustness and stability (Chubar et al., 2017). Standard methods for estimating the error of Fourier series truncating are expressed in terms of derivatives of the decomposed function. In X-ray optics, we deal with functions digitized on a mesh on a plane, and estimating the truncating error may be more difficult.

Usually, when the Fourier method in X-ray optics is applied, the paraxial equation is solved,

where E = exp(ikx)A. When calculating wave propagation in air, a simpler equation with constant coefficients is solved,

whose solution can be written analytically in Fourier series, due to which one can apply FFT. At the same time, the effect of lenses on wave propagation is calculated based on the simplified equation

It is easy to find that the solution of equation (21) in the case of kinoforms contains discontinuities in the variables y and z. Thus, the boundary condition for equation (20) will contain discontinuities. It is known from the Fourier series theory that, in the case of a function with discontinuities of the first kind, the terms of the Fourier series decrease slowly, as 1/k_m, where k_m is the wavenumber of the harmonic m. It is also known that the formulas for the truncated series error for this case are poorly developed, since this error depends on subtle details of the function's behaviour, and it is difficult to write a general effective formula. Thus, although the Fourier method converges in the case of kinoforms, and the solution can in principle be obtained, nevertheless, in the case under consideration, the method converges slowly and also there may be difficulties in controlling the approximate solution error.

Thus, in the case of kinoforms, the described method of oriented Gaussian beams may have some advantage, because the influence of discontinuity points on the solution is easily controlled. Only a relatively small number of Gaussian beams effectively intersect with discontinuities, while the image is created by all Gaussian beams together. With a large number of Gaussian beams, possible inaccuracies in the parameters of Gaussian beams intersecting with discontinuities are not of great importance, because there are few such beams relative to their total number. Gaussian beams intersecting with discontinuities can even be ignored, which is shown when the formula (11) is constructed and investigated.

Moreover, in the case of kinoform lenses, the mathematical interpretation of equation (1) causes certain difficulties due to the presence of a jump change in the coefficient. When the wave falls on a discontinuity approximately perpendicular to the discontinuity surface, we use the condition of the electric field continuity at the discontinuity. When the angle between the direction of wave propagation and the discontinuity surface is small, then – physically – diffraction effects should arise. At the same time, some clear principles should be implemented mathematically. We can use an approximate solution in the form of a sum of Gaussian beams to define the solution for this case. If we discard in the sum the Gaussian beams effectively intersecting with discontinuities and go to the limit of an infinite number of Gaussian beams, then we obtain a physically justified definition of the solution for this case. Expressions for Gaussian beams are specified more exactly by Maslov (1994).

To solve the X-ray wave propagation problem, Frechet or Kirchhoff integral representations can be used in the paraxial approximation. The integrands of these integrals are products of the electric field function on the surface and a very fast oscillating function (in the paraxial case, this is the Kirchhoff propagator). In simple cases, these integrals can be calculated analytically. Due to fast-oscillating functions in the integrands, these integrals often cannot be calculated without using special numerical methods. Using the theorem that the integral of the product of functions is equal to the integral of the product of the Fourier transforms of the functions is a popular way to calculate these integrals numerically (Chubar & Elleaume, 2013). The use of this theorem essentially returns the question of these integrals to the question of the Fourier series. The Fourier series has to be truncated, and control of the truncated error is desirable. This issue has already been discussed above. In Kshevetskii & Wojda (2015), an alternative way of calculating the Kirchhoff integral is proposed, which does not use the Fourier transform but is a version of Philon's well known numerical method for calculating the integrals of fast-oscillating functions (Levin, 1982; Iserles, 2004).

Sums (3) and (5) over Gaussian beams may be considered as integral sums, which in the limit h → 0 give integral representations for a solution like Frechet or Kirchhoff integrals. However, a detailed investigation of this issue is beyond the scope of this work.

The development of a theoretical apparatus for solving various problems of X-ray wave propagation and investigation of the behaviour of these waves on the medium parameter jumps is the purpose of this study. Over the course of the research, a computer program for solving the problems of interest has been developed. The computer program is a research tool but was not the purpose of this research, and therefore, according to the authors, it is not the main result of this work.

5.3. Results of calculating the X-ray wave propagation through kinoforms and ordinary lenses

The theory of calculating the X-ray wave propagation based on the Helmholtz equation using the Gaussian beam technique is described by Wojda & Kshevetskii (2019). Since this work has some specifics, for example we consider kinoforms and not ordinary lenses, the proposed methods require some development and generalization, and this is done in this work.

A stepwise change in lens thickness results in diffraction effects caused by this step-change in thickness. The study of these diffraction effects is interesting from a physical point of view. But we also draw attention to some theoretical problems. Derivatives in (9) do not exist where a lens thickness jump occurs. Therefore, formulas (9) require modification or some special justification. Perhaps some development of the Helmholtz equation theory is required to make the theory applicable to the considered non-standard case.

Diffraction effects arising from jumps in the kinoform-lens thickness are more pronounced in the case of planar lenses, as compared with three-dimensional kinoform lenses. Therefore, planar lenses are the most interesting for theoretical research, and we will pay the greatest attention to such lenses.

When an X-ray wave propagates through a biconcave lens, we can observe diffraction caused only by the edge of the lens aperture (Wojda & Kshevetskii, 2019), and the effect is not very strong. If an X-ray wave propagates through a kinoform lens, then spherical waves appear in places where there is a sharp change in the lens thickness. They appear as spots and are a consequence of X-ray diffraction. Fig. 4 shows the result of calculating the X-ray wave propagation through the kinoform lens; the figure shows the detailed structure of the electric field module at a distance of 1 m from the optical object. Calculations are performed for the X-ray wavelength λ₀ = 10⁻¹⁰ m, β = 3.1801 × 10⁻¹⁰ and δ = 2.2156 × 10⁻⁶. For comparison, in Fig. 4, there is a red dashed line that shows the results of calculations for a conventional single beryllium lens. The lens has a curvature of R = 20 000/m. The kinoform lens has the same surface curvature as the conventional lens; the calculations for the kinoform lens are shown with the blue line.

Figure 4 Calculation results for |E| at distances of 1, 3, 20 and 40 m after one conventional lens (red dashed line) and one kinoform lens (blue line).

The performed calculations show that a kinoform lens does not focus X-rays as well as a perfect lens does. Fig. 5 demonstrates this. The reason a kinoform lens is not as efficient at focusing radiation as a perfect lens is because a step change in the kinoform width causes some of the radiation to scatter. The distance at which this scattering occurs (less than 1 m from the kinoform) is much shorter than the distance at which the lens focuses the X-ray waves (11 m).

Figure 5 Results for |E| before the focus spot (at the distance 11 m after one lens) obtained with 1 kinoform lens (blue line) and 1 ideal lens (red dashed line).

The cause of X-ray radiation scattering is an abrupt change in the wave electric field phase due to a jump-like change in the kinoform lens thickness. From these breakpoints, X-ray waves propagate in the form of spherical waves. At the distance of 20 m and 40 m, we can observe a strong influence of stepwise changes in the thickness of the kinoform lens on the shape of the graph of the electric field intensity; see plots in Fig. 4.

However, calculations show that by increasing the number of kinoform lenses used, the X-ray scattering effect can be significantly reduced. Focusing X-rays on a large number of kinoform lenses is significantly more efficient because more lenses focus the radiation at a shorter distance.

To take into account in detail the challenging shape of the kinoform lens, when calculating the propagation of X-rays through the lens one should use Gaussian beams with a very small width. The narrower the Gaussian beams used, the smaller the error in the obtained approximate solution arising from inaccurate accounting for the stepwise change in the kinoform lens thickness. The formulas for the electric wavefield contain many Gaussian beams. The smaller the widths of the Gaussian beams used, the more Gaussian beams are used and the less the influence of each Gaussian beam on the final wavefield. Therefore, jumps in the thickness of kinoform lenses have a limited effect on the wavefield, and the influence of jumps on the wavefield may be better calculated. Below we will consider in more detail the influence of jumps in the lens thickness on the wave electric field.

Given the calculation results presented, it would appear that kinoform lenses are not suitable for focusing X-rays. Let us see, however, if increasing the number of kinoform lenses used will improve the focusing of X-rays.

The graphs in Fig. 6 show a comparison of the absolute value of the electric field obtained with 10 ideal lenses and 10 kinoform lenses, respectively. The main differences are as follows. The maximum electric field value with ideal lenses was obtained at a distance of 1.126 m, and the highest value with kinoform lenses was obtained at a distance of 1.133 m. Also, when using kinoform lenses, the electric field intensity at the best focus distance is not exactly concise with the Gaussian law. However, additional peaks are noticeable, the size of which is only several times smaller than the maximum value.

Figure 6 Focus spot obtained with 10 kinoform lenses (blue line) and 10 ideal lenses (red dashed line).

We should note that an increase in the number of kinoform lenses used significantly improves the quality of focusing, and, with a significant number of kinoform lenses, focusing is practically not inferior to that of conventional ideal lenses. Fig. 7 shows for comparison the maximum value of the electric wavefield, which was obtained with 30 ideal lenses (red) and 30 kinoform lenses (blue). One can see that using 30 kinoform lenses allows very good focusing of X-ray waves.

Figure 7 Focus spot obtained using 30 perfect lenses (red dashed line) and 30 kinoforms (blue line).

Optical aberration occurs in an optical system that uses a kinoform lens. Increasing the number of kinoform lenses allows us to reduce the defects of the optical system, and also enables better focusing of X-rays.

Although the optical system of many kinoform lenses is good for focusing X-rays, the calculations performed show that kinoform lenses are probably not suitable for obtaining high-quality magnified images of objects. Calculations show that a stepwise change in the thickness of a kinoform lens causes significant deviations of the X-ray radiation intensity from the Gaussian shape both before the focusing point and after the focusing point. We can see this in the graphs in Fig. 8. The diffraction effects visible in the graphs in Fig. 8 will undoubtedly manifest themselves during attempts to obtain images of objects.

Figure 8 Results of calculation of electric field obtained for 30 kinoform lenses (blue line) and 30 ideal lenses (red dashed line) 0.166 m after optical objects and 1.5 m after optical objects.

Figs. 9 and 10 show the focal spot of one 3D kinoform lens. Diffraction effects due to jumps in the kinoform lens thickness in the 3D case are weaker than in the planar case, but these effects certainly take place and they are significant.

Figure 9 Distribution of the electric field modulus in a focal spot created by one 3D kinoform lens (size in millimetres).

Figure 10 Distribution of the electric field modulus in a focal spot created by one 3D kinoform lens (size in millimetres). The field magnitude is shown with the density plot.

Figs. 11, 12 and 13 show grating images obtained from the comparison between one kinoform lens and one ideal lens. The calculations were performed for diamond lenses under the assumption that the X-ray wavelength is equal to 1 Å. For diamond and an X-ray wavelength of 1 Å, the complex refractive index values are: δ = 4.7423 × 10⁻⁶ and β = 2.8366 × 10⁻⁹. The distances between the grating and the optical system were different in the calculations to obtain 1, 2, 4 and 8 magnification factors. The distances between the optical system and the image of the grating behind the optical system were determined using the well known thin lens formula

where the focal length f is found from the formula f = 1/(2δR). Thus, in all figures, these are the clearest images possible within the framework of the optical systems used.

Figure 11 Image of the grating directly after the grating.

Figure 12 Gating image at a distance of two, three, five and nine focus distances from the kinoform lens. The field magnitude is shown with the density plot.

Figure 13 Grating image at a distance of two, three, five and nine focus distances from the ideal lens. The field magnitude is shown with the density plot.

6. Some generalizations: formulas for lenses with defects

Below we generalize the previous formulas for the case when lenses contain internal defects (cavities, inclusions) or external defects (surface irregularities, surface roughness). This item, in the physical sense, is close to kinoform lenses since the main thing that distinguishes both cases is a possible non-smooth, abrupt behaviour of the phase of an X-ray wave passing through such lenses.

Consider formulas (9): these contain partial derivatives of functions describing lens thickness. Differentiable functions are a mathematical abstraction, so the differentiability of functions describing real physical objects cannot be verified experimentally. Obviously, formulas (8) and (9) are for ideal­ized lenses, whose thickness may contain some abrupt changes. On the whole, the thickness is described with smooth functions, that need some adjustments to make them suitable for cases of non-smooth behaviour of the parameters.

Below are some modifications of formulas (8) and (9) that allow taking into account the irregularities in the lens shape and even internal lens defects,

where

In (25), x₁ and x₂ are the coordinates of the beginning and end of the lens plate. Based on physical intuition, differentiability of the function Φ(y, z) with respect to the variables y, z is not assumed. [The case of a piecewise-differentiable function Φ(y, z), as a simpler one, has already been considered above.] It is assumed that the function Φ(y, z) is piecewise continuous. The assumption that the function Φ(y, z) is piecewise continuous is not burdensome and seems quite physically justified. Recall that the differentiability of a function implies its continuity, but the converse is generally speaking not true. [Weierstrass gave a wonderful example of a function that is continuous at each point, but not differentiable at any point. This is the function (Weierstrass, 1894).]

The Gaussian beam technique, developed in Wojda & Kshevetskii (2019) and in this paper, assumes that at each step of the procedure used we start from the field E(x₀, y, z), given on some plane x = x₀, and the function E(x₀, y, z) is piecewise differentiable with respect to the variables y and z, or smoother. Under this condition, the approximate solution, constructed within the framework of the Gaussian beam technique, converges to the exact one as . Accordingly, within the framework of this method, the function Φ(y, z), as it relates to the function E(x₀, y, z), has to be at least piecewise-differentiable. Therefore, the need exists for removing the existing limitations of the method to allow similarly piecewise continuous functions. Below, such a fundamental possibility will be shown.

Let us consider the function

This function Ψ(y, z, a) is differentiable with respect to both variables y and z, and as at all points of continuity of the function Φ(y, z). Let us replace the function Φ(y, z) with the function Ψ(y, z, a) everywhere in formulas (23) and (24), considering the obtained expressions as approximate solutions, while the exact solution is given by the limit . This then implies that the complete solution of the problem, for the electromagnetic wavefield, is . This is a general scheme for obtaining the solution for the considered case of non-smooth function n(x, y, z). Then from (26) the approximate relations follow,

In the formulas (27) and (24), symbols of partial derivatives are understood not as ordinary derivatives but as generalized derivatives, or derivatives in a weak sense (Ladyzenskaya, 1985). The solution obtained in the limit is understood as a generalized or weak solution of the equation (1) under study (Ladyzenskaya, 1985).

The formulas (27) are similar to the formulas (10). In practical calculations, the choice of the parameter a is important. When deriving formulas (8) and (9), it was assumed that the beams are narrow and that the incident phase can locally be approximated with a linear function. So it makes sense to use a > σ.

A few words about the importance of studying the effect of small-scale defects in lenses on focusing and imaging. In this paper, kinoform lenses, which roughly correspond to real lenses of this type, are considered. However, we can also imagine an extremely thin kinoform lens containing a very large number of concentric rings. Such a lens, although extremely thin, is likely to be capable of focusing X-rays, albeit with poor quality. We can think of such a very thin kinoform lens as a model of the influence of lens defects on focusing and imaging. The effect may be significant. Therefore, the question requires research.

7. Conclusion

The article contains a description of a new mathematical apparatus for high-precision calculation of X-ray wave propagation for cases when the medium is highly inhomogeneous or the optical element shape is sharply changing. The paper examines the propagation of X-ray waves through kinoform lenses and analyzes the use of these lenses for focusing and imaging. The thickness jumps in kinoform lenses produce diffraction effects that merit interest and are difficult to calculate, making them a worthwhile field of study. In the research, preference is given to the study of planar lenses, since in this case the diffraction effects are greater than in three-dimensional lenses. These diffraction effects manifest themselves in the form of the electric field oscillations of significant amplitudes, but at very fine spatial scales and with a very fine structure. These fine diffraction structures disappear as the distance from the lens increases, but the optical system properties change.

To calculate the X-ray wave propagation, the method of oriented Gaussian beams is developed and applied. Within the framework of this approach, the X-ray wave is represented as a superposition of oriented Gaussian beams. Explicit formulas are written describing the Gaussian beam propagation through a lens or a sample.

Let us first consider the results obtained in calculations for planar lenses. Lens thickness jumps lead to diffraction effects when X-ray waves propagate through the lens. As a result, on the plane at the focal length, we can observe one main peak of the X-ray wave intensity and peaks of smaller but significant amplitudes. These satellite peaks prevent the formation of quality images with kinoform lenses.

However, when many kinoform lenses are used, the number of peaks on the plane of the focal length decreases, and the value of the main peak increases. Thus, the focal spot obtained with kinoform lenses is not very different from that obtained with ideal lenses. Therefore, the use of a large number of kinoforms can have a positive effect on imaging. Nevertheless, although the focal spot is quite good, diffraction effects due to lens thickness jumps are observed at distances both shorter and longer than the focal length. This effect should degrade imaging quality.

Compared with diffraction effects on thickness jumps of planar lenses, 3D kinoform lenses generate wavy diffraction effects. These effects have a lower amplitude than in the case of planar lenses, but their characteristics (number of waves) are the same (compare Fig. 5 and Fig. 9). Hence, we can conclude that the study of diffraction effects at thickness jumps of planar lenses gives a good notion of this kind of diffraction phenomena.

To illustrate the capabilities of the Gaussian beam technique, the results of calculating the image of a grating placed in front of either an ideal lens or a kinoform lens are presented. The image is presented at various distances behind the lens, according to the thin lens formula. One can see that in the case of kinoform lenses the image has many defects.

The use of kinoform lenses for focusing X-rays was previously experimentally investigated by Gorelick et al. (2019), and it has been argued that kinoform lenses may have limitations in focusing X-rays. In this study, we showed that kinoform lenses can also have imaging limitations. However, kinoform lenses have been considered here rather for studying diffraction effects and with the aim of developing theoretical methods for their calculation, because one can expect similar diffraction effects on defects (cavities, inclusions) of real lenses and when X-rays propagate through porous materials. Thus, the development of theoretical methods for calculating these diffraction effects in X-rays optics may be useful.