Robust inversion of time-resolved data via forward-optimisation in a trajectory basis

An inversion method for time-resolved data from ultrafast experiments is introduced, based on forward-optimisation in a trajectory basis. The method is applied to experimental data from x-ray scattering of the photochemical ring-opening reaction of 1,3-cyclohexadiene and electron diffraction of the photodissociation of CS2. In each case, inversion yields a model that reproduces the experimental data, identifies the main dynamic motifs, and agrees with independent experimental observations. Notably, the method explicitly accounts for continuity constraints and is robust for noisy data.


Introduction
Inverse problems are important across science and engineering, [1][2][3] and have a long history at the borderline between physics and mathematics. 4 They concern the determination of key target characteristics from an observed set of outputs, for instance determining the shape of a vibrating membrane from the spectrum or reconstructing an object from a photograph. Important applications appear in medical imaging and material characterisation. [1][2][3]5,6 In contrast to the for-ward, also known as the direct, problem, inverse problems are often mathematically illposed, meaning they may exhibit instability, be underdetermined, or lack unique solutions. These challenges must be overcome via new computational schemes, regularization techniques, objective functionals, and experimental procedures.
For molecular structure determination the inverse problem is central. The structure of a molecule is generally underdetermined by the experimental data and the molecular model is optimised (bootstrapped/refined) subject to auxiliary constraints. Examples include the fitting of spectroscopic data to a Morse oscillator or Dunham coefficients in a diatomic, 7 bond-distance, bond-angle, and torsional constraints during x-ray crystallography or NMR refinement, 8 and the use of supplementary experimental or computational data in gas-phase electron diffraction structure determination. 9 In recent years, ultrafast imaging of photoexcited molecules have developed rapidly, with experiments increasingly capable of tracking molecular dynamics on fundamental timescales, observing phenomena such as vibrations, bond breaking, or charge transfer. [10][11][12] The techniques include spectroscopy (e.g. TRPES), [13][14][15][16] ultrafast electron diffraction (UED), 17,18 ultrafast x-ray scattering (UXS), 19,20 Coulomb explosion imaging (CEI), 21 and others. Beyond fitting data to heuristic rate models, extracting a more detailed molecular model is often challenging. From the perspective of structure determination, for instance, a minimal structural model would consist of a sequence of structures fulfilling at least some degree of temporal continuity. A complete model would conceivably consist of the full temporal evolution of the molecular wavefunction with all associated time-dependent nuclear and electronic distributions and populations.
Part of the issue is that most experiments only provide partial information of complex and multidimensional dynamics, which makes the complete problem severely underdetermined. One obvious strategy is to combine data from complementary experiments. For instance, spectra are sensitive to electronic populations while scattering experiments primarily, but not exclusively, detect the molecular geometry. 12,[22][23][24][25] Another important strategy is comparison to theoretical results. Least but not last, inversion algorithms that account for at least some of the physics stand to play an increasingly important role.
The challenges involved in inversion in this context are so significant that the default option is to rely on comparison between theory and experiment only. The more reliable variety of this approach makes the comparison in the observable space, with observables forward mapped (i.e. predicted) from simulations. However, such comparison only provides a qualitative understanding of the data and leaves little recourse if agreement between experiment and theory is not achieved. Inversion algorithms capable of providing a rapid, first-order, assessment of the observed dynamics would carry great value for experimental progress, even if by simply determining the range of molecular models that are, in principle, commensurate with the observations.
Therefore, efforts to tackle the inverse problem for ultrafast dynamics are intensifying. One recently developed approach used for the interpretation of ultrafast x-ray scattering data samples a large pool of randomly gen-erated molecular structures to find the best fit, 12,26,27 and genetic algorithms have been proposed to invert electron diffraction data by exploring molecular structure via consecutive in silico mutations. [28][29][30][31][32] These approaches include phase retrieval algorithms for the structural determination, but this requires that the molecules are at least partially aligned. [33][34][35][36][37][38] The preceding approaches are applicable to both static and dynamic data, but methods specifically aimed at dynamic data have also been attempted. This includes machine learning approaches that use a variational recurrent neural network trained on temporally correlated frames 39 and an approach that systematically perturbs the molecular structure at each step, starting with the well-known initial structure at time zero. 40 In this paper, we present a detailed discussion of an approach for modelling timedependent data which explicitly includes the time-evolution, thereby ensuring that continuity, inclusion of known initial or final structures and any other known physics. It proceeds by optimising the weights of semiclassical trajectories from quantum molecular dynamics simulations against experimental data. In essence, this biases the simulations towards experimental observations. The approach is general, and provides a platform for merging data from several complementary experiments, with initial focus on the analysis of experimental data from ultrafast x-ray and electron scattering experiments. [41][42][43] The aim of the paper is to provide the first unified and general presentation of this methodology, and to critically evaluate its performance, establish best-practice, and explore avenues for improving the methodology further.

Forward optimisation
The time evolution of any observable X(t) can be calculated from the molecular wavefunc-tion Ψ(t) via a forward mapping M , where the index j identifies the type of measurement, which could be anything from photoelectron spectra to electron diffraction or Coulomb explosion imaging. The observable X j (t) may be resolved with respect to several implicit variables, for instance photoelectron kinetic energy and angular distribution. In contrast to the inverse problem, the forward mapping is mathematically well-conditioned and does not suffer stability or underdetermination issues. We therefore proceed to tackle the inverse as a forward optimisation problem, [41][42][43] where our goal is to find a model molecular wavefunction Ψ mod (t) that yields predicted observablesX j mod (t) that reproduce the experimental observables X j exp (t). The calculation of the observable consists of two steps, (2) where the first step is the forward mapping M j from eq 1, which produces the theoretically predicted signal, and a second apparatus mapping step S j which replicates the effect of the measurement apparatus on the data, for instance due to limited time resolution. Such distortions are unavoidable despite that the experimental data X j exp (t) will have been through extensive preprocessing to remove known artefacts and distortions.
The optimisation proceeds by modifying the model wavefunction until close agreement between the experimental and predicted data is achieved, as measured by the target function F , where the index j runs over the different types of experiments included. In practice it is common for the integration over time t to be replaced by a summation over a temporal grid {t i }. The α j is a regularisation factor that must be included when data from several different types of experiments is considered. The factor is determined from the numerical profile of each of the different data sets, and scales data so that it can be combined in a balanced manner. The reminder of the paper is organised such that Section 2.2 presents the reference molecular wavefunction, Section 2.3 presents the parameterised model wavefunction and the target function that result from the model function, and Section 2.4 discusses the forward mapping with an emphasis on scattering experiments. In Section 3, the range of applicable numerical optimisation techniques are discussed, as well as the S j mapping to match the predicted observables to the experiment. The preprocessing of the experimental data is touched upon in Section 4 on data treatment and in the Supplementary Information (SI). In the Results, Section 5, two applications to recent ultrafast electron diffraction and ultrafast x-ray scattering data are examined in detail and the convergence and the resulting interpretation are assessed.

