I think this is the usual argument to derive electrostatic field energy. It's not that you consider just the free motion of the charges, for which total energy and momentum are of course conserved. BTW from this you can also much more elegantly and generally derive the field energy for the general case, using Maxwell's equations and the Lorentz force between charge-current distributions.
The thought experiment is as follows: You consider point charges first. In the thought experiment all the point charges are kept at infinity. Then you put charge ##Q_1## to its place ##\vec{x}_1##. For this you don't need any work, because there's no force acting on it. Then you put the 2nd charge at its place ##\vec{x}_2##, keeping the 1st charge fixed at its position (!). On this 2nd charge the Coulomb force is acting, i.e.,
$$\vec{F}_2=-Q_2 \nabla \Phi_1, \quad \Phi_1=-\frac{Q_1 Q_2}{4 \pi \epsilon_0 |\vec{x}-\vec{x}_1|}.$$
The total work is independent of the path, because the electrostatic force is conservative. Thus the work done finally is
$$E_{\text{mech}}=-\frac{Q_1 Q_2}{4 \pi \epsilon_0 |\vec{x}_2-\vec{x}_1|}.$$
So if ##Q_1 Q_2>0## you have to use energy (because of the repulsive force) and if ##Q_1 Q_2## you get energy back (because of the attractive force). That means the change in total field energy is
$$E_{\text{field}}^{(2)}=\frac{Q_1 Q_2}{4 \pi \epsilon_0 |\vec{x}_1-\vec{x}_2|}.$$
Now you keep the two charges fixed and move the 3rd one at its position ##\vec{x}_3##. Now you have the superposition of the Coulomb forces from particle 1 and 2 and thus for three particles you get
$$E_{\text{field}}^{(3)}= \frac{Q_1 Q_2}{4 \pi \epsilon_0 |\vec{x}_1-\vec{x}_2|} + \frac{Q_1 Q_3}{4 \pi \epsilon_0 |\vec{x}_1-\vec{x}_3|}+ \frac{Q_2 Q_3}{4 \pi \epsilon_0 |\vec{x}_2-\vec{x}_3|},$$
where the first term is from bringing the 2nd charge at its position in presence of the charge 1 as discussed before and the other 2 terms from bringing the 3rd charge at its position in the presence of charges 1 and 2. This argument you can repeat. So the total change of field energy is
$$E_{\text{field}}^{(N)}=\sum_{j=1}^N \sum_{k<j} \frac{Q_j Q_k}{4 \pi \epsilon_0 |\vec{x}_j-\vec{x}_k|} = \frac{1}{2} \sum_{j \neq k} \frac{Q_j Q_k}{4 \pi \epsilon_0 |\vec{x}_j-\vec{x}_k|}.$$
It is important to note that we did not count the "self-energy" of each point charge, i.e., not the energy contained in the electric field of charge 1 in the first step of the gedanken experiment nor the self-energies of the other charges. The reason is that these self-energies are infinite and also unimportant for the energy balance between the field an charges, since only energy differences between different configurations are physically observable but not absolute energies.
That's also the reason that for continuous charge distributions ##\rho## you can simply write
$$E_{\text{field}}=\frac{1}{2} \int_{\mathbb{R}^3} \mathrm{d}^3 r \int_{\mathbb{R}^3} \mathrm{d}^3 r' \frac{\rho(\vec{r}) \rho(\vec{r}')}{4 \pi \epsilon_0 |\vec{r}-\vec{r}'|},$$
where now the "self-energy" is contained in this expression, but it fulfills the purpose for the energy balance as well and is not divergent, because the singularities at ##|\vec{r}-\vec{r}'|=0## are integrable.
This expresses the field energy from the point of view of the charge distribution making up this field, which is of course due to our "mechanistic" approach defining it. As usual in electrodynamics you can express such quantities also using the fields (in macroscopic electrodynamics you can to a certain degree even arbitrarily shuffle parts of the quantities between "matter" and "field" descriptions, but that's another story). To see this we use the fact that the electrostatic potential is given by
$$\Phi(\vec{x})=\int_{\mathbb{R}^3} \mathrm{d}^3 r' \frac{\rho(\vec{r}')}{4 \pi \epsilon_0 |\vec{r}-\vec{r}'|}.$$
Plugging this in our expression for the field energy we get
$$E_{\text{field}}=\frac{1}{2} \int_{\mathbb{R}^3} \mathrm{d}^3 r \rho(\vec{r}) \Phi(\vec{r}).$$
On the other hand from Gauss's Law we have
$$\rho(\vec{r})=\epsilon_0 \vec{\nabla} \cdot \vec{E}(\vec{r}).$$
Plugging this into the above integral and integrating by parts (assuming that the charge distribution vanishes outside a finite region), we get
$$E_{\text{field}}=\frac{\epsilon_0}{2} \int_{\mathbb{R}^3} \mathrm{d}^3 r \vec{E}^2(\vec{r}).$$
This suggests that
$$u_{\text{field}}(\vec{x})=\frac{\epsilon_0}{2} \vec{E}^2(\vec{x})$$
is the energy density of an electrostatic field, and this can indeed also justified by using the local form of energy conservation to derive the general expressions for energy density and energy-current density from Maxwell's equations and the Lorentz-force Law, but that's another story.
