Finite Difference Time Domain (FDTD) Method

Finite Difference Time Domain (FDTD) Method February 14, 2014

Yee Cell (3-D) • Kane Yee (1966) suggested the following cell for evaluating discretized Maxwell Equations:

Yee Cell (2-D)

Finite Difference Equations • Maxwell’s Equations are continuous. • For all but the simplest geometry, closed form solutions of Maxwell’s Equations are not able to be obtained. • We must convert the continuous equations (with infinitesimal differences) to discrete equations with finite differences

Discretize and Conquer ((Ey(i+1, j)-Ey(i, j)) – (Ex(i, j+1) – Ex(i, j)))*(1/dx) = -(mu_0)*(Hz(i, j)now - Hz(i, j)before)/dt Hz(i,j) = Hz(i,j) – (dt/mu_0)*(Ey(i,j)-Ey(i+1,j) + Ex(i,j+1) - Ex(i,j))*(1/dx) • (Hz(i, j)-Hz(i-1, j))*(1/dx) = (eps_0)*(Ey(i,j)now - Ey(i,j)before) / dt Ey(i,j) = Ey(i,j) - (dt/eps_0)*(Hz(i-1,j)-Hz(i,j))*(1/dx)

Convergence • When we go from continuous to discrete, does the solution become “incorrect”? • If you keep making the spatial step smaller and smaller, eventually the results converge to the correct result • Rule of thumb: make Dx equal to 1/100 of the wavelength of the incident light • The Courant Condition specifies a value for Dt:

Programming • I wrote examples of 1-D and 2-D FDTD in Fortran 90, C and Python. • Maxim uses Fortran 90 for everythingand therefore so should you. • Ask him for the Absoft Fortran 90 editor. If you try to compile and run with this editor it will give you trouble. I suggest using the cluster.

Fortran 90 • Variables must be declared (but not initialized) at the very beginning • Constants must be declared and initialized at the beginning • After variables are declared and constants are declared and initialized, variables must be set • Then you can start writing your code. • End your code with the word “end”.

Code Example implicit none ! !~~~ Fundamental Constants ~~~! ! double precision, parameter :: pi=3.141592653589793D0,c=299792458.0D0 double precision, parameter :: mu0=4.0D-7*pi,eps0=1.0D0/(c*c*mu0) ! !~~~ Grid, Time Steps, Unit Cells ~~~! ! double precision, parameter :: dx=1.0d-9,dy=1.0d-9 integer, parameter :: Nt=5000 integer, parameter :: Nx=601

! !~~~ Source ~~~! ! integer, parameter :: ms=10 !m position of the incident source double precision, parameter :: E0=1.0D0 !incident pulse peak amplitude double precision pulse(Nt) double precision, parameter :: omega=2.0*pi*2.4180D14*5.00 !laser frequency in eV ! !~~~ Time Step ~~~! ! double precision, parameter :: dt=dx/(2.0D0*c) ! !~~~ Dielectric ~~~! ! double precision eps_r(Nx) double precision, parameter :: eps_rel_dielectric=1.0!4.0 double precision, parameter :: eps_rel_vacuum=1.0 integer, parameter :: xnd = 300 !Grid point where dielectric starts

! !~~~ EM field components ~~~! ! double precision Ex(Nx),Hz(Nx-1) ! !~~~ Boundary condition ~~~! ! double precision ex_low_m2 double precision ex_low_m1 double precision ex_high_m2 double precision ex_high_m1 ! !~~~ Other parameters and variables ~~~! ! integer n,i double precision t double precision, parameter :: dt_eps0=dt/eps0,dt_mu0=dt/mu0

! !~~~ Fire off the source ~~~! ! !Sinewave do n=1,Nt t=dt*n pulse(n)=E0*SIN(omega*t) enddo ! !~~~ Structure ~~~! ! do i=1,Nx if((i<xnd)) then eps_r(i)=eps_rel_vacuum else eps_r(i)=eps_rel_dielectric endif enddo

! !~~~ Initialize remaining variables ~~~! ! Ex=0.0 Hz=0.0 ex_low_m2=0.0 ex_low_m1=0.0 ex_high_m2=0.0 ex_high_m1=0.0 do n=1,Nt

!Hz Update do i=1,Nx-1 Hz(i)=Hz(i)-dt_mu0*(Ex(i+1)-Ex(i))/dx enddo !Ex Update do i=2,Nx-1 if(i==ms) then !Excitation Ex(i)=Ex(i)-(1/eps_r(i))*dt_eps0*(Hz(i)-Hz(i-1))/dx+ & pulse(n) else Ex(i)=Ex(i)-(1/eps_r(i))*dt_eps0*(Hz(i)-Hz(i-1))/dx endif enddo

! IncBoundary Condition Ex(1) = ex_low_m2 ex_low_m2 = ex_low_m1 ex_low_m1 = Ex(2) Ex(Nx) = ex_high_m2 ex_high_m2 = ex_high_m1 ex_high_m1 = Ex(Nx-1) enddo !Nt

Import into Matlab and type Plot(Ex): !Record Ex OPEN(file='Ex.dat',unit=32) do i=1,Nx write(32,*) Ex(i) enddo CLOSE(unit=32) end

Results • As one is learning FDTD, it is helpful to plot electric field over the entire grid

Results • Graphs of the T and R coefficients (obtained from the Poynting vector) are measured in experiments

Next Steps • Look at a single-processor 2-D code with a dielectric cylinder (ask me or Maxim for it) • The absorbing boundary conditions (ABC’s) are implemented via convolutional perfectly matched layers (CPML) • It’s there and it works. In general, there’s no need to mess with it.

Later On • The best metals for SPP’s are silver and gold • Metal is represented by the Drude Model • The interaction of the two-level emitters with the classical EM field is handled by the Liouville-von Neumann equation • Longer simulations benefit from message passing interface (MPI); the simulation is split up between several processors • Ask me or Maxim for the MPI slits code • Rather than using point sources, we can define a Total Field / Scattered Field boundary to introduce a plane wave into the system

Finite Difference Time Domain (FDTD) Method