Excited Electronic States in Total Isotropic Scattering From Molecules

Ultrafast x-ray scattering experiments are routinely analyzed in terms of the isotropic scattering component. Here we present an analytical method for calculating total isotropic scattering directly from ab initio two-electron densities of ground and excited electronic states. The method is generalized to compute isotropic elastic, inelastic, and coherent mixed scattering. The computational results focus on the potential for differentiating between electronic states and on the composition of the total scattering in terms of elastic and inelastic scattering. By studying the umbrella motion in the first excited state of ammonia, we show that the associated electron density redistribution leaves a comparably constant fingerprint in the total signal that is similar in magnitude to the contribution from the changes in molecular geometry.


Introduction
The X-ray Free-Electron Laser (XFEL) facilities that have emerged around the world in the last decade provide coherent and ultrashort x-ray pulses, whose peak brightness is more than ten orders of magnitude larger than for synchrotron radiation. With a pulse duration that rivals that of optical lasers, XFELs have greatly enriched the palette of experimental techniques used to study the most fundamental aspects of chemistry -how molecules move, 1 how chemical bonds are made or broken, 2 and how electrons rearrange after in-teraction with light. 3 One such powerful technique is non-resonant ultrafast x-ray scattering from gasphase samples. 1,[4][5][6][7][8][9] In a pump-probe fashion, an ensemble of molecules is pumped by an optical laser to an excited electronic state, and the resulting photodynamics is probed via hard x-rays with varying delay time. Due to the fast nature of the dynamics and the small number of scattering molecules, these experiments require ultrashort pulse duration and large photon numbers currently only provided by XFELs.
It is not always recognized that gas-phase xray scattering has fundamental differences from xray crystallography. On account of the large average separation between molecules and the absence of a regular lattice, the gas-phase scattering is free from intermolecular interferences for all but the smallest values of the scattering vector. 10 The signal should thus be understood as an incoherent sum of scattering intensities from isolated molecules. It is on this single-molecule scale that quantum effects are most easily observed. In contrast, the periodicity of crystals means that the signal at the Bragg peaks is strongly dominated by elastic scattering, 11 which, unlike the total scattering in gas phase, is a one-electron property. 12,13 It is furthermore worth noting that xray crystallography is traditionally concerned with molecules in their thermal ground states, while the laser-induced dynamics in pump-probe experiments evolves on multiple electronic states, each characterized by its own distinctive electron distribution. It follows that the theoretical tools devel-oped for and successfully applied in x-ray crystallography for more than a century are not always best suited for ultrafast x-ray scattering. [14][15][16][17][18][19][20][21][22][23][24] Indeed, gas-phase scattering from ground state molecules initially, and, more recently, ultrafast xray scattering, have prompted the development of a number of algorithms that aim the accurate prediction of x-ray scattering starting directly from the ab initio electronic structure of molecules. As gas-phase samples in thermal equilibrium are isotropic, a central question in these methodologies is that of rotational averaging. Wang and Smith have first suggested a direct analytical method for evaluating isotropic scattering intensities. 25 Various permutations of this approach have been devised since with the ultimate goal to reduce computational effort. [26][27][28] Alternative strategies that relay on numerical rotational averaging 21,23,29 or grid-density methods 24 have also been suggested. Although it might not be immediately obvious, isotropic scattering is directly applicable to gas-phase ultrafast x-ray scattering despite the anisotropy imposed on the sample by the pump laser's polarisation. Separation of the isotropic from the anisotropic part of the experimental scattering signal can be achieved by means of a Legendre decomposition, [30][31][32] with the isotropic component equivalent to the outcome if the ensemble was fully isotropic. 33,34 In this context, one of the challenges that ultrafast x-ray scattering is facing is how to accurately predict scattering intensities from the huge conformational space explored by the molecule during photochemical dynamics. This issue does not only hold the key to accurate simulations of ultrafast scattering but also couples to the problem of inversion of experimental data to molecular geometries via iterative procedures as done in x-ray crystallography. Currently available algorithms are often computationally too demanding for such a high-throughput task. In this article we extend the method developed by Crittenden and Bernard 28 for calculating isotropic scattering as a sum of spherical Bessel functions. We demonstrate the existence of a recursive relationship between the expansion coefficients that allows for a significant speed-up of the calculation, in addition to generalising the approach to an arbitrary angular momentum. We illustrate how the method scales with the level of theory and the basis set used. We show that it is applicable not only to elastic and total scattering but also to inelastic and coherent mixed terms between different electronic states. 22,35 We make use of our methodology to investigate the different components comprising the total scattering by a model multielectron two-state system, namely the ground and first excited states of ammonia along the umbrella normal mode. We demonstrate that in this particular case, the change in the scattering signal upon optical excitation cannot be simply attributed to changes in the molecular conformation only, and is in fact strongly influenced by the electron density redistribution in the excited state due to excitation.

