Signatures of electronic and nuclear coherences in ultrafast molecular x-ray and electron diffraction

Femtosecond x-ray and electron diffraction hold promise to image the evolving structures of single molecules. We present a unified quantum-electrodynamical formulation of diffraction signals, based on the exact many-body nuclear + electronic wavefunction that can be extracted from quantum chemistry simulations. This gives a framework for analyzing various approximate molecular dynamics simulations. We show that the complete description of ultrafast diffraction signals contains interesting contributions involving mixed elastic and inelastic scattered photons that are usually masked by other larger contributions and are neglected. These terms include overlaps of nuclear wavepackets between different electronic states that provide an electronic decoherence mechanism and are important for the time-resolved imaging of conical intersections.


I. INTRODUCTION
Over the past century, x-ray and electron diffraction signals have become the main tool for exploring the structure of matter. 2 The stationary x-ray diffraction (XRD) pattern in crystals is given by the square of the Fourier transform of the ground state electronic charge density. Electrons, in contrast, are diffracted from the total (electronic þ nuclear) charge density. Diffraction is usually dominated by the 2-molecule (or more generally 2-scatterer) coherent response that is responsible for the formation of Bragg peaks. For single molecule diffraction or diffraction in the absence of order, e.g., from liquids, the structure factor vanishes due to the absence of a long-range order and a continuous pattern is observed in momentum space. 21,28 In contrast to the 2-scatterer case, this 1-molecule diffraction pattern may not be written as a modulus square of a scattering amplitude. 6 Stationary x-ray and electron diffraction signals directly probe the charge densities. The advent of x-ray Free Electron Lasers (XFELs) has permitted the measurement of diffraction patterns from nanocrystals and time-evolving structures 8,13 with subfemtosecond resolution. 7,9,12,22,23,34 It is tempting to assume that for a time-evolving charge density, these signals simply give instantaneous snapshots of the charge density. This picture is, however, incomplete and the signals depend on coherences and inelastic scattering that go beyond the charge density. 26 For an excited molecule prepared in a time-evolving superposition of many-electron states, the signal may not be expressed in terms of the instantaneous charge density alone.
Here, we derive general expressions for diffraction signals in molecules that are ready to employ ab initio electronic structure calculations of nonadiabatic dynamics in the joint electronic þ nuclear space. Specific signatures of conical intersections are identified. A correct description of the signals requires treating the charge density as an operator and taking into account its different matrix elements. Since the diagonal matrix elements scale as the number of electrons in the system and off diagonal elements have contributions from a few electrons involved in the relevant electronic transitions, the ultrafast XRD or UED (ultrafast electron diffraction) signals are dominated by the former. The off diagonal elements that contribute to inelastic signals are thus often neglected. The independent atoms model (IAM) is widely used to describe the diffraction signals, 17,24 which assumes that each atom contributes to the signal with its own structure factor. This approximation fails to describe valence electrons involved in chemical bonding. This is justified for stationary diffraction but misses bond forming and breaking. Ab initio calculations are widely used to recover the experimental XRD and UED signals. 30,32,33 There exist different approximations for the XRD and UED signals, owing to the different levels of treatment of the molecular wavefunction. Phenomenological descriptions are often sufficient to reproduce specific experiments, but they miss certain features of the diffraction signal.
We present a rigorous formulation of time-resolved diffraction signals in the joint electronic þ nuclear space that contains all contributions to the signals. Contributions that mix diagonal and off diagonal matrix elements of the charge density are highlighted. These involve nuclear wavepackets in different electronic surfaces and are of importance in the study of decoherence effects in conical intersection. The electronic and nuclear charge densities are single-body operators in their respective subspaces but many-body effects enter through the many-body wavefunction used to compute their matrix elements. We show that such contributions can in principle be extracted by the separate detection of elastic, inelastic, and non-frequency resolved terms. Covariance techniques that make use of the stochastic properties of FEL SASE light can be adopted to that end. 3 Electron diffraction probes the total charge density. 10 Interestingly, the terms often missed in XRD (mixed terms of diagonal and off diagonal electronic charge densities) also enter in UED through mixed nuclear-electronic charge densities (diagonal nuclear and off diagonal electronic densities).

