Real-time (RT) electron density propagation with time-dependent density functional theory (TDDFT) or Hartree–Fock (TDHF) is one of the most popular methods to model the charge transfer in molecules and materials. However, both RT-TDHF and RT-TDDFT within the adiabatic approximation are known to produce inaccurate evolution of the electron density away from the ground state in model systems, leading to large errors in charge transfer and erroneous shifting of peaks in absorption spectra. Given the poor performance of these methods with small model systems and the widespread use of the methods with larger molecular and material systems, here we bridge the gap in our understanding of these methods and examine the size-dependence of errors in RT density propagation. We analyze the performance of RT density propagation for systems of increasing size during the application of a continuous resonant field to induce Rabi-like oscillations, during charge-transfer dynamics, and for peak shifting in simulated absorption spectra. We find that the errors in the electron dynamics are indeed size dependent for these phenomena, with the largest system producing the results most aligned with those expected from linear response theory. The results suggest that although the RT-TDHF and RT-TDDFT methods may produce severe errors for model systems, the errors in charge transfer and resonantly driven electron dynamics may be much less significant for more realistic, large-scale molecules and materials.

Simulating electronic dynamics is essential for developing an improved understanding of electronic processes in molecules and materials, including charge transfer and the evolution of excited states. Real-time time-dependent density-functional theory (TDDFT) [and time-dependent Hartree–Fock theory (TDHF)] methods can be effectively utilized in such cases as a means of simulating electronic dynamics at relatively affordable computational costs. Real-time TDDFT has been used to study electron (and nuclear) dynamics in a myriad of cases: multinucleon transfer reactions through molecular and atomic collisions,1,2 molecules in oscillating electromagnetic fields of varying strengths,3–7 high-harmonic generation,8,9 resonant excitation dynamics (e.g., charge transfer,10–14 excitation-energy transfer,15 strong-field ionization,16 core excitations,17–19 plasmonic excitations20), perturbations in organic,21 bio-molecules,22 chiral molecules,23,24 metallic25–27 systems, periodic28,29 systems, semiconductor materials,30 optical cavities,31 electronic stopping,32 etc. Runge and Gross33 proved that there exists a one-to-one mapping between the time-dependent density of a system and the external potential, justifying the use of TDDFT to simulate time-dependent electronic phenomena. This was further shown by van Leeuwen to be a special case of a more fundamental theorem, which states that the time-dependent density of a many-particle system can be reproduced by a unique external potential in another many-particle system with a differing two-particle interaction.34 

Although formally exact, an extremely common approximation made when using TDDFT in practice is the adiabatic approximation. One of the primary conditions within the formalism of TDDFT is that the evolving density depends on the initial state of the system and the history of its time-evolution. Within the Kohn–Sham (KS) formalism35 of non-interacting particles, this dependence is incorporated into the exchange-correlation (XC) potential of the time-dependent Kohn–Sham Hamiltonian. Assuming that this potential only depends on the instantaneous density and not the density at previous points in time or the initial wave function or initial Kohn–Sham state is the widely employed adiabatic approximation. This dependence on the initial state and the electron density at previous points in time is termed the “memory” of the potential, and in practice, incorporation of memory is rarely satisfied.36–41 

The adiabatic approximation leads to qualitative errors in electron dynamics:38,41 incorrect charge transfer dynamics,42 peak shifting in time-resolved absorption spectra43,44 due to violation of the resonance condition,45 and incorrect electron dynamics when driven at resonant frequencies.46 The errors associated with the adiabatic approximation are often referred to as memory effects; they arise from the use of only the instantaneous density in the exchange-correlation potential. The TDHF method, which can be considered a special case within the generalized Kohn–Sham formalism of TDDFT,47 suffers from similar issues when used for real-time electron density propagation due to the dependence of the Fock operator on the electronic density, creating a similar set of nonlinear equations.48,49

It has been shown in one-dimensional model systems (such as the asymmetric Hubbard dimer and electron–hydrogen atom scattering)42,44,45,50 as well as three-dimensional systems43,51,52 that adiabatic approximations made to the time-dependent exchange-correlation potential produce qualitative errors when used to simulate coherent electron dynamics, especially in cases involving charge transfer. In violation of one of the exact conditions that the time-dependent counterpart of the KS formalism (TDKS) XC potential must follow, some model systems44 and small molecules43 show frequency shifts in peaks in their absorption spectra upon being resonantly driven from stationary states. The degree and direction of this shift are dependent on the amount of population driven from the ground state and the transition frequencies between the higher-lying states.43 Because larger systems will have a smaller relative change in the total electron density upon excitation of a single electron, the reference electron density input into the potential will be closer to that of the ground state compared to small model systems. It is unclear how such large systems would behave when driven resonantly starting in stationary states, one of the common applications for real-time TDDFT methods. Given the popularity of TDDFT methods for large-scale chemical applications, understanding and mitigating any associated systemic size-dependence (including errors due to memory effects) is an important undertaking.

In this study, we aim to analyze the system size dependence of errors due to the adiabatic approximation. We restrict our analysis to electron dynamics generated with real-time TDHF and time-dependent configuration interaction singles (TDCIS) methods. For both methods, the time-independent, ground-state wave functions are exactly the same—the Hartree–Fock wave function. However, the TDCIS method is a wave function based method that does not suffer from errors due to the dependence of the potential on the density that is found with TDHF and TDDFT within the adiabatic approximation. By comparing TDHF and TDCIS, any consideration of electron correlation effects is removed from our analysis. For three molecules [(CH3)2N–(CH=CH)nH, n = 1, 2, 3] with increasing π-conjugation, we show that size-dependent effects do indeed exist due to the dependence of the potential on the time-dependent density, specifically in the behavior with a resonant field and in the shifting of peaks in the absorption spectrum. We examine the size-dependence of this qualitatively unphysical behavior to determine if TDHF and TDDFT methods could be more confidently applied to larger chemical systems.

Non-relativistic electronic evolution is governed by the time-dependent Schrödinger equation (TDSE). For the kind of external perturbations considered in the present case, the TDSE, within the Born–Oppenheimer approximation, is given by (in atomic units, used throughout the rest of this text),
itΨe(x,t)=Ĥe(r,r,t)Ψe(x,t),
(1)
where Ĥe is the electronic Hamiltonian made up of the one- and two-body electronic operators, Ψe is the anti-symmetric many-body electronic wave function (parametrically dependent on the nuclear coordinates), and x{ri,σi} represents the collective spin and spatial coordinates of all the electrons in the system. In the present case, we employ stationary nuclei, which is common for simulations of electronic processes on the femtosecond timescale.
The electronic Hamiltonian in Eq. (1), referred to as the exact Hamiltonian henceforth, is given as
Ĥe(r,r,t)=r22AZA|rRA|+V̂ext(r,t)+rr1|rr|,
(2)
where ZA is the atomic number and RA is the spatial coordinates of nucleus A, r and r′ are the electronic spatial coordinates, and V̂ext is an external, time-dependent potential applied to the system that couples with the electronic motion.

The exact solution of TDSE is often impractical except for systems with a small number of electrons, and one can proceed toward an approximate solution by assuming that Ψe can be described by a single Slater determinant built from a set of N one-electron orbitals in an N-electron system (this is the case in the Kohn–Sham and the Hartree–Fock formalisms).53,54 This approximation may be inadequate to describe the effects of strong electron correlation and can be rectified by incorporating multiple determinants into the wave function (as in the configuration interaction formalism).55,56

The methods used in the present study, TDCIS and real-time TDHF (RT-TDHF), share the restricted Hartree–Fock (RHF) wave function as their ground state. The two methods have been shown to exhibit similar time-dependent behavior in some cases.57 However, there are some key differences between the two methods in terms of how they incorporate the time-dependence of electronic motion, as outlined in Secs. II A and II B.

One of the simplest ways to solve for an electronic excited state is to include single-electron transitions between the occupied ({i, j, k}) and unoccupied ({a, b, c}) orbitals, generating a configuration interaction singles (CIS) wave function,
ΨCISM(x)=cHFMΦHF(x)+i,aciaMΦia(x).
(3)
Here ΨCISM is the CIS ansatz58,59 for the wave function corresponding to the stationary state M of the system under consideration. Using the superposition principle, we can further construct the time-dependent electronic wave function as a linear combination of the CIS wave functions,60 
ΨTDCIS(x,t)=MaM(t)ΨCISM(x).
(4)
The equation of motion for the electronic state-vector of the system is given by
ita(t)=He(t)a(t),
(5)
where a is the time-dependent TDCIS wave function represented as a state-vector, and He is the Hamiltonian matrix, both on the basis of the stationary CIS states. Our solution to Eq. (5) is obtained by using the first order term of the Magnus expansion.61–63 The set of reference orbitals obtained upon solving the HF equations in the present case and the subsequent CI states built from the corresponding Slater determinants remain unchanged during the TDCI propagation of the system’s electronic wave function for a given molecular geometry. The explicit time-dependence is incorporated entirely through the time-dependent expansion coefficients in the many-body TDCI wave function.

