Fourier–Hankel–Abel Nyquist-limited tomography: A spherical harmonic basis function approach to tomographic velocity-map image reconstruction

Charged-particle imaging techniques are among the most popular tools for studying gas-phase molecular photochemistry and photophysics.^1–5 In particular, the velocity-map imaging (VMI) variant introduced by Eppink and Parker⁶ provides a route to the measurement of high-resolution three-dimensional (3D) velocity and angular distributions of photoproducts following an optical interaction. In the simplest implementation, the entire 3D distribution is projected onto a two-dimensional (2D) detector using a set of annular electrostatic lenses. As discussed elsewhere,^7,8 this projection is robust (i.e., free from significant distortions) provided the kinetic energy gained in the accelerating electric field by the charged-particles is significantly greater than the initial kinetic energy the particles have when created—a situation that is common in such measurements. The detector is typically a pair of microchannel plates (MCP) coupled to a phosphor screen, and the resulting image is then captured with a digital camera. To retrieve the original 3D distribution from the 2D projection, some form of data processing is then required. As an alternative, more sophisticated experiments have also been designed to remove the need for this processing step and instead aim to measure the entire 3D distribution (or a portion of it) directly.^9–18 Although there is clearly great interest in performing these more complex measurements, many research groups still rely on the numerical reconstruction of 3D distributions from 2D projection images. This approach does, however, impose certain restrictions on the distributions (and, consequently, the types of photochemical phenomena) that can be easily studied using 2D imaging alone. In the case of unimolecular photoionization/photodissociation experiments involving linearly polarized light, the distributions of photoproducts produced will exhibit cylindrical symmetry about the laser polarization axis. If this symmetry axis lies parallel to the plane of the imaging detector, then the entire 3D distribution can be uniquely determined from a single 2D projection using the inverse Abel transform.¹⁹ This transform may also be applied when circular laser polarizations are used, with the symmetry axis now being the optical propagation direction. The Abel transform is extensively used for processing data from all sorts of VMI-based experiments, and numerous different software implementations of this transform are in circulation throughout the charged-particle imaging community.^7,8,19–30

With the increasing use of exotic combinations of different laser polarizations in photoexcitation/photoionization schemes, however, methods to accurately recover non-cylindrically symmetric distributions are also highly desirable. We have recently demonstrated how artificial neural networks can be used to predict non-cylindrically symmetric 3D distributions from single 2D images.^31,32 This relies on being able to simulate large volumes of training data, consisting of 3D distributions and their corresponding 2D projections that approximate real experimental scenarios. To build confidence in such novel approaches, however, having alternative benchmarking methods for reliably measuring non-cylindrically symmetric distributions is still important. Fortunately, it is also possible to reconstruct general distributions analytically, although now multiple different projections are required. With multiple images, a tomographic approach based around the Radon transform (of which the Abel transform is a special case) can be performed, and distributions with different or no symmetry constraints can be fully recovered.³³ These reconstruction approaches have some advantages over strategies employing slice-imaging^9–11 or 3D-VMI^13,16–18 spectrometers mentioned earlier, particularly in regard to photoelectron-based measurements. Directly recording 3D photoelectron angular distributions (PADs) requires a more complex multielectrode spectrometer, along with a shot-to-shot analysis with sub-ns timing resolution for the arrival of each electron. This imposes an inherent limitation on the number of multi-hit events per laser shot to maintain a reliable association between particle position and arrival time.^13,16–18 It also potentially hinders future developments of these 3D spectrometers for high-repetition rate laser sources. In contrast, tomographic measurements can be conducted with a relatively straightforward detection system, and no additional equipment or acquisition software is required. Furthermore, by relying on only 2D projection measurements, no multi-hit restrictions are placed on the rate of data acquisition.

