A new method for the fully generalized reconstruction of three-dimensional (3D) photoproduct distributions from velocity-map imaging (VMI) projection data is presented. This approach, dubbed Fourier–Hankel–Abel Nyquist-limited TOMography (FHANTOM), builds on recent previous work in tomographic image reconstruction [C. Sparling and D. Townsend, J. Chem. Phys. **157**, 114201 (2022)] and takes advantage of the fact that the distributions produced in typical VMI experiments can be simply described as a sum over a small number of spherical harmonic functions. Knowing the solution is constrained in this way dramatically simplifies the reconstruction process and leads to a considerable reduction in the number of projections required for robust tomographic analysis. Our new method significantly extends basis set expansion approaches previously developed for the reconstruction of photoproduct distributions possessing an axis of cylindrical symmetry. FHANTOM, however, can be applied generally to any distribution—cylindrically symmetric or otherwise—that can be suitably described by an expansion in spherical harmonics. Using both simulated and real experimental data, this new approach is tested and benchmarked against other tomographic reconstruction strategies. In particular, the reconstruction of photoelectron angular distributions recorded in a strong-field ionization regime—marked by their extensive expansion in terms of spherical harmonics—serves as a key test of the FHANTOM methodology. With the increasing use of exotic optical polarization geometries in photoionization experiments, it is anticipated that FHANTOM and related reconstruction techniques will provide an easily accessible and relatively low-cost alternative to more advanced 3D-VMI spectrometers.

## I. INTRODUCTION

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^{6} 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.^{19} 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.^{33} 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)^{33} are commonly used, as are matrix algebra-based methods such as the Algebraic Reconstruction Technique (ART),^{34} Simultaneous Iterative Reconstruction Technique (SIRT),^{35} and Simultaneous Algebraic Reconstruction Technique (SART).^{36}

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.*^{44} 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.

^{35}algorithm and demonstrated, for the first time, how it can be applied to photoelectron imaging measurements and outperform FBP when presented with the same input data.

^{49}This publication also focused on a second (and much less well-known) tomographic reconstruction method based on the Hankel transform. The Hankel Transform Reconstruction (HTR)

^{50}approach drastically simplifies the tomographic retrieval of 3D objects that are best described within a

*spherical*coordinate system. Photofragment angular distributions fall into this category and can usually be characterized by an expansion in spherical harmonic functions,

^{51}

^{,}

*l*and

*m*are zero or positive integers (with

*l*≥

*m*), $Plm(cos\theta )$ are the associated Legendre polynomials, and

*θ*and

*ϕ*are the polar and azimuthal spherical angles in the

*yz*- and

*xy*-planes, respectively (see Fig. 1). The

*B*

_{lm}(

*r*) coefficients (known as

*anisotropy parameters*) encode the relative weightings associated with each spherical harmonic in the expansion for each velocity/energy. A form of this equation was first derived by Yang for describing angular distributions in nuclear reactions,

^{52}and in the context of photoionization/dissociation or molecular scattering experiments, it is commonly known as a

*partial-wave expansion*. For many experimental scenarios, the upper bounds on

*l*(denoted

*l*

_{max}) and

*m*(determined by

*l*and any additional symmetry constraints, as

*m*= 0 for cylindrical symmetry) are typically small, and so the expansion can be truncated to just a few terms. As shown in our previous uses of the HTR approach for photoelectron imaging applications,

^{32,49}these upper bounds are the

*only*critical factor in determining how many distinct angular projections are required to uniquely recover the 3D distribution using tomography. In photoionization experiments, for example, the maximum value that

*l*can take in the lab-frame description of a PAD is 2

*N*, where

*N*is the total number of photons absorbed in the ionization process. Truncating the partial-wave expansion in this way means that the highest frequency oscillations in the PAD will be described by terms in cos 2

*Nθ*(for the

*l*expansion) and/or cos 2

*Nϕ*(for the

*m*expansion). As discussed in the original HTR publication of Higgins and Munson,

^{50}this highest frequency dependence in the source distribution is mapped directly into the projection data, leading to a maximum oscillation of cos 2

*Nα*, where

*α*is the angle between the principal polarization axis and the imaging plane (here, the

*yz*-plane). To accurately measure the highest frequency features in the projection data and, therefore, reconstruct the highest frequency features in the full 3D distribution, just 2

