Chemistry - Hartree-Fock dipole moment
Solution 1:
I originally wrote this as a comment, but I will post it as an answer as I'm fairly certain this is most of the answer.
In general, there are two distinct types of electron correlation which can be described. Exchange correlation, which results from the indistinguishability of electrons, and Coulomb correlation which correlates the motions of electrons by considering individual repulsions between electrons.
Hartree-Fock (HF) theory accounts for the exchange interaction, but does not have Coulomb correlation as HF is a mean-field theory. That is, each electron experiences the average potential of the other electrons.
So, there is a sort of imbalance here from the beginning. It is known that the presence of this exchange interaction, which only happens for electrons of the same spin, results in two phenomena called Fermi holes and Fermi heaps. The references in that Wikipedia page are quite good. Essentially, a Fermi heap results when you have electrons of opposite spin. Because they must have a spatially symmetric wavefunction, the ends up being a greater probability of finding both electrons at the same location. This sounds silly, but read the references in the wikipedia page and you'll see what I mean. A fermi hole results when you have two electrons of the same spin. These must have an antisymmetric spatial wavefunction, so there is no probability of finding them in the same location.
If we think about these two effects in terms of electrical repulsion, the Fermi hole is a lower energy configuration than the Fermi heap. Because we almost always have both spin-up and spin-down electrons in multi-electron systems, both of these phenomena will take place, and the overall electron density will be some kind of average of all of these. (This is only really a conceptual framework.) If we take a simple case of a homonuclear diatomic, then these heaps and holes will all average out and we will get a symmetric electron density. If, on the other hand, we take a heteronuclear diatomic, then we see that there will be a tendency for these Fermi heaps to form over the atom with the larger nuclear charge just based on Coulombic attraction.
From here, I believe it is fairly clear that Coulomb correlation will tend to make the effect of these Fermi heaps smaller, and hence the overall separation of charge will decrease rather than increase, which explains why HF tends to overestimate dipole moments. By the same reasoning, the electron density is already polarized due to this heaping effect, so the system is not likely to be as polarizable as if this electron density were allowed to relax as is the case when Coulomb correlation is included.
This is all a very conceptual way of thinking about this, but I don't know that the math describing HF will be all that much more useful.
Solution 2:
jheindel has given an interesting and valid perspective on the issue, but I believe that one does not necessarily have to go holes and heaps to understand the situation if configuration interaction (CI) will do.
The explanation that I have been given comes down to: restricted HF yields wave-functions that have a too high portion of ionic character. This ionic character also manifests itself in the dipole moments. The overly ionic character is known to cause the absurd dissociation curves of RHF. Note that unrestricted HF can partially repair this, at the cost of introducing spin contamination. At typical bond distances, UHF will collapse to the RHF solution anyway.
For an extensive discussion, I suggest: F. Jensen, Introduction to Computational Chemistry, Wiley, 2nd Edition (2007), Chapter 4.3: Illustrating how [Configuration Interaction] Accounts for Electron Correlation, and the RHF Dissociation Problem. Of that discussion, I will write up the key points here.
Let's assume the following model system: $\ce{H2}$, minimal basis set. RHF will yield the following ground state wave-function (and since it's RHF; this is true for all nucleus-nucleus distances): $$ \Phi_0 = [\chi_A(1) + \chi_B(1)][\chi_A(2) + \chi_B(2)] \\ = \underbrace{\chi_A(1)\chi_A(2) + \chi_B(1)\chi_B(2)}_{\mathrm{Ionic \ Part}} + \underbrace{\chi_A(1)\chi_B(2) + \chi_B(1)\chi_A(2)}_{\mathrm{Covalent \ Part}} $$ $\chi$ are the AOs. Atom/AO labels $A,B$. Electron labels $1,2$, where given and hopefully obvious otherwise. This is after the spin functions have been integrated out/away and the LCAO ansatz has been used. The ionic part refers to a $\ce{H+}\ce{H-}/\ce{H-}\ce{H+}$ situation, whereas the covalent part looks like two hydrogen atoms next to each other.
Several other (excited) configurations exist, but only one of them is relevant (due to symmetry and spin considerations), which is: $$ \Phi_1 = [\chi_A(1) - \chi_B(1)][\chi_A(2) - \chi_B(2)] \\ = {\chi_A\chi_A + \chi_B\chi_B} - {\chi_A\chi_B - \chi_B\chi_A} $$ They can be combined by linear combination (obeying normalization) as in: $$ \Phi_{\mathrm{CI}} = c_0 \Phi_0 + c_1 \Phi_1 $$ which allows for the down-scaling of the ionic term as appropriate for the nucleus-nucleus distance.
Jensen goes on to note that the too-high dissociation energy of RHF manifests itself also in short bond lengths and high vibrational frequencies because the potential hypersurface is distorted. The too-high ionic character finally leads to high dipole moments and atomic charges (as in Mulliken charges etc.).
Note: I am not bored enough for quantitative analysis yet.