Molecular wavefunction
The time evolution of the molecular wavefunction Ψ(t) is governed by the time-dependent Schrödinger equation, whereĤ is the molecular Hamiltonian and r and R the electronic and nuclear coordinates. The molecular wavefunction can be expanded in the Born-Huang form as, 44 where χ k (t, R) are time-dependent nuclear wavepackets which propagate on electronic eigenstates ψ k (r; R) that depend parametrically on the nuclear coordinates R. 1 In prac-tice, the expansion in eq 5 is truncated to include only the N s electronic states visited during the dynamics of interest. A wide range of numerical techniques to solve eq 4 exist. Accurate methods such as numerical grid propagators 45,46 and multiconfigurational time-dependent Hartree (MCTDH) 47 require precalculated potential energy surfaces which are not feasible for most molecules of interest. An alternative is direct dynamics (dd) methods that expand the molecular wavefunction by classical or semiclassical trajectories. Examples include surface hopping (SH), ab initio multiple spawning (AIMS), ab initio multiconfigurational Ehrenfest (AIMCE), and direct dynamics variational multiconfigurational Gaussians (dd-vMCG). 48-51 In a general form, the direct dynamics wavefunction Ψ dd (t) can be expressed as, with c n (t) the expansion coefficient for each of the N TBF trajectory basis functions (TBFs) given by, where the parenthesis contains the electronic states |ψ k and their populations |a n k (t) | 2 . The nuclear basis functions χ n (t) follow phasespace trajectories (R n (t), P n (t)) where R n (t) and P n (t) are the nuclear positions and momenta, respectively. The equations of motion, which govern the trajectories, populations, and auxiliary coefficients such as for instance wavepacket width coefficients, are different for each method. Thus, the details of each specific dd wavefunction will vary. For instance, the nuclear wavepacket is a Gaussian 52 in all methods except SH where it is a δfunction, for SH and AIMS only one electronic state is occupied by each TBF at any given time, while Ehrenfest methods such as AIMCE occupy several states simultaneously, in AIMS (AIMCE) spawning (cloning 53 ) increases the number of TBFs as the simulation progresses, etc.

Model wavefunction
We adapt the wavefunction in eq 6 as our parameterised model wavefunction by rescaling the expansion of the TBFs. In effect, the TBFs can be thought of physically reasonable constraints on a system far from equilibrium, which in addition automatically fulfill continuity requirements. The resulting model function is given by, where w n (t) are the weights for the TBFs which are adjusted to bias the theoretical model towards the experiment, subject to a normalisation condition 1 = Ψ mod (t) Ψ mod (t) . If observables are taken to only depend on the nuclear coordinates and are calculated in the diagonal zeroth-order bracket-averaged Taylor expansion (BAT) approximation (see ref 54 ), as we will do here, the coefficients w n (t) and c n (t) can be taken to be timeindependent with the normalisation straightforwardly given by n w n = 1 for c n ≡ 1. For SH wavefunctions, this is always the case. The target function in eq 3 then becomes a function of the time-independent weights w = (w 1 . . . w N TBF ).
As a final aside, we note that the current approach can also be used when the target is best described by a density matrix, by augmenting the model with additional model wavefunctions weighted by their population factors.

Forward mapping
The forward optimisation exploits that the direct problem, i.e. the forward mapping, is mathematically well-defined and has stable solutions. For each type of experimental observable the mapping M j will be different and based on different theoretical approximations and computational techniques. For instance, there is an extensive body of work on the prediction of time-resolved photoelectron spectra 55,56 with available techniques ranging from approximate Dyson orbitals calculations 25,57 to the highly accurate R-matrix method. 58 The method used to calculate the x-ray scattering and electron diffraction signals in the current paper will be discussed in Section 3.
(b) Dissociation of CS 2 yields singlet or triplet sulphur atoms. Figure 1: Schematic of the two photochemical processes probed in the experiments for which data is inverted in this paper, via a basis of trajectory basis functions (TBFs).

1,3-cyclohexadiene ring-opening
The ring-opening reaction of 1,3-cyclohexadiene (CHD) to 1,3,5-hexatriene (HT), shown schematically in Figure 1a, is a prototypical Woodward-Hoffman photoinduced electrocyclic reaction. 59 It has been the target for a large number of pioneering time-resolved experiments that include UXS, UED, and timeresolved spectroscopies. 23,25,41,42,[59][60][61][62][63][64][65][66][67] Upon absorption of a 267 nm photon the molecule undergoes a π → π * transition to the steeply sloped 1B electronic state, glancing the conical intersection with the 2A electronic state while staying on the adiabatic potential energy surface. Passage through the conical intersection to the electronic ground state returns the molecule to the original ring-closed CHD or breaks a C−C bond to yield the ring-open HT. A distribution of various cis(Z)/trans(E) HT isomers are observed given the high internal energy of the system. Semiclassical trajectories, TBFs, for the dynamics are calculated using the ab initio Ehrenfest (AIMCE) method. 41,42 The electronic structure calculations use the ab initio package MOLPRO, 68 which supplies the forces and the nonadiabatic couplings at 3SA-CAS(6,4)-SCF/cc-pVDZ level of theory for the ground and excited states. A set of 100 TBFs, propagated for 200 fs, with initial conditions sampled from a Wigner distribution in the Franck-Condon region of equilibrium ground state CHD provide the basis for the trajectoryfitting procedure.
Semiclassical trajectories (TBFs) to be used in the forward optimisation are calculated using the SHARC surface-hopping package. 80,81 The trajectories start on the optically bright 1 Σ + u state with a total energy corresponding excitation by a 200 nm pulse and initial coordinates sampled from the ground state Wigner distribution. The forces and nonadiabatic couplings are calculated using SA8-CASSCF(10,8)/SVP electronic structure theory using MOLPRO 68 and 197 trajectories are propagated for 1 ps. Further information can be found in SI Section ??.

