Simulation of Bragg coherent diffraction imaging

The arrangement of atoms within a crystal and information on deviations from the ideal lattice is encoded in the diffraction pattern obtained from an appropriately conducted Bragg coherent diffraction imaging (BCDI) experiment. A foreknowledge of how specific displacements of atoms within the unit cell alter the BCDI diffraction pattern and the subsequent real-space image is often useful for interpretation and can provide valuable insight for materials design. Here we report on an atomistic approach to efficiently simulate BCDI diffraction patterns by factorising and eliminating certain redundancies in the conventional approach. Our method is able to reduce the computation time by several orders of magnitude without compromising the recovered phase information and therefore enables feasible atomistic simulations on nanoscale crystals with arbitrary lattice distortions.


Introduction
Displacements of atoms from their ideal position in the crystal structure of a quantum material can result in dramatic changes in a range of properties and the emergence of entirely new phases of the material that are often of considerable technological utility. Novel electronic and magnetic phases can arise as a result of subtle atomic displacements of a single atom in each unit cell and as a result, there is considerable interest in developing robust means to control displacements for use in technological applications that require on demand switching between phases [1][2][3][4][5][6].
Bragg coherent x-ray diffraction imaging (BCDI) is a lensless imaging technique that encodes information on the displacement of atoms with Angstrom sensitivity in the associated far-field diffraction pattern [7][8][9][10][11]. BCDI is performed by illuminating an isolated crystal with a spatially coherent x-ray source that exceeds the dimensions of the crystal. A crystal of the order of a few hundred nanometres is usually required for this condition to be met. The scattered light from the entire volume of the crystal then interferes in the far-field to produce a diffraction pattern. The diffraction pattern is then recorded using a pixelated photon counting detector that records the scattering intensity. The third dimension of the diffraction pattern is then obtained by rotating the crystal through the Bragg condition while maintaining a largely fixed incident and reflected wave vector. Information on atomic displacements is then obtained in three dimensions at the surface and in the bulk of the material by computationally inverting the three-dimensional diffraction pattern using iterative projection or machine learned neural network phase retrieval techniques [12][13][14].
The resolution and accuracy of the recovered three-dimensional image is highly dependent on the BCDI experimental geometry (oversampling ratio), the presence of noise in the measurement and to some extent, the sample crystalline quality. A means therefore to calculate the diffraction pattern of a single nanocrystal with arbitrary atomic displacements can provide useful insight into the experimentally observed diffraction pattern and can also facilitate in recovery of the associated morphology and phase information.
To date, practical attempts to simulate diffraction from a single crystal with arbitrary atomic displacements have either exclusively relied on approximating scattering from the atomic structure to the far-field imaging plane with a Fourier transform or employed random sampling methods [15][16][17][18]. In such a framework, there is generally no direct means to specify atomic displacements and the associated phase information is instead used to approximate atomic displacements. As a result, potential difficulties arise with such approaches when attempting to correlate specific atomic displacements with the phase information.
The most general approach within the kinematic approximation consists of a summation of each scattering event from each atom to arrive at the scattering intensity at a single location on the imaging plane. This can result in a total number of summation of the order of 10 9 for a pixelated imaging plane. Although this approach is not always prohibitive on modern accelerated Graphical Processing Unit (GPU) computing facilities or High Performance Computing (HPC) facilities, there are certain redundancies (outlined below) in the standard general approach that when removed, further and significantly increases the efficiency of the calculation and enable its use in more demanding applications such as training data for machine learned phase retrieval optimisation [19,20].
In this report a highly accurate and efficient means to simulate diffraction of coherent x-rays from a single crystal with arbitrary lattice distortions is presented. By doing so, we show that existing implementations based on summation of scattering intensities can benefit in efficiency when utilising the approach presented herein.

Model description
Our method utilises the general scattering amplitude equation within the kinematic approximation to generate a diffraction speckle pattern by individually considering the contribution of each atom within an isolated crystal [21,22]. Refinements are then made to this standard scheme. The illumination source is assumed to be larger than the crystal and fully coherent. Partial coherence effects however are readily incorporated, as previously demonstrated for a subset of special cases [23,24].

