LaunchKernels[];

Electromagnetic finite-difference time-domain (FDTD) in Mathematica

Dr. John Kohoutek*
*National Institute of Standards and Technology/University of Maryland
jack.kohoutek@gmail.com

Abstract

Maxwell’s equations in differential form are mapped onto two dimensions and discretized.  Mathematica is used to solve the discretized equations for a fine mesh, with each mesh point being represented by a given epsilon (dielectric constant).  The dielectric constant-dependent geometry is built using the new Region[] functionality of Mathematica 10.  Utilizing this approach, complex, abritrary geometries can be realized and solved for.  Results of various computations will be shown.  Highlights include movies of interesting solved geometries, field intensity profiles, and different boundary conditions, including periodic, perfect metals, and PML.  Far field behavior is analyzed using a given method.

Maxwell’s Equations

Maxwell’s equations are given below in differential form:

∇·E=(ρ/(Subscript[ϵ, 0]))
∇·B=0
∇\[Cross]E=-((∂B)/(∂t))
∇\[Cross]B=(Subscript[μ, o])J+(Subscript[μ, 0])(Subscript[ϵ, 0])((∂E)/(∂t))

In a sourceless, two dimensional region containing a tranverse wave, the above vector equations reduce to:

((∂(Subscript[E, x]))/(∂x))+((∂(Subscript[E, y]))/(∂y)) = 0
((∂(Subscript[B, z]))/(∂z)) = 0
((∂(Subscript[E, y]))/(∂x))-((∂(Subscript[E, x]))/(∂y)) = -((∂(Subscript[B, z]))/(∂t))
((∂(Subscript[B, z]))/(∂y))=(Subscript[μ, 0])(Subscript[ϵ, 0])((∂(Subscript[E, x]))/(∂t))
-((∂(Subscript[B, z]))/(∂x))=(Subscript[μ, 0])(Subscript[ϵ, 0])((∂(Subscript[E, y]))/(∂t))

Discretization

Maxwell’s equations are discretized and solved for in a time-dependent manner, assuming a 3x3 grid per point, in a point-by-point manner.

F={(Subscript[E, x]),(Subscript[E, y]),(Subscript[B, z])}

Defining Geometry

Regions

Regions of different epsilon (dielectric constant) are defined using the new Region[] functionality of Mathematica version 10, defining the mesh and plotting the final structures:

offset={300*^-9,400*^-9};
R1=Disk[offset-{0,0},250*^-9];
R2=Rectangle[offset-{50*^-9,0},offset-{-50*^-9,-2000*^-9}];
R3=Triangle[{offset-{-50*^-9,-2000*^-9},offset-{50*^-9,-2000*^-9},offset-{0,-2100*^-9}}];
{RegionPlot[BooleanRegion[Or,{R1,R2,R3}],AspectRatio->Automatic,ImageSize->Medium],Show[{Graphics[R1],Graphics[R2],Graphics[R3]},ImageSize->Medium]}

Structure Definition

A total structure region is defined using Boolean predicates with a combination of the ImplicitRegion[] and BooleanRegion[] commands.  The region must be inclusive and not just be a “boundary region”.  Therefore, when defining the region, we should use “<”, “>”, etc, instead of “=”, if using algebraic notation, or use “Disc” instead of “Circle”, for instance, if using Geometric primitives.  Structures of different epsilon can be defined as different final regions.

Mesh Generation

