You can solve the divergence equation in spherical coordinates to produce a ## \rho_p ## that is constant everywhere, where ## -\nabla \cdot P=\rho_p ##, but polarization still creates an electrically neutral system, so that the surface polarization charge density of ## \sigma_p=P \cdot \hat{n} ## will be such that it neutralizes any positive ## \rho_p ## in the interior. Solving the divergence equation, you get ## P_r(r)=-(\rho_p/3)r ##. ## \\ ## ## \int \rho_p dv=+(4/3) \pi R^3 \rho_p. ## ## \\ ## Now ## P_r=-(\rho_p/3)R ## at ## r=R ##, so that ## \sigma_p=-(\rho_p/3)R ##. ## \\ ## ## \int \sigma_p dA=-(\rho_p/3)R(4 \pi R^2) ## precisely neutralizing ## \int \rho_p dv ##.