Electrostatic interaction energy I have a confusion about in which part of a Hamiltonian an electrostatic interaction energy can be found. Consider a system of electrons interacting through electric field $\mathbf{E}$ with the Hamiltonian
$$
H=H_\mathrm{el}+H_\mathrm{field},
$$
where the Hamiltonians of electrons and of the field are
$$
H_\mathrm{el}=\sum_i\frac{\mathbf{p}_i^2}{2m}+\sum_ie\varphi(\mathbf{r}_i),\qquad H_\mathrm{field}=\int\frac{\mathbf{E}^2}{8\pi}dV
$$
(I omit magnetic field and vector potential). Using $\mathbf{E}=-\nabla\varphi$ and $\mathrm{div}\,\mathbf{E}=4\pi\rho$, we can integrate $H_\mathrm{field}$ by parts
$$
H_\mathrm{field}=\int\frac{(\nabla\varphi)^2}{8\pi}dV=\oint\frac{\varphi\nabla\varphi}{8\pi}\cdot d\mathbf{S}-\int\frac{\varphi\nabla^2\varphi}{8\pi}dV=\frac12\int\rho\varphi\:dV,
$$
neglecting, as usual, the surface integral.
Now, assuming each electron creates its own electrostatic potential $\varphi_i$, so that $\varphi=\sum_i\varphi_i$, and writing the charge density as $\rho(\mathbf{r})=\sum_je\delta(\mathbf{r}-\mathbf{r}_j)$, we arrive at
$$
H_\mathrm{el}=\sum_i\frac{\mathbf{p}_i^2}{2m}+\sum_{ij}e\varphi_i(\mathbf{r}_j),\qquad H_\mathrm{field}=\frac12\sum_{ij}e\varphi_i(\mathbf{r}_j).
$$
As we see, both $H_\mathrm{el}$ and $H_\mathrm{field}$ contain electron self-energies $\frac12e\varphi_i(\mathbf{r}_i)$ (which are divergent and are usually neglected) and pairwise Coulomb interaction energies $e\varphi_i(\mathbf{r}_j)=e^2/|\mathbf{r}_i-\mathbf{r}_j|$. 
But why these energies occur in the Hamiltonian three times $-$ once in $H_\mathrm{field}$ and twice in $H_\mathrm{el}$? Is something incorrect in my derivations or interpretations?
UPDATE
As said in the answers, the correct Hamiltonian indeed seems to be
$$
H_\mathrm{correct}=\sum_i\frac{\mathbf{p}_i^2}{2m}+\frac12\sum_{ij}e\varphi_i(\mathbf{r}_j).
$$
Still I don't understand why my Hamiltonian $H=H_\mathrm{el}+H_\mathrm{field}$ is incorrect.
Let's start from the relativistic action for particles and electromagnetic field (I think everybody agree it is correct):
$$
S=-mc\sum_i\int ds_i-\frac{e}{c}\sum_i\int A_\mu(x_i) dx^\mu_i-\frac1{16\pi c}\int F_{\mu\nu}F^{\mu\nu}cdtdV.
$$
In the non-relativistic limit, neglecting magnetic field and vector potential, we have:
$$
-mc\sum_i\int ds_i=\sum_i\int\frac{m\mathbf{v}_i^2}{2}dt,
$$
$$
-\frac{e}{c}\sum_i\int A_\mu(x_i) dx^\mu_i=-\sum_i\int e\varphi(\mathbf{r}_i)dt,
$$
$$
-\frac1{16\pi c}\int F_{\mu\nu}F^{\mu\nu}cdtdV=\int\frac{\mathbf{E}^2}{8\pi}dVdt.
$$
Writing the action as a time integrated Lagrange function $S=\int Ldt$, we get
$$
L=\sum_i\frac{m\mathbf{v}_i^2}{2}-\sum_ie\varphi(\mathbf{r}_i)+\int\frac{\mathbf{E}^2}{8\pi}dV.
$$
Now we can identify the particle momenta $\mathbf{p}_i=\partial L/\partial\mathbf{v}_i=m\mathbf{v}_i$ and scalar potential "momenta" $\Pi_k=\partial L/\partial(\partial_k\varphi)=\partial_k\varphi/4\pi$ ($k=x,y,z$) associated with its "velocities" $\partial_k\varphi$. Finally, the Hamiltonian is
$$
H=\sum_i\mathbf{p}_i\mathbf{v}_i+\sum_k\Pi_k\partial_k\varphi-L=
$$
$$
=\sum_i\frac{\mathbf{p}_i^2}{2m}+\sum_ie\varphi(\mathbf{r}_i)+\int\frac{\mathbf{E}^2}{8\pi}dV.
$$
We see the Hamiltonian $H$ can be derived from the first principles, so why it is wrong?
 A: There are two problems with your derivation of the Hamiltonian
