Time dependent adaptive configuration interaction applied to attosecond charge migration

Recent developments of high-intensity attosecond laser pulses¹ can provide valuable insight into the phenomena of charge migration, defined as the electronic motion, e.g., following ionization, that occurs before a nuclear response.^2–5 Charge migration has been used to explain site-selective reactivity in electronically excited peptides^6,7 and can potentially be used to direct chemical reactions into normally inaccessible pathways. While initial work on experimental measurement and control of charge migration is promising,⁸ theoretical techniques are required to understand specific electron dynamics pathways and to begin answering questions on the greater feasibility of charge migration controlled chemistry.

Cederbaum and Zoberly first introduced the idea that electron correlation in populated excited cationic states drives attosecond charge migration following ionization.^9–13 Following this work, numerous studies using the non-Dyson intermediate-state representation of the third or fourth order algebraic diagrammatic construction [ADC(3), ADC(4)],^5,14–23 time-dependent density functional theory (TD-DFT),^24–26 and time-dependent density matrix renormalization group (TD-DMRG)^27,28 have been employed to further understand pure electron dynamics following ionization. The major theoretical challenge of describing coherent electron dynamics following ionization is accurately characterizing the populated cationic states. ADC(3), ADC(4), and TD-DFT are affordable options to study charge migration, but they inherently depend on a single-reference description of correlation effects. Failure to characterize strong correlation in excited states can lead to a qualitatively incorrect description of the dynamics and is likely to occur when numerous near-degenerate cationic states with many coupled excitations become populated. TD-DMRG can treat an electronic state beyond the traditional excitation level hierarchy and is thus well-positioned for computing complicated dynamics, despite challenges associated with excited-state DMRG algorithms and propagating a matrix product state wave function.^27–30 

In this work, we generalize the adaptive configuration interaction (ACI)^31–33 method to a time-dependent version (TD-ACI) for simulating electron dynamics. ACI is based on the very old idea of selected CI^34–48 and is one of numerous new manifestations of these techniques.^49–62 TD-ACI can be viewed as an approximation to time-dependent complete active space (CAS) techniques, which have been successfully applied to photoionization processes,^63–68 where determinant screening in ACI enables the use of much larger CAS spaces. Due to its systematic improvability and demonstrated ability to treat many strongly correlated electrons in ground and excited states, TD-ACI will be essential for benchmarking time dependent methods and for providing insight into the role of electron correlation in attosecond electron dynamics.

Studying ultrafast dynamics of electrons requires the solution to the time-dependent Schrödinger equation (TDSE), in atomic units, $i∂∂tΨ(t)=ĤΨ(t)$⁠, which dictates how a wave function, Ψ(t), evolves in time according to the Hamiltonian (Ĥ). For a time-independent Ĥ, the evolution of a wave function can be written exactly as

where Δt is the time step over which the wave function evolves. To avoid the combinatorial complexity of the exact full configuration interaction (FCI) solution to the TDSE, many approximate time-dependent methods in electronic structure theory have been proposed to study pure electron dynamics, including real-time versions of Hartree–Fock (HF) theory,^69,70 configuration interaction,^67,71–76 multiconfigurational self consistent field wave functions,^77–79 density functional theory,^24,80–82 and coupled-cluster theories.^83–87 To study charge migration following ionization, Eq. (1) needs to be accurately approximated by (i) generating an accurate initial ionized state, (ii) building a tractably sized cationic Hamiltonian that well-describes all populated cationic states throughout propagation, and (iii) avoiding evaluation of the exponential using conventional numerical integrators.

We use ACI to build a compact determinantal representation of Ĥ and a well-defined initial ionized state. ACI uses an iterative screening algorithm to build a model space, M^N, composed of N-electron Slater determinants (Φ^N), such that the resultant wave function, $ΨN=∑ΦμN∈MNΦμNcμ$⁠, produces an energy error approximately equal to a user-defined parameter, σ (|E_ACI − E_FCI| ≈ σ).³¹ We have reported several algorithms to compute excited states with ACI, including a state-averaged ACI (SA-ACI) which optimizes a single model space by defining determinant importance as the contribution to the correlation energy averaged over several roots.³² This algorithm is very useful in formulating a time-dependent theory as it produces a compact Hamiltonian with a controllable average energy error over any number of roots.

