research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775

Oriented Gaussian beams for high-accuracy computation with accuracy control of X-ray propagation through a multi-lens system

aGdansk University of Technology, Faculty of Applied Physics and Mathematics, Gdansk, Poland, bCenter for Functionalized Magnetic Materials (FunMagMa), Immanuel Kant Baltic Federal University, Kaliningrad, Russia, and cImmanuel Kant Baltic Federal University, Kaliningrad, Russia
*Correspondence e-mail: pawwojda@pg.edu.pl

Edited by G. Grübel, HASYLAB at DESY, Germany (Received 18 July 2018; accepted 7 December 2018; online 21 January 2019)

A highly accurate method for calculating X-ray propagation is developed. Within this approach, the propagating wave is represented as a superposition of oriented Gaussian beams. The direction of wave propagation in each Gaussian beam agrees with the local direction of propagation of the X-ray wavefront. When calculating the propagation of X-ray waves through lenses, the thin lens approximation is applied. In this approximation, the wave parameters change discontinuously when the wave passes through a lens; the corresponding explicit formulae are derived. The theory is applied to highly accurate calculation of the focusing of X-rays by a system of many beryllium lenses. Fine structure of the wave electric field on the focal plane is revealed and studied. The fine structure is formed due to the diffraction of waves at the edges of the lens apertures. Tools for controlling the calculation accuracy are proposed. The amplitude of the electric field on the focal plane and the focal spot width are shown to be very sensitive to the quality of the calculation, while the best focus position can be obtained even from simple calculations.

1. Introduction

At present, X-ray optics are developing rapidly due to the possibilities of visualization of very small objects that they provide. To visualize the structure of a very small object, it is first necessary to focus the X-rays. Currently, commonly used focusing elements are bent crystals, multilayers or mirrors, capillaries, waveguides, refractive lenses and diffraction elements, such as Fresnel zone plates or Laue multilayer lenses (Snigirev et al., 1996[Snigirev, A., Kohn, V., Snigireva, I. & Lengeler, B. (1996). Nature, 384, 49-51.]; Roth et al., 2017[Roth, T., Alianelli, L., Lengeler, D., Snigirev, A. & Seiboth, F. (2017). MRS Bull. 42, 430-436.]).

In this paper, focusing in X-ray optics using X-ray refraction mechanisms will be considered. The first successes in this field were presented by Snigirev et al. in 1996, where a composite refractive lens was used (Snigirev et al., 1996[Snigirev, A., Kohn, V., Snigireva, I. & Lengeler, B. (1996). Nature, 384, 49-51.]). Since then, many experiments have been performed using focused X-rays. The use of lenses with a spherical or parabolic concave profile is one of the most popular methods of X-ray focusing nowadays. Recently, researchers have also used diamond lenses for X-ray focusing (Roth et al., 2017[Roth, T., Alianelli, L., Lengeler, D., Snigirev, A. & Seiboth, F. (2017). MRS Bull. 42, 430-436.]).

The materials used for X-ray focusing have some common features. First of all, they are characterized by a complex refractive index n = 1 − δ + iβ, where δ is responsible for the refraction of X-rays in lenses and β for the attenuation of the X-ray intensity. The ideal material for making X-ray lenses should have a δ value as large as possible and an absorption coefficient β as small as possible. In addition, the material should allow the technological possibility of manufacturing lenses with a small curvature radius R. Knowing the material properties, it is possible to calculate with high accuracy the focal length f of an optical system consisting of N lenses with a given surface curvature R using the formula f = [{R}/({N\delta})]. In addition to the listed properties, the materials must be strong and inexpensive. Unfortunately, it is not easy to satisfy the latter two conditions. Beryllium is one of the most popular materials currently used for focusing X-rays; the cost of manufacturing a beryllium lens is about 5000 euros. Unfortunately, while the lens is in use the lens material is subject to oxidation, and after a while the lens has to be replaced. In order to effectively focus X-rays at relatively small distances and for achieving a large magnification, use of complex optical systems consisting of several dozen lenses is necessary. Therefore, the experiments are associated with significant costs. A way to reduce this cost would be to find or develop cheaper and more effective materials for lens manufacturing. Performing more extensive theoretical studies based on mathematical modeling is another way. This work is devoted to the development of a more advanced mathematical apparatus for theoretical studies and the application of advanced methods for calculating the X-ray focusing by a system of many lenses.

At present, in theoretical studies in X-ray optics, the most popular way to solve the paraxial wave equation is by fast Fourier transform (FFT) (Goodman, 1996[Goodman, J. W. (1996). Introduction to Fourier Optics, 2nd ed. New York: McGraw-Hill.]). FFT is a very effective calculation method. However, the exact solution of the paraxial equation has the form of an infinite Fourier series, and this series has to be truncated when we move to the FFT, which leads to error. The magnitude of the error obviously depends on the circumstances, but the issue of accuracy and reliability of calculation has been poorly investigated. The methods for estimating the Fourier series truncating error are well developed for functions specified analytically, and the truncation error is expressed in terms of derivatives of the decomposed function. In calculations in X-ray optics, we deal with functions digitized on a mesh, and classical methods of estimating the Fourier series truncation error are difficult to apply to functions specified in such a form. The Nyquist principle considers the necessary digitizing frequency of a function, but this is not an estimate of the Fourier series truncating error. In addition, studies show that various focal spot characteristics have significantly different sensitivities to the quality of the calculation (Kshevetskii et al., 2016[Kshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305-1314.]). Therefore, considering the accuracy of calculating specific parameters instead of the calculation accuracy on the whole might be more appropriate. The function describing the electric field of the X-ray wave passing through a system of several tens of lenses becomes rapidly oscillating and requires use of a very dense numerical mesh for its digitization and many Fourier series terms in calculations. In addition, it has been proved that, for a larger number of lenses, the faster the oscillating wavefield after the lenses, and the more Fourier series terms have to be taken into account for an acceptable description (Kshevetskii et al., 2016[Kshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305-1314.]). These circumstances show the limited possibilities of Fourier methods and stimulate the development of new methods for solving the problems of X-ray optics. New methods should provide high accuracy of calculations in complex cases and should be equipped with tools for controlling the accuracy and reliability of calculations (Kshevetskii et al., 2016[Kshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305-1314.]; Kshevetskii & Wojda, 2015[Kshevetskii, S. P. & Wojda, P. (2015). Math. Appl. 43, 193-206.]).