Scattering observables
The forward mapping required to compare the model to the UXS and UED data is summarised below. Although scattering signals can be calculated from first principles for gasphase molecules [82][83][84][85][86][87][88][89][90] and recent x-ray scattering measurements have demonstrated that given sufficient accuracy different electronic states can be resolved, 10,91-93 the experimental data considered in this paper can be modelled using the significantly simpler independent atom model (IAM), which approximates the scattering as a summation over the coherent scattering from isolated atoms.
The instantaneous energy-integrated total scattering cross section into the solid angle dΩ is then, 54,94,95 where q = k 0 − k 1 is the scattering vector expressed in terms of the wave vectors of the incoming and outgoing photon or electron. 2 For x-ray scattering, the prefactor (dσ/dΩ) x is the differential Thomson scattering cross section (dσ/dΩ) Th , here made to include the |e 0 · e * 1 | polarisation factor, while for electron scattering it is the Rutherford cross section (dσ/dΩ) Rh which here includes the scaling factor s −4 . 96,97 Note that the expression above does not account for the duration of the x-ray pulse, which is accounted for via the temporal convolution discussed in Section 3.3.
Using the diagonal BAT approximation and assuming time-independent expansion coefficients w = 1, the scattering intensity I(q) can be calculated as, w n I n q, R n (t) . (11) More general expressions that account for the full wavefunction in eq 6, including the nonlocal nature of the individual TBF nuclear wavepackets, have been derived previously. 54 Continuing with the current simplified form, sufficient for our present needs, yields the scattering intensity as, I n (q, R n (t)) = |f (q, R n (t))| 2 + S inel (q), (12) where S inel (q) is the inelastic scattering, which is independent of molecular geometry and given by an incoherent summation over the atomic contributions, with N at the number of atoms. The corresponding elastic contribution is given by the form factor f (q, R n (t)), where f A (q) are the atomic form factors and R nA (t) the position vector for atom A in trajectory n. Both f A (q) and S A (q) are tabulated. 98 The form factors for electron scattering are where Z A is the atomic number and f x A the x-ray scattering form factor. 96,97 For high energy electron scattering, notably MeV-UED, it is sometimes necessary to use form factors with relativistic corrections. 99,100 When the target is a gas of anisotropic molecules, rotational averaging of eq 14 results in, 101 between atoms A and B in trajectory n.

Apparatus mapping
In this section we discuss the apparatus mapping S j in eq 2 that is required to match the forward mapped signal X j mod (t) from eq 9 to the experimentally observed signal.

Temporal alignment and convolution
In the scattering experiments, the time-zero is roughly calibrated by the instrument, with the exact time-zero inferred from the observed data. The absence of independent validation means that one must check the alignment of the temporal axes in the experiment and the model. 3 We define the relationship between the experimental and model time axes as t = t + t 0 , where t is the experiment, t the model, and t 0 the time shift. The temporal alignment t 0 is one of the global parameters optimised. A Gaussian convolution in the temporal domain is included as, where G(t) = b c exp (−a c t 2 ) mimics the instrument response function, with the normalisation constant given as b c = 4 ln 2/(τ 2 c π) and τ c the full-width half-maximum (FWHM). The convolution mainly equates to the crosscorrelation of the pump and the probe, effectively compensating for the δ(t)-pulse excitation approximation used to generate the TBFs. In practice, the instrument response function also accounts for other limits on the temporal resolution, such as temporal jitter.
Finally, we note that the experimental resolution with respect to the momentum transfer q is such that no convolution of the model is required. The amount of structural information in the signal is limited by the q-range, q ∈ [q min , q max ], measured in the experiment. 42,95

Percent difference signal
The experimental signal is considered in the percent difference form to minimise systematic multiplicative errors, where I on (q, t) is the optically pumped 'laseron' signal and I off q, t j ≤ t 0 is the static 'laser-off' reference signal measured at delay times t t 0 . The theoretical equivalent is calculated from the model wavefunction Ψ mod (t) and defined as, where the excitation fraction γ scales the intensity according to the implicit degree of excitation. The theoretical 'laser-off' signal I th off (q) is calculated from a suitable reference geometry, or more accurately using a Wigner distribution of the system in its ground state at the equilibrium geometry. In some cases, it may be necessary to modify the definition in eq 18 to scale the signal by the ratio of the integrated intensity of the 'laser-on' and 'laseroff' signal, or replace the uniform excitation fraction γ with a q dependent excitation fraction γ (q). Such modifications are discussed in the SI.

Target function and confidence matrix
The experimental percentage difference signal %∆I exp (q, t) and the signal predicted from the model %∆I mod (q, t, w) give the target function introduced in eq 3 the following specific form, where w = (w 1 , w 2 , . . . w N TBF ) are the normalised trajectory weights w = 1, and c consists of additional global parameters to be optimised, such as e.g. the excitation fraction γ and the time-shift t 0 . In principle, more complicated global parameters designed to off-set shortcomings in the quality of the TBFs could be included, such as time-warping to offset inaccuracies in the kinetic energy. The double sum runs over all experimental data points, identified by q i momentum transfer and the temporal coordinates t j = t j + t 0 .
To avoid overfitting of inherently noisy experimental data, the target function includes a confidence matrix. The matrix p conf (q i , t j ) weights data points, identified by their value of the momentum transfer q i and time bin t j , according to the experimental confidence in the accuracy of that data point. For instance, data points subject to poor statistics or systematic errors are given smaller weight while points that are known or expected to be accurate are weighted accordingly. The exact form of p conf (q i , t j ) depends on the data set, as discussed in the SI.
If needed, one can define a confidence threshold p min conf , such that all points p conf q i , t j ≤ p min conf are zero, thus excluding them from the optimisation. The advantage of excluding poor quality points from the optimisation must be balanced against the reduced size of the experimental data set.

