Derivation of the Transfer Function and Impulse Response for Fresnel

Diffraction
Maxwell showed that the electric field must satisfy the wave equation:
1 ∂2E
∇2 E − =0
v2 ∂t2
which means that the individual components (“polarizations”) also must satisfy this equation:
n2 ∂ 2 Ej
∇2 Ej − =0
c2 ∂t2
where j = x, y, z specifies the component and v ≡ nc , where n is the refractive index. The complete solution
for the electric field is the sum of the three polarizations over the integral of fields for all values of k and
ω weighted by (generally complex-valued) weight function F [k, ω] = F [kx , ky , kz , ω], which is analogous to
the frequency spectrum in the Fourier transform:
3
X Z +∞ Z +∞ Z +∞ Z +∞
E [x1 , y1 , z1 , t] = âj dω dkx dky dkz (F [kx , ky , kz , ω] exp [i (kx x + ky y + kz z − ωt)])
j=1 −∞ −∞ −∞ −∞

where the summation over j allows us to assign any polarization to the electric field. From this point
forward we ignore the summation over polarizations, but the linearity of the equation allows us to add
any additional contributions at the end if needed. If we apply the well-known trial solution E [x, y, z, t] =
E0 exp [i (k • r − ωt)] to the wave equation in vacuum, we obtain two expressions that are related by the
wave equation:

∇2 E = ∇2 (E0 exp [i (k • r − ωt)])

µ 2 ¶
∂ ∂2 ∂2
= (E0 exp [−iωt]) + + exp [i (kx x + ky y + kz z)]
∂x2 ∂y 2 ∂z 2
¡ ¢
= (E0 exp [−iωt]) · kx2 + ky2 + kz2 exp [i (kx x + ky y + kz z)]
2
= i2 |k| E = −k 2 E
∂2E
= (−iω)2 E = −ω 2 E
∂t2
Thus the wave equation has the form yields the identity:
1 ∂2E ω2 ³ ω ´2
∇2 E − = −k 2
E − E = 0 =⇒ k 2
= =⇒ k = ω
c
v2 ∂t2 c2 c
This relation £ applies¤ a constraint to the general expression that may be expressed as the 1-D Dirac delta
function δ k − ωc . Thus the constrained solution for the electric field is:
Z +∞ Z +∞ Z +∞ Z +∞ h ωi
E [x1 , y1 , z1 , t] = dω dkx dky dkz (F [kx , ky , kz , ω] exp [i (kx x + ky y + kz z − ωt)]) δ k −
−∞ −∞ −∞ −∞ c

If we assume that the light is monochromatic with frequency ω 0 (wavelength λ0 = 2πc

ω0 ), then the integral over
ω may be evaluated immediately by setting the temporal frequency dependence to a Dirac delta function:

F [kx , ky , kz , ω] = F [kx , ky , kz ] · δ [ω − ω 0 ]

which leads to the expression for the electric field::

E [x1 , y1 , z1 , t]
Z +∞ Z +∞ Z µZ h ¶
+∞ +∞
ωi
= dkx dky dkz
dω F [kx , ky , kz ] · δ [ω − ω 0 ] δ k − exp [i (kx x + ky y + kz z − ωt)]
−∞ −∞ −∞ −∞ c
Z +∞ Z +∞ Z +∞ h ω0 i
= exp [−iω 0 t] dkx dky dkz (F [kx , ky , kz ] exp [i (kx x + ky y + kz z)]) δ k −
−∞ −∞ −∞ c

1
where the sifting property of the Dirac delta function has been used. This remaining 1-D Dirac delta function
reduces the 3-D integral to 2-D. This means that the 3-D representation of E for one wavelength may be
evaluated from a 2-D function. To evaluate ththat function, we need to integrate over one of the components
of k. We’ll select the integral of kz , because we will define the source
h¡ distribution ¢in1 the 2-D
i x − y plane.
£ ω0
¤ 2 2 2 2 ω0
Therefore, we must put the 1-D Dirac delta function δ k − c = δ kx + ky + kz − c in the form of
δ [kz − (kz )0 ], where (kz )0 is the positive number that satisfies:

