Extreme ultraviolet (EUV) and soft X-ray microscopy techniques, which can exploit element-specific absorption contrast to form high-resolution chemically resolved imaging, have garnered wide applications for nano-imaging^[1,2]. Furthermore, many modern EUV sources, featuring ultrashort pulses, are highly suitable for investigating transient phenomena on the femtosecond and attosecond timescales^[3]. Due to the challenging fabrication of high-quality lenses in the EUV or X-ray regimes, coherent diffraction imaging (CDI), which is a lensless imaging technique, has been proposed and widely applied in synchrotron facilities and free-electron lasers^[4–6]. CDI can recover the sample image from the measured diffraction pattern via phase-retrieval algorithms. Recently, with the development of tabletop X-ray and EUV sources, X-ray CDI has been shown to be feasible on laboratory-scale sources^[7].

Due to its high coherent photon flux, high-harmonic generation (HHG) has been successfully used for ptychography^[2,8–10], a recent extension of CDI. In ptychography, a set of diffraction patterns through scanning the sample area by area ensures overlap between adjacent scans. The redundant information provided by scanning with overlap enables ptychography to quickly converge and obtain information on complex objects and illumination probes. Thanks to the robustness of ptychography, ptychographic information multiplexing^[11] (PIM) has been proposed to recover the spectral response of a sample at different wavelengths simultaneously; when combined with high-harmonic comb sources, it can be used to obtain images with elemental contrasts^[12]. Building on PIM, a universal ptychographic algorithm utilizing automatic differentiation significantly advances spectral multiplexing with EUV sources, enabling high-resolution and chemically specific imaging without wavelength scanning^[13]. However, the scanning requirement of ptychography limits its ability for in situ observation of dynamic samples, especially when studying the phase transition dynamics of electrode materials^[1]. Hence, CDI with single-frame and multispectral imaging capabilities is still indispensable. Several broadband CDI algorithms such as multiwavelength CDI^[14], polyCDI^[15], and numerical monochromatization^[16–18] have been developed for polychromatic data. These methods mainly solve the problem of low temporal coherence through numerical algorithms, thereby improving the photon flux and reducing exposure time instead of obtaining images of each wavelength. Recently, a multiwavelength phase-retrieval algorithm has demonstrated that CDI can realize reconstruction with multiwavelength illumination through numerical simulations^[19], which utilizes the wavelength dependence of the support to separate and recover each spectral component. However, to achieve this method on tabletop EUV sources, further optimization of the algorithm remains necessary.

In this paper, we propose a single-frame multiwavelength CDI (mw-CDI) method with a novel iterative projection model that can reconstruct each spectral response of the sample with multiple harmonics well, so that we can perform chemically resolved imaging of dynamically changing samples. Our algorithm is shown to provide better image quality for each wavelength compared to traditional methods. The proof-of-principle experiment was carried out with an HHG EUV source containing three harmonic peaks (namely, 34.3, 39.6, and 46.8 nm).

A basic CDI algorithm conducts iterations between two planes: a sample plane and a detector plane. Our mw-CDI follows a similar flow chart, as outlined in Fig. 1. A coherent mode decomposition is adopted to form a modified phase-retrieval algorithm, allowing exit waves at different wavelengths to be recovered simultaneously. When sources contain several distinct wavelengths (${\lambda }_n$, $n=1,2,\dots ,N$), there will be $N$ coherent modes with different wavelengths propagating from the sample plane to the detector plane. The multiwavelength diffraction intensity can be described by the incoherent sum of far-field propagated waves ${\psi }_n$, $$I(\mathit{q})=\sum _{n=1}^N{S}_n{|{\psi }_n(\mathit{q})|}^2,$$(1)where $\mathit{q}$ is a position vector in reciprocal space and $S_n$ are the spectral weightings measured beforehand.