Global optimisation
In the ideal scenario, the experimental signal has sufficient quality (high temporal resolution, large q-range, and excellent signalto-noise) that simultaneous optimisation of all model parameters is possible. Since the number of global parameters is small, the most straightforward approach is to optimise the trajectory weights w for fixed global parameters c, and to determine from a scan of reasonable values of c the best overall solution.
For each particular c, a pool of N init initial weights {w} init are Monte Carlo sampled with the values of the individual weights in each initial set w i ∈ [0, 1]. The target function given in eq 19 is then minimised utilising a non-linear trust-region reflective leastsquares algorithm for each initial set. Among the optimisations that converge to a local minimum, within the tolerance constraints, the best minimum is considered a candidate for the global minimum. The number of initial conditions N init generated in the weight space is increased until apparent convergence is achieved. The convergence of the optimisation with respect to N init is discussed in SI, Section ??. In challenging cases, an iterative approach using more targeted sampling can be used, which is discussed in the SI.

Two-step optimisation
When the experimental data is of lower quality, more stable solutions are found via a twostep optimisation procedure. The first step ensures that global parameters, such as t 0 and τ c , are determined as accurately as possible before, in the second step, optimisation of the trajectory weights is attempted.

Step 1: Global parameters
First, we identify the strongest features in the data by inspection of the confidence matrix p conf (q i , t j ). This allows us to minimise the negative effects of noise on the optimisation. Selecting the strongest feature, we define the net integrated percentage difference signal as, where q conf min and q conf max are the bounds on the section of highest confidence. The global parameters c are then optimised against this integrated signal, %∆I int exp (t). Since this step must be performed independently from the optimisation of the trajectory weights w, we assume that the t ≈ t 0 wavefunction can be approximated by an equally weighted sum of TBFs. This is reasonable given the limited dispersion and dephasing in the wavepacket at early times.
The t 0 is inherently linked to τ c and γ since the pulse width τ c affects the onset of the signal and γ scales the strength of the signal. Therefore, for different combinations of t 0 and τ c , the sum of square error between the experimental and model integrated difference signal in eq 20 is minimised with the γ excitation fraction as the free parameter.

Step 2: TBF weights
The second step in the optimisation takes the best values of (t 0 , τ c , γ) from the previous step with the goal of identifying the optimal trajectory weights w. Since the γ in step one is determined on the assumption of equally weighted TBFs, the parameter γ is reoptimised alongside the weights w in this second step. If the final value of γ diverges significantly from its initial value, this may indicate that the parameters determined in the first step are not optimal. In order to ensure self-consistency, it is advisable to repeat the second step optimisation for several different sets of global parameter values (t 0 , τ c , γ) which correspond to good fits in the first step. The procedure is complete once the best sets of these values from step one and two agree, and when the target function F (w, c) has converged.

Metrics of fit quality
In addition to the target function F (w, c), it is useful to have other measures of the quality of the fit. The improvement in overall agreement between the unweighted and weighted model function can be quantified by the relative absolute error (RAE), where N t is the number of time steps, N q the number of points in q, and %∆I th (q, t) the theoretical signal calculated from the theoretical wavefunction not subjected to any bias or optimisation. The RAE measure, as defined above, is independent of the size of the N t ×N q grid. Another helpful metric is the root mean squared error (RMSE) defined as, Particular solutions can also be characterised using the variance of the N TBF weights from a single optimisation, i.e. the spread of weights resulting from one of the N init initial conditions (typically the initial condition that yields the best solution for a specific set of global parameters), defined as, where w is the mean of the weights. In addition, to compare different solutions we define the distance D b from the overall best set of optimised weights, which describes the distance for a particular solution b from the best set of optimised weights.

Experimental data 4.1 Ultrafast x-ray scattering
The ultrafast x-ray scattering (UXS) data for the ring-opening reaction of CHD shown in Figure 2 is taken from ref 42 (further details in the SI, Section ??). The confidence matrix p conf in eq 3 in this case is based on the number of photon hits per frame (SI, eq ??). Due to a long interaction region in these experiments, we include a q-dependent excitation fraction γ (q) which is indirectly optimised by allowing it to be uniformly scaled by a factor x such that it yields the scaled excitation fraction γ x (q) (for more details see the SI). Accounting for the convolution of the signal as in eq 16 and the subsequent temporal binning into bins of size ∆t = 25 fs, the total length of the signal used becomes 275 fs. The quality of the experimental data allows a one-step global optimisation in which we scan over a range of t 0 values and x scaling factors, where t 0 ∈ [−38, −14] and x ∈ [0.7, 1.3]. The duration of the pump and probe pulse were measured from experiment as 60 fs and 30 fs respectively, which fixes the value of τ c .

Ultrafast electron diffraction
The ultrafast electron diffraction (UED) data for the photodissociation of CS 2 is taken from ref 43 (further details in the SI, Section ??). The data does not support the real-space pair distribution function (PDF) analysis obtained by a sine-transform of the modified scattering signal ∆sM that is common in UED (see SI, Section ??). Instead, the signal is evaluated in the percent difference form as per eq 17.
The experimental signal %∆I exp (s, t) in Figure 3 displays a strong enhancement band (red) next to a strong depletion band (blue) in the range 3.5 < s < 6.0Å −1 . Also note less intense enhancement/depletion features for s < 2Å −1 which appear at later times and which correlate with the onset of strong dissociation.
Given the noise and limited temporal resolution in the data, we carry out a two-step optimisation. In the first step, we fit t 0 using eq 20 with bounds [s min , s max ] = [2.8, 4.2] A −1 . This procedure is repeated in the range t 0 ∈ [− 16,83] fs and τ c ∈ [150, 250] fs. The best global parameters are then used in the determination of the weights w. Initial conditions for w are generated both using unbiased Monte Carlo sampling and using the biased iterative sampling procedure described in SI, Section. ??. In either case, the confidence matrix p conf is based on estimated experimental standard deviations (see SI, eq ??).

