$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\erf}{\operatorname{erf}}$ $\newcommand{\dag}{\dagger}$ $\newcommand{\const}{\mathrm{const}}$ $\newcommand{\arcsinh}{\operatorname{arcsinh}}$
Let $n\ge 3$. Consider spherically symmetric density $f$ (thus depending on $r=(x_1^2+x_2^2+\ldots+x_n^2)^{\frac{1}{2}}$ only). Then it creates a density which is also spherically symmetric.
Assume first that $f=0$ in $B(0,R)$ and consider $u$ in $B(0,R)$. Here $u$ must be harmonic and then due to mean-value theorem $u=u(0)$ is constant (in particular $\nabla u=0$). More precisely \begin{equation} u= -\frac{1}{n-2} \int_R^\infty \rho f(\rho)\,d\rho \label{eq-7.3.1} \end{equation} where we replaced lower limit $0$ by $R$ since $f(\rho)=0$ as $r< R$.
Assume now that $f=0$ in $\{r\ge R\}$ and consider $u$ there. Then $u=r^{2-n} A+B$.
Exercise 1. Prove that if $f$ is spherically symmetric and $f=0$ in $\{r\ge R\}$ then $u=r^{2-n} A+B$ there.
To have potential at infinity equal to $0$ we must take $B=0$. Further plug it into (7.2.2) with $\Omega=B(0,R)$. The left-hand expression becomes the total charge $\int_{B(0,R)} f(x)\,dV$ while the right-hand expression becomes $-(n-2)\sigma_n A$ and therefore \begin{equation} u= -\frac{1}{(n-2)\sigma_n}r^{2-n}\int _{B(0,R)} f\,dx = -\frac{1}{n-2}r^{2-n}\int _0^R \rho^{n-1} f(\rho)\,d\rho. \label{eq-7.3.2} \end{equation} Then in the general case we to calculate $u$ we break $f=f_1+f_2$ where $f_1=0$ for $r\ge R$ and $f_2=0$ for $r\le R$ and then calculate $u$ as $r=R$ and finally set $R=r$: \begin{equation} u= -\frac{1}{n-2} \int_r^\infty \rho f(\rho)\,d\rho-\frac{1}{n-2}r^{2-n}\int _0^r f(\rho)\,d\rho. \label{eq-7.3.3} \end{equation}
In particular, if
Theorem 1. Let $f$ be spherically symmetric. Then
Remark 1. Statement (b) is often expressed by a catch phrase "There is no gravity in the cavity". In particular if $f(r)=\const$ as $R_1\le r\le R_2$ and $f(r)=0$ as $r\le R_1$ or $r\ge R_2$ , then all calculations are easy. Amazingly, the same is true for a shell between two proportional ellipsoids (Ivory theorem).
\begin{multline} u(y)=\int_{\Omega} G^0(x,y)\Delta u(x)\,dV \\ + \int_\Sigma \bigl(-u(x)\frac{\partial G^0}{\partial \nu_x}(x,y)+ G^0(x,y)\frac{\partial u}{\partial \nu}(x)\bigr)\, dS \qquad \tag{7.2.7}\label{eq-7.3.26-7} \end{multline} with \begin{equation} G^0(x,y)=\left\{\begin{aligned} -&\frac{1}{(n-2)\sigma_n}|x-y|^{2-n}&& n\ne 2,\\ -&\frac{1}{4\pi}|x-y|^{-1} &&n=3,\\ &\frac{1}{2\pi}\log |x-y| &&n=2,\\ &\frac{1}{2} |x-y| &&n=1 \end{aligned}\right. \tag{7.2.8}\label{eq-7.3.26-8} \end{equation} where we changed notation $G^0$ instead of $G$ as we will redefine $G$ later.
Let $g(x,y)$ be solution to problem \begin{align} &\Delta_x g(x,y)=0\qquad&&\text{in }\ \Omega,\label{eq-7.3.4}\\[3pt] &g(x,y)=-G^0(x,y) &&\text{as }\ x\in\Sigma\label{eq-7.3.5} \end{align} and condition $u\to 0$ as $|x|\to \infty$ if $\Omega$ is unbounded. In virtue of (\ref{eq-7.3.4}) and (7.2.5) \begin{multline} 0=\int_{\Omega} g(x,y)\Delta u(x)\,dV + \\ \int_\Sigma \bigl(-u(x)\frac{\partial g}{\partial \nu_x}(x,y)+ g(x,y)\frac{\partial u}{\partial \nu}(x)\bigr)\, dS. \qquad \label{eq-7.3.6} \end{multline} Adding to (\ref{eq-7.3.26-7}) we get \begin{multline} u(y)=\int_{\Omega} G(x,y)\Delta u(x)\,dV + \\ \int_\Sigma \bigl(-u(x)\frac{\partial G}{\partial \nu_x}(x,y)+ G(x,y)\frac{\partial u}{\partial \nu}(x)\bigr)\, dS \qquad \label{eq-7.3.7} \end{multline} with \begin{equation} G(x,y):=G^0(x,y)+g(x,y). \label{eq-7.3.8} \end{equation} So far we have not used (\ref{eq-7.3.5}) which is equivalent to $G(x,y)=0$ as $x\in \Sigma$. But then \begin{equation*} u(y)=\int_{\Omega} G(x,y)\Delta u(x)\,dV - \int_\Sigma \frac{\partial G}{\partial \nu_x}(x,y)u(x)\, dS \end{equation*} and therefore \begin{equation} u(y)=\int_{\Omega} G(x,y)f(x)\,dV - \int_\Sigma \frac{\partial G}{\partial \nu_x}(x,y)\phi(x)\, dS \label{eq-7.3.9} \end{equation} is a solution to \begin{align} &\Delta u=f \qquad&&\text{in }\ \Omega,\label{eq-7.3.10}\\[3pt] &u=\phi &&\text{on }\ \Sigma.\label{eq-7.3.11} \end{align}
Similarly if we replace (\ref{eq-7.3.5}) by Robin boundary condition \begin{equation} \bigl(\frac{\partial g}{\partial \nu_x}-\alpha g\bigr)(x,y)= -\bigl(\frac{\partial G^0}{\partial \nu_x}-\alpha G^0\bigr) (x,y) \qquad \text{as }\ x\in\Sigma \label{eq-7.3.12} \end{equation} we get $\bigl(\frac{\partial G}{\partial \nu_x}-\alpha G\bigr)(x,y)=0$ as $x\in \Sigma$ and rewriting (\ref{eq-7.3.7}) as \begin{multline} u(y)=\int_{\Omega} G(x,y)\Delta u(x)\,dV + \\ \int_\Sigma \bigl(-u(x) \bigl[ \frac{\partial G}{\partial \nu_x}- \alpha G \bigr] (x,y)+ G(x,y)\bigl[ \frac{\partial u}{\partial \nu}- \alpha u\bigr]\bigr)\, dS \qquad \label{eq-7.3.13} \end{multline} we get \begin{equation} u(y)=\int_{\Omega} G(x,y)f(x)\,dV + \int_\Sigma G(x,y)\psi(x)\, dS \label{eq-7.3.14} \end{equation} for solution of the (\ref{eq-7.3.10}) with the boundary condition \begin{equation} \frac{\partial u}{\partial \nu}-\alpha u=\psi \qquad \text{on }\ \Sigma. \label{eq-7.3.15} \end{equation}
Finally, if we take $g$ satisfying mixed boundary condition \begin{align} & g(x,y)=-G^0(x,y) && \text{as }\ x\in \Sigma',\label{eq-7.3.16}\\ & \bigl(\frac{\partial g}{\partial \nu_x}-\alpha g\bigr)(x,y)= -\bigl(\frac{\partial G^0}{\partial \nu_x}-\alpha G^0\bigr) (x,y) &&\text{as }\ x\in\Sigma'' \qquad\label{eq-7.3.17} \end{align} we get \begin{multline} u(y)=\int_{\Omega} G(x,y)f(x)\,dV + \int_{\Sigma'} G(x,y)\psi(x)\, dS - \\ \int_{\Sigma''} \frac{\partial G(x,y)}{\partial \nu_x}\phi(x)\, dS\qquad \label{eq-7.3.18} \end{multline} for solution of the (\ref{eq-7.3.10}) with the boundary condition \begin{align} &u=\phi && \text{on }\ \Sigma',\label{eq-7.3.19}\\ &\frac{\partial u}{\partial \nu}-\alpha u=\psi && \text{on }\ \Sigma'', \label{eq-7.3.20} \end{align} where $\Sigma=\Sigma'\cup\Sigma''$ and $\Sigma'\cap\Sigma''=\emptyset$.
Definition 1. $G(x,y)$ (defined for corresponding boundary problem) is called Green's function.
Remark 2. It is similar (by usage) to heat kernel $\frac{1}{2\sqrt{kt}}e^{-\frac{|x-y|^2}{4kt}}$.
One can prove
Theorem 2. \begin{equation}G(x,y)=G(y,x). \label{eq-7.3.21} \end{equation}
Consider now purely Neumann problem in the connected domain. We cannot solve $\Delta g=0$ with the boundary condition $\frac{\partial g}{\partial \nu_x}=-\frac{\partial G^0}{\partial \nu_x}$ as such problem requires one solvability condition to the right-hand expression and boundary value \begin{equation} \int_\Omega f\, dV+ \int _\Sigma \psi\,dS =0 \label{eq-7.3.22} \end{equation} and this condition fails as $\int_\Sigma \frac{\partial G^0}{\partial \nu_x}\, dS\ne 0$ (one can prove it).
To fix it we consider $g$ satisfying \begin{align} &\Delta_x g(x,y)=c\qquad&&\text{in }\ \Omega,\label{eq-7.3.23}\\[3pt] &\frac{\partial g}{\partial \nu_x}(x,y)= -\frac{\partial G^0}{\partial \nu_x}(x,y) &&\text{as }\ x\in\Sigma \label{eq-7.3.24} \end{align} with unknown constant $c$; chosing it correctly one can satisfy solvability condition. Then \begin{equation*} u(y)=\int_{\Omega} G(x,y)f(x)\,dV + \int_{\Sigma} G(x,y)\psi(x)\, dS + c \int _V u\,dx \end{equation*} and therefore \begin{equation} u(y)=\int_{\Omega} G(x,y)f(x)\,dV + \int_{\Sigma} G(x,y)\psi(x)\, dS +C \label{eq-7.3.25} \end{equation} gives us solution if it exists (and it exists under solvability condition (\ref{eq-7.3.24})). This solution is defined up to a constant.