3.029 Spring 2022
Lecture 24 - 05/02/2022

Diffraction and Imaging

Last lecture, we saw had some fun approximating periodic functions using Fourier Series

More rigorously, we say we can decompose any periodic signal into a sum of (spatial) frequencies, each with their own amplitude and phase

The Fourier Transform is a continuous analog of the same concept, where one takes the period to infinity

Here, we will investigate a very important application of the Fourier Transform in Materials Science

Namely diffraction and microscopy

Fraunhofer diffraction

Positive (Converging) Lens

For reasons which will become apparent soon, we’ll start by investigating optics using a positive (converging) lens

We could code this up ourselves, but there’s this wonderful manipulate on the Wolfram Function Repository

["ShowEquations"->False,Ticks->None]

How exactly does this relate to diffraction and imaging?

Or more importantly to a Fourier Transform?

First, note that rays from a parallel beam, all get converged (hens the lens’ name) to a single point on the focal plane

if the rays forming the plane wave hit the lens “face on”, the focused point lies on the  optical axis (horizontal x-axis above)

if the rays arrive at an angle however, they will get focused on a point displaced from the optical axis, but still at the same focal plane

As such, the lens “transforms” the plane wave’s incident angles to displacement positions on the focal plane

Import["https://upload.wikimedia.org/wikipedia/commons/9/9a/Lens_and_wavefronts_rotated.gif","Animation",Paneled->False,AnimationRunning->False]

Second, recall that a monochromatic light beam can be represented as a sum of many plane waves traveling at different angles

A complicated wave made up of many plane wave components hitting the lense will thus by focused on various different displacement spots on the focal plane

These are the light wave’s “spatial frequencies” and are effectively the components of the 2-D Fourier transform

Aperture Diffraction

More generally, when a beam’s (not necessarily made of light) path is blocked by an object, the wave gets scattered around the object

these scattered wave-fronts interfere constructively/destructively to form bright/dark bands respectively

this is known as a diffraction pattern

Since we’ve seen diffraction patterns can be modeled using Fourier Transforms, let’s investigate diffraction off some common objects, or apertures

To do this, we’ll use the built-in Discrete Fourier Transform (DFT) in Mathematica, Fourier

Fun fact: this formula is also important in science and engineering, a particularly-clever implementation for this algorithm is called the Fast Fourier Transform

with one popular library to compute it called FFTW (Fastest Fourier Transform in the West), and was developed by Prof. Steven G. Johnson at MIT

Fourier however places the largest frequencies at the beginning of a 1D diffraction pattern, or equivalently at the corners of a 2D diffraction pattern

as such we need to center the diffraction pattern by shifting it

fftshift[data_?ArrayQ]:=RotateRight[data,Quotient[Dimensions[data],2]]
ifftshift[data_?ArrayQ]:=RotateRight[data,Ceiling[Dimensions[data]/2]]

Clear[fftshift,ifftshift]

Note the two functions above are equivalent for even dimensions

Finally, note that Fourier returns a complex array

the diffracted waves indeed have both amplitude and phase

However, in typical diffraction/imaging setups we are only sensitive to the intensity of the diffracted wave

which is given by the squared absolute amplitude

diffractionPattern[data_?ArrayQ]:=((Abs[fftshift[Fourier[data]]])^2)

Circular aperture

With[{aperture=DiskMatrix[16,{256,256}]},
Multicolumn[{
ArrayPlot[aperture],
ArrayPlot[((diffractionPattern[aperture])^0.5),ColorFunction->"SunsetColors"]},2]
]

Note we used the same ‘power-law scaling’ trick we used in L13 to increase our dynamical range

Since the diffraction pattern is radially-symmetric, let’s look at a 1D “radial” cut of the diffraction pattern

With[{aperture=DiskMatrix[16,{256,256}]},
ListLinePlot[((diffractionPattern[aperture][[128]])^0.5),FrameTicks->None,Frame->True,LabelStyle->Directive[16,Black],ImageSize->350,PlotStyle->Red,PlotRange->All]
]

