This commit is contained in:
Pierre-Francois Loos 2020-07-22 14:20:55 +02:00
parent fbdf708629
commit ef219c6753
2 changed files with 203 additions and 197 deletions

View File

@ -1,13 +1,21 @@
%% This BibTeX bibliography file was created using BibDesk.
%% http://bibdesk.sourceforge.net/
%% Created for Pierre-Francois Loos at 2020-06-26 09:45:06 +0200
%% Created for Pierre-Francois Loos at 2020-07-22 10:40:20 +0200
%% Saved with string encoding Unicode (UTF-8)
@article{Loos_2020f,
Author = {P. F. Loos and J. Authier},
Date-Added = {2020-07-22 10:38:59 +0200},
Date-Modified = {2020-07-22 10:40:20 +0200},
Journal = {in preparation},
Title = {Dynamical Kernels for Optical Excitations},
Year = {2020}}
@article{Petersilka_1996,
Author = {M. Petersilka and U. J. Gossmann and and E. K. U. Gross},
Date-Added = {2020-06-26 09:43:33 +0200},
@ -17,7 +25,8 @@
Pages = {1212},
Title = {Excitation Energies From Time-Dependent Density-Functional Theory},
Volume = {76},
Year = {1996}}
Year = {1996},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevLett.76.1212}}
@article{Nielsen_1980,
Author = {Egon S. Nielsen and Poul Jorgensen},
@ -14426,14 +14435,12 @@
Bdsk-Url-2 = {https://doi.org/10.1103/PhysRevB.93.235113}}
@article{Boulanger_2014,
author = {Boulanger, Paul and Jacquemin, Denis and Duchemin, Ivan and Blase, Xavier},
title = {Fast and Accurate Electronic Excitations in Cyanines with the Many-Body Bethe-Salpeter Approach},
journal = {J. Chem. Theory Comput.},
volume = {10},
number = {3},
pages = {1212--1218},
year = {2014},
doi = {10.1021/ct401101u},
}
Author = {Boulanger, Paul and Jacquemin, Denis and Duchemin, Ivan and Blase, Xavier},
Doi = {10.1021/ct401101u},
Journal = {J. Chem. Theory Comput.},
Number = {3},
Pages = {1212--1218},
Title = {Fast and Accurate Electronic Excitations in Cyanines with the Many-Body Bethe-Salpeter Approach},
Volume = {10},
Year = {2014},
Bdsk-Url-1 = {https://doi.org/10.1021/ct401101u}}

View File

@ -224,7 +224,7 @@ Our calculations are benchmarked against high-level (coupled-cluster) calculatio
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%
The Bethe-Salpeter equation (BSE) formalism \cite{Salpeter_1951,Strinati_1988} is to the $GW$ approximation \cite{Hedin_1965,Golze_2019} of many-body perturbation theory (MBPT) \cite{Martin_2016} what time-dependent density-functional theory (TD-DFT) \cite{Runge_1984,Casida_1995} is to Kohn-Sham density-functional theory (KS-DFT), \cite{Hohenberg_1964,Kohn_1965} an affordable way of computing the neutral excitations of a given electronic system.
The Bethe-Salpeter equation (BSE) formalism \cite{Salpeter_1951,Strinati_1988} is to the $GW$ approximation \cite{Hedin_1965,Golze_2019} of many-body perturbation theory (MBPT) \cite{Martin_2016} what time-dependent density-functional theory (TD-DFT) \cite{Runge_1984,Casida_1995} is to Kohn-Sham density-functional theory (KS-DFT), \cite{Hohenberg_1964,Kohn_1965} an affordable way of computing the neutral (or optical) excitations of a given electronic system.
In recent years, it has been shown to be a valuable tool for computational chemists with a large number of systematic benchmark studies on large families of molecular systems appearing in the literature \cite{Boulanger_2014,Jacquemin_2015a,Bruneval_2015,Jacquemin_2015b,Hirose_2015,Jacquemin_2017a,Jacquemin_2017b,Rangel_2017,Krause_2017,Gui_2018} (see Ref.~\onlinecite{Blase_2018} for a recent review).
Taking the optical gap (\ie, the lowest optical excitation energy) as an example, BSE builds on top of a $GW$ calculation by adding up excitonic effects $\EB$ to the $GW$ HOMO-LUMO gap
@ -243,8 +243,8 @@ where
\begin{equation} \label{eq:Egfun}
\EgFun = I^N - A^N
\end{equation}
is the the fundamental gap, \cite{Bredas_2014} $I^N = E_0^{N-1} - E_0^N$ and $A^N = E_0^{N+1} - E_0^N$ being the ionization potential and the electron affinity of the $N$-electron system.
Here, $E_s^{N}$ is the total energy of the $s$th excited state of the $N$-electron system, and $E_0^N$ corresponds to its ground-state energy.
is the the fundamental gap, \cite{Bredas_2014} $I^N = E_0^{N-1} - E_0^N$ and $A^N = E_0^{N} - E_0^{N+1}$ being the ionization potential and the electron affinity of the $N$-electron system, respectively.
Here, $E_S^{N}$ is the total energy of the $S$th excited state of the $N$-electron system, and $E_0^N$ corresponds to its ground-state energy.
Because the excitonic effect corresponds physically to the stabilization implied by the attraction of the excited electron and its hole left behind, we have $\EgOpt < \EgFun$.
Due to the smaller amount of screening in molecules as compared to solids, a faithful description of excitonic effects is paramount in molecular systems.
@ -269,7 +269,7 @@ In these two latter studies, they also followed a (non-self-consistent) perturba
It is important to note that, although all the studies mentioned above are clearly going beyond the static approximation of BSE, they are not able to recover additional excitations as the perturbative treatment makes ultimately the BSE kernel static.
However, it does permit to recover, for transitions with a dominant single-excitation character, additional relaxation effects coming from higher excitations.
These higher excitations would be explicitly present in the BSE Hamiltonian by ``unfolding'' the dynamical BSE kernel, and one would recover a linear eigenvalue problem with, nonetheless, a much larger dimension.
These higher excitations would be explicitly present in the BSE Hamiltonian by ``unfolding'' the dynamical BSE kernel, and one would recover a linear eigenvalue problem with, nonetheless, a much larger dimension. \cite{Loos_2020f}
Based on a simple two-level model which permits to analytically solve the dynamical equations, Romaniello and coworkers \cite{Romaniello_2009b,Sangalli_2011} evidenced that one can genuinely access additional excitations by solving the non-linear, frequency-dependent eigenvalue problem.
For this particular system, it was shown that a BSE kernel based on the random-phase approximation (RPA) produces indeed double excitations but also unphysical excitations. \cite{Romaniello_2009b}
@ -284,7 +284,7 @@ In this regard, MBPT provides key insights about what is missing in adiabatic TD
In the present study, we extend the work of Rohlfing and coworkers \cite{Rohlfing_2000,Ma_2009a,Ma_2009b,Baumeier_2012b} by proposing a renormalized first-order perturbative correction to the static BSE excitation energies.
Importantly, our correction goes beyond the plasmon-pole approximation as the dynamical screening of the Coulomb interaction is computed exactly.
Moreover, we investigate quantitatively the effect of the TDA by computing both the resonant and anti-resonant dynamical corrections to the BSE excitation energies.
%Moreover, we investigate quantitatively the effect of the TDA by computing both the resonant and anti-resonant dynamical corrections to the BSE excitation energies.
Unless otherwise stated, atomic units are used.
%%%%%%%%%%%%%%%%%%%%%%%%
@ -293,7 +293,7 @@ Unless otherwise stated, atomic units are used.
%%%%%%%%%%%%%%%%%%%%%%%%
In this Section, following Strinati's seminal work, \cite{Strinati_1988} we first derive in some details the theoretical foundations leading to the dynamical BSE.
Additional details about this derivation are provided as {\SI}.
%Additional details about this derivation are provided as {\SI}.
We present, in a second step, the perturbative implementation of the dynamical correction as compared to the standard static approximation.
%================================
@ -306,7 +306,7 @@ The two-body correlation function $L(1,2; 1',2')$ --- a central quantity in the
\end{equation}
where, \eg, $1 \equiv (\bx_1 t_1)$ is a space-spin plus time composite variable.
The relation between $G$ and the one-body charge density $\rho(1) = -i G(1,1^+)$ provides a direct connection with the density-density susceptibility $\chi(1,2) = L(1,2;1^+,2^+)$ at the core of TD-DFT.
(The notation $1^+$ means that the time $t_1$ is taken at $t_1^{+} = t_1 + 0^+$, where $0^+$ is a small positive infinitesimal.)
(The notation $1^+$ means that the time $t_1$ is taken at $t_1^{+} = t_1 + 0^+$, where $0^+$ is a positive infinitesimal.)
The two-body correlation function $L$ satisfies the self-consistent BSE \cite{Strinati_1988}
\begin{multline} \label{eq:BSE}
@ -345,7 +345,7 @@ is the BSE kernel that takes into account the self-consistent variation of the H
[where $\delta$ is Dirac's delta function and $v$ is the bare Coulomb operator] and the xc self-energy $ \Sigma_\text{xc}$ with respect to the variation of $G$.
In Eqs.~\eqref{eq:G1} and \eqref{eq:G2}, the field operators $\Hat{\psi}(\bx t)$ and $\Hat{\psi}^{\dagger}(\bx't')$ remove and add (respectively) an electron to the $N$-electron ground state $\ket{N}$ in space-spin-time positions ($\bx t$) and ($\bx't'$), while $T$ is the time-ordering operator.
The resolution of the dynamical BSE starts with the expansion of $L_0$ and $L$ [see Eqs.~\eqref{eq:L0} and \eqref{eq:L}] over the complete orthonormalized set of $N$-electron excited states $\ket{N,s}$ (with $\ket{N,0} \equiv \ket{N}$). \cite{Strinati_1988}
The resolution of the dynamical BSE starts with the expansion of $L_0$ and $L$ [see Eqs.~\eqref{eq:L0} and \eqref{eq:L}] over the complete orthonormalized set of $N$-electron excited states $\ket{N,S}$ (with $\ket{N,0} \equiv \ket{N}$). \cite{Strinati_1988}
In the optical limit of instantaneous electron-hole creation and destruction, imposing $t_{2'} = t_2^+$ and $t_{1'} = t_1^+$, and using the relation between the field operators in their time-dependent (Heisenberg) and time-independent (Schr\"{o}dinger) representations, \eg,
\begin{equation} \label{Eisenberg}
\hpsi(1) = e^{ i \hH t_1 } \hpsi(\bx_1) e^{-i \hH t_1 },
@ -354,32 +354,32 @@ In the optical limit of instantaneous electron-hole creation and destruction, im
\begin{equation}
\begin{split}
iL(1,2; 1',2')
& = \theta(+\tau_{12}) \sum_{s > 0} \chi_s(\bx_1,\bx_{1'}) \tchi_s(\bx_2,\bx_{2'}) e^{ - i \Om{s}{} \tau_{12} }
& = \theta(+\tau_{12}) \sum_{s > 0} \chi_S(\bx_1,\bx_{1'}) \tchi_S(\bx_2,\bx_{2'}) e^{ - i \Om{S}{} \tau_{12} }
\\
& - \theta(-\tau_{12}) \sum_{s > 0} \chi_s(\bx_2,\bx_{2'}) \tchi_s(\bx_1,\bx_{1'}) e^{ + i \Om{s}{} \tau_{12} },
& - \theta(-\tau_{12}) \sum_{s > 0} \chi_S(\bx_2,\bx_{2'}) \tchi_S(\bx_1,\bx_{1'}) e^{ + i \Om{S}{} \tau_{12} },
\end{split}
\end{equation}
where $\tau_{12} = t_1 - t_2$, $\theta$ is the Heaviside step function, and
\begin{subequations}
\begin{align}
\chi_s(\bx_1,\bx_{2}) & = \mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{2})] }{N,s},
\chi_S(\bx_1,\bx_{2}) & = \mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{2})] }{N,S},
\\
\tchi_s(\bx_1,\bx_{2}) & = \mel{N,s}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{2})] }{N}.
\tchi_S(\bx_1,\bx_{2}) & = \mel{N,S}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{2})] }{N}.
\end{align}
\end{subequations}
The $\Om{s}{}$'s are the neutral excitation energies of interest.
The $\Om{S}{}$'s are the neutral excitation energies of interest.
Picking up the $e^{+i \Om{s}{} t_2 }$ component of both $L(1,2; 1',2')$ and $L(6,2;5,2')$, simplifying further by $\tchi_s(\bx_2,\bx_{2'})$ on both sides of the BSE [see Eq.~\eqref{eq:BSE}], we seek the $e^{-i \Om{s}{} t_1 }$ Fourier component associated with the right-hand side of a modified dynamical BSE, which reads
Picking up the $e^{+i \Om{S}{} t_2 }$ component of both $L(1,2; 1',2')$ and $L(6,2;5,2')$, simplifying further by $\tchi_S(\bx_2,\bx_{2'})$ on both sides of the BSE [see Eq.~\eqref{eq:BSE}], we seek the $e^{-i \Om{S}{} t_1 }$ Fourier component associated with the right-hand side of a modified dynamical BSE, which reads
\begin{multline} \label{eq:BSE_2}
\mel{N}{T [ \hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}') ] } {N,s} e^{ - i \Om{s}{} t_1 }
\mel{N}{T [ \hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}') ] } {N,S} e^{ - i \Om{S}{} t_1 }
\theta ( \tau_{12} )
\\
= \int d3456 \, L_0(1,4;1',3) \Xi(3,5;4,6)
\\
\times \mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}
\times \mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,S}
\theta [\min(t_5,t_6) - t_2].
\end{multline}
For the lowest neutral excitation energies falling in the fundamental gap of the system (\ie, $\Om{s}{} < \EgFun$ due to excitonic effects), $L_0(1,2;1',2')$ cannot contribute to the $e^{-i \Om{s}{} t_1 }$ response term since its lowest excitation energy is precisely the fundamental gap [see Eq.~\eqref{eq:Egfun}].
For the neutral excitation energies falling in the fundamental gap of the system (\ie, $\Om{S}{} < \EgFun$ due to excitonic effects), $L_0(1,2;1',2')$ cannot contribute to the $e^{-i \Om{S}{} t_1 }$ response term since its lowest excitation energy is precisely the fundamental gap [see Eq.~\eqref{eq:Egfun}].
Consequently, special care has to be taken for high-lying excited states (like core or Rydberg excitations) where additional terms have to be taken into account (see Refs.~\onlinecite{Strinati_1982,Strinati_1984}).
Dropping the space/spin variables, the Fourier components with respect to $t_1$ of $L_0(1,4;1',3)$ reads
@ -395,29 +395,29 @@ We now adopt the Lehman representation of the one-body Green's function in the q
\end{equation}
where $\eta$ is a positive infinitesimal and $\mu$ is the chemical potential.
The $\e{p}$'s in Eq.~\eqref{eq:G-Lehman} are quasiparticle energies (\ie, proper addition/removal energies) and the $\MO{p}(\bx)$'s are their associated one-body (spin)orbitals.
In the following, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, while $p$ and $q$ indicate arbitrary orbitals.
Projecting the Fourier component $L_0(\bx_1,4;\bx_{1'},3; \omega_1 = \Om{s}{} )$ onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$ yields
In the following, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, while $p$, $q$, $r$, and $s$ indicate arbitrary orbitals.
Projecting the Fourier component $L_0(\bx_1,4;\bx_{1'},3; \omega_1 = \Om{S}{} )$ onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$ yields
\begin{multline} \label{eq:iL0bis}
\iint d\bx_1 d\bx_{1'} \, \MO{a}^*(\bx_1) \MO{i}(\bx_{1'}) L_0(\bx_1,4;\bx_{1'},3; \Om{s}{})
\iint d\bx_1 d\bx_{1'} \, \MO{a}^*(\bx_1) \MO{i}(\bx_{1'}) L_0(\bx_1,4;\bx_{1'},3; \Om{S}{})
\\
=
\frac{ \MO{a}^*(\bx_3) \MO{i}(\bx_4) e^{i \Om{s}{} t^{34} }} { \Om{s}{} - ( \e{a} - \e{i} ) + i \eta }
\qty[ \theta( \tau_{34} ) e^{i \qty( \e{i} + \frac{\Om{s}{}}{2}) \tau_{34} } + \theta( - \tau_{34} ) e^{i \qty(\e{a} - \frac{\Om{s}{}}{2}) \tau_{34} } ].
\frac{ \MO{a}^*(\bx_3) \MO{i}(\bx_4) e^{i \Om{S}{} t^{34} }} { \Om{S}{} - ( \e{a} - \e{i} ) + i \eta }
\qty[ \theta( \tau_{34} ) e^{i \qty( \e{i} + \frac{\Om{S}{}}{2}) \tau_{34} } + \theta( - \tau_{34} ) e^{i \qty(\e{a} - \frac{\Om{S}{}}{2}) \tau_{34} } ].
\end{multline}
with $t^{34} = (t_3 + t_4)/2$ and $\tau_{34} = t_3 -t_4$.
More details are provided in the Appendix.
As a final step, we express the terms $\mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}')] }{N,s}$ and $\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}$ from Eq.~\eqref{eq:BSE_2} in the standard electron-hole product (or single-excitation) space.
%with $t^{34} = (t_3 + t_4)/2$ and $\tau_{34} = t_3 -t_4$.
More details are provided in Appendix \ref{app:A}.
As a final step, we express the terms $\mel{N}{T [\hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}')] }{N,S}$ and $\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,S}$ from Eq.~\eqref{eq:BSE_2} in the standard electron-hole product (or single-excitation) space.
This is done by expanding the field operators over a complete orbital basis of creation/destruction operators.
For example, we have (see derivation in the Appendix)
For example, we have (see derivation in Appendix \ref{app:B})
\begin{multline} \label{eq:spectral65}
\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}
\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,S}
\\
= - \qty( e^{ -i \Omega_s t^{65} } ) \sum_{pq} \MO{p}(\bx_6) \MO{q}^*(\bx_5)
\mel{N}{\ha_q^{\dagger} \ha_p}{N,s}
= - \qty( e^{ -i \Om{S}{} t^{65} } ) \sum_{pq} \MO{p}(\bx_6) \MO{q}^*(\bx_5)
\mel{N}{\ha_q^{\dagger} \ha_p}{N,S}
\\
\times \qty[ \theta( \tau_{65} ) e^{- i \qty( \e{p} - \frac{\Om{s}{}}{2} ) \tau_{65} } + \theta( - \tau_{65} ) e^{ - i \qty( \e{q} + \frac{\Om{s}{}}{2}) \tau_{65} } ],
\times \qty[ \theta( \tau_{65} ) e^{- i \qty( \e{p} - \frac{\Om{S}{}}{2} ) \tau_{65} } + \theta( - \tau_{65} ) e^{ - i \qty( \e{q} + \frac{\Om{S}{}}{2}) \tau_{65} } ],
\end{multline}
with $t^{65} = (t_5 + t_6)/2$ and $\tau_{65} = t_6 -t_5$. The $\mel{N}{\ha_q^{\dagger} \ha_p}{N,s}$ are the unknown particle-hole amplitudes.% in the orbital product basis.
with $t^{65} = (t_5 + t_6)/2$ and $\tau_{65} = t_6 -t_5$. The $\mel{N}{\ha_q^{\dagger} \ha_p}{N,S}$ are the unknown particle-hole amplitudes.% in the orbital product basis.
%================================
@ -435,35 +435,35 @@ leads to the following simplified BSE kernel
where $W$ is the dynamically-screened Coulomb operator.
The $GW$ quasiparticle energies $\eGW{p}$ are usually good approximations to the removal/addition energies $\e{p}$ introduced in Eq.~\eqref{eq:G-Lehman}.
Substituting Eqs.~\eqref{eq:iL0bis}, \eqref{eq:spectral65}, and \eqref{eq:Xi_GW} into Eq.~\eqref{eq:BSE_2}, and projecting onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$, one gets after a few tedious manipulations (see {\SI}) the dynamical BSE:
Substituting Eqs.~\eqref{eq:iL0bis}, \eqref{eq:spectral65}, and \eqref{eq:Xi_GW} into Eq.~\eqref{eq:BSE_2}, and projecting onto $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$, one gets after a few tedious manipulations the dynamical BSE:
\begin{equation} \label{eq:BSE-final}
\begin{split}
( \eGW{a} - \eGW{i} - \Om{s}{} ) X_{ia,s}
& + \sum_{jb} \qty[ \kappa \ERI{ia}{jb} - \widetilde{W}_{ij,ab}(\Om{s}{}) ] X_{jb,s} \\
& + \sum_{jb} \qty[ \kappa \ERI{ia}{bj} - \widetilde{W}_{ib,aj}(\Om{s}{}) ] Y_{jb,s}
( \eGW{a} - \eGW{i} - \Om{S}{} ) X_{ia,S}
& + \sum_{jb} \qty[ \kappa \ERI{ia}{jb} - \widetilde{W}_{ij,ab}(\Om{S}{}) ] X_{jb,S} \\
& + \sum_{jb} \qty[ \kappa \ERI{ia}{bj} - \widetilde{W}_{ib,aj}(\Om{S}{}) ] Y_{jb,S}
= 0,
\end{split}
\end{equation}
with $X_{jb,s} = \mel{N}{\ha_j^{\dagger} \ha_b}{N,s}$ and $Y_{jb,s} = \mel{N}{\ha_b^{\dagger} \ha_j}{N,s}$, and where $\kappa = 2 $ or $0$ for singlet and triplet excited states (respectively).
Neglecting the anti-resonant terms, $Y_{jb,s}$, in the dynamical BSE, which are (usually) much smaller than their resonant counterparts, $X_{jb,s}$, leads to the well-known Tamm-Dancoff approximation (TDA).
with $X_{jb,S} = \mel{N}{\ha_j^{\dagger} \ha_b}{N,S}$ and $Y_{jb,S} = \mel{N}{\ha_b^{\dagger} \ha_j}{N,S}$, and where $\kappa = 2 $ or $0$ for singlet and triplet excited states (respectively).
Neglecting the anti-resonant terms, $Y_{jb,S}$, in the dynamical BSE, which are (usually) much smaller than their resonant counterparts, $X_{jb,S}$, leads to the well-known TDA.
In Eq.~\eqref{eq:BSE-final},
\begin{equation}
\ERI{ia}{jb} = \iint d\br d\br' \, \MO{i}^*(\br) \MO{a}(\br) v(\br -\br') \MO{j}^*(\br') \MO{b}(\br'),
\ERI{pq}{rs} = \iint d\br d\br' \, \MO{p}^*(\br) \MO{q}(\br) v(\br -\br') \MO{r}^*(\br') \MO{s}(\br'),
\end{equation}
are the bare two-electron integrals in the spatial orbital basis $\lbrace \MO{p}(\br{}) \rbrace$, and
\begin{multline} \label{eq:wtilde}
\widetilde{W}_{ij,ab}(\Om{s}{})
= \frac{ i }{ 2 \pi} \int d\omega \; e^{-i \omega 0^+ } W_{ij,ab}(\omega)
\widetilde{W}_{pq,rs}(\Om{S}{})
= \frac{ i }{ 2 \pi} \int d\omega \; e^{-i \omega 0^+ } W_{pq,rs}(\omega)
\\
\times \qty[ \frac{1}{ \Omega_{ib}^s - \omega + i \eta } + \frac{1}{ \Omega_{ja}^{s} + \omega + i\eta } ],
\times \qty[ \frac{1}{ \Om{ps}{S} - \omega + i \eta } + \frac{1}{ \Om{qr}{S} + \omega + i\eta } ],
\end{multline}
is an effective dynamically-screened Coulomb potential, \cite{Romaniello_2009b} where $\Om{pq}{s} = \Om{s}{} - ( \eGW{q} - \eGW{p} )$ and
is an effective dynamically-screened Coulomb potential, \cite{Romaniello_2009b} where $\Om{pq}{S} = \Om{S}{} - ( \eGW{q} - \eGW{p} )$ and
\begin{equation}
W_{ij,ab}({\omega})
= \iint d\br d\br' \, \MO{i}(\br) \MO{j}^*(\br) W(\br ,\br'; \omega) \MO{a}^*(\br') \MO{b}(\br').
W_{pq,rs}({\omega})
= \iint d\br d\br' \, \MO{p}(\br) \MO{q}^*(\br) W(\br ,\br'; \omega) \MO{r}^*(\br') \MO{s}(\br').
\end{equation}
%\xavier{\sout{ A second coupled equation for the $(X_{ia}^{s}, Y_{ia}^{s} )$ vector can be obtained by projecting $\mel{N}{T [ \hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}') ] } {N,s}$ and $L_0(\bx_1,4;\bx_{1'},3; \Om{s}{})$ onto $\MO{i}^*(\bx_1) \MO{a}(\bx_{1'})$ instead of $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$. } }
%\xavier{\sout{ A second coupled equation for the $(X_{ia}^{S}, Y_{ia}^{S} )$ vector can be obtained by projecting $\mel{N}{T [ \hpsi(\bx_1) \hpsi^{\dagger}(\bx_{1}') ] } {N,S}$ and $L_0(\bx_1,4;\bx_{1'},3; \Om{S}{})$ onto $\MO{i}^*(\bx_1) \MO{a}(\bx_{1'})$ instead of $\MO{a}^*(\bx_1) \MO{i}(\bx_{1'})$. } }
%================================
@ -519,38 +519,42 @@ The RPA matrices $\bA{\RPA}$ and $\bB{\RPA}$ in Eq.~\eqref{eq:LR-RPA} are of siz
The analysis of the poles of the integrand in Eq.~\eqref{eq:wtilde} yields
\begin{multline}
\widetilde{W}_{ij,ab}( \Om{s}{} )
\widetilde{W}_{ij,ab}( \Om{S}{} )
= \ERI{ij}{ab} + 2 \sum_m \sERI{ij}{m} \sERI{ab}{m}
\\
\times \qty[ \frac{1}{\Om{ib}{s} - \Om{m}{\RPA} + i\eta} + \frac{1}{\Om{ja}{s} - \Om{m}{\RPA} + i\eta} ].
\times \qty[ \frac{1}{\Om{ib}{S} - \Om{m}{\RPA} + i\eta} + \frac{1}{\Om{ja}{S} - \Om{m}{\RPA} + i\eta} ].
\end{multline}
One can verify that, in the static limit where $\Om{m}{\RPA} \to \infty$, the matrix elements $\widetilde{W}_{ij,ab}$ correctly reduce to the static ${W}_{ij,ab}$ ones, evidencing that the standard static BSE problem is recovered from the present dynamical formalism in this limit.
Due to excitonic effects, the lowest BSE excitation energy, $\Om{1}{}$, stands lower than the lowest RPA excitation energy, $\Om{1}{\RPA}$, so that, $\Om{ib}{s} - \Om{m}{\RPA} < 0 $ and $\widetilde{W}_{ij,ab}(\Om{s}{})$ has no resonances.
Furthermore, $\Om{ib}{s}$ and $\Om{ja}{s}$ are necessarily negative quantities for in-gap low-lying BSE excitations.
Thus, we have $\abs*{\Omega_{ib}^{s} - \Om{m}{\RPA}} > \Omega_m^{\RPA}$.
One can verify that, in the static limit where $\Om{m}{\RPA} \to \infty$, the matrix elements $\widetilde{W}_{ij,ab}$ correctly reduce to their static expression
\begin{equation}
\label{eq:Wstat}
\W{ij,ab}{\text{stat}}
\equiv W_{ij,ab}(\omega = 0)
= \ERI{ij}{ab} - 4 \sum_m \frac{\sERI{ij}{m} \sERI{ab}{m}}{\OmRPA{m}{} - i \eta},
\end{equation}
evidencing that the standard static BSE problem is recovered from the present dynamical formalism in this limit.
Due to excitonic effects, the lowest BSE excitation energy, $\Om{1}{}$, stands lower than the lowest RPA excitation energy, $\Om{1}{\RPA}$, so that, $\Om{ib}{S} - \Om{m}{\RPA} < 0 $ and $\widetilde{W}_{ij,ab}(\Om{S}{})$ has no resonances.
Furthermore, $\Om{ib}{S}$ and $\Om{ja}{S}$ are necessarily negative quantities for in-gap low-lying BSE excitations.
Thus, we have $\abs*{\Om{ib}{S} - \Om{m}{\RPA}} > \Om{m}{\RPA}$.
As a consequence, we observe a reduction of the electron-hole screening, \ie, an enhancement of electron-hole binding energy, as compared to the standard static BSE, and consequently smaller (red-shifted) excitation energies.
This will be numerically illustrated in Sec.~\ref{sec:resdis}.
%================================
\subsection{Tamm-Dancoff (TDA) approximation }
\xavier{
The analysis of the (off-diagonal) screened Coulomb potential matrix elements multiplying the $Y_{jb,s}$ coefficients
\begin{multline}
\widetilde{W}_{ib,aj}( \Om{s}{} )
%================================================
\subsection{Dynamical Tamm-Dancoff approximation}
%================================================
The analysis of the (off-diagonal) screened Coulomb potential matrix elements multiplying the $Y_{jb,S}$ coefficients in Eq.~\eqref{eq:BSE-final}, \ie,
\begin{multline} \label{eq:W-Y}
\widetilde{W}_{ib,aj}(\Om{S}{})
= \ERI{ib}{aj} + 2 \sum_m \sERI{ib}{m} \sERI{aj}{m}
\\
\times \qty[ \frac{1}{\Om{ij}{s} - \Om{m}{\RPA} + i\eta} + \frac{1}{\Om{ba}{s} - \Om{m}{\RPA} + i\eta} ].
\times \qty[ \frac{1}{\Om{ij}{S} - \Om{m}{\RPA} + i\eta} + \frac{1}{\Om{ba}{S} - \Om{m}{\RPA} + i\eta} ],
\end{multline}
reveals on the contrary strong divergences even for low-lying excitations with, \eg,
$$
\Om{ba}{s} - \Om{m}{\RPA} = \Om{s}{} - \Om{m}{\RPA} - ( \eGW{a} - \eGW{b} )
$$
Since $( \eGW{a} - \eGW{b})$ can take small to large positive or negative values, $(a,b)$ indexes such that $(\Om{ba}{s} - \Om{m}{\RPA})$ cancels can always occur, even for low lying $\Om{s}{}$ excitations, namely negative $( \Om{s}{} - \Om{m}{\RPA} )$ energies. Such divergences may explain that in previous calculations dynamical effects were only accounted for at the TDA level.
Going beyond the TDA is beyond the present study.
}
reveals strong divergences even for low-lying excitations with, for example, $\Om{ba}{S} - \Om{m}{\RPA} = \Om{S}{} - \Om{m}{\RPA} - ( \eGW{a} - \eGW{b} ) \approx 0$.
Such divergences may explain that in previous studies dynamical effects were only accounted for at the TDA level. \cite{Strinati_1988,Rohlfing_2000,Ma_2009a,Ma_2009b,Romaniello_2009b,Sangalli_2011,Zhang_2013,Rebolini_2016}
To avoid confusions here, enforcing the TDA for the dynamical correction (which corresponds to neglecting the dynamical correction originating from the anti-resonant part of the BSE Hamiltonian) will be labeled as dTDA in the following.
Going beyond the dTDA is outside the scope of the present study but shall be addressed eventually.
%=================================
%================================
\subsection{Perturbative dynamical correction}
%=================================
@ -558,44 +562,44 @@ From a more practical point of view, Eq.~\eqref{eq:BSE-final} can be recast as a
\begin{equation}
\label{eq:LR-dyn}
\begin{pmatrix}
\bA{}(\Om{s}{}) & \bB{}(\Om{s}{})
\bA{}(\Om{S}{}) & \bB{}(\Om{S}{})
\\
-\bB{}(-\Om{s}{}) & -\bA{}(-\Om{s}{})
-\bB{}(-\Om{S}{}) & -\bA{}(-\Om{S}{})
\\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{s}{} \\
\bY{s}{} \\
\bX{S}{} \\
\bY{S}{} \\
\end{pmatrix}
=
\Om{s}{}
\Om{S}{}
\begin{pmatrix}
\bX{s}{} \\
\bY{s}{} \\
\bX{S}{} \\
\bY{S}{} \\
\end{pmatrix},
\end{equation}
where the dynamical matrices $\bA{}$ and $\bB{}$ have the same $\Nocc \Nvir \times \Nocc \Nvir$ size than their RPA counterparts.
Same comment applies to the eigenvectors $\bX{s}{}$, and $\bY{s}{}$ of length $\Nocc \Nvir$.
Note that, due to its non-linear nature, Eq.~\eqref{eq:LR-dyn} may provide more than one solution for each value of $s$. \cite{Romaniello_2009b,Sangalli_2011,Martin_2016}
where the dynamical matrices $\bA{}$ and $\bB{}$ have the same $\Nocc \Nvir \times \Nocc \Nvir$ size than their RPA counterparts, and we assume real quantities from hereon.
Same comment applies to the eigenvectors $\bX{S}{}$, and $\bY{S}{}$ of length $\Nocc \Nvir$.
Note that, due to its non-linear nature, Eq.~\eqref{eq:LR-dyn} may provide more than one solution for each value of $S$. \cite{Romaniello_2009b,Sangalli_2011,Martin_2016}
Accordingly to Eq.~\eqref{eq:BSE-final}, the BSE matrix elements in Eq.~\eqref{eq:LR-dyn} read
\begin{subequations}
\begin{align}
\label{eq:BSE-Adyn}
\A{ia,jb}{}(\Om{s}{}) & = \delta_{ij} \delta_{ab} (\eGW{a} - \eGW{i}) + \kappa \ERI{ia}{jb} - \tW{ij,ab}{}(\Om{s}{}),
\A{ia,jb}{}(\Om{S}{}) & = \delta_{ij} \delta_{ab} (\eGW{a} - \eGW{i}) + \kappa \ERI{ia}{jb} - \tW{ij,ab}{}(\Om{S}{}),
\\
\label{eq:BSE-Bdyn}
\B{ia,jb}{}(\Om{s}{}) & = \kappa \ERI{ia}{bj} - \tW{ib,aj}{}(\Om{s}{}).
\B{ia,jb}{}(\Om{S}{}) & = \kappa \ERI{ia}{bj} - \tW{ib,aj}{}(\Om{S}{}).
\end{align}
\end{subequations}
Now, let us decompose, using basic perturbation theory, the non-linear eigenproblem \eqref{eq:LR-dyn} as a zeroth-order static (hence linear) reference and a first-order dynamic (hence non-linear) perturbation, such that
Now, let us decompose, using basic Rayleigh-Schr\"odinger perturbation theory, the non-linear eigenproblem \eqref{eq:LR-dyn} as a zeroth-order static (hence linear) reference and a first-order dynamic (hence non-linear) perturbation, such that
\begin{multline}
\label{eq:LR-PT}
\begin{pmatrix}
\bA{}(\Om{s}{}) & \bB{}(\Om{s}{}) \\
-\bB{}(-\Om{s}{}) & -\bA{}(-\Om{s}{}) \\
\bA{}(\Om{S}{}) & \bB{}(\Om{S}{}) \\
-\bB{}(-\Om{S}{}) & -\bA{}(-\Om{S}{}) \\
\end{pmatrix}
\\
=
@ -607,8 +611,8 @@ Now, let us decompose, using basic perturbation theory, the non-linear eigenprob
\end{pmatrix}
+
\begin{pmatrix}
\bA{(1)}(\Om{s}{}) & \bB{(1)}(\Om{s}{}) \\
-\bB{(1)}(-\Om{s}{}) & -\bA{(1)}(-\Om{s}{}) \\
\bA{(1)}(\Om{S}{}) & \bB{(1)}(\Om{S}{}) \\
-\bB{(1)}(-\Om{S}{}) & -\bA{(1)}(-\Om{S}{}) \\
\end{pmatrix},
\end{multline}
with
@ -625,37 +629,31 @@ and
\begin{subequations}
\begin{align}
\label{eq:BSE-A1}
\A{ia,jb}{(1)}(\Om{s}{}) & = - \tW{ij,ab}{}(\Om{s}{}) + \W{ij,ab}{\text{stat}},
\A{ia,jb}{(1)}(\Om{S}{}) & = - \tW{ij,ab}{}(\Om{S}{}) + \W{ij,ab}{\text{stat}},
\\
\label{eq:BSE-B1}
\B{ia,jb}{(1)}(\Om{s}{}) & = - \tW{ib,aj}{}(\Om{s}{}) + \W{ib,aj}{\text{stat}},
\B{ia,jb}{(1)}(\Om{S}{}) & = - \tW{ib,aj}{}(\Om{S}{}) + \W{ib,aj}{\text{stat}}.
\end{align}
\end{subequations}
where we have defined the static version of the screened Coulomb potential [see Eq.~\eqref{eq:W-RPA}]
\begin{equation}
\label{eq:Wstat}
\W{ij,ab}{\text{stat}}
\equiv W_{ij,ab}(\omega = 0)
= \ERI{ij}{ab} - 4 \sum_m \frac{\sERI{ij}{m} \sERI{ab}{m}}{\OmRPA{m}{} - i \eta}.
\end{equation}
According to perturbation theory, the $s$th BSE excitation energy and its corresponding eigenvector can then expanded as
According to perturbation theory, the $S$th BSE excitation energy and its corresponding eigenvector can then expanded as
\begin{subequations}
\begin{gather}
\Om{s}{} = \Om{s}{(0)} + \Om{s}{(1)} + \ldots,
\Om{S}{} = \Om{S}{(0)} + \Om{S}{(1)} + \ldots,
\\
\begin{pmatrix}
\bX{s}{} \\
\bY{s}{} \\
\bX{S}{} \\
\bY{S}{} \\
\end{pmatrix}
=
\begin{pmatrix}
\bX{s}{(0)} \\
\bY{s}{(0)} \\
\bX{S}{(0)} \\
\bY{S}{(0)} \\
\end{pmatrix}
+
\begin{pmatrix}
\bX{s}{(1)} \\
\bY{s}{(1)} \\
\bX{S}{(1)} \\
\bY{S}{(1)} \\
\end{pmatrix}
+ \ldots.
\end{gather}
@ -669,51 +667,51 @@ Solving the zeroth-order static problem
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{s}{(0)} \\
\bY{s}{(0)} \\
\bX{S}{(0)} \\
\bY{S}{(0)} \\
\end{pmatrix}
=
\Om{s}{(0)}
\Om{S}{(0)}
\begin{pmatrix}
\bX{s}{(0)} \\
\bY{s}{(0)} \\
\bX{S}{(0)} \\
\bY{S}{(0)} \\
\end{pmatrix},
\end{equation}
yields the zeroth-order (static) $\Om{s}{(0)}$ excitation energies and their corresponding eigenvectors $\bX{s}{(0)}$ and $\bY{s}{(0)}$.
Thanks to first-order perturbation theory, the first-order correction to the $s$th excitation energy is
yields the zeroth-order (static) $\Om{S}{(0)}$ excitation energies and their corresponding eigenvectors $\bX{S}{(0)}$ and $\bY{S}{(0)}$.
Thanks to first-order perturbation theory, the first-order correction to the $S$th excitation energy is
\begin{equation}
\label{eq:Om1}
\Om{s}{(1)} =
\Om{S}{(1)} =
\T{\begin{pmatrix}
\bX{s}{(0)} \\
\bY{s}{(0)} \\
\bX{S}{(0)} \\
\bY{S}{(0)} \\
\end{pmatrix}}
\cdot
\begin{pmatrix}
\bA{(1)}(\Om{s}{(0)}) & \bB{(1)}(\Om{s}{(0)}) \\
-\bB{(1)}(-\Om{s}{(0)}) & -\bA{(1)}(-\Om{s}{(0)}) \\
\bA{(1)}(\Om{S}{(0)}) & \bB{(1)}(\Om{S}{(0)}) \\
-\bB{(1)}(-\Om{S}{(0)}) & -\bA{(1)}(-\Om{S}{(0)}) \\
\end{pmatrix}
\cdot
\begin{pmatrix}
\bX{s}{(0)} \\
\bY{s}{(0)} \\
\bX{S}{(0)} \\
\bY{S}{(0)} \\
\end{pmatrix}.
\end{equation}
From a practical point of view, if one enforces the TDA for the dynamical correction (which we label dTDA in the following), we obtain the very simple expression
From a practical point of view, if one enforces the dTDA, we obtain the very simple expression
\begin{equation}
\label{eq:Om1-TDA}
\Om{s}{(1)} = \T{(\bX{s}{(0)})} \cdot \bA{(1)}(\Om{s}{(0)}) \cdot \bX{s}{(0)}.
\Om{S}{(1)} = \T{(\bX{S}{(0)})} \cdot \bA{(1)}(\Om{S}{(0)}) \cdot \bX{S}{(0)}.
\end{equation}
This correction can be renormalized by computing, at basically no extra cost, the renormalization factor which reads, in the dTDA,
\begin{equation}
\label{eq:Z}
Z_{s} = \qty[ 1 - \T{(\bX{s}{(0)})} \cdot \left. \pdv{\bA{(1)}(\Om{s}{})}{\Om{s}{}} \right|_{\Om{s}{} = \Om{s}{(0)}} \cdot \bX{s}{(0)} ]^{-1}.
Z_{S} = \qty[ 1 - \T{(\bX{S}{(0)})} \cdot \left. \pdv{\bA{(1)}(\Om{S}{})}{\Om{S}{}} \right|_{\Om{S}{} = \Om{S}{(0)}} \cdot \bX{S}{(0)} ]^{-1}.
\end{equation}
This finally yields
\begin{equation}
\Om{s}{\text{dyn}} = \Om{s}{\text{stat}} + \Delta\Om{s}{\text{dyn}} = \Om{s}{(0)} + Z_{s} \Om{s}{(1)}.
\Om{S}{\text{dyn}} = \Om{S}{\text{stat}} + \Delta\Om{S}{\text{dyn}} = \Om{S}{(0)} + Z_{S} \Om{S}{(1)}.
\end{equation}
with $\Om{s}{\text{stat}} \equiv \Om{s}{(0)}$ and $\Delta\Om{s}{\text{dyn}} = Z_{s} \Om{s}{(1)}$.
with $\Om{S}{\text{stat}} \equiv \Om{S}{(0)}$ and $\Delta\Om{S}{\text{dyn}} = Z_{S} \Om{S}{(1)}$.
This is our final expression.
%%% FIG 1 %%%
@ -726,7 +724,7 @@ This is our final expression.
%%% %%% %%%
In terms of computational cost, if one decides to compute the dynamical correction of the $M$ lowest excitation energies, one must perform, first, a conventional (static) BSE calculation and extract the $M$ lowest eigenvalues and their corresponding eigenvectors [see Eq.~\eqref{eq:LR-BSE-stat}].
These are then used to compute the first-order correction from Eq.~\eqref{eq:Om1} or Eq.~\eqref{eq:Om1-TDA}, which also require to construct and evaluate the dynamical part of the BSE Hamiltonian for each excitation one wants to dynamically correct.
These are then used to compute the first-order correction from Eq.~\eqref{eq:Om1-TDA}, which also require to construct and evaluate the dynamical part of the BSE Hamiltonian for each excitation one wants to dynamically correct.
The static BSE Hamiltonian is computed once during the static BSE calculation and does not dependent on the targeted excitation.
Searching iteratively for the lowest eigenstates, via Davidson's algorithm for instance, can be performed in $\order*{\Norb^4}$ computational cost.
@ -737,7 +735,7 @@ Although it might be reduced to $\order*{\Norb^4}$ operations with standard reso
\section{Computational details}
\label{sec:compdet}
%%%%%%%%%%%%%%%%%%%%%%%%
All systems under investigation have close-shell singlet ground states.
All systems under investigation have closed-shell singlet ground states.
We then adopt a restricted formalism throughout this work.
The $GW$ calculations performed to obtain the screened Coulomb operator and the quasiparticle energies are done using a (restricted) HF starting point.
Perturbative $GW$ (or {\GOWO}) \cite{Hybertsen_1985a, Hybertsen_1986} quasiparticle energies are employed as starting points to compute the BSE neutral excitations.
@ -768,62 +766,62 @@ All the static and dynamic BSE calculations have been performed with the softwar
& \mc{3}{c}{cc-pVTZ ($\Eg^{\GW} = 20.21$ eV)}
& \mc{3}{c}{cc-pVQZ ($\Eg^{\GW} = 20.05$ eV)} \\
\cline{3-5} \cline{6-8} \cline{9-11}
State & Nature & \tabc{$\Om{s}{\stat}$} & \tabc{$\Delta\Om{s}{\dyn}$(dTDA)} & \tabc{$\Delta\Om{s}{\dyn}$}
& \tabc{$\Om{s}{\stat}$} & \tabc{$\Delta\Om{s}{\dyn}$(dTDA)} & \tabc{$\Delta\Om{s}{\dyn}$}
& \tabc{$\Om{s}{\stat}$} & \tabc{$\Delta\Om{s}{\dyn}$(dTDA)} & \tabc{$\Delta\Om{s}{\dyn}$} \\
State & Nature & \tabc{$\Om{S}{\stat}$} & \tabc{$\Om{S}{\dyn}$} & \tabc{$\Delta\Om{S}{\dyn}$}
& \tabc{$\Om{S}{\stat}$} & \tabc{$\Om{S}{\dyn}$} & \tabc{$\Delta\Om{S}{\dyn}$}
& \tabc{$\Om{S}{\stat}$} & \tabc{$\Om{S}{\dyn}$} & \tabc{$\Delta\Om{S}{\dyn}$} \\
\hline
$^1\Pi_g(n \ra \pis)$ & Val. & 9.90 & -0.32 & -0.33 & 9.92 & -0.40 & -0.42 & 10.01 & -0.42 & -0.42 \\
$^1\Sigma_u^-(\pi \ra \pis)$ & Val. & 9.70 & -0.33 & -0.31 & 9.61 & -0.42 & -0.36 & 9.69 & -0.44 & -0.40 \\
$^1\Delta_u(\pi \ra \pis)$ & Val. & 10.37 & -0.31 & -0.32 & 10.27 & -0.39 & -0.42 & 10.34 & -0.41 & -0.41 \\
$^1\Sigma_g^+$(R) & Ryd. & 15.67 & -0.17 & -0.59 & 15.04 & -0.21 & -0.46 & 14.72 & -0.21 & -0.33 \\
$^1\Pi_u$(R) & Ryd. & 15.00 & -0.21 & -0.26 & 14.75 & -0.27 & -0.30 & 14.72 & -0.29 & -0.30 \\
$^1\Sigma_u^+$(R) & Ryd. & 22.88\fnm[1] & -0.15 & -0.20 & 19.03 & -0.08 & -0.07 & 16.78 & -0.06 & -0.07 \\
$^1\Pi_u$(R) & Ryd. & 23.62\fnm[1] & -0.11 & -0.10 & 19.15 & -0.11 & -0.10 & 16.93 & -0.09 & -0.08 \\
$^1\Pi_g(n \ra \pis)$ & Val. & 9.90 & 9.58 & -0.32 & 9.92 & 9.53 & -0.40 & 10.01 & 9.59 & -0.42 \\
$^1\Sigma_u^-(\pi \ra \pis)$ & Val. & 9.70 & 9.37 & -0.33 & 9.61 & 9.19 & -0.42 & 9.69 & 9.25 & -0.44 \\
$^1\Delta_u(\pi \ra \pis)$ & Val. & 10.37 & 10.05 & -0.31 & 10.27 & 9.88 & -0.39 & 10.34 & 9.93 & -0.41 \\
$^1\Sigma_g^+$ & Ryd. & 15.67 & 15.50 & -0.17 & 15.04 & 14.84 & -0.21 & 14.72 & 14.43 & -0.21 \\
$^1\Pi_u$ & Ryd. & 15.00 & 14.79 & -0.21 & 14.75 & 14.48 & -0.27 & 14.80 & 14.59 & -0.29 \\
$^1\Sigma_u^+$ & Ryd. & 22.88\fnm[1] & 22.73 & -0.15 & 19.03 & 18.95 & -0.08 & 16.78 & 16.71 & -0.06 \\
$^1\Pi_u$ & Ryd. & 23.62\fnm[1] & 23.51 & -0.11 & 19.15 & 19.04 & -0.11 & 16.93 & 16.85 & -0.09 \\
\\
$^3\Sigma_u^+(\pi \ra \pis)$ & Val. & 7.39 & -0.48 & -0.69 & 7.46 & -0.59 & -0.41 & 7.59 & -0.62 & -0.50 \\
$^3\Pi_g(n \ra \pis)$ & Val. & 8.07 & -0.42 & -0.43 & 8.14 & -0.52 & -0.48 & 8.24 & -0.54 & -0.48 \\
$^3\Delta_u(\pi \ra \pis)$ & Val. & 8.56 & -0.41 & -0.40 & 8.52 & -0.52 & -0.40 & 8.62 & -0.55 & -0.41 \\
$^3\Sigma_u^-(\pi \ra \pis)$ & Val. & 9.70 & -0.33 & -0.31 & 9.61 & -0.42 & -0.36 & 9.69 & -0.44 & -0.40 \\
$^3\Sigma_u^+(\pi \ra \pis)$ & Val. & 7.39 & 6.91 & -0.48 & 7.46 & 6.87 & -0.59 & 7.59 & 6.97 & -0.62 \\
$^3\Pi_g(n \ra \pis)$ & Val. & 8.07 & 7.65 & -0.42 & 8.14 & 7.62 & -0.52 & 8.24 & 7.70 & -0.54 \\
$^3\Delta_u(\pi \ra \pis)$ & Val. & 8.56 & 8.15 & -0.41 & 8.52 & 8.00 & -0.52 & 8.62 & 8.07 & -0.55 \\
$^3\Sigma_u^-(\pi \ra \pis)$ & Val. & 9.70 & 9.37 & -0.33 & 9.61 & 9.19 & -0.42 & 9.69 & 9.25 & -0.44 \\
\hline
\\
& & \mc{3}{c}{aug-cc-pVDZ ($\Eg^{\GW} = 19.49$ eV)}
& \mc{3}{c}{aug-cc-pVTZ ($\Eg^{\GW} = 19.20$ eV)}
& \mc{3}{c}{aug-cc-pVQZ ($\Eg^{\GW} = 19.00$ eV)} \\
\cline{3-5} \cline{6-8} \cline{9-11}
State & Nature & \tabc{$\Om{s}{\stat}$} & \tabc{$\Delta\Om{s}{\dyn}$(dTDA)} & \tabc{$\Delta\Om{s}{\dyn}$}
& \tabc{$\Om{s}{\stat}$} & \tabc{$\Delta\Om{s}{\dyn}$(dTDA)} & \tabc{$\Delta\Om{s}{\dyn}$}
& \tabc{$\Om{s}{\stat}$} & \tabc{$\Delta\Om{s}{\dyn}$(dTDA)} & \tabc{$\Delta\Om{s}{\dyn}$} \\
State & Nature & \tabc{$\Om{S}{\stat}$} & \tabc{$\Om{S}{\dyn}$} & \tabc{$\Delta\Om{S}{\dyn}$}
& \tabc{$\Om{S}{\stat}$} & \tabc{$\Om{S}{\dyn}$} & \tabc{$\Delta\Om{S}{\dyn}$}
& \tabc{$\Om{S}{\stat}$} & \tabc{$\Om{S}{\dyn}$} & \tabc{$\Delta\Om{S}{\dyn}$} \\
\hline
$^1\Pi_g(n \ra \pis)$ & Val. & 10.18 & -0.41 & -0.43 & 10.42 & -0.42 & -0.41 & 10.52 & -0.43 & -0.29 \\
$^1\Sigma_u^-(\pi \ra \pis)$ & Val. & 9.95 & -0.44 & -0.41 & 10.11 & -0.45 & -0.42 & 10.20 & -0.45 & -0.41 \\
$^1\Delta_u(\pi \ra \pis)$ & Val. & 10.57 & -0.41 & -0.41 & 10.75 & -0.42 & -0.45 & 10.85 & -0.42 & -0.48 \\
$^1\Sigma_g^+$ & Ryd. & 13.72 & -0.04 & -0.05 & 13.60 & -0.03 & -0.03 & 13.55 & -0.02 & -0.03 \\
$^1\Pi_u$ & Ryd. & 14.07 & -0.05 & -0.05 & 13.98 & -0.04 & -0.04 & 13.96 & -0.03 & -0.03 \\
$^1\Sigma_u^+$ & Ryd. & 13.80 & -0.08 & -0.10 & 13.98 & -0.07 & -0.10 & 14.08 & -0.06 & -0.08 \\
$^1\Pi_u$ & Ryd. & 14.22 & -0.04 & -0.03 & 14.24 & -0.03 & -0.03 & 14.26 & -0.03 & -0.02 \\
$^1\Pi_g(n \ra \pis)$ & Val. & 10.18 & 9.77 & -0.41 & 10.42 & 9.99 & -0.42 & 10.52 & 10.09 & -0.43 \\
$^1\Sigma_u^-(\pi \ra \pis)$ & Val. & 9.95 & 9.51 & -0.44 & 10.11 & 9.66 & -0.45 & 10.20 & 9.75 & -0.45 \\
$^1\Delta_u(\pi \ra \pis)$ & Val. & 10.57 & 10.16 & -0.41 & 10.75 & 10.33 & -0.42 & 10.85 & 10.42 & -0.42 \\
$^1\Sigma_g^+$ & Ryd. & 13.72 & 13.68 & -0.04 & 13.60 & 13.57 & -0.03 & 13.54 & 13.52 & -0.02 \\
$^1\Pi_u$ & Ryd. & 14.07 & 14.02 & -0.05 & 13.98 & 13.94 & -0.04 & 13.96 & 13.93 & -0.03 \\
$^1\Sigma_u^+$ & Ryd. & 13.80 & 13.72 & -0.08 & 13.98 & 13.91 & -0.07 & 14.08 & 14.03 & -0.06 \\
$^1\Pi_u$ & Ryd. & 14.22 & 14.19 & -0.04 & 14.24 & 14.21 & -0.03 & 14.26 & 14.23 & -0.03 \\
\\
$^3\Sigma_u^+(\pi \ra \pis)$ & Val. & 7.75 & -0.63 & -2.42 & 8.02 & -0.64 & -0.45 & 8.12 & -0.64 & -0.53 \\
$^3\Pi_g(n \ra \pis)$ & Val. & 8.42 & -0.54 & -0.50 & 8.66 & -0.56 & -0.79 & 8.75 & -0.56 & -0.47 \\
$^3\Delta_u(\pi \ra \pis)$ & Val. & 8.86 & -0.54 & -0.47 & 9.04 & -0.56 & -0.52 & 9.14 & -0.56 & -0.51 \\
$^3\Sigma_u^-(\pi \ra \pis)$ & Val. & 9.95 & -0.44 & -0.41 & 10.11 & -0.45 & -0.42 & 10.20 & -0.45 & -0.41 \\
$^3\Sigma_u^+(\pi \ra \pis)$ & Val. & 7.75 & 7.12 & -0.63 & 8.02 & 7.38 & -0.64 & 8.12 & 7.48 & -0.64 \\
$^3\Pi_g(n \ra \pis)$ & Val. & 8.42 & 7.88 & -0.54 & 8.66 & 8.10 & -0.56 & 8.75 & 8.20 & -0.56 \\
$^3\Delta_u(\pi \ra \pis)$ & Val. & 8.86 & 8.32 & -0.54 & 9.04 & 8.48 & -0.56 & 9.14 & 8.57 & -0.56 \\
$^3\Sigma_u^-(\pi \ra \pis)$ & Val. & 9.95 & 9.51 & -0.44 & 10.11 & 9.66 & -0.45 & 10.20 & 9.75 & -0.45 \\
\end{tabular}
\end{ruledtabular}
\fnt[1]{Excitation energy larger than the fundamental gap.}
\end{table*}
\end{squeezetable}
First, we investigate the basis set dependency of the dynamical correction as well as the validity of the dTDA (which corresponds to neglecting the dynamical correction originating from the anti-resonant part of the BSE Hamiltonian).
First, we investigate the basis set dependency of the dynamical correction.
Note that, in the present calculations, the zeroth-order Hamiltonian is always the ``full'' BSE static Hamiltonian, \ie, without TDA.
The singlet and triplet excitation energies of the nitrogen molecule \ce{N2} computed at the BSE@{\GOWO}@HF level for the aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets are reported in Table \ref{tab:N2}, where we also report the $GW$ gap, $\Eg^{\GW}$, to show that each corrected transition is well below this gap.
The singlet and triplet excitation energies of the nitrogen molecule \ce{N2} computed at the BSE@{\GOWO}@HF level for the aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ basis sets are reported in Table \ref{tab:N2}, where we also report the $GW$ gap, $\Eg^{\GW}$, to show that corrected transitions are usually well below this gap.
The \ce{N2} molecule is a very convenient example as it contains $n \ra \pis$ and $\pi \ra \pis$ valence excitations as well as Rydberg transitions.
As we shall further illustrate below, the magnitude of the dynamical correction is characteristic of the type of transitions.
One key result of the present investigation is that the dynamical correction is quite basis set insensitive with a maximum variation of $0.03$ eV between in smallest (aug-cc-pVDZ) and largest (aug-cc-pVQZ) basis sets.
This is quite a nice feature as one does not need to compute this more expensive correction in a very large basis.
The second important observation extracted from the results gathered in Table \ref{tab:N2} is that the dTDA is a rather satisfactory approximation, especially for the singlet states where one observes a maximum discrepancy of $0.03$ eV between the ``full'' and dTDA excitation energies.
The story is different for the triplet states for which deviations of the order of $0.3$ eV is observed between the two sets, the dTDA of excitation energies being, as we shall see later on, more accurate.
Indeed, although the dynamical correction systematically red-shift the excitation energies (as anticipated in Sec.~\ref{sec:dynW}), taking into account the coupling between the resonant and anti-resonant parts of the BSE Hamiltonian [see Eq.~\eqref{eq:BSE-final}] yields a systematic blue-shift of the correction, the overall correction remaining negative but by a smaller amount.
This outcome is similar to the conclusions of several benchmark studies \cite{Jacquemin_2017b,Rangel_2017} which clearly concluded that static BSE triplet excitations are notably too low in energy and that the use of the TDA is able to partly reduce this error.
In accordance with the success of the dTDA, the remaining calculations of the present study are performed within this approximation.
%The second important observation extracted from the results gathered in Table \ref{tab:N2} is that the dTDA is a rather satisfactory approximation, especially for the singlet states where one observes a maximum discrepancy of $0.03$ eV between the ``full'' and dTDA excitation energies.
%The story is different for the triplet states for which deviations of the order of $0.3$ eV is observed between the two sets, the dTDA of excitation energies being, as we shall see later on, more accurate.
%Indeed, although the dynamical correction systematically red-shift the excitation energies (as anticipated in Sec.~\ref{sec:dynW}), taking into account the coupling between the resonant and anti-resonant parts of the BSE Hamiltonian [see Eq.~\eqref{eq:BSE-final}] yields a systematic blue-shift of the correction, the overall correction remaining negative but by a smaller amount.
%This outcome is similar to the conclusions of several benchmark studies \cite{Jacquemin_2017b,Rangel_2017} which clearly concluded that static BSE triplet excitations are notably too low in energy and that the use of the TDA is able to partly reduce this error.
%In accordance with the success of the dTDA, the remaining calculations of the present study are performed within this approximation.
%%% TABLE I %%%
\begin{squeezetable}
@ -838,7 +836,7 @@ In accordance with the success of the dTDA, the remaining calculations of the pr
\begin{tabular}{llcdddddddddd}
& & & \mc{5}{c}{BSE@{\GOWO}@HF} & \mc{5}{c}{Wave function-based methods} \\
\cline{4-8} \cline{9-13}
Mol. & State & Nature & \tabc{$\Eg^{\GW}$} & \tabc{$\Om{s}{\stat}$} & \tabc{$\Om{s}{\dyn}$} & \tabc{$\Delta\Om{s}{\dyn}$} & \tabc{$Z_{s}$}
Mol. & State & Nature & \tabc{$\Eg^{\GW}$} & \tabc{$\Om{S}{\stat}$} & \tabc{$\Om{S}{\dyn}$} & \tabc{$\Delta\Om{S}{\dyn}$} & \tabc{$Z_{S}$}
& \tabc{CIS(D)} & \tabc{ADC(2)} & \tabc{CCSD} & \tabc{CC2} & \tabc{TBE} \\
\hline
\ce{HCl} & $^1\Pi$ & CT & 13.43 & 8.30 & 8.19 & -0.11 & 1.009
@ -932,7 +930,7 @@ In accordance with the success of the dTDA, the remaining calculations of the pr
\begin{tabular}{llcdddddddddd}
& & & \mc{5}{c}{BSE@{\GOWO}@HF} & \mc{5}{c}{Wave function-based methods} \\
\cline{4-8} \cline{9-13}
Mol. & State & Nature & \tabc{$\Eg^{\GW}$} & \tabc{$\Om{s}{\stat}$} & \tabc{$\Om{s}{\dyn}$} & \tabc{$\Delta\Om{s}{\dyn}$} & \tabc{$Z_{s}$}
Mol. & State & Nature & \tabc{$\Eg^{\GW}$} & \tabc{$\Om{S}{\stat}$} & \tabc{$\Om{S}{\dyn}$} & \tabc{$\Delta\Om{S}{\dyn}$} & \tabc{$Z_{S}$}
& \tabc{CIS(D)} & \tabc{ADC(2)} & \tabc{CCSD} & \tabc{CC2} & \tabc{TBE} \\
\hline
\ce{H2O} & $^3B_1(n \ra 3s)$ & Ryd. & 13.58 & 7.62 & 7.48 & -0.14 & 1.009
@ -1000,8 +998,8 @@ In accordance with the success of the dTDA, the remaining calculations of the pr
Tables \ref{tab:BigTabSi} and \ref{tab:BigTabTr} report, respectively, singlet and triplet excitation energies for various molecules computed at the BSE@{\GOWO}@HF level and with the aug-cc-pVTZ basis set.
For comparative purposes, excitation energies obtained with the same basis set and several second-order wave function methods [CIS(D), ADC(2), CCSD, and CC2] are also reported.
Finally, the highly-accurate TBEs of Refs.~\onlinecite{Loos_2018a,Loos_2019,Loos_2020b} will serve us as reference.
For each excitation, we report the static and dynamic excitation energies, $\Om{s}{\stat}$ and $\Om{s}{\dyn}$, as well as the value of the renormalization factor $Z_s$ defined in Eq.~\eqref{eq:Z}.
As one can see in Tables \ref{tab:BigTabSi} and \ref{tab:BigTabTr}, the value of $Z_s$ is always quite close to unity which shows that the perturbative expansion behaves nicely, and that a first-order correction is probably quite a good estimate of the non-perturbative result.
For each excitation, we report the static and dynamic excitation energies, $\Om{S}{\stat}$ and $\Om{S}{\dyn}$, as well as the value of the renormalization factor $Z_S$ defined in Eq.~\eqref{eq:Z}.
As one can see in Tables \ref{tab:BigTabSi} and \ref{tab:BigTabTr}, the value of $Z_S$ is always quite close to unity which shows that the perturbative expansion behaves nicely, and that a first-order correction is probably quite a good estimate of the non-perturbative result.
Moreover, we have observed that an iterative, self-consistent resolution [where the dynamically-corrected excitation energies are re-injected in Eq.~\eqref{eq:Om1}] yields basically the same results as its (cheaper) renormalized version.
%%% TABLE I %%%
@ -1016,7 +1014,7 @@ Moreover, we have observed that an iterative, self-consistent resolution [where
\begin{tabular}{llcdddddd}
& & & \mc{5}{c}{BSE@{\GOWO}@HF} \\
\cline{4-8}
Molecule & State & Nature & \tabc{$\Eg^{\GW}$} & \tabc{$\Om{s}{\stat}$} & \tabc{$\Om{s}{\dyn}$} & \tabc{$\Delta\Om{s}{\dyn}$} & \tabc{$Z_{s}$} & \tabc{CC3} \\
Molecule & State & Nature & \tabc{$\Eg^{\GW}$} & \tabc{$\Om{S}{\stat}$} & \tabc{$\Om{S}{\dyn}$} & \tabc{$\Delta\Om{S}{\dyn}$} & \tabc{$Z_{S}$} & \tabc{CC3} \\
\hline
acrolein & $^1A''(n \ra \pis)$ & Val. & 11.67 & 4.62 & 4.28 & -0.35 & 1.030 & 3.77 \\
& $^1A'(n \ra \pis)$ & Val. & & 6.86 & 6.70 & -0.16 & 1.023 & 6.67 \\
@ -1063,9 +1061,9 @@ Moreover, we have observed that an iterative, self-consistent resolution [where
This is the conclusion
%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Supplementary material}
%\section*{Supplementary material}
%%%%%%%%%%%%%%%%%%%%%%%%
See the {\SI} for a detailed derivation of the dynamical BSE equation.
%See the {\SI} for a detailed derivation of the dynamical BSE equation.
%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgements{
@ -1073,7 +1071,7 @@ See the {\SI} for a detailed derivation of the dynamical BSE equation.
PFL thanks the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481) for financial support.
This work was performed using HPC resources from GENCI-TGCC (Grant No.~2019-A0060801738) and CALMIP (Toulouse) under allocation 2020-18005.
Funding from the \textit{``Centre National de la Recherche Scientifique''} is acknowledged.
This work has also been supported through the EUR grant NanoX ANR-17-EURE-0009 in the framework of the \textit{``Programme des Investissements d'Avenir''.}}
This study has been (partially) supported through the EUR grant NanoX No.~ANR-17-EURE-0009 in the framework of the \textit{``Programme des Investissements d'Avenir''.}
%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Data availability}
@ -1083,7 +1081,7 @@ The data that support the findings of this study are available within the articl
\appendix
\section{$L_0(1,3; 1',4)$ $(t_1)$-time Fourier transform}
\label{appendixA}
\label{app:A}
In this Appendix, we derive Eqs.~\eqref{eq:iL0} to \eqref{eq:iL0bis}.
Defining the $t_1$-time Fourier transform of $L_0(1,3;4,1')$ with
@ -1129,30 +1127,31 @@ With the change of variable $\omega \to \omega + {\omega_1}/2$ one obtains readi
\end{equation}
where (pp) and (hh) labels particle-particle and hole-hole channels neglected here.
Projecting onto $\phi_a^*(\bx_1) \phi_i(\bx_{1'})$ selects the first line of the RHS, leading to Eq.~\eqref{eq:iL0bis}
with $ (\omega_1 \to \Omega_s )$.
with $ (\Om{1}{} \to \Om{s}{} )$.
\section{ $\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}$ in the electron/hole product basis }
\section{ $\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,S}$ in the electron/hole product basis }
\label{app:B}
We now derive in some more details Eq.~\eqref{eq:spectral65}.
Starting with
\begin{equation}
\begin{split}
\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,s}
& = \theta(+\tau_{65}) \mel{N}{ \hpsi(6) \hpsi^{\dagger}(5) }{N,s}
\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)] }{N,S}
& = \theta(+\tau_{65}) \mel{N}{ \hpsi(6) \hpsi^{\dagger}(5) }{N,S}
\\
& - \theta(-\tau_{65}) \mel{N}{ \hpsi^{\dagger}(5) \hpsi(6) }{N,s}
& - \theta(-\tau_{65}) \mel{N}{ \hpsi^{\dagger}(5) \hpsi(6) }{N,S}
\end{split}
\end{equation}
we use the relation between operators in their Heisenberg and Schr\"{o}dinger representations [see Eq.~\eqref{Eisenberg}] to obtain
\begin{equation}
\begin{split}
& \mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)]}{N,s} = \\
& + \theta(\tau_{65}) \mel{N}{ \hpsi(\bx_6) e^{-i\hH \tau_{65}} \hpsi^{\dagger}(\bx_5) }{N,s} e^{ i E^N_0 t_6 } e^{ - i E^N_s t_5 }
& \mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)]}{N,S} = \\
& + \theta(\tau_{65}) \mel{N}{ \hpsi(\bx_6) e^{-i\hH \tau_{65}} \hpsi^{\dagger}(\bx_5) }{N,S} e^{ i E^N_0 t_6 } e^{ - i E^N_S t_5 }
\\
& - \theta(-\tau_{65}) \mel{N}{ \hpsi^{\dagger}(\bx_5) e^{ i\hH \tau_{65}} \hpsi(\bx_6) }{N,s} e^{ i E^N_0 t_5 } e^{ - i E^N_s t_6 }
& - \theta(-\tau_{65}) \mel{N}{ \hpsi^{\dagger}(\bx_5) e^{ i\hH \tau_{65}} \hpsi(\bx_6) }{N,S} e^{ i E^N_0 t_5 } e^{ - i E^N_S t_6 }
\end{split}
\end{equation}
with $E^N_0$ the $N$-electron ground-state energy and $E^N_s$ the energy of the $s$th excited state $\ket{N,s}$.
with $E^N_0$ the $N$-electron ground-state energy and $E^N_S$ the energy of the $S$th excited state $\ket{N,S}$.
Expanding now the field operators with creation/destruction operators in the orbital basis
\begin{subequations}
\begin{align}
@ -1164,12 +1163,12 @@ Expanding now the field operators with creation/destruction operators in the orb
one gets
\begin{equation}
\begin{split}
\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)]}{N,s}
\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)]}{N,S}
\\
= \sum_{pq} \phi_p(\bx_6) \phi_q^{*}(\bx_5)
[ & \theta(+\tau_{65}) \mel{N}{ \ha_p e^{-i \hH \tau_{65}} \ha^{\dagger}_q }{N,s} e^{ i E^N_0 t_6 } e^{ - i E^N_s t_5 }
[ & \theta(+\tau_{65}) \mel{N}{ \ha_p e^{-i \hH \tau_{65}} \ha^{\dagger}_q }{N,S} e^{ i E^N_0 t_6 } e^{ - i E^N_S t_5 }
\\
- & \theta(-\tau_{65}) \mel{N}{ \ha^{\dagger}_q e^{ i \hH \tau_{65}} \ha_p }{N,s} e^{ i E^N_0 t_5 } e^{ - i E^N_s t_6 } ]
- & \theta(-\tau_{65}) \mel{N}{ \ha^{\dagger}_q e^{ i \hH \tau_{65}} \ha_p }{N,S} e^{ i E^N_0 t_5 } e^{ - i E^N_S t_6 } ]
\end{split}
\end{equation}
We now act on the $N$-electron ground-state with
@ -1183,18 +1182,18 @@ We now act on the $N$-electron ground-state with
\end{align}
\end{subequations}
where $\lbrace \e{p/q} \rbrace$ are quasiparticle energies, such as the $GW$ ones, namely proper addition/removal energies.
Taking the associated bras that we plug into the orbital product basis expansion of $\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)]}{N,s}$ one obtains:
Taking the associated bras that we plug into the orbital product basis expansion of $\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)]}{N,S}$ one obtains:
\begin{equation}
\begin{split}
\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)]}{N,s}
\mel{N}{T [\hpsi(6) \hpsi^{\dagger}(5)]}{N,S}
\\
= \sum_{pq} \phi_p(\bx_6) \phi_q^{*}(\bx_5)
[ & \theta(+ \tau_{65}) \mel{N}{ \ha_p \ha^{\dagger}_q }{N,s} e^{ -i \e{p} \tau_{65} } e^{ - i \Omega_s t_5 }
[ & \theta(+ \tau_{65}) \mel{N}{ \ha_p \ha^{\dagger}_q }{N,S} e^{ -i \e{p} \tau_{65} } e^{ - i \Om{S}{} t_5 }
\\
- & \theta(-\tau_{65}) \mel{N}{ \ha^{\dagger}_q \ha_p }{N,s} e^{ -i \e{q} \tau_{65} } e^{ - i \Omega_s t_6 } ]
- & \theta(-\tau_{65}) \mel{N}{ \ha^{\dagger}_q \ha_p }{N,S} e^{ -i \e{q} \tau_{65} } e^{ - i \Om{S}{} t_6 } ]
\end{split}
\end{equation}
leading to Eq.~\eqref{eq:spectral65} with $\Omega_s = (E^N_s - E^N_0)$, $t_6 = \tau_{65}/2 + t^{65}$ and $t_5 = - \tau_{65}/2 + t^{65}$.
leading to Eq.~\eqref{eq:spectral65} with $\Om{S}{} = (E^N_S - E^N_0)$, $t_6 = \tau_{65}/2 + t^{65}$ and $t_5 = - \tau_{65}/2 + t^{65}$.
\bibliography{BSEdyn}