Thin Lens - Fourier Optics     H. Q. Le  copyrighted

ECE 5368/6358                              han q le - copyrighted
Use solely for students registered for UH ECE 6358/5368 during courses - DO NOT distributed (copyrighted materials)

Segment 3 - Grating devices

2. Application to grating spectrometer

2.0 Introduction

http://scienceblogs.com/startswithabang/2012/10/24/what-is-the-sun-made-out-of/

2.1 Input half

2.1.1 Impulse response function (delta function input).

If we have infinitely narrow slit, it is a delta function. Then, at the grating, apply focal plane-to-focal-plane theorem, the output field incident on the grating is:  
                                    (Subscript[U, G])[x,y] =  k ((e^(2i k f))/(2 π i f))(F^(-1))[δ][k(x/f),k(y/f)]                   (2.1.1a)
                                    (Subscript[U, G])[x] =  k((e^(2i k f))/(2 π i f)) (1/(Sqrt[2 π]))(∫_(-∞)^∞)δ[s-(Subscript[s, 0])](e^(-i k s x/f))ds          
                                                    = k((e^(2i k f))/(2 π i f)) (1/(Sqrt[2 π]))(e^(-i k (Subscript[s, 0]) x/f))ds                  (2.1.1b)
                                                    = k((e^(2i k f))/(2 π i f)) (1/(Sqrt[2 π]))(e^(-i k Tan[α] x))ds               (2.1.1c)   
This shows a plane wave making a angle α with the axis, Tan[α]≡((Subscript[s, 0])/f)

2.1.2 Finite slit swidth

If the slit is finite size a, then it can be treated as a convolution with a delta function:
                              (Subscript[U, input])[s] = (∫_(-∞)^∞)δ[s-(Subscript[s, 0])](Subscript[S, A])[(Subscript[s, 0])] ds                              (2.1.2a)
The purpose of this is that we can apply convolution, for example let the slit be from -a/2 to a/2, then:
                                       (Subscript[U, 3])[u] ~ (∫_(-a/2)^(a/2))(Subscript[U, 3;δ])[s] ds                                      (2.1.2b)

One caveat is that the theorem works only if the lens is infinite. Since δ function beam spread more than any lens, hence, the result is not valid. However, in practice, realistic beam has finite beam width, and if the entire beam of light is within the lens and becomes zero near the edge of the lens, then this approximation is quite valid.

2.2 Output half

Now that we have the beam above incident on a grating (which is the front focal plane of the second lens), what is the field in its back focal plane? 
We can use the result above again:
                                    (Subscript[U, 3])[u,v] =  k ((e^(2i k f))/(2 π i f))(F^(-1))[T][k(u/f),k(v/f)]      (1.4.6b) (2.2.1)
except now, T is the grating multiplied by the incident field (Subscript[U, G])[x] above.
Let's assume that the transmission function is periodic function of some basic function g[x]:
                                    T[x]=g[x-m d]   for    x∈[m d, (m+1) d]    (2.2.2)
where g[x] is defined over one period {0,d}, and is 0 outside {0, N d}.

http://lightsmyth.com/products/nanopatterned_diffraction_gratings/index.php

Then:
        (F^(-1))[T (Subscript[U, G])][k(u/f)]=(1/(Sqrt[2 π])) k((e^(2i k f))/(2 π i f)) (1/(Sqrt[2 π]))(∑_(m=0)^N)(∫_(m d)^((m+1) d))g[x-m d](e^(-i k x u/f))(e^(-i k (Subscript[s, 0]) x/f)) dx
                                  = k((e^(2i k f))/((((2 π))^2) i f)) (∑_(m=0)^N)(∫_(m d)^((m+1) d))g[x-m d](e^(-i  k x  (((u+(Subscript[s, 0]))/f)))) dx
                                   = k((e^(2i k f))/((((2 π))^2) i f))(∑_(m=0)^N)(∫_0^d)g[t](e^(-i k (t+m d)  (((u+(Subscript[s, 0]))/f)))) dt
                                   = k((e^(2i k f))/((((2 π))^(3/2)) i f))(∑_(m=0)^N)(e^(-i k (t+m d)  (((u+(Subscript[s, 0]))/f))))(1/(Sqrt[2 π]))(∫_0^d)g[t](e^(-i k t (((u+(Subscript[s, 0]))/f)))) dt     (2.2.3)
