Real-Time Extension of TAO-DFT

Thermally assisted occupation density functional theory (TAO-DFT) has been an efficient electronic structure method for studying the ground-state properties of large electronic systems with multi-reference character over the past few years. To explore the time-dependent (TD) properties of electronic systems (e.g., subject to an intense laser pulse), in this work, we propose a real-time (RT) extension of TAO-DFT, denoted as RT-TAO-DFT. Moreover, we employ RT-TAO-DFT to study the high-order harmonic generation (HHG) spectra and related TD properties of molecular hydrogen H2 at the equilibrium and stretched geometries, aligned along the polarization of an intense linearly polarized laser pulse. The TD properties obtained with RT-TAO-DFT are compared with those obtained with the widely used time-dependent Kohn–Sham (TDKS) method. In addition, issues related to the possible spin-symmetry breaking effects in the TD properties are discussed.


Introduction
Over the last thirty years, Kohn-Sham density functional theory (KS-DFT) [1] has been a popular electronic structure method for the ground-state (GS) properties of physical systems in the presence of static external potentials at zero electronic temperature (θ el = 0) due to its low computational cost and reasonable accuracy [2][3][4][5].Conventional timedependent density functional theory (TD-DFT) [6] (also called the time-dependent Kohn-Sham (TDKS) method, real-time TD-DFT (RT-TD-DFT), or real-time density functional theory (RT-DFT)), which is the time-dependent (TD) extension of KS-DFT, has been recently applied to explore the TD and excited-state properties of electronic systems under the influence of TD external potentials [7][8][9].Recently, a frequency-domain formulation of linear-response TD-DFT (LR-TD-DFT) [10] has also been adopted to obtain excitation energies (i.e., limited to the weak-field perturbative regime), owing to its computational efficiency and reasonable accuracy [7][8][9].Nevertheless, for the study of TD phenomena or excitation energies beyond the linear response, conventional TD-DFT [6], which involves propagating the TDKS equation in the time domain without any restriction to the TD external potentials, remains a promising method.
In KS-DFT [1], since the exact exchange-correlation (xc) energy functional E xc [ρ], in terms of the GS density ρ(r), has not been discovered, it remains necessary to adopt density functional approximations (DFAs) for E xc [ρ] to perform practical calculations [2][3][4][5].The xc energy functionals based on the frequently adopted DFAs, such as the LDA (local density approximation) [11,12] and GGAs (generalized gradient approximations) [13], are computationally efficient for the study of large systems.However, the DFA xc energy functionals have a few intrinsic shortcomings [2-5] and can yield the following qualitative errors: the self-interaction error (SIE), non-covalent interaction error (NCIE), and static correlation error (SCE).Since conventional TD-DFT [6], which usually takes the GS of a physical system as the initial state and often employs the GS xc potential (i.e., the functional derivative of E xc [ρ]) evaluated at the instantaneous density ρ(r, t) in the so-called adiabatic approximation [7][8][9], the qualitative errors of E xc [ρ] can also degrade the accuracy of conventional TD-DFT results [14 -17].
These qualitative errors can generally be reduced with the modification of the DFA functionals.For example, the SIE can be reduced by mixing the Hartree-Fock (HF) exchange energy into the parent DFA functionals (commonly called hybrid functionals) [18][19][20].The NCIE can be reduced by combining the parent DFA functionals with the dispersion energy correction (also known as dispersion-corrected functionals) [21,22] or with the secondorder Møller-Plesset (MP2) correlation energy (often called double-hybrid functionals) [23].The SCE can be reduced by incorporating a fully nonlocal correlation energy component, such as the RPA (random phase approximation) correlation energy [24,25], into the parent DFA functionals.Nonetheless, the DFA, dispersion-corrected, hybrid, and double-hybrid functionals fail to resolve the SCE problem, while the RPA and related functionals are very demanding in computational expense and hence are impractical for large systems.
To circumvent the SCE problem at low computational cost, thermally assisted occupation density functional theory (TAO-DFT) [26] (i.e., a density functional theory with fractional orbital occupations) has been recently developed.Note that TAO-DFT is an electronic structure method for the GS properties of physical systems at zero electronic temperature (θ el = 0), even though it adopts a reference system of noninteracting electrons at some fictitious temperature θ.The xc energy functionals developed in KS-DFT can also be used in TAO-DFT [26][27][28][29].Nonetheless, in strong contrast to KS-DFT, TAO-DFT, even with the commonly used DFA, dispersion-corrected, and hybrid functionals, can approximately describe strong static correlation effects, especially when an appropriate value of θ is chosen [26][27][28].Consequently, TAO-DFT is very promising for studying the GS properties of large systems with strong static correlation effects [30][31][32][33][34][35][36][37].Other TAO-DFT extensions include the schemes that determine the system-independent [38] and system-dependent [39] values of θ, TAO-DFT-based ab initio molecular dynamics (for equilibrium thermodynamic and dynamical properties) [40], and TAO-DFT-based polarizable continuum model (for solvation effects) [41].
Within the framework of TAO-DFT, Yeh et al. have recently proposed a frequencydomain formulation of linear-response time-dependent TAO-DFT [42], denoted as TDTAO-DFT (or, more precisely, LR-TDTAO-DFT by its inherent linear-response (LR) nature), allowing excitation energy calculations in the frequency domain (i.e., using Casida's formulation [10]).In TDTAO-DFT, the TD effective one-electron potential (see Equations ( 6) and (B6) of Ref. [42]) is defined with the TD pure state Ψ TAO (t) of a noninteracting reference system (also see Appendix B1 of Ref. [42]).However, the TD density ρ(r, t) (see Equation (5) of Ref. [42]) in TDTAO-DFT is generally not associated with a TD noninteracting pure state Ψ TAO (t) but associated with a TD noninteracting ensemble (which should be described by a TD density operator, as will be discussed later).For example, in TDTAO-DFT (with θ = 0), at the initial time t 0 , the initial density ρ(r, t 0 ) is simply the TAO-DFT GS density ρ(r) (see Equation (1) of Ref. [42]), which should be associated with a thermal ensemble [26] (i.e., not associated with a pure state) of noninteracting electrons at a nonvanishing fictitious temperature (θ = 0).Therefore, the underlying assumption of TDTAO-DFT (i.e., that the TD density ρ(r, t) is assumed to be associated with the TD pure state Ψ TAO (t) of a noninteracting reference system) is generally incorrect, except only for the θ = 0 case (wherein TDTAO-DFT reduces to conventional TD-DFT [6] or, more precisely, LR-TD-DFT [10] by its inherent LR nature).
The rest of this paper is organized as follows.In Section 2, we review TAO-DFT and discuss closely related electronic structure methods.The formulation of RT-TAO-DFT is presented in Section 3. In Section 4, we describe the details of our RT-TAO-DFT calculations for the HHG spectra and related TD properties of molecular hydrogen H 2 at the equilibrium and stretched geometries, aligned along the polarization of an intense linearly polarized laser pulse.The TD properties computed using RT-TAO-DFT are discussed and compared with the results of conventional TD-DFT [6].Moreover, issues related to the possible spinsymmetry breaking effects in the TD properties are also discussed.Our conclusions are provided in Section 5.

Overview of TAO-DFT
Consider a physical system of N interacting electrons moving in an external potential v ext (r) at zero electronic temperature (θ el = 0).In TAO-DFT [26], the GS density ρ(r) of the physical system is represented by the thermal equilibrium density of a reference system (called the thermally assisted occupation (TAO) reference system) of noninteracting electrons in the presence of a local potential v TAO (r) (called the TAO potential) at some fictitious temperature θ (i.e., the temperature of the TAO reference system).In other words, ρ(r) is represented by the TAO orbitals {φ j (r)} and TAO orbital occupation numbers (TOONs) { f j } (atomic units (a.u.) are adopted throughout this work): ( Here, f j is the occupation number of the j-th TAO orbital φ j (r), given by the Fermi-Dirac (FD) distribution function where 0 ≤ f j ≤ 1, j is the energy of the j-th TAO orbital φ j (r) and µ is the chemical potential chosen to conserve N (i.e., the number of electrons): On the basis of the Hohenberg-Kohn (HK) theorems [61] for the physical system at θ el = 0 and the Mermin theorems [62] for the TAO reference system at the fictitious temperature θ, a set of self-consistent equations (i.e., the TAO equations) that determine the TAO orbitals {φ j (r)}, the TAO orbital energies { j }, and hence the TOONs { f j } (see Equations ( 2) and ( 3)) and the GS density ρ(r) (see Equation ( 1)) are given by [26] ĥTAO (r)φ j (r) = j φ j (r).
Here, ĥTAO (r) is the TAO effective one-electron Hamiltonian: with the TAO potential (i.e., the TAO effective one-electron potential) where v ext (r) is the external potential of the physical system, v is the Hartree potential (i.e., the functional derivative of the Hartree energy functional δρ(r) is the xcθ potential, which is the functional derivative of the xcθ energy functional being the xc energy functional (as defined in KS-DFT [1]) and E θ [ρ] being the θ-dependent energy functional (e.g., see Equation (14) of Ref. [26]).
In TAO-DFT [26], to obtain the GS density ρ(r) (i.e., represented by Equation (1) with the TAO orbitals {φ j (r)} and TOONs { f j }) of the physical system, Equations ( 1) to ( 6) should be solved self-consistently.After the self-consistency is achieved, the GS energy E[ρ] of the physical system (at θ el = 0) is given by where the first term is the external potential energy, E H [ρ] and E xcθ [ρ] are the Hartree and xcθ energy functionals, respectively, and A θ s [{ f j , φ j }] is the noninteracting kinetic free energy at the fictitious temperature θ: i.e., the sum of the kinetic energy and entropy contribution of noninteracting electrons at the fictitious temperature θ, which can be exactly computed using the TAO orbitals {φ j (r)} and TOONs { f j }.Note that, for the special case of θ = 0, E θ=0 [ρ] = 0 and TAO-DFT (with E xcθ [ρ]) [26] reduce to KS-DFT (with E xc [ρ]) [1].

Density Representation in TAO-DFT
The GS density ρ(r) of a physical system is interacting v-representable (I-VR) as the exact ρ(r) belongs to a GS wavefunction of an interacting N-electron Hamiltonian for some external potential v ext (r), which can be exactly computed using the full configuration interaction (FCI) method at the complete basis set limit [63].Moreover, the exact ρ(r) can be represented by the natural orbitals {χ j (r)} and natural orbital occupation numbers (NOONs) {n j } [64]: where the NOONs {n j } satisfy the following conditions, Nevertheless, in KS-DFT [1], ρ(r) is assumed to be noninteracting pure-state vrepresentable (NI-PS-VR) as it belongs to a one-determinantal GS wavefunction of a noninteracting N-electron Hamiltonian (i.e., the Kohn-Sham (KS) Hamiltonian) for some local potential (i.e., the KS potential) [65][66][67].Accordingly, in KS-DFT, ρ(r) is represented by the occupied KS orbitals {φ KS i (r)}: As has been shown in a number of studies [66][67][68][69][70][71], there are some reasonable GS densities (e.g., the GS densities of some electronic systems with strong static correlation effects) that are not NI-PS-VR.Apparently, these GS densities cannot be obtained with KS-DFT even adopting the exact xc energy functional E xc [ρ].
On the other hand, in TAO-DFT [26], ρ(r) (given by Equation ( 1)) is assumed to be noninteracting thermal ensemble v-representable (NI-TE-VR) as it belongs to a thermal ensemble of a reference system of noninteracting electrons in the presence of a local potential (i.e., the TAO potential) at some fictitious temperature θ.Consequently, in TAO-DFT, ρ(r) is represented by the TAO orbitals {φ j (r)} and TOONs { f j }: where the TOONs { f j } (given by the FD distribution function) satisfy the following conditions, Owing to the similar expressions of the TAO-DFT GS density ρ TAO (r) (see Equation ( 14)) and the exact GS density ρ FCI (r) (see Equation (11)), the fictitious temperature θ in TAO-DFT can be so chosen that the distribution of TOONs is close to the distribution of the exact NOONs, which is closely related to the stability (i.e., the single-reference (SR)/multireference (MR) character) of the GS of an electronic system [26].Accordingly, the exact GS density is more likely to be NI-TE-VR with this choice of θ.In contrast to KS-DFT (i.e., TAO-DFT with θ = 0), TAO-DFT has an extra degree of freedom in choosing the θ value to improve the GS density representability.

Approximate Energy Functionals and Fictitious Temperatures in TAO-DFT
Since the exact xcθ energy functional E xcθ [ρ] (i.e., one of the key ingredients in TAO-DFT), in terms of the GS density ρ(r), has not been known, it remains necessary to employ DFAs for E xcθ [ρ] to perform practical calculations using TAO-DFT.Conventional DFAs, such as the LDA and GGAs, for E xcθ [ρ] (i.e., the DFA xcθ energy functional E DFA xcθ [ρ]) can be adopted [26,27].In addition to TAO-DFA (i.e., TAO-DFT with the DFA functional E DFA xcθ [ρ]), TAO-DFT with the exact exchange [28] and related hybrid functionals [28,29] may also be employed.
For the GS of an electronic system, the fictitious temperature θ of a given energy functional in TAO-DFT should be so selected that the distribution of TOONs simulates the distribution of the exact NOONs.In this situation, the static correlation associated with the electronic GS can be properly captured by the entropy contribution (see Equation ( 10)) in TAO-DFT [26].In other words, the optimal θ should be closely related to the SR/MR character of the electronic GS.For systems with electronic ground states possessing SR character (i.e., SR systems), all the NOONs should be close to either 0 (fully empty) or 1 (fully occupied), and, hence, the optimal θ values in TAO-DFT should be sufficiently small (but nonvanishing for real electronic systems [38]).However, for systems with electronic ground states possessing MR character (i.e., MR systems), the distributions of NOONs (and hence the optimal θ values) can be highly system-dependent.While it remains very challenging to devise a scheme that always yields the best θ of each system for a given energy functional in TAO-DFT, some progress has been achieved in recent years.
To improve the optimal system-independent θ scheme, a self-consistent scheme that determines the optimal θ values of electronic systems has been recently proposed [39].

Comparison of KS-DFT, TAO-DFT, and FT-DFT
Here, we compare three generally different electronic structure methods (see Table 1): KS-DFT [1], TAO-DFT [26], and FT-DFT (finite-temperature density functional theory, also called the Mermin-Kohn-Sham (MKS) method) [1,62], each of which employs a reference system of noninteracting electrons in the presence of a local potential at some fictitious temperature θ.

KS-DFT TAO-DFT FT-DFT
Both KS-DFT and TAO-DFT are electronic structure methods for the GS properties of physical systems at zero electronic temperature (θ el = 0).Note that θ ≡ θ el = 0 is assumed in KS-DFT, while the fictitious temperature (θ ≥ 0) can be different from the electronic temperature (θ el = 0) in TAO-DFT.Accordingly, the GS density ρ(r) of a physical system at θ el = 0 is assumed to be NI-PS-VR in KS-DFT and NI-TE-VR in TAO-DFT.Moreover, in KS-DFT, the HK universal functional (i.e., the sum of the interacting kinetic energy and the electron-electron repulsion energy at θ el = 0) [61], which is a functional of the GS density ρ(r), is given by where T s [{φ KS i }] (exactly computed using the occupied KS orbitals {φ KS i (r)}) is the noninteracting kinetic energy at zero fictitious temperature (θ = 0), and E xc [ρ] is the xc energy functional, which needs to be approximated for practical KS-DFT calculations.By contrast, in TAO-DFT, the HK universal functional [61], which is a functional of the GS density ρ(r), is expressed as where A θ s [{ f j , φ j }] (exactly computed using the TAO orbitals {φ j (r)} and TOONs { f j }) is the noninteracting kinetic free energy at the fictitious temperature θ, and E xcθ [ρ] is the xcθ energy functional, which needs to be approximated for practical TAO-DFT calculations.Note that TAO-DFT (with θ = 0) reduces to KS-DFT.
On the other hand, FT-DFT (i.e., the MKS method) is an electronic structure method for the thermal equilibrium properties of physical systems at finite electronic temperatures (θ el ≥ 0), and θ ≡ θ el is assumed in FT-DFT.Therefore, the thermal equilibrium density ρ θ el (r) of a physical system at θ el is assumed to be NI-TE-VR in FT-DFT.Moreover, in FT-DFT, the Mermin (M) universal functional (i.e., the sum of the interacting kinetic free energy and the electron-electron repulsion energy at θ el ) [62], which is a θ el -dependent functional of the thermal equilibrium density ρ θ el (r), is given by where Consequently, for the GS properties of physical systems at θ el = 0, FT-DFT reduces to KS-DFT, while TAO-DFT (with θ = 0) can be very different from KS-DFT (especially for MR systems) [26,38,41].

TAO-DFT (with
As mentioned previously, TAO-DFT is an electronic structure method for the GS properties of physical systems at zero electronic temperature (θ el = 0).In TAO-DFT, the xcθ energy functional . At a sufficiently small fictitious temperature (θ ≈ 0), the magnitude of E θ≈0 [ρ] should remain small compared to that of E xc [ρ], and hence the approximation ) is an approximate TAO-DFT method (good for θ ≈ 0), which may be adopted to describe the strong static correlation effects of some GS systems (wherever θ ≈ 0 can be an appropriate fictitious temperature) at θ el = 0.
On the other hand, FT-DFT is an electronic structure method for the thermal equilibrium properties of physical systems at finite electronic temperatures (θ el ≥ 0), wherein θ ≡ θ el is assumed.In FT-DFT, at zero electronic temperature (θ el = 0), F is an approximate FT-DFT method (good for θ ≡ θ el ≈ 0), which may be used to study the temperature effects of thermal equilibrium systems at θ el ≈ 0 [72,73].
According to their mathematical expressions, TAO-DFT with 62].However, owing to their distinctly different physical meanings, one can easily distinguish the two approximate methods simply based on the electronic properties computed.For example, for the GS properties of physical systems at θ el = 0, FT-DFT with F [26,38,41].Therefore, a number of recent studies on the GS properties of physical systems at θ el = 0 [74][75][76][77][78][79][80][81] have actually been performed using TAO-DFT with

KS-DFA with the rTAO Energy Correction
On the basis of the TAO-DFA (with some fictitious temperature θ) energy expression [26,27], Yeh, Yang, and Hsu have recently proposed a post-KS energy correction, called the rTAO energy correction [82], which is a θ-dependent energy correction evaluated with the KS-DFA (i.e., KS-DFT with the DFA xc energy functional) orbitals.Owing to the post-KS nature, for clarity, we denote this method as the KS-DFA+rTAO method.
For a small fictitious temperature (θ 40 mhartree), the KS-DFA+rTAO method has been shown to approximately reproduce the TAO-DFA results for some selected electronic properties [82].This may in part be due to the post-KS nature (i.e., good for small θ) and a very limited amount of test data.For a considerably large θ (e.g., the optimal θ of TAO-DFT with the exact exchange for the dissociation of molecular hydrogen H 2 [28,38]) or for other electronic properties (e.g., atomization energies), the results obtained with this method can be very different from those obtained with TAO-DFT.
More importantly, the GS density obtained with the KS-DFA+rTAO method (with any value of θ) is the same as the KS-DFA GS density (i.e., NI-PS-VR).In other words, some reasonable non-NI-PS-VR GS densities (e.g., the GS densities of some MR systems) cannot be obtained with KS-DFA and the KS-DFA+rTAO method (with any value of θ) [66][67][68][69][70][71].
In particular, whenever the spin-symmetry constraint [3,4,26,83] on the singlet GS density of an electronic system is violated with KS-DFA (which can commonly happen for MR systems), it must also be violated with the KS-DFA+rTAO method (with any value of θ).In such a situation, the spin-unrestricted KS-DFA/KS-DFA+rTAO results can be very different from the corresponding spin-restricted KS-DFA/KS-DFA+rTAO results, yielding unphysical spin-symmetry breaking effects in the spin-unrestricted KS-DFA/KS-DFA+rTAO calculations.By contrast, the spin-symmetry breaking issues can be greatly resolved by TAO-DFA (with an appropriate θ) [26][27][28]30,33,[38][39][40][41], highlighting the significance of the GS density representation in TAO-DFT.

RT-TAO Equation
Consider a physical system of N interacting electrons moving in a TD external potential v ext (r, t).The Hamiltonian operator of the physical system is given by containing the operators of the kinetic energy the electron-electron interaction and the TD external potential Let Ψ(t) be the TD state of the physical system.For most TD cases of physical interest, in this work, the time propagation is assumed to start from the GS Ψ GS (i.e., the lowest energy eigenstate of Ĥ(t 0 ), which is a stationary state) of the unperturbed physical system at time t = t 0 , and Ψ GS is assumed to be non-degenerate.Accordingly, the TD state Ψ(t) of the physical system is a solution of the TD Schrödinger equation (TDSE): with the initial state From the TD state Ψ(t) , the TD density ρ(r, t) of the physical system can be determined by where ρ(r) = ∑ N i=1 δ(r − r i ) is the number density operator.In particular, the initial density ρ(r, t 0 ) of the physical system is given by the GS density ρ GS (r) of the unperturbed physical system: According to the Runge-Gross (RG) theorems for the physical system (i.e., consisting of a TD pure state) [6], the TD state Ψ(t) is a functional of the TD density ρ(r, t) (i.e., formally depending on the density ρ(r, t ) at all previous times t ≤ t) and the initial state Ψ(t 0 ) , i.e., Ψ(t) = Ψ[ρ, Ψ(t 0 ) ](t) .In this work, the initial state Ψ(t 0 ) = Ψ GS is a functional of the initial density ρ(r, t 0 ) = ρ GS (r) based on the HK theorems [61].Since the dependence of initial state Ψ(t 0 ) is implicitly included in the TD density ρ(r, t), for brevity, Ψ[ρ, Ψ(t 0 ) ](t) is denoted as Ψ[ρ](t) hereafter.Now, we define the action functional of the physical system: where B[ρ] is a universal functional of the TD density ρ(r, t): Note that the action functional A[ρ] has a stationary point at the exact TD density ρ(r, t) of the physical system, given by the Euler equation: In order to develop an RT method compatible with TAO-DFT [26], we introduce the RT-TAO reference system, consisting of an ensemble of noninteracting electrons moving in a TD local potential v s (r, t).The RT-TAO reference system can interchange electrons with its environment, and, hence, the electron number N e in the RT-TAO reference system can be varied from 0 to ∞.The Hamiltonian operator of the RT-TAO reference system is given by Ĥs (t) = Ts + vs (t), (30) containing the operators of the kinetic energy Ts and the TD local potential vs (t).
At time t = t 0 , the time propagation starts from the grand canonical ensemble (i.e., a stationary ensemble) of the TAO reference system at some fictitious temperature θ, obtained with TAO-DFT [26].Therefore, the initial density operator Γs (t 0 ) of the RT-TAO reference system is given by the TAO density operator ΓTAO (i.e., the grand canonical density operator [84] of the TAO reference system at the fictitious temperature θ): with the equilibrium statistical weights satisfying the following conditions, Here, Φ 0 N e ,n denotes the n-th N e -electron eigenstate of Ĥs (t 0 ) = ĤTAO (i.e., the Hamilto- nian operator of the TAO reference system) and E N e ,n the corresponding energy eigenvalue.
For t > t 0 , the noninteracting ensemble of the RT-TAO reference system is generally non-stationary due to the presence of the TD local potential v s (r, t).According to the TDSE with the initial state Φ N e ,n (t 0 ) = Φ 0 N e ,n , we know how the TD state Φ N e ,n (t) evolves in time.Therefore, the TD noninteracting ensemble of the RT-TAO reference system can be properly defined by the following TD density operator Γs (t): where the statistical weights w N e ,n , which are assumed to be time-independent, are given by Equation (32), i.e., the initial statistical weights.Note that the TD density operator Γs (t) of the RT-TAO reference system is a solution of the Liouville-von Neumann equation: with the initial density operator Γs (t 0 ) = ΓTAO (given by Equation ( 31)).From the TD density operator Γs (t), the TD density ρ s (r, t) of the RT-TAO reference system can be determined by where ρs (r) is the number density operator [84].In particular, the initial density ρ s (r, t 0 ) of the RT-TAO reference system is given by the thermal equilibrium density ρ TAO (r) (see Equation ( 14)) of the TAO reference system [26], which is the same as the GS density ρ GS (r) of the unperturbed physical system (note that ρ GS (r) is assumed to be NI-TE-VR with this θ) and hence the initial density ρ(r, t 0 ) of the physical system (see Equation ( 26)): Here, we seek v s (r, t) that yields the same solution ρ s (r, t) = ρ(r, t) for t ≥ t 0 .According to the Li-Tong (LT) theorems for the RT-TAO reference system (i.e., consisting of a TD noninteracting ensemble) [85], the TD density operator Γs (t) is a functional of the TD density ρ s (r, t) (i.e., formally depending on the density ρ s (r, t ) at all previous times t ≤ t) and the initial density operator Γs (t 0 ), i.e., Γs (t) = Γs [ρ s , Γs (t 0 )](t).In this work, the initial density operator Γs (t 0 ) = ΓTAO (see Equation (31)) is a functional of the initial density ρ s (r, t 0 ) = ρ TAO (r) (see Equation ( 14)) based on the Mermin theorems [62].Since the dependence of initial density operator Γs (t 0 ) is implicitly included in the TD density ρ s (r, t), for brevity, Γs [ρ s , Γs (t 0 )](t) is denoted as Γs [ρ s ](t) hereafter.Now, we define the action functional of the RT-TAO reference system: where B s [ρ s ] is a universal functional of the TD density ρ s (r, t): Note that the action functional A s [ρ s ] has a stationary point at the exact TD density ρ s (r, t) of the RT-TAO reference system, given by the Euler equation: In RT-TAO-DFT, the universal functional B[ρ] (given by Equation ( 28)) is partitioned as where the universal functional B s [ρ] is given by Equation ( 40), A H [ρ] is the Hartree action functional: and A xcθ [ρ] is the xcθ action functional: which is a universal functional of the TD density ρ(r, t).Applying Equation (42) to Equation ( 29), we obtain Comparing Equation (41) with Equation (45) shows that the same solution ρ s (r, t) = ρ(r, t) can be obtained if we choose the TD effective one-electron potential v s (r, t) (up to a purely TD function C(t)) of the RT-TAO reference system as where v ext (r, t) is the TD external potential of the physical system, v H (r, t) = δA H [ρ]  δρ(r,t) = dr ρ(r ,t) |r−r | is the TD Hartree potential, and v xcθ (r, t) = δA xcθ [ρ]  δρ(r,t) is the TD xcθ potential.
Since the RT-TAO reference system consists of a TD noninteracting ensemble, the Hamiltonian Ĥs (t) (see Equation ( 30)) is separable, with the RT-TAO effective one-electron Hamiltonian ĥs (r, where the TD effective one-electron potential v s (r, t) (denoted as the RT-TAO potential) is given by Equation ( 46).The RT-TAO orbitals {φ j (r, t)} evolve in time according to the effective one-electron TDSE (denoted as the RT-TAO equation): with i.e., the initial j-th RT-TAO orbital φ j (r, t 0 ) is given by the j-th TAO orbital φ 0 j (r) (the j-th energy eigenfunction of ĥs (r, t 0 ) = ĥTAO (r) (see Equation ( 5)), associated with the GS of the unperturbed physical system) [26].The TD density ρ(r, t) of the physical system, which is the same as the TD density ρ s (r, t) (given by Equation ( 37)) of the RT-TAO reference system, can be computed using [84] where the occupation number f j of the j-th RT-TAO orbital φ j (r, t), which is time-independent, is given by Equation (2), i.e., its initial occupation number [26], also satisfying the conditions 0 ≤ f j ≤ 1 and ∑ j f j = N.
For the special case of θ = 0, RT-TAO-DFT (with the xcθ action functional A xcθ [ρ]) reduces to conventional TD-DFT (with the xc action functional A xc [ρ]) [6], providing that, at time t = t 0 , the initial state of the physical system is the non-degenerate GS of the unperturbed physical system.
Here, we discuss the representation of the TD density ρ(r, t) of a physical system in conventional TD-DFT [6] and RT-TAO-DFT.In conventional TD-DFT, ρ(r, t) is assumed to be TD noninteracting pure-state v-representable (TD-NI-PS-VR) as it belongs to a TD one-determinantal wavefunction of a noninteracting N-electron Hamiltonian for some TD local potential.By contrast, in RT-TAO-DFT, ρ(r, t) (given by Equation ( 50)) is assumed to be TD noninteracting ensemble v-representable (TD-NI-E-VR) as it belongs to a TD noninteracting ensemble (described by a TD density operator; e.g., see Equation ( 35)) in the presence of a TD local potential (i.e., the RT-TAO potential).
In RT-TAO-DFT, since we specify the initial state Ψ(t 0 ) = Ψ GS of the physical system and the initial density operator Γs (t 0 ) = ΓTAO of the RT-TAO reference system, the two conditions (i) the same initial density ρ(r, t 0 ) = ρ s (r, t 0 ) and (ii) the same initial time derivative of the density ∂ ∂t ρ(r, t) t=t 0 = ∂ ∂t ρ s (r, t) t=t 0 = 0 can be satisfied for the physical and RT-TAO reference systems, providing that the GS density ρ GS (r) of the unperturbed physical system is NI-TE-VR with a given value of θ (see Equation (38)).Note that conditions (i) and (ii), which ensure the existence of TD-NI-PS-VR densities [86], may also be the conditions for the existence of TD-NI-E-VR densities [87].

Matrix Representation
In RT-TAO-DFT, the j-th RT-TAO orbital φ j (r, t) can be expanded in the orthonormal one-electron basis {φ 0 p (r)}, spanned by the GS TAO orbitals (i.e., the TAO orbitals associated with the GS of the unperturbed physical system) [26]: where {C pj (t)} are the TD expansion coefficients.Accordingly, the TD density ρ(r, t) (see Equation ( 50)) can be expressed as where P(t) is the one-electron density matrix at time t, with matrix elements Moreover, at time t, F(t) is the RT-TAO matrix (commonly known as the Fock matrix), which is the matrix representation of the RT-TAO effective one-electron Hamiltonian ĥs (r, t) (see Equation ( 47)), with matrix elements In the orthonormal one-electron basis {φ 0 p (r)}, the RT-TAO equation (given by Equation ( 48)) can be reformulated in terms of the TD one-electron density matrix P(t) [8]: As the time propagation is assumed to start from the GS of the unperturbed physical system at time t = 0 (without loss of generality, the initial time t 0 ≡ 0 is assigned hereafter), the initial one-electron density matrix is given by and the initial RT-TAO matrix is given by where f p and p are the occupation number and energy, respectively, of the p-th GS TAO orbital φ 0 p (r) [26].The formal solution of RT-TAO equation (see Equation ( 55)) for the TD one-electron density matrix P(t) is given by [60,88,89] Here, U(t b , t a ) is a unitary time propagator from t a to t b : where T denotes time ordering.However, since F(t a ) does not necessarily commute with F(t b ) for t a = t b (see Equation ( 59)), it remains challenging to obtain P(t) directly from P(0) (see Equation ( 58)) for a long time interval [0, t] in RT-TAO-DFT.
In practical calculations, to reduce the error in the propagation for a long time interval [0, t], U(t, 0) is commonly split into a product of multiple time propagators, each corresponding to a small time step ∆t: where t n = n∆t denotes the value of t at the n-th time step, noting that t 0 = 0 and t m = m∆t = t.U(t n+1 , t n ) is the time propagator from t n to t n+1 = t n + ∆t, given by which takes P(t n ) to P(t n+1 ): We denote U n = U(t n+1 , t n ) and P n = P(t n ) for brevity and apply Equation ( 60) to Equation ( 58).Accordingly, the density matrix P m = P(t m ) = P(t) can be obtained from the initial density matrix P 0 = P(t 0 ) = P(0) via the following expression: For a sufficiently small time step ∆t, F(t n ) remains nearly commutative with F(t n + ∆t), and, hence, U n = U(t n+1 , t n ) (see Equation ( 61)) can be computed without considering the time ordering.Note, however, that the exact time-ordered propagator can only be obtained in the limit of an infinitesimal time step (i.e., ∆t → 0).Recently, various algorithms [60,88,89] have been developed for the numerical construction of the time propagation of TDKS equation [6], which may also be adopted for the time propagation of RT-TAO equation.
In short, it takes the following key steps to run an RT-TAO-DFT calculation for describing the time evolution of the electron density following a perturbation: • Construct the initial one-electron density matrix P(0) (see Equation ( 56)) and the initial RT-TAO matrix F(0) (see Equation ( 57)) for the GS of the unperturbed physical system at time t = 0 using TAO-DFT (i.e., the respective GS theory).

•
Apply the TD field to the physical system for t > 0, and propagate the one-electron density matrix P(t) and the RT-TAO matrix F(t) in the time domain, according to the RT-TAO equation (given by the matrix representation, e.g., see Equation ( 55)).

HHG Spectra from RT-TAO-DFT
HHG from an electronic system (e.g., an atom or molecule) is a nonlinear optical process driven by an intense laser pulse, wherein the laser frequency can be converted into its integer multiples [43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60].HHG has recently attracted much attention since it can be used to explore the structure and dynamics of electronic systems and chemical reactions on a femtosecond timescale.In addition, HHG can be employed to generate attosecond pulse trains as well as individual attosecond pulses [45,46].HHG can be qualitatively described by the semiclassical three-step model [43,44].First, an electron tunnels out from an electronic system in an intense laser field (i.e., tunnel ionization).Second, the electron is driven away from or back to the parent ion by the laser field.Finally, the electron recombines with the parent ion, emitting a high-energy photon.
In this work, we perform RT-TAO-DFT calculations to explicitly obtain the HHG spectra and related TD properties of molecular hydrogen H 2 at the equilibrium and stretched geometries: • H 2 with an equilibrium bond length of 1.45 bohr (≈0.767Å). • H 2 with a stretched bond length of 3.78 bohr (≈2.00 Å).
Here, the nuclei of H 2 are positioned along the z-axis (i.e., the laser polarization) with the center of mass being located at the origin.
To obtain the HHG spectrum, the electronic system H 2 , which starts from the GS at time t = 0, experiences an intense laser pulse for t > 0. In order to mimic the commonly used Ti:sapphire laser [90], the strong-field interaction is generated by a laser pulse with an oscillating electric field linearly polarized along the z-axis (see Figure 1): Here, the interaction with the electric field is treated in the dipole approximation and the length gauge [91].The electric-field amplitude of the laser pulse A 0 = 0.0534 a.u.
(corresponding to the peak intensity I 0 ≈ 1 × 10 14 W/cm 2 ), the laser frequency (also called the fundamental frequency) ω 0 = 1.5498 eV (corresponding to the wavelength λ 0 ≈ 800 nm), and σ p = 500 a.u.(≈12.1 fs) are adopted.In the HHG process, the electron released by tunnel ionization can travel far away from the center of the electronic system H 2 .To capture strong-field ionization process and to remove artificial reflections induced by the finite extent of Gaussian basis set (which will be adopted to describe the TAO/RT-TAO orbitals), a complex absorbing potential (CAP) [60], −iv CAP (r), is also employed for t > 0. For an electronic system consisting of N A atoms, the CAP function v CAP (r) is defined as the minimum of the values of the atom-centered spherical absorbing potentials: v CAP (r) = min{g 1 (r), ..., g N A (r)}. (65) Here, g I (r) is a spherical absorbing potential around the I-th nucleus: for I = 1, ..., N A , where R I is the position of the I-th nucleus.The cutoff radius r 0 should be small enough to interact with the electron density (because the space extended by the Gaussian basis set is finite), while it should be large enough to not overly perturb the original electronic system.Here, we adopt the cutoff radius r 0 = 9.524 bohr (≈5.040Å), the curvature η = 4.0 hartree/bohr 2 , and the maximum potential value V max = 10 hartree.
To propagate the one-electron density matrix P(t) in the time domain, we adopt a time step of ∆t = 0.02 a.u.(≈ 0.484 as) and a total propagation time of τ = 1000 a.u.(≈24.1 fs), which corresponds to a total of 5 × 10 4 time steps.The modified mid-point unitary transformation (MMUT) algorithm [88,89] is employed for the numerical time propagation of P(t).As the use of the CAP breaks the conservation of the norm of the RT-TAO orbitals, the time propagation is no longer unitary [60].
Following sufficient time propagation, the TD one-electron density matrix P(t) is determined, and the TD density ρ(r, t) is given by Equation ( 52), yielding various TD properties.While the RT-TAO orbital occupation numbers, which are the same as the GS TOONs { f j }, are time-independent, the norm of the RT-TAO orbitals can decrease with time t due to the CAP.To describe electron ionization, the number of bound electrons is computed using Moreover, the induced dipole moment along the laser polarization (i.e., the z-axis) is calculated by Accordingly, the HHG spectrum can be computed using [9,92] where the HHG spectrum has been smoothed using the Hamming window function w H (t) = 0.54 − 0.46 cos( 2πt τ ) to reduce the numerical noise.In the HHG spectrum, the harmonic order is defined as ω/ω 0 , with ω 0 being the fundamental frequency (see Equation ( 64)).
Here, we present the approximations made in the RT-TAO-DFT calculations, the computational details, and the numerical results.As the exact TD xcθ potential v xcθ (r, t) = δA xcθ [ρ]   δρ(r,t) (see Equation (46)) remains unknown, approximations to v xcθ (r, t) are necessary for practical RT-TAO-DFT calculations.While the exact v xcθ (r, t) formally depends on the density ρ(r, t ) at all previous times t ≤ t, in this study, we make the adiabatic approximation: where the TD xcθ potential v xcθ (r, t) is approximated by the GS xcθ potential δE xcθ [ρ] δρ(r) (see Equation ( 6)) evaluated at the instantaneous density ρ(r, t).In the adiabatic approximation, since the exact xcθ energy functional E xcθ [ρ] also remains unknown, a DFA to E xcθ [ρ] should be made as well.In this work, we adopt the LDA (i.e., the simplest DFA) xcθ being the LDA xc energy functional [11,12] and E LDA θ [ρ] being the LDA θ-dependent energy functional [26].For brevity, RT-TAO-DFT with the adiabatic LDA xcθ potential is denoted as RT-TAO-ALDA.At time t = 0, the initial state (i.e., the GS of the unperturbed H 2 with a given bond length) is obtained with the underlying GS theory, TAO-LDA (i.e., TAO-DFT with the LDA xcθ energy functional E LDA xcθ [ρ]) [26].Note that RT-TAO-ALDA (with θ = 0) corresponds to TD-ALDA (i.e., conventional TD-DFT with the adiabatic LDA xc potential) [7][8][9], and its underlying GS theory, TAO-LDA (with θ = 0), corresponds to KS-LDA (i.e., KS-DFT with the LDA xc energy functional E LDA xc [ρ]) [11,12].In this study, we also investigate the possible spin-symmetry breaking effects in the TD properties, in analogy to the GS counterparts [3,4,26,83].At time t = 0, the initial state (i.e., the GS of the unperturbed H 2 with a given bond length) is a singlet state, and, hence, the initial up-spin and down-spin densities obtained with an exact theory must be the same based on the spin-symmetry constraint [3,4,26,83].Moreover, for t > 0, since the TD external potential adopted is spin-independent (see Equation ( 67)), the up-spin and down-spin densities, which are equally propagated in the time domain, must be the same at any subsequent time t [93].Therefore, the TD properties (which depend on the TD spin densities) of H 2 obtained with the spin-unrestricted formalism must be identical to those obtained with the spin-restricted formalism.
To examine whether this spin-symmetry constraint can be satisfied by RT-TAO-ALDA, we perform spin-restricted and spin-unrestricted RT-TAO-ALDA (with θ = 0, 7, 20, and 40 mhartree) calculations for the TD properties, such as the number of bound electrons, induced dipole moment, and HHG spectrum, of H 2 at the equilibrium and stretched geometries (aligned along the polarization of an intense linearly polarized laser pulse) using the d-aug-cc-pVTZ basis set and a high-quality grid EML(99,590) containing 99 Euler-Maclaurin radial grid points and 590 Lebedev angular grid points.We note that the choice of basis set can significantly affect the HHG spectrum [94].For the special case of θ = 0, RT-TAO-ALDA reduces to TD-ALDA.All numerical results are obtained with a development version of Q-Chem 5.4 [95].
Since the GS of the unperturbed H 2 with an equilibrium bond length of 1.45 bohr exhibits mainly SR character, the spin-symmetry constraint can be satisfied by spinunrestricted TAO-LDA (with θ = 0, 7, 20, and 40 mhartree) [26], producing the same up-spin and down-spin densities at time t = 0.In addition, for t > 0, owing to the use of a TD spin-independent external potential, the up-spin and down-spin densities, which are equally propagated in the time domain, should be the same at any subsequent time t [93].Therefore, the TD properties, such as the number of bound electrons (see Figure 2), induced dipole moment (see Figure 3), and HHG spectrum (see Figure 4), of H 2 with an equilibrium bond length of 1.45 bohr, obtained with spin-restricted and spin-unrestricted RT-TAO-ALDA (with θ = 0, 7, 20, and 40 mhartree), are essentially the same (i.e., within the numerical precision considered).
On the other hand, the GS of the unperturbed H 2 with a stretched bond length of 3.78 bohr exhibits a noticeable MR character [26], and, hence, the spin-symmetry constraint is violated with spin-unrestricted KS-LDA (i.e., TAO-LDA with θ = 0), producing symmetry-broken spin densities at time t = 0.In this situation, even when a TD spinindependent external potential is applied to the initial state (i.e., a spin-symmetry-broken GS) for t > 0, the TD effective one-electron potentials can be spin-dependent, and, hence, the up-spin and down-spin densities, which are unequally propagated in the time domain, can be very different at any subsequent time t.As shown, the TD properties, such as the number of bound electrons (see Figure 5), induced dipole moment (see Figure 6), and HHG spectrum (see Figure 7), of H 2 with a stretched bond length of 3.78 bohr, obtained with spin-restricted and spin-unrestricted TD-ALDA (i.e., RT-TAO-ALDA with θ = 0), are distinctly different, yielding unphysical spin-symmetry breaking effects in all the TD properties examined.Such an unphysical spin-symmetry breaking feature of spin-unrestricted TD-ALDA is apparently undesirable for RT simulations.By contrast, the spin-symmetry breaking effects in the TD properties obtained with RT-TAO-ALDA are shown to be reducible with the increase in θ at essentially no additional computational cost.In particular, the TD properties obtained with spin-restricted and spin-unrestricted RT-TAO-ALDA (with θ = 40 mhartree) are essentially the same, yielding essentially no unphysical spin-symmetry breaking effects in all the TD properties examined.This desirable feature can be attributed to the satisfaction of spin-symmetry constraint on the singlet GS density of the stretched H 2 by spin-unrestricted TAO-LDA (with θ = 40 mhartree) [26].Here, the θ = 0 case corresponds to TD-ALDA.

Conclusions
In conclusion, we have developed RT-TAO-DFT (i.e., an RT extension of TAO-DFT [26]), allowing the study of TD properties of both SR and MR systems.By resorting to an ensemble formalism, RT-TAO-DFT has resolved the aforementioned inconsistency of TDTAO-DFT (especially for θ = 0) [42].Since the assumption of a weak perturbation is not required in RT-TAO-DFT, spin-restricted and spin-unrestricted RT-TAO-DFT (with various θ) calculations have been performed to explore the TD properties (e.g., the number of bound electrons, induced dipole moment, and HHG spectrum) of H 2 at the equilibrium and stretched geometries, aligned along the polarization of an intense linearly polarized laser pulse.The TD properties obtained with RT-TAO-DFT (with various θ) have been compared with those obtained with conventional TD-DFT [6], which corresponds to RT-TAO-DFT (with θ = 0).Moreover, issues related to the possible spin-symmetry breaking effects in the TD properties are also discussed.
exactly computed using the MKS orbitals {φ MKS k (r)} and MKS orbital occupation numbers { f MKS k }) are the noninteracting kinetic free energy at the fictitious temperature θ ≡ θ el , and F θ el xc [ρ θ el ] is the xc free energy functional, which needs to be approximated for practical FT-DFT calculations.Note that FT-DFT (with θ el = 0) reduces to KS-DFT.

Figure 1 .
Figure 1.Electric field of the pulse adopted.

Figure 6 .
Figure6.Induced dipole moment for H 2 with a stretched bond length of 3.78 bohr, with spin-restricted and spin-unrestricted RT-TAO-ALDA (with various θ).Here, the θ = 0 case corresponds to TD-ALDA.