It appears that the first explicit statement about the symmetry-breaking of Unrestricted Hartree-Fock (UHF) was by P. O. Löwdin in Discussion on The Hartree-Fock Approximation, P. Lykos and G. W. Pratt,
Rev. Mod. Phys. 35, 496 – Published 1 July 1963.
Löwdin writes:
"I would like to comment on some peculiarities with respect to the symmetry properties. ...
Confusion may arise from the fact that the exact eigenfunction $\Psi$ and the approximate eigenfunction in the form of a Slater determinant may have rather different properties. For instance,
if $\Lambda$ is a normal constant of motion satisying the relation $H\Lambda = \Lambda H$,
then every eigenfunction
to $H$ is automatically an eigenfunction to $\Lambda$ or (in the case of a degenerate energy level) may be
chosen in that way, so that
$H\Psi = E\Psi$
$\Lambda\Psi = \lambda\Psi$
... On the other hand, if one drops the symmetry constraint and considers only the relation
$\delta \langle D\vert H\vert D\rangle=0$
one obtains a, nonrestricted Hartree Fock scheme, and the solution $D$ corresponding to the
absolute minimum has now usually lost its eigenvalue property with respect to $\Lambda$, i.e., the corresponding Hartree-Fock functions are no longer symmetry-adapted....
In my opinion, the Hartree-Fock scheme based on a single Slater determinant $D$ is in a dilemma
with respect to the symmetry properties and the normal constants of motion $\Lambda$. The assumption
that $D$ should be symmetry-adapted or an eigenfunction to $H$ leads to an energy $\langle H\rangle $ high above
the absolute minimum, and the energy difference amounts to at least 1eV per electron pair and
more. In the sense of Eckart's criterion [C. Eckart, Phys. Rev. 36, 877 (1930); B. A. Lengyel,
J. Math. Analysis Appl. 5, 451 (1962)], the absolute minimum of $\langle D\vert H\vert D\rangle$ leads certainly to a
better wave function $D$, but the symmetry properties are now lost and the determinant is a
"mixture" of components of various symmetry types."