Imaging of complex density in silver nanocubes by coherent x-ray diffraction

When using coherent x-rays to perform lensless imaging, it is the complex wave field exiting the sample or, in the case of the Bragg geometry, the deformed electron density distribution of a crystal, that is being sought. For most samples, to some extent, the image will be complex, containing both an amplitude and phase variation across the sample. We have developed versions of the hybrid input–output (HIO) and error reduction (ER) algorithms that are very robust for the inversion to complex objects from three-dimensional (3D) coherent x-ray diffraction (CXD) data measured around a Bragg spot of a small crystal. The development and behavior of these algorithms will be discussed in the context of inverting a 3D CXD pattern measured around a (111) Bragg spot of a silver nanocube.


Introduction
It has recently been realized that coherent x-ray diffraction (CXD) in the Bragg geometry presents a unique capability with sensitivity to lattice distortion caused by strain within a small crystal [1]- [3]. The technique can image, in three dimensions, not only the shape of the object, but also a projection of the distortion of the crystal lattice on to the Q vector of the Bragg spot measured [4]. The lattice distortion will be manifest in the measured diffraction pattern as a local asymmetry around a given Bragg peak. Upon inversion to direct space, the resulting image of the object will be complex, where the amplitude represents the density of the object and the phase is a simple projection of the displacement of the atomic positions of the crystal from an ideal lattice. The inversion step requires the recovery of the reciprocal space phases lost in the measurement of the diffraction intensities [5,6]. In this paper, we explore recent improvements in the computational algorithms that are appropriate for inverting diffraction patterns with this unique type of phase contrast.
In a more traditional incarnation, phase contrast imaging is widely used in optical microscopy, while more recently it has also become an x-ray technique [7,8]. With varying degrees of sophistication, they all share some method of converting phase contrast, usually due to a complex refractive index, into a measurable intensity variation. Diffraction-based microscopy seeks to directly exploit the asymmetry inherent in diffraction from complex objects to image the phase structure of the object. The measurements are carried out in a variety of ways. Most recently, curved beam illumination has been demonstrated to enhance the recovery of the reciprocal space phases lost in the measurement, while at the same time providing a lower resolution hologram for use in the diffraction inversion [9]- [11]. The known phase structure of the incident illumination is exploited in the reconstruction algorithm to give very rapid convergence. A curved beam measurement has not yet been exploited to enhance the recovery of phases in the Bragg CXD geometry. In our case, we illuminate the sample with a plane wave and measure the far-field coherent diffraction around a Bragg peak of the crystal. To enable the inversion to direct space, the diffraction is over sampled at greater than the Nyquist frequency associated with the size of the object under consideration. Algorithms using known information in direct space, such as the finite size of the object and reasonable limits on the phases, can then be used to invert the measured intensities to an image [6].
The computational methods are based on image-processing techniques pioneered by Fienup in the 1970s and 1980s that calculate the object as the Fourier transform of the (phased) diffraction pattern [12]. As originally found, the methods tend to stagnate with a solution that mixes the true object and its 'twin', obtained by reversing the directions of all the three coordinates. Both the object and the twin give the same diffraction pattern, but the superpositions do not; however, the calculations are frequently found to stagnate. Fienup proposed the hybrid input-output (HIO) method that worked well to break the stagnation. Variations of HIO are still widely used today, including in our current work. But the choice of 'direct-space constraint' that is applied to update the object in direct space on every iteration within HIO is still the subject of experimentation. We present results here for a new direct-space constraint that is particularly effective for inverting the diffraction patterns of sub-micrometer crystalline objects and modifications of the reciprocal space operations of the algorithms that improve our results.
The samples studied here were chemically synthesized 'nanocubes' of silver, obtained by precipitation from silver nitrate solution in the presence of a PVP surfactant [13]. By controlling its concentration, the PVP limits the size of the precipitate to nanometer dimensions and also influences its shape. Cube-shaped crystals with about 200 nm edge dimension were used in these experiments, similar to those shown in figure 1. A sample similar to this has recently been imaged at 3 nm resolution in two dimensions (2D) using coherent diffraction imaging in smallangle geometry [14], where the data were manually symmeterized to result in a real image. This is likely an accurate representation of the true object given that the sample is almost a perfect crystal. The structures of these and other nanocrystals is generally interesting because of their potential applications as an industrial material. Because of the small size, a strong influence of the surfaces, edges and vertices can be expected in their physical and chemical properties. At this size, the crystal lattice structure is expected to remain FCC, although for much smaller nanocrystals of gold, this is known to break down and be replaced by decahedral and icosahedral phases [15]. When the crystal size starts to approach this limit, deviations in the structure are expected to become more pronounced. These would be visible as strain (or complex phase) in our direct-space images. In effect, the surface, edge and vertex contributions would be seen to encroach upon and eventually dominate the FCC core.
Since the growth conditions are deliberately far from equilibrium, other entirely new phenomena may be present that could be detected from phase contrast microscopy. Both silver and copper films grown at 100 K are known to be highly defective with 2% vacancy concentration [16] and up to 1% strain [17]. When the growth conditions are sufficiently extreme, the crystal remains FCC but has so many imperfections that the lattice constant and density are visibly altered. Within the chemically synthesized Ag crystals, there is the distinct possibility of finding density variations (at the 2% level) and/or strains up to 1%. The size-limiting effect of the capping effect of the PVP surfactant will mean that the rates of growth of the center and the outside regions of the crystals will be very different. Note that an inhomogeneous 1% strain of a 200 nm Ag crystal measured at its 111 Bragg peak will cause a phase shift of 10π or five complete wraps of the phase. The sensitivity of the CXD method is below 0.1 rad. In the present case, we observe a total phase range within the cube image of about 0.7 rad as shown in figure 2 We previously measured CXD diffraction patterns from Ag nanocubes similar to those used in the present study [18]. The advance reported at that time was the use of Kirkpatrick-Baez (KB) focusing optics to concentrate the incident x-ray beam once it had been made coherent with a small slit. In the experiments reported here, 200 nm cubes are illuminated by KB optics with a focus size around 1.5 µm; in this case the object can be assumed to sit within a plane wave illumination.