II. THE ELECTRONIC CHARGE DENSITY OPERATOR
The electronic charge density is a single body operator given by where n is the number of electrons and w † ðr 0 Þ and wðr 0 Þ are the fermionic field creation and annihilation operators expanded in a singleelectron basis / a , wðrÞ ¼ X a / a ðrÞc a : Substituting Eqs. (2) and (3) in Eq. (1) gives This expression can be recast as where c ð1Þ ðr; r 0 Þ ¼ w † ðr 0 ÞwðrÞ is the one-electron density operator. X-ray diffraction measures the Fourier transform of the charge density operator where q is the momentum transfer vector. XRD carried out on a single molecule measures the expectation value of the product of two charge density operators In real space, this gives r E ðrÞr E ðr 0 Þ ¼ e 2 w † ðrÞwðrÞw † ðr 0 Þwðr 0 Þ ¼ c ð2Þ ðr 0 ; r; r 0 rÞ þ c ð1Þ ðr; rÞdðr À r 0 Þ ; where c ð2Þ ðr 1 ; r 2 ; r 0 1 ; r 0 2 Þ ¼ w † ðr 0 2 Þw † ðr 0 1 Þwðr 1 Þwðr 2 Þ is the twoelectron density operator.

III. THE X-RAY DIFFRACTION SIGNAL FROM A SINGLE MOLECULE
The XRD signal is defined in Appendix A as the integrated rate of photon number change during the scattering process. Starting from Eq. (A5) and assuming a very short (impulsive) x-ray pulse, the single molecule x-ray diffraction signal can be recast as where jWðTÞi is the molecular many-body wavefunction in the joint electronic and nuclear space. Note that although r E ðrÞ is a single body electronic operator, its expectation value is calculated in the joint nuclei þ electrons space and depends on nuclear many-body wavefunction overlaps. In the following, we expand jWðTÞi in the adiabatic basis set consisting of products of nuclear and electronic wavefunctions where ju i i are the many-body electronic states which may depend parametrically on nuclear coordinates and jv i ðTÞi are the nuclear many-body wavefunction in the state i. jv i ðTÞi in Eq. (10) are not normalized and include the amplitude in-state i. Alternatively, normalized nuclear wavefunctions c i ðTÞjṽ i ðTÞi ¼ jv i ðTÞi can be introduced, where c i ðTÞ is the amplitude in-state i. This wavefunction acts in the jr; Ri space where r and R are the electronic and nuclear coordinates, respectively, The adiabatic basis set provides a convenient representation of the exact molecular wavefunction; note that we are not making the adiabatic approximation. Electronic structure simulations are usually performed in the adiabatic basis under the Born-Oppenheimer approximation, which states that the electron and nuclear motions can be separated due to their different timescales. The resulting adiabatic electronic states u i ðr; RÞ depend parametrically on the nuclear configuration. Exact quantum dynamical simulations can then be performed on these surfaces by including non-adiabatic couplings that account for conical intersection regions, where the motion of electrons and nuclei becomes strongly coupled and non-Born-Oppenheimer effects kick in. This yields a time-dependent nuclear wavepacket v i ðR; TÞ in each adiabatic state, and the exact many-body wavefunction jWðTÞi, Eq. (10). Other basis sets can be chosen, e.g., the diabatic or exact factorization of the molecular wavefunction. 1 This affects the wavefunction and the formulation of the diffraction signal, but the terms and contributions discussed here are incorporated there in a similar fashion. Different techniques and approximations have been developed to describe the non-adiabatic dynamics in large molecules. These include surface hopping, 4,25,29 frozen Gaussian approaches of various Structural Dynamics ARTICLE scitation.org/journal/sdy types, 5,18,31 or the multi-configurational time-dependent Hartree. 20 These approximations of the molecular wavefunction miss certain contributions to the diffraction signal. The present formulation is based on the exact electronic and nuclear wavefunction, and thus a offers a benchmark for other approximate diffraction simulations.
Inserting the many-body wavefunction, Eq. (10), into Eq. (9) yields the various contributions to the 1 molecule XRD signal depicted in Fig. 1, When i ¼ k or k ¼ j, the matrix elements r ik or r kj are the ordinary electronic charge densities that involve all the electrons. On the other hand, when i 6 ¼ k, the matrix elements represent transition charge densities that only involve the orbitals involved in a particular transition. As such, they are more sensitive to charge migration between the ground and valence excited states that typically only involve a few electrons.
The various contributions to the signal are conveniently represented by loop diagrams that represent the time evolution of the bra and ket of the molecular many-body wavefunction and the perturbative interactions with the incoming x-ray or electron beams. 19 If the electronic wavefunction is independent of the nuclear coordinates (the crude-adiabatic basis 16 ), then the electronic operators in Eq. (12) can be factorized into electronic and nuclear parts and the signal then depends on an overlap of the nuclear wavepackets For a two electronic state model with a ground state g and an excited state e, Eq. (12) The first two terms are the elastic contributions from the ground and excited states, whereas the following two terms represent inelastic Stokes and anti-Stokes processes. Finally, the last two terms represent mixed elastic/inelastic contributions. As can be read from the last two diagrams in Fig. 2, the mixed terms involve an elastic scattering involving r gg or r ee and an inelastic one with r eg or r ge . The two scattered photon amplitudes have frequencies centered around x X and x X 6 x eg , where x X is the incoming x-ray frequency and x eg is the transition frequency between states e and g. In order to generate a signal, a population must be created on the detector, which is possible only if these two amplitudes have a frequency within the detector bandwidth. The mixed terms can be observed only by a broadband detector with a bandwidth larger than x eg . This is usually the case for XRD detectors since in most cases x eg ( x X . This quantity depends on the molecular properties, i.e., the energy splitting of the involved electronic states. The latter is usually larger in the Franck-Condon region (between one and several eV), and smaller between higher excited states or in the vicinity of conical intersections.
XRD signals are usually computed either within the Independent Atom Model (IAM) or by ab initio techniques. The IAM completely misses the motion of valence electrons (bond forming and breaking). 24 Even with ab initio calculations, the last two contributions in Eq. (14) are often not included in the analysis of experimental data 11,17 but have been pointed out in former studies. 6,14,15,27 These are prominent at conical intersections, where the overlap of nuclear wavepackets provides a decoherence mechanism, and thus carries valuable information about the non-adiabatic passage.

IV. TWO-MOLECULE X-RAY DIFFRACTION SIGNALS
We now repeat the above discussion for the 2-molecule coherent diffraction. This signal arises from orientationally ordered molecular

FIG. 2. Loop diagrams for one-molecule XRD [Eq. (14)] and UED [Eq. (25)] expanded in electronic eigenstates for a two-state model. XRD is given by diagrams (i)-(vi), while all diagrams contribute to UED. The total XRD signal is the sum of electronic elastic [(i) and (ii)], inelastic electronic Stokes and anti-Stokes [(iii) and (iv)], and mixed electronic elastic/inelastic terms [(v) and (vi)]. UED has additional contributions from nuclear elastic [(vii) and (viii)], the mixed nuclear/electronic elastic [(ix) and (x)], and mixed nuclear/electronic elastic/inelastic [(xi) and (xii)] scattering.
Structural Dynamics ARTICLE scitation.org/journal/sdy assemblies such as crystals. In contrast to the one-molecule signal [Eq. (14)] described by two-point correlation functions hwðTÞjrrjwðTÞi, the two-molecule one does not contain products of the density operators, but absolute squares of single operator contributions. 6,19 This stems from the fact that the interaction with each of the two charge densities occurs on different molecules and the two-point correlation function can be factorized into a product of two single point ones. The signal, Eq. (A5), can then be recast as S XRD ðq; TÞ / jhr E ðq; TÞij 2 : By using the expansion of Eq. (10), we obtain For a two electronic states model, see Fig. 3, this gives S XRD ðq; TÞ / jhv g ðTÞjr E gg ðqÞjv g ðTÞi þ hv e ðTÞjr E ee ðqÞjv e ðTÞi þ hv e ðTÞjr E eg ðqÞjv g ðTÞi þ hv g ðTÞjr E ge ðqÞjv e ðTÞij 2 : The first two terms represent elastic scattering from states g and e, and the remaining two terms are the Stokes and anti-Stokes which represent inelastic scattering. When the detection has no frequency resolution, the last two terms become complex-conjugates and can be combined.
Expanding the square in Eq. (17) where each contribution can be read from Fig. 2,

VI. FREQUENCY-RESOLVED XRD
If the scattered light is frequency-dispersed, it is possible to separate the elastic and the inelastic (Stokes and anti-Stokes) contributions to the signal. We have three contributions S antiÀStokes XRD ðq; T; x s ¼ x x þ x eg Þ / <ðhv e ðTÞjr E ge ðÀqÞr E ge ðqÞjv e ðTÞiÞ: Without frequency resolution, one recovers Eq. (14). Equations (29)-(31) assume that one can separate the contributions. In realistic cases, one has to compute the contribution from the different states within the bandwidth of the pulses, see the diagram in Fig. 4, where a X ðx x Þ is the incident x-ray pulse envelope, Dx ¼ x s À x x andx ij is the transition frequency between states i and j. Importantly, x ij is an operator in nuclear space and contains inelastic vibrational and vibronic contributions to the inelastic scattering. Thus, the resonant factor is not factorized out of the expectation value since it also needs to be averaged over the nuclear wavepacket. Equation (32) is derived in Appendix B.
In order to extract the contribution including nuclear overlaps, one can detect the frequency-integrated signal, Eq. (14), that contains all contributions and subtract from it the ones measured with frequency-resolution, Eqs. (29)- (31).

VII. CONCLUSIONS
We have presented a rigorous formulation of ultrafast XRD and UED signals suitable for ab initio simulations. The main expressions are Eqs. (12), (14), and (16), (17) for the single molecule and the coherent two molecules signals, respectively. We have paid special attention to terms that are usually ignored in diffraction simulations arising from a superposition of states.
Contributions that mix elastic and inelastic scattering require a broadband detector. They usually make only a minor correction to the signal, which explains why they have not been studied so far. Targeted efforts are needed to extract them from the data. Subtraction of the elastic and inelastic diffraction will be measured separately with narrowband detection from the nonfrequency-resolved signal that could single out these terms. Alternatively, correlation measurements using stochastic light also have the potential to separate the various contributions. 3 These terms involve nuclear wavepackets in different electronic states and are thus highly sensitive to vibronic coherences. They can offer a novel window for conical intersections.
where r E is the molecular electronic charge density operator. The vector potential operator is given by where A x ðr; tÞ is the classical incoming x-ray field. In the last line, we have explicitly written the interaction picture density matrix propagator, where H int;À ðt 0 Þ is the Liouville space operator version of H int ðt 0 Þ given in Eq. (A1). Since the scattered photon mode s is initially in the vacuum state, this expression vanishes at 0th order in the probe field and a first order expansion must be done to recover the x-ray diffraction signal When the field amplitudes A x ðr; tÞ do not depend on the r (i.e., the incoming fields are plane waves), the integrals over the spatial coordinates correspond to Fourier transforms and we get Â hr E ðÀq; tÞr E ðq; t 0 ÞiA x ðtÞA Ã x ðt 0 Þe ixsðtÀt 0 Þ : (A6) Finally, assuming that the pulses are impulsive compared to the matter dynamics and incident at a delay T, A x ðtÞ ¼ A x dðt À TÞ, we can set t ¼ t 0 ¼ T and recover Eq. (9). Similarly, Eq. (15) can be recovered by assuming that the two-point correlation function in Eq. (A5) can be factorized over molecular pairs. 6

APPENDIX B: FREQUENCY-DISPERSED SIGNALS
Without frequency resolution, the signal is given by With the change of variable s ¼ t À t 0 , the integral becomes S XRD ðq; T; x s Þ / < ð dx X dx 0 X dtdse ixss hWðt 0 ÞjG † ðt À t 0 Þ Â r E ðÀqÞGðsÞr E ðqÞGðt À s À t 0 ÞjWðt 0 Þi Âa x ðx x Þa Ã x ðx 0 x Þe ixxðtÀTÞ e Àix 0 x ðtÀsÀTÞ ; where G(t) is the molecular propagator that contains the interaction with the actinic pulse. This propagation can be treated numerically in order to include possible non-adiabatic dynamics. Assuming that the pulse is short compared to the dynamics of the molecule, we can make the approximation WðtÞ ¼ Wðt À sÞ ¼ WðTÞ. The integral over t can then be computed and gives a Dirac delta function dðx x À x 0 x Þ, The integral over s can be computed assuming that GðsÞ is the field-free propagator ð þ1 0 dse iðxsÀxxþiÞs e ÀiL0s ¼ 1 Inserting this in Eq. (B6) and expanding over states lead to Eq. (32) in the main text.

DATA AVAILABILITY
The data that supports the findings of this study are available within the article.