Filming movies of attosecond charge migration with high harmonic spectroscopy

Electron migration in molecules is the progenitor of chemical reactions and biological functions after light-matter interaction. Following this ultrafast dynamics, however, has been an enduring endeavor. Recently, it has been sug-gested that high-harmonic spectroscopy (HHS) is able to probe dynamics with attosecond temporal and sub-˚angstr¨om spatial resolution. Still, real-time visu-alization of single-molecule dynamics continues to be a great challenge because experimental harmonic spectra are due to the coherent averages of light emission from individual molecules of diﬀerent alignments. Here we demonstrate that the uniting of machine learning (ML) algorithm and HHS in two-color laser pulses enables us to retrieve the complex amplitude and phase of harmonics from single ﬁxed-in-space molecule. From the complex single-molecule dipoles for diﬀerent harmonics, we construct a movie of electron migration after tun-nel ionization of the molecules at time steps of 50 attoseconds. Moreover, the angular dependence of molecular charge migration is fully resolved. By exam-ining the movie, we observe that electron holes do not just “migrate” along the laser direction, but may “swirl” around the atom centers. Our ML-based HHS establishes a general reconstruction scheme for studying ultrafast charge migration in molecules, paving a way for further advance in tracing and controlling photochemical reactions by femtosecond lasers.


I. MAIN
Charge migration in molecules plays a fundamental role in numerous chemical and biological processes [1][2][3][4]. Following the ultrafast electron dynamics is of crucial importance for the understanding, and ultimately, the steering of complex chemical and biological reactions. However, such studies have not been possible in the last century since experimental tools reaching the typical natural temporal resolutions of few-femtoseconds or attoseconds are not available [5,6]. At the dawn of the 21st century, the advent of attosecond pulses based on an extremely nonlinear process known as high-order harmonic generation (HHG) breaks the femtosecond barriers, inaugurating a new era for time-resolved metrology [7][8][9]. By using attosecond pulse in combination with a femtosecond infrared (IR) pulse in a pump-probe conguration, various experimental techniques have been implemented to partially monitor ultrafast electron dynamics in atoms and molecules, such as attosecond extreme-ultraviolet (XUV) transient absorption spectroscopy [10,11], attosecond photoemission [12] and photo-fragmentation [13] spectroscopy, etc. With such conguration, charge motion was reported on a bio-relevant molecule (phenylalanine) in experiment [14] with temporal resolution of a few femtoseconds.
High-order harmonics from molecules also contains information on the structure and dynamics of molecules and High harmonic spectroscopy (HHS) recently has been established as a burgeoning technique for ultrafast molecular detection [15][16][17][18][19][20][21][22]. HHS relies on an "internal" built-in pumpprobe process in high-order harmonic generation (HHG) known as a three-step model [23,24]: ionization (pump), acceleration, and recombination (probe) [see Fig. 1(a)]. In this process, the electron dynamics is triggered at the ionization step, and probed at the recombination instant by the returning electron with information encoded in the harmonic spectrum. The temporal resolution of this unique pump-probe process arises from the intrinsic frequency chirp of harmonic emission [25], i.e., different harmonic orders are associated with well-defined ionization-recombination delays spent by the electron in the continuum [see Fig. 1(b)]. The typical resolution for a commonly-used 800 nm driving laser is ∼100 attoseconds. The spatial resolution of HHS comes from the de Broglie wavelength of the recombining electron and can reach sub-ångström [19][20][21][22]. With this technique, measurement of ultrafast nuclear wave packets created by ionization in isotopic molecules was realized [15,16]. Very recently, the HHS technique has been further nicely extended to extract attosecond charge migration in iodoacetylene molecule by Kraus et al. [18]. However, so far in the literature all the single-molecule dynamics has been directly taken from the experimental data averaged over an ensemble of angular distributed molecules. This practice has severely impacted 3 the accuracy of the extracted single-molecule dynamics.
As demonstrated in [18], to reconstruct electron dynamics of single molecules from experimental harmonic spectra, one has to address many challenges. First, reconstruction is an inverse scattering problem, which is generally solved by iterative methods. Reconstruction is possible only if the scattering theory is on firm ground. Second, since molecules can not be aligned perfectly in experiment, the measured harmonics are the result of a coherent superposition of individual radiation weighted by the angular distribution of the molecules [26][27][28]. Thus, the critical challenge is to extract fixed-in-space single-molecule harmonic amplitudes and phases. This second step has not been properly treated in all of existing publications using HHS techniques. By taking angular averaged harmonic spectra from the experimental data, the phases of harmonics are severely compromised (see Supplementary Note 3), which would incur the loss of dynamic information of electron wave packet that is at the heart of charge migration. Besides, to extract charge migration, extra parameters (i.e., the population coefficients of different cation states) that describe electron wave packet of a single molecule should be obtained as well. Thus, it is also necessary to evaluate if adequate experimental data are available to make sure that accurate retrieval of charge migration is possible.
Here we demonstrate that, by resolving the above three issues, we have succeeded in filming movies of charge migration in the molecular ion using HHS to reveal substantial sub-ångström electron migration. Since electron dynamics runs on the timescale of about 20 attoseconds (one atomic unit of time is 24 attoseconds), while the emission time difference between two odd harmonics by an 800-nm laser is about 100 attoseconds, we shorten this interval by half, by generating both even and odd harmonics with a two-color laser field. With harmonic spectra from both single-color and two-color lasers in standard pump-probe measurements, together with the powerful modern machine learning (ML) algorithm, we have abundant data to retrieve the complex dipole amplitude and phase of each fixed-in-space molecule to construct charge migration movies.
The alignment-angle-resolved molecular charge migration has been fully characterized for the first time.

