\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,txfonts}
\usepackage[version=4]{mhchem}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts}

\usepackage[
	colorlinks=true,
    citecolor=blue,
    breaklinks=true
	]{hyperref}
\urlstyle{same}

\begin{document}	

\title{Spin-Conserved and Spin-Flip Optical Excitations From the Bethe-Salpeter Equation Formalism}

\author{Enzo \surname{Monino}}
	\affiliation{\LCPQ}	 
\author{Pierre-Fran\c{c}ois \surname{Loos}}
	\email{loos@irsamc.ups-tlse.fr}
	\affiliation{\LCPQ}

\begin{abstract}
\alert{Here comes the abstract.}
%\bigskip
%\begin{center}
%	\boxed{\includegraphics[width=0.5\linewidth]{TOC}}
%\end{center}
%\bigskip
\end{abstract}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\alert{Here comes the introduction.}
Unless otherwise stated, atomic units are used, and we assume real quantities throughout this manuscript.
In the following, we consider systems with collinear spins and a spin-independent hamiltonian without contributions such as spin-orbit interaction.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Unrestricted $GW$ formalism}
\label{sec:UGW}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 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$.
Let us denote as $\MO{p_\sig}$ the $p$th orbital of spin $\sig$ (where $\sig =$ $\up$ or $\dw$) and $\e{p_\sig}{}$ its one-electron energy.
In the present context these orbitals can originate from a HF or KS calculation.
In the following, $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.
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 a spin-flip excitation, the hole and particle states, $\MO{i_\sig}$ and $\MO{a_\bsig}$, have opposite spins, $\sig$ and $\bsig$.

%================================
\subsection{The dynamical screening}
%================================
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}
\begin{equation}
		G^{\sig}(\br_1,\br_2;\omega) 
		= \sum_i \frac{\MO{i_\sig}(\br_1) \MO{i_\sig}(\br_2)}{\omega - \e{i_\sig}{} - i\eta}	
		+ \sum_a \frac{\MO{a_\sig}(\br_1) \MO{a_\sig}(\br_2)}{\omega - \e{a_\sig}{} + i\eta}	
\end{equation}
where $\eta$ is a positive infinitesimal.
Based on the spin-up and spin-down components of $G$, one can easily compute the non-interacting polarizability (which is a sum over spins)
\begin{equation}
	\chi_0(\br_1,\br_2;\omega) = - \frac{i}{2\pi} \sum_\sig \int G^{\sig}(\br_1,\br_2;\omega+\omega') G^{\sig}(\br_1,\br_2;\omega') d\omega'
\end{equation}
and subsequently the dielectric function
\begin{equation}
	\epsilon(\br_1,\br-2;\omega) = \delta(\br_1 - \br_2) - \int \frac{\chi_0(\br_1,\br_3;\omega) }{\abs{\br_2 - \br_3}} d\br_3
\end{equation}
where $\delta(\br_1 - \br_2)$ is the Dirac function.
Based on this latter ingredient, one can access the dynamically-screened Coulomb potential
\begin{equation}
	W(\br_1,\br_2;\omega) =  \int \frac{\epsilon^{-1}(\br_1,\br_3;\omega) }{\abs{\br_2 - \br_3}} d\br_3
\end{equation}
which is spin independent as the bare Coulomb interaction $\abs{\br_1 - \br_2}^{-1}$ does not depend on spin coordinates.

Within the $GW$ formalism, the is computed at the RPA level by considering only the manifold of the spin-conserved neutral excitation.
In the orbital basis, the spectral representation of $W$ reads
\begin{multline}
\label{eq:W}
	W_{p_\sig q_\sig,r_\sigp s_\sigp}(\omega) = \ERI{p_\sig q_\sig}{r_\sigp s_\sigp} 
	+ \sum_m \ERI{p_\sig q_\sig}{m}\ERI{r_\sigp s_\sigp}{m} 
	\\
	\times \qty[ \frac{1}{\omega - \Om{m}{\spc,\RPA} + i \eta} - \frac{1}{\omega + \Om{m}{\spc,\RPA} - i \eta} ]