Scattering from a nanocrystal supercell ensemble
The description begins by considering a monochromatic beam of light with intensity Φ 0 (photons/m 2 /s) and wavelength λ impinging on a single nanocrystal. This nanocrystal is located at the origin O which coincides with the origin in the laboratory frame of reference (figure 1). The primary coherent x-ray beam is oriented in the positive z-direction in the global co-ordinate system and has no component in the x or y-directions. The unit vector for the primary beam is denoted byk i . The position of each atom in the crystal is given by  detector in the general directionk P f, . We define the wavevector transfer at the centroid as Q = k f − k i and the general wavevector transfer to a point P as q = k f,P − k i .
The instantaneous scattering amplitude A at a point P, a distance R from the nanocrystal, is then in general given by: where r 0 is the Thompson scattering length and f m (q) is the atomic form factor. When there is no variation in atomic arrangements within the unit cell, equation (1) is often divided into contributions from the lattice and unit cell structure factor S(Q): If atom m is allowed to displace by an amount u m (R n ), equation (2) is no longer so easily separated and the scattering amplitude becomes: The scattering intensity is then given by: In general, to obtain a slice of the three-dimensional diffraction pattern (the third dimension and rocking curve procedure are discussed below), equations (3) and (4) are computed for each detector pixel by varying the reflected scattering vector k f,P . The number of photons scattered into a detector pixel of length d x,y is then given by (3) is often prohibitively challenging due to the summation over each atom in each unit cell in the crystal which can extend beyond 10 9 terms. To overcome this challenge, the summation over each unit cell is instead replaced with a summation over groups of unit cells (or supercells) of a fixed size that are adjacent within the crystal. To accomplish this, lattice vector displacements n i are replaced by n i d i , where d i are positive constant integer values for each index i. We then define a supercell lattice vector R R r nd m nd m = + and R nd = n 1 d 1 a 1 + n 2 d 2 a 2 + n 3 d 3 a 3 . Substituting into the scattering amplitude equation and simplifying yields:

Computation of equation
whereũ m is the average displacement of atom m in the supercell at position R nd . We will demonstrate below that grouping unit cells into supercells within the scattering amplitude summation is justified provided the total size of each supercell remains below the resolution of the measurement which is determined by the experimental geometry.

Geometric considerations
In order to determine the appropriate supercell size, geometric considerations largely related to the sample-todetector distance and the rocking curve angular width are considered.

Detector scattering
Consider a detector of size X×Y pixels with pixel width of d x,y . If the centroid is pointing at the centre of the detector at a distance of R, the angles α i , α j subtending pixel indices i, j (in the x and y-directions respectively) are given by: Reflected vector k f,P is then rotated onto each pixel such that:

Sample to detector distance
The angle α subtended by a single pixel is given by α = d x,y /R. The fringe spacing w f (1/m) for a given real-space direction is approximated as w f ≈ 2π/(N 1/3 |a|), where a is a lattice vector and N is the total number of unit cells. To sample the data at an oversampling ratio of β, each fringe must at least cover 2β pixels such that w f /(2β) = α|k|. The optimal sample to detector distance at the centroid is then given by: x y x y f x y , The third dimension of the diffraction pattern is obtained by incrementally rotating the nanocrystal by an angle Δθ while recording the scattering intensity at the detector. This procedure is equivalent to rotating the Q-vector relative to the nanocrystal. The rocking curve increment Δθ is the angle subtending w f /(2β), normal to the scattering plane (aboutˆk k i f ) and is given by: The rocking curve width max q is the angular range within which the object crystal must be rotated in order to acquire the third dimension of diffraction. A good estimate of this can be obtained by considering the distribution of diffraction intensity as a sinc function with a maximum intensity of ( ) = F and a fringe width of w f . For an exposure time of τ seconds, the average decay in intensity at a point P with wavevector k f,P can be described by the envelope function: and is included for practical purposes. We can then evaluate q q max = when ( ) E q max takes its minimum value of ( ) E q n d x y max 0 , 2 = , where n 0 = 1 is the minimum photon count at that pixel. Using equation (10), it is then possible to obtain q max and hence obtain the total rocking width max q as: The resolution ΔX at which the diffraction data is obtained then becomes in units of metres.

