## Scalar Diffraction – Fourier

In an earlier post we had a look at various diffraction formalisms that either originated in the Huygens-Fresnel principle or led to essentially the same results. The principle modeled a field distribution in a aperture (or on a surface) as a source of infinitely many spherical waves whose amplitudes and phases were prescribed by the aperture field. However, this spherical wave decomposition is not the only way to look at the problem. Provided that the aperture / surface is two-dimensional (flat), we can use Fourier analysis to decompose any field on this surface into a number of plane waves, which can then be propagated to any other surface where they interfere (like in the Huygens-Fresnel principle).

#### Going Fourier

Similar to the frequency decomposition of a time signal using the Fourier transform, one can use the two-dimensional Fourier transform to decompose a two-dimensional field in some aperture at $z=0$ into its spatial frequencies:

$$\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 spatial frequencies $\omega_x$ and $\omega_y$ are somewhat abstract quantities. This is just Fourier analysis formalism. Since $x$ and $y$ (and $z$) are position quantities, the corresponding Fourier transform quantities must actually be (angular) wavenumbers [1]. Hence we can use the wavevector,

$$\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 equation (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 corresponding Fourier synthesis equation 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 written as a superposition of plane waves with amplitudes (and phases) $\tilde A(\mathbf{k})$ in the direction of $\mathbf{k}$. Since for plane waves we have

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

we can easily determine the missing $\omega_z$ for any given $\omega_x$ and $\omega_y$, and the latter two suffice to uniquely identify the direction $\mathbf{k}$.

#### Angular Spectrum

Since $\mathbf{k}$ basically describes a direction in 3D space (and a wavelength of course), the Fourier transform $\tilde A(\mathbf{k})$ of the initial aperture field distribution $A(\mathbf{p})$ is called the angular spectrum (instead of the frequency spectrum). We have to be a bit careful here since the angular frequency spectrum usually also refers to the Fourier transform of a time signal, in terms of the angular frequencies $\omega$ instead of $\nu$.

Figure 1 shows schematically the aperture and its 2D Fourier transform. The amount of padding determines the resolution in $k$-space of the transform while the sampling resolution of the source field determines the range of the transform, see [2]. When using a source resolution of 2 samples per $\lambda$, the range of the Fourier transform includes all forward propagating waves (as shown in the figure). For illustration, orange circles correspond to specific angles (relative to the surface normal) of propagation of the plane waves. We will use these later.

Figure 1: Aperture in the source plane (top) and resulting angular spectrum (logarithmic magnitude) for an aperture size 20λ. Orange circles correspond to specific angles (relative to the surface normal) in the spectrum. The amplitudes of non-propagating waves (those having an imaginary propagation constant and/or angle greater than π/2) have been set to zero (black).

So, once we have determined the angular spectrum $\tilde A(\mathbf{k})$ using (3), we can determine the field $A(\mathbf{q})$ at any position $\mathbf{q}$ in 3D space by simply propagating to and then superposing the plane waves at $\mathbf{q}$. The necessary propagation distance (and thus phase increment) is given simply by $\mathbf{k} \cdot \mathbf{q}$, as shown in Figure 2. While this part is much easier for plane waves than it is for the spherical waves of Huygens and Fresnel, it is more complicated to obtain the initial amplitudes of the waves in the Fourier approach. In any case, we have to solve integrals. Fortunately, the much quicker (discrete) numerical Fourier transform can be used if doing so carefully to avoid subsampling and/or aliasing.

Figure 2: Illustration of the contribution of a single angular spectrum component k to the resulting field at location q on a screen.

With (5) and using the definition of the scalar vector product, we can write the exponent 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 subtended by the vectors $\mathbf{q}$ and $\mathbf{k}$ and $\bigl|\mathbf{q}\bigr|$ is the distance of the observation point from the origin. Looking at the field distribution on a sphere centered on the origin so that $\bigl|\mathbf{q}\bigr|$ is constant and equal to $q_{0}$, the expression (6) is very sensitive to even very small changes in $\alpha$ for $q_{0} \gg \lambda$. It will oscillate rapidly if $\mathbf{q}$ and $\mathbf{k}$ are not parallel. As a result of the near-uniform distribution of the resulting phases, these contributions to $A(\mathbf{q})$ in (4) tend to average to zero and can be neglected. Hence, the resulting field on a sphere of a large radius $q_{0}$ is the two-dimensional Fourier transform of the source field with

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

The projection of the field on a planar screen is then only a trivial task. If $\mathbf{q}$ is very far from the source and the involved angles are very small (the so-called paraxial approximation) so that $\omega_z \gg \omega_{x,y}$, we can neglect the curvature of the sphere near the optical axis and write (10) directly 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} \biggr) \tag{8}$$

Hence, at a large distance from the diffracting aperture the resulting field distribution on a plane is described by the Fourier transform of the aperture field. This is the same result as that which Fraunhofer obtained with his approximation to the Fresnel-Kirchhoff integral for large distances between object and screen!

Figure 3 shows the evolution of the field magnitude at a point on the optical axis in various distances from the source $S$ when increasing the portion of the angular spectrum used to compute it. Clearly, at larger distances only plane waves propagating at small angles to the direction of the point on screen contribute to the field at that point, while close to the source the whole angular spectrum plays a role.

Figure 3: Convergence behavior for the on-axis field in various distances from the source. The angle parameter describes the largest angle of the plane waves to be included in the field computation (see also Figure 1).

