Ultrafast X-ray scattering from molecules

We present a theoretical framework for the analysis of ultrafast X-ray scattering experiments using nonadiabatic quantum molecular dynamics simulations of photochemical dynamics in molecules. A detailed simulation of a pump-probe experiment in ethylene is used to examine the sensitivity of the scattering signal to simulation parameters. The results are robust with respect to the number of wavepackets included in the total expansion of the molecular wavefunction, but the degree of delocalization of the nuclear wavepackets is shown to put a fundamental limit on the resolution attainable. Overall, the calculated scattering signals correlate closely with the dynamics of the molecule.

Scattering-based techniques are complementary to spectroscopic techniques because they probe structure and structural dynamics in a more direct manner. 17 The first successful time-resolved X-ray scattering experiments were performed at synchrotrons using mechanical choppers to create picosecond X-ray pulses. [17][18][19] XFELs are capable of delivering significantly shorter, even less than 30 fs duration, pulses of coherent X-rays with an intensity orders of magnitude higher than that of synchrotrons. [5][6][7] Importantly, the intensity is sufficient to compensate for small X-ray scattering cross-sections and thus enable time-resolved gas-phase X-ray scattering experiments, 7,20 allowing direct comparison with results from other gasphase ultrafast techniques such as time-resolved photoelectron spectroscopy. 21 Theoretical and computational analysis plays an important role in the interpretation of ultrafast experiments. However, although reactive scattering calculations are quite sophisticated 22 and complex processes such as ionization, 23 dissociation, 24,25 and the competition between the two 26,27 can be modeled accurately in very small molecules, dynamics in even slightly larger molecules remains difficult. Ultimately, the reason for this can be traced to the nonlocal nature of quantum mechanics and the exponential scaling with the number of degrees of freedom. 15 The most promising methods 15,16,28 involve solving the time-dependent Schrödinger equation with potential energy surfaces and nonadiabatic couplings calculated onthe-fly. Such calculations compare directly with the observations in time-resolved experiments, and are made feasible by the fairly localized nature of the nuclear wavefunction at short times, reducing the size of the phase-space. 29 In contrast to classical surface-hopping, 30 these methods attempt to capture the full quantum propagation of the nonadiabatic nuclear wavefunction. Similar methods have been used to simulate ultrafast electron dynamics in intense laser fields. 31,32 Progress in the understanding of photochemistry will require close integration of theory and experiments. Accurate quantum molecular dynamics simulations that directly predict the experimental signal can provide a detailed interpretation of experimental data and allow critical scrutiny of the reaction mechanism. In this paper, which presents a simulation of an ultrafast X-ray scattering pumpprobe experiment, we examine the fundamental equations describing X-ray scattering in the context of the ab-initio multiconfigurational Ehrenfest (AI-MCE) method 15 for quantum molecular dynamics. Related methods such as ab-initio multiple spawning (AIMS) 28 and variational multiconfigurational Gaussians (v-MCG) 16 share many characteristics with AI-MCE, and derivations for these methods will follow a similar path. Importantly, we also derive simplified forms of the Xray scattering expressions, suitable for the situation when elastic scattering is dominant, and extend the elastic treatment to include the independent atom model with and without rotational averaging. Computational results are shown for the case of elastic ultrafast X-ray scattering from photoexcited ethylene, and we investigate the sensitivity of the X-ray scattering signal to the specifics of the quantum simulation and the approximations made.