Theory
For ultrafast x-ray scattering from photoexcited molecules, the differential scattering cross section per solid angle Ω) can be shown to take the following form, 36 e ιq q q·(r r r i −r r r j ) |Ψ(t) , where Ψ(t) is the time-dependent wavefunction of the molecule, q q q is the momentum transfer vector, and r i and r j are the position vectors of electrons i and j, respectively. The detected signal is proportional to (dσ/dΩ) Th , which is the differential Thomson scattering cross-section for a free electron which includes the polarization factor |e 1 · e 2 | of the incoming and scattered x-rays), 33,37 and W(∆ω), the window function with a detection window defined by ∆ω. The braket notation implies integration over both the electronic and nuclear coordinates. Eq. (1) is valid in the limit of a large detection window, i.e. in the absence of energy resolution on the detector, as discussed in Ref. 35. Note that, in deriving this equation, the high photon energies of the x-rays compared to the energy spectrum of a typical molecule allow us to apply the Waller-Hartree approximation 38 and disregard the comparatively small changes in the photon en-ergy. Upon photoexcitation the molecular wave function, |Ψ(t) , is described by a Born-Huang expansion in the basis of the N electronic eigenstates ψ I (r r r;R R R) accessed during the dynamics. The electronic eigenstates are functions of the electronic coordinatesr r r of the N e electrons in the molecule and depend parametrically on the molecular frame nuclear coordinatesR R R, In this expansion, each electronic state I is multiplied by the corresponding time-dependent rovibrational nuclear wave packet |χ rv I (t) , which depends on the internal nuclear coordinatesR R R and on the three Euler angles αβγ, which relate the molecular and laboratory frames. Using Eq. (2), the differential scattering cross section becomes, The key quantity in Eq. (3) is I I J (q q q,R R R), which is the two-electron scattering matrix element given by, Since terms with i = j in Eq. (4) reduce to the Kronecker delta δ I J , the two-electron scattering matrix element can be written as, where I I J (q q q,R R R) is the pure two-electron part of I I J (q q q,R R R) with i j. Going further, using the sifting property of the Dirac delta function, exp[ιq q q · (r r r i − r r r j )] can be expressed in an integral form: e ιq q q·(r r r i −r r r j ) = dr r r 1 dr r r 2 e ιq q q·(r r r 1 −r r r 2 ) δ(r r r 1 −r r r i )δ(r r r 2 −r r r j ), so that the integral over the electronic coordinates in Eq. (4) becomes, I I J (q q q,R R R) = N e δ I J + 2 dr r r 1 dr r r 2 ρ (2) I J (r r r 1 , r r r 2 ,R R R)e ιq q q·(r r r 1 −r r r 2 ) , where ρ (2) I J (r r r 1 , r r r 2 ,R R R) is the expectation value of two-electron density operatorρ(r r r 1 , r r r 2 ) = (1/2) N e i N e j i δ(r r r 1 −r r r i )δ(r r r 2 −r r r j ). 39 At that point, it is prudent to differentiate between the diagonal elements with respect to the electronic states, I = J, and off-diagonal (mixed) terms, with I J.
In the former case, ρ (2) II (r r r 1 , r r r 2 ,R R R) gives the probability of finding one of the electrons of the system in state I at r r r 1 , while another electron is at r r r 2 . This term can be further separated into two contributions by expanding the two-electron density as a sum of products over one-electron density functions, ρ (1) IK (r r r,R R R), which follows from insertion of the resolution of the identity in the basis of the electronic states. The contribution to the total scattering from the first term in Eq. (8) is the elastic scattering, while the contribution from the second term is the inelastic scattering. When I J, the quantity ρ (2) I J (r r r 1 , r r r 2 ,R R R) is referred to as the two-electron transition density function (i.e. the diagonal part of the density matrix) in order to differentiate it from the case of I = J, which is simply known as the two-electron density function. The I J scattering terms play a critical role in coherent mixed scattering, which appears when there is a coherence between two electronic states. 35,40,41 This concludes the presentation of the fundamental theory.
The main focus of this article is an efficient methodology to calculate isotropic scattering signals. This is motivated by the central role that the isotropic scattering signal plays in the interpretation of experiments, irrespective of the degree of alignment in the sample. 33,34 To appreciate this, we must consider the standard approach in gasphase scattering experiments to decompose the ob-served signal in the basis of orthogonal Legendre polynomials, P α (cos θ q ), 30 dσ dΩ (q q q) = α P α (cos θ q )S α (q), (9) where cos θ q is the component of the unit scattering vector in direction of the laser polarization axis shown in Fig. 1. The angle θ q is related to the detector angles by cos θ q = sin (θ d /2) cos δ − cos (θ d /2) cos φ d sin δ, where δ is angle between the directions of the laser polarisation and the X-ray beam propagation. In the case of a perpendicular pump-probe arrangement, i.e. the case discussed here, when the laser polarisation axis is perpendicular to the X-ray propagation, this reduces to cos θ q = − cos θ d 2 cos φ d . It should be emphasised that the Legendre polynomials are functions of θ q and not of the detector angles, θ d and φ d . The advantage of this approach is that the contribution from the internal and external molecular degrees of freedom can be separated out, as will be shown below. It is important to point out that in many pump-probe gas phase scattering experiments, the detected signal is rarely fully isotropic on account of the preferential excitation of the molecules whose transition dipole moments align with the polarization axis of the linearly polarized pump laser. 8 A remarkable property of the decomposition in Eq. (9) is that the entire information for the internal degrees of freedom can be extracted from any of the components, S α (q), or a combination of them. 42 While in the high-order terms, the internal dynamics is mingled with information about the rotational wavepacket, the zeroth order term in the Legendre expansion can be analysed solely from the point of view of the internal molecular degrees of freedom. It is exactly equivalent to the scattering from a fully isotropic ensemble and consistent with the Debye scattering formalism (see Supplementary Information for derivation). Using the orthogonality of the Legendre polynomials and that P 0 (cos θ q ) = 1, the zeroth-order term, which is henceforth referred to as isotropic scattering, is Figure 1: Illustration of the geometrical relations in X-ray scattering. The laboratory frame coordinate system XYZ is defined so that the incoming X-ray beam wave vector, k k k 0 , is aligned with the Y-axis, while the direction of the optical pumplaser polarisation points in the Z direction. After interaction with the sample, radiation is scattered in the direction k k k. The detector angles θ d and φ d are defined as the polar and the azimuthal angles of k k k with respect to the Y axis and the ZX plane, respectively. The momentum transfer vector, q q q = k k k 0 − k k k, forms the polar angle θ q with the positive Z axis, and azimuthal the angle φ q upon projection onto the XY plane.
given by, where the average is taken over the angular coordinates of the momentum transfer vector. With reference to the full expression for the differential scattering cross-section in Eq. (3), only the twoelectron scattering matrix element, I I J (q q q,R R R), depends on the scattering vector. Furthermore, the average renders the scattering signal independent of the orientation of the molecule in space, meaning that I I J (q,R R R) = I I J (q q q,R R R) θ q φ q does not depend on the Euler angles (detailed derivation included in Supplementary Information). That allows for a separate integration of the rotational wave packets, resulting in a simple scaling factor, which is equal to one for I = J, where the vibrational, |χ v , and rotational, |Θ r , components of the rovibrational wave packet have been made explicit, and their dependence on time have been omitted for brevity. It follows that the isotropic component of the experimental scattering signal can be analyzed, as if the entire ensemble showed a fully isotropic rotational distribution.

