Scalar Diffraction - Huygens, Fresnel, Fraunhofer

One can work years in fiber optic com­mu­ni­ca­tions and not real­ly have to care about dif­frac­tion (it did come up in a lec­ture or two), but when work­ing with free-space light prop­a­ga­tion and/or illu­mi­nat­ed sur­faces it’s bet­ter to know the dif­fer­ence between your Fres­nel and your Fraun­hofer dif­frac­tion, and what the heck the Huy­gens-Fres­nel prin­ci­ple is all about. Luck­i­ly, dif­frac­tion is real­ly very straight­for­ward once one grasps said Huy­gens-Fres­nel prin­ci­ple; the rest is just math­e­mat­i­cal sleights of hand as we’ll see.

A word of warn­ing: While most of the fol­low­ing is extract­ed from books and oth­er pub­lic sources, I make no claim on every­thing being total­ly cor­rect. This is basi­cal­ly a sum­ma­ry of my own under­stand­ing of dif­frac­tion, which may pos­si­bly be faulty (though I sure hope it isn’t).

Huygens-Fresnel Principle

In Chris­ti­aan Huy­gens’ orig­i­nal prin­ci­ple, a wave­front was prop­a­gat­ed by gen­er­at­ing so-called wavelets at each point of the old wave­front to con­struct the new wave­front. The prin­ci­ple is shown in Fig. 1. The wavelets had a radius of the wave­length, and the process was iter­a­tive, wave­length-by-wave­length.1 Augustin-Jean Fres­nel lat­er improved on the prin­ci­ple by replac­ing the spher­i­cal wavelets with spher­i­cal waves and pos­tu­lat­ing that the waves inter­fere with each oth­er. This way he was able to prop­a­gate the known field at some loca­tion to some oth­er loca­tion in a sin­gle step (pro­vid­ed there were no obsta­cles between the two).

Fig­ure 1: Illus­tra­tion of Huy­gens’ con­struc­tion of a prop­a­gat­ing wave. A spher­i­cal wave impinges on an aper­ture (blue). The evo­lu­tion of the result­ing wave­front (red) can then be con­struct­ed iter­a­tive­ly using wavelets (gray) with the radius of a wave­length.

Accord­ing to Wikipedia the prin­ci­ple dry­ly states “that each point of a medi­um (dis­turbed by a pass­ing wave) becomes a source of dis­tur­bance which prop­a­gates from this point in all direc­tions indis­crim­i­nate­ly.” [1] For a field that is known on some sur­face, each point of the sur­face becomes a source of a spher­i­cal wave with the inten­si­ty and phase that the field has at this point. The waves from all the points are then super­posed to obtain the field at oth­er points in space. This can be nice­ly illus­trat­ed with Young’s dou­ble-slit exper­i­ment in which the field in the plane of the aper­ture is zero except inside the slits:

Fig­ure 2: Illus­tra­tion of the Huy­gens-Fres­nel con­struc­tion with the dou­ble slit. Top: thin black lines rep­re­sent wave­fronts sep­a­rat­ed by one wave­length. The incom­ing mono­chro­mat­ic plane wave is par­tial­ly blocked at an aper­ture with two nar­row slits (width is much small­er than the wave­length). Cylin­dri­cal waves emanate from the holes in the aper­ture and inter­fere. Blue dashed lines show areas of con­struc­tive inter­fer­ence at which the total inten­si­ty of the com­pound wave will be high. Bot­tom: same con­struc­tion shown with col­ors rep­re­sent­ing the field ampli­tude (red for pos­i­tive val­ues, green for neg­a­tive). Blue dashed lines again show the direc­tions of high­est inten­si­ty.

The beau­ty of this prin­ci­ple is that, despite being so sim­ple, it is quite uni­ver­sal, at least for scalar fields (when polar­iza­tion effects can be neglect­ed) and can explain reflec­tion, refrac­tion and dif­frac­tion effects.

Diffraction Computations