Quantum Molecular Dynamics
Our starting point is quantum molecular dynamics simulations of photoexcited dynamics. We base our representation of the nonstationary wavefunction of the photoexcited molecule on the AI-MCE method, 15,33,34 which is closely related to similar methods such as AIMS 28 and v-MCG. 16 Collectively, these nonadiabatic trajectory-based meth-ods amount to the most successful attempt to date to treat complex photochemical dynamics in comparatively large molecules. It should also be noted that the results obtained here could be adapted to surface-hopping trajectory simulations 30,35,36 in a straightforward manner.
In AI-MCE, the molecular wavefunction at time t is expanded as a sum of N Ehrenfest wavepackets with complex coefficients D k (t), The wavepackets |ψ k (t) , also referred to as Ehrenfest configurations, form a basis set for the calculation. Each consists of electronic and nuclear components, where the nuclear Gaussian wavepacket |g k (t) is shared across N s electronic states, |φ i k . The electronic wavefunctions are multiplied by a complex amplitude, a i k (t), such that |a i k (t)| 2 corresponds to the population on that electronic state. This ansatz has shortcomings, but if a sufficient number of wavepackets are included in the expansion in Eq.
(1), good convergence is attained. 34,37,38 Convergence can be improved by inclusion of a numerical mechanism for branching of the wavepackets when the constituent electronic states have very different gradients. 37 The time evolution of the total molecular wavefunction, |Ψ(t) in Eq. (1), is given by the propagation of each individual wavepacket, |ψ k (t) , corresponding to a semiclassical trajectory given by phase-space coordinatesQ k (t),P k (t) (see below), together with the time dependence of the Ehrenfest coefficients, a i k (t), and the Ehrenfest configuration amplitudes, D k (t). The amplitudes D k (t) are due to the coupling between different trajectories. The propagation uses electronic potential energies and nonadiabatic couplings obtained from ab initio electronic structure calculations. Readers interested in further details are directed to Refs. 15,33,34 For our present purposes, the specific form of the wavepacket in Eq. (2) is of interest. The nuclear wavepacket |g k = |g k (Q k (t),P k (t)) = |g k 1 (Q k 1 (t), P k 1 (t)) . . . |g k N at (Q k N at (t), P k N at (t)) is a product of three-dimensional Gaussian wavepackets (coherent states) for each atom where N at is the number of atoms. The labels Q k α (t) and P k α (t) are the coordinates and the momenta of a trajectoryguided 3D Gaussian coherent state (CS) with the width γ α which is chosen individually for each atom according to the prescription developed in Ref. 39 The component for the atom with index α in coordinate representation is then where R α = (R αx , R αy , R αz ) are the Cartesian nuclear coordinates in three dimensions. We use a bar inR,Q k (t), andP k (t) to refer to the coordinates of all atoms. The probability distribution of a single nuclear wavepacket is thus given by, It is trivial to confirm that a wavepacket is normalized such that 1 = |v k (R)| 2 dR, while the overlap of two nuclear wavepackets is, wherez k = γ 1/2Q k + ıγ −1/2 −1P k , with analogous definition ofz l , allows for a compact representation of phase-space coordinates. 40 To describe electronically nonadiabatic effects the electronic wavefunction in coordinate representation is commonly written as dependent parametrically on nuclear coordinates so that Eq. (2) becomes, where the electronic wavefunction φ i (r;R) depends directly on the electronic coordinatesr and parametrically on the positions of the nucleiR. The overlap matrix for Ehrenfest configurations is similar to Eq. (5), but must include the overlap of the electronic components, where we have removed the superfluous subscripts l and k on the electronic states |φ i . The final equality in Eq. (7) makes use of the orthonormality of the electronic states φ i (r;R)|φ j (r;R) r = δ i j .
In the next section, the identity operator for the Ehrenfest basis will be used. The identity operator must account for the nonorthonormality of the CS basis in which the nuclear wavefunction is expressed, and is given by,

X-ray scattering
We now introduce the equations required to calculate X-ray scattering from the nonstationary molecular wavefunction excited by the pump pulse. We take as a starting point the doubledifferential scattering cross-section given by the expression, 41 where dΩ is the scattering angle and dω k 1 is proportional to the energy of the scattered radiation. The theory of X-ray scattering outlined here follows closely the development by K. B. Møller and N. Henriksen. [41][42][43] The scattering operatorL in Eq. (9) is defined as,L = e 2 2m e j e ıqr j , where e and m e are the charge and mass of an electron, and r j is the position of electron j. Formally the sum runs over all charged particles, but since electrons are much lighter than nuclei and the mass appears in the denominator, only the electrons are included. The momentum transfer vector is defined as q = k 0 − k 1 , where k 0 and k 1 are the incoming and outgoing wave vectors for the Xrays. The HamiltonianĤ M in Eq. (9) corresponds to the material system, in our case the photoexcited molecule. The pre-factor α is, with the term P 2 accounting for the polarization of the X-rays, c the speed of light, 0 the permittivity of vacuum, and ω = kc, where k = k.
The electric field of the incoming X-rays in Eq.
, with t p the delay time between pump and probe for a pump pulse centered at t = 0. For simplicity we consider a Gaussian X-ray pulse profile, (t − t p ) = e (t−t p ) 2 /2γ 2 d , which has the full width at half maximum intensity (FWHM) duration d the X-ray pulse intensity profile and C P (δ) = √ (δ)e ıω k 0 δ the normalized X-ray pulse coherence function.
Rewriting the double-differential scattering cross-section in Eq. (9) using these definitions results in, with the scattering from the material system given by W(τ, δ),

