Is it always true that the Hartree-Fock slater determinant is the minimum energy single slater determinant solution to a given molecular Hamiltonian $H$?
By construction, any determinant $|\chi_{i_1}\ldots \chi_{i_N}\rangle$ where the $\chi_{i_j}$ are distinct solutions of the Hartree-Fock equation $f(x) \chi(x) = \epsilon_i \chi(x)$ will have stationary energy $\langle \chi_{i_1}\ldots \chi_{i_N} | H |\chi_{i_1}\ldots \chi_{i_N}\rangle$ (where stationarity is with respect to perturbing the orbitals). The Hartree-Fock slater determinant $|\Phi_0\rangle =|\chi_1\ldots\chi_N\rangle$ is constructed by choosing the lowest $N$ energy solutions of the Hartree-Fock equation.
Now, $|\chi_{i_1}\ldots \chi_{i_N}\rangle$ is an eigenstate of the Hartree-Fock Hamiltonian $H^{HF}$ with energy $\sum_{i=1}^N \epsilon_{i_j}$, so constructing $|\Phi_0\rangle$ from the $N$ lowest energy solutions of the HF equation clearly minimizes the energy of $H^{HF}$. But the molecular Hamiltonian $H$ is given by $$H = H^{HF} + \sum_{i=1}^N\sum_{j>i}^N \frac{1}{r_{ij}} - \sum_{i=1}^N v^{HF}(i)$$ where $v^{HF}$ is the HF potential. This satisfies $$\langle{\Phi_0}|H|\Phi_0\rangle = \sum_{i=1}^N \epsilon_i - \frac{1}{2}\sum_{i=1}^N \sum_{j=1}^N \langle ij|\,|ij\rangle $$ where $$\langle ij|\,|ij\rangle:=\int dx_1 dx_2 \chi_i^*(x_1)\chi_j^*(x_2)\frac{1}{r_{12}}(\chi_i(x_1)\chi_j(x_2) - \chi_j(x_1)\chi_i(x_2))$$ Due to the $\langle ij|\,|ij\rangle$ term, it is not clear to me that constructing $\Phi_0$ from the $N$ lowest energy solutions of the HF equations still results in the lowest energy single slater determinant for $H$. An excited determinant would have a higher energy contribution to the $\sum_i \epsilon_i$ term, but could this is offset by a greater deduction from the $\langle ij|\,|ij\rangle$ term?