Problems using Green's Second Theorem to derive the Kirchhoff-Helmholtz equation Edit - I've tried to answer my own questions but I think there could be some flaws. Can someone help me please. I think I'm pretty close to it just need someone to give me the extra push or give the answer a counterexample. I've been struggling on this for a long time now. I've also seen that for some reason people are taking the volume integral without the point of the receiver, which does indeed follow through. However this relationship needs the volume integral of the receiver to to be equal with the volume integral formed in the question.
I'm really struggling with one line when it comes through the derivation of the Kirchhoff-Helmholtz equation used in Theoretical Acoustics and Wave Scattering from Rough Surfaces. To me the missing line seems pretty critical in showing the required rigor needed to convince oneself that this holds. 
Effectively you begin with the Inhomogenous Helmholtz equation and it's coupled Green's function solution;
$$ (\Delta + k^2)u = -f(\mathbf{r}),\quad (\Delta + k^2)G(\mathbf{r,r_0}) = -\delta(\mathbf{r-r_0}).  $$
Multiply the first equation by $G(\mathbf{r,r_0})$ and the second equation by $u$ and then take a volume integral with respect to $x_0,y_0,z_0$ (which I assume is the co-ordinates of $\delta$ and $G$?)  of the whole volume occupied by the medium (as well as swapping $r$ and $r_0$ due to reciprocity with Green's and delta) which simply gives;
$$ \int\int\int [G(\mathbf{r,r_0})\Delta u(\mathbf{r_0})-u(\mathbf{r_0})\Delta G(\mathbf{r,r_0})] dv_0 = \int\int\int u(\mathbf{r_0})\delta(\mathbf{r-r_0}) dv_0 - \int\int\int f(\mathbf{r_0})G(\mathbf{r,r_0}) dv_0.$$
The above equation is just;
$$ u(\mathbf{r}) = \int\int\int f(\mathbf{r_0})G(\mathbf{r,r_0}) dv_0 +  \int\int\int [G(\mathbf{r,r_0})\Delta u(\mathbf{r_0})-u(\mathbf{r_0})\Delta G(\mathbf{r,r_0})] dv_0 . $$
From here this is where the problems begin. The next line in textbooks will just apply Green's second theorem to the second term on the right which turns it into a surface integral. However I think this is non-trivial. Surely there is something happening at the point $\mathbf{r_0 = r} $ which contradicts the fact that everywhere in the volume must be non-singular and differentiable (the conditions to apply Green's Second). 
The natural next step I have been trying is to modify the volume so there is a ball of radius tending to $\mathbf{r}$ or splitting the surface integrals up to have a surface integral of that ball as well as the other surface. But then there is a surface integral centered at a pole which should be non-zero using Residue Theorem (Edit - trying to use the residue theorem has been debunked by the "answer"). However the next step in text books is just:
$$ u(\mathbf{r}) = \int\int\int f(\mathbf{r_0})G(\mathbf{r,r_0}) dv_0 +  \int\int\int [G(\mathbf{r,r_0})\frac{\partial}{\partial n} u(\mathbf{r_0})-u(\mathbf{r_0})\frac{\partial}{\partial n} G(\mathbf{r,r_0})] dS_0  $$ 
I feel like there is something critical I am missing. But for the life of me I just cannot see what is happening. Any help would be of immense help! Thank you in advance!
Edit: Thought I had a lead with Sommerfeld's Optics Vol 4 but that also in no way assists the line jump.
Edit 2: possibly there's an approach by making two volumes as shown in the "answer" and then showing the ball volume tends to zero and the surface of that ball also tends to 0 as the radius tends to zero. I can't see how to do this. Any help would be greatly appreciated.
Edit 3: the tending to zero volume argument seems to be the way forward. But I can't evaluate the integral. Intuitively I feel like there should be a volume of sphere multiplied by something but I can't find any theory to assist that. 
Edit 4: I think there is something in the picture below. The text is in Russian so I am not exactly sure how the variables have been changed for context we are looking at the bottom of the first page and the start of the next page. PLEASE can someone help me with some understanding. Thank you.
Edit 5: so even the Russian text drops the volume integral of f multiplied by the greens theorem. But I can't find anything which allows that to happen. The only way I saw that is by considering the scattered field. But that makes the methodology in the Russian text questionable. 
Edit 6: I heard there is something along the lines of using Green's second theorem specifically for distributions. I haven't found anything for it though.
 A: Given that you've written