Some­times we don’t want to sit down for a day or two and draw lit­tle cir­cles to fig­ure out what the field ema­nat­ing from some object is going to be. Also, if the aper­ture and/or the inci­dent field is a lit­tle more com­pli­cat­ed, this gets unwieldy fast. If we had some form of math­e­mat­i­cal descrip­tion we could sit back, relax, and let the com­put­er do all the nec­es­sary cal­cu­la­tions. The good news is that in order to do that we don’t need much more than the for­mu­la­tion of a spher­i­cal 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} \big­gr)\\
&= \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 pos­si­ble ini­tial phase of the field is incor­po­rat­ed into the phase of com­plex $A_0$. Due to ener­gy con­ser­va­tion the ener­gy den­si­ty falls off as $4\pi \bigl|\mathbf{r}-\mathbf{r}_0\bigr|^2$ and thus the field ampli­tude scales inverse­ly with increas­ing dis­tance from the point of ori­gin of the spher­i­cal wave $\mathbf{r}_0$. Com­pare this post for use of a sim­i­lar for­mal­ism (I am try­ing to keep a con­sis­tent nota­tion), except that here we observe at a fixed point in time while in the oth­er post we held the posi­tion con­stant. As in the oth­er post, we can com­fort­ably work with only the expo­nen­tial argu­ment of the real part func­tion $\Re\bigl\lbrace\cdot\bigr\rbrace$ instead of the full field $E$ as long as our oper­a­tions are lin­ear (such as is inte­gra­tion). To dis­tin­guish the com­plex field from the real-val­ued 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 Huy­gens-Fres­nel prin­ci­ple

$$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 intim­i­dat­ing, but is quite straight­for­ward, real­ly, when we look at the illus­tra­tion in Fig. 3. Here, $\mathbf{p}$ is the posi­tion of the ele­ment $dS$ on sur­face $S$ and $A(\mathbf{p})$ is pro­por­tion­al2 to the field at $\mathbf{p}$. The vec­tor $\mathbf{q}$ describes the loca­tion on the screen, the field on which we are inter­est­ed in. The inte­gral then describes the super­po­si­tion of spher­i­cal waves orig­i­nat­ing on each and every point on $S$ at the point $\mathbf{q}$. Final­ly, $\kap­pa$ is a “real­i­ty fac­tor” to unite the­o­ry and obser­va­tion, which we need to dis­cuss a bit more in detail.

Fig­ure 3: Geom­e­try of the dif­frac­tion prob­lem, as described by (3). Top fig­ure shows side view, bot­tom fig­ure shows the view through the screen onto the sur­face S (irreg­u­lar shape). The vari­able ρ describes the length of the pro­jec­tion of the red arrow onto the x-y plane and is used in the Fres­nel approx­i­ma­tion below.

The inverse depen­dence of the spher­i­cal wave ampli­tude on $\bigl|\mathbf{q}-\mathbf{p}\bigr|$ means that as $\mathbf{q}$ approach­es $\mathbf{p}$ we will obtain extreme­ly high field ampli­tudes, which is unphys­i­cal – the area for which $\bigl|\mathbf{q}-\mathbf{p}\bigr| \lt \lamb­da$ is not shown in Fig. 2 for this very rea­son. Hence the Huy­gens-Fres­nel con­struc­tion should not be used to deter­mine the field less that a few wave­lengths from the aper­ture.

A Plethora of Coefficients

Turns out that to make (3) work for e.g. a plane unob­struct­ed wave, $\kap­pa$ must be pro­por­tion­al to $1/i\lambda$, $\lamb­da$ being the wave­length. This means that an addi­tion­al 90° phase shift exists between the oscil­la­tion of the field on $S$ and the screen. Fur­ther­more, to pre­vent the (math­e­mat­i­cal) prop­a­ga­tion of back­ward waves – which would be pos­si­ble with Huy­gens’ con­struc­tion – Fres­nel incor­po­rat­ed an incli­na­tion or obliq­ui­ty fac­tor into $\kap­pa$. Fres­nel mere­ly pre­scribed that it depends on the angle $\alpha$ between the sur­face nor­mal of $S$ at $\mathbf{p}$ and the vec­tor from $\mathbf{p}$ to $\mathbf{q}$ in the fol­low­ing way: for small angles it should be approx­i­mate­ly uni­ty and fall off (at first slow­ly then fast) to zero as the angle approach­es $\pi/2$. This way the wave would prop­a­gate main­ly in the direc­tion of the sur­face nor­mal at $\mathbf{p}$. For more detail on how to arrive at these expres­sions, see e.g. [3].