Experimental details
A suspension of single-crystal silver nanocubes measuring about 200 nm was deposited on a glass slide and allowed to dry at room temperature. The sample was then mounted on the diffractometer at Sector 34ID-C of the Advanced Photon Source. A front-illuminated direct-detection CCD camera was mounted on the goniometer arm and oriented at the (111) Bragg angle of the silver reciprocal lattice. The 9 keV coherent x-ray beam was slitted to 50 µm × 2 µm (V × H ) and focused on to the sample with a pair of KB mirrors to give approximately 1 µm beam size. The sample was then scanned across the beam with an incidence angle of 1 • until diffraction from a randomly oriented individual nanocrystal was seen on the CCD camera. While thousands of crystals are likely in the beam at any given position of the sample substrate, only one or two will be oriented at the (111) Bragg angle of the silver lattice. Frequently, samples of this nature are not stable enough in the x-ray beam to acquire data for a very long time. This could be due to charging effects of the particles, which then repel each other and move out of the beam. Upon alignment of this crystal in the beam, a 3D CXD pattern was recorded by rotating the crystal in 0.01 • steps through a total of 0.3 • , similar to a traditional rocking curve measurement in crystallography. Part way through the repeat scan, the particle was seen to leave the beam. Due to the limited dynamic range of the CCD, the individual acquisitions at a given angle were accumulated in five exposures that were 2 s each. A selection of the recorded frames is shown in figure 3 along with a 3D representation of the data at a single intensity contour level of 400 analogue digital units (ADU) just above the background. The asymmetry visible in the 3D diffraction pattern is a manifestation of distortion caused by strain within the lattice of the crystal [4]. The readout error from the CCD is visible in figure 3(e). The known shape of this object has made it an ideal playground for the development and testing 6 of phasing algorithms for the recovery of complex objects in the Bragg CXD geometry. The effect of a given algorithm or parameter can be immediately deduced by inspection of the object. That which results in something with flatter faces or sharper edges is deemed a better result.