A schematic of a typical tomographic imaging experiment, along with the associated coordinate system and overall reconstruction procedure, is shown in Fig. 1. By rotating the dissociating/ionizing laser polarization in angular increments of α with respect to the (detection) yz-plane, various projections of the 3D recoiling photofragment distribution may be individually recorded. These projection data can then be reconstituted into a set of sinograms, each of which corresponds to the Radon transform of a particular slice of the distribution parallel to the xy-plane. By performing an inverse Radon transform on each sinogram, in turn, the full 3D distribution can be recovered. As is the case with the Abel transform, there are numerous ways to perform this specific step. Within the medical imaging community (where tomographic reconstruction was originally developed), analytical solutions based around the Fourier-slice theorem and filtered back-projection (FBP)³³ are commonly used, as are matrix algebra-based methods such as the Algebraic Reconstruction Technique (ART),³⁴ Simultaneous Iterative Reconstruction Technique (SIRT),³⁵ and Simultaneous Algebraic Reconstruction Technique (SART).³⁶ 

Within the photoelectron imaging community, the FBP approach has proved the most popular to date.^37–43 This is likely due to its implementation as a standard function in many software packages such as MATLAB and Python. One of the initial demonstrations of photoelectron tomography by Wollenhaupt et al.⁴⁴ made use of the Fourier-slice theorem approach to reconstruction, and this has also been adopted by others.^45–48 These methods do, however, rely on a large volume of projection data, with α typically varied in small angular steps between 1° and 5° to produce high-resolution sinograms (resulting in around 50-100 projections). The data acquisition time of a tomographic imaging experiment using these approaches is, therefore, around two orders of magnitude longer when compared with typical single image experiments.

Going further, it is still possible to improve the quality of tomographic reconstructions beyond that of HTR while still operating at the Nyquist sampling limit. Only information pertaining to m, the order of the spherical harmonic [and, consequently, only the ϕ–dependent information; see Eq. (2b)], is used to constrain the HTR reconstruction procedure. The information encoded in the degree of the harmonic terms, l (encoded in the θ-dependence), remains unharvested, and yet, this is something that is also typically known and constrained in 3D photoproduct distributions. Instead of creating a single sinogram for each z-position of the distribution and proceeding slice-by-slice (as depicted in Fig. 1), an alternative reconstruction approach that operates over all z-values at once would be able to capture both ϕ- and θ-dependent behavior simultaneously. Building on this idea, a new form of tomographic reconstruction procedure is presented here that we refer to as Fourier–Hankel–Abel Nyquist-limited TOMography, or FHANTOM. This process makes explicit use of the partial-wave description of photofragment distributions to determine which spherical harmonic terms contribute to the projection data. From these extracted partial-wave expansion coefficients, the original distribution may be fully recovered. Limiting the output 3D distribution to only contain smoothly varying spherical harmonic basis functions produces extremely high-quality tomographic reconstructions even with the sole use of the Nyquist-limited minimum number of projections, as in the case of the HTR method. The following section describes the mathematical details of FHANTOM and its implementation. An illustrative walk-through example of the reconstruction procedure is then provided, and using simulated data, the ability of FHANTOM to deal with varying levels of realistic image noise is assessed. Finally, the algorithm is demonstrated using real experimental data recorded from the 2 + 1 multiphoton ionization of (R)-camphor with elliptically polarized 400 nm light, and the strong-field ionization of ethanol with 1030 nm linearly polarized laser pulses. In all cases, FHANTOM is shown to be robust and dependable, with similar performance to other reconstruction approaches on high-quality data and more impressive performance on poor data. An additional, highly compelling advantage of FHANTOM is the automatic retrieval of the spherical harmonic expansion coefficients throughout the reconstruction process, thus eliminating the necessity for additional time-consuming fitting steps (as required for HTR, FBP, etc.).