X-ray scattering and dynamics
This section brings together the equations for the AI-MCE wavefunction from Section 2.1 and the X-ray scattering from Section 2.2. The scattering signal is given by Eq. (12), with W(τ, δ) the key quantity to calculate for the material system. A formal solution can be obtained by inserting the identity operator from Eq. (8) into the definition of W(τ, δ) in Eq. (13). However, considering the simplifying approximations we wish to make, a more convenient ansatz is found by first rewriting W(τ, δ) using the propagation operator, U(t, t 0 ) = e −ıH M (t−t 0 ) , such that, which translates into, Eq. (15) is given in matrix form with D the column vector containing the wavefunction expansion coefficients, D k (t) in Eq. (1), and withΩ −1 the inverse of the wavepacket overlap matrix from Eq. (7), both of which come from the AI-MCE simulations. The actual scattering matrix elements are contained in the matrix W on the right-hand-side of Eq. (15), and are given by, where |ψ l (τ) and |ψ k (τ) are Ehrenfest wavepackets given by Eq. (2), and the electronic states on the second line are written as |n and |m . Based on the different time-scales for nuclear and electronic motion, one may proceed to simplify these matrix elements, where the first approximation is obtained by inserting an electronic-state identity operator 1 = i |i i|, introducing the electronic scattering matrix elements L nm = n|L|m , and assuming an adiabatic Born-Oppenheimer Hamiltonian such that The second approximation is obtained by assuming that the time-scale for nuclear motion is significantly slower than electronic motion.

The elastic scattering limit
The expression in Eq. (17) can be simplified further if one assumes that the electronic states are well separated, so that the oscillating phase term, exp ıδ(V n − 2V l + V m )/2 , cancels all nondiagonal contributions to the integral, W lk ≈ n g l |a n * l a n k |L nn | 2 |g k .
The diagonal L nn terms that remain correspond to elastic scattering, and can be rewritten using the elastic scattering form factors, f 0 (q; n,R), for each electronic state |n(r;R) , It is straightforward, for instance by partial integration of all but one electron coordinates in the bracket, to show that the form factor is the Fourier transform of the electron density, ρ n (r;R), These molecular form factors can be calculated directly from ab initio electronic wavefunctions. 44 The difference in form factor, or equivalently in electron density, between electronic states is comparatively small. 44 Neglecting these differences, the ground state electron density can be approximated further by a sum over single-atom fragments, where ρ α (r) is the electron density of each atom centered at its position R α . A Fourier transform of the approximate electron density in Eq. (21) yields a sum of atomic form factors, f 0 α (q), multiplied by phase-factors corresponding to the positions of the atoms, 46 where R α is the position of atom α, and q the momentum transfer vector with q = |q|. This approximation, known as the independent atom model (IAM), affords significant computational savings since the atomic form factors are known and tabulated. 45 In many circumstances, not the least in traditional X-ray diffraction from thermal samples, the IAM is a very reasonable approximation, despite that it fails to account for the distortion of the electron density due to chemical bonding and does not distinguish between different electronic states. 44 At this point we should point out the restrictions on the elastic approximation when X-ray scattering at XFELs is considered. The elastic approximation is valid when the duration of the coherent X-ray pulses is short compared to the time-scale for nuclear motion, yet sufficiently long that the scattering cross-terms between different electronic states average out. [41][42][43] When these conditions are not fulfilled, then the full scattering must be taken into account. 47

