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.

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 Parker6 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-imaging9–11 or 3D-VMI13,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 

FIG. 1.

Schematic of a typical VMI tomographic imaging procedure along with a definition of the coordinate system used throughout this publication. Starting from the top and following the arrows clockwise, I(x, y, z) is the 3D distribution to be reconstructed. A hemispherical view of the full distribution is shown here. Different projections, Pα(y, z), are recorded with the distribution rotated by different angles α. The angles α and ϕ (the azimuthal angle used along with θ to define the spherical harmonic functions) both lie in the xy-plane. These projection images can be rearranged into sinograms Pz(y, α) at each z-location of the distribution. Using an inverse Radon Transform (iRT), each sinogram can be reconstructed to give the corresponding slice of the 3D distribution in the xy-plane. Assembling these sequential planar cuts together then yields the fully reconstructed 3D distribution. Each of the different tomographic reconstruction approaches currently used within the VMI community only differ in the iRT step of this cycle.

FIG. 1.

Schematic of a typical VMI tomographic imaging procedure along with a definition of the coordinate system used throughout this publication. Starting from the top and following the arrows clockwise, I(x, y, z) is the 3D distribution to be reconstructed. A hemispherical view of the full distribution is shown here. Different projections, Pα(y, z), are recorded with the distribution rotated by different angles α. The angles α and ϕ (the azimuthal angle used along with θ to define the spherical harmonic functions) both lie in the xy-plane. These projection images can be rearranged into sinograms Pz(y, α) at each z-location of the distribution. Using an inverse Radon Transform (iRT), each sinogram can be reconstructed to give the corresponding slice of the 3D distribution in the xy-plane. Assembling these sequential planar cuts together then yields the fully reconstructed 3D distribution. Each of the different tomographic reconstruction approaches currently used within the VMI community only differ in the iRT step of this cycle.

Close modal

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.

In a recent publication, we made use of the linear algebra-based SIRT35 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,
Ir,θ,ϕ=l=0lmaxm=0lBlmrYlmθ,ϕ.
(1)
Here, the real cosine spherical harmonic functions, Ylmθ,ϕ, are defined as
Yl0θ,ϕ=2l+14π×Pl0cosθ,
(2a)
Ylmθ,ϕ=(1)m2l+12πlm!l+m!×Plmcosθ×cosmϕ,
(2b)
where l and m are zero or positive integers (with lm), Plm(cosθ) 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 Blm(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 lmax) 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 2N, 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 (for the l expansion) and/or cos 2 (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, 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 2N + 1 distinct projections need to be recorded over the interval 0° ≤ α ≤ 180°. This is a consequence of the Nyquist sampling theorem53 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 2N + 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.).

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

FIG. 2.

Illustration of the generalized Fourier–Hankel (FH) transform—a key concept behind FHANTOM. Starting in the top panel, the zeroth-order FH-transform (exactly equivalent to the Abel transform) can be used to recover a cylindrically symmetric distribution (here, the Y20 spherical harmonic) from a single 2D projection. The same is true for any 3D distribution that has a cos() angular dependence; only now an mth-order FH-transform is required. The center and lower panels show examples of the m = 2 and 4 cases of this relationship (specifically the Y42 and Y44 spherical harmonics), but it holds generally for any integer value of m and is independent of l or any other θ-dependence.

FIG. 2.

Illustration of the generalized Fourier–Hankel (FH) transform—a key concept behind FHANTOM. Starting in the top panel, the zeroth-order FH-transform (exactly equivalent to the Abel transform) can be used to recover a cylindrically symmetric distribution (here, the Y20 spherical harmonic) from a single 2D projection. The same is true for any 3D distribution that has a cos() angular dependence; only now an mth-order FH-transform is required. The center and lower panels show examples of the m = 2 and 4 cases of this relationship (specifically the Y42 and Y44 spherical harmonics), but it holds generally for any integer value of m and is independent of l or any other θ-dependence.

Close modal

The HTR algorithm demonstrated in our previous publication49 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 pBASEX23 (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 transform3,19 and BASEX7). 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.

