Temporal and spatially resolved imaging of the correlated nuclear-electronic dynamics and of the ionized photoelectron in a coherently electronically highly excited vibrating LiH molecule

Few-cycle IR pulses are becoming available for exciting, ionizing, and probing the photoinduced subsequent dynamics in molecules.^1–3 They provide a versatile tool for investigating the correlation between electronic and early time nuclear dynamics both in the neutral and in the cation. Such ultrafast pulses may also enable the steering of nuclear motion in molecules toward a nonstatistical outcome⁴ that is typically not possible for nanosecond or even picosecond multiphoton excitation⁵ of either vibrational^6,7 or electronic^8,9 states.

Ultrafast electronic dynamics and charge migration have been experimentally evidenced using attosecond or few femtosecond pulses in various pump-probe setups¹ based on XUV attosecond and few femtosecond NIR pulses,^10–13 two XUV attosecond pulses,¹⁴ and high-harmonic generation,^15–18 or more recently using attosecond transient absorption spectroscopy of core electrons.^19–21 In addition, few-cycle carrier envelope phase controlled pulses^22–24 and deep UV attopulses²⁵ have been used to steer the outcome of chemical reactions.

Understanding the molecular response to such short pulses requires a quantum dynamical description of the correlated electronic and nuclear motions. When higher excited electronic states are accessed, this remains a challenge for quantum chemistry. We investigate the coherent electron-nuclear dynamics induced by few-cycle IR pulses at the quantum level all the way up to ionization. To integrate the time-dependent nuclear Schrödinger equation (TDSE), we used the partitioning technique^26,27 for the neutral and the cationic states coupled by their interaction with the electric field of the pulse. In addition, the nonadiabatic coupling (NAC) between the electronic states of the neutral driven by the nuclear motion is fully taken into account during and after the pulse. We show that such a level of description is necessary to capture in space and in time the phase correlations that shape the vibronic wave packets in the neutral and cationic spaces.

The molecular system studied is the four electron diatomic molecule LiH. LiH is a heteronuclear diatomic molecule with a dense manifold of low-lying excited electronic states of different polarities in the UV range that fragment into chemically different species. Its electronic structure and in particular the Coulomb ionic character Li⁺H⁻ of its ground electronic state has been extensively studied since the early days of quantum chemistry.^28–34 More recently, high level multiconfiguration electronic structure computations have been carried out.^35–39 The different polarity of the excited electronic states arises from the transfer of the ionic character to higher electronic states as the internuclear distances increases, which leads to a rich charge migration dynamics when a superposition of electronic states of different polarities is built by short femtosecond pulses. The LiH molecule is therefore an ideal prototype for investigating charge migration and dissociation control at the fully quantum level using short few femtosecond pulses.

Dynamical simulations of the photoexcitation of the LiH molecule using ultrashort pulses have been investigated previously by others^40–42 and in our group. We investigated the purely electronic quantum dynamical response to short femtosecond pulses in the neutral^43–46 or including the coupling to the photoionization continuum in pump probe schemes.^47–49 Coupled dynamics of electrons and nuclei in the photoexcited neutral was simulated using a nuclear grid.^38,50–52 The method that we describe in this paper takes the next step by simulating the coupled dynamics of electrons and nuclei, including photoionization and photoexcitation, thereby getting a better description of the field-induced dynamics including the amplitude exchange between the neutral and the continuum states when the molecule interacts with strong IR pulses.

Up to now, taking into account the role of nuclear motion in describing the photoionization process has been mostly limited to simulations relying on the sudden ionization approximation at a given (frozen) geometry of the nuclei.^53–60 The full quantum dynamics of molecular photoionization coupled with nuclear motion have been investigated in the case of single and two-electron diatomic molecules. The dynamics of the H₂⁺ molecular ion has been computed by solving the full 3-body TDSE.^61–67 The dynamics of the photoionization of two- and multi-electron diatomic molecules have been computed in the single active electron (SAE) approximation for H₂^25,68 and alkaline molecules.^69,70 Complete reviews can be found in Refs. 1 and 71–73. Our method allows describing the dynamics of the dipole interaction of a multi-electron molecule with an intense IR pulse, including the photoexcitation to several excited electronic states, their coupling to the ionization continuum in the single active electron approximation and the nonadiabatic interactions driven by the nuclear motion during and after the pulse. Furthermore, we propose an approach for computing the photoionization matrix elements at each nuclear geometry efficiently, which could otherwise represent a rate limiting step when several electronic states are coupled to the ionization continua. This approach can be extended to more than one nuclear degree of freedom.

In our numerical implementation, the multiphoton excitation and ionization processes induced by moderately strong few femtosecond pulses (with a peak intensity of about 10¹³ W/cm²) are treated nonperturbatively. The partitioning technique is used to describe the two orthogonal subspaces and the neutral and the ionized subspaces,^49,68,74 where we extend the methodology to two kinds of grids: The nuclear wave functions of the neutral and cationic states of LiH are represented on a 1 dimensional (1D) nuclear grid, and a 3D electronic grid is used to describe the wave function of the photoelectron. A summary of our computational scheme is given in Fig. 1 and in the Appendixes. Our methodology allows visualizing the spatial and temporal localization of the photoelectron and of the electronic densities and coherences during the pulse.

Solving the TDSE simultaneously on a 1D nuclear and a 3D electronic grid for coupled electronic states of the neutral and the cation becomes rapidly computationally very heavy because the photoionization matrix elements need to be computed at each nuclear geometry for a dense enough 3 dimensional representation of the electronic coordinate of the quasicontinuum, which substantially increases the dimensionality of the simulation. Therefore, several approximations are required to be able to account for the dynamics in the coupled neutral and cationic subspaces. Photoionization is described assuming a single active electron (SAE) so that the total wave function of the ionized subspace is an antisymmetrized product of the electronic wave function of the cation and the orthogonalized plane waves.^49,74–77 Within this approximation, the photoionization matrix elements reduce to computing the dipole between a Dyson orbital and the wave function of the photoelectron. Using orthogonalized plane waves as the basis set for the photoelectron allowed us to develop an efficient numerical approach for computing the dipole photoionization matrix elements in the reciprocal space, exploiting analytical properties of Fourier transform;^78–81 see  Appendix B. Numerical efficiency is critical here because these matrix elements need to be computed at each nuclear grid point.

