Assume that we've "solved" the RHF Hartree-Fock equations in some finite basis, that is, we have arrived at a self-consistent set of coefficients $c_{\mu{}i}$ such that each spatial molecular orbital is expanded as $\chi_i = \sum_{\mu=1}^K c_{\mu{}i}\varphi_i$. By construction, this set of orbitals is orthonormal, and forms a collection of eigenfunctions for the Fock operators in the given basis, so $f_i |\chi_i\rangle = \varepsilon_i |\chi_i\rangle$ where $\varepsilon_i$ is the $i$th orbital energy.
Question: what, if anything, can we say about these orbitals in terms of their relationship to the "true" orbitals in the basis set limit, denoted $\chi_i^\mathrm{BS}$? Here, I mean those orbitals that would be obtained if the calculation were to be extended to an (infinite) complete basis set. (I realise that the assumed existence of a single unique set of true orbitals may not be justified...)
For example, is it possible to estimate, bound, or otherwise describe the overlap $\langle \chi_i | \chi_i^\mathrm{BS} \rangle$ between an approximate solution and the true solution? How about matrix elements such as $\langle \chi_i | \mathcal{H} | \chi_i^\mathrm{BS} \rangle$, where $\mathcal{H}$ is the electronic Hamiltonian, or indeed any self-adjoint operator?
I suspect/hope that some kind of a statement towards the first part might be possible, given that the Roothan-Hall equations are effectively a collection of Galerkin equations. One of the standard hand-waving statements about solutions obtained using Galerkin methods is that they are in some sense projections of the true solution onto finite-dimensional subspaces, and have residual vectors that are "orthogonal" to that subspace(for a potentially complicated definition of "orthogonal").
(I'm also interested in answers that would apply in any of the HF settings, i.e. UHF, ROHF, etc.)