Normally, CDI can start the reconstruction by multiplying the initial support guess and object guess. In mw-CDI, the support domains, $A_n^k$, for each coherent mode are wavelength-dependent scaling. Therefore, a primary wavelength component will be selected, and its corresponding support domain will be scaled in pixel size to generate support domains for other components while maintaining the same number of computed pixels. Hence, the spatial constraint using the support domains at the $k$th iteration is given by $$\begin{gathered} {\phi }_n^k(\mathit{r})={\widehat{\phi }}_n^{k-1}(\mathit{r})A_n^k(\mathit{r})+\beta [{\widehat{\phi }}_n^{k-1}(\mathit{r})-{\phi }_n^{k-1}(\mathit{r})] \\ ·[1-A_n^k(\mathit{r})]·\mathit{a}, \end{gathered}$$(2)and $A_n^{k=1}(\mathit{r})$ is the circular-shaped support function, $$A_n^{k=1}(\mathit{r})=\{\begin{array}{cc}1& \mathit{r}\in \text{support region}\\ 0& \text{others}\end{array},$$(3)where $\mathit{r}$ is a sample space coordinate. ${\phi }_n^k(\mathit{r})$ is the complex scalar wavefield at the $n$th wavelength and ${\widehat{\phi }}_n^{k-1}(\mathit{r})$ denotes the revised estimate wave obtained in the ($k-1$)th iteration. $\beta$ is an update weight factor that is varied within [0,1], which can control how fast the values outside the support region are set to zero. In this work, a value of $\beta =0.5$ was used. $\mathit{a}$ serves as a switch between the error reduction (ER) and hybrid-output (HIO) algorithms. When its value is 0, the ER algorithm^[20] is employed, while the HIO algorithm^[21] is used for 1. The ER and HIO algorithms are alternately iterated 20 and 50 times, respectively.

Then, the exit waves of the sample propagate to the detector in the far field so that the estimated wavefront at the detector plane can be given by $${\psi }_n^k(\mathit{q})=\mathcal{F}[2{\phi }_n^k(\mathit{r})-{\phi }_n^{k-1}(\mathit{r})],$$(4)where $\mathcal{F}$ represents the Fourier transform operator. The difference map (DM) algorithm^[22] is used and every 15 iterations ${\phi }_n^{k-1}(\mathit{r})$ is replaced by ${\phi }_n^k(\mathit{r})$ to avoid stagnation. The new iterative projection model of mw-CDI, a combination of ER, HIO, and DM algorithms, significantly improves the convergence.

The primary distinction between the mw-CDI and the monochromatic method lies in the necessity of incoherent superposition computed diffraction patterns from different wavelengths when applying the modulus constraint. Therefore, the revised diffracted wave ${\widehat{\psi }}_n^k(\mathit{q})$ can be obtained by $${\widehat{\psi }}_n^k(\mathit{q})={\psi }_n^k(\mathit{q})\frac{\sqrt{I_m}}{\sqrt{{\sum }_{n=1}^N{S}_n^k{|{\psi }_n^k(\mathit{q})|}^2}+\epsilon },$$(5)where $I_m$ is the measured multiwavelength diffraction intensity and $\epsilon$ is a minute value employed to prevent division by zero. The accurate values of $S_n$ are important and can be modified by the properties of the sample. Using the known spectral weight as a starting point, we update $S_n^k$ after the modulus constraint from the intensity percentage of each spectral mode, $${\widehat{S}}_n^k=\frac{{\sum }_{\mathit{q}}{|{\widehat{\psi }}_n^k(\mathit{q})|}^2}{{\sum }_n[{\sum }_{\mathit{q}}{|{\widehat{\psi }}_n^k(\mathit{q})|}^2]}.$$(6)

Then, we propagate the revised wavefront to the sample plane and jump to the spatial constraint step to start the next round of iteration. The support domains will be updated using a multiwavelength version of the shrink wrap method^[23] before applying the spatial constraint. The support domain corresponding to the primary wavelength is obtained by thresholding a Gaussian-smoothed version of ${\widehat{\phi }}_n^k(\mathit{r})$, followed by scaling this support with the ratio of wavelengths to obtain the remaining support domains. In our experiment, the Gaussian kernel array was $151\mathrm{pixel}\times 151\mathrm{pixel}$ in size with a standard deviation of 0.2 and a threshold of 0.08.

Our experimental setup for EUV CDI is shown in Fig. 2, where an HHG source, containing three harmonics of 34.3, 39.6, and 46.8 nm, was used as the light source [the spectrum is shown in Fig. 2(b)]. To produce the higher flux of harmonics around 26 eV, the 0.45 mJ, 180 fs laser pulses with a central wavelength of 515 nm are obtained by bunching a commercial femtosecond laser with a center wavelength of 1030 nm (Light Conversion, Pharos) to about 5 mm and passing through a 0.5 mm thick BBO crystal. The beam then passes through two dichroic mirrors to remove the original 1030 nm laser, and the remaining pure 515 nm driving light source is spatially filtered by a half-wave plate and an iris, in order to approach the highest HHG conversion efficiency and the highest EUV spot quality. The intensity characteristics of the three harmonic beams are shown in Fig. 2(c). The spot diameter (FWHM) measured at 1 m is $\sim 2.5\mathrm{mm}$, thus determining that the full angle of HHG beam divergence is about 2.5 mrad. The laser pulses are focused by a lens with a focal length of 300 mm into a gas cell with a hole of 200 µm, resulting in an estimated laser peak power of $2.93\times {10}^{14}\mathrm{W}/{\mathrm{cm}}^2$. Our gas cell is a tube with an inner diameter of 2.0 mm, which can get the highest HHG conversion efficiency with our test parameters. Subsequently, two 250 nm thick aluminum films are used to remove the residual driving laser mixed in the high-order harmonics. Hence, the tabletop EUV source we utilized is compact and cost-effective, which can facilitate the widespread application of EUV CDI.

Our imaging experimental schematic is illustrated in Fig. 2(a). Before illuminating the sample, the EUV beam is focused by a toroidal mirror and spatially filtered through a pinhole to create an illumination probe. Because of the poor convergence of CDI, the pinhole surface is typically placed in complete contact with the sample plane, thereby generating a sharp probe with a diameter of 30 µm on the sample plane to facilitate convergence. An optical microscope image of the sample is shown in Fig. 3(a). The sample featuring a maze pattern was fabricated through two-photon polymerization 3D printing techniques (Nanoscribe), while the material composition is a photosensitive resin known as IP-Dip. It exhibits structures with a linewidth of 1.5 µm, gaps with a linewidth of 2.5 µm, and a thickness of 5 µm, while being completely drilled through. The detector with $11\mathrm{\mu m}\times 11\mathrm{\mu m}$ in pixel size (XF95, Tucsen) was placed downstream of the sample at a distance of 40 mm. To sufficiently capture high-frequency information, a detector area of $1024\mathrm{pixel}\times 1024\mathrm{pixel}$ was used in our experiment. The obtained diffraction pattern, with an exposure time of only 8 ms, is shown in Fig. 3(b). According to the experimental setup, our oversampling ratio reached 30, which facilitates obtaining high-quality reconstructed images^[24].

In the processing of multiwavelength experimental data, the iterative projection algorithm combining ER, HIO, and DM played a significant role in facilitating the convergence of mw-CDI. To compare the advantages of our new method, we first employed a multiwavelength CDI using a combined ER and HIO algorithm to process the data, which was previously used to process simulated and experimental data of dual-wavelength EUV^[19,25]. Figures 3(c) and 3(e) display the reconstruction corresponding to wavelengths of 34.3, 39.6, and 46.8 nm, respectively. All wavelength components failed to achieve satisfactory reconstruction; only some lines can be retrieved. However, when the new combined approach was applied to mw-CDI, the amplitude of sample at all three wavelengths was successfully reconstructed [shown in Figs. 3(f)–3(h)]. The reconstructed sample structure and linewidths closely match the measurement result obtained using optical microscopy. Furthermore, by calculating the mean squared error (MSE) between the diffraction patterns, our method achieved an MSE value of 0.003, whereas the traditional method had a value of 0.14, demonstrating the superior quality of our approach. The detailed comparison of the diffraction patterns is presented in Fig. S1 in the Supplement 1. While theoretically shorter wavelengths result in higher reconstruction resolution, the reconstruction quality for the wavelength of 46.8 nm is superior among the three harmonics, as it possesses the highest energy proportion. In contrast, the harmonic at 34.3 nm, which has a minimal energy proportion, exhibits the poorest reconstruction quality. This indicates that in multiwavelength reconstruction, the energy proportion of each wavelength will have a significant impact on the quality of the results. To enhance the reconstruction quality, we also employed a multi-exposure acquisition strategy to obtain diffraction patterns with a higher dynamic range. The reconstruction is shown in Fig. S2 in the Supplement 1, compared to single-frame, the multi-exposure results exhibit significant improvements in quality. The detailed descriptions of the multi-exposure experiment are given in the Supplement 1. We also incorporated several commonly used iterative projection algorithms into mw-CDI. Comparative studies using both experimental and simulated data show our novel method offers superior performance; details can be found in the Supplement 1.

To demonstrate our mw-CDI method combined with HHG EUV can achieve chemical-specific imaging, we conducted simulation using a model metal mixture at a thickness of 40 nm as the sample. The EUV harmonics, whose wavelengths are 23.95 and 25.12 nm, are used as illumination to create the polychromatic diffraction pattern [Fig. 4(a)]. Notably, the absorption characteristics of the sample resulted in a modification of the spectral weights, shifting from [0.5, 0.5] before the sample to [0.522, 0.478] at the exit surface. Here, we used a photon count of $2\times {10}^9$ and introduced Poisson noise into the diffraction pattern. The simulated experimental geometry and parameters are essentially consistent with those used in our experiment above, while only the probe size changed to 15 µm to satisfy the oversampling ratio requirement. The model metal mixture sample consisted of Mg and Al, whose refractive indices and extinction coefficients according to Refs. [26,27] are shown in Figs. 4(b) and 4(c). The refractive indices of the two metals are relatively similar; however, there is a notable difference in their extinction coefficients: Mg exhibits an absorption edge in the wavelength of 24.7 nm. The distribution of the two metals is shown in Fig. 4(d), with the dark regions representing Mg, and the light regions representing Al. The absorption rates of the two elements at a wavelength of 25.12 nm are similar, making it difficult to distinguish between them [Fig. 4(e)]. However, at a wavelength of 23.95 nm, due to the significant difference in absorption rates, the amplitude image exhibits high contrast, clearly revealing the distribution of the two elements [Fig. 4(d)]. We conducted mw-CDI reconstruction on this simulated sample and successfully recovered the results from the polychromatic diffraction pattern, which matched the sample very well as shown in Figs. 4(f) and 4(g). Furthermore, the spectral weights of the sample exit surface recovered during our reconstruction process were found to be [0.526, 0.474], which is consistent with the theoretical values above. Therefore, our method enables element-specific imaging without wavelength scanning.

In summary, we have proposed a mw-CDI method that enables simultaneous reconstruction at multiple wavelengths, and we have achieved single-frame multiwavelength EUV CDI using HHG for the first time, to the best of our knowledge. Our method incorporates a combination algorithm of ER, HIO, and DM to process the multiwavelength diffraction data. Compared to previous EUV dual-wavelength CDI, our method exhibits a significant improvement in reconstruction quality at three harmonics. Additionally, we showed a simulation with two harmonics around the Mg and Al absorption edge, to demonstrate the feasibility of single-frame chemical mapping. We believe that our method, which can enable in situ reconstruction of element-specific information, will further promote the widespread application of tabletop EUV sources in the imaging of materials.