Deciphering charge migration with HHS
We first explain how ultrafast charge migration is encoded in and how it can be identified from the molecular HHG signals. After strong-field tunneling ionization, the molecule in general is left in the ground as well as several nearby electronically excited states of the ion [1-4, [29][30][31]. Once the electron is ejected, the occupation amplitudes of these electronic states would change with   time under the influence of the external laser field, creating a time-dependent many-electron wave packet, or equivalently, an electron hole wave packet. The time dependence of the modulus square of the wave packet is called electron charge (or hole) migration. Such hole dynamics is encoded in the harmonic spectra of each single molecule [29][30][31].
To decipher charge migration from the harmonic spectra, we choose the two most widely studied molecules N 2 and CO 2 as the candidates to demonstrate our scheme. First, we consider N 2 . As illustrated in Fig. 1(a), the two channels (Channel X and A) relevant to HHG in N 2 are associated with the groundX 2 Σ g and first excitedÃ 2 Π u states of N + 2 ion. In molecular orbital picture, these two channels correspond to ionization of an electron from the two highest-occupied molecular orbitals, HOMO and HOMO-1, respectively. Owing to the different symmetries of these orbitals, HHG from different ionization channels depends differently on molecular alignment. Figure 1 to laser polarization, the signals of lower order harmonics (e.g., H21 and H25) are strong (weak) because they are generated from the HOMO orbital. For higher orders, like H31 and H33, the signals receive substantial contribution from the HOMO-1 orbital and thus becomes quite large when the molecules are anti-aligned. These properties of the harmonic spectra of N 2 have been well studied experimentally [32] and can be calculated with quantitative rescattering (QRS) theory [33][34][35][36]. Similar results are also observed in HHG from CO 2 molecule (see Supplementary Note 2). The different alignment dependence of the contribution to HHG from multiple molecular orbitals allows us to access their complex population amplitudes in the complex wave packet of the molecular ion, from which ultrafast charge migration can be uncovered.

Retrieval of single-molecule dipole
For a parallel configuration of the pump and probe polarizations, the time-dependent harmonic signal is given by [26,27] where θ is the alignment angle of the molecule, τ is the pump-probe delay, D(ω, θ) is the fixedangle-dependent total dipole moment from the single-molecule response, and ρ(θ, τ ) is the angular distribution of molecules at delay τ . In our experiment, ρ(θ, τ ) is determined with method in [37].
Note that perfect molecular alignment can not be achieved in experiment. Our simulations show   ML to this problem. The modern ML algorithm has been demonstrated to have remarkable abilities in characterizing complex sets of data with a high degree of accuracy, and has been widely utilized in genetics [38], condensed-matter physics [39], and material science [40]. We show here that it can effectively deal with the complicated decoding in HHS (see Methods and Supplementary Note 4).
The second step of our reconstruction is to disentangle the multichannel contributions from the total single-molecule dipole moment D(ω, θ) obtained in the first step. Driven by a one-color 800-nm probe pulse, only one set of D(ω, θ) can be obtained for each harmonic order which is insufficient for the reconstruction. To overcome this problem, we employed an additional parallel two-color laser field to generate harmonics. This two-color field consists of an intense 800-nm fundamental field and a weak second-harmonic (SH) field. The SH field is weak (∼2 × 10 −3 of the fundamental field) such that it barely alters the electron dynamics of the molecular ion (see Supplementary Note 5). Measurement of harmonics at different relative phases between the fundamental and SH fields replenishes the additional data set needed to retrieve multichannel contributions in the second step. Figures 2(a,b) show the time-dependent signals of H22 and H27 from N 2 measured as a function of the relative phase of the two-color laser fields. One can see that the modulation of HHG intensity depends sensitively on the relative phase, and different harmonic orders present different dependences. Applying the ML-based reconstruction procedure to each harmonic order, we can obtain the time-dependent complex mixing coefficients of the multiple orbitals of the molecular cation, for each fixed-in-space angle, at the excursion time when the recombination of that harmonic order occurs, as illustrated in Fig. 1(b).
Figure 2(c) shows the reconstructed population amplitudes of the groundX state of the N + 2 ion for four specific alignment angles versus the excursion time. The circles indicate the data extracted from the experiment for harmonics from H15 to H27 (including even and odd orders). Figure 2(d) shows the reconstructed relative phase between theX andÃ states. With these parameters, a complex-valued time-dependent wave packet from the two holesX andÃ can be constructed for We have also carried out the measurements for CO 2 molecules. Harmonic spectra of CO 2 9 involves three electronic states of the ion, the groundX (HOMO), the firstÃ (HOMO-1) and secondB (HOMO-2) excited states [29,30]. Using the same procedure as for N 2 , the reconstructed populations and relative phases of these three ion states are shown in Figs. 2(e,f) and Figs. 2(g,h), respectively. With these parameters, a complex-valued time-dependent wave packet from the three holesX,Ã andB can be constructed for CO + 2 . To evaluate our reconstructions, we have carried out calculations based on the time-dependent density functional theory (TDDFT) to simulate the parameters we obtained. The theoretical calculations are in "reasonable" agreement with the reconstructions (see Supplementary Note 5).

Filming attosecond charge migration
To visualize the hole dynamics, we have calculated the modulus square of the wave packet versus time for N + 2 and CO + 2 with the data in Figs. 2(c-h). To quantify attosecond charge migration processes, we first examine hole dynamics of CO + 2 when the fixed-in-space molecule is aligned parallel to the laser polarization axis. In Fig. 3(c), we plot the relevant molecular orbitals to help understanding. Figure 3(a) presents hole densities extracted at the recombination times of harmonics H15 to H26. Here, t=0 is when electron is born. In the figures, the frames are separated by steps of about 50 attoseconds. At a quick glance, we can see substantial hole migration from +x side to the -x side over the 0.5 fs duration. To take a closer look at how the hole migrates, in Fig. 3(b) we zoom in the region of x=[-1.2 a.u., 1.2 a.u.]. In the first 0.25 fs, we can see that the hole density coalesces in the middle region along the x-axis slowly from the +x side to the -x side, but in the second half, the pace turns faster where at later time we can see some density is moving back to the +x side. In addition, the migration seems to deplete density from the 4th→2nd→3rd→1st→2nd quadrants. In fact, its time evolution is more like a counterclockwise swirling than a direct migration. To observe the continuous evolution of hole density, we have constructed movies at steps of 10 attoseconds by interpolating the mixing coefficients in the complex wave packets shown in Figs. 2(e-h) (see Supplementary movie 1).
Since quantum mechanics can predict single hole behavior only statistically, it is inappropriate to associate such "movement" with trajectory or speed. It is more proper to think such movement like that of a flock of migrating birds. Each bird has its own mind and there is nothing to govern the trajectory of each bird, yet the flock will move on. In the presence of an external disturbance, the movement of the flock will change. So is the migration of an electron or a hole in the presence of the laser field. To describe the "movement" of charge density of a quantum system, consider the equation of continuity in quantum mechanics, ∂ρe ∂t +∇· J=0, where ρ e is the probability density and J is the probability current density. By calculating the latter quantum mechanically from the complex hole wave packet, we can then calculate the total flux crossing the x=0 a.u. plane at each time instant, as shown by the red line in Fig. 3(d). We have confirmed that the total flux change is equal to the rate of change of the total charge density. Likewise, the total flux crossing the y=0 a.u. plane can be calculated and shown as the green line. The total flux varies widely versus time, implying the statistical nature of the probability current density.
Our reconstruction method for each harmonic obtains the occupation amplitude and phase for all angles of each fixed-in-space molecules. In Figs. 4 (a-c), Comparing to results previously reported by Kraus et al. [18], we are able to make movies of charge migration at steps of 10 attoseconds by interpolating the experimental data, for molecules at all alignment angles. This is because we were able to retrieve single molecule complex dipole at each fixed-in-space angle and probe charge migration at recombination times for both even and odd harmonics at every 50 attoseconds, instead of 100 attoseconds if only odd harmonics are measured.
With twice as many data points within the same interval, it is more accurate to interpolate. Figure   3(b) has shown that such 10 attoseconds resolution is needed in order to follow the charge migration in CO 2 .

II. DISCUSSION AND CONCLUSION
In summary, we have used a ML-based HHS method to construct movies of charge migration in a molecular ion at its natural timescale of 10 attoseconds to follow the change of hole density during the time between electron is removed till it has recombined with the ion to emit harmonics. We presented charge migration at the most fundamental level for each single fixed-in-space molecule at all alignment angles from the experimental HHG spectra. The method presented here in principle can be extended to other molecules, even to more complex molecules. In reality, the method is limited by the complexity of the electron wave packet of the cations. If the electron wave packet composes of many molecular orbitals, more mixing parameters among them should be retrieved which in turn would require more experimental data. Thus, even though iodoacetylene studied in Kraus et al. [18] is a larger molecule, only two molecular orbitals are active in harmonic generation as in N 2 . For CO 2 , the wave packet is more complicated, thus charge migration as shown in Fig.   3, is more complicated and interesting.
In contrast to the conventional belief that attosecond electron dynamics should be probed with attosecond-pump-attosecond-probe scheme, HHS is the simplest route to probing attosecond electron dynamics involving valence electrons of molecules. Using HHS, a unique one-to-one timeenergy relation is established through constructive interference of harmonic emissions from a multicycle long driving laser. Both the technology and methods of analysis for studying attosecond charge migration are now available. Looking ahead, it is highly desirable to extend our method to other molecules, to explore how the charge migration depends on the structure of molecules, or more specifically, the spatial properties of active molecular cation orbitals that underlie the wave packet that governs attosecond electron dynamics.

III. METHODS
Experimental methods. Our experiment is carried out by using a commercial Ti:sapphire laser system (Legend Elite-Duo, Coherent, Inc.), which delivers 35-fs, 800-nm laser pulses at a repetition rate of 1 kHz. The output laser is split into a pump and a probe pulse. The pump pulse with moderate intensity is used to induce nonadiabatic alignment of molecules along its polarization.
The intense probe pulse has been used either directly (one-color scheme) or remoulded to a parallel two-color laser field (see Supplementary Note 1) to interact with the aligned molecules to generate high-order harmonics. The pump and probe pulses are parallel in the polarization. A motorized delay line is installed in the arm of the pump pulse to adjust the pump-probe delay. These two pulses are focused to a supersonic gas jet by a spherical mirror (f=250 mm). The gas jet is placed 2 mm after the laser focus to ensure good phase matching for short-trajectory harmonics, and the backing pressure is maintained at 0.8 bar. The generated high harmonics are detected by a homemade flat-field soft x-ray spectrometer. In the one-color experiment, the harmonic signals measured at different pump-probe delays [ Fig. 1(d) and Supplementary Fig. 2(e)] are used to identify the multiple orbitals effect in HHG. In the two-color experiment, the time-dependent harmonic signals are measured at various relative phases of the two-color laser field, providing a two-dimensional data set for decomposing the multiple orbital contributions in HHG process.
Reconstruction methods. Our reconstruction procedure has two steps. The first step is to retrieve single-molecule dipole moment from the measured time-dependent HHG signals. HHG from aligned molecular ensemble is expressed as Eq. (1). For a given harmonic order, Eq. (1) can be expanded as [41,42] S(τ ) = π 0 D(θ)ρ(θ, τ )sinθdθ * π 0 D(θ)ρ(θ, τ )sinθdθ = π 0 π 0 D * (θ 1 )D(θ 2 )ρ(θ 1 , τ )ρ(θ 2 , τ )sinθ 1 sinθ 2 dθ 1 dθ 2 . Let and discretize θ 1 , θ 2 in the range of [0, π] with a step of 0.01 rad, the S(τ ) then can be rewritten as Note that in Eq. (3), the imaginary part of D * (θ 1 )D(θ 2 ) is omitted due to its asymmetry upon the exchange of θ 1 and θ 2 , which will vanish after the convolution. In Eqs. (2) and (4), the 13 molecular axis distribution ρ(θ, τ ) and therefore ρ(θ 1 , θ 2 , τ ) can be determined from the one-color experiment with the method in [37]. To retrieve single-molecule dipole moment D(θ), we first solve the matrix R(θ 1 , θ 2 ) according to Eq. (5). Here, we utilize a widely-used ML algorithm, sparse representation [43], to do the retrieval. In the reconstruction, we first build a dictionary matrix for S(τ ) by expanding the R matrix with a series of two-dimensional Legendre polynomial The second step of the reconstruction is to decompose the multichannel contributions from the total dipole moment D(θ; α) obtained in the first step. In molecular HHG, the total dipole moment D(ω, θ) is a coherent superposition of the induced dipole moment of each emission channel, i.e., where D i (ω, θ) is the induced dipole moment of each emission channel. The subscript i denotes the involved molecular orbitals (molecular ion states) during the HHG process. Consider the initially populated molecular ion state i, D i (ω, θ) is given by where d i ion (ω, θ) and d i rec (ω, θ) are the transition dipole moments related to the ionization and recombination steps, respectively. a acc (ω) denotes the propagation amplitude of the electron wave packet in the continuum. d i ion (ω, θ) can be expressed as Here Ψ 0 (N ) is the ground state of the N-electron molecule, Φ i (N − 1) is the ground state or excited state of the (N-1)-electron molecular ion, where a molecular orbital ψ i has been removed to a continuum state χ k . Both many-electron wavefunctions Ψ 0 (N ) and Φ i (N − 1)χ k are properly antisymmetrized andD is the dipole operator from all the electrons. In obtaining Eq. 8(b), we assume that all the molecular orbitals are in the neutral and the ion does not change before and after ionization.
To calculate the recombination transition dipole, we take into account that the molecular ion has been modified by the laser field during the time interval between ionization and recombination Thus, the recombination dipole Inserting Eq. (10) to Eq. (7), we obtain Eq. (11) implies that hole hopping occurs before recombination. The degree of hopping depends on laser parameter and the alignment angle θ of the molecule. With Eq. (11), the total dipole moment for HHG can be expressed as whereD ij (ω, θ) = d i ion (ω, θ)a acc (ω)d j rec (ω, θ).
In our reconstruction, only the most relevant molecular orbitals ψ j are considered (that is,X and A states for N 2 , andX,Ã, andB states for CO 2 ). The dipole momentD ij (ω, θ) is calculated under the experimental laser conditions. The complex-valued coefficients C ij (θ) are directly associated to the electron dynamics in the molecular ion. To determine C ij (θ), we have performed two-color experiment. Since the SH field in the two-color experiment is weak enough and hardly alters the laser-induced electron dynamics, the coefficients C ij (θ) can be assumed to be independent of the relative phase α. The coefficients C ij (θ) are then retrieved from the dipole moment D(θ; α) obtained at different relative phases α by solving Eq. (12) with the genetic algorithm.
With the coefficients C ij (θ) retrieved, we can further calculate the wave function of the molecular ion at the recombination instant as 15 and the population coefficient of the molecular orbital ψ j as [18], where γ i (θ) is the initial population of the molecular ion state i, which is related to the alignmentangle-dependent ionization rate η i (θ) by |γ i (θ)| 2 = η i (θ)/Σ i η i (θ). In our work, the alignmentangle-dependent ionization rates of each molecular orbital are simulated by the MO-ADK theory [45,46]. Repeating the above procedure for different harmonic orders, we can then construct the hole dynamics in the ionized molecular ion according to time-frequency mapping underlying the HHG process.