A table of corresponding epsilons is made in the same order as a table of the regions defined above.  Next, we generate the mesh by asking each point in the mesh area whether it is a member of each region defined above (using RegionMember[], and if it is, give it the corresponding epsilon.  If it doesn’t belong to any of the regions, assign it epsilon of free space.  
For now, we use a uniform mesh spacing, dx, and dy, but that can change in the future, based on the local refractive index.  The mesh resolution is determined by changing the Δx and Δy parameters.

regionsizex=600*^-9;
regionsizey=2700*^-9;
Δx=Δy=40*^-9;
jtot=Round[regionsizex/Δx];
itot=Round[regionsizey/Δy];
e0=8.85*^-12;
erlist={4,6,7};
regionlist={R1,R2,R3};(*erlist and regionlist must be same length*)
epsmat1=ParallelTable[Piecewise[Append[Transpose[{erlist*e0,RegionMember[#,{j*Δx,i*Δy}]&/@regionlist}],{e0,True}]],{i,itot},{j,jtot}];
ArrayPlot[epsmat1,ColorFunction->"Rainbow"]

Boundary Conditions

Periodic

Reducing Maxwell’s equations to one dimension, discretizing them, setting all physical constants to 1, and using a unit mesh, one is left with the below.  In the below example we start with a unit field pulse at the middle of the 1-D space and allow it to move based simply on the rules of the space.  Notice the wraparound.

Ex=Table[0,{i,30}];
Ex[[15]]=1;
Bz=Table[0,{i,30}];
Bz[[15]]=1;
totaltime=15;
differencekernel={-1,0,1};
result=Table[{Ex,Bz}={Ex,Bz}+{ListConvolve[differencekernel,Bz,1],ListConvolve[differencekernel,Ex,1]},{i,totaltime}];
Manipulate[ListLinePlot[result[[i,1]]],{i,1,totaltime,1}]

PEC = Perfect Electric Conductor

PEC is effectively implemented considering a Drude model with an infinite number of free electrons.  This yields perfect reflection anywhere there is a PEC placed.

ϵ=1-(((Subscript[ω, p])^2)/(ω^2))
(Subscript[ω, p])=(Sqrt[((n (e^2))/(m (Subscript[ϵ, o])))])
n->∞

PML

PML is considered using an auxilliary differential equation (ADE) method, in two dimensions.  This leads to a high quenching of the fields with small back-reflection normal to the PML.

Results of numerical computation

Results

Simulation structure below has a refractive index of n=1.4.  The structure was chosen because the circular region focuses light like a lens.  Regions of higher index tend to “confine” light.  Wavelength of source is 500 nm.

From left to right: 1) Actual mesh of dielectric region used in simulation (derived from region).  Fine mesh is used., 2)  ArrayPlot[] of electromagnetic source,  3) PML grading.

Movie of result.

Above is showing the summed FFT of every point in the simulation region.

Above is showing the spatial map of the FFT at the index of the maximum point in the above plot.

Typical application of nanosphere setup.

The above plot shows the result of an expression for the total energy in the simulation region, showing an intial divergence followed by a smooth convergence as the energy leaves the simulation region due to PML boundaries.

Results of numerical computation cont’d

Because the higher refractive index confines light, and the lens focused it, we can use the straight region as a waveguide after being focused by the lens:

From left to right: 1) Actual mesh of dielectric region used in simulation (derived from region).  Fine mesh is used., 2)  ArrayPlot[] of electromagnetic source,  3) PML grading.

Movie 2 above.

Above is showing the summed FFT of every point in the simulation region.

Above is showing the spatial map of the FFT at the index of the maximum point in the above plot.

The above plot shows the result of an expression for the total energy in the simulation region, showing an intial divergence followed by a smooth convergence as the energy leaves the simulation region due to PML boundaries.

Far field!

In the below I have taken a field pulse with a wavelength in the radio regime and bounced it off of an antenna with a length of 25 cm, corresponding to a common dipole antenna.  I have rotated the dipole antenna and analyzed the reflection coefficient by seeing how much power crosses the red region after bouncing off of the antenna.

The below plot is the reflection coefficient for different angles of the antenna, at different input frequencies - corresponding to the far-field radiation pattern of the dipole antenna.

Future work...

Measuring molecules via EELS LSPR peak shift (gas/molecular sensing)...

Making an electric monopole via simulation!

((∂(Subscript[E, x]))/(∂x))+((∂(Subscript[E, y]))/(∂y)) =(ρ/(Subscript[ϵ, 0]))
((∂(Subscript[B, z]))/(∂z)) = 0
((∂(Subscript[E, y]))/(∂x))-((∂(Subscript[E, x]))/(∂y)) = -((∂(Subscript[B, z]))/(∂t))
((∂(Subscript[B, z]))/(∂y))=(Subscript[μ, 0])(Subscript[ϵ, 0])((∂(Subscript[E, x]))/(∂t))
-((∂(Subscript[B, z]))/(∂x))=(Subscript[μ, 0])(Subscript[ϵ, 0])((∂(Subscript[E, y]))/(∂t))