PhaseT3M: 3D imaging at 1.6 Å resolution via electron cryo-tomography with nonlinear phase retrieval

Domain: physics
Published: 2025-12-17
Editorial outcome: major_revision

Abstract
Electron cryo-tomography (cryo-ET) enables 3D imaging of complex, radiation-sensitive structures with molecular detail. However, image contrast from the interference of scattered electrons is nonlinear with atomic density and multiple scattering further complicates interpretation. These effects degrade resolution, particularly in conventional reconstruction algorithms, which assume linearity. Particle averaging can reduce such issues but is unsuitable for heterogeneous or dynamic samples ubiquitous in biology, chemistry, and materials sciences. Here, we develop a phase retrieval-based cryo-ET method, PhaseT3M. We experimentally demonstrate its application to an approximately 7 nm Co3O4 nanoparticle on an approximately 30 nm carbon substrate, achieving a maximum resolution of 1.6 Å, surpassing conventional limits using standard cryo-TEM equipment. PhaseT3M uses a multislice model for multiple scattering and Bayesian optimization for alignment and computational aberration correction, with a positivity constraint to recover ‘missing wedge’ information. Applied directly to biological particles, it enhances reconstruction quality and reduces artifacts, establishing a standard for routine 3D imaging with phase contrast.

Introduction
Transmission electron microscopy (TEM) enables the formation of two-dimensional (2D) structural images from atomic-to-micrometer scales, making it a valuable characterization tool for materials of nearly any composition. However, critical structural details such as defects, interfaces, complex architectures, and weak electron signals are hidden by the projection of the three-dimensional (3D) material potential to 2D1–5. Electron tomography (ET) addresses this limitation by combining images from multiple angles to reconstruct 3D volumes from 2D projections. This approach reveals 3D features that are inaccessible with 2D imaging alone1,6–11, yet are essential for determining material properties. Recent advances in ET have pushed its resolution to the atomic scale through the use of aberration-corrected annular dark-field (ADF) scanning transmission electron microscopy (STEM)1,5,7,10,12–17. Despite the advances, ADF-STEM suffers from inherent limitations, including nonlinear contrast due to multiple scattering, a restricted depth of focus, poor sensitivity to light elements, and high electron dose requirements due to slow scanning time. These drawbacks make it less effective for 3D imaging of thick or radiation-sensitive materials18–24.

In contrast to scanning techniques, conventional TEM (CTEM) uses a parallel electron beam, allowing for rapid, low-dose imaging of dynamic processes and radiation-sensitive materials8,11. CTEM image contrast arises from the interference of scattered electrons from the electrostatic potential. Similar to STEM, stronger electron-atom interactions and increased sample thickness result in more multiple scattering, leading to more nonlinear interference contrast compared to atomic density. Moreover, multiple scattering in samples thicker than approximately twice the mean free electron path25 degrades beam coherence and therefore limits resolution. Additional contrast modulation is created by lenses in the imaging system. These effects hinder accurate structural reconstruction using conventional tomography techniques that assume linear image contrast8,11,23,24, except in special cases26,27 where electrons scatter only once from a single atom type. The special cases can achieve atomic resolution but require thousands of aberration-corrected high-resolution TEM (HRTEM) images and high electron doses, and are limited to tiny nanoparticles, making them impractical for larger or beam-sensitive samples26,27.

Weak electron-atom interactions in organic materials allow their structures to be determined from single-scattering reconstructions8,11, but their radiation sensitivity prevents atomic details from being resolved in individual particles, and their image contrast is not strictly linear. Advanced 3D imaging techniques, such as single-particle cryo-electron microscopy (cryo-EM), reconstruct the average 3D atomic positions of particles as small as protons with near-atomic resolution28–33, overcoming radiation sensitivity limits. However, these techniques are limited to materials with thousands of nearly identical copies that can be averaged together, which precludes their application to heterogeneous materials.

Interpretable 3D imaging of arbitrary materials requires a physical model to solve the electron-atom scattering problem and instrumental contrast transfer modulation to overcome their limitations, while also being robust to noise arising from sensitive or thick samples. Emerging 3D phase retrieval algorithms based on the multislice method address these limitations by solving the inverse scattering problems and enabling reliable 3D reconstruction of thick samples under low electron doses23,24,34–38. STEM-based 3D phase retrieval has recently made major advances in retrieving 3D phase information with improved resolution and sensitivity to light elements39–42. However, implementing these methods broadly in practice remains challenging, especially for radiation-sensitive materials, due to beam damage from the intense probe, long acquisition times for large areas, and distortions from scan instabilities. In contrast, HRTEM in CTEM enables fast imaging over large fields of view. Combined with fractional imaging and motion correction, it can achieve high-resolution without significant loss, even for radiation-sensitive materials43–45.

Here, we experimentally demonstrate an HRTEM tomography method based on phase retrieval, PhaseT3M, showing that 3D phase-contrast imaging can achieve a maximum resolution of 1.6 Å with a standard cryo-electron microscope without physical aberration correction. We reconstruct a 7 nm Co3O4 nanoparticle supported on 30 nm of carbon using the multislice model to account for multiple scattering and Bayesian optimization to computationally correct image and tilt alignment, contrast transfer, and microscope aberrations. Using the same method, we also reconstruct human immunodeficiency virus (HIV-1) particles from the publically available EMPIAR-10164 dataset30, which represents higher Fourier ring correlation (FRC) values and 20–50% R-factor improvement compared to conventional reconstruction methods. Our method is applicable to both organic and inorganic materials, vitrified in ice or under vacuum. Autonomous data acquisition and reconstruction scripts are widely available and easily implemented on standard commercial microscopes, making this a powerful scientific tool that spans materials, biological, physical, chemical, and geological sciences.