CHD ring-opening (UXS)
For the CHD reaction, the data is of sufficient quality that a global one-step optimisation is feasible. Optimisation yields the best fit parameters [t 0 , x] = [−38 fs, 1.3], with a RAE of 0.775. Recall, the scaling factor x is optimised as to uniformly scale the q dependent excitation fraction such that, γ x (q) = γ (q) x. The values of γ (q) and the resulting optimised γ x (q) over the available q range, can be seen in SI Figure ??. The convergence with respect to x and the t 0 shift in Figure 4 shows that more negative t 0 shifts and larger x result in lower values of the target function F (w, c) and better RAE's (with the trend more pronounced for higher x). Generally, F (w, c) is more sensitive to x than t 0 .
The change in the fraction of ring-open trajectories is small as t 0 changes, which is unsurprising given that the majority of the ringopening occurs in a concerted fashion within the first 140 fs. 59 Current literature values are found in the range 40-60%. 59 We note that higher values of x redistribute some of the ring-open weights to ring-closed weights. Since one does not physically expect the fraction of open to closed trajectories to change as a function of the excitation fraction, this is an artefact which stems from the fact that the percent difference signal is stronger for ring-open molecules than ring-closed, due to the reference equilibrium CHD structure being ringclosed. In principle, this could be addressed by a different choice of reference structure.  Table 1. These account for almost all the weights. The trajectories can be visualised in terms of their characteristic C 1 −C 6 bond distance as shown in Figure 5. Three main families of trajectories are observed: direct ring-opening, a slower indirect ring-opening which undergoes several C 1 −C 6 stretches before breaking the bond, and finally ring-closed paths with initially strong oscillations in the C 1 −C 6 bond which are damped out as the energy disperses across all motions. In total, the ring-opening and closed trajectories have a weight of 52 % and 48 % respectively.
The ground state HT has several cis-/trans-(Z) / (E) isomers. Out of the 48% HT product, we predominately observe the presence of cEc-HT (8.7%) and cZc-HT (34.5%) isomers. Due to the length of the 275 fs temporal window used in the fit it is no surprise we can not clearly detect the tZt-HT isomer. However, we note that we observe a small fraction (3.4%) with a configuration between cZt and tZt. With data available over a longer temporal range, one could better refine the ground state dynamics and eventually observe as the system settles into a thermal equilibrium between different HT isomers. In Figure 6, one can see there is little variation from the best set of weights if one examines solutions with more negative t 0 shifts or larger x. The cluster of very good optimisations are all found within a small radius of D 2 b from the optimal solution, and have RAE's , which is a consequence of the optimisation redistributing weight between ringclosed trajectories and selecting a different set of ring-opening trajectories. In addition, we note that the better optimisations (as measured by their RAE) tend to exhibit a larger variance in their set of weights, this can be seen in SI Figure ??.

CS 2 photodissociation (UED)
We now turn to the more difficult case of CS 2 , where the resolution forces us to employ the two-step optimisation procedure. The first step of optimisation yields the optimal global parameters [t 0 , τ c , γ] = −83 fs, 230 fs, 3% .   This emphasises the importance of a robust method for fitting, inverting, and interpreting the experimental data.
For the second step of the optimisation, we take the global parameters determined in step one as the starting point for the unconstrained determination of the TBF weights w. We also repeat the procedure for three different values of t 0 to check that the optimisations are consistent. During optimisation, the pre-determined value of the global parameter γ is allowed to readjust. The value should change little but if it does change significantly this may indicate that the values from step one are suboptimal.
The weighting of the diffraction signal with the confidence matrix, p conf s i , t j , is  crucial to achieve a sensible fit, with significantly poorer results obtained without it. To investigate this effect, the optimisation was repeated for several values of the confidence matrix cut-off threshold, p min conf ∈ [0, 0.45, 0.50, 0.55, 0.60, 0.65]. As seen in SI Figure ??, increasing the threshold corresponds to hard filtering of the data and lower thresh-olds correspond to including more data. For example, p min conf = 0 includes all data, while p min conf = 0.65 only includes the main feature in Figure 3. Examining the RAE in Figure 9d, it is evident that at all values of t 0 , high p min conf is detrimental to the optimisation and p min conf = 0 yields the best results. This suggests that including all data is beneficial as long as the data is given appropriate confidence weights. This is further confirmed by examining the physical observables resulting from the inversion, as discussed next.
Figures 9a-c show the convergence of three physical observables as a function of the convergence threshold p min conf for t 0 = [17, −33, −83] fs. The observables are not explicitly part of the fit, but are indirectly a function of the optimised weights. The branching ratio between singlet and triplet dissociation product varies significantly over the range of times and confidence thresholds, and convergence to the spectroscopy-derived literature value (≈ 3) is only obtained for t 0 = 83 fs and then only at low values of confidence threshold p min conf → 0. In contrast, the convergence of the fraction of bound trajectories at 1 ps is more uniform across the different t 0 values. Values in the range of 20-30%, congruent with spectroscopic results, are observed for the lower confidence thresholds. A similar trend can be observed for the excitation fraction γ, which is comparatively stable across all values of t 0 and p min conf . This indicates that a suitable choice of τ c was made and is also reflective of the fact that the value of γ is most governed by the strongest feature in the experimental data (which is also recorded with the lowest noise).
Overall, the best optimisation yields the final values t 0 , τ c , γ, p min conf = −83 fs, 230 fs, 3.4%, 0 , with the convergence of F (w, c) shown in Figure ?? in the SI for reference. The initial guess of γ = 3.0% from step one is reoptimised to 3.4% in the second step, with the comparatively small change suggestive of a stable solution. Examination of the optimised weights show that eight dominant trajectories account for 99.9% of the total weight. Two of these are bound trajectories that make up a total of 28%, three correspond to singlet dissociation,   Figure ??. From this we attribute the onset of the peak below 2Å −1 in the experimental signal %∆I exp (s, t) to dissociation. In Figure 10 we see that the set of very good solutions, as judged by their RAE value, clus- ter around the best solution as measured by the distance D b (eq 24). Broadly, this indicates that the optimisation procedure succeeds in locating robust global optima. The clustering trend is especially evident when p min conf ≤ 0.5 and t 0 = [−83, −33]. As the value of p min conf increases and t 0 is shifted towards earlier times, the distances increase. These deviations at high values of p min conf are explained by the fact that the bound population at 1 ps is overestimated due to the omission of a distinct feature just below 2Å −1 with a delay in the onset of the characteristic of dissociation (see Figure 3). As with CHD, we see that the better optimisations (as measured by their RAE) also exhibit a larger variance in their set of weights, this can be seen in SI Figure ??.
It is worth noting that for the best t 0 and for each confidence threshold, the optimisation converges to some combination of the same subset of 16 trajectories. While only 8 of these are significant in the best optimisation where (p min conf ) = 0, generally the weights of these trajectories show little variation as a function of the confidence threshold. Some are re-weighted significantly, for example the dominant triplet trajectory with a weight of nearly 45% only appears in the set of trajectories once the confidence threshold is lowered and more data included. The same is true of a singlet trajectory which sees its weight increase 0% → 4.5% as the confidence threshold is decreased. The opposite effect is seen for a triplet trajectory and a bound trajectory. As more data is included, these trajectories see their weight decrease 20% → 7.5% and 32% → 18%, respectively. Thus the contribution of some trajectories can be over-estimated and the contribution of others under-estimated when the optimisation is based on high confidence thresholds, i.e. smaller data sets. A closer examination of the distribution of trajectory weights as a function of t 0 shows that the relative distribution between singlet and triplet trajectories is sensitive to t 0 . This is not surprising given that singlet and triplet dissociation is separated temporally, with the singlet dissociation appearing earlier. Shifts in the temporal alignment of experiment and model thus forces the optimisation to try to compensate by changing the relative composition of the singlet vs triplet dissociation.
Finally, using the iterative weight sampling procedure outlined in the SI, we can confirm that convergence to the global minimum is achieved for the parameters [t 0 , τ c , γ, p min conf ] = −83 fs, 230 fs, 3.4%, 0 (for further details, see the SI).

