Why is the total electrostatic energy bigger than the sum of separate energies? I don't understand why, if we add the two energies $\,\boxed{1}\,$ and $\,\boxed{2}\,$ from two electric fields, the component $\,\boxed{3}\,$ is added and so it's different from their sum
\begin{equation}
E=\int\tfrac12\epsilon_0\left(\mathbf E_1+\mathbf E_2\right)^2\mathrm dV=\underbrace{\int\tfrac12\epsilon_0\mathbf E_1^2\,\mathrm dV}_{\boxed{1}}+\underbrace{\int\tfrac12\epsilon_0\mathbf E_2^2\,\mathrm dV}_{\boxed{2}}+\underbrace{\int\epsilon_0\mathbf E_1\cdot\mathbf E_2\,\mathrm dV}_{\boxed{3}}  
\tag{01}\label{01}
\end{equation}
 A: From the formal point of view, the charging process going from zero-field to the total field ${\bf E}_1 + {\bf E}_2$ always gives an energy density
$$
\cal E = \varepsilon_0\int_0^{{\bf E}_1 + {\bf E}_2} {\bf E}\cdot {\mathrm d}{\bf E}=\frac{\varepsilon_0}{2}\left(  {\bf E}_1 + {\bf E}_2  \right)^2.
$$
If we split the process into a charging from zero to ${\bf E}_1$ and from ${\bf E}_1$ to ${\bf E}_1 + {\bf E}_2 $, the first part of the process gives an energy density proportional to ${\bf E}^2_1$, but the second is proportional to $\left(  {\bf E}_1 + {\bf E}_2  \right)^2 - {\bf E}^2_1$, differing from
${\bf E}^2_2$ by a term $2 {\bf E}_1 \cdot {\bf E}_2$.
What about intuition?
Apart from the bare rewording of the result of integration, there is an important conceptual reason for the non-additivity of the energy density: electrostatic energies originate from long-range coulomb interactions between point-like charges. Generally, at least for macroscopic samples, the additivity of energies is only justified if the interactions between microscopic constituents are short-range. In the presence of long-range interactions, it is not possible to justify that the total energy of a system made by a part with energy $\cal E_1$ and another with energy $\cal E_2$ is simply
$\cal E_1+ \cal E_2$.
A: In electrostatics:
$$\iiint_{\tau}\epsilon_0\vec{E}_{1} \cdot \vec{E}_{2} d\tau$$
$$=\iiint \epsilon_0\vec{E}_1\cdot(-\nabla V_{2})\,d\tau.$$
Rearranging the vector field identity
$$\nabla \cdot (V_{2}\vec{E}_1) = V_{2} \nabla \cdot \vec{E}_1 + \vec{E}_{1} \cdot \nabla V_{2},$$
gives
$$\vec{E}_{1} \cdot (-\nabla V_{2})= -\nabla\cdot(V_{2}\vec{E}_1) + V_{2} \nabla \cdot\vec{E}_1.$$
Substituting this in the integral gives
$$\iiint\epsilon_0\left[-\nabla \cdot (V_{2}\vec{E}_1) + V_{2} \nabla \cdot \vec{E}_1\right]d\tau$$
$$=-\iiint\epsilon_0 \nabla \cdot (V_{2}\vec{E}_1)\,d\tau+\iiint \epsilon_0V_{2} \nabla \cdot \vec{E}_{1}\,d\tau.$$
Invoking Stokes' Theorem on the first integral gives
$$-\iint\epsilon_0 (V_{2}\vec{E}_1) \cdot d\vec{S}+
\iiint\epsilon_0 V_{2} \nabla \cdot \vec{E}_1\,d\tau.$$
With the integration volume as all of space, the surface integral evaluates to zero for localized sources, since at large distances $(r\rightarrow\infty)$,
the integrand falls off as $V_{2}\vec{E}_{1}\sim r^{-3}$. This leaves
$$\iiint\epsilon_0V_{2} \nabla \cdot \vec{E}_1\,d\tau.$$
Invoking Gauss' Law,
$$\iiint\epsilon_0 V_{2} \frac{\rho_1}{\epsilon_0}\,d\tau$$
$$=\iiint V_{2} \rho_1\,d\tau$$
This expression is more intuitive, each $V_{2} [\rho_{1} d\tau]$ element, represents an infinitesimal amount of work required to build up a charge $\rho_{1} d\tau$ in the presence of the potential from distribution 2,$[V_{2}]$
[as qV = work needed to  move an object from infinity to its location, instead we are building up a charge dq, and summing all of these elements in an integral]
This expression represents the potential energy between the distributions. Work is required to build up distribution 2  in the presence of distribution 1, as distribution 2 now needs to overcome the electrostatic forces of distributions 1, as aswell as its own field.
A: Consider not just the fields, but the charge and current distributions required to produce them. The first term would be the energy required to assemble the charge and current distribution required for the first field by itself. Similarly the second term is the energy required to assemble the distribution for the second field by itself. However, if the first distribution is already present, then the presence of its field can make assembling the second distribution more or less difficult than it was to assemble it by itself. That difference is the third term.
