Singular Spectrum Analysis of Time-series Data from Time-dependent density-functional theory in Real-time

In this section, we outline the real-time TDDFT calculation procedure. The total energy of the ground state, based on density-functional theory (DFT) with the local density approximation, is derived from the Kohn–Sham (KS) equation. While DFT is inadequate for accurately describing optical responses and absorption spectra involving electronically excited states, this limitation was addressed by Runge and Gross through the introduction of time-dependent evolution of the DFT equation [1]. The equations of motion of TDDFT coupled with pseudopotentials can be written as follows:

where $\psi_{j}$ represents the $j$ th wave function. Hamiltonian $H$ is the KS Hamiltonian $H_{KS}$ to which a perturbation $V_{ext}({\bf r},t)$ is added, i.e., $H=H_{KS}+V_{ext}({\bf r},t)$. $V_{ion}^{ps}$ represents the ionic pseudopotential, $V_{H}$ is the Hartree potential, and $V_{XC}$ is the exchange-correlation potential. The Hartree potential and the exchange-correlation potential are expressed by the electronic charge density $\rho({\bf r},t)=\sum_{j}|\psi_{j}({\bf r},t)|^{2}$. The summation is performed over all occupied states $j$, and the Hartree potential is determined by $\nabla^{2}V_{H}=-4\pi\rho$. For simplicity, all calculations in this study are performed using the atomic unit system.

Prior to the calculation of the optical response, we first determine the stationary state using conventional DFT to optimize the electronic structure. Next, we added the external potential $V_{ext}$ as a perturbation to the system and tracked the linear response of the system in real-time. In our calculations, we employ the real-time, real-space approach to solve Eq.(3) using the finite difference method. This approach is advantageous due to its suitability for parallel computation. In this study, a uniform grid is used for simplicity.

The external electric field $V_{ext}(t)$ is applied as a perturbation in the form $V_{ext}(t)=E(t)\xi$ in the $\xi$-direction at $t=0$. When a very weak electric field $E(t)$ is applied, a dipole moment $\mu_{\xi}(t)$ is generated as the system’s response. The polarizability, $\alpha_{\xi}(\omega)$, which characterizes the linear response, is expressed as follows. Here, $\xi$ represents the $x,y,z$ directions.

When an external stimulus in the form of a delta function, $V_{ext}(t)=-k\xi\delta(t)$, is applied, the polarizability $\alpha_{\xi}(\omega)$ can be directly obtained from FT of the dipole moment $\mu_{\xi}(t)$

where $k$ represents the strength of the external perturbation in the $\xi$ direction. Furthermore, the imaginary part of the polarizability is used to compute the total oscillator strength $S(\omega)$

where $\alpha=(\alpha_{x}+\alpha_{y}+\alpha_{z})/3$ is the average polarizability. From Eq.(3), the time-dependent wave function $\psi_{j}({\bf r},t)$ is expressed as follows:

The initial wave function ${\tilde{\psi}}_{j}$ for this perturbation at $t=0$ is obtained from the following equation:

It becomes a very simple expression:

The time evolution of the dipole moment, $\mu_{\xi}(t)$, is numerically computed up to a total time $T$ using discrete time steps of $\Delta t$, where $T=N\Delta t$. The spectral resolution of $S(\omega)$, obtained from Eq.(6), is approximately $O(1/T)$. Therefore, insufficient data points reduce spectral resolution, requiring a sufficiently large $T$ for a clear spectrum. To address this limitation, we propose applying SSA to the time-series $\mu_{\xi}(t)$ to enable clear spectral analysis even with limited data.

The procedure for applying SSA to the time-series data $\mu_{\xi}(t)$ is described below. SSA is a method for decomposing time-series data and extracting components based on eigenvalues. A trajectory matrix is created from the time-series data, followed by Singular Value Decomposition (SVD). The resulting components were then reconstructed to identify the fundamental oscillation components in the data, which were analyzed using FT [7, 8].

