Simulating dark-field X-ray microscopy images with wavefront propagation techniques

Carlsen, M.; Detlefs, C.; Yildirim, C.; Ræder, T.; Simons, H.

doi:10.1107/S205327332200866X

research papers

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 78| Part 6| November 2022| Pages 482-490

https://doi.org/10.1107/S205327332200866X

Open

access

Simulating dark-field X-ray microscopy images with wavefront propagation techniques

Mads Carlsen,^a ^* Carsten Detlefs,^b Can Yildirim,^b Trygve Ræder ^a and Hugh Simons ^a

^aDepartment of Physics, Technical University of Denmark (DTU), Fysikvej, Building 311, 2800 Kgs Lyngby, Denmark, and ^bEuropean Synchrotron Radiation Facility, 72 avenue des Martyrs, 38043 Grenoble, France
^*Correspondence e-mail: madsac@dtu.dk

Edited by I. A. Vartaniants, Deutsches Electronen-Synchrotron, Germany (Received 1 March 2022; accepted 30 August 2022; online 10 October 2022)

Dark-field X-ray microscopy is a diffraction-based synchrotron imaging technique capable of imaging defects in the bulk of extended crystalline samples. Numerical simulations are presented of image formation in such a microscope using numerical integration of the dynamical Takagi–Taupin equations and wavefront propagation. The approach is validated by comparing simulated images with experimental data from a near-perfect single crystal of diamond containing a single stacking-fault defect in the illuminated volume.

Keywords: diffraction imaging; dynamical diffraction; simulation.

1. Introduction

Dark-field X-ray microscopy (Simons et al., 2015 ) (DFXM) is a full-field X-ray imaging technique similar to X-ray topography (Berg, 1931 ; Lang, 1957 ). However, unlike the latter, DFXM utilizes an X-ray objective lens placed between the sample and the detector to create a magnified image and can therefore achieve a spatial resolution better than the detector pixel size. The spatial bandwidth is limited by the small numerical aperture (NA ≃ 10⁻³) of this objective lens. Compared with classical X-ray topography, this provides angular resolution of the scattered beam direction (Poulsen et al., 2017 ), which makes it possible to quantitatively measure strains and rotations of the crystal lattice by sequentially collecting images while rotating the sample and moving the objective lens.

Traditionally, the quantitative analysis of DFXM data and the theoretical description (Poulsen et al., 2017, 2021 ) of the method have relied on strong approximations. Most important is the kinematical approximation, which is to omit multiple scattering effects, and which holds for small and for highly deformed crystals. Another approximation, which is built into the geometric optics treatment of Bragg scattering, is that infinitesimal sub-volumes in the sample scatter according to the Bragg law for a perfect infinite crystal, and that the intensities scattered from such sub-volumes add together incoherently. In some cases (e.g. for near-perfect crystallites and small defect structures), however, such approximations are not valid.

Here we present a method for simulating DFXM images based on wavefront propagation combined with a framework that treats multiple scattering events (known as the dynamical theory of X-ray diffraction) (Authier, 2001 $[Authier, A. (2001). Dynamical Theory of X-ray Diffraction, pp. 534-551. Dordrecht: Springer Netherlands.]$ ). This method is able to handle the effects of coherence, dynamical scattering, aberrations of the objective lens and detector imperfections, and as such it is a more realistic model of DFXM for near-perfect crystals than the geometric optics model. This will be useful for understanding the type of contrast observed in DFXM images and can aid in experimental planning and data analysis.

Furthermore, a number of advanced approaches to DFXM have been suggested in the literature and tested by various authors [e.g. magnified topo-tomography (Jakobsen et al., 2019 ), confocal Bragg microscopy and tele-ptychography (Pedersen et al., 2020 ), Fourier ptychographic DFXM (Carlsen et al., 2022 )]. To date, the theoretical models used in these examples rely on idealized models of the underlying physics. Again, a more realistic forward model would be useful for testing the limits of their applicability.

In this paper we discuss the steps that are required for such a simulation and present an implementation and some initial results. We compare simulated and experimental data for a near-perfect single-crystal diamond containing only one single stacking fault in the imaged field of view (FOV).

2. Method

Recent developments in DFXM are moving towards studying the dynamics of isolated defects, such as dislocations, domain walls and acoustic waves in near-perfect single crystals (Dresselhaus-Marais et al., 2021 ; Holstad et al., 2021 ). Here, dynamical effects cannot be ignored. In these cases, the analysis of DFXM data has often relied on the weak beam approximation, where it is assumed that even highly perfect crystals scatter approximately kinematically when the sample is rotated to the tails of the rocking curve. At this position, the bulk of the crystal does not scatter strongly. Instead, only small strained volumes near defects and surfaces contribute to the scattered intensity.

Much has been published on the subject of simulating X-ray topography images (Taupin, 1964 ; Authier, 1968 ; Epelboin & Soyer, 1985 ). These methods are applicable in our case, but a few extra precautions must be taken due to our need for high quantitative accuracy, and the added complication caused by the objective lens. There exist useful solutions (Chubar et al., 2013 ; Pedersen et al., 2018 ; Celestre et al., 2020 ) for simulating synchrotron sources and X-ray optical components of the kinds applied in a DFXM instrument with coherent wavefront methods. A full simulation based on coherent wavefronts is therefore within reach by combining these established methods.