As with all other tomographic reconstruction procedures, FHANTOM initially takes a series of images P_α(y, z), each of which is recorded at a specified projection angle, α. Unlike other approaches, however, the FHANTOM algorithm does not then generate a set of sinograms from these data. Instead, it evaluates the Fourier transform of this set of projections along the α-direction for each individual pixel to filter the contributions from the different m components. It follows from the HTR approach that the mth Fourier-profile corresponds directly to the projection of the part of the distribution that exhibits a cos mϕ angular dependence.^49,50 Considering first the m = 0 (i.e., constant term) Fourier-profile, this image represents the projection of the cylindrically symmetric component of the full 3D distribution. It is already well established in the imaging community that the central slice through a cylindrically symmetric distribution (here, in the yz-plane; see Fig. 1), may be recovered from its projection via the inverse Abel transform. This mathematical transform may alternatively be performed by taking a Fourier transform, followed by a zeroth-order Hankel transform of the image data. In fact, this is how the inverse Abel transform was calculated in the seminal first demonstration of ion imaging.³ A schematic of this is shown in the top row of Fig. 2. A concept central to the FHANTOM reconstruction approach is the realization that this is just the simplest case of a more general relationship. For any 3D distribution with a ϕ-dependence given by cos mϕ, the central slice can be perfectly recovered from a single projection by performing, in succession, a Fourier and mth-order Hankel transform on the image. This result holds for any θ-dependence in the 3D distribution (i.e., any value or combination of values for l). Having this central slice and knowing how the distribution evolves as a function of ϕ then provides the complete information needed to reconstruct the entire 3D distribution. Example cases when m = 2 and 4 are shown in the middle and lower panels of Fig. 2, respectively. This fact alone, however, is not particularly useful since the photoproduct distributions recorded under symmetry-breaking conditions are typically composed of several different cos mϕ contributions. It only becomes useful once multiple projections have been recorded, such that these different m components may be separated out in an initial Fourier transform filtering step.

The HTR algorithm demonstrated in our previous publication⁴⁹ operates on the same basic principle, but instead of using the Fourier-profile decomposition images explicitly in the analysis, a separate sinogram is created for each z-pixel coordinate, as is performed with all other conventional tomographic reconstruction techniques (see Fig. 1). This approach is fast and simple to implement, but it does come with some disadvantages. As is already well-recognized within the photoion/photoelectron imaging community, numerically calculating the Hankel transform of real experimental data can lead to significant noise being directed along the center axis of the reconstruction (here, the z-axis).^19,22 This may be detrimental for extracting angular anisotropy parameters from reconstructed photoproduct distributions. This issue is typically encountered whenever an inverse Abel transform is calculated in Cartesian coordinates.^{7,19,20,27,29}

For the conventional Abel transform, this center-line noise problem is overcome by operating in polar coordinates instead of using a line-by-line Cartesian approach.^23,24,26 Notably, the pBASEX²³ (polar BAsis Set EXpansion) algorithm makes use of annular basis functions with a radial Gaussian intensity profile and a certain degree of angular anisotropy, described using the m = 0 spherical harmonics. The forward Abel transform of each basis function is numerically calculated before expanding the measured projection data as a linear combination of these projected basis functions. The central slice through the distribution may then be found by expanding the original (i.e., non-Abel transformed) basis functions using the same calculated expansion coefficients. This approach concentrates any noise present in the data toward a single spot at the image center, leading to significantly higher-quality reconstructions. The use of polar basis sets further enhances the reconstruction quality by limiting the output to only consist of smooth functions and fits to the entire projection image at once rather than line-by-line (as is performed in Cartesian methods such as the Fourier–Hankel transform^3,19 and BASEX⁷). This process has the final additional advantage that the angle-integrated radial distribution and any angular distribution parameters are automatically determined analytically from the basis set expansion coefficients. Therefore, they do not need to be calculated in a separate fitting step.

All of the basis set creation outlined earlier, along with the matrix inversion steps, may be performed in advance of any data processing and then be used repeatedly on any data recorded under the same experimental resolution and set of limiting l and m values (which is typically dictated by the photon order and polarization geometry, as expanded upon in Sec. V). Figure 4 shows some typical total basis set creation times (including both the write and inversion steps) for different image resolutions and values of l_max. Timing measurements were carried out on a 2015 MacBook Pro with a 2.5 GHz Intel Core i7 processor and 16 GB of RAM running a MATLAB R2017b implementation of FHANTOM.

For many experimental applications, it is possible to construct a basis set in just a few seconds. Even the highest resolution examples shown in Fig. 4 were completed in less than one hour, significantly faster than the typical write times for pBASEX basis sets.²⁵ For example, the l_max = 4 basis sets would be suitable for reconstructing data in time-resolved photoelectron circular dichroism (TR-PECD) measurements, which utilize a 1 + 1′ linear-pump circular-probe ionization scheme.^55–59 Such measurements have, so far, not made use of tomographic reconstruction approaches and have instead made simplifying cylindrical symmetry assumptions when reconstructing the 3D-PADs. This is likely because of the daunting prospect of recording >50 projections for each pump–probe time delay, as would be required using FBP, for example. With FHANTOM, however, high-quality PAD reconstructions along with all the dynamical information encoded in the time evolution of the B_lm parameters may be reconstructed from just five projections per time step (or, within certain limiting cases, only 3, as will be discussed further shortly). This is a significant improvement over what is possible using standard tomographic reconstruction techniques and will make otherwise prohibitively time-consuming experiments feasible. The same basic approach can also be used, for example, when interrogating vector correlation properties in ion photofragment distributions, which regularly make use of multicolor crossed-linear polarization geometries.⁶⁰ With the basis set written, the actual reconstruction time is fast (typically a few seconds) and comparable to other tomographic algorithms.