Conclusions
In this paper, we have presented a forward optimisation method for the inversion of timeresolved data and evaluated its performance on data from actual experiments. The method matches experimental data to a model molecular wavefunction by optimising the weights of trajectory basis functions (TBFs). The TBFs provide appropriate constraints on molecular systems far from equilibrium and ensure that continuity relations are fulfilled. The method can be applied to any type of experiment for which observables can be calculated from the model wavefunction and to any molecular system for which a basis of trajectories can be generated.
Notably, the TBFs ensure a physically sensible solution even when the inversion is underdetermined by the available data, which is the common situation in ultrafast experiments where the observable often only illuminates one specific aspect of the complex dynamics, or when the data has limited temporal, spatial, or energy resolution (for instance a limited q-range in scattering experiments). Importantly, we demonstrate that our approach is robust for noisy data and show that using a confidence matrix for the experimental data is far superior to excluding noisy data from the inversion. For very noisy data, we use a twostep optimisation which first optimises global parameters and then the model wavefunction.
We apply the method to ultrafast scattering data for the molecules CS 2 and CHD. Good agreement with the experimental data is found, and the model performs well also on additional quality measures introduced to evaluate the model wavefunction and the stability of the solution. The appraisal of the solution includes comparison to spectroscopic data not included for in the inversion, with the optimised model reproducing key physical properties such as the branching ratio of ring-open to ring-closed product molecules for CHD and the ratio of singlet to triplet dissociation products in CS 2 .
The procedure allows key dynamic motifs that contribute to the signal to be identified. For CHD, the method disentangles three pathways involved in the ring-opening and for CS 2 it identifies a key marker of dissociation in the signal.
Ultrafast experiments tend to include theory and simulations as part of their analysis. The current approach fits naturally into this workflow by allowing the discrepancies between theory and experiment to be identified while providing a robust interpretation of the experimental data. The forward optimisation analysis can be performed at small additional cost and provides a bridge between the theory and the experiment, and introduces a systematic strategy for addressing potential shortcomings in simulations.
Looking ahead, forward optimisation should be applied to other observables and, more importantly, to combined data sets (e.g. photoelectron spectroscopy and scattering data). For scattering experiments specifically, higher quality data will potentially allow for the inversion to account for state-specific scattering, 83,87 inelastic effects, 102,103 possibly coherent mixed scattering, 104 and, finally, alignment effects. 105 Sophisticated analysis of multidimensional, low-noise data is a prerequisite for unambiguous identification of subtle yet important effects in photochemical dynamics, including interferences or passage through conical intersections, and inversion methods are expected to play an increasingly central role in ultrafast imaging.
Acknowledgement KA acknowledges an EP-SRC doctoral studentship and a public engagement scholarship from the University of Edinburgh. AK acknowledges funding from the EPSRC (EP/V006819 and EP/V049240) and the Leverhulme Trust (RPG-2020-208). This work was also supported by the Department of Energy, Office of Science, Basic Energy Sciences, under award number DE-SC0020276. Finally, AK acknowledges a Fellowship at the Swedish Collegium for Advanced Studies with financial support from the Erling-Persson Family Foundation and the Knut and Alice Wallenberg Foundation.

Graphical TOC Entry
A new inversion method for data from ultrafast experiments, mapping photochemical reactions.

Limitations of ∆PDF(R,t)
It is common to represent time-dependent electron scattering signals in terms of the modified scattering ∆sM(s,t), defined as, ∆sM(s,t) ≈ s I on (s,t) − I(s) off I at(s) . (S1) Here we use the electron scattering convention of denoting the momentum transfer vector as s.
I on (s,t) and I off (s) refer to the signal at time t and the t 0 signal recorded without a pump pulse.
In the denominator, we find the background that results from the incoherent summation of the scattered intensity from individual atoms, which is denoted I at (s). While the signal represented by defined eq S1 allows one to sine-transform the signal into real space, the division by the theoretically calculated I at (s) also results in an absence of error cancellation in comparison to the signal represented as %∆I exp (s,t). See Figure S1 for lineouts of the ∆sM signal. Thus, we note that for data with poor signal-to-noise ratio, it is often beneficial to perform the inversion procedure in the reciprocal space of %∆I exp (s,t).
The sine-transform that allows for the direct inversion of ∆sM(s,t) to real space can be defined as, The resulting pair distribution function, ∆PDF(R,t), reveals the change in the distribution of interatomic distances R with time. Ideally, the limits of integration in eq S2 should extend from s min = 0 to s max = ∞. In reality, for current MeV-UED experiments, one can expect an upper limit s max ≈ 12Å. Moreover, the reliability of the higher range of s values is limited due to the s −4 , 1 hence a dampening factor of e −αs 2 is included to limit the effect of large s noise on the transform. Care must be taken in the tuning of the value of α, there is a balance between limiting the detrimental effect of noise at high s and ensuring that the range of s included is sufficiently large to accurately reproduce the distribution of interatomic distances. At the same time, the signal is usually not recorded in the very low range of s due to centre-shifting and streaking resulting from 2 plasma effects in the target. In the event that only a small amount of data is missing from the low range, linear extrapolation to s = 0 can be attempted. However, when s min is too large, this can result in inaccurate or erroneous PDFs. In brief, the lack of data at low s amounts to the exclusion of low frequencies in the sine transform, resulting in the PDF being artificially shifted by a near constant term. The experimentally determined PDF will then exhibit interatomic distances below that of the acceptable physical limit, as well as an unexpected increase in intermediate distances at later times. If the size of the transform window is truncated to exclude the region below 1Å, then the theoretical PDF can be made to match the experimentally determined PDF, as in Figure S2.
While one could envision performing our procedure directly in real space for very high quality UED data in the future, due to these reasons, we perform the fitting procedure in reciprocal space on the signal represented as %∆I exp (s,t). Overall, this is a much more stable alternative.