Data preparation
The raw data frames were cropped to isolate the diffraction from the cube. The individual frames were further binned since the raw data were oversampled by more than a factor of 20 in the plane of the detector. The oversampling of the diffraction pattern is measured by essentially counting the pixels covered by a given fringe in the pattern. The Nyquist frequency, or critical sampling rate of the diffraction pattern, is one half the reciprocal of the size of the object. The Nyquist frequency will sample the complex amplitudes at both the peaks and the troughs of the fringes. This sampling rate in reciprocal space, however, is undersampling the intensities. Upon loss of the phase information in the diffracted wave, critically sampling the diffraction complex amplitudes will not show the zero of intensity between the fringes. To sample the intensities critically, one must measure the diffraction pattern at least twice as finely as the complex amplitudes, or oversample the diffraction amplitudes by a factor of two. The intensities can also be thought of as the Fourier transform of the autocorrelation of the object. The autocorrelation is twice the size of the object in every dimension, so the critical sampling rate in reciprocal space will be twice that of the object [6,19]. In our case, this oversampling ratio is much higher and is controlled by the distance between the detector and the sample. Approximately, 1 m distance from the sample to the detector resulted in fringes covering about 20 of the camera's 20 µm pixels in each direction parallel to the detector. Our sample rate perpendicular to the detector is given by the incremental rotation of the sample. In this case, the oversampling rate is about 14 and was not further binned.
The data were then shifted around in the 3D array until the central bright spot of the autocorrelation had an approximately uniform phase of zero, indicating that the amplitudes are reasonably centered in the reciprocal space array. The centered array was then symmetrically cropped and zero padded so that the extent of the measured diffraction pattern was about 50% the total size of each array dimension. The reason for zero padding the final data array is to prevent aliasing from the fast Fourier transform interfering with the recovery of the phases of the measured data [19]. Aliasing cannot be avoided when discrete Fourier transforms are implemented. A compact object in direct space inherently requires infinite Fourier components to represent it in reciprocal space. When a discrete Fourier transform is used, the higher order harmonics will wrap across the boundary of the computational array and be added to the complex amplitudes on the other side. By zero padding the array, the aliased complex amplitudes can be prevented from overlapping and corrupting the phases recovered for the measured moduli.
A single-photon event in a pixel of the CCD results in approximately 300 ADU. A measured dark image (shutter closed) averaging about 1400 ADU in each pixel was subtracted from each frame prior to use. There is significant noise inherent to the CCD and can be seen in the isosurface through the data in figure 3(e) as streaks to the right of the greatest intensities recorded. This noise is dealt with in the algorithms described below.
The ith element of the reciprocal space array after the modulus constraint is applied The square root of the ith measured intensity The amplitudes and phases of the Fourier transform of u