From the time-series data $F=(f_{1},f_{2},\dots,f_{N})$, a trajectory matrix $Y$ of size $\tau\times n$ is created using window width $\tau(1<\tau<N/2)$, where $n=N-\tau+1$. This matrix is skew-symmetric and captures the correlation in the original time-series data:

The column vectors of $Y$ are expressed as ${\bf y}^{(1)},{\bf y}^{(2)},\dots,{\bf y}^{(n)}$. SVD is applied to this matrix, decomposing it into $Y=U\Sigma V^{T}$, where $U$ is a $\tau\times\tau$ matrix with column vectors $({\bf u}_{1},\dots,{\bf u}_{\tau})^{T}$. $\Sigma$ is a diagonal matrix of size $\tau\times n$, and each diagonal component is given by $(\sigma_{1},\dots,\sigma_{\tau})$. $V^{T}$ is also a $n\times n$ orthogonal matrix with row vectors $({\bf v}_{1},\dots,{\bf v}_{\tau})^{T}$, respectively. To remove high-frequency noise, low-rank approximations are used in SVD. Based on the singular values, $Y$ is decomposed into its cor-responding components $Y_{i}$:

where $r$ is the adopted rank. For simplicity, the product of $U_{r}$ and $\Sigma_{r}$ is defined as ${U^{\prime}}_{r}$. The $i$ th component $Y_{i}$ can be expressed as

In general, this matrix $Y_{i}$ is not skew-symmetric like the original matrix $Y$. Using the average as shown in Eq.(20), we reconstruct the following time-series data ${\tilde{F}}_{i}=({\tilde{f}}_{i,1},{\tilde{f}}_{i,2},\dots,{\tilde{f}}_{i,N})$, which corresponds to each singular value $\sigma_{i}$:

Consequently, we can obtain a matrix ${\tilde{Y}}_{i}$

The reconstructed matrix ${\tilde{Y}}_{i}$ retains symmetry similar to that of the original trajectory matrix $Y$. Now we also place the column vector ${\tilde{Y}}_{i}$ as ${\tilde{\bf y}}_{i}^{(1)},{\tilde{\bf y}}_{i}^{(2)},\dots,{\tilde{\bf y}}_{i}^{(n)}$. The original time-series can be approximated as $F\approx\Sigma_{i}{\tilde{F}}_{i}$.

In some cases, multiple $\tilde{F}_{i}$ components may exhibit oscillations of similar frequency and amplitude. These components are grouped into the same cluster and treated as a single oscillation. To classify these components into clusters, the similarity of their oscillation behavior was evaluated using the W-correlation matrix

This section explains the forecasting of the time-series data $\tilde{F}_{i}$ using SSA [9]. As mentioned above, since the time-series data $\tilde{F}_{i}$, extracted as specific oscillation components, is significantly simpler than the original $F$, it allows for stable forecasting. Forecasting $\tilde{F}_{i,N+1}$ using $\tilde{F}_{i}=(\tilde{f}_{i,1},\tilde{f}_{i,2},\dots,\tilde{f}_{i,N})$ involves adding a new column vector of $\tilde{Y}_{i}$ as shown Eq.(25).

Since only ${\tilde{f}}_{i,N+1}$ is unknown, we use the first to the $\tau-1$ of ${\tilde{Y}}_{i}$, and obtain the following linear combination of the components ${\tilde{u}}_{i,jk}$ with coefficients $h_{k}$ as follows:

where ${\tilde{u}}_{j,k}~{}(j=1,\dots,\tau-1)$ are the components of ${U^{\prime}}_{r}^{(i)}$, obtained from the decomposition as per Eq.(15). The first term in Eq.(32) is the average of each row $f_{i,a}~{}(a=1,\dots,\tau-1)$ of ${\tilde{Y}}_{i}$, and the second term is the deviation of the component. The new point, ${\tilde{f}}_{i,N+1}$, is obtained by solving Eq.(32) for $h_{k}$ using the least-squares method. The forecast is made by

