# Interpretation Fourier Diffraction

I got familiar with interference and now want to understand it with the help of Fourier optics, as it is called.

For the derivation of the diffraction integral the following setup is to consider:

Relating to the picture (respectively the lecture based on it) the $$\textsf {Electric field interference pattern}$$ $$E(x,y)$$, that forms when light is diffracted at an object, described by an $$\textsf{aperture function}$$ $$g(u,v)$$ (here a 2-D-rectangle) can be calculated simply by:

$$\displaystyle E(x,y) = c\,\int\int g(u,v)\,\dfrac{e^{i\,k\,r}}{r} \,\mathrm{du\,dv}$$ where $$r$$ is the distance of a pointsource to a point on the diffraction pattern:

$$r = \sqrt{(x-u)^2+(y-v)^2+d^2}$$.

Intuitively this makes sense: by integrating you're summing up all the pointy waves phases modified by the aperture function and wonder, what the outcome looks like.

By some parabolic approximations (large distance) this transfers over to $$\textsf{Fraunhofer Diffraction}$$:

$$\displaystyle E(x,y) = c\,\dfrac{e^{i\,k\,d}}{d}\,e^{i\,\frac{k}{2\,d}(x^2+y^2)}\,\int\int g(u,v)\,e^{-i\,\frac{k}{d}(x\,u+y\,v)}\,\mathrm{du\,dv}$$

Now, here I have to ask: where does a Fourier transform come into play? The Integral could be considered as such, but what is it transforming? I just see a function of position $$(u,v)$$ staying a function of position $$(x,y)$$

Here's a more precise example, what I can't associate though. It's 1-D diffraction on a single slit:

In this specific case the aperture function is set by $$g(x) = \begin{cases} 1 & |x| < \frac{a}{2} \\ 0 & |x| > \frac{d}{2}\end{cases}$$ as you expect a slit to act. Now, without any consideration taking the Fourier transform shifting over to the wave vector in x-direction:

$$\displaystyle E(k_x) = C\,\int_{-a/2}^{a/2} 1\,e^{-i\,k\,x}\,\mathrm{dx} = C\,\dfrac{2\,\sin(k_x\,\frac{a}{2})}{k_x}$$ and using $$k_x = k\,\sin(\alpha)$$:

$$E(\alpha) = C\,\dfrac{2\,\sin(k\,\sin(\alpha)\,\frac{a}{2})}{k\,\sin(\alpha)}$$

Now this is exactly the Electric field pattern lectures are talking about. However how is this connected to the $$\textsf{Fraunhofer Approximation}$$? To what kind of space the Fourier Transform was shunting. Position to wave vector and then to angle?

The thing is I only know Fourier from time to frequency domain. (Even if it still troubles me). Maybe you can help me out for this extra example.

You start out with the Helmholtz equation that holds true in a homogeneous lossless medium for any rectangular component of both the $$E$$ and $$H$$ fields: $$\nabla^2 U + \kappa^2 U = 0 \tag{1}\label{1}$$ That is $$U$$ may stand for any of the $$E_x, E_y,...,H_z$$. A general solution of $$\eqref{1}$$ is the plane wave $$e^{\mathfrak {j} \kappa \hat {\mathbf{k}} \cdot \mathbf{r}}$$ where $$\hat {\mathbf{k}}$$ is the unit vector in the direction of propagation and $$\kappa = \omega/c =2\pi/\lambda$$.

Now assume that a wave, say $$U(x,y,z)$$ propagates *mostly* along the $$z$$ axis (paraxial approximation) and is incident upon the $$z=0$$ plane. Decompose this field at an arbitrary $$z$$ into plane wave components using the 2D Fourier transform: $$U(x,y,z)=\int \int dudv A_z(u,v) e^{\mathfrak j \kappa(xu+yv)}\tag{2}\label{2}$$ Here $$dU(x,y,z)= A_z(u,v) e^{\mathfrak j \kappa (xu+yv)} dudv$$ is a differential plane wave moving in the direction represented by the direction cosines $$u,v,w$$ where $$u^2+v^2+w^2=1$$, i.e., $$w=\sqrt{1-u^2-v^2}$$.