Algorithms
Many algorithms are employed for inversion of coherent diffraction. Fienup's HIO and error reduction (ER) have been implemented with a variety of constraints that depend on the problem being solved [12]. The basic principle behind all of these algorithms is to iterate between direct and reciprocal space while imposing constraints of some type on the data in each realm. The goal is to modify the complex amplitudes, in both spaces, in some way that drives the algorithm towards a self-consistent solution to the phase problem. Table 1 gives the definition of the reciprocal space error (RSE) metric we use to monitor the progress of the algorithm toward convergence, as well as the definitions of symbols used in descriptions of the algorithms below.
In the case of HIO the term constraint is perhaps a bit severe. By its nature HIO tends toward a solution in a rather indirect way [20]. The rules are formulated in a manner that inhibits direct space amplitudes from attaining an appreciable magnitude in regions of the computational array inconsistent with the known properties of the sample. In the case of a support constraint, the amplitudes are discouraged simply based on the position of the given volume element (voxel). Those outside the compact support region, identified with the object position, are inhibited from growing larger through the feedback component of the algorithm, the lower part in equation (1). Another rule that can be implemented in direct space is on the phase of the given voxel. If upon entering the direct space portion of the algorithm a particular voxel within the current direct space estimate of the object is found to have a phase outside of some specified range, the HIO feedback mechanism will be used to inhibit or perhaps reduce in magnitude such an amplitude. This algorithm, referred to as phase-constrained HIO, is summarized in equation (1). In each expression, i indicates that the operations are done on a per voxel basis, i being the index of a given voxel in the computational array. The only free parameter in the algorithm is β and is typically set between 0.7 and 0.9 [12]. It essentially controls the convergence rate of the iterations through the feedback component (lower condition in equation (1)) of the algorithm.
Frequently, HIO is combined with iterations of the ER algorithm to solve stagnation problems with HIO. It has been found that HIO will get stuck in orbit around one or more local minima in the solution space. Iterations of ER can be used to collapse the current iterate to one of them [20]; ER, with the particular constraints of equation (2), can be shown to be identical to a steepest descent or gradient search algorithm [12]. In the case of the data obtained from this cube, the ER iterations tend to destroy the object recovered by HIO. We believe that this is due to the rather high noise level of the data, which is over-fit by a steepest descent algorithm.
To solve both the stagnation of HIO and the over fitting of ER, we have implemented a new form of ER. Termed 'phase-only ER', it has been designed to directly deal with a particular source of stagnation of HIO in the recovery of complex objects and to avoid the pitfalls of ER's steepest descent nature. The algorithm is typically run after a number of iterations of phaseconstrained HIO have already recovered an image of the object. It operates solely on the phases of the current iterate. The idea is to leave the object within the support region untouched and just rotate the complex amplitudes outside of the support region to the nearest branch of the real axis at every iteration through direct space. This is represented by equation (3), where the '±' indicates that the complex amplitudes are rotated to the nearest branch of the real axis. It will be shown later that the stagnation of HIO in our case is directly associated with phase features outside of the support region. By forcing the complex amplitudes to be real outside of the support, the object recovered by HIO is left largely intact, which is not the case with ER in a steepest descent mode, and also gives a relatively low error metric.
Up to this point, only the direct space operations of the algorithms have been discussed. The reciprocal space constraint we use during HIO is significantly more straightforward. The current iterate is Fourier transformed and the amplitudes of the Fourier transform are replaced with those of the measurement, leaving the calculated phases intact. The only addendum to this recipe in reciprocal space is the use of an amplitude threshold. We only assign the experimental amplitudes to the calculated phases where those amplitudes exceed a threshold value; the remaining amplitudes and phases are left at their calculated values. Equation (4) expresses this operation. The threshold of 400 in intensity or 20 in amplitude was used for all subsequent work with this data set.
Yet another variant of ER has proven, in our case, to give good results in the face of noisy data, particularly in the case of a noisy background. This variant of ER is illustrated in equation (5) by the reciprocal space operation. The algorithm is simply ER with the usual support constraint in direct space. By setting the amplitudes in reciprocal space to zero, where the experimental amplitudes are less than a given threshold, this version of ER avoids the formation of the strong phase features in direct space that leads to a poor RSE. Doing the same in the reciprocal space operation of HIO has not worked for these data, however; an image of the cube is not recovered.

