4.1 Sample selection and initial characterization
We sampled a ~ 2 kg block of Bishop Tuff welded ignimbrite from unit Ig1Eb according to the stratigraphic definitions of Wilson and Hildreth (1997). Oxygen isotopic data from this unit suggests a lack of post-emplacement hydrothermal alteration and probable preservation of igneous Fe-oxides (Holt and Taylor 1998). After breaking the block into smaller volumes using non-magnetic BeCu tools, we conducted further crushing using alumina grinding vessels in a Spex Shatterbox® at the M.I.T. Radiogenic Isotope Laboratory. Subsequent manual picking yielded approximately 40 zircons. We previously analyzed these zircons for their paleomagnetic record using the SQUID Microscope, finding paleointensities with < 15% uncertainties that are concordant with bulk sample results (Gee et al. 2010; Fu et al. 2017). Two zircons (RF1 and RF2) that were not part of the paleointensity heating experiments were selected for ptycho-tomography analysis (Figs. S1a, d). These crystals were mounted in a block of non-magnetic EPO-TEK 301 epoxy and given a 0.4 T isothermal remanent magnetization in the out-of-plane direction. We imaged these zircons using a QDM at Harvard University, which took data in projective magnetic microscopy (PMM) mode using a 0.9 mT bias field canceled during measurement to < 1 µT (Glenn et al. 2017). The effective spatial resolution of the QDM maps (Figs. S1b, e) was 2.4 µm. The isothermal remanent magnetic moment of each crystal was obtained using an upward continuation/dipole fitting approach, yielding 1.0 x10-10 and 7.6 x 10–11 Am2 for RF1 and RF2, respectively. A fit focussing on just the strong magnetic signal observed in the upper left of RF1 (Figs. S1b) yielded 8.0 x 10–11 Am2 with an estimated 20% uncertainty.
Following the QDM mapping, the two crystals were carefully extracted from the epoxy and mounted on sharp tungsten carbide tips using ultraviolet (UV) setting glue (Figs. S1c, f). Secondary electron (SE) imaging (Figs. S2a, c) and energy-dispersive X-ray (EDX) spectroscopy (Figs. S2b, d) of part of the surface of each zircon crystal was performed using a Quanta-650F scanning electron microscope (SEM) at the Department of Earth Sciences, University of Cambridge. Imaging was performed in low-vacuum mode to prevent charging of the uncoated samples at an accelerating voltage of 20 kV (RF1) and 30 kV (RF2). Chemical maps were calculated using the open source programme HyperSpy by integrating X-ray counts for the Zr Lα, Si Kα, Ca Kα and Fe Kα peaks. Integration was performed without background subtraction over an integration window equal to twice the full-width at half maximum of the peak. Chemical maps for Fe, Ca and Zr were normalised to their maximum counts and used as the red, green and blue channels, respectively, of the RGB images shown in Figs. S2b, d.
4.2 X-ray ptycho-tomography and X-ray fluorescence measurements
All experiments were performed using the I13-1 coherence branchline of the I13 long beamline for imaging and coherence at the Diamond Synchrotron Source, Didcot, United Kingdom. I13-1 is ideally suited to conduct multi-scale and multi-modal ptychography and ptycho-tomography (Rau 2017; Kuppili et al. 2017; Sala et al. 2018; Batey et al. 2018; Cipiccia et al. 2019; Cipiccia et al. 2021; Batey et al. 2022). We employed a hybrid, dual-mode experimental geometry with an X-ray energy of 15 keV. A focusing optics ensemble, comprising a 400 µm diameter blazed Fresnel zone plate, a 60 µm diameter central beam stop and a 25 µm diameter order-sorting aperture was employed to obtain the optimum X-ray illumination function (Kuppili 2020). The Merlin detector was positioned 1.29 m from the sample to record the far-field diffraction patterns in high-resolution mode (HR-mode). This geometry allowed us to obtain 5 nm pixel size with 384x384 detector cropping. The Excalibur detector was positioned at 14.6m to conduct ptycho-tomography in the large FOV survey mode (S-mode). The larger area of the Excalibur detector enabled imaging at multiple resolutions depending on the choice of detector cropping. Pixel sizes of 85 nm and 21 nm were achieved with 256x256 and 1024x1024 detector cropping, respectively. We employed optimum overlap parameters and imaged the entirety of each crystal in each case. To satisfy the over-sampling criterion we employed a probe of ~ 1 µm extent in HR-mode and ~ 10 µm extent in S-mode. In each case we acquired the optimum number of ptychographic angular projections. In S-mode we calculated this value depending on the moderate resolution corresponding to 256x256 cropping. The Merlin detector was positioned on top of a long-range linear stage, and could be moved in and out of beam, thus giving us the ability to shift between S- and HR-modes seamlessly. We employed a Vortex X-ray Fluorescence (XRF) detector positioned perpendicular to the beam propagation to acquire chemical information.
4.3 X-ray ptychography processing and reconstruction
Ptychographic reconstructions were carried out using ‘Ptycho’, a ptychographic reconstruction code based on the Difference Map algorithm coded in the Python programming language (Thibault et al. 2009; Enders and Thibault 2016). The code was parallelised using mpi4py (Dalcín et a. 2007) and can be accelerated by employing multiple computer cores for a given single ptychographic projection. As a single ptycho-tomography dataset may contain 600–1000 ptychographic projections, reconstruction is inherently computationally intensive. A sample ptychographic projection was put through multiple rounds of the algorithmic image reconstruction process to optimise the reconstruction parameters (e.g., number of iterations, Fourier relaxation factor, and X-ray illumination modes) (Kuppili 2020). The optimized parameters were then used for all the ptychographic projections of a given ptycho-tomography dataset. As this process becomes near impractical with conventional computational resources, we employed the ARCHER supercomputing cluster to bring down the processing times from months/weeks to days/hours. Even here, we faced challenges in reconstructing the S-mode datasets with 1024x1024 cropping owing to the large RAM requirement. In this regard, we were successful in exploiting the information from QDM mapping to identify sub-regions of interest within the larger field-of-view. Only the most important regions were reconstructed at the highest obtainable resolution (21 nm) in S-mode. This step reduced the computational challenge and made the reconstructions feasible in practical time frames.
Owing to the ultra-large field-of-view, combined with the dense nature of samples, the ptychographic reconstruction process was particularly challenging in this study. We developed two innovative constraint mechanisms to mitigate these issues. First, the high density of the sample can be mitigated by using higher beam energy (15 keV in this case). Imaging samples at high energies (> 12 keV) opens up exciting possibilities, such as the ability to image thick, dense 3D samples. However, at higher energies there is a possibility that the amplitude and phase-shift maps diverge in their appearance (Supplementary Text). This divergence introduces artifacts primarily in the amplitude component of the object’s transmission function during the ptychographic reconstruction process, ultimately leading to sub-optimal ptychographic reconstructions that are prone to artifacts and a sub-optimal reconstruction of the illumination function. We introduced a Gaussian blurring step that blurs the amplitude of the object’s transmission function before every iteration step. This new constraint improved the ptychographic reconstructions, whereby artifacts corrupting the reconstructions were removed (Supplementary Text). Second, while the aforementioned constraint improved ptychographic reconstructions, the reconstructed illumination function still appeared sub-optimal in some of the reconstructions. The reconstructed illumination function acts as a counterweight in assuring the quality of a reconstruction (Kuppili 2020) and one needs to diagnose/mitigate the factors resulting in sub-optimal illumination functions. We developed a novel diagnosis technique based on Principal Component Analysis (PCA) that is able to analyze and correct for the root causes of sub-optimal probe function (Supplementary Text). In ptychography the far-field diffraction pattern is formed by interaction between object’s transmission function and a partially coherent probe function. A partially coherent probe can be thought of as the resultant of 𝑁 mutually incoherent, coherent illumination function modes with different amounts of power associated with each mode (Thibault and Menzel, 2013, Kuppili 2020). Usually, an initial illumination guess is made, whereby the primary probe mode (depending on experimental geometry) possesses the majority of probe power, while the remaining modes possess random noise with nominal power. Such an initial guess has too many degrees of freedom to absorb unwanted structure into the illumination function. Over multiple iterations, this aspect corrupts the probe ultimately corrupting the object’s transmission function. We have therefore devised a new constraining methodology in which one starts with a guess where all modes are identical containing identical amount of power. This prevents the probe from getting corrupted, thus improving the ptychographic reconstructions (Supplementary Text).
Tomographic reconstructions were performed using the phase component of the ptychographic projections, which has far greater contrast than the absorption component in the hard X-ray regime. Prior to reconstruction, the phase components underwent multiple pre-reconstruction steps. The datasets were corrected for phase ramp using a combination of automated and manual phase-ramp correction techniques (Kuppili 2020). The projections were aligned in the axial direction, employing automated and manual edge-detection methods (Kuppili 2020). The projections were corrected in the radial direction employing a center-of-mass alignment method (Kuppili 2020). Finally, 3D reconstructions were obtained by employing a modified ‘filtered back projection’ method on the derivatives of the aligned ptychographic phase projections, which circumvents the need to unwrap the phase images prior to reconstruction (Guizar-Sicarios et al. 2011; Pfeiffer et al. 2007).
4.4 Phantom calculation
To explore the limits of detectability for magnetite particles using X-ray ptychography, and how the size of the reconstructed particles may be affected by the phase-retrieval and reconstruction processing, a ‘phantom’ dataset containing a mixture of magnetite, apatite and void space embedded in zircon was created and processed in a manner that replicates as closely as possible the treatment of the experimental data. The phantom for simulations was a modified 3D Shepp-Logan phantom (Kak and Slaney 1988) containing 1024x1024x1024 pixels with 85 nm pixel size (similar to the S-mode data with 256x256 detector cropping). The total extent of the phantom is 87 µm, similar to the sample dimensions imaged in this study. The original phantom was modified by adding ellipsoids of various sizes with random orientations, embedded within, or lying on the edges of, the existing primary ellipsoids. Each ellipsoid is assigned a “refractive index” value based on the chemical phase these ellipsoids are expected to represent.
Ptychography datasets at a given rotation angle (φ) were generated by finding the 2D complex transmission function from the rotated 3D refractive index phantom (Kuppili 2020). The 2D transmission function was multiplied by a complex illumination function (borrowed from one of the experimental datasets) with optimum overlap. As each scan point, the intensity of far-field diffraction (assuming a sample to detector distance of 14.6 m) pattern was saved (256x256 cropping) along with the position information of the scan point. This process was then repeated for 1600 equally spaced angles in the − 90° to 90° angular range. It should be noted that the generated datasets were quite ideal in theory and something similar will be very hard to achieve in a real experiment.
The ptychographic phase projection at a given angle was reconstructed using the same algorithmic framework as the experimental datasets. Same process of optimization was employed and all datasets were reconstructed on the ARCHER supercomputer. We employed the same constraint mechanisms in reconstructing simulated datasets and also tested the efficacy of our newly developed constrained mechanisms in terms of quality/resolution of the obtained reconstructions.
The reconstructed phase projections were taken through tomographic reconstruction routine same as the experimental datasets. Automated phase ramp correction was employed, a manual edge detection was used to correct for drifts/shifts in the axial direction and correction in radial direction was done using center of mass alignment method. 3D reconstructions were obtained by employing a modified “filtered back projection” method on the derivatives of the aligned ptychographic phase projections same as the experimental datasets.
4.5 Simulation studies on factors effecting ptychographic resolution
In order to fully understand the limits of resolution especially in the HR mode and evolve a criterion for optimum flux required to resolve features we have carried out extensive simulations. The first step involved engineering illumination functions to accurately depict the experimental conditions. We have achieved this by borrowing the probe functions from the ptychographic reconstruction showed in Fig. 10c. The amplitude of the borrowed probe was scaled by a given factor Pfsc such that one gets increasing or decreasing photon values for various Pfsc. These engineered probes were then used in conjunction with the reconstructed object functions to generate the forward modeled ptychographic data sets. These data sets were then reconstructed using the same algorithmic parameters as the reconstruction shown in Fig. 10c. The reconstructed objects were then compared to obtain parity between the forward model involving engineered illumination functions and the synthetically scaled experimental data sets (Fig. S10). Once parity was established, we used the same engineered illumination functions to generate ptychographic data sets with resolution phantoms. The resolution phantom contains a series of (6) magnetite spheres embedded midway in a 100 um thick zircon matrix. The diameters are 5000, 1000, 500, 300, 200 and 100 nm arranged in a spiral fashion (Fig. S5). The scaling of the illumination functions varied between 0.1 to 100 in order to fully gauge the effects of photons on the obtainable resolution. In order to further push the understanding on resolution limits we have carried similar studies on a finer phantom with 7 magnetite spheres 5000, 1000, 500, 300, 200, 100 and 50 nm. The scaling of the illumination function was limited to 100 (Fig. S11).
4.6 Segmentation and visualization
Tomographic reconstructions were analysed using the software Dragonfly Pro by ORS Inc. A non-local means (NLM) and/or median filter was applied to the data to reduce noise prior to segmentation. Segmentation was performed using intensity thresholding to highlight individual phases and microstructural features of interest, followed by manual clean up. Smoothing of the segmented surfaces was applied to individual objects to remove outlying voxels and improve the quality of the visualizations, taking care not to visibly compromise their morphological features. Visualizations presented throughout the paper show the zircon host (grey), apatite (green), magnetite (red) and pores (black). All features discussed in detail are identified as magnetite due to their association with Fe signals observed in XRF scans and magnetic signals observed in the QDM. However, in the absence of direct diffraction evidence it is not possible to be definitive regarding this phase identification, so we refer to all particles as putative magnetite.