Quantum state tomography of molecules by ultrafast diffraction

Ultrafast electron diffraction and time-resolved serial crystallography are the basis of the ongoing revolution in capturing at the atomic level of detail the structural dynamics of molecules. However, most experiments capture only the probability density of the nuclear wavepackets to determine the time-dependent molecular structures, while the full quantum state has not been accessed. Here, we introduce a framework for the preparation and ultrafast coherent diffraction from rotational wave packets of molecules, and we establish a new variant of quantum state tomography for ultrafast electron diffraction to characterize the molecular quantum states. The ability to reconstruct the density matrix, which encodes the amplitude and phase of the wavepacket, for molecules of arbitrary degrees of freedom, will enable the reconstruction of a quantum molecular movie from experimental x-ray or electron diffraction data.


Eq. 3 in the main text
Modification by experimental data Pr m1−m2 (θ, t) Eq. 6-Eq. 11 Converged Not converged Supplementary Figure 1. Schematic flow chart for imposing constraints to the wavepacket probability distribution. The internal procedure for the "constraints of density matrix" is separately elaborated in Supplementary Fig. 2. The superscript n represents n-th iteration.
Supplementary Figure 3. Simulated probability distribution and diffraction pattern of rotational wavepacket. The first row shows the initial angular probability for N 2 molecules prepared at a rotational temperature of 30 K and the expectation values of cos 2 θ of the time evolving wavepacket for N 2 molecules after laser pulse [2]. The alignment laser pulse is linearly polarized with a Gaussian envelope of duration τ L = 50 fs and 10 13 W/cm 2 peak intensity, and θ is the polar angle between the polarization and the molecular axes. The duration is much shorter than the characteristic rotational time τ L T . The second and third rows show the angular probability distribution changes from aligned to anti-aligned, and the difference of their diffraction intensity with respect to t = 0. The X-ray photon energy is assumed to be 20 keV. ) for λ ranging from 10 −5 to 10 8 . The Tikhonov regularization procedure minimizes I − KPr 2 2 + λ Pr 2 2 . The black point marked on the curve is the turning point corresponding to λ ≈ 10 4 . The yellow area starting from log λ = 1 and ending at log λ = 4 illustrates the admissible range of regularization parameter λ. Figure 5. Quantum tomography result of numerical trial with initial guess of correct diagonal elements of density matrix. The modulus of density matrix elements are shown in the upper panel, where J = |m|, |m| + 1, · · · , J max within each m-block. The phases of all density matrix elements are zero at t = 0. The lower panel shows angular probability distribution, the recovered modulus and phases of density matrix elements faithfully reproduce the reference Pr(θ, t), which is cylindrically symmetric in azimuthal direction of φ. Error functions of density matrix and probability distribution are shown in the rightmost column. Figure 6. Quantum tomography result of numerical trial with random initial guess of density matrix. Only the measured probability distribution and general properties of density matrix (being Hermitian, positive semidefinite and with unity trace) are imposed as constraints during the iteration algorithm. The density matrix to be recovered and its probability distribution are identical to that in Fig. 5. The modulus of density matrix elements are shown in the upper panel, where J = |m|, |m| + 1, · · · , J max within each m-block. The phases of all density matrix elements are zero at t = 0. The lower panel shows angular probability distribution, the recovered modulus and phases of density matrix elements faithfully reproduce the reference Pr(θ, t) (cylindrically symmetric in azimuthal direction of φ). Error functions of density matrix and probability distribution are shown in the rightmost column. initial guess ⟨n 1 n 2 · · · n N |ρ|m 1 m 2 · · · m N ⟩ Pr ∆1,∆2,··· ,∆N (x 1 , x 2 , · · · , x N )