Phasing
The support was initially determined from the size of the autocorrelation of the data. The autocorrelation should be twice the size of the object in every direction, so taking an initial support that is greater than one half of this size will be sufficient to contain the object. After a couple of cycles of phasing and adjusting the support, one was found that just contains the object. A phase range of ±π/2 was used as the phase constraint in HIO for the initial inversions. Figure 4 shows the object obtained with 20 iterations of phase-constrained HIO. The translucent box in figure 4(a) shows the cube in the support obtained after a few adjustments of the support.
In general, if the support box is too small, as in figure 4(b), the recovered object will look distorted and will have flat sides corresponding to the sides of the support. If the support is qualitatively too big, as in figure 4(c), something relatively compact will be visible in the support box, but it is usually pinned to one wall or the other. It is relatively easy to estimate the amount by which the support should be reduced in this case. We are also developing algorithms based on the familiar Shrink Wrap methods for use in 3D phasing of Bragg geometry CXD [21]. The RSE as a function of iteration during phase-constrained HIO typically falls rapidly to some minimum and then climbs to usually between 0.1 and 0.2 for reasonable-looking objects from these data. This effect is illustrated in the plots of RSE in figure 5. Figure 6 shows the object recovered in 100 iterations of phase-constrained HIO as a function of the magnitude of Figure 5. The RSE as a function of iteration for each of the listed phase constraints in HIO. In each case, phase-only ER was run for 20 iterations after 100 iterations of HIO. The support is the one that was visibly determined to be the best fit to the recovered object with a ±π/2 phase constraint shown in figure 4. All runs were identically initialized with the support filled with a constant real amplitude. the phase constraint. In each case, the HIO phase constraint was between positive and negative of the given value. The algorithm is initialized with an object that is the support box filled with a constant uniform real density in all cases. Since our goal is to test individual components of the phasing procedure, this seems to be a justifiable starting point for all runs. Since no iterations of any variety of ER have been run, all of the results shown in figure 6 have a relatively high error metric. The RSE of each is between 0.12 and 0.15. We have found that the proper phase range for a given data set can be determined in a similar fashion to the support. After establishing a support size with a large phase range in HIO, the phase of the object was inspected and the phase range adjusted to better fit the observation. In the case of this cube, the initial inversions done with a phase range of ±π/2 indicated that a phase range of 0.87 would be adequate for the object. This value was used in HIO for subsequent work on the data set. If the phase range were made too small, we see obvious holes or a broken density in the recovered object, as in figure 6(a), with the phase at the boundaries of the density running right up to the phase limits that were imposed. This effect can be seen in the object recovered with the ±π/8 phase range in figure 6. The entire object is badly deformed. When the algorithm was given too large a phase range, the recovered object has a distorted shape as seen in figures 6(d)-(h). Phase-constrained HIO probably improves the phasing procedure for complex objects by inhibiting the formation of the centrosymmetric twin. In this case, the twin image would not only be inverted around the origin but also be the complex conjugate of the other. If the actual distortion of the crystalline lattice is less than a total of π in phase, and the phase range admitted by HIO is adjusted accordingly, the phase constraint will prevent the twin from forming in the reconstruction process.
Clearly HIO with a phase constraint gives a good-looking solution, as seen in figure 6(c). The recovered object looks like the expected cube, with relatively sharp corners and flat sides. It is common for a very good-looking reconstruction, regardless of algorithm parameters, coming out of HIO to have an RSE around 0.1, whose Fourier transform is even visibly a poor fit to the data. Figure 7 shows a line through the center of the 3D diffraction data and the identical line through the Fourier transform of the object after HIO. The spikes in the line through the amplitudes of the Fourier transform of the object must be due to features in direct space that are spaced quite far apart, since small features in reciprocal space are due to long-range correlations in direct space. Since the feedback mechanism in HIO limits the magnitude of the density outside of the support, the strong features in direct space giving rise to the oscillations in reciprocal space must be coming from something other than the amplitude of the density. We have identified phase vortices in the region of low amplitude outside the support. Phase vortices are well known to be a source of stagnation in phasing methods [2]. Usually these vortices occur around pixels of low amplitude in reciprocal space and give rise to features similar to what we are seeing, but in the direct space object. In our case, since the density is complex and has low amplitude outside of the support region, we are seeing vortices forming in the direct space image leading to strong features in reciprocal space. Figure 8 shows the cube at the center of the array with a slice plane through the phase of the complex density. The strong phase features leading to the oscillations in reciprocal space are evident outside the support region. Figure 9 shows the RSE as a function of iteration during the HIO execution. Plotted with the RSE is the number of direct space voxels around which a phase vortex was found. The two graphs clearly track each other well. In principle, the ER algorithm should overcome this sort of stagnation by just setting the problem amplitudes in direct space, outside of the support, to zero in both amplitude and phase. The result after 20 iterations of HIO, figure 10(a), followed by ten iterations of ER with just a support constraint is shown in figure 10(b). As mentioned earlier, ER with support as the only constraint in direct space is equivalent to a steepest descent search. It seems the noise of the data is over-fit, which leads to the object recovered in HIO  being destroyed. Essentially we think that HIO is doing the bulk of the phasing work, but the direct space vortices are artificially corrupting the RSE measure of the fit to the data.
To deal with the vortices in direct space, in the region outside of the support box, we have developed a couple of variants of the ER algorithm that are effective. In the first, chronologically, the phases of the amplitudes outside of the support region are rotated to the nearest branch of the real axis. Coined phase-only ER, this operation in direct space will directly remove the Figure 11. (a) Solid isosurface (top) at the 50% amplitude in the result obtained from 20 iterations of HIO with the phase constraint of ±0.87 rad. Below, the isosurface is rendered semitransparent with a slice plane through the phase of the object. The RSE was 0.14. (b) The result of ten iterations of phase-only ER following the HIO done to produce (a) with an identical slice through the phase of the object. The RSE of (b) is 0.0009. The phase scale bar is in radians.
phase vortices outside of the support. As can be seen in figure 11(b), the object recovered by HIO, shown in figure 11(a), is left largely intact in both amplitude and phase by the phase-only ER iterations. At the completion of 20 iterations of HIO, with a phase constraint of ±0.87 rad, the RSE was 0.14. At the conclusion of ten iterations of phase-only ER, the RSE had fallen to 0.0009. The HIO iterations were initialized with a constant real amplitude within the support region. A recently published work of Newton et al applied this algorithm to successfully reconstruct images from multiple Bragg peaks of a zinc oxide microcrystal. The multiple reconstructions were then combined to obtain the full 3D vector displacement field of the crystal and hence the six independent components of the strain tensor of the sample [22].
Another method to reduce vortices in direct space, and improve the RSE of a given solution, is through a modification of the reciprocal space operation of ER. We have found that ER with nothing more than a support constraint in direct space combined with setting amplitudes in reciprocal space to zero, at points below the given threshold, works almost as well as phaseonly ER for this data set. Figure 12 shows that the shape of the cube was left intact, and the phase structure within the cube was also maintained close to the values obtained by HIO. The RSE fell to 0.005, not quite as low but similar to that obtained from phase-only ER iterations following HIO. If the low amplitudes were zeroed in the HIO iterations, the results are less reliable. Since the same HIO result is used to test both variants of ER, the objects recovered by  figure 11; the 3D image is identical. (b) Ten iterations of ER, following the HIO of (a), with a support constraint in direct space and amplitudes less then the threshold in reciprocal space set to zero. The RSE of (b) is 0.005. The phase scale bar is in radians.
both are very similar, although it is clear from the differences visible in figure 12 that zeroing the amplitudes in reciprocal space has a small impact on the phase structure of the solution.

Summary
In summary, we have found that a phase-constrained HIO algorithm combined with a variant of the ER algorithm is ideal for the phasing of a diffraction pattern from a compact complex object such as the silver nanocube presented here.