Methods
We now turn our attention to the analytical evaluation of the isotropic differential scattering cross section. We will seek a solution for the rotationalaverage two-electron scattering matrix elements, which are a prerequisite for a more detailed description later on, that involves the role of nuclear motion. From Eq. (7), we will need to evaulate the expression: dr r r 1 dr r r 2 ρ (2) I J (r r r 1 , r r r 2 ,R R R)e ιq q q·(r r r 1 −r r r 2 ) , In the absence of energy resolution, i.e. the standard set-up for current time-resolved scattering experiments, the separation of the total scattering into elastic and inelastic components is not needed, but it is instructive to show that our methodology is applicable to all four cases: total, elastic, inelastic and coherent mixed. Thus, we can split the expression for total scattering (I = J) using Eq. (8): where terms with K = I are elastic, while K I are inelastic. The key quantities in Eq. (12) and Eq. (13) are the one-and two-electron density functions, ρ (1) I J (r r r 1 ) and ρ (2) I J (r r r 1 , r r r 2 ), which can be expressed as weighted products of molecular orbitals (dropping the parametric dependence on the nuclear coordinates), 39 and where the indices run over all occupied molecular orbitals, N MO . The terms D II i j and d II i jkr are referred to as the elements of the one-and twoelectron reduced density matrix (1-and 2-RDM), respectively. If I J, they are known as one-and two-electron reduced transition density matrix elements.
Introducing Eq. (14) and Eq. (15) into the expression for the isotropic scattering, Eq. (12) and Eq. (13), shows that the efficient ab initio solution to the isotropic scattering problem requires the evaluation of integrals of the type, where A common strategy in molecular electronic structure theory is to expand the orbitals as a weighted sum of N bf primitive Cartesian Gaussian-type orbitals (GTOs), where A A A i i i is the centre of the i th primitive Gaussian and M (a) i the molecular orbital expansion coefficient. In the case of contracted Gaussian functions, i is premultiplied by a contraction coefficient. The sums of l i , m i and n i specify the orbital angular momentum. The integral K I J then takes the form, with the two-electron charge density r , and where we have labeled the Fourier integrals over r r r 1 and r r r 2 as J i j (q q q) and J * kr (q q q), respectively, J i j (q q q) = dr r r 1 g i (r r r 1 )g j (r r r 1 )e ιq q q·r r r 1 (20) and J * kr (q q q) = dr r r 2 g k (r r r 2 )g r (r r r 2 )e −iq q q·r r r 2 .
Eq. (19) reveals that the calculation of the ab initio isotropic scattering requires the evaluation of N 4 bf angular integrals. This is a formidable computational challenge, even for the smallest molecules with an adequate basis set. However, schemes for the efficient evaluation of such integrals have been proposed before 25,28 and make use of the properties of the Gaussian functions and their analytical Fourier transforms. The first step is to express the angular momentum properties of the product of two primitive Gaussians, g i and g j in a derivative form, 43 Π i j (r r r) = g i (r r r)g j (r r r) where we do not explicitly show the dependence of the McMurchie-Davidson coefficients, Λ , on the Cartesian components of the centres of the Gaussian functions and Gaussian exponents. The expression also exploits that the product of two Gaussians functions is a new Gaussian function, leading to the following definitions that appear in Eq. (22), The McMurchie-Davidson expansion is widely used in computational chemistry programme packages to express integrals of Gaussian functions with high angular momentum as derivatives of stype Gaussian integrals. Their utility stems from the existence of a recursive relationship between the coefficients, which enables their rapid evaluation. Substituting Π i j (r r r) into the expression for J i j (q q q) and taking the derivative in front of the Fourier integral, results in, J i j (q q q) = dr r r Π i j (r r r)e ιq q q·r r r dr r re −γ P (r r r−P P P) 2 e ιq q q·r r r .
The Fourier Transform of the s-type Gaussian function can be evaluated using the Fourier shift property, dr r re −γ P (r r r−P P P) 2 e ιq q q·r r r = π γ P 3/2 e −q 2 /4γ P e ιq q q·P P P . (25) Having analytically performed the Fourier Transform of the product of two arbitrarily GTOs, the full expression for J i j (q q q)J * kr (q q q) prior to integrating out the angular dependence becomes, J i j (q q q)J * kr (q q q) = dr r r 1 Π i j (r r r 1 )e ιq q q·r r r 1 dr r r 2 Π i j (r r r 2 )e −iq q q·r r r 2 where e ιq q q·(P P P−Q Q Q) .
Introducing H H H = P P P − Q Q Q and the combined angular momentum quantum numbers, L = L 1 + L 2 , M = M 1 + M 2 and N = N 1 + N 2 , Eq. (27) can be written as, It should be recognized that, if there was no need to perform the rotational average calculation, i.e. for scattering in the molecular frame, the derivatives in Eq. (28) trivially evaluate to However, even in the case of isotropic scattering the quantity J i j (q q q)J * kr (q q q) θ q φ q expressed in its current form has a relatively simple analytic solution. Resolving the angular integrals, which only affects e ιq q q·H H H results in, Solutions to the the equation above are discussed by Wang et al., 25 where it is given as a four dimensional sum over trigonometric functions scaled by precalculated numerical factors. Here, we follow more closely the approach suggested by Crittenden et al., 28 who calculated the result analytically for a limited number of angular momenta as a sum of spherical Bessel functions. In contrast to their approach, we recognize the existence of a recursive relationship between the expansion coefficients, which allows for a fast calculation and handling of arbitrarily large angular momenta. In the simple case when H < cut , the exponential in Eq. (28) is approximately unity and Eq. (29) takes the form, where B LMN = sin L+M θ q cos N θ q cos L φ q sin M φ q θ q φ q (31) is a numerical constant. As discussed above, in the case when H ≥ cut , the evaluation of the derivatives relies on the properties of the spherical Bessel functions, j β (qH), where β = min[(L + M + N − p − s − t)/2] + p + s + t, with min[] denoting the least integer greater than or equal to the quantity in the brackets. The coefficients a p L , b s M and c t N are related to the Hermite polynomials and obey the following recursive re-lation (here given for a p L ), The implementation of a sensible cut-off value, cut , is essential for the numerical stability of the algorithm, which might otherwise be affected by prohibitively large values of the (q/H) β factor in Eq. (32). The efficient procedure for evaluating Eq. (29) makes use of the recursive relationship for the coefficients a p L , b s M and c t N as well as of the recursive formula for the spherical Bessel functions. In addition to that, it is of paramount importance for the computational efficiency to take into account the symmetries J i j (q q q)J * kr (q q q) θ q φ q = J kr (q q q)J * i j (q q q) θ q φ q and J i j (q q q) = J ji (q q q), which together result in a speed-up factor of approximately 8. Another important simplification stems from careful consideration of the contraction scheme of the basis set used. If a given primitive GTO is a part of multiple contractions, the corresponding integral should be performed only once and the constant Z I J i jkr needs to be modified to reflect the combined contribution of this primitive to the molecular orbitals. The final trick for improving the computational performance is to treat together all GTOs whose centres and exponents are the same. Careful examination of these cases shows that they ultimately lead to the same values of H H H and can only differ in their angular momentum numbers l, m and n. Treating them together results in a family of F LMN (q q q, H H H) θ q φ q integrals, for which most of the terms in Eq. (32) are shared. In addition, a global cut-off linked to the relative size of z I J abcd could significantly speed up the calculation, at the expense of an effective decrease of the total number of electrons integrated. We have noticed that allowing for 0.1% electron density loss can lead to an approximate speed-up factor of two in most molecules explored without significant effect on the results.

