Electric fields due to $Q_1$ and $Q_2$:
$$\vec{E}_1=\frac{Q_1}{4\pi\epsilon_0}\frac{x\hat{x}+y\hat{y}+z\hat{z}}{\left(x^2+y^2+z^2\right)^{3/2}},$$
$$\vec{E}_2=\frac{Q_2}{4\pi\epsilon_0}\frac{x\hat{x}+y\hat{y}+(z-R)\hat{z}}{\left(x^2+y^2+(z-R)^2\right)^{3/2}}.$$
Dot product of electric fields in spherical coordinates:
$$\begin{align}
\vec{E}_1\cdot\vec{E}_2
&=\frac{Q_1Q_2}{16\pi^2\epsilon_0^2}\frac{x^2+y^2+z^2-Rz}{\left(x^2+y^2+z^2\right)^{3/2}\left(x^2+y^2+(z-R)^2\right)^{3/2}}\\
&=\frac{Q_1Q_2}{16\pi^2\epsilon_0^2}\frac{r^2-Rr\cos{\theta}}{r^3\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}\\
&=\frac{Q_1Q_2}{16\pi^2\epsilon_0^2}\frac{r-R\cos{\theta}}{r^2\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}.\\
\end{align}$$
Setting up the integral:
$$\begin{align}
U
&=\epsilon_0\int_{V}\vec{E}_1\cdot\vec{E}_2\,\mathrm{d}V\\
&=\epsilon_0\int_{0}^{2\pi}\int_{0}^{\pi}\int_{0}^{\infty}\vec{E}_1\cdot\vec{E}_2\,r^2\sin{\theta}\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\varphi\\
&=\epsilon_0\int_{0}^{2\pi}\int_{0}^{\pi}\int_{0}^{\infty}\frac{Q_1Q_2}{16\pi^2\epsilon_0^2}\frac{r-R\cos{\theta}}{r^2\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}\,r^2\sin{\theta}\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\varphi\\
&=\frac{Q_1Q_2}{16\pi^2\epsilon_0}\int_{0}^{2\pi}\int_{0}^{\pi}\int_{0}^{\infty}\frac{r-R\cos{\theta}}{\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}\,\sin{\theta}\,\mathrm{d}r\,\mathrm{d}\theta\,\mathrm{d}\varphi\\
&=\frac{Q_1Q_2}{8\pi\epsilon_0}\int_{0}^{\pi}\int_{0}^{\infty}\frac{r-R\cos{\theta}}{\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}\,\sin{\theta}\,\mathrm{d}r\,\mathrm{d}\theta\\
\end{align}$$
Substituting $r=Rx$ and $t=\cos{\theta}$,
$$\begin{align}
U
&=\frac{Q_1Q_2}{8\pi\epsilon_0}\int_{0}^{\pi}\int_{0}^{\infty}\frac{r-R\cos{\theta}}{\left(r^2-2Rr\cos{\theta}+R^2\right)^{3/2}}\,\sin{\theta}\,\mathrm{d}r\,\mathrm{d}\theta\\
&=\frac{Q_1Q_2}{8\pi\epsilon_0}\int_{0}^{\pi}\int_{0}^{\infty}\frac{Rx-R\cos{\theta}}{\left(R^2x^2-2R^2x\cos{\theta}+R^2\right)^{3/2}}\,\sin{\theta}\,R\,\mathrm{d}x\,\mathrm{d}\theta\\
&=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\pi}\int_{0}^{\infty}\frac{x-\cos{\theta}}{\left(x^2-2x\cos{\theta}+1\right)^{3/2}}\,\sin{\theta}\,\mathrm{d}x\,\mathrm{d}\theta\\
&=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\infty}\int_{0}^{\pi}\frac{x-\cos{\theta}}{\left(x^2-2x\cos{\theta}+1\right)^{3/2}}\,\sin{\theta}\,\mathrm{d}\theta\,\mathrm{d}x\\
&=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\infty}\int_{-1}^{1}\frac{x-t}{\left(x^2-2xt+1\right)^{3/2}}\,\mathrm{d}t\,\mathrm{d}x\\
&=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\infty}\left[\frac{x-1}{x^2\sqrt{x^2-2x+1}}+\frac{x+1}{x^2\sqrt{x^2+2x+1}}\right]\,\mathrm{d}x\\
&=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\infty}\left[\frac{x-1}{x^2|x-1|}+\frac{x+1}{x^2|x+1|}\right]\,\mathrm{d}x\\
&=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{\infty}\left[\frac{\operatorname{sgn}{(x-1)}}{x^2}+\frac{1}{x^2}\right]\,\mathrm{d}x\\
&=\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{0}^{1}\left[\frac{-1}{x^2}+\frac{1}{x^2}\right]\,\mathrm{d}x+\frac{Q_1Q_2}{8\pi\epsilon_0\,R}\int_{1}^{\infty}\left[\frac{1}{x^2}+\frac{1}{x^2}\right]\,\mathrm{d}x\\
&=\frac{Q_1Q_2}{4\pi\epsilon_0\,R}\int_{1}^{\infty}\frac{\mathrm{d}x}{x^2}\\
&=\frac{Q_1Q_2}{4\pi\epsilon_0\,R},
\end{align}$$
as expected.
Jacob, you probably won't like this, but the easiest way I see to do the calculation rigorously is to use differential forms. (One reference that's somewhat accessible is my textbook Multivariable Mathematics ..., but you can find plenty of others.)
If $S$ is a closed (oriented) surface in $\Bbb R^3$ not containing the origin and $S^2$ is the unit sphere, then we consider the mapping $f\colon S\to S^2$ given by $f(\mathbf x)=\dfrac{\mathbf x}{\|\mathbf x\|}=\dfrac{\mathbf x}r$. Let $\omega = x\,dy\wedge dz + y\,dz\wedge dx + z\,dx\wedge dy$ be the area $2$-form on $S^2$. (We can also think of $\omega$ as the restriction to the unit sphere of the $2$-form $$\eta =\frac{x\,dy\wedge dz + y\,dz\wedge dx + z\,dx\wedge dy}{r^3}$$ on $\Bbb R^3-\{0\}$. This is the $2$-form corresponding to the electric field of a point charge with strength $1$ at the origin.) It follows from the change of variables theorem that
$$\int_S f^*\omega = \int_{S^2}\omega = 4\pi.$$
(This takes care of the subtle cancellation issues you were worrying about when $f$ is not one-to-one. Because $S$ is a smooth surface, the projection map $f$ has degree $1$.)
The crucial calculation (which embodies the surface area statement in your question) is this: $f^*\omega = \eta$. To verify this, you'll need to know that $f^*(\mathbf x) = f(\mathbf x) = \dfrac{\mathbf x}r$, and $f^*(d\mathbf x) = df = \dfrac{d\mathbf x}r - \dfrac{\mathbf x}{r^2}dr$. So (take a deep breath), remembering that $\wedge$ is skew-symmetric:
\begin{align*}f^*\omega &= \frac xr\left(\frac{dy}r-\frac y{r^2}dr\right)\wedge\left(\frac{dz}r-\frac z{r^2}dr\right)+\frac yr\left(\frac{dz}r-\frac z{r^2}dr\right)\wedge\left(\frac{dx}r-\frac x{r^2}dr\right)+\\& \hspace{1.5in}\frac zr\left(\frac{dx}r-\frac x{r^2}dr\right)\wedge\left(\frac{dy}r-\frac y{r^2}dr\right)\\
&= \eta - \frac1{r^4}\left(xz\,dy\wedge dr+xy\,dr\wedge dz + xy\,dz\wedge dr + yz\,dr\wedge dx + yz\,dx\wedge dr + xz\,dr\wedge dy\right) \\
&= \eta,\end{align*}
as required.
For further motivation, in spherical coordinates $\eta = \sin\phi\,d\phi\wedge d\theta$ (using mathematicians' convention that $\phi$ is the angle from the positive $z$-axis). But you still ultimately need to use my remark above that the degree of $f$ is $1$.
Best Answer
The multinomial theorem, for the case $n=3$, says that $$(x_1 + x_2 + x_3)^n = \sum_{k_1+k_2+k_3 = n} \frac{n!}{k_1!k_2!k_3!} x_1^{k_1} x_2^{k_2} x_3^{k_3} $$ This generalizes to possibly non-integer real exponents $r$ by writing $r! = \Gamma(r+1)$ and $$ \sum_{k_1,k_2,k_3 \in \Bbb Z} \frac{r!}{k_1!k_2!k_3!(r-k_1-k_2-k_3)!} x_1^{k_1} x_2^{k_2} x_3^{k_3} $$
(Notice that in the case of positive integer $r$ this series ends because a $(-1)!$ starts appearing in the denominator, and that is infinite leading to the terms being zero.)
In this case, $(-\frac12)! = \Gamma(\frac12)$ is $\sqrt{\pi}$ and the factors of $\sqrt{\pi}$ will cancel.
For the expression at hand, the lowest order terms in an expansion for $a << r$ are: $$ \left[ 1+ \frac{2a}{r}\cos\theta + \left(\frac{a}{r} \right)^2 \right]^{-\frac12}= 1 - \frac{a}{2r} + \left( \frac38-\cos\theta\right)\frac{a^2}{r^2} - \left( \frac5{16}-\frac32\cos\theta\right)\frac{a^3}{r^3}+ \cdots $$