Helmholtz - Kirchhoff Integral Theorem

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.

The standard method of proof of the HKIT exploits a general relation between two functions (Green's second identity) to deduce a property of a general function ($U$) from the properties of a particular one ($G$). In this case the "general relation" is simplified (its LHS becomes zero) because both functions are solutions of the Helmholtz equation ($U$ by hypothesis, and $G$ by construction). The "duct tape" impression arises because it is counterintuitive to deduce a property of one function from the properties of another.

Using the HKIT, or otherwise,* one can show that the wave function in a region containing no sources, due to "primary" sources outside that region, is as if the primary sources were replaced by a certain distribution of "secondary" sources on the boundary of the region. But the secondary sources turn out to have a certain directionality — which is not captured by your last proposed integral, which is therefore not equivalent to the HKIT.

* Here's my attempt at doing it "otherwise", whereby I get a form of the HKIT from the secondary sources, instead of the other way around: "Exact derivation of Kirchhoff's integral theorem and diffraction formula using high-school math". So, unless I screwed up big-time, the HKIT can be obtained without Green's identity (and without assuming that the surface of integration is planar).