Benchmarking and scaling
Our method is extensively tested and validated in a series of ab initio calculations for the ammonia molecule, NH 3 , varying the basis set and size of the active space. All ab initio calculations are performed with the MOLPRO electronic structure software package. 44,45 We chose ammonia because it has previously been used by Hoffmeyer et al. 29 to illustrate the importance of multiconfigurational wavefunctions in total x-ray and electron scattering. Although, as shown above, our methodology also encompasses individual inelastic and coherent mixed transitions, we focus on total and elastic scattering here. The discussion of inelastic scattering is limited to the cumulative sum of all inelastic transitions, which is given by the difference between total and elastic scattering. The richness of the information encoded in the inelastic and coherent mixed terms will be the target of a follow-up publication. Fig. 2 compares total and elastic scattering at different levels of theory, namely Hartree-Fock (HF) and Complete Active Space Self-Consistent Field (CASSCF) with 6 and 8 active orbitals and all electrons active, i.e. CASSCF (10,6) and CASSCF (10,8), respectively.
The basis sets include Pople's and Dunning's correlationconsistent basis sets with double-zeta, doublezeta plus diffuse functions and triple-zeta plus diffuse functions. Calculations with STO-3G minimal basis set are also performed and are included in SI Table 1 and SI Table 2, which summarise the benchmarking results. The geometry used in all calculations is optimized at the CASSCF(10,8)/aug-cc-pVTZ level of theory. The results are first calculated as the fractional signal change (commonly expressed in percent), We choose to use the fractional signal change (in percentage), %∆S (q), to ensure that the difference at each momentum transfer vector is relative to the absolute value of the intensity at this point. It is worth noting that the signals in the limit of large q for the total and elastic scattering differ by the number of the electrons in the molecule, which results in a larger value in the denominator in Eq. (34) and a smaller integral in Eq. (35) for total scattering. It is clear that the two families of basis functions are comparable at a given level of theory. Interestingly, the Pople's basis sets with the smaller active space and Hartree-Fock seem to be closer to the reference CASSCF(10,8)/aug-cc-pVTZ compared to the correlation-consistent basis sets. It is possible that the split-valence basis sets provide a better description in cases where the active space is not sufficient to capture the static electron correlation adequately. In addition, we found that STO-3G is largely unsuitable for scattering calculations (SI Table 1 and SI Table 2). In fact, it is comparable to the Independent Atom Model (IAM), 19 which gives E = 17.2 Å −1 and E = 32.8 Å −1 for total and elastic scattering. Overall, in both Pople's and Dunning's basis sets, the convergence is smooth with increasing the number of basis functions. The most significant effect is seen from going from the smaller CASSCF (10,6) to the larger CASSCF (10,8) active space, which can be attributed to the effects of electron correlation on the electron density of the system. We note that the total scattering is two-electron property and will be significantly affected by electron correlation. However, it is clear that elastic scattering convergence is equally affected suggesting that the electron density relaxation associated with static correlation, implicit in the CASSCF calculations, is a major factor in that case. For accurate total scattering calculation, ab initio methods that capture dynamic correlation are desirable. Fig. 3 compares the computational time required for the calculations discussed above. Generally, the calculations scale with the fourth power of the number of basis functions, N 4 bf , as all distinctive Fourier transforms over four basis functions need to be considered if they are not related by permutational symmetry. As discussed above, grouping together integrals with equivalent origin in Cartesian space and exponents partially offsets this scaling. Practically, as the number of atoms increases, so does the number of unique centres so that the speed up for large molecules could be small. In the case of elastic scattering, the two-electron charge density Z I J i jkr can be expressed as two independent pairs Z I J i j and Z I J kr . After precomputing these results, they can easily be accessed in the calculation with no added computational cost. Hence, the elastic scattering is largely independent of the level of theory or active space used. For the total scattering, the two-electron charge density remains a function of all four GTOs, so that its on-the-fly calculation increases the computational time compared to the elastic scattering in a manner that scales with the number of active orbitals. The scalings for specific inelastic and coherent mixed terms are equivalent to elastic and total scattering, respectively. Overall, we note that the computational time required per scattering calculation with the largest basis set used in this work is three orders of magnitude more expensive than a simple calculation with a minimal basis set. This scaling is of great practical importance for choosing an optimal method when considering scattering from a wide range of molecular geometries along a reaction coordinate or for the purpose of iterative inversion of experimental data. Given the comparatively small overhead cost of total scattering, it is sensible to adopt this type of calculation in lieu of the common approach of calculating total scattering as a sum of elastic scattering and tabulated inelastic corrections. 9 In Fig. 4 we show total, elastic and inelastic scattering curves of NH 3 at the reference CASSCF(10,8)/aug-cc-pVDZ level. The comparison with previous MR-SDCI (Multi Reference Single and Doubles Configuration Interaction) calculations by Hoffmeyer et al. 29 shows rather good agreement given the differences in the methods, levels of theory and basis sets. Their approach relies on numerical integration, whereas our result is strictly analytical. Our best calculation employs (a) NH 3 total scattering (b) NH 3 elastic scattering Figure 2: Convergence of the total and elastic X-ray scattering signals of NH 3 computed with HF, CASSCF(10,6) and CASSCF (10,8) using various basis sets. The height of the bars represent the integral of the percentage intensity change with respect to CASSCF(10,8)/aug-cc-pVTZ total and elastic scattering for the range 0 ≤ q ≤ 11.34 Å −1 . Figure 3: Comparison of the computational time required for the calculation of the isotropic elastic and total x-ray scattering signals in NH 3 . Different ab initio methods are considered as well as basis sets. Note that the computational scaling is expressed as a logarithm of the CPU time in milliseconds. The solid part of the bars represent the elastic scattering, whereas the shaded area at the top shows the extra time required to compute the total scattering with the same level of theory and basis set.
aug-cc-pVDZ, while their work reports a smaller double-zeta basis set with polarization and diffuse functions, [5s3p2d/3s2p]. Unlike MR-SDCI, CASSCF calculations account only for static electron correlation. The dynamic electron correlation has a smaller effect on the elastic component of the scattering signal as seen by comparing to the MR-SDCI results. Its influence increases when total and inelastic scattering are considered, demonstrating the importance of electron correlation for these quantities. A systematic study of the effect of electron correlation is critically important for fully understanding gas-phase X-ray scattering experiments and will be addressed in subsequent work.