IAM scattering matrix elements 2.4.1 Ansatz for W IAM
The wide use of the independent atom model and the significant computational efficiencies it offers, motivate us to explore this model further. We can use the IAM approximation to simplify the matrix W in Eq. (16). Starting with the final line of Eq. (17), we set the ground state population to unity, a (1) k = 1, and assume that the ground state electron density is well described by the independent atom model, such that L 11 ≈ (e 2 /2m e ) f IAM , where f IAM is the form factor according to IAM. We thus get with the IAM scattering matrix, W IAM , defined as, In its square amplitude form Eq. (22) becomes, where R αβ = R α − R β is the vector between two atoms, and indices α and β run over all atoms in the molecule. The atomic contribution, in Eq. (25) carries no structural information. The comparatively simple expression that results from the independent atom model makes it possible to find near-analytical forms for the matrix elements of W IAM .

Diagonal elements of W IAM
We first consider the diagonal elements, W IAM kk , of the IAM scattering matrix defined in Eq. (24). These matrix elements correspond to a convolution of | f IAM | 2 from Eq. (25) by the nuclear probability distribution for a wavepacket given by Eq. (4), which yields, The resulting expression is similar to Eq. (25) with the scattering interference between pairs of atoms proportional to their separation, but with the nuclear coordinates (R 1 , . . . , R N at ) replaced by the wavepacket coordinates (Q k 1 , . . . , Q k N at ), and an additional damping term exp (−q 2 /2γ αβ ) proportional to the combined coherent state widths for the two atoms, γ α and γ β (see Eq. (3)). The new factor γ αβ is defined as which reduces to a constant prefactor exp (−q 2 /2γ) when all {γ α } are identical. This damping is proportional to the degree of localization of the nuclear wavepacket. In the limit of strong localization (γ → ∞) the traditional IAM formula is recovered, while, conversely, the scattering interference responsible for structural information is lost when the atoms are strongly delocalized (γ → 0).

Off-diagonal elements of W IAM
Derivation of the off-diagonal matrix elements of W IAM lk (l k) is tedious but straightforward. The re-sult is, where the overlap Ω lk is defined by Eq. (5), and Ω αβ lk = g l α |g k α g l β |g k β is the partial overlap of particles α and β in wavepackets l and k respectively. Furthermore, and with the general definition of difference vectors Q lk αα and P lk αα given by and, finally, Overall, Eq. (29) for the off-diagonal elements shares much the same structure as Eq. (27) for the diagonal elements, including the damping factor exp (−q 2 /2γ αβ ), but is weighted by the overlap between the nuclear wavepackets, via Ω lk , and additional damping terms proportional to the phasespace distance between identical atoms via lk α . Finally, in the off-diagonal matrix elements the interference between pairs of atoms is given by the average of the distance between atoms in each of the two wavepackets, via the imaginary part ofρ.

Rotational averaging of W IAM
Rotational averaging becomes important in the absence of rotational alignment or orientation of the molecules and yields an isotropic signal that only depends on the magnitude of the momentum transfer vector q. In IAM, rotational averaging can be carried out in a straightforward fashion since the atomic form factors, f 0 α (q), are spherically sym-metric. We proceed by integrating over all directions of the vector q. For the IAM form factor in Eq. (25) this yields the well-known formula, 46 | f IAM (q)| 2 = I at +2 N at where R αβ = |R α − R β | is the distance between two atoms. The same rotational averaging for the diagonal matrix elements of W IAM in Eq. (27) yields, where W IAM rot denotes the rotationally averaged matrix, and Q kk αβ = |Q kk αβ | = |Q k α − Q k β | as previously defined in Eq. (32). Again, it is worthwhile noting the similarity to the original IAM formula in Eq. (35).
Averaging the off-diagonal elements of W IAM is more intricate due to the vectorρ having both real and imaginary components, but amounts to the replacement of the "2 cosh (ρq)"-term in Eq. (29) by a function F rot (ρ, q). This gives the rotationally averaged off-diagonal elements as, with F rot (ρ, q) defined as, F rot (ρ, q) = 1 2 π 0 e ρ ri q cos θ J 0 (ϕ(θ)) sin θdθ . . .

