Stochastic Real-Time Second-Order Green’s Function Theory for Neutral Excitations in Molecules and Nanostructures

We consider the electronic Hamiltonian of a finite system interacting with an explicit electric field. In second quantization the Hamiltonian is given by

where $i$, $j$, $k$, and $l$ denote indexes of a general basis, $\hat{a}^{\dagger}_{i}$ ($\hat{a}_{i}$) is the creation (annihilation) operator for an electron in orbital $\chi_{i}$, and $h_{ij}$ and $v_{ijkl}$ are the matrix elements of the one-body and two-body interactions, respectively. The two-body terms are given by the 4-index electron repulsion integral (ERI):

where we have assumed that the basis set is real. We use atomic units throughout the manuscript, where the electron charge $e=1$, the electron mass $m_{e}=1$, $\hbar=1$, the Bohr radius $a_{0}=1$, and $4\pi\epsilon_{0}=1$.

The last term in Eq. (1) is a time-dependent perturbation. Here, to describe neutral excitation, we couple the system to an external electric field, $E(t)$, within the dipole approximation, where $\Delta_{ij}(t)=E(t)\cdot{\bf\mu}_{ij}$ and

We choose to explicitly include this term in the Hamiltonian (rather than introducing a linear-response perturbation in the initial wave function) because we make no assumption that the field is weak.

Following Ref. 6, the equations of motion for the single-particle lesser Green’s function, ${\bf G}^{<}(t_{1},t_{2})$, are given by the KB equations:

and

where $t_{1}$ and $t_{2}$ are projections onto the real-time branch, $\rho(t)=-i{\bf G}^{<}(t,t)$ is the time-dependent density matrix, and ${\bf F}[\rho(t)]$ is the Fock operator, with matrix elements given by

In the above, the Hartree and exchange potentials are given by $v^{H}_{ij}[\rho]=\sum_{kl}v_{ijkl}\rho_{kl}$ and $v_{ij}^{x}[\rho]=\sum_{kl}v_{ikjl}\rho_{kl}$, respectively.

The last terms in Eqs. (4) and (5) are the collision integrals, given by[6]

and

respectively. In the above equations, ${\bf\Sigma}$ is the self-energy (which encodes all many-body interactions) and the superscript "$R/A$" denotes retarded/advanced components.

To obtain an approximate expression for the self-energies, we use the second-order Born approximation, where the self-energy is expanded to second-order in the Coulomb interaction. The resulting retarded component can be written in terms of the retarded and greater screened Coulomb integrals ($\delta W^{R/>}$)[6]

where

and

A particular feature of the self-energy in Eq. (9) is the inclusion of dynamical exchange as diagrammatically illustrated in Fig. 1.

The equations of motion (4) and (5) for the GFs together with the expression for the self-energy (Eq. (9)) form a close set of equations, but depend on two times, $t_{1}$ and $t_{2}$. To further simplify the time evolution of the GF, we assume that the retarded self-energy responds instantaneously to the application of external driving forces (e.g. the adiabatic limit)[6]

while the lesser self-energy is assumed to be negligible[6]

In the above, $\tilde{\bf\Sigma}^{\rm ad}[t]$ is the so-called adiabatic GF2 self-energy with matrix elements[6]

where

is the Fourier transform of the screened Coulomb interaction, $f(\varepsilon)$ is the Fermi-Dirac distribution, $\eta$ is a small positive regularization parameter, and $\varepsilon_{k}$ are the quasiparticle energies obtained using a stochastic GF2 for charge excitations (see Ref. 34 for more information on how to calculate the quasiparticle energies using GF2). Using Eqs. (12) and (13) for the self-energy, the time evolution of the GFs given by Eqs. (4) and (5) can be reduced to a simpler form for the equal time ($t_{1}=t_{2}\equiv t$) GF (we assume an orthogonal basis from now on)[6]

where, as before, $\rho(t)=-i{\bf G}^{<}(t,t)$ and $[A,B]=AB-BA$. Excitation energies obtained using Eq. (16) will be referred to as TD-GF2 (or TD-G0F2 when the quasi-particle energies are corrected using a single-shot, non-self-consistent GF2 [34]).

In TD-GF2, the computational limiting step is the calculation of the self-energy at time $t$, $\tilde{\bf\Sigma}^{ad}(t)$. The formal computational cost scales as $O(N_{e}^{5})$ with system size, limiting the application of TD-GF2 to small system sizes. To reduce the number of self-energy evaluations, we adopt the dynamic mode decomposition (DMD) method to describe the long-time limit of the density matrix, $\rho(t)$, as described in the next subsection. In addition, we develop a stochastic approach that reduces the scaling of computing the self-energy to $O(N_{e}^{3})$ at the account of introducing a controlled statistical error. This is described in the next section.

The dynamic mode decomposition method allows the extrapolation of the density matrix dynamics to long times without the need to further solve the equation of motion. As developed in Ref. 42, the DMD method is a data-driven model order reduction procedure used to predict the long-time nonlinear dynamics of high-dimensional systems. The method is based on Koopman’s theory [45, 46] for reduced order modeling. The general strategy is to find a few ($r$) modes $\phi_{\ell}^{ij}$ with associated frequencies $\omega^{\ell}_{ij}$ to approximate the density matrix dynamics as