Accurate electron dynamics can be obtained by solving the TDCI equation of motion for various CI ansatzes beyond CIS.64–69 Although the scaling of TDCI methods is factorial in the number of electrons when using full CI, truncated CI ansatzes can be used to save computational cost.70,71 Including excitations beyond singles in the CI expansion adds to electron correlation effects and will lead to improved energies compared to CIS, which tends to overestimate excitation energies. However, we expect that the results presented here would be consistent with any CI ansatz, and we choose CIS as the reference ground state Hartree–Fock wave function is identical to that used with TDHF, providing an identical starting electron density and related properties.

The time-dependent counterpart of the KS formalism (TDKS), following from Eq. (1) and the Runge–Gross theorem,33 is described by the following equation:72,73
itϕiKS(x,t)=ĥS[n,{ϕjKS}](r,t)ϕiKS(x,t),
(6)
where ĥS is the KS Hamiltonian for a single electron, {ϕiKS} are the KS orbitals, and n is the one-electron total density of the restricted, closed-shell system under consideration,
n(r,t)=2iocc|ϕiKS(r,t)|2.
(7)
The KS Hamiltonian is a functional of the time-dependent one-electron density and the KS orbitals and is given by
ĥS[n,{ϕjKS}](r,t)=r22AZA|rRA|+V̂ext(r,t)+drn(r,t)|rr|+vXC[nt,nt<t,Ψ0,Φ0](r,t).
(8)

The fourth term in Eq. (8), known as the Hartree potential, corresponds to the classical electron–electron interaction energy contribution from the Coulomb operator, whereas vXC is the unknown exchange-correlation (XC) potential defined as the functional derivative of the XC energy with respect to the one-electron density. The XC potential formally depends on the initial state of the interacting system, Ψ0, and the Kohn–Sham initial state, Φ0, as well as the intermediate TDKS densities, nt′<t.38 

When using the real-time TDKS (RT-TDKS) method, the initial state and intermediate TDKS densities are often not taken into consideration, an assumption leading to what is known as the “adiabatic approximation,” where only the instantaneous electron density nt is input into the XC potential just as is performed with the ground-state density functional theory. The adiabatic approximation appears to work well in some instances, such as ionization,37,39,74,75 where electron removal may keep the system closer to the neutral reference state compared to electron excitation. There are several proposed XC functionals and time-dependent formulations that incorporate memory effects, at least partially.76–78 However, in none of these cases is the resonance condition, an important consideration when modeling coherent electron dynamics, satisfied.79 

Time-dependent Hartree–Fock (TDHF) theory has developed independently80 of Kohn–Sham TDDFT but can be thought of as a special case of the generalized Kohn–Sham formalism.47,81 The XC potential in Eq. (8) is then the exchange operator from HF theory. This substitution results in a Liouville–von Neumann equation involving the time-dependent Fock operator,82 
itP(t)=[F(t),P(t)],
(9)
where F′ is the Fock matrix, P′ is the density matrix, and the primes indicate that the basis in which the matrices are represented is orthogonal.
If the time-dependent HF molecular orbitals (MOs), {ϕi}, are expanded as linear combinations of a set of basis functions, ϕk(r, t) = ∑μcμk(t)χμ(r) ({cνi} being known as MO expansion coefficients), then we can define the corresponding density-matrix and Fock matrix elements to be
Pμν(t)=2iocccμi(t)cνi*(t),
(10)
Fμν(t)=drχμ*(r)r22AZA|rRA|+V̂ext(r,t)χν(r)+drdrχμ*(r)λ,σχλ*(r)|rr|[Pλσ(t)χσ(r)χν(r)Pλν(t)χν(r)χσ(r)].
(11)

Equation (9) can then be solved numerically, here using the modified mid-point unitary transformation (MMUT) algorithm, which is restarted every 100 time-steps to remedy energy drifts.6,82–86

The RT-TDHF method, owing to the dependence of the Fock operator on the electronic density that parallels the dependence of the Kohn–Sham Hamiltonian on the electronic density, suffers from similar issues as the RT-TDDFT electron dynamics within the adiabatic approximation, namely incorrect Rabi oscillations, peak shifting, and a poor description of charge transfer.48 The TDHF theory does not formally account for memory effects; therefore, there is no path forward for improving these deficiencies in the TDHF theory. However, it is important to understand the limitations and applicability of TDHF theory in comparison to TDCIS, a method similar to TDHF in terms of spatial nonlocality but one that has no density-dependence in its potential.

Note that within the matrix formulation of TDDFT, solving for the eigenvalues and eigenvectors of the TDDFT equations within the Tamm–Dancoff approximation (TDA)87 leads to a set of matrix equations nearly identical to the CIS equations. Both the TDA and CIS tend to lead to an overestimation of excitation energies compared to their full matrix equivalents, TDDFT (without TDA) and TDHF.88–91 

The set of molecules chosen for the study are π-conjugated systems that increase systematically in size with 40, 54, and 68 electrons. The systems are characterized by a dimethyl amide donor moiety, (CH3)2N, attached to a series of unsaturated, π-conjugated hydrocarbon chains increasing in length, (CH=CH)nH (n = 1, 213), acting as the electron-accepting moieties (Fig. 1). The amide group acts as an electron donor, overall increasing the excited state charge transfer character compared to simple π-conjugated hydrocarbons and leading to a larger change in the dipole moment between the ground state and excited state.

FIG. 1.

Set of π-conjugated molecules with increasing size.

FIG. 1.

Set of π-conjugated molecules with increasing size.

Close modal

Molecular geometries were optimized at the B3LYP/6-31+G* level of theory.92–94 The characterization of the excited states was performed with linear response TDHF/STO-3G and CIS/STO-3G calculations using the GAUSSIAN quantum chemistry program (Development Version i14+).95 The STO-3G basis set is chosen for the sake of ease of analysis of the composition of the excited states and electron dynamics, with some results shown for the larger 6-31G basis set in the supplementary material. The main conclusions of this study are not affected by the choice of basis set. We additionally show in the supplementary material the results of using an unrestricted wave function for propagation, which we find are identical to the closed-shell results presented here.

To obtain an estimate of the energy of the doubly excited state, the S0, S1, and S2 states were computed using the state-averaged complete active space self-consistent field (SA-CASSCF) method,96 with active spaces of (4,3), (6,5), and (8,7) for Systems 1, 2, and 3, respectively. The active spaces were chosen to include the frontier n- and π-type orbitals involved in the S0 → S1 and S1 → S2 transitions. The state-averaging was carried out by assigning equal weights to the S0, S1, and S2 states. Energies are discussed in the text and given in the supplementary material.

The external perturbation chosen is the light–matter interaction described by a sinusoidal electric field, coupled with the electrons under the dipole approximation as72,73
V̂ext(r,t)=Er=Emaxsin(ωt)r,
(12)
where Emax=Emax,xî+Emax,yĵ+Emax,zk̂ is the amplitude vector and ω is the frequency of the oscillatory electric field. The molecular axes align best with the x axis, and the S0 → S1 transition possesses a large transition dipole moment along this axis; hence, the field is chosen to be linearly polarized along the x axis. The dipole moment matrices in the STO-3G basis, (μx, μy, μz) corresponding to the dipole moment operators (x̂,ŷ,ẑ), are obtained from the GAUSSIAN program. The time-dependent dipole moment along each axis is given by
μζ{x,y,z}(t)=Tr(Pχ(t)μζχ),
(13)
where χ is the basis set used for the linear expansion of the MOs. The contribution of V̂ext to the Hamiltonian matrix is given by
Vextχ(t)=Emax,xsin(ωt)μxχ.
(14)

We use a propagation step size of 0.002 fs (∼0.082 68 a.u.). Trajectories with an applied sinusoidal external field have been generated by perturbing the molecules from their ground states. Trajectories with a delta-kick perturbation are generated by applying a constant electric field to the ground state self-consistent field calculation before beginning electron propagation. The TDHF trajectories have been generated using the GAUSSIAN program, and the TDCIS trajectories have been generated using an in-house Python code that uses the GAUSSIAN generated CIS energies and transition dipoles as input parameters.