The SAE approximation has been compared with correlated methods for describing photoelectron spectra resulting in either above threshold ionization^82–84 or strong field ionization.^85–88 In our treatment of photoionization, the static electron correlation between the electrons of the neutral and the cation is included via the Dyson orbitals computed from multideterminantal electronic wave functions which can comprise multiple excitation configurations of the neutral and the cation.^{47,76,77,89,90} Since both the ground state and excited state of the neutral and of the cation are multideterminantal wave functions, all the electronic states of the neutral can have a photoionization matrix element with the states of the cation.

For the moderate field strengths investigated here, the value of the Keldysh parameter,^91,92 $γ=Ip/2Up$⁠, where U_p = E²/4ω_L is the ponderomotive energy of the laser pulse and I_p is the ionization potential of the ground state of the system (in atomic unit), is larger than 1. In this regime, the photoionization process is nonadiabatic with respect to the lowering and raising of the ionization potential due to the interaction of the electric field and the photoionization matrix elements depend on the laser frequency. We therefore neglect the tunneling ionization process through the lowered Coulomb barrier in the computations reported below. This process becomes important for values of the Keldysh parameter <1, when the leaving electron follows adiabatically the oscillations of the electric field and the photoionization matrix elements no longer depend on the laser frequency.⁹³ 

Coulomb interactions between the departing photoelectron and cation lead to a breakdown of the SAE. They become important when the photoelectron is localized close to the cation core. They play a significant role for values of the Keldysh parameter <1 and can induce coupling between the open ionization channels built on different states of the cation.^56,94,95 Coulomb interactions are also responsible for autoionizing resonances embedded in the continua above the ionization potential and the Auger processes.^96–99 In static computations of photoionization cross sections, Coulomb interactions are included at various levels of theory using several kinds of basis sets for the wave functions of the cation and of the photoelectron.^1,100 Earlier approaches based on variational scattering theory¹⁰¹ and algebraic diagrammatic construction (ADC)^102,103 were developed for weak laser intensities and long picosecond and nanosecond pulses or synchrotron radiation. Recent approaches are based on describing the photoelectron by Coulomb waves basis sets,⁸⁹ B-spline,¹⁰⁴ B-spline combined with Gaussian orbitals,¹⁰⁵ B-spline combined with ADC,^106–108 and the R matrix theory.^109–111 For dynamical studies at frozen nuclei, mixed B spline-Gaussian basis set,¹⁰⁰ discrete variable representation,¹¹² and finite element representation of the radial coordinates¹¹³ have been used. When a representation of the photoelectron wave function is not needed, one can use very diffuse sets of atomic orbitals to describe the electronic structure of the neutral and absorbing potentials to compute the rate of photoionization.^114–117 One can also describe strong field ionization classically or semiclassically.^{114,118–123}

In the electronic structure computations, a large basis set of AO augmented with Rydberg orbitals is used which allows the description of electron correlation effects for the electronic states of the neutral, including high Rydberg states, and of the cation and at the level of the Dyson orbitals for the photoionization process. In the single active electron approximation used in the simulations below, the Coulomb interactions between the cation and the leaving photoelectron are not explicitly taken into account. The first excited state of the LiH cation is 12 eV above the ground state of the cation. Since the carrier frequency of the pulse is in the IR and the peak intensity is 10¹³ W/cm², a single cation state can be accessed and the interchannel coupling induced by Coulomb interactions is not expected to play a major role. The interaction between the photoelectron and the electric field is included by rotating the plane waves to Volkov states which makes the Hamiltonian diagonal in the continuous subspace.⁴⁹ Channel mixing due to the dipole interaction is thereby fully accounted for.

The paper is organized as follows. Section II provides details on the potential energy curves, dipoles, and nonadiabatic couplings of LiH, computed using a complete active space (CAS) average description. In Sec. III, we summarize the partitioning approach used to integrate the time-dependent Schrodinger equation with many more details given in  Appendix A. The derivation of the expression of the dipole photoionization matrix elements is given in  Appendix B. In Sec. IV, we show that the carrier envelope phases (CEPs) of a few-cycle intense IR pulse can be used to initiate different nonequilibrium electronic and vibrational dynamics in the neutral, leading to different asymptotic yields. We also show that photoionization acts as a filter to ionize preferentially certain photoexcited states, thereby providing an extra knob on the control. Because the dynamics is coherent, the phase of the pulse is imprinted in the photoelectron wave function as well as on the wave function in the bound states even when photoionization and nonadiabatic coupling during the pulse are allowed in the computation. In particular, we show that the value of the CEP governs the phase of the electronic coherences which strongly affects the amplitude exchanges when the vibronic wave packet reaches regions of nonadiabatic coupling. Guiding the dynamics through nonadiabatic interactions remains a challenging goal.^124–126 We report on snapshots of the electronic density as the vibronic wave packet and coherence go through the NAC region, thereby contributing to a better understanding of the role of the nonequilibrium electronic density in their control.

The potential energy, permanent and transition dipoles, and nonadiabatic couplings (NACs) are computed as a function of the interatomic distance, R, at the state averaged CAS-SCF (4, 20) level with the 6-311++G(3df, 3dp) Gaussian basis set augmented by 2 S, 3 P, and 3 D Rydberg orbitals centered on the H atom (with ζ = 0.010 95, 0.002 737 5, 0.046 875, 0.117 2, 0.00 293, 0.033 333, 0.011 111, and 0.037 00) for a band of 10 ∑ states, 4 doubly degenerate ∏ states, 1 Δ state for the neutral, and 4 ∑ states for the cation. As can be seen from Fig. 2(a), the potential energy curves possess an attractive Coulomb-like character, that is transferred from the GS to the higher excited electronic states as the internuclear distance increases. This feature of the potential energy curves is due to the nature of the electronic wave function that includes a Li⁺H⁻ electron transfer configuration.^28,30 Including Rydberg orbitals is essential for describing correctly the shift of the charge transfer character by NAC from the GS to high ∑ states as the internuclear distance increases, as can be seen from Fig. 2(b). The quantum chemistry software MOLPRO¹²⁷ has been used throughout. The subprogram MULTI is used for averaged CAS SCF.^128,129 The ground state of the cation and its first excited were computed at the same CAS level. The first excited state of the cation is 12 eV above the GS of the cation and therefore is not expected to play a role in the simulations described in Sec. IV.

