## Scalar Diffraction - Fourier

In an ear­li­er post we had a look at var­i­ous dif­frac­tion for­malisms that either orig­i­nat­ed in the Huy­gens-Fres­nel prin­ci­ple or led to essen­tial­ly the same results. The prin­ci­ple mod­eled a field dis­tri­b­u­tion in a aper­ture (or on a sur­face) as a source of infi­nite­ly many spher­i­cal waves whose ampli­tudes and phas­es were pre­scribed by the aper­ture field. How­ev­er, this spher­i­cal wave decom­po­si­tion is not the only way to look at the prob­lem. Pro­vid­ed that the aper­ture / sur­face is two-dimen­sion­al (flat), we can use Fouri­er analy­sis to decom­pose any field on this sur­face into a num­ber of plane waves, which can then be prop­a­gat­ed to any oth­er sur­face where they inter­fere (like in the Huy­gens-Fres­nel prin­ci­ple).

#### Going Fourier

Sim­i­lar to the fre­quen­cy decom­po­si­tion of a time sig­nal using the Fouri­er trans­form, one can use the two-dimen­sion­al Fouri­er trans­form to decom­pose a two-dimen­sion­al field in some aper­ture at $z=0$ into its spa­tial fre­quen­cies:

$$\tilde A(\omega_x, \omega_y) = \intop^\infty \!\!\!\!\!\!\!\! \intop_{-\infty} A(\mathbf{p}) \exp\bigl( i \omega_x x + i \omega_y y\bigl) dx\,dy \tag{1a}$$

with

$$\mathbf{p} = \begin{pmatrix} x \\ y \\ 0 \end{pmatrix} \tag{1b}$$

At this point the spa­tial fre­quen­cies $\omega_x$ and $\omega_y$ are some­what abstract quan­ti­ties. This is just Fouri­er analy­sis for­mal­ism. Since $x$ and $y$ (and $z$) are posi­tion quan­ti­ties, the cor­re­spond­ing Fouri­er trans­form quan­ti­ties must actu­al­ly be (angu­lar) wavenum­bers [1]. Hence we can use the wavevec­tor,

$$\mathbf{k} = \begin{pmatrix} k_x \\ k_y \\ k_z \end{pmatrix} = \begin{pmatrix} \omega_x \\ \omega_y \\ \omega_z \end{pmatrix} \tag{2}$$

and equa­tion (1b) to rewrite (1a) as

$$\tilde A(\mathbf{k}) = \intop^\infty \!\!\!\!\!\!\!\! \intop_{-\infty} A(\mathbf{p}) \exp\bigl( i \mathbf{k} \cdot \mathbf{p} \bigr) dx\,dy \tag{3}$$

The cor­re­spond­ing Fouri­er syn­the­sis equa­tion is the inverse of (3):

$$A(\mathbf{q}) = \intop^\infty \!\!\!\!\!\!\!\! \intop_{-\infty} \tilde A(\mathbf{k}) \exp\bigl( -i \mathbf{k} \cdot \mathbf{q} \bigl) d\omega_x\,d\omega_y \tag{4}$$

That is, the field in the plane $z=0$ can be writ­ten as a super­po­si­tion of plane waves with ampli­tudes (and phas­es) $\tilde A(\mathbf{k})$ in the direc­tion of $\mathbf{k}$. Since for plane waves we have

$$\bigl|\mathbf{k}\bigr| = \frac{2\pi}{\lambda} \tag{5}$$

we can eas­i­ly deter­mine the miss­ing $\omega_z$ for any giv­en $\omega_x$ and $\omega_y$, and the lat­ter two suf­fice to unique­ly iden­ti­fy the direc­tion $\mathbf{k}$.

#### Angular Spectrum

Since $\mathbf{k}$ basi­cal­ly describes a direc­tion in 3D space (and a wave­length of course), the Fouri­er trans­form $\tilde A(\mathbf{k})$ of the ini­tial aper­ture field dis­tri­b­u­tion $A(\mathbf{p})$ is called the angu­lar spec­trum (instead of the fre­quen­cy spec­trum). We have to be a bit care­ful here since the angu­lar fre­quen­cy spec­trum usu­al­ly also refers to the Fouri­er trans­form of a time sig­nal, in terms of the angu­lar fre­quen­cies $\omega$ instead of $\nu$.

Fig­ure 1 shows schemat­i­cal­ly the aper­ture and its 2D Fouri­er trans­form. The amount of padding deter­mines the res­o­lu­tion in $k$-space of the trans­form while the sam­pling res­o­lu­tion of the source field deter­mines the range of the trans­form, see [2]. When using a source res­o­lu­tion of 2 sam­ples per $\lamb­da$, the range of the Fouri­er trans­form includes all for­ward prop­a­gat­ing waves (as shown in the fig­ure). For illus­tra­tion, orange cir­cles cor­re­spond to spe­cif­ic angles (rel­a­tive to the sur­face nor­mal) of prop­a­ga­tion of the plane waves. We will use these lat­er.

