One can work years in fiber optic communications and not really have to care about diffraction (it did come up in a lecture or two), but when working with free-space light propagation and/or illuminated surfaces it’s better to know the difference between your Fresnel and your Fraunhofer diffraction, and what the heck the Huygens-Fresnel principle is all about. Luckily, diffraction is really very straightforward once one grasps said Huygens-Fresnel principle; the rest is just mathematical sleights of hand as we’ll see.

A word of warning: While most of the following is extracted from books and other public sources, I make no claim on everything being totally correct. This is basically a summary of my own understanding of diffraction, which may possibly be faulty (though I sure hope it isn’t).

#### Huygens-Fresnel Principle

In Christiaan Huygens’ original principle, a wavefront was propagated by generating so-called *wavelets* at each point of the old wavefront to construct the new wavefront. The principle is shown in Fig. 1. The wavelets had a radius of the wavelength, and the process was iterative, wavelength-by-wavelength.^{1} Augustin-Jean Fresnel later improved on the principle by replacing the spherical wavelets with spherical waves and postulating that the waves interfere with each other. This way he was able to propagate the known field at some location to some other location in a single step (provided there were no obstacles between the two).

According to Wikipedia the principle dryly states “that each point of a medium (disturbed by a passing wave) becomes a source of disturbance which propagates from this point in all directions indiscriminately.” [1] For a field that is known on some surface, each point of the surface becomes a source of a spherical wave with the intensity and phase that the field has at this point. The waves from all the points are then superposed to obtain the field at other points in space. This can be nicely illustrated with Young’s double-slit experiment in which the field in the plane of the aperture is zero except inside the slits:

The beauty of this principle is that, despite being so simple, it is quite universal, at least for scalar fields (when polarization effects can be neglected) and can explain reflection, refraction and diffraction effects.

#### Diffraction Computations

Sometimes we don’t want to sit down for a day or two and draw little circles to figure out what the field emanating from some object is going to be. Also, if the aperture and/or the incident field is a little more complicated, this gets unwieldy fast. If we had some form of mathematical description we could sit back, relax, and let the computer do all the necessary calculations. The good news is that in order to do that we don’t need much more than the formulation of a spherical wave, which is [2]

$$\begin{align} E(\mathbf{r}) &= \frac{1}{2}\biggl( \frac{A_0}{\sqrt{4\pi}\bigl|\mathbf{r}-\mathbf{r}_0\bigr|} \exp \bigl(- i k \bigl|\mathbf{r}-\mathbf{r}_0\bigr| \bigr) + \mathrm{c.c} \biggr)\\

&= \Re\biggl\lbrace \frac{A_0}{\sqrt{4\pi}\bigl|\mathbf{r}-\mathbf{r}_0\bigr|} \exp \bigl(- i k \bigl|\mathbf{r}-\mathbf{r}_0\bigr| \bigr) \biggr\rbrace \end{align}\tag{1}$$

where a possible initial phase of the field is incorporated into the phase of complex $A_0$. Due to energy conservation the energy density falls off as $4\pi \bigl|\mathbf{r}-\mathbf{r}_0\bigr|^2$ and thus the field amplitude scales inversely with increasing distance from the point of origin of the spherical wave $\mathbf{r}_0$. Compare this post for use of a similar formalism (I am trying to keep a consistent notation), except that here we observe at a fixed point in time while in the other post we held the position constant. As in the other post, we can comfortably work with only the exponential argument of the real part function $\Re\bigl\lbrace\cdot\bigr\rbrace$ instead of the full field $E$ as long as our operations are linear (such as is integration). To distinguish the complex field from the real-valued one, we write

$$A(\mathbf{r}) = \frac{A_0}{\sqrt{4\pi}\bigl|\mathbf{r}-\mathbf{r}_0\bigr|} \exp \bigl(- i k \bigl|\mathbf{r}-\mathbf{r}_0\bigr| \bigr)\tag{2}$$

With (2) we are ready to write down the Huygens-Fresnel principle

$$A(\mathbf{q}) = \intop_S \kappa(\mathbf{p},\mathbf{q}) \frac{A(\mathbf{p})}{\bigl|\mathbf{q}-\mathbf{p}\bigr|} \exp \bigl(- i k \bigl|\mathbf{q}-\mathbf{p}\bigr| \bigr) \, dS\tag{3}$$