Even as we go very far out (where we expect the spatial frequencies to “damp” out), we see noticeable intensity oscillations

This is called Gibbs ringing phenomenon, and arises whenever we try to approximate a discontinuity or sharp edge using Fourier methods

in-fact, we saw this last time with the square wave too

A simple solution to avoid this then, is to use a smoothed-out, or blurred aperture

With[{aperture=GaussianFilter[DiskMatrix[16,{256,256}],5]},
Multicolumn[{
ArrayPlot[aperture],
ArrayPlot[((diffractionPattern[aperture])^0.25),ColorFunction->"SunsetColors"]},2]
]

Transmission Electron Microscope

Electro-magnetic Lenses

Similar to an optical microscope, which uses visible light and optical lenses to magnify objects

A parallel electron beam scatters off a sample, and can be focused and magnified using electro-magnetic lenses

However, since the wavelength of electrons is much smaller than that of visible light, the resolution limit is much better than with an optical microscope

Let’s consider a simple “Gaussian” magnetic lens, i.e. a magnetic field with distribution

where (Subscript[z, 0]) is the center of the distribution and σ is the distribution’s Full-Width-Half-Maximum (FWHM) divided by 2.355

An electron beam passing through this magnetic field will get deflected off its optical axis according to the path-equation

where  is the accelerating voltage,  is the electron charge, and  is the electron’s rest mass

This is a simple ODE we can solve numerically

We’ll consider a magnetic field of strength (Subscript[B, 0])=1.8T, with a FWHM=2mm

and start our electron source a distance 5mm before the center of the lens