Fig­ure 1: Aper­ture in the source plane (top) and result­ing angu­lar spec­trum (log­a­rith­mic mag­ni­tude) for an aper­ture size 20λ. Orange cir­cles cor­re­spond to spe­cif­ic angles (rel­a­tive to the sur­face nor­mal) in the spec­trum. The ampli­tudes of non-prop­a­gat­ing waves (those hav­ing an imag­i­nary prop­a­ga­tion con­stant and/or angle greater than π/2) have been set to zero (black).

So, once we have deter­mined the angu­lar spec­trum $\tilde A(\mathbf{k})$ using (3), we can deter­mine the field $A(\mathbf{q})$ at any posi­tion $\mathbf{q}$ in 3D space by sim­ply prop­a­gat­ing to and then super­pos­ing the plane waves at $\mathbf{q}$. The nec­es­sary prop­a­ga­tion dis­tance (and thus phase incre­ment) is giv­en sim­ply by $\mathbf{k} \cdot \mathbf{q}$, as shown in Fig­ure 2. While this part is much eas­i­er for plane waves than it is for the spher­i­cal waves of Huy­gens and Fres­nel, it is more com­pli­cat­ed to obtain the ini­tial ampli­tudes of the waves in the Fouri­er approach. In any case, we have to solve inte­grals. For­tu­nate­ly, the much quick­er (dis­crete) numer­i­cal Fouri­er trans­form can be used if doing so care­ful­ly to avoid sub­sam­pling and/or alias­ing.

Fig­ure 2: Illus­tra­tion of the con­tri­bu­tion of a sin­gle angu­lar spec­trum com­po­nent k to the result­ing field at loca­tion q on a screen.

With (5) and using the def­i­n­i­tion of the scalar vec­tor prod­uct, we can write the expo­nent in (4) as

$$-i \mathbf{k} \cdot \mathbf{q} = -i \bigl|\mathbf{k}\bigr| \bigl|\mathbf{q}\bigr| \cos \alpha = -\frac{i 2 \pi}{\lambda} \bigl|\mathbf{q}\bigr| \cos \alpha \tag{6}$$

where $\alpha$ is the angle sub­tend­ed by the vec­tors $\mathbf{q}$ and $\mathbf{k}$ and $\bigl|\mathbf{q}\bigr|$ is the dis­tance of the obser­va­tion point from the ori­gin. Look­ing at the field dis­tri­b­u­tion on a sphere cen­tered on the ori­gin so that $\bigl|\mathbf{q}\bigr|$ is con­stant and equal to $q_{0}$, the expres­sion (6) is very sen­si­tive to even very small changes in $\alpha$ for $q_{0} \gg \lamb­da$. It will oscil­late rapid­ly if $\mathbf{q}$ and $\mathbf{k}$ are not par­al­lel. As a result of the near-uni­form dis­tri­b­u­tion of the result­ing phas­es, these con­tri­bu­tions to $A(\mathbf{q})$ in (4) tend to aver­age to zero and can be neglect­ed. Hence, the result­ing field on a sphere of a large radius $q_{0}$ is the two-dimen­sion­al Fouri­er trans­form of the source field with

$$A(\mathbf{q}) = \tilde A (\mathbf{q}) = \tilde A \biggl(\frac{q_0}{\bigl|\mathbf{k}\bigr|} \mathbf{k} \big­gr) = \tilde A \biggl(\frac{\lambda q_0}{2 \pi} \mathbf{k} \big­gr) \tag{7}$$

The pro­jec­tion of the field on a pla­nar screen is then only a triv­ial task. If $\mathbf{q}$ is very far from the source and the involved angles are very small (the so-called parax­i­al approx­i­ma­tion) so that $\omega_z \gg \omega_{x,y}$, we can neglect the cur­va­ture of the sphere near the opti­cal axis and write (10) direct­ly as

$$A(\mathbf{q}) = A\begin{pmatrix} x \\ y \\ z_0 \end{pmatrix} \approx \tilde A \biggl(\frac{\lambda z_0}{2 \pi} \mathbf{k} \big­gr) \tag{8}$$

Hence, at a large dis­tance from the dif­fract­ing aper­ture the result­ing field dis­tri­b­u­tion on a plane is described by the Fouri­er trans­form of the aper­ture field. This is the same result as that which Fraun­hofer obtained with his approx­i­ma­tion to the Fres­nel-Kirch­hoff inte­gral for large dis­tances between object and screen!