Up to this point, the FHANTOM treatment considered so far assumes a “worst-case” scenario for determining which m values contribute to the 3D photoproduct distribution. As was mentioned briefly in our HTR publication,⁴⁹ it is often the case that there is no difference between “up” and “down” in imaging experiments, and consequently, the top half of a VMI image will mirror the bottom half. This is the result of a plane of reflection symmetry in the original 3D distribution. Enforcing this top-bottom symmetry means that only even m components are present in the partial-wave expansion, which further reduces the data required for an accurate tomographic reconstruction. With this assumption built in, now only projections in the range 0° ≤ α ≤ 90° need to be considered. These are then duplicated to fill in the remaining 90° ≤ α ≤ 180° range required for the reconstruction, effectively halving the total number of projections (and, consequently, the acquisition time) needed. This symmetry restriction has already been utilized in other photoelectron tomography studies,^37,38 but these do not make use of the Nyquist-limited aspect available only through the HTR and FHANTOM approaches. This idea will be developed further in Sec. V, where examples of typical multiphoton ionization schemes and their corresponding symmetry and reconstruction requirements will be discussed.

To showcase FHANTOM, its performance is evaluated and compared to HTR on some noisy datasets. In our previous publication,⁴⁹ we provided extensive comparisons between HTR and other alternative reconstruction approaches and concluded that HTR always performs as well as or better than FBP and SIRT on the same data. Therefore, these alternatives will not be discussed further in this section. The first column of Fig. 5 shows simulated projection data (only the α = 0° images are shown as a representative example) with varying levels of experimental noise. This was quantified by the average number of signal counts-per-pixel (cpp). The noisy projections were simulated by first generating a smooth and noise-free 3D distribution on a 256 × 256 × 256 grid, with l_max = 4. This form of distribution is representative of photofragment data that may be recorded in a 1 + 1′ crossed-linear polarization geometry, for example.⁶¹ Projections were then obtained by numerically integrating a copy of the distribution rotated through five different angles—the Nyquist-limiting number of projections—specifically, α = −90°, −45°, 0°, 45°, and 90°. Note that the Nyquist-limit does not dictate that these five projections be evenly spaced in α, and in principle, any five unique projections contain the required information for tomographic reconstruction. This is not discussed further, however, since this requires the use of a non-uniform fast-Fourier transform step, which is not currently implemented in FHANTOM.

Random samples were then generated from these noise-free projections using the Poissonian imnoise filter in MATLAB, simulating the stochastic arrival of individual ion/electron strikes. By varying the cpp in each simulation—0.1, 1, 10, and 100 cpp are considered here—the effects of realistic experimental noise can be evaluated for both FHANTOM and HTR. The full sets of simulated projection data were then fed into the HTR and FHANTOM algorithms to compare their respective outputs. In this example, the basis used in the FHANTOM reconstruction consisted of functions with Δr_k = 2, σ = 2, and l_max = 4.