Gus­tav Kirch­hoff showed that the same result could be obtained rig­or­ous­ly from the wave equa­tion by choos­ing appro­pri­ate inte­gra­tion con­tours, with­out resort­ing to wavelets or arbi­trar­i­ly choos­ing some coef­fi­cient $\kap­pa$. Kirchhoff’s solu­tion suf­fi­cient­ly far away from the sur­face $S$ can be writ­ten in the form (3) by using

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

Lord Rayleigh and Arnold Som­mer­feld also derived two more solu­tions to the dif­frac­tion prob­lem direct­ly from the wave equa­tion. The Rayleigh-Som­mer­feld solu­tions are char­ac­ter­ized by

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

Inter­est­ing­ly, the Kirch­hoff solu­tion is the arith­metic mean of the Rayleigh-Som­mer­feld solu­tions. The first Rayleigh-Som­mer­feld solu­tion is often used used to for­mu­late the Huy­gens-Fres­nel prin­ci­ple. While the coef­fi­cient $1/i\lambda$ can be phys­i­cal­ly explained [3], the $\cos \alpha$ depen­dence can­not, it seems.3

Since the coef­fi­cients affect the ampli­tude of the expand­ing spher­i­cal wave (dif­fer­ent­ly) depend­ing on $\alpha$, they should cor­rect­ly also con­tain a nor­mal­iza­tion fac­tor to sat­is­fy ener­gy con­ser­va­tion laws, see foot­note 2. Most pop­u­lar lit­er­a­ture on the top­ic seems to avoid a dis­cus­sion of this.

Fresnel’s Approximation

For the (sim­ple) case that all $\mathbf{p}$ lie on the $z=0$ plane and all points $\mathbf{p}$ lie on a pla­nar screen at some arbi­trary $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 com­pared to the dimen­sions of the sur­face $S$ and the screen, all the coef­fi­cients from the pre­vi­ous sec­tion approach the same val­ue and we can write

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

The prob­lem now in solv­ing the inte­gral (3) in a gen­er­al way are the $\bigl|\mathbf{q}-\mathbf{p}\bigr|$ terms. It is espe­cial­ly crit­i­cal in the expo­nent, where it is mul­ti­plied by $k$ which makes the phase of this expres­sion very sen­si­tive to even small changes in $\mathbf{p}$ or $\mathbf{q}$. The denom­i­na­tor is rather relaxed – it hard­ly changes dur­ing inte­gra­tion 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 denom­i­na­tor only.

The expo­nen­tial is a bit trick­i­er, and this is where Fresnel’s approx­i­ma­tion hap­pens. We can write

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


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

was illus­trat­ed in Fig. 3. Fres­nel then used the series expan­sion for the square root expres­sion [5]4:

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

If the con­tri­bu­tion of the third term (and thus all high­er order terms) to the expo­nen­tial in (3) were neg­li­gi­ble, i.e. for the result­ing 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 approx­i­mate $\bigl|\mathbf{q}-\mathbf{p}\bigr|$. This holds for $\rho \ll z$ as long as $\lamb­da$ 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 \big­gl(- i k \frac{\rho^2 }{2z}\biggr) \tag{10}$$

and by using (10) in (3) we obtain the Fres­nel dif­frac­tion inte­gral:

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

Since $z$ is fixed, we only need to fig­ure out $\rho$ when cal­cu­lat­ing the inte­gral. When writ­ing the sur­face inte­gral over $S$ in polar coor­di­nates cen­tered on $(x_q, y_q)$, this becomes par­tic­u­lar­ly easy.