To study charge migration following ionization, we first perform a SA-ACI computation on the lowest one or few states of the neutral molecule of interest to build the ground state wave function, $Ψ0N$⁠. We invoke the sudden ionization approximation,^5,9,27 where the initial wave function for the simulation [Ψ^N−1(t = 0)] is defined immediately after the ionization process by annihilating an electron in spin orbital ϕ_i and normalized appropriately as

The sudden ionization approximation provides a well-defined initial state suitable for benchmarking the propagation in TD-ACI field-free with a time-independent Hamiltonian. In practice, TD-ACI simulations can include the ground electronic state interacting with an ionizing pulse, but invoking the sudden ionization approximation allows us to effectively decouple errors associated with TD-ACI and those connected to our choice of initial conditions. While computationally convenient, the sudden ionization approximation will inevitably introduce discrepancies with the experiment, even when modeling quasi-field-free dynamics. Since ACI uses a linear expansion of Slater determinants, the initial state can be defined by either a localized hole or a superposition of holes without complications possible in other approaches.^26,88,89

We then define the basis for the cationic Hamiltonian, M^N−1, as the set of cationic determinants, {Φ^N−1}, generated by applying a single annihilation over all spin orbitals to all determinants in the original basis M^N, $MN−1={ΦIN−1:ΦIN−1=âiΦμN,∀i∈Aand∀ΦμN∈MN}$⁠, where A denotes the set of all active orbitals. Finally, we define the TD-ACI wave function with the time-dependence entirely encoded in the expansion coefficients [c_I(t)] as $ΨN−1(t)=∑ΦIN−1∈MN−1cI(t)ΦIN−1$⁠. Starting from the initial condition [Eq. (2)], the wave function Ψ^N−1(t) can be determined at any time t by integration of the TDSE. Exact propagation via Eq. (1) requires complete diagonalization of the cationic Hamiltonian and is unfeasible for realistic simulations. Instead, we use the fourth-order Runge–Kutta (RK4) algorithm with fixed time steps,^27,30,90,91 which for a time-independent Ĥ is equivalent to approximating the exact exponential propagator using a Taylor series truncated to fourth order. For the relatively short time steps (0.02–0.05 as) and total propagation times (2 fs) used in this work, we do not anticipate significant errors related to the lack of norm and energy conservation in the RK4 integrator. Moreover, we prefer RK4 to the norm and energy conserving Krylov subspace time evolution as both methods perform similarly at short time-scales, and RK4 avoids numerical instabilities associated with linear dependencies in the Krylov subspace. We monitor the migration of the ionized hole throughout the molecule using the hole occupation number, n_i for an orbital i, defined as the difference between the occupation of orbital i in the ground state and the ionized state at time t, $ni(t)=⟨Ψ0N|âi†âi|Ψ0N⟩−⟨ΨN−1(t)|âi†âi|ΨN−1(t)⟩$⁠.

We first study dynamics triggered by valence ionization in benzene to understand the effects of using a truncated cationic model space and an approximate time propagator. The dynamics of the ionized state is characterized by hole migration within the π/π^* manifold, with weak hole-occupation of a nearby σ-bonding orbital. The benzene geometry was optimized using DFT with a B3LYP functional and the cc-pVDZ basis set using the PSI4 program.⁹² ACI and CASCI computations use a CAS(8,8) containing the π/π^* valence space and the energetically nearest σ/σ^* bonding/antibonding pair, and they employ a restricted HF reference computed with the cc-pVDZ basis and density fitted integrals using the cc-pVDZ-JKFIT basis.⁹³ The computation was run using C₁ symmetry, but we refer to using corresponding D_2h labels for clarity, all plotted in Fig. S1. In the TD-ACI simulation, we use a time step of 0.05 as for a total time of 2 fs, and the initial state is prepared by annihilating the alpha 1b_1u spin orbital.

Upon ionization of the 1b_1u orbital, the hole migrates to a superposition of the degenerate 1b_2g and 1b_3g orbitals and back smoothly with a computed frequency of roughly 750 as using CASCI and the exponential propagator. The oscillation is faster by about 200 as compared to previously reported ADC(3) results because of the minimal CAS(8,8) employed in this work.²² This reduced model is nonetheless an effective test of our theory because the oscillating hole occupations require determinants different in character from those that define the ground cationic state. We test three different schemes to build M^N−1 from a SA-ACI computation of M^N optimized with respect to one, two, or four roots of the neutral species. These time-independent computations are run with values of σ chosen to produce similar dimensions of M^N−1 to facilitate comparisons.