Results
Tomography experiment setup
Figure 1a–c shows the experiment setup of our 3D phase-contrast tomography, PhaseT3M. We collected a tomographic tilt series of Co3O4 nanoparticles supported by a carbon grid, by acquiring HRTEM images at tilt angles ranging from −65∘ to +65∘ in 1∘ increments (see Supplementary Figs. 1–3). At each tilt angle, three images were acquired at different defocus values. The electron dose per image was 21.5 e−/Å2, resulting in a total accumulated dose of 8,460 e−/Å2 across the entire series. As shown in Fig. 1d–f, we performed image preprocessing for 3D reconstruction. These steps include motion correction, focal and tilt alignment, estimation of microscopy parameters, and denoising (see “Methods”). From the preprocessed tilt series, we reconstructed a 3D potential of the sample, which includes an approximately 7 nm diameter Co3O4 nanoparticle and a 30 nm thick amorphous carbon support. To retrieve 3D phase structural information and account for multiple scattering, we applied a forward and inverse multislice approach36,46 (see “Methods”), as illustrated in Fig. 1g. We iteratively reconstructed the 3D potential using the gradient descent method, minimizing the differences between the calculated projections from the estimated 3D potential and the experimentally measured projections.

Experimental tomography via phase retrieval
The reconstructed volume in Fig. 2a represents a 3D potential of the 7 nm Co3O4 nanoparticle partially embedded in the 30 nm carbon support. The intensity of the reconstructed potential is approximately proportional to the atomic number: lighter elements like carbon and oxygen produce lower potential signal, while heavier elements like cobalt show a higher potential magnitude. As seen in Fig. 2a, the Co3O4 nanoparticle and the carbon support are clearly distinguishable.

To estimate the quality and resolution of the 3D reconstruction, we reconstructed a 3D potential using a cropped tilt series focusing on the Co3O4 nanoparticle with a pixel size of 0.52 Å (see “Methods”). We then applied a 3D mask to the reconstructed volume in Fig. 2b to remove the intensity of the carbon support film. As shown in Fig. 2c, the central slices of the 3D reconstructed volume along different axes reveal bright atomic intensities and lattice lines of the nanoparticle. The 2D Fourier transform of each central slice in Fig. 2d shows the {400} diffraction peak at a resolution of 2 Å, except for the slice where the missing wedge region overlaps with one of the {400} diffraction peaks. The 3D Fourier transform of the 3D reconstructed image in Fig. 2e also exhibits diffraction peaks at 2 Å resolution, consistent with the 2D Fourier transform results. Figure 2f shows a 3D isosurface image with a lower threshold than that in Fig. 2e, presented from a different viewpoint to highlight the missing wedge region. Interestingly, the diffraction peaks of 4.6 Å {111} and 2.8 Å {220} resolution in the missing wedge region were restored. This improvement is attributed to the application of a physical constraint, positivity (i.e., enforcing non-negative potential values), which helps suppress missing wedge artifacts and the elongation effect, ultimately aiding in the recovery of diffraction peaks. As shown in Supplementary Fig. 4, only the reconstruction with the positivity constraint exhibits the recovery of diffraction peaks within the missing wedge. Additionally, {422} diffraction peaks at 1.6 Å resolution in Fig. 2e, f are visible, indicating the high-resolution capability of our reconstruction. This suggests that the maximum resolution of our 3D reconstruction is 1.6 Å. However, due to the limited tilt range and single, fixed rotation axis, the resolution is anisotropic, varying depending on direction. We also plotted the power spectrum in Fig. 2g to cross-check the resolution estimation, and the results are consistent with the previous Fourier transform analysis.

To investigate the effect of multiple scattering correction using the inverse multislice method, we compared slices from volumes reconstructed with the multislice and single-slice phase retrieval approaches. As seen in Supplementary Fig. 5, multiple scattering correction via the multislice method improves the resolution. Furthermore, Supplementary Fig. 6 shows that the multislice-based method yields a lower reconstruction error compared to the single-slice approach. Additionally, as shown in Supplementary Fig. 7, the multislice method can account for defocus gradients arising from variations in sample height and tilting.

To evaluate the resolution capability of our reconstruction method under different dose conditions, we leveraged our raw image stacks, which consist of 24 subframes per tilt and defocus, allowing us to analyze varying total doses by selecting subsets of these subframes. Reconstructions were performed at total doses of 8460 (24 subframes), 4230 (12 subframes), 705 (2 subframes), 353 (1 subframe), and 118 (one defocus image with one subframe) e−/Å2. As shown in Supplementary Fig. 8, the maximum achievable resolution systematically decreases as the total dose decreases. Notably, at the lowest dose of 118 e−/Å2, some diffraction peaks (corresponding to 5.8 Å in the 2D Fourier transform and 2.8 Å in the 3D Fourier transform) remain visible, but the lattice structure cannot be resolved. These results suggest that while low-dose conditions can retain limited high-resolution information, achieving near-atomic resolution beyond lattice-level resolution requires higher electron doses or averaging across multiple particles to increase the effective signal-to-noise ratio.

Simulated tomography via phase retrieval
To verify the resolution estimated from our experimental reconstruction, we conducted simulations to reproduce our phase retrieval-based tomography experiment. By comparing the experimental and simulated results, we confirmed the validity of our resolution estimation. As shown in Fig. 3a, we first created 3D atomic positions of a Co3O4 nanoparticle and an amorphous carbon membrane (see “Methods”). Next, we performed TEM simulations to generate the simulated tilt and focal series under the same experimental conditions. We applied Poisson noise to match the electron dose level used in the experiment (see “Methods”). Finally, we reconstructed a 3D potential from the simulated tilt series using the same reconstruction process as in the experiment (see “Methods”). We then applied a 3D mask to the reconstructed volume to extract the nanoparticle intensity. Additionally, for using to provide a reference for the reconstruction, we generated 3D atomic potentials from the generated atomic model (see “Methods”).

Following the same resolution analysis as performed on the experimental results, we evaluated the simulated reconstruction. Figure 3b shows the 3D density map of the reconstructed volume of the Co3O4 nanoparticle, revealing clear lattice fringes. The 3D isosurface plots of the Fourier transform images in Fig. 3c exhibit diffraction peaks at 1.6 Å resolution and recover diffraction peaks within the missing wedge region, consistent with the experimental results. Additionally, faint peaks are observed at 1.4 Å resolution, which is beyond the maximum resolution achieved in the experiment.

