1. 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 (
), due to its low computational cost and reasonable accuracy [
2,
3,
4,
5]. Conventional time-dependent 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
, in terms of the GS density
, has not been discovered, it remains necessary to adopt density functional approximations (DFAs) for
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,
3,
4,
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
) evaluated at the instantaneous density
in the so-called adiabatic approximation [
7,
8,
9], the qualitative errors of
can also degrade the accuracy of conventional TD-DFT results [
14,
15,
16,
17,
18,
19].
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) [
20,
21,
22,
23]. The NCIE can be reduced by combining the parent DFA functionals with the dispersion energy correction (also known as dispersion-corrected functionals) [
24,
25] or with the second-order Møller-Plesset (MP2) correlation energy (often called double-hybrid functionals) [
23,
26]. The SCE can be reduced by incorporating a fully nonlocal correlation energy component, such as the RPA (random phase approximation) correlation energy [
27,
28], 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) [
29] (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 (
), 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 [
29,
30,
31,
32]. 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 [
29,
30,
31]. Consequently, TAO-DFT is very promising for studying the GS properties of large systems with strong static correlation effects [
33,
34,
35,
36,
37,
38,
39,
40,
41,
42,
43,
44,
45]. Other TAO-DFT extensions include the schemes that determine the system-independent [
46] and system-dependent [
47] values of
, TAO-DFT-based
ab initio molecular dynamics (for equilibrium thermodynamic and dynamical properties) [
48], and TAO-DFT-based polarizable continuum model (for solvation effects) [
49].
Within the framework of TAO-DFT, Yeh
et al. have recently proposed a frequency-domain formulation of linear-response time-dependent TAO-DFT [
50], 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 Eq. (6) and Eq. (B6) of Ref. [
50]) is defined with the TD pure state
of a noninteracting reference system (also see Appendix B1 of Ref. [
50]). However, the TD density
(see Eq. (5) of Ref. [
50]) in TDTAO-DFT is generally not associated with a TD noninteracting pure state
, 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
), at the initial time
, the initial density
is simply the TAO-DFT GS density
(see Eq. (1) of Ref. [
50]), which should be associated with a thermal ensemble [
29] (i.e., not associated with a pure state) of noninteracting electrons at a nonvanishing fictitious temperature (
). Therefore, the underlying assumption of TDTAO-DFT (i.e., that the TD density
is assumed to be associated with the TD pure state
of a noninteracting reference system) is generally incorrect, except only for the
case (wherein TDTAO-DFT reduces to conventional TD-DFT [
6] or more precisely, LR-TD-DFT [
10] by its inherent LR nature).
To resolve the aforementioned inconsistency of TDTAO-DFT (especially for
) [
50], in the work, we reformulate the TD extension of TAO-DFT by introducing a new reference system, consisting of an ensemble of noninteracting electrons moving in a TD local potential. This real-time (RT) extension of TAO-DFT is denoted as RT-TAO-DFT. Besides, since the assumption of a weak perturbation is not required in RT-TAO-DFT, we also employ RT-TAO-DFT to study strong-field electron dynamics in molecules as well as high-order harmonic generation (HHG) [
51,
52,
53,
54,
55,
56,
57,
58,
59,
60,
61,
62,
63,
64,
65,
66,
67,
68].
The rest of this paper is organized as follows. In Section II, we review TAO-DFT and discuss closely related electronic structure methods. The formulation of RT-TAO-DFT is presented in Section III. In Section IV, we describe the details of our RT-TAO-DFT calculations for the HHG spectra and related TD properties of molecular hydrogen
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 spin-symmetry breaking effects in the TD properties are also discussed. Our conclusions are provided in Section V.
4. 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 [
51,
52,
53,
54,
55,
56,
57,
58,
59,
60,
61,
62,
63,
64,
65,
66,
67,
68]. 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 [
53,
54].
HHG can be qualitatively described by the semiclassical three-step model [
51,
52]. 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 at the equilibrium and stretched geometries:
Here, the nuclei of 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
, which starts from the GS at time
, experiences an intense laser pulse for
. In order to mimic the commonly used Ti:sapphire laser [
95], 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 [
96]. The electric-field amplitude of the laser pulse
= 0.0534 a.u. (corresponding to the peak intensity
W/cm
), the laser frequency (also called the fundamental frequency)
= 1.5498 eV (corresponding to the wavelength
800 nm), and
= 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
. 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) [
68],
, is also employed for
. For an electronic system consisting of
atoms, the CAP function
is defined as the minimum of the values of the atom-centered spherical absorbing potentials:
Here,
is a spherical absorbing potential around the
I-th nucleus:
for
I = 1, ...,
, where
is the position of the
I-th nucleus. The cutoff radius
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
= 9.524 bohr (≈ 5.040 Å), the curvature
= 4.0 hartree/bohr
, and the maximum potential value
= 10 hartree.
To sum up, in the present HHG study, the RT-TAO effective one-electron Hamiltonian
(see Eq. (
47)) is given by
where
is the field-free RT-TAO effective one-electron Hamiltonian,
(see Eq. (
64)) is the strong-field interaction, and
(see Eq. (
65)) is the CAP.
To propagate the one-electron density matrix
in the time domain, we adopt a time step of
= 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
time steps. The modified mid-point unitary transformation (MMUT) algorithm [
93,
94] is employed for the numerical time propagation of
. As the use of the CAP breaks the conservation of the norm of the RT-TAO orbitals, the time propagation is no longer unitary [
68].
Following sufficient time propagation, the TD one-electron density matrix
is determined, and the TD density
is given by Eq. (
52), yielding various TD properties. While the RT-TAO orbital occupation numbers, which are the same as the GS TOONs
, 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
Besides, 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,
97]
where the HHG spectrum has been smoothed using the Hamming window function
to reduce the numerical noise. In the HHG spectrum, the harmonic order is defined as
, with
being the fundamental frequency (see Eq. (
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
(see Eq. (
46)) remains unknown, approximations to
are necessary for practical RT-TAO-DFT calculations. While the exact
formally depends on the density
at all previous times
, in this study, we make the adiabatic approximation:
where the TD xc
potential
is approximated by the GS xc
potential
(see Eq. (
6)) evaluated at the instantaneous density
. In the adiabatic approximation, since the exact xc
energy functional
also remains unknown, a DFA to
should be made as well. In this work, we adopt the LDA (i.e., the simplest DFA) xc
energy functional
, with
being the LDA xc energy functional [
11,
12] and
being the LDA
-dependent energy functional [
29]. For brevity, RT-TAO-DFT with the adiabatic LDA xc
potential is denoted as RT-TAO-ALDA. At time
, the initial state (i.e., the GS of the unperturbed
with a given bond length) is obtained with the underlying GS theory, TAO-LDA (i.e., TAO-DFT with the LDA xc
energy functional
) [
29]. Note that RT-TAO-ALDA (with
) 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
) corresponds to KS-LDA (i.e., KS-DFT with the LDA xc energy functional
) [
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,
29,
88]. At time
, the initial state (i.e., the GS of the unperturbed
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,
29,
88]. Besides, for
, since the TD external potential adopted is spin-independent (see Eq. (
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 [
98]. Therefore, the TD properties (which depend on the TD spin densities) of
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
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. For the special case of
, RT-TAO-ALDA reduces to TD-ALDA. All numerical results are obtained with a development version of Q-Chem 5.4 [
99].
Since the GS of the unperturbed
with an equilibrium bond length of 1.45 bohr exhibits mainly SR character, the spin-symmetry constraint can be satisfied by spin-unrestricted TAO-LDA (with
= 0, 7, 20, and 40 mhartree) [
29], producing the same up-spin and down-spin densities at time
. In addition, for
, 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 [
98]. 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
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
with a stretched bond length of 3.78 bohr exhibits a noticeable MR character [
29], and hence, the spin-symmetry constraint is violated with spin-unrestricted KS-LDA (i.e., TAO-LDA with
), producing symmetry-broken spin densities at time
. In this situation, even when a TD spin-independent external potential is applied to the initial state (i.e., a spin-symmetry-broken GS) for
, 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
with a stretched bond length of 3.78 bohr, obtained with spin-restricted and spin-unrestricted TD-ALDA (i.e., RT-TAO-ALDA with
) 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 of
, 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
by spin-unrestricted TAO-LDA (with
= 40 mhartree) [
29].