The nonadiabatic coupling (NAC) matrix elements [Fig. 2(b)], τ_ij(R), between the adiabatic electronic states are computed using the DDR program^130–132 of MOLPRO.¹³³ The signs of the permanent [Fig. 2(c)] and transition dipole [Fig. 2(d)] moments and of the NAC were corrected by maximizing the overlap between the electronic wave functions between adjacent grid points along R. The electronic structure was computed for 512 grid points from which curves were interpolated for integrating the time dependent Schrödinger equation (TDSE). The molecular frame orientation is shown as an inset in Fig. 2 where the molecular axis is along the z axis with the Li atom in the +z direction. We report only on NAC and dipole coupling curves between the ∑ states because in all the simulations below, we use optical pulses polarized along the LiH molecular axis, which by the selection rules, can only photoexcite neutral ∑ states. The ∏ states could be accessed through the interaction between the bound states and the continuum, which are not subject to symmetry selection rules. However, for the field strength and the short pulses investigated below, the amplitudes in the ∏ states due to the photoinduced interactions with the continuum are negligible.

The NAC between $∑3$ and $∑4$ and to a lesser extent that between $∑2$ and $∑3$ are already large in the FC region and of opposite sign. There are also two regions of strong interaction between these two pairs of states at larger R values, between 8 and 14 Å. The permanent dipoles of the four lowest states [Fig. 2(c)] are of alternating polarity and large in the FC region: $∑3$ being Li^δ+H^δ− as the $∑1$⁠, which is opposite to the polarity Li^δ+H^δ− of the GS and $∑2$ and $∑4$⁠. Partial charge analysis of the low lying electronic states of LiH based on valence bond electronic structure computations can be found in Refs. 134–136. A natural bond orbital (NBO) computation of the partial charges of the Li and H atoms at the CAM-B3LYP/6-311++G(3df,3pd) density functional theory level for the ground electronic state gives −0.817 |e| on H and +0.817 |e| on the Li and −0.846 |e| sur H and +0.846 |e| sur Li at the HF level. Consequently, the GS has a large enough equilibrium permanent dipole [+2.24 a.u., see Fig. 2(c)] for LiH to be spatially oriented experimentally.^137–139 The signs of the transition dipoles between the lowest excited states reflect the polarity of the states they connect. The transition dipoles are also large in the FC region, see Fig. 2(d), which leads to a rich transient dynamics during the pulse.

The nuclear TDSE is solved numerically using the partitioning technique for the 10 ∑ states below the IP, including the NAC and dipole coupling for photoexcitation within the neutral ∑ manifold and also photoionization from these states to the ground state of the cation. A complete description of the methodology can be found in  Appendixes A and  B. The basis set for integrating the TSDE on a finite grid^140,141 is defined as follows.^52,142 Nuclear wave functions are orthogonal door functions centered on grid points, θ(R_g). The electronic wave functions of the neutral are adiabatic N electron wave functions, $Ψineutr;Rg$⁠, where i is the index of the electronic state and r stands for the 3N electronic coordinates. The total wave function of the neutral are products, $Ψineutr;RgθRg$⁠. The wave function of the ionized states is taken as an antisymmetrized product of an N − 1 electron wave function of the cation, $ΨGScatr;Rg$⁠, where r stands here for the 3(N − 1) electronic coordinates, and the wave function of the ionized electron, $ϕk⊥elecrl$ indexed by l, with momentum k. We discretize the continuum using a basis of orthogonalized plane waves, $ϕk⊥elecrl$⁠, Eq. (B3). This basis set allows for efficient computation of the photoionization matrix elements, as is described in the Appendixes. The plane waves are orthogonalized to all the molecular orbitals (MOs) included in the CAS active space of the neutral^47,49 so that we can define two orthogonal subspaces the bound subspace: $Q=∑i,gNneutNgΨineutr;RgθRgΨineutr;RgθRg$ and the ionization subspace $P=∑j,g,kNcat,Ng,NkΨjcatr;RgθRgϕk⊥elecrlΨjcatr;RgθRgϕk⊥elecrl$⁠, where r_l is the coordinate of the single active electron. Using these two subspaces, the partitioning of the nuclear TDSE [see  Appendix A, Eqs. (A1)–(A3)] leads to two sets of coupled equations,

with the total wave function projected on two subspaces,

where $cigneutt$ is the amplitude of the bound electronic state at grid point g and $cjgkneutt$ is the amplitude of an electronic state j of the cation at grid point g with a momentum k of the photoelectron. In Eqs. (1) and (2), N_neut is the number of electronic states of the neutral, N_cat is the number of electronic states of the cation, N_k is the number of plane waves, and N_g is the number of grid points along R. In the simulations below, N_neut = 10, N_cat = 1, N_k = 28 672 and N_g = 256. The number of nuclear grid points, N_g, and k values for discretizing the plane waves, N_k, for each grid point in R has been checked for convergence. For each R value, 56 values of momentum |k| (ranging linearly from 0.0001 to 1.5 a.u.) and 512 angular distributions for each value of momentum (θ and ϕ values are randomly sampled within a single quadrant of the unit sphere, and the distribution are symmetrized around the origin to ensure the symmetry of the photoionization coupling elements) are included. In total, we have 7 342 592 basis set functions. See also Fig. 1 above.

The nuclear kinetic energy is computed using the finite difference approximation [See Eq. (A6) in  Appendix A] which allows us to maintain realistic numerical precision at a reasonable computer cost.^52,143 Detailed expressions of the Hamiltonian matrix elements in Eq. (1) are given in  Appendixes A and  B, and its structure is schematically shown in Fig. 1(b).

From the electronic structure computation at each R value, one can compute isocontours of the electronic densities and of the transition densities between them. Such a plot is shown in Fig. 3 for the NAC region between $∑2$ and $∑3$ just outside of the FC region. Because these two states are also coupled to $∑4$⁠, see Fig. 2(b), the picture is more complex than a simple two electronic state interaction, but, nevertheless, one can clearly see the charge transfer switching between the two states occurring between the two sides of the avoided crossing. At 2.2 Å, the two states still have a parentage with their character at $Req:∑2$ has a Li^δ+H^δ− character, while $∑3$ has the opposite one, Li^δ−H^δ+. At 2.8 Å, where the NAC is maximum, the two states have a mixed character. At the exit of the NAC region, at 3.3 Å, a charge transfer between the two adiabatic states has occurred: the density of $∑2$ has an excess of charge on the Li, as that of $∑3$ before the crossing and $∑3$ has more density on the H atom.