Supplementary
Eq. 39 constraints of probability with measured Eq. 40 constraints of density matrix Supplementary Figure 7. Quantum tomography of vibrational state. The iterative transform is again between the spaces of density matrix and the blockwise probability distribution The details of the keV UED setup and experimental conditions for nitrogen alignment experiment have been previously introduced in [3,4]. We use a tilted infrared laser pulse to excite the rotational wave packet of the nitrogen molecule ensemble with a laser pulse duration of 60 fs, a spot size of 190 um (horizontal) × 260 um (vertical), and pulse energy of 1mJ. The tilted angle is about 60 degrees, which is designed to remove the group velocity mismatch due to the lower speed (0.526c, where c is the speed of light) of the electron pulse. The probe electron pulse is generated by shinning a 266 nm UV laser onto a copper cathode, which is accelerated by a 90 keV DC voltage and then compressed by a 3GHz RF electric field to minimize the temporal pulse duration on the sample. The electron beam is truncated using a platinum aperture with a diameter of 100 um to deliver a beam current of 8 pA, corresponding to 10,000 electrons per pulse. A de Laval nozzle with an inner diameter of 30 um is used to deliver the nitrogen molecules to the interaction as a supersonic molecular beam with a diameter of 200 um, and the nozzle backing pressure is 1200 mbar of nitrogen. The instrument response time was determined to be 240 fs by fitting the experimental anisotropy to its corresponding simulation. The timing jitter was 50 fs rms over several hours [4]. The electron diffraction patterns are recorded by an electron-multiplying chargecoupled device (EMCCD) camera, and the time delay between the pump and probe is controlled by an optical stage. Here the step of time delay is 100 fs.

SUPPLEMENTARY NOTE 2. DIFFRACTION PATTERN TREATMENT
The details of how to retrieve the angular distribution from the measured diffraction patterns have been explained in [4]. Briefly, the diffraction difference pattern for each image is calculated with ∆I(s, t) = I(s, t) − I(s, t < −1ps) to remove the background of atomic scattering, and then are averaged over the four quadrants using its symmetry. The simulated random molecular scattering with a rescaling factor of 0.35, which is obtained by fitting the experimental anisotropy evolution and its corresponding simulation, is added to ∆I(s, t) to recover molecular diffraction intensity I(s, t). The modified pair distribution function (MPDF) [4] is calculated by applying the inverse Fourier transform of I(s, t), followed by an Abel inversion, giving the information of angular distribution Pr(θ, φ, t).
The angular distribution retrieved from experimental data covers the initial alignment through the revivals up to about 7 ps, which is deconvolved using the algorithm in [5][6][7]. The instrumental point spread function (PSF) is assumed to be a one-dimensional Gaussian function with a full width at half maximum of 0.28 ps, corresponding to the general temporal resolution of the setup, for the deconvolution, which eliminates the blurring due to the limit temporal resolution of the setup. The general temporal resolution 0.28 ps is obtained by comparing the deconvoluted angular distribution evolution with the simulation and minimizing their difference. The temporal evolution of Pr(θ, φ, t) is extended to obtain the data up to 11ps by a reflection of the angular distribution evolution from 6.1ps to 1.2 ps to approximate the data from 6.1 ps to 11 ps according to the approximate symmetry based on the simulation.