The top row of Fig. 5 shows data simulated and reconstructed at an average count rate of 100 cpp. As can be seen from the example projection image (left), this level of signal corresponds to extremely high-quality data. From a visual inspection of the output reconstructions (center and right columns of Fig. 5), in this regime, FHANTOM and HTR perform equally well. Only a small amount of center line (HTR) and center spot (FHANTOM) noise is present in the reconstructed data. As the signal level decreases to 10 cpp (corresponding to a more typical high-quality VMI data acquisition regime), both reconstruction methods still perform very well, although the center noise becomes more pronounced in both cases. For FHANTOM, however, this remains very small and localized in the center of the distribution, whereas for HTR, the noise is spread across the (z axis) center line. At 1 cpp, the center line noise in the HTR reconstruction is now of similar overall intensity to the true signal level of the 3D data—something not seen with FHANTOM until extremely sparse average count rates (0.1 cpp) are reached. Using a Cartesian line-by-line reconstruction approach (as in HTR) means that neighboring slices parallel to the xy-plane are reconstructed independently, and no extra information on the shape of the 3D distribution in perpendicular planes (i.e., the yz-plane) can be used to improve the overall output quality. This factor becomes increasingly prominent as signal levels decrease and the Poissonian nature of the image noise becomes significant. The FHANTOM approach instead effectively couples together every single voxel in the 3D distribution to solve for the reconstruction that best matches with all the sampled projection data. The performance of both approaches can be compared further by examining the integrated intensity distributions and angular parameters produced in each case. It should be noted again that these values are returned automatically by FHANTOM [using Eq. (9)] during the reconstruction procedure, while a secondary lengthy fitting step is required to retrieve these same values from the HTR reconstructions (details of which may be found in our earlier publication³²).

Figure 6 shows the angle-integrated intensity distributions extracted from the noisy projection data using FHANTOM (top panel) and HTR (bottom panel). For both approaches at all noise levels, the intensity plots agree very well, although a small amount of oscillatory noise is present for count rates below 1 cpp. For such measurements, there is, therefore, no significant advantage to the FHANTOM approach over HTR besides the time saved in determining the velocity distribution in parallel with the reconstruction process. Similar observations are made when comparing Cartesian and polar Abel inversion strategies; the integral nature of the velocity distribution is extremely robust to even high levels of noise.⁶² Instead, it is in the angular distributions where the major benefits of FHANTOM can be found. Figure 7 shows bar plots of the spherical harmonic coefficients extracted using both reconstruction methods, averaged over each of the two peaks in the velocity distribution (Region 1 and Region 2 in Fig. 6), and normalized by the average intensity. The original simulated anisotropy parameter values are also presented. FHANTOM (shown in the two left-hand panels) performs extremely well here, accurately extracting the coefficients for all noise levels. Even in the extremely low signal regime of 0.1 cpp, the calculated values remain in excellent agreement with the original simulated input, although now with a larger associated error. This is, however, a very extreme case, and when operating at more realistic signal levels, FHANTOM is capable of accurately retrieving the angular coefficients with a high degree of precision. HTR (right-hand panels) still manages to reliably extract angular parameters in good agreement with the simulated values, but not with the same high accuracy or precision possible with FHANTOM. The effects of uniform random background noise to emulate detector dark counts have very little effect on either method and are, therefore, not discussed in detail here. This is because they are effectively filtered away during the initial Fourier transform step,^63,64 which is a further benefit of adopting a Fourier-based tomographic reconstruction approach.

It may even be feasible to further improve the reconstruction quality possible with FHANTOM at still lower signal levels by taking directly into account the Poisson statistics of the image noise. Similar approaches have been applied to situations where the inverse Abel transform is used for image analysis of cylindrically symmetric distributions.^21,28,30 Instead of finding the linear-least squares solution to the problem (as is performed in most reconstruction approaches), the maximum-likelihood solution can instead be found, which allows for statistical noise models to be incorporated directly into the reconstruction process. Alternatively, we have shown previously how artificial neural networks may be used to effectively enhance sparse-count VMI data and produce high-resolution “denoised” images.⁶⁵ Both of these approaches, however, are beyond the scope of this publication and are not considered further.

Before re-examining these data using FHANTOM, it is first important to confirm that the modulation of the polarization direction with respect to the static extracting electric field of the VMI spectrometer does not have any influence on the total ionization cross section. As has recently been highlighted by Kalaitzis et al. for the case of ionization very close to threshold,⁶⁷ if the static field does exert an influence, the critical assumption that the sole role of the VMI instrument is to project the 3D-PAD onto the detector will be violated. This can be confirmed by measuring the total yield of photoelectron intensity at each projection angle. These data are presented in Fig. 8 for the case of ionization using linear polarization. As can be clearly seen, there is minimal variation in the total photoelectron yield at each projection angle, and this is well within the error bounds of the measurement. Similar results are found across all laser ellipticities. Therefore, tomographic reconstruction approaches may be reliably applied under these experimental conditions.

