Phase retrieval by coherent modulation imaging

Phase retrieval is a long-standing problem in imaging when only the intensity of the wavefield can be recorded. Coherent diffraction imaging is a lensless technique that uses iterative algorithms to recover amplitude and phase contrast images from diffraction intensity data. For general samples, phase retrieval from a single-diffraction pattern has been an algorithmic and experimental challenge. Here we report a method of phase retrieval that uses a known modulation of the sample exit wave. This coherent modulation imaging method removes inherent ambiguities of coherent diffraction imaging and uses a reliable, rapidly converging iterative algorithm involving three planes. It works for extended samples, does not require tight support for convergence and relaxes dynamic range requirements on the detector. Coherent modulation imaging provides a robust method for imaging in materials and biological science, while its single-shot capability will benefit the investigation of dynamical processes with pulsed sources, such as X-ray free-electron lasers.

P hase retrieval arises in fields such as electron microscopy 1 , crystallography 2 and astronomy 3 , and is a core aspect of coherent diffraction imaging (CDI) 1,[4][5][6] . Iterative projection algorithms 7,8 commonly used in CDI search for a solution by iteratively using projections and constraints in two separated planes. Since the first demonstration of CDI with X-rays 4 , the high brightness available at third-generation synchrotrons and X-ray free-electron lasers (XFELs) 9,10 along with the development of more sophisticated algorithms and experimental set-ups [11][12][13][14] , has led to impressive results in materials and biological science, such as imaging of nano-particles 15 , cells 16 and viruses 17 .
The current challenges for CDI include the high sensitivity of reconstruction convergence to noise or missing experimental data [17][18][19] , the limited ability to retrieve images for extended objects or complex-valued samples from single measurements 20 , and the need to collect data with a high dynamic range 17,19,21 . Various modifications to the illumination [22][23][24] have been proposed to address those issues for which the object appears as a perturbation on the illumination. A substantial a priori knowledge of the exit wave reduces the possibility that the algorithm becomes stagnated. Modifications to the sample exit wave have also been suggested [25][26][27] in which multiple measurements have been a key requirement for the convergence of reconstruction. We note that single-shot imaging schemes that use aperture arrays 28,29 or structured light to illuminate a sample 24 have recently been proposed, but so far only simulations and visible-light demonstrations have been reported.
The phase problem in CDI is usually formulated as recovering an image from its Fourier-magnitude data. In this paper we show that breaking the Fourier transform relation by inserting a wavefront modulator between the sample and the detector can greatly help to solve the phase problem with a single measurement. This coherent modulation imaging (CMI) method was initially demonstrated with visible light 30 . Here we extend its application to the X-ray regime using a substantially improved algorithm and show that CMI addresses many of the challenges that affect conventional CDI.