Iterative Weight Generation Procedure
In cases where convergence is not easily achieved through the standard Monte-Carlo approach to weight generation on the interval [0, 1], it may be desirable to sample the total weight configuration space iteratively through a series of configuration space samples increasing in size. This involves defining a series of bounding limits (L bound ,U bound ) on the allowed values of the weights. One starts the procedure with a set of bounds that highly constrain the sample space so that the generated initial weights are narrowly distributed around that of the equally weighted case, i.e. where w = 1/N TBF . The target function is then minimised with respect to each of the N init initial weights in the pool. If required, the number of sets of initial weights N init can be selectively tuned so that convergence of the target function to the minima within the allowed bounds is observed. From this initial sample of the weight space, the sets of initial conditions that give rise to the lowest value of the target function are selected and entered into the subsequent pool of initial conditions that are generated by slightly increasing the bounds (L bound ,U bound ) on the allowed values of the weights.
This weight sampling procedure is then repeated for a series of bounds, each of which correspond to relieving some of the constraint on the allowed values and hence amount to increasing the size of the weight configuration space which is sampled. At each step, the best set of initial conditions are selected and inserted into the next iterations pool. One must be careful that the number of initial conditions N init at each iteration is sufficient to allow convergence. Generally, N should increase as the the bounds are lifted as the size of the weight sampling space is increasing. This is repeated until the bounds are such that the weights are sampled on the interval [0, 1], which corresponds to the whole weight configuration space. By projecting through a series of different sized weight sampling spaces we provide a way of ensuring the sampling of the space is sufficient enough for convergence to correct global minima to be achieved. It also allows for the refinement of the number of sample sets N init so that it is large enough for global convergence without having to blindly select a huge value which may result in a large computational overhead.
where β (t, γ) is given by, Here Q off and Q γ (t) refer to the integrated signal on the detector with the laser off and on respectively. Notice the excitation fraction γ (q) now has a dependence on q, this can be seen in Figure S14. The purpose of the explicit dependence of γ (q) on q is to account for the non-uniform nature of the scattered intensity across the measured q range, this arises from the long and extended interaction region. Note, we do not optimise γ (q) directly with respect to q. Rather, we optimise a scaling factor x that acts to uniformly scale γ (q) across the available q range such that it yields a scaled optimised excitation fraction γ x (q), which is related to the initial estimate of the q dependent excitation fraction γ (q) as, γ x (q) = γ (q) x. The initial guess of γ (q) is calculated from Beer-Lambert law. Clearly Q γ (t) is excitation fraction and q dependant. The laser-on case is defined as, where F instr (q) is the instrument response function that accounts for variation in detector sensitivity with q. The limits of integration q 0 and q 1 refer to the radial limits on the detector. The laser off intensity is given as, In the event that the intensity is equal for both the on and off pulses, β reduces to zero and eq S3 reduces to the generalised expression in the main manuscript. As outlined in the main text, the percent difference signal is calculated from a series of trajectories assuming the IAM and using the x-ray scattering form factors f x A (q). 3 In the case of UXS, the beam is linearly polarised, and so we scale %∆I mod (q,t, γ) by a polarisation factor 1 2 1 + cos 2 θ . Next, the signal %∆I mod (q,t, γ) is convoluted to mimic the temporal resolution observed in experiment, as described in the main manuscript. The duration of the UV pump laser and the x-ray pulse were measured as 60 fs and 30 fs respectively, thus the value of τ c is fixed to the product of two Gaussian's with the corresponding width as to mimic the instrument response function. Finally, the signal is binned along the temporal axis into bins of size ∆t = 25 fs. Accounting for all of this, the length of the subset of the signal used in the optimisation becomes t + 275 fs.
To limit the effect of data points with lower statistics on the minimisation of F (w, c), we define the confidence matrix p conf based on the number of photon hits per frame as, where N(t j ) and N max are the number of hits per frame and the maximum number of hits over all frames respectively. The relative accuracy of each point in q is given by p q (q j ).

