Smooth sphere data have a harmonic replacement

Translate $a$ to zero. Put $q=\mathrm{\mid }x-y\mathrm{\mid }$, $A=R^2-\mathrm{\mid }x{\mathrm{\mid }}^2$, and $P(x,y)=A{q}^{-n}\mathrm{/}(R{\omega }_{n-1})$, where $\mathrm{\mid }y\mathrm{\mid }=R$. The sphere has finite positive measure ${\omega }_{n-1}{R}^{n-1}$. On compact interior sets $q$ is bounded away from zero; all $x$-derivatives of the kernel times bounded $g$ have a constant integrable majorant. The mean value theorem bounds their difference quotients likewise. Dominated convergence therefore permits differentiation of every order under the sphere integral and proves continuity of those derivatives.

Cartesian differentiation gives $\mathrm{\nabla }{q}^{-n}=-n(x-y)q^{-n-2}$, $\mathrm{\Delta }{q}^{-n}=2n{q}^{-n-2}$ and $\mathrm{\Delta }A=-2n$. The product rule gives $\mathrm{\Delta }(A{q}^{-n})=2n{q}^{-n-2}(A-q^2+2x\cdot (x-y))=2n{q}^{-n-2}(R^2-\mathrm{\mid }y{\mathrm{\mid }}^2)=0$. Thus both $h$ and $I(x)=\int P(x,y)dS(y)$ are smooth harmonic functions.

Rotation invariance of surface measure and the kernel makes $I$ constant on each sphere centered at zero. Its spherical mean equals $I(0)$ by harmonic mean values; since it is already constant on that sphere, $I(x)=I(0)$. At zero, $P(0,y)=1\mathrm{/}({\omega }_{n-1}{R}^{n-1})$, so $I(0)=1$. Also $P>0$ in the ball.

Fix $p\in \mathrm{\partial }{B}_R$. The mass identity gives $h(x)-g(p)=\int P(x,y)(g(y)-g(p))dS(y)$. For a given $\epsilon >0$, take $\delta >0$ so that $\mathrm{\mid }g(y)-g(p)\mathrm{\mid }<\epsilon$ when $\mathrm{\mid }y-p\mathrm{\mid }<\delta$. That part of the integral has absolute value at most $\epsilon$. On the remaining sphere, if $\mathrm{\mid }x-p\mathrm{\mid }<\delta \mathrm{/}2$, then $\mathrm{\mid }x-y\mathrm{\mid }\ge \delta \mathrm{/}2$, and $P(x,y)\le (R^2-\mathrm{\mid }x{\mathrm{\mid }}^2)\mathrm{/}(R{\omega }_{n-1}(\delta \mathrm{/}2{)}^n)$ uniformly tends to zero. The remaining integral is bounded by this supremum times $2\mathrm{\parallel }g{\mathrm{\parallel }}_{\mathrm{\infty }}\mathrm{\mid }\mathrm{\partial }{B}_R\mathrm{\mid }$. Thus $h(x)\to g(p)$, proving continuity on the closed ball.

Any two such harmonic replacements have identical boundary data on a bounded ball and are continuous on its closure, so classical Dirichlet uniqueness makes them equal. Translation restores the stated formula at $a$.