*N*+ 1 distinct projections need to be recorded over the interval 0° ≤

*α*≤ 180°. This is a consequence of the Nyquist sampling theorem

^{53}and drastically reduces the data burden that has been previously associated with, for example, FBP and Fourier-slice reconstruction. HTR has so far remained a relatively unexplored approach for the more typical medical applications of tomography, simply because this angular sampling rule does not in any way simplify the reconstruction of medical images. Human organs cannot be effectively described with a small number of sine/cosine functions, and so many projection images are still required, and techniques such as FBP remain the most suitable. When the 2

*N*+ 1 rule is appropriate, however, HTR allows for more ambitious experiments to be carried out that would otherwise be prohibitively time consuming to perform with more standard tomographic techniques, as demonstrated, for example, in a recent publication investigating photoelectron elliptical dichroism effects in the enantiomers of camphor.

^{32}

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.).

## II. FOURIER–HANKEL–ABEL NYQUIST-LIMITED TOMOGRAPHY (FHANTOM)

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 *m*th 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.^{3} 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 *m*th-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^{49} 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^{23} (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^{7}). 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.

*m*is no longer constrained to be zero). Using the established partial-wave expansion introduced in Eq. (1), the 3D distribution is assumed to take the form

*f*

_{klm}(

*r*,

*θ*,

*ϕ*) is given by the equation

*k*is an index for each radial Gaussian basis function,

*r*

_{k}and

*σ*are the pixel radius and width of the

*k*th basis function, and

*Y*

_{lm}(

*θ*,

*ϕ*) is a real cosine spherical harmonic function described by Eq. (2). Using exclusively cosine-like spherical harmonics provides a simpler derivation of the FHANTOM algorithm, but additional

*sine-like*spherical harmonic functions (required for describing more general 3D distributions that are rotated with respect to the lab reference frame) can easily be integrated into the methodology. This is discussed further in the Appendix. Examples of basis functions and associated projections are shown in Fig. 2. The case of

*m*= 0 is, by definition,

*ϕ*-independent, and so these elements reduce to the conventional (and cylindrically symmetric) Legendre polynomial basis functions used in the likes of pBASEX (to within a normalization factor). The projected set of basis functions

*g*

_{klm}(

*r*,

*θ*) may then be defined by integrating Eq. (4) along the entire time-of-flight axis (here, the

*x*axis),

*r*

_{k},

*l*and

*m*. It also results in a very large basis set, the size of which (and the associated write time) scales with the

*cube*of the number of radial basis functions,

*k*

_{max}. This can, therefore, place heavy demands on local computer memory. Far more efficient basis generation can be realized by invoking the

*forward*version of the Fourier–Hankel transform. In this way, the

*yz*-planar slice through the basis function distributions can be easily constructed numerically, and then the full projection can be simulated by using a forward

*m*th-order Fourier–Hankel transform (see Fig. 2). Now, the basis set consists only of 2D images, so the size is proportional to the

*square*of

*k*

_{max}. With the basis constructed, all that remains is to expand each Fourier-profile of the projection data $Pmy,z$ in terms of the relevant associated Legendre polynomial basis functions,

**P**

_{m}is a vector containing each element of the $Pmy,z$ image,

**C**

_{m}is a vector containing the unknown

*c*

_{klm}expansion coefficients, and

**G**

_{m}is a matrix containing each of the

*g*

_{klm}(

*r*,

*θ*) functions for the given

*m*value. Inverting this equation leads to an expression for

**C**

_{m},

**G**

_{m}. Even though

**G**

_{m}is potentially a very large matrix, its inverse can be calculated efficiently using, for example, the pinv function in MATLAB. Successively applying Eq. (8) for each different

*P*

_{m}(

*y*,

*z*) Fourier-profile leads to the complete set of

*c*

_{klm}expansion coefficients, allowing for the entire original distribution to be reconstructed using Eq. (3). The angular anisotropy parameters

*B*

_{lm}(

*r*) can also be calculated from

*c*

_{klm}by expanding the angularly isotropic (i.e., Gaussian)

*f*

_{k00}basis functions,

*B*

_{lm}(

*r*) in a separate fitting step after the 3D distribution has been recovered, as is required when using other tomographic approaches.

^{38,39}This makes FHANTOM a uniquely streamlined tomographic reconstruction technique for VMI applications. A schematic of the reconstruction process is shown in Fig. 3. The