The diagonal and transition densities are key ingredients that define the nonequilibrium electronic density in the neutral built by the pulse $ρ(r,t;Rg)$ and its time evolution as the nuclear dynamics unfolds. Of special interest is the coherence between two electronic states in going through a NAC region. The charge transfer character of the electronic coherence is governed by the transition densities. The time-dependence of $ρ(r,t;Rg)$ is governed by $cigneutt$ computed by integration of the TDSE. The expression of the nonequilibrium electronic density is given in Eq. (A13). We focus in Sec. IV on the coherence between $∑2$ and $∑3$ shown in Fig. 3.

In this section, we examine interplay between nuclear motion and the photoexcitation and photoionization dynamics during the pulse. We follow the evolution of the vibronic wave packet built by the pulse during the subsequent dynamics that involves nonadiabatic coupling. These couplings can both change the phase relation between the components of the vibronic wave packet and the asymptotic population in each adiabatic electronic state.

The GS of the aligned LiH molecule is excited by two CEP controlled few-cycle IR pulses linearly polarized along the molecular axis, z; see the inset in Fig. 2. The pulses have a carrier wavelength of 720 nm (ω = 0.063 a.u. = 1.17 eV), a FWHM of 3.5 fs (σ = 62 a.u.), and a peak intensity of 1.14 10¹³ W/cm² (f₀ = 0.017 a.u.) [see Eq. (A5) in  Appendix A]. The time profile of the electric field has one major maximum and two secondary extrema of opposite sign in its envelope; see Fig. 4 for the electric field profiles. As in previous works,^3,38,45,46 we use CEP values differing by π to steer the motion of the electronic density and the subsequent vibronic dynamics. For a CEP value of 0 in Eq. (A5) in  Appendix A, at its maximal strength, the electric field points in the +z direction, toward the Li nucleus, while for the CEP value of π, the electrical field at its maximum points in the −z direction toward the H nucleus; see Fig. 4.

Due to its short duration, the pulse has a large bandwidth of 1.03 eV. Because of this large bandwidth, the pulse coherently excites a superposition of a large number of vibrational states in all the excited electronic states that are accessed. The potential energy curves shown in Fig. 2(a) are very anharmonic. The harmonic vibrational period of the GS is 26 fs. The one of the Σ₁ state is 90 fs. The higher excited electronic states are even more anharmonic and shallow and exhibit secondary minima. The Σ₁ state [see Fig. 2(a)] can be accessed from the ground state by a two photon transition. The continuum can be reached sequentially through the photoexcitation of higher states (Σ₂–Σ₉) that have large transition dipole moments with Σ₁ (and also between themselves) and large photoionization cross sections. The Keldysh parameter of the ground and five lowest excited states is above 1, suggesting that multiphoton photoionization should be the dominant process.

Since the dynamics is coherent in the coupled bound and continuous subspaces, the CEP control is maintained in the presence of significant photoionization, photoexcitation, and nuclear motion during the pulse. In particular, as shown in Fig. 4, the two CEP values allow achieving different superpositions of bound electronic states at the end of the pulses, thereby building different nonequilibrium electronic and vibrational densities that evolve to different asymptotic yields; see Table I for the electronic states of the fragments.

As can be seen from Fig. 4, photoexcitation occurs throughout the duration of the pulse, while significant photoionization only begins to take place when the electric field reaches its major maximum. Overall, the two CEP values lead to about 20%–25% of ionization and this yield is slightly higher for moving nuclei. By the end of the pulse, the yields in the bound excited electronic states are of a few percent. One can see that the populations in the excited sates computed with and without nuclear motion differ significantly. The role of the nuclear motion during the pulse is twofold. The spatial extension of the nuclear wave function of the GS combined with the variation of the potentials and permanent and transition dipoles over the Franck Condon region leads to a different photoexcitation and photoionization dynamics. In addition, in the second half of the pulse, the field free states are coupled by both the dipole coupling and the NAC, which leads to complex pattern of the transient amplitude transfer. The strong NAC between field free $∑3$ and $∑4$ states and the onset of the NAC between $∑3$ and $∑2$ at the exit of the FC region [see Fig. 2(b)] lead to a significant reorganization of the amplitudes and populations in the bound states toward the end of the pulse.

During the first half of the pulse, the CEP value controls the extent of photoexcitation of the low ∑ states and the population dynamics computed for moving and frozen nuclei are similar. The $∑1$ and $∑3$ states have a negative polarity (negative permanent dipole moment) in the FC region and are primarily accessed when the electric field is negative, while the GS, $∑2$⁠, and $∑4$ have a positive one, see Fig. 2(c), and are primarily accessed when the field is positive. This polarity control leads to a significant excitation of $∑1$ during the first negative extremum before the positive maximum of the electric field for CEP = 0, both for moving and frozen nuclei, as shown in Figs. 4(a) and 4(b). In contrast, $∑1$ is only populated when the field reaches its negative maximum for CEP = π, as shown in Figs. 4(c) and 4(d). As a consequence, for the CEP = 0 dynamics, the extensive photoionization of $∑1$ that takes place at the positive maximum of the electric field prevents the photoexcitation of the higher ∑ states in the second half of the pulse which are all below 3%, $∑2,∑3$⁠, and $∑4$ being populated by similar amounts. On the other hand, for the CEP = π dynamics, the yield in $∑2$ reaches about 8.3% at the end of the pulse and is larger than the population in $∑1$ (4.3%) and of $∑9$ (2.2%), while all the other excited states have populations below 1%. The populations of these states are strongly affected by the NAC induced by nuclear motion at the exit of the Franck-Condon region. Therefore, the two CEP values lead to significantly different yields in the low excited states $∑2,∑3$⁠, and $∑4$ at the end of the pulse (at 15 fs for CEP = 0: $∑2=0.7%,∑3=2.0%$⁠, and $∑4=2.3%$ and for CEP = π $∑2=8%,∑3=0.4%$⁠, and $∑4=0.4%$⁠).

