Suppose we have the following Slater determinant: \begin{equation} | \Psi \rangle = \prod \limits_{i,i'} a^+_{i\alpha} a^+_{i'\beta} | \rangle \end{equation} where $a^+_{i\alpha}$ creates an electron in state $i$ with spin $\alpha$ and, in general, $i \neq i'$. I want to evaluate $\langle \Psi | S^2 | \Psi \rangle$ using second quantization. We can express the $S^2$ operator as \begin{equation} S^2 = S_- S_+ + S_z (S_z +1) \end{equation} with \begin{equation} S_- = \sum_p a^+_{p\beta} a_{p \alpha} \hspace{1.5cm} S_+ = \sum_p a^+_{p\alpha} a_{p \beta}. \end{equation} Since $|\Psi \rangle$ is an eigenfunction of $S_z$, evaluating $\langle \Psi | S_z | \Psi \rangle$ terms becomes trivial and the problem reduces to the evaluation of $\langle \Psi | S_- S_+ | \Psi \rangle$. In the standard restricted Hartree-Fock method, $i = i'$ and it is easy to show that $\langle \Psi | S_- S_+ | \Psi \rangle = 0$ using the canonical anticommutation relations. When $i \neq i'$ (unrestricted Hartree-Fock) we must have that $\langle \Psi | S_- S_+ | \Psi \rangle = N_\beta - \text{Tr}(PQ)$ (this I verified in books) where $N_\beta$ is the number of $\beta$ electrons and $P_{ij} = \langle \Psi | a^+_{j\alpha} a_{i\alpha} | \Psi \rangle$ and $Q_{ij} = \langle \Psi | a^+_{j\beta} a_{i\beta} | \Psi \rangle$.
My question is then, specifically, how do we obtain $\langle \Psi | S_- S_+ | \Psi \rangle = N_\beta - \text{Tr}(PQ)$ using the anticommutation relations of second quantization? If I try to use these anticommutation relations just in the way they are written in textbooks I do not get the right answer. Clearly, I am getting something wrong or some special consideration needs to be taken into account when $i \neq i'$. If someone can show me how to correctly use second quantization to get the right answer here, I would appreciate it a lot.