The classical electromagnetic field Lagrangian density can be written as $$-\frac{1}{4} F^{\mu \nu} F_{\mu \nu} - \frac{1}{c} j^\mu A_{\mu} = \frac{1}{2}\vec{E}(\vec{x})^{2} -\frac{1}{2}\vec{B}(\vec{x})^{2} - \rho (\vec{x}) \phi (\vec{x}) + \vec{\jmath}(\vec{x})\cdot\vec{A}(\vec{x}).$$ However, the electric field is subject to the constraint $\vec{\nabla} \cdot \vec{E} = \rho$, and we can write $\vec{E}= - \partial_t \vec{A} - \vec{\nabla}\phi$, essentially splitting it into its transverse (divergence-free) and longitudinal (rotation-free) components. The integral over the field energy of the electric field then yields (because the fields are orthogonal in each Fourier mode): \begin{align} \int d^3x \, \left( \vec{E}_{t}\!{}^{2} + \vec{E}_{l}{}^{2} \right) &= \int d^3x \, \left[ \vec{E}_{t}\!{}^{2} + (\vec{\nabla} \phi )^2 \right] \\ &= \int d^3 x \, \left(\vec{E}_{t}\!{}^{2}- \phi \vec{\nabla}{}^{2} \phi \right) = \int d^3 x \, \left(\vec{E}_{t}\!{}^{2} - \phi \rho \right). \end{align} I have used one of Green’s identities here, and afterwards applied Gauss' law (in the Coulomb gauge).
So the field energy contains (after one applies the constraint (Gauss' law, which is one of the Euler-Lagrange equations)) already the interaction with a charge distribution. What's going on here? Does that mean the other term is counted too much, and is not necessary? Or should I simply not worry, because the Euler-Lagrange equations do follow from this Lagrangian, and can't be applied beforehand? If so, are there any implications connected with this?
If one does the Legendre transform, and writes down the Hamiltonian, one still has 2 terms with the interaction energy, after applying the equations of motion.