The time-dependent MO occupations, Ni(t), plotted for RT-TDHF and TDCIS are obtained by projecting the time-dependent MOs and wave functions onto the initial set of Hartree–Fock MOs and the Hartree–Fock wave function,
NiRT-TDHF(t)=Ci(t)P(t)Ci(t),
(15)
NiTDCIS(t)=|ai(t)|2×NiCIS(0).
(16)
Here Ci(t) is the time-dependent MO coefficient vector for orbital i, and NiCIS(0) is the occupation number of the ith orbital calculated as the sum of its occupations over all the determinants used for a given static CIS calculation.
Linear absorption spectra are generated using the time-dependent electron density matrices and Eq. (13). The total time-dependent electronic dipole moment, μ(t)=μx2(t)+μy2(t)+μz2(t), is used to calculate the dipole moment time-correlation function, which is Fourier-transformed to obtain the absorption cross-section at 0 K in the frequency domain as97 
α(ω)=2πω3dtμ(0)μ(t)eiωt.
(17)

To determine the stationary excited states of the molecules, CIS and linear response (LR)-TDHF calculations were performed. The resulting excitation energies and x-axis components of the excited state and transition dipole moments are listed in Table I for state S1 (total dipole moments are given in the supplementary material). Given the same ground state, the S0 dipole moments for the three molecules are the same across the two methods. The large difference in ground S0 and excited S1 state dipole moments (consistently >1 D across the set of molecules) indicates that S1 is a charge-transfer state. The magnitude of the dipole moment of S1 increases as the system size increases. This increase is, however, larger for excited state dipole moments calculated using CIS compared to TDHF. The difference in the discrepancy between excited state dipole moments decreases as the system size increases. The charge-transfer excitation energies calculated using the linear response methods decrease in value with increasing system size, and the corresponding S0 → S1 density-difference plots show an increase in electron density in the π-cloud and a decrease around the N-center in the molecules (Fig. S1).

TABLE I.

Transition energies and x-axis components of state and transition dipole moments calculated using linear response TDHF and CIS methods for Systems 1, 2, and 3.

ΔE (eV)μx (D)
LR-TDHFS0 → S1S0S1S0 → S1
System 1 8.95 −1.12 −2.84 −4.13 
System 2 7.05 −1.94 −4.29 −6.60 
System 3 5.99 −2.58 −5.47 −8.65 
ΔE (eV)μx (D)
CISS0 → S1S0S1S0 → S1
System 1 9.47 −1.12 −3.82 −4.52 
System 2 7.51 −1.94 −5.09 −7.04 
System 3 6.38 −2.58 −6.13 −9.33 
ΔE (eV)μx (D)
LR-TDHFS0 → S1S0S1S0 → S1
System 1 8.95 −1.12 −2.84 −4.13 
System 2 7.05 −1.94 −4.29 −6.60 
System 3 5.99 −2.58 −5.47 −8.65 
ΔE (eV)μx (D)
CISS0 → S1S0S1S0 → S1
System 1 9.47 −1.12 −3.82 −4.52 
System 2 7.51 −1.94 −5.09 −7.04 
System 3 6.38 −2.58 −6.13 −9.33 

Table S4 lists the composition of the S0 → S1 transition for the three molecules calculated in terms of one-electron transitions between orbitals. For System 1, the transition is primarily characterized by a HOMO → LUMO transition, with small contributions from other one-electron transitions. With increasing system size, the one-electron transitions involving the frontier orbitals beyond the HOMO and LUMO contribute to slightly higher extents, but the electronic character of the transition across the three molecules is primarily that of a one-electron HOMO → LUMO transition for both excited state methods. Corresponding density-difference plots show a similar nature of the re-distribution of the electron density upon excitation from the ground state S0 to S1.

Previous studies by one of the authors showed that higher lying excited states that are not accessible by linear response TDHF can contribute to the electron density evolved using the RT-TDHF method and that the resonance energy is affected by the superposition of the ground and excited states; the linear response of a superposition state composed of the S0 and S2 states has been shown to yield a resonant frequency that is an average of the S0 → S1 and S2 → S1 transitions when using the real-time TDHF method in H2 and HeH+ systems. Therefore, the S2 state is implicitly accounted for in real-time TDHF, and the gap between the S1 and S2 states can determine both the direction and the degree of the peak shifts.43,48 To determine the direction of the peak shifts, the excitation energies of the S0 → S1 and S1 → S2 transitions were calculated using the SA-CASSCF method, where the S2 state corresponds to a doubly excited state with the nearly double occupation of the LUMO. The resulting S0 → S1 excitation energies show a similar trend as with CIS and LR-TDHF: decreasing excitation energies with increasing system size (ΔES0S1=9.24,7.23,5.92 eV). The energies for the next transition are smaller: ΔES1S2=3.97,1.18,1.29 eV. Given the energetically lower-lying S1 → S2 transition compared to the S0 → S1 transition, the resonance peaks are expected to shift toward lower energies as the population of the LUMO increases and the contribution of the doubly excited configuration increases.

Moving to the time-domain, we can use both the RT-TDHF and TDCIS methods to propagate the electrons in the presence of an applied field. For RT-TDHF, the linear absorption spectra of the three systems perturbed with a δ-pulse show peaks at expected positions, corresponding to energies calculated with the linear response methods. To force the evolving TDHF electronic structure away from the regime where the adiabatic approximation holds, the molecules are significantly perturbed from their initial ground states using applied fields with field frequencies resonant with the linear response S0 → S1 transitions. We first excite the systems using a continuous field with frequency resonant with the S0 → S1 transition for each molecule to induce Rabi-like oscillations, and then later examine population and peak shifting trends with finite-time applied fields.

1. Continuous resonant field: Rabi oscillation

The external field with resonant frequency can be used to drive the system significantly away from its reference state (the Hartree–Fock state in the current study). With the application of a resonant field to a two-level system, the population of the S0 and S1 states should undergo Rabi oscillation, inverting at a frequency directly proportional to the field amplitude and the transition dipole moment between the two states, calculated as Ωcalc.=|μx,S0S1|×Emax,x. Previous studies on two-level systems have shown that the Rabi cycle-like oscillatory behavior of the electron density propagated using RT-TDHF can deviate significantly from the exact electron dynamics.49,98 The difference in electron dynamics is attributed to the detuning of the resonant frequency due to the dependence of the XC potential on the instantaneous electron density.98 This effect can be significant depending on the fraction of electron density rearranged with respect to the total electron density.38 Understanding the potential errors of a given methodology for modeling a system in the presence of a resonant field is important for simulating many varieties of pump-probe and nonlinear spectroscopy.

To examine the differences between the electron dynamics calculated using TDCIS and RT-TDHF methods for the three systems of interest in the presence of a resonant field, we apply a sinusoidal driving field for the entire duration of propagation with a field frequency equal to the S0 → S1 transition energy at maximum field strengths of 0.001, 0.003, and 0.005 a.u. The molecules are not formally two-level systems, so some deviation from ideal Rabi oscillation behavior is expected. As the populations of electronic states cannot be directly quantified with TDHF, we use the occupation of the virtual orbital space as a quantifier of the deviation of TDHF electron dynamics from ideal behavior. Because the S0 → S1 transition is primarily HOMO → LUMO in character, we monitor the occupation of the LUMO as a metric for the population of the S1 state. See Fig. 2 for population oscillation with an applied field of Emax = 0.003 a.u. Plots with Emax = 0.001 and 0.005 a.u. are shown in the supplementary material, Figs. S2(a) and S2(b).

FIG. 2.

Time-dependent LUMO occupations for Systems 1–3 with the field turned on, with applied field frequencies resonant with the S0 → S1 linear response transition energies, and a field amplitude of 0.003 a.u. The degree of electron density transferred to the LUMO changes with system size.

FIG. 2.

Time-dependent LUMO occupations for Systems 1–3 with the field turned on, with applied field frequencies resonant with the S0 → S1 linear response transition energies, and a field amplitude of 0.003 a.u. The degree of electron density transferred to the LUMO changes with system size.

Close modal

The TDCIS behavior is as expected, with population oscillating between the S0 and S1 states as seen by the near single occupation of the LUMO for each system; the maximum population of the LUMO is 0.95, 0.95, and 0.94 e for Systems 1, 2, and 3, respectively. The maximum value remains consistent throughout the duration of the applied field over many Rabi cycles. Additionally, the frequency of the population oscillation is well within the numerical error [O(10−5 a.u.)] of the expected Rabi frequency for all three systems.

