The equations are often rewritten in terms of a so-called density matrix, which is defined as $$ \begin{equation} \rho_{\gamma\delta}=\sum_{i=1}^{N}\langle\gamma|i\rangle\langle i|\delta\rangle = \sum_{i=1}^{N}C_{i\gamma}C^*_{i\delta}. \tag{18} \end{equation} $$ It means that we can rewrite the Hartree-Fock Hamiltonian as $$ \hat{h}_{\alpha\beta}^{HF}=\epsilon_{\alpha}\delta_{\alpha,\beta}+ \sum_{\gamma\delta} \rho_{\gamma\delta}\langle \alpha\gamma|V|\beta\delta\rangle_{AS}. $$ It is convenient to use the density matrix since we can precalculate in every iteration the product of two eigenvector components \( C \).