Let us consider an electronic system consisting of $n = n_\up+ n_\dw$ electrons (where $n_\up$ and $n_\dw$ are the number of spin-up and spin-down electrons, respectively) and $N$ one-electron basis functions.
The number of spin-up and spin-down occupied orbitals are $O_\up= n_\up$ and $O_\dw= n_\dw$, respectively, and, assuming no linear dependencies in the one-electron basis set, there is $V_\up= N - O_\up$ and $V_\dw= N - O_\dw$ spin-up and spin-down virtual (\ie, unoccupied) orbitals.
The number of spin-conserved single excitations is then $S^\spc= S_{\up\up}^\spc+ S_{\dw\dw}^\spc= O_\up V_\up+ O_\dw V_\dw$, while the number of spin-flip excitations is $S^\spf= S_{\up\dw}^\spf+ S_{\dw\up}^\spf= O_\up V_\dw+ O_\dw V_\up$.
It is important to understand that, in a spin-conserved excitation the hole orbital $\MO{i_\sig}$ and particle orbital $\MO{a_\sig}$ have the same spin $\sig$.
In the following, we assume real quantities throughout this manuscript, $i$ and $j$ are occupied orbitals, $a$ and $b$ are unoccupied orbitals, $p$, $q$, $r$, and $s$ indicate arbitrary orbitals, and $m$ labels single excitations.
Moreover, we consider systems with collinear spins and a spin-independent hamiltonian without contributions such as spin-orbit interaction.
The pillar of Green's function many-body perturbation theory is the (time-ordered) one-body Green's function, which has poles at the charged excitations (i.e., ionization potentials and electron affinities) of the system. \cite{ReiningBook}
The spin-$\sig$ component of the one-body Green's function reads \cite{ReiningBook,Bruneval_2016a}
Within the $GW$ formalism, the dynamical screening is computed at the random-phase approximation (RPA) level by considering only the manifold of the spin-conserved neutral excitations.
In the orbital basis, the spectral representation of $W$ is
In Eqs.~\eqref{eq:W_spectral} and \eqref{eq:sERI}, the RPA spin-conserved neutral excitations $\Om{m}{\spc,\RPA}$ and their corresponding eigenvectors, $\bX{m}{\spc,\RPA}$ and $\bY{m}{\spc,\RPA}$, are obtained by solving a linear response system of the form
In the absence of instabilities, the linear eigenvalue problem \eqref{eq:LR-RPA} has particle-hole symmetry which means that the eigenvalues are obtained by pairs $\pm\Om{m}{}$.
In such a case, $(\bA{}{}-\bB{}{})^{1/2}$ is positive definite, and Eq.~\eqref{eq:LR-RPA} can be recast as a Hermitian problem of half the dimension
Within the Tamm-Dancoff approximation (TDA), the coupling terms between the resonant and anti-resonant parts, $\bA{}{}$ and $-\bA{}{}$, are neglected, which consists in setting $\bB{}{}=\bO$.
In such a case, Eq.~\eqref{eq:LR-RPA} reduces to straightforward Hermitian problem of the form:
where $G_{\KS}^{\sig}$ is the Kohn-Sham Green's function built with Kohn-Sham orbitals $\MO{p_\sig}^{\KS}(\br)$ and one-electron energies $\e{p_\sig}^{\KS}$ according to Eq.~\eqref{eq:G} and $v^{\xc}(\br)$ is the Kohn-Sham local exchange-correlation potential.
The target quantities here are the quasiparticle energies $\eGW{p_\sig}$, \ie, the poles of $G$ [see Eq.~\eqref{eq:G}], which correspond to well-defined addition/removal energies (unlike the KS orbital energies).
Because the exchange-correlation part of the self-energy is, itself, constructed with the Green's function [see Eq.~\eqref{eq:Sig}], the present process is, by nature, self-consistent.
The same comment applies to the dynamically-screened Coulomb potential $W$ entering the definition of $\Sig{}{\xc}$ [see Eq.~\eqref{eq:Sig}] which is also constructed from $G$ [see Eqs.~\eqref{eq:chi0}, \eqref{eq:eps}, and \eqref{eq:W}].
In its simplest perturbative (\ie, one-shot) version, known as {\GOWO}, a single iteration is performed, and the quasiparticle energies $\eGOWO{p_\sig}$ are obtained by solving the frequency-dependent quasiparticle equation
Because, from a practical point of view, one is usually interested by the so-called quasiparticle solution (or peak), the quasiparticle equation \eqref{eq:QP-eq} is often linearized around $\omega=\e{p_\sig}^{\KS}$, yielding
is a renormalization factor which also represents the spectral weight of the quasiparticle solution.
In addition to the principal quasiparticle peak which, in a well-behaved case, contains most of the spectral weight, the frequency-dependent quasiparticle equation \eqref{eq:QP-eq} generates a finite number of satellite resonances with smaller weights.
Within the ``eigenvalue'' self-consistent $GW$ scheme (known as ev$GW$), several iterations are performed during which only the one-electron energies entering the definition of the Green's function [see Eq.~\eqref{eq:G}] are updated by the quasiparticle energies obtained at the previous iteration (the corresponding orbitals remain evaluated at the Kohn-Sham level).
Finally, within the quasiparticle self-consistent $GW$ scheme (qs$GW$), both the one-electron energies and the orbitals are updated until convergence is reached.
These are obtained via the diagonalization of an effective Fock matrix which includes explicitly a frequency-independent and hermitian self-energy defined as
Like its TD-DFT cousin, BSE deals with the calculation of (neutral) optical excitations as measured by absorption spectroscopy. \cite{Salpeter_1951,Strinati_1988}
Using the BSE formalism, one can access the spin-conserved and spin-flip excitations.
The Dyson equation that links the generalized four-point susceptibility $L^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega)$ and the BSE kernel $\Xi^{\sig\sigp}(\br_3,\br_5;\br_4,\br_6;\omega)$ is
where, as usual, we have not considered the higher-order terms in $W$ by neglecting the derivative $\partial W/\partial G$. \cite{Hanke_1980, Strinati_1982, Strinati_1984, Strinati_1988}
Within the static approximation which consists in neglecting the frequency dependence of the dynamically-screened Coulomb potential, the spin-conserved and spin-flip optical excitation at the BSE level are obtained by solving a similar linear response problem
The dynamical correction to the static BSE kernel is defined in the Tamm-Dancoff approximation as \cite{Strinati_1988,Rohlfing_2000,Ma_2009a,Ma_2009b,Romaniello_2009b,Sangalli_2011,Loos_2020e}
The explicit expressions of $\Delta\expval{S^2}_m^{\spc}$ and $\Delta\expval{S^2}_m^{\spf}$ can be found in the Appendix of Ref.~\onlinecite{Li_2010} for spin-conserved and spin-flip excitations, and are functions of the $\bX{m}{}$ and $\bY{m}{}$ vectors and the orbital overlaps.
As explained in Ref.~\onlinecite{Casanova_2020}, there are two sources of spin contamination: i) spin contamination of the reference, and ii) spin-contamination of the excited states due to the spin incompleteness.
For the systems under investigation here, we consider either an open-shell doublet or triplet reference state.
We then adopt the unrestricted formalism throughout this work.
The $GW$ calculations performed to obtain the screened Coulomb operator and the quasiparticle energies are done using a (unrestricted) UHF starting point.
Perturbative $GW$ (or {\GOWO}) \cite{Hybertsen_1985a,Hybertsen_1986,vanSetten_2013} quasiparticle energies are employed as starting points to compute the BSE neutral excitations.
These quasiparticle energies are obtained by linearizing the frequency-dependent quasiparticle equation, and the entire set of orbitals is corrected.
Further details about our implementation of {\GOWO} can be found in Refs.~\onlinecite{Loos_2018b,Veril_2018,Loos_2020,Loos_2020e}.
%Note that, for the present (small) molecular systems, {\GOWO}@UHF and ev$GW$@UHF yield similar quasiparticle energies and fundamental gap.
%Moreover, {\GOWO} allows to avoid rather laborious iterations as well as the significant additional computational effort of ev$GW$.
%In the present study, the zeroth-order Hamiltonian [see Eq.~\eqref{eq:LR-PT}] is always the ``full'' BSE static Hamiltonian, \ie, without TDA.
The dynamical correction is computed in the TDA throughout.
As one-electron basis sets, we employ the Dunning families cc-pVXZ and aug-cc-pVXZ (X = D, T, and Q) defined with cartesian Gaussian functions.
Finally, the infinitesimal $\eta$ is set to $100$ meV for all calculations.
%It is important to mention that the small molecular systems considered here are particularly challenging for the BSE formalism, \cite{Hirose_2015,Loos_2018b} which is known to work best for larger systems where the amount of screening is more important. \cite{Jacquemin_2017b,Rangel_2017}
%For comparison purposes, we employ the theoretical best estimates (TBEs) and geometries of Refs.~\onlinecite{Loos_2018a,Loos_2019,Loos_2020b} from which CIS(D), \cite{Head-Gordon_1994,Head-Gordon_1995} ADC(2), \cite{Trofimov_1997,Dreuw_2015} CC2, \cite{Christiansen_1995a} CCSD, \cite{Purvis_1982} and CC3 \cite{Christiansen_1995b} excitation energies are also extracted.
%Various statistical quantities are reported in the following: the mean signed error (MSE), mean absolute error (MAE), root-mean-square error (RMSE), and the maximum positive [Max($+$)] and maximum negative [Max($-$)] errors.
All the static and dynamic BSE calculations have been performed with the software \texttt{QuAcK}, \cite{QuAcK} freely available on \texttt{github}.
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481).}
