  Spin-flip methods applied to excited-state approaches like the Bethe-Salpeter
Equation allow access to the excitation energies of open-shell systems, such as
molecules and defects in solids. The eigenstates of these solutions, however,
are generally not eigenstates of the spin operator $\hat{S}^2$. Even for simple
cases where the excitation vector is expected to be, for example, a triplet
state, the value of $\langle \hat{S}^2 \rangle$ may be found to differ from
2.00; this difference is called ``spin contamination.'' The expectation values
$\langle \hat{S}^2 \rangle$ must be computed for each excitation vector, to
assist with the characterization of the particular excitation and to determine
the amount of spin contamination of the state. Our aim is to provide for the
first time in the spin-flip methods literature a comprehensive resource on the
derivation of the formulas for $\langle \hat{S}^2 \rangle$ as well as its
computational implementation. After a brief discussion of the theory of the
Spin-Flip Bethe-Salpeter Equation and some examples further illustrating the
need for calculating $\langle \hat{S}^2 \rangle$, we present the derivation for
the general equation for computing $\langle \hat{S}^2 \rangle$ with the
eigenvectors from an SF-BSE calculation, how it is implemented in a Python
script, and timing information on how this calculation scales with the size of
the SF-BSE Hamiltonian.