As can be seen from Fig. 5, the populations in the $∑2,∑3$⁠, and $∑4$ states are also extensively reorganized by the NACs at later times. For CEP = 0, $∑2,∑3$ strongly interact in the range 20–40 fs where they undergo several amplitude exchanges. In the range 80–120 fs, the interaction is smaller and the populations in $∑2$ and $∑3$ oscillate but do not cross; see Fig. 5(a). For a CEP = π, there is only one strong interaction between $∑2$ and $∑3$ taking place around 30 fs, as shown in Fig. 5(b), accompanied by 3 oscillations in the populations without crossing. Note also that there is no significant population in $∑4$ after the pulse for CEP = π, in contrast to CEP = 0 and that the yields in $∑1$ and $∑9$ are similar for both CEP values and essentially unaffected by nonadiabatic coupling after the pulse. The different asymptotes of the excited states of LiH correspond to different excited states of the Li atom (see Table I) and can be measured. In the case of the CEP = 0 dynamics, we get that $Σ1$ (3.4%) is about equal to $Σ4$ (3.2%) and larger than $Σ3$ (1.1%) and $Σ2$ (0.6%). The relative yields are essentially reversed for the CEP = π dynamics [$Σ2$ (8%) > $Σ1$ (4%) > $Σ3$ and $Σ4$ (0.6%)].

The two time ranges of population oscillations between $∑2,∑3$⁠, and $∑4$ correspond to the NAC occurring in the regions 2–4 Å and 8–14 Å in Fig. 2(b); see also Fig. 3 that shows isocontours of the electronic diagonal and transition densities in the 2–4 Å NAC region for $∑2$ and $∑3$⁠. Because the two values of the CEP lead to a different localization and spreading of the wave packets on the different electronic states at the end of the pulse, these are affected differently when they reach the NAC regions. This effect can be seen from the maps of the $∑2−∑3$ coherence plotted as a function of the interatomic distance (x axis) and of time (y axis) in Fig. 6 where the two NAC regions are marked by a green square. The mean value of the position of the wave packet on $∑2$ is shown as a plain green line, while that on $∑3$ is shown in dashes. Vertical arrows connect the same time and same R values on the two maps. One can see that at the end of the pulse (first vertical arrow), the $∑2−∑3$ coherence for the CEP = 0 pulse has an opposite phase compared that of the CEP = π one. In the NAC regions, the wave packet exhibits different branches depending on its spreading in R (and correspondingly in momentum, not shown) which is different for the two CEP values. These branches subsequently acquire different phases. After the first NAC region, the $∑2−∑3$ coherence exhibits two branches for both CEP values. For the top branch (second vertical arrow), the two coherences are now in phase, while for the bottom one (third vertical arrow), they remain out of phase. After the second NAC region, the branching pattern is more complex for the CEP = π dynamics than for the CEP = 0 one. The two main branches within each coherence are out of phase, the lowest branch of the coherence remains out of phase for the CEP = 0 and the CEP = π dynamics, and the most upper branch remains in phase.

One also sees that while the mean R values on $∑2$ and $∑3$ are similar at short times, they reach larger values for the dynamics with the CEP = 0 than the CEP = π. Because of the similar localization of the wave packet at the early stage (20–40 fs) of both dynamics, the first NAC region plays an important role for both [see Figs. 5(a) and 5(b) for the population exchange]. In the 80–120 fs region, the effect of the NAC is small for the two CEP values and smaller for the CEP = π dynamics than for the CEP = 0 one. The reason is that the wave packets on the two potential curves are less localized in the second NAC region than in the first one and less for the CEP = π dynamics than for the CEP = 0 one, which leads to very small population exchange as shown in Fig. 5. Another reason contributing to the small effect of the second NAC interaction in the case of the CEP = π dynamics is that the population in $∑3$ is considerably smaller than in $∑2$ in the CEP = π dynamics, while the two populations are more commensurate for the CEP = 0 one. These results show that the CEP value affects the population and the phase at the end of the pulse, but also the ensuing nonadiabatic dynamics because it determines the initial conditions for the motion of the nuclear wave packets on each potential.

We next illustrate in Fig. 7 how the nonequilibrium electronic density, $ρ(r,t;Rg)$ [Eq. (A13)], evolves when the $∑2−∑3$ coherence passes through the first NAC region for the CEP = 0 dynamics [lowest green square, Fig. 6(a)]. To understand the evolution of the nonequilibrium densities, the static picture shown in Fig. 3 needs to be convoluted with the fact that the time-dependent electronic densities are a superposition of electronic states with amplitudes that are determined by the solutions of the TDSE [Eq. (1)]. In Fig. 7, we show isocontours of the electronic density computed at four values of (indicated by yellow dots) along the mean R values of wavepackets traveling on and the curves. The isocontours of the electronic density are plotted on a heatmap of the coherence corresponding to a zoom in the first green square of Fig. 6(a). Four states $∑1,∑2,∑3$⁠, and $∑9$ contribute to the electronic density [see Fig. 5(b)], but only $∑2$ and $∑3$ contribute to the coherent beating shown in the heatmap so that a two electronic state picture provides a good semiquantitative description. Such a two-state description has been used before to provide understanding on the control of charge migration in neutrals and cations and its control with the pulse parameters.^{23,24,38,43,44,67,124,144–149} Using the simple picture of a coherent superposition of only $∑2$ and $∑3$⁠, at a given time and R value, the electronic density takes the form of a sum of two terms, a population term and an interference term [see Eq. (3) below]. The population term is the sum of the electronic densities at this R value weighted by the populations at this time and R value. The interference term depends on the electronic transition density matrix, $ρ∑2−∑3r;R$⁠, which is weighted by time-dependent amplitude of the coherence between $∑2$ and $∑3$⁠, that we can write $2cg∑2tcg∑3tcos Δϕt$ with $Δϕt=ϕg∑2t−ϕg∑3t$