So, while not being applic­a­ble to as many dif­frac­tion prob­lems as the Huy­gens-Fres­nel method (e.g. non-pla­nar sur­faces, large diver­gence angles, etc.), the Fres­nel dif­frac­tion inte­gral can be used to solve some dif­frac­tion prob­lems ana­lyt­i­cal­ly. The require­ment of small angles $\alpha$, as described by $\rho \ll z$, is equiv­a­lent to the parax­i­al approx­i­ma­tion, since both aper­ture and screen are lim­it­ed to regions near the opti­cal axis.

Fraunhofer’s Approximation

Joseph von Fraun­hofer fig­ured out that under cer­tain restric­tions the image on the screen would be the Fouri­er trans­form of the field on (the still pla­nar) sur­face $S$. We start by expand­ing expres­sion (7b) for the dis­tance $\rho$

\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}

with $\mathbf{p}$ and $\mathbf{q}$ as defined above.
Now the sec­ond term only depends on the dimen­sions of $S$, or rather the part of $S$ (e.g. an aper­ture or a sur­face illu­mi­nat­ed by a laser spot) that has a non-neg­li­gi­ble field $A(\mathbf{p})$. When such an aper­ture is much, much small­er than the dis­tance to the screen, the con­tri­bu­tion of this term to the phase in the Fres­nel inte­gral (11) is5

$$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 neglect­ed. The first term in (12) does not depend on the point $\mathrm{p}$ on $S$ and we can fac­tor it out­side the dif­frac­tion inte­gral such that

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


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

Equa­tion (14a) describes the two-dimen­sion­al Fouri­er trans­form of the field on $S$ for any point $\mathbf{q}$ on the screen. The restric­tion (13), how­ev­er, is quite severe – e.g. at a wave­length of 633nm and for an aper­ture of 1mm, the dis­tance $z$ to the screen must already be much larg­er than a meter

Inter­est­ing­ly, a cor­rect­ly placed, reg­u­lar (qua­drat­ic) lens will exact­ly remove the term in (12) that we approx­i­mat­ed by zero in (13). Hence, the Fraun­hofer approx­i­ma­tion becomes the exact solu­tion of the field on the screen after the lens, and is one way to show that a lens can act as a opti­cal Fouri­er trans­former.

Last Words

When com­put­ing the dif­frac­tion pat­tern of a field on some sur­face, we have basi­cal­ly three options: a) com­pute the Huy­gens-Fres­nel dif­frac­tion inte­gral (3) – this works with arbi­trar­i­ly shaped sur­faces and geome­tries, as long as the dis­tance to the screen is larg­er than a few wave­lengths (oth­er­wise we may enjoy a rig­or­ous treat­ment of Maxwell’s equa­tions); b) com­pute the sim­pler Fres­nel inte­gral (11) if $z \gg \rho$, i.e. the involved angles are not too large, and the sur­face $S$ is on a plane; c) com­pute the Fraun­hofer inte­gral (14) if $z$ is real­ly huge or if there is a lens involved (see [3]).


There is still anoth­er way to fig­ure out the dif­frac­tion pat­tern from an arbi­trary (pla­nar) field dis­tri­b­u­tion, which uses the two-dimen­sion­al Fouri­er trans­form to cal­cu­late an angu­lar spec­trum of plane waves from the source dis­tri­b­u­tion. These can then be prop­a­gat­ed to any sur­face in order to deter­mine the result­ing field on that sur­face. This will be the top­ic of anoth­er post, though.


I post­ed a bunch of images cal­cu­lat­ed from the var­i­ous inte­grals above and dif­fer­ent aper­tures here. These are prob­a­bly a bit more illus­tra­tive than just the math by itself.

1 Inter­est­ing­ly, the “Off­set Path” func­tion (with round joints) in graph­ics pro­grams such as Adobe Illus­tra­tor seems to have the same result as a Huy­gens con­struc­tion would have were the path a wave­front.

