The (electrostatic) force on an extended object It is well known that, if I have a system of $N$ particles acted upon only by conservative internal and external forces, then I can obtain the force on the $\mathrm{i^{th}}$ particle as
$$\textbf{F}_i = -\nabla_iU$$
where the potential energy of the entire system $U= U_{int} + U_{ext}$ is the sum of internal (due to inter-particle interactions) and external (due to externally impressed forces) potential energies.
Now suppose I want the net force $\textbf{F}$ on the system of particles. My first question is: I suppose this is $\textit{defined}$ as the sum of forces from external sources on the individual particles? If so, then I should be able to proceed, assuming that the internal forces obey Newton's third law, with
$$\sum_{i}-\nabla_i U=\sum_{i}-\nabla_i U_{ext} = \sum_{i}\textbf{F}_{ext} \equiv \textbf{F},$$
where in the first equality I've observed that internal forces cancel pairwise given the obeyance of Newton's third law.
Now this entire setup has been for the following question, where my goal is to obtain the Coulomb force as a theorem from an assumed form for the Coulomb potential energy.
Suppose that a finite charge distribution $\rho$ is of the form $\rho(\textbf{r}) = \rho_2(\textbf{r}-\textbf{R})$ where $\textbf{R}$ is a parameter "locating" the body and, per this answer, would usefully be placed at the centre of mass of the distribution if we were interested in the dynamics of the system. Plunging forward, we ask about the force on $\rho$ due to an external electrostatic potential $\varphi_1$, where we suppose that the potential energy $U$ of the system is given by (an integral with no bounds is assumed to be "over all space")
$$U \equiv \int d^3r \ \rho(\textbf{r})\varphi_1(\textbf{r})=\int d^3r \ \rho_2(\textbf{r}-\textbf{R})\varphi_1(\textbf{r}).$$
Now I want to prove the following theorem for this particular electrostatic system: $\textbf{F} = -\nabla_\textbf{R}U$ where, importantly, $\textbf{F}$ takes the typical Coulomb form. Is the following correct?
$$-\nabla_\textbf{R}U = -\nabla_\textbf{R}\int d^3r \ \rho_2(\textbf{r}-\textbf{R})\varphi_1(\textbf{r})=\int d^3r \ (-\nabla_\textbf{R})\rho_2(\textbf{r}-\textbf{R})\varphi_1(\textbf{r}) \stackrel{(1)}{=} \int d^3r \ \nabla_\textbf{r}\rho_2(\textbf{r}-\textbf{R})\varphi_1(\textbf{r}) \stackrel{(2)}{=} \\-\int d^3r \ \rho_2(\textbf{r}-\textbf{R})\nabla_\textbf{r}\varphi_1(\textbf{r} )\stackrel{(3)}{=} \int d^3r \ \rho_2(\textbf{r}-\textbf{R})\textbf{E}_1(\textbf{r}) = \int d^3r \ \rho(\textbf{r})\textbf{E}_1(\textbf{r}) \stackrel{(4)}{=} \textbf{F}$$
where in (1) I've changed the variables with respect to which the gradient is taken, in (2) I've used an integration by parts and that the integrand vanishes on a boundary at infinity (since the charge distribution is finite), in (3) I've used the definition of the electric field as the gradient of its potential, and in (4) I've used the identification of $\textbf{F}$ as the sum of external forces (in the continuum limit) as I asked earlier in the question.
Now it seems like what I've done intimately requires that the charge distribution be finite (in particular, in the aforementioned integration by parts step), whereas the second equation in this post seems to suggest that my result should hold in the framework of classical mechanics no matter the nature of $\rho$. Have I made a mistake somewhere then in either analysis?
 A: First off a mathematical loophole is that typically, $\phi_1=O(r^{-1})$. In that case, you just need $\rho_2=O(r^{-k})$ with $k>1$ for the boundary term to vanish. A compact support is not necessary.
The general answer is the conservation of charge. When it is not conserved in the bulk, you’ll need to add a boundary term to the potential taking into account the additional energy due to charge accumulation at the boundary.
Things are clearer in a compact domain, but can be extrapolated to an infinite one. I’ll note $\Omega$ the domain of interest and I’ll start by the force to derive the potential via the principle of virtual work.
Using a general displacement, $\delta s$ possibly depending on $r$, I have:
$$
\begin{align}
\delta W &= \int_\Omega d^3r \rho_2 E_1 \cdot \delta s \\
&= -\int_\Omega d^3r \rho_2 \nabla \phi_1 \cdot \delta s \\
&= -\int_{\partial\Omega} d^2r \cdot (\phi_1 \rho_2\delta s)+\int_\Omega d^3r \phi_1 \nabla\cdot (\rho_2\delta s) \\
&= -\int_{\partial\Omega} d^2r (\phi_1 \delta \sigma_2)-\int_\Omega d^3r \phi_1 \delta\rho_2\\
&= -\delta\left[\int_{\partial\Omega} d^2r (\phi_1 \sigma_2)+\int_\Omega d^3r \phi_1 \rho_2\right]
\end{align}
$$
I have set $\rho_2$ the volume charge density in the bulk and $\sigma_2$ the surface charge density. From the conservation of charge, in the bulk there is the continuity equation:
$$
\delta \rho_2 +\nabla\cdot(\rho_2\delta s)=0
$$
and at the surface an analogous relation due to the accumulation of charge when there is a non zero flux with $n$ the outward normal vector at the boundary:
$$
\delta\sigma_2=n\cdot(\rho_2\delta s)
$$
This gives the full expression of the potential energy:
$$
V=\int_{\partial\Omega} d^2r (\phi_1 \sigma_2)+\int_\Omega d^3r \phi_1 \rho_2
$$
In general, gauge invariance is a good sanity check. Here, it’s just that $\phi_1$ is defined up to an additive constant. This gauge transformation is consistent with $V$ (adding the same constant) iff you have conservation of charge, more specifically $0$ net charge:
$$
\int_{\partial\Omega} d^2r \sigma_2+\int_\Omega d^3r \rho_2=0
$$
Back to the infinite case, your boundary term when doing the integration by parts can be important if there is a net flux to “infinty.” This can be the case from a simple translation of charge if charge density does not decay fast enough.
Hope this helps.