with coefficients $\lambda^{\ell}_{ij}$. This model is constructed from the short-time nonlinear dynamics of the density matrix and can be seen as a finite-dimensional linear approximation to the dynamics.

We define a stochastic orbital $\theta$ as a vector in the Hilbert space of the system with random elements $\pm 1$. The average of the outer product of the stochastic vectors

represents the identity matrix, referred to as the stochastic resolution of the identity.[40] Here, $\langle\theta\otimes\theta^{T}\rangle_{N_{s}}\equiv\frac{1}{N_{s}}\sum_{\xi=1}^{N_{s}}\theta_{\xi}\otimes\theta_{\xi}^{T}$ is an average over the set $\{\theta_{\xi}\}$ of uncorrelated stochastic orbitals $\theta_{\xi}$, with $\xi=1,2,\dots,N_{s}$.

Analogous to the deterministic resolution of the identity (also known as density fitting), in which 3-index $(ij|A)$ and 2-index $V_{AB}=(A|B)$ ERIs are used to approximate the 4-index ERI as

where $A(B)$ is an auxiliary basis of dimension $N_{\text{aux}}$, the stochastic resolution of the identity can be used as a resolution basis to approximate the 4-index ERIs as [40]

where $R_{\alpha\beta}=\sum_{A}^{N_{aux}}(\alpha\beta|A)\sum_{B}^{N_{aux}}V_{AB}^{-1/2}\theta_{B}$. One advantage of using this approximation is that the indexes $ij$ and $kl$ are decoupled, allowing to perform tensor contractions and reduce the computational scaling.[40, 47]

The use of the stochastic resolution of the identity to approximate the ERIs introduces a controllable statistical error that can be tuned by changing the number of stochastic orbitals, with a convergence rate proportional to $1/\sqrt{N_{s}}$. An alternative for controlling the error is to use the range-separated stochastic resolution of the identity in which the largest contributions to the ERIs are treated deterministically, while the remaining terms are treated stochastically. Specifically, as proposed in Ref. 39, we first identify large contributions (denoted by the superscript $L$) to the 3-index ERIs with respect to a preset threshold,

where $\varepsilon^{\prime}$ is a parameter in the range $[0,N_{\rm e}]$. The factor $\frac{\varepsilon^{\prime}}{N}$ guarantees that the total non-vanishing elements in $(ij|A)^{L}$ scales as $O_{\rm e}^{2})$. Then, we define the large K-tensors

and keep only their larger elements, according to a second threshold

in which $\varepsilon$ is a parameter in the range $[0,1]$. We then define large and small (denoted by the superscript $S$) R-tensors as

and

where $R_{ij}$ is defined as in Eq. (20). Using these expressions, a range-separated 4-index ERI can be written as

To derive a stochastic expression for the self-energy we insert Eq. (20) (or Eq. (26) for range-separated computations) into Eqs. (14) and (15), to yield

where $\delta\rho(t)=\rho(t)-\rho(t_{0})$. In the above equation, the "prime" superscript denotes that a different set of stochastic orbitals is used to construct the $R^{\prime}$-tensors. Next, we rearrange the equation of motion for the density matrix (cf., Eq. (16)) as:

where $\Delta H=\Sigma(t_{0})$ is the GF2 (or G0F2) quasiparticle energy correction. Excitation energies obtained using Eq. (28) in combination with Eq. (27) will be referred to as sTD-GF2 (or sTD-G0F2).

Fig. 2 shows the stochastic (sTD-G0F2) and deterministic (TD-G0F2) induced dipole dynamics of a hydrogen dimer chain (H₂₀, containing ten H₂ dimers with bond length of 0.74 Å  and intermolecular distance of 1.26 Å, align along the $z$-axis). The equations of motion for the density matrix were propagated using Eq. (28) with stochastic (Eq. (27)) and deterministic (Eq. (14)) self energies, respectively. In both cases, a Gaussian-pulse centered at $t_{0}=1$ fs, was used to represent the electric field, with an amplitude $\gamma E_{0}=0.02$ V/Å and a variance of 0.005 fs; the regularization parameter appearing in Eq. (15) $\eta=0.01$ and the inverse temperature is set to $\beta=50$ in all computations. For the stochastic approach, averages were computed using $N_{s}=80$ stochastic orbitals. In all cases, the minimal basis set STO-3G was used.

Fig. 2(a) exemplifies how the stochastic approach reproduces the deterministic dynamics, by comparing the TD-G0F2 and sTD-G0F2 induced dipole dynamics for the H₂₀ chain. The shaded region in red is the standard deviation (SD) obtained from $6$ independent runs. The statistical error can be reduced by increasing the number of stochastic orbitals (with a convergence rate proportional to $1/\sqrt{N_{s}}$) or by changing the range separation parameters $\varepsilon$ and $\varepsilon^{\prime}$, as is further discussed in Sec. IV.3 below.

