Briefly speaking, from Helmholtz equation you can deduce the Rayleigh-Sommerfeld diffraction formula. Further, with paraxial approximation you deduce Fresnel diffraction formula.
In math, for wave $u_0(\mathbf{r})$ propagating distance $z$ as $u_z(\mathbf{r})$, by Rayleigh-Sommerfeld equation:
$$u_z(\mathbf{r}) = -\frac{1}{2\pi}u_0(\mathbf{r}) \star \frac{1}{2\pi}\frac{z e^{j kR}}{R^3}(1 - jkR),$$
where $\mathbf{r} = (x,y)$, $R = \sqrt{x^2 + y^2 + z^2}$, $k$ the wave number, $\star$ denotes convolution.
To get Fresnel formula, you need:
With these two aproximations, you yield the Frenel formula:
$$u_z(\mathbf{r}) = u_0(\mathbf{r}) \star \frac{e^{j kz}}{j\lambda z} e^{j \frac{k \|\mathbf{r}\|^2}{2z}}.$$
You start out with the Helmholtz equation that holds true in a homogeneous lossless medium for any rectangular component of both the $E$ and $H$ fields:
$$\nabla^2 U + \kappa^2 U = 0 \tag{1}\label{1}$$
That is $U$ may stand for any of the $E_x, E_y,...,H_z$. A general solution of $\eqref{1}$ is the plane wave $e^{\mathfrak {j} \kappa \hat {\mathbf{k}} \cdot \mathbf{r}}$ where $\hat {\mathbf{k}}$ is the unit vector in the direction of propagation and $\kappa = \omega/c =2\pi/\lambda$.
Now assume that a wave, say $U(x,y,z)$ propagates *mostly* along the $z$ axis (paraxial approximation) and is incident upon the $z=0$ plane. Decompose this field at an arbitrary $z$ into plane wave components using the 2D Fourier transform:
$$U(x,y,z)=\int \int dudv A_z(u,v) e^{\mathfrak j \kappa(xu+yv)}\tag{2}\label{2}$$
Here $$dU(x,y,z)= A_z(u,v) e^{\mathfrak j \kappa (xu+yv)} dudv$$ is a differential plane wave moving in the direction represented by the direction cosines $u,v,w$ where $u^2+v^2+w^2=1$, i.e., $w=\sqrt{1-u^2-v^2}$.
Next we wish to find the relationship between the field (component) $U(x,y,z)$ at $z=0$ and at arbitrary $z$. From linearity $ A_z(u,v) = A_0(u,v) e^{\mathfrak{j}\kappa w}$. The inverse Fourier transform gives:
$$A_0(x,y)=b_1\int \int dx'dy' U(x',y',0) e^{-\mathfrak j \kappa(x'u+y'v)}\tag{3}\label{3}$$ for some constant $b$ (in fact $b_1=1/(2\pi)^2$
Again because of the linear shift invariant nature of hte optical (EM) system the *output* field at $z$ is a convolution of the *input* field at $z=0$ and kernel, say, $G(x,y)$, that is
$$U(x,y,z) =\int \int dx'dy' U(x',y',0) G(x-x',y-y')\tag{4}\label{4}$$
and
$$G(x,y) = \int \int du dv e^{\mathfrak j \kappa w} e^{\mathfrak j \kappa (ux+vy)} \tag{5}\label{5}$$ where $w=\sqrt{1-u^2-v^2}$.
When polar coordinates are introduced $u=R cos\theta$ and $v=R sin \theta$ after much pain we get (see the Erdelyi Higher Transcendental Functions, 4.52)
$$G(R cos\theta,R sin \theta) = \frac{1}{\mathfrak j \lambda}\frac{e^{\mathfrak j \kappa \sqrt{R^2+z^2}}}{\sqrt{R^2 + z^2}} \frac{z}{\sqrt{R^2 + z^2}}\left( 1-\frac{1}{\mathfrak j \kappa \sqrt{R^2 + z^2}}\right) \tag{6}\label{6}$$
Now in the far field (Fraunhofer) we have
- when $z \gg R$ the term in parentheses is approximated as $1$
- Define $cos\phi = \frac{z}{\sqrt{R^2+z^2}}$, this angle $\phi$ is the direction the elementary wavelet has with the axis $z$ as it propagates, and $cos\phi$ is the obliquity factor introduced first by Fresnel. Taking $\phi \approx 0$ that is $cos\phi \approx 1$ is the paraxial approximation.
- The term $\sqrt{R^2+z^2}$ is the distance between the point $x',y'$ in the $z=0$ plane and the field point $x,y,z$, and in the paraxial regime this can be approximated as $z$.
Therefore
$$G(R cos\theta,R sin \theta)\approx \tilde G(R) = \frac{1}{\mathfrak j \lambda}\frac{e^{\mathfrak j \kappa \sqrt{R^2+z^2}}}{z} \tag{7}\label{7}$$
This gives us the field at $z$ approximately by calcuating the convolution integral as
$$ U(x,y,z)\approx \frac{b_2}{z}\int \int dx'dy'U(x',y',0)e^{\mathfrak j \kappa \sqrt{R^2+z^2}} \tag{8}\label{8}$$
Finally, in the Fraunhofer far-field we also have the approximation
$$\sqrt{R^2+z^2} = \sqrt{z^2+(x-x')^2+(y-y')^2} \\
\approx z+ \frac{1}{2z}(x-x')^2+\frac{1}{2z}(y-y')^2,$$ and thus
$$ U(x,y,z)\approx \frac{b_2 e^{\mathfrak j \kappa z}}{z}\int \int dx'dy'U(x',y',0)e^{\mathfrak j \kappa (\frac{1}{2z}(x-x')^2+\frac{1}{2z}(y-y')^2)} \tag{9}\label{9}$$
The result $\eqref{9}$ is the diffraction integral, the input (incident) field $U(x,y,0)$is at $z=0$ represents any component of $E$ or $H$ and is supposed to be known. It propagates in the free half-space $z>0$ and via this convolution integral one can calculate the same field component at any point $x,y,z$.
If you wish so you can also write the $\eqref{9}$ convolution integral as a Fourier transform by expanding the quadratic terms in the exponent:
$$ U(x,y,z)\approx \frac{b_2 e^{\mathfrak j \kappa (z+x^2/(2z) +y^2/(2z))}}{z}\int \int dx'dy' e^{\mathfrak j \tfrac{\kappa}{2z} (x'^2 +y'^2)} U(x',y',0)e^{-\mathfrak j \kappa (xx'+yy')/z} \tag{10}\label{10}$$
but notice that $\eqref{10}$ is not the Fourier transform of $U(x,y,0)$ but instead there is a quadratic phase factor also involved, so that aside from some unit modulus phase factor outside the integral this is the Fourier transform of $e^{\mathfrak j \kappa(x^2 +y^2)/(2z)} U(x,y,0)$. Now if the extent of the diffracting object at $z=0$ is such that $\tfrac{\kappa}{2z}(x^2 +y^2) \ll 1 $ or $x^2 +y^2 \ll z \lambda $ then the quadratic phase factor can be neglected and we do have the Fourier transform of the incident field.
[1] Stark: Application of Optical Fourier Transforms, Academic Press 1982
Best Answer
Sorry for my poor english. My native language is french.
I have no time to look at it in detail but it seems to me that you have to take into account the phase factors which multiply the Fourier transform. They are usually forgotten because they disappear in the calculation of the intensity. But they are not constant and here you have to keep them in the calculation of the second Fourier transform?
In Goodman, "introduction to Fourier optics" : $E(x_1,y_1 )=\exp(jkz)\exp(jk(x_1^2+y_1^2 )/2z) FT(E(x_0,y_0 ))$