Results
Description of the experiment. Our CMI experiments were conducted at the cSAXS (coherent small-angle X-ray scattering) beamline of the Swiss Light Source. The experimental set-up is illustrated in Fig. 1a. Coherent X-rays of energy 6.2 keV passed through a 2.2 mm diameter pinhole, forming an illumination probe with a photon flux of 3.6 Â 10 6 photons per second on an object located 1 mm downstream. The X-rays transmitted by the object then passed through a wavefront modulator 2.5 mm further downstream before propagating 7.2 m to a Pilatus 2M detector with square pixels with a side length of 172 mm (ref. 31). The central 512 Â 512 pixels of the recorded pattern were used for reconstruction. The effective image pixel size at the planes of the modulator and the object was 16.4 nm.
Placing a wavefront modulator downstream of the sample is essential in this method. In our experiments, the modulator consisted of randomly distributed tungsten pillars with a nominal width and spacing of 200 nm. A tilted view of the modulator structure obtained using a scanning electron microscope is shown in Fig. 1b. The transmission of the modulator was characterized using ptychography [32][33][34] before measurements of the sample using CMI (Methods). The modulator transmission functions are shown in Fig. 1c,d, and these were used as a priori information in the CMI phase retrieval algorithm. In our experiments, about 70% of the X-rays were transmitted through the phase plate modulator and used for reconstruction. The fluence on the modulator was calculated to be 0.24 J cm À 2 , which is below the damage threshold of 0.5 J cm À 2 found for Fresnel zone plates made of tungsten and is well below the threshold of 59 J cm À 2 for Fresnel zone plates made of diamond when exposed to freeelectron laser pulses 35 . No beam-induced damage to the modulator was observed during the course of the experiments.
A flowchart of the CMI phase retrieval algorithm is shown in Fig. 1e. The support plane is one where the wavefield is known to have non-zero values within a finite region. It may correspond to the sample plane, but this is not a general requirement. Accurate knowledge of the axial separation between the sample and the modulator is not critical. Numerical wave propagation after reconstruction can be performed to bring the retrieved sample exit wave to focus when necessary. Reconstruction in CMI is performed by iteratively propagating a wavefield estimate between the support and detector planes, via the modulator. The applied constraints and the necessary inputs in each plane are indicated in Fig. 1e. The steps of modulation and wave propagation between the modulator and the support plane yield a rapid algorithmic convergence due to the removal of the twin image and spatial shift ambiguities whose competition with the correct solution is a potential cause of stagnation in conventional algorithms that use two planes 36,37 . Ambiguity removal in CMI is discussed in more detail in the Supplementary Note 1 and Supplementary Fig. 1. Details of the newly developed three-plane phase retrieval algorithm are given in the Methods section.
Single-shot reconstructions with a relaxed support. An example of CMI data recorded with a 3 s exposure time is shown on a logarithmic scale in Fig. 2a. The dashed square indicates the extent of the diffraction data used for the simulation of missing data shown in Fig. 3a-c. The requirement for high dynamic range detection and the consequent need for a beam stop are reduced in CMI ( Supplementary Fig. 2). The illumination wavefield shown in Fig. 2b was retrieved from a single-shot CMI measurement after 150 computational iterations with the loose support of a 3.2 mm-diameter disk. The amplitude is displayed as image brightness and the phase as image hue. The amplitude profile on the right side of Fig. 2b is taken along the central vertical line. It shows that the illumination oscillates on either side of the central maximum and has slowly tapered edges. The overall probe width is around 2.4 mm, measured to where the edge taper falls to half the value of the central maximum. The support requirement on the object exit wave can be significantly relaxed in CMI and this allows robust reconstruction even in the case of a tapered wavefield that would pose difficulties for conventional phase retrieval algorithms 38 . In our experiment, the modulator centre was approximately aligned with the illumination beam. A loose support centrally located around the sample exit wavefield was found sufficient to recover an image from the recorded data. In practice, the lateral object position relative to the modulator only needs to be known roughly to provide an initial estimate of the support position. The exact lateral position of the object wave relative to the modulator can be retrieved automatically by the phasing algorithm as long as a significant part of the wavefield falls inside the specified support. This is shown in Supplementary  Fig. 3b, where a part of the wave can still be reconstructed when the chosen support is laterally displaced. The resulting reconstruction gives a clear indication of where to reposition the support. In our experiments, the detector was in the far field of the modulator so the lateral position of the modulator relative to the detector did not affect the behaviour of the reconstruction algorithm. We have found the new CMI algorithm very robust to changes of the modulator's lateral position.
The CMI reconstruction of part of a test pattern nanofabricated in a 1.5 mm thick tungsten film by electron beam lithography and reactive ion etching is shown in Fig. 2c. The amplitude variations within the probe cause some parts of the sample to be illuminated with fewer photons, and thus affect the reconstruction quality of those areas. A flat-top illumination would therefore be preferable. A 3.2 mm diameter disk defined a loose support for the exit wavefield, and during the single-shot reconstruction the sample features started to appear after only a few tens of iterations from a random initial guess (Supplementary Video 1). A reconstruction using the looser support of a 3.7 mm diameter disk is shown in Supplementary Fig. 3a and is of comparable quality to Fig. 2c. A sequence of single-shot reconstructions covering a larger region of the test pattern is shown in Supplementary Video 2, where the phase of the illumination wavefield has been subtracted to reveal the phase change induced by the object.
We estimated the reconstruction resolution from the range of spatial frequencies of the retrieved wavefields that were stable for different initial guesses. In CMI the recorded data are related to the object via a modulator whose transmission function might not be known accurately. Resolution metrics that compare the computed estimate directly with the recorded data, such as the phase retrieval transfer function 21,39 , are not readily suitable. The red curves in Fig. 2d show the azimuthally integrated spectral amplitude (AISA) for 100 individually retrieved wavefields with different random starts. These wavefields were also averaged after their arbitrary phase offsets have been removed, and the blue curve shows its corresponding azimuthally integrated spectral amplitude. The blue curve deviates from the red ones at a spatial frequency, where the data noise starts to dominate the fine structure in the image. A 5% deviation indicated by the vertical line gives a resolution estimate of 37 nm.
In our experiments, the allowed wave extent on the modulator was 4.2 mm to fulfil the Nyquist sampling condition of the diffraction intensity in each dimension. This sets a limit on the sample area that can be imaged from a single measurement. In our experiments, however, a small fraction of the light was spread outside this area due to wave propagation effects from the pinhole to the modulator plane. To reduce the possible artefacts arising from this spreading, we performed reconstructions with a larger object window and applied the modulus constraint with suitable down-sampling and up-sampling operators (Methods).
Extension of the field of view. CMI can retrieve the sample exit wave from a single measurement or the illumination probe when no sample is present. If multiple measurements are taken from were determined retrospectively by aligning the single-shot reconstructions 40 , which avoided the positioning errors of the stepper motors. Some artefacts bearing the amplitude profile of the illumination probe can be seen, and these are due to small changes in the probe between each single-shot measurement. When the probe function is stable over a number of shots, a measurement of the probe before the measurements on the sample can be used to reduce such artefacts by dividing out the illumination from the exit wave to provide the sample transmission, provided that the probe amplitude is well above the noise level.
Tolerance to missing data. In CDI, the diffraction data missing from a measurement will result in unconstrained spatial modes in the phase retrieval that hinder convergence 17,18,21 .
The robustness of CMI to missing data was demonstrated by reconstructing measured diffraction data with its central region of 40 pixels in diameter being deliberately masked out, as shown in Fig. 3b. The speckles omitted from the centre were well retrieved as shown in Fig. 3c and are comparable to the recorded data in Fig. 3a. The quality of the retrieved image in Fig. 3d is comparable to that in Fig. 2c.
Reconstructions of a zone plate. In addition to the strongly scattering test sample, we performed CMI on a weaker sample of a zone-doubled Fresnel zone plate with buried structures 41 . The zone plate was made of hydrogen silsesquioxane resist on a Si 3 N 4 membrane, patterned by electron beam lithography and coated with a uniform and conformal thin film of iridium by atomic layer deposition. The zone plate has a diameter of 100 mm and an outermost zone width of 25 nm. A stitched phase image from a 39 Â 8 array of single-shot reconstructions of the zone plate sample, covering an area of 56 Â 11 mm 2 , is shown in Fig. 4a (the absorption image is shown in Supplementary Fig. 4). As one example, the single-shot reconstruction of the left circled region obtained using the loose support of a 3.2 mm disc is shown in Fig. 4b. Figure 4c is a close-up view of the boxed area in Fig. 4a. The iridium coating with a nominal thickness of 30 nm is clearly revealed on the steep sidewall of the photoresist, as shown in Fig. 4d. The buried hydrogen silsesquioxane resist pillars shown in the cross-section in Fig. 4d cannot be visualized by surfacesensitive techniques such as scanning electron microscopy, demonstrating one advantage of transmission X-ray imaging.