#### Imagery

Similar to this post for the Huygens-Fresnel diffraction and its descendants, we can numerically calculate the intensity resulting from diffraction at an aperture in various distances from the aperture. Since most of these plots look very similar to the Huygens-Fresnel results, only a small selection of patterns is shown.

Figure 4: Diffraction pattern for a planar square aperture of 20λ edge length in a distance of 10λ from the aperture, calculated from the angular spectrum. Bottom compares intensity to Huygens-Fresnel solution (exact).

Figure 5: Diffraction pattern for a planar square aperture of 20λ edge length in a distance of 1000λ from the aperture, calculated from the angular spectrum. Bottom compares intensity to Huygens-Fresnel solution (exact).

If we look closely in the diffraction pattern, and more clearly in the cross section, there are noticeable fluctuations within the outer fringes at $z = 1000\lambda$. These are a result of the rapidly oscillating phase that the plane waves at large angles to the direction of $\mathbf{q}$ have. If these oscillations are not sampled sufficiently densely, such aliasing occurs. This can be avoided by either sampling the angular spectrum more densely (by increasing the padding in the aperture field description) or by disregarding the outer plane waves (those with small $k_z$) in the field computation.

#### Discussion

Unfortunately, the Fourier method only works for planar apertures / surfaces. We could define a three-dimensional Fourier transform by not letting $z=0$ in (1b), but then all three Cartesian components of $\mathbf{k}$ would already be defined by (3) and there would be no remaining degree of freedom to ensure that (5) can be fulfilled.

We could also approximate an arbitrary surface by many locally planar surfaces for which the Fourier decomposition is possible. However, in the limit of infinitesimal surface elements $dS$ the Fourier transform $\tilde{A}$ would be a constant for all directions $\mathbf{k}$ and thus would be a (half) spherical wave – together with the required integration over all surface elements we obtain the Fresnel-Kirchhoff diffraction formalism as limiting case and don’t even have to worry about the coefficients. Wow.

#### Computation

Like before I’ll post the code I have used to calculate the intensity distributions. The cross section method used in the other post can be applied here, for the square aperture, also. Arbitrary apertures are more complicated, but quite possible. Since the numerical integration necessary for the Huygens-Fresnel approach becomes much slower for arbitrary boundaries, the angular spectrum approach may be faster in that case – I haven’t checked.
As always I take no responsibility for the correctness of the code (or anything here, really), but at least the results are consistent.

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 = 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] Spatial frequency [Wikipedia]
[2] DFT Cheat Sheet [Nerd Rage]

last posts in diffraction:

• Anonymous wrote:

Thanks for the wonderfully illustrative post! I’ve a question about your kvector – why have you not included the 2*pi/lambda multiplier?

• Marcus wrote:

I’m not sure what you mean. Eq. (5) shows the relation between k and λ. Did I miss something?

• Hallo Markus,
ich habe eine Frage zu Figure 1; Du sagst ja, dass du die Samplingfrequenz gerade so wählst, dass der Range des Winkelspektrum gerade alle Vorwärtsrichtungen (also +/-pi/2) abdeckt. Das deckt sich auch mit deinen orangen Kreisen. Allerdings erklären sich mir nicht deine Achsenbeschriftung von -2pi – 2pi?

• Allgemein gilt wie in (5) beschrieben
$$\lambda \sqrt{ k_x^2 + k_y^2 + k_z^2 } = \lambda \left|\mathbf{k}\right| = 2\pi$$
Das heißt $\lambda \, k_x \le 2\pi$ (und ebenso für $k_y$ und $k_z$). An dieser Stelle sieht man, dass die Achsenbeschriftungen in der Tat falsch sind: es müsste $\lambda \, k_x$ anstatt $k_x / \lambda$ heißen. Ich habe die Graphen inzwischen korrigiert. Die Grenzen $-2\pi$ und $2\pi$ sind aber richtig – vielleicht wäre es einleuchtender gewesen, $k_x / \left|\mathbf{k}\right|$ zu schreiben und die Achsen von $-1$ bis $1$ laufen zu lassen. Letztendlich ging es nur darum, die Wellenlänge herauszunormieren…
Der Winkel $\theta$ zur Oberflächennormalen, der auch mit den orangen Kreisen dargestellt ist, wird hier durch $\cos \theta = \lambda \, k_z$ definiert, beschreibt also erst mal etwas ganz anderes.

• Dave wrote:

Dear Markus,
thank you for your fine example!
I am already familiar with the angular wave spectrum method, like Goodman described it. But I am really interested in your Fraunhofer approximation by angular wavespectrum.
I am having a hard time to understand the following sentence:
“As a result of the near-uniform dis­tri­b­u­tion of the result­ing phases, these con­tri­bu­tions to A(q) in (4) tend to aver­age to zero and can be neglected.”
– What are the “resulting” phases?
– Why do they follow a near-uniform distribution?
– And how do they average out?
Esspecially why do they only average out for $q_0>>\lambda$, if they are uniform distributed shouldn’t they also average out for closer distances?
The backround of my question is, that I have multi slit problem (slits stacked on top of each other), I can easily calculate the field propagation between these slits, but for the last slit the field should be viewed from a greater distance. So I am looking for a way to examine when the Fraunhofer approximation is valid, if I only know the complex field distribution in the upper aperture and the distance!
Thank you!
Dave