$F={\frac {\mu ep}{4\pi }}\int \int \int {\frac {1}{PQ}}{\frac {d^{2}}{dx^{2}}}{\frac {1}{\rho }}dx\ dy\ dz$;
Now ${\frac {d^{2}}{dx^{2}}}{\frac {1}{\rho }}={\frac {Y_{2}}{\rho ^{3}}}$, where Y_{2} is a surface harmonic of the second order. And when ρ > R,
${\frac {1}{PQ}}={\frac {1}{\rho }}+{\frac {R}{\rho ^{2}}}Q_{1}+{\frac {R^{2}}{\rho ^{3}}}Q_{2}+\dots$;
and when ρ < R,
${\frac {1}{PQ}}={\frac {1}{R}}+{\frac {\rho }{R^{2}}}Q_{1}+{\frac {\rho ^{2}}{R^{3}}}Q_{2}+\dots$;
where Q_{1}, Q_{2}, &c. are zonal harmonics of the first and second orders respectively referred to OP as axis.
Let Y'_{2} denote the value of Y_{2} along OP. Then, since $\int Y_{n}Q_{m}ds$, integrated over a sphere of unit radius, is zero when n and m are different, and ${\frac {4\pi }{2n+1}}Y_{n}^{'}$ when n=m, Y'_{n} being the value of Y_{n} at the pole of Q_{n}, and since there is no electric displacement within the sphere,
$F={\frac {\mu ep}{4\pi }}\times {\frac {2\pi Y_{2}^{'}}{5}}\left\{\int _{R}^{\infty }{\frac {R^{2}}{\rho ^{4}}}d\rho +\int _{a}^{R}{\frac {\rho \ d\rho }{R^{3}}}\right\}$
$={\frac {\mu ep}{5}}Y_{2}^{'}\left({\frac {5}{6R}}{\frac {a^{2}}{2R^{3}}}\right)$,


or, as it is more convenient to write it,
$={\frac {\mu ep}{5}}\left({\frac {5R^{2}}{6}}{\frac {a^{2}}{2}}\right){\frac {d^{2}}{dx^{2}}}{\frac {1}{R}}$.
By symmetry, the corresponding values of G and H are
$G={\frac {\mu ep}{5}}\left({\frac {5R^{2}}{6}}{\frac {a^{2}}{2}}\right){\frac {d^{2}}{dx\ dy}}{\frac {1}{R}}$,
$H={\frac {\mu ep}{5}}\left({\frac {5R^{2}}{6}}{\frac {a^{2}}{2}}\right){\frac {d^{2}}{dx\ dz}}{\frac {1}{R}}$.


These values, however, do not satisfy the condition
${\frac {dF}{dx}}+{\frac {dG}{dy}}+{\frac {dH}{dz}}=0$.
If, however, we add to F the term ${\frac {2\mu ep}{3R}}$, this condition will be satisfied; while, since the term satisfies Laplace's equation, the other conditions will not be affected: thus we have finally