³ ´ 12 ω0 h 2 ¡ ¢i 12
ω
kx2 + ky2 + (kz )20 − = 0 =⇒ (kz )0 = c20 − kx2 + ky2
c

In other words, the argument of the 1-D Dirac delta function is a function of kz , and we must recall that
expression. In this case, we have:
h ω0 i h¡ ¢ 1 ω0 i
δ k− = δ kx2 + ky2 + kz2 2 − = δ [g [kz ]]
c c
¡ ¢ 1 ω0
where g [kz ] ≡ kx2 + ky2 + kz2 2 −
c

Dirac Delta Function of Functional Argument

Recall that if g [x0 = 0], then:
δ [x − x0 ]
δ [g [x]] = ¯¯ ¯¯
dg
¯ dx ¯
x=x0
dg
so we need to evaluate the derivative dkz using the chain rule:

dg 1¡ 2 ¢− 1 kz
= kx + ky2 + kz2 2 · 2kz = ¡ ¢1
dkz 2 kx2 + ky2 + kz2 2

At the root (kz )0 , this derivative takes the value:

¯
dg ¯¯ (kz )0
¯ = ³ ´ 12
dkz kz =(kz )
0 kx2 + ky2 + (kz )20
h 2 ¡ ¢i 12
ω0 2 2
c2 − k x + k y
= ³ h 2 ¡
ω ¢i´ 12
kx2 + ky2 + c20 − kx2 + ky2
h 2 ¡ ¢i 12
ω0 2 2
c2 − k x + k y
= ³ 2 ´ 12
ω0
c2
h ¡ ¢i 12
ω20 ∙ ¸1
c2 − kx2 + ky2 c2 ¡ 2 ¢ 2
¡ ω0 ¢ 2
= = 1 − 2 kx + ky
c
ω0

Thus the expression for the 1-D Dirac delta function in terms of (kz )0 is:
h ω0 i δ [k − (kz )0 ]
δ k− = ¯ z¯ ¯
c ¯ dg ¯ ¯
¯ ¯
¯ dkz ¯kz =(kz ) ¯
0
" µ ¶ 12 #
1 ω 20 ¡ 2 ¢
= h δ kz − − kx + ky2
c2
¡ ¢i 12 c2
1− ω20
kx2 + ky2

2
Evaluate the Integral over kz
When evaluating the 1-D integral over kz , the 1-D Dirac delta function allows us to substitute (kz )0 =
³ 2 ¡ ¢´ 12
ω0 2 2
c2 − kx + ky at every appearance of kz :
Z h
+∞
ω0 i
dkz (F [kx , ky , kz ] exp [i (kx x + ky y + kz z)]) δ k −
−∞ c
⎛ ⎞
Z " µ ¶ 12 #
+∞ ¡ ¢
⎜ 1 ω 20 ⎟
= dkz (F [kx , ky , kz ] exp [i (kx x + ky y + kz z)]) · ⎝ h δ kz − − kx2 + ky2 ⎠
−∞ 2 ¡ ¢i 12 c2
1 − ωc 2 kx2 + ky2
0
⎛ ⎞
" Ã µ 2 ¶ 1 !#
⎜ 1 ⎟ ω0 ¡ 2 2
¢ 2
= ⎝F [kx , ky , (kz )0 ] h exp i kx x + ky y + − kx + ky z
2 ¡ ¢i 12 ⎠ c2
1 − ωc 2 kx2 + ky2
0
⎛ ⎞
" µ ¶ 12 #
⎜ 1 ⎟ ω 20 ¡ 2 ¢
= ⎝F [kx , ky , (kz )0 ] h exp [i (kx x + ky y)] exp i − kx + ky2 z
c2
¡ ¢i 12 ⎠ c2
2
1 − ω 2 kx + ky 2
0