$$
\sum_i\frac{\mathbf{p}_i^2}{2m}+\sum_ie\varphi(\mathbf{r}_i)+\int\frac{\mathbf{E}^2}{8\pi}dV,
$$
one in the final step, the other one in the basic assumptions.
The first problem is in the last step where you do Legendre's transformation to get $H$ from $L$. It should be the case that the term $\int\frac{\mathbf{E}^2}{8\pi}dV$ should appear with negative sign in the resulting Hamiltonian so the terms $1/2e\varphi_j(\mathbf r_i)$ inside it get subtracted from terms $e\varphi_j(\mathbf r_i)$ originating in the particle term and the resulting Hamiltonian contains only terms $1/2e\varphi_j(\mathbf r_i)$.
I think you got the sign wrong because you added the terms $\Pi_k \partial_k\varphi$ which shouldn't be there.
EDIT
True, the field has its degrees of freedom. But when passing from Lorentz-covariant Lagrangian description that is treating all spatial coordinates in the same way to customary Hamiltonian description that is treating $t$ in a special way, one must first alter the Lagrangian description to a different one, where there are only 4 Lagr. field velocities (derivatives of 4 field coordinates with respect to time), as opposed to 16 velocities in the original description (all possible derivatives of 4 field coordinates). Then, $L$ depend only on field velocities $\partial_t A_1,\partial_t A_2,\partial_t A_3$, but it does not depend on $\partial_t\varphi$. So if we ignore the magnetic field terms, there is no contribution due to field "pv" terms in the definition
$$
H=\sum_i\mathbf{p}_i\mathbf{v}_i + \int \sum_\mu \frac{\partial \mathscr{L} }{\partial (\partial_t A_\mu)}\partial_t A_\mu \,d^3\mathbf x - L
$$
we just have
$$
H=\sum_i\mathbf{p}_i\mathbf{v}_i - L.
$$
This expression leads to the correct result
$$
H_\mathrm{correct}=\sum_i\frac{\mathbf{p}_i^2}{2m}+\frac12\sum_{ij}e\varphi_i(\mathbf{r}_j).
$$
provided the correct Lagrangian is used (see the one in the links I gave).
END OF EDIT
The other problem is that (although the result is correct) this whole derivation is still invalid because the main assumption - the action integral
$$
S=-mc\sum_i\int ds_i-\frac{e}{c}\sum_i\int A_\mu(x_i) dx^\mu_i-\frac1{16\pi c}\int F_{\mu\nu}F^{\mu\nu}cdtdV.
$$
and the implicated Lagrangian function
$$
L = \sum_i -mc^2\sqrt{1-v_i^2/c^2} - \int \sum_{i}  j_i^\mu/c\, A_\mu(x_i^\mu) dV - \int \frac{1}{16\pi}F_{\mu\nu} F^{\mu\nu} dV
$$
are mathematically invalid. One indication of this is the fact that the interaction integral is not mathematically defined and the field integral is infinite. Another indication is that this Lagrangian invariably lead people to self-interaction models (Lorentz-Abraham-Dirac equation and others) that nobody was able to make consistent.
All these problems stem from the mistaken assumption that one can interpret Poynting's formulae as field energy quantities even if particles are points. This is not possible, because point particles have field that diverges in the points where particles are and so the usual derivation of the work-energy theorem for EM field and matter based on manipulating the expression
$$
\int \mathbf E \cdot \mathbf j \,dV
$$
fails.
Luckily, there is a consistent theory of charged point particles, due to Frenkel (I think he was the first to formulate the problem and the solution in a general way). See for example
J. Frenkel, Zur Elektrodynamik punktfoermiger Elektronen, Zeits. f. Phys., 32, (1925), p. 518-534. http://dx.doi.org/10.1007/BF01331692
For a shorter, more easy-to-read account, see also
R. C. Stabler, A Possible Modification of Classical Electrodynamics, Physics Letters, 8, 3, (1964), p. 185-187. http://dx.doi.org/10.1016/S0031-9163(64)91989-4
I gave some details on this approach in these answers:
Complete classical description of two interacting charges
About Self-Energy/Self-Potential Energy
A: I would say that you are counting twice since the beginning. When you write:
$$H=H_{el}+H_{field}$$
your writing of $H_{el}$ suppose that $H_{field}$ is due to an external field. But you are including the potential of the electrons in it, potential you already included in your electron hamiltonian:
$$\sum_{i,j} e \varphi_i(r_j) $$ 
So I think the idea is to write the hamiltonian as the kinetic part of your electron, and an interaction part that include the interaction of your charge with all the potential you are considering (external field and other electrons).
I don't know anything like electron self energy (your electron doesn't interact with himself or with it's own electric field), your sum should be for $i $ different than $j$, as you are writing the potential of one electron in the field of the other electrons. 
So the proper formula is:
$$ H = H_{kin} + H_{int}$$
$$ H_{int} = \frac{1}{2}\sum_{i,j,i\neq j}e \varphi_i(r_j) + \sum_i eV_{ext}(r_i)$$
A: I'm thinking along similar lines as denis. The Hamiltonian for the system is simply composed of a kinetic part and a potential part, so
$$H_{el}=H_{kin}+H_{pot}$$For $H_{kin}$ we can write
$$H_{kin}=\sum_i\frac{\mathbf{p}_i}{2m}$$
 and for $H_{pot}$
$$ H_{pot} = \frac 1 2\sum_{i,j,i\neq j}e \varphi_i(r_j)$$
The factor $\frac 1 2$ appears because every pair electrons is taken twice in the summation

.
So the total Hamiltonian $H_{tot}$ becomes the sum of those two:
$$H_{tot}=H_{kin}=\sum_i\frac{\mathbf{p}_i}{2m}+ \frac 1 2\sum_{i,j,i\neq j}e \varphi_i(r_j)$$
$H_{field}$ is already contained in this expression (unless there was an external field penetrating the collection of electrons).