It's a bit cumbersome to keep (((u+(Subscript[s, 0]))/f)) around. We will call this (OverTilde[u])
                                   =k((e^(2i k f))/((((2 π))^(3/2)) i f))G[k (OverTilde[u])](∑_(m=0)^N)(e^(-i k (t+m d)  (OverTilde[u])))                         (2.2.4a)
 where we define:                       G[k(OverTilde[u])]=(1/(Sqrt[2 π]))(∫_0^d)g[t](e^(-i k t (OverTilde[u]))) dt                  (2.2.4b)
 Then:
                      (F^(-1))[T (Subscript[U, G])][k(u/f)]=k((e^(2i k f))/((((2 π))^(3/2)) i f))G[k (OverTilde[u])](e^(-i (N/2)k d (OverTilde[u])))((Sin[(((N+1))/2) k d (OverTilde[u])])/(Sin[(1/2) k d (OverTilde[u])]))     (2.2.5)
 The field at the back focal plane is:
                              (Subscript[U, 3])[u,v] =  k ((e^(2i k f))/(2 π i f))(F^(-1))[T (Subscript[U, G])][k(u/f),k(v/f)]                                    (2.2.6a)
                              (Subscript[U, 3])[u] =  k ((e^(2i k f))/(2 π i f))k((e^(2i k f))/((((2 π))^(3/2)) i f))G[k (OverTilde[u])](e^(-i (N/2)k d (OverTilde[u])))((Sin[(((N+1))/2) k d (OverTilde[u])])/(Sin[(1/2) k d (OverTilde[u])]))          (2.2.6b)
We don't really care the constant phase and amplitude term k ((e^(2i k f))/(2 π i f))k((e^(2i k f))/((((2 π))^(3/2)) i f)), hence, we write:
                                          (Subscript[U, 3])[u] ~ G[k (OverTilde[u])](e^(-i (N/2)k d (OverTilde[u])))((Sin[(((N+1))/2) k d (OverTilde[u])])/(Sin[(1/2) k d (OverTilde[u])]))                  (2.2.7)
 Its intensity is:
                                 (((|(Subscript[U, 3])[u,v,(Subscript[z, 2])]|))^2) =(((|G[k (OverTilde[u])]|))^2)(((((Sin[(((N+1))/2) k d (OverTilde[u])])/(Sin[(1/2) k d (OverTilde[u])]))))^2)          (2.2.8)
 We see that it is a product of 2 functions, (((|G[k (OverTilde[u])]|))^2) which is the diffraction pattern from a single period (envelope), and (((((Sin[(((N+1))/2) k d (OverTilde[u])])/(Sin[(1/2) k d (OverTilde[u])]))))^2) which is a very sharp function with peaks at:
                                                                (1/2) k d ((u+(Subscript[s, 0]))/f)=m π             (2.2.9a)
                                                                 (π/λ) d ((u+(Subscript[s, 0]))/f)=m π
                                                                  u=f m (λ/d)-(Subscript[s, 0])               (2.2.9b)

What is N in practical term? It is the number of grating teeth seen by the beam. Usually, the beam is smaller than the grating. If the beam is larger, than N is the number of grating teeth of the entire grating.
Example: f=10 cm, m=1, d=1.4 μ, the beam= 1 cm, grating = 5 cm Then:

d=1.4 ; Beam=1 ; GW=5 ;
Ns=Floor[((Min[Beam,GW] (10^4))/d)]

Source: Grating spectrometer basic calculation

Block[{GW,Ns,spectrometer},

Manipulate[
If[RUN== 1,
spectrometer=Show[, ImageSize-> {400,200}] ; 
Column[{
TextCell["Intensity pattern (one side only)"
,Black,20,Italic,FontFamily->"Calibri",Background->Hue[0.15,1,1]]
,
Manipulate[
GW=5 ;
f=10.;  Ns=Floor[((Min[Beam,GW] (10^4))/d)] ; 
Manipulate[
Grid[{{spectrometer
,Show[Plot[(((((Sin[(Ns+1) π d (u+s)/(λ f)])/((Ns+1) Sin[π d (u+s)/(λ f)]))))^2),{u,-0.25,10.}, PlotRange->{0,1.2}
, Filling-> Axis, PlotStyle-> Red, FillingStyle->Orange]
,Graphics[{Blue,PointSize[0.04],Point[{s,1.1}]} ] 
, ImageSize-> {400,200}]}}]
,{λ,0.5,1},{{s,0.},-0.2,0.2}]
, {{d,1.4},0.5,3}, {Beam,0.01,0.5}
,{{f,10.},2,20}]
}]

,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.7]]]
,Button["GRATING CALCULATION",BaseStyle->{24,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]
]