2 Why do we cite pro­por­tion­al­i­ty here instead of keep­ing the $1/\sqrt{4\pi}$ coef­fi­cient? In the sec­tion titled A Pletho­ra of Coef­fi­cients the sim­ple spher­i­cal wave is mod­i­fied to account for the direc­tion­al­i­ty between $\mathbf{p}$ and $\mathbf{q}$ with respect to the opti­cal axis by intro­duc­ing the obliq­ui­ty fac­tor $\kap­pa$. While it would be pos­si­ble to incor­po­rate $\kap­pa$ into the ener­gy con­ser­va­tion cal­cu­la­tion, it would be too off-top­ic for this post (and seem­ing­ly many books on this top­ic, either). After all, we are main­ly inter­est­ed in a qual­i­ta­tive descrip­tion of the dif­fract­ed field and not the quan­ti­ta­tive cal­cu­la­tion of cor­rect field ampli­tudes.

3 Ini­tial­ly, I thought that this would sim­ply account for the spher­i­cal pro­jec­tion of the sur­face ele­ment $dS$. This can be visu­al­ized more eas­i­ly if we look at the reverse prop­a­ga­tion from $\mathbf{q}$ to $\mathbf{p}$. The spher­i­cal wave ema­nat­ing from $\mathbf{q}$, hav­ing a uni­form inten­si­ty dis­tri­b­u­tion, illu­mi­nates $dS$. The sol­id angle (the rel­e­vant quan­ti­ty when deal­ing with spher­i­cal waves) sub­tend­ed by $dS$ is then $\cos \alpha \, dS$. But, and that’s a big but, this is valid for the opti­cal inten­si­ty, not the field, of the spher­i­cal wave, so I think this would account only for a fac­tor of $\sqrt{\cos \alpha}$ or some­thing more com­pli­cat­ed. I real­ly wish I would have paid more atten­tion in my math cours­es.
Then I thought, maybe each $dS$ can be mod­eled as a Lam­bert­ian scat­ter­er [4] (with­out an under­ly­ing phys­i­cal rea­son, though) which also has a $\cos \alpha$ char­ac­ter­is­tic, though again regard­ing inten­si­ty, not the field. Maybe it’s both and the square roots com­bine to give the right $\kap­pa$, but some­one smarter has prob­a­bly already fig­ured all this out and not told me.

4 This is one of my most often used series expan­sions and very handy. You’ll need this e.g. when try­ing to fig­ure out the spec­trum (hence, the field) for an inten­si­ty-mod­u­lat­ed sig­nal.

5 This con­di­tion, writ­ten as

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

being the aper­ture dimen­sion, looks almost like the sec­ond term in (8). How­ev­er, they are not the same since $\rho$ also con­tains the coor­di­nates of $\mathbf{q}$. Neglect­ing the sec­ond term in the series expan­sion would cause the field on the screen to be the aver­age of the field on the sur­face $S$ every­where – cer­tain­ly not the cor­rect dif­frac­tion pat­tern.

[1] http://​sci​ence​world​.wol​fram​.com/​p​h​y​s​i​c​s​/​S​p​h​e​r​i​c​a​l​W​a​v​e​.​h​tml
[2] R.W. Ditch­burn, Light, Inter­science Pub­lish­ers, New York, 1961.
[3] J.W. Good­man, Intro­duc­tion to Fouri­er Optics, 3rd ed., Roberts & Com­pa­ny Pub­lish­ers, 2005.
[5] http://​en​.wikipedia​.org/​w​i​k​i​/​S​q​u​a​r​e​_​r​o​o​t​#​P​r​o​p​e​r​t​ies

last posts in diffraction:

One Comment

  • […] The best book I’ve found on Fouri­er Optics numer­i­cal sim­u­la­tion is « Com­pu­ta­tion­al Fouri­er Optics » by D. Voeltz. It is far bet­ter than anoth­er one my col­leagues refer to, « Numer­i­cal Sim­u­la­tion of Opti­cal Wave Prop­a­ga­tion«   by JD Schmidt. I found an inter­est­ing dis­cus­sion on Nerd’s Rage blog (espe­cial­ly the arti­cle about on scalar dif­frac­tion) […]

Post a Comment

Your email is never shared.