Note that the only dependence on z appears in the last term. To simplify (and shorten!) this expression, we
rename part of the integrand:
" µ 2 ¶1 #
ω0 ¡ 2 2
¢ 2 1
A [kx , ky ; z] ≡ F kx , ky , − kx + ky h exp [i (kz )0 z]
c2 2 ¡ ¢i 12
1 − ωc 2 2
kx + ky 2
0
" µ 2 ¶ 12 # " µ ¶1 #
ω0 ¡ 2 ¢ 1 ω 20 ¡ 2 ¢ 2
= F kx , ky , − kx + ky2 h exp i 2
− kx + ky z
c2
2 ¡ ¢i 12 c2
1 − ωc 2 kx2 + ky2
0

The same function evaluated at z = 0 is:

" µ ¶ 12 #
ω 20 ¡ 2 ¢ 1
A [kx , ky ; z = 0] = F kx , ky , − kx + ky2 h
c2 c2
¡ ¢i 12
1− ω 20
kx2 + ky2

The entire integral now may be written as:

µZ +∞ Z +∞ ¶
E [x1 , y1 , z1 , t] = exp [−iω 0 t] · dkx dky A [kx , ky ; z] exp [i (kx x + ky y)]
−∞ −∞
ÃZ Z Ã " µ ¶ 1 #! !
+∞ +∞
ω 20 ¡ 2 2
¢ 2
= exp [−iω 0 t] · dkx dky A [kx , ky ; 0] exp i − kx + ky z exp [i (kx x + ky y)]
−∞ −∞ c2

ω0 2π 2π
which has the form of an inverse 2-D Fourier transform if we identify c = λ0 , k = |k| = λ , kx = 2πξ, and
ky = 2πη. This leads to:
" µ ¶1 #
ω 20 ¡ 2 2
¢ 2
exp [i (kz )0 z] = exp i − kx + ky z
c2
∙ ¸
2π ¡ 2
¡¡ 2 2
¢¢¢ 12
= exp i 1−λ ξ +η z
λ0

and:
exp [i (kx x + ky y)] = exp [+2πi (ξx + ηy)]

3
The exact expression for the integral becomes:
Z +∞ Z +∞ µ ∙ ¸¶
2π ¡ ¡ ¢¢ 1
E [x1 , y1 , z1 , t] = (2π)2 exp [−iω 0 t] dξ dη A [ξ, η; 0] · exp i 1 − λ20 ξ 2 + η2 2 z exp [+2πi (ξx + ηy)]
−∞ −∞ λ0

We will ignore the temporal phase term exp [−iω 0 t] from this point forward, because it contributes no spatial
dependence. The corresponding expression in the frequency domain is:
µ ∙ ¸¶
2π ¡ ¡ ¢¢ 1
F2 {E [x1 , y1 , z1 , t]} = A [ξ, η; 0] · exp i 1 − λ20 ξ 2 + η2 2 z
λ0

Approximating the Dependence on z

The task now is to make approximations that allow the integral to be h evaluated more readily. iThe obvious
2π
¡ 2
¡ 2 ¢¢ 1
candidate for this treatment is to remove the square root from exp i λ0 1 − λ0 ξ + η2 2 z . Recall the
binomial theorem:
n n (n − 1) 2
(1 + ε) = 1 + nε + ε + ···
2!
If |ε| << 1, then this simplifies to:
(1 + ε)n ∼= 1 + nε
£ ¤1
If the radial spatial frequency ρ = ξ 2 + η2 2 is small compared to the reciprocal of the wavelength (ρ <<
¡ 2 ¢
λ−1 1
0 ), then we can substitute n = 2 and ε = −λ
2
ξ + η2 :
∙ ¸ ∙ µ ¶ ¸
2πi ¡ 2
¡ 2
¢¢ 12
∼ 2πi 1 2¡ 2 2
¢
exp 1 − λ0 ξ + η z = exp 1 − λ0 ξ + η z
λ0 λ0 2
∙ ¸
2πz £ ¡ ¢¤
= exp i · exp −iπλ0 z ξ 2 + η2
λ0