BAT approximation
In a recent publication Makhov et al. introduced a braket-averaged Taylor expansion (BAT) as a means of efficient calculation of matrix elements which depend on the nuclear coordinates. 37 The simplest form of BAT calculates potential energy and nonadiabatic coupling matrix elements, which are used to couple the wavepackets in Eq.
(1) and to propagate the equations for the coefficients D k (t), as an average of the matrix elements for individual Ehrenfest trajectories (wavepackets). BAT economises greatly the computational cost of quantum propagation of the wavefunction in Eq. (1) as the number of the matrix elements calculations scales as N, the number of Ehrenfest wavepackets, instead of N 2 as would normally be required without BAT. Also, BAT allows Ehrenfest trajectories to be calculated independently from each other. Running trajectories on-the-fly is the most expensive part of the quantum propagation but can be done in parallel. Thus, although the total CPU time for a basis of hundreds of trajectories can be several months or even years, parallelization makes these calculation tractable. Once the trajectories are stored on disk, BAT makes it possible to propagate the coupled equations for the coefficients D k (t) without additional electronic structure calculations, reusing potential energy and nonadiabatic coupling matrix elements already calculated for the individual trajectories.
The BAT approach is easier to derive with the help of a different form of the Ehrenfest configuration, (39) in which unlike in Eq. (6) the electronic wavefunction is "attached" toQ k (t), the center of the nuclear wavepacket, instead of φ i (r;R) as in Eq. (6). If the electronic wavefunction is smooth and does not change much at the scale of the nuclear wavepacket, the two approaches that rely on Eqs. (6) and (39) are dynamically equivalent. This imposes certain limits on the parameter γ, which describes the width of the wavepacket, namely the Gaussian in Eq. (3) must be sufficiently narrow and γ cannot be too small. BAT has been tested in Ref. 37 and produces results in quantitative agreement with the benchmark. This is not surprising as BAT is expected to fail only in regions where the potential energy and/or nonadiabatic coupling matrix elements change rapidly between coupled trajectories. This may occur near conical intersections, especially if the conical intersection is between the trajectories, but is rare in practice.
BAT also simplifies the calculation of matrix elements required in the current theory. Specifically, in the 0th order BAT approximation the matrix elements of W IAM become, with | f IAM (q;Q l )| 2 defined as in Eq. (25). The removal of theR-dependence significantly simplifies the evaluation of the matrix elements compared to Eq. (24).
Similarly, the corresponding rotationally averaged form, W IAM,BAT rot , is obtained by using the rotationally averaged form of f IAM in Eq. (35), where Q ll αβ = |Q ll αβ | following the definition in Eq. (32). The effect of the nuclear wavepacket when evaluating the matrix elements using the BAT approximation amounts to a simple weighting factor g l |g k since the smearing-out effect of nuclear averaging is no longer included.

Computational details 3.1 Quantum molecular dynamics
The nonadiabatic photodynamics of the ethylene molecule (H 2 C=CH 2 ) upon excitation by a pump laser pulse was simulated using the AI-MCE quantum molecular dynamics method, and the nonstationary molecular wavefunction thus obtained provided the input for the calculation of the timeresolved X-ray scattering.
The simulations of the photoexcited ethylene molecule followed the protocol developed in Ref. 34 The electronic potential energies, their gradients, and the nonadiabatic couplings were calculated on-the-fly at the three-state-averaged complete active space self-consistent field (SA3-CASSCF) level using the MOLPRO quantum chemistry package. 48 A small but balanced CAS(2,2) active space was used, known to describe the lowest two electronic excited states semi-quantitatively, 28,49 with the Dunning's cc-ppVDZ (correlation consistent, polarised valence, double zeta) basis set. 50 The Ehrenfest trajectories were sampled in the Franck-Condon region using a Wigner distribution, 51 with initial population completely localized on the first excited S 1 ππ * state.