In our previous publication, the 3D-PEELDs were reconstructed using the HTR technique with seven equally spaced angular projections between α = 0° and 180°.³² This implicitly assumes the “worst case” scenario for three-photon ionization and allows (in principle) for contributions from odd m harmonics in the 3D-PEELD. These terms, however, should not be present in the data due to the up-down symmetry of the light–matter interaction, and only even m components will be non-zero. Therefore, just four projections on the interval α = 0°–90° contain sufficient information to perfectly reconstruct the corresponding 3D-PAD/PEELD, and that scenario will now be considered here.

The top row of Fig. 9 shows examples of 3D-PEELD distributions recorded at different laser ellipticities (quantified by the third Stokes parameter, S₃^68,69) reconstructed using HTR with just four equally spaced projections between α = 0° and 90°. The bottom row shows 3D-PEELD distributions reconstructed from the same projection data, but these have instead been processed using FHANTOM. The same radial basis functions (with Δr_k = 2 and σ = 2) described in the previous section were used, but now for the angular part, l_max = 6 because of the overall 3-photon REMPI process. Purely by visual inspection, the outputs reconstructed using FHANTOM are smoother than the HTR results and exhibit less noise in their structure. In the $S3=0.50$ case (and to a lesser extent, $S3=0.77$⁠), it is also possible to observe the break in cylindrical symmetry about the z-axis in the FHANTOM reconstructions. This is much harder to visualize in the HTR examples, where background artifacts make discerning subtle structural changes in the 3D-PEELD by eye far more challenging.

The spherical harmonic expansion coefficients B_lm extracted from each 3D-PEELD can be plotted as a function of ellipticity. As discussed previously,³² the most significant result of our previous work considering these data was the presence of previously unreported spherical harmonic terms with m ≠ 0 in the 3D-PEELD distribution. In particular, the B₃₂ and B₅₂ expansion coefficients were found to be of comparable magnitude to the m = 0 contributions and varied as a function of $S3×1−S32$⁠. This trend is attributed to alignment introduced in the initial two-photon resonant excitation step of the 2 + 1 REMPI process.

The top row of Fig. 10 shows the plot of B_lm vs S₃ obtained using the HTR method. The equivalent outputs returned from FHANTOM [directly from the c_klm coefficients using Eq. (9)] are shown in the lower row. Both reconstruction approaches produce results that are in excellent quantitative agreement with each other. The error bars in Fig. 9 are, however, slightly smaller for the FHANTOM method. This is likely to be a consequence of the reconstruction noise being more centrally localized, corrupting the extraction of the angular anisotropy parameters to a lesser extent. It should also be noted that camphor is a popular system for study with gas-phase imaging techniques because it yields such a large and structured asymmetry in the 3D-PAD. We anticipate that for other more challenging systems (i.e., those with smaller and more subtle overall asymmetries, or for molecules exhibiting low vapor-pressures or ionization cross-sections), the advantages of using FHANTOM will be far more evident.

The high-quality reconstructions available through FHANTOM may even be obtained in scenarios where the limits of l_max are not necessarily dictated by theory (unlike the case of multiphoton ionization). In such instances, it is not immediately clear if a small set of tomographic projections can be effectively extended to describe the entire original distribution. In the case of strong-field ionization, l_max is unbounded in the description of the laboratory frame PAD. Consequently, an infinite number of m terms (and, consequently, projections) are formally needed to be able to perfectly reconstruct the distribution. Of course, this is not practically feasible. Instead, a Fourier analysis of the projection image data can be used to empirically estimate a suitable value for l_max, above which the distribution will not contain any significant contributions. As a walk-through example, the strong-field (3 × 10¹³ W cm⁻²) ionization of ethanol with linearly polarized light at 1030 nm at 1 MHz repetition rate is considered. The instrument and procedure used for this form of measurement have been described in detail previously.^63,64 Briefly, the alignment of the linear polarization with respect to an imaging detector was varied continuously using a motorized half-waveplate moving at 45° s⁻¹, and the corresponding 2D projection of the 3D-PAD was recorded on a CMOS camera with a 50 ms exposure time. Through this process, the 3D-PAD was rotated a total of 40 times, sampling 150 projection images per revolution. This acquisition method, along with the high-repetition rate laser, allows for many high-quality tomographic projections to be recorded in rapid succession. In previous publications using this setup, the images were Fourier transformed pixel-by-pixel (along the waveplate angle direction) to filter out any background noise in the data,^63,64 as discussed briefly earlier. Here, we show how the same Fourier-filtered data can be used in conjunction with FHANTOM to produce high-quality tomographic reconstructions.