In contrast, deviation from Rabi oscillation behavior is observed for TDHF propagation. The maximum LUMO population is 0.31, 0.58, and 0.76 for System 1 for fields with strengths of 0.001, 0.003, and 0.005 a.u. applied for 100, 60, and 20 fs, respectively. Systems 2 and 3 have a closer to the full population of the LUMO, with maximum occupation values of 0.84 and 0.99, respectively, for a field strength of 0.003 a.u. However, all systems also show an additional high frequency oscillation of the population when the LUMO is maximally occupied, indicative of a mixed state rather than being purely in the S1 state,46,48 and System 3 also no longer undergoes full population inversion after three Rabi cycles, which is presumably due to the occupation of higher lying states as we see the increased occupation of the virtual space (see the discussion below).

Another obvious discrepancy in the TDHF electron dynamics from the ideal Rabi behavior is in the frequency of oscillation, and here we find a strong size-dependent trend. We can quantify the trend by computing the Rabi frequencies using linear response transition dipole moments (Table I) and comparing this calculated Rabi oscillation frequency to those deduced from the oscillations of time-dependent LUMO occupations for applied electric fields with differing amplitudes, where we extract the Rabi oscillation time by taking the time for LUMO occupation to reach a minimum within the electron dynamics. A comparison of the LR calculated and RT observed Rabi oscillation frequencies is reported in Table II. For System 1, in the presence of the weakest field applied here, Emax = 0.001 a.u., the observed population oscillation has a frequency of 4.6 × 10−3, whereas that calculated is 1.6 × 10−3, yielding a deviation of ∼180%. The deviation decreases to ∼75% with a stronger field of Emax = 0.005 a.u. The agreement in the oscillation frequency is improved with increasing system size, going to deviations of ∼125% and ∼100% for Systems 2 and 3 with Emax = 0.001 a.u. and ∼35% and ∼25% for Systems 2 and 3 with Emax = 0.005 a.u. Overall, the trend suggests that the errors in the transfer of electron density between S0 and S1 in the presence of a resonant field tend to decrease with increasing system size, implying potentially a similar increase in accuracy in charge-transfer dynamics.

TABLE II.

Rabi frequencies (in ×10−3 a.u.) for field resonant with the S0 → S1 transition, calculated using linear response values for the transition dipole moment (Ωcalc.) and observed from the real-time electron dynamics (Ωobs.) for Systems 1–3, for field-amplitudes of Emax = 0.001, 0.003, and 0.005 a.u.

0.001 a.u.0.003 a.u.0.005 a.u.
TDHFΩcalc.Ωobs.Ωcalc.Ωobs.Ωcalc.Ωobs.
System 1 1.6 4.6 4.9 9.9 8.1 14.2 
System 2 2.6 5.9 7.8 12.4 13.0 17.6 
System 3 3.4 6.8 10.2 14.5 17.0 21.1 
0.001 a.u.0.003 a.u.0.005 a.u.
TDHFΩcalc.Ωobs.Ωcalc.Ωobs.Ωcalc.Ωobs.
System 1 1.6 4.6 4.9 9.9 8.1 14.2 
System 2 2.6 5.9 7.8 12.4 13.0 17.6 
System 3 3.4 6.8 10.2 14.5 17.0 21.1 

For a two-level system, increasing the intensity of the applied resonant field increases the frequency of population oscillation. With the three systems propagated by the TDCIS method, we indeed see the population driven at a faster rate between S0 and S1, with a little population in other states. To directly compare TDCIS and TDHF maximum populations, we next scan the time for which the external resonant sinusoidal field is turned on, as well as scan the field strength. The range of field strengths is between 0.005 and 0.075 a.u., and the number of field cycles is between 1 and 10. In Fig. 3, the map of maximum virtual MO occupations during the resonant field application is shown for a range of field strengths and a number of field cycles. The TDCIS method produces a nearly uniform maximum occupation of one electron as long as the field is on long enough and is intense enough.

FIG. 3.

Maximum virtual MO occupations (sum of LUMO, LUMO+1, and LUMO+2) within the duration of the applied field obtained from TDCIS (top) and RT-TDHF (bottom) propagation methods, plotted as a function of field-on time and field amplitude for Systems 1–3 perturbed with a sinusoidal electric field using CIS and linear response TDHF S0 → S1 transition-resonant frequencies.

FIG. 3.

Maximum virtual MO occupations (sum of LUMO, LUMO+1, and LUMO+2) within the duration of the applied field obtained from TDCIS (top) and RT-TDHF (bottom) propagation methods, plotted as a function of field-on time and field amplitude for Systems 1–3 perturbed with a sinusoidal electric field using CIS and linear response TDHF S0 → S1 transition-resonant frequencies.

Close modal

However, this behavior does not hold for the TDHF electron density propagation, where we instead see much different maximum virtual populations and much higher populations of the virtual orbital space. For all three systems, we see that a more intense field leads to a higher occupation of the virtual space, well beyond the single electron occupation of the LUMO. With the field parameters explored, System 1 approaches nearly double occupation of the LUMO, System 2 has three electrons in the virtual MOs, and System 3 has five to six electrons in the virtual space for the more intense fields. The greater population of System 3 could be expected, given that there are many more excited states that could be occupied by a larger system. Overall, this population map suggests that for TDHF electron propagation, applying a field resonant with the S0 → S1 transition populates states much higher lying than the S1 state and that the resonant field does not lead to well-behaved population oscillation between S0 and S1.

2. Single electron occupation of LUMO: Charge transfer

In this section and Sec. IV B 3, we analyze the behavior of TDHF electron dynamics when the electron occupation is close to that expected for the S1 state. In order to drive the system away from the initial S0 reference state, the molecules considered were again perturbed with fields using the linear response resonant frequencies. Field strengths and durations were chosen by picking trajectories from the set used to plot Fig. 3 that, when the field was removed, had average HOMO and LUMO occupations close to one for a “clean” S0 → S1 transition (the HOMO → LUMO transition being the major contributor to the S1 state, cf. Table S4). Then the field parameters corresponding to the most stationary evolution [constant MO occupation(s) with time] were chosen to induce the transition. In Fig. 4, we show the TDHF MO occupations with the chosen field parameters. For the three molecules, the final occupation is close to one electron in the HOMO and one electron in the LUMO, as expected for the S1 state (see Fig. S4 for the corresponding TDCIS MO occupation and S1 state population plot).

FIG. 4.

Time-dependent occupations obtained by projecting evolving molecular orbitals (MOs) onto the initial set of MOs, plotted for different perturbations with systems initialized in S0. The amplitudes of the perturbing fields are 0.01, 0.005, and 0.005 a.u., and the fields are turned on for 8, 8, and 5 cycles (gray area in the plot) for Systems 1, 2 and 3, respectively.

FIG. 4.

Time-dependent occupations obtained by projecting evolving molecular orbitals (MOs) onto the initial set of MOs, plotted for different perturbations with systems initialized in S0. The amplitudes of the perturbing fields are 0.01, 0.005, and 0.005 a.u., and the fields are turned on for 8, 8, and 5 cycles (gray area in the plot) for Systems 1, 2 and 3, respectively.

Close modal

We next analyze the degree of charge transfer during this population change by monitoring the time-dependent dipole moment for each molecule. Based on the stationary state characterization in Table I, we expect the dipole moment to increase as the population is driven from S0 to S1, with the change in dipole moment increasing with increasing system size. This is indeed what we see with both the TDCIS and TDHF methods, as shown in Fig. 5. For all systems, the initial dipole moment is identical for both methods, as both TDCIS and RT-TDHF have the same initial HF ground state. The stationary state S1 dipoles are shown with dashed lines, and we see that the TDCIS time-dependent dipole moment approaches this value and remains centered around the S1 dipole moment. There is very little oscillation for Systems 1 and 2, and a small amount of oscillation for System 3.

FIG. 5.

Moving-averaged time-dependent dipole moments (along the x axis) of the three molecules obtained from the resonant-frequency field perturbed RT-TDHF and TDCIS trajectories. The systems are initialized in the state S0. The dipole moments plotted are averaged over a window of 1 fs, or 500 time-steps, as a moving average.

FIG. 5.

Moving-averaged time-dependent dipole moments (along the x axis) of the three molecules obtained from the resonant-frequency field perturbed RT-TDHF and TDCIS trajectories. The systems are initialized in the state S0. The dipole moments plotted are averaged over a window of 1 fs, or 500 time-steps, as a moving average.

Close modal