X-ray Scattering
The time-resolved X-ray scattering signal for the photoexcited ethylene was calculated using the IAM approximation. Only elastic scattering was considered, with the assumption that ω k 1 = ω k 0 (k 0 =k 1 ) and C P (δ) → δ(0), where δ(0) is the Dirac delta function. In this situation, the scattering signal in Eq. (12) from the excited molecular wavefunction, |Ψ(τ) , reduces to where W(τ, q) is taken from Eq. (15) withΩ ≈ Ω and W ≈ (e 2 /2m e ) 2 W IAM (or W IAM rot for rotational averaging). When calculating the scattering matrix W IAM , it is numerically convenient to separate out the constant atomic component I at such that, W IAM = Ω • I at + W IAM,mol , where • denotes the element-wise Hadamard matrix product. The time t p corresponds to the delay time between pump and probe, for a pump pulse centered at time τ = 0.
The experimental observations are normally represented by a 'laser on' -'laser off' difference sig-nal in the following form, 20 where γ excit is the fraction of excited molecules, dS (q, t p ) is the signal corresponding to the excited molecular wavefunction, and dS off (q) is the background signal from unpumped molecules. Inelastic scattering corrections commensurate with the present approximations 7,20 have been tabulated using Waller-Hartree theory. 45,52,53 In the calculations, 13.8 keV X-ray photons corresponding to 0 ≤ q ≤ 14 Å −1 were used, with an X-ray pulse duration of 25 fs and an excitation fraction for the molecules of γ excit = 9%. The atomic form factors f 0 α (q) were taken from tables in Ref. 45

Quantum molecular dynamics
The AI-MCE simulations reveal that following photoexcitation into the Franck-Condon region of the ππ * electronic S 1 state, the ethylene molecule undergoes cis-trans isomerization around the C 1 =C 2 double bond. The molecule then decays via two competing processes. The first is nonradiative decay through a twisted or pyramidalized conical intersection, and the second is H-atom migration to form ethylidene (CH 3 CH), which then decays through a different conical intersection. In terms of populations, the population of S 1 changes slowly for the first 30 fs after the photoexcitation. At this point, the excited molecular wavefunction reaches a region where the gap between S 1 and S 0 is sufficiently small to allow for efficient population transfer, and the population then decays exponentially with an approximate lifetime of τ ≈ 112 fs.
The time evolution of the nuclear wavefunction is shown in Fig. 1 in terms of a probability density contour plot for the twist and pyramidalization angles, which together with the C 1 -C 2 distance discussed below can be used to characterize the photodynamics of ethylene. The two angles are defined mathematically in the caption of Fig. 1. The twist angle corresponds to a twist around the C 1 =C 2 double bond, and the pyramidalization angle reflects the degree of deviation of the two carbons from sp 2 hybridization. The nuclear wavefunction initially moves ballistically along the twisting coordinate with a period of ∼40 fs (the twisting of ethylene on the electronic ground state S 0 has a period of ∼33 fs). During the first twisting cycle, the wavefunction is almost totally located on the S 1 state but at later times the population starts transferring to the S 0 state. This picture is consistent with the results of Ref. 28 obtained using the AIMS method at the same level of electronic structure theory. From 50 fs and onwards, the wavefunction becomes quite dispersed, as evidenced by the bond angle distributions in Fig. 1 and the bond-length distributions in Fig. 2.
A second characteristic of the ethylene photodynamics is an early oscillation in the C 1 -C 2 bond distance as shown in Fig. 2. The graph shows the median C 1 -C 2 bond distance, as well as the first and third quartile (dashed lines) which enclose the bond distances for 50% of all the trajectories. At shorter times there is a coherent oscillation in the C 1 -C 2 distance, which leaves a distinct signature in the X-ray scattering pattern as we will see in the next section.