Looks a bit intimidating, but is quite straightforward, really, when we look at the illustration in Fig. 3. Here, $\mathbf{p}$ is the position of the element $dS$ on surface $S$ and $A(\mathbf{p})$ is proportional^{2} to the field at $\mathbf{p}$. The vector $\mathbf{q}$ describes the location on the screen, the field on which we are interested in. The integral then describes the superposition of spherical waves originating on each and every point on $S$ at the point $\mathbf{q}$. Finally, $\kappa$ is a “reality factor” to unite theory and observation, which we need to discuss a bit more in detail.

The inverse dependence of the spherical wave amplitude on $\bigl|\mathbf{q}-\mathbf{p}\bigr|$ means that as $\mathbf{q}$ approaches $\mathbf{p}$ we will obtain extremely high field amplitudes, which is unphysical – the area for which $\bigl|\mathbf{q}-\mathbf{p}\bigr| \lt \lambda$ is not shown in Fig. 2 for this very reason. Hence the Huygens-Fresnel construction should not be used to determine the field less that a few wavelengths from the aperture.

#### A Plethora of Coefficients

Turns out that to make (3) work for e.g. a plane unobstructed wave, $\kappa$ must be proportional to $1/i\lambda$, $\lambda$ being the wavelength. This means that an additional 90° phase shift exists between the oscillation of the field on $S$ and the screen. Furthermore, to prevent the (mathematical) propagation of backward waves – which would be possible with Huygens’ construction – Fresnel incorporated an *inclination* or *obliquity factor* into $\kappa$. Fresnel merely prescribed that it depends on the angle $\alpha$ between the surface normal of $S$ at $\mathbf{p}$ and the vector from $\mathbf{p}$ to $\mathbf{q}$ in the following way: for small angles it should be approximately unity and fall off (at first slowly then fast) to zero as the angle approaches $\pi/2$. This way the wave would propagate mainly in the direction of the surface normal at $\mathbf{p}$. For more detail on how to arrive at these expressions, see e.g. [3].

Gustav Kirchhoff showed that the same result could be obtained rigorously from the wave equation by choosing appropriate integration contours, without resorting to wavelets or arbitrarily choosing some coefficient $\kappa$. Kirchhoff’s solution sufficiently far away from the surface $S$ can be written in the form (3) by using

$$\kappa_\mathrm{K} = \frac{1-\cos \alpha}{2 i\lambda} \tag{4a}$$

Lord Rayleigh and Arnold Sommerfeld also derived two more solutions to the diffraction problem directly from the wave equation. The Rayleigh-Sommerfeld solutions are characterized by

$$\kappa_\mathrm{RS1} = \frac{\cos \alpha}{i\lambda} \tag{4b}$$

and

$$\kappa_\mathrm{RS2} = \frac{1}{i\lambda} \tag{4c}$$

Interestingly, the Kirchhoff solution is the arithmetic mean of the Rayleigh-Sommerfeld solutions. The first Rayleigh-Sommerfeld solution is often used used to formulate the Huygens-Fresnel principle. While the coefficient $1/i\lambda$ can be physically explained [3], the $\cos \alpha$ dependence cannot, it seems.^{3}

Since the coefficients affect the amplitude of the expanding spherical wave (differently) depending on $\alpha$, they should correctly also contain a normalization factor to satisfy energy conservation laws, see footnote 2. Most popular literature on the topic seems to avoid a discussion of this.

#### Fresnel’s Approximation

For the (simple) case that all $\mathbf{p}$ lie on the $z=0$ plane and all points $\mathbf{p}$ lie on a planar screen at some arbitrary $z$, we can write

$$\mathbf{p} = \begin{pmatrix} x_p \\ y_p \\ 0 \end{pmatrix} \quad \text{and} \quad \mathbf{q} = \begin{pmatrix} x_q \\ y_q \\ z \end{pmatrix}$$

If $z$ is large compared to the dimensions of the surface $S$ and the screen, all the coefficients from the previous section approach the same value and we can write

$$\kappa_\mathrm{K} \approx \kappa_\mathrm{RS1} \approx \kappa_\mathrm{RS2} = \frac{1}{i\lambda}\tag{6}$$

The problem now in solving the integral (3) in a general way are the $\bigl|\mathbf{q}-\mathbf{p}\bigr|$ terms. It is especially critical in the exponent, where it is multiplied by $k$ which makes the phase of this expression *very* sensitive to even small changes in $\mathbf{p}$ or $\mathbf{q}$. The denominator is rather relaxed – it hardly changes during integration when $z$ is large. Since we already assumed that this be true, we can set $\bigl|\mathbf{q}-\mathbf{p}\bigr| \approx z$ *in the denominator only*.

