Optical Spectrum Analysis of Real-Time TDDFT Using the Maximum Entropy Method

In the calculation of time-dependent density-functional theory in real time, we apply an external field to perturb the optimized electronic structure, and follow the time evolution of the dipole moment to calculate the oscillator strength distribution. We solve the time-dependent equation of motion, keeping track of the dipole moment as time-series data. We adopt Burg's maximum entropy method (MEM) to compute the spectrum of the oscillator strength, and apply this technique to several molecules. We find that MEM provides the oscillator strength distribution at high resolution even with a half of the evolution time of a simple FFT of the dynamic dipole moment. In this paper we show the effectiveness and efficiency of MEM in comparison with that of FFT. Not only the total number of time steps, but also the length of the autocorrelation, the lag, plays an important role in improving the resolution of the spectrum.


Introduction
Time-dependent density functional theory (TDDFT) is one of the most prominent and widely used methods for calculatiing excited states of medium-to-large molecules, and it is recognized as a powerful tool for studying electronic transitions of molecules [1,2]. TDDFT is, in principle, capable of analyzing electronic behaviors if the exact local exchange-correlation (xc) functional is known. In practical calculations, for valence-excited states with the excitation energy well below the ionization potential, accurate results have been obtained for an approximate xc functional [2].
We used TDDFT to study the optical properties of several materials within this framework where the electron-hole interaction is at the LDA level. Though the excitonic effects may be involved by solving the Bethe-Salpeter equation for the electron-hole Green's function [3], we neglect them for the sake of simplicity.
We focused on the maximum entropy method (MEM) to derive the electronic transition spectrum from the time-dependent dipole moment given by TDDFT, and evaluated the efficiency of the method. In our TDDFT, we employed a spatial grid with a fixed domain to express the physical system and solved the equation of motion by the finite difference approach, i.e., the real-time and real-space technique [4], without using explicit bases such as plane waves or Gaussian. The time evolution of the dipole moment of the system was tracked by solving the time-dependent equation of motion.
To determine the effectiveness of the new method, we used MEM, instead of using the FFT of the dynamic dipole moment to calculate the absorption spectrum. For smaller molecules, we found MEM provided spectra of fairly high resolution with a relatively small number of time-series data. MEM is attractive because a significant portion of the computational time is spent for calculations of time-series data of dipole moment. The theory with our numerical details is given in the next section, and the observed characteristics of MEM follow along with our conclusions.

Theory and numerical details 2.1. Time-dependent density-functional theory
In this section, we briefly describe our TDDFT calculation procedure [5]. The basis is the density functional theory (DFT) [6] with the local density approximation (LDA). The total energy of the ground state can be derived from the Kohn-Sham equation (KS) [7]. DFT is much less successful in describing optical responses and absorption spectra where electronic excited states are involved. However, this difficulty is, in principle, overcome by the extension of DFT to its time-dependent version, TDDFT, which was established by Runge and Gross [1].
In analogy to the time-independent case, the TDDFT equation of motion coupled with pseudopotentials is given by where V ps ion is the ionic pseudopotential, V H is the Hartree potential, and V XC is the xc potential. Since the exact time-dependent xc potential is not known, the originally non-local time-dependent xc potential is replaced with a time-independent local one. This is reasonable when the density varies slowly with time. This approximation allows the use of a standard local ground-state xc potential in the TDDFT framework. The Hartree and xc potentials can be determined from the electronic charge density, ρ(r, t) = j |ψ j (r, t)| 2 . The summation is over all the occupied states j. The Hartree potential is determined by ∇ 2 V H = −4πρ, and the usual local density approximation (LDA) is used as the xc potential V XC in our study. For the ionic potential, we employed the pseudopotential V ps ion in the separable form so that only the valence electrons were considered [8]. Prior to the calculation of optical responses, we first solved the usual, time-independent formulation of the pseudopotential-DFT method [9] to obtain the optimized electronic structure [10]. We then applied an external field V ext to the system as a perturbation and followed the linear responses of the system in real time. The xc functional we adopted is originally for the electronic ground states. This method of TDDFT has been effectively used, in cases where the potential was time-dependent, to study the behaviors of electrons in oscillating electric and magnetic fields and hence excited-state reactions [11]. In our calculations, the real-time and real space technique was adopted in solving Eq. (1) by means of the finite difference approach [4]. The uniform grid was used in our study for simplicity.
The time-dependent wave function is given by where H is the Hamiltonian of the system and k is a small wave number of the external perturbation in the z direction. In the linear response, the time-dependent polarizability is proportional to the dipole matrix element μ ξ (t) = ψ(t)|ξ|ψ(t) , where ξ represents x, y, and z. The polarizability α ξ (ω) is numerically calculated by the fast Fourier transform (FFT) of μ ξ (t) and averaged as α = (α x + α y + α z )/3. The oscillator strength distribution S(ω) is related to the imaginary part of the polarizability, In our FFT calculation for C 60 as shown in Fig. 1 was a certain amount of space between the edge and any atom and the absorption boundary condition was applicable. For the time evolution, the time step was Δt = 0.001 eV −1 and the total number of the steps N was varied. The calculated oscillator strength distribution of C 60 for N = 5000 shows a clear separation between the absorption peaks in comparison with that for N = 2000. The width of each peak was inversely proportional to the total evolution time, T = N Δt. In the spectrum for N = 5000, the low energy absorption peaks at 3.5, 4.4, and 5.6 eV agreed well with the experimentally observed peaks at 3.8, 4.8, and 5.8 eV [10]. The real-space algorithm offers the computational complexity of O(N 2 ) and is attractive for the calculation of large systems. Also of potential interest is the intuitive feature of the method: the real-space approach allows us to grasp a clear physical image of the wave functions, which is not the case when the wave functions are expanded to the explicit basis sets. With the realtime approach, we can follow the evolution of the wave functions much more directly. These advantages helped our coding procedures.

