I am not able to understand the math of this.
After taking the Slater determinant as a wavefunction in Hartree-Fock, does the procedure to find wavefunctions remain the same as in the Hartree method i.e. the self-consistent method, or does it change?
Both Hartree and Hartree-Fock have a self-consistent field method, as does density functional theory, too; the difference is just that Hartree theory uses a bosonic wave function (symmetric with interchange of particles), while Hartree-Fock uses a fermionic wave function i.e. a Slater determinant that is the correct one for electrons.
In both cases, the self-consistent field method is obtained by the calculus of variations, which gives you a condition for achieving an extremum of the energy functional; e.g. in finite-basis Hartree-Fock, this is the Roothaan equation FC = SCE that gives you the one-particle states. (Note that the extremum does not have to be a minimum, it can also be a saddle point!)
Because the Hamiltonian is a two-particle operator, the potential the electrons move in is dependent on the movement of the electrons, or F = F(C). This means that you have to find a self-consistent solution. It turns out that the Roothaan procedure of solving FC=SCE for C and recomputing F=F(C) usually does not converge, and you have to stabilize the iterations e.g. by damping or extrapolation.
This is discussed in several textbooks, as well as our recent overview.
Edit: as is pointed out by André below, the Hartree method actually is inconsistent: the proper bosonic wave function is indistinguishable, whereas the Hartree method places specific particles on specific orbitals. A proper bosonic model would have a Fock matrix similar to Hartree-Fock; the difference is just that the exchange term would have a plus sign instead of a minus sign as in Hartree-Fock. For bosonic systems one often uses the Gross-Pitaevskii model where all the particles occupy the lowest orbital (exclusion principle is only for fermions).
The procedure to find the Hartree equations and the Hartree-Fock equations is very similar, we have to minimize the expectation value of the Hamiltonian under the orthonormalization constraint. However, both methods differ in the form of the wavefunction.
In the Hartree method, the total wavefunction is a Hartree product: $$ \Phi(\mathbf{r}_1,\mathbf{r}_2,\dots,\mathbf{r}_N)=\prod_{i=1}^{N}\phi_i(\mathbf{r}_i) $$ The constraint is the normalization of every orbital, so then, we can construct the functional: $$ L[\{\phi_i,\phi_i^*\}]=\langle\hat{H}\rangle-\sum_{i=1}^{N}\varepsilon_i\left(\langle\phi_i|\phi_i\rangle-1\right) $$ where $\varepsilon_i$ are Lagrange multipliers. Now, the expectation value of the Hamiltonian is $$ \begin{aligned} \langle\hat{H}\rangle&=\sum_{i=1}^{N}\langle\phi_i|\frac{-\nabla_i^2}{2}-\frac{Z}{r_i}|\phi_i\rangle\left(\prod_{j\neq i}^{N}\langle\phi_j|\phi_j\rangle\right)+\frac{1}{2}\sum_{i=1}^{N}\sum_{j\neq i}^{N}\langle\phi_i\phi_j|r_{ij}^{-1}|\phi_i\phi_j\rangle\left(\prod_{k\neq i,j}^{N}\langle\phi_k|\phi_k\rangle\right) \\ &=\sum_{i=1}^{N}h_{ii}+\frac{1}{2}\sum_{i=1}^{N}\sum_{j\neq i}^{N}J_{ij} \end{aligned} $$ The functional $L$ depends on the set of all the $\phi_i$ and its complex conjugates, so we can vary either of them arbitrarily. For instance, a variation in a $\phi_k^*$ can be written as $$ L[\phi_1,\phi_1^*,\dots,\phi_k,\phi_k^*+\lambda\delta\phi_k^*,\dots,\phi_N,\phi_N^*]\equiv L[\phi_k^*+\lambda\delta\phi_k^*] $$ and the difference between the original functional and this is $$ L[\phi_k^*+\lambda\delta\phi_k^*]-L[\{\phi_i,\phi_i^*\}]=\lambda\left(\langle\delta\phi_k|\frac{-\nabla_k^2}{2}-\frac{Z}{r_k}|\phi_k\rangle+\sum_{j\neq i}^{N}\langle\delta\phi_k\phi_j|r_{kj}^{-1}|\phi_k\phi_j\rangle-\varepsilon_k\langle\delta\phi_k|\phi_k\rangle\right) $$ Dividing by $\lambda$ and taking the limit where $\lambda\to 0$, we have that an extremum of the functional gives the Hartree equations (in summary, we look for $\delta L/\delta\phi_k^*=0$): $$ \left(\frac{-\nabla_k^2}{2}-\frac{Z}{r_k}+\sum_{j\neq k}^{N}\int\frac{|\phi_j(\mathbf{r}_j)|^2}{r_{jk}}d\mathbf{r}_j\right)\phi_k(\mathbf{r}_k)=\varepsilon_k\phi_k(\mathbf{r}_k) $$ In short, the total energy can be written in terms of the orbital energies $\varepsilon_i$ as $$ E_H=\sum_{i=1}^{N}\varepsilon_i-\frac{1}{2}\sum_{i=1}^{N}\sum_{j\neq i}^{N}J_{ij} $$ The problem of the Hartree method is that the wavefunction is not antisymmetric with respect to the exchange of particles and, more fundamentally, it doesn't consider the electrons as indistinguishable particles. Also, in the original formulation, it doesn't consider spin.
We can fix these problems by using an Slater determinant as a wavefunction. Now, we have $$ \Psi(\mathbf{x}_1,\dots,\mathbf{x}_N)=\frac{1}{\sqrt{N!}}\sum_{n=1}^{N!}(-1)^{n_p}\hat{P}_n(\psi_i(\mathbf{x}_1),\psi_j(\mathbf{x}_2)\dots\psi_k(\mathbf{x}_N)), $$ where $\hat{P}_n$ is the permutation operator, $n_p$ is the number of required transpositions to get an specific permutation, and $\mathbf{x}_i$ are spatial and spin coordinates.
Now, we construct the functional for this situation, with $\langle\hat{H}\rangle=\sum_{i=1}^{N}h_{ii}+\frac{1}{2}\sum_{i=1}^{N}\sum_{j\neq i}^{N}J_{ij}-K_{ij}$ and the orthonormality condition of the orbitals. $$ L[\{\psi_i,\psi_i^*\}]=\langle\hat{H}\rangle-\sum_{i=1}^{N}\sum_{j=1}^{N}\varepsilon_{ij}\left(\langle\psi_i|\psi_j\rangle-\delta_{ij}\right) $$ Briefly, using the same method we get to the Hartree-Fock equations $$ \left(\frac{-\nabla_k^2}{2}-\frac{Z}{r_k}+\sum_{j\neq k}^{N}\int\frac{|\psi_j(\mathbf{x}_j)|^2}{r_{jk}}d\mathbf{x}_j-\int\frac{\psi_j^*(\mathbf{x}_j)\hat{P}_{kj}\psi_j(\mathbf{x}_j)}{r_{kj}}d\mathbf{x}_j\right)\psi_k(\mathbf{x}_k)=\sum_{j=1}^{N}\varepsilon_{kj}\psi_j(\mathbf{x}_k) $$ Finally, doing an unitary transformation we get the canonical form of the HF equations and the total energy in terms of the spin-orbitals is $$ E_{HF}=\sum_{i=1}^{N}\varepsilon_i-\frac{1}{2}\sum_{i=1}^{N}\sum_{j\neq i}^{N}J_{ij}-K_{ij} $$
As you see, the method for finding an extremum of the energy with respect of the orbitals with the constraint of orthonormalization is similar. However, considering an Slater determinant instead of a Hartree product improves a lot the results in terms of physical significance and numerical accuracy. An Slater determinant considers that the electrons are antisymmetric with respect to exchange of coordinates and they are indistinguishable particles. Moreover, we have another term involved in the energy: the exchange integral, which is always positive and lowers the energy when it's not zero. In a pseudo-classical interpretation of the determinantal energies (as Szabo & Ostlund called in their book) we see that exchange terms arise only when we have electrons with the same spin, but, obviously, in different spin-orbitals.