X-ray scattering 4.2.1 Overall characteristics
The simulated dynamics of the photoexcited ethylene molecule was used to calculate elastic X-ray scattering following the procedure outlined in Section 3.2. The scattering signal is shown as the percentage difference signal, ∆dS (q, t p ), defined in Eq. (43). In the following we discuss the properties of this signal, and its dependence on factors such as the size of the nuclear basis in the simulations and the wavepacket width parameters.
The rotationally averaged elastic X-ray scattering difference signal calculated from the full simulation is shown in Fig. 3a as a function of time. The signal is dominated by the two carbon atoms, which contribute 12 out of 16 electrons in ethylene, with comparatively minor contributions from the four hydrogen atoms. The variation in the signal in Fig. 3a becomes progressively smaller with time as the increasing dispersion and associated   delocalization of the nuclear wavefunction averages out specific motions. This increase in nuclear dispersion with time was already evident in Figs. 1 and 2. At times t < 50 fs, however, the coherent stretch of the C 1 -C 2 bond, clearly visible in Fig. 2, results in a strong signature in the scattering signal across the whole range of the momentum transfer q.
The sensitivity of the calculated scattering signal to the nuclear basis used in the simulations, i.e. the number of Ehrenfest trajectories included when calculating ∆dS (q, t p ), can be seen when comparing Fig. 3a, which shows the signal calculated for 1000 trajectories, and Fig. 3b, which shows the signal for a randomly selected subset of 20 trajectories. Although the exact appearance of the signal from the smaller set will depend on which 20 trajectories are included, a subset of 20 is sufficiently large to make these differences rather small. Interestingly, even the small subset correctly depicts the main features of the scattering signal. This is not entirely surprising, since in a recent analysis of ultrafast X-ray experiments, 7 we demonstrated that a judicious choice of trajectories allows an accurate representation of the experimental signal with only a small number of trajectories. However, the smaller subset does underestimate the dispersion of the nuclear wavepacket at longer times, as can be seen from the features present at long times in the small set in Fig. 3b, but absent from the full set in Fig. 3a.
We end this section with a comparison between the rotationally averaged signal and the signal from perfectly aligned molecules, shown in  Figure 4: Comparison of the elastic scattering difference signal, ∆dS (q, t p ), at time t p = 25 fs for rotationally averaged (4a) and perfectly aligned (4b) molecules. The signal is calculated for 200 trajectories without convolution, and the contour plots are drawn to correspond to detector images for the incoming X-ray pulse perpendicular to the C 1 -C 2 molecular axis and the initial plane of the molecule. The maximum radius shown corresponds to momentum transfer q = 10 Å −1 . X-ray polarization is not taken into account.
to fully represent the scattering signal as a function of the amplitude of the momentum transfer q, as done in Fig. 3 for instance. On the other hand, the scattering signal from aligned molecules, shown in Fig. 4b, reveals a great deal more detail. The scattering image is distinctly asymmetric due to the loss of symmetry in the molecule by time t p = 25 fs (see also scattering patterns in Ref. 44 ) but the contour lines over-emphasize somewhat the difference between the left and the right side of the image. Although the fully aligned signal in Fig. 4b represents an idealized scenario that cannot be achieved experimentally, even for partial alignment it is clear that the information content in the scattering signal would increase dramatically, something illustrated by a number of recent time-independent X-ray 54 and electron 55 scattering experiments. Note that the attenuation due to the X-ray polarization is not included, which will mask parts of the image depending on the orientation and character of the polarization. 46   We begin by examining the effect of the width of the individual wavepackets. The widths are determined by the factors γ of each coherent state as shown in Eq. (3). In Fig. 5 we compare the scattering difference signal, ∆dS (q, t p ), calculated for the default values of γ with the same scattering signals re-calculated using modified values of γ. Specifically, the comparison is to γ/3, corresponding to a more delocalized wavepacket, and to 3γ, corresponding to a more localized wavepacket. The results show that the more delocalized the wavepacket is, the weaker the signal becomes for large values of q. Even the comparatively modest reduction to γ/3 is sufficient to eliminate the signal for q > 8 Å −1 . The trends observed in Fig. 5 were anticipated from the expressions derived earlier, see the discussion below Eq. (28). The observed damping is commensurate with Eq. (27) in which the strength of the damping increases with q. As discussed in Section 2.4.5, the ab-initio MCE quantum dynamics approach is robust only if sufficiently narrow Gaussians are used as a basis for the nuclear dynamics. Therefore complete washing out of the signal only occurs in the smallγ regime where the ab-initio MCE method is no longer reliable.

