In-plane wavevector distribution in partially coherent X-ray propagation

Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Zhangheng Road 239, Pudong District, Shanghai 201800, People’s Republic of China, University of Chinese Academy of Sciences, Yuquan Road 19, Shijingshan District, Beijing 100049, People’s Republic of China, Shanghai Synchrotron Radiation Facility, Shanghai Advanced Research Institute, Zhangheng Road 239, Shanghai 201204, People’s Republic of China, and Advanced Photon Source, Argonne National Laboratory, 9700 South Cass Avenue, Argonne, IL 60439, USA. *Correspondence e-mail: wangyong@sinap.ac.cn, xshi@anl.gov, tairenzhong@sinap.ac.cn


Introduction
The rapid development of synchrotron radiation technologies has led to orders-of-magnitude increases in the brilliance and spatial coherence of X-ray sources, which in turn have led to the development of many experimental techniques using radiation coherence, such as coherent diffraction imaging (Whitehead et al., 2009;Liang et al., 2015), X-ray lithography (Zhang et al., 2014;Achenbach et al., 2018), X-ray holography (Schaffert et al., 2013;Robisch et al., 2016) and X-ray photon correlation spectroscopy (Stephenson et al., 2009). The design and construction of coherent beamlines requires simulation tools to analyze the coherence properties of the beam. Several software packages have been established to simulate the propagation of a partially coherent beam through different beamline optics. SHADOW3 stores the phase for each ray and sums the coherent rays to simulate the propagation of partially coherent light (Sanchez del Rio et al., 2011). SRW calculates the radiation from a single electron and its propagation through a beamline and sums the intensities incoherently from all electrons (Chubar et al., 2013;Chubar, 2014;Samoylova et al., 2011). PHASE separates a short pulse into many temporal slices, calculates the propagation of every fully coherent slice and sums their intensities (Bahrdt et al., 2014(Bahrdt et al., , 2011. The HYBRID method combines ray tracing and wavefront propagation and performs an intensity convolution (Shi et al., 2014). COMSYL propagates the cross spectral density function along the beamline optics by means of coherent mode decomposition (Glass & Sanchez del Rio, 2017). We developed the MOI code, based on statistical optics, to calculate the propagation of the mutual optical intensity through beamline optics (Meng et al., 2015). Combined with local ray tracing, the ISSN 1600-5775 MOI code was further developed to address non-ideal mirrors (Meng et al., 2017). The main advantage of the MOI code is that it can generate and store mutual optical intensities, including information on the intensity, degree of coherence and statistical phase difference between every pair of points at a previously defined beamline position.
In the original MOI code it was assumed that, within each element of the calculation grid, the wavefield is fully coherent with a constant amplitude and phase. This assumption is easy to satisfy as long as the element is much smaller than the beam and the coherence length. However, achieving a constant phase is difficult when the X-ray wavelength is very short, unless a large number of elements are used. In this paper, we have updated the MOI code by considering the phase distribution inside each element and calculating the onedimensional propagation of mutual optical intensity. The improved code provides accurate results with higher efficiency. Furthermore, the stored in-plane wavevector allows detailed analysis of the phase distribution in a partially coherent wavefield. The MOI code is free to use and can be downloaded at http://www.moixray.cn. The current version is V1.01, which can calculate the one-dimensional propagation of mutual optical intensity through a synchrotron radiation beamline.

Mutual optical intensity propagation
The version of the MOI code developed in our previous work (Meng et al., 2015(Meng et al., , 2017 uses mutual optical intensity to describe a partially coherent wavefield. The wavefield transverse plane is divided into many small elements, within which the wavefield is assumed to be fully coherent with a constant amplitude. However, the phase distribution within each element is not necessarily constant, which is important for the accurate propagation of the mutual optical intensity. In this work, we describe the phase distribution in three steps. Firstly, the phase distribution within each element is assumed to be constant (Meng et al., 2015(Meng et al., , 2017. Here we repeat some of the key equations to help describe the improved code. The one-dimensional mutual optical intensity (Cittert, 1934;Zernike, 1938) between two points (x 1 + x) and x 2 in the transverse plane at a particular propagation position z of the beamline for a single photon energy can be written as where x 1 and x 2 represent the central points of any two small elements, u(x 1 ) and u(x 2 ) correspond to their complex amplitudes, respectively, and x is the relative coordinate and limited inside the element x 1 , thus u(x 1 + x) is equal to u(x 1 ) according to the assumption of full coherence and constant amplitude. The mutual optical intensity propagation through free space can be expressed as (Meng et al., 2015) where ( 1 ) and ( 2 ) are the tilt factors with tilt angles 1 and 2 , is the X-ray wavelength, w 1 and w 2 are the central points of any two elements in the image plane, and r 1 and r 2 are the distances from the object plane to the image plane, x 1 w 1 and x 2 w 2 , respectively, as shown in Fig. 1. When the transverse plane is separated into many small elements, equation (2) can be written as where m and n are the indices of the elements. Owing to the assumption of full coherence and constant complex amplitude within each element, J(x 1 , x 2 ) can be removed from the integration. Equation (3) becomes or where A(x 1 , w 1 ) is the integral within the x 1 element. J xx is an MÂM matrix of the mutual optical intensity for all elements, where M is the number of elements in the source plane. J ww is an NÂN matrix of the mutual optical intensity in the image plane. A xw is an MÂN matrix representing the integral in every element in the source plane with every element in the image plane. A(x 1 , w 1 ) can be calculated numerically within the Fresnel or Fraunhofer approximation if the element size is small enough, or A schematic diagram of mutual optical intensity propagation through free space.
where r 01 is the distance between the central points in the element x 1 in the source plane and element w 1 in the image plane. Secondly, we assume the phase distribution within each element is a linear function, and therefore the in-plane wavevector given by the phase gradient, k(x) = d'/dx, is constant. The in-plane wavevector is a local and one-point function which depends on only one point. It exists in any partially coherent wavefield since the local area can always be regarded as fully coherent. The integral of the in-plane wavevector between two points describes a statistical average phase difference since the two points may be partially coherent and not have a stable phase difference. The mutual optical intensity can also provide the statistical average phase difference between any two points. However, the mutual optical intensity is a non-local and two-point function which depends on two independent points in the transverse plane and needs much more memory to describe the phase distribution information than the in-plane wavevector. We further describe the mutual optical intensity propagation by considering a constant in-plane wavevector within each element. The phase difference between x 1 and (x 1 + x) is k(x 1 )x, where k(x 1 ) is the constant in-plane wavevector in the element x 1 . k(x 1 ) can also be written as k(x 1 ) = (2/) sin[(x 1 )], where (x 1 ) represents the propagation angle between the wavevector and the normal of the transverse plane. Equation (1) can be modified as Equation (5) needs to be rewritten as Thirdly, the in-plane wavevector in the whole transverse plane is step distributed and discontinuous because of the assumption of a constant in-plane wavevector within each element. Linear interpolation is used between each pair of adjacent elements. The in-plane wavevector within an element can be written as k(x) = k 0 (x 1 ) + (2/)a 1 x = 2/ {sin [ 0 (x 1 )] + a 1 x}, where k 0 (x 1 ) is the in-plane wavevector at the central point of the element x 1 , 0 (x 1 ) is the propagation angle between the wavevector at the central point and the element normal, and a 1 is the linear interpolation coefficient. There are two different a 1 in the left and right halves of the element coming from the linear interpolation with the left and right adjacent elements. We describe them both as a 1 for simplicity. The mutual optical intensity propagation is modified with the consideration of the linear interpolation of the in-plane wavevector. The phase difference between x 1 and ( Equation (6) needs to be modified as The integral A(x 1 , w 1 ) in equation (7) can be rewritten as Equations (4), (8) and (9) together are the basic equations for the propagation of the mutual optical intensities with the in-plane wavevectors in the source plane.

In-plane wavevector calculation
The in-plane wavevectors at the image plane must also be calculated and stored to perform the next propagation of the mutual optical intensity. The mutual optical intensity between two points (w 1 + dw 1 ) and w 2 at the image plane can be expressed as The derivative of the mutual optical intensity is Thus, the in-plane wavevector can be written as research papers Using equation (4a), the derivative of the mutual optical intensity can be obtained from Substituting equation (13) into equation (12), we can obtain the in-plane wavevectors in the image plane.

Mutual optical intensity propagation through free space 3.1. Propagation through a fully open aperture
The propagation of partially coherent light through free space is now analyzed using the improved MOI code. The light source is assumed to obey the Gaussian Schell model (GSM) (Starikov & Wolf, 1982;Gori et al., 2001) for simplicity. The mutual optical intensity in the source plane is where J 12 is the mutual optical intensity between two points x 1 and x 2 , represents the coherence length and represents the root-mean-square (r.m.s.) beam size. Note that the GSM is not necessary for the MOI code, and any mutual optical intensity distribution can be used for the calculation. The propagation of partially coherent light through a beamdefining aperture (BDA) is calculated. The optical setup is shown in Fig. 2, where the distance from the source to the BDA is r 12 = 118 m, and the distance from the BDA to the image plane is r 23 = 100 m. The optical parameters of the source are = 2 nm, = 100 mm and = 100 mm. The source is a secondary source defined by a slit and its length is cut to be 2 = 200 mm. A limited source is more common than an unlimited source for demonstration. The  The optical setup for the propagation of a beam through a beam-defining aperture (BDA).  not only with a real source but also with a secondary source defined by a slit in the beamline.
The two propagation steps (from source to BDA and then from BDA to image) of the mutual optical intensity through a fully open BDA are carried out to check the validity of the improved MOI code. The BDA opening is chosen to be 10000 mm in the calculation to ensure that the greatest portion of the beam can pass through the BDA, on which the FWHM of the spot is 1340 mm.
The mutual optical intensity at the image plane is calculated using both the original and improved MOI codes with different numbers of elements, N. The source plane, the BDA plane and the image plane are separated into the same number of elements for each calculation. The distributions of the intensity and degree of coherence relative to the central point at the image plane are shown in Fig. 3. For comparison, we also calculate the propagation directly from the source to the image plane, where only one propagation step is used with N = 200. The one-step propagation result serves as a reference, while the two-step propagation with fully open BDA should give the same results. The intensity profiles and degree of coherence calculated using the original MOI code with different numbers of elements are shown in Figs. 3(a) and 3(b), respectively. The profiles become smoother when the number of elements increases, meaning the calculation accuracy improves. A large number of elements (N ! 500) is necessary for the original MOI code to give reasonable results. The calculation results from the improved MOI code are shown in Figs. 3(c) and 3(d). The intensity and coherence degree profiles are relatively smooth for N = 200. When the number of elements continues to increase, the calculation results show negligible changes and agree well with the reference case.
The r.m.s. spot sizes at the image plane for different numbers of elements are extracted and summarized in Fig. 4(a). The spot sizes obtained from the original MOI code fluctuate dramatically for low numbers of elements. When the number of elements increases to greater than 600, the spot size starts to converge to the reference value, but still with a noticeable deviation (1.8%). This calculation error is the result of the assumption of a constant phase inside each element in the original MOI code, which ignores the beam's local divergence and thus generates a smaller spot size. On the other hand, the results from the improved MOI code coincide well with the reference when the number of elements is greater than 200, with a deviation of only 0.02%, or even less for larger numbers of elements.
The computation times for the two MOI codes with different numbers of elements are shown in Fig. 4(b). With the same numbers of elements, the computation time for the improved MOI code is about twice that of the original code, owing to the additional computation of the in-plane wavevector. However, considering the result in Fig. 4(a) that fewer elements are needed in the improved code, the total computation time is much less to achieve the same or even better accuracy. For example, the computation time is 20 s with a deviation of 0.02% from the reference value for the improved MOI code with N = 200, while these values are 300 s and 0.4%, respectively, for the original MOI code with N = 1800. The improved MOI code can provide results that are more accurate with less computation time, representing an improvement in both accuracy and efficiency.
The in-plane wavevector obtained from the improved MOI code also enables detailed analysis of the in-plane phase distribution. Fig. 5 shows results from the improved MOI code with high spatial resolution using N = 1000. The in-plane wavevector distribution at the image plane with a fully open BDA (black line) and without a BDA (red line) are shown in Fig. 5(a). The two curves are mostly overlapping, with only small deviations near the edges which are caused by the beam diffraction from the BDA. The central part of the in-plane wavevector is linearly proportional to the x coordinate with a slope of 1.44 Â 10 À5 . The corresponding spherical wave radius is 218.1 m, which agrees well with the propagation distance of 218 m. The in-plane wavevector profile with the BDA is shown again in Fig. 5(b) together with the intensity profile for further comparisons. Intensity oscillations with low frequencies originate from the diffraction of the limited source. Thanks to the high accuracy of the improved MOI code, we can even observe the high-frequency oscillations in both intensity and wavevector profiles near the BDA boundaries, as shown in the   inset of Fig. 5(b). These high-frequency oscillations appear prominent only near the edges of the BDA, not in the center, indicating that they are caused by near-field diffraction at the BDA edges. The phase distribution from the k x integral and MOI at the image plane with the BDA are shown in Fig. 5(c). The phase at the center of the transverse plane is assumed to be zero. The phase distribution from MOI is calculated by summing the phase difference between every pair of neighboring elements from MOI. It can be seen that the phase distributions from the k x integral and MOI coincide well and they are parabolic since the k x distribution is nearly linear.

Propagation through a limited BDA size
Propagation of the mutual optical intensity through a limited BDA size is studied using the improved MOI code. The optical setup is the same as that in Section 3.1 (cf. Fig. 2). The numbers of elements at the source plane, the BDA plane and the image plane are all the same (N = 300). The intensity and in-plane wavevector profiles at the image plane for various BDA sizes are shown in Fig. 6.
The results with a BDA size of 30 mm are shown in Figs. 6(a) and 6(b). The intensity distribution [cf. Fig. 6(a)] demonstrates a Fraunhofer diffraction profile, since the BDA size is much smaller than the propagation distance. The in-plane wavevector profile [cf. Fig. 6(b)] is linear, agreeing with a spherical wave distribution. The corresponding radius of the spherical wave is 100.004 m, which matches the propagation distance r 23 = 100 m from the BDA to the image plane. This result indicates that the propagated wavefield is similar to that of a point source when the BDA size is very small.
Figs. 6(c) and 6(d) show the intensity and in-plane wavevector profiles, respectively, for a BDA size of 300 mm. The BDA acts not like a point source but like a limited line source. The in-plane wavevector profile shows oscillations at the positions of the intensity valleys owing to the diffraction of the BDA. The corresponding spherical wave radius is 106 m, only slightly larger than r 23 = 100 m, also indicating that the wavefield at the image plane is determined by the diffraction of the BDA.
Figs. 6(e) and 6( f ) show the results for a BDA size of 800 mm. The wavevector profile clearly shows two different slopes. The corresponding spherical wave radius for the center region is 207 m, which is similar to the distance r 13 = 218 m from the source to the image plane. The corresponding spherical wave radius for the two edge regions is 90.6 m, close to r 23 . The double-peak feature of the intensity profile is a combination of the Fresnel diffraction of the BDA and the Fraunhofer diffraction of the limited source.
When the BDA size is increased further to 2000 mm [cf. Figs. 6(g) and 6(h)], the spherical wave radii for the center and edge regions are 218 and 89.5 m, respectively. The center region, where most of the beam intensity resides, becomes larger. The intensity distribution is dominated by the Fraunhofer diffraction of the source, with a small contribution from the BDA diffraction.
The BDA size in Figs. 6(i) and 6( j) is 5000 mm, which is much larger than the spot size at the BDA plane ( = 570 mm). The intensity profile on the image plane can be regarded as the result of direct propagation from the source, which is similar to the case in Section 3.1. The in-plane wavevector profile still shows two slopes. In the center area where the source can be directly 'seen,' the corresponding spherical wave radius is  219 m, close to r 13 . For the 'invisible' edge regions, the corresponding spherical wave radius is 98.4 m, close to r 23 . We note that the in-plane wavevector is independent of the intensity. Even though the intensity within the edge regions is very weak, the in-plane wavevector still shows clear highfrequency oscillations caused by the BDA diffraction.
The corresponding spherical wave radii for the central region, where most of the beam intensity is located, at the image plane for different BDA sizes are summarized and shown in Fig. 7. The radius of the central region is a direct measure of the location of the effective source. One can clearly see the transition of the effective source location moving from the BDA plane to the source plane with increasing size of the BDA. All this information can easily be obtained by analyzing the in-plane wavevector.  (a, c, e, g, i) Intensity and (b, d, f, h, j) in-plane wavevector distributions at the image plane for various BDA sizes, which are (a, b) 30 mm, (c, d) 300 mm, (e, f ) 800 mm, (g, h) 2000 mm and (i, j) 5000 mm.

Analysis of the in-plane wavevector with different source coherence lengths
We further analyze the in-plane wavevector with different source coherence lengths. A simple optical setup is considered in which the beam propagates from the source to the image plane over a distance of 218 m. The size of the source is apertured to 200 mm. The source is assumed to obey the GSM with an r.m.s. size of = 100 mm. The intensity and in-plane wavevector profiles at the image plane under various coherence conditions (different coherence lengths ) are shown in Fig. 8.
When the source is fully coherent ( = 1), the in-plane wavevector shows jumps to near-zero values at the positions of those intensity valleys [cf. Fig. 8(a)] where there is a shift in the phase (Nye & Berry, 1974). When the source is partially coherent [e.g. = 3000 mm in Fig. 8(b)], the in-plane wavevector still has those jumps but is not close to zero as in Fig. 8(a). When the source coherence is even worse ( = 500 mm), the jumps disappear completely, as shown in Fig. 8(c). This phenomenon occurs because the phase at a given position is a statistical result. The effect of the phase shift is smoothed out by the partial coherence of the beam. The lower the beam coherence, the less profound are the diffraction effects that can be observed. The improved MOI code provides accurate calculation of the in-plane wavevector, which can be useful for detailed partial coherence analysis of the beamline.

Mutual optical intensity propagation through a non-ideal mirror
We calculate the mutual optical intensity propagation through a non-ideal mirror using both the original and improved MOI codes. The optical setup is shown in Fig. 9. The partially coherent source is focused on the focal plane by an elliptical cylinder mirror. The entrance and exit arm lengths of the mirror are 34 m and 8 m, respectively. Intensity (black curves) and in-plane wavevector profiles (red curves) at the image plane calculated using the improved MOI code with different source coherence lengths of (a) = 1, (b) = 3000 mm and (c) = 500 mm.

Figure 9
The optical setup for the propagation of a partially coherent beam through an elliptical cylinder mirror.

Figure 7
The corresponding spherical wave radii for the central region at the image plane for various BDA sizes.
The GSM source has parameters of = 153 mm and = 109.6 mm, and is limited to a size of 1530 mm. The X-ray wavelength is = 2.55 nm. The incident grazing angle is 1.5 . The elliptical cylindrical mirror has an r.m.s. slope error of 0.5 mrad, with the height error profile generated following the methodology of Shi et al. (2016) and shown in Fig. 10.
The intensity profiles and degree of coherence at the focal plane calculated using the original and improved MOI codes with different numbers of elements are shown in Fig. 11. As shown in Figs. 11(a) and 11(b), the intensity profiles and degree of coherence derived from the original MOI code are significantly distorted when the number of elements N < 200. The calculation results from the improved MOI code demonstrate a more stable performance with smaller N, as shown in Figs. 11(c) and 11(d). The in-plane wavevector is not zero at the off-axis. The maximum in-plane wavevector is 0.11 mm À1 in the propagation, and the corresponding divergence angle between the wavevector and element normal is 0.035 nrad. The original MOI model ignores the divergence angle within an element and assumes it to be zero. The improved MOI model considers the divergence angle within an element in the propagation and shows more accurate results. This means the divergence angle in an element is important to get a more accurate mutual optical intensity propagation with a low number of elements. The r.m.s. spot sizes at the focal plane using the original and improved MOI codes with different numbers of elements are shown in Fig. 12. Similar to the results for free-space propagation [cf. Fig. 4(a)], the improved MOI code provides much more stable calculations while requiring smaller numbers of elements. Again, improvements in both calculation accuracy and efficiency can be obtained using the improved MOI code.

Conclusions
The capability of the MOI code is extended by considering the phase distribution inside each element through the calculation of the in-plane wavevector. The developed formulae are The surface-height error profile of the elliptical cylinder mirror with an r.m.s. slope error of 0.5 mrad.  general and can be applied to any optics. The improved MOI code is compared with the original code by the propagation simulation of the mutual optical intensity through free space and non-ideal mirrors. The improved code is demonstrated to provide both higher accuracy and higher efficiency. This is an important step towards the future development of the 2D MOI code.
Knowledge of the in-plane wavevector also allows detailed analysis of the beam wavefield. Since there is no singlecoordinate phase distribution for partially coherent light, the phase information stored in the mutual optical intensity is described as the phase difference of every pair of coordinates, which is a statistical distribution. The in-plane wavevector is a single-coordinate function that describes the phase information of partially coherent light, which is more convenient and accurate than the phase difference.