Ultrafast electron diffraction (UED)
Data was recorded over the range of momentum transfer 0.76 ≤ s ≤ 12.25Å −1 and over a temporal range of −1.26 ≤ t e ≤ 1.73 ps. 4 An initial value of t = −120 fs was selected. In the CS 2 UED case, we define the percentage difference signal in a somewhat simpler fashion as, We use relativistic electron scattering form factors f e A (s) calculated using ELSEPA. 5 Note that due to instabilities in the ab-initio electronic structure calculations at very large nuclear separations, some dissociative trajectories crash before they reach t = 1 ps. These trajectorties are extended beyond the point of crashing using a harmonic model of the vibrational motion of the CS fragment remaining. The dissociated atom is extended linearly according to its last known velocity. We define dissociation as the point at which a trajectory reaches a C−S bond length greater than 3.4Å, which is the point of no return in the simulations.
Due to a lack of statistics on the electron beam, we define a confidence matrix p conf as heuristically based on estimated experimental standard deviations from a bootstrapping procedure 4 , σ (q i , q j ). The elements of the confidence matrix are normalised as follows, 6 The matrices p conf (q i ,t j ) that result from setting p conf (q i ,t j ) ≤ p min conf to zero, where p min conf ∈ [0, 0.45, 0.50, 0.55, 0.60, can be visualised in Figure S9. The case where p min conf = 0 amounts to including all the data points, but weighting them according to the values in p conf (q i ,t j ). The higher the value of p min conf , the more data points are excluded from the optimisation. The largest value of p min conf = 0.65 amounts to including only a narrow band of the main enhancement seen in the signal.

Initial Conditions
Given the Monte-Carlo nature of the sampling of initial trajectory weights w, the choice in the number of initial conditions N init is crucial. Too small a number of N init risks that the sampling density of the weight configuration space is inadequate for ensuring that the global minimum is found. In Figure S7, we see that the target function converges rapidly as a function of the number of initial conditions N init used for the optimisations. We observe the lowest value of F (w, c) when N init = 250. However, from N init = 15 onwards, we see that F (w, c) only changes by 1e−10 between each successive increase in N init . For N init = 100, which is the value of N init we use in unconstrained optimisations, we see that 46 % of the initial conditions within the sample pool converge to a solution within 1e−3 of the best optimisation.
For the best optimisation, the convergence of the value of the target function and the step size at each iteration in the minimisation, can be seen in Figure S8. The value of F (w, c) decreases rapidly within the first few iterations as the step size increases. By iteration 19, F (w, c) reaches a value of 12.18, before converging to 12.12 as the size of the step decreases and the global minimum is located.

Iterative Weight Sampling -CS 2
Due to the flatness of the target function and the presence of many local minima, we determine that the best optimisation (t 0 , τ c , γ, p min conf ) = (83 fs, 230 fs, 3.4%, 0) is the true global minimum through the weight generation procedure outlined in Section 2. We define a set of 35 upper bounding limits U bound , where U bound ∈ [0.075, 1]. Note that the set of limits is non-linear and the corresponding values can be seen in Figure S10. For each U bound we perform a constrained optimisation that starts by generating N init = 100 initial values of w. At each successive iteration we take the best three sets of weights from the previous iteration and add them to the pool of new initial conditions. Generally, the same set of trajectories are picked out in all the optimisations, although the lower the value of U bound , the more trajectories are given an appreciable weight. However, these trajectories are very narrowly distributed. This is apparent in the lower panel of Figure S11, where one can see how small the variance v 2 is for the set of weights that result in the lowest values of F (w, c) corresponding to local minima in close proximity to the global minimum. This can be explained as by forcing a broader distribution in the number of trajectories selected through placing significant constraint on the optimisation, the trajectories are all given a weight close to the limit of U bound . As the constraint is lifted, the variance within the weights w increases. This effect can also be seen in Figure S13, which shows the trajectory weights as a function of U bound . Once the value of U bound = 0.45 is reached, the weights converge and the optimisation distills out the key dynamics from the broad distribution of trajectories. This is also reflected in the upper panel of Figure S11, which shows the variance of the weights for each U bound from the best optimisation where (t 0 , τ c , γ, p min conf ) = (83 fs, 230 fs, 3.4%, 0). Figure S12 shows the convergence of key parameters as a function of U bound . The lack of sensitivity to the value of γ is apparent, indicating that the initial guess of γ as determined in the initial t 0 fit is reasonable. The convergence of the bound fraction appears to take a bit longer, although values remain in the range of 20-30% for the most part. In comparison, the convergence of the branching ratio is significantly slower, the predicted ratio starting off around 1:2 which is in reasonable agreement with literature. However, this quickly gets worse as some constraint on 8 the optimisation are lifted. The effect is due to over-estimation of the contribution of the singlet population, and it is not until there is sufficient freedom that some of this weight is redistribution to several triplet trajectories. The failure of the constrained optimisations to pick out the finer details is not surprising, given the low resolution in the experimental data and that there is some linear dependence between the singlet and triplet contributions to dissociation.

R( 8 A)
" PDF  Figure S2 Lineouts of experimental and theoretical ∆PDF. The experimental ∆PDF exp and the theoretical ∆PDF th (truncated), which results from truncating the transformation range to exclude the low s region, show good agreement. The absence of the low frequencies from small s yields a vertical shift in the ∆PDFs and an unphysical shoulder at low atomic separations. As a consequence, these two curves indicate the presence of a very short and an intermediate (≈ 2.5Å) bond length that have no physical basis in CS 2 , especially not at large times t. In contrast, the ∆PDF calculated from the theoretical ∆sM which includes the low s data does not suffer these artifacts.

Figure S3
The optimised model CS 2 signal %∆I mod (s,t) decomposed into its contributions from a) dissociative singlet, b) dissociative triplet, and c) bound trajectories. Dashed black lines represent the dissociation times of the individual trajectories from which the signal is composed. Note the presence of the peak just below 2Å −1 that is delayed in onset and the shift in the main enhancement, these are both correlated with dissociation. Figure S4 The value of F (w, c) with respect to confidence threshold p min conf at three t 0 shifts. Note that comparisons between different values of p min conf hold no meaning as thresholding data points will always lead to a lower value of F (w, c). The general trend that t 0 = −83 fs gives the lowest value is true across all values of p min conf . This includes p min conf = 0 which is selected as optimal when the convergence of the statistical measures and physical constants are taken into account, as discussed in the main manuscript.

Figure S5
The variance v 2 w of the best set of weights for a series of CS 2 optimisations for which p conf min and t 0 are varied. For each value of p min conf , the varying opacity represents a different t 0 shift, ranging from -83, -33 to +17 fs from dark to light respectively. The variance of the weights increases for the better optimisations with a lower relative absolute error (RAE).

Figure S6
The variance v 2 w of the best set of weights for a series of CHD optimisations for which the scaling factor x and t 0 are varied. For each value of x, the varying opacity represents a different t 0 ranging from -38 to -14 fs (the darker the shade, the earlier the t 0 shift). The variance of the weights tends to increase for the better optimisations with a lower relative absolute error (RAE).

Figure S7
Convergence with respect to the number of initial conditions N init in the Monte-Carlo sampling procedure for the CS 2 reaction. The figure shows the convergence of the minimum value of the target function F (w, c) (top panel) and also the difference between F (w, c) at each successive sampling density governed by N init (bottom panel). Both are plotted on a logarithmic scale.

Figure S8
The convergence of the best selected optimisation for CS 2 , with respect to the number of iterations. The top panel shows the convergence of F (w, c), and the bottom panel shows the convergence of the step size taken along the target function surface at each iteration in the optimisation.

Figure S9
The confidence matrices p conf (q i ,t j ) for a series of thresholds p min conf used in the optimisation process of CS 2 . Each point below the chosen value of p min conf corresponds to thresholding out the data by setting to zero.   Highly constrained optimisations result in the weight being distributed more evenly between a broader range of trajectories. As constraint is lifted, the archetypical trajectories are filtered out.

Figure S14
The q dependent excitation fraction γ (q) used in the CHD optimisation. This mitigates the effect of different contributions of γ (q) across q due to the large interaction region. This is scaled uniformly across q in the target function.