The central slices of the 3D volume in Fig. 3d exhibit bright atomic intensities with clear lattice fringes. The simulated atomic structure in Fig. 3d, e aligns well with the intensity maxima in the central slice, with small deviations arising from experimental limitations such as restricted electron dose, elongation from the missing wedge, the thick carbon background, and residual alignment errors. For reference, the ground truth 3D atomic potential in Fig. 3e, g shows perfect agreement with the atomic positions, along with very sharp diffraction peaks up to the Nyquist limit. As seen in Fig. 3f, the {400} diffraction peaks at a resolution of 2 Å are visible in the Fourier transform images, except within the missing wedge, similar to the experimental results (Fig. 2d). However, simulation results show much sharper diffraction peaks and additional peaks of 1.8 Å beyond the 2.0 Å resolution limit of the experiment. This difference is due to imperfections in alignment, parameter estimation, and noise modeling, as well as limitations of the tomography method, which cannot fully account for all experimental effects, such as inelastic scattering. It is well known that simulation results generally outperform experimental ones because they do not consider all sources of systematic errors and artifacts present in experimental data47. Although the simulation performs slightly better, as expected, the simulated 3D reconstruction, which mimics our experimental conditions, demonstrates results comparable to the experimental data in terms of resolution and diffraction peak retrieval within the missing wedge region.

Additionally, to test the validity of the phase retrieval method, we examined its performance at different electron doses: 8460 e−/Å2 (used in the experiment), 84,600 e−/Å2, and infinite dose. For simplicity, we used a clean model of Co3O4 nanoparticle without alignment errors and carbon background, including only the missing wedge and electron dose effects. As shown in Supplementary Fig. 9, the reconstruction at infinite dose is in agreement with the ground truth, while the reconstruction quality progressively degrades as the dose decreases. Note that applying a low-pass filter of 0.8 Å−1, together with missing wedge effects, hinders the capability to fully resolve the oxygen intensities in the reconstructions. This demonstrates the validity of our method, but also highlights its practical limit under experimental conditions.

Comparative analysis of phase retrieval and conventional tomography
To compare our approach with conventional tomography, we reconstructed the 3D volume of the Co3O4 nanoparticle and its carbon support using the projection-based tomography algorithm, specifically the simultaneous iterative reconstruction technique (SIRT) (see “Methods”)48. As shown in Fig. 4a, e and Supplementary Fig. 5a, d, both the multislice and single-slice phase retrieval methods can retrieve not only the nanoparticle but also the carbon support, whereas conventional tomography fails to reconstruct the carbon support. This suggests that the ability to solve the inverse scattering problem is the critical factor for reconstructing the carbon support, rather than multiple scattering correction. In Supplementary Fig. 10c, d, the carbon intensity is even lower than the vacuum area, and ring-shaped artifacts appear prominently at the boundaries due to nonlinear complex interference contrast, which violates the tomography projection rule.

To further evaluate the performance of the two tomography methods, we compared cropped regions close to the nanoparticle. As shown in Fig. 4b, f, both methods enable the reconstruction of the nanoparticle structures. However, phase retrieval provides higher overall quality and resolution than conventional tomography. This is supported by the Fourier images in Fig. 4d, h and Supplementary Fig. 11c, f, which show higher index diffraction peaks for phase-contrast imaging. Specifically, Fig. 2 demonstrates that our phase retrieval imaging achieves a resolution of 1.6–2.0 Å, exceeding the typical limits of the conventional tomographic method and the diffraction-limited resolution under our experimental conditions. In comparison, Supplementary Fig. 12c–g shows that SIRT reconstruction yields a lower resolution of 2.0–2.4 Å, which is similar to the first zero-crossing point of 2.4 Å indicated by the 2D CTF curve in Supplementary Fig. 13). Additionally, the brighter intensities observed at the nanoparticle boundaries in Fig. 4g and Supplementary Fig. 11e indicate a deficiency in low-frequency information in the SIRT reconstruction. In contrast, the phase retrieval-based reconstruction compensates for these anisotropic intensities by effectively retrieving the missing frequency information through the lens transfer characteristics. This capability can be demonstrated by simulations of the spectral signal-to-noise ratio (SSNR) in 2D focal series reconstruction, presented in Supplementary Fig. 14. The SSNR indicates that information can be retrieved even at the zero-crossing points of the CTF in raw HRTEM images, where contrast transfer is otherwise absent (see “Methods”).

HIV-1 particles reconstruction
A major advantage of our method is that it can be directly applied to typical HRTEM tomographic tilt data, including open-sourced datasets such as EMPIAR. Notably, while the previous example used tomographic tilt data acquired with multiple defocus values, our method is also applicable to tilt series collected at a single defocus, enabling the retrieval of the full 3D phase volume.

To demonstrate the flexibility of our 3D method, we applied it to a biological sample, HIV-1 particles, using the publicly available EMPIAR-10164 dataset30. We first prepared a single HIV-1 tomographic tilt-series dataset (Supplementary Fig. 15). We then reconstructed the 3D volumes of HIV-1 particles using the phase retrieval and conventional tomography methods for comparison. Figure 5 presents the 3D volumes obtained from our phase retrieval method and the SIRT algorithm48, clearly showing the ring-shaped structures of HIV-1 particles and bright fiducial gold markers. As observed with the Co3O4 nanoparticle reconstructions (Fig. 4e), the conventional tomography reconstructions (Fig. 5b and Supplementary Fig. 16c–h) exhibit noticeable ring-shaped artifacts near the edges because the interference contrast does not strictly satisfy the assumptions of the tomographic projection theorem. In the enlarged views of the reconstructions (Fig. 5e–h and Supplementary Fig. 16), the overall structural features appear similar across different methods. However, conventional tomography reconstructions exhibit pronounced edge artifacts. Notably, iterative-based methods such as SIRT and real space iterative reconstruction (RESIRE) produce smoother intensities due to the suppression of high-frequency noise through iterative updates. In contrast, non-iterative reconstructions, such as filtered back projection (FBP) and weighted back projection (WBP), display sharper intensity features and reduced edge artifacts. However, this apparent improvement can be misleading, as it mainly results from the suppression of low-frequency components by the reconstruction filter, leading to an incomplete recovery of the true structure despite the visually enhanced high-frequency details. This interpretation is supported by our iterative method results using the same reconstruction filter as in FBP and WBP, which produced similar sharpening effects and low-frequency artifact reduction, as shown in Supplementary Fig. 17.