We find that for a fixed number of stochastic orbitals and for fixed values of $\varepsilon$ and $\varepsilon^{\prime}$, the statistical error increases with the propagation time, as shown in Fig. 2(b). This is consistent with our previous finding for the stochastic time-dependent density functional theory [48] and for the stochastic BSE approach.[49] Since the induced dipole decays rather rapidly in time, the increase in the statistical error at long times does not affect the absorption spectrum in any significant way. Nonetheless, to mitigate the divergence of the dynamics at long times, we have multiplied the induced dipole by a damping function $e^{-10t/t_{\text{max}}}$, where $t_{\text{max}}$ corresponds to the total propagation time and plays a similar role as the regularization parameter, $\eta$.

The long-time dynamics of the density matrix and the resultant time-dependent induced dipole were obtained using the DMD technique outlined above. Fig. 2(c) shows a comparison between the extrapolated DMD dynamics and the dynamics obtained by solving Eq. (28) for both the deterministic and stochastic methods. The shaded region (first $6$ fs in Fig. 2(c)) indicates the portion of the dynamics that was used to train the reduced DMD model (defined as DMD window), while the remaining $34$ fs (of which 14 fs are shown in Fig. 2(c)) corresponds to the extrapolated dynamics. We find that the DMD technique accurately captures the main dynamical features, even for the noisy stochastic data. Naturally, the time scale of the events that can be captured by the reduced DMD model depends on the DMD window length. Below, we analyze the accuracy of the DMD approach in reproducing the absorption spectra.

In Fig. 3 panels (a) and (b) we compare the absorption spectra obtained from the stochastic and deterministic real-time dynamics and the reference deterministic frequency-domain linear-response approach (G0F2-BSE), for two representative hydrogen dimer chains. The absorption spectra obtained from the three different approaches (vertically shifted for clarity) are numerically identical, demonstrating that the real-time implementations are consistent with the frequency domain reference methods (in the weak coupling-linear response limit) and, in particular, that the stochastic approach can reproduce the benchmark results with only $80$ stochastic orbitals.

The frequency resolution of the absorption spectra can be improved by propagating the density matrix dynamics to longer times using the DMD technique. Fig. 3 panels (c) and (d) show the corresponding absorption spectra for a $40$ fs extrapolated trajectory (sTD-G0F2+DMD) considering three different DMD window lengths. Even a short ($2$fs) window length provides a reasonable description of the low-excitation features (main absorption peak at $\sim$15 eV), but the quality of the spectra improves with increasing DMD windows lengths, especially for the higher-excitation peaks. Specifically, for the sTD-G0F2+DMD method, we observed that the average DMD spectral error is proportional to $1/\sqrt{t_{\text{DMD}}}$, with $t_{\text{DMD}}$ being the DMD window length.

The variation of the range-separated threshold parameters, $\varepsilon$ and $\varepsilon^{\prime}$ (see Eqs. (21) and (23)), allows us to control the ratio of deterministic to stochastic Coulomb tensor elements. As $\varepsilon^{\prime}\to N$ or $\varepsilon\to 1$ the approach reduces to the fully stochastic limit. By contrast, when $\varepsilon^{\prime}\to 0$ or $\varepsilon\to 0$ the approach is fully deterministic. Fig. 4(a) shows the dependence of the statistical error on $\varepsilon^{\prime}$ and $\varepsilon$ for H₂₀. The average error was estimated using $n=6$ independent stochastic runs as

where $N_{\omega}$ is the number of frequencies used in the range $E=10-30$ eV. As $\varepsilon$ increases the statistical error increases and approaches the fully stochastic limit (dotted line in Fig. 4(a)). For the case with the lowest statistical error in Fig. 4(a) ($\varepsilon^{\prime}=0.002$, $\varepsilon=0.001$), the amount of ERI elements computed deterministically corresponds to $\approx 10\%$ for H₂₀, resulting in an error reduction of almost $2$ orders of magnitude compared to the fully stochastic limit.

The main advantage of using the stochastic formulation of GF2 in the real-time domain is the reduction in the computational complexity and scaling. Formally, GF2 in the frequency-domain scales as $O(N_{\rm e}^{6})$ with the system size ($N$) while the real-time deterministic implementation scales as $O_{\rm e}^{5})$. By contrast, when the stochastic resolution of identity is used in the time-domain, the computational scaling is further reduced to $O_{\rm e}^{3})$, as long as the number of stochastic orbitals does not increase with system size to achieve a similar statistical error (which is the case for the systems studied here). The computational limiting step in the sTD-GF2 method is the computation of the self-energy (Eq. (27)), with a formal scaling of $O(N_{s}N_{\rm e}^{3})$ when appropriate tensor contractions are used. Figure 4(b) shows the computational cost associated with the stochastic and deterministic real-time methods and the equivalent frequency-domain linear-response implementation for hydrogen dimer chains with varying lengths. The lowest scaling corresponds to the stochastic real-time implementation, sTD-GF2, which exhibits an $O_{\rm e}^{3})$ behavior, with a large pre-factor. For the current target statistical error, the stochastic approach is computationally more efficient than the deterministic approach for system sizes that exceed $N\approx 200$ basis functions.