Where does the label ‘Dq’ to denote the field split in coordination compounds come from?
I will expand my previous comment into an answer.
In the $O_h$ case the metal ion is surrounded by $6$ ligand ions, each with charge $q=-Ze$ located at $\langle\pm a,0,0\rangle$, $\langle0, \pm a,0\rangle$, and $\langle0,0,\pm a\rangle$.
The electrostatic potential of a point charge $q$ is
$$V=\frac{kq}R$$
where
$$k=\frac1{4\pi\epsilon_0}$$
is the Coulomb's Law constant and $\epsilon_0$ is the permittivity of free space and $R=\left\|\vec{r_f}-\vec{r_s}\right\|$ is the distance from the source point $\vec{r_s}$ to the field point $\vec{r_f}$.
Working for the moment with the potential $V_1$ due to the charge at $\langle a,0,0\rangle$,
$$R=\left\|\langle x,y,z\rangle-\langle a,0,0\rangle\right\|=\left\|\langle x-a,y,z\rangle\right\|=\sqrt{(x-a)^2+y^2+z^2}=\sqrt{a^2-2ax+x^2+y^2+z^2}=a\sqrt{1-\frac{2x}a+\frac{r^2}{a^2}}$$
Where $r^2=x^2+y^2+z^2$ is the squared distance from the origin. Then
$$V_1=-\frac{Ze}{4\pi\epsilon_0a}\left(1-\frac{2x}a+\frac{r^2}{a^2}\right)^{-1/2}$$
Now we need to expand $(1+x)^{-1/2}$ to fourth order in $x$. The Taylor series yields
$$(1+x)^{-1/2}=1-\frac12x+\frac38x^2-\frac5{16}x^3+\frac{35}{128}x^4+O\left(x^5\right)$$
We could have obtained this result via algebra alone by assuming that $(1+x)^{-1/2}=1+ax+bx^2+cx^3+dx^4+O\left(x^5\right)$, then squaring and multiplying by $(1+x)$ so that
$$\left(1+ax+bx^2+cx^3+dx^4\right)^2(1+x)=1+O\left(x^5\right)$$
and then finding the unknowns $a$, $b$, $c$, and $d$ from the conditions that the coefficients of $x$, $x^2$, $x^3$, and $x^4$ should be zero. Then we have
$$\begin{align}V_1= &amp; -\frac{Ze}{4\pi\epsilon_0a}\left\{1-\frac12\left[-\frac{2x}a+\frac{r^2}{a^2}\right]+\frac38\left[\frac{4x^2}{a^2}-\frac{4xr^2}{a^3}+\frac{r^4}{a^4}\right]\right. \\
&amp; \left.-\frac5{16}\left[-\frac{8x^3}{a^3}+\frac{12x^2r^2}{a^4}\right]+\frac{35}{128}\frac{16x^4}{a^4}\right\}+O\left(\left(\frac ra\right)^5\right)\end{align}$$
This is complicated, but now we simplify by combining the potentials due to all $6$ ligands. We want to find group operations from $O_h$ that map the first ligand to the others, and fortunately the powers of $S_6$, the rotoreflection by angle $\frac{\pi}3$ about the axis through the origin and $\langle1,1,1\rangle$ will do the trick.
Its powers cycle the original ligand point through all the others. We just need to find the mappings of the powers of $x$ by the powers of $S_6$ and sum over all powers of $S_6$:
$$\begin{array}{c|cccccc|c}
f(x,y,z)&amp;Ef&amp;S_6f&amp;S_6^2f&amp;S_6^3f&amp;S_6^4f&amp;S_6^5f&amp;\sum\limits_{i=0}^5S_6^if \\
\hline r^2&amp;r^2&amp;r^2&amp;r^2&amp;r^2&amp;r^2&amp;r^2&amp;6r^2 \\
x&amp;x&amp;-z&amp;y&amp;-x&amp;z&amp;-y&amp;0 \\
x^2&amp;x^2&amp;z^2&amp;y^2&amp;x^2&amp;z^2&amp;y^2&amp;2r^2 \\
x^3&amp;x^3&amp;-z^3&amp;y^3&amp;-x^3&amp;z^3&amp;-y^3&amp;0 \\
x^4&amp;x^4&amp;z^4&amp;y^4&amp;x^4&amp;z^4&amp;y^4&amp;2\left(x^4+y^4+z^4\right)
\end{array}$$
So now
$$\begin{align}V_{O_h}= &amp; \sum\limits_{i=0}^5S_6^iV_1 \\
&amp; =-\frac{Ze}{4\pi\epsilon_0a}\left\{6-\frac12\left[0+\frac{6r^2}{a^2}\right]+\frac38\left[\frac{8r^2}{a^2}-0+\frac{6r^4}{a^4}\right]\right. \\
&amp; \left.-\frac5{16}\left[0+\frac{24r^4}{a^4}\right]+\frac{35}{4}\frac{x^4+y^4+z^4}{a^4}\right\}+O\left(\left(\frac ra\right)^5\right) \\
&amp; =-\frac{Ze}{4\pi\epsilon_0a}\left[6-\frac{21r^4}{4a^4}+\frac{35\left(x^4+y^4+z^4\right)}{4a^4}\right]+O\left(\left(\frac ra\right)^5\right) \\
&amp; =-\frac{6Ze}{4\pi\epsilon_0a}-\frac{35Ze}{16\pi\epsilon_0a^5}\left(x^4+y^4+z^4-\frac35r^4\right)+O\left(\left(\frac ra\right)^5\right)\end{align}$$
All this mathjax is exhausting. I can't finish in one sitting, but at least I have gotten to
$$D=\frac{35Ze}{16\pi\epsilon_0a^5}$$
EDIT: Now that the potential for $O_h$ geometry has been developed, the same will be done for the $T_d$ case where the metal ion is surrounded by $4$ ligand ions, each with charge $q=-Ze$, this time located at
$\left\langle\frac a{\sqrt3},\frac a{\sqrt3},\frac a{\sqrt3}\right\rangle$, $\left\langle\frac a{\sqrt3},-\frac a{\sqrt3},-\frac a{\sqrt3}\right\rangle$, $\left\langle-\frac a{\sqrt3},\frac a{\sqrt3},-\frac a{\sqrt3}\right\rangle$, and $\left\langle-\frac a{\sqrt3},-\frac a{\sqrt3},\frac a{\sqrt3}\right\rangle$.
The distance from the source point at $\vec{r_s}=\left\langle\frac a{\sqrt3},\frac a{\sqrt3},\frac a{\sqrt3}\right\rangle$ to the field point at $\vec{r_f}=\langle x,y,z\rangle$ is
$$\begin{align}R= &amp; \left\|\vec{r_f}-\vec{r_s}\right\|=\left\|\left\langle x-\frac a{\sqrt3},y-\frac a{\sqrt3},z-\frac a{\sqrt3}\right\rangle\right\| \\
&amp; =\sqrt{\left(x-\frac a{\sqrt3}\right)^2+\left(y-\frac a{\sqrt3}\right)^2+\left(z-\frac a{\sqrt3}\right)^2} \\
&amp; =\sqrt{a^2-\frac{2a}{\sqrt3}(x+y+z)+x^2+y^2+z^2} \\
&amp; =a\sqrt{1-\frac{2(x+y+z)}{\sqrt3a}+\frac{r^2}{a^2}}\end{align}$$
Then the potential due to this charge is
$$\begin{align}V_2 &amp; =\frac{kq}R=-\frac{Ze}{4\pi\epsilon_0a}\left(1-\frac{2(x+y+z)}{\sqrt3a}+\frac{r^2}{a^2}\right)^{-1/2} \\
&amp; =-\frac{Ze}{4\pi\epsilon_0a}\left\{1-\frac12\left[-\frac{2(x+y+z)}{\sqrt3a}+\frac{r^2}{a^2}\right]+\frac38\left[\frac4{3a^2}(x+y+z)^2-\frac{4(x+y+z)r^2}{\sqrt3a^3}+\frac{r^4}{a^4}\right]\right. \\
&amp; \left.-\frac5{16}\left[-\frac8{3\sqrt3a^3}(x+y+z)^3+\frac{4(x+y+z)^2r^2}{a^4}\right]+\frac{35}{128}\left[\frac{16}{9a^4}(x+y+z)^4\right]\right\}+O\left(\left(\frac ra\right)^5\right)\end{align}$$
This time the group operations that map the model point to all the others will be the members of the $D_2$ point group: the identity $E$, and the rotations by $180°$ about the x-axis $C_{2x}$, the y-axis $C_{2y}$, and the z-axis $C_{2z}$.
Again a table for combining the contributions from all point charges:
$$\begin{array}{c|cccc|c}
f(x,y,z)&amp;R_1f=Ef&amp;R_2f=C_{2x}f&amp;R_3f=C_{2y}f&amp;R_4f=C_{2z}f&amp;\sum\limits_{i=1}^4R_if \\
\hline r^2&amp;r^2&amp;r^2&amp;r^2&amp;r^2&amp;4r^2 \\
x&amp;x&amp;x&amp;-x&amp;-x&amp;0 \\
xy&amp;xy&amp;-xy&amp;-xy&amp;xy&amp;0 \\
xyz&amp;xyz&amp;xyz&amp;xyz&amp;xyz&amp;4xyz
\end{array}$$
And the reader can readily verify that also
$$\sum\limits_{i=1}^4R_iy=\sum\limits_{i=1}^4R_iz=\sum\limits_{i=1}^4R_ixz=\sum\limits_{i=1}^4R_iyz=0$$
So the rule is that terms in the expression for $V_2$ with odd powers of any of the coordinates cancels except for $xyz$, which, along with all terms with only even powers of coordinates, adds.
Expanding powers of $(x+y+z)$ and summing over the four sources,
$$\begin{align}V_{T_d} &amp; =\sum\limits_{i=1}^4R_iV_2=-\frac{Ze}{4\pi\epsilon_0a}\left\{4-\frac12\left[0+\frac{4r^2}{a^2}\right]+\frac38\left[\frac4{3a^2}(4r^2)-0+\frac{4r^4}{a^4}\right]\right. \\
&amp; \left.-\frac5{16}\left[-\frac8{3\sqrt3a^3}(6)(4xyz)+\frac{4(4r^2)r^2}{a^4}\right]\right. \\
&amp; \left.+\frac{35}{128}\frac{16}{9a^4}\left[4\left(x^4+y^4+z^4\right)+6(4)\left(x^2y^2+x^2z^2+y^2z^2\right)\right]\right\}+O\left(\left(\frac ra\right)^5\right) \\
&amp; = -\frac{Ze}{4\pi\epsilon_0a}\left\{4+\frac{3r^4}{2a^4}+\frac{20}{\sqrt3a^3}xyz-\frac{5r^4}{a^4}\right. \\
&amp; \left.+\frac{35}{72a^4}\left[4\left(x^4+y^4+z^4\right)+12\left(r^4-x^4-y^4-z^4\right)\right]\right\}+O\left(\left(\frac ra\right)^5\right) \\
&amp; = -\frac{Ze}{4\pi\epsilon_0a}\left[4+\frac{20}{\sqrt3a^3}xyz-\frac{35}{9a^4}\left(x^4+y^4+z^5-\frac35r^4\right)\right]+O\left(\left(\frac ra\right)^5\right) \\
&amp; = -\frac{4Ze}{4\pi\epsilon_0a}-\frac{5Ze}{\sqrt3\pi\epsilon_0a^4}xyz+\frac{35Ze}{36\pi\epsilon_0a^5}\left(x^4+y^4+z^4-\frac35r^4\right)+O\left(\left(\frac ra\right)^5\right)\end{align}$$
On comparison with $V_{O_h}$, the first term is the potential due to the net charge of all the ligands uniformly spread over the surface of a sphere of radius $a$.
The second term has no counterpart in $V_{O_h}$ and causes no splitting but does create mixing between the $d_{xy}$ and $p_z$ orbitals and would lower the energy of the former and raise the energy of the latter.
I don't know how significant that effect is, though. The splitting term is $\frac49D\left(x^4+y^4+z^4-\frac35r^4\right)$ and is in fact $-\frac49$ the value of the splitting term for $O_h$ geometry.
Still not done, but it's getting late.  