In dark-field transmission electron microscopy (DFTEM), simulating images with dynamical diffraction effects is done routinely. Although superficially similar to DFXM, the methods cannot be transferred from one technique to the other. The most common methods in DFTEM are either Bloch-wave methods that rely on the column approximation which cannot be applied to DFXM where the scattering angles often exceed 10°, or multi-slice methods that require the atomic structure of the sample to be well sampled, which is not feasible for DFXM where sample sizes often exceed 100 µm (Kirkland, 1998 ).

X-ray topography methods, on the other hand, rely on the two-beam approximation (Taupin, 1964), which cannot be applied to electron microscopy where potentially hundreds of reflections contribute significantly to the images.

Here we apply the formalism of the Takagi–Taupin equations, which make use of the two-beam approximation to simplify the scattering problem. This allows us to treat Bragg scattering in an averaged way, removing the requirement to over-sample the unit cell. Instead, the sampling is limited by the divergence of the incident X-rays and the size of features in the sample. Propagation from the sample to the detector is handled by established paraxial Fourier optics methods.

We split the simulation flow into five steps as follows:

(i) Beam: calculate the amplitudes of the modulated waves of the beam incident on the sample.

(ii) Sample: calculate the complex scattering function throughout the sample crystal.

(iii) Integration: numerically integrate the Takagi–Taupin equations to get the complex amplitudes of the scattered beam on the exit surface of the sample.

(iv) Propagation: propagate the scattered complex amplitudes through the imaging optics to the detector.

(v) Detector: interpolate the propagated field to the detector pixels and account for detector characteristics.

A schematic drawing of the steps and flow of the simulation is shown in Fig. 1.

Figure 1
Schematic drawing illustrating the different steps of the simulation flow.

2.1. Defining the beam

We work in a paraxial wave-optics formalism and use a coherent mode decomposition to describe the state of the X-ray beam at a given plane. This is given by a number of modes, labeled with p: E^(p)(x, y; z), where E is the amplitude function of a modulated plane wave: $[{\bf E} = \Re[E({\bf r})\exp(i{\bf k}_{0}\cdot{\bf r}-i\omega t)]\, \hat{{\bf p}}]$ , where r = (x, y, z) is the spatial position, k₀ is the modes' wavevector, ℏω is the photon energy and $[\hat{{\bf p}}]$ is the polarization vector.

If a mode is known on a plane z = 0, then the mode's amplitude on another plane can be found by applying a (linear) coherent wavefront propagator:

$[E^{(p)}(x,y\semi z) = {\cal P}_{0\rightarrow z}[E^{(p)}(x,y\semi 0)].\eqno(1)]$

We assume that the radiation is beam-like, i.e. that it has a limited extent in two dimensions in real space and in all dimensions in reciprocal space. The reciprocal-space distribution is around a central wavevector k₀, which defines the nominal direction and wavelength of the incident beam.

2.2. Defining the sample

In the Takagi–Taupin (Takagi, 1962 , 1969 ; Taupin, 1964) approach to X-ray scattering in the two-beam approximation, the only quantities of interest are the average electric susceptibility of the crystal, χ₀, and the two scattering constants χ_h and $[\chi_{\overline{h}}]$ that describe the cross section for scattering and back-scattering, respectively. These may be spatially modified by a displacement field u(r) of the strained crystal:

$[\chi^{\prime}_{h}({\bf r}) = \exp[-i{\bf Q}\cdot{\bf u}({\bf r})]\chi_{h}\eqno(2)]$

$[\chi^{\prime}_{\overline{h}}({\bf r}) = \exp[i{\bf Q}\cdot {\bf u}({\bf r})]\chi_{\overline{h}}.\eqno(3)]$

Here Q is the reciprocal-lattice vector of the given reflection in an undeformed reference lattice. This dependence on the displacement field explains the high sensitivity to small strains. χ_h and $[\chi_{\overline{h}}]$ can often be considered constants, but in samples with twin boundaries there may be discontinuous jumps in the values of this function, which explains contrast observed at presumably strain-free inversion twin domain boundaries in polar materials (Klapper, 1987 ).

In this paper we focus on slab-shaped single crystals, i.e. crystals with two parallel surfaces and infinite size in the orthogonal directions (see Fig. 2). The normal of the exit surface is called $[\hat{{\bf n}}]$ .

Figure 2
Scattering geometry and simulated box in the case where the scattering plane is normal to $[\hat{{\bf y}}]$ . The two angles α₀ and α_h fully specify the scattering geometry in this case.

We utilize a discrete representation of the sample structure on an orthogonal grid defined by the three directions $[\hat{{\bf x}}]$ , $[\hat{{\bf y}}]$ and $[\hat{{\bf z}} = \hat{{\bf n}}]$ , and corresponding step sizes d_x, d_y and d_z. It is necessary that the surface normal of the crystal slab is parallel to the $[\hat{{\bf z}}]$ axis, but we make no requirement on the last free rotation. The number of grid points in each direction will be labeled N_x, N_y and N_z, giving a total size of the simulated domain of L_x = d_x(N_x − 1), L_y = d_y(N_y − 1), L = d_z(N_z − 1). L is the thickness of the simulated crystal slab.