Using the dipole moments along the x axis, we see that the RT-TDHF time-dependent dipole moment for System 1 approaches the LR-TDHF dipole moment of −2.84 D but does not complete the full charge transfer. We can quantify the extent of charge transferred as (μx,avgRTμx,S0HF)/(μx,S1LRμx,S0HF), where μx,S0HF is the HF ground state dipole moment, μx,S1LR is the LR-TDHF S1 state dipole moment, and μx,avgRT is the average RT-TDHF dipole moment (see Fig. S3). We find that System 1 completes ∼48% of the charge transfer when starting at −1.12 D and ending at an average value of ∼−1.97 D. The RT-TDHF dipole moment for System 2 also approaches the LR-TDHF S1 dipole moment but again falls short, obtaining an average value of ∼−3.44 D and completing only ∼64% of the expected charge transfer. Additionally, System 2 shows significant dipole oscillations that are larger in magnitude than the total change in dipole moment from the ground state. Similarly, System 3 does not quite get to the LR-TDHF S1 dipole moment, but the average dipole moment of ∼−4.57 D yields ∼69% of the charge transfer, so it completes the most of all three molecules. Similar to System 2, the dipole oscillations of System 3 are very large. These large dipole oscillations suggest a mixed state for RT-TDHF, whereas nearly identical MO occupations with TDCIS produce very little dipole oscillation, suggesting a near-stationary state.

Overall, driving the RT-TDHF MO occupation to what would be expected for the S1 state produces an electron density indicative of some charge transfer but not as much as would be expected from the LR-TDHF S1 dipole moments. In all cases, the RT-TDHF average dipole moment is smaller than the corresponding LR-TDHF value, showing that the charge transfer is incomplete. There is a size-dependent trend in the degree of charge transfer, with System 1 showing the largest disagreement in the dipole moment charge compared to linear response theory, with some improvement with increasing system size. However, the improvement in the average dipole moment change for Systems 2 and 3 comes with a large dipole oscillation, supporting that the systems are in a mixed state. A mixed state at the equal population of the HOMO and LUMO agrees with the finding in some previous work by one of the authors that suggests that single electron occupation of the HOMO and LUMO with single determinant RT-TDHF may be more representative of a mixed state composed of half the S0 state and half the doubly excited S2 state than the S1 state.

3. Single electron occupation of LUMO: Peak shifting in absorption spectra

Peak shifting in the absorption spectra when simulating time-dependent electron dynamics is a violation of the resonance condition. We physically expect that peaks in the spectrum will have changes in intensity as the population evolves, but that the energy of the peaks will remain constant rather than shifting. The unphysical peak shifting phenomenon has been observed when using real-time TDDFT methods with the adiabatic approximation as the density is perturbed far from the reference ground state.38,41,43,44 To determine if there are size-dependent trends in peak shifting, the chosen molecules are driven from their initial ground states using frequencies resonant with S0 → S1 transition energies calculated from LR-TDHF as described in Sec. IV B 2, resulting in a single occupation of the HOMO and LUMO. The system is then propagated with the field off for 40 fs, and the resulting field-off dipole moment is Fourier-transformed to produce the absorption spectrum for each system.

The resulting spectra are shown in Fig. 6, with a comparison given to the spectra computed from a trajectory with the electron density perturbed by a delta-kick electric field. The LR-TDHF S0 → S1 excitation energies are also shown in the spectra for reference, with the LR-TDHF energies being in excellent agreement with the delta-kick absorption spectra peaks. For all three molecules, the spectra generated from the single electron occupation show significant shifts from those obtained with the delta-kick pulse, in violation of the resonance condition. The peak shifts to lower energies by 0.85, 0.62, and 0.57 eV for Systems 1, 2, and 3, respectively, showing a decreasing shift with a larger system size. If the amount of the peak shift is determined as a percentage of the linear response excitation energy, then both Systems 1 and 3 show a shift of 9.5% and System 2 a shift of 8.8%, suggesting no size-dependent trend when examining percentage (See Table S5 in the supplementary material).

FIG. 6.

Peak shifts observed in the linear absorption spectra, calculated from time-dependent dipole moments obtained using RT-TDHF propagation, obtained from the resonant-frequency-perturbed systems in comparison to those obtained from the weak delta-pulse-perturbed systems. The frequency used for perturbing the molecules is resonant with the LR-TDHF S0 → S1 transition.

FIG. 6.

Peak shifts observed in the linear absorption spectra, calculated from time-dependent dipole moments obtained using RT-TDHF propagation, obtained from the resonant-frequency-perturbed systems in comparison to those obtained from the weak delta-pulse-perturbed systems. The frequency used for perturbing the molecules is resonant with the LR-TDHF S0 → S1 transition.

Close modal

In previous studies of TDDFT peak shifting in small model molecular systems, the direction and magnitude of the peak shift were found to be directly related to the S1 → S2 transition energy, where the S0 → S2 transition was dark and the S2 state was predominantly a doubly occupied LUMO. As the LUMO became populated, the total electron density became a mixture of the S0 and the S2 states, with the peak frequency corresponding to the relative mixture of the S0 → S1 transition and the S2 → S1 transition. Therefore, the S1 → S2 peak became coupled to the S0 → S1 peak as the LUMO occupation grew. With the relatively larger systems considered here, it is difficult to characterize a similar S2 state, but with an estimate of the energies given by the SA-CASSCF method, we know that the energy of the S1 → S2 transition is likely smaller than the energy of the S0 → S1 transition, causing a shift to lower energies. Additionally, just based on the relative energies of the SA-CASSCF S0 → S1 and S1 → S2 transitions, we would also expect a decrease in the degree of peak shifting with increasing system size. Overall, we find that the magnitude of the peak shift is significant, but it does indeed decrease as the system size increases, which may be due to the relative energies in the chosen molecules or to larger systems being less subject to errors due to the adiabatic approximation.

In this study, we simulate the electron dynamics for conjugated molecules of increasing size with 40, 54, and 68 electrons, which exhibit charge-transfer nature in the S1 state, using both TDCIS and RT-TDHF methods. Both methods have the same ground state, with TDCIS having correct physical behavior during field driven processes, whereas the RT-TDHF method has a density-dependent potential, as in RT-TDDFT, and exhibits similar errors as RT-TDDFT within the adiabatic approximation as the electron density is driven from the ground state. We here demonstrate for the first time that if the electron density is driven beyond the ground state using real-time electron density propagation, significant size-dependent errors are observed.

We resonantly drive electron dynamics to induce Rabi-like oscillatory behavior in the three molecules, showing that as the system size increases, we observe better agreement between the linear response Rabi frequency and the Rabi frequency calculated from real-time TDHF dynamics. For the case of resonant excitation from the S0 state to a population consistent with the S1 state, we find that the shifts in peaks of the absorption spectra and the degree of charge transfer obtained from RT-TDHF dynamics show better agreement with the expected LR result with increasing system size.

Overall, the RT-TDHF propagation of the electron density for the largest of the three molecules shows the smallest errors from the expected behavior. However, the errors are still significant for a system with 68 electrons, and it is unclear if the errors will continue to decrease with increasing system size. Although the larger system shows more significant state mixing, as seen in the large dipole oscillations and large population of the virtual orbital space, it may be that for charge transfer and resonantly driven electron dynamics, there are indeed smaller errors in real-time electron dynamics for large molecular and material systems. This finding provides key insight and guidance when performing simulations of the electron transfer pathways in molecules and materials.

Furthermore, details presented in the supplementary material include the Cartesian coordinates of the three systems studied here: ground and excited state dipole moments, MO transition coefficients, density differences, and time-dependent MO occupations.

This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award No. DE-SC0020203. We acknowledge the computational time on the MERCED cluster (funded by NSF ACI-1429783) and on the Pinnacles cluster (funded by NSF ACI-2019144). The authors appreciate helpful discussions with Dr. Hrant P. Hratchian.

The authors have no conflicts to disclose.