WITH FIXED PROJECTION QUANTUM NUMBERS
We extend the treatment in Ref. [8] to show that the density matrix element J 1 m 1 |ρ|J 2 m 2 in the (m 1 , m 2 )-block subspace can be solved analytically, once the blockwise probability density Pr m 1 ,m 2 (θ, t) of given projection quantum numbers m 1 , m 2 is determined. We expand the blockwise probability density with eigenbasis, where the energy level difference is and I is the moment of inertia of the rotating molecule. For the sake of convenience, we define normalized associated Legendre polynomials with orthonormal relations We use the orthogonal relations of Legendre polynomials and exponential functions in the integral transformation [8]. Firstly, consider the motion along rotational polar coordinate θ. The product of two associated Legendre polynomials occur in Eq. 1 can be expanded by single associated Legendre polynomials Thus, integrate over θ, Let T = 2πI, which is the rotational revival period [9], and integrate over t, where the range of integration 2T is the minimal period of e i[β(α+1)/2I−∆ω]t because which is the integral multiple of 2π. The range of α and β is set to be |∆J| ≤ |β| ≤ α ≤ J, where β and ∆J are of the same sign. If β(α + 1) has unique integer factorization, the only term remaining in the sum satisfying is β = ∆J and α = J. The corresponding density matrix element can be derived as If the factorization of β(α + 1) is not unique, we calculate all integrations I m 1 m 2 (α , β ) where β(α + 1) = β (α + 1). For example, when β = 0, all of the ∆J = 0 terms remain. When changing the value of α, all these I m 1 m 2 and corresponding density matrix elements constitute a set of linear algebraic equations (where α = 2J can only be which has unique solution because all diagonal terms of the upper triangular matrix are nonzero.

SUPPLEMENTARY NOTE 4. LASER ALIGNMENT OF ROTATING MOLECULE
The effective Hamiltonian of rotating molecule-laser interaction is [2] where J is the rotational angular momentum, (t) is the electric field of the laser pulse, B = 1/2I is the rotational constant, α and α ⊥ are the components of the static polarizability, parallel and perpendicular to the molecular axes. The molecule is assumed to be in the vibrational and electronic ground state. An initial rotational eigenstate |J 0 M 0 evolves to a pendular state [2] where J and J 0 are of the same parity. The coupling coefficients d J 0 m 0 J is induced by laser field, satisfying selection rules ∆m = 0 and ∆J = 0, ±2. d J 0 m 0 J is invariant after the laser pulse, and the evolution of rotational angular distribution originates from interference of each dynamical phase. The coherence of the created quantum state can be maintained for several revival periods, and the alignment is reconstructed at predetermined times and survives for a perfectly controllable period [2], the sufficiently long coherence time makes the time evolution measurement of quantum state tomography feasible.
The initial system in thermal equilibrium can be characterized by the following density operator where ω J 0 is the Boltzmann statistical factor determined by the rotational temperature. The density operator of the laser-aligned system iŝ So the partial trace of m subspace with odd (or even) J is invariant in the dynamics of laser alignment, since it is a general property of laser-molecule interaction, where we used the normalization property of coefficients d J 0 M J (t) in Eq. 14.
Notice that density matrix of opposite magnetic quantum number m and −m is symmetric forρ ini , which also remains symmetric for transition matrix element induced by laser interaction H eff (t). From Eq. 13, taking into account selection rule ∆M = 0, In this section we show the detailed procedure for making an arbitrary density matrix and probability distribution to satisfy the physical constraints given in the main text. Most physical constraints are given in the summation form. For example, from Eq. 18, From the measured probability distribution and the constraint can be expressed as They can be sataisfied by scaling with a common factor Pr m 1 ,m 2 (θ, t) → β(θ, t)Pr m 1 ,m 2 (θ, t) , β = Pr m 1 −m 2 (θ, t) The constraints in probability space is given by Eq. 26, and illustrated with flow chart in Fig. 1.
Further constraints in density matrix space include being Hermitian, positive semidefinite and having invariant partial traces (the procedure is presented with the flow chart in Fig. 2).
As a general rule to guarantee the completeness of constraint conditions, we can firstly analyse the time-dependent molecular probability distribution Pr(θ, φ, t) can be obtained by solving the Fredholm integral equation of the first kind using Tikhonov regularization procedure [11]. We assume τ = − cos θ and replace the integral by Riemann summation, at each instant, where ∆φ = 2π a , ∆τ = 2 b , i is ranging from 1 to a, j is ranging from 1 to b, k is ranging from 1 to c, and l is ranging from 1 to d. φ and θ are the azimuthal and levitation angles of the linear molecular rotor, Θ and Φ are the scattering angle of the X-ray photon in the lab system (as is shown in Fig. 1 in the main text). We can write the total diffraction intensity in the matrix . . .
To avoid singular matrix inversion, we use Tikhonov regularization to get the rotational probability distribution, where E is identity matrix of size (c × d) and K T is the transpose of matrix K.
The Tikhonov regularization performs excellently in dealing with experimental data with measurement errors and preventing overfitting, and can faithfully recover the probability density distribution. To validate the faithfulness of the obtained probability distribution Pr(θ, φ), we define the condition number where A 2 = i A 2 i is the L 2 Euclid norm. The condition number characterizes the degree of variation of the solution Pr(θ, φ) with respect to the input data of measured diffraction intensity I(s), its value provides a measure for the sensitivity of the solution with respect to the measurement error and choice of regularization parameters. From Fig. 4, we can estimate that λ ≥ 10 is required to ensure cond ≤ 10, and subsequently to ensure the reliability of the solution.
Quantum tomography of the rotational wavepacket gives the result shown in Fig. 3 in the main text. After 50 iterations, both density matrix and probability distribution are precisely recovered.

MATRIX AND INITIAL GUESS
We have verified the new quantum tomographic method by the rotational wavepacket of a laseraligned molecule. We also illustrate the power of the new method by applying it to a randomly chosen density matrix rather than that in the laser-aligned case. The iterative QT algorithm also converges after about 20 iterations and density matrix is recovered with considerable accuracy.
We impose the error functions of density matrix and probability distribution to measure the accuracy of iteration results, which are defined by where the subscript n represents the result of n-th iteration, and 0 represents the correct result.
In Fig. 5 we show the result of identical algorithm given in Fig. 1 and Fig. 2, only with smaller order J max of density matrix to be recovered. The initial state is given by correct diagonal elements of density matrix. The iteration converged to the expected result with error 20 (ρ) = 3.5 × 10 −3 and 20 (Pr) = 1.7 × 10 −3 .
Especially, we show with the proof-of-principle example that this iterative QT algorithm is insensitive with the initial guess of density matrix. The rotational temperature which provides much information such as initial guess and partial trace, is actually not indispensable to the QT method. Assume we are dealing with a pure QT problem without any additional knowledge to the density matrix to be recovered. As is shown in Fig. 6, a random initial guess will also lead to a converged result after about 30 iterations with error 30 (ρ) = 3.9 × 10 −2 and 30 (Pr) = 9.0 × 10 −3 .

SUPPLEMENTARY NOTE 8. VIBRATIONAL QUANTUM TOMOGRAPHY
Vibrational quantum tomography recovers the density matrix of N vibrational modes from the probability distribution evolution Pr(x 1 , x 2 , · · · , x N , t) where φ n i (x i ) is the harmonic oscillator wavefunction of the i-th vibrational mode with energy eigenvalue (n i + 1 2 )ω i . The dimension problem arises naturally. Here the probability is (N + 1)dimensional and density matrix is 2N -dimensional, which is inadmissible for analytical solutions when N > 1. In conventional QT method that is based on integral transform, the orthogonal properties cancel out one summation by integrating over one parameter. For example, where T = 2π ω 0 . f mn (x) is the sampling function [12] defined by where φ m (x) and ϕ n (x) are respectively regular and irregular wavefunctions of harmonic oscillator. The bi-orthogonal properties of sampling function is under frequency constraints m − n = m − n .
Our theory, based on the following two procedures, fully utilizes the above orthogonal properties and imposes constraints for lack of dimension. First, we set up the transformation between probability and density matrix in a subspace Second, starting from an initial guess, effective physical constraints can be imposed by iterative projection method to get the converged result. For example, the priori knowledge of density matrix of being Hermitian, positive semidefinite and normalized. The algorithm of vibrational state QT and an example of 2D vibrational quantum tomography is shown in Fig. 7 and Fig. 8. The initial guess is given randomly, and only the probability distribution and general properties of density matrix are imposed as constraints during the iteration algorithm.
Similar to rotational QT, the dimension problem can be reflected by the fact that for unless only one single combination of {∆ i } satisfies N i=1 ∆ i r i = k, there is no direct way to obtain Pr ∆ 1 ,∆ 2 ,··· ,∆ N (x 1 , x 2 , · · · , x N ) from the measured wavepacket density distribution, only their sum can be available through Fourier transform of the measured probability distribution evolution where we assume ω i = r i ω 0 (r i are integers and T = 2π/ω 0 , r i 's are the set of smallest integers to represent the measured frequencies). In the new iterative QT method for N -dimensional vibrational system, we do not need infinitely long time of measurement anymore, which used to be indispensable to fill the whole space of N -dimensional phases [13] while physically infeasible.
Besides, in the new iterative QT method, the ratio of frequencies does not have to be irrational, which is important because in reality N -dimensional vibrational systems with commensurable frequencies are ubiquitous.
The pattern function can be approximated around x = 0 as [14] f nn ∼ − 2 π sin[−π(n + 1/2) In order to resolve a period of the oscillation of the pattern function that arises in the convolution (Eq. 40), the required spatial resolution for reconstructing vibrational density matrix up to Nth order has to be better than δx ≤ π/2 √ 2N + 1. The maximal order of the desired density matrix also sets demand on the temporal resolution. Suppose d time intervals are measured for a half period T /2 = π/ω 0 . From Eq. 42, we have a phase resolution of kπ/d for the Fourier transformation of probability distribution function. The aliasing phenomena defines the maximal order of density matrix we can access to be N = d/k − 1, thus the required temporal resolution is The quantum tomography procedure presented above can be easily generalized to systems with strong anharmonicity , or when coupling among different vibrational modes exist. In general case, the Hamiltonian [15] whereĥ i is the separable part for i-th vibrational mode with eigenstate φ n i (x i ), and V (x 1 , x 2 , · · · , x N ) is coupling potential among N vibrational modes. The eigenstate is a linear combination of product 1D wavefunctions assigned with quantum numbers I = {I 1 , I 2 , · · · , I N } with energy The wavefunction ansatz and the self-consistent field state interaction (SCF-SI) method can be extended to the treatment of molecular vibrations with anharmonic potential and strong coupling.
To solve the vibrational Hamiltonian efficiently for this general case, the separable part Hamiltonianĥ i for the anharmonic vibrational modes can be replaced by, for example, Morse potential model and corresponding pattern function [16] can be used. The iterative projection algorithm for quantum tomography should be set up based on the transformation between probability and density matrix in a subspace where the frequency constraint of sampling function requires ω iα − ω jα = ω ∆α (α = 1, 2, · · · , N ).
The density matrix element can be solved from the linear equation of 48. If there are n basis eigenstate for i-th uncoupled vibrational mode φ n i (x i ), the coupled density matrix can be recovered to the order of (2n) N/2 . Similarly, the procedure starts from an initial guess and imposes constraints to both density matrix space and probability space. Besides basic properties of density matrix and probability distribution, the subspace probability should also satisfy where ω I and ω J are energy eigenvalues of the coupled Hamiltonian, T is the common period for all vibrational frequency intervals.
To enhance the convergence of iterative QT procedure for vibrational states, physical constraints can be imposed on the diagonal matrix elements of the density matrix, which is experimentally accessible, e.g. through photoelectron spectra and absorption spectra, which can directly provide constraints on diagonal density matrix elements of basis states with eigenenergy E [17].
As a final remark, for vibrational QT, it is sometimes neccessary to use the velocities of nuclei as constraining physical conditions, in the case that the basis states of density matrix is energetically degenerate. For example, given the ratio of two vibrational frequencies r 1 /r 2 = 1/2, consider a mixed state consisting of |20 and |10 (the pure state is a special case of it), their density matrix is The probability distribution Pr(x 1 , x 2 , t) = ρ 11 φ 2 2 (x 1 )φ 2 0 (x 2 ) + ρ 22 φ 2 0 (x 1 )φ 2 1 (x 2 ) (51) could not reflect the imaginary part of the off-diagonal density matrix elements because the degeneracy of the two basis states smears out the temporal evolution of the probability distribution.
If |20 and |01 belong to the same symmetry representation, their coupling will lead to Fermi resonance and the degeneracy can be lifted. In the case that |20 and |01 are exactly degenerate, additional constraints must be imposed. Because with the ultrafast diffraction method, the velocity of nuclei and thus their momenta can be extracted experimentally, we can naturally construct physical constraints through products of momenta, such as p 2 contains information of imaginary part of non-diagonal density matrix elements Im[ρ 12 ] = −Im[ρ 21 ], with which we can effectively determine the imaginary part of the off-diagonal density matrix elements between exactly degenerate basis states, by using the products of velocities as physical constraints in the iterative QT procedure.
Throughout the paper, we focus on recovering the density matrix, which is interconnected with the Wigner function W (q, p) via the overlapping formula, where WÔ(q, p) = (1/2π) dx exp(−ipx) q − x 2 |Ô|q + x 2 . Especially, the Wigner function can be expressed in terms of the density operatorρ as W (q, p) = Wρ(q, p).