The complex value of these scattering functions needs to be known with high resolution. For simple test cases, where the displacement field is given by an analytical expression, this is not a point of concern. If, however, the displacement field is generated by a numerical simulation, then it needs to be computed with sufficient resolution to at least match the resolution of the experiment (≃50 nm) throughout the volume of the sample. The size of the sample can be up to several hundred micrometres.

For highly deformed samples, the scattering function contains a phase factor of the shape exp(iQ · u). In order to limit the phase variation between adjacent voxels to less than 2π, step sizes must be below |∇u|/(2π|Q|), which means small steps must be used for highly deformed samples. In samples with nanometre-sized domains or other small structures, all structural features must be resolved.

2.3. Integrating the Takagi–Taupin equations

Scattering of X-rays by a deformed crystal is treated in the formalism of the Takagi–Taupin equations (TTEs). In formulating the TTEs, there is an arbitrary choice of the vector k₀, which leads to slightly different versions of the final equations. The ones chosen here are referred to as the symmetrical TTEs (Vartanyants & Robinson, 2001 ) that arise when k₀ is taken to be the vacuum wavevector of the incident beam as opposed to the more conventional choice of the refracted wavevector used in the original publications:

$[\eqalignno{ 2\left({\bf k}_{0}\cdot{\boldnabla}\right)E_{ 0}^{(p)}& = -ik^{2}\chi_{0}E_{0}^{(p)}-ik^{2}C\chi_{\overline{h }}^{\prime}E_{h}^{{\prime}(p)}&\cr 2\left({\bf k}_{h}\cdot{\boldnabla}\right)E_{h}^{(p)}& = -ik^{2}(\chi_{0}+\beta)E_{h}^{(p)}-ik^{2}C\chi_{h}^{\prime}E _{0}^{(p)},&(4)}]$

where $[\beta = 2\sin(2\theta)\phi]$ is a measure of the misorientation away from the vacuum Bragg condition, ϕ is the rocking angle and C is the polarization factor. E₀^(p) are the modes of the plane-wave decomposition of the incident beam and E_p^(p) are the corresponding modes of the scattered beam, which is given relative to the wavevector k_h = k₀ + Q.

The TTEs are typically solved using a finite-difference integration scheme. A number of different algorithms have been used in the literature with the pioneering work done by Taupin (1964), Authier (1968), Epelboin & Soyer (1985). In this paper we use a novel, more flexible scheme which is documented by Carlsen & Simons (2021 ). In contrast to other established methods, this method is able to operate on an orthogonal grid representation of the sample structure, which is achieved by implicitly utilizing Fourier interpolation. This, however, causes a less efficient use of computer resources.

The geometry of the problem is set by the shape of the incident beam, the vectors k₀ and Q, as well as the choice of a computational grid. The plane spanned by the two vectors k₀ and Q is called the scattering plane.

When the scattering plane is normal to the $[\hat{{\bf y}}]$ axis, the TTEs decouple into a set of 2D problems that can be solved slice-by-slice. If we are free to choose the orientation of the computational grid, we can always choose a geometry where this becomes the case.

When Q is parallel to the surface of the crystal, we say that we have a symmetric Laue geometry. In this case we can choose an orientation of the computational grid where $[{\bf Q}||\hat{{\bf x}}]$ , where the TTEs take a particularly simple form.

In order to solve the TTEs, we impose zero Dirichlet boundary conditions in the two transverse dimensions, x and y. These boundary conditions require the sample grid to be large enough to fit the entire Borrmann triangle extending from every point where the incident beam amplitude is non-zero. This is fulfilled if the non-zero part of the amplitude is fully contained in the rectangle defined by (see Fig. 2)

$[\eqalignno{&\max(0, Lk_{0x}/k,Lk_{hx}/k)\,\lt\,x\,\lt&\cr &\quad\min(L_{x},L_{x}+Lk_{0x}/k,L_{x}+Lk_{hx}/k)&\cr &\max(0, Lk_{0y}/k,Lk_{hy}/k)\,\lt\,y\,\lt&\cr &\quad\min(L_{y},L_{y}+Lk_{0y}/k,L_{y}+Lk_{hy}/k),&(5)}]$

where L_x, L_y are the lengths of the simulated domain in the x and y directions, and L is the thickness of the crystal. For a given incident beam and Q, this sets a minimum on the size of the simulated domain in the two transverse directions.

With the sample structure and the boundary conditions given, the TTEs constitute an initial value problem, where the initial value is the amplitude of the incident X-rays on the z = 0 surface. This can be solved by an appropriate finite-difference scheme to yield the amplitudes of both the transmitted and scattered beams on the exit surface z = L.

2.4. Propagating through the optics

In DFXM, the scattered beam at the exit surface of the sample is imaged onto a detector using an objective lens in a magnifying geometry. The propagation through the lens is challenging to simulate due to the inherent thick-lens behavior of the compound refractive lens (CRL) which is typically used as the objective lens. For this, we use a computational approach where each lens in the CRL is treated as a thin lens and a paraxial FFT (fast Fourier transform) propagator is used to propagate the wavefront between each lens.

To do this, we make use of a method for propagating wavefronts with rapidly varying quadratic phases that come from the transmission function of the lenses that treats the quadratic component analytically (Ozaktas et al., 1996 ; Chubar & Celestre, 2019 ). To simulate an aberrated CRL, a separate aberration function is included at each lenslet in the CRL. These methods have previously been used to model thick CRLs like the one used in this study (Pedersen et al., 2018; Celestre et al., 2020).