EDIT:To get the $q$ part of $Dq$ we are going to have to do some integrals. I am going to use the notation
$$\oint f(x,y,z)d^2\Omega$$ to represent the integral of $f(x,y,z)$ over the surface of the unit sphere.
One integral we know right off the bat is
$$\oint 1d^2\Omega=4\pi$$
using the mensuration formula from geometry for the area of a sphere.
Actually, if we didn't know this formula, it would be no problem since if the answer were an unknown, say $\xi$, all the $\xi$'s would cancel out in the final results.
If the integrand is a monomial with an odd power of $z$, then
$$\begin{align}\oint x^ay^bz^{2c+1}d^2\Omega &amp; =\oint \sigma_{hz}\left(x^ay^bz^{2c+1}\right)d^2\Omega=(-1)^{2c+1}\oint x^ay^bz^{2c+1}d^2\Omega \\
&amp; =\frac{1+(-1)^{2c+1}}2\oint x^ay^bz^{2c+1}d^2\Omega=0\end{align}$$
Where $\sigma_{hz}$ is a reflection though the $xy$-plane. Clearly any proper rotation of the integrand about an axis through the origin doesn't affect the value of the integral because it just corresponds to viewing the same function from a different angle.
Also inversion of the integrand through the origin doesn't change the integral and any improper rotation including reflection can be written as inversion times a proper rotation. Thus it can be seen that if
$$\oint x^ay^bz^cd^2\Omega \ne 0$$
Then $a$, $b$, and $c$ must all be even. Also we can rotate coordinates so that we only need compute
$$\oint x^{2a}y^{2b}z^{2c}d^2\Omega$$
for $a \ge b \ge c$.
We can calculate
$$\oint r^2d^2\Omega=\oint(x^2+y^2+z^2)d^2\Omega=3\oint x^2d^2\Omega=4\pi$$
where we have rotated individual terms and then used $r^2=1$ on the surface of the unit sphere. Now,
$$\oint x^2d^2\Omega=\frac{4\pi}3=\oint x^2r^2d^2\Omega=\oint \left(x^4+x^2y^2+x^2z^2\right)d^2\Omega=\oint x^4d^2\Omega+2\oint x^2y^2d^2\Omega$$
Rotating by $\frac{\pi}4$ about the $z$-axis,
$$\begin{align}\oint x^4d^2\Omega &amp; =\oint \left(\frac{x+y}{\sqrt2}\right)^4d^2\Omega=\frac14\oint\left(x^4+4x^3y+6x^2y^2+4xy^3+y^4\right)d^2\Omega \\
&amp; =\frac12\oint x^4d^2\Omega+\frac32\oint x^2y^2d^2\Omega\end{align}$$
This gives us two eqautions in two unknowns and we can solve using Excel or Matlab or whatever to get
$$\oint x^4d^2\Omega=\frac{4\pi}5\text{ and }\oint x^2y^2d^2\Omega=\frac{4\pi}{15}$$
Again we have
$$\oint x^2d^2\Omega=\frac{4\pi}3=\oint x^2r^4d^2\Omega=\oint x^6d^2\Omega+6\oint x^4y^2d^2\Omega+2\oint x^2y^2z^2d^2\Omega$$
$$\oint x^4d^2\Omega=\frac{4\pi}5=\oint x^4r^2d^2\Omega=\oint x^6d^2\Omega+2\oint x^4y^2d^2\Omega$$
$$\oint x^6d^2\Omega=\oint \left(\frac{x+y}{\sqrt2}\right)^6d^2\Omega=\frac14\oint x^6d^2\Omega+\frac{15}4\oint x^4y^2d^2\Omega$$
Solving the last $3$ equations in $3$ unknowns we get
$$\oint x^6d^2\Omega=\frac{4\pi}7\text{, }\oint x^4y^2d^2\Omega=\frac{4\pi}{35}\text{, and }\oint x^2y^2z^2d^2\Omega=\frac{4\pi}{105}$$
Almost there :)
$$\begin{align}\oint x^2d^2\Omega &amp; =\frac{4\pi}3=\oint x^2r^6d^2\Omega \\
&amp; =\oint x^8d^2\Omega+8\oint x^6y^2d^2\Omega+6\oint x^4y^4d^2\Omega+12\oint x^4y^2z^2d^2\Omega\end{align}$$
$$\begin{align}\oint x^4d^2\Omega &amp; =\frac{4\pi}5=\oint x^4r^4d^2\Omega \\
&amp; =\oint x^8d^2\Omega+4\oint x^6y^2d^2\Omega+2\oint x^4y^4d^2\Omega+2\oint x^4y^2z^2d^2\Omega\end{align}$$
$$\oint x^6d^2\Omega=\frac{4\pi}7=\oint x^6r^2d^2\Omega=\oint x^8d^2\Omega+2\oint x^6y^2d^2\Omega$$
$$\oint x^8d^2\Omega=\oint \left(\frac{x+y}{\sqrt2}\right)^8d^2\Omega=\frac18\oint x^8d^2\Omega+\frac72\oint x^6y^2d^2\Omega+\frac{35}8\oint x^4y^4d^2\Omega$$
Solving these $4$ equations in $4$ unknowns,
$$\oint x^8d^2\Omega=\frac{4\pi}{9}\text{, }\oint x^6y^2d^2\Omega=\frac{4\pi}{63}\text{, }\oint x^4y^4d^2\Omega=\frac{4\pi}{105}\text{, and }\oint x^4y^2z^2d^2\Omega=\frac{4\pi}{315}$$
With these integrals in hand we want to normalize the polynomials $xy$ and $x^2-y^2$.
$$\oint (xy)^2d^2\Omega=\oint x^2y^2d^2\Omega=\frac{4\pi}{15}$$
So a normalized spherical harmonic that transforms as the $xy$ partner of the $T_{2g}$ representation of $O_h$ and as the $z$ partner the $T_2$ representation of $T_d$ is
$$Y_{xy}^2(x,y,z)=\left(\frac{15}{4\pi}\right)^{1/2}\frac{xy}{r^2}$$
And we have
$$\psi_{xy}(x,y,z)=Y_{xy}^2(x,y,z)f_{xy}(r)$$
To find the electrostatic potential energy of an electron with such a wave function in the $O_h$ potential we integrate over its charge density to get
$$\begin{align}U_{O_h,xy} &amp; =\int V_{O_h}(x,y,z)\rho_{xy}(x,y,z)d^3volume \\
&amp; =\int\left[-\frac{6Ze}{4\pi\epsilon_0a}-D\left(x^4+y^4+z^4-\frac35r^4\right)\right]\times \\
&amp; (-e)\left[\left(\frac{15}{4\pi}\right)^{1/2}\frac{xy}{r^2}f_{xy}(r)\right]^*\left[\left(\frac{15}{4\pi}\right)^{1/2}\frac{xy}{r^2}f_{xy}(r)\right]d^3volume \\
&amp; =\frac{6Ze^2}{4\pi\epsilon_0a}\oint\left|\left(\frac{15}{4\pi}\right)^{1/2}\frac{xy}{r^2}\right|^2d^2\Omega\int_0^{\infty}\left|f_{xy}(r)\right|^2r^2dr \\
&amp; +De\left(\frac{15}{4\pi}\right)\oint \left(\frac{xy}{r^2}\right)^2\frac{x^4+y^4+z^4-\frac35r^4}{r^4}\int_0^{\infty}r^4\left|f_{xy}(r)\right|^2r^2dr \\
&amp; =\frac{6Ze^2}{4\pi\epsilon_0a}(1)(1)+\frac{15De}{4\pi}\left(\frac{4\pi}{63}+\frac{4\pi}{63}+\frac{4\pi}{315}-\frac35\frac{4\pi}{15}\right)\left\langle r^4\right\rangle \\
&amp; =\frac{6Ze^2}{4\pi\epsilon_0a}-\frac{8De}{105}\left\langle r^4\right\rangle\end{align}$$
Where we have used the separate normalizations of the angular part of the wave function $Y_{xy}^2(x,y,z)$ and the radial part $f_{xy}(r)$ and the average value of $r^4$ for the radial distribution, $\left\langle r^4\right\rangle$.
Now we calculate
$$\oint \left(x^2-y^2\right)^2d^2\Omega=2\oint x^4d^2\Omega-2\oint x^2y^2d^2\Omega=2\left(\frac{4\pi}{5}\right)-2\left(\frac{4\pi}{15}\right)=\frac{16\pi}{15}$$
So a normalized spherical harmonic that transforms as the $x^2-y^2$ partner of the $E_g$ representation of $O_h$ and the $E$ representation of $T_d$ is
$$Y_{x^2-y^2}^2(x,y,z)=\left(\frac{15}{16\pi}\right)^{1/2}\left(x^2-y^2\right)$$
We could also have obtained this normalization by rotating $Y_{xy}^2(x,y,z)$ by $-90°$ about the $z$-axis.
Now we can get the electrostatic potential energy of an electron with wave function
$$\psi_{x^2-y^2}(x,y,z)=Y_{x^2-y^2}^2(x,y,z)f_{x^2-y^2}(r)$$
in the $O_h$ potential as
$$\begin{align}U_{O_h,x^2-y^2} &amp; =\int V_{O_h}(x,y,z)\rho_{x^2-y^2}(x,y,z)d^3volume \\
&amp; =\int\left[-\frac{6Ze}{4\pi\epsilon_0a}-D\left(x^4+y^4+z^4-\frac35r^4\right)\right]\times \\
&amp; (-e)\left[\left(\frac{15}{16\pi}\right)^{1/2}\frac{x^2-y^2}{r^2}f_{x^2-y^2}(r)\right]^*\left[\left(\frac{15}{16\pi}\right)^{1/2}\frac{x^2-y^2}{r^2}f_{x^2-y^2}(r)\right]d^3volume \\
&amp; =\frac{6Ze^2}{4\pi\epsilon_0a}\oint\left|\left(\frac{15}{16\pi}\right)^{1/2}\frac{x^2-y^2}{r^2}\right|^2d^2\Omega\int_0^{\infty}\left|f_{x^2-y^2}(r)\right|^2r^2dr \\
&amp; +De\left(\frac{15}{16\pi}\right)\oint \left(\frac{x^2-y^2}{r^2}\right)^2\frac{x^4+y^4+z^4-\frac35r^4}{r^4}\int_0^{\infty}r^4\left|f_{x^2-y^2}(r)\right|^2r^2dr \\
&amp; =\frac{6Ze^2}{4\pi\epsilon_0a}(1)(1)+\frac{15De}{16\pi}\left((2)\frac{4\pi}{9}+(4)\frac{4\pi}{105}-(4)\frac{4\pi}{63}-(2)\frac{4\pi}{315}-\frac65\frac{4\pi}{5}+\frac65\frac{4\pi}{15}\right)\left\langle r^4\right\rangle \\
&amp; =\frac{6Ze^2}{4\pi\epsilon_0a}+\frac{4De}{35}\left\langle r^4\right\rangle\end{align}$$
Where we have assumed the average value of $r^4$ is the same as for $\psi_{xy}(x,y,z)$ because they both started out as $d$ orbitals and we haven't taken into account higher-order effects.
Since there are $3$ $T_{2g}$ orbitals and $2$ $E_g$ orbitals, the average potential energy over all orbitals is
$$U_{O_h,ave}=\frac15\left[3\left(\frac{6Ze^2}{4\pi\epsilon_0a}-\frac{8De}{105}\left\langle r^4\right\rangle\right)+2\left(\frac{6Ze^2}{4\pi\epsilon_0a}+\frac{4De}{35}\left\langle r^4\right\rangle\right)\right]=\frac{6Ze^2}{4\pi\epsilon_0a}$$
The splitting in potential energies of the orbitals is
$$\begin{align}U_{O_h,x^2-y^2}-U_{O_h,xy} &amp; =\left(\frac{6Ze^2}{4\pi\epsilon_0a}+\frac{4De}{35}\left\langle r^4\right\rangle\right)-\left(\frac{6Ze^2}{4\pi\epsilon_0a}-\frac{8De}{105}\left\langle r^4\right\rangle\right) \\
&amp; =\frac{4De}{21}\left\langle r^4\right\rangle\equiv10Dq\end{align}$$
Then we get
$$q=\frac{2e}{105}\left\langle r^4\right\rangle$$
and so
$$U_{O_h,x^2-y^2}=\left(\frac{6Ze^2}{4\pi\epsilon_0a}\right)+6Dq$$and$$U_{O_h,xy}=\left(\frac{6Ze^2}{4\pi\epsilon_0a}\right)-4Dq$$
This derivation was a little bloody because I chose to do everything with only algebra.
There is a much more elegant exposition possible that requires higher level math.
Is it worth adding to the end of this answer or posting as another answer?