*B*

_{lm}(

*r*) coefficients extracted using FHANTOM will always correspond to spherical harmonic functions defined in the coordinate system of Fig. 1 (i.e., with the angle

*θ*lying in the

*yz*-plane and

*ϕ*in the

*xy*-plane). In some cases, however, it may also make physical sense to define them with respect to an alternative coordinate system. For example, PADs produced in experiments involving linearly polarized light are typically described and analyzed in a coordinate system defined with respect to the laser polarization vector rather than the optical propagation direction. Fortunately, it is relatively straightforward to redefine

*B*

_{lm}(

*r*) in any chosen rotated coordinate reference frame using Wigner

*D*-matrices.

^{54}Details on this procedure may be found in the Appendix.

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.^{25} 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.^{60} 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,^{49} 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.

## III. RESILIENCE TO NOISE

To showcase FHANTOM, its performance is evaluated and compared to HTR on some noisy datasets. In our previous publication,^{49} 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.^{61} 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^{32}).

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.^{62} 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.^{65} Both of these approaches, however, are beyond the scope of this publication and are not considered further.

## IV. EXPERIMENTAL DEMONSTRATIONS

### A. Multiphoton photoelectron elliptical dichroism of (*R*)-camphor

^{32}This process generates 3D-PADs that do not, in general, exhibit cylindrical symmetry about any axis. The overall form for this type of 3D-PAD is given by

*p*indicates the helicity of the elliptical polarization (with extreme values of $+1$ for left-circularly polarized, $\u22121$ for right-circularly polarized, $0$ being linearly polarized, and elliptical polarizations taking on any intermediate values). PEELD presents itself as a forward-backward asymmetry in the PAD (with respect to the laser propagation direction, or

*z*-axis) due to contributions from odd-degree (i.e., odd values of

*l*) spherical harmonics. The sign of this asymmetry (and that of the corresponding odd-degree

*B*

_{lm}coefficients) then flips upon reversing the handedness of the elliptical polarization (i.e., the sign of

*p*). 3D-PEELD distributions (containing only the odd

*l*components) were obtained by ionizing (

*R*)-camphor with both left- and right-elliptically polarized light, tomographically reconstructing the resulting PADs, and subtracting one from the other to remove the contributions from even

*l*terms (which do not change with the sign of

*p*),

^{47}

^{,}

^{55–59,66}Through an in-depth tomographic analysis using the HTR method, however, it was shown that symmetry-breaking (i.e.,

*m*≠ 0 spherical harmonic) terms can, in fact, make a significant contribution to the 3D-PEELD distributions, being of similar magnitude to the cylindrically symmetric (i.e.,

*m*= 0) terms. An in-depth discussion regarding this may be found in the original publication.

^{32}

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,^{67} 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°.^{32} 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*_{3}^{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,^{32} 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*_{32} and *B*_{52} expansion coefficients were found to be of comparable magnitude to the *m* = 0 contributions and varied as a function of $S3\xd71\u2212S32$. 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*_{3} 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.

### B. Strong-field ionization of ethanol

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^{13} W cm^{−2}) 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^{−1}, 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.

## V. A USER GUIDE FOR COMMON POLARIZATION GEOMETRIES

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.

Polarization geometry (notes) . | l values
. | m values
. | No. of projections . |
---|---|---|---|

N × ↑ (one-color multiphoton ionization) | 0, 2, …, 2N | 0 | 1 |

N × ↺ (with chiral molecules) | 0, 1, 2, …, 2N | 0 | 1 |

N × ↺ (elliptical multiphoton ionization with chiral | 0, 1, 2, …, 2N | 0, 2, …, l | N + 1 (0°–90°) |

molecules) | |||

$N\xd7\u2191+M\xd7\u2191$ (parallel polarizations e.g., two-color | 0, 2, …, 2(N + M) | 0 | 1 |

pump–probe measurements) | |||

$N\xd7\u2191+(M\xd7\u2192)$ (perpendicular polarizations) | 0, 2, …, 2(N + M) | 0, 2, …, l | N + M + 1 (0°–90°) |

$N\xd7\u2191+(M\xd7\u21ba)$ (with chiral molecules) | 0, 1, 2, …, 2(N + M) | 0, 2, …, l | N + M + 1 (0°–90°) |