Demo: Grating spectrometer basic calculation

Source code: Grating spectrometer spectral resolution

Block[{GW,Ns,uRange,spectrometer},
Manipulate[
If[RUN== 1,
spectrometer=Show[, ImageSize-> {400,200}] ; 
Manipulate[
GW=5 ;
Ns=Floor[((Min[Beam,GW] (10^4))/d)] ; 
uRange={ 0.499,0.501}*f/d ;
Column[{
TextCell["Intensity pattern (1st order) - λ-resolution"
,Black,20,Italic,FontFamily->"Calibri",Background->Hue[0.15,1,1]]
,Manipulate[
Grid[{{spectrometer
,Show[Plot[{(((Sin[(Ns+1.) π d u/(λ f)])/((Ns+1.) Sin[π d u/(λ f)])))^2,(((((Sin[(Ns+1.) π d u/( (λ+dλ) f)])/((Ns+1.) Sin[π d u/( (λ+dλ)  f)]))))^2)
,(((((Sin[(Ns+1.) π d u/(λ f)])/((Ns+1.) Sin[π d u/(λ f)]))))^2)+(((((Sin[(Ns+1.) π d u/( (λ+dλ) f)])/((Ns+1.) Sin[π d u/( (λ+dλ)  f)]))))^2)}
,{u,uRange[[1]],uRange[[2]]}, PlotStyle-> {{Thick,Red}, {Thick,Hue[0.35,1,0.7]},{Thick,Dashed,Hue[0.8,1,0.7]}},PlotRange->{0,2}]
, ImageSize-> {400,200}]}}]

,{{λ,0.5,"(Subscript[λ, 0])"},0.499,0.501},{dλ,-0.0005,0.0005}]
}]
, {{d,1.4},0.5,3},{ {Beam,1.},0.05,3.}
,{{f,10.},2,20}]

,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.7]]]
,Button["GRATING SPECTROMETER RESOLUTION",BaseStyle->{24,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]
]

Demo: Grating spectrometer spectral resolution

2.3 Finite slit width

If we have a finite slit width, open from -a/2 to a/2, then the output is a convolution:
                      (Subscript[U, 3])[u] ~  (∫_(-a/2)^(a/2))G[k (OverTilde[u])](e^(-i (N/2)k d (OverTilde[u])))((Sin[(((N+1))/2) k d (OverTilde[u])])/(Sin[(1/2) k d (OverTilde[u])]))d(Subscript[s, 0])                               (2.3.1a)
                         (Subscript[U, 3])[u] ~  (∫_(-∞)^∞)G[k ((u+s)/f)](e^(-i (N/2)k d ((u+s)/f)))((Sin[(((N+1))/2) k d ((u+s)/f)])/(Sin[(1/2) k d ((u+s)/f)]))Π[-s]ds           (2.3.1b)
where Π[-s] is a unit box function from -a/2 to a/2.
                         (Subscript[U, 3])[u] ~  (∫_(-∞)^∞)G[k ((u-(-s))/f)](e^(-i (N/2)k d ((u-(-s))/f)))((Sin[(((N+1))/2) k d ((u-(-s))/f)])/(Sin[(1/2) k d ((u-(-s))/f)]))Π[-s]ds    (2.3.2)
Assuming that the slit is symmetric:                         
                         (Subscript[U, 3])[u] ~  (∫_(-∞)^∞)G[k ((u-s)/f)](e^(-i (N/2)k d ((u-s)/f)))((Sin[(((N+1))/2) k d ((u-s)/f)])/(Sin[(1/2) k d ((u-s)/f)]))Π[s]ds                  (2.3.3)