By incorporating this value into ${\tilde{Y}}_{i}$ as a new component to ${\bf\tilde{y}}_{i}^{(n+1)}$, an updated ${\tilde{Y}}_{i}$ is obtained, expanded by one column. SVD is then reapplied to the updated ${\tilde{Y}}_{i}$, allowing the operations described in Eq.(32) and (33) to be repeated iteratively, enabling forecasts to be extended to the desired number of points.

To evaluate the effectiveness of SSA, we applied it to the dipole moments obtained from TDDFT calculations for ethylene, a molecule with well-characterized electronic states. The initial external stimulus was applied along the C-C bond direction, resulting in a slowly varying dipole moment $\mu(t)$. Spectral analysis of $\mu(t)$ was performed using FT as described in Eq.(5), focusing on the band-edges associated with optical absorption. Figure 1 shows the dipole moments and the resulting spectrum $S(\omega)$ for two cases: (1) sufficient data ($N\simeq 20000$) and (2) insufficient data ($N\simeq 5000$), with a time step$\Delta t=0.002$ [1/eV]. When sufficient time evolution is available, the band-edge spectra are clear, with distinct peaks around $7.5$ eV. In contrast, with insufficient data points ($N\simeq 5000$), the spectrum becomes weak and broad, making it challenging to discern whether it is a single peak or composed of multiple oscillations. Identifying peaks associated with the band-edges also becomes difficult. We then analyzed the spectra using the shorter time-evolved time-series $F=(\mu_{1},\mu_{2},\dots,\mu_{5000})$ up to $N=5000$ steps. It is necessary to select the appropriate oscillation components related to the band-edge among several obtained by decomposing the time-series $F$. To identify the major oscillation components in $F$, a trajectory matrix $Y$ is constructed with the window width $\tau=1000$ according to Eq.(14). The results are obtained in the order of the largest singular values. Using $Y_{1},Y_{2},\dots,Y_{r}$ according to Eq.(16), and by averaging according to Eq.(20), each time-series data ${\tilde{F}}_{1},{\tilde{F}}_{2},\dots$ corresponding to each singular value is obtained.

Figure 2 shows the first six oscillations. Pairs such as (${\tilde{F}}_{1}$,${\tilde{F}}_{2}$), (${\tilde{F}}_{3}$,${\tilde{F}}_{4}$), and (${\tilde{F}}_{5}$,${\tilde{F}}_{6}$) exhibit similar frequencies and amplitudes, indicating they are related oscillations. Each pair was clustered and denoted as $({\tilde{F}}_{i},{\tilde{F}}_{j})$. For these reconstructed oscillations, we counted the wavenumbers, excluding the first and the last approximately $100$ steps to avoid boundary effects. The spectra derived from this wavenumber analysis showed that the main oscillations correspond to peaks around $7.5$ eV for $({\tilde{F}}_{1},{\tilde{F}}_{2})$, $11.8$ eV for $({\tilde{F}}_{3},{\tilde{F}}_{4})$ and $18.4$ eV for $({\tilde{F}}_{5},{\tilde{F}}_{6})$. Comparing these results with Fig.1(b) confirms that these peaks are included.

The peaks near the corresponding energies are notable. In particular, $({\tilde{F}}_{1},{\tilde{F}}_{2})$ represent the lowest energy peaks, which are associated with the band-edge. These oscillations exhibit minimal beating and large amplitudes, indicating that they are relatively simple signals. In contrast, $({\tilde{F}}_{3},{\tilde{F}}_{4})$ and $({\tilde{F}}_{5},{\tilde{F}}_{6})$ show noticeable beats and are mixed with adjacent oscillations. This indicates that SSA does not always achieve a perfect decomposition into single-frequency components.