pathEquationSolution=With[{σ=0.002/2.355,B0=1.8,z0=0,maxLength=0.01,
e=QuantityMagnitude["ElementaryCharge","SIBase"],
m0=QuantityMagnitude["ElectronMass","SIBase"]
},
ParametricNDSolveValue[{r''[z]==-(e/(8m0 U))(((B0 Exp[-((((z-z0))^2)/(σ^2))]))^2)r[z],r[-0.005]==aperture,r'[-0.005]==initialVelocity},r,{z,-0.005,maxLength},{U,initialVelocity,aperture}]
];

We can then plot our electron paths for various accelerating voltages

Manipulate[
With[{
cols=ResourceFunction["ColorBrewerData"]["Blues",9],
paths=pathEquationSolution[U 1000,#,aperture]&/@Subdivide[-1,1,8]},
Plot[Evaluate[Through[paths[z]]],{z,-0.005,0.01},PlotRange->{-1,1}0.006,PlotLabel->StringTemplate["Accelerating Voltage, U= `1`kV"][U],BaseStyle->14,ImageSize->500,PlotRange->All,Frame->True,LabelStyle->Directive[Black,16],PlotStyle->cols]
],{{U,100,"Accelerating Voltage"},100,300},{{aperture,0.,"Aperture Position"},-0.005,0.005},Paneled->False,SaveDefinitions->True]

Imaging Modes

You’ll notice we also had the capability of starting our electrons off-the optic axis

and we called this parameter “aperture” above

we can understand what this is doing, by looking at different imaging modes

These are typically introduced using “ray-diagrams” like these

In essence, as we saw above, each lens in the TEM column performs a Fourier transform

If 2 (or 4) lenses are used, you get an image

If 3 (or 1) lenses are used, you get the Fourier transform of the image, which we said is called the diffraction pattern

The above is true at the exactly at the focal plane after the last lens

In practice, one uses the strength of the magnetic lenses and a combination of apertures and detectors at different focal planes to switch between imaging and diffraction modes

We will model this simply using multiple gaussian lenses in series

lensPrefactor[acceleratingVoltage_,z_][lensStrength_,lensFWHM_,position_]:=With[{
e=QuantityMagnitude["ElementaryCharge","SIBase"],
m0=QuantityMagnitude["ElectronMass","SIBase"]
},
(e/(8m0 acceleratingVoltage))(((lensStrength Exp[-((((z-position))^2)/(((lensFWHM/2.355))^2))]))^2)
]

multiLensPathEquation[acceleratingVoltage_,lensesDetails_,rayPosition_,rayVelocity_]:=
multiLensPathEquation[acceleratingVoltage,lensesDetails,rayPosition,rayVelocity]=With[{prefactor=Plus@@lensPrefactor[acceleratingVoltage,z]@@@lensesDetails},
NDSolveValue[{r''[z]+prefactor r[z]==0,r[-0.1]==rayPosition,r'[-0.1]==rayVelocity},r,{z,-0.1,0.5}]
]

Finally, we note that at these high speeds electrons are actually relativistic - and thus use a modified accelerating voltage

relativisticElectrons[acceleratingVoltage_]=With[{
m0=QuantityMagnitude["ElectronMass","SIBase"],
c=QuantityMagnitude["SpeedOfLight","SIBase"],
eV=QuantityMagnitude["Electronvolts","SIBase"]
},
(1+((acceleratingVoltage eV )/(m0 (c^2))))acceleratingVoltage];

Putting everything together we have

axisFlip=#/.{x_Line|x_GraphicsComplex:>MapAt[#~Reverse~2&,x,1],x:(PlotRange->_):>x~Reverse~2}&;

rayDiagramTEM[lensStrengths_,lensFWHMs_,flip_:True,lensPositions_:{0.0,0.1,0.2,0.25,0.35},rayPositions_:{0.01,-0.01},acceleratingVoltage_:200 1000]:=
Block[{plot,rayPaths=multiLensPathEquation[relativisticElectrons[acceleratingVoltage],(({lensStrengths,lensFWHMs,lensPositions})^\[Transpose]),#,0]&/@rayPositions},
plot=Plot[Evaluate[Through[rayPaths[z]]],{z,-0.1,0.5},PlotRange->{-1,1}0.025,BaseStyle->14,PlotRange->All,Frame->True,AspectRatio->If[flip,2.4,1/2.4],FrameTicks->None,ImageSize->If[flip,{Automatic,400},{400,Automatic}],Prolog->If[flip,{{{Gray,Opacity[0.25],Disk[{0,-#},{0.0125,0.01}]},
{Gray,Dashed,Line[{{-#,-0.025},{-#,0.025}}]}}&/@lensPositions},{}],FrameStyle->Directive[Black,Thick],PlotStyle->{Red,Blue},ScalingFunctions->{If[flip,"Reverse",Automatic],Automatic}];
If[flip,plot//axisFlip,plot]]

Manipulate[
rayDiagramTEM[{C1,C2,OL,IL,PL},{C1σ,C2σ,OLσ,ILσ,PLσ},True],
{{C1,.15,"First Condenser Lens, C1 Strength [T]"},.05,5},{{C2,1.9725,"Second Condenser Lens, C2 Strength [T]"},.05,5},{{OL,2,"Objective Lens, OL Strength [T]"},.05,5},{{IL,2.0825,"Intermediate Lens, IL Strength [T]"},2.05,2.1},{{PL,2.0015,"Projector Lens, PL Strength [T]"},.05,5},
Delimiter,
{{C1σ,0.025,"First Condenser Lens, C1 FWHM [mm]"},0.002,0.05},{{C2σ,0.01,"Second Condenser Lens, C2 FWHM [mm]"},0.002,0.05},{{OLσ,0.01,"Objective Lens, OL FWHM [mm]"},0.002,0.05},{{ILσ,0.01,"Intermediate Lens, IL FWHM [mm]"},0.002,0.05},{{PLσ,0.01,"Projector Lens, PL FWHM [mm]"},0.002,0.05}
,Paneled->False,SaveDefinitions->True]

Play around with the lens values above to switch between imaging mode and diffraction modes at the focal plane at the bottom of the diagram!

<|"Bright Field Imaging"-><|"Schematic"->rayDiagramTEM[{0.15,1.9725,2,2.0825,2.0015},{2.5,1,1,1,1}/100,False],
"C1"-><|"strength"->Quantity[0.15,"Teslas"],"FWHM"->Quantity[2.5,"Centimeters"]|>,
"C2"-><|"strength"->Quantity[1.9725,"Teslas"],"FWHM"->Quantity[1,"Centimeters"]|>,
"OL"-><|"strength"->Quantity[2,"Teslas"],"FWHM"->Quantity[1,"Centimeters"]|>,
"IL"-><|"strength"->Quantity[2.0825,"Teslas"],"FWHM"->Quantity[1,"Centimeters"]|>,"PL"-><|"strength"->Quantity[2.0015,"Teslas"],"FWHM"->Quantity[1,"Centimeters"]|>|>,
"Selected Area Diffraction"-><|"Schematic"->rayDiagramTEM[{0.15,1.9725,2,2.05925,2.0095},{2.5,1,1,1,1}/100,False],"C1"-><|"strength"->Quantity[0.15,"Teslas"],"FWHM"->Quantity[2.5,"Centimeters"]|>,
"C2"-><|"strength"->Quantity[1.9725,"Teslas"],"FWHM"->Quantity[1,"Centimeters"]|>,
"OL"-><|"strength"->Quantity[2,"Teslas"],"FWHM"->Quantity[1,"Centimeters"]|>,
"IL"-><|"strength"->Quantity[2.05925,"Teslas"],"FWHM"->Quantity[1,"Centimeters"]|>,"PL"-><|"strength"->Quantity[2.0095,"Teslas"],"FWHM"->Quantity[1,"Centimeters"]|>|>|>//Transpose//Dataset//Rasterize

TEM Simulator!

Finally, let’s have some with a simple ‘TEM simulator’

Source Code

TEMSimulator[image_Image,OptionsPattern[]]:=
If[Not[TrueQ[OptionValue["Interactive"]]],
Module[{fft,imageData,originalDiffractionImage,numRows,numCols,filteredFourierData,uv, fourier=Fourier,inverseFourier=InverseFourier,magnificationStatic=1,centerSpotStatic=True,
radiusStatic=If[TrueQ[OptionValue["Interactive"]],5,Last[OptionValue["Interactive"]]["Radius"]],
posStatic=If[TrueQ[OptionValue["Interactive"]],ImageDimensions[image]/2,Last[OptionValue["Interactive"]]["Location"]]},
imageData = ImageData[ColorConvert[RemoveAlphaChannel[image],"Grayscale"],"Real"];
fft = fourier[imageData];
originalDiffractionImage=Image[(((Abs@centerDiffractionData[fft]))^2),ImageSize->400];
{numCols,numRows}= Dimensions[imageData];
filteredFourierData = aperture[magnificationStatic radiusStatic,Dimensions[fft],posStatic-{ Round[numRows/2], Round[numCols/2]},centerSpotStatic]fft;
Return[<|"realSpace"-> image,
"diffractionPattern"->Show[originalDiffractionImage,PlotRange-> {{numCols/2 (1-magnificationStatic),numCols/2 (1+magnificationStatic)},{numRows/2 (1-magnificationStatic),numRows/2 (1+magnificationStatic)}}],
"apertureDiffraction"-> apertureImage[Image[((Abs[centerDiffractionData@filteredFourierData])^2)],
magnificationStatic radiusStatic,{numCols-posStatic[[2]],posStatic[[1]]-numRows}],
"reconstructedImage"->ImageAdjust[Image[Abs[inverseFourier[filteredFourierData]]]]|>]
],
 DynamicModule[
{fft,imageData,originalDiffractionImage,numRows,numCols,filteredFourierData,uv, fourier=Fourier,inverseFourier=InverseFourier},
imageData = ImageData[ColorConvert[RemoveAlphaChannel[image],"Grayscale"],"Real"];
fft = fourier[imageData];
originalDiffractionImage=Image[(((Abs@centerDiffractionData[fft]))^2),ImageSize->400];
{numCols,numRows}= Dimensions[imageData];

Manipulate[
uv=Round[u];
filteredFourierData = aperture[magnification radius,Dimensions[fft],uv-{ Round[numRows/2], Round[numCols/2]},centerSpot]fft;
Deploy[
Show[originalDiffractionImage,PlotRange-> {{numCols/2 (1-magnification),numCols/2 (1+magnification)},{numRows/2 (1-magnification),numRows/2 (1+magnification)}}
]
],
Dynamic[GraphicsRow[
{image,
apertureImage[Image[((Abs[centerDiffractionData@filteredFourierData])^2)],
magnification radius,{numCols-uv[[2]],uv[[1]]-numRows}],ImageAdjust@Image[Abs[inverseFourier[filteredFourierData]]]},
ImageSize-> Large
]
],
Dynamic[{ -Round[numRows/2], -Round[numCols/2]}+uv],
{{u,{numCols,numRows}/2},{1,1},{numRows,numCols},ControlType->Locator,Appearance-> Graphics[{Dashed,Yellow,Thickness[0.005],Circle[{0,0},Scaled[ radius/numCols]]}]},

{{radius,5,"Aperture Radius"},1,Round@Ceiling[Max[{numRows,numCols}]/2],1},
{{centerSpot,False,"Two-Beam-Condition"},{False,True}},
{{magnification,1, "k-space viewport"},0.001,1},
SynchronousUpdating->True,Paneled->False
]
]
]

Options[TEMSimulator] ={"Interactive"->True};

centerDiffractionData[matrix_]:=
 Block[{columns,rows},
{columns,rows}= Dimensions[matrix];
RotateRight[matrix,{Round[columns/2],Round[rows/2]}]]

aperture[radius_,matrixDims:{ncols_,nrows_}, pos:{x_,y_},includeCenterBeam_:False]:=
If[includeCenterBeam,
Clip[RotateRight[DiskMatrix[2,matrixDims],Round[matrixDims/2 ]]+RotateRight[DiskMatrix[radius,matrixDims],Round[matrixDims/2+{-y,x}]]],
RotateRight[DiskMatrix[radius,matrixDims],Round[matrixDims/2+{-y,x}]]]

apertureImage[image_,radius_,u:{x_,y_}]:= 
Block[{imageDimensions,rows,columns},
imageDimensions= {rows,columns}= ImageDimensions[image];
Colorize[ImageResize[ImageTake[image,Round[{x-radius,x + radius}],Round[{y-radius,y+ radius}]],imageDimensions],ColorFunction-> ColorData["StarryNightColors"]]]

Lattice Functions

1D Lattice Functions

angledCosinePlanes[θ_][x_,y_]=1/2+1/2Cos[{x,y}.AngleVector[θ]];
imageData=N[Table[angledCosinePlanes[-10°][x,y],{x,-30,30,0.25},{y,-30,30,0.25}]];
image=Image[imageData];
TEMSimulator[image]

2D Lattice Functions - Hexagonal Lattice

hexagonalLattice[a_,phase_][x_,y_]=1/2+1/2Sum[Cos[((2π)/(a(Sqrt[3])/2)) {Cos[θ-phase],Sin[θ-phase]}.{x,y}],{θ,π/6,π,π/3}];
imageData=N[Table[hexagonalLattice[5,0.][x,y],{x,-30,30,0.25},{y,-30,30,0.25}]];
image=Image[imageData];
TEMSimulator[image]

1D Lattice Functions - Polycrystalline material

TEMSimulator[]

Experimental Images

Experimental Image - Twin Boundary

TEMSimulator[]

Quasicrystals

Quasi Particles - 5-fold symmetry

quasiCrystal[a_][x_,y_]:=1/2+1/2Sum[Cos[((2π)/a){Cos[θ],Sin[θ]}.{x,y}],{θ,π/10,π,π/5}]
imageData=N[Table[quasiCrystal[5][x,y],{x,-30,30,0.25},{y,-30,30,0.25}]];
image=Image[imageData];
TEMSimulator[image]