To solve the center-line noise problem encountered with HTR, FHANTOM adopts a similar approach to pBASEX but extends the basis set concept beyond the projection of cylindrically symmetric distributions (i.e., where 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
Ix,y,zIr,θ,ϕ=k=1kmaxl=0lmaxm=0lcklmfklmr,θ,ϕ.
(3)
Here, each basis function fklm(r, θ, ϕ) is given by the equation
fklmr,θ,ϕ=errkσ2×Ylmθ,ϕ,
(4)
where k is an index for each radial Gaussian basis function, rk and σ are the pixel radius and width of the kth basis function, and Ylm(θ, ϕ) 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 gklm(r, θ) may then be defined by integrating Eq. (4) along the entire time-of-flight axis (here, the x axis),
gklmr,θ=fklmr,θ,ϕdx.
(5)
Implementing this approach is, however, very computationally expensive and requires constructing and numerically integrating a 3D distribution for each combination of rk, 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, kmax. 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 mth-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 kmax. 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,
Pm(y,z)Pmr,θ=k=1kmaxl=mlmaxcklmgklmr,θ.
(6)
This can be rewritten using matrix notation as
Pm=CmGm.
(7)
Here, Pm is a vector containing each element of the Pmy,z image, Cm is a vector containing the unknown cklm expansion coefficients, and Gm is a matrix containing each of the gklm(r, θ) functions for the given m value. Inverting this equation leads to an expression for Cm,
Cm=PmGm1,
(8)
where Gm1 is the Moore–Penrose pseudo-inverse of Gm. Even though Gm 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 Pm(y, z) Fourier-profile leads to the complete set of cklm expansion coefficients, allowing for the entire original distribution to be reconstructed using Eq. (3). The angular anisotropy parameters Blm(r) can also be calculated from cklm by expanding the angularly isotropic (i.e., Gaussian) fk00 basis functions,
Blm(r)=k=1kmaxcklmfk00r.
(9)
This removes the need for determining each Blm(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 Blm(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 Blm(r) in any chosen rotated coordinate reference frame using Wigner D-matrices.54 Details on this procedure may be found in the  Appendix.
FIG. 3.

A schematic of the FHANTOM procedure. Starting from the top and following the arrows clockwise, I(x, y, z) is the 3D distribution to be reconstructed (shown here in hemispherical view), with lmax = 4 in this example. Instead of rearranging the Pα(y, z) projection data to generate sinograms, as is performed using more conventional tomography techniques (see Fig. 1), the projections are Fourier transformed along the α-direction to separate out the projection data into the different m contributions. The total number of profiles for any given application may vary, but this will be determined by the experimental geometry or by the number of projections recorded. Using a basis set expansion approach, the cklm coefficients are determined, and the mth-order Fourier–Hankel transform is calculated for each Fourier-profile to give Im(y, z). A cos() dependence can then be applied to sweep the 2D profile into 3D, yielding Im(x, y, z). The final 3D reconstruction is then simply the sum over all the different m values.

FIG. 3.

A schematic of the FHANTOM procedure. Starting from the top and following the arrows clockwise, I(x, y, z) is the 3D distribution to be reconstructed (shown here in hemispherical view), with lmax = 4 in this example. Instead of rearranging the Pα(y, z) projection data to generate sinograms, as is performed using more conventional tomography techniques (see Fig. 1), the projections are Fourier transformed along the α-direction to separate out the projection data into the different m contributions. The total number of profiles for any given application may vary, but this will be determined by the experimental geometry or by the number of projections recorded. Using a basis set expansion approach, the cklm coefficients are determined, and the mth-order Fourier–Hankel transform is calculated for each Fourier-profile to give Im(y, z). A cos() dependence can then be applied to sweep the 2D profile into 3D, yielding Im(x, y, z). The final 3D reconstruction is then simply the sum over all the different m values.

Close modal

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

FIG. 4.

Typical times for the writing and matrix inversion of some representative FHANTOM basis sets. Timings are measured for pixel resolutions of 128 × 128, 256 × 256, and 512 × 512 pixels at different values of lmax. In each case, the spacing between adjacent basis functions Δrk = 2.

FIG. 4.

Typical times for the writing and matrix inversion of some representative FHANTOM basis sets. Timings are measured for pixel resolutions of 128 × 128, 256 × 256, and 512 × 512 pixels at different values of lmax. In each case, the spacing between adjacent basis functions Δrk = 2.

Close modal

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 lmax = 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 Blm 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.

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 lmax = 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.

FIG. 5.

Noisy simulated VMI projection data (left) and corresponding hemispherical 3D reconstructions (center and right) showcasing the performance of FHANTOM and HTR with different signal levels (measured in average counts-per-pixel or cpp). The simulated distributions consist of spherical shells with different anisotropy parameters defined with respect to the laser propagation (horizontal) axis (B20 = −0.4, B42 = 0.2, B44 = 0.1, and B40 = B22 = 0 for the inner shell, and B20 = −0.4, B40 = 0.1, B22 = 0.8, and B42 = B44 = 0 for the outer shell). The axes and labels present in similar figures are omitted here to better show the center line/spot noise. See the text for a full discussion.

FIG. 5.

Noisy simulated VMI projection data (left) and corresponding hemispherical 3D reconstructions (center and right) showcasing the performance of FHANTOM and HTR with different signal levels (measured in average counts-per-pixel or cpp). The simulated distributions consist of spherical shells with different anisotropy parameters defined with respect to the laser propagation (horizontal) axis (B20 = −0.4, B42 = 0.2, B44 = 0.1, and B40 = B22 = 0 for the inner shell, and B20 = −0.4, B40 = 0.1, B22 = 0.8, and B42 = B44 = 0 for the outer shell). The axes and labels present in similar figures are omitted here to better show the center line/spot noise. See the text for a full discussion.

Close modal

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 Δrk = 2, σ = 2, and lmax = 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 publication32).

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.

FIG. 6.

Angle-integrated intensity plots for each of the reconstructions shown in Fig. 5. Across all noise levels (from 0.1 to 100 cpp on average), the intensity distributions appear broadly similar when using HTR and FHANTOM. Even for the extremely sparse images with an average of just 0.1 cpp, clear peaks are still easily identified. Regions 1 and 2 (shaded in green and blue, respectively) are used in Fig. 7 for averaging the Blm parameters. See the main text for a full discussion.

FIG. 6.

Angle-integrated intensity plots for each of the reconstructions shown in Fig. 5. Across all noise levels (from 0.1 to 100 cpp on average), the intensity distributions appear broadly similar when using HTR and FHANTOM. Even for the extremely sparse images with an average of just 0.1 cpp, clear peaks are still easily identified. Regions 1 and 2 (shaded in green and blue, respectively) are used in Fig. 7 for averaging the Blm parameters. See the main text for a full discussion.

Close modal
FIG. 7.

Extracted angular anisotropy parameters for each of the simulated reconstructions shown in Fig. 5. Regions 1 and 2 correspond to the shaded regions of Fig. 6. All parameters are averaged over their respective regions, and the error bars denote the 2σ uncertainty in the mean value. The FHANTOM parameters (shown on the right) are all in excellent agreement with the original simulated values (black bars), with the small exception of the extremely noisy 0.1 cpp reconstructions. For the HTR reconstructions (right), some small discrepancies are seen between the reconstructed and simulated values. In addition, the error bars are larger for the HTR reconstructions compared with FHANTOM. See the main text for a full discussion.

FIG. 7.

Extracted angular anisotropy parameters for each of the simulated reconstructions shown in Fig. 5. Regions 1 and 2 correspond to the shaded regions of Fig. 6. All parameters are averaged over their respective regions, and the error bars denote the 2σ uncertainty in the mean value. The FHANTOM parameters (shown on the right) are all in excellent agreement with the original simulated values (black bars), with the small exception of the extremely noisy 0.1 cpp reconstructions. For the HTR reconstructions (right), some small discrepancies are seen between the reconstructed and simulated values. In addition, the error bars are larger for the HTR reconstructions compared with FHANTOM. See the main text for a full discussion.

Close modal

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.

To validate the performance of FHANTOM on real experimental data, we first revisit a recent photoelectron elliptical dichroism (PEELD) study on the enantiomers of the chiral camphor molecule using a 2 + 1 resonance-enhanced multiphoton ionization (REMPI) scheme at 400 nm.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
Ipr,θ,ϕ=l=02N=6m=0lBlmprYlmθ,ϕ,
(10)
where the superscript p indicates the helicity of the elliptical polarization (with extreme values of +1 for left-circularly polarized, 1 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 Blm 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,
3DPEELDr,θ,ϕ=I+r,θ,ϕIr,θ,ϕ.
(11)
In some previous reports, deviations from cylindrical symmetry about the propagation axis in 3D-PEELD distributions (as well as in examples using mixed linear/circular laser polarizations) have been assumed to be negligibly small. Working under such a supposition, VMI data were then processed in an approximate manner using a standard inverse Abel transform.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.

FIG. 8.

Measurements of the total photoelectron signal following multiphoton ionization of (R)-camphor with linearly polarized 400 nm laser pulses over a range of different projection angles α. The mean intensity (averaged over four experiment cycles) at each projection angle is shown (blue line) along with the mean standard deviation of the intensity (gray shaded area). Since the intensity remains effectively constant for all projection angles, the polarization orientation with respect to the VMI lens system does not noticeably influence the photoionization cross section, and a tomographic reconstruction approach can be used. See the text for additional details.

FIG. 8.

Measurements of the total photoelectron signal following multiphoton ionization of (R)-camphor with linearly polarized 400 nm laser pulses over a range of different projection angles α. The mean intensity (averaged over four experiment cycles) at each projection angle is shown (blue line) along with the mean standard deviation of the intensity (gray shaded area). Since the intensity remains effectively constant for all projection angles, the polarization orientation with respect to the VMI lens system does not noticeably influence the photoionization cross section, and a tomographic reconstruction approach can be used. See the text for additional details.

Close modal

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, S368,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 Δrk = 2 and σ = 2) described in the previous section were used, but now for the angular part, lmax = 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.

FIG. 9.

Reconstructed hemispherical 3D-PEELD distributions following the 2 + 1 ionization of (R)-camphor with 400 nm elliptically polarized laser pulses using the HTR (top) and FHANTOM (bottom) approaches. Just four projections (with α = 0°, 30°, 60°, and 90°) are used to produce these reconstructions. The color scale represents the asymmetry in photoelectron intensity relative to the averaged intensity from the left- and right-elliptically polarized PADs.

FIG. 9.

Reconstructed hemispherical 3D-PEELD distributions following the 2 + 1 ionization of (R)-camphor with 400 nm elliptically polarized laser pulses using the HTR (top) and FHANTOM (bottom) approaches. Just four projections (with α = 0°, 30°, 60°, and 90°) are used to produce these reconstructions. The color scale represents the asymmetry in photoelectron intensity relative to the averaged intensity from the left- and right-elliptically polarized PADs.

Close modal

The spherical harmonic expansion coefficients Blm 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 B32 and B52 expansion coefficients were found to be of comparable magnitude to the m = 0 contributions and varied as a function of S3×1S32. 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 Blm vs S3 obtained using the HTR method. The equivalent outputs returned from FHANTOM [directly from the cklm 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.

FIG. 10.

Selected spherical harmonic expansion coefficients for the 3D-PEELD from the 2 + 1 ionization of (R)-camphor with 400 nm laser pulses. The top panel shows B32 and B52 reconstructed and extracted using HTR and a secondary fitting step. For the lower panel, these expansion coefficients were returned directly from the FHANTOM algorithm. Error bars denote 2σ uncertainty. The full contextualization and interpretation of these results may be found elsewhere.32 

FIG. 10.

Selected spherical harmonic expansion coefficients for the 3D-PEELD from the 2 + 1 ionization of (R)-camphor with 400 nm laser pulses. The top panel shows B32 and B52 reconstructed and extracted using HTR and a secondary fitting step. For the lower panel, these expansion coefficients were returned directly from the FHANTOM algorithm. Error bars denote 2σ uncertainty. The full contextualization and interpretation of these results may be found elsewhere.32 

Close modal

The high-quality reconstructions available through FHANTOM may even be obtained in scenarios where the limits of lmax 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, lmax 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 lmax, above which the distribution will not contain any significant contributions. As a walk-through example, the strong-field (3 × 1013 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.

FIG. 11.

Measurements and analysis of the total photoelectron signal following the strong-field ionization of ethanol at 1030 nm. The top panel shows the total signal as the waveplate is continuously rotated. The lower panel shows the mean intensity averaged over multiple revolutions (blue line) and the mean standard deviation of the intensity (gray shaded area). The modulation in signal intensity seen in the top panel is mostly due to random fluctuations. Only a small oscillation in the mean signal is observed, which is attributed to detector inhomogeneity. As such, the polarization orientation with respect to the VMI lens system does not noticeably influence the photoionization cross section. See the main text for additional details.

FIG. 11.

Measurements and analysis of the total photoelectron signal following the strong-field ionization of ethanol at 1030 nm. The top panel shows the total signal as the waveplate is continuously rotated. The lower panel shows the mean intensity averaged over multiple revolutions (blue line) and the mean standard deviation of the intensity (gray shaded area). The modulation in signal intensity seen in the top panel is mostly due to random fluctuations. Only a small oscillation in the mean signal is observed, which is attributed to detector inhomogeneity. As such, the polarization orientation with respect to the VMI lens system does not noticeably influence the photoionization cross section. See the main text for additional details.

Close modal

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 Δrk = 2 and lmax = 14) was written and inverted in around 3 min on the same MacBook Pro used in Sec. II.

FIG. 12.

Fourier analysis of projection data to determine the appropriate basis set for FHANTOM when reconstructing images recorded from the strong-field ionization of ethanol at 1030 nm. 6000 projection images are recorded by continuously rotating the half-waveplate. In the center panel, a pixel-by-pixel Fourier transform reveals which m components contribute significantly to the projection data. The m = 0, 2, 4, and 6 Fourier-profiles are shown inset as examples. The cumulative spectrum (lower panel) shows what fraction of the projection data is accounted for by considering the first successive m terms. In this example, the projection data are described perfectly with terms up to and including m = 14. See the main text for details.

FIG. 12.

Fourier analysis of projection data to determine the appropriate basis set for FHANTOM when reconstructing images recorded from the strong-field ionization of ethanol at 1030 nm. 6000 projection images are recorded by continuously rotating the half-waveplate. In the center panel, a pixel-by-pixel Fourier transform reveals which m components contribute significantly to the projection data. The m = 0, 2, 4, and 6 Fourier-profiles are shown inset as examples. The cumulative spectrum (lower panel) shows what fraction of the projection data is accounted for by considering the first successive m terms. In this example, the projection data are described perfectly with terms up to and including m = 14. See the main text for details.

Close modal

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.

FIG. 13.

Top row: 3D-PADs recorded from the strong-field ionization of ethanol with 1030 nm laser pulses reconstructed using FHANTOM, HTR, and FBP. Bottom row: The residual between FHANTOM and the HTR and FBP approaches. In both cases, the residuals are incredibly small (note the ×100 multiplication factor for the color scale), indicating the use of a basis set during the reconstruction accurately captures all significant structure in the 3D-PAD.

FIG. 13.

Top row: 3D-PADs recorded from the strong-field ionization of ethanol with 1030 nm laser pulses reconstructed using FHANTOM, HTR, and FBP. Bottom row: The residual between FHANTOM and the HTR and FBP approaches. In both cases, the residuals are incredibly small (note the ×100 multiplication factor for the color scale), indicating the use of a basis set during the reconstruction accurately captures all significant structure in the 3D-PAD.

Close modal

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.

TABLE I.

Some common multiphoton ionization schemes and the corresponding symmetry constraints that the resulting photoproduct angular distributions must obey. These symmetries directly determine how many projection images are required to be able to reconstruct the full 3D distribution and in what range of α these must be sampled. Increasingly more complex polarization geometries are considered before reaching the most general angular distribution, described by Yang’s theorem [Eq. (1)].

Polarization geometry (notes)l valuesm valuesNo. of projections
N × ↑ (one-color multiphoton ionization) 0, 2, …, 2N 
N × ↺ (with chiral molecules) 0, 1, 2, …, 2N 
N × ↺ (elliptical multiphoton ionization with chiral 0, 1, 2, …, 2N 0, 2, …, l N + 1 (0°–90°) 
molecules) 
N×+M× (parallel polarizations e.g., two-color 0, 2, …, 2(N + M
pump–probe measurements) 
N×+(M×) (perpendicular polarizations) 0, 2, …, 2(N + M0, 2, …, l N + M + 1 (0°–90°) 
N×+(M×) (with chiral molecules) 0, 1, 2, …, 2(N + M0, 2, …, l N + M + 1 (0°–90°) 
N×+(M×) (i.e., angle combinations other than  0, 2, …, 2(N + Ml, …, −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 valuesm valuesNo. of projections
N × ↑ (one-color multiphoton ionization) 0, 2, …, 2N 
N × ↺ (with chiral molecules) 0, 1, 2, …, 2N 
N × ↺ (elliptical multiphoton ionization with chiral 0, 1, 2, …, 2N 0, 2, …, l N + 1 (0°–90°) 
molecules) 
N×+M× (parallel polarizations e.g., two-color 0, 2, …, 2(N + M
pump–probe measurements) 
N×+(M×) (perpendicular polarizations) 0, 2, …, 2(N + M0, 2, …, l N + M + 1 (0°–90°) 
N×+(M×) (with chiral molecules) 0, 1, 2, …, 2(N + M0, 2, …, l N + M + 1 (0°–90°) 
N×+(M×) (i.e., angle combinations other than  0, 2, …, 2(N + Ml, …, −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

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

The derivation of the FHANTOM reconstruction approach presented in this article has made the initial simplifying assumption that the original coordinate system of the 3D distribution aligns perfectly with the laboratory measurement frame. When this is true, the distribution will be composed entirely of real 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
Ylmθ,ϕ=(1)m2l+12πlm!l+m!×Plmcosθ×sinmϕ,
(A1)
where 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 Blm expansion coefficients from the imaginary part of the Fourier transformed projections. This is accomplished using the same method as described earlier and returns Blm 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.
For FHANTOM to be able to do this and operate correctly, however, the spherical harmonic basis functions must be defined and generated in the specific laboratory frame of reference introduced in Fig. 1. 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. It is, therefore, desirable to be able to transform the extracted Blm coefficients into different coordinate reference frames. For complex spherical harmonics,
Ylmθ,ϕ=1m2l+14πlm!l+m!Plmcosθeimϕ,
(A2)
this problem is easily solved using the Wigner D-matrices, Dlα,β,γ, 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 Dmml(α,β,γ) as weighting coefficients:
Ylmθ,ϕ=m=llDmmlα,β,γYlmθ,ϕ.
(A3)
Since each real spherical harmonic can be described in terms of the more general complex harmonics,
Ylmθ,ϕ=i2Ylm1mYlmifm<0,Yl0ifm=0,12Ylm+1mYlmifm>0,
(A4)
it is still possible to make use of the Wigner 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,
UB=12000000i000i000i0i0000020000010100010001000000Bl2Bl1Bl0Bl1Bl2=12iBl2iBl2iBl1+iBl12Bl0Bl1Bl1Bl2+Bl2=Bl2Bl1Bl0Bl1Bl2=BC,
(A5)
where 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
BC=Dlα,β,γBC,
(A6)
and, finally, can then be recast once more into real spherical harmonics using the inverse of U,
BR=U1BC.
(A7)
Combining each of these steps together, a rotation operator R(α,β,γ) can be defined as
Rα,β,γ=U1Dlα,β,γU,
(A8)
which directly converts the real expansion coefficients into the new rotated reference,
B=Rα,β,γB.
(A9)
By performing this matrix operation for each 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.
1.
D. W.
Chandler
,
P. L.
Houston
, and
D. H.
Parker
,
J. Chem. Phys.
147
,
013601
(
2017
).
2.
M. N. R.
Ashfold
,
D.
Murdock
, and
T. A. A.
Oliver
,
Annu. Rev. Phys. Chem.
68
,
63
(
2017
).
3.
D. W.
Chandler
and
P. L.
Houston
,
J. Chem. Phys.
87
,
1445
(
1987
).
4.
M. N. R.
Ashfold
,
N. H.
Nahler
,
A. J.
Orr-Ewing
,
O. P. J.
Vieuxmaire
,
R. L.
Toomes
,
T. N.
Kitsopoulos
,
I. A.
Garcia
,
D. A.
Chestakov
,
S.-M.
Wu
, and
D. H.
Parker
,
Phys. Chem. Chem. Phys.
8
,
26
(
2006
).
5.
A. G.
Suits
,
Rev. Sci. Instrum.
89
,
111101
(
2018
).
6.
A. T. J. B.
Eppink
and
D. H.
Parker
,
Rev. Sci. Instrum.
68
,
3477
(
1997
).
7.
V.
Dribinski
,
A.
Ossadtchi
,
V. A.
Mandelshtam
, and
H.
Reisler
,
Rev. Sci. Instrum.
73
,
2634
(
2002
).
8.
C.
Bordas
,
F.
Paulig
,
H.
Helm
, and
D. L.
Huestis
,
Rev. Sci. Instrum.
67
,
2257
(
1996
).
9.
C. R.
Gebhardt
,
T. P.
Rakitzis
,
P. C.
Samartzis
,
V.
Ladopoulos
, and
T. N.
Kitsopoulos
,
Rev. Sci. Instrum.
72
,
3848
(
2001
).
10.
D.
Townsend
,
M. P.
Minitti
, and
A. G.
Suits
,
Rev. Sci. Instrum.
74
,
2530
(
2003
).
11.
J. J.
Lin
,
J.
Zhou
,
W.
Shiu
, and
K.
Liu
,
Rev. Sci. Instrum.
74
,
2495
(
2003
).
12.
M.
Ryazanov
and
H.
Reisler
,
J. Chem. Phys.
138
,
144201
(
2013
).
13.
S. K.
Lee
,
Y. F.
Lin
,
S.
Lingenfelter
,
L.
Fan
,
A. H.
Winney
, and
W.
Li
,
J. Chem. Phys.
141
,
221101
(
2014
).
14.
K.
Mizuse
,
R.
Fujimoto
, and
Y.
Ohshima
,
Rev. Sci. Instrum.
90
,
103107
(
2019
).
15.
C.
Vallance
,
M.
Brouard
,
A.
Lauer
,
C. S.
Slater
,
E.
Halford
,
B.
Winter
,
S. J.
King
,
J. W. L.
Lee
,
D. E.
Pooley
,
I.
Sedgwick
,
R.
Turchetta
,
A.
Nomerotski
,
J.
John
, and
L.
Hill
,
Phys. Chem. Chem. Phys.
16
,
383
(
2013
).
16.
E. S.
Goundreau
,
A. E.
Boguslavskiy
,
D. J.
Moffat
,
V.
Makhija
,
M.
Hemsworth
,
R.
Lausten
,
C.
Marceau
,
I.
Wilkinson
, and
A.
Stolow
,
Rev. Sci. Instrum.
94
,
063002
(
2023
).
17.
Y.
Ranathunga
,
T. A.
Olowolafe
,
S. K.
Lee
, and
W.
Li
,
Rev. Sci. Instrum.
94
,
063303
(
2023
).
18.
Y.
Ranathunga
,
T. A.
Olowolafe
,
E.
Orunesajo
,
H.
Musah
,
S. K.
Lee
, and
W.
Li
,
J. Chem. Phys.
158
,
191104
(
2023
).
19.
B. J.
Whitaker
,
Imaging in Molecular Dynamics
(
Cambridge University Press
,
Cambridge
,
2003
).
20.
E. W.
Hansen
and
P. L.
Law
,
J. Opt. Soc. Am. A
2
,
510
(
1985
).
21.
B.
Dick
,
Phys. Chem. Chem. Phys.
16
,
570
(
2013
).
22.
D. D.
Hickstein
,
S. T.
Gibson
,
R.
Yurchak
,
D. D.
Das
, and
M.
Ryazanov
,
Rev. Sci. Instrum.
90
,
065115
(
2019
).
23.
G. A.
Garcia
,
L.
Nahon
, and
I.
Powis
,
Rev. Sci. Instrum.
75
,
4989
(
2004
).
24.
G. M.
Roberts
,
J. L.
Nixon
,
J.
Lecointre
,
E.
Wrede
, and
J. R. R.
Verlet
,
Rev. Sci. Instrum.
80
,
053104
(
2009
).
25.
G. R.
Harrison
,
J. C.
Vaughan
,
B.
Hidle
, and
G. M.
Laurent
,
J. Chem. Phys.
148
,
194101
(
2018
).
26.
J. O. F.
Thompson
,
C.
Amarasinghe
,
C. D.
Foley
, and
A. G.
Suits
,
J. Chem. Phys.
147
,
013913
(
2017
).
27.
J. O. F.
Thompson
,
C.
Amarasinghe
,
C. D.
Foley
,
N.
Rombes
,
Z.
Gao
,
S. N.
Vogels
,
S. Y. T.
van de Meerakker
, and
A. G.
Suits
,
J. Chem. Phys.
147
,
074201
(
2017
).
28.
B.
Dick
,
Phys. Chem. Chem. Phys.
21
,
19499
(
2019
).
29.
R. A.
Livingstone
,
J. O. F.
Thompson
,
M.
Iljina
,
R. J.
Donaldson
,
B. J.
Sussman
,
M. J.
Paterson
, and
D.
Townsend
,
J. Chem. Phys.
137
,
184304
(
2012
).
30.
F.
Renth
,
J.
Riedel
, and
F.
Temps
,
Rev. Sci. Instrum.
77
,
033103
(
2006
).
31.
C.
Sparling
,
A.
Ruget
,
J.
Leach
, and
D.
Townsend
,
Rev. Sci. Instrum.
93
,
023303
(
2022
).
32.
C.
Sparling
,
A.
Ruget
,
L.
Ireland
,
N.
Kotsina
,
O.
Ghafur
,
J.
Leach
, and
D.
Townsend
,
J. Chem. Phys.
159
,
214301
(
2023
).
33.
A. C.
Kak
and
M.
Slaney
,
Principles of Computerized Tomographic Imaging
(
IEEE Press
,
1999
).
34.
R.
Gordon
,
R.
Bender
, and
G. T.
Herman
,
J. Theor. Biol.
29
,
471
(
1970
).
36.
A. H.
Andersen
and
A. C.
Kak
,
Ultrason. Imaging
6
,
81
(
1984
).
37.
C.
Smeenk
,
L.
Arissian
,
A.
Staudte
,
D. M.
Villeneuve
, and
P. B.
Corkum
,
J. Phys. B: At., Mol. Opt. Phys.
42
,
185402
(
2009
).
38.
P.
Hockett
,
C.
Lux
,
M.
Wollenhaupt
, and
T.
Baumert
,
Phys. Rev. A
92
,
013412
(
2015
).
39.
P.
Hockett
,
M.
Staniforth
, and
K. L.
Reid
,
Mol. Phys.
108
,
1045
(
2010
).
40.
M.
Eklund
,
H.
Hultgren
,
I.
Kiyan
,
H.
Helm
, and
D.
Hanstorp
,
Phys. Rev. A
102
,
023114
(
2020
).
41.
J.-H.
Oelmann
,
T.
Heldt
,
L.
Guth
,
J.
Nauta
,
N.
Lackmann
,
V.
Wössner
,
S.
Kokh
,
T.
Pfeifer
, and
J. R. C.
López-Urrutia
,
Rev. Sci. Instrum.
93
,
123303
(
2022
).
42.
M.
Weger
,
J.
Maurer
,
A.
Ludwig
,
L.
Gallmann
, and
U.
Keller
,
Opt. Express
21
,
21981
(
2013
).
43.
C. A.
Mancuso
,
D. D.
Hickstein
,
P.
Grychtol
,
R.
Knut
,
O.
Kfir
,
X.-M.
Tong
,
F.
Dollar
,
D.
Zusin
,
M.
Gopalakrishnan
,
C.
Gentry
,
E.
Turgut
,
J. L.
Ellis
,
M.-C.
Chen
,
A.
Fleischer
,
O.
Cohen
,
H. C.
Kapteyn
, and
M. M.
Murnane
,
Phys. Rev. A
91
,
031402(R)
(
2015
).
44.
M.
Wollenhaupt
,
M.
Krug
,
J.
Köhler
,
T.
Bayer
,
C.
Sarpe-Tudoran
, and
T.
Baumert
,
Appl. Phys. B
95
,
647
(
2009
).
45.
M.
Wollenhaupt
,
C.
Lux
,
M.
Krug
, and
T.
Baumert
,
ChemPhysChem
14
,
1341
(
2013
).
46.
C.
Lux
,
M.
Wollenhaupt
,
C.
Sarpe
, and
T.
Baumert
,
ChemPhysChem
16
,
115
(
2014
).
47.
A.
Comby
,
E.
Bloch
,
C. M. M.
Bond
,
D.
Descamps
,
J.
Miles
,
S.
Petit
,
S.
Rozen
,
J. B.
Greenwood
,
V.
Blanchet
, and
Y.
Mairesse
,
Nat. Commun.
9
,
5212
(
2018
).
48.
E.
Bloch
,
S.
Larroque
,
S.
Rozen
,
S.
Beaulieu
,
A.
Comby
,
S.
Beauvarlet
,
D.
Descamps
,
B.
Fabre
,
S.
Petit
,
R.
Taïeb
,
A. J.
Uzan
,
V.
Blanchet
,
N.
Dudovich
,
B.
Pons
, and
Y.
Mairesse
,
Phys. Rev. X
11
,
041056
(
2021
).
49.
C.
Sparling
and
D.
Townsend
,
J. Chem. Phys.
157
,
114201
(
2022
).
50.
W. E.
Higgins
and
D. C.
Munson
,
IEEE Trans. Med. Imaging
7
,
59
(
1988
).
53.
M. M.
Woolfson
and
M. S.
Woolfson
,
Mathematics for Physics
(
Oxford University Press
,
Oxford
,
2007
).
54.
E. P.
Wigner
,
Group Theory: And its Application to the Quantum Mechanics of Atomic Spectra
(
Academic Press, Inc.
,
New York
,
1959
).
55.
A.
Comby
,
S.
Beaulieu
,
M.
Boggio-Pasqua
,
D.
Descamps
,
F.
Légaré
,
L.
Nahon
,
S.
Petit
,
B.
Pons
,
B.
Fabre
,
Y.
Mairesse
, and
V.
Blanchet
,
J. Phys. Chem. Lett.
7
,
4514
(
2016
).
56.
S.
Beaulieu
,
A.
Comby
,
D.
Descamps
,
B.
Fabre
,
G. A.
Garcia
,
R.
Géneaux
,
A. G.
Harvey
,
F.
Légaré
,
Z.
Mašín
,
L.
Nahon
,
A. F.
Ordonez
,
S.
Petit
,
B.
Pons
,
Y.
Mairesse
,
O.
Smirnova
, and
V.
Blanchet
,
Nat. Phys.
14
,
484
(
2018
).
57.
V.
Blanchet
,
D.
Descamps
,
S.
Petit
,
Y.
Mairesse
,
B.
Pons
, and
B.
Fabre
,
Phys. Chem. Chem. Phys.
23
,
25612
(
2021
).
58.
V.
Svoboda
,
N. B.
Ram
,
D.
Baykusheva
,
D.
Zindel
,
M. D. J.
Waters
,
B.
Spenger
,
M.
Ochsner
,
H.
Herburger
,
J.
Stohner
, and
H. J.
Wörner
,
Sci. Adv.
8
,
eabq2811
(
2022
).
59.
V.
Wanie
,
E.
Bloch
,
E. P.
Månsson
,
L.
Colaizzi
,
S.
Riabchuk
,
K.
Saraswathula
,
A. F.
Ordonez
,
D.
Ayuso
,
O.
Smirnova
,
A.
Trabattoni
,
V.
Blanchet
,
N. B.
Amor
,
M.-C.
Heitz
,
Y.
Mairesse
,
B.
Pons
, and
F.
Calegari
, arXiv:2301.02002 (
2023
).
60.
A. G.
Suits
and
O. S.
Vasyutinskii
,
Chem. Rev.
108
,
3706
(
2008
).
61.
K. L.
Reid
,
D. J.
Leahy
, and
R. N.
Zare
,
J. Chem. Phys.
95
,
1746
(
1991
).
62.
C.
Sparling
, Ph.D. thesis (
Heriot-Watt University
,
2023
).
63.
D.
Rajak
, Ph.D. thesis (
Université de Bordeaux
,
2023
).
64.
D.
Rajak
,
S.
Beauvarlet
,
O.
Kneller
,
A.
Comby
,
R.
Cireasa
,
D.
Descamps
,
B.
Fabre
,
J. D.
Gorfinkiel
,
J.
Higuet
,
S.
Petit
,
S.
Rozen
,
H.
Ruf
,
N.
Thiré
,
V.
Blanchet
,
N.
Dudovich
,
B.
Pons
, and
Y.
Mairesse
,
Phys. Rev. X
14
,
011015
(
2024
).
65.
C.
Sparling
,
A.
Ruget
,
N.
Kotsina
,
J.
Leach
, and
D.
Townsend
,
ChemPhysChem
22
,
76
(
2021
).
66.
S.
Beauvarlet
,
E.
Bloch
,
D.
Rajak
,
D.
Descamps
,
B.
Fabre
,
S.
Petit
,
B.
Pons
,
Y.
Mairesse
, and
V.
Blanchet
,
Phys. Chem. Chem. Phys.
24
,
6415
(
2022
).
67.
P.
Kalaitzis
,
S.
Danakas
,
C.
Bordas
, and
S.
Cohen
,
Phys. Rev. A
108
,
013106
(
2023
).
68.
G.
Stokes
,
Trans. Cambridge Philos. Soc.
9
,
399
(
1851
).
69.
B.
Schaefer
,
E.
Collett
,
R.
Smyth
,
D.
Barrett
, and
B.
Fraher
,
Am. J. Phys.
75
,
163
(
2007
).
70.
A. P.
Clark
,
M.
Brouard
,
F.
Quadrini
, and
C.
Vallance
,
Phys. Chem. Chem. Phys.
8
,
5591
(
2006
).
Published open access through an agreement with JISC Collections