The nuclear wavefunction
Next, we examine the influence of the offdiagonal elements in the W IAM and the W IAM rot matrices on the X-ray scattering signal. The inclu-sion of the off-diagonal elements increases computational overheads significantly by changing the scaling from N to N 2 , where N is the number of Ehrenfest trajectories included in the wavefunction expansion. Fig. 6 compares the calculated signal, ∆dS (q, t p ), at time t p = 25 fs for 100 trajectories, with and without the off-diagonal elements included. The difference is the greatest at large values of q. This is to be expected, since the interference between trajectories, mediated by the off-diagonal elements, is strongest when the distance in phase space between the trajectories is small. The greater the number of trajectories and the more delocalized they are (i.e. small γ), the greater the effect on the signal from including the off-diagonal elements, although the latter factor (delocalization) will be counteracted by a corresponding increase in damping at large q.
An important point is that if the quantum molecular dynamics simulations are converged (and if the off-diagonal elements are included), the dependence of the scattering signal on the width of the individual wavepackets γ will vanish, since the nuclear wavefunction will have the correct shape given by the quantum mechanics of the situation.

BAT approximation
Having examined the effect of the nuclear wavefunction on the X-ray scattering, it is appropriate to compare the calculations so far to the far simpler BAT approximation discussed in Section 2.4.5. The BAT approximation is computationally faster by strongly reducing the dependence of the matrix elements on the nuclear coordinates. To a degree, the BAT approximation emphasises the 'classical' aspect of each trajectory by removing the nuclear wavepacket from the matrix elements. As the comparison in Fig. 7 shows, this changes the scattering signal for large values of q significantly, with a marked difference between the BAT approximation and the conventional calculation for q > 4 Å −1 . This is unsurprising, since for large values of q the spatial resolution in the scattering is sufficient to detect the comparatively fine differences in the nuclear wavefunction as described in the conventional propagation versus the BAT approximation.
On the other hand, for small values of q the BAT Figure 7: Contour plot of the absolute difference between the scattering signals, ∆dS (q, t p ), calculated using the standard method (see Fig. 3a) and the BAT approximation, shown as a function of q in Å −1 and t p in fs. The difference between the two calculations is negligible for q < 4, but becomes significant for q > 4.
approximation agrees well with the signal, which supports previous work using the diagonal BAT approximation to interpret ultrafast X-ray scattering experimental data. 7 In that case, the maximum q range of the experimental data was approximately 4.3 Å −1 , comfortably within the validity range of the approximation.

Conclusions
In this paper, we have constructed a theoretical and computational framework for the analysis of ultrafast X-ray scattering experiments using nonadiabatic quantum molecular dynamics simulations. The fundamental expressions for X-ray scattering have been put in a form suitable for state-of-the-art on-the-fly quantum molecular dynamics methods. Particular attention has been given to the development of a hierarchy of approaches for the pre-diction of elastic X-ray scattering, some of which have already been used for the analysis of recent experiments. 7 At one extreme, the BAT approximation combined with IAM form factors and a comparatively small number of trajectories is sufficient to produce good results for relatively small values of the momentum transfer. This can be improved upon systematically by the inclusion of the full dependence of the wavefunction on nuclear coordinates, and also by the replacement of the IAM form factors by ab-initio molecular form factors. 44,56,57 In recognition of the potential of molecular alignment in scattering experiments, 54 the equations in this paper are given both in a rotationally averaged form and in a general form more appropriate for full or partial alignment. Although the presented theoretical framework allows for the calculation of inelastic effects in the scattering of coherent X-rays, the focus of the present calculations has been on elastic scattering. A detailed simulation of a pump-probe experiment in ethylene has been used to examine the sensitivity of the predicted scattering to the parameters of the simulation. We have found that the results are robust with respect to the number of wavepackets included in the total expansion of the molecular wavefunction. The results are however quite sensitive to the degree of delocalization of the wavepackets. The damping of the scattering signal with q implies that there might be limited rewards for measuring signals for large values of q in time-resolved scattering experiments. On the upside, scattering should provide a sensitive probe of the actual shape and dispersion of the nuclear wavefunction.
Ultrafast X-ray scattering has been long anticipated 58 and is now in a period of rapid development. The first experimental results have appeared in the literature, [5][6][7] and with the ongoing construction of new XFEL facilities we expect intense interactions between experiments and theory over the coming years.