Electron Self-Energy Suppose an electron is bound to a nucleus, leading to an overall charge density
$$
\rho(x) = Ze\delta(x) - e \left|\psi(x) \right|^2.
$$
If I interpret this charge density just as I would in the context of classical electrostatics, I could easily calculate the electric field arising from it. I routinely think about such an electric field, especially in chemistry, where I often informally consider the electron orbitals to be interacting with each other via Coulomb repulsion. This electric field contains an energy given by
$$
U = \frac{1}{8\pi} \int \left| \vec{E(x)} \right|^2 d^3x.
$$
Using Gauss' Law, this energy can be rewritten in terms of the charge density that gives rise to it as
$$
U = \frac{1}{2}\int \frac{\rho(x_1)\rho(x_2)}{|x_{1}-x_{2}|} d^3x_1 d^3x_2.
$$
Substituting my expression for $\rho(x)$ gives three terms:
$$
U = \frac{(Ze)^2}{2} \int \frac{\delta(x_1)\delta(x_2)}{|x_{1}-x_{2}|} d^3x_1 d^3x_2 \\ 
  - Ze^2 \int \frac{|\psi(x)|^2}{|x|} d^3x \\
  + \frac{e^2}{2} \int  \frac{|\psi(x_1)|^2|\psi(x_2)|^2}{|x_{1}-x_{2}|} d^3x_1 d^3x_2.
$$
The first and second terms are totally understandable. The first is the electrostatic self-energy of the nucleus. If one remembers that the charge density of a nucleus is not actually a delta function, but rather some distribution with a length scale on the order of the nuclear radius, it can be understood to have a finite value that gives an upper bound on, for example, the energy released in some process like nuclear fission. Because the nuclear radius typically doesn't change in processes I am interested in, I regard this as a constant and ignore it.
The second term just gives the Coulomb attraction between the nucleus and electron. This is the term routinely used in the Hamiltonian to calculate e.g. the energy of a hydrogen-like atom.
The third term I have little understanding of. The only length scale in the integrand is the Bohr radius $a$, and so the third term will be on the order of $\frac{e^2}{a}$, which is on the order of the binding energy of the nucleus and electron. The fine-structure constant doesn't appear, and so this cannot be regarded as a self-energy interaction stemming from some QED correction (I am treating the electric field classically anyway, so that's expected).
So what is it? Does it have any physical significance, and why is it ignored when one usually calculates the energy of a hydrogen-like atom? This term should increase the atomic radius and cause the electronic charge density to 'repel itself'. I think my result is incorrect and the third term should be absent, but I am not sure what assumption is wrong.
 A: The third term is the self-energy of the electron charge distribution, aka the classical work needed to assemble the electron. Of course the electron is not actually a smeared charge around the nucleus, but that’s a problem with the semi-classical analysis.
A: 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.
A: This question is related to the initial Schrödinger's interpretation of $|\psi|^2$ as charge density (up to a constant factor). Taking into account the self-energy seemed to yield incorrect results for the energies of the atomic levels.
A.O. Barut wrote a lot about this problem (see, e.g., Ann. Physik Leipzig 45 (1988) 1, 31-36 and references there).
Abstract:

Radiation theory is reformulated by reviving the interpretation of
$\psi$-waves as realistic material waves of electronic matter and the
relation of this formulation to second quantized quantum field theory
is discussed.

I did not check his derivations in detail, but he seems to say that taking into account the self-field gives correct results after some renormalization:

The contribution of the self energy can then be calculated by an
iterative procedure. Infinite terms separate and can be absorbed into
energy and mass terms on the left hand side of eq. (10). A more
careful analysis [14] shows that the infrared and ultra-violet
divergences actually do not appear if the intergrals are properly
defined by their Mellin transformations in the same way as the
infinite terms in classical electrodynamics disappear by the use of
Riesz potential [15]. In these last treatments no renormalization
seems to be necessary, the parameters of the theory $e$ and $m$, are
already the physical charge and mass.

Let me also note that I recently proposed a plasma-like modification of Schrödinger's interpretation (Entropy 2022, 24(2), 261 ; Quantum Rep. 2022, 4(4), 486-508)