When Δϕ(t) is 0 or 2nπ, the interference term is added to the population term, while when Δϕ(t) is (2n + 1)π, it is subtracted. Note that because of the NAC, the electronic state and transition densities change character when going through the avoided crossing region, as shown in Fig. 3. Starting from the left in Fig. 7, at the first yellow (R, t) point, the populations in $∑2$ and $∑3$ are 0.5% and 2%, respectively. The coherence term is negative, meaning that the interference term was subtracted to the stationary part in Eq. (3). The interference term adds density on the top of the Li and between H and Li. At the next point, the coherence is positive (red color) and the position is that at which the NAC is maximum (R_g = 2.66 Å). The populations in the two states are about equal (0.8% and 0.6%, respectively). The electronic density is essentially an average between the two stationary densities of the adiabatic states shown in Fig. 3. At the third point, the populations in $∑2$ and $∑3$ are about equal (1%) and the interference term needs to be subtracted, which leads to an excess of density at the H atom and a contour closer to that of $∑3$ after the crossing. At the fourth point, the population in $∑2$ is 0.5% and smaller than in $∑3$ (1%). The interference term which is positive contributes to reinforce the $∑2$ character, meaning that the polarity of the electronic density has been inverted as the coherence went through the NAC region.

The CEP control achieved in the presence of photoionization and NAC is imprinted in the oscillations of the dipole moment in the neutral bound states [Eq. (A10) of  Appendix A] plotted in Fig. 8. The oscillations of the dipole moment reflect the coherent charge migration and the localization of the electronic density on the Li and the H nuclei⁴⁶ and could be probed experimentally by time resolved photoelectron angular distributions^{14,47,60,147,150–155} or transient absorption.^50,51 In Fig. 8, the fast oscillation of one femtosecond corresponds to the period of the $GS-∑1$⁠. One clearly sees from the inset that this coherence oscillates with a phase difference of π for the two CEP dynamics during the first 60 fs. This fast oscillation is dumped as the $∑1$ wave packet travels on the shallow $∑1$ potential and stops to overlap with the wave packet on the tight GS potential [see Fig. 2(a)]. The vibrational period on $∑1$ is 90 fs, which leads to a revival of the $GS−∑1$ coherence at 90 fs and 180 fs. At these two revivals, the oscillations of the $GS−∑1$ coherences are in phase. For the CEP = 0 case, one can also clearly see in the inset the beatings of the 3 coherences (⁠$∑2−∑3,∑3−∑4$⁠, and $∑2−∑4$⁠) which take place on a 5–8 fs time scale and a faster oscillation of 2 fs that corresponds to the $Σ1−Σ2$ coherence. Such beatings are not present for the CEP = π dynamics at longer times because the populations in $∑3$ and $∑4$ are much lower than that in $∑2$⁠. For the CEP = π dynamics, after 50 fs, the dipole value is dominated by the contributions of the GS and of $∑2$ which have positive permanent dipoles, which leads to comparatively higher values than for the CEP = 0 dynamics. All the $Σ$ states have very shallow wells and are not purely repulsive, see Fig. 2(a), which leads to very long lived vibronic coherences that persist as long as the wave packets on two different electronic states overlap, as shown in the case of the $Σ2−Σ3$ above. This case is different from the situation of a superposition between a bound and a purely repulsive state for which the overlap is fading as the wave packet is moving on the repulsive state, or for the coherent excitation of two purely repulsive states. Such a situation has been reported is for H₂⁺ and T₂⁺ in Ref. 156. For LiH, such a case will occur when the polarization of the pulse allows accessing both $Σ$ and $Π$ states.⁵¹ The vibrational periods that reflect the motion of the vibrational wave packet on each electronic state can be probed by transient absorption spectroscopy.^50,51

We compare in Fig. 9 the dipole moment of the neutral during the pulse to that of the velocity dipole of the photoelectron [Eq. (A12) of  Appendix A]. In Fig. 9, the dipole of the photoelectron is computed in the velocity gauge and the dipole of the neutral molecule in the length gauge [Eq. (A10)]. It has been shown that for complete basis sets, the velocity gauge dipole and the length gauge dipole are equivalent.^157–160 The velocity dipole is used because it can be computed analytically [see Eq. (A12)], which leads to a more accurate value than a numerical integration as is required for the length dipole. For the moderate peak power of the pulse used in the simulation (10¹³ W/cm²), the value of the dipole at the end of the pulse is close to zero. In this case, it was shown that the dipoles computed in the two gauges lead to similar spectra.^158–160 In agreement with the short time inset in Fig. 8, Fig. 9 shows that in the first half of the pulse, when photoexcitation dominates, the dipole moment in the neutral follows the electric field essentially adiabatically (with a phase difference in the response of about π/2 as expected) and that the momentum of the photoelectron is essentially 0. At the maximum of the pulse, extensive photoionization takes place, which leads to a rise in the population of the cation [see Figs. 4(a) and 4(c)]. Note first that the dipole moment of the neutral and that of the photoelectron are in phase and that the dipole moment of the photoelectron computed at frozen nuclei is systematically smaller than that of the photoelectron because the ionization yield is smaller in that case [see Figs. 4(a) and 4(c)]. After the maximum of the pulse, the dipole moment of the neutral decreases, because of the rise of the cation yield, and is less adiabatic with respect to the field because of the onset of coherent beatings between the electronic states of the neutral analyzed above. On the other hand, the dipole of the photoelectron remains essentially adiabatic with respect to the field but exhibits small subhalf cycle oscillations. Because of the CEP difference between the two dynamics, ionization occurs from different sides of the molecule. For CEP = 0, there is a high yield of photoexcitation to $∑1$ that ionizes before the maximum of the pulse. At the maximum, the population is transiently high in $∑2$ which like the GS has an excess of electrons at the H side (see Fig. 10) which leads to a maximum of the dipoles of both the neutral and the photoelectron in the +z direction. Note that the dipole in the bound states then starts to oscillate nonadiabatically with respect to the field because other ∑ states are transiently populated. For the CEP = π dynamics, the excess of electronic density (see Fig. 10) at the maximum of the field is on the Li side (in agreement with the fact that there is a high population in $∑1$⁠), and the maximum of the dipole moments occurs in the −z direction. At the end of the pulse, the dipole moment of the photoelectron including nuclei motion for the CEP = 0 dynamics is negative, while it is slightly positive for the CEP = π dynamics.

