Three-dimensional momentum imaging of dissociation in flight of metastable molecules

We investigate dissociation in flight of metastable molecular dications formed by ultrashort, intense laser pulses using the cold target recoil ion momentum spectroscopy technique. A method for retrieving the lifetime(s) of the transient metastable state(s) as well as the complete three-dimensional momenta of the dissociating fragments is presented. Specifically, we demonstrate and discuss this approach by focusing on dissociation in flight of the ethylene dication going to the deprotonation channel. Two lifetimes are found to be associated with this process, C2H 4 2 + → C2H3+ + H+: τ 1 = 202 ± 10 ns and τ 2 = 916 ± 40 ns. For the corresponding channel in deuterated ethylene, lifetimes of τ 1 = 269 ± 29 ns and τ 2 = 956 ± 83 ns are obtained.


Introduction
Many processes in molecules happen on fast timescales. For example, rotations and vibrations typically occur on picosecond and femtosecond timescales, respectively. Molecular bond rearrangement and fragmentation often proceed on similarly short timescales, as shown in [1][2][3][4][5] and many others. Hence, femtosecond laser pulses, possessing temporal durations shorter than these timescales, can be used to influence and shed light on molecular dynamics [6][7][8][9].
Not all processes in molecules, however, proceed so swiftly. Multiply charged molecular ions can exist in metastable states that lead to fragmentation happening on long timescales that range from picoseconds to even seconds [10][11][12]. The lifetimes of these transient systems are governed by the relevant potential energy landscape and the mechanisms responsible for decay, which can include tunneling, predissociation, and radiative decay to repulsive states. Investigating the formation, properties, and decay of these metastable molecular ions experimentally and theoretically has been a prominent field of research (see review papers [10][11][12] and [13][14][15][16][17][18][19][20][21][22][23][24], for example).
We study the decay dynamics of metastable molecules by employing coincidence three-dimensional (3D) momentum imaging, which provides the complete 3D momenta of the fragments and therefore their kinetic energy release (KER) and angular distributions. This information can in turn facilitate understanding of the dissociation mechanism(s), demonstrated for example in [17,25,26]. Hence, the 3D momentum imaging technique has been a powerful tool in studies of molecular fragmentation following ionization by ultrashort intense laser pulses, x-ray (or extreme ultraviolet) photons, or fast ion impact [27,28], as long as the breakup is prompt. Prompt breakup happens on a sub-picosecond timescale, much shorter than the flight times of the fragments to the detector. The ionization processes in such experiments can readily form multiply charged metastable molecular ions, seen for example in [17,22,23,[29][30][31][32][33][34].
Importantly, in coincidence measurements, a metastable molecular ion may survive beyond the interaction region, traveling through the spectrometer for a non-negligible time before undergoing dissociation in flight. In this unimolecular fragmentation process, which is a subset of delayed dissociation, the survival time of a fragmenting metastable molecule, t d , is a significant fraction of the time of flight (TOF) of intact metastable ions, t m . In the present experiments, t d is on the order of hundreds of nanoseconds to a few microseconds. In general, the observable range of t d may differ depending on the specific metastable system, as well as the conditions of the experiment.
Field and Eland developed a method to extract the lifetimes of metastable states decaying in flight by fitting Monte-Carlo simulated time-of-flight-difference distributions N (t t 2 1 -) to the corresponding measured timedifference spectrum. Here, t 1 and t 2 are the times of flight of the first and second fragments, respectively. They have demonstrated their technique for a vast array of molecules [36]. Subsequent studies have implemented this method of lifetime retrieval for other metastable molecules [39,[41][42][43][44]. Recently, making some simplifying assumptions, Larimian et al [47] calculated t d kinematically and retrieved the lifetime for deprotonation of the metastable ethylene dication. They also retrieved the momentum distribution of the fragments using Abel inversion [49] of the position image and discussed possible decay routes. Our aim in this work is to present a more direct approach for extracting information about dissociation in flight of metastable molecules from coincidence momentum imaging measurements, employing basic principles. This method takes advantage of the known symmetries regarding the fragmentation process and allows one to extract not only the lifetime(s) of the metastable molecule but also the momenta, KER, and angular distributions. This information can enable pinpointing of the likely metastable state(s) dissociating in flight, as well as their dissociation mechanisms. Furthermore, this technique is general and can be applied to many different systems that undergo dissociation in flight. While our method is versatile and can provide a wealth of information to deepen understanding of dissociation processes, the focus of this manuscript is the analysis method of retrieving this information from the measurement.

Experimental method
To demonstrate our approach, we examine the same dissociation-in-flight channel as Larimian et al [47]. Namely, we look at deprotonation of metastable ethylene dications, C 2 H 4 2  + C 2 H 3 + +H + (as well as the deuterated equivalent), using the cold target recoil ion momentum spectroscopy (COLTRIMS) technique [27,28]. Laser pulses with central wavelength of 790 nm, 23 fs duration (full width at half maximum (FWHM) in intensity), and peak intensity of about 3×10 14 W cm −2 are used to doubly ionize the ethylene molecules, introduced as a supersonic jet. The pulse duration was measured using second harmonic frequency-resolved optical gating (SHG FROG) [50], and the intensity was evaluated using the kink in the photoelectron spectrum of neon associated with 2U p (where U p is the average quiver energy of the free electron), which represents the transition from direct to rescattered electrons. To determine the 2U p point, we measured the momenta of Ne + recoil ions at low extraction field, following the method detailed in [51,52]. As shown in figure 1(a), we observe the well-known signature of dissociation in flight, a long stripe in the CTOF spectrum. In figure 1(b), we note another signature in a density plot of the ion yield as a function of the measured TOF and position, N(TOF, x), where two 'halos' extend from the light and heavy fragments to the small C 2 D 4 2+ spot. Notably, in both the CTOF and N(TOF, x) spectra, the distributions of the two fragments converge to that of the intact metastable dication. Furthermore, as highlighted in figure 1(a), the predictions of Newton's equations of motion for dissociation in flight, which are marked by the open triangles, agree well with the measured data.

Analysis method and results
To accomplish our goal of directly extracting information about dissociation in flight of a metastable molecule from our measurements, we start from the kinematic equations. The coordinate system utilized is depicted in figure 2(a). We employ principles similar to those that have been applied for years to image collision-and laserinduced prompt molecular dissociation in CTOF measurements [27,28]. The equations of motion become slightly different from those for prompt dissociation to account for the survival time of the dication, t d . For twobody dissociation in flight of a generic metastable dication, AB 2  + A + + B + , we have the following equations in the x direction, which in our case is along the laser beam propagation and transverse to the spectrometer axis: Here, x j is the measured position of fragment j, and x 0 is the initial position of the metastable dication. v 0x is the average initial x-component velocity of the dication in the laboratory frame (v 0 x 0  for a cold jet), v jx ¢ is the xcomponent dissociation velocity in the AB 2+ center-of-mass frame, t j is fragment jʼs time of flight, and m j is its mass. Clearly, t d is needed to properly calculate the transverse momenta of the fragments. Note that the ycomponent equations (transverse to the spectrometer axis and along the jet flow) are similar to the x equations, except that v 0y , the supersonic jet velocity, is not negligible.
At this point, it is interesting to contrast momentum imaging of dissociation in flight to that of prompt breakup, for which t 0 d  . In the case of prompt fragmentation, one can readily see that the transverse equations of motion are decoupled from motion along the z direction, which is parallel to the spectrometer axis. For the problem of dissociation in flight, however, this is not the case, as t d leaves us with more unknowns than equations in the transverse directions. Thus, we need to determine t d from the z-component kinematic equations first.

t d and lifetime determination
How exactly can we retrieve t d from the z-direction motion? First, we write the z-component equations of motion: Here, ℓ is the ion flight distance, a m is the AB 2+ acceleration, v jz ¢ is the z-component dissociation velocity of the jth fragment in the AB 2+ center-of-mass frame, and a j is its acceleration. We proceed to write the equations of motion in a more convenient dimensionless form. To that end, we multiply equation (2) by a 2 m : we replace the first term on the left-hand side of equation (3) with t m 2 . Further, dividing both sides by t m 2 leads us to dimensionless z-component equations of motion. We also write the equation for momentum conservation in the AB 2+ center-of-mass frame: , the velocity of the dication.
The equations above suggest that we can solve for t dm (and hence t d ), as we have three equations and three unknowns, t dm , v z 1 ¢ , and v z 2 ¢ . Combining the above equations of motion and the equation for momentum conservation, we eliminate v jz ¢ by substitution, resulting in an equation that can be solved for t dm . Several subsequent algebraic steps lead us to the following quadratic equation: is the ratio of the mass of the light fragment to that of the heavy fragment. While the simplicity of our derivation and resulting equations makes the calculation of t d and subsequent momentum imaging look quite straightforward, the problem is more convoluted than it initially seems. The quadratic equation for t d has two solutions, which we denote as t d + ( ) and t d -( ) , where the superscripts correspond to the sign that is chosen in the quadratic formula.
For the vast majority of events, it is not clear which solution is correct, as both are physical based on obvious criteria: t d must be real and t t 0 ). This dilemma also thwarts retrieval of v z 1 ¢ and the z-component momenta, for which one needs to properly evaluate t d .
As a noteworthy aside, this dual-solution t d retrieval problem belongs to an extensive family of inverse problems, in which one is trying to retrieve initial conditions from observable parameters. This is an important problem faced in a wide array of fields in science and mathematics, such as medical imaging, x-ray crystallography, optics, geology, acoustics, and many others [53][54][55]. Frequently, the solution to an inverse problem is not unique, as is the case for dissociation in flight.
To address the inverse problem at hand, we utilize symmetry concepts. For any given t d , reflection symmetry of v z 1 ¢ about 0 is expected, as the light fragment is equally likely to be ejected in either the forward or backward directions. Thus, one could select t d + ( ) or t d -( ) randomly with equal probability and fit an exponential decay function to the resultant N(t d ) distribution to get a lifetime τ.
In the case of dissociation in flight, however, the symmetry of the overall v z 1 ¢ distribution should be broken to some degree due to the lifetime of the metastable molecule. That is, if τ is the lifetime of the dication, the ratio of the number of molecules that survive for times t d Thus, one could correct for the symmetry breaking using this factor. The extent of this correction is determined by the magnitude of t t . Since t d + ( ) and t d -( ) are typically not dramatically different compared to the lifetime τ, this correction is small. Thus, we use a self-consistent approach in which we start as suggested above, by choosing t d randomly with equal likelihood for each event. Recall that this choice is exact in the limit t t 1 is then fitted to the resulting N(t d ) distribution to retrieve the lifetime. This lifetime allows computation of the aforementioned factor N N + -( ) ( ) , given in equation (6), which is then used to weight the choice of t d + ( ) or t d -( ) in the next iteration. The obtained N(t d ) distribution is again fit with an exponential decay function to retrieve a more accurate lifetime, again allowing calculation of a new weighting factor for the choice of t d This process is repeated until the lifetime τ converges. Note that for a given iteration, the choice of t d + ( ) or t d -( ) and the fitting procedure is repeated for multiple trials to account for the finite sample size of our data. Also, the lifetime τ used to compute the weighting factor for the subsequent iteration is the mean value of those obtained in the multiple trials. For more details about our iterative approach, visit appendix B.
When this self-consistent method is applied to our data, we find good convergence within just a few iterations. We performed the analysis on the C 2 H 4 and C 2 D 4 data to explore the possibility of isotopic effects.
As mentioned, the two-term exponential fit indicates that at least two metastable states are responsible for the observed dissociation in flight. The single lifetime for C 2 H 4 2+ reported by Larimian et al, 498±12 ns [47], lies between the two lifetimes that we have measured. This discrepancy could be due to a number of reasons, such as differences in the laser pulse parameters, the method used to compute t d , or the t d range chosen for the exponential decay fit. Note that when we perform a single exponential decay fit, shown in figure 3(a), we obtain a lifetime of 491±19 ns, consistent with the previous measurement.
We also note that our measurements suggest a possible small isotopic effect in the shorter lifetime, 1 t . The difference between the two shorter lifetimes is on the level of 2.2σ, while the longer lifetimes are the same within the measurement uncertainty. As dissociation in flight of ethylene dications is a low-rate channel, we expect that higher statistics data would make the presence or absence of an isotopic effect in the lifetimes more clear cut.
Since we currently lack the good quality electronic structure information on these molecular ions needed to understand this isotopic effect and also to keep the focus on the method, we limit this discussion to highlighting the rich information afforded by our technique. Note that in both plots, the statistical error bars are smaller than the symbols. Two-term exponential decay fits, plotted in blue, agree nicely with the data. The gray dashed line shown in (a), which doesn't agree with the data, is a single-term exponential decay fit. will be assigned incorrectly. Therefore, this method of computing t d cannot be used for momentum imaging, as it is done on an event-by-event basis. As such, we need a single value of t d for each event, even if it is approximate.

Momentum imaging
To approximate t d , we neglect v z 1 ¢ in equations (4), as the term containing this quantity is typically on the order of a few percent compared to the other terms, as further detailed in appendix C. Moreover, as shown in the same appendix, the error that this approximation introduces in the retrieved t d and the transverse momenta is also estimated to be at most a few percent. Having neglected v z 1 ¢ , the equation for t d becomes linear and thus has a single solution for each event, given explicitly by t t t t t Employing this approximation of t d , we have computed the transverse momentum of the D + fragments from C 2 D 4 2+ , shown in figure 4. The distribution of the radial momentum, p p p r x y 1 1 2 1 2 = + , shown in figure 4(b), agrees well with the functional form for the projection of an isotropic distribution onto a plane [56]. To verify that the momentum distribution is isotropic, the complete 3D momentum distribution is needed. Therefore, the p 1z component should also be measured.
Is there any way to retrieve the z-direction momentum? Recall our initial fundamental problem in evaluating v z 1 ¢ that is associated with the t d + ( ) and t d -( ) solutions. Now, we have gone even further and neglected v z 1 ¢ entirely, eliminating the possibility of recovering p 1z directly from the measurement. It is important to note that the polarization is typically aligned along the spectrometer axis (z direction) in COLTRIMS measurements in order to reduce losses of fast fragments, which are usually ejected along the laser field. This choice leads to equivalent x and y momentum components due to the axial symmetry about the laser field and prevents the direct determination of p 1z . To retrieve the missing information, i.e., the p 1z momentum component along the laser field, we take advantage of this axial symmetry and align the laser polarization along the y axis. Under these conditions, the measured p 1y distribution is along the laser polarization, while p 1x is transverse. Moreover, the 'lost' p 1z distribution can be recovered from the measured p 1x distribution by taking advantage of the axial symmetry about the laser polarization. Under ideal conditions, this measurement is sufficient to retrieve the complete 3D momentum distributions of the fragments. In many cases, however, imperfections like spatial non-uniformities in the detector response may bias the results.
To circumvent this issue and verify that the momentum distribution is isotropic, we performed two measurements with the polarization along the z and y directions, as illustrated in figures 5(a) and (b), respectively. Note that while the angular distributions drawn in this figure do not resemble the isotropic distribution we measure in the present experiment, they help to better convey the difference between the two measurement schemes. In the first measurement with the polarization along z, the momentum distribution along the polarization cannot be retrieved because v z 1 ¢ has been neglected. In the second measurement with the polarization along y, the momentum distribution parallel to the laser polarization can be retrieved directly, while the complete transverse momentum distribution can be recovered by using the axial symmetry about the laser field as discussed above.
Let us define j as a rotation angle in the xy-plane in both measurements. This angle is sketched on figure 4(a). We denote N ( 1 j ) as the distribution obtained in the first measurement (polarization along z) and N ( 2 j ) as the distribution found in the second measurement (polarization along y). Computing the ratio of these distributions yields the result shown in figure 5(c). Note that the position-dependent detection efficiency cancels out in this ratio, thus eliminating the impact of detector imperfections. The ratio shown in figure 5(c) is rather flat, directly demonstrating that dissociation in flight yields an isotropic momentum distribution. This distribution is likely the result of t d being much longer than the rotational timescale of the molecule. Thus, information about the initial alignment of the molecule with respect to the laser polarization is lost, and the resulting distribution is isotropic.

Kinetic energy release
Finally, just as accessing the z-component momenta is problematic, so too is retrieving the KER on an event-byevent basis. To obtain a KER distribution, we utilize a method based on the onion peeling technique, which has been widely used to analyze photofragment images [57][58][59][60]. The transverse momentum distribution, which was obtained using the method described in the previous section, serves as our projected 'onion,' which we slice along the p x direction. Some sample slices for dissociation in flight of C 2 D 4 2+ are shown in figure 6(a). As we carry out iterative onion peeling subtraction on the slices, the counts that are 'peeled' away are allocated into the appropriate KER bins to accumulate a distribution. For a given iteration being performed on a particular slice, the KER is found using the edges of the slice, as only events with the maximum KER reach this region. We employ this technique for measurements in which the laser polarization is along the y axis, even though the isotropic nature of the distribution would allow one to use data from either of the previously discussed measurement schemes. Further details about this KER retrieval method will be described in a forthcoming publication about the dissociation in flight of metastable carbon dioxide dications [61].
The KER distribution obtained using our 'sliced' onion peeling technique for dissociation in flight of C 2 D 4 2+ is presented in figure 6(b). We estimate the uncertainty in the obtained KER to be about 0.3 eV 2 . The centroid of the KER, at about 4.2 eV, is in good agreement with that obtained by Larimian et al [47]. Finally, as an alternative method of retrieving the KER, we compare the measured N p r 1 ( ) with several simulated N p r 1 ( ) distributions corresponding to Gaussian KER distributions with different centroids and widths, as illustrated in figure 6(c). As Figure 5. (a) and (b) Schematic angular distributions. Note that these do not reflect our measured angular distributions, but they help better illustrate the concept of our two-measurement scheme. (a) First measurement: the laser polarization is aligned along the z direction (spectrometer axis). As v z 1 ¢ is neglected, momentum information parallel to the polarization is not accessible in this configuration. (b) Second measurement: the laser polarization is parallel to the y axis (jet direction), and thus, the momentum distribution along the laser polarization can now be retrieved. Additionally, exploiting the azimuthal symmetry about the laser polarization allows for retrieval of the p z distribution in this scheme, as N(p 1z ) N = (p 1x ). (c) Ratio of the N(j) distributions for the two measurement schemes. N 1 j ( ) corresponds to the j distribution measured with the laser polarization parallel to the time-of-flight axis, z, and N 2 j ( ) is measured after rotating the polarization to be along the jet direction, y. The dotted red line corresponds to the average value of this ratio.
can be seen, the simulated N p r 1 ( ) distribution with a KER centroid of 4.15 eV and a FWHM of 0.5 eV agrees well with the measured N p r 1 ( ). Furthermore, this result supports that obtained using the 'sliced' onion peeling approach. As mentioned, the KER can supply a great deal of insight into dissociation pathways, e.g., [17,25,26]. This pursuit, however, is beyond the scope of this paper, which focuses on the method.

Summary and outlook
In summary, we have developed a method to study dissociation in flight of metastable molecular ions using coincidence momentum imaging measurements. Our approach, which supplies valuable information about the relevant metastable states, including the lifetime(s) and momentum distributions of the dissociating fragments, has been realized through the application and symmetries of the relevant kinematic equations.
Encountered hurdles such as the inverse problem of choosing t d + ( ) or t d -( ) and the related problem of retrieving v z 1 ¢ have been addressed by exploiting symmetries of the fragmentation. The readily expressed forward-backward symmetry breaking in the v z 1 ¢ distribution was used in a self-consistent manner to obtain the N(t d ) distribution of the sample and hence the lifetimes of the metastable states dissociating in flight. This analysis allowed us to find two lifetimes in the deprotonation of metastable ethylene dications and a possible isotopic effect in the shorter lifetime.
The necessity of a single t d value for each event to obtain the momenta was fulfilled by neglecting v z 1 ¢ , an assumption which we have shown to be on solid ground, as it introduces minimal error in the calculation of t d and the momenta. Furthermore, the azimuthal symmetry about the laser polarization was exploited to obtain all the components of the momentum rendered unretrievable by the t d inverse problem.
Finally, while we have demonstrated this method for the specific case of deprotonation of metastable ethylene dications formed by intense femtosecond laser pulses, this technique is applicable to coincidence measurements on a variety of metastable molecular systems dissociating in flight, which could also be formed via other means, such as x-ray photoabsorption or fast charged particle impact. The long, curved dissociation-in-flight stripe in the CTOF spectrum fairly closely follows a third-order polynomial dependence as a function of t 1 . To more effectively select true events while suppressing the contribution of random pairs (i.e., improve the signal to 'noise' ratio), we use the coefficients of a third-order polynomial fit to the curved coincidence stripe to straighten it and then apply a simple rectangular gate. More specifically, t ⊥ , the perpendicular distance from each (t 1 , t 2 ) data point to the polynomial fit is calculated. A Gaussian function is fit to the N t( ) distribution, and a±3σ gate is applied to the data, as shown in figure A1(a). Figure A1(b) shows the straightened stripe and gate around it. Note that for C 2 H 4 , the left and right gating bounds for t 1 were chosen to be 825 ns and 2400 ns, respectively, and for C 2 D 4 , they were 1200 ns and 2400 ns, respectively. These data selection schemes allow one to avoid the prompt breakup region, as well as the area of the CTOF spectrum near the end of the dissociation-in-flight stripe, where 'noise' and other channels dominate. As we are applying this analysis to experimental data of a relatively weak channel for which statistics are limited, one may worry about the robustness of the proposed method that relies on random number generation. To address this issue, as mentioned in section 3.1, each iteration of the lifetime determination procedure consists of multiple trials. That is, for each iteration, the analysis is simply repeated multiple times (each time with a randomly selected seed). Each trial uses the same solution choice weighting scheme. An exponential decay function is fit to the resulting N t d ( ) distributions to retrieve lifetimes for each trial. At the end of an iteration, the amplitudes N 0 To obtain a single value of t d needed for evaluating the momentum of each event, an approximation is necessary. We start with the first expression in equations (4) and divide both sides by t t jm dm to obtain Figure A1. Illustration of the procedure used for selecting the CTOF stripe. (a) N t( ) for C 2 H 4 , where t ⊥ is the perpendicular distance of a measured (t 1 ,t 2 ) pair to a third-order polynomial fit to the CTOF stripe. Shown in red, a Gaussian function is fit to N t( ). As indicated by the purple dashed lines, a±3σ gate about the centroid of this Gaussian is then set on the data. (b) Density plot of ion yield as a function of t 1 and t ⊥ , N t t , 1( ). In addition to the±3σ gating, lower and upper bounds are chosen for t 1 to minimize the contribution of noise and channels other than the dissociation-in-flight channel. The surviving data inside the gate (dashed purple box) is then used to retrieve the lifetimes and momenta.
This equation is solved for t dm after neglecting the second term on the right-hand side that contains v jz ¢ . How valid is neglecting this term? To explore this question, we perform simulations using a few typical input random functional distributions, including a single-term exponential decay N t d ( ) distribution with 900 t = ns, a Gaussian KER distribution centered at 4 eV with a 1 eV FWHM, and an isotropic angular distribution. Let us denote the terms on the right-hand side of equation (9) as '1st term,' '2nd term,' and '3rd term,' in left-to-right order. We evaluate the quantity R ≡|2nd term|/(1st term + 3rd term), shown in figure C1(a). For 99% of the events, R 10 < %. Furthermore, the error that neglecting the second term introduces into the recovered t d and p 1x , shown in figures C1(b) and (c), is also reasonably small. For t d , about 77% of the simulated events lie below the 5% error level, and for p 1x , about 98% of the events lie below the same error level.
Since the validity of neglecting v z 1 ¢ depends on the magnitude of v z 1 ¢ , and for our simulations we have assumed values of this quantity approximately matching the measured ones, it is reasonable to explore how large v z 1 ¢ can be before the approximation breaks down. Thus, we performed simulations with a few larger KER values (and hence larger maximum values of v z 1 ¢ ). Even for a high KER of 20 eV (Gaussian distribution with 1 eV FWHM and the same lifetime and angular distribution as before), 99% of the simulated events have R 20 < %, as shown in figure C2(a). Moreover, as shown in figures C2(b) and (c), 75% of the events have <20% error in the  retrieved t d , and about 98% of the events are below the same error level in p 1x . Given this extreme example, we are thus assured that neglecting v z 1 ¢ for momentum computation is a reasonable approximation for our case. It is noteworthy that in general it is not merely the KER (and hence the maximum v z 1 ¢ ) that is important for consideration but the ratio v v jz m ¢ , as can be readily seen in equation (9). Recall that v m is the velocity of the dication. Therefore, in certain cases, it may also be desirable to increase the spectrometer voltage in the experiment to increase v m and thereby improve the validity of this approximation.
In solving for t 1 and t 2 , we choose the positive root in the quadratic formula because choosing the negative root yields negative t 1 and t 2 values or t t d 1 < , which is unacceptable, in contrast to the t d equation for which both roots can make physical sense. The values of t 1 and t 2 are truncated to varying levels of precision in ns, simulating the possible digitizer accuracy. Then a random fraction is added to the truncated number, as is done for the measured (digitized) data, to match the original number of digits. For example, if t i is a time-of-flight value truncated to n decimal place(s), the new TOF after adding the random fraction is t t r 10  While our simulations have proven instructive in identifying influential error sources such as those related to the measured time-of-flight values, note that the errors presented in the body of this manuscript reflect those evaluated by statistical means.