It is useful to introduce a new optical axis ( $[\hat{{\bf z}}_{i}]$ ) aligned with the average wavevector of the scattered beam, k_h. To bridge the gap between these two coordinate systems, we project the values from the exit surface of the crystal to the z = 0 plane of the new coordinate system.

The inherent near-field nature of the imaging geometry, which is determined by the fact that the FOV is as large as the aperture of the objective lens, is handled by first multiplying the field with the near-field phase factors:

$[\phi_{{\rm n.f.}} = \exp\left[{{i\pi(x_{i}^{2}+y_{i}^{2})} \over {\lambda d_{1}}}\right].\eqno(6)]$

The need to over-sample this function can effectively limit the size of the FOV that is possible to simulate for a given pixel size.

2.5. Detector characteristics

With the mode amplitudes on the detector plane given, we now need to interpolate these values to the detector pixels and incoherently sum over the modes of the incident beam.

If the purpose of the simulation is to estimate the resolution or to create a data set to be used to test the data analysis procedures, one should remember to include non-ideal behavior introduced at the detector. The most important effects are the incoherent point spread of the detector, background signal, counting noise and non-linear response.

The detector used here is an indirect detector composed of a scintillator crystal (25 µm-thick gadolinium gallium garnet), optical microscope (Mitutoyu M Plan Apo 10×, NA = 0.28 and tube lens) and pco.edge 5.5 sCMOS camera with pixel size 6.5 µm. Contributions to the incoherent point spread function arise from the diffraction limit due to the finite NA of the optical microscope, the finite thickness of the scintillator (especially when it is larger than the depth-of-focus of the optical microscope), scattering within the scintillator and the optical microscope, and aberrations.

Counting statistics/readout noise can be the critical factor when imaging small grains or weak reflections. The non-linear response (saturation) might be quite important for perfect crystals, as the images have interesting features over a very large dynamic range.

3. Comparison with experiment

To test our approach, we simulate a section-topography-type experiment with a near-perfect single-crystal diamond slab of thickness 300 µm containing a single stacking-fault defect. The sample is imaged in a symmetric Laue geometry in a {111} reflection with [110] entrance and exit surfaces.

The investigated defect is a stacking fault, which arises by the addition or removal of a single close-packed plane of atoms in the face-centered cubic (f.c.c.) parent lattice of the diamond. The fault vector is of the family $[{\bf b}_{{\rm s.f.}} = ({1} / {3})\langle 111\rangle]$ which is not a translational symmetry of the f.c.c. lattice. These planar defects are bounded by the surfaces of the crystal and Frank-type partial dislocations (Frank, 1951 ; Kowalski et al., 1989 ). In the described experiment, the edges of this defect lie outside the Borrmann triangle and the effects of the strain originating at the edges can be ignored. This allows us to treat the stacking fault as an infinite planar defect [on one side u(r) = 0, on the other u(r) = b_s.f.]. In the Takagi–Taupin description, the stacking fault thus becomes a discrete jump in the phase of the scattering function of magnitude 2π[hkℓ] × b_s.f. = ±2/3π when imaged using a [111] reflection not orthogonal to the stacking-fault normal (Klapper, 1987). A sketch of the sample geometry is shown in Fig. 3.

Figure 3
(a) 3D sketch of the simulated geometry where the blue sheet represents the incident/transmitted X-ray beam, the red represents the scattered rays and the yellow plane represents the stacking fault. (b) 2D slice of the simulated geometry showing the same features as (a) with definitions of useful angles and distances.

This constitutes a good test sample, as the defect (aside from the 12 possible orientations of the defect) only has a single continuous degree of freedom: the position of the stacking fault along the normal. Furthermore, diamond is one of a few materials where macroscopic crystals with very low defect density are available. Provided our method captures all the relevant physics, we should therefore be able to perfectly recreate the experimental data.

The dynamical scattering patterns produced by isolated stacking faults in diamonds have previously been studied in detail by classical X-ray topography (Kowalski et al., 1989).

Experiments were carried out at the ESRF dark-field X-ray microscopy beamline, ID06-HXM (Kutsal et al., 2019 ), after the EBS (Extremely Brilliant Source) upgrade. A Si(111) Bragg–Bragg double-crystal monochromator selected X-rays with photon energy 17 keV from the undulator source. The spot size of the beam on the condenser lens is limited by a slit in the vertical direction to 0.2 mm. Based on parameters of the X-ray source given by the ESRF, we estimate a vertical coherence length of 340 µm at the distance of 50 m from the source to the condenser lens, which is significantly larger than the aperture of the condenser lens. The horizontal coherence length is estimated to be 40 µm, which is smaller than the FOV but much larger than the extent of the coherent point spread function of the objective lens. The objective lens consists of 70 individual biconcave Be lenses of apex curvature 50 µm. The CRL is further modified with a 100 µm square aperture after the last lens

To describe the incident X-rays, we use only one coherent mode for each wavelength and 41 different wavelengths covering a relative energy spread of 6 × 10⁻⁴ in total. For each energy component, the amplitude at the sample position is the Fourier transform of the complex transmission function of an aberration-free 1D condenser lens with a Gaussian-shaped absorbing aperture and a hard cut-off at 140 µrad. Consequently, the angular spectrum in the transverse direction is a top-hat profile with a small Gaussian smoothing at the edges, whereas the spatial profile is a diffraction-limited focal spot with some ringing due to the hard cut-off (see Fig. 4). We omit the use of several transverse mutually incoherent modes. The effects of including these would be small, since the effective source size is small compared with the NA of the objective lens.