Maximum entropy method
As shown in the previous section, the resolution of the FFT absorption spectrum S(ω) is ∼ T −1 , corresponding to the finite spectral bandwidth of the limited time series of μ(t), where T = N Δt is the total evolution time. MEM, in contrast, estimates the full spectrum from the limited data based on information theory. The entropy in information theory has been recognized as a measure of uncertainty [12,13]. Any inferences made from incomplete information should use with the probability distribution which maximizes the entropy under the constraints of available information [14]. The present study is based on Burg's method [15], where the dipole moment μ is assumed to be a random variable and the samples are given by TDDFT, μ m = μ(mΔt).
The calculated autocorrelation C m at the time lag m is For N → ∞, the Fourier transformation of the autocorrelation function is the power spectrum P (ω) [16] which is directly comparable to S(ω). When each random variable μ n obeys the Gaussian distribution, is maximized, and thus h is taken as the entropy. Conversely, under the constraints of Eq. (6) with the given value of Eq. (4), we find the valueP (ω) which maximizes Eq. (7). The solution isP where M (≤ N ) is the maximum number of |m| in Eq. (4). The parameters a k are the Lagrange multipliers which are the solution of the Yule-Walker equation or equivalently, the linear Toeplitz matrix equation [16], The Yule-Waker equation Eq. (9) is also derived by the autoregressive model of where e m 's are the white noise errors [17]. In this context, the important number for M is where FPE is Akaike's final prediction error [18,19].

Optical spectra of MEM
In this section, we discuss the calculated absorption spectra of MEM compared with those of FFT.

Ethylene
A simple molecule ethylene was chosen to confirm the efficiency of MEM. The molecular structure was based on the ground state. The calculation of the time evolution was the same as the realtime TDDFT procedure described in Sec. 2.1. Figure 2   The results of MEM in Fig. 2 (f), (g) and (h) are similar to those of FFT. The resolution in Fig. 2 (h) is comparable with that in Fig. 2 (e) but the number of time steps is one half, which is the greatest advantage of MEM. MEM becomes more efficient if a suitable number M is chosen for the lag of autocorrelation. In the cases Fig. 2 (f) and (g), we used Akaike's FPE (M A = 48). In the case of Fig. 2 (h), M A was 38 but we adopted M = 1000, which was almost 25 times larger, to make the spectrum similar to that of FFT. Use of more data are in the autocorrelation calculation with larger lags and the estimation of lower-frequency components gives greater reliability. However, there are cases in which too large a value of M provides unphysical results such as pseudo-peaks and splits.

Fluorene
Poly (9,9-dialkyl-fluorene) and their substituted derivatives are expected to be the basic materials for blue emitting polymer diodes, and their electronic structures have been extensively studied [3,20,21,22,23]. We employed the oligomer of fluorene (FL) with n = 8, and the results are shown in Fig. 3. Figure 3 (a) shows the dynamic dipole moment with the total time steps N = 20000. The spectrum of oligo-FL in Fig. 3 (b) was calculated by MEM with N = 10000 and M = 1000. As in the calculation of ethylene, the absorption spectrum is comparable with that of the FFT calculation in Fig. 3 (c), where FFT took N = 20000 steps to obtain the spectrum. The broad peak at 490 nm for MEM and FFT may correspond to the experimentally observed one at 348 nm for poly-FL. This discrepancy is due to the inherent problem of DFT [5]. In this MEM calculation, we employed M = 1000, which was almost forty times larger than Akaike's FPE (M A = 26). We did not observe a significant change in the spectrum with further increased M .  Figure 4 (a) shows the difference between the MEM and FFT spectra obtained in the same time steps N = 5000 for C 60 . In this calculation of MEM, Akaike's FPE (M A = 166) was employed. As noted in Section 2.1, the FFT and thus MEM spectra reasonably simulated the experimental absorption peaks at 3.8, 4.8, and 5.8 eV.

C60
The MEM spectrum has a different intensity profile from that of FFT: the MEM intensities of the low-energy region increased and those of the high-energy region decreased in comparison with the FFT intensities. As shown in Fig. 4 (b), in the case of N = 5000 (M A = 166), there are three peaks in the region from 3 to 7 eV, whereas in the case of N = 3000 (M A = 62) and 4000 (M A = 62), there are no separated peaks there.

Concluding remarks
We confirmed that MEM is an efficient method for obtaining absorption spectra from TDDFT. In our calculation, in comparison with the conventional FFT, only a half as many steps were required to obtain the same spectral resolution, and the peak positions were as reliable as in FFT.
Since the FFT or MEM procedure is carried out after the TDDFT calculation, the reduction in computational time afforded by MEM is a remarkable advantage, for example, in screening processes to find candidate materials from thousands of trial materials. The total number of steps N and the lag of autocorrelation M are quite sensitive to the spectrum.