We recall the convolution theorem for Fourier transform.We recall that:
                G[k (OverTilde[u])](e^(-i (N/2)k d (OverTilde[u])))((Sin[(((N+1))/2) k d (OverTilde[u])])/(Sin[(1/2) k d (OverTilde[u])]))=(F^(-1))[T (Subscript[U, G])][k(u/f)]
                                              =(1/(Sqrt[2 π]))∫T[x] (Subscript[U, G])[x] (e^(-i x ((k/f)u))) dx
                                              =(1/(Sqrt[2 π]))(f/k)∫T[(f/k)q] (Subscript[U, G])[(f/k)q] (e^(-i q u)) dq               (2.3.4)
neglecting some amplitude and phase constant. The Fourier transform of Π[s] is ((2 sin(((a q)/2)))/q). Thus, the convolution is the Fourier transform of their product:
                 (f/k)T[(f/k)q] (Subscript[U, G])[(f/k)q]((2 sin(((a q)/2)))/q) =(f/k)g[(f/k)q](e^(-i q (Subscript[s, 0]) ))((2 sin(((a q)/2)))/q)                (2.3.5)
This can be used for numerical approach:
                           (Subscript[U, 3])[u]~(1/(Sqrt[2 π]))∫(f/k)g[(f/k)q](e^(-i q (Subscript[s, 0]) ))((2 sin(((a q)/2)))/q) (e^(-i q u)) dq     
                           (Subscript[U, 3])[u]~(1/(Sqrt[2 π]))∫(f/k)g[(f/k)q]((2 sin(((a q)/2)))/q)(e^(-i q (u+(Subscript[s, 0])) ))  dq                 (2.3.6)
if we know g[x]. 
For an approximation, if we neglect G[k (OverTilde[u])], which is just a broad envelop function (insensitive to wavelength) as the resolution is determined by the term ((Sin[(((N+1))/2) k d (OverTilde[u])])/(Sin[(1/2) k d (OverTilde[u])])), which acts as a delta function, then, the output linewidth and resolution is related to the slit width with simple relation:
                                                       u=f m (λ/d)                        (2.3.7a)
                                                      Δu=f m ((Δλ)/d)                     (2.3.7b)
                                                      Δλ=Δu (d/(m f))                      (2.3.7c)

m=1 ;d=1.4 ; Beam=1 ; GW=5 ;
f=10.;
Manipulate[Plot[ 10^3*s d/(m f*10^4) ,{s,0.,200}],{f,1,100},{d,1.,5},{m,1,2,1}]

We see how the resolution in nm as a function of slit width.

SpecRes[u_,SpectParm_,λ_]:=Module[{d,GW,f,Beam ,a, Ns},
{d,GW,f,Beam ,a}=SpectParm ;
Ns=Floor[((Min[Beam,GW] (10^4))/d)] ;
NIntegrate[((Sin[(Ns+1) π d ((u-s)/(λ f))])/(Sin[π d ((u-s)/(λ f))])),{s,-(a*10^-4)/2,(a*10^-4)/2}]
]

SpectParm={d,GW,f,Beam ,a};
d=1.4; GW=5. ; f=10 ; Beam=0.5; a=100. ;
Plot[{SpecRes[u,{d,GW,f,Beam ,10},0.5]^2
,SpecRes[u,{d,GW,f,Beam ,50},0.5]^2
,SpecRes[u,{d,GW,f,Beam ,100},0.5]^2
,SpecRes[u,{d,GW,f,Beam ,200},0.5]^2}
,{u,3.5,3.6}, PlotRange-> All]

Including both input and output slit, what we have is a convolution of two slits.

SpectParm={d,GW,f,Beam ,a};
d=1.4; GW=5. ; f=10 ; Beam=1; a=100. ;
W1=Table[SpecRes[u,{d,GW,f,Beam ,100},0.5]^2,{u,3.5,3.6,0.0005}]  ;
W2=Table[UnitBox[u/(0.01)],{u,-0.05,0.05,0.0005}]  ;

ListPlot[{W1,15*W2},Joined-> True, PlotRange-> All]
Lineout=ListConvolve[W2,W1,1] ;
ListPlot[Lineout/Length[Lineout],Joined-> True, PlotRange-> All]

2.4 Realistic beam profile and apodization

In practice, we rarely see the diffraction fringes:

for 2 reasons:
1- they tend to be bunching more tightly than resolution: smaller than the slit width to be seen
2- The actual beam shape does not see a sharp edge cut-off. The oscillatory fringes are associated with a sharp edge. The beam profile usually drops off gradually within the grating or the lens/mirror.

Plot[{Exp[-4*x^2],UnitBox[x]},{x,-2,2}]