Figure 4
Vertical profile in (a) reciprocal and (b) real space of the beam used in the simulation of this paper.

In the experimental realization we make the observations in a laboratory frame defined by $[\hat{{\bf x}}_{{\rm lab}}]$ , $[\hat{{\bf y}}_{{\rm lab}}]$ and $[\hat{{\bf z}}_{{\rm lab}}]$ , where $[\hat{{\bf z}}]$ is parallel to k₀. In this experiment the scattering is vertical in the laboratory frame, meaning that $[\hat{{\bf y}}_{{\rm lab}}]$ is normal to the scattering plane. We choose the integration grid such that $[\hat{{\bf y}} = \hat{{\bf y}}_{{\rm lab}}]$ .

The simulation used a grid of 2048 × 2048 × 3001 points with step sizes of 60, 60 and 100 nm, respectively. This matches the 300 µm thickness of the sample and the ≃100 µm FOV in the y direction. The large size in the x direction was necessary to satisfy the constraints of equation (5). Step sizes in x and y directions are chosen to be larger than the resolution of the final images, and the step size in z is chosen such that the integration error of the finite-difference scheme is lower than the noise level of the final images. The execution time is dominated by the integration of the TTEs, which took 3.5 h per energy point running on a single core. The simulation was parallelized over the modes.

The computation time was dominated by the integration of the TTEs, which involved ≃10 000 1D Fourier transforms of the 2048 × 2048 arrays used in this example compared with just ≃150 2D Fourier transforms for the simulation of the optics. The use of computer resources scales linearly with the number of modes, so if a transversely incoherent model was to be used, the needed resources would increase by a factor of the number of modes. In the study presented here computation time was not an issue, but if high-throughput simulation was needed, large performance gains could likely be achieved by applying a more efficient algorithm for the integration of the TTEs (Epelboin & Soyer, 1985). The integration of the TTEs can also be parallelized into separate 2D problems for each vertical slice of the sample and should easily scale with available computer resources. If the thick-lens effects are not critically important, the CRL may well be approximated by a single convolution step (Pedersen et al., 2018). Furthermore, the large number of energy modes used in this example was only necessary to ensure proper sampling of the high-frequency features used as lens aberration and avoid aliasing problems in Fig. 9. If these are not important, a much smaller number of energy points would be sufficient.

3.1. Near-field measurements

The microscope used for the experiment provides the possibility to additionally carry out traditional X-ray topography by placing a detector closely behind the sample (40 mm in the examples shown here) without using the objective lens. This is useful for alignment and for low-resolution characterization of the sample.

For a single coherent mode [Fig. 5(b)] this propagation distance causes recognizable Fresnel-diffraction fringes around sharp features in the scattered field. In a polychromatic simulation [Fig. 5(c)], these fringes are blurred out. The difference in wavelength (which causes a small difference in the free-space propagator) is not sufficient to explain this blurring. Rather, the broadening is caused by the slight difference in scattering angle of the different energy components in the incident beam [Fig. 5(d)]. The features in the simulated image, however, are not as wide as in the measured data [Fig. 5(a)]. This is likely explained by the incoherent point spread of the detector.

Figure 5
X-ray section topography of the single diamond (110) sample containing a single stacking-fault defect in the imaged volume, imaged at a propagation distance of 40 mm. (a) Experimental image, (b) simulated image with one coherent mode, (c) polychromatic simulated image. The dotted white lines in (a)–(c) mark the edge of the Borrmann triangle. (d) Sketch of the scattering geometry in reciprocal space in the laboratory frame. The incident beam consists of a continuum of rays of different k₀ vectors (blue rectangle); due to the divergence and finite bandwidth of the incident beam, this smears out the Ewald's sphere to a shell. For a given rotation of the crystal only points on the dotted line (the bisection of the Q vector) satisfy the Bragg condition. Different energy components scatter at slightly different angles.

The vertical divergence of the condensed line beam is large compared with the intrinsic width of the dynamical rocking curve of the diamond sample. The crystal therefore acts as an analyzer when rocked in the condensed line beam (see Fig. 6). The width of the rocking curve is largely determined by the divergence of the incident beam, but the finite-energy bandwidth blurs the sharp cut-off caused by the condenser lens' physical aperture. The asymmetry of the measured rocking curve [Fig. 6(c)] reveals a misalignment of the condenser lens.

Figure 6
(a) Simulated spatially resolved rocking curve of a slice of the diamond crystal away from the stacking fault using an aberrated model of the condenser lens. (b) Measured spatially resolved rocking curve of a slice of the diamond crystal away from the stacking fault. (c) Measured and simulated spatially integrated rocking curves of a diamond single crystal in the region containing a single stacking fault. (d) Sketch of the scattering geometry in reciprocal space in the laboratory frame. The incident beam consists of a continuum of rays of different k₀ vectors (blue rectangle); due to the divergence and finite bandwidth of the incident beam, this smears out the Ewald's sphere to a shell. For a given rotation of the crystal only points on the dotted line (the bisection of the Q vector) satisfy the Bragg condition. When the crystal is rotated by an angle ϕ (thereby rotating Q into Q′), this line is moved and a different part of the incident spectrum satisfies the Bragg condition. This leads to scattering in a different direction in the laboratory frame (comparing red and pink arrows).

