Validity of Fresnel diffraction integral for arbitrary field The Fresnel diffraction integral is used to calculate the electric field after it has been propagated over a distance $L$. Usually, the validity of the Fresnel diffraction integral is given by an upper limit of the Fresnel number. This number is easy to calculate for specific applications like a homogenously illuminated square or circular aperture. However, I like to numerically propagate arbitrary field distributions. How do I make sure that the Fresnel approximation is still valid?
 A: The Fresnel approximation is basically the same as the paraxial approximation. The latter makes some assumptions about the angular spectrum of the beam that you want to propagate. Let's denote this by $F(\mathbf{a})$, where $\mathbf{a}$ is the two-dimensional spatial frequency vector (the Fourier domain coordinates). To get the angular spectrum, you can simply compute the two-dimensional Fourier transform of your optical beam 
$$ F(\mathbf{a})={\cal F}\{f(\mathbf{x})\} , $$ 
where $f(\mathbf{x})$ is the complex function in a plane perpendicular to the propagation direction that represents the cross-section of your beam. Then one computes the modulus square of this angular spectrum $|F(\mathbf{a})|^2$. This is basically the spatial power spectral density of your beam. If your beam satisfies the requirement for paraxial propagation, then the bulk of the optical power of this power spectral density 
$$ {\rm optical~power} = \iint_A |F(\mathbf{a})|^2\ {\rm d}^2 a , $$
(integrated power spectral density over a specific area $A$) would sit inside a radius that is much smaller than (at most one tenth of) the maximum radius for propagating wave, which is given by $1/\lambda$. Stated differently, if one were to compute the variance
$$\sigma^2 = \iint_{-\infty}^{\infty} |\mathbf{a}|^2 |F(\mathbf{a})|^2\ {\rm d}^2 a . $$
(Here we assumed the first moment is zero.) Then the result would give one a scale for the size of the spectrum on the Fourier plane. If this scale is much smaller than the radius of the propagating waves on the Fourier plane
$$ \sigma \ll \frac{1}{\lambda} , $$
then the beam is paraxial.
For a more comprehensive understanding of this, one needs to study Fourier optics.
A: To answer this question, let's derive the Fresnel integral and see what we need to assume to make it valid.
You begin with the Helmholtz equation in a homogeneous medium $(\nabla^2 + k^2)\psi = 0$. All transverse Cartesian components of the electromagnetic field vectors (electric and magnetic fields) fulfill this equation in a homogeneous medium, as do all the Cartesian/time components of the potentials in Lorenz gauge.
If the field comprises only plane waves in the positive $z$ direction then we can represent the diffraction of any scalar field on any transverse (of the form $z=c$) plane by:
$$\begin{array}{lcl}\psi(x,y,z) &=& \frac{1}{2\pi}\int_{\mathbb{R}^2} \left[\exp\left(i \left(k_x x + k_y y\right)\right) \exp\left(i \left(k-\sqrt{k^2 - k_x^2-k_y^2}\right) z\right)\,\Psi(k_x,k_y)\right]{\rm d} k_x {\rm d} k_y\\
\Psi(k_x,k_y)&=&\frac{1}{2\pi}\int_{\mathbb{R}^2} \exp\left(-i \left(k_x u + k_y v\right)\right)\,\psi(x,y,0)\,{\rm d} u\, {\rm d} v\end{array}\tag{1}$$
To understand this, let's put carefully into words the algorithmic steps encoded in these two equations:


*

*Take the Fourier transform of the scalar field over a transverse plane to express it as a superposition of scalar plane waves $\psi_{k_x,k_y}(x,y,0) = \exp\left(i \left(k_x x + k_y y\right)\right)$ with superposition weights $\Psi(k_x,k_y)$. We can do this, and be assured our decomposition is unique, with reasonable assumptions about $\psi$, e.g. that $\psi$ should be a Tempered Distribution is a sufficient condition for existence and uniqueness of such a superposition;

*Note that plane waves propagating in the $+z$ direction fulfilling the Helmholtz equation vary as $\psi_{k_x,k_y}(x,y,z) = \exp\left(i \left(k_x x + k_y y\right)\right) \exp\left(i \left(k-\sqrt{k^2 - k_x^2-k_y^2}\right) z\right)$ (this is where the Helmholtz equation enters our argument);

*Propagate each such plane wave from the $z=0$ plane to the general $z$ plane using the plane wave solution noted in step 2;

*Inverse Fourier transform the propagated waves to reassemble the field at the general $z$ plane.


Now we look at (1) and see how we can approximate it. Indeed, you should be able to see that you will get Fresnel's integral if we can assume:
$$k-\sqrt{k^2 - k_x^2-k_y^2}\approx \frac{1}{2\,k}(k_x^2+k_y^2)\tag{2}$$
i.e. that the first order Taylor approximation of the exact quantity on the left of (2) is valid. This is what Flippiefanus's answer means by "paraxial approximation". 
