How do I express the interaction energy between two charged spheres?

Consider two identical insulating spheres each with radius $R$ and uniform charge $Q$ through their volume. They are separated from their centers by a distance of $d>2R$.

Here is my general equation for potential energy

$$U_E = \frac{1}{2}\int_V \rho(r)\phi(r)\mathrm{d}\tau$$

How would I express the interaction energy between the two spheres? With this, would I use $\rho$ from sphere 2 and $\phi$ from sphere 1, evaluated at 2?

Is this going in the right direction? I found the potential energy inside the sphere using and to come up with but I believe this is for a single sphere, and am unsure how the second sphere would affect my results. Do I ignore the second sphere altogether?

You have an expression for the total energy

$$U_E = \frac{1}{2}\int_V \rho(\vec{r})\phi(\vec{r})\mathrm{d}^3x.$$

And now you can break the charge density into two parts $\rho=\rho_1+\rho_2$ and the potential into two parts $\phi=\phi_1+\phi_2$, where each is due to just one sphere (so the additions hold by linearity of charge and superposition of fields/potential). Put them in and you get

$$U_E = \frac{1}{2}\int_V (\rho_1(\vec{r})+\rho_2(\vec{r}))(\phi_1(\vec{r})+\phi_2(\vec{r}))\mathrm{d}^3x.$$

Which equals

$$U_E = \frac{1}{2}\int_V \rho_1(\vec{r})\phi_1(\vec{r})+\rho_2(\vec{r})\phi_2(\vec{r})+ \rho_2(\vec{r})\phi_1(\vec{r})+ \rho_1(\vec{r})\phi_2(\vec{r})\mathrm{d}^3x.$$

The first two terms integrate out to be the energy of the spheres in isolation. So the last two terms are the interaction energy:

$$U_{int} = \frac{1}{2}\int_V \rho_2(\vec{r})\phi_1(\vec{r})+ \rho_1(\vec{r})\phi_2(\vec{r})\mathrm{d}^3x.$$

Now, each term integrates to tell you the total work done by one sphere on the other sphere as you (hypothetically) bring the isolated objects from super far away to their current positions. Since the forces were equal and opposite, the work done by one on the other is equal to the work done by the other on the one.

So the remaining two terms actually give equal integrals.

$$U_{int} = \frac{1}{2}\int_V \rho_2(\vec{r})\phi_1(\vec{r})+ \rho_1(\vec{r})\phi_2(\vec{r})\mathrm{d}^3x= \int_V \rho_1(\vec{r})\phi_2(\vec{r})\mathrm{d}^3x= \int_V \rho_2(\vec{r})\phi_1(\vec{r})\mathrm{d}^3x.$$

If what you want is total electrostatic energy of the system, then take $\rho,\phi$ to be total charge density and total electrostatic potential (due to both spheres) and calculate the integral as you wrote it. If you want just the potential energy of the sphere 2 in the (fixed) field of the sphere 1, then take $\rho$ to be charge density of the sphere 2 only and $\phi$ due to sphere 1 only, and calculate the integral

$$U_{21} = \int_V \rho_2 \phi_1 d\tau$$

without the factor $1/2$.