In their landmark paper of 1964, Inhomogeneous Electron Gas, Hohenberg and Kohn wrote down:
We shall be considering a collection of an arbitrary number of electrons, enclosed in a large box and moving under the influence of an external potential $v(\mathbf{r})$ and the mutual Coulomb repulsion. The Hamiltonian has the form
$$H=T+V+U \tag{1}$$
where[...atomic units are used.]
\begin{eqnarray*} T &\equiv& \frac{1}{2} \int \nabla \psi^*(\mathbf{r}) \ \nabla\psi(\mathbf{r}) \ d \mathbf{r} \tag{2} \\ V &\equiv& \int v(\mathbf{r}) \ \psi^*(\mathbf{r}) \ \psi(\mathbf{r}) \ d\mathbf{r} \tag{3} \\ U &=& \frac{1}{2} \int \frac{1}{|\mathbf{r}-\mathbf{r}^{\prime}|} \ \psi^\ast(\mathbf{r}) \ \psi^*(\mathbf{r}^{\prime}) \ \psi(\mathbf{r}^{\prime}) \ \psi(\mathbf{r}) \ d \mathbf{r} \ d \mathbf{r}^{\prime} \tag{4} \end{eqnarray*}
After applying the Born-Oppenheimer approximation, one obtains the following $N$-electron Hamiltonian,
$$ \hat{H} = - \frac{1}{2} \sum_{i=1}^N \nabla_i^2 + \sum_{i=1}^N v(\hat{\mathbf{r}_i}) + \frac{1}{2} \sum_{i \neq j} \frac{1}{|\hat{\mathbf{r}}_i-\hat{\mathbf{r}}_j|} $$
and the time-independent Schrodinger equation,
$$ \hat{H} \ \psi(\mathbf{r}_1, \mathbf{r}_2, ..., \mathbf{r}_N) = E \ \psi(\mathbf{r}_1, \mathbf{r}_2, ..., \mathbf{r}_N) $$
My question is how to obtain Eq. (4) given the above Hamiltonian and Schrodinger equation?
I know Eq. (4) represents the classical Hartree energy and it can be deduced from Electrodynamics arguments (of course not in terms of the wavefunction form depicted above but in terms of the charge density of a charge distribution), BUT how to obtain it from the expectation value of the above Hamiltonian, i.e.,
$$ U = \int_{\mathbb{R}^3} \ldots \int_{\mathbb{R}^3} \ d \mathbf{r}_1 \ldots d \mathbf{r}_N \ \psi^*(\mathbf{r}_1, ..., \mathbf{r}_N) \ \left(\frac{1}{2} \sum_{i \neq j} \frac{1}{|\hat{\mathbf{r}}_i-\hat{\mathbf{r}}_j|}\right) \ \psi(\mathbf{r}_1, ..., \mathbf{r}_N) $$