If we have a gradual Gaussian beam falling on a grating, what does the diffraction pattern lood like?

Plot[{2Exp[-0.01*x^2],(1+Sign[Cos[2 π x] ])/2},{x,-20,20}, Filling-> {2->  Axis}, AspectRatio-> 0.2]

FullSimplify[(∫_(m d)^((m+1/2) d))(E^(-I k x u/f)) \[DifferentialD]x]

-((I f)/(k u)) (E^(-(m+(1/4)) I d k ((  u)/f))) ((E^((1/4)I d k ((  u)/f)))-(E^(-(1/4)I d k ((  u)/f))))

-((I f)/(k u))(E^(-(m+(1/4)) I d k ((  u)/f)))   ExpToTrig[((E^((1/4)I d k ((  u)/f)))-(E^(-(1/4)I d k ((  u)/f))))]

We see that it is now approximately:
                                     (∑_(m=-N/2)^(N/2))(e^(-(((m d))^2)/(2 (σ^2) )))(∫_(m d)^((m+1/2) d))(e^(-i k x u/f)) dx                  (2.4.1a)
or:                               (1/2)(e^(-((i d k u)/(4 f))))((sin(((d k u)/(4 f))))/((k u/(4 f)))))(∑_(m=-N/2)^(N/2))(e^(-(((m d))^2)/(2 (σ^2) )))(e^(-i m d k(( u)/f)))         (2.4.1b)
Without the factor (e^(-(((m d))^2)/(2 (σ^2) ))), the sum (∑_(m=-N/2)^(N/2))(e^(-i m d k(( u)/f))) would give us the function ((Sin[(((N+1))/2) k d (OverTilde[u])])/(Sin[(1/2) k d (OverTilde[u])])) with fringes. The effect of gradual beam is to soften it.

Source code: Gaussian apodization

Block[{d, f},
Manipulate[
If[RUN== 1,d=1.4 ; f=10. ;
Column[{
Manipulate[{Plot[(∑_(m=-3Ceiling[σ/d])^(3 Ceiling[σ/d]))(E^(-(((m d))^2)/(2 (σ^2) )))(E^(-I m (d/λ) 2 π(( u)/f))),{u,0.,4}, PlotRange->{{3,4},All}, PlotPoints-> 200]
,Plot[Log[Abs[(∑_(m=-3Ceiling[σ/d])^(3 Ceiling[σ/d]))(E^(-(((m d))^2)/(2 (σ^2) )))(E^(-I m (d/λ) 2 π(( u)/f)))]],{u,0.,4}, PlotRange->{{3,4},All}, PlotPoints-> 200]}, {λ,0.45,0.55},{σ,20,200}]}]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.7]]]
,Button["GAUSSIAN APODIZATION",BaseStyle->{24,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]
]

Demo: Gaussian apodization

The output is Gaussian-like. If the beam is not Gaussian-like, it is quite similar that the beam shape Fourier transform is transfer into the output slit. Essentially, the function ((Sin[(((N+1))/2) k d (OverTilde[u])])/(Sin[(1/2) k d (OverTilde[u])])) acts like a delta-function with respect to the beam shape.
In general, people can design to have specific output shape for some expected input beam shape, this is known as apodization (like apodizing a filter). Thus, one can have a box-like output, obtained by engineering the grating itself.

2.5 Etendue of a system

Given a system, we cannot arbitrarily send a beam in with any divergence and any aperture.

A system can only accept a finite A×Ω product, defined as the etendue of the system.
A beam also has its etendue: A Gaussian beam has beam width w will have divergence:
                                                                      Tan[θ]~(( λ)/(π w)) 
Its one dimensional etendue is:          ~2 w 2 θ = ((4 λ)/π) 
The 2D etendue:                                   π ((Tan[θ])^2) π (w^2)~π (λ^2) 
This is diffraction-limited etendue, in other words, it cannot be smaller than this because of fundamental diffraction effect. Hence, a very frequent concept is "diffraction-limited" or non diffraction-limited to characterize a beam property.
A system etendue  A×Ω cannot be smaller than the diffraction value, and of course, is often to be designed as large as possible to accommodate applications. As system is designed so that it has a consistent etendue through out for effient light throughput.

Etendue of lighting system

Etendue of fiber

2.6 Advanced topic: spatial frequency - optical transfer function

See section 3.