To evaluate the performance of our approach, we conducted a leave-one-out cross-validation49 and calculated the R-factor16,17. For the leave-one-out test, we computed FRC curves between the experimental zero-degree projection and the calculated projection from the reconstructed volume, excluding the zero-degree tilt image (see “Methods”). This approach is particularly effective for assessing biological specimens where the number of tilt projections is limited. As shown in Fig. 5c and Supplementary Fig. 18, both the multislice and single-slice phase retrieval methods achieve higher correlations compared to conventional tomography algorithms. Notably, the low correlation observed in the low-frequency range of the FRC for FBP and WBP explain their misleading appearance of sharpness, as these methods tend to preserve high-frequency components while suppressing low-frequency information during reconstruction. We also calculated the R-factor between the measured projections and the corresponding forward projections from the reconstructed volumes (see “Methods”). The multislice phase retrieval method achieved the lowest average R-factor of 3.5%, compared to 4.7% for SIRT, 5.8% for RESIRE, 10.1% for FBP, 10.1% for WBP, 5.3% for simultaneous algebraic reconstruction technique (SART), and 5.2% for algebraic reconstruction technique (ART) (see Supplementary Fig. 19). These results indicate that our phase retrieval method demonstrates an enhanced capability to recover missing structural information and provides superior agreement with the experimental projections.

We note that these HIV-1 results are based on a single tilt series, leading to a significantly lower resolution than the previously reported values achieved through subtomogram averaging30. Furthermore, the tomographic tilt series was acquired with a total electron dose of approximately 100 e−/Å2, which inherently limits the achievable resolution compared to the 1.6 Å resolution obtained for the Co3O4 nanoparticle, where two orders of magnitude higher dose was used. It is also important to consider that organic materials like HIV-1 are composed of light elements, which generally exhibit weak electron-atom scattering. This is consistent with our observation that the multislice and single-slice approaches produce nearly identical FRC curves, reconstruction errors, and R-factor, with the multislice results showing a very slight improvement, as shown in Supplementary Figs. 19 and 20. Consequently, conventional reconstruction methods, including computationally efficient, non-iterative approaches such as back projection-based methods, can yield reasonably accurate results for such samples, particularly when combined with particle averaging techniques. Nevertheless, the objective of this work is to demonstrate that our physics-based phase retrieval framework can be directly applied to biological datasets and still provide improved resolution and artifact suppression, even from a single tomographic tilt series. Importantly, we anticipate that the benefits of our method will be more pronounced for thicker or more structurally heterogeneous biological specimens, where multiple scattering and interference effects degrade the performance of conventional reconstruction algorithms.

Methods
Sample preparation
Co3O4 nanoparticles were synthesized using the sol-gel method50. In a 100 mL round-bottom flask, a solution of 0.073 g cobalt (II) perchlorate hexahydrate (Aldrich) in 15 mL 1-octanol (Aldrich) was reacted with 2.67 g oleylamine (Acros) and heated to 120 ∘C in air. As the mixture heated to 80 ∘C, 0.7 mL of distilled water was introduced, and the reaction was maintained at 120 ∘C for 6 h with stirring. The reaction was then cooled to room temperature. To purify the nanoparticles, excess acetone was used initially, followed by a second purification where the nanoparticles were dispersed in hexane and precipitated with ethanol. After each purification, the nanoparticles were collected via centrifugation at > 200×g. The supernatant was discarded, and the nanoparticles were redispersed in toluene for drop-casting onto a TEM grid with a carbon substrate.

HRTEM data acquisition
A tomographic tilt and focal series of a Co3O4 nanoparticle supported on carbon was acquired in electron counting mode using a Titan Krios G3i microscope (Thermo Fisher Scientific) operating at an acceleration voltage of 300 kV, equipped with a Gatan K3 direct electron detector operating in correlated-double sampling mode51. The data was collected automatically using SerialEM. The tomography experiment was performed under cryogenic conditions without ice encapsulation, utilizing the super-resolution mode of the K3 detector. The super-resolution pixel size was calibrated to be 0.52 Å. Tilt angles ranged from −65∘ to +65∘ in 1∘ increments, with a total electron dose of 8460 e− Å−2. HRTEM images were recorded as a stack of 24 subframes per image. At each tilt angle, a focal series was acquired at initial defocus values of 100, 250, and 900 nm. Note that the final estimated defocus values differed from the initial settings due to changes in defocus resulting from tilting. The estimated defocus values are plotted in Fig. 1d.

A tomographic tilt series of HIV-1 particles (EMPIAR-10164) was obtained from the open-source EMPIAR database30. We used only the first of the 60 tilt series included in the EMPIAR-10164 dataset for our reconstruction, without using subtomogram averaging. As reported in ref. 30, the tilt range spans from −57∘ to +60∘ in 3∘ increments, with a total electron dose of approximately 100 e− Å−2. HRTEM images were recorded as a stack of 8 subframes per image.

Preprocessing for 3D reconstruction
We applied a series of preprocessing steps to the measured tilt series images, including motion correction, defocus estimation, alignment of focal and tilt series, microscope parameters optimization, and BM3D denoising:

Motion correction. To compensate for sample movement during measurement, we employed the motion correction algorithm from Ref. 43. Each image in our HRTEM tilt dataset consists of 24 subframes, and for each image, the motion-corrected subframes were summed to produce a single corrected 2D image.

Defocus estimation. To estimate the defocus values of the HRTEM images, we used CTFFIND4 software52, which applies Thon ring fitting. Accurate defocus estimation is crucial for the subsequent focal and tilt alignment steps. Therefore, the fitted defocus values from CTFFIND4 are also later refined during the microscope parameter optimization step, described below.

Focal series alignment. During the HRTEM acquisition, undesirable image shifts may be introduced in each image along the x and y-axes. Correcting these shifts is crucial to achieving accurate, high-resolution reconstructions. We first corrected the alignment within each focal series collected at a given tilt angle, where each focal series consists of three images captured at different defocus values.

To align the images within each focal series, we first performed a grid search over possible 2D shift values (Δx and Δy) between the first and second defocus images. The optimal shift minimizes the reconstruction error of a reduced dataset containing only the first defocus image and the shifted version of the second defocus image. We then applied the same approach to align the first and third defocus images in the same manner, thereby completing the alignment correction for the focal series. This alignment procedure was repeated independently for each tilt angle in our tomographic dataset.