Diffraction pattern coordinate transformation
Diffraction amplitude measurements are in general recorded in the coordinate system of the experimental setup which is typically not a rectilinear coordinate system. As a result, the amplitude and phase information reconstructed from diffraction amplitude measurements will also not reside in a rectilinear coordinate system. It is however possible to transform the experimental coordinate system into a rectilinear coordinate system. We demonstrate this using the following procedure for a θ rocking curve measurement. Analogous procedures apply for a f rocking curve measurement. Ifẑ,ŷ andẑ are unit vectors in a rectilinear coordinate system, unit vectors for the incident k i and reflected k f wavevectors are generally given by: where θ and f are angular rotations about theŷ andẑ axes respectively. This system of rotations is chosen so that thek i unit vector is aligned along thek vector direction in the absence of rotation. Angles θ 0 and f 0 are additional angular rotations about theŷ andẑ axes respectively needed for the k f vector to adopt the Bragg condition for diffraction. Angle It is then possible to define reciprocal lattice vectors * a j in a rectilinear coordinate system as follows: where Δθ are angular increments in the rocking-curve measurement and Δθ 0 , Δf 0 are angular displacements of thek f vector. The coordinates generated will in general differ for each speckle pattern that is obtained from a unique Bragg reflection. In order to generate coordinate points common to a number of speckle patterns, it is in general necessary to interpolate onto a regular grid.

Effect of supercell size on speckle pattern
To demonstrate the efficacy of our framework, coherent x-ray scattering from a range of nanocrystal material systems was calculated first for a range of supercell sizes, in the absence of lattice distortions. The detector is assumed capable of single photon counting so that one photon detected corresponds to one count on the detector. Each diffraction speckle pattern was computed using a Nvidia V100 Graphics Processing Unit (GPU) as provided by the Iridis High Performance Computing Facility at the University of Southampton. The total compute time varied from approximately 2 hours (without supercell grouping) down to a few seconds, inversely proportional to the supercell grouping size. Poisson noise was then added to each diffraction pattern and the χ 2 error, which has a possible range of [0, ∞ ], was then calculated using equation (15). In equation (15), I d is the scattering intensity for a supercell of size d (unit cells) and I 0 is the scattering intensity in the absence of supercell grouping (i.e. 1 unit cell per supercell, which is the ground truth). The resulting χ 2 values are plotted in figure 2 using the relevant parameters in table 1. Results are shown for a single Bragg reflection of each material only as similar results were obtained for other reflections (See supporting material is available stacks.iop.org/JPCO/6/055003/mmedia). In each case, a relatively small χ 2 error (< 0.01) is seen to persist for supercell sizes below the resolution (ΔX 3 in units of nm 3 ) and rapidly increase at an exponential rate beyond the resolution limit. The normalised compute time is also plotted against each supercell size in figure 2. This is shown to rapidly reduce from the initial values (9, 11 and 7 minutes, respectively) down to a few seconds with an increase in supercell size.
Once a suitable supercell size was established below the resolution threshold, coherent x-ray scattering was then calculated for the same range of nanocrystal material systems with the addition of predetermined atomic lattice distortions. Diffraction patterns from nanocrystals of zinc oxide (ZnO), barium titanate (BaTiO 3 or BTO) and lithium cobalt oxide (LiCoO 2 or LCO) were obtained using the parameters shown in table 1. In each case, the size of the supercell remained below the resolution of the measurement (ΔX 3 ).
The nanocrystal electron density and phase information was subsequently reconstructed from the resulting scattering patterns with the inclusion Poisson noise. Standard iterative phase retrieval techniques were used and are available in the Interactive Phase Retrieval Suite software package [25]. This typically began with 1000 iterations of the HIO algorithm using a fixed solvent support region and a mask to eliminate amplitude information below the single photon (per pixel) threshold. This was immediately followed by 5,000 iterations of the Shrink Wrap algorithm with the HIO Mask constraint, as before. Finally 100 iterations of the Error Reduction algorithm were performed. Figure 4(a) shows the diffraction pattern of the (101) reflection of a ZnO nanocrystal with nominal dimensions of 95 × 95 × 200 nm. ZnO has Wurtzite crystal structure (space group P6 3 mc) with lattice parameters a = b = 0.325 nm and c = 0.52 nm. Zn ions are known to migrate to either in-plane or out-of-plane interstitial sites, as shown in figure 3(a). In this crystal phase information in the reconstructed image results from sinusoidal displacement of the Zn ion, at the origin of each Wurtzite unit cell, normal to the (002) direction and along the (110) direction. The period of the displacement was commensurate with the length of the nanocrystal with a peak-to-peak amplitude of 20%, in Wynckoff fractional units. The total compute time was less than 10 minutes. The diffraction pattern ( figure 4(a), left) clearly shows six fold rotational symmetry as expected from a hexagonal crystal with some distortions due to the aforementioned atomic deviation from the ideal lattice. Iterative phase retrieval as performed on the diffraction pattern successfully recovered the morphology of the crystal with the correct dimensions and also the phase information which is visible in figure 4(a), upper right. A sinusoidal variation in the phase between ± 0.3 radians is observed along the length of the nanocrystal while uniform phase is observed in transverse planes, in good agreement with the calculated ground truth ( figure 4(a), lower right).