\end{multline}
where the bare two-electron integrals are \cite{Gill_1994}
\begin{equation}
\label{eq:sERI}
		\ERI{p_\sig q_\tau}{r_\sigp s_\taup} = \int  \frac{\MO{p_\sig}(\br_1) \MO{q_\tau}(\br_1) \MO{r_\sigp}(\br_2) \MO{s_\taup}(\br_2)}{\abs{\br_1 - \br_2}} d\br_1 d\br_2
\end{equation}
and the screened two-electron integrals (or spectral weights) are explicitly given by
\begin{equation}
		\ERI{p_\sig q_\sig}{m} = \sum_{ia\sigp} \ERI{p_\sig q_\sig}{r_\sigp s_\sigp} (\bX{m}{\spc,\RPA}+\bY{m}{\spc,\RPA})_{i_\sigp a_\sigp}
\end{equation}
In Eqs.~\eqref{eq:W} and \eqref{eq:sERI}, the RPA spin-conserved neutral excitations $\Om{m}{\spc,\RPA}$ and their corresponding eigenvectors $(\bX{m}{\spc,\RPA}+\bY{m}{\spc,\RPA})$ are obtained by solving a linear response system of the form
\begin{equation}
\label{eq:LR-RPA}
	\begin{pmatrix}
		\bA{}{}	&	\bB{}{}	\\
		-\bB{}{}	&	-\bA{}{}	\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{m}{}	\\
		\bY{m}{}	\\
	\end{pmatrix}
	=
	\Om{m}{}
	\begin{pmatrix}
		\bX{m}{}	\\
		\bY{m}{}	\\
	\end{pmatrix}
\end{equation}
where the expressions of the matrix elements of $\bA{}{}$ and $\bB{}{}$ are specific of the method and of the spin manifold. 
The spin structure of these matrices, though, is general 
\begin{align}
\label{eq:LR-RPA-AB}
	\bA{}{\spc} & = \begin{pmatrix}
		\bA{}{\upup,\upup}	&	\bA{}{\upup,\dwdw}	\\
		\bA{}{\dwdw,\upup}	&	\bA{}{\dwdw,\dwdw}	\\
	\end{pmatrix}
	&
	\bB{}{\spc} & = \begin{pmatrix}
		\bB{}{\upup,\upup}	&	\bB{}{\upup,\dwdw}	\\
		\bB{}{\dwdw,\upup}	&	\bB{}{\dwdw,\dwdw}	\\
	\end{pmatrix}
\\
\label{eq:LR-RPA-AB}
	\bA{}{\spf} & = \begin{pmatrix}
		\bA{}{\updw,\updw}	&	\bO	\\
		\bO					&	\bA{}{\dwup,\dwup}	\\
	\end{pmatrix}
	&
	\bB{}{\spf} & = \begin{pmatrix}
		\bO					&	\bB{}{\updw,\dwup}	\\
		\bB{}{\dwup,\updw}	&	\bO	\\
	\end{pmatrix}
\end{align}
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 
\begin{equation}
\label{eq:small-LR}
	(\bA{}{} - \bB{}{})^{1/2} \cdot (\bA{}{} + \bB{}{}) \cdot (\bA{}{} - \bB{}{})^{1/2} \cdot \bZ{}{} = \bOm{2} \cdot \bZ{}{}
\end{equation}
where the excitation amplitudes are
\begin{equation}
	\bX{}{} + \bY{}{} = \bOm{-1/2} \cdot (\bA{}{} - \bB{}{})^{1/2} \cdot \bZ{}{}
\end{equation}
Within the Tamm-Dancoff approximation (TDA), the coupling terms between the resonant and anti-resonant parts, $\bA{}{}$ and $-\bA{}{}$, are neglected, which consist in setting $\bB{}{} = \bO$.
In such a case, Eq.~\eqref{eq:LR-RPA} reduces to $\bA{}{} \cdot \bX{m}{} = \Om{m}{} \bX{m}{}$.

At the RPA level, the matrix elements of $\bA{}{}$ and $\bB{}{}$ are
\begin{subequations}
\begin{align}
	\label{eq:LR_RPA-A}
	\A{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} & = \delta_{ij} \delta_{ab} \delta_{\sig \sigp}  \delta_{\tau \taup} (\e{a_\tau} - \e{i_\sig}) + \ERI{i_\sig a_\tau}{b_\sigp j_\taup}
	\\ 
	\label{eq:LR_RPA-B}
	\B{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} & = \ERI{i_\sig a_\tau}{j_\sigp b_\taup}