First, using the procedure introduced in the previous section, the total photoelectron signal recorded at each waveplate position is plotted in the top panel of Fig. 11. Although the signal appears to oscillate over the waveplate rotation, only a very small modulation in the photoelectron yield can be seen when the intensity is averaged for a single rotation (lower panel of Fig. 11). We attribute this to mild detection inhomogeneities present in the MCP used to record these data, particularly toward the center of the detector. Since the 3D-PAD is elongated along the ionizing polarization direction, the α = 0° projection images are mostly unaffected by this since the photoelectron intensity is dispersed across the detector face. For the α = ±90° projections, however, most of the photoelectrons in the 3D-PAD arrive over a small area toward the center of the MCP, where detection efficiency is lowest. This explains the shape of the measured total photoelectron signal as a function of α. More generally, this can occur when recording images containing pixels with high peak intensities, and so we stress that caution should be taken to avoid any artifacts caused by saturating the imaging camera when using any form of tomography for VMI applications.

With this verification step complete, performing a pixel-by-pixel Fourier transform on the imaging data is then able to separate out the different m components of the projected 3D-PAD. This procedure is shown schematically in Fig. 12. A power spectrum indicating the relative strengths of the different m contributions can be produced by plotting the normalized root-mean-square of the pixel intensities of each Fourier-profile. This reveals clear and prominent peaks at even values of m (as anticipated from the up-down symmetry of the distribution), with only a negligible level of background signal present otherwise. In fact, the meaningless non-integer frequency values can be used to establish a threshold below which integer m terms may be ignored. This processed spectrum is shown in the center row of Fig. 12, with the m = 0, 2, 4, and 6 Fourier-profiles shown inset next to their corresponding amplitude peak. Enough projection data are available in this instance to measure the contribution of m terms up to and including m = 74. From the processed spectrum, however, no terms above m = 14 are non-zero. A cumulative spectrum can be calculated by summing together the individual contributions from each successive m, revealing what fraction of the distribution is encompassed by the first m terms. This is shown in the final panel of Fig. 12. From this plot, for example, terms up to and including m = 12 are found to make up 99.4% of the projection data. This would be a sensible threshold, but since 100% of the distribution is described by including m = 14, this value was used for generating a FHANTOM basis set for this specific reconstruction problem. This basis set (with a resolution of 256 × 256 pixels, σ = 2 Gaussian basis functions, with Δr_k = 2 and l_max = 14) was written and inverted in around 3 min on the same MacBook Pro used in Sec. II.

The basis set expansion step of FHANTOM can be performed directly on each of the Fourier profiles below the m = 14 threshold to calculate their contribution to the strong-field 3D-PAD. Alternatively, the entire volume of projection data can be processed using FBP or HTR, and comparisons can be made between the different approaches. Figure 13 shows the reconstructions produced via each of these distinct methods, along with the residual difference between FHANTOM and the HTR or FBP reconstructions. It can be clearly seen that the distributions appear almost identical to each other, and the residuals are extremely small, with a maximum absolute value of ∼1%. This indicates that the basis set approach central to FHANTOM does not compromise the reconstruction of 3D-PADs, even in cases where the partial-wave expansions are not a priori truncated at small values of l and m.

Given that the ionization potential of ethanol is around 10.5 eV, a minimum of nine photons at 1030 nm would be required to drive a multiphoton ionization process. Therefore, terms up to and including m = 2 × 9 = 18 might be expected to contribute to the overall 3D-PAD. Further analysis of the reconstructed distributions also reveals the presence of photoelectrons with a maximum kinetic energy of around 25 eV, which would correspond to the absorption of at least 30 photons within a multiphoton picture, resulting in an even higher value of m = 60. The absence of any angular terms beyond m = 14 in the reconstructed data leads to the conclusion that the angular dependency of the photoelectron lies within the expectation of a multiphoton ionization framework, whereas the kinetic energy distribution clearly reveals the strong-field ionization regime.

