\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.
The spin-$\sig$ component of the one-body Green's function reads
\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 it, 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 functions.
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.
Therefore, within the $GW$ formalism, the dynamical screening $W(\omega)$ is computed at the RPA level using the spin-conserved neutral excitations.
In the orbital basis, the spectral representation of $W(\omega)$ 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 the following linear response system  
\begin{equation}
\label{eq:LR-RPA}
	\begin{pmatrix}
		\bA{}{\spc,\RPA}	&	\bB{}{\spc,\RPA}	\\
		-\bB{}{\spc,\RPA}	&	-\bA{}{\spc,\RPA}	\\
	\end{pmatrix}
	\cdot
	\begin{pmatrix}
		\bX{m}{\spc,\RPA}	\\
		\bY{m}{\spc,\RPA}	\\
	\end{pmatrix}
	=
	\Om{m}{\spc,\RPA}
	\begin{pmatrix}
		\bX{m}{\spc,\RPA}	\\
		\bY{m}{\spc,\RPA}	\\
	\end{pmatrix}
\end{equation}

The spin structure of matrices $\bA{}{}$ and $\bB{}{}$ 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}
and does not only apply to the RPA but also to RPAx (\ie, RPA with exchange), BSE and TD-DFT.
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, 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 read
\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 exchange (x) and correlation (c) parts for convenience.
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.

%================================
\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.
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 counterpart 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}

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}
%================================

\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}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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{sf-BSE}
%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}