The corresponding localization of the photoelectron [Eq. (A11)] and of the electronic density of the neutral [Eq. (A13)] during the pulse is shown in Fig. 10, panels (a) and (c) for the photoelectron and (b) and (d) for the electronic density, for the two values of the CEP. To highlight the variations in the localization of the electronic density, Figs. 10(b) and 10(d) show the difference between the total electronic density at time t and the electronic density of the GS. CEP = 0 is shown in the top row and CEP = π in the bottom one. When the field is negative, the electronic density is localized on the Li and when the field is positive, on the H. This can be seen from the isocontours of the electronic density at time t = 9.4 fs, for which electrons still follow the pulse essentially adiabatically. For CEP = 0 [panel (b)], the electric field is negative and there is an excess of electronic (yellow color) density on the Li and a defect on the H nucleus compared to the GS of the neutral in agreement with the large population in $∑1$ [Fig. 4(c)]. For CEP = π [panel (d)], the field is positive and there is a defect of density compared to the neutral (blue color) on the Li and an excess on the H nucleus (corresponding to higher population in $∑2$⁠). The localization of the photoelectron exhibits the same trend even though there has not been yet a significant photoionization. At time t = 9.4 fs, there is more density localized close to the Li for CEP = 0 [panel (a)] than for CEP = π [panel (c)].

The time t = 10.3 fs falls just after the maximum of the field. For CEP = 0, at the maximum of the electric field, a large fraction of the population photoexcited to $∑1$ has been ionized and $∑2$ is transiently populated [see Fig. 4(a)]. $∑2$ has the same polarity than the GS with an excess of electronic density on the H atom. The isocontour of the difference with the GS shown in panel (c) reflects the excited configurations present in $∑2$⁠. Since for CEP = 0, the maximum of field is in the +z direction which corresponds to an ionization from the H side, there is a higher localization of the photoelectron at the H side just after the maximum [panel (a)]. For CEP = π, the field is negative at its maximum and the ionization occurs primarily from the Li side, in agreement with a larger localization of the photoelectron probability density at the Li [panel (c)] and an excess of electron for the electronic density on the Li side [panel (d)]. The localization of the photoelectron on the H side for a positive extremum of the field is also very clear at t = 11.5 fs for CEP = 0 [panel (a)], while there is much less probability density at the H atom at the same time for CEP = π since the field has now a negative extremum. Overall, in the second half of the pulse, there is more density probability on the H side for the CEP = 0 dynamics than for the CEP = π one, and one clearly sees the photoelectron leaving the cation core, in agreement with the slightly positive and negative values of the dipole of the photoelectron shown in Fig. 9. One can also discern in Figs. 10(a) and 10(c) subhalf cycle features in the localization of the photoelectron; see, in particular, Fig. 10(a). These features translate in very small amplitude subhalf cycle oscillations in Fig. 9. They reflect the partial return of the photoelectron to Li or H side. Similar subcycle electron localization during a strong field process have been reported in the case of H₂⁺.¹⁶¹ Note that the isocontour of the bound density difference in the neutral at 11.5 fs in panel (d) looks similar to that at 10.3 fs in panel (c) and also reflects the high population in the $∑2$ state; see Fig. 4(c).

We reported on the coupled electronic-nuclear dynamics in the neutral and cation of the LiH molecule driven by few-cycle IR pulses of controlled CEP. Our methodology is fully quantum and includes photoexcitation, photoionization, and the NAC driven by nuclear motion during and after the pulse. To get an efficient time propagation, we developed a method of computation of photoionization coupling elements at every grid point in R based on the properties of the Fourier transforms of atomic orbitals. Our approach relies on the single active electron approximation which makes it tractable for including nuclear motion. The SAE neglects dynamical electron-electron correlation between the departing electron and the electrons remaining on the cation. This approximation could affect the value of the relative populations in the ionization continuum and in the excited electronic states at the end of the pulse but does not affect our main conclusion. Even in the presence of photoexcitation, photoionization, and NAC during the pulse, the value of the CEP provides a control on the dynamics in the neutral. The nonequilibrium vibronic wave packets built at the end of the CEP = 0 and CEP = π pulses in the neutral have different populations on the excited ∑ states and electronic coherences of opposite phases. They reach the NAC regions at larger internuclear distances with different spatial and momentum distributions, which leads to different amplitude transfers and eventually different asymptotic fragment yields. Therefore, such pump-probe schemes could be used to steer the charge migration of ionized molecules as well as to control the relative product yields. Our methodology to solve the TDSE is based on the partitioning technique and provides a wave function for both the states of the neutral and the states of the cation, for the 1D nuclear degree of freedom as well as for the 3D electronic wave function of the photoelectron. This allowed us to shed light on the correlations in space and in time of the entangled nuclear and electronic degrees of freedom during the pulse and also after the pulse when the wave packets passed through NAC regions.

This work was supported by the AMOS program within the Chemical Sciences, Geosciences, and Biosciences Division of the Office of Basic Energy Sciences, Office of Science, US Department of Energy, Award No. DE-SC0012628, and by the F.R.S.-FNRS research grants (Nos. T.0132.16 and J.0012.18). F.R. acknowledges support from the Fonds National de la Recherche Scientifique, Belgium (F.R.S.-FNRS). S.v.d.W. acknowledges the support of a Ph.D. fellowship from FRIA-F.R.S.-FNRS, and B.M. acknowledges the support from the University of Liège for a post doctoral fellowship. Computational resources have been provided by the Consortium des Equipements de Calcul Intensif (CECI), funded by the F.R.S.-FNRS under Grant No. 2.5020.11.

The set of coupled equation for the amplitudes of the electronic states at each grid point given in Eq. (1) can be derived from the TDSE for the total wave function |Φ(R, r, t)⟩,

where R is the nuclear coordinate and r stands for the electronic coordinates.

Using the projection operators defined in Sec. III and using the closure relation P + Q = 1, we can partition the TDSE into its neutral and cation components,²⁷ 

As explained in the main text, the total wave function Φ(R, r, t) is expanded in a basis of product of electronic states for the cation and on a grid for the nuclear coordinates,

Projecting the first and the second equation in (A2) on each product basis function leads to Eqs. (1) and (2) of the main text.

In Eq. (1) of the main text, the matrix elements in the bound subspace, $Hig,i′g′t$⁠, take into account the dipole coupling leading to photoexcitation and the NAC,

where μ is the reduced mass of LiH. The electric field of the pulse, E(t), is defined from the time-derivative of the vector potential to ensure that the time integral of the electric field is zero,⁶⁶ 

where ω is the carrier frequency, f₀ is the field strength, σ is the pulse duration (FWMH = 2.35σ), ε is the polarization vector, and ϕ is the carrier envelope phase (CEP), i.e., the phase between the electric field of the pulse and the Gaussian envelope.

