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
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:
$z \gg \|\mathbf{r}\|$. And thus: $R = \sqrt{\|\mathbf{r}\|^2 + z^2} \approx z + \frac{\|\mathbf{r}\|^2}{2z}$.
$z \gg \lambda$. This leads to: $kR \gg 1$.
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}}.$$