Ferroelectric phase transition in BTO
Barium Titanate (BaTiO 3 or BTO) is a ferroelectric material that undergoes a phase transition at T c = 393 K [26][27][28]. Below this temperature, ferroelectricity arises from the non-centrosymmetric properties of the lattice as the Ti cations and O anions are displaced from their equilibrium positions in opposite directions along the caxis giving rise to a spontaneous polarization ( figure 3 (b)). A 180°ferroelectric domain structure is considered for a cubic nanocrystal of nominal dimensions 400 × 400 × 400 nm. The Ti ions are displaced by 5% in Wyckoff fractional units while the oxygen octahedron ions are displaced in the opposite direction; the O ion at Wyckoff position (0.5 0.05 0) is displaced by 11% and the remaining O ions are displaced by 7%. The domain structure is incorporated by using a square wave function for the displacement of the ions along the (110) direction with a period of 400 nm. Figure 4(b) shows the simulated diffraction pattern of the (101) reflection with the four-fold symmetry expected of a cubic crystal.
Reconstruction of the diffraction pattern using iterative phase retrieval has recovered the correct morphology of the nanocrystal as well as the expected domain structure (4(b), upper right). A sudden shift in the phase from positive to negative clearly exhibits the 180°shift in the polarization. A region in the centre of the crystal of zero phase is clearly visible indicative of the formation of a planar domain wall structure.

Lithium ion migration in LCO
Lithium Cobalt Oxide (LiCoO 2 or LCO) has a layered structure consisting of O atoms bonding with Co more strongly than with the Li ion. When an electric field is applied to the LCO crystals within a battery, charging occurs as the Li ions begin to migrate from their original positions towards the anode ( [29,30]). However, the crystallographic structures can present dislocations in the material. These can prove beneficial for cycling, as the dislocations could facilitate an energetically improved diffusion path for the Li ions to be transported between their lattice site within the crystal structure and the electrolyte or electrode ( [31,32]). Yet, dislocations are by definition defects within the lattice, that can eventually initiate mechanical degradations ( [33]). To simulate dislocations, an LCO nanocrystal of 350 × 350 × 200nm nominal dimensions was generated with a slip taking place along the1120 á ñ lattice directions in the {0001} planes. The resulting diffraction pattern is shown in figure 4(c) (left) for the (104) Bragg reflection. Iterative phase retrieval algorithms have successfully reconstructed the nanocrystal's morphology, including the phase information ( figure 4(c) top right). The morphology of the reconstruction resembles the expected morphology (figure 4(c) bottom right), but the observed additional surface roughness of the amplitude is likely due to the high aspect ratio unit cell, resulting in a relatively wider variation in the fringe spacing, and hence in the resolution, both along and normal to the high symmetry axis of the diffraction pattern. Regarding the phase information within, it is predominantly neutral, except for the created dislocation along the1120 á ñ direction.

Discussions
This work demonstrates an atomistic approach to simulating Bragg coherent diffraction from a single nanocrystal. We have shown that provided that the unit cells are grouped into supercells with dimensions below the resolution of the measurement, it is possible to reconstruct the amplitude and phase information from the scattering intensity alone. Our measure of the resolution is an approximation based on scattering from a small and perfect nanocrystal. Bragg CDI experiments performed in practice will likely result in lower resolution measurements due to inherent defects present in most crystals that can reduce visibility in the speckle pattern.
In addition it should be noted that although our framework correctly assumes that the scattering intensity recorded by each pixel of the area detector results from scattering into a solid angle subtending that pixel, the phase reconstruction process instead considers each pixel as having recorded intensity from an infinitesimally small region of reciprocal space. Coupled with the finite boundary of the acquired scattering data, this can result in aliasing in the reconstruction process. A means to mitigate this effect is to expands the spatial frequency range by padding the data.

Conclusions
In summary, we have demonstrated the feasibility of atomistic computation of Bragg coherent x-ray diffraction scattering from nanoscale crystals and three-dimensional image recovery that is limited in resolution primarily by the diffraction experimental geometry. By grouping unit cells into supercells of dimensions below the BCDI experiment resolution, it was demonstrated that our approach can accurately simulate the coherent x-ray diffraction pattern of a single nanocrystal on a significantly reduced time scale (typically two orders of magnitude) and subsequently recover the correct electron density and phase information. Our approach is therefore of considerable utility for use in BCDI diffraction pattern simulation, data interpretation and materials design.