Discussion
We have shown a practical solution to the phase problem with a single exposure of the sample that applies to any form of coherent radiation and to samples of different diffracting strengths and physical sizes. Using a wavefront modulator we break the Fourier transform relationship between the object and the measured data. The method overcomes several challenges with conventional single-shot CDI. The CMI phase retrieval algorithm converges rapidly with only a loose support. In addition, we show the method can robustly recover significant areas of missing diffraction data, which can occur experimentally when beam stops are used or when the detector has a segmented structure with gaps between the segments. In addition, the use of a wavefront modulator can substantially reduce the intensity of the central beam and spread the photons more evenly over the detector, thus reducing the dynamic range requirement on the detector and consequently the need for a beam stop. The reduction of the data dynamic range is especially important for experiments on XFELs, where non-counting detectors of limited dynamic range, such as charge-coupled devices, are commonly used. Diffractive optics made from diamond and tungsten have been successfully tested and found to survive XFEL pulses 35,42 , and these materials could be used to form a radiation-resistant wavefront modulator for XFEL CMI applications. The single-shot capability of CMI can potentially become an important asset in circumventing the resolution limits imposed due to radiation damage 43,44 in X-ray microscopy, in the same way that the imaging of single biomolecules was first proposed 45 and later demonstrated 46 using XFELs. The development of single-shot diffraction imaging techniques in combination with the ultrashort pulses of highly brilliant XFEL sources is a promising method to record the diffraction pattern of biological samples before the radiation damage reveals itself. On one hand, it is complementary to current efforts to add cryogenic capabilities to existing ptychographic coherent diffraction imaging 47 and X-ray microscopy systems 48,49 to preserve the state of radiation-sensitive samples during data acquisition. On the other hand, the single-shot capability offers the possibility to study ultrafast dynamical processes with a temporal resolution limit only determined by the pulse structure of the light source and the rate at which the detector can acquire successive diffraction patterns.