Next we wish to find the relationship between the field (component) $$U(x,y,z)$$ at $$z=0$$ and at arbitrary $$z$$. From linearity $$A_z(u,v) = A_0(u,v) e^{\mathfrak{j}\kappa w}$$. The inverse Fourier transform gives: $$A_0(x,y)=b_1\int \int dx'dy' U(x',y',0) e^{-\mathfrak j \kappa(x'u+y'v)}\tag{3}\label{3}$$ for some constant $$b$$ (in fact $$b_1=1/(2\pi)^2$$ Again because of the linear shift invariant nature of hte optical (EM) system the *output* field at $$z$$ is a convolution of the *input* field at $$z=0$$ and kernel, say, $$G(x,y)$$, that is $$U(x,y,z) =\int \int dx'dy' U(x',y',0) G(x-x',y-y')\tag{4}\label{4}$$ and $$G(x,y) = \int \int du dv e^{\mathfrak j \kappa w} e^{\mathfrak j \kappa (ux+vy)} \tag{5}\label{5}$$ where $$w=\sqrt{1-u^2-v^2}$$. When polar coordinates are introduced $$u=R cos\theta$$ and $$v=R sin \theta$$ after much pain we get (see the Erdelyi Higher Transcendental Functions, 4.52) $$G(R cos\theta,R sin \theta) = \frac{1}{\mathfrak j \lambda}\frac{e^{\mathfrak j \kappa \sqrt{R^2+z^2}}}{\sqrt{R^2 + z^2}} \frac{z}{\sqrt{R^2 + z^2}}\left( 1-\frac{1}{\mathfrak j \kappa \sqrt{R^2 + z^2}}\right) \tag{6}\label{6}$$

Now in the far field (Fraunhofer) we have

1. when $$z \gg R$$ the term in parentheses is approximated as $$1$$
2. Define $$cos\phi = \frac{z}{\sqrt{R^2+z^2}}$$, this angle $$\phi$$ is the direction the elementary wavelet has with the axis $$z$$ as it propagates, and $$cos\phi$$ is the obliquity factor introduced first by Fresnel. Taking $$\phi \approx 0$$ that is $$cos\phi \approx 1$$ is the paraxial approximation.
3. The term $$\sqrt{R^2+z^2}$$ is the distance between the point $$x',y'$$ in the $$z=0$$ plane and the field point $$x,y,z$$, and in the paraxial regime this can be approximated as $$z$$.

Therefore $$G(R cos\theta,R sin \theta)\approx \tilde G(R) = \frac{1}{\mathfrak j \lambda}\frac{e^{\mathfrak j \kappa \sqrt{R^2+z^2}}}{z} \tag{7}\label{7}$$ This gives us the field at $$z$$ approximately by calcuating the convolution integral as $$U(x,y,z)\approx \frac{b_2}{z}\int \int dx'dy'U(x',y',0)e^{\mathfrak j \kappa \sqrt{R^2+z^2}} \tag{8}\label{8}$$

Finally, in the Fraunhofer far-field we also have the approximation $$\sqrt{R^2+z^2} = \sqrt{z^2+(x-x')^2+(y-y')^2} \\ \approx z+ \frac{1}{2z}(x-x')^2+\frac{1}{2z}(y-y')^2,$$ and thus

$$U(x,y,z)\approx \frac{b_2 e^{\mathfrak j \kappa z}}{z}\int \int dx'dy'U(x',y',0)e^{\mathfrak j \kappa (\frac{1}{2z}(x-x')^2+\frac{1}{2z}(y-y')^2)} \tag{9}\label{9}$$

The result $$\eqref{9}$$ is the diffraction integral, the input (incident) field $$U(x,y,0)$$is at $$z=0$$ represents any component of $$E$$ or $$H$$ and is supposed to be known. It propagates in the free half-space $$z>0$$ and via this convolution integral one can calculate the same field component at any point $$x,y,z$$.

If you wish so you can also write the $$\eqref{9}$$ convolution integral as a Fourier transform by expanding the quadratic terms in the exponent:

$$U(x,y,z)\approx \frac{b_2 e^{\mathfrak j \kappa (z+x^2/(2z) +y^2/(2z))}}{z}\int \int dx'dy' e^{\mathfrak j \tfrac{\kappa}{2z} (x'^2 +y'^2)} U(x',y',0)e^{-\mathfrak j \kappa (xx'+yy')/z} \tag{10}\label{10}$$

but notice that $$\eqref{10}$$ is not the Fourier transform of $$U(x,y,0)$$ but instead there is a quadratic phase factor also involved, so that aside from some unit modulus phase factor outside the integral this is the Fourier transform of $$e^{\mathfrak j \kappa(x^2 +y^2)/(2z)} U(x,y,0)$$. Now if the extent of the diffracting object at $$z=0$$ is such that $$\tfrac{\kappa}{2z}(x^2 +y^2) \ll 1$$ or $$x^2 +y^2 \ll z \lambda$$ then the quadratic phase factor can be neglected and we do have the Fourier transform of the incident field.

[1] Stark: Application of Optical Fourier Transforms, Academic Press 1982

• wow, really detailed derivation. To be honest some parts of it I don't understand at first, but it gives me the feeling there is much more behind it to back up. After I had made some additional research I figured $x = k_x\,d/k$ what explains how to convert from position to angle (for small ones). But as said, your huge elaboration provides the whole background story how to truly get to the diffraction integral itself. Hats off.
– Leon
Commented Aug 23, 2021 at 12:19
• This spectral decomposition derivation is preferred by engineers rather than physicists who almost always start with Kirchhoff's integral formula relating the surface values and derivatives of the incident field to the enclosed field but that has the problem how to assign the surface field for the values and their derivatives are not independent of each other and we have no idea what they are. Commented Aug 23, 2021 at 14:33
• The Kirchhoff (Sommerfeld, Runge, Debye, etc.) method, see en.wikipedia.org/wiki/Kirchhoff%27s_diffraction_formula, is the standard but this (Duffieux, Clemmow, Michelson (?)) is faster and is easier to understand the various approximation steps physically, but you need to know some amount of Fourier analysis before (linear systems, convolution, etc.). Commented Aug 23, 2021 at 14:34
• Oh, Kirchhoff-Diffraction. There it is originating from. For me the near-field approximations are sufficient. Concerning those I've got an additional question: It's being said the field pattern $U(x,y)$ is beging calculated like $B \int \int U(x',y')\,e^{\textstyle i\,\frac{k}{2\,z}\,\left[(x-x')^2+(y-y')^2\right]}$. In fact this is a convolution: $U(x,y) = B\,U(x,y)\ast h(x,y)$ where $h(x,y) = e^{\textstyle i\,\frac{k}{2\,z}\,\left[x +y\right]}$. Now this should be much easier to compute in Fourier space.
– Leon
Commented Aug 23, 2021 at 20:20
• Your $h(x,y)$ needs "squares" in its exponent, with that it is a convolution $U\otimes h$ but the way to calculate it is by expanding the exponential (Eq.9) and that step converts it automatically into a true Fourier transform (Eq.10) but with the quadratic phase modulation on $U$. In the Fresnel region besides the phase modulation you probably want to keep the denominator $\sqrt{R^2+z^2}$ as well not its approximant $z$. At any rate, it being a true Fourier transform after the exponent is expanded you can calculate it via FFT. Commented Aug 23, 2021 at 20:35