High-frequency defects in the condenser lens are visible in the spatially resolved rocking curve of perfect parts of the crystal [see Fig. 6(b)] as vertical stripes. The stripes are blurred along the rocking-angle direction due to the finite bandwidth of the incident X-rays – an effect which is confirmed by the simulations shown in Fig. 6(a).

When the crystal is rocked, a different part of the spectrum of the incident beam will be in the Bragg condition and therefore the sample will scatter in slightly different directions as a function of the rocking angle ϕ [sketched in Fig. 6(d)]. Due to the finite propagation distance from the sample to the detector, this change in angle translates to a change in position of the measured image. This mix of position and angular information, which is avoided by the use of an objective lens, is unavoidable in X-ray topography methods due to the finite propagation length, but it is exaggerated in this study due to the relatively large propagation distance and large vertical divergence compared with more usual topography experiments.

3.2. DFXM images

To simulate the DFXM images, we use a model of an ideal CRL. The physical aperture of the CRL is defined by a 0.1 × 0.1 mm square absorbing slit placed at the exit of the lens. The only fitted parameters for the whole simulation are the relative positions of the sample, lens and detector as well as the noise level of the detector. The intensity of the simulated images is scaled to match the intensity of the measured images. The same scaling parameter is used for all images. The experimental images are overexposed (saturated) at the direct image of the stacking fault; therefore we here choose a color map that clips the highest intensities in the simulated images.

As can be seen in Fig. 7, the DFXM simulation qualitatively recreates the features of the experimental images. However, there are a number of deviations:

Figure 7
Measured (a), (c) and simulated (b), (d) DFXM images of a stacking-fault defect in a diamond single crystal. (a), (b) show images where the objective lens is placed in the center of the diffracted beam. (c), (d) show images where the lens is displaced from the diffracted beam leaving the bottom half of the FOV in the dark field. The dotted white line marks the edge of the Borrmann triangle.

(i) We underestimate the magnification of the imaging setup by about 5%, which leads to an incorrect scaling of the images. This is most likely due to a small deviation of the apex radius of curvature from the nominal value of 50 µm in the individual lenses of the CRL used as objective lens.

(ii) The simulated images contain a smaller number of Borrmann fringes (the horizontal stripe features) than the measured data. We attribute this to the known high sensitivity of the spacing of Borrmann fringes to small macroscopic strain gradients (Rodriguez-Fernandez et al., 2021 ). Alterntively, a slight miscalibration of the photon energy or sample thickness could also lead to this difference.

(iii) The simulated images have a regular pattern of vertical streaks close to the right-hand side of the images. These are due to Fresnel diffraction from the hard edge of the square aperture in the objective lens. This is likely an artifact of the assumption of perfect transverse coherence in the horizontal direction or of the perfectly sharp edges of the aperture that are somewhat jagged in practice.

(iv) The measured images contain noise with the appearance of vertical streaks and speckle-like features close to the brightest features. This can be explained either by the aberrations in the condenser lens or in the objective lens, as will be discussed later.

