I will treat the electric field semiclassically as a function of the position of the electron, so that $$ \langle \mathbf{x_1}|\hat{\mathbf{E}}(\mathbf{x})|\mathbf{x_2} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} -e\frac{\mathbf{x}-\mathbf{x_1}}{\lVert \mathbf{x}-\mathbf{x_1} \rVert}^3 \right). $$$$ \langle \mathbf{x_1}|\hat{\mathbf{E}}(\mathbf{x})|\mathbf{x_2} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} -e\frac{\mathbf{x}-\mathbf{x_1}}{\lVert \mathbf{x}-\mathbf{x_1} \rVert^3} \right). $$ As before, it is implied that a charge $Ze$ of infinite mass lies at the origin. Similarly for the charge density, \begin{equation} \langle {\mathbf{x_1}|\hat{\rho}(\mathbf{x})|\mathbf{x_2}} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze\delta^{3}(\mathbf{x}) - e\delta^{3}(\mathbf{x}-\mathbf{x_1}) \right). \tag{1} \end{equation} The average electric field of an arbitrary state $\Psi$ is then given by $$ \langle \mathbf{E}(\mathbf{x}) \rangle = \langle {\Psi|\hat{\mathbf{E}}(\mathbf{x})|\Psi} \rangle = Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} + \int -e\lvert \Psi(\mathbf{x^\prime}) \rvert^2 \frac{\mathbf{x}-\mathbf{x^\prime}}{{\lVert \mathbf{x}-\mathbf{x^\prime} \rVert}^3} \text{d}^3\text{x}^\prime $$ which is precisely the electric field created by a charge density of $Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2$. In fact, we also have $$ \langle \rho(\mathbf{x}) \rangle = \langle {\Psi|\hat{\rho}(\mathbf{x})|\Psi} \rangle = Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2. $$ Using these relations, the expecation values of the electric field and charge density can be shown to obey the Ehrenfest Theorem-esque relationship $$ \nabla \cdot \langle \mathbf{E}(\mathbf{x}) \rangle = 4\pi\langle\rho(\mathbf{x})\rangle. $$ This immediately suggests that the error in the question is caused by the simple fact that $$ U = \frac{1}{8\pi} \int \langle {\lVert \mathbf{E}(\mathbf{x}) \rVert}^2 \rangle \text{d}^3\text{x} \neq \frac{1}{8\pi} \int {\lVert \langle \mathbf{E}(\mathbf{x}) \rangle \rVert}^2 \text{d}^3\text{x}, $$ or alternatively \begin{equation} U = \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1})\hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1\text{d}^3\text{x}_2 \tag{2} \end{equation} \begin{equation} \neq \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1}) \rangle \langle \hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2. \tag{3} \end{equation} In the question I have implicitly replaced the average of the product with the product of the averages, neglecting the covariance. Eq.3 was calculated rather than Eq.2. Using Eq.1, Eq.2 can be rewritten as $$ U = \frac{1}{2} \int \frac{ (Ze\delta^3(\mathbf{x_1})-e\delta^3(\mathbf{x_1}-\mathbf{x})) (Ze\delta^3(\mathbf{x_2})-e\delta^3(\mathbf{x_2}-\mathbf{x})) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \lvert \Psi(\mathbf{x}) \rvert^2 \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x} $$ $$ = \frac{(Ze)^2}{2} \int \frac{ \delta^3(\mathbf{x_1})\delta^3(\mathbf{x_2})}{{\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 - Ze^2 \int \frac{ \lvert \Psi(\mathbf{x}) \rvert^2 }{ \lVert \mathbf{x} \rVert } \text{d}^3\text{x} $$ $$ + \hspace{2mm} \frac{e^2}{2} \int \lvert \Psi(\mathbf{x}) \rvert^2 \frac{ \delta^3(\mathbf{x_1}-\mathbf{x})\delta^3(\mathbf{x_2}-\mathbf{x}) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x}. $$ The subsitutions $\mathbf{x_1}^\prime=\mathbf{x_1}-\mathbf{x}$ and $\mathbf{x_2}^\prime=\mathbf{x_2}-\mathbf{x}$ allow the integral over $\mathbf{x}$ to be done independently of the integrals over $\mathbf{x_1}$ and $\mathbf{x_2}$, and the third term simplifies to $$ \frac{e^2}{2} \int \frac{ \delta^3(\mathbf{x_1}^\prime)\delta^3(\mathbf{x_2}^\prime) }{ {\lVert \mathbf{x_1}^\prime-\mathbf{x_2}^\prime \rVert} } \text{d}^3\text{x}^\prime_1 \text{d}^3\text{x}^\prime_2. $$ This expression is exactly the (formally infinite) electrostatic self-energy of a point particle with charge $e$, as expected.
I will treat the electric field semiclassically as a function of the position of the electron, so that $$ \langle \mathbf{x_1}|\hat{\mathbf{E}}(\mathbf{x})|\mathbf{x_2} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} -e\frac{\mathbf{x}-\mathbf{x_1}}{\lVert \mathbf{x}-\mathbf{x_1} \rVert}^3 \right). $$ As before, it is implied that a charge $Ze$ of infinite mass lies at the origin. Similarly for the charge density, \begin{equation} \langle {\mathbf{x_1}|\hat{\rho}(\mathbf{x})|\mathbf{x_2}} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze\delta^{3}(\mathbf{x}) - e\delta^{3}(\mathbf{x}-\mathbf{x_1}) \right). \tag{1} \end{equation} The average electric field of an arbitrary state $\Psi$ is then given by $$ \langle \mathbf{E}(\mathbf{x}) \rangle = \langle {\Psi|\hat{\mathbf{E}}(\mathbf{x})|\Psi} \rangle = Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} + \int -e\lvert \Psi(\mathbf{x^\prime}) \rvert^2 \frac{\mathbf{x}-\mathbf{x^\prime}}{{\lVert \mathbf{x}-\mathbf{x^\prime} \rVert}^3} \text{d}^3\text{x}^\prime $$ which is precisely the electric field created by a charge density of $Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2$. In fact, we also have $$ \langle \rho(\mathbf{x}) \rangle = \langle {\Psi|\hat{\rho}(\mathbf{x})|\Psi} \rangle = Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2. $$ Using these relations, the expecation values of the electric field and charge density can be shown to obey the Ehrenfest Theorem-esque relationship $$ \nabla \cdot \langle \mathbf{E}(\mathbf{x}) \rangle = 4\pi\langle\rho(\mathbf{x})\rangle. $$ This immediately suggests that the error in the question is caused by the simple fact that $$ U = \frac{1}{8\pi} \int \langle {\lVert \mathbf{E}(\mathbf{x}) \rVert}^2 \rangle \text{d}^3\text{x} \neq \frac{1}{8\pi} \int {\lVert \langle \mathbf{E}(\mathbf{x}) \rangle \rVert}^2 \text{d}^3\text{x}, $$ or alternatively \begin{equation} U = \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1})\hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1\text{d}^3\text{x}_2 \tag{2} \end{equation} \begin{equation} \neq \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1}) \rangle \langle \hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2. \tag{3} \end{equation} In the question I have implicitly replaced the average of the product with the product of the averages, neglecting the covariance. Eq.3 was calculated rather than Eq.2. Using Eq.1, Eq.2 can be rewritten as $$ U = \frac{1}{2} \int \frac{ (Ze\delta^3(\mathbf{x_1})-e\delta^3(\mathbf{x_1}-\mathbf{x})) (Ze\delta^3(\mathbf{x_2})-e\delta^3(\mathbf{x_2}-\mathbf{x})) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \lvert \Psi(\mathbf{x}) \rvert^2 \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x} $$ $$ = \frac{(Ze)^2}{2} \int \frac{ \delta^3(\mathbf{x_1})\delta^3(\mathbf{x_2})}{{\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 - Ze^2 \int \frac{ \lvert \Psi(\mathbf{x}) \rvert^2 }{ \lVert \mathbf{x} \rVert } \text{d}^3\text{x} $$ $$ + \hspace{2mm} \frac{e^2}{2} \int \lvert \Psi(\mathbf{x}) \rvert^2 \frac{ \delta^3(\mathbf{x_1}-\mathbf{x})\delta^3(\mathbf{x_2}-\mathbf{x}) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x}. $$ The subsitutions $\mathbf{x_1}^\prime=\mathbf{x_1}-\mathbf{x}$ and $\mathbf{x_2}^\prime=\mathbf{x_2}-\mathbf{x}$ allow the integral over $\mathbf{x}$ to be done independently of the integrals over $\mathbf{x_1}$ and $\mathbf{x_2}$, and the third term simplifies to $$ \frac{e^2}{2} \int \frac{ \delta^3(\mathbf{x_1}^\prime)\delta^3(\mathbf{x_2}^\prime) }{ {\lVert \mathbf{x_1}^\prime-\mathbf{x_2}^\prime \rVert} } \text{d}^3\text{x}^\prime_1 \text{d}^3\text{x}^\prime_2. $$ This expression is exactly the (formally infinite) electrostatic self-energy of a point particle with charge $e$, as expected.
I will treat the electric field semiclassically as a function of the position of the electron, so that $$ \langle \mathbf{x_1}|\hat{\mathbf{E}}(\mathbf{x})|\mathbf{x_2} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} -e\frac{\mathbf{x}-\mathbf{x_1}}{\lVert \mathbf{x}-\mathbf{x_1} \rVert^3} \right). $$ As before, it is implied that a charge $Ze$ of infinite mass lies at the origin. Similarly for the charge density, \begin{equation} \langle {\mathbf{x_1}|\hat{\rho}(\mathbf{x})|\mathbf{x_2}} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze\delta^{3}(\mathbf{x}) - e\delta^{3}(\mathbf{x}-\mathbf{x_1}) \right). \tag{1} \end{equation} The average electric field of an arbitrary state $\Psi$ is then given by $$ \langle \mathbf{E}(\mathbf{x}) \rangle = \langle {\Psi|\hat{\mathbf{E}}(\mathbf{x})|\Psi} \rangle = Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} + \int -e\lvert \Psi(\mathbf{x^\prime}) \rvert^2 \frac{\mathbf{x}-\mathbf{x^\prime}}{{\lVert \mathbf{x}-\mathbf{x^\prime} \rVert}^3} \text{d}^3\text{x}^\prime $$ which is precisely the electric field created by a charge density of $Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2$. In fact, we also have $$ \langle \rho(\mathbf{x}) \rangle = \langle {\Psi|\hat{\rho}(\mathbf{x})|\Psi} \rangle = Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2. $$ Using these relations, the expecation values of the electric field and charge density can be shown to obey the Ehrenfest Theorem-esque relationship $$ \nabla \cdot \langle \mathbf{E}(\mathbf{x}) \rangle = 4\pi\langle\rho(\mathbf{x})\rangle. $$ This immediately suggests that the error in the question is caused by the simple fact that $$ U = \frac{1}{8\pi} \int \langle {\lVert \mathbf{E}(\mathbf{x}) \rVert}^2 \rangle \text{d}^3\text{x} \neq \frac{1}{8\pi} \int {\lVert \langle \mathbf{E}(\mathbf{x}) \rangle \rVert}^2 \text{d}^3\text{x}, $$ or alternatively \begin{equation} U = \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1})\hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1\text{d}^3\text{x}_2 \tag{2} \end{equation} \begin{equation} \neq \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1}) \rangle \langle \hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2. \tag{3} \end{equation} In the question I have implicitly replaced the average of the product with the product of the averages, neglecting the covariance. Eq.3 was calculated rather than Eq.2. Using Eq.1, Eq.2 can be rewritten as $$ U = \frac{1}{2} \int \frac{ (Ze\delta^3(\mathbf{x_1})-e\delta^3(\mathbf{x_1}-\mathbf{x})) (Ze\delta^3(\mathbf{x_2})-e\delta^3(\mathbf{x_2}-\mathbf{x})) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \lvert \Psi(\mathbf{x}) \rvert^2 \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x} $$ $$ = \frac{(Ze)^2}{2} \int \frac{ \delta^3(\mathbf{x_1})\delta^3(\mathbf{x_2})}{{\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 - Ze^2 \int \frac{ \lvert \Psi(\mathbf{x}) \rvert^2 }{ \lVert \mathbf{x} \rVert } \text{d}^3\text{x} $$ $$ + \hspace{2mm} \frac{e^2}{2} \int \lvert \Psi(\mathbf{x}) \rvert^2 \frac{ \delta^3(\mathbf{x_1}-\mathbf{x})\delta^3(\mathbf{x_2}-\mathbf{x}) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x}. $$ The subsitutions $\mathbf{x_1}^\prime=\mathbf{x_1}-\mathbf{x}$ and $\mathbf{x_2}^\prime=\mathbf{x_2}-\mathbf{x}$ allow the integral over $\mathbf{x}$ to be done independently of the integrals over $\mathbf{x_1}$ and $\mathbf{x_2}$, and the third term simplifies to $$ \frac{e^2}{2} \int \frac{ \delta^3(\mathbf{x_1}^\prime)\delta^3(\mathbf{x_2}^\prime) }{ {\lVert \mathbf{x_1}^\prime-\mathbf{x_2}^\prime \rVert} } \text{d}^3\text{x}^\prime_1 \text{d}^3\text{x}^\prime_2. $$ This expression is exactly the (formally infinite) electrostatic self-energy of a point particle with charge $e$, as expected.
I will treat the electric field semiclassically as a function of the position of the electron, so that $$ \langle \mathbf{x_1}|\hat{\mathbf{E}}(\mathbf{x})|\mathbf{x_2} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} -e\frac{\mathbf{x}-\mathbf{x_1}}{\lVert \mathbf{x}-\mathbf{x_1} \rVert}^3 \right). $$ As before, it is implied that a charge $Ze$ of infinite mass lies at the origin. Similarly for the charge density, \begin{equation} \langle {\mathbf{x_1}|\hat{\rho}(\mathbf{x})|\mathbf{x_2}} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze\delta^{3}(\mathbf{x}) - e\delta^{3}(\mathbf{x}-\mathbf{x_1}) \right). \tag{1} \end{equation} The average electric field of an arbitrary state $\Psi$ is then given by $$ \langle \mathbf{E}(\mathbf{x}) \rangle = \langle {\Psi|\hat{\mathbf{E}}(\mathbf{x})|\Psi} \rangle = Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} + \int -e\lvert \Psi(\mathbf{x^\prime}) \rvert^2 \frac{\mathbf{x}-\mathbf{x^\prime}}{{\lVert \mathbf{x}-\mathbf{x^\prime} \rVert}^3} \text{d}^3\text{x}^\prime $$ which is precisely the electric field created by a charge density of $Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2$. In fact, we also have $$ \langle \rho(\mathbf{x}) \rangle = \langle {\Psi|\hat{\rho}(\mathbf{x})|\Psi} \rangle = Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2. $$ Using these relations, the expecation values of the electric field and charge density can be shown to obey the Ehrenfest Theorem-esque relationship $$ \nabla \cdot \langle \mathbf{E}(\mathbf{x}) \rangle = 4\pi\langle\rho(\mathbf{x})\rangle. $$ This immediately suggests that the error in the question is caused by the simple fact that $$ U = \frac{1}{8\pi} \int \langle {\lVert \mathbf{E}(\mathbf{x}) \rVert}^2 \rangle \text{d}^3\text{x} \neq \frac{1}{8\pi} \int {\lVert \langle \mathbf{E}(\mathbf{x}) \rangle \rVert}^2 \text{d}^3\text{x}, $$ or alternatively \begin{equation} U = \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1})\hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1\text{d}^3\text{x}_2\text{d}^3\text{x} \tag{2} \end{equation}\begin{equation} U = \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1})\hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1\text{d}^3\text{x}_2 \tag{2} \end{equation} \begin{equation} \neq \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1}) \rangle \langle \hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x}. \tag{3} \end{equation}\begin{equation} \neq \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1}) \rangle \langle \hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2. \tag{3} \end{equation} In the question I have implicitly replaced the average of the product with the product of the averages, neglecting the covariance. Eq.3 was calculated rather than Eq.2. Using Eq.1, Eq.2 can be rewritten as $$ U = \frac{1}{2} \int \frac{ (Ze\delta^3(\mathbf{x_1})-e\delta^3(\mathbf{x_1}-\mathbf{x})) (Ze\delta^3(\mathbf{x_2})-e\delta^3(\mathbf{x_2}-\mathbf{x})) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \lvert \Psi(\mathbf{x}) \rvert^2 \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x} $$ $$ = \frac{(Ze)^2}{2} \int \frac{ \delta^3(\mathbf{x_1})\delta^3(\mathbf{x_2})}{{\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 - Ze^2 \int \frac{ \lvert \Psi(\mathbf{x}) \rvert^2 }{ \lVert \mathbf{x} \rVert } \text{d}^3\text{x} $$ $$ + \hspace{2mm} \frac{e^2}{2} \int \lvert \Psi(\mathbf{x}) \rvert^2 \frac{ \delta^3(\mathbf{x_1}-\mathbf{x})\delta^3(\mathbf{x_2}-\mathbf{x}) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x}. $$ The subsitutions $\mathbf{x_1}^\prime=\mathbf{x_1}-\mathbf{x}$ and $\mathbf{x_2}^\prime=\mathbf{x_2}-\mathbf{x}$ allow the integral over $\mathbf{x}$ to be done independently of the integrals over $\mathbf{x_1}$ and $\mathbf{x_2}$, and the third term simplifies to $$ \frac{e^2}{2} \int \frac{ \delta^3(\mathbf{x_1}^\prime)\delta^3(\mathbf{x_2}^\prime) }{ {\lVert \mathbf{x_1}^\prime-\mathbf{x_2}^\prime \rVert} } \text{d}^3\text{x}^\prime_1 \text{d}^3\text{x}^\prime_2. $$ This expression is exactly the (formally infinite) electrostatic self-energy of a point particle with charge $e$, as expected.
I will treat the electric field semiclassically as a function of the position of the electron, so that $$ \langle \mathbf{x_1}|\hat{\mathbf{E}}(\mathbf{x})|\mathbf{x_2} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} -e\frac{\mathbf{x}-\mathbf{x_1}}{\lVert \mathbf{x}-\mathbf{x_1} \rVert}^3 \right). $$ As before, it is implied that a charge $Ze$ of infinite mass lies at the origin. Similarly for the charge density, \begin{equation} \langle {\mathbf{x_1}|\hat{\rho}(\mathbf{x})|\mathbf{x_2}} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze\delta^{3}(\mathbf{x}) - e\delta^{3}(\mathbf{x}-\mathbf{x_1}) \right). \tag{1} \end{equation} The average electric field of an arbitrary state $\Psi$ is then given by $$ \langle \mathbf{E}(\mathbf{x}) \rangle = \langle {\Psi|\hat{\mathbf{E}}(\mathbf{x})|\Psi} \rangle = Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} + \int -e\lvert \Psi(\mathbf{x^\prime}) \rvert^2 \frac{\mathbf{x}-\mathbf{x^\prime}}{{\lVert \mathbf{x}-\mathbf{x^\prime} \rVert}^3} \text{d}^3\text{x}^\prime $$ which is precisely the electric field created by a charge density of $Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2$. In fact, we also have $$ \langle \rho(\mathbf{x}) \rangle = \langle {\Psi|\hat{\rho}(\mathbf{x})|\Psi} \rangle = Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2. $$ Using these relations, the expecation values of the electric field and charge density can be shown to obey the Ehrenfest Theorem-esque relationship $$ \nabla \cdot \langle \mathbf{E}(\mathbf{x}) \rangle = 4\pi\langle\rho(\mathbf{x})\rangle. $$ This immediately suggests that the error in the question is caused by the simple fact that $$ U = \frac{1}{8\pi} \int \langle {\lVert \mathbf{E}(\mathbf{x}) \rVert}^2 \rangle \text{d}^3\text{x} \neq \frac{1}{8\pi} \int {\lVert \langle \mathbf{E}(\mathbf{x}) \rangle \rVert}^2 \text{d}^3\text{x}, $$ or alternatively \begin{equation} U = \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1})\hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1\text{d}^3\text{x}_2\text{d}^3\text{x} \tag{2} \end{equation} \begin{equation} \neq \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1}) \rangle \langle \hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x}. \tag{3} \end{equation} In the question I have implicitly replaced the average of the product with the product of the averages, neglecting the covariance. Eq.3 was calculated rather than Eq.2. Using Eq.1, Eq.2 can be rewritten as $$ U = \frac{1}{2} \int \frac{ (Ze\delta^3(\mathbf{x_1})-e\delta^3(\mathbf{x_1}-\mathbf{x})) (Ze\delta^3(\mathbf{x_2})-e\delta^3(\mathbf{x_2}-\mathbf{x})) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \lvert \Psi(\mathbf{x}) \rvert^2 \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x} $$ $$ = \frac{(Ze)^2}{2} \int \frac{ \delta^3(\mathbf{x_1})\delta^3(\mathbf{x_2})}{{\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 - Ze^2 \int \frac{ \lvert \Psi(\mathbf{x}) \rvert^2 }{ \lVert \mathbf{x} \rVert } \text{d}^3\text{x} $$ $$ + \hspace{2mm} \frac{e^2}{2} \int \lvert \Psi(\mathbf{x}) \rvert^2 \frac{ \delta^3(\mathbf{x_1}-\mathbf{x})\delta^3(\mathbf{x_2}-\mathbf{x}) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x}. $$ The subsitutions $\mathbf{x_1}^\prime=\mathbf{x_1}-\mathbf{x}$ and $\mathbf{x_2}^\prime=\mathbf{x_2}-\mathbf{x}$ allow the integral over $\mathbf{x}$ to be done independently of the integrals over $\mathbf{x_1}$ and $\mathbf{x_2}$, and the third term simplifies to $$ \frac{e^2}{2} \int \frac{ \delta^3(\mathbf{x_1}^\prime)\delta^3(\mathbf{x_2}^\prime) }{ {\lVert \mathbf{x_1}^\prime-\mathbf{x_2}^\prime \rVert} } \text{d}^3\text{x}^\prime_1 \text{d}^3\text{x}^\prime_2. $$ This expression is exactly the (formally infinite) electrostatic self-energy of a point particle with charge $e$, as expected.
I will treat the electric field semiclassically as a function of the position of the electron, so that $$ \langle \mathbf{x_1}|\hat{\mathbf{E}}(\mathbf{x})|\mathbf{x_2} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} -e\frac{\mathbf{x}-\mathbf{x_1}}{\lVert \mathbf{x}-\mathbf{x_1} \rVert}^3 \right). $$ As before, it is implied that a charge $Ze$ of infinite mass lies at the origin. Similarly for the charge density, \begin{equation} \langle {\mathbf{x_1}|\hat{\rho}(\mathbf{x})|\mathbf{x_2}} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze\delta^{3}(\mathbf{x}) - e\delta^{3}(\mathbf{x}-\mathbf{x_1}) \right). \tag{1} \end{equation} The average electric field of an arbitrary state $\Psi$ is then given by $$ \langle \mathbf{E}(\mathbf{x}) \rangle = \langle {\Psi|\hat{\mathbf{E}}(\mathbf{x})|\Psi} \rangle = Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} + \int -e\lvert \Psi(\mathbf{x^\prime}) \rvert^2 \frac{\mathbf{x}-\mathbf{x^\prime}}{{\lVert \mathbf{x}-\mathbf{x^\prime} \rVert}^3} \text{d}^3\text{x}^\prime $$ which is precisely the electric field created by a charge density of $Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2$. In fact, we also have $$ \langle \rho(\mathbf{x}) \rangle = \langle {\Psi|\hat{\rho}(\mathbf{x})|\Psi} \rangle = Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2. $$ Using these relations, the expecation values of the electric field and charge density can be shown to obey the Ehrenfest Theorem-esque relationship $$ \nabla \cdot \langle \mathbf{E}(\mathbf{x}) \rangle = 4\pi\langle\rho(\mathbf{x})\rangle. $$ This immediately suggests that the error in the question is caused by the simple fact that $$ U = \frac{1}{8\pi} \int \langle {\lVert \mathbf{E}(\mathbf{x}) \rVert}^2 \rangle \text{d}^3\text{x} \neq \frac{1}{8\pi} \int {\lVert \langle \mathbf{E}(\mathbf{x}) \rangle \rVert}^2 \text{d}^3\text{x}, $$ or alternatively \begin{equation} U = \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1})\hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1\text{d}^3\text{x}_2 \tag{2} \end{equation} \begin{equation} \neq \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1}) \rangle \langle \hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2. \tag{3} \end{equation} In the question I have implicitly replaced the average of the product with the product of the averages, neglecting the covariance. Eq.3 was calculated rather than Eq.2. Using Eq.1, Eq.2 can be rewritten as $$ U = \frac{1}{2} \int \frac{ (Ze\delta^3(\mathbf{x_1})-e\delta^3(\mathbf{x_1}-\mathbf{x})) (Ze\delta^3(\mathbf{x_2})-e\delta^3(\mathbf{x_2}-\mathbf{x})) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \lvert \Psi(\mathbf{x}) \rvert^2 \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x} $$ $$ = \frac{(Ze)^2}{2} \int \frac{ \delta^3(\mathbf{x_1})\delta^3(\mathbf{x_2})}{{\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 - Ze^2 \int \frac{ \lvert \Psi(\mathbf{x}) \rvert^2 }{ \lVert \mathbf{x} \rVert } \text{d}^3\text{x} $$ $$ + \hspace{2mm} \frac{e^2}{2} \int \lvert \Psi(\mathbf{x}) \rvert^2 \frac{ \delta^3(\mathbf{x_1}-\mathbf{x})\delta^3(\mathbf{x_2}-\mathbf{x}) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x}. $$ The subsitutions $\mathbf{x_1}^\prime=\mathbf{x_1}-\mathbf{x}$ and $\mathbf{x_2}^\prime=\mathbf{x_2}-\mathbf{x}$ allow the integral over $\mathbf{x}$ to be done independently of the integrals over $\mathbf{x_1}$ and $\mathbf{x_2}$, and the third term simplifies to $$ \frac{e^2}{2} \int \frac{ \delta^3(\mathbf{x_1}^\prime)\delta^3(\mathbf{x_2}^\prime) }{ {\lVert \mathbf{x_1}^\prime-\mathbf{x_2}^\prime \rVert} } \text{d}^3\text{x}^\prime_1 \text{d}^3\text{x}^\prime_2. $$ This expression is exactly the (formally infinite) electrostatic self-energy of a point particle with charge $e$, as expected.
I will treat the electric field semiclassically as a function of the position of the electron, so that $$ \langle \mathbf{x_1}|\hat{\mathbf{E}}(\mathbf{x})|\mathbf{x_2} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} -e\frac{\mathbf{x}-\mathbf{x_1}}{\lVert \mathbf{x}-\mathbf{x_1} \rVert}^3 \right). $$ As before, it is implied that a charge $Ze$ of infinite mass lies at the origin. Similarly for the charge density, \begin{equation} \langle {\mathbf{x_1}|\hat{\rho}(\mathbf{x})|\mathbf{x_2}} \rangle = \delta^{3}(\mathbf{x_1}-\mathbf{x_2}) \left( Ze\delta^{3}(\mathbf{x}) - e\delta^{3}(\mathbf{x}-\mathbf{x_1}) \right). \tag{1} \end{equation} The average electric field of an arbitrary state $\Psi$ is then given by $$ \langle \mathbf{E}(\mathbf{x}) \rangle = \langle {\Psi|\hat{\mathbf{E}}(\mathbf{x})|\Psi} \rangle = Ze \frac{\mathbf{x}}{\lVert \mathbf{x} \rVert^3} + \int -e\lvert \Psi(\mathbf{x^\prime}) \rvert^2 \frac{\mathbf{x}-\mathbf{x^\prime}}{{\lVert \mathbf{x}-\mathbf{x^\prime} \rVert}^3} \text{d}^3\text{x}^\prime $$ which is precisely the electric field created by a charge density of $Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2$. In fact, we also have $$ \langle \rho(\mathbf{x}) \rangle = \langle {\Psi|\hat{\rho}(\mathbf{x})|\Psi} \rangle = Ze\delta^3(\mathbf{x})-e\lvert \Psi(\mathbf{x}) \rvert^2. $$ Using these relations, the expecation values of the electric field and charge density can be shown to obey the Ehrenfest Theorem-esque relationship $$ \nabla \cdot \langle \mathbf{E}(\mathbf{x}) \rangle = 4\pi\langle\rho(\mathbf{x})\rangle. $$ This immediately suggests that the error in the question is caused by the simple fact that $$ U = \frac{1}{8\pi} \int \langle {\lVert \mathbf{E}(\mathbf{x}) \rVert}^2 \rangle \text{d}^3\text{x} \neq \frac{1}{8\pi} \int {\lVert \langle \mathbf{E}(\mathbf{x}) \rangle \rVert}^2 \text{d}^3\text{x}, $$ or alternatively \begin{equation} U = \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1})\hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1\text{d}^3\text{x}_2\text{d}^3\text{x} \tag{2} \end{equation} \begin{equation} \neq \frac{1}{2} \int \frac{ \langle \hat{\rho}(\mathbf{x_1}) \rangle \langle \hat{\rho}(\mathbf{x_2}) \rangle }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x}. \tag{3} \end{equation} In the question I have implicitly replaced the average of the product with the product of the averages, neglecting the covariance. Eq.3 was calculated rather than Eq.2. Using Eq.1, Eq.2 can be rewritten as $$ U = \frac{1}{2} \int \frac{ (Ze\delta^3(\mathbf{x_1})-e\delta^3(\mathbf{x_1}-\mathbf{x})) (Ze\delta^3(\mathbf{x_2})-e\delta^3(\mathbf{x_2}-\mathbf{x})) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \lvert \Psi(\mathbf{x}) \rvert^2 \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x} $$ $$ = \frac{(Ze)^2}{2} \int \frac{ \delta^3(\mathbf{x_1})\delta^3(\mathbf{x_2})}{{\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 - Ze^2 \int \frac{ \lvert \Psi(\mathbf{x}) \rvert^2 }{ \lVert \mathbf{x} \rVert } \text{d}^3\text{x} $$ $$ + \hspace{2mm} \frac{e^2}{2} \int \lvert \Psi(\mathbf{x}) \rvert^2 \frac{ \delta^3(\mathbf{x_1}-\mathbf{x})\delta^3(\mathbf{x_2}-\mathbf{x}) }{ {\lVert \mathbf{x_1}-\mathbf{x_2} \rVert} } \text{d}^3\text{x}_1 \text{d}^3\text{x}_2 \text{d}^3\text{x}. $$ The subsitutions $\mathbf{x_1}^\prime=\mathbf{x_1}-\mathbf{x}$ and $\mathbf{x_2}^\prime=\mathbf{x_2}-\mathbf{x}$ allow the integral over $\mathbf{x}$ to be done independently of the integrals over $\mathbf{x_1}$ and $\mathbf{x_2}$, and the third term simplifies to $$ \frac{e^2}{2} \int \frac{ \delta^3(\mathbf{x_1}^\prime)\delta^3(\mathbf{x_2}^\prime) }{ {\lVert \mathbf{x_1}^\prime-\mathbf{x_2}^\prime \rVert} } \text{d}^3\text{x}^\prime_1 \text{d}^3\text{x}^\prime_2. $$ This expression is exactly the (formally infinite) electrostatic self-energy of a point particle with charge $e$, as expected.