In this paper, a method for calculating the propagation and focusing of X-ray waves based on the use of oriented Gaussian beams is presented. The idea of using Gaussian beams in optics is not new. The technique of using beams in optics has been given elsewhere (Kogelnik, 1965[Kogelnik, H. (1965). Appl. Opt. 4, 1562-1569.]; Keller & Streifer, 1971[Keller, J. B. & Streifer, W. (1971). J. Opt. Soc. Am. 61, 40-43.]; Deschamps, 1972[Deschamps, G. (1972). Proc. IEEE, 60, 1022-1035.]; Deschamps et al., 1983[Deschamps, J., Courjon, D. & Bulabois, J. (1983). J. Opt. Soc. Am. 73, 256.]). In later works (Chabory et al., 2005[Chabory, A., Sokoloff, J., Bolioli, S. & Combes, P. F. (2005). C. R. Phys. 6, 654-662.], 2010[Chabory, A., Sokoloff, J., Bolioli, S. & Elis, K. (2010). Proceedings of the 4th European Conference on Antennas and Propagation (EuCAP'2010), 12-16 April 2010, Barcelona, Spain, pp. 1-5.]; Ghannoum et al., 2009[Ghannoum, I., Letrou, C. & Beauquet, G. (2009). Proceedings of RADAR 2009 - International Radar Conference `Surveillance for a Safer World', 12-16 October, 2009, Bordeaux, France.]), Gaussian beams (a Gaussian shooting algorithm) were used in radio-location for finding the location of an object. X-ray waves are characterized by an exceptionally high wave frequency; therefore, the method using Gaussian beams needs to be adapted to solve the problems of X-ray optics.

Here, we set the incident wave in the form of a simple coherent wavefront, and we do not touch on the problem of real sources of X-ray waves, partially coherent and non-strictly monochromatic. However, Gaussian beams can also be used to study partially coherent radiation beams, as exemplified by the Gaussian Shell Model (Deschamps, 1972[Deschamps, G. (1972). Proc. IEEE, 60, 1022-1035.]; Deschamps et al., 1983[Deschamps, J., Courjon, D. & Bulabois, J. (1983). J. Opt. Soc. Am. 73, 256.]). There exists a numerical method based on the analytical stationary phase method (Bahrdt, 2007[Bahrdt, J. (2007). Phys. Rev. ST Accel. Beams, 10, 060701.]), which resembles to some extent the method presented in this paper. It does not use Gaussian beams but allows simulation of electric field propagation `along geometrical rays'.

In X-ray optics, the paraxial wave equation (Goodman, 1996[Goodman, J. W. (1996). Introduction to Fourier Optics, 2nd ed. New York: McGraw-Hill.]; Kohn, 2012[Kohn, V. G. (2012). J. Synchrotron Rad. 19, 84-92.]) is now widely used for calculating the electric wavefield. The paraxial equation is approximately derived from the Helmholtz equation (Levy, 2000[Levy, M. (2000). Parabolic Equation Methods for Electromagnetic Wave Propagation. Londo: The Institution of Electrical Engineers.]; Ishimaru, 1991[Ishimaru, A. (1991). Electromagnetic Wave Propagation, Radiation and Scattering. Englewood Cliffs: Prentice Hall.]; Babich & Buldyrev, 1991a[Babich, V. M. & Buldyrev, V. S. (1991a). Asymptotic Methods in Short-Wavelength Diffraction Theory. Oxford: Alpha Science International.],b[Babich, V. M. & Buldyrev, V. S. (1991b). Short-Wavelength Diffraction Theory: Asymptotic Methods. Berlin: Springer Verlag.]). Gaussian beams are exact solutions of the paraxial equation. Since the widths of the beams are much larger than the X-ray wavelength, Gaussian beams can also be treated as almost exact solutions of the Helmholtz equation. As we shall see, practically all solutions of the problems of X-ray optics can be constructed in the form of specially selected sums of Gaussian beams. We use oriented Gaussian beams, that is, Gaussian beams propagating at certain angles to the optical axis of the optical system. The condition that the propagation directions of the beams have to be orthogonal to the phase front of the simulated wave determines the inclination angles.

The developed approach combines properties of geometric optics (the Gaussian beams have clear directions of wave propagation and are localized in transverse directions) and of wave optics (the Gaussian beams are almost exact solutions of the Helmholtz wave equation). The method is well suited for solving problems with many lenses, and allows us to highly accurately calculate not only the focusing but also the fine-structure effects on the focal plane that result from X-ray diffraction at the edges of lens apertures.

Also, the paper proposes tools for controlling the calculation accuracy and shows the possibility of calculating various focusing parameters with the required accuracy. It is revealed that the half-width of the focal spot and, in particular, the amplitude of the electric field on the focal plane are very sensitive to the calculation quality. The position of the best focus can be obtained even with simple calculations.

2. Basic equations and starting points

The propagation of a monochromatic electromagnetic wave with frequency ω0 in a medium with complex refractive index n = 1 − δ + iβ is described by the Helmholtz equation (Levy, 2000[Levy, M. (2000). Parabolic Equation Methods for Electromagnetic Wave Propagation. Londo: The Institution of Electrical Engineers.]; Ishimaru, 1991[Ishimaru, A. (1991). Electromagnetic Wave Propagation, Radiation and Scattering. Englewood Cliffs: Prentice Hall.]; Babich & Buldyrev, 1991a[Babich, V. M. & Buldyrev, V. S. (1991a). Asymptotic Methods in Short-Wavelength Diffraction Theory. Oxford: Alpha Science International.],b[Babich, V. M. & Buldyrev, V. S. (1991b). Short-Wavelength Diffraction Theory: Asymptotic Methods. Berlin: Springer Verlag.]),

[\Delta E({\bf{r}})+k_{0}^{2}\,n^{2}E({\bf{r}})=0,\qquad k_{0} = {{\omega_{0}} / {c}}, \eqno(1)]

where E is the electric field, r = (x, y, z), k0 is the wavenumber, c is the speed of light, β is an absorption coefficient and δ is a refractive decrement (β and δ are non-negative).

More convenient for practical applications, a simplified wave equation can be derived from the Helmholtz equation (1)[link] (Leontovich, 1944[Leontovich, M. A. (1944). Izv. Akad. Nauk. SSSR Ser. Phys. 8, 16.]). Let us consider the case when the wave propagates along the x-axis, and the characteristic scales ly, lz of the wave along the axes y and z are much larger than the characteristic scale lx of the wave along the x-axis: [l_{x} \ll l_{y}], [l_{x} \ll l_{z}]. In this case, equation (1)[link] can be greatly simplified. If we substitute

[E\left({\bf{r}}\right) = A\left({\bf{r}}\right)\exp\left(ik_{0}x\right), \eqno(2)]

into equation (1)[link] we obtain the equation for the function A(r). We suppose that A(r) varies only slowly with variable r. Neglecting the small term of a second derivative of the function A(r) with respect to the x-variable, we arrive at the paraxial equation

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

Equation (3)[link] was first proposed by Leontovich in 1944 (Leontovich, 1944[Leontovich, M. A. (1944). Izv. Akad. Nauk. SSSR Ser. Phys. 8, 16.]) for describing the propagation of a monochromatic electromagnetic wave within the paraxial approach. Equation (3)[link] is often called the paraxial equation.

Equations (1)[link] and (3)[link] describe the propagation of monochromatic waves with frequency ω0. Considering quasimonochromatic waves characterized by a set of frequencies ω close to ω0, we can construct the solution of such a problem in the form of a superposition of waves with different ω. Thus, a solution of the general problem of propagation of nonmonochromatic waves still leads to equations (1)[link] or (3)[link].

3. The propagation of X-rays in a vacuum

Let us consider equation (3)[link] for the case of a vacuum, that is, for complex refractive index n = 1, then equation (3)[link] takes the following form,

[{{ {{\partial}} A} \over { {{\partial}} x}} = {{i} \over {2k_{0}}}\left({{ {{\partial}} ^{2}A} \over { {{\partial}} y^{2}}}+{{ {{\partial}} ^{2}A} \over { {{\partial}} z^{2}}}\right). \eqno(4)]

This equation has important applications in optics in general (Goodman, 1996[Goodman, J. W. (1996). Introduction to Fourier Optics, 2nd ed. New York: McGraw-Hill.]; Levy, 2000[Levy, M. (2000). Parabolic Equation Methods for Electromagnetic Wave Propagation. Londo: The Institution of Electrical Engineers.]), where it describes the propagation of electromagnetic waves, and in acoustics (Babich & Buldyrev, 1991a[Babich, V. M. & Buldyrev, V. S. (1991a). Asymptotic Methods in Short-Wavelength Diffraction Theory. Oxford: Alpha Science International.],b[Babich, V. M. & Buldyrev, V. S. (1991b). Short-Wavelength Diffraction Theory: Asymptotic Methods. Berlin: Springer Verlag.]). Equation (4)[link] has an exact partial solution in the form of a Gaussian beam (Kogelnik, 1965[Kogelnik, H. (1965). Appl. Opt. 4, 1562-1569.]),

[\eqalignno{ H(x,y,z,x_{0},y_{0},z_{0}) = {}& {{1} \over {1+i\,\left(x-x_{0}\right)/(k_{0}\sigma^{2})}} \cr& \times \exp\left\{-{{(y-y_{0})^{2}+(z-z_{0})^{2}} \over {2\left[\sigma^{2}+i\,\left(x-x_{0}\right)/k_{0}\right]}}\right\}. &(5)}]

Here σ is a real parameter, which controls the Gaussian beam width. At x = x0, the solution (5)[link] takes the form of a usual Gauss function,

[H(x_{0},y,z,x_{0},y_{0},z_{0}) = \exp\left\{ -{{1}\over{2\sigma^{2}}} \left[\left(y-y_{0}\right)^{2}\,+\,\,\left(z-z_{0}\right)^{2}\right]\right\}. \eqno(6)]

We can consider (6)[link] as a boundary condition given at the plane x = 0 to the equation (4)[link]. In this case, the formula (5)[link] gives a solution to equation (4)[link] that meets the boundary condition (6)[link]. Numerical examples of solutions of equation (4)[link] with other boundary conditions can be found elsewhere (Kshevetskii et al., 2016[Kshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305-1314.]; Kshevetskii & Wojda, 2015[Kshevetskii, S. P. & Wojda, P. (2015). Math. Appl. 43, 193-206.]).

Applying the formula (5)[link] to the formula (2)[link], we obtain an approximate solution

[\eqalignno{ G\left({\bf{r}},x_{0},y_{0},z_{0}\right) = {}& \exp\left[ik_{0}\left(x-x_{0}\right)\right] {{1}\over{1+i\,\left(x-x_{0}\right)/\left(k_{0}\sigma^{2}\right)}} \cr& \times \exp\left\{-{{\left(y-y_{0}\right)^{2}\,+\,\,\left(z-z_{0}\right)^{2}} \over {2\left[\sigma^{2}+i\,\left(x-x_{0}\right)/k_{0}\right]}}\right\} &(7)}]

of the Helmholtz equation (1)[link]. In our investigation, the condition (σk0)2 ≃ 109 ≫ 1 is satisfied, and this means that the function (7)[link] satisfies the Helmholtz equation (1)[link] with very high accuracy; this is an almost exact solution of equation (1)[link]. Therefore, we shall often name (7)[link] as a solution of the Helmholtz equation, for convenience's sake. We also shall call for convenience the solution (7)[link] of the Helmholtz equation a Gaussian beam. The solution (7)[link] of the Helmholtz equation (1)[link] describes the wave propagation along the OX-axis; we assume that this axis coincides with the optical axis of the optical system.

Waves in X-ray optics can propagate at small angles to the OX-axis. Therefore, it is desirable to generalize the particular solution (7)[link] of the Helmholtz equation (1)[link] to the case when the wave propagates at an arbitrary angle to the OX-axis. Let the vector [{\bf{e}}_{1}] indicate the direction of wave propagation and let the vector [{\bf{j}}] be directed along the OZ-axis. We denote [{\bf{e}}_{2}] = [{\bf{j}} \times {\bf{e}}_{1}], [{\bf{e}}_{3}] = [{\bf{e}}_{1} \times {\bf{e}}_{2}]. Then the generalization of the formula (7)[link] to the case when the wave propagates along the direction indicated by the vector [{\bf{e}}_{1}] is written as follows,

[\eqalignno{ G\left({\bf{r}},{\bf{r}}_{0},{\bf{e}}_{1}\right) = {}& \exp\left[i\,{\bf{k}}\left({\bf{r}}-{\bf{r}}_{0}\right)\right] {{ 1 }\over{ 1+i\,\left[\left({\bf{r}}-{\bf{r}}_{0}\right)\,{\bf{e}}_{1}\right]/\left(|{\bf{k}}|\sigma^{2}\right)}} \cr& \times \exp\left(-{{\left[\left({\bf{r}}-{\bf{r}}_{0}\right)\,{\bf{e}}_{2}\right]^{2}\,+\,\,\left[\left({\bf{r}}-{\bf{r}}_{0}\right)\,{\bf{e}}_{3}\right]^{2}} \over {2\left\{\sigma^{2}+i\,\left[\left({\bf{r}}-{\bf{r}}_{0}\right)\,{\bf{e}}_{1}\right]/|{\bf{k}}|\right\}}}\right). &(8)}]

Here r0 = (x0, y0, z0), [{\bf{k}}] = [k_ {0}{\bf{e}}_{1}]. Since we are only interested in Gaussian beams propagating at small angles to the OX-axis, then with high accuracy kxk0. In the general case, however, ky ≠ 0 and kz ≠ 0. As before, σ is the parameter that determines the Gaussian beam width.

Since the formula (8)[link] yields an almost exact particular solution of the Helmholtz equation (1)[link] in the case (k0σ)2 ≫ 1, there is a temptation to construct a general solution of equation (8)[link] of the problem of X-ray propagation in the form of a superposition of particular solutions of the kind (8)[link], with appropriately chosen parameters.

We introduce in space a system of lines given by the equations (y = ym, z = zl), where ym and zl are real numbers and ym+1 = ym + h, zl+1 = zl + h, and h is the step of the mesh κ. Also, we denote each of the introduced lines by the pair (ym, zl).

Suppose we know the digitized electric wave field E on the plane x = x0. That is, on the plane x = x0, we have the mesh κ of points (ym, zl) where we know the values of the field E. We denote E0ml = E(x0, ym, zl). The problem is to find out the field E in the whole space. In principle, at each mesh point (ym, zl) we can specify its own direction [{\bf{e}}_{1,ml}] of the Gaussian beam propagation determined from physical considerations. The general solution of the X-ray wave propagation problem

[E_{\rm{approx}}\left({\bf{r}}\right) = \sum\limits_{m,l} C_{ml}G\left({\bf{r}},{\bf{r}}_{ml},{\bf{e}}_{1,ml}\right). \eqno(9)]

Here [{\bf{r}}_{ml}] = (x0,ym,zl), [{\bf{k}}_{ml}] is a wavevector of the (m,l)-beam, Cml is the factor that determines the (m,l)-beam amplitude, and [{\bf{e}}_{1,ml}] is a unit vector determining the propagation direction of waves in the (m,l)-beam.

We have to reasonably determine the wavevectors [{\bf{k}}_{ml}], the multipliers Cml, the unit vectors [{\bf{e}}_{1, ml}] and the parameter σ of Gaussian beams. First of all, we define the wavevectors [{\bf{k}}_{ml}] of the Gaussian beams used. We represent the X-ray wavefield in the plane x = x0 in the neighborhood of the point [{\bf{r}}_{ml}] = (x0,ym,zl) in the form E = [E_{0ml}\exp[ik_{0}(x-x_{0})]\exp(i\phi)], where ϕ is a real local phase of the wave. Expanding the wave phase ϕ in a Taylor series in the neighborhood of the point [ {\bf{r}}_{ml}], we obtain

[\eqalignno{ E &= E_{0ml} \exp\left[ik_{0}\left(x-x_{0}\right)\right]\exp\left(i\phi\right) \cr& \simeq E_{0ml} \exp\Bigg\{i\Bigg[\phi\left({\bf{r}}_{ml}\right) + \left(k_{0}+\left.{{ {{\partial}}\phi}\over{{{\partial}}x}}\right|_{{\bf{r}} ={\bf{r}}_{ml}}\right) \left(x-x_{0}\right) \cr& \quad + \left.{{{\partial}\phi}\over{{\partial}y}}\right|_{{\bf{r}}={\bf{r}}_{ml}} \left(y-y_{m}\right) + \left.{{{\partial}\phi}\over{{\partial}z}}\right|_{{\bf{r}}={\bf{r}}_{ml}} \left(z-z_{l}\right)\Bigg]\Bigg\}. &(10)}]

We see that the definition of the local wavevector [{\bf{k}}_{ml}] by the relations

[k_{ml,y}= \left.{{{\partial}\phi}\over{{\partial}y}}\right|_{{\bf{r}}_{ml}}, \quad k_{ml,z}= \left.{{ {\partial}\phi}\over{{\partial}z}}\right|_{{\bf{r}}_{ml}}, \quad k_{ml,x}=k_{0}+ \left.{{{\partial}\phi}\over{{\partial}x}}\right|_{{\bf{r}}_{ml}} \eqno(11)]

is reasonable. From equation (4)[link] it follows that kml,x can be represented as

[k_{ml,x} = k_{0} - {{k_{ml,y}^2 +k_{ml,z}^2}\over{2k_0}}.]

Obviously, the unit vectors [{\bf{e}}_{1,ml}] in (9)[link] are directed along the [{\bf{k}}_{ml}] vector, that is, [{\bf{e}}_{1,ml}] = [{{{\bf{k}}_{ml}}/{|{\bf{k}}_{ml}|}}]. It is also obvious that [{\bf{e}}_{2,ml}] = [{\bf{j}}\times{\bf{e}}_{1,ml}], [{\bf{e}}_{3,ml}] = [{\bf{e}}_{1,ml}\times{\bf{e}}_{2,ml}]. On the κ grid, the components of the vector [{\bf{k}}_{ml}] can be approximated by the finite differences kml,y[({\phi_{m+1,l}-\phi_{m-1,l}})/2h], kml,z[({\phi_{m,l+1}-\phi_{m,l-1}})/2h], where h is the mesh step.

In our work, the X-rays propagate at very small angles to the OX-axis and for all points of the plane x = x0 the approximate relations [ {\bf{k}}_{ml}({\bf{r}}-{\bf{r}}_{ml})] ≃ 0 are valid because the vectors [{\bf{k}}_{ml}] are almost orthogonal to the plane x = x0; respectively, we obtain

[E\left(x=x_{0},y,z\right) \simeq \sum\limits_{m,l} C_{ml}\exp\left[-{{\left(y-y_{m}\right)^{2}\,+\,\left(z-z_{l}\right){}^{2}} \over {2\sigma^{2}}}\right]. \eqno(12)]

We consider on the plane x = x0 an arbitrary closed star domain Ω containing the point [(m^{\prime}, l^{\,\prime})] of the grid and we find the coefficients Cml from the condition

[\eqalignno{ \int\limits_{\Omega}E&\left(x=x_{0},y,z\right)\,{\rm{d}}y\,{\rm{d}}z \, \simeq \cr& \int\limits_{\Omega}\sum\limits_{m,l}C_{ml}\exp\left[-{{\left(y-y_{m}\right)^{2}\,+\,\,\left(z-z_{l}\right)^{2}}\over{2\sigma^{2}}}\right]\,{\rm{d}}y\,{\rm{d}}z. &(13)}]

The relation (13)[link] is assumed to be valid for any domain Ω whose diameter is much greater than both h and σ. From (13)[link] then follows

[\sum\limits_{\left({m,l}\right)\,\in\,\Omega} \!\!\! E_{0ml}h^{2} \,\,\simeq \sum\limits_{\left({m,l}\right)\,\in\,\Omega} \!\!\! 2\pi\sigma^{2}C_{ml}. \eqno(14)]

From (14)[link] it follows that CmlE0mlh2/(2πσ2), where E0mlE(x0, ym, zl). The left-hand side of the formula (14)[link] is obtained by applying the usual quadrature formula to (13)[link] for calculating the integral over the region Ω. When calculating the integrals on the right side of (14)[link], assuming the beam width σ to be small compared with the diameter of the Ω region, we used the fact that the Gaussian functions decrease quickly with increase of [(yym)2 + (zzl)2], and we replaced the integrals over Ω of the Gaussian functions with integrals over the entire plane R2.

In the derivation (14)[link] we used the simplifying approximation [{\bf{k}}_{ml} \parallel OX]. In a more general case, the wave can propagate at a small angle to the OX-axis. Understanding the principle of derivation of the formula for the factors Cml, we below formulate the reasoning that does not use the assumption [{\bf{k}}_{ml} \parallel OX]. We draw through the point [(x_ {0}, y_ {m^ \prime}, z_ {l^{\,\prime}})] a plane [\Pi] perpendicular to the vector [{\bf{k}}_{m^{\prime} l^{\,\prime}}]. We introduce a special local coordinate system. Let the [OX ^ \prime]-axis be along the direction indicated by the vector [{\bf{e}}_{1, m^{\prime} l^{\,\prime}}], the [OY ^ \prime]-axis be along [{\bf{e}}_{2, m^{\prime} l^{\,\prime}}] = [{\bf{j}} \times {\bf{e}}_{1, m^{\prime} l^{\,\prime}}] and the [OZ ^ \prime]-axis be along [{\bf{e}}_{3, m^{\prime} l^{\,\prime}}] = [{\bf{e}}_{1, m^{\prime} l^{\,\prime}} \times {\bf{e}}_{2, m^{\prime} l^{\,\prime}}]. We have taken into account here that the vector [{\bf{k}}_{m^{\prime} l^{\,\prime}}] may be nonparallel to the optical axis of the system. We consider an arbitrary closed star-shape domain [\Omega ^ \prime] in the plane [\Pi] containing the point [(m^{\prime}, l^{\,\prime})] of the plane [\Pi]. We assume that the domain [\Omega ^ \prime] is not very large; in this case, [{\bf{k}}_{ml}][{\bf{k}}_{m^{\prime} l^{\,\prime} }] is fulfilled for all points (m,l) of the mesh inside [\Omega^\prime].

We repeat the above derivation for Ω′ and obtain for Cml a formula that coincides with the previous one,

[C_{ml} = E_{0ml}h^{2}/\left(2\pi\sigma^{2}\right), \eqno(15)]

but the new derivation of the formula (15)[link] is more accurate and takes into account that the vectors [{\bf{k}}_{ml}] may be nonparallel to the optical axis, but can form small angles with the optical axis.

In Fig. 1[link] for the planar case, an example of approximation of the function E(x0,y) = [E_0\exp\{\,i\,[k_0(x-x_0)]\}\exp(-y^2/2\alpha^2)] (E0 = 1.6 × 107 V m−1, α = 100 µm, k0 = ω0/c, ω0 = 6π × 1018 Hz) on the line x = x0, using a linear combination of the functions

[\eqalignno{ \exp&\left[\,i\,{\bf{k}}_{m}\left({\bf{r}}-{\bf{r}}_{m}\right)\right] {{ 1 }\over{ 1+i\,\left[\left({\bf{r}}-{\bf{r}}_{m}\right) \,{\bf{e}}_{1,m}\right] /\left(|{\bf{k}}_{m}|\sigma^{2}\right) }} \cr& \,\,\times \exp \left( -{{ \left[\left({\bf{r}}-{\bf{r}}_{m} \right) \,{\bf{e}}_{2,m}\right]^{2} }\over{ 2 \left\{\sigma^{2}+ i\,\left[\left({\bf{r}}-{\bf{r}}_{m}\right)\,{\bf{e}}_{1,m}\right]/ |{\bf{k}}_{m}|\right\} }} \right), &(16) }]

is shown. The coefficients of the linear superposition are determined by the formula Cm = E0mh/[(2π)1/2σ], the components of the wavevector km are determined by the formula (11)[link], rm = (x0, ym), km = |km|, [{\bf{e}}_{1,m}] = [{\bf{k}}_m/k_m]. Some difference between these formulas and the above formulae (9)[link] and (15)[link] is explained by the fact that in this example we have only one transverse coordinate. We see that the system of functions (16)[link] used allows us to approximate a given function very accurately on the line x = x0. Since we apply a superposition of almost exact solutions of the Helmholtz equation, which satisfy the boundary condition with high accuracy, this method allows solving the problem with high accuracy.

[Figure 1]
Figure 1
Approximation of the wave electric field on the line x = x0 by a linear superposition of functions of the form (16)[link]. The blue dots denote the grid points, and they are also the centers of the Gaussian functions. h = σ = 20 µm. The exact function and the approximating function practically merge.

4. Oriented Gaussian beams for describing the X-ray propagation through optical elements

Now consider the propagation of X-rays through lenses. To begin with, the consideration of the propagation of X-rays through a single lens will suffice. X-rays penetrate any material and propagate in any material with very little reflection and absorption. This allows the boundary condition of the continuity of the function E to be used at the lens–air boundaries (Kohn, 2002[Kohn, V. G. (2002). J. Environ. Toxicol. Publ. Heal. Lett. 76, 701-704.]; Lengeler et al., 1999[Lengeler, B., Schroer, C., Tümmler, J., Benner, B., Richwin, M., Snigirev, A., Snigireva, I. & Drakopoulos, M. (1999). J. Synchrotron Rad. 6, 1153-1167.]). Materials used in X-ray optics for lens manufacturing have such small parameters β and δ that equation (3)[link] in the case of propagation of an X-ray wave inside the lens can be represented in the following form,

[{{ {{\partial}} A} \over { {{\partial}} x}}-ik_{0}(-\delta+i\beta) \, A = {{i} \over {2k_{0}}} \left({{ {{\partial}} ^{2}A} \over { {{\partial}} y^{2}}}+{{ {{\partial}} ^{2}A} \over { {{\partial}} z^{2}}}\right). \eqno(17)]

The substitution

[A(x,y,z) = B(x,y,z)\exp\big[ik_{0}(-\delta+i\beta)\,x\big] \eqno(18)]

transforms equation (17)[link] into a standard paraxial equation,

[{{ {{\partial}} B} \over { {{\partial}} x}} = {{i} \over {2k_{0}}}\left({{ {{\partial}} ^{2}B} \over { {{\partial}} y^{2}}}+{{ {{\partial}} ^{2}B} \over { {{\partial}} z^{2}}}\right). \eqno(19)]

The method of solving equation (19)[link], and also the Helmholtz equation (1)[link], is described in Section 2. We can use the formulas (18)[link], (5)[link], (7)[link], (10)[link] and (11)[link] for obtaining the solution to equation (1)[link] to describe the propagation of an X-ray wave through a lens.

The formulas (18)[link] and (19)[link] show that if an X-ray wave encounters a lens with a local thickness d(y, z) then at the exit from the lens at the lens surface the X-ray wave obtains an additional phase ϕ(y, z) = k0(−δ + iβ)d(y, z), which is also equivalent to the fact that the function describing the electric field obtains the factor [\exp[ik_{0}(-\delta+i\beta)\,d(y,z)]]. The scheme explaining the local thickness of the lens is shown in Fig. 2[link].

[Figure 2]
Figure 2
Schematic representation of the lens and its local thickness. The functions xL(y, z), xR(y, z) describe the surface of the lens closer to the light source and the lens surface farther from the light source. For the parabolic biconcave lens xL(y,z) = max{ xC(1/2R) (y2+z2)(Wsm/2), xC-(W/2) }, xR(y, z) = min{xC + (1/2R)(y2+z2) + (Wsm/2), xC + (W/2)}, where R denotes the curvature radius of the lens, W is the maximum thickness of the lens and Wsm is the minimum thickness of the lens; xC is a position of the geometric center of the lens.

Since the lens thickness is much less than the focal length, using the thin lens approximation is permissible and convenient. Within this approximation we take the lens center at the coordinate xC of the lens on the OX-axis and assume that the lens lies on the plane x = xC. Let the notations (xC − 0) and (xC + 0) call the limit on the left and the limit on the right to the point xC; these notations will be used for writing the wave fallen on the thin lens from the left and the wave passed through the thin lens to the right. When an X-ray wave passes through a thin lens, the wave phase (and, correspondingly, the electric field) changes abruptly. We suppose that a narrow Gaussian beam G(r, xC − 0, ym, zl) described by the formula (7)[link] falls onto the lens from the left. Also, let the width σ of a Gaussian beam be much smaller than the lens aperture. At the lens output at the lens surface, we obtain the wave

[\eqalignno{ \exp&\Bigg\{ i k_{0}(-\delta+i\beta)\Bigg[d\left(y_{m},z_{l}\right)+{{ {{\partial}} d\left(y_{m,}z_{l}\right)} \over { {{\partial}} y}}\left(y-y_{m}\right)\cr& \quad+{{ {{\partial}} d\left(y_{m,}z_{l}\right)} \over { {{\partial}} z}}\left(z-z_{l}\right)\Bigg]\Bigg\} \, G\left({\bf{r}},x_{\rm{C}}+0,y_{m},z_{l}\right). &(20) }]

We have expanded the lens thickness in a Taylor series in the neighborhood of (ym, zl) and have taken into account linear terms of the series. Since the lens is considered as a thin one, the wave parameters changed a lot when the wave passed through the lens.

Now we consider the same question in terms of oriented Gaussian beams (8)[link]. The X-ray wave [G({\bf{r}},{\bf{r}}_{ml},{\bf{e}}_{1,ml})], having rml = (xC, ym, , zl) and the unit vector e1, ml indicating the direction of wave propagation parallel to the OX-axis, fell on the lens from the left. At the output, after the lens, we obtained the wave [\exp[ik_{0}(-\delta+i\beta)\,d(y_{m},z_{l})]\,G({\bf{r}},{\bf{r}}_{ml},{\bf{e}}^{\,\prime}_{1,ml})]. The parameters of this wave are determined by the formula (10)[link], which leads to

[\eqalign{ k_{ml,y}^{\,\prime} & = k_{ml,y}- k_{0}\,\delta {{ {{\partial}} d\left(y_{m},z_{l}\right) }\over{ {{\partial}}y}}, \cr k^{\,\prime}_{ml,z} & = k_{ml,z}-k_{0}\,\delta {{ {{\partial}} d\left(y_{m},z_{l}\right) }\over{ {{\partial}} z}} , \cr k_{ml,x}^{\,\prime} & = k_{0}- {{ k^{\,\prime\,2}_{ml,y}+k^{\,\prime\,2}_{ml,z} }\over{ 2k_0}} .} \eqno(21)]

Thus, the wave has changed its propagation direction and phase, but otherwise has remained unchanged. The addition to the wavevector [{\bf{k}}_{ml}] in (21)[link] has been generated by the factor before G(r, xC + 0, ym, zl) in (20)[link].

Now we generalize the obtained result to the case when the wave in the form of a superposition of narrow Gaussian beams [E\left({\bf{r}}\right) = \textstyle\sum_{m,l}C_{ml}G({\bf{r}},{\bf{r}}_{ml},{\bf{e}}_{1,ml})] falls from the left on the plane x = xC containing a lens. Then after the lens we obtain the wave

[E\left({\bf{r}}\right) = \sum_{m,l} C^{\,\prime}_{ml} G\left({\bf{r}},{\bf{r}}_{ml},{\bf{e}}^{\prime}_{1,ml} \right). \eqno(22)]

This wave is a superposition of oriented Gaussian beams [G({\bf{r}},{\bf{r}}_{ml},{\bf{e}}^{\prime}_{1,ml})] in which the parameters [ {\bf{k}}^{\prime}_{ml}] and [{\bf{e}}^{\prime}_{1,ml}] are determined by the formula (21)[link] and

[C^{\,\prime}_{ml} = \exp\big[ ik_{0}(-\delta+i\beta) \, d\left(y_{m},z_{l}\right)\big] C_{ml}. \eqno(23)]

In the approach used, the formula (22)[link] describes very accurately the propagation of X-ray waves in the air after the lens. High accuracy of the formulas (22)[link], (23)[link] and (21)[link] is due to the lens symmetry, so that the lens center is clear. Also, the explicit formulas for very narrow Gaussian X-ray beams make it possible to calculate explicitly the phase shifts with allowance for fine details.

In the derivation of the formulas (21)[link] and (23)[link], we assumed that the X-ray wave propagates in the lens along the optical axis. This is a justified approximation, since X-ray waves have a great penetrating power and practically do not change their propagation direction. Nevertheless, if we use a lot of lenses, about a hundred or more, the focal length of the optical system is relatively small. In this case, when calculating the local phase shifts, it is necessary to take into account the local slopes of the wave propagation direction with respect to the optical axis. Corresponding amendments can be taken into account; they will somewhat complicate the formulas. If we use only a few dozen lenses, the local slopes of the direction of wave propagation inside the lens are very small and practically do not affect the result.

We hope that the calculation method presented here will be convenient and adequate for studying the effects of lens defects (foreign inclusions, oxides, caverns) on focusing and imaging. A theoretical study of the influence of lens defects on focusing and imaging may yield more results than an experimental one, since we are unlikely to be able to obtain from the experiments the internal lens structure in detail. We believe that modeling the internal structure of lenses as random would be preferable. Obviously, the effect of one small defect in the lens will be small, but the effects of many defects in the lens can accumulate, especially if we apply a lot of lenses. Therefore, in calculations, high accuracy is required. When examining the impact of lens defects on the focusing and imaging, it is reasonable to not consider lenses as thin ones but to cut the lenses into layers. This will allow the introduction of a local wavevector inside each lens, and a better understanding of the wave behavior when the wave will meet a defect or multiple defects.

5. An estimate of the acceptable step h and the method error

When (k0σ)2 ≃ 109 ≫ 1, the Gaussian beams (8)[link] satisfy well the Helmholtz equation (1)[link]. Therefore, we can consider them as practically exact solutions of the Helmholtz equation (1)[link]. A linear superposition of exact solutions (12)[link] for xx0 gives an exact solution (9)[link] of the Helmholtz equation (1)[link]. Thus, the error of the approximate solution is determined by the error in calculating [{\bf{k}}_{ml}] and the error of approximating the boundary condition with a series (12)[link]. From the formula (10)[link] it follows that the error in calculating the local wavevectors is O(h2). This means that the error decreases in proportion to h2 when h decreases. As for the approximations (12)[link]–(15)[link], the analysis shows that the error in this case also decreases in proportion to h2. Below we will apply the Runge rule (28)[link] to check the calculation accuracy. The derivation of the Runge rule is given by Kshevetskii et al. (2016[Kshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305-1314.]) and Kraus & Langer (2007[Kraus, J. & Langer, U. (2007). Editors. Lectures on Advanced Computational Methods in Mechanics. p. 164. Berlin: de Gruyter.]). We have to substitute k = 2 into the Runge rule formula (28)[link] because the approximation order is equal to two.

The problem of focusing X-rays is solved step by step. First we calculate the wavefield in the plane in front of the lenses. Then we compute the X-ray propagation through the lenses. After that we calculate the wavefield after the lenses on the plane perpendicular to the optical axis. We finally use the wave electric field on the plane after the lenses for calculating the focusing. As shown in the work of Kshevetskii et al. (2016[Kshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305-1314.]), in the case of many lenses the wavefield on the plane after the lenses is fast-oscillating. The calculation quality is determined by how accurately we are able to describe this fast-oscillating field.

To obtain an estimate of the acceptable step h, we assume that the phase of the wave incident on the lens system is approximately constant. The wave phase is changed significantly after the X-ray wave has passed through the lenses. The formula

[d(y,z) = {{y{}^{2}+z{}^{2}} \over {R}}+W_{\rm sm} \eqno(24)]

describes the dependence on the transverse coordinates of the thickness of each lens within the aperture. Here Wsm is the smallest distance between the parabolic surfaces of the lens and R is the curvature radius of the lens surfaces. According to the formulas (18)[link] and (19)[link], when the X-ray wave has propagated through a system of N lenses, the X-ray wave has the phase

[\phi \simeq k_{0}\delta N \left({{y^{2}+z^{2}}\over{R}}+W_{\rm sm}\right), \eqno(25)]

dependent on the transverse coordinates. We have neglected the wave attenuation.

The approximation (9)[link] assumes that the function (25)[link] locally in some neighborhood of each point and at distances much larger the grid steps h can be approximated by a plane. Then we can introduce the local wavevector [{\bf{k}}\left(x,y,z\right)] and the approximation (9)[link] is effective. Obviously, the function (25)[link] can be approximated in the neighborhood of each mesh point with a plane if within the aperture the conditions

[\left|{{{{\partial}}^{2}\varphi}\over{{{\partial}}y^{2}}}\,h^{2}\right| \ll 1, \qquad \left|{{{{\partial}}^{2}\varphi}\over{{{\partial}}z^{2}}}\,h^{2}\right| \ll 1\eqno(26)]

are satisfied. The conditions (26)[link] give a restriction on the step h of the mesh κ,

[h \ll \left( {{ cR }\over{ 2\omega_{0}\delta{N} }} \right)^{1/2}. \eqno(27)]

Actually we use [h \simeq [(0.05cR)/(\omega_0\delta{N})]^{1/2}] in our calculations.

In the work of Kshevetskii et al. (2016[Kshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305-1314.]), the question of a fundamental restriction to the step h was considered, and the condition [h \ll ({{cR}/{\omega_0}})\,\delta{N}a] was derived. Here a is the lens aperture. The condition (27)[link] does not depend on the lens aperture and in practice it turns out to be substantially softer than the condition described by Kshevetskii et al. (2016[Kshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305-1314.]). It allows the grid step h to be increased by about an order of magnitude. This possibility is explained by the fact that we use in (9)[link] inclined Gaussian beams, which, due to their slopes, take into account fast oscillations on the OYZ-plane of the wave electric field. In the derivation of (27)[link], the lenses are assumed to be ideal, having a parabolic shape. When calculating the X-ray propagation through non-ideal-shaped lenses, the h step may have to be reduced for catching and investigating the effects of non-ideal lenses.

A special issue concerns the influence of the edges of lens apertures on focusing and imaging. The derivatives of the refractive index n change abruptly at the edges of the lens apertures, and the wave phase φ also changes abruptly. Therefore, the condition (26)[link] is violated at the edges of lens apertures, and even [\left|{{ {{\partial}} ^{2}\varphi}/{ {{\partial}} x^{2}}}\right|], [\left|{{ {{\partial}} ^{2}\varphi}/{ {{\partial}} z^{2}}}\right|] tend to infinity. Let us consider all points (m, l) in the sum (9)[link], which are in a thin ring of width containing the contribution of the edge of the lens aperture. The contribution of this thin ring into the total wavefield does not exceed E*/σ2. Here E* is the maximum value of |E| in this ring. Since the width of the considered ring can be taken as arbitrarily small, the contribution of the edge of the lens aperture in the formula (9)[link] to the general wave picture can be neglected. This conclusion can also be justified by the fact that the beam of X-ray waves incident on the lens usually has an effective width that is smaller than the lens aperture size.

Nevertheless, the influence of the edges of lens apertures on the wave propagation and focusing exists. Formula (9)[link] allows calculation with high accuracy of this influence on the wavefield in the focal plane; this will be demonstrated later.

5.1. Volume of calculations

The computational complexity (number of operations) of the presented method using Gaussian beams is estimated as MIMF. Here MI is the number of points digitizing the electric field after the lenses. MF 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 MI. The number MF is determined only by the desirable resolution of the picture of interest.

At present, Fourier methods are often used for calculations in X-ray optics. It is known that the computational complexity of the FFT is estimated as Mlog2(M), where M is the number of points used in the FFT method (Cooley et al., 1967[Cooley, J., Lewis, P. & Welch, P. (1967). IEEE Trans. Audio Electroacoust. 15, 76-79.]). The FFT method is very fast. However, the FFT technology based on FFT is somewhat lacking in numerical robustness and stability (Chubar et al., 2017[Chubar, O., Rakitin, M., Chen-Wiegart, Y., Chu, Y. S., Fluerasu, A., Hidas, D. & Wiegartr, L. (2017). Proc. SPIE, 10388, 1038805.]). Currently, work is underway to develop software that allows automatic control of the accuracy of calculations using FFT (Chubar et al., 2017[Chubar, O., Rakitin, M., Chen-Wiegart, Y., Chu, Y. S., Fluerasu, A., Hidas, D. & Wiegartr, L. (2017). Proc. SPIE, 10388, 1038805.]). Therefore, the question of calculation accuracy control is an important part of modern research related to the computational problems of X-ray optics. When the number of lenses increases, the frequency of oscillations of the wavefield in the plane behind the lenses perpendicular to optical axis increases fast. Therefore, the number of Fourier terms that have to be taken into account also increases fast. Standard methods for estimating the error of truncating Fourier series are applied to functions specified analytically, and they 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 applying standard formulas for estimating the truncating error is difficult. In addition, studies on the calculation accuracy of X-ray focusing have revealed that the sensitivity of various parameters characterizing the focusing to the calculation quality is very different (Kshevetskii et al., 2016[Kshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305-1314.]). Therefore, evaluating the accuracy of calculations is a serious problem and the calculation accuracy issue should be formulated for specific characteristics of focusing, and not as a general question.

Thus, the method developed in this paper has some advantages in cases where we need to guarantee the reliability and accuracy of calculations. At the same time, the developed method is quite universal and uncomplicated.

6. Technique for modeling the propagation and focusing of X-ray waves in terms of oriented Gaussian beams

The propagation and focusing of X-ray waves is investigated for the following conditions. We consider the situation that some source of X-ray waves with wave frequency ω0 = 6π × 1018 Hz produces in the plane x = x0 (x0 = 0 m) a Gaussian distribution of electric wavefield intensity with a width [\varrho] = 100 µm. For beryllium lenses and given frequency of X-ray waves, β = 3.1801 × 10−10 and δ = 2.2156 × 10−6. Then the X-rays propagate further to a distance of 40 m, and then the waves propagate through a system of 30 parabolic lenses made of beryllium with a curvature radius of 50 µm and smallest width of 30 µm. After the X-rays have propagated through the lenses, the X-rays are focused, and we study them in the focal plane.

The general scheme of calculating the wave propagation is shown in Fig. 3[link]. The problem is solved as follows. The specified wave electric field from a source of monochromatic X-ray waves is digitized on the mesh on the plane x = x0, and then the X-ray propagation to the region x > x0 is calculated using formulas (9)[link], (8)[link], (10)[link] and (15)[link]. The resulting wave solution is digitized again on the mesh on the plane x = x1. The wave propagation through the first lens is calculated using the formulas (22)[link], (21)[link] and (23)[link]. The resulting wavefield after the first lens is digitized on the mesh on the plane x = x2. The wave propagation through the second lens is calculated using formulas (22)[link], (21)[link] and (23)[link]. The obtained wavefield after the second lens is digitized on the mesh on the plane x = x3. The formulas (22)[link], (21)[link] and (23)[link] are used for calculating the wave propagation through the third lens, and so on. We repeat this until the wavefield is calculated after the last lens corresponding to the coordinate x30. The resulting wavefield after the last lens is used to find the position of the best focus. It is digitized on the mesh on the plane spaced at a distance corresponding to the position of the best focus, which we denote x31. Finally, we investigate the behavior of the digitized electric field on the focal plane x = x31.

[Figure 3]
Figure 3
Optical system and general scheme of calculations.

6.1. Results of calculations

The calculations are performed for a one-dimensional case. We consider the case where the wave electric field of the X-ray beam from the synchrotron far in front of the lenses in the plane x = x0 is described by the formula E(0,y) = [E_0\exp(-y^2/2\alpha^2)] (where E0 = 1.6 × 107 V m−1 and α = 100 µm). Table 1[link] presents the calculation results obtained by approximating the wave electric field with the help of 1000, 2000 and 4000 oriented Gaussian beams, respectively. The widths σ = h of these Gaussian beams are equal to 1 µm, 0.5 µm and 0.25 µm, respectively. The |E| = |A| function takes the largest value at the distance of 0.36658 m after the lens system. Table 1[link] lists the values of |E| and the FWHM at the distance equal to the focal length (x = x31).

Table 1
Results of the computation of the maximum value of |E| and FWHM at the focal distance (FD: 0.36658 m after the symmetry center of last lens) for different values of h; R = 20000 m−1 and a system of 30 beryllium lenses is used

Number of points h (µm) Maximum of |E| at FD FWHM at FD (nm)
1000 1 5.2063 × 108 191.6
2000 0.5 5.6027 × 108 189.5
4000 0.25 5.6299 × 108 189.0

Using the data in Table 1[link] and applying the Runge rule, we estimate the calculation errors of the results given in Table 1[link]. The Runge rule (Kraus & Langer, 2007[Kraus, J. & Langer, U. (2007). Editors. Lectures on Advanced Computational Methods in Mechanics. p. 164. Berlin: de Gruyter.]) allows to estimate computational errors of values independent of coordinates, such as the FWHM and the focal length. The electric field depends on the coordinates, but the maximum value of |E| in the focal spot is independent of the coordinates. Therefore, we can control the accuracy of calculating this physical quantity using the Runge rule, too. The error Q of computing the value of Z within the Runge rule can be estimated using the approximate equality

[Q \simeq \left|{{Z_{h_{1}}-\,Z_{h_{2}}}\over{\left({{h_{2}}/{h_{1}}}\right)^{k}\,-\,1}}\right|, \eqno(28)]

where Zh1 and Zh2 are approximate values of Z obtained with the space steps h1 and h2, respectively, and [h_{2}\,\lt\,h_{1}]. The derivation of the formula (28)[link] is given by Kshevetskii et al. (2016[Kshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305-1314.]).

According to the analysis performed in Section 5[link], the oriented Gaussian beam method is a method of the second order of accuracy with respect to the step h. Therefore, we have to take k = 2 in (28)[link]. Applying the Runge rule (28)[link] to the results of computations for the case h1 = σ1 = 1 µm and and for the case h2 = σ2 = 0.5 µm, we have found out that the calculation error of the FWHM is approximately equal to 1.5% of its value. The error of computation of a maximal value of the wave electric field in the focal spot is 9% of its value. These results show that the amount of M = 1000 Gaussian beams with h = σ = 1 µm is not sufficient for a high-precision calculation of the maximum |E| in the focal spot. The results of the calculation for the case of h1 = σ1 = 0.5 µm and h2 = σ2 = 0.25 µm show that the FWHM calculation error is approximately 0.4%, and the error in calculating the maximum value of the wave electric field is 0.6%. Therefore, we see that the choice of h = σ = 0.5 µm gives highly accurate results.

Figs. 4[link] and 5[link] shows the behavior of the function |E| at the focal plane x = x31. The curve shape in the center of Fig. 4[link] is close to the Gaussian function, and this is a so-called focal spot. Fig. 5[link] shows the fine structure of the function |E|, formed in some neighborhood of the focal spot. This fine structure consists of a set of Gaussian-like peaks. The maximum value of the function |E| in this fine structure is at least two orders of magnitude less than the maximum value of the function |E| in the focal spot. The values of the local maxims |E| in the fine structure decrease with increasing distance from the center. To make sure that this fine structure is not a computational artifact, but is a physical effect, in Fig. 5[link] the results of calculations with the ordinary and double steps h are pictured. The results of calculations with ordinary and double steps coincide with high accuracy. This indicates a physical nature of the calculated fine structure.

[Figure 4]
Figure 4
Focus spot for 2000 points.
[Figure 5]
Figure 5
Detailed focus spot for 2000 (the red line) and 4000 points (the blue line).

Comparing the results obtained with different widths α of the initial Gaussian distribution of the electric wavefield at x = x0, we have found that the amplitudes of the Gaussian-like peaks of the fine structure in the focal plane x = x31 correlate with |E| at the edges of the lens apertures. This observation shows that the fine structure of the wave electric field in the focal plane has arisen as a result of X-ray diffraction at the edges of the lens apertures.

7. Conclusions

A high-precision method for solving X-ray optics problems, based on the use of superposition of oriented Gaussian beams, is developed. A Gaussian beam is an exact solution to the paraxial equation of optics, and for appropriate beam parameters (the wavelength is much smaller than the beam width) it can be interpreted as an almost exact solution of the Helmholtz equation. In an oriented Gaussian beam, the waves propagate at a given angle to the optical axis. Some relations that determine the optimal choice of Gaussian beam parameters and formulas for linear superposition coefficients are obtained. Also, the formulas that determine the variation of the parameters of oriented Gaussian beams and the superposition coefficients when the X-ray wave propagates through the lens are derived. The propagation of a Gaussian beam through a lens leads to changes in the propagation direction of the wave and the complex amplitude of the electric field in the beam.

For the method presented, tools for controlling the accuracy and reliability of calculations are developed. It is shown that different characteristics of the focal spot have a significantly different sensitivity to the computation quality. Among the parameters analyzed ([\max\left(|E\,|\right)], FWHM, focus length), the maximum amplitude of the electric field in the focal spot is the most sensitive and exacting to the computation quality. A clear formula for the mesh step (the number of Gaussian beams used) required for high-precision calculations is derived.

As an application, the problem of focusing of X-rays by a system of 30 beryllium lenses in an X-ray microscope was solved. The high accuracy of the presented method made possible calculating not only the X-rays focusing but also the fine structure of the wave electric field on the focal plane. The fine structure arises from the diffraction of X-rays at the edges of the lens apertures.

The developed high-accuracy modeling method is quite universal and is designed to solve various problems of X-ray optics and, in particular, to study the effect of various defects in lenses (cavities, foreign inclusions, oxides) on focusing and imaging. The presented method allows solving problems with complex lens shapes or calculating the propagation of X-ray waves through a sample with a complex internal structure, and with control of the accuracy of calculations.

References

First citationBabich, V. M. & Buldyrev, V. S. (1991a). Asymptotic Methods in Short-Wavelength Diffraction Theory. Oxford: Alpha Science International.  Google Scholar
First citationBabich, V. M. & Buldyrev, V. S. (1991b). Short-Wavelength Diffraction Theory: Asymptotic Methods. Berlin: Springer Verlag.  Google Scholar
First citationBahrdt, J. (2007). Phys. Rev. ST Accel. Beams, 10, 060701.  CrossRef Google Scholar
First citationChabory, A., Sokoloff, J., Bolioli, S. & Combes, P. F. (2005). C. R. Phys. 6, 654–662.  CrossRef CAS Google Scholar
First citationChabory, A., Sokoloff, J., Bolioli, S. & Elis, K. (2010). Proceedings of the 4th European Conference on Antennas and Propagation (EuCAP'2010), 12–16 April 2010, Barcelona, Spain, pp. 1–5.  Google Scholar
First citationChubar, O., Rakitin, M., Chen-Wiegart, Y., Chu, Y. S., Fluerasu, A., Hidas, D. & Wiegartr, L. (2017). Proc. SPIE, 10388, 1038805.  Google Scholar
First citationCooley, J., Lewis, P. & Welch, P. (1967). IEEE Trans. Audio Electroacoust. 15, 76–79.  CrossRef Google Scholar
First citationDeschamps, G. (1972). Proc. IEEE, 60, 1022–1035.  CrossRef Google Scholar
First citationDeschamps, J., Courjon, D. & Bulabois, J. (1983). J. Opt. Soc. Am. 73, 256.  CrossRef Google Scholar
First citationGhannoum, I., Letrou, C. & Beauquet, G. (2009). Proceedings of RADAR 2009 – International Radar Conference `Surveillance for a Safer World', 12–16 October, 2009, Bordeaux, France.  Google Scholar
First citationGoodman, J. W. (1996). Introduction to Fourier Optics, 2nd ed. New York: McGraw-Hill.  Google Scholar
First citationIshimaru, A. (1991). Electromagnetic Wave Propagation, Radiation and Scattering. Englewood Cliffs: Prentice Hall.  Google Scholar
First citationKeller, J. B. & Streifer, W. (1971). J. Opt. Soc. Am. 61, 40–43.  CrossRef Google Scholar
First citationKogelnik, H. (1965). Appl. Opt. 4, 1562–1569.  CrossRef Google Scholar
First citationKohn, V. G. (2002). J. Environ. Toxicol. Publ. Heal. Lett. 76, 701–704.  Google Scholar
First citationKohn, V. G. (2012). J. Synchrotron Rad. 19, 84–92.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationKraus, J. & Langer, U. (2007). Editors. Lectures on Advanced Computational Methods in Mechanics. p. 164. Berlin: de Gruyter.  Google Scholar
First citationKshevetskii, S., Wojda, P. & Maximov, V. (2016). J. Synchrotron Rad. 23, 1305–1314.  CrossRef CAS IUCr Journals Google Scholar
First citationKshevetskii, S. P. & Wojda, P. (2015). Math. Appl. 43, 193–206.  Google Scholar
First citationLengeler, 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
First citationLeontovich, M. A. (1944). Izv. Akad. Nauk. SSSR Ser. Phys. 8, 16.  Google Scholar
First citationLevy, M. (2000). Parabolic Equation Methods for Electromagnetic Wave Propagation. Londo: The Institution of Electrical Engineers.  Google Scholar
First citationRoth, T., Alianelli, L., Lengeler, D., Snigirev, A. & Seiboth, F. (2017). MRS Bull. 42, 430–436.  CrossRef Google Scholar
First citationSnigirev, A., Kohn, V., Snigireva, I. & Lengeler, B. (1996). Nature, 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.

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775
Follow J. Synchrotron Rad.
Sign up for e-alerts
Follow J. Synchrotron Rad. on Twitter
Follow us on facebook
Sign up for RSS feeds