= Solution
One normalization of the <Papkovich–Neuber representation> uses a harmonic vector $\boldsymbol\Phi$ and harmonic scalar $\psi$:
$$
\boxed{\mathbf u=\boldsymbol\Phi-\frac12\nabla(\mathbf x\cdot\boldsymbol\Phi+\psi),\qquad
p=-\mu\nabla\cdot\boldsymbol\Phi,\qquad
\nabla^2\boldsymbol\Phi=0,\quad\nabla^2\psi=0.}
$$
These formulae give $\nabla\cdot\mathbf u=0$ because $\nabla^2(\mathbf x\cdot\boldsymbol\Phi)=2\nabla\cdot\boldsymbol\Phi$, and give $\mu\nabla^2\mathbf u=\nabla p$. Rescaling the potentials gives the other common normalizations of the representation.
A point <force> requires a decaying, rotationally covariant field linear in $\mathbf F$, with no length scale. Thus choose the harmonic monopole $\boldsymbol\Phi=c\mathbf F/r$, $\psi=0$, for $r>0$. Differentiating gives
$$
\mathbf u=\frac c2\left(\frac{\mathbf F}{r}+\frac{(\mathbf F\cdot\mathbf x)\mathbf x}{r^3}\right),\qquad
p=\mu c\frac{\mathbf F\cdot\mathbf x}{r^3}.
$$
The <Newtonian fluid stress tensor> is $\sigma_{ij}=-p\delta_{ij}+\mu(\partial_i u_j+\partial_j u_i)$, which simplifies to $\sigma_{ij}=-3\mu c(\mathbf F\cdot\mathbf x)x_ix_j/r^5$. Its traction resultant on a <sphere>, using an outward normal from the origin, is $-4\pi\mu c\mathbf F$. Momentum balance with the <force> $\mathbf F\delta(\mathbf x)$ therefore fixes $c=1/(4\pi\mu)$. The resulting <Stokeslet> is
$$
\boxed{u_i=\frac1{8\pi\mu}\left(\frac{F_i}{r}+\frac{(\mathbf F\cdot\mathbf x)x_i}{r^3}\right),\qquad
p=\frac{\mathbf F\cdot\mathbf x}{4\pi r^3},\qquad
\sigma_{ij}=-\frac3{4\pi}\frac{(\mathbf F\cdot\mathbf x)x_ix_j}{r^5}.}
$$
An arbitrary constant <pressure> may be added. The stress has purely radial traction on concentric <spheres>, and its integrated traction $-\mathbf F$ supplies the required distributional point-force normalization.
Back to article page