With the various M^N−1 bases, we propagate the initial wave functions using the RK4 algorithm with a 0.05 as time step, and we plot the error in their overlaps with respect to CASCI in Fig. 1. Due to the short time step used, the observed errors are largely resultant from truncation alone, and we see it is not necessarily correlated with the number of determinants in each computation. For example, the propagations with M^N−1 generated from a single state [1(a)] show errors in the overlap up to 0.3 at 2 fs when the cationic space contains about 200 determinants. With a similar dimension, the two-state variant already shows an approximate 10-fold increase in accuracy. This increase in accuracy is due to the influence of two-hole/one-particle states and simple hole-excited states in the hole dynamics.²² To describe all of these states, M^N−1 needs determinants with single, double, or higher excitation character with respect to ionized ground and excited states. Our results indicate that such a determinantal makeup is most effectively achieved by ionizing determinants from an SA-ACI computation done with respect to several excited states.

The hole migration from the 1b_1u orbital to a superposition of 1b_2g and 1b_3g orbitals is also shown in Fig. 2. In the left plots, we show hole occupations computed with TD-ACI using the exact propagator with either one [2(a)] or four [2(b)] optimized roots in defining M^N−1. The plots on the right show the same data, but instead using the RK4 propagator. Our first observation is that the TD-ACI dynamics are indistinguishable when using either exact or RK4 propagatiors, showing that our expected propagation errors are negligible even with determinant selection so long as an appropriate time step is chosen. When M^N−1 is optimized with respect to one root, the hole occupation oscillates too quickly and does not transfer from the 1b_1u orbital (red) to the 1b_2g and 1b_3g orbitals (blue and green) with enough magnitude. Interestingly, the hole occupations of the 1b_2g and 1b_3g orbitals are correctly always identical. The four-state procedure to build a similarly sized M^N−1 gives occupations nearly indistinguishable from the exact result. These results indicate that TD-ACI is a viable technique in studying charge migration resultant from ionization, particularly if the basis for the cationic Hamiltonian is truncated using importance criteria that can consider multiple roots. This use of an SA-ACI computation to build the basis enables a faster convergence to the σ = 0 limit and does not bias the determinantal makeup of cationic Hamiltonian toward the ground cationic state.

Iodoacetylene is a valuable model of charge migration for theorists and experimentalists alike. Recently, Kraus et al. have used high-harmonic generation to achieve 100 as resolution and control over charge migration following ionization.^4,8 The dynamics following ionization was also investigated theoretically in a TD-DMRG study using a CAS(16,36) active space.²⁷ In iodoacetylene, two cationic states drive the dynamics, one characterized by a hole in the 5p-like orbitals of iodine perpendicular to the molecular axis, and the other having a hole in the two π-bonding orbitals in the acetylene triple bond. These states are near degenerate, and a multiconfigurational time-dependent approach may be required. The quasi-field-free dynamics of ionized iodoacetylene involve the hole migrating between the iodine and acetylene groups with an experimentally measured frequency of about 0.93 fs.^4,8 Our final test of TD-ACI is to reproduce this migration frequency using a truncated basis of cationic determinants.

We optimize the geometry with DFT using a PBE0⁹⁴ functional and def2-SVP basis set.⁹⁵ This basis was used in all iodoacetylene computations in addition to an effective core potential used to remove 28 core electrons. Within the remaining 59 orbitals, we use a CAS(16,22) for the neutral species, which includes all molecular orbitals generated from the 2p-like orbitals in the acetylene group and 4d and 5p-like molecular orbitals from the iodine that are subsequently split-localized. In the dynamics simulations, we use a time step of 0.02 as for a total simulation time of 2 fs. The initial state is prepared by annihilating a valence 5p-like orbital on the iodine atom, and all computations are run in C₁ symmetry with split-localized orbitals in the active space.