To encourage the adoption of FHANTOM within the charged-particle imaging community, a specific guide for users employing the approach for the first time is presented in this concluding section. Examples of some common polarization geometries that are regularly used in N or N + M multiphoton-ionization schemes are listed in Table I, and for each case, the limiting values of l and m, along with the number and range of angular projections required for a robust reconstruction, are given.

The multiphoton-PEELD example considered in Sec. IV A made use of the case shown in the third row of the table for elliptical polarizations. It is also always possible to empirically determine the limits of l and m (and, hence, the minimum number of required projections) using a Fourier analysis of the recorded projection data, as used in Sec. IV B for the case of strong-field ionization. As an aside, we also note that similar sampling requirements must also apply to slice-imaging experiments.^9–11 For full 3D reconstruction using a slice-imaging approach, multiple slices must be recorded at different points along the projection direction (by carefully gating the detector at different times) or instead by rotating the distribution (as is performed in tomographic measurements) and slicing through different angles. To the authors’ knowledge, the latter approach has not yet been adopted. In either case, the Nyquist-limit sampling argument that is applicable to the projection data (discussed throughout this article) will also be applicable to the sliced image data. Therefore, measurements of ϕ-dependent photoproduct angular distributions require the same minimum volume of imaging data, regardless of whether slicing is used or not. With this in mind, we stress that tomography (which, at present, is typically only used in photoelectron imaging experiments) may well be a simple and low-cost alternative to slice-imaging and is particularly attractive for reconstructing ϕ-dependent photofragment angular distributions, such as those encountered when studying photofragment angular momentum polarization.^60,70

The partial-wave expansion model routinely used throughout many gas-phase imaging experiments allows for the simple analysis of 3D photoproduct distributions. The general angular form and symmetry properties of the spherical harmonic functions used in the expansion make them uniquely suitable for tomographic reconstruction using extremely angle-sparse datasets. In this report, we have shown how the FHANTOM reconstruction process can produce 3D photoproduct distributions from projection data by using the partial-wave model to simplify and make intuitive sense of the reconstruction process. We have also demonstrated how FHANTOM produces the best quality reconstructions from the minimum volume of projection image data and, furthermore, also returns the important B_lm coefficients without resorting to extra fitting steps. We anticipate that any experiment aiming to investigate photoproduct distributions lacking cylindrical symmetry within the VMI community will be able to exploit the clear benefit of FHANTOM. Moreover, we feel that advanced tomographic image reconstruction approaches are currently an under-appreciated alternative to more complicated 3D-VMI schemes, both for imaging ions and electrons.

This work was supported by Leverhulme Trust Research Project Grant No. RPG-2012-735, Carnegie Trust Research Incentive Grant No. 70264, EPSRC Platform Grant No. EP/P001459, Agence Nationale de la Recherche (ANR)—Shotime (Grant No. ANR-21-CE30-038-01), and the European Research Council (ERC) under the European Union’s Horizon 2022 research and innovation Program No. 682978—EXCITERS. Heriot-Watt University (HWU) is acknowledged for providing Chris Sparling with Ph.D. funding. Finally, we acknowledge Lewis Ireland (HWU) for fruitful discussions.

The authors have no conflicts to disclose.

Chris Sparling: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Writing – original draft (lead); Writing – review & editing (equal). Debobrata Rajak: Investigation (supporting); Methodology (supporting); Writing – review & editing (supporting). Valérie Blanchet: Funding acquisition (supporting); Investigation (supporting); Methodology (supporting); Supervision (supporting); Writing – review & editing (supporting). Yann Mairesse: Conceptualization (supporting); Formal analysis (supporting);Funding acquisition (supporting); Investigation (supporting);Methodology (supporting); Supervision (supporting); Writing– review & editing (supporting). Dave Townsend: Conceptualization(supporting); Funding acquisition (lead); Methodology (supporting);Project administration (lead); Supervision (lead); Writing –original draft (supporting); Writing – review & editing (equal).

To assist in the adoption of FHANTOM within the VMI community, a MATLAB implementation is available from the corresponding authors upon reasonable request.