In Figs. 7(c), 7(d) we simulate an image where the objective lens is displaced from the center of the scattered beam such that rays that are specularly reflected fall outside the aperture of the objective lens in the bottom part of the displayed ROI (region of interest). In that region, only diffusely scattered X-rays will contribute to the image. This results in the disappearance of the dynamical features, while the direct image of the stacking fault can still be seen. In visible light microscopy, this is referred to as `dark-field contrast'.

In the geometric optics treatment, weak beam contrast is explained by the presence of small regions where the lattice is strained and rotated away from the average lattice. In these regions, rays are scattered if they satisfy the exact Bragg condition for the deformed lattice. A stacking fault is in principle a perfectly sharp defect with no spatial extent so no such region exists. The appearance of weak beam contrast therefore illustrates the inability of the geometric model to handle diffraction effects that are important when describing scattering from small structures with a characteristic size on the order of the wavelength or smaller.

Since the stacking fault is thought to be a perfectly sharp defect, the width of the image of this defect can serve as a rough estimate of the resolution of the instrument. The stacking fault is a 2D feature, and therefore the width of the image is not only determined by the resolution of the imaging optics, but also by the projection of the part of the stacking fault illuminated by the sheet beam along the scattered beam direction.

In Fig. 8 (inset) we compare the width of this feature in the experimental images with that in the simulated images. We see that the polychromaticity does not contribute significantly to the width of the feature in the simulations. Previous studies of the chromatic aberrations in CRLs, using the same computational approach as we apply here, also find that the chromatic aberration only adds a small part to the point spread of the imaging optics (Pedersen et al., 2018).

Figure 8
Comparison of the width of the dark-field image of the stacking fault as seen in Figs. 7

(c), 7

(d). The x axis refers to distances in the sample plane coordinates.

The experimental image has an FWHM of 1.3 µm, about 0.5 µm wider in the demagnified sample plane coordinates than that predicted by the simulations. We believe that the resolution of the experiment is degraded by aberrations in the CRL lenses.

3.3. Aberrated lenses

So far, we have ignored the effect of the aberrations in the lenses. In transmission images (Lyatun et al., 2020 ) and Bragg images taken without the condenser lens, short-wavelength aberrations are known to cause strong speckle-like noise in the final images. The apparent absence of this noise in DFXM images is surprising at first. However, as previously observed, this noise is averaged out when the imaged field is only partially coherent (Falch et al., 2019 ; Carlsen et al., 2022). Normally we think of the dynamically scattered X-rays as highly coherent as the Bragg scattering effectively collimates the incident radiation, but this argument does not consider the polychromaticity of the incident radiation.

While CRLs have been shown to be nearly achromatic over the bandwidth of the monochromated beam (Pedersen et al., 2018), Bragg scattering is not: a higher-energy component of the incident beam scatters at a smaller angle and vice versa, as sketched in Fig. 5(d). Since the incident beam has a large vertical divergence (compared with both the energy bandwidth and the Darwin width) set by the aperture of the condenser lens, the integration over energies corresponds to integrating over a small spread of angles of the scattered beam. This integration averages out the high-frequency parts of the aberrations in the vertical direction, leaving features elongated in the vertical direction. A relative energy spread of 1.0 × 10⁻⁴ corresponds to an angular difference of 1.0 × 10⁻⁴ rad × tan θ which gives 8 µm at the lens plane – comparable with the average grain size (15 µm) of the O30H-grade beryllium (Lyatun et al., 2020) used in our CRL.

This effect is demonstrated in Figs. 9(a), 9(b), where an aberrated lens is constructed by multiplying the wavefront by an aberration function at the position of the first lens and the last lens in the CRL. The aberration functions used here are pure phase objects and are made by randomly placing a number of circles of random size. The amplitude and number of circles are chosen to make the simulated and measured images similar. The partially averaged speckle noise has the appearance of vertical stripes and is qualitatively similar to that observed in the real images [Fig. 7(a)]. In a typical experiment, we do not acquire sufficient data to uniquely determine the aberrations, but an effective aberration function can be recovered using Fourier ptychography (Carlsen et al., 2022). The qualitative similarity between simulated and measured images confirms that the vertical stripe artifacts seen in DFXM images of highly perfect crystals can be explained by high-frequency errors in the objective lens, which we know to be present.

Figure 9
(a) Monochromatic and (b) polychromatic simulation using an aberrated objective lens. (c) Monochromatic and (d) polychromatic simulation using an aberrated condenser lens. The scale bar refers to distances on the detector.

In Figs. 9(c), 9(d) we investigate the effect of similar aberrations, to those used in the objective lens, in the condenser lens. Once again, averaging over the energy bandwidth significantly reduces the strength of the noise and results in vertical stripes. In this case the stripes are unbroken and can be followed from the top to the bottom of the image, in contrast to the noise observed in the real images and with an aberrated objective lens.

4. Conclusion

DFXM is based on well known physics and we can predict the images it will produce – if we know the structure of the sample. It is possible to simulate the full FOV of the prototypical DFXM instrument at ID06-HXM at the ESRF.

Comparison of our simulations and experimental findings from a near-perfect single-crystal diamond suggests that most deviations between our simulations and the observed images can be explained by non-ideal behavior of the lenses. This suggests that the performance of DFXM instruments is critically limited by the quality of the objective lens.

In general, we do not have a sufficiently accurate model of the sample structure to do full simulations of the DFXM experiment, and the data collected in a typical experiment are not enough to build a full 3D model of the sample at sufficiently high resolution. Nevertheless, simulations like the ones shown here should prove useful for evaluating possible upgrades of the instrumentation and to qualitatively study the type of contrast observed from different types of defects in near-perfect single crystals, such as isolated dislocations, twin boundaries and point defects.

More deformed crystals are difficult to simulate, both because models of the displacement field in such crystals are not easily obtainable and because the large strain would require impractical, small step sizes to ensure proper sampling of the scattering function. In these samples, dynamical diffraction effects are not thought to be important and the speckle-like noise due to lens aberrations should also be less strong, as the scattering is more diffuse. So a wavefront-based simulation approach is less appropriate in this type of sample. It may however be interesting to investigate the transition from the dynamical patterns from near-perfect crystallites to kinematical scattering from deformed crystals using our new approach.

Acknowledgements

The authors acknowledge the provision of beamtime from the European Synchrotron.

Funding information

MC and HS acknowledge funding from ERC Starting Grant No. 804665.

References

Authier, A. (1968). Phys. Status Solidi B, 27, 77–93. CrossRef Web of Science Google Scholar
Authier, A. (2001). Dynamical Theory of X-ray Diffraction, pp. 534–551. Dordrecht: Springer Netherlands. Google Scholar
Berg, W. (1931). Naturwissenschaften, 19, 391–396. CrossRef CAS Google Scholar
Carlsen, M., Raeder, T. M., Yildirim, C., Rodriguez-Lamas, R., Detlefs, C. & Simons, H. (2022). Opt. Express, 30, 2949–2962. Web of Science CrossRef PubMed Google Scholar
Carlsen, M. & Simons, H. (2021). arXiv:2106.12412. Google Scholar
Celestre, R., Chubar, O., Roth, T., del Rio, M. S. & Barrett, R. (2020). Advances in Computational Methods for X-ray Optics V, edited by O. Chubar & K. Sawhney, Vol. 11493, pp. 88–101. International Society for Optics and Photonics, SPIE. Google Scholar
Chubar, O. & Celestre, R. (2019). Opt. Express, 27, 28750–28759. Web of Science CrossRef PubMed Google Scholar
Chubar, O., Fluerasu, A., Berman, L., Kaznatcheev, K. & Wiegart, L. (2013). J. Phys. Conf. Ser. 425, 162001. CrossRef Google Scholar
Dresselhaus-Marais, L. E., Winther, G., Howard, M., Gonzalez, A., Breckling, S. R., Yildirim, C., Cook, P. K., Kutsal, M., Simons, H., Detlefs, C., Eggert, J. H. & Poulsen, H. F. (2021). Sci. Adv. 7, eabe8311. Web of Science PubMed Google Scholar
Epelboin, Y. & Soyer, A. (1985). Acta Cryst. A41, 67–72. CrossRef CAS Web of Science IUCr Journals Google Scholar
Falch, K. V., Detlefs, C., Christensen, M. S., Paganin, D. & Mathiesen, R. (2019). Opt. Express, 27, 20311–20322. Web of Science CrossRef CAS PubMed Google Scholar
Frank, F. C. (1951). London, Edinb. Dubl. Philos. Mag. J. Sci. 42, 809–819. Google Scholar
Holstad, T. S., Raeder, T. M., Carlsen, M. A., Knudsen, E. B., Dresselhaus-Marais, L., Haldrup, K., Simons, H., Nielsen, M. M. & Poulsen, H. F. (2021). arXiv:2111.01545. Google Scholar
Jakobsen, A. C., Simons, H., Ludwig, W., Yildirim, C., Leemreize, H., Porz, L., Detlefs, C. & Poulsen, H. F. (2019). J. Appl. Cryst. 52, 122–132. Web of Science CrossRef CAS IUCr Journals Google Scholar
Kirkland, J. (1998). Advanced Computing in Electron Microscopy. New York: Springer US. Google Scholar
Klapper, H. (1987). Prog. Cryst. Growth Charact. 14, 367–401. CrossRef CAS Web of Science Google Scholar
Kowalski, G., Lang, A. R., Makepeace, A. P. W. & Moore, M. (1989). J. Appl. Cryst. 22, 410–430. CrossRef CAS Web of Science IUCr Journals Google Scholar
Kutsal, M., Bernard, P., Berruyer, G., Cook, P. K., Hino, R., Jakobsen, A. C., Ludwig, W., Ormstrup, J., Roth, T., Simons, H., Smets, K., Sierra, J. X., Wade, J., Wattecamps, P., Yildirim, C., Poulsen, H. F. & Detlefs, C. (2019). IOP Conf. Ser.: Mater. Sci. Eng. 580, 012007. Google Scholar
Lang, A. R. (1957). Acta Metall. 5, 358–364. CrossRef CAS Web of Science Google Scholar
Lyatun, I., Ershov, P., Snigireva, I. & Snigirev, A. (2020). J. Synchrotron Rad. 27, 44–50. Web of Science CrossRef CAS IUCr Journals Google Scholar
Ozaktas, H., Arikan, O., Kutay, M. & Bozdagt, G. (1996). IEEE Trans. Signal Process. 44, 2141–2150. CrossRef Web of Science Google Scholar
Pedersen, A. F., Chamard, V. & Poulsen, H. F. (2020). Opt. Express, 28, 15770–15782. Web of Science CrossRef PubMed Google Scholar
Pedersen, A. F., Simons, H., Detlefs, C. & Poulsen, H. F. (2018). J. Synchrotron Rad. 25, 717–728. Web of Science CrossRef IUCr Journals Google Scholar
Poulsen, H. F., Dresselhaus-Marais, L. E., Carlsen, M. A., Detlefs, C. & Winther, G. (2021). J. Appl. Cryst. 54, 1555–1571. Web of Science CrossRef CAS IUCr Journals Google Scholar
Poulsen, H. F., Jakobsen, A. C., Simons, H., Ahl, S. R., Cook, P. K. & Detlefs, C. (2017). J. Appl. Cryst. 50, 1441–1456. Web of Science CrossRef CAS IUCr Journals Google Scholar
Rodriguez-Fernandez, A., Diaz, A., Iyer, A. H. S., Verezhak, M., Wakonig, K., Colliander, M. H. & Carbone, D. (2021). Phys. Rev. Lett. 127, 157402. Web of Science PubMed Google Scholar
Simons, H., King, A., Ludwig, W., Detlefs, C., Pantleon, W., Schmidt, S., Stöhr, F., Snigireva, I., Snigirev, A. & Poulsen, H. F. (2015). Nat. Commun. 6, 6098. Google Scholar
Takagi, S. (1962). Acta Cryst. 15, 1311–1312. CrossRef CAS IUCr Journals Web of Science Google Scholar
Takagi, S. (1969). J. Phys. Soc. Jpn, 26, 1239–1253. CrossRef CAS Web of Science Google Scholar
Taupin, D. (1964). Bull. Soc. Fr. Mineral. Cristallogr. 87, 469–511. CAS Google Scholar
Vartanyants, I. A. & Robinson, I. K. (2001). J. Phys. Condens. Matter, 13, 10593–10611. Web of Science CrossRef CAS Google Scholar