Solution (source code)

= Solution

The <electric polarization> $P$ is electric dipole moment per unit volume. A polarized body acts as the bound volume and <surface charge density>
$$
\rho_{\rm b}=-\nabla\mathbin{\cdot}P,
\qquad \sigma_{\rm b}=P\mathbin{\cdot}n.
$$
Indeed, summing the <scalar potential> of its dipoles and integrating by parts gives
$$
\phi(x)=\frac1{4\pi\epsilon_0}
\int_V P(x')\mathbin{\cdot}\nabla'\frac1{|x-x'|}\,d^3x'
$$
$$
=\frac1{4\pi\epsilon_0}
\left[
\int_V\frac{-\nabla'\mathbin{\cdot}P(x')}{|x-x'|}\,d^3x'
+\int_S\frac{P(x')\mathbin{\cdot}n}{|x-x'|}\,dS'
\right].
$$
This is exactly the potential of $\rho_{\rm b}$ and $\sigma_{\rm b}$.

For the central free charge, <Gauss's law> for the <electric displacement field> gives
$$
D(r)=\frac{q}{4\pi r^2}\hat r,
\qquad E=\frac D\epsilon,
\qquad P=D-\epsilon_0E
=\left(1-\frac{\epsilon_0}{\epsilon}\right)
\frac{q}{4\pi r^2}\hat r.
$$
Thus on the sphere,
$$
\boxed{\sigma_{\rm b}
=\left(1-\frac{\epsilon_0}{\epsilon}\right)
\frac{q}{4\pi R^2}}.
$$