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.

## I. INTRODUCTION

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 excitations^{20}), perturbations in organic,^{21} bio-molecules,^{22} chiral molecules,^{23,24} metallic^{25–27} systems, periodic^{28,29} systems, semiconductor materials,^{30} optical cavities,^{31} electronic stopping,^{32} etc. Runge and Gross^{33} 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) formalism^{35} 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 spectra^{43,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 systems^{43,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 systems^{44} and small molecules^{43} 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 [(CH_{3})_{2}N–(CH=CH)_{n}H, *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.

## II. THEORY AND METHODS

_{e}is the anti-symmetric many-body electronic wave function (parametrically dependent on the nuclear coordinates), and $x\u20d7\u2261{ri,\sigma 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.

*Z*

_{A}is the atomic number and

**R**

_{A}is the spatial coordinates of nucleus

*A*,

**r**and

**r**′ are the electronic spatial coordinates, and $V\u0302ext$ 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.

### A. Time-dependent configuration interaction singles (TDCIS)

*i*,

*j*,

*k*}) and unoccupied ({

*a*,

*b*,

*c*}) orbitals, generating a configuration interaction singles (CIS) wave function,

^{58,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}

**a**is the time-dependent TDCIS wave function represented as a state-vector, and

**H**

_{e}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.

### B. Real-time time-dependent density functional theory (RT-TDDFT) and time-dependent Hartree–Fock (RT-TDHF)

^{33}is described by the following equation:

^{72,73}

*n*is the one-electron total density of the restricted, closed-shell system under consideration,

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 *v*_{XC} 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, *n*_{t′<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 *n*_{t} 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}

^{80}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}

**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.

*ϕ*

_{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

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}

## III. COMPUTATIONAL DETAILS

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)_{n}H (*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.

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 S_{0}, S_{1}, and S_{2} 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 S_{0} → S_{1} and S_{1} → S_{2} transitions. The state-averaging was carried out by assigning equal weights to the S_{0}, S_{1}, and S_{2} states. Energies are discussed in the text and given in the supplementary material.

^{72,73}

*ω*is the frequency of the oscillatory electric field. The molecular axes align best with the

*x*axis, and the S

_{0}→ S

_{1}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\u0302,y\u0302,z\u0302)$, are obtained from the GAUSSIAN program. The time-dependent dipole moment along each axis is given by

*χ*is the basis set used for the linear expansion of the MOs. The contribution of $V\u0302ext$ to the Hamiltonian matrix is given by

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.

*N*

_{i}(

*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,

**C**

_{i}(

*t*) is the time-dependent MO coefficient vector for orbital

*i*, and $NiCIS(0)$ is the occupation number of the

*i*th orbital calculated as the sum of its occupations over all the determinants used for a given static CIS calculation.

^{97}

## IV. RESULTS AND DISCUSSION

### A. Stationary state characterization

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 S_{1} (total dipole moments are given in the supplementary material). Given the same ground state, the S_{0} dipole moments for the three molecules are the same across the two methods. The large difference in ground S_{0} and excited S_{1} state dipole moments (consistently $>1$ D across the set of molecules) indicates that S_{1} is a charge-transfer state. The magnitude of the dipole moment of S_{1} 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 S_{0} → S_{1} density-difference plots show an increase in electron density in the *π*-cloud and a decrease around the N-center in the molecules (Fig. S1).

ΔE (eV)
. | μ_{x} (D)
. | |||
---|---|---|---|---|

LR-TDHF . | S_{0} → S_{1}
. | S_{0}
. | S_{1}
. | S_{0} → S_{1}
. |

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)
. | ||

CIS . | S_{0} → S_{1}
. | S_{0}
. | S_{1}
. | S_{0} → S_{1}
. |

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-TDHF . | S_{0} → S_{1}
. | S_{0}
. | S_{1}
. | S_{0} → S_{1}
. |

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)
. | ||

CIS . | S_{0} → S_{1}
. | S_{0}
. | S_{1}
. | S_{0} → S_{1}
. |

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 S_{0} → S_{1} 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 S_{0} to S_{1}.

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 S_{0} and S_{2} states has been shown to yield a resonant frequency that is an average of the S_{0} → S_{1} and S_{2} → S_{1} transitions when using the real-time TDHF method in H_{2} and HeH^{+} systems. Therefore, the S_{2} state is implicitly accounted for in real-time TDHF, and the gap between the S_{1} and S_{2} 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 S_{0} → S_{1} and S_{1} → S_{2} transitions were calculated using the SA-CASSCF method, where the S_{2} state corresponds to a doubly excited state with the nearly double occupation of the LUMO. The resulting S_{0} → S_{1} excitation energies show a similar trend as with CIS and LR-TDHF: decreasing excitation energies with increasing system size ($\Delta ES0\u2192S1=9.24,7.23,5.92$ eV). The energies for the next transition are smaller: $\Delta ES1\u2192S2=3.97,1.18,1.29$ eV. Given the energetically lower-lying S_{1} → S_{2} transition compared to the S_{0} → S_{1} 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.