To efficiently determine optimal shifts with sub-pixel accuracy, each grid search was conducted in two stages: an initial coarse search covering a 25×25 pixel range with a step size of 2 pixels, followed by a finer search covering a 2×2 pixel range with a step size of 0.5 pixels. To reduce the likelihood of being trapped in local minima, we applied a low-pass filter to the defocus images during each search, with a cutoff of 0.3 Å−1 for the coarse search and 0.7 Å−1 for the fine search.

Tilt series alignment. After correcting the alignment within each focal series, we corrected the alignment between the focal series at different tilts. First, we reconstructed an initial 3D volume using only the focal series at the zero-degree tilt. Next, we performed a grid search over the possible 2D shift values (Δx and Δy) for the focal series at the +1∘ tilt angle. The optimal shift was determined by minimizing the reconstruction error of the +1∘ focal series, using the zero-degree 3D reconstruction as the initial guess.

We then repeated this procedure to obtain the optimal shift for the focal series at the −1∘ tilt angle, this time using a new initial 3D volume reconstructed from both the zero-degree and previously aligned +1∘ tilt series. This alignment process continued sequentially for all remaining tilt angles, sorted by absolute value and alternating between positive and negative tilts (e.g., +2∘, −2∘, +3∘, −3∘, and so forth). For each tilt angle θ, the optimal shift was identified by minimizing the reconstruction error of the corresponding focal series, with the initial 3D volume guess obtained from the reconstruction of all previously tilt-aligned focal series.

Similar to the focal series alignment, we divided each grid search into a coarse search and fine search, with low-pass filter cutoffs of 0.3 Å−1 and 0.7 Å−1, respectively.

Microscope parameters optimization. We optimized the defocus values, C value, and pixel size using Bayesian optimization to minimize the reconstruction error (using scikit-optimize Python package for Bayesian optimization). Because defocus estimation through Thon ring fitting can be imperfect for thick samples, we refined the three defocus values at each tilt by minimizing the reconstruction error of the focal series at that tilt. These defocus values were further refined using the gradient descent method during the final 3D reconstruction to enhance accuracy. Similarly, we applied Bayesian optimization to determine the optimal C value and pixel size, which were determined to be 2.3 mm and 0.52 Å, respectively. The estimated defocus and C values from this step were used as initial parameters for the 3D reconstruction. During phase retrieval-based reconstruction, these values were further refined using a gradient descent method.

Image denoising. Before 3D reconstruction, we applied a denoising procedure to each HRTEM image as follows. The HRTEM experimental images are primarily affected by Poisson noise, which can be described as:
Y=αP(n), (1)
where Y is the intensity at each pixel, α is the gain parameter, and P(n) is the Poisson distribution of n electrons. We estimated the gain parameter α through statistical analysis of the 24 subframes across each image, following previous works5,16,17. The images were denoised by the block-matching and 3D filtering (BM3D) algorithm53. To meet the BM3D algorithm’s requirement for Gaussian noise, we applied the Anscombe transformation54 to convert the Poisson noise into pure Gaussian noise. After denoising with BM3D, we applied the inverse Anscombe transformation.

PhaseT3M, phase retrieval algorithm
By using the HRTEM multislice method as a forward model, we formulate an inverse problem to retrieve the 3D electrostatic potential V from 2D HRTEM images acquired at various tilt angles and focal planes, as illustrated in Fig. 1b, c. This inverse problem is posed as an optimization task that minimizes the error function E2, which quantifies the discrepancy between the experimentally measured HRTEM images and the calculated HRTEM images from the 3D potential V:
V=argmin V E2(θ,Δf)=argmin V ∥Iθ,Δf(r;V)−Iθ,Δf(r)∥2 2, (2)
where ∥⋅∥2 represents the L2 norm, and Iθ,Δf(r;V) and Iθ,Δf(r) are the calculated and measured HRTEM image at each tilt angle θ and each defocus value Δf, respectively. r denotes 2D spatial coordinates (x, y) in real space.

Our reconstruction algorithm for solving the inverse problem has three main components: (1) the forward multislice model to calculate HRTEM images, (2) the inverse multislice model to retrieve a 3D electrostatic potential V, and (3) regularization to account for experimental noise and enhance the quality of the reconstructed volume. These three steps are iteratively applied across all tilt angles until the error converges sufficiently, yielding the final 3D electrostatic potential V. For further detail, see the pseudocode presented in the Supplementary Information.

Forward multislice model. The forward multislice model uses the current estimate of the 3D electrostatic potential V to calculate HRTEM images, denoted as I. To describe the scattering interactions between the electron beam and the sample, we begin by considering the Schrödinger equation for fast electrons within the framework of the forward model23:
∂ψ(x,y,z)/∂z = iλ/4π ∇2 ψ(x,y,z)+iσV(x,y,z)ψ(x,y,z), (3)
where ψ(x, y, z) is the electron wave function, λ is the relativistic electron wavelength, σ is the relativistic interaction constant. To efficiently solve the Schrödinger equation for fast electrons, we use the multislice method, which provides an effective numerical approximation for the beam-sample interaction55 by discretizing the full 3D potential into a stack of 2D projected potentials. This approach facilitates the calculation of simulated HRTEM images by sequentially propagating the electron wave function through each potential slice, accounting for both phase shifts and scattering events.

Specifically, the 3D potential V is discretized into Nz slices along the z-direction, each with a thickness Δz, yielding a 2D projected potential:
V2D(x,y)=∫zm+Δz zm V(x,y,z)dz (4)
where m=1 to Nz represent the slice index.

Next, the evolution of the electron wave function is computed by sequentially applying the free-space propagation function and transmission function for each slice of the material. Specifically, ψm+1, an electron wave function after transmitting through the m-th slice, can be expressed using ψm as
ψm+1,rp(r)=F−1[P(q)F[tm(r)ψm,rp(r)]], (5)
where F and F−1 represent the forward and inverse Fourier transform, respectively.