Fig­ure 3 shows the evo­lu­tion of the field mag­ni­tude at a point on the opti­cal axis in var­i­ous dis­tances from the source $S$ when increas­ing the por­tion of the angu­lar spec­trum used to com­pute it. Clear­ly, at larg­er dis­tances only plane waves prop­a­gat­ing at small angles to the direc­tion of the point on screen con­tribute to the field at that point, while close to the source the whole angu­lar spec­trum plays a role.

Fig­ure 3: Con­ver­gence behav­ior for the on-axis field in var­i­ous dis­tances from the source. The angle para­me­ter describes the largest angle of the plane waves to be includ­ed in the field com­pu­ta­tion (see also Fig­ure 1).

#### Imagery

Sim­i­lar to this post for the Huy­gens-Fres­nel dif­frac­tion and its descen­dants, we can numer­i­cal­ly cal­cu­late the inten­si­ty result­ing from dif­frac­tion at an aper­ture in var­i­ous dis­tances from the aper­ture. Since most of these plots look very sim­i­lar to the Huy­gens-Fres­nel results, only a small selec­tion of pat­terns is shown.

Fig­ure 4: Dif­frac­tion pat­tern for a pla­nar square aper­ture of 20λ edge length in a dis­tance of 10λ from the aper­ture, cal­cu­lat­ed from the angu­lar spec­trum. Bot­tom com­pares inten­si­ty to Huy­gens-Fres­nel solu­tion (exact).

Fig­ure 5: Dif­frac­tion pat­tern for a pla­nar square aper­ture of 20λ edge length in a dis­tance of 1000λ from the aper­ture, cal­cu­lat­ed from the angu­lar spec­trum. Bot­tom com­pares inten­si­ty to Huy­gens-Fres­nel solu­tion (exact).

If we look close­ly in the dif­frac­tion pat­tern, and more clear­ly in the cross sec­tion, there are notice­able fluc­tu­a­tions with­in the out­er fringes at $z = 1000\lambda$. These are a result of the rapid­ly oscil­lat­ing phase that the plane waves at large angles to the direc­tion of $\mathbf{q}$ have. If these oscil­la­tions are not sam­pled suf­fi­cient­ly dense­ly, such alias­ing occurs. This can be avoid­ed by either sam­pling the angu­lar spec­trum more dense­ly (by increas­ing the padding in the aper­ture field descrip­tion) or by dis­re­gard­ing the out­er plane waves (those with small $k_z$) in the field com­pu­ta­tion.

#### Discussion

Unfor­tu­nate­ly, the Fouri­er method only works for pla­nar aper­tures / sur­faces. We could define a three-dimen­sion­al Fouri­er trans­form by not let­ting $z=0$ in (1b), but then all three Carte­sian com­po­nents of $\mathbf{k}$ would already be defined by (3) and there would be no remain­ing degree of free­dom to ensure that (5) can be ful­filled.

We could also approx­i­mate an arbi­trary sur­face by many local­ly pla­nar sur­faces for which the Fouri­er decom­po­si­tion is pos­si­ble. How­ev­er, in the lim­it of infin­i­tes­i­mal sur­face ele­ments $dS$ the Fouri­er trans­form $\tilde{A}$ would be a con­stant for all direc­tions $\mathbf{k}$ and thus would be a (half) spher­i­cal wave – togeth­er with the required inte­gra­tion over all sur­face ele­ments we obtain the Fres­nel-Kirch­hoff dif­frac­tion for­mal­ism as lim­it­ing case and don’t even have to wor­ry about the coef­fi­cients. Wow.

#### Computation

Like before I’ll post the code I have used to cal­cu­late the inten­si­ty dis­tri­b­u­tions. The cross sec­tion method used in the oth­er post can be applied here, for the square aper­ture, also. Arbi­trary aper­tures are more com­pli­cat­ed, but quite pos­si­ble. Since the numer­i­cal inte­gra­tion nec­es­sary for the Huy­gens-Fres­nel approach becomes much slow­er for arbi­trary bound­aries, the angu­lar spec­trum approach may be faster in that case – I haven’t checked.
As always I take no respon­si­bil­i­ty for the cor­rect­ness of the code (or any­thing here, real­ly), but at least the results are con­sis­tent.

aperture = 20; % in multiples of lambda
padding = 100; % in multiples of lambda
resolution = 2; % samples per lambda

% construct square aperture
fftsize = (2*padding + aperture) * resolution;
Ap = zeros(fftsize);
Ap((padding * resolution + 1):((padding + aperture) * resolution),...
(padding * resolution + 1):((padding + aperture) * resolution)) = 1;
Ap = ifftshift(Ap); % necessary for DFT

% generate \mathbf{k}
kvector = linspace(-resolution/2, resolution/2, fftsize + 1);
kvector = kvector(1:length(kvector)-1); % remove last element