Excited state total scattering
Here, we consider a simple model that illustrates both the utility of our methodology and the nature of the signal detected in ultrafast x-ray scattering. For our purpose, a suitable candidate is the photoexcitation of ammonia to the first excited singlet state, whose initial dynamics follows an umbrella motion. 46 The goal is to track the changes and dominant contributions to the total scattering as the geometry changes. In order to achieve that, we first optimize the ground state ion geometry at the CASSCF(9,8)/6-31+G* level as an approximation for the first excited state, which has Rydberg 3s character. The normal modes are calculated. The molecular geometry is then displaced in a series of steps along the umbrella mode  29 Subfigure 4a shows the total intensity, while subfigure 4b shows the difference.  1,5,8,9), where I tot exc (q,R R R) and I tot gs (q,R R R 0 ) are the total scattering intensities for the excited state and the ground state, respectively. This expression gives the change of the signal for the excited state at a specific geometry,R R R, relative to the scattering from the ground state at its equilibrium geometry, R R R 0 , under the assumption that there is no geometry change in the ground state upon excitation. In order to investigate its underlying contributions, it is conceptually useful to rewrite the expression in Eq. (36) as sum of two contributions, ∆S tot (q,R R R) = ∆S tot elec (q,R R R) + ∆S tot nucl (q,R R R) = I tot exc (q,R R R) − I tot gs (q,R R R) I tot gs (q,R R R 0 ) The first term, henceforth called electronic, shows the difference solely due to the electronic redistribution at any given geometry. The second nuclear term indicates the contribution due to structural changes, and is defined with respect to the electronic structure of the ground state only. Furthermore, making use of the fact that the total scattering is a sum of elastic and inelastic scattering, each of these terms can be split into two contributions, The breakdown of the total signal in terms of these four components is given in Fig. 5. The umbrella motion is tracked from 0 to 33.0 degrees in the pyramidization angle formed between the plane of the hydrogen atoms and one of the N−H bonds. The displacement along the normal mode is accompanied by a N−H bond elongation from 1.03 Å to 1.23 Å, which is seen to be the dominant factor for the variability in the nuclear part of the signal. The ground state equilibrium geometry used for this work has a pyramidization angle of 23.5 degrees and an N−H bond length of 1.02 Å. This implies that the planar geometry of the excited state is more similar to the ground state geometry, with referene to the N−H distance that dominates the scattering. This is reflected by the small magnitude of the nuclear scattering in Fig. 5a for this geometry. As the pyramidization angle increases, so does the bond length, which ultimately results in a maximum amplitude of %∆S of about 14%.
Meanwhile, the elastic part of the electronic contribution in Fig. 5b shows much less variability with the change of the geometry. The overall shape seems to change very little and the maximum amplitude of the central peak ranges only from 7% to 13%. It is of paramount importance to point out that the magnitude of the electronic component is comparable to that attributed to the nuclear motion. While at large nuclear pyramidization angles, the nuclear contribution is approximately two times larger, at near-planar geometries the elastic signal is almost exclusively attributed to the redistribution of the electronic density in the Rydberg state. Given the range of momentum transfer vectors for which the dip in Fig. 5b occurs, the observation can be explained by the effective loss of electron density in the molecular core associated with the delocalization of the Rydberg electron in the excited state. The small changes in the electronic component along the umbrella mode align with the fact that Rydberg electrons are not strongly affected by structural evolution of the ionlike core.
As expected, inspection of Figs. 5c and 5d reveals that the magnitude of the inelastic scattering is generally about 10 times smaller than the corresponding elastic contribution. Nonetheless, the size of the inelastic contribution is clearly large enough to have a tangible effect on interpretation of the experimental data. Here, the patterns parallel those seen in the elastic scattering. The electronic contribution is relatively constant with a magnitude comparable to the nuclear effect. The largest pyramidization angle has the largest nuclear component and the smallest electronic, with about twofold difference between the two. The planar NH 3 is dominated by the electronic scattering. Interestingly, the geometry dependence of the inelastic scattering is rarely accounted for in experiments. Specifically, the inelastic scattering is often approximated as an incoherent sum of stationary inelastic Compton factors for individual atoms. This is clearly a poor approximation in this case and most likely in general when considering small effects in the fractional signal change in time-resolved gas-phase experiments. 23 Given that the real experimental observable is the total scattering that stems from the Fourier transform of the two-electron electron density, it seems natural to consider the inelastic effect on an equal footing with elastic scattering. As seen here, inelastic scattering can account for up to 10% of the fractional signal change partially attributed to the difference between the two states considered, and exhibits a geometry-dependence similar to elastic scattering.
When pumped into the excited state, ammonia undergoes fast umbrella motion. The picture that emerges from this work is that the observed signal will alternate between two extrema driven by the elongation and contraction of the N−H distance. However, as shown in Fig. 5e, the baseline for this oscillation is set by the shape of the relatively constant difference in the scattering signals between the Rydberg state and the ground state. The latter is attributed to the electron-density hole brought about by the promotion of an electron from the HOMO to the diffuse 3s Rydberg orbital. The inelastic scattering has smaller but far from negligible effect on the scattering signal.