Figure 3 shows the correlation matrix $W_{ij}$ for ${\tilde{F}}_{1}-{\tilde{F}}_{10}$, calculated using Eq.(26). The diagonal components are, as expected, equal to $1$, as they represent self-correlation. If different time-series are completely separated, their off-diagonal correlations vanish. However, Fig.(3) shows high correlations within the pairs $({\tilde{F}}_{1},{\tilde{F}}_{2})$,$({\tilde{F}}_{3},{\tilde{F}}_{4})$, and $({\tilde{F}}_{5},{\tilde{F}}_{6})$, confirming their relationship. The effectiveness of this decomposition depends on the window width $\tau$, which must be optimized for proper separation. In SSA, the outputs are ordered from the highest singular value to the lowest, which means that the main oscillations may appear in the high-energy region, depending on the direction of the dipole moment oscillation. To analyze the band-edges effectively, it is useful to rearrange the decomposed and reconstructed time-series data in ascending order of wavenumber, renaming them as ${\tilde{F}}_{1},{\tilde{F}}_{2},\dots$, for clarity.

Figure 4 demonstrates the analysis and forecasting of oscillations using the method described in Section 2.3. In Fig.4(a), the time-series data $F=(\mu_{1},\mu_{2},\dots,\mu_{5000})$ is shown for the dynamic dipole moment $\mu(t)$ up to $N=5000$ steps. Fig.4(b) displays the time-series data ${\tilde{F}}={\tilde{F}}_{1}+{\tilde{F}}_{2}$, and the forecasted oscillation data using the prediction method from Section 2.3. In the figure, ${\tilde{\mu}}(t)$ represents ${\tilde{F}}=({\tilde{\mu}}_{1},{\tilde{\mu}}_{2},\dots,{\tilde{\mu}}_{5000})$. For this ${\tilde{F}}$(blue), the prediction is repeated up to $N=20000$ steps, and the result is shown in orange. The reconstructed oscillation remains almost simple and stable over a long period of time. Figure.4(c) shows the spectrum obtained by applying FT directly to the original time-series $F$ from Fig.4(a). Since time-series up to $N=5000$ steps are used, the total time $T$ is not sufficient, and the spectrum is very broad with low resolution. Focusing on the low-energy peak, the spectrum in the energy region from $6$ eV to $9$ eV is shown in Fig.4(d), with each peak normalized to the magnitude obtained from the results at $N=2000$ steps (orange). The spectra compare the FT results for ${\tilde{F}}$ using data from $N=5000$ steps (blue), $N=10000$ steps (green), and $N=20000$ steps (orange). It is evident that the signal becomes stronger and sharper as the number of steps increases.

In this section, we confirm the effect of the SSA forecast by comparing our results with those of TDDFT up to the same number of steps, using small molecules. For molecules with simple electronic structures, such as benzene, naphthalene, anthracene, and tetracene, we examined the effect of forecasting on the spectra near the band-edge, as in Section 3.1 It is well-known that the energy of the peak positions decreases as the molecular size increases.

Figure 5 compares spectra derived from the dynamic dipole moments calculated by TDDFT with those extended and forecasted using SSA. Fig.5(a) shows the spectra near the band-edge, based on the dynamic dipole moment $\mu(t)$ obtained by TDDFT up to the step number ($N=2000$). The results reveal a shift of the band-edge to lower energy as the number of benzene rings and atomic size increase. However, the spectrum in this energy region appears ambiguous, influenced by other nearby peaks.

To clarify the band-edge spectrum, we applied SSA to extract ${\tilde{\mu}}(t)$ from each time-series data at $N=2000$ steps, isolating the oscillations associated with the band-edge. The extracted fundamental oscillation was then forecasted and extended to $N=8000$ steps, as shown in Fig.5(b), with a window width of $\tau=500$. Since SSA effectively isolates the oscillation components near the band-edge, the influence of other peaks is relatively reduced. On the other hand, Fig.5(c) shows spectra obtained from the dipole moments calculated up to the step number ($N=8000$) simply by the TDDFT calculation.

Comparing the peak shapes of these spectra, it is found that they are nearly identical and equivalent.