$N\xd7\u2191+(M\xd7\u2197)$ (i.e., angle combinations other than | 0, 2, …, 2(N + M) | −l, …, −2, 0, 2, …, l | N + M + 1 (0°–90°) |

parallel or perpendicular) | |||

“Worst case” N-photon process (i.e., Yang’s theorem) | 0, 1, 2, …, 2N | −l < m < l | 2N + 1 (0°–180°) |

Polarization geometry (notes) . | l values
. | m values
. | No. of projections . |
---|---|---|---|

N × ↑ (one-color multiphoton ionization) | 0, 2, …, 2N | 0 | 1 |

N × ↺ (with chiral molecules) | 0, 1, 2, …, 2N | 0 | 1 |

N × ↺ (elliptical multiphoton ionization with chiral | 0, 1, 2, …, 2N | 0, 2, …, l | N + 1 (0°–90°) |

molecules) | |||

$N\xd7\u2191+M\xd7\u2191$ (parallel polarizations e.g., two-color | 0, 2, …, 2(N + M) | 0 | 1 |

pump–probe measurements) | |||

$N\xd7\u2191+(M\xd7\u2192)$ (perpendicular polarizations) | 0, 2, …, 2(N + M) | 0, 2, …, l | N + M + 1 (0°–90°) |

$N\xd7\u2191+(M\xd7\u21ba)$ (with chiral molecules) | 0, 1, 2, …, 2(N + M) | 0, 2, …, l | N + M + 1 (0°–90°) |

$N\xd7\u2191+(M\xd7\u2197)$ (i.e., angle combinations other than | 0, 2, …, 2(N + M) | −l, …, −2, 0, 2, …, l | N + M + 1 (0°–90°) |

parallel or perpendicular) | |||

“Worst case” N-photon process (i.e., Yang’s theorem) | 0, 1, 2, …, 2N | −l < m < l | 2N + 1 (0°–180°) |

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}

## VI. SUMMARY AND CONCLUSION

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX: ROTATIONS AND CHANGING COORDINATE REFERENCE FRAMES

*cosine*-like spherical harmonic functions [see Eqs. (1) and (2)], and the Fourier transform of the projection data will also be entirely real. If instead the distribution is rotated with respect to the lab frame—either by some systematic error or through a real physical effect—additional

*sine*-like (but still

*real*) spherical harmonics will contribute to the expansion

*m*is now a

*negative*integer, by convention. This will result in a

*complex*Fourier transform of the projection data. Accounting for this rotation is straightforward and is achieved in FHANTOM by calculating the sine-like

*B*

_{lm}expansion coefficients from the

*imaginary*part of the Fourier transformed projections. This is accomplished using the same method as described earlier and returns

*B*

_{lm}for both positive and negative

*m*values that are defined in the laboratory frame. Since the real spherical harmonic functions form a complete set of orthonormal functions, these values are sufficient to, in principle, reconstruct

*any*real 3D distribution.

*B*

_{lm}coefficients into different coordinate reference frames. For

*complex*spherical harmonics,

*Wigner D-matrices*, $Dl\alpha ,\beta ,\gamma $, where

*α*,

*β*, and

*γ*are the Euler angles that transform between the original reference frame (described with the spherical coordinate angles

*θ*and

*ϕ*) and the new rotated frame (denoted

*θ*′ and

*ϕ*′).

^{54}For each spherical harmonic, the rotation can be described as an expansion in other spherical harmonics of the same degree, with the Wigner

*D*-matrix elements $Dm\u2032ml(\alpha ,\beta ,\gamma )$ as weighting coefficients:

*D*-matrices for the rotation of real expansions, provided they are first re-expressed using complex functions. Another (square) matrix,

**U**, may be constructed for each

*l*to do this,

**B**is a vector containing each of the

*real*expansion coefficients for a certain value of

*l*, and $BC$ is the vector of corresponding

*complex*coefficients. These are then rotated into the new reference frame (denoted with primes) using the Wigner

*D*-matrix

**U**,

**B**with different

*l*values, it is simple to transform between coordinate systems, adapting the output of FHANTOM to suit any relevant experimental context without any extensive reanalysis.

## REFERENCES

*Imaging in Molecular Dynamics*

*Principles of Computerized Tomographic Imaging*

*Mathematics for Physics*

*Group Theory: And its Application to the Quantum Mechanics of Atomic Spectra*