Karnamohit Ranka: Data curation (lead); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Christine M. Isborn: Conceptualization (lead); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
H.
Flocard
,
S. E.
Koonin
, and
M. S.
Weiss
, “
Three-dimensional time-dependent Hartree-Fock calculations: Application to 16O + 16O collisions
,”
Phys. Rev. C
17
(
5
),
1682
(
1978
).
2.
K.
Sekizawa
, “
TDHF theory and its extensions for the multinucleon transfer reaction: A mini review
,”
Front. Phys.
7
,
20
(
2019
).
3.
H.
Eshuis
,
G. G.
Balint-Kurti
, and
F. R.
Manby
, “
Dynamics of molecules in strong oscillating electric fields using time-dependent Hartree–Fock theory
,”
J. Chem. Phys.
128
(
11
),
114113
(
2008
).
4.
C. M.
Isborn
and
X.
Li
, “
Singlet–triplet transitions in real-time time-dependent Hartree–Fock/density functional theory
,”
J. Chem. Theory Comput.
5
(
9
),
2415
(
2009
).
5.
U.
De Giovannini
,
G.
Brunetto
,
A.
Castro
,
J.
Walkenhorst
, and
A.
Rubio
, “
Simulating pump–probe photoelectron and absorption spectroscopy on the attosecond timescale with time-dependent density functional theory
,”
ChemPhysChem
14
(
7
),
1363
(
2013
).
6.
M.
Wibowo
,
T. J. P.
Irons
, and
A. M.
Teale
, “
Modeling ultrafast electron dynamics in strong magnetic fields using real-time time-dependent electronic structure methods
,”
J. Chem. Theory Comput.
17
(
4
),
2137
(
2021
).
7.
A.
Ghosal
and
A. K.
Roy
, “
A real-time TDDFT scheme for strong-field interaction in Cartesian coordinate grid
,”
Chem. Phys. Lett.
796
,
139562
(
2022
).
8.
J.
Heslar
,
D. A.
Telnov
, and
S.-I.
Chu
, “
Subcycle dynamics of high-harmonic generation in valence-shell and virtual states of Ar atoms: A self-interaction-free time-dependent density-functional-theory approach
,”
Phys. Rev. A
91
(
2
),
023420
(
2015
).
9.
E.
Coccia
and
E.
Luppi
, “
Time-dependent ab initio approaches for high-harmonic generation spectroscopy
,”
J. Phys.: Condens. Matter
34
(
7
),
073001
(
2021
).
10.
K.
Lopata
and
N.
Govind
, “
Modeling fast electron dynamics with real-time time-dependent density functional theory: Application to small molecules and chromophores
,”
J. Chem. Theory Comput.
7
(
5
),
1344
(
2011
).
11.
C.
Andrea Rozzi
,
S.
Maria Falke
,
N.
Spallanzani
,
A.
Rubio
,
E.
Molinari
,
D.
Brida
,
M.
Maiuri
,
G.
Cerullo
,
H.
Schramm
,
J.
Christoffers
, and
C.
Lienau
, “
Quantum coherence controls the charge separation in a prototypical artificial light-harvesting system
,”
Nat. Commun.
4
(
1
),
1602
(
2013
).
12.
S. M.
Falke
,
C. A.
Rozzi
,
D.
Brida
,
M.
Maiuri
,
M.
Amato
,
E.
Sommer
,
A.
De Sio
,
A.
Rubio
,
G.
Cerullo
,
E.
Molinari
, and
C.
Lienau
, “
Coherent ultrafast charge transfer in an organic photovoltaic blend
,”
Science
344
(
6187
),
1001
(
2014
).
13.
S. A.
Fischer
,
C. J.
Cramer
, and
N.
Govind
, “
Excited state absorption from real-time time-dependent density functional theory
,”
J. Chem. Theory Comput.
11
(
9
),
4294
(
2015
).
14.
A.
Bruner
,
S.
Hernandez
,
F.
Mauger
,
P. M.
Abanador
,
D. J.
LaMaster
,
M. B.
Gaarde
,
K. J.
Schafer
, and
K.
Lopata
, “
Attosecond charge migration with TDDFT: Accurate dynamics from a well-defined initial state
,”
J. Phys. Chem. Lett.
8
(
17
),
3991
(
2017
).
15.
T.
Trepl
,
I.
Schelter
, and
S.
Kümmel
, “
Analyzing excitation-energy transfer based on the time-dependent density functional theory in real time
,”
J. Chem. Theory Comput.
18
(
11
),
6577
(
2022
).
16.
A.
Sissay
,
P.
Abanador
,
F.
Mauger
,
M.
Gaarde
,
K. J.
Schafer
, and
K.
Lopata
, “
Angle-dependent strong-field molecular ionization rates with tuned range-separated time-dependent density functional theory
,”
J. Chem. Phys.
145
(
9
),
094105
(
2016
).
17.
T.
Akama
and
H.
Nakai
, “
Short-time Fourier transform analysis of real-time time-dependent Hartree–Fock and time-dependent density functional theory calculations with Gaussian basis functions
,”
J. Chem. Phys.
132
(
5
),
054104
(
2010
).
18.
A.
Bruner
,
S. M.
Cavaletto
,
N.
Govind
, and
S.
Mukamel
, “
Resonant X-ray sum-frequency-generation spectroscopy of K-edges in acetyl fluoride
,”
J. Chem. Theory Comput.
15
(
12
),
6832
(
2019
).
19.
M.
Chen
and
K.
Lopata
, “
First-principles simulations of X-ray transient absorption for probing attosecond electron dynamics
,”
J. Chem. Theory Comput.
16
(
7
),
4470
(
2020
).
20.
T. P.
Rossi
,
M.
Kuisma
,
M. J.
Puska
,
R. M.
Nieminen
, and
P.
Erhart
, “
Kohn–Sham decomposition in real-time time-dependent density-functional theory: An efficient tool for analyzing plasmonic excitations
,”
J. Chem. Theory Comput.
13
(
10
),
4779
(
2017
).
21.
J.
Krumland
,
A. M.
Valencia
,
S.
Pittalis
,
C. A.
Rozzi
, and
C.
Cocchi
, “
Understanding real-time time-dependent density-functional theory simulations of ultrafast laser-induced dynamics in organic molecules
,”
J. Chem. Phys.
153
(
5
),
054106
(
2020
).
22.
A.
Dreuw
, “
Quantum chemical methods for the investigation of photoinitiated processes in biological systems: Theory and applications
,”
ChemPhysChem
7
(
11
),
2259
(
2006
).
23.
J.
Mattiat
and
S.
Luber
, “
Recent progress in the simulation of chiral systems with real time propagation methods
,”
Helv. Chim. Acta
104
(
12
),
e2100154
(
2021
).
24.
J.
Mattiat
and
S.
Luber
, “
Comparison of length, velocity, and symmetric gauges for the calculation of absorption and electric circular dichroism spectra with real-time time-dependent density functional theory
,”
J. Chem. Theory Comput.
18
(
9
),
5513
(
2022
).
25.
E. W.
Draeger
,
X.
Andrade
,
J. A.
Gunnels
,
A.
Bhatele
,
A.
Schleife
, and
A. A.
Correa
, “
Massively parallel first-principles simulation of electron dynamics in materials
,”
J. Parallel Distrib. Comput.
106
,
205
(
2017
).
26.
G. U.
Kuda-Singappulige
and
C. M.
Aikens
, “
Excited-state absorption in silver nanoclusters
,”
J. Phys. Chem. C
125
(
45
),
24996
(
2021
).
27.
L.
Bhan
,
C.
Covington
, and
K.
Varga
, “
Signatures of atomic structure in subfemtosecond laser-driven electron dynamics in nanogaps
,”
Phys. Rev. B
105
(
8
),
085416
(
2022
).
28.
C.
Shepard
,
R.
Zhou
,
D. C.
Yost
,
Y.
Yao
, and
Y.
Kanai
, “
Simulating electronic excitation and dynamics with real-time propagation approach to TDDFT within plane-wave pseudopotential formulation
,”
J. Chem. Phys.
155
(
10
),
100901
(
2021
).
29.
A.
Kononov
,
C.-W.
Lee
,
T. P.
dos Santos
,
B.
Robinson
,
Y.
Yao
,
Y.
Yao
,
X.
Andrade
,
A. D.
Baczewski
,
E.
Constantinescu
,
A. A.
Correa
,
Y.
Kanai
,
N.
Modine
, and
A.
Schleife
, “
Electron dynamics in extended systems within real-time time-dependent density-functional theory
,”
MRS Commun.
12
(
6
),
1002
(
2022
).
30.
S.
Yamada
,
M.
Noda
,
K.
Nobusada
, and
K.
Yabana
, “
Time-dependent density functional theory for interaction of ultrashort light pulse with thin materials
,”
Phys. Rev. B
98
(
24
),
245147
(
2018
).
31.
T. S.
Haugland
,
C.
Schäfer
,
E.
Ronca
,
A.
Rubio
, and
H.
Koch
, “
Intermolecular interactions in optical cavities: An ab initio QED study
,”
J. Chem. Phys.
154
(
9
),
094113
(
2021
).
32.
C.-K.
Li
,
J.-m.
Xue
, and
F.-S.
Zhang
, “
Channeling electronic stopping power of lithium ions in diamond: Contribution of projectile inner-shell electrons
,”
Phys. Rev. A
106
(
2
),
022807
(
2022
).
33.
E.
Runge
and
E. K. U.
Gross
, “
Density-functional theory for time-dependent systems
,”
Phys. Rev. Lett.
52
(
12
),
997
(
1984
).
34.
R.
van Leeuwen
, “
Mapping from densities to potentials in time-dependent density-functional theory
,”
Phys. Rev. Lett.
82
(
19
),
3863
(
1999
).
35.
W.
Kohn
and
L. J.
Sham
, “
Self-consistent equations including exchange and correlation effects
,”
Phys. Rev.
140
(
4A
),
A1133
(
1965
).
36.
N. T.
Maitra
,
K.
Burke
, and
C.
Woodward
, “
Memory in time-dependent density functional theory
,”
Phys. Rev. Lett.
89
(
2
),
023002
(
2002
).
37.
L.
Mancini
,
J. D.
Ramsden
,
M. J. P.
Hodgson
, and
R. W.
Godby
, “
Adiabatic and local approximations for the Kohn-Sham potential in time-dependent Hubbard chains
,”
Phys. Rev. B
89
(
19
),
195114
(
2014
).
38.
N. T.
Maitra
, “
Perspective: Fundamental aspects of time-dependent density functional theory
,”
J. Chem. Phys.
144
(
22
),
220901
(
2016
).
39.
M. J. P.
Hodgson
and
J.
Wetherell
, “
Accurate real-time evolution of electron densities and ground-state properties from generalized Kohn-Sham theory
,”
Phys. Rev. A
101
(
3
),
032502
(
2020
).
40.
M. T.
Entwistle
and
R. W.
Godby
, “
Exact nonadiabatic part of the Kohn-Sham potential and its fluidic approximation
,”
Phys. Rev. Mater.
4
(
3
),
035002
(
2020
).
41.
D.
Dar
,
L.
Lacombe
, and
N. T.
Maitra
, “
The exact exchange–correlation potential in time-dependent density functional theory: Choreographing electrons with steps and peaks
,”
Chem. Phys. Rev.
3
(
3
),
031307
(
2022
).
42.
N. T.
Maitra
, “
Charge transfer in time-dependent density functisonal theory
,”
J. Phys.: Condens. Matter
29
(
42
),
423001
(
2017
).
43.
M. R.
Provorse
,
B. F.
Habenicht
, and
C. M.
Isborn
, “
Peak-shifting in real-time time-dependent density functional theory
,”
J. Chem. Theory Comput.
11
(
10
),
4791
(
2015
).
44.
K.
Luo
,
J. I.
Fuks
, and
N. T.
Maitra
, “
Studies of spuriously shifting resonances in time-dependent density functional theory
,”
J. Chem. Phys.
145
(
4
),
044101
(
2016
).
45.
J. I.
Fuks
,
K.
Luo
,
E. D.
Sandoval
, and
N. T.
Maitra
, “
Time-resolved spectroscopy in time-dependent density functional theory: An exact condition
,”
Phys. Rev. Lett.
114
(
18
),
183002
(
2015
).
46.
B. F.
Habenicht
,
N. P.
Tani
,
M. R.
Provorse
, and
C. M.
Isborn
, “
Two-electron Rabi oscillations in real-time time-dependent density-functional theory
,”
J. Chem. Phys.
141
(
18
),
184112
(
2014
).
47.
A.
Seidl
,
A.
Görling
,
P.
Vogl
,
J. A.
Majewski
, and
M.
Levy
, “
Generalized Kohn-Sham schemes and the band-gap problem
,”
Phys. Rev. B
53
(
7
),
3764
(
1996
).
48.
C. M.
Isborn
and
X.
Li
, “
Modeling the doubly excited state with time-dependent Hartree–Fock and density functional theories
,”
J. Chem. Phys.
129
(
20
),
204107
(
2008
).
49.
M.
Ruggenthaler
and
D.
Bauer
, “
Rabi oscillations and few-level approximations in time-dependent density functional theory
,”
Phys. Rev. Lett.
102
(
23
),
233001
(
2009
).
50.
J. I.
Fuks
,
P.
Elliott
,
A.
Rubio
, and
N. T.
Maitra
, “
Dynamics of charge-transfer processes with time-dependent density functional theory
,”
J. Phys. Chem. Lett.
4
(
5
),
735
(
2013
).
51.
S.
Raghunathan
and
M.
Nest
, “
The lack of resonance problem in coherent control with real-time time-dependent density functional theory
,”
J. Chem. Theory Comput.
8
(
3
),
806
(
2012
).
52.
D.
Dar
,
L.
Lacombe
,
J.
Feist
, and
N. T.
Maitra
, “
Exact time-dependent density-functional theory for nonperturbative dynamics of the helium atom
,”
Phys. Rev. A
104
(
3
),
032821
(
2021
).
53.
R.
McWeeny
,
An Overview of Molecular Quantum Mechanics
(
Springer US
,
1992
), pp.
3
17
.
54.
D. A.
McQuarrie
and
J. D.
Simon
,
Physical Chemistry: A Molecular Approach
(
University Science Books
,
Sausalito, CA
,
1997
).
55.
Ideas of Quantum Chemistry
, edited by
L.
Piela
(
Elsevier
,
2007
).
56.
R. J.
Bartlett
and
M.
Musiał
, “
Coupled-cluster theory in quantum chemistry
,”
Rev. Mod. Phys.
79
(
1
),
291
(
2007
).
57.
P.
Hoerner
,
M. K.
Lee
, and
H. B.
Schlegel
, “
Angular dependence of strong field ionization of N2 by time-dependent configuration interaction using density functional theory and the Tamm-Dancoff approximation
,”
J. Chem. Phys.
151
(
5
),
054102
(
2019
).
58.
J. B.
Foresman
,
M.
Head-Gordon
,
J. A.
Pople
, and
M. J.
Frisch
, “
Toward a systematic molecular orbital theory for excited states
,”
J. Phys. Chem.
96
(
1
),
135
(
1992
).
59.
C.
David Sherrill
and
H. F.
Schaefer
III
, “
The configuration interaction method: Advances in highly correlated approaches
,”
Adv. Quantum Chem.
34
,
143
269
(
1999
).
60.
T.
Klamroth
, “
Laser-driven electron transfer through metal-insulator-metal contacts: Time-dependent configuration interaction singles calculations for a jellium model
,”
Phys. Rev. B
68
(
24
),
245421
(
2003
).
61.
F.
Casas
and
A.
Iserles
, “
Explicit Magnus expansions for nonlinear equations
,”
J. Phys. A: Math. Gen.
39
(
19
),
5445
(
2006
).
62.
S.
Blanes
,
F.
Casas
,
J. A.
Oteo
, and
J.
Ros
, “
The Magnus expansion and some of its applications
,”
Phys. Rep.
470
(
5–6
),
151
(
2009
).
63.
D.
Tannor
,
Introduction to Quantum Mechanics
,
1st ed.
(
University Science Books
,
2008
).
64.
P.
Krause
,
T.
Klamroth
, and
P.
Saalfrank
, “
Time-dependent configuration-interaction calculations of laser-pulse-driven many-electron dynamics: Controlled dipole switching in lithium cyanide
,”
J. Chem. Phys.
123
(
7
),
074105
(
2005
).
65.
P.
Krause
,
T.
Klamroth
, and
P.
Saalfrank
, “
Molecular response properties from explicitly time-dependent configuration interaction methods
,”
J. Chem. Phys.
127
(
3
),
034107
(
2007
).
66.
H. B.
Schlegel
,
S. M.
Smith
, and
X.
Li
, “
Electronic optical response of molecules in intense fields: Comparison of TD-HF, TD-CIS, and TD-CIS(D) approaches
,”
J. Chem. Phys.
126
(
24
),
244110
(
2007
).
67.
J. A.
Sonk
,
M.
Caricato
, and
H. B.
Schlegel
, “
TD-CI simulation of the electronic optical response of molecules in intense fields: Comparison of RPA, CIS, CIS(D), and EOM-CCSD
,”
J. Phys. Chem. A
115
(
18
),
4678
(
2011
).
68.
W.-T.
Peng
,
B. S.
Fales
, and
B. G.
Levine
, “
Simulating electron dynamics of complex molecules with time-dependent complete active space configuration interaction
,”
J. Chem. Theory Comput.
14
(
8
),
4129
(
2018
).
69.
R.
Ramakrishnan
, “
Charge-transfer selectivity and quantum interference in real-time electron dynamics: Gaining insights from time-dependent configuration interaction simulations
,”
J. Chem. Phys.
152
(
19
),
194111
(
2020
).
70.
I. S.
Ulusoy
,
Z.
Stewart
, and
A. K.
Wilson
, “
The role of the CI expansion length in time-dependent studies
,”
J. Chem. Phys.
148
(
1
),
014107
(
2018
).
71.
V.
Kochetov
and
S. I.
Bokarev
, “
RhoDyn: A ρ-TD-RASCI framework to study ultrafast electron dynamics in molecules
,”
J. Chem. Theory Comput.
18
(
1
),
46
(
2022
).
72.
C. A.
Ullrich
,
Time-Dependent Density-Functional Theory
(
Oxford University Press
,
2011
).
73.
Fundamentals of Time-Dependent Density Functional Theory
, edited by
M. A. L.
Marques
,
N. T.
Maitra
,
F.
Nogueira
,
E. K. U.
Gross
, and
A.
Rubio
(
Springer
,
2012
).
74.
R.
Baer
, “
Prevalence of the adiabatic exchange-correlation potential approximation in time-dependent density functional theory
,”
J. Mol. Struct.: THEOCHEM
914
(
1–3
),
19
(
2009
).
75.
I. S.
Wahyutama
,
D. D.
Jayasinghe
,
F.
Mauger
,
K.
Lopata
,
M. B.
Gaarde
, and
K. J.
Schafer
, “
All-electron, density-functional-based method for angle-resolved tunneling ionization in the adiabatic regime
,”
Phys. Rev. A
106
(
5
),
052211
(
2022
).
76.
C. A.
Ullrich
and
K.
Burke
, “
Excitation energies from time-dependent density-functional theory beyond the adiabatic approximation
,”
J. Chem. Phys.
121
(
1
),
28
(
2004
).
77.
H. O.
Wijewardane
and
C. A.
Ullrich
, “
Time-dependent Kohn-Sham theory with memory
,”
Phys. Rev. Lett.
95
(
8
),
086401
(
2005
).
78.
H. O.
Wijewardane
and
C. A.
Ullrich
, “
Real-time electron dynamics with exact-exchange time-dependent density-functional theory
,”
Phys. Rev. Lett.
100
(
5
),
056404
(
2008
).
79.
J. I.
Fuks
,
L.
Lacombe
,
S. E. B.
Nielsen
, and
N. T.
Maitra
, “
Exploring non-adiabatic approximations to the exchange-correlation functional of TDDFT
,”
Phys. Chem. Chem. Phys.
20
(
41
),
26145
(
2018
).
80.
A. D.
McLachlan
and
M. A.
Ball
, “
Time-dependent Hartree–Fock theory for molecules
,”
Rev. Mod. Phys.
36
(
3
),
844
(
1964
).
81.
R.
Baer
and
L.
Kronik
, “
Time-dependent generalized Kohn–Sham theory
,”
Eur. Phys. J. B
91
(
7
),
170
(
2018
).
82.
X.
Li
,
S. M.
Smith
,
A. N.
Markevitch
,
D. A.
Romanov
,
R. J.
Levis
, and
H. B.
Schlegel
, “
A time-dependent Hartree–Fock approach for studying the electronic optical response of molecules in intense fields
,”
Phys. Chem. Chem. Phys.
7
(
2
),
233
(
2005
).
83.
W.
Liang
,
C. T.
Chapman
, and
X.
Li
, “
Efficient first-principles electronic dynamics
,”
J. Chem. Phys.
134
(
18
),
184102
(
2011
).
84.
A.
Gómez Pueyo
,
M. A. L.
Marques
,
A.
Rubio
, and
A.
Castro
, “
Propagators for the time-dependent Kohn–Sham equations: Multistep, Runge–Kutta, exponential Runge–Kutta, and commutator free Magnus methods
,”
J. Chem. Theory Comput.
14
(
6
),
3040
(
2018
).
85.
Y.
Zhu
and
J. M.
Herbert
, “
Self-consistent predictor/corrector algorithms for stable and efficient integration of the time-dependent Kohn-Sham equation
,”
J. Chem. Phys.
148
(
4
),
044117
(
2018
).
86.
D. B.
Williams-Young
,
A.
Petrone
,
S.
Sun
,
T. F.
Stetina
,
P.
Lestrange
,
C. E.
Hoyer
,
D. R.
Nascimento
,
L.
Koulias
,
A.
Wildman
,
J.
Kasper
,
J. J.
Goings
,
F.
Ding
,
A.
Eugene DePrince
,
E. F.
Valeev
, and
X.
Li
, “
The Chronus Quantum software package
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci. Sci.
10
(
2
),
e1436
(
2019
).
87.
S.
Hirata
and
M.
Head-Gordon
, “
Time-dependent density functional theory within the Tamm–Dancoff approximation
,”
Chem. Phys. Lett.
314
(
3–4
),
291
(
1999
).
88.
F.
Cordova
,
L. J.
Doriol
,
A.
Ipatov
,
M. E.
Casida
,
C.
Filippi
, and
A.
Vela
, “
Troubleshooting time-dependent density-functional theory for photochemical applications: Oxirane
,”
J. Chem. Phys.
127
(
16
),
164111
(
2007
).
89.
C. M.
Isborn
,
N.
Luehr
,
I. S.
Ufimtsev
, and
T. J.
Martínez
, “
Excited-state electronic structure with configuration interaction singles and Tamm–Dancoff time-dependent density functional theory on graphical processing units
,”
J. Chem. Theory Comput.
7
(
6
),
1814
(
2011
).
90.
J. E.
Subotnik
, “
Communication: Configuration interaction singles has a large systematic bias against charge-transfer states
,”
J. Chem. Phys.
135
(
7
),
071104
(
2011
).
91.
S.
Grimme
, “
A simplified Tamm-Dancoff density functional approach for the electronic excitation spectra of very large molecules
,”
J. Chem. Phys.
138
(
24
),
244104
(
2013
).
92.
A. D.
Becke
, “
Density-functional thermochemistry. III. The role of exact exchange
,”
J. Chem. Phys.
98
(
7
),
5648
(
1993
).
93.
P. J.
Stephens
,
F. J.
Devlin
,
C. F.
Chabalowski
, and
M. J.
Frisch
, “
Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields
,”
J. Phys. Chem.
98
(
45
),
11623
(
1994
).
94.
B. P.
Pritchard
,
D.
Altarawy
,
B.
Didier
,
T. D.
Gibson
, and
T. L.
Windus
, “
New basis set exchange: An open, up-to-date resource for the molecular sciences community
,”
J. Chem. Inf. Model.
59
(
11
),
4814
(
2019
).
95.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
G. A.
Petersson
,
H.
Nakatsuji
,
X.
Li
,
M.
Caricato
,
A. V.
Marenich
,
J.
Bloino
,
B. G.
Janesko
,
R.
Gomperts
,
B.
Mennucci
,
H. P.
Hratchian
,
J. V.
Ortiz
,
A. F.
Izmaylov
,
J. L.
Sonnenberg
,
D.
Williams-Young
,
F.
Ding
,
F.
Lipparini
,
F.
Egidi
,
J.
Goings
,
B.
Peng
,
A.
Petrone
,
T.
Henderson
,
D.
Ranasinghe
,
V. G.
Zakrzewski
,
J.
Gao
,
N.
Rega
,
G.
Zheng
,
W.
Liang
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
K.
Throssell
,
J. A.
Montgomery
, Jr.
,
J. E.
Peralta
,
F.
Ogliaro
,
M. J.
Bearpark
,
J. J.
Heyd
,
E. N.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
T. A.
Keith
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A. P.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
J. M.
Millam
,
M.
Klene
,
C.
Adamo
,
R.
Cammi
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
O.
Farkas
,
J. B.
Foresman
, and
D. J.
Fox
, Gaussian Development Version Revision I.14+,
Gaussian, Inc.
,
Wallingford, CT
,
2018
.
96.
M. W.
Schmidt
and
M. S.
Gordon
, “
The construction and interpretation of MCSCF wavefunctions
,”
Annu. Rev. Phys. Chem.
49
(
1
),
233
(
1998
).
97.
J. L.
McHale
,
Molecular Spectroscopy
(
CRC Press
,
2017
).
98.
J. I.
Fuks
,
N.
Helbig
,
I. V.
Tokatly
, and
A.
Rubio
, “
Nonlinear phenomena in time-dependent density-functional theory: What Rabi oscillations can teach us
,”
Phys. Rev. B
84
(
7
),
075107
(
2011
).

Supplementary Material