A more high-level answer this time. In this approach we first find out the form of the term in the potential that will split the $d$ orbitals.
Recall that the character of a proper rotation by angle $\alpha$ for the set of spherical harmonics with $L^2Y_l^m(\theta,\phi)=\hbar^2l(l+1)Y_l^m(\theta,\phi)$ is
$$\chi_l(R(\alpha))=\sum_{m=-l}^le^{im\alpha}=\frac{\sin(m+\frac12)\alpha}{\sin\frac12\alpha}$$
for $\alpha \neq n\pi$. $\chi_l(R(0))=2l+1$ because all the diagonal elements are $1$ and $\chi_l(R(\pi))=(-1)^l$ because the diagonal elements alternate between $\pm1$ and there is an odd number of them so their sum is one of their boundary values.
Starting from the character table for $O$, we can decompose the characters of the multipole moments via character orthogonality.
$$\begin{array}{c|ccccc|c}
O&E&3C_{2h}&8C_3&6C_4&6C_{2d}&\Gamma \\
\hline A_1&1&1&1&1&1& \\
A_2&1&1&1&-1&-1& \\
E&2&2&-1&0&0& \\
T_1&3&-1&0&1&-1 \\
T_2&3&-1&0&-1&1 \\
\hline S&1&1&1&1&1&A_{1g} \\
P&3&-1&0&1&-1&T_{1u} \\
D&5&1&-1&-1&1&E_g+T_{2g} \\
F&7&-1&1&-1&-1&A_{2u}+T_{1u}+T_{2u} \\
G&9&1&0&1&1&A_{1g}+E_g+T_{1g}+T_{2g} \\
H&11&-1&-1&1&-1&E_u+2T_{1u}+T_{2u} \\
I&13&1&1&-1&1&A_{1g}+A_{2g}+E_g+T_{1g}+2T_{2g}\end{array}$$
We have expanded this all the way out to hexacontatetrapole moments just in case it were desired to investigate the splitting of $f$ orbitals, but the only multipole that matters for $d$ orbitals is the $A_{1g}$ hexadecapole moment.
This is because the product of a spherical harmonic $Y_l^m(\theta,\phi)$ with $Y_l^{m^{\prime}}(\theta,\phi)$ contains only spherical harmonics $Y_{l^{\prime\prime}}^{m^{\prime\prime}}(\theta,\phi)$ with even $l^{\prime\prime}$, $0\le l^{\prime\prime}\le2l$.
We have used the fact that the parity of electric multipole moments is even for even $l$ and odd for odd $l$ to determine the characters over the full $O_h$ group. The reverse is the case for magnetic multipole moments.
We use the projection operator
$$P_{jk}^{\lambda}f(x,y,z)=\frac{\chi^{\lambda}(E)}{|G|}\sum_{R\in G}D_{kj}^{\lambda}\left(R^{-1}\right)Rf(x,y,z)$$
Where $\chi^{\lambda}(E)$ is the number of partners of irreducible representation $\lambda$, $|G|$ is the order of the group $G$, the sum is taken over all elements $R$ of the group, and $D_{kj}^{\lambda}\left(R^{-1}\right)$ is the matrix element of the matrix representing $R^{-1}$.
If $f(x,y,z)$ has any character of the $k^{th}$ partner of the $\lambda^{th}$ irreducible representation of $G$, this will be nonzero and transform as the $j^{th}$ partner of the $\lambda^{th}$ irreducible representation.
As an example, with $G=T_d$, $\lambda=T_2$, and $f(x,y,z)=x+xz+z\left(2z^2-3x^2-3y^2\right)$ the different $P_{jk}^{T_2}f(x,y,z)$ can spit out $9$ different functions.
But for our purposes today we only need to project into $A_{1g}$ because any multipole moments must have the full symmetry of the charge configuration that creates them.
Then all $D_{11}^{A_{1g}}\left(R^{-1}\right)=1$ so
$$P_{11}^{A_{1g}}f(x,y,z)=\frac1{|G|}\sum_{R\in G}Rf(x,y,z)$$
We can simplify the sum by finding the stabilizer $H$ of $F(x,y,z)$ in $G$: $Sf(x,y,z)=f(x,y,z)$ for all $S\in H$.
Then given a set $C$ of representatives of $H$ in $G$, we have
$$P_{11}^{A_{1g}}f(x,y,z)=\frac1{|G|}\sum_{T\in C}\sum_{S\in H}TSf(x,y,z)=\frac{|H|}{|G|}\sum_{T\in C}Tf(x,y,z)$$
Starting with a fourth order polynomial $f(x,y,z)=z^4$, its stabilizer in $G=O_h$ is $H=D_{4h}$, the dihedral group of order $16$ with $z$ axis as principal axis of symmetry.
The cosets can be written as the powers of $C_{3[111]}$, a rotation of $120°$ about the axis through the origin and $\langle1,1,1\rangle$.
$$P_{11}^{A_{1g}}f(x,y,z)=\frac{16}{48}\left(Ez^4+C_{3[111]}z^4+C_{3[111]}^2z^4\right)=\frac13\left(z^4+x^4+y^4\right)$$
We use the Gram-Schmidt process to turn this into a normalized spherical harmonic.
Looking back up at the character table we see that also monopole moments transform as $A_{1g}$, so we start out with $g(x,y,z)=1$ and find its magnitude:
$$\oint\left|1\right|^2d^2\Omega=\frac{2\Gamma\left(\frac12\right)\Gamma\left(\frac12\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac32\right)}=\frac{2\pi^{3/2}}{\frac12\pi^{1/2}}=4\pi$$
Where we have used the formula
$$\oint x^{2a}y^{2b}z^{2c}d^2\Omega=\frac{2\Gamma\left(a+\frac12\right)\Gamma\left(b+\frac12\right)\Gamma\left(c+\frac12\right)}{\Gamma\left(a+b+c+\frac32\right)}$$
to evaluate the integrals over the unit sphere of a monomial with even powers of $x$, $y$, and $z$. If any exponent is odd, the integral will be zero.
Also $\Gamma\left(\frac12\right)=\sqrt{\pi}$ and for $x>0$, $\Gamma(x+1)=x\Gamma(x)$.
Hopefully a reviewer can find a thread like 'What Every Chemistry Student Needs to Know About the Gamma and Beta Functions' and post a link in comments.
So with the normalized spherical harmonic
$$Y_{A_{1g}1}^0(\theta,\phi)=\left(\frac1{4\pi}\right)^{1/2}1$$
We can compute the unnormalized spherical harmonic
$$\begin{align}cY_{A_{1g}1}^4(\theta,\phi) & =\frac1{r^4}P_{11}^{A_{1g}}f(x,y,z)-Y_{A_{1g}1}^0(\theta,\phi)\oint \left[Y_{A_{1g}1}^0(\theta^{\prime},\phi^{\prime})\right]^*P_{11}^{A_{1g}}f(x^{\prime},y^{\prime},z^{\prime})d^2\Omega^{\prime} \\
& =\frac1{3r^4}\left(x^4+y^4+z^4\right)-\frac1{12\pi}\oint\left((x^{\prime})^4+(y^{\prime})^4+(z^{\prime})^4\right)d^2\Omega^{\prime} \\
& =\frac1{3r^4}\left(x^4+y^4+z^4\right)-\frac1{12\pi}(3)\frac{2\Gamma\left(\frac52\right)\Gamma\left(\frac12\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac72\right)} \\
& =\frac1{3r^4}\left(x^4+y^4+z^4\right)-\frac1{4\pi}\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac32\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)} \\
& =\frac1{3r^4}\left(x^4+y^4+z^4\right)-\frac15=\frac1{3r^4}\left(x^4+y^4+z^4-\frac35r^4\right)\end{align}$$
This has magnitude
$$\begin{align}|c|^2 & =\oint \left[cY_{A_{1g}1}^4(\theta,\phi)\right]^*cY_{A_{1g}1}^4(\theta,\phi)d^2\Omega=\frac19\oint\left(x^4+y^4+z^4-\frac35r^4\right)^2d^2\Omega \\
& =\frac19\left[(3)\frac{2\Gamma\left(\frac92\right)\Gamma\left(\frac12\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac{11}2\right)}+(6)\frac{\Gamma\left(\frac52\right)\Gamma\left(\frac52\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac{11}2\right)}\right. \\
& \left.-\frac35(6)\frac{\Gamma\left(\frac52\right)\Gamma\left(\frac12\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac72\right)}+\frac9{25}\frac{\Gamma\left(\frac12\right)\Gamma\left(\frac12\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac32\right)}\right] \\
& =\frac19\left[(3)\frac{4\pi}{9}+(6)\frac{4\pi}{105}-\frac35(6)\frac{4\pi}5+\frac9{25}(4\pi)\right]=\frac{64\pi}{4725}\end{align}$$
So
$$Y_{A_{1g}1}^4(\theta,\phi)=\left(\frac{525}{64\pi}\right)^{1/2}\frac1{r^4}\left(x^4+y^4+z^4-\frac35r^4\right)$$
Given this spherical harmonic, we can find the electrostatic potential readily.
To find the potential at a field point $\vec{r_f}=\left\langle x_f,y_f,z_f\right\rangle$ we add up the contribution of the charge $dq=\rho(\vec{r_s})d^3volume$ at source point $\vec{r_s}=\left\langle x_s.y_s,z_s\right\rangle$ and sum over all source points:
$$V_{O_h}\left(\vec{r_f}\right)=\int\frac{\rho\left(\vec{r_s}\right)}{4\pi\epsilon_0\left\|\vec{r_f}-\vec{r_s}\right\|}d^3volume$$
Here
$$\begin{align}\left\|\vec{r_f}-\vec{r_s}\right\| & =\sqrt{\left(\vec{r_f}-\vec{r_s}\right)\cdot\left(\vec{r_f}-\vec{r_s}\right)}=\sqrt{r_s^2-2\vec{r_s}\cdot\vec{r_f}+r_f^2} \\
& =r_s\sqrt{1-2\frac{\vec{r_s}\cdot\vec{r_f}}{r_s^2}+\frac{r_f^2}{r_s^2}}=r_s\sqrt{1-2\frac{r_f}{r_s}\cos\gamma+\frac{r_f^2}{r_s^2}}\end{align}$$
Where $\cos\gamma=\frac{\vec{r_s}\cdot\vec{r_f}}{r_sr_f}$ is cosine of the angle between the source point and the field point as seen at the origin.
We can use the generating function for the Legendre polynomials and the spherical harmonic addition theorem to expand
$$\begin{align}\frac1{r_s\sqrt{1-2\frac{r_f}{r_s}\cos\gamma+\frac{r_f^2}{r_s^2}}} & =\frac1{r_s}\sum_{l=0}^{\infty}\left(\frac{r_f}{r_s}\right)^lP_l(\cos\gamma) \\
& =\frac1{r_s}\sum_{l=0}^{\infty}\left(\frac{r_f}{r_s}\right)^l\frac{4\pi}{2l+1}\sum_{m=-l}^l\left[Y_l^m\left(\theta_s,\phi_s\right)\right]^*Y_l^m\left(\theta_f,\phi_f\right)\end{align}$$
Valid for $r_f\lt r_s$. If $r_f>r_s$, just replace $\frac{r_f^l}{r_s^{l+1}}$ with $\frac{r_s^l}{r_f^{l+1}}$.
This is OK for spherical symmetry, but in chemistry we want a different basis set adapted to the point group of the complex.
We choose an orthonormal basis set $\left\{Z_l^m(\theta,\phi),m=-l,l\right\}$ so we can write the old basis set as
$$Y_l^m(\theta,\phi)=\sum_{m=-l}^lA_{m^{\prime}m}^lZ_l^{m^{\prime}}(\theta,\phi)$$
Since both bases are orthormal sets,
$$\begin{align}\delta_{mm^{\prime}} & =\oint\left[Y_l^{m^\prime}(\theta,\phi)\right]^*Y_l^m(\theta,\phi)d^2\Omega \\
& =\oint\sum_{m_1=-l}^l\left[A_{m_1m^{\prime}}^lZ_l^{m_1}(\theta,\phi)\right]^*\sum_{m_2=-l}^lA_{m_2m}^lZ_l^{m_2}(\theta,\phi)d^2\Omega \\
& =\sum_{m_1=-l}^l\left(A_{m_1m^{\prime}}^l\right)^*\sum_{m_2=-l}^lA_{m_2m}^l\delta_{m_1m_2}=\sum_{m_1=-l}^l\left(A_{m_1m^{\prime}}^l\right)^*A_{m_1m}^l\end{align}$$
Which says that as matrices, $(A^l)^{\dagger}A^l=I$. Inverses commute because their product being nonsingular requires that both $(A^l)^{\dagger}$ and $A^l$ must be nonsingular and so their product must have an inverse, $\left(A^l(A^l)^{\dagger}\right)^{-1}$. Then
$$\begin{align}A^l(A^l)^{\dagger} & =I\left(A^l(A^l)^{\dagger}\right)=\left[\left(A^l(A^l)^{\dagger}\right)^{-1}\left(A^l(A^l)^{\dagger}\right)\right]\left(A^l(A^l)^{\dagger}\right) \\
& =\left(A^l(A^l)^{\dagger}\right)^{-1}\left[\left(A^l(A^l)^{\dagger}\right)\left(A^l(A^l)^{\dagger}\right)\right]=\left(A^l(A^l)^{\dagger}\right)^{-1}\left\{A^l\left[(A^l)^{\dagger}\left(A^l(A^l)^{\dagger}\right)\right]\right\} \\
& =\left(A^l(A^l)^{\dagger}\right)^{-1}\left\{A^l\left[\left((A^l)^{\dagger}A^l\right)(A^l)^{\dagger}\right]\right\}=\left(A^l(A^l)^{\dagger}\right)^{-1}\left[A^l\left(I(A^l)^{\dagger}\right)\right] \\
& =\left(A^l(A^l)^{\dagger}\right)^{-1}\left(A^l(A^l)^{\dagger}\right)=I\end{align}$$
As matrix elements, this is
$$\sum_{m_1=-l}^lA_{m^{\prime}m_1}^l\left(A_{mm_1}^l\right)^*=\delta_{mm^{\prime}}$$
Then
$$\begin{align}\sum_{m=-l}^l\left[Y_l^m\left(\theta_s,\phi_s\right)\right]^*Y_l^m\left(\theta_f,\phi_f\right) & =\sum_{m=-l}^l\sum_{m_1=-l}^l\left[A_{m_1m}^lZ_l^{m_1}\left(\theta_s,\phi_s\right)\right]^*\sum_{m_2=-l}^lA_{m_2m}^lZ_l^{m_2}\left(\theta_f,\phi_f\right) \\
& =\sum_{m_1=-l}^l\left[Z_l^{m_1}\left(\theta_s,\phi_s\right)\right]^*\sum_{m_2=-l}^lZ_l^{m_2}\left(\theta_f,\phi_f\right)\delta_{m_2m_1} \\
& =\sum_{m_1=-l}^l\left[Z_l^{m_1}\left(\theta_s,\phi_s\right)\right]^*Z_l^{m_1}\left(\theta_f,\phi_f\right)\end{align}$$
This is good news because it means that we can use the $Y_l^m$ basis to prove the spherical harmonic addition theorem, but then substitute a more convenient basis when it comes time to apply it.
Now we have
$$\begin{align}V_{O_h}\left(\vec{r_f}\right) & =\int_0^{r_f}\int_0^{\pi}\int_0^{2\pi}\frac{\rho\left(\vec{r_s}\right)}{4\pi\epsilon_0}\sum_{l=0}^{\infty}\frac{r_s^l}{r_f^{l+1}}\frac{4\pi}{2l+1}\sum_{m=-l}^l\left[Y_l^m\left(\theta_s,\phi_s\right)\right]^*Y_l^m\left(\theta_f,\phi_f\right)d\phi_s\sin\theta_sd\theta_sr_s^2dr_s \\
& +\int_{r_f}^{\infty}\int_0^{\pi}\int_0^{2\pi}\frac{\rho\left(\vec{r_s}\right)}{4\pi\epsilon_0}\sum_{l=0}^{\infty}\frac{r_f^l}{r_s^{l+1}}\frac{4\pi}{2l+1}\sum_{m=-l}^l\left[Y_l^m\left(\theta_s,\phi_s\right)\right]^*Y_l^m\left(\theta_f,\phi_f\right)d\phi_s\sin\theta_sd\theta_sr_s^2dr_s\end{align}$$
We confine our interest to the region near the origin so $r_s>r_f$ and the first integral above will be zero.
Also we are just trying to derive a formula for the splitting of the $d$ orbitals, so we cut our sums off at $l=4$.
Only the $l=0$ and $l=4$ irreducible representations of the full rotation group have a partner that transforms as $A_{1g}$ in $O_h$ and we have those partners to hand at this point.
Given the locations of the $6$ charges $-Ze$ at distance $a$ along the coordinate axes we can write $\rho\left(\vec{r_s}\right)=-Ze\sum\limits_{j=0}^5\delta^{(3)}\left(\vec{r_s}-S_6^j\langle0,0,a\rangle\right)$ where $S_6$ is a $60°$ rotoreflection about the $[111]$ axis, so
$$\begin{align}V_{O_h}\left(\vec{r_f}\right) & =\int_{r_f}^{\infty}\int_0^{\pi}\int_0^{2\pi}\frac{-Ze}{\epsilon_0}\frac1{r_s}\sum\limits_{j=0}^5\delta^{(3)}\left(\vec{r_s}-S_6^j\langle0,0,a\rangle\right)\left(\frac1{4\pi}\right)d\phi_s\sin\theta_sd\theta_sr_s^2dr_s \\
& +\int_{r_f}^{\infty}\int_0^{\pi}\int_0^{2\pi}\frac{-Ze}{\epsilon_0}\frac{r_f^4}{9r_s^5}\sum\limits_{j=0}^5\delta^{(3)}\left(\vec{r_s}-S_6^j\langle0,0,a\rangle\right)\left(\frac{525}{64\pi}\right)\times \\
& \frac1{r_s^4}\left(x_s^4+y_s^4+z_s^4-\frac35r_s^4\right)\frac1{r_f^4}\left(x_f^4+y_f^4+z_f^4-\frac35r_f^4\right)d\phi_s\sin\theta_sd\theta_sr_s^2dr_s \\
& =\frac{-Ze}{\epsilon_0a}(6)\left(\frac1{4\pi}\right)-\frac{Zer_f^4}{9\epsilon_0a^5}(6)\left(\frac{525}{64\pi}\right)\frac25\frac1{r_f^4}\left(x_f^4+y_f^4+z_f^4-\frac35r_f^4\right) \\
& =-\frac{6Ze}{4\pi\epsilon_0a}-\frac{35Ze}{16\pi\epsilon_0a^5}\left(x_f^4+y_f^4+z_f^4-\frac35r_f^4\right)\end{align}$$
In the above we used the $3$ dimensional Dirac delta functions to do all the integrals.
An $A_{1u}$ representation in $O_h$ is $A+1$ for the $T_d$ case, so we might want to consider such a multipole moment there.
Looking for the stabilizer of $xyz$ in $T_d$, we see that it's the whole group, so it already transforms as $A_1$ in $T_d$.
It's also already orthogonal to $1$ so $xyz=cY_{A_11}^3$. Finding the magnitude,
$$\begin{align}|c|^2 & =\oint\left|cY_{A_11}^3\right|^2d^2\Omega=\oint x^2y^2z^2d^2\Omega=\frac{2\Gamma\left(\frac32\right)\Gamma\left(\frac32\right)\Gamma\left(\frac32\right)}{\Gamma\left(\frac92\right)} \\
& =\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac12\right)\left(\frac12\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)\left(\frac72\right)}=\frac{4\pi}{105}\end{align}$$
So our normalized spherical harmonic is
$$Y_{A_11}^3(\theta,\phi)=\left(\frac{105}{4\pi}\right)^{1/2}\frac{xyz}{r^3}$$
The charge distribution for $T_d$ geometry is
$$\rho\left(\vec{r_s}\right)=-Ze\sum_{j=0}^3\delta^{(3)}\left(\vec{r_s}-S_4^j\left\langle\frac a{\sqrt3},\frac a{\sqrt3},\frac a{\sqrt3}\right\rangle\right)$$
Where $S_4$ is a $90°$ rotoreflection about the $z$ axis. Making the same simplifications as for the $O_h$ case, we have
$$\begin{align}V_{T_d}\left(\vec{r_f}\right) & =\int_{r_f}^{\infty}\int_0^{\pi}\int_0^{2\pi}\frac{-Ze}{\epsilon_0}\frac1{r_s}\sum_{j=0}^3\delta^{(3)}\left(\vec{r_s}-S_4^j\left\langle\frac a{\sqrt3},\frac a{\sqrt3},\frac a{\sqrt3}\right\rangle\right)\left(\frac1{4\pi}\right)d\phi_s\sin\theta_sd\theta_sr_s^2dr_s \\
& =\int_{r_f}^{\infty}\int_0^{\pi}\int_0^{2\pi}\frac{-Ze}{\epsilon_0}\frac{r_f^3}{7r_s^4}\sum_{j=0}^3\delta^{(3)}\left(\vec{r_s}-S_4^j\left\langle\frac a{\sqrt3},\frac a{\sqrt3},\frac a{\sqrt3}\right\rangle\right)\left(\frac{105}{4\pi}\right)\times \\
& \frac{x_sy_sz_s}{r_s^3}\frac{x_fy_fz_f}{r_f^3}d\phi_s\sin\theta_sd\theta_sr_s^2dr_s \\
& +\int_{r_f}^{\infty}\int_0^{\pi}\int_0^{2\pi}\frac{-Ze}{\epsilon_0}\frac{r_f^4}{9r_s^5}\sum_{j=0}^3\delta^{(3)}\left(\vec{r_s}-S_4^j\left\langle\frac a{\sqrt3},\frac a{\sqrt3},\frac a{\sqrt3}\right\rangle\right)\left(\frac{525}{64\pi}\right)\times \\
& \frac1{r_s^4}\left(x_s^4+y_s^4+z_s^4-\frac35r_s^4\right)\frac1{r_f^4}\left(x_f^4+y_f^4+z_f^4-\frac35r_f^4\right)d\phi_s\sin\theta_sd\theta_sr_s^2dr_s \\
& =\frac{-Ze}{\epsilon_0a}(4)\left(\frac1{4\pi}\right)-\frac{Zer_f^3}{7\epsilon_0a^4}(4)\left(\frac{105}{4\pi}\right)\left(\frac1{3\sqrt3}\right)\frac{x_fy_fz_f}{r_f^3} \\
& -\frac{Zer_f^4}{9\epsilon_0a^5}(4)\left(\frac{525}{64\pi}\right)\left(\frac{-4}{15}\right)\frac1{r_f^4}\left(x_f^4+y_f^4+z_f^4-\frac35r_f^4\right) \\
& =-\frac{4Ze}{4\pi\epsilon_0a}-\frac{5Z}{\sqrt3\pi\epsilon_0a^4e}x_fy_fz_f+\frac{35Ze}{36\pi\epsilon_0a^5}\left(x_f^4+y_f^4+z_f^4-\frac35r_f^4\right)\end{align}$$
So now we are ready to get our $d$ orbitals and compute the splitting. Starting with the $d_{xy}$ orbital:
$$\oint x^2y^2d^2\Omega = \frac{2\Gamma\left(\frac32\right)\Gamma\left(\frac32\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac72\right)}=\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac12\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)}=\frac{4\pi}{15}$$
So the normalized $d_{xy}$ orbital is
$$\psi_{xy}\left(\vec{r}\right)=\left(\frac{15}{4\pi}\right)^{1/2}\frac{xy}{r^2}f(r)$$
Where f(r) is assumed normalized to unity. We can compute
$$\begin{align}U_{O_hxy} &=-e\int\left[\psi_{xy}\left(\vec{r}\right)\right]^*V_{O_h}\left(\vec{R}\right)\psi_{xy}\left(\vec{r}\right)d^3volume \\
& =-e\oint\int_0^{\infty}\left[-\frac{6Ze}{4\pi\epsilon_0a}-\frac{35Ze}{16\pi\epsilon_0a^5}\frac1{r^4}\times\right. \\
& \left.\left(x^4+y^4+z^4-\frac35r^4\right)r^4\right]\left(\frac{15}{4\pi}\right)\frac{x^2y^2}{r^4}|f(r)|^2r^2drd^2\Omega \\
& =\frac{6Ze^2}{4\pi\epsilon_0a}+\frac{35Ze^2}{16\pi\epsilon_0a^5}\left\langle r^4\right\rangle\left(\frac{15}{4\pi}\right)\left[(2)\frac{2\Gamma\left(\frac72\right)\Gamma\left(\frac32\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac{11}2\right)}\right. \\
& \left.+\frac{2\Gamma\left(\frac52\right)\Gamma\left(\frac32\right)\Gamma\left(\frac32\right)}{\Gamma\left(\frac{11}2\right)}-\frac35\frac{2\Gamma\left(\frac32\right)\Gamma\left(\frac32\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac72\right)}\right] \\
& =\frac{6Ze^2}{4\pi\epsilon_0a}+\frac{525Ze^2}{64\pi^2\epsilon_0a^5}\left\langle r^4\right\rangle\left[(2)\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)\left(\frac12\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)\left(\frac72\right)\left(\frac92\right)}\right. \\
& \left.+\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac32\right)\left(\frac12\right)\left(\frac12\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)\left(\frac72\right)\left(\frac92\right)}-\frac35\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac12\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)}\right] \\
& =\frac{6Ze^2}{4\pi\epsilon_0a}-\frac{Ze^2}{6\pi\epsilon_0a^5}\left\langle r^4\right\rangle\end{align}$$
On to the $d_{x^2-y^2}$ orbital:
$$\begin{align}\oint \left(x^2-y^2\right)^2d^2\Omega &= (2)\frac{2\Gamma\left(\frac52\right)\Gamma\left(\frac12\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac72\right)}-(2)\frac{2\Gamma\left(\frac32\right)\Gamma\left(\frac32\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac72\right)} \\
& =(2)\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac32\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)}-(2)\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac12\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)}=\frac{16\pi}{15}\end{align}$$
So the normalized $d_{x^2-y^2}$ orbital is
$$\psi_{x^2-y^2}\left(\vec{r}\right)=\left(\frac{15}{16\pi}\right)^{1/2}\frac{x^2-y^2}{r^2}f(r)$$
We have
$$\begin{align}U_{O_hx^2-y^2} &=-e\int\left[\psi_{x^2-y^2}\left(\vec{r}\right)\right]^*V_{O_h}\left(\vec{R}\right)\psi_{x^2-y^2}\left(\vec{r}\right)d^3volume \\
& =-e\oint\int_0^{\infty}\left[-\frac{6Ze}{4\pi\epsilon_0a}-\frac{35Ze}{16\pi\epsilon_0a^5}\frac1{r^4}\times\right. \\
& \left.\left(x^4+y^4+z^4-\frac35r^4\right)r^4\right]\left(\frac{15}{16\pi}\right)\frac{\left(x^2-y^2\right)^2}{r^4}|f(r)|^2r^2drd^2\Omega \\
& =\frac{6Ze^2}{4\pi\epsilon_0a}+\frac{35Ze^2}{16\pi\epsilon_0a^5}\left\langle r^4\right\rangle\left(\frac{15}{16\pi}\right)\left[(2)\frac{2\Gamma\left(\frac92\right)\Gamma\left(\frac12\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac{11}2\right)}-(4)\frac{2\Gamma\left(\frac72\right)\Gamma\left(\frac32\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac{11}2\right)}\right. \\
& \left.+(4)\frac{2\Gamma\left(\frac52\right)\Gamma\left(\frac52\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac{11}2\right)}-(2)\frac{2\Gamma\left(\frac52\right)\Gamma\left(\frac32\right)\Gamma\left(\frac32\right)}{\Gamma\left(\frac{11}2\right)}\right. \\
& \left.-\frac35(2)\frac{2\Gamma\left(\frac52\right)\Gamma\left(\frac12\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac72\right)}+\frac35(2)\frac{2\Gamma\left(\frac32\right)\Gamma\left(\frac32\right)\Gamma\left(\frac12\right)}{\Gamma\left(\frac72\right)}\right] \\
& =\frac{6Ze^2}{4\pi\epsilon_0a}+\frac{525Ze^2}{256\pi^2\epsilon_0a^5}\left\langle r^4\right\rangle\left[(2)\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)\left(\frac72\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)\left(\frac72\right)\left(\frac92\right)}\right. \\
& \left.-(4)\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)\left(\frac12\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)\left(\frac72\right)\left(\frac92\right)}+(4)\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac32\right)\left(\frac12\right)\left(\frac32\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)\left(\frac72\right)\left(\frac92\right)}\right. \\
& \left.-(2)\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac32\right)\left(\frac12\right)\left(\frac12\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)\left(\frac72\right)\left(\frac92\right)}-\frac35(2)\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac32\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)}\right. \\
& \left.+\frac35(2)\frac{2\pi^{3/2}\left(\frac12\right)\left(\frac12\right)}{\pi^{1/2}\left(\frac12\right)\left(\frac32\right)\left(\frac52\right)}\right] \\
& =\frac{6Ze^2}{4\pi\epsilon_0a}+\frac{Ze^2}{4\pi\epsilon_0a^5}\left\langle r^4\right\rangle\end{align}$$
With
$$Dq=\left(\frac{35Ze}{16\pi\epsilon_0a^5}\right)\left(\frac{2e}{105}\right)\left\langle r^4\right\rangle=\frac{Ze^2}{24\pi\epsilon_0a^5}\left\langle r^4\right\rangle$$
We see that
$$U_{O_hxy}=\frac{6Ze^2}{4\pi\epsilon_0a}-4Dq$$
and
$$U_{O_hx^2-y^2}=\frac{6Ze^2}{4\pi\epsilon_0a}+6Dq$$
For the $T_d$ case, the $Y_{A_11}^3$ multipole moment won't contribute to $U_{T_dxy}$ or $U_{T_dx^2-y^2}$ because it has odd parity and we have done the other hard integrals already, so
$$U_{T_dxy}=\frac{4Ze^2}{4\pi\epsilon_0a}-4\left(-\frac49\right)Dq=\frac{4Ze^2}{4\pi\epsilon_0a}+\frac{16}9Dq$$
and
$$U_{T_dx^2-y^2}=\frac{4Ze^2}{4\pi\epsilon_0a}+6\left(-\frac49\right)Dq=\frac{4Ze^2}{4\pi\epsilon_0a}-\frac83Dq$$