To analyze the charge migration dynamics, we compute a representation of the density of the ionized hole (ρ_h) by scaling the HF orbitals by the corresponding hole occupation,^9,27 ρ_h = ∑_i|ϕ_i|²n_i. Our first simulation uses a two-state SA-ACI computation with σ = 10 mE_h to generate a cationic basis containing 244 361 determinants, where the cationic Hilbert space contains roughly 5 × 10¹⁰ determinants. We plot the hole occupation numbers for the ionized 5p-like orbital on the iodine and the in-plane π-bonding orbital of the acetylene in Fig. 3 for the 2 fs simulation. The hole migrates from the iodine to the acetylene and back in about 0.91 fs, agreeing closely with the experimental frequency of 0.93 fs. The migration still leaves some degree of hole occupation on each moiety, which we show Fig. 3. The relatively good agreement between our computed frequency of hole migration and the experimental value suggests that the ACI procedure is building an adequate space of determinants to describe the relevant cationic states.

To test the sensitivity of the dynamics with respect to M^N−1, we show hole migration times with TD-ACI using cationic spaces built from one, two, and three root SA-ACI computations in Table I. The acetylene migration time is defined as the time required for the ionized hole, initially on the iodine, to maximally populate on the acetylene group. The iodine migration time is the time elapsed from initial ionization to repopulation of the hole density on iodine after it has migrated. For a nearly constant size of M^N−1, we see that increasing the number of states used to generate the cationic basis has a negligible effect on the hole migration times. Only when the total number of cationic determinants is increased, even when optimized for the neutral ground state, do we see the migration times approach the experimental and TD-DMRG values. This result depends solely on the determinantal makeup of the cationic states most important in forming the evolving wave packet. For iodoacetylene, these states are simple 1-hole cationic states with the hole located on the iodine or acetylene moieties, and no coupled electronic excitations are as significant as they were for our previous study on benzene. As a result, the dominant contribution for both hole states comes from different annihilations of the ground state and not from annihilations of excited states.

Our correlation treatment provides accurate dynamics despite the neglect of dynamical correlation. While fortuitous cancellation of error is possibly present and some dynamical correlation may be recovered within our active spaces, the success of the active-space treatment suggests that dynamical correlation is relatively unimportant in accurately defining the relevant cationic states for propagation. While dynamical correlation is likely necessary for accurately computing other time-dependent properties, TD-ACI shows great promise in computing time-dependent reference wave functions. A very appealing property of our TD-ACI approach is that we can effectively optimize M^N−1 regardless of the character of the relevant cationic states, whether simple or complex.

In this work, we have extended the ACI approach to simulate ultrafast electron dynamics. TD-ACI is able to effectively model the real-time propagation of the exact wave function using a truncated space of determinants selected by the time-independent ACI algorithm. We apply this methodology to charge migration dynamics that follow ultrafast ionization. We find that ACI is well-suited to find determinants relevant to the initial state, and determinants important at later points in time can be initially identified in excited state computations to reduce systematic increases in error as the wave function evolves. We also find that the single-state ACI computation can recover these determinants if a sufficiently large value of σ is used. Due to the short time scale of charge migration, we did not encounter significant issues in using approximations to the exact time propagator although these effects could be important in longer simulations. Propagation and truncation errors were studied using a charge migration model in benzene, where we found that cationic spaces optimized for multiple roots are needed to efficiently capture dynamics in which two-hole one-particle states are relevant. This type of state is not always relevant to the dynamics in propagating our initial state as we saw in our application to hole migration in iodoacetylene. Using various schemes to build our cationic basis, we were able to compute hole migration times between acetylene and iodine groups that match experimental results.

We anticipate numerous future directions of study. While the SA-ACI procedure was successful in building a determinantal space for the entire dynamics, we envision an even more efficient scheme where the selection and removal of determinants can be done during the simulation itself. For example, one can imagine an algorithm that estimates the importance of a determinant at a future time step using low-order approximate propagation. Additionally, integration of the current TD-ACI scheme with a time-dependent Hamiltonian would enable ab initio studies of molecules interacting with fluctuating electric or magnetic fields. Finally, we envision including dynamical correlation effects beyond our active space treatment by combining an effective Hamiltonian theory with our time-dependent reference wave functions.

See the supplementary material for the following: (i) discussion on the RK4 propagator, plot of active space used for benzene (PDF), and (ii) all geometries in the xyz format (ZIP).

This work was supported by the U.S. National Science Foundation under Grant No. CHE 1900532 and the Camille Dreyfus Teacher–Scholar Award (No. TC-18-045). J.B.S. would like to thank Dr. Tianyuan Zhang for helpful discussions on approximate propagation techniques.