\end{align}
\end{subequations}
from which we obtain the following expressions
\begin{subequations}
\begin{align}
	\label{eq:LR_RPA-Asc}
	\A{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} & = \delta_{ij} \delta_{ab} \delta_{\sig \sigp} (\e{a_\sig} - \e{i_\sig}) + \ERI{i_\sig a_\sig}{b_\sigp j_\sigp}
	\\ 
	\label{eq:LR_RPA-Bsc}
	\B{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} & = \ERI{i_\sig a_\sig}{j_\sigp b_\sigp}
\end{align}
\end{subequations}
for the spin-conserved excitations and 
\begin{subequations}
\begin{align}
	\label{eq:LR_RPA-Asf}
	\A{i_\sig a_\bsig,j_\sig b_\bsig}{\spf,\RPA} & = \delta_{ij} \delta_{ab} (\e{a_\bsig} - \e{i_\sig})
	\\ 
	\label{eq:LR_RPA-Bsf}
	\B{i_\sig a_\bsig,j_\bsig b_\sig}{\spf,\RPA} & = 0
\end{align}
\end{subequations}
for the spin-flip excitations.

%================================
\subsection{The $GW$ self-energy}
%================================
Within the acclaimed $GW$ approximation, \cite{Hedin_1965,Golze_2019} the exchange-correlation (xc) part of the self-energy
\begin{equation}
\begin{split}
	\Sig{}^{\text{xc},\sig}(\br_1,\br_2;\omega) 
	& = \Sig{}^{\text{x},\sig}(\br_1,\br_2) + \Sig{}^{\text{c},\sig}(\br_1,\br_2;\omega) 
	\\
	& = \frac{i}{2\pi} \int G^{\sig}(\br_1,\br_2;\omega+\omega') W(\br_1,\br_2;\omega') e^{i \eta \omega'} d\omega'
\end{split}
\end{equation}
is, like the one-body Green's function, spin-diagonal, and its spectral representation reads
\begin{gather}
	\SigX{p_\sig q_\sig}
	= - \frac{1}{2} \sum_{i\sigp} \ERI{p_\sig i_\sigp}{i_\sigp q_\sig}	 
	\\
\begin{split}
	\SigC{p_\sig q_\sig}(\omega) 
	& = \sum_{im} \frac{\ERI{p_\sig i_\sig}{m} \ERI{q_\sig i_\sig}{m}}{\omega - \e{i_\sig} + \Om{m}{\spc,\RPA} - i \eta}
	\\
	& + \sum_{am} \frac{\ERI{p_\sig a_\sig}{m} \ERI{q_\sig a_\sig}{m}}{\omega - \e{a_\sig} - \Om{m}{\spc,\RPA} + i \eta}			 
\end{split}
\end{gather}
which has been split in its exchange (x) and correlation (c) contributions.
The Dyson equation linking the Green's function and the self-energy holds separately for each spin component, and the quasiparticle energies $\eGW{p_\sig}$ are obtained by solving the frequency-dependent quasiparticle equation
\begin{equation}
	\omega = \e{p_\sig}{} - V_{p_\sigma}^{\xc} + \SigX{p\sigma} + \SigC{p\sigma}(\omega)
\end{equation}
with 
\begin{equation}
	V_{p_\sigma}^{\xc} = \int \MO{p_\sig}(\br) v^{\xc}(\br) \MO{p_\sig}(\br) d\br 
\end{equation}
where $v^{\xc}(\br)$ the Kohn-Sham exchange-correlation.
\alert{Adding the Dyson equation? Introduce linearization of the quasiparticle equation and different degree of self-consistency.}

%================================
\subsection{The Bethe-Salpeter equation formalism}
%================================
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
\begin{multline}
	L^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega) 
	= L_{0}^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega)
	\\
	+ \int L_{0}^{\sig\sigp}(\br_1,\br_4;\br_1',\br_3;\omega) 
	\Xi^{\sig\sigp}(\br_3,\br_5;\br_4,\br_6;\omega)
	\\
	\times L^{\sig\sigp}(\br_6,\br_2;\br_5,\br_2';\omega)
	d\br_3 d\br_4 d\br_5 d\br_6
\end{multline}
where 
\begin{multline}
	L_{0}^{\sig\sigp}(\br_1,\br_2;\br_1',\br_2';\omega) 
	\\
	= \frac{1}{2\pi} \int G^{\sig}(\br_1,\br_2';\omega+\omega') G^{\sig}(\br_1',\br_2;\omega') d\omega'
\end{multline}
is the non-interacting analog of the two-particle correlation function $L$.

Within the $GW$ approximation, the BSE kernel is
\begin{multline}
	i \Xi^{\sig\sigp}(\br_3,\br_5;\br_4,\br_6;\omega) 
	= \frac{\delta(\br_3 - \br_4) \delta(\br_5 - \br_6) }{\abs{\br_3-\br_6}} 
	\\
	- \delta_{\sig\sigp} W(\br_3,\br_4;\omega) \delta(\br_3 - \br_6) \delta(\br_4 - \br_6) 
\end{multline}
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
\begin{equation}
\label{eq:LR-BSE}
	\begin{pmatrix}
		\bA{}{\BSE}		&	\bB{}{\BSE}	\\
		-\bB{}{\BSE}	&	-\bA{}{\BSE}	\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{m}{\BSE}	\\
		\bY{m}{\BSE}	\\
	\end{pmatrix}
	=
	\Om{m}{\BSE}
	\begin{pmatrix}
		\bX{m}{\BSE}	\\
		\bY{m}{\BSE}	\\
	\end{pmatrix}
\end{equation}
Defining $W^{\stat}_{p_\sig q_\sig,r_\sigp s_\sigp} = W_{p_\sig q_\sig,r_\sigp s_\sigp}(\omega = 0)$, we have
\begin{subequations}
\begin{align}
	\label{eq:LR_BSE-A}
	\A{i_\sig a_\tau,j_\sigp b_\taup}{\BSE} & = \A{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig j_\sigp,b_\taup a_\tau}
	\\ 
	\label{eq:LR_BSE-B}
	\B{i_\sig a_\tau,j_\sigp b_\taup}{\BSE} & = \B{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig b_\taup,j_\sigp a_\tau}
\end{align}
\end{subequations}
from which we obtain, at the BSE level, the following expressions
\begin{subequations}
\begin{align}
	\label{eq:LR_BSE-Asc}
	\A{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\BSE} & = \A{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig j_\sigp,b_\sigp a_\sig}
	\\ 
	\label{eq:LR_BSE-Bsc}
	\B{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\BSE} & = \B{i_\sig a_\sig,j_\sigp b_\sigp}{\spc,\RPA} - \delta_{\sig \sigp} W^{\stat}_{i_\sig b_\sigp,j_\sigp a_\sig}
\end{align}
\end{subequations}
for the spin-conserved excitations and 
\begin{subequations}
\begin{align}
	\label{eq:LR_BSE-Asf}
	\A{i_\sig a_\bsig,j_\sig b_\bsig}{\spf,\BSE} & = \A{i_\sig a_\bsig,j_\sig b_\bsig}{\spf,\RPA} - W^{\stat}_{i_\sig j_\sig,b_\bsig a_\bsig}
	\\ 
	\label{eq:LR_BSE-Bsf}
	\B{i_\sig a_\bsig,j_\bsig b_\sig}{\spf,\BSE} & = - W^{\stat}_{i_\sig b_\sig,j_\bsig a_\bsig}
\end{align}
\end{subequations}
for the spin-flip excitations.

%================================
\subsection{Dynamical correction}
%================================
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}

\begin{multline}
	\widetilde{W}_{p_\sig q_\sig,r_\sigp s_\sigp}(\omega) = \ERI{p_\sig q_\sig}{r_\sigp s_\sigp} 
	+ \sum_m \ERI{p_\sig q_\sig}{m}\ERI{r_\sigp s_\sigp}{m} 
	\\
	\times \qty[ \frac{1}{\omega - (\e{s_\sigp}{} - \e{q_\sig}{}) - \Om{m}{\spc,\RPA} + i \eta} + \frac{1}{\omega - (\e{r_\sigp}{} - \e{p_\sig}{}) - \Om{m}{\spc,\RPA} + i \eta} ]
\end{multline}

\begin{equation}
\label{eq:LR-dyn}
	\begin{pmatrix}
		\bA{}{\dBSE}(\omega)		&	\bB{}{\dBSE}(\omega)	
		\\
		-\bB{}{\dBSE}(-\omega)	&	-\bA{}{\dBSE}(-\omega)	
		\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{m}{\dBSE}	\\
		\bY{m}{\dBSE}	\\
	\end{pmatrix}
	=
	\Om{m}{\dBSE}
	\begin{pmatrix}
		\bX{m}{\dBSE}	\\
		\bY{m}{\dBSE}	\\
	\end{pmatrix}
\end{equation}

\begin{subequations}
\begin{align}
	\label{eq:LR_dBSE-A}
	\A{i_\sig a_\tau,j_\sigp b_\taup}{\dBSE}(\omega) & = \A{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} \widetilde{W}_{i_\sig j_\sigp,b_\taup a_\tau}(\omega) 
	\\ 
	\label{eq:LR_dBSE-B}
	\B{i_\sig a_\tau,j_\sigp b_\taup}{\dBSE}(\omega) & = \B{i_\sig a_\tau,j_\sigp b_\taup}{\RPA} - \delta_{\sig \sigp} \widetilde{W}_{i_\sig b_\taup,j_\sigp a_\tau}(\omega) 
\end{align}
\end{subequations}


\begin{multline}
\label{eq:LR-PT}
	\begin{pmatrix}
		\bA{}{\dBSE}(\omega)		&	\bB{}{\dBSE}(\omega)	\\
		-\bB{}{\dBSE}(-\omega)	&	-\bA{}{\dBSE}(-\omega)	\\
	\end{pmatrix}
	\\
	=
	\begin{pmatrix}
		\bA{}{(0)}	&	\bB{}{(0)}	
		\\
		-\bB{}{(0)}	&	-\bA{}{(0)}	
		\\
	\end{pmatrix}
	+
	\begin{pmatrix}
		\bA{}{(1)}(\omega)		&	\bB{}{(1)}(\omega)	\\
		-\bB{}{(1)}(-\omega)	&	-\bA{}{(1)}(-\omega)	\\
	\end{pmatrix}
\end{multline}
with
\begin{subequations}
\begin{align}
	\label{eq:BSE-A0}
	\A{i_\sig a_\tau,j_\sigp b_\taup}{(0)} & = \A{i_\sig a_\tau,j_\sigp b_\taup}{\BSE}
	\\ 
	\label{eq:BSE-B0}
	\B{i_\sig a_\tau,j_\sigp b_\taup}{(0)} & = \B{i_\sig a_\tau,j_\sigp b_\taup}{\BSE}
\end{align}
\end{subequations}
and 
\begin{subequations}
\begin{align}
	\label{eq:BSE-A1}
	\A{i_\sig a_\tau,j_\sigp b_\taup}{(1)}(\omega) & = - \delta_{\sig \sigp} \widetilde{W}_{i_\sig j_\sigp,b_\taup a_\tau}(\omega) + \delta_{\sig \sigp} W^{\stat}_{i_\sig j_\sigp,b_\taup a_\tau}
	\\ 
	\label{eq:BSE-B1}
	\B{i_\sig a_\tau,j_\sigp b_\taup}{(1)}(\omega) & = - \delta_{\sig \sigp} \widetilde{W}_{i_\sig b_\taup,j_\sigp a_\tau}(\omega) + \delta_{\sig \sigp} W^{\stat}_{i_\sig b_\taup,j_\sigp a_\tau}
\end{align}
\end{subequations}

\begin{subequations}
\begin{gather}
	\Om{m}{\dBSE} = \Om{m}{(0)} + \Om{m}{(1)} + \ldots
	\\
	\begin{pmatrix}
		\bX{m}{\dBSE}	\\
		\bY{m}{\dBSE}	\\
	\end{pmatrix}
	= 
	\begin{pmatrix}
		\bX{m}{(0)}	\\
		\bY{m}{(0)}	\\
	\end{pmatrix}
	+
	\begin{pmatrix}
		\bX{m}{(1)}	\\
		\bY{m}{(1)}	\\
	\end{pmatrix}
	+ \ldots
\end{gather}
\end{subequations}

\begin{equation}
\label{eq:LR-BSE-stat}
	\begin{pmatrix}
		\bA{}{(0)}	&	\bB{}{(0)}	\\
		-\bB{}{(0)}	&	-\bA{}{(0)}	\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{S}{(0)}	\\
		\bY{S}{(0)}	\\
	\end{pmatrix}
	=
	\Om{m}{(0)}
	\begin{pmatrix}
		\bX{m}{(0)}	\\
		\bY{m}{(0)}	\\
	\end{pmatrix}
\end{equation}

\begin{equation}
\label{eq:Om1}
	\Om{m}{(1)} =
	\T{\begin{pmatrix}
		\bX{m}{(0)}	\\
		\bY{m}{(0)}	\\
	\end{pmatrix}}
	\cdot
	\begin{pmatrix}
		\bA{}{(1)}(\Om{m}{(0)})	&	\bB{}{(1)}(\Om{m}{(0)})	\\
		-\bB{}{(1)}(-\Om{m}{(0)})	&	-\bA{}{(1)}(-\Om{m}{(0)})	\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{m}{(0)}	\\
		\bY{m}{(0)}	\\
	\end{pmatrix}
\end{equation}

\begin{equation}
\label{eq:Om1-TDA}
	\Om{S}{(1)} = \T{(\bX{m}{(0)})} \cdot \bA{}{(1)}(\Om{m}{(0)}) \cdot \bX{m}{(0)}
\end{equation}

\begin{equation}
\label{eq:Z}
	Z_{m} = \qty[ 1 - \T{(\bX{m}{(0)})} \cdot \left. \pdv{\bA{}{(1)}(\Om{m}{})}{\Om{S}{}} \right|_{\Om{m}{} = \Om{m}{(0)}} \cdot \bX{m}{(0)} ]^{-1}
\end{equation}

\begin{equation}
	\Om{m}{\dBSE} = \Om{m}{(0)} + Z_{m} \Om{m}{(1)}
\end{equation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Oscillator strengths}
\label{sec:os}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
For the spin-conserved transition, the $x$ component of the transition dipole moment is
\begin{equation}
	\mu_{x,m}^{\spc} = \sum_{ia\sig} (i_\sig|x|a_\sig)(\bX{m}{\spc}+\bY{m}{\spc})_{i_\sig a_\sig}
\end{equation}
with 
\begin{equation}
	(p_\sig|x|q_\sigp) = \int \MO{p_\sig}(\br) \, x \, \MO{q_\sigp}(\br) d\br
\end{equation}
and the total oscillator strength is given by
\begin{equation}
	f_{m}^{\spc} = \frac{2}{3} \Om{m}{\spc} \qty[ \qty(\mu_{x,m}^{\spc})^2 + \qty(\mu_{x,m}^{\spc})^2 + \qty(\mu_{x,m}^{\spc})^2 ]
\end{equation}
For spin-flip transitions, we have $f_{m}^{\spf} = 0$ as the transition matrix elements $(i_\sig|x|a_\bsig)$ vanish.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Spin contamination}
\label{sec:spin}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{equation}
	\expval{S^2}_m = \expval{S^2}_0 + \Delta \expval{S^2}_m
\end{equation}

\begin{equation}
	\expval{S^2}_{0} 
	= \frac{n_{\up} - n_{\dw}}{2} \qty( \frac{n_{\up} - n_{\dw}}{2} + 1 )
	+ n_{\dw} - \sum_p (p_{\up}|p_{\dw})^2
\end{equation}
where
\begin{equation}
	(p_\sig|q_\sigp) = \int \MO{p_\sig}(\br) \MO{q_\sigp}(\br) d\br
\end{equation}
is the overlap between spin-up and spin-down orbitals.

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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational details}
\label{sec:compdet}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results}
\label{sec:res}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% TABLE I %%%
%\begin{table}
%
%\end{table}
%%% %%% %%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
\label{sec:ccl}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgements{
We would like to thank Xavier Blase and Denis Jacquemin for insightful discussions. 
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).}
%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Data availability statement}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The data that supports the findings of this study are available within the article and its supplementary material.

%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{sfBSE}
%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}