$$
(\nabla^2 + k^2)G(\mathbf{r},\mathbf{r}_0) = -\delta(\mathbf{r}-\mathbf{r}_0),
$$
you must be using distributional derivatives. Otherwise, how else did you end up with a distribution from your differential operator? The distributional derivative has its own divergence theorem, and the functions seen here satisfy its conditions, so passing to the boundary works fine.
If you really want to do it using standard derivatives, you're going to have to work with the equation in integral form. It's not pretty, but it goes like this...
Let $V$ be the volume we're solving the equation in. For all volumes $W\subseteq V$ whose boundary does not contain $\mathbf{r}_0$, $G$ satisfies the following equation:
$$
\oint_{\partial W} \boldsymbol\nabla G(\mathbf{r},\mathbf{r}_0)\cdot d^2\mathbf{A} + \int_{W} k^2 G(\mathbf{r},\mathbf{r}_0)d^3\mathbf{r} = -\chi_{W}(\mathbf{r}_0),
$$
where $\chi_W$ is the indicator function of the volume $W$. Note that $G$ is differentiable everywhere on $\partial W$, so this is well-defined. We're going to use this to show that the equation you derived for $u(\mathbf{r})$ holds when integrated over any $W\subseteq V$, and thus holds at every point of $V$. We start by multiplying the equation for $G$ by  $u(\mathbf{r}_0)$ and integrating $\mathbf{r}_0$ over the whole volume $V$:
$$
\int_V u(\mathbf{r}_0)\left[\oint_{\partial W} \boldsymbol\nabla G(\mathbf{r},\mathbf{r}_0)\cdot d^2\mathbf{A}\right]d^3\,\mathbf{r}_0 +  \int_V\int_{W} k^2 u(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0)d^3\mathbf{r}d^3\mathbf{r}_0 = -\int_Wu(\mathbf{r}_0)d^3\mathbf{r}_0.
$$
At this point we use Helmholtz's equation to write $k^2 u(\mathbf{r}_0)$ in terms of $f$ and $\nabla^2 u(\mathbf{r}_0)$, rename the dummy variable in the $u$ integral, and rearrange a little:
\begin{multline}
\int_V  u(\mathbf{r}_0)\left[\oint_{\partial W}\boldsymbol\nabla G(\mathbf{r},\mathbf{r}_0)\cdot d^2\mathbf{A}\right]d^3\,\mathbf{r}_0 \\ = \int_W\left[\int_V\nabla_0^2 u(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0)d^3\mathbf{r}_0+ \int_V f(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0)d^3\mathbf{r}_0-u(\mathbf{r})\right]d^3\mathbf{r}.
\end{multline}
Now we need to turn that surface integral on $\partial W$ into a volume integral on $W$. We can't appeal to the divergence theorem directly because we can't differentiate $G$ over the whole volume. Instead, we use the fact that $\boldsymbol\nabla G(\mathbf{r},\mathbf{r}_0) = \boldsymbol\nabla_0 G(\mathbf{r},\mathbf{r}_0)$ to apply the divergence theorem to the $\mathbf{r}_0$ integral and get rid of the $G$ gradients:
\begin{multline}
\int_Vu(\mathbf{r}_0) \left[\oint_{\partial W} \boldsymbol\nabla G(\mathbf{r},\mathbf{r}_0)\cdot d^2\mathbf{A}\right]d^3\,\mathbf{r}_0 = \int_V u(\mathbf{r}_0)\left[\oint_{\partial W} \boldsymbol\nabla_0 G(\mathbf{r},\mathbf{r}_0)\cdot d^2\mathbf{A}\right]d^3\,\mathbf{r}_0 \\= \int_V u(\mathbf{r}_0)\boldsymbol\nabla_0\cdot\left[\oint_{\partial W}  G(\mathbf{r},\mathbf{r}_0) d^2\mathbf{A}\right]d^3\,\mathbf{r}_0 \\= \oint_{\partial V} \oint_{\partial W}  u(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0) d^2\mathbf{A}\cdot d^2\mathbf{A}_0 - \int_V \boldsymbol\nabla_0u(\mathbf{r}_0)\cdot\left[\oint_{\partial W}  G(\mathbf{r},\mathbf{r}_0) d^2\mathbf{A}\right]d^3\,\mathbf{r}_0.
\end{multline}
Now with no more $G$ gradients, we can interchange the order of integration and apply the divergence theorem to get rid of all the surface integrals in $W$:
\begin{multline}
= \int_W \boldsymbol\nabla\cdot\left[\oint_{\partial V}u(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0) d^2\mathbf{A}_0\right]d^3\mathbf{r} - \int_{W}\boldsymbol \nabla\cdot\left[\int_V G(\mathbf{r},\mathbf{r}_0)\boldsymbol\nabla_0 u(\mathbf{r}_0)d^3\mathbf{r}_0 \right]d^3\mathbf{r}
\\ = \int_W\left[\oint_{\partial V} u(\mathbf{r}_0)\boldsymbol \nabla_0 G(\mathbf{r},\mathbf{r_0})\cdot d^2\mathbf{A}_0-\boldsymbol \nabla\cdot\int_V G(\mathbf{r},\mathbf{r}_0)\boldsymbol\nabla_0 u(\mathbf{r}_0)d^3\mathbf{r}_0 \right]d^3\mathbf{r}.
\end{multline}
Note that differentiating under the integral sign is justified here because $G$ is differentiable at all points of $\partial V$, and I have again used $\boldsymbol\nabla G(\mathbf{r},\mathbf{r}_0) = \boldsymbol\nabla_0 G(\mathbf{r},\mathbf{r}_0)$. So we have
\begin{multline}
\int_W\left[\oint_{\partial V} u(\mathbf{r}_0)\boldsymbol \nabla_0 G(\mathbf{r},\mathbf{r_0})\cdot d^2\mathbf{A}_0-\boldsymbol \nabla\cdot\int_V G(\mathbf{r},\mathbf{r}_0)\boldsymbol\nabla_0 u(\mathbf{r}_0)d^3\mathbf{r}_0\right]d^3\mathbf{r}\\ = \int_W\left[\int_V\nabla_0^2 u(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0)d^3\mathbf{r_0}+\int_V f(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0)d^3\mathbf{r}_0-u(\mathbf{r})\right]d^3\mathbf{r}.
\end{multline}
Since this is true for all volumes $W\subseteq V$, it must be true at all points in $V$, and thus we have for all $\mathbf{r}\in V$,
\begin{multline}
u(\mathbf{r}) = \int_V f(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0)d^3\mathbf{r}_0 - \oint_{\partial V} u(\mathbf{r}_0)\boldsymbol \nabla_0 G(\mathbf{r},\mathbf{r_0})\cdot d^2\mathbf{A}_0\\+\boldsymbol \nabla\cdot\int_V G(\mathbf{r},\mathbf{r}_0)\boldsymbol\nabla_0 u(\mathbf{r}_0)d^3\mathbf{r}_0+\int_V\nabla_0^2 u(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0)d^3\mathbf{r}_0.
\end{multline}
OK, at this point I can't remember the proper way to show that that last line is the other half of the boundary term. If you allow me the leeway of differentiating under the integral sign even though $G$ is not differentiable at $\mathbf{r} = \mathbf{r}_0$, then it follows from
\begin{multline}
\boldsymbol \nabla\cdot\int_V G(\mathbf{r},\mathbf{r}_0)\boldsymbol\nabla_0 u(\mathbf{r}_0)d^3\mathbf{r}_0+\int_V\nabla_0^2 u(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0)d^3\mathbf{r}_0\\ = 
\int_V \left[\boldsymbol\nabla_0 G(\mathbf{r},\mathbf{r}_0)\cdot\boldsymbol\nabla_0 u(\mathbf{r}_0)d^3\mathbf{r}_0+\nabla_0^2 u(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0)\right]d^3\mathbf{r}_0 \\= \int_V \boldsymbol\nabla_0\cdot\left[G(\mathbf{r},\mathbf{r}_0)\boldsymbol\nabla_0 u(\mathbf{r}_0)\right]d^3\mathbf{r}_0 = \oint_{\partial V} G(\mathbf{r},\mathbf{r}_0)\boldsymbol\nabla_0 u(\mathbf{r}_0)\cdot d^2\mathbf{A}_0
\end{multline}
which gives
$$
u(\mathbf{r}) = \int_V f(\mathbf{r}_0)G(\mathbf{r},\mathbf{r}_0)d^3\mathbf{r}_0 + \oint_{\partial V} \left[ G(\mathbf{r},\mathbf{r}_0)\boldsymbol\nabla_0 u(\mathbf{r}_0)-u(\mathbf{r}_0)\boldsymbol \nabla_0 G(\mathbf{r},\mathbf{r_0})\right]\cdot d^2\mathbf{A}_0
$$
as desired.
A: It's a function in 3 variables not 2, I can't see how apply Residue Theorem from complex variable.
The green function for Helmholtz equation in $\mathbb{R}³$ should be 
$$ G(x,y) = \frac{e^{ik|x-y|}}{4\pi|x-y|}$$
For find the green function. Just solve de Helmholtz homogeneous equation $\Delta G + k²G = -\delta $ using separation of variables and solve de Bessel ODE which appears when we apply that technique.
Now note 
$$ \nabla_{y} G(x,y) = \left(\frac{1}{|x-y|}-ik\right) \frac{e^{ik|x-y|}}{4\pi|x-y|}n(y),$$
where $n$ is the normal vector of surface $\partial D$ directed to the exterior.
Consider $x\in D$ and the sphere $S(x,\rho)=\{x \in \mathbb{R}³:|x-y|=\rho\}$ circumscribed in $D$ and $D_\rho= \{y \in D :|x-y|>\rho\}$ in $S(x,\rho)\cup \partial D$ by Green's formula results:
$$\int_{\partial D \cup S(x,\rho)} \partial_nu(y)G(x,y)-\partial_nG(x,y)u(y) ds = 
\int_{D_p}(\Delta u(y)+k²u(y))G(x,y)dV $$
Then by a mean value theorem you can prove that
$$\lim_{\rho \to 0} \int_{\partial D \cup S(x,\rho)} \partial_nu(y)G(x,y)-\partial_nG(x,y)u(y) ds =u(x).$$
Note here we assume that de function $u$ is $C²(\bar{D})$ for have normal derivative continuous.
A: Okay so I've mostly figured it out. First we split the volume integral into the volume without a ball of radius $\epsilon$ and a volume of a ball centered at $x_0$ of radius $\epsilon$ as follows;
$$ \int_{V_0}G(r,r_0)\Delta u(r_0) - u(r_0)\Delta G(r,r_0)dV_0 = \int_{V_\epsilon}G(r,r_0)\Delta u(r_0) - u(r_0)\Delta G(r,r_0) dV_\epsilon + \int_{V'}G(r,r_0)\Delta u(r_0) - u(r_0)\Delta G(r,r_0) dV'$$ 
Now the volume integral V' contains no discontinuities and therefore can be transformed into three surface integrals. One will be the boundary of the exterior volume (which tends to 0 by Sommerfield radiation condition so we can omit that - wave decays tending to infinity so $u$ and it's normal derivative will be 0). As well as the surface that we are scattering and the surface integral of the ball;
$$ \int_{V'}G(r,r_0)\Delta u(r_0) - u(r_0)\Delta G(r,r_0) dV'= \int_SG(r,r_0)\frac{\partial}{\partial n }u(r_0) - u(r_0)\frac{\partial}{\partial n}G(r,r_0)dS_0 + \int_{S_{\epsilon}} G(r,r_0)\frac{\partial}{\partial n }u(r_0) - u(r_0)\frac{\partial}{\partial n}G(r,r_0)dS_\epsilon  $$
So, now we have two terms we need to cancel out. We do this by tending the radius $\epsilon$ to zero therefore $G(r,r_0) = \frac{e^{ik\epsilon}}{\epsilon}$ in $V_\epsilon$ and on $S_\epsilon$
Firstly, the $V_\epsilon$ term:
$$lim_{\epsilon\rightarrow 0}\int_{V_\epsilon} G(r,r_0)\Delta u(r_0) - u(r_0)\Delta G(r,r_0)dV_\epsilon = lim_{\epsilon\rightarrow 0}\int_{V_\epsilon} G(r,r_0)\Delta u(r_0) + u(r_0)k^2G(r,r_0)dV_\epsilon $$
This is because we are not at the point $r_0 = r$ therefore the delta function is zero. This leaves:
$$ lim_{\epsilon\rightarrow 0}\int_{V_\epsilon} G(r,r_0)f(r_0)dV_\epsilon  $$ as this is the solution to the Helmholtz equation (there could potentially be a problem here but it doesn't seem like it.
Now all is left to do is to parameterise $dV_\epsilon = \epsilon^2 \sin(\phi) d\epsilon d\theta d\phi$
$$ lim_{\epsilon\rightarrow 0}\int_{V_\epsilon} G(r,r_0)f(r_0)\epsilon^2 \sin(\phi) d\epsilon d\theta d\phi $$
$$ lim_{\epsilon\rightarrow 0}\int_{V_\epsilon} \frac{e^{ik\epsilon}}{\epsilon}f(r_0)\epsilon^2 \sin(\phi) d\epsilon d\theta d\phi $$
Now G is just the free field Green's function so has $\frac{1}{\epsilon}$ term which just leaves the volume integral tending to zero the same way.
The only issue is;
$$\lim_{\epsilon\rightarrow 0}\int_{S_{\epsilon}} G(r,r_0)\frac{\partial}{\partial n }u(r_0) - u(r_0)\frac{\partial}{\partial n}G(r,r_0)dS_\epsilon \rightarrow 0 $$
I feel like this could be done by the same logic. Please can someone pick up this one. Possibly by just substituting the value of $dS_\epsilon$ by the same way in equation 9 of http://math.mit.edu/~jorloff/suppnotes/suppnotes02/v9.pdf so
$$ \lim_{\epsilon\rightarrow 0}\int_{S_{\epsilon}} G(r,r_0)\frac{\partial}{\partial n }u(r_0) - u(r_0)\frac{\partial}{\partial n}G(r,r_0)\epsilon^2sin\phi d\phi d\theta $$ So clearly 
$$ \lim_{\epsilon\rightarrow 0}\int_{S_{\epsilon}} G(r,r_0)\frac{\partial}{\partial n }u(r_0)\epsilon^2sin\phi d\phi d\theta \rightarrow 0 $$
Now what to do with;
$$ \lim_{\epsilon\rightarrow 0}\int_{S_{\epsilon}} u(r_0)\frac{\partial}{\partial n}G(r,r_0)\epsilon^2sin\phi d\phi d\theta $$
$$  = \lim_{\epsilon\rightarrow 0}\int_{S_{\epsilon}} u(r_0)(\frac{1}{\epsilon}-ik)\frac{e^{ik\epsilon}}{\epsilon}n(r_0)  \epsilon^2sin\phi d\phi d\theta $$
n(r_0) is normal to the surface. I can't see how this tends to zero to remain the equation unchanged. can someone jump in here