Thus we have the expression:

µ ∙ ¸ ¶
2πz £ ¡ ¢¤
F2 {E [x1 , y1 , z1 , t]} ∼
= A [2πξ, 2πη; 0] · exp i · exp −iπλ0 z ξ 2 + η 2 if ξ 2 + η2 << λ−2
0
λ0

Transfer Function of Light Propagation

This expression has the form of a modulation in the frequency domain, so we can identify the multiplicative
term as a transfer function:
h i £ ¡ ¢¤
H [ξ, η; z1 ] = exp i 2πz
λ0
1
· exp −iπλ0 z1 ξ 2 + η2

which describes the effect of light propagation on the field evaluated at z = 0 for ρ << λ−1
0 , which defines
in the Fresnel diffraction region.

4
Impulse Response of Light Propagation
We can then use the known Fourier transform to evaluate the corresponding impulse response in the Fresnel
diffraction region.
∙ ¸
2πz1 © £ ¡ ¢¤ª
h [x, y; z1 ] = exp i · F2 exp −iπλ0 z1 ξ 2 + η2
λ0
∙ ¸
2πz1 © £ ¤ª © £ ¤ª
= exp i · F1 exp −iπλ0 z1 ξ 2 · F1 exp −iπλ0 z1 η2
λ0
∙ ¸ µ h πi ∙ ¸¶ µ h πi ∙ ¸¶
2πz1 1 x2 1 y2
= exp i · √ exp −i exp +iπ · √ exp −i exp +iπ
λ0 λ0 z1 4 λ0 z1 λ0 z1 4 λ0 z1
∙ ¸ µ ¶2 ³ " ¡ ¢ #
2πz1 1 h i´
π 2 x2 + y 2
= exp i · √ exp −i exp +iπ
λ0 λ0 z1 4 λ0 z1
∙ ¸ " ¡ 2 ¢ #
2πz1 1 x + y2
= exp i · (−i) exp +iπ
λ0 λ0 z1 λ0 z1

This the impulse response of light propagation in the Fresnel diffraction region:
h i ∙ ¸
(x2 +y2 )
h [x, y; z1 ] = exp i 2πz
λ0
1
· 1
iλ0 z1 exp +iπ λ0 z1

In other words, the electric field at the plane z = z1 , which we will call g [x, y; z1 ], can be derived from the
field at the plane z = 0, which we will call f [x, y; 0] via the 2-D convolution:
µ ∙ ¸¶ à " ¡ 2 2
¢ #!
1 2πz 1 x + y
g [x, y; z1 ] ∼
= exp i f [x, y; 0] ∗ exp +iπ
iλ0 z1 λ0 λ0 z1
³ h i´ µ ∙ ¸¶
(x2 +y2 )
g [x, y; z1 ] ∼
= λ01z1 exp i 2πzλ0
1
− i π
2 f [x, y; 0] ∗ exp +iπ λ0 z1

Interpretation of Fresnel Diffraction

Since we can evaluate the equation for light propagation from a 2-D distribution of monochromatic light
in the plane z = 0 as a convolution, then the process of light propagation in this approximation must be
linear and shift invariant. The impulse response is proportional to a quadratic-phase factor, and thus is an
“allpass” filter. The leading factor of z1−1 is due to the inverse square
h lawifor light propagating “straight
down” the z-axis from z = 0 to z = z1 . The constant-phase factor exp i 2πz λ0
1
represents the change in phase
of light along this same path. The approximation assumes that the effect of the inverse square law for light
that travels at any angle relative to the z-axis is too small to be considered. The spherical waves emitted by
the source(s) in the plane [x, y; 0] propagate sufficiently far so that the variation in amplitude with transverse
distance [x, y] is ignored (i.e., the inverse square law only and the spherical wavefront is approximated by a
paraboloidal wavefront.
The Fourier transform of the quadratic-phase factor demonstrates that the transver function is a “downchirp,”
a quadratic-phase factor with a negative sign.