Solution (source code)

= Solution

Away from $x=y$ the Newton kernel is harmonic, while integrating its normal <derivative> over a small sphere gives one, proving $\nabla^2\hat G=\delta$. Green’s second identity is $\int_V(u\nabla^2v-v\nabla^2u)=\int_{\partial V}(u\partial_nv-v\partial_nu)$. For the half-space, images give $G(x,y)=-(4\pi|x-y|)^{-1}+(4\pi|x-y^*|)^{-1}$. Substitution in the boundary formula yields the Poisson kernel $z/[2\pi((x_1-y_1)^2+(x_2-y_2)^2+z^2)^{3/2}]$; polar integration against the Gaussian gives exactly the stated one-dimensional <integral>.