not specifically answering the question, but i think it's usefull to understand why this integral represents potential energy between 2 charge distributions. Another way to look at this is 

$$\iiint \epsilon_0\vec{E_1}.\vec{E_2}d\tau$$
$$\iiint \epsilon_0\vec{E_1}.(-\nabla V_{2})d\tau$$

Using 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}$$


$$\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( -\nabla \cdot (V_{2}\vec{E_1}) + V_{2} \nabla \cdot \vec{E_1}  )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$$

Envoking stokes theorem on the first integral

$$\iint-\epsilon_0 (V_{2}\vec{E_1}) \cdot d\vec{S}+
\iiint\epsilon_0 V_{2} \nabla \cdot \vec{E_1}  d\tau$$

Lim V goes to infinity the surface integral evaluates to zero

$$\iiint\epsilon_0V_{2} \nabla \cdot \vec{E_1}  d\tau$$

Envoking gauss law

$$\iiint\epsilon_0 V_{2} \frac{\rho_1}{\epsilon_0}  d\tau$$


$$\iiint V_{2} \rho_1 d\tau$$


In this form it is clear to see why it represents potential energy. As you are building up  charge distribution 1 in the presence of potential 2