Methods
Design considerations for the modulator. There are two considerations in the design of the transmission function of the modulator: to reduce the data dynamic range and to resolve ambiguities during phase retrieval. The exact form of modulator function is not critical and it is not restricted to a particular form such as to resemble a Fresnel zone plate. The resolution of the reconstructed image is not constrained by the feature size of the modulator. For ease of fabrication by electron beam lithography, a metal film with a binary pattern of open areas and phase-shifting areas was used in our experiment. The density of phase-shifting areas was designed to minimize the intensity of the zero-order component and hence to reduce the data's dynamic range. The angular spectrum of the light scattered by the modulator should fall within the angular acceptance of the detector. For the removal of ambiguities, a broader angular spectrum and a larger variation of the modulation function would discriminate the real solution from the ambiguous ones more effectively (Supplementary Fig. 1). In our simulation and experiments with various samples, the convergence rate of a reconstruction and the quality of the resulting images showed little dependence on the wavefront (sample) structures when the modulator phase varies more strongly than the wavefield that was retrieved. Our set-up and the design of the modulator have been found to work for general samples.
In our experiment, the modulator consisted of a silicon nitride support membrane coated with a tungsten film into which a randomly distributed pattern of 200 nm wide spaces was etched using electron beam lithography. Figure 1b shows an scanning electron microscope image of the modulator taken at a tilted angle to reveal the fabrication quality. The design thickness of the tungsten coating was 1.26 mm to induce a p radians phase shift at an X-ray energy of 6.2 keV. The pillar density was chosen to ensure their contribution to the zero-order diffracted amplitude was comparable to the contribution from the open spaces. The relative phase difference of these contributions was designed to be p radians, resulting in destructive interference that reduces the strength of the zero-order beam. In this way, the need for a beam stop was reduced and the dynamic range requirement of the detector was relaxed.
Characterization of the modulator. The phase and amplitude of the modulator transmission function were used as a priori information in our calculation and were measured by ptychography 32-34 before the CMI single-shot measurements. A total of 177 diffraction patterns over an area of 12 mm Â 12 mm of the modulator were recorded as the modulator was scanned across the X-ray beam using a two-axis piezo stage. The scan positions lay at the vertices of concentric polygons, starting with a pentagon. There was a radial increment of 0.8 mm and the number of vertices increased by five between successive rings, with an exposure time of 3 s at each scan point. The extended ptychographical iterative engine algorithm with position correction 50 was run on the recorded data for 120 iterations, yielding the results shown in Fig. 1c,d. The phase shift due to the tungsten pillars was measured to be 4.1 radians, 33% larger than the design value.
Reconstruction algorithm. The algorithm used here is an improved version of the one described in ref. 30. A new de-modulating approach has been used, which is suitable for a modulator with a strong absorption variation. The raised-power magnitude constraint 30 has been replaced by the approach described below in step (4) to prevent the algorithm from being trapped in local minima. The magnitude constraint in the detector plane has been revised to include up-sampling and downsampling at the detector plane as described below in step (5). The results presented in this paper were only possible with these improvements to the algorithm.
The wave estimates at the support, modulator and detector planes are specified by their coordinate subscript S, m or M and D, respectively. Normally we start with a random initial guess at the support plane and proceed as follows for the jth iteration, (i) Apply the support constraint to drive towards zero all values outside the support c j ðr S Þ ¼ĉ j À 1 ðr S ÞSðr S Þ þ bĉ j À 1 ðr S Þ À c j À 1 ðr S Þ 1 À Sðr S Þ ð Þ ð 1Þ where S(r S ) denotes the support; c j À 1 (r S ) is the previous estimate; and c j À 1 ðr S Þ is the estimate after revising by the measured data in the (j À 1)th iteration, as shown in Fig. 1e. The estimate outside the support is gradually driven towards zero at a rate controlled by the constant b; throughout the results reported here b ¼ 0.5 was used. (ii) Propagate c j (r S ) to the modulator, which yields c j (r m ). For the results presented in this paper, the angular spectrum method was used for propagation between the support plane and the modulator. (iii) Apply modulation where T (r m ) is the complex transmission function of the modulator. The subscripts m and M represent planes directly upstream or downstream of the modulator. (iv) Generate the estimate at the detector plane according to where P z { Á } represents the wave propagator for a distance z from the modulator to the detector and was a Fourier transform in our case. Every N iterations, where N ¼ 40 in the paper, c j À 1 (r M ) is set to equal c j (r M ), which appeared to help to avoid stagnation. (v) Apply the magnitude constraint To accommodate some energy outside of the object space area defined by the Nyquist sampling condition at the detector, we use estimates of the wavefield with a larger array size than the recorded diffraction pattern. A larger array size in the calculation accounts for some degree of under-sampling in the diffraction data and helps reduce artefacts. This has been found beneficial in improving the image quality of our results. After calculating the diffracted wave c j (r D ) at the detector, we down-sampled its amplitude to match the sampling interval of the diffraction pattern. First, the amplitude of the calculated wave c j (r D ) is convolved with a detector pixel form function and then down-sampled to the array size of the recorded data, yielding |c j (r D )| ds . A constant function of 2 Â 2-pixel wide was used to approximate the detector pixel form function for our results. Next, we calculated the weighting function b ds (r D ) of the measured data with respect to the down-sampled amplitude, according to b ds ðr D Þ ¼ ffiffiffiffiffiffiffiffiffiffi ffi and then up-sampled the weighting function, here using nearest-neighbour interpolation, yielding b (r D ); where I(r D ) is the measured intensity in the detector plane. Finally, the revised wavefield at the detector is calculated aŝ where a is a constant. For the results in this paper a ¼ 1. (viii) Back-propagateĉ j ðr m Þ to the support plane, yielding the revised wave field c j ðr S Þ.
Steps 1-8 are repeated until a termination condition is met.
The algorithm for wavefield propagation in step (2) between the support plane and the modulator, and in step (4) between the modulator and detector, should be adapted to the experiment geometry. For near-field propagation, the angular spectrum method is used; the Fresnel propagator is more suitable for the case of intermediate-range propagation 51 .
Data availability. All the relevant data are available from the authors on request.