Conclusions
The mathematical framework and computational approach presented in this article allows for the efficient calculation of the isotropic total, elastic, inelastic and coherent mixed scattering signals. In order to perform the integration over the Euler angles needed to achieve spherical averaging, we consider ab initio wavefunctions expressed in a basis of Gaussian-type orbitals. Analytic solutions to the Fourier transform from real to reciprocal space results in a series of products of Gaussian and spherical Bessel functions, which are relatively easy to evaluate computationally. The approach is benchmarked against previous numerical calculations in the case of ammonia. We demonstrate the scaling of the algorithms with the basis sets and levels of theory used.
The methodology described in the paper is utilized to investigate a simplified model of the photoexcitation of ammonia to its 3s Rydberg state.
The observed elastic scattering signal shows a strong signature of the shift in the electron density associated with the promotion of an electron from the HOMO to the 3s Rydberg orbital. The magnitude of this purely electronic effect is comparable to the geometry dependent part of the signal. In addition, ammonia shows strong change in the inelastic scattering upon excitation, which is driven by changes in the electronic structure. It was furthermore shown that changes in geometry also play a role in inelastic scattering.
As ultrafast gas-phase x-ray scattering experiments are becoming more and more successful in obtaining high-quality data, it is of paramount im-portance to have the right tools to analyze the results. The internal dynamics, both nuclear and electronic, is encoded in the isotropic part of the signal, which can be extracted by means of Legendre decomposition of the detector signal. The isotropic signal should be understood as the spherical average of the Fourier transform of the correlated two-electron density of the molecule. The approach presented here allows this signal to be calculated efficiently for ground and excited states and can be used to aid the interpretation of pumpprobe ultrafast x-ray scattering experiments.
Going further, a similar mathematical apparatus can be applied to the case of static molecules or arbitrary high-order terms in the Legendre decomposition of the fully dimensional signal. The former can be achieved readily by forgoing the Spherical Bessel function expansion, while the latter can be achieved, albeit in a less straightforward manner, by using higher-order spherical Bessel functions. Equally important is the question of the ab initio level of theory used to calculate the molecular wavefunctions. Given that the signal is related to the two-electron density, it can be expected that the total scattering shows high sensitivity to electron correlation. Hence, thorough investigation of the impact of the post-Hartee-Fock methods is urgently needed. (e) Total scattering, %∆S tot (q) Figure 5: Breakdown of the total percentage fractional intensity change along the umbrella mode in the 3s Rydberg state of ammonia relative to the ground state equilibrium geometry. The displacement along the normal mode is labeled in terms of the pyramidization angle between an N−H bond and the plane of the three hydrogen atoms. The nuclear elastic (a) and nuclear inelastic (c) terms indicate the changes due to the geometry evolution in the elastic and inelastic scattering, respectively. Similarly, the electronic elastic (b) and inelastic (d) components show the changes in scattering as a result of the difference between the electronic structure of the ground and the excited state at each geometry. The total scattering, i.e. the sum of (a)-(d), is given in (e).