The exponential is a bit trickier, and this is where Fresnel’s approximation happens. We can write

$$\bigl|\mathbf{q}-\mathbf{p}\bigr| = \sqrt{\rho^2 + z^2} = z \sqrt{\frac{\rho^2}{z^2}+1}\tag{7a}$$

where

$$\rho^2 = \bigl(x_q – x_p\bigr)^2 + \bigl(y_q – y_p\bigr)^2\tag{7b}$$

was illustrated in Fig. 3. Fresnel then used the series expansion for the square root expression [5]^{4}:

$$\bigl|\mathbf{q}-\mathbf{p}\bigr| = z + \frac{\rho^2}{2z} – \frac{\rho^4}{8z^3} + …\tag{8}$$

If the contribution of the third term (and thus all higher order terms) to the exponential in (3) were negligible, i.e. for the resulting phase change in (3)

$$k \frac{\rho^4}{8z^3} = \frac{2\pi \rho^4}{8\lambda z^3} \ll 2\pi \tag{9}$$

we may use only the first two terms to approximate $\bigl|\mathbf{q}-\mathbf{p}\bigr|$. This holds for $\rho \ll z$ as long as $\lambda$ does not become too small. We can then write

$$\exp \bigl(- i k \bigl|\mathbf{q}-\mathbf{p}\bigr| \bigr) \approx \exp \bigl(- i k z \bigr) \exp \biggl(- i k \frac{\rho^2 }{2z}\biggr) \tag{10}$$

and by using (10) in (3) we obtain the Fresnel diffraction integral:

$$A(\mathbf{q}) = \frac{\exp \bigl(- i k z \bigr)}{i\lambda z} \intop_S A(\mathbf{p}) \exp \biggl(- i k \frac{\rho^2}{2z}\biggr) \, dS\tag{11}$$

Since $z$ is fixed, we only need to figure out $\rho$ when calculating the integral. When writing the surface integral over $S$ in polar coordinates centered on $(x_q, y_q)$, this becomes particularly easy.

So, while not being applicable to as many diffraction problems as the Huygens-Fresnel method (e.g. non-planar surfaces, large divergence angles, etc.), the Fresnel diffraction integral can be used to solve some diffraction problems analytically. The requirement of small angles $\alpha$, as described by $\rho \ll z$, is equivalent to the *paraxial approximation*, since both aperture and screen are limited to regions near the optical axis.

#### Fraunhofer’s Approximation

Joseph von Fraunhofer figured out that under certain restrictions the image on the screen would be the Fourier transform of the field on (the still planar) surface $S$. We start by expanding expression (7b) for the distance $\rho$

$$\begin{aligned}

\rho^2 &= \bigl(x_q – x_p\bigr)^2 + \bigl(y_q – y_p\bigr)^2\\

&= x_q^2 + x_p^2 – 2 x_q x_p + y_q^2 + y_p^2 – 2 y_q y_p\\

&= \bigl(x_q^2 + y_q^2 \bigr) + \bigl(x_p^2 + y_p^2\bigr) – 2 \mathbf{q} \cdot \mathbf{p} \tag{12}

\end{aligned}$$

with $\mathbf{p}$ and $\mathbf{q}$ as defined above.

Now the second term only depends on the dimensions of $S$, or rather the part of $S$ (e.g. an aperture or a surface illuminated by a laser spot) that has a non-negligible field $A(\mathbf{p})$. When such an aperture is much, much smaller than the distance to the screen, the contribution of this term to the phase in the Fresnel integral (11) is^{5}

$$k \frac{x_p^2 + y_p^2}{2z} = 2 \pi \frac{x_p^2 + y_p^2}{2\lambda z} \ll 2\pi \tag{13}$$

and can be neglected. The first term in (12) does not depend on the point $\mathrm{p}$ on $S$ and we can factor it outside the diffraction integral such that

$$A(\mathbf{q}) = a \intop_S A(\mathbf{p}) \exp \biggl(- i k \frac{\mathbf{q} \cdot \mathbf{p}}{z}\biggr) \, dS\tag{14a}$$

with

$$a = \frac{\exp \bigl(- i k z \bigr)}{i\lambda z} \exp \biggl(- i k \frac{x_q^2 + y_q^2}{2z}\biggr)\tag{14b}$$

Equation (14a) describes the two-dimensional Fourier transform of the field on $S$ for any point $\mathbf{q}$ on the screen. The restriction (13), however, is quite severe – e.g. at a wavelength of 633nm and for an aperture of 1mm, the distance $z$ to the screen must already be much larger than a meter