In the nuclear kinetic energy term and in the nonadiabatic coupling term, the first and second derivatives were computed using the finite difference approximation to the fourth order,

which allows us to keep numerical precision at a reasonable computer cost.⁵² 

The Hamiltonian of the cation is written as separable between one electron that ionizes and the remaining cation, neglecting the electronic correlation between the N − 1 electrons of the cation and the active electron l,

As discussed in the Introduction, this approximation is warranted in the case when the departing electron is far from the cation core and the field strength and the carrier frequency of the pulse are such that photoionization dominates over strong field and tunnel ionization, which is the case for the simulations reported in this work. The electronic states $Ψjcatr,Rg$ are adiabatic states of $Ĥcatelec$⁠, which leads for the matrix elements, $Hjgk,j′g′k′t$⁠, in the ionization subspace to

where the orthogonalized plane waves have been rotated so as to diagonalize the dipole coupling in the Hamiltonian of the photoelectron.^47,49,74 Note that both the permanent and the transition electronic dipole coupling elements are included in Eqs. (A4) and (A8). The permanent dipole terms for the neutral and the cation result in the AC shift of the IP with the oscillations of the electric field, sometimes called the AC or the dynamical Stark shift. Since the TDSE is solved numerically, the effect is included to all orders. The first excited state of the cation is 12 eV above the GS of the cation and cannot be reached by multiphoton IR ionization. Since there is only one state of the cation included in the basis set, we do not explicitly write the NAC between electronic states for the P subspace.

For orthogonalized planes waves,^47,49,77 the off diagonal matrix elements, $Hig,jg′kt$⁠, leading to photoionization are given by

where $ϕijDrl;Rg$ is the Dyson orbital between the electronic state of the neutral i and the state of the cation j. Since these matrix elements need to be computed at each grid point Rg, we developed an original method, based on the analytical evaluation of the Fourier transform of the Dyson orbitals. As the Dyson orbitals can be expressed as linear combination of atomic orbitals, their dipole moment with a basis of orthogonalized plane waves can be computed readily in Fourier space. All the details of the numerical implementation of this approach are given  Appendix B.

Using Eqs. (1) and (A3), the time-dependent dipole moment of the neutral is given by

The probability density of the photoelectron along the molecular axis is defined as

where ρ(|k|) = k² is the density of states. The velocity dipole of the photoelectron is only nonzero along the molecular axis for a pulse linearly polarized along z and an oriented molecule,

The time- and geometry-dependent electronic densities ρ(r; R_g, t) shown in Figs. 6 and 9 are computed as the total electronic densities at a given grid R_g and instant t. The densities are defined for the superposition of adiabatic electronic states at this (R_g, t) point,

where the diagonal terms, ρ_II(r; R_g), are the one electron density matrix elements of the adiabatic electronic states of the neutral and the off diagonal terms, ρ_IJ(r; R_g), are the transition density matrix elements between the adiabatic electronic states I and J. Isocontours of ρ_II(r; R_g) and ρ_IJ(r; R_g) across the first NAC are shown in Fig. 3 for the states $Σ2$ and $Σ3$⁠.

The diagonal and transition density matrix elements are obtained from the electronic structure computation, fully accounting for the multideterminantal nature of the electronic wave functions.

Using plane waves as a basis set for photoelectron brings the photoionization matrix elements [Eq. (A9)] to Fourier transform (FT) which allows for numerically efficient numerical implementations of their computations. This is crucial because they need to be computed for each position of the nuclei. In this section, we detail the main steps of our implementation. The photoionization matrix elements need to be computed for each value of the nuclear coordinate, R_g. In the equations below, we do not write explicitly this dependence and r stands for the coordinate of the photoelectron.

In Eq. (A9), the integral is the transition dipole between a Dyson orbital and an orthogonalized plane wave,

The Dyson orbital, $ϕijDr$⁠, is the overlap between the N electron many determinantal wave function of the neutral and the N − 1 one of the cation computed at the nuclear geometry,

$ϕk⊥elecr$ are orthogonalized by the Gram-Schmidt procedure and normalized to the volume (L³ for a cubic box) which are given in direct space by

where $ϕmMOr$ is computed at each grid point R_g and $ϕ̃mMOk$ is its Fourier transform. The Gram-Schmidt procedure conserves the norm of the wave function, provided that the size of the box, L, is large enough.

The orthogonalization of the plane waves implies that we have a new set of orthogonalized plane waves and Dyson orbitals at each position of the nuclei, R_g. Using Eq. (B3), and writing the transition dipole, Eq. (B1), in reciprocal space, we get

where $ϕ̃ijDk$ is the FT of the Dyson orbital. The Dyson orbital can be expanded into MOs,

which leads to

where $μmnMO$ are the transition dipoles of the MOs at the grid point R_g.

The next step is to use the expansion of the MOs in terms of the atomic orbitals, AOs,

Using this expansion, we arrive at an expression for the photoionization transition moment, Eq. (B6), that depends only on LCAO and Dyson orbital coefficients and the FT of the AOs, $χ̃sAOk$⁠, and the gradient of their FT,

Since AO basis sets are usually expressed in terms of Gaussian Type Orbitals (GTOs), written as linear combinations of contractions¹⁶² as follows

where r_0r is the position vector of the nucleus to which the $χsAO$ basis function is attached. $Ylsmsr−r0sr−r0s$ are the real spherical harmonic used in Gaussian AO basis sets.¹⁶² α_ps are normalized contraction coefficients that are tabulated for each type of Gaussian basis set.

Equation (B8) can be rewritten in terms of the FT of the GTOs, whose analytical expressions can be derived using the spherical harmonics expansion of plane waves. The Fourier transform of the contractions in (B9) is equivalent to the integral,

where $k^=kk$ and j_ls(kr) are the first kind spherical Bessel functions. This integral exists and is finite, and it can be written as a power series converging to the radial factor; see Eqs. 10.1.1, 9.1.10, 6.1.12, and 7.4.4 of Ref. 163,

In Eq. (B11), we define a radial factor $Klspk=−iklsexp−k2/4ζp2ζp23+ls$⁠. Using this notation, the analytical expressions for the three components of gradient with respect to the coordinate k expressed in spherical coordinates |k|, θ, φ take the form