kx = repmat(kvector, [fftsize, 1]);
ky = repmat(kvector.', [1, fftsize]);
kz = sqrt(1 - kx.^2 - ky.^2);

% compute angular spectrum
spectrum = fftshift(fft2(Ap));
phaseshift = exp(-1i * kvector * pi/2);
phaseshift = repmat(phaseshift, [fftsize, 1]) .* repmat(phaseshift.', [1, fftsize]);
spectrum = spectrum .* phaseshift; % compensate the symmetry of Ap for the DFT;
% otherwise Aq will not be centered

spectrum(imag(kz) ~= 0) = 0; % disregard evanescent waves
spectrum = spectrum / sqrt(sum(sum(abs(spectrum).^2)) / sum(sum(Ap.^2))) / fftsize;
% power normalization (and Parseval)

% field on screen
z = 1000; % in multiples of lambda

qx = linspace(-200, 200, 501);
qy = linspace(-200, 200, 501);

for ii = 1:length(qx)
for jj = 1:length(qy)

Aq(ii,jj) = sum(sum(spectrum .* exp(1i * 2*pi * (kx * qx(ii) ...
+ ky * qy(jj) + kz * z))));

end
end


[1] Spa­tial fre­quen­cy [Wikipedia]
[2] DFT Cheat Sheet [Nerd Rage]

last posts in diffraction:

• Anonymous wrote:

Thanks for the won­der­ful­ly illus­tra­tive post! I’ve a ques­tion about your kvec­tor - why have you not includ­ed the 2*pi/lambda mul­ti­pli­er?

• Marcus wrote:

I’m not sure what you mean. Eq. (5) shows the rela­tion between k and λ. Did I miss some­thing?

• Hal­lo Markus,
ich habe eine Frage zu Fig­ure 1; Du sagst ja, dass du die Sam­plingfre­quenz ger­ade so wählst, dass der Range des Winkel­spek­trum ger­ade alle Vor­wärt­srich­tun­gen (also +/-pi/2) abdeckt. Das deckt sich auch mit deinen orangen Kreisen. Allerd­ings erk­lären sich mir nicht deine Achsenbeschrif­tung von -2pi - 2pi?

• All­ge­mein gilt wie in (5) beschrieben
$$\lamb­da \sqrt{ k_x^2 + k_y^2 + k_z^2 } = \lamb­da \left|\mathbf{k}\right| = 2\pi$$
Das heißt $\lamb­da \, k_x \le 2\pi$ (und eben­so für $k_y$ und $k_z$). An dieser Stelle sieht man, dass die Achsenbeschrif­tun­gen in der Tat falsch sind: es müsste $\lamb­da \, k_x$ anstatt $k_x / \lamb­da$ heißen. Ich habe die Graphen inzwis­chen kor­rigiert. Die Gren­zen $-2\pi$ und $2\pi$ sind aber richtig - vielle­icht wäre es ein­leuch­t­en­der gewe­sen, $k_x / \left|\mathbf{k}\right|$ zu schreiben und die Achsen von $-1$ bis $1$ laufen zu lassen. Let­z­tendlich ging es nur darum, die Wellen­länge her­auszunormieren…
Der Winkel $\theta$ zur Ober­flächen­nor­malen, der auch mit den orangen Kreisen dargestellt ist, wird hier durch $\cos \theta = \lamb­da \, k_z$ definiert, beschreibt also erst mal etwas ganz anderes.

• Dave wrote:

Dear Markus,
thank you for your fine exam­ple!
I am already famil­iar with the angu­lar wave spec­trum method, like Good­man described it. But I am real­ly inter­est­ed in your Fraun­hofer approx­i­ma­tion by angu­lar wavespec­trum.
I am hav­ing a hard time to under­stand the fol­low­ing sen­tence:
“As a result of the near-uni­form dis­tri­b­u­tion of the result­ing phas­es, these con­tri­bu­tions to A(q) in (4) tend to aver­age to zero and can be neglect­ed.”
- What are the “result­ing” phas­es?
- Why do they fol­low a near-uni­form dis­tri­b­u­tion?
- And how do they aver­age out?
Esspe­cial­ly why do they only aver­age out for $q_0»\lambda$, if they are uni­form dis­trib­uted shouldn’t they also aver­age out for clos­er dis­tances?
The back­round of my ques­tion is, that I have mul­ti slit prob­lem (slits stacked on top of each oth­er), I can eas­i­ly cal­cu­late the field prop­a­ga­tion between these slits, but for the last slit the field should be viewed from a greater dis­tance. So I am look­ing for a way to exam­ine when the Fraun­hofer approx­i­ma­tion is valid, if I only know the com­plex field dis­tri­b­u­tion in the upper aper­ture and the dis­tance!
Thank you!
Dave