### B. Electron dynamics

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 S_{0} → S_{1} transitions. We first excite the systems using a continuous field with frequency resonant with the S_{0} → S_{1} 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 S_{0} and S_{1} 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 $\Omega calc.=|\mu x,S0\u2192S1|\xd7Emax,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 S_{0} → S_{1} 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 S_{0} → S_{1} transition is primarily HOMO → LUMO in character, we monitor the occupation of the LUMO as a metric for the population of the S_{1} state. See Fig. 2 for population oscillation with an applied field of *E*_{max} = 0.003 a.u. Plots with *E*_{max} = 0.001 and 0.005 a.u. are shown in the supplementary material, Figs. S2(a) and S2(b).

The TDCIS behavior is as expected, with population oscillating between the S_{0} and S_{1} 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 S_{1} 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, *E*_{max} = 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 *E*_{max} = 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 *E*_{max} = 0.001 a.u. and ∼35% and ∼25% for Systems 2 and 3 with *E*_{max} = 0.005 a.u. Overall, the trend suggests that the errors in the transfer of electron density between S_{0} and S_{1} 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.

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 S_{0} and S_{1}, 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.

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 S_{0} → S_{1} transition populates states much higher lying than the S_{1} state and that the resonant field does not lead to well-behaved population oscillation between S_{0} and S_{1}.

#### 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 S_{1} state. In order to drive the system away from the initial S_{0} 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” S_{0} → S_{1} transition (the HOMO → LUMO transition being the major contributor to the S_{1} 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 S_{1} state (see Fig. S4 for the corresponding TDCIS MO occupation and S_{1} state population plot).

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 S_{0} to S_{1}, 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 S_{1} dipoles are shown with dashed lines, and we see that the TDCIS time-dependent dipole moment approaches this value and remains centered around the S_{1} dipole moment. There is very little oscillation for Systems 1 and 2, and a small amount of oscillation for System 3.

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 $(\mu x,avgRT\u2212\mu x,S0HF)/(\mu x,S1LR\u2212\mu x,S0HF)$, where $\mu x,S0HF$ is the HF ground state dipole moment, $\mu x,S1LR$ is the LR-TDHF S_{1} state dipole moment, and $\mu 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 S_{1} 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 S_{1} 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 S_{1} state produces an electron density indicative of some charge transfer but not as much as would be expected from the LR-TDHF S_{1} 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 S_{0} state and half the doubly excited S_{2} state than the S_{1} 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 S_{0} → S_{1} 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 S_{0} → S_{1} 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).

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 S_{1} → S_{2} transition energy, where the S_{0} → S_{2} transition was dark and the S_{2} state was predominantly a doubly occupied LUMO. As the LUMO became populated, the total electron density became a mixture of the S_{0} and the S_{2} states, with the peak frequency corresponding to the relative mixture of the S_{0} → S_{1} transition and the S_{2} → S_{1} transition. Therefore, the S_{1} → S_{2} peak became coupled to the S_{0} → S_{1} peak as the LUMO occupation grew. With the relatively larger systems considered here, it is difficult to characterize a similar S_{2} state, but with an estimate of the energies given by the SA-CASSCF method, we know that the energy of the S_{1} → S_{2} transition is likely smaller than the energy of the S_{0} → S_{1} transition, causing a shift to lower energies. Additionally, just based on the relative energies of the SA-CASSCF S_{0} → S_{1} and S_{1} → S_{2} 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.

## V. CONCLUSIONS

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 S_{1} 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 S_{0} state to a population consistent with the S_{1} 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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

^{16}O +

^{16}O collisions

*ab initio*approaches for high-harmonic generation spectroscopy

*ab initio*QED study

*Physical Chemistry: A Molecular Approach*

_{2}by time-dependent configuration interaction using density functional theory and the Tamm-Dancoff approximation

*Introduction to Quantum Mechanics*

*Time-Dependent Density-Functional Theory*

*Fundamentals of Time-Dependent Density Functional Theory*

*Ab initio*calculation of vibrational absorption and circular dichroism spectra using density functional force fields