The free-space propagation function P(q) can be written as P(q)=exp(−iπΔzλq2), with q=√qx2+qy2, qx and qy representing the 2D spatial frequency components in the x and y directions, respectively. The transmission function tm(r) is given as tm(r)=exp(iσV2D,m(x,y)). While the slice thickness Δz can be allowed to vary between slices to accommodate different material properties or experimental conditions, we opted to use a uniform Δz thickness across all slices in our calculation to simplify the computational model.

Finally, the incident parallel electron wave function, ψ0, is initialized to unity in the HRTEM forward model. The exit wave function, ψexit (also denoted as ψNz+1), represents the electron wave after transmission through the sample. To obtain ψexit, we recursively apply Eq. (5) starting from the incident wave function ψ0 and propagate it through each of the Nz slices of the sample, where Nz is the total number of slices. To incorporate defocus effects and other lens aberrations, the final wave function is expressed as:
ψfinal,θ,Δf(r)=F−1[ψexit,θ(q)exp(−iχ(q))], (6)
where χ(q) is the aberration function in Krivanek notation37,56:
χ(q,ϕ)=2π/λ ∑n,n0 (λn+1 qn+1)/(n+1) [Cnn0,a cos(n0ϕ)+Cnn0,b sin(n0ϕ)], (7)
Using the derived gradient formula, the target quantity U, such as V and Cnn0 is iteratively updated as:
U=U−αU ∇U E2(θ), (10)
where αU is the step size controlling the strength of gradient descent. We use Eq. (10) to update the 3D potential V by optimizing image shift, defocus values, and higher-order aberration coefficients across the entire tilt series.

Regularization. Regularization is incorporated to mitigate experimental noise and enhance the quality of the reconstructed volume. Based on the physical assumption that the 3D electrostatic potential is non-negative and real23, we apply a positivity constraint, setting any negative potential values to zero at each iteration. Additionally, to suppress high-frequency noise generated during iterative reconstruction, a 3D low-pass filter is applied to the reconstructed volume at each iteration.

3D phase retrieval-based reconstruction of a Co3O4 nanoparticle
Using the 3D phase retrieval algorithm described above, we reconstructed the 3D potential of a Co3O4 nanoparticle supported on carbon. We used Nz = 600 slices for the multislice calculation, 100 iterations, a positivity constraint, and a 3D low-pass filter cutoff of 0.8 Å−1 to reduce high-frequency noise. Note that in Supplementary Fig. 4b, the only difference is that the positivity constraint was not applied. In Supplementary Fig. 5d–f, a single-slice calculation was used with Nz = 1. The resulting size of the reconstructed 3D volume is 600×600×600 voxels, with a voxel size of 0.52 Å. The error curves for both multislice and single-slice reconstructions converge, as shown in Supplementary Fig. 6. In addition to reconstructing the 3D potential, we further refined the image shifts, defocus values, and high-order aberrations (up to C32) during the process. For the larger field of view reconstruction, as shown in Figs. 1h, 2a, and 4a, the 3D volume was reconstructed with pixel binning of 4, yielding a voxel size of 2.08 Å and a volume size of 225×225×225 voxels.

The iterative phase retrieval problem is inherently non-convex, and convergence to a global minimum cannot be guaranteed in theory. To mitigate the risk of becoming trapped in local minima, we implemented several practical strategies to ensure reliable convergence: (1) validation on simulated datasets with known ground truth, (2) monitoring smooth error reduction over iterations, (3) random shuffling of tilt projections to avoid bias, and (4) using physics-based constraints and the Adam optimizer to enhance stability. Together, these strategies provide strong empirical evidence that our reconstructions converge reliably without being trapped in poor local minima.

Simulation to reproduce the phase retrieval-based tomography experiment
To reproduce the phase retrieval-based tomography experiment, we first created a 3D atomic model of a Co3O4 nanoparticle on a carbon support. The Co3O4 nanoparticle was modeled as a spherical, 7 nm-diameter crystal with a perfect Co3O4 structure. For the carbon support, we generated an amorphous structure by randomly placing atoms within a 60×60×30 nm3 volume, maintaining an atomic density of 100 atoms per nm3. As shown in Fig. 3a, the nanoparticle was placed on the amorphous carbon support, positioned halfway embedded in the carbon structure. Atomic positions of carbon within the overlapping region of the nanoparticle and the support were removed.

Using the generated atomic structure, we simulated a tilt and focal series using the multislice-based HRTEM simulation implemented in the abTEM python package57. We used tilt angles ranging from −65∘ to +65∘ in 1∘ increments, and we applied the same refined defocus values we obtained from our experimental dataset. To mimic experimental conditions, Poisson noise was applied to the images based on our experimental electron dose, and random image shifts were introduced following a Gaussian distribution with a mean of 0 pixels and a standard deviation of 0.5 pixels.

Using this generated tilt and focal series, we reconstructed the simulated 3D potential. The reconstruction parameters matched those used in the experiment, including 600 slices for the multislice calculation, 100 iterations, and a 3D low-pass filter cutoff of 0.8 Å−1. The resulting full volume size was 600×600×600 voxels, with a voxel size of 0.52 Å. A 3D mask was then applied to the reconstructed volume to isolate only the Co3O4 nanoparticle intensity, as shown in Fig. 3b.

As a simulated reconstruction reference, we generated a very sharp 3D potential from the generated atomic model using abTEM57. We then applied Gaussian smoothing to this 3D potential with a standard deviation of 0.65 pixels.

SIRT reconstruction of a Co3O4 nanoparticle
We reconstructed a 3D volume from the tilt and focal series after the preprocessing step using the SIRT algorithm48 with 200 iterations. To suppress high-frequency noise, a 3D low-pass filter with a cutoff frequency of 0.8 Å−1 was applied. The final reconstruction volume size was 600×600×600 voxels, with a voxel size of 0.52 Å. For a larger field of view, as shown in Figs. 2a and 4a and Supplementary Fig. 12a, the 3D volume was reconstructed with pixel binning of 4, resulting in a voxel size of 2.08 Å and a volume size of 225×225×225 voxels.

