differential geometry - Helmholtz - Kirchhoff Integral Theorem

I'm reading Introduction to Fourier Optics - J. Goodman and got to the Helmholtz-Kirchhoff Integral Theorem (HKIT). My problem with it is that in this material and all of the ones I looked into (Principles of Optics - M. Born, E. Wolf and some lecture notes found googling), the following reasoning is used.

1. 
   

   We consider Green's second identity where $U(P)$ is viewed as the disturbance made by the field at some point $P$ $$\iiint_V U \nabla^2G-G\nabla^2Udv=\iint_{\partial V} U\frac{\partial G}{\partial n} - G\frac{\partial U}{\partial n}ds$$ $U$ also satisfies the Helmholtz equation.

   
   


2. 
   
   

   We take a setup on which we will use Green's second identity that is given in the following image.

   
   

3. We use Kirchhoff's choice for the auxiliary function $G$ as follows, where $r_{01}$ is the distance from point $P_0$ to some point $P_1$ that is on the surface $S+S_{\epsilon}$ $$G=\frac{exp(ikr_{01})}{r_{01}}$$ As far as I understoond, the reason for taking $G$ of this form is such that the LHS of Green's identity will become $0$.


4. We plug $G$ from point $3$ in Green's second identity, do some computations and get to HKIT that is

$$U(P_0)=\frac{1}{4\pi}\iint_{\partial V} \frac{\partial U}{\partial n}\Bigr[\frac{exp(ikr_{01})}{r_{01}}\Bigr] -U\frac{\partial}{\partial n}\Bigr[\frac{exp(ikr_{01})}{r_{01}}\Bigr]ds$$

Ok, I kind of understand how this goes, although to me it looks like sticking some formulas with duct tape. What I consider not very clear is the choice on $G$. First time I read this chapter, I considered that $G$ describes a point source at $P_0$. This later made no sense, since we were interested in computing $U(P_0)$. So, if $G$ described what happens at $P_0$ there would have been no need to compute $U(P_0)$.

Later I thought that given the formula for $G$ I can look at it in a reversed manner, that is it gives the amplitude at $P_0$ from a point source on $S$. This made a lot more sense to me, but then it was not clear why can't we just integrate the contribution from all the sources on $S$ and write

$$U(P_0)=\iint_S \frac{exp(ikr_{01})}{r_{01}} ds$$

Then I considered the following scenario. Suppose I have a point source at the center of the reference system that is described by $\psi(r) = A_0\frac{exp(ikr)}{r}$. Now consider a sphere of radius $R$ that has in its center the point source. The existence of the source allows based on Huygens principle to onsider that every point on the sphere acts like a source of amplitude $\psi(R)=A_0\frac{exp(ikR)}{R}$, and also, by considering the contributions on those secondary sources we must be able to retrieve the equation for the first point source by computing

$$\frac{1}{4\pi}\iint_S \psi(R) \frac{exp(-ikR)}{R} ds = A_0$$

However HKIT does not look like the result I got above.

Tl;dr Does $G$ have a physical meaning in HKIT, and if so where can I find it explained. For my attempt at doing the same thing as HKIT for a particular case, is the result correct? If so, why can't HKIT be expressed without the use of Green's second identity? If not, where did I mess up?

Answer

Not to distract from the greatness of the work that these people did way back or in any way diminish the value of Goodman's book on this topic, one can use a different approach to treat diffraction, which is equally rigorous and in my view much more understandable. It does not need any duct tape. By explaining this approach, I will hopefully address your questions.

The idea is to treat free space as a system. This system is linear and shift invariant. The input to the system is the complex function that represents the scalar optical field in the input plane $g(x,y,0)$ and the output is a complex optical field in the output plane $g(x,y,z)$, located at some distance $z$ from the input plane. Based on the linearity and shift-invariance, one can immediate write down a general expression for the output in terms of the input $$ g(x,y,z)=\int g(x',y',0) G(x-x',y-y';z) dx dy . $$ You may notice that this is a convolution integral and $G(x,y;z)$ is a Green function (also called the impulse response in linear systems theory).

How does one get the expression for this Green function? One can proceed as follows, using a standard approach. The idea is to expand the input function in terms of the eigenfunctions of the system and then reconstruct the output function from these eigenfunctions after one allowed them to pass through the system.

What are the eigenfunctions for free space propagation? They would be the solutions of the Helmholtz equation; for our purposes here we can use the plane waves. What happens to plane waves when they propagate over a distance $z$? They pick up a phase factor that depends on the propagation vector and the propagation distance. In particular, $$ \exp[-i (x k_x + y k_y)] \rightarrow \exp[-i (x k_x + y k_y + z k_z)] , $$ where $$ k_z = \sqrt{k^2-k_x^2-k_y^2} . $$

Now to expand the input function in terms of plane waves, one simply needs to perform a 2D Fourier transform of the input function. The result is the angular spectrum $F(k_x,k_y)$. Then one would multiply this angular spectrum with the propagation phase factor $\exp(-i z k_z)$, and perform the inverse Fourier transform to obtain the output function.

The whole process can be expressed as $$ g(x,y,z)=\int \int g(x',y',0) e^{i (x' k_x + y' k_y)} dx dy\ e^{-i (x k_x + y k_y + z k_z)} d k_x d k_y . $$ The integral over $x,y$ represents the first Fourier transform and the integral over $k_x,k_y$ is the final inverse Fourier transform.

One can now interchange the order of integration and write the integral as $$ g(x,y,z)=\int g(x',y',0) \int e^{-i [(x-x') k_x + (y-y') k_y + z k_z]} d k_x d k_y dx dy . $$ Then one can first evaluate the integrals over $k_x,k_y$ to get $$ g(x,y,z)=\int g(x',y',0) G(x-x',y-y';z) dx dy , $$ where $$ G(x,y;z)=\int e^{-i (x k_x + y k_y + z k_z)} d k_x d k_y . $$

So we see that we do indeed get a convolution integral, and the Green function is simply the integral over all plane waves with $k_z$ expressed in terms of $k_x$ and $k_y$. However, the latter integral is not so easy to evaluate. In the paraxial limit, it gives the Fresnel kernel. If one evaluates the integral for the Green function numerically one can see that it resembles the well-known Green function for free space propagation $$ G(r) = \frac{\exp(i k r)}{r} . $$ However, it is not exactly the same.

Hopefully this gives you a more intuitively appealing way to understand free-space propagation.