Interestingly, a correctly placed, regular (quadratic) lens will exactly remove the term in (12) that we approximated by zero in (13). Hence, the Fraunhofer approximation becomes the exact solution of the field on the screen after the lens, and is one way to show that a lens can act as a optical Fourier transformer.

#### Last Words

When computing the diffraction pattern of a field on some surface, we have basically three options: a) compute the Huygens-Fresnel diffraction integral (3) – this works with arbitrarily shaped surfaces and geometries, as long as the distance to the screen is larger than a few wavelengths (otherwise we may enjoy a rigorous treatment of Maxwell’s equations); b) compute the simpler Fresnel integral (11) if $z \gg \rho$, i.e. the involved angles are not too large, and the surface $S$ is on a plane; c) compute the Fraunhofer integral (14) if $z$ is really huge or if there is a lens involved (see [3]).

There is still another way to figure out the diffraction pattern from an arbitrary (planar) field distribution, which uses the two-dimensional Fourier transform to calculate an angular spectrum of plane waves from the source distribution. These can then be propagated to any surface in order to determine the resulting field on that surface. This will be the topic of another post, though.

#### Update

I posted a bunch of images calculated from the various integrals above and different apertures here. These are probably a bit more illustrative than just the math by itself.

^{1} Interestingly, the “Offset Path” function (with round joints) in graphics programs such as Adobe Illustrator seems to have the same result as a Huygens construction would have were the path a wavefront.

^{2} Why do we cite proportionality here instead of keeping the $1/\sqrt{4\pi}$ coefficient? In the section titled *A Plethora of Coefficients* the simple spherical wave is modified to account for the directionality between $\mathbf{p}$ and $\mathbf{q}$ with respect to the optical axis by introducing the obliquity factor $\kappa$. While it would be possible to incorporate $\kappa$ into the energy conservation calculation, it would be too off-topic for this post (and seemingly many books on this topic, either). After all, we are mainly interested in a qualitative description of the diffracted field and not the quantitative calculation of correct field amplitudes.

^{3} Initially, I thought that this would simply account for the spherical projection of the surface element $dS$. This can be visualized more easily if we look at the reverse propagation from $\mathbf{q}$ to $\mathbf{p}$. The spherical wave emanating from $\mathbf{q}$, having a uniform intensity distribution, illuminates $dS$. The solid angle (the relevant quantity when dealing with spherical waves) subtended by $dS$ is then $\cos \alpha \, dS$. But, and that’s a big but, this is valid for the optical intensity, not the field, of the spherical wave, so I think this would account only for a factor of $\sqrt{\cos \alpha}$ or something more complicated. I really wish I would have paid more attention in my math courses.

Then I thought, maybe each $dS$ can be modeled as a Lambertian scatterer [4] (without an underlying physical reason, though) which also has a $\cos \alpha$ characteristic, though again regarding intensity, not the field. Maybe it’s both and the square roots combine to give the right $\kappa$, but someone smarter has probably already figured all this out and not told me.

^{4} This is one of my most often used series expansions and very handy. You’ll need this e.g. when trying to figure out the spectrum (hence, the field) for an intensity-modulated signal.

^{5} This condition, written as

$$\frac{a^2}{2z} \ll 1 \quad \text{with} \quad a^2 = x_p^2 + y_p^2$$

being the aperture dimension, looks almost like the second term in (8). However, they are not the same since $\rho$ also contains the coordinates of $\mathbf{q}$. Neglecting the second term in the series expansion would cause the field on the screen to be the average of the field on the surface $S$ *everywhere* – certainly not the correct diffraction pattern.

[1] http://en.wikipedia.org/wiki/Huygens–Fresnel_principle

[1] http://scienceworld.wolfram.com/physics/SphericalWave.html

[2] R.W. Ditchburn, *Light*, Interscience Publishers, New York, 1961.

[3] J.W. Goodman, Introduction to Fourier Optics, 3rd ed., Roberts & Company Publishers, 2005.

[4] http://en.wikipedia.org/wiki/Lambert’s_cosine_law

[5] http://en.wikipedia.org/wiki/Square_root#Properties

last posts in diffraction:

[…] The best book I’ve found on Fourier Optics numerical simulation is « Computational Fourier Optics » by D. Voeltz. It is far better than another one my colleagues refer to, « Numerical Simulation of Optical Wave Propagation« by JD Schmidt. I found an interesting discussion on Nerd’s Rage blog (especially the article about on scalar diffraction) […]