Spectral signal-to-noise ratio (SSNR) calculation
To calculate the SSNR, we first generated a 2D white noise potential with a size of 100×100 pixels. Using the abTEM Python package57, we simulated HRTEM focal series from the white noise potential with defocus values of 100, 250, and 900 nm. Poisson noise, corresponding to electron dose, was applied to the simulated images, and 50 ensemble focal series were generated. We then reconstructed 50 potentials from these noisy focal series using PhaseT3M. Finally, the SSNR was computed from the reconstructed potentials following the method described in Ref. 58.

HIV-1 reconstruction
A tomographic tilt series of HIV-1 particles (EMPIAR-10164), as described in the HRTEM data acquisition section of “Methods”, was used for reconstruction. Motion correction was performed using the algorithm from ref. 43, applied to the 8 subframes of each image in the tilt series. The motion-corrected subframes were then summed to generate a single corrected image per tilt. Defocus values were estimated using CTFFIND4 (version 4.1.14)52. Tilt alignment was carried out in IMOD (version 5.1)59,60, using a cross-correlation-based coarse alignment followed by fiducial marker-based fine alignment. Before reconstruction, we applied binning by a factor of 9, taking the pixel size from 1.35 Å to 12.15 Å. Various reconstruction methods were independently applied to the binned tilt series data, including phase retrieval, SIRT48, RESIRE61, FBP, WBP62, SART63, and ART64. For the reconstructions, we used the tomopy package (version 1.15.2) for SIRT and ART, the ASTRA toolbox (version 2.3.1) for FBP, WBP, and SART, and the RESIRE MATLAB code61 for the RESIRE reconstruction. The error curves for both multislice and single-slice phase retrieval reconstructions converge, as shown in Supplementary Fig. 20. The final reconstructed volume size was 360×360×360 voxels, with a voxel size of 12.15 Å. We performed 100 iterations for our phase retrieval method and 500 iterations for the SIRT, RESIRE, and SART reconstructions. For the ART reconstruction, 20 iterations were used, yielding the clearest image. Note that all reconstructions of HIV-1 particles were performed without applying a positivity constraint, as variations in vitreous ice thickness can lead to regions of negative potential.

Leave-one-out method
To avoid overfitting and to objectively evaluate the reconstruction quality, we adopted a leave-one-out strategy. In this approach, each projection is systematically excluded from the reconstruction process, and a forward projection is generated from the resulting volume for comparison with the omitted experimental image. To calculate the leave-one-out FRC, we first prepared one measured projection and the corresponding calculated linear projection generated from the 3D volume reconstructed using all tilt projections except for the measured projection, ensuring an unbiased comparison. To minimize Fourier artifacts, the prepared projections were padded using the calculated mean intensity as a constant value to match the background. The FRC was then computed as:
FRC(r)=∑i F1(ri)F2(ri)* / √∑i |F1(ri)|2 ∑i |F2(ri)|2, (11)
where F1 and F2 are the Fourier transforms of the two 2D images, and F2* is the complex conjugate of F2. The index ri refers to an individual pixel at radius r in Fourier space. The FRCs shown in Fig. 5c and Supplementary Fig. 18 were calculated using the zero-degree projection.

R-factor calculation
The R-factor is a metric used to quantify the difference between two images. To evaluate the similarity between the measured projection and the calculated projection from the reconstructed volume, we computed the R-factor as:
R=∑i |Imeas,i−Icalc,i| / ∑i |Imeas,i|, (12)
where Imeas,i and Icalc,i represent the measured and calculated intensities at pixel i, respectively. The index i runs over all pixels in real space. Note that calculated projections Icalc,i were generated from a forward multislice or single-slice model for PhaseT3M, and a linear projection model for all other conventional reconstruction methods. Supplementary Fig. 19 presents the R-factors calculated for all tilt angles across the different reconstruction methods.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The research paper by Lee et al. reports PhaseT∃M, a 3D reconstruction method from images of a tilt series. The novelty of the method is the use of a multislice model to account for multiple scattering and Bayesian optimization for alignment and aberration correction. When applied to radiation-hardy nanoparticles with total electron flounce of about 8000 electrons/Å2 (this is about ~100 times permitted for biological samples), PhaseT∃M attains a maximum resolution of 1.6 Å for a 7 nm Co₃O₄ nanoparticle on a 30 nm carbon substrate, which is very impressive. When compared to classical SIRT reconstruction, PhaseT3M, improves resolution by 20-30% (to 9.0 nm) for HIV-1 particle. As shown, PhaseT3M reduces artifacts and enabling applications across materials science and biology. Despite higher computational demands, PhaseT∃M's compatibility with standard microscopes and open-source datasets highlights its potential for broad scientific impact. There are a few issues that the author should address for convince serious practitioners of cryogenic electron tomography (cryoET) to adapt PhaseT3M, particularly for thick biological samples where multi-slice might offer some advantages.\n\nMajor issues\n1. Figures illustrating reconstructed 3D features are hard to see meaningful features due to the high level of noise in the tomograms. Authors should display three orthogonal density slices, x-y, y-z, z-x, reviews for the tomograms compared in all figures comparing reconstructed volumes.\n2. For the HIV tomogram comparisons, the authors claim that they compared PhaseT3M with a number of methods. Yet, only SIRT result is shown next to the PhaseT3M result in Figure 5. Results from the other reconstruction method should be displayed, particular that from RESIRE, which made similar claim for reconstructing nanoparticles before.\n3. In biological cryoET, SIRT is rarely used nowadays due to poor performance. For example, for reconstruction without subtomogram averaging, IsoNet is most commonly used to obtain directly interpretable tomograms. It would be desirable to see how PhaseT3M compares to IsoNet result, at least for the biological application.\n4. Claim about high resolution for cryoET might be misleading. The problem of biological cryoET is not computational rather, is the inherent limitation of radiation damage electrons have on biological, frozen-hydrated samples. The dose used in the current study is three orders of magnitude higher than typical fluance tolerated by biological samples. In fact, the\n5. While the paper demonstrates PhaseT∃M’s application to HIV-1 particles (Page 6, Lines 294–304), achieving a resolution of 9.0 nm, it notes that this is lower than sub-tomogram averaging results (Page 6, Line 301). The discussion of limitations for biological samples, particularly those with low electron-atom scattering, is brief and could be expanded (Page 7, Lines 250–253).\n6. Include a brief statement in the discussion or methods section about the method’s adaptability to other microscope configurations or any preliminary tests conducted on alternative setups. This would reassure readers of its broad applicability.\n\nMinor issues:\n1. Line 94: “on” to “with” or “using”\n2. Line 336: delete “the” in front of program CTFFIND4", "reviewer_2": "PhaseTEM introduces a novel approach to three-dimensional (3D) phase-contrast imaging by combining high-resolution TEM with nonlinear phase retrieval. Building upon the physical modeling of multiple electron scattering, this method leverages the multislice formalism and incorporates Bayesian optimization to iteratively reconstruct the 3D electrostatic potential of complex samples. In this work, Lee et al. aim to establish PhaseTEM as a generalizable technique for high- fidelity, low-dose 3D imaging, with applications to vitrified biological specimens in cryo-EM. The authors conducted experiments on a 7 nm Co₃O₄ nanoparticle embedded in 30 nm of carbon, without ice encapsulation, using cryo- tomography. They claimed a 3D spatial resolution of 1.6 Å, surpassing conventional limits in standard cryo-EM. However, standard cryo-EM instruments are typically dedicated to the investigation of biological samples that are vitrified by freezing. Even without considering the embedded biological sample itself, the surrounding vitrified ice presents a major challenge for cryo-EM to achieve high resolution due to strict dose limitations and beam-induced motion. Although this work was performed on a cryogenic TEM instrument, the sample condition used to demonstrate the 1.6 Å resolution is effectively that of a standard materials science TEM experiment, and does not realistically mimic the requirements of cryo-EM for biological imaging. It is plausible that the method could be applied to cryo-EM, but the demonstrated resolution cannot be considered a new benchmark for cryo-EM unless the sample is vitrified. Therefore, I feel that the achievement, while technically sound, is not sufficiently significant to the related fields to justify publication in Nature Communications in its current form.\n\nIn addition, several points need to be addressed:\n• Since the manuscript claims that a resolution of 1.6 Å surpasses conventional limits using standard cryo-EM instruments, the authors should provide a summary of current cryo-EM resolution and dose conditions used with appropriate references.\n• The total accumulated dose reported is 8460 e⁻/Å². Can the authors re-examine their experiment and assess the achievable resolution at a dose level realistically tolerated by a cryo sample?\n• How does the method account for defocus variations across the sample thickness, particularly at high tilt angles during phase retrieval?\n• In Figure 3c, only 2D slices from the final reconstruction of the particle are shown. How do these compare directly with the original 3D atomic model? A side-by-side comparison should be provided to help assess the reconstruction accuracy.\n• How did the error metrics evolve with the number of iterations? Can the authors provide a plot showing the convergence behavior?\n• How did the authors determine whether the reconstruction had converged to a global minimum, rather than being trapped in a local minimum?"}, "author_response": "We thank the editor and the reviewers for taking the time and effort to provide very constructive comments, which have been very helpful in making improvements to the manuscript. We hope that the depth and quality of our responses indicate a corresponding dedication to the highest possible quality of a published manuscript resulting from this work.\nPlease find below our response to the reviewers’ comments on our manuscript “PhaseT∃M: 3D Imaging at 1.6 Å Resolution via Electron Cryo-Tomography with Nonlinear Phase Retrieval”. We have carefully addressed all the reviewers' criticisms and suggestions and trust that the manuscript is now acceptable for publication in Nature Communications.\nResponse to reviewers’ comments\nReviewer #1 (Remarks to the Author):\n[...]"}, {"round": 2, "reviewer_comments": {"reviewer_2": "The manuscript provides more information and analysis in this revised version. I am satisfied with most of the responses to my questions.\nHowever, regarding the response to my Question 4, I remain concerned about the validity of the methodology, even for the simulation data. It is highly subjective to state, as in Revision Fig. 5, that “the atomic positions of the simulated Co₃O₄ nanoparticle are consistent with the intensity maxima of the reconstructed volume, although small deviations between the atomic positions and reconstructed intensities are observed.”\nIn the figure (see the circled area as an example, or any other areas), there is no general trend supporting this startment and the statement in the main text (lines 128–131): “The intensity of the reconstructed potential is approximately proportional to the atomic number: lighter elements like carbon and oxygen produce lower potential signal, while heavier elements like cobalt show higher potential magnitude.” In fact, the intensities on Co and O sites appear to fluctuate randomly, and in some locations are even reversed, although the resolution is insufficient to resolve individual atoms.\nSince this work focuses on phase retrieval, the reconstructed intensities should correspond to the correct atomic positions to demonstrate that the reconstruction is free from contrast reversal or multiple-scattering artefacts. As a constructive next step, I suggest first validating the methodology at a substantially higher electron dose, where signal-to-noise limitations are minimal, to establish the best achievable accuracy of the algorithm. The dose can then be progressively reduced to determine the threshold at which reconstruction fidelity begins to degrade. This would provide a clearer separation between algorithmic performance and experimental constraints, and strengthen confidence in the results."}, "author_response": "We sincerely appreciate the reviewer’s positive feedback on the revised manuscript. Also, we fully understand the concern, and we have carefully addressed it in this revision.\nUpon re-examination of the reconstruction, we found a mismatch between atomic position and reconstructed volume, which originated from a slight off-centering during 3D rotation. This has now been corrected, the updated figures (main Fig. 3) in the revised manuscript reflect this improvement.\nIn addition, following the reviewer’s constructive suggestion, we systematically tested the performance and validity of the phase retrieval approach at different electron doses. Specifically, we performed simulations with 8,460 e/Å² (used in experiment), 84,600 e/Å², and infinite dose. To focus the effect of electron dose, we used a model of Co O nanoparticles without carbon background. The simulation includes only missing wedge effects and dose effects. Revision Figure 1 shows these doses tests. As expected, the reconstructions at infinite and high doses are agreement in the ground truth, while progressive reduction in reconstructed resolution and contrast is observed as the dose is reduced. Note that applying a low pass filter of 0.8 Å-1, together with missing wedge effects, hinders full resolution of all the oxygens intensities in the reconstructions. This demonstrates both the validity of our method and its practical limitation under experimental constraints. We have included this discussion (pages 6, lines 220-230) and added the Revision Fig. 1 (Supplementary Fig. 9) in the revised manuscript."}]

Decision letter