1 FDTD method Finite difference time domain EMLAB

1 FDTD method (Finite difference time domain) EMLAB

Contents 2 EMLAB

Approximating the Time Derivative (1 of 3) 3 An intuitive first guess at approximating the time derivatives in Maxwell’s equations is: This is an unstable formulation. EMLAB

Approximating the Time Derivative (2 of 3) 4 We adjust the finite‐difference equations so that each term exists at the same point in time. These equations will get messy if we include interpolations. Is there a simpler approach? EMLAB

Approximating the Time Derivative (3 of 3) 5 We stagger E and H in time so that E exists at integer time steps (0, Δt, 2 Δ t, …) and H exists at half time steps (0, Δ t/2, 3Δ t/2, 4Δ t/2 …). We will handle the spatial derivatives in curl next lecture in a very similar manner. EMLAB

Derivation of the Update Equations 6 • The “update equations” are the equations used inside the main FDTD loop to calculate the field values at the next time step. • They are derived by solving our finite‐difference equations for the fields at the future time values. EMLAB

Anatomy of the FDTD Update Equation Update coefficient Field at the future time step. Field at the previous time step. 7 (To speed simulation, we calculate these before iteration. ) Curl of the “other” field at an intermediate time step EMLAB

The FDTD Algorithm 8 EMLAB

9 Yee Grid Scheme EMLAB

Representing Functions on a Grid Example physical (continuous) 2 D function A grid is constructed by dividing space into discrete cells Function is known only at discrete points 10 Representation of what is actually stored in memory EMLAB

Grid Unit Cell 11 Whole Grid A function value is assigned to a specific point within the grid unit cell. EMLAB

3 D Grids 12 A three‐dimensional grid looks like this: A unit cell from the grid looks like this: EMLAB

Collocated Grid 13 Within the unit cell, we need to place the field components Ex, Ey, Ez, Hx, Hy, and Hz. A straightforward approach would be to locate all of the field components within in a grid cell at the origin of the cell. EMLAB

Yee Grid 14 Instead, we are going to stagger the position of each field component within the grid cells. EMLAB

Reasons to Use the Yee Grid Scheme 1. Divergence‐free 15 3. Elegant arrangement to approximate Maxwell’s curl equations 2. Physical boundary conditions are naturally satisfied EMLAB

Yee Cell for 1 D, 2 D, and 3 D Grids 1 D Yee Grid 2 D Yee Grids 16 3 D Yee Grid EMLAB

Consequences of the Yee Grid 17 • Field components are in physically different locations • Field components may reside in different materials even if they are in the same unit cell • Field components will be out of phase • Recall the field components are also staggered in time. EMLAB

Visualizing Extended Yee Grids 2 X 2 X 2 Grid 18 4 X 4 Grid (E Mode) EMLAB

19 Finite‐Difference Approximation to Maxwell’s Equations EMLAB

Normalize the Magnetic Field 20 We satisfied the divergence equations by adopting the Yee grid scheme. We now only have to deal with the curl equations. The E and H fields are related through the impedance of the material they are in, so they are roughly three orders of magnitude different. This will cause rounding errors in your simulation and it is always good practice to normalize your parameters so they are all the same order of magnitude. Here we choose to normalize the magnetic field. EMLAB

Curl Equations with Normalized Magnetic Field 21 Using the normalized magnetic field, the curl equations become Proof Note: EMLAB

Expand the Curl Equations 22 EMLAB

Assume Only Diagonal Tensors 23 EMLAB

Final Analytical Equations 24 These are the final form of Maxwell’s equations from which we will formulate the FDTD method. Next, we will approximate these equations with finite‐differences in the Yee grid. EMLAB

Finite‐Difference Equation for Hx 25 EMLAB

Finite‐Difference Equation for Hy 26 EMLAB

Finite‐Difference Equation for Hz 27 EMLAB

Finite‐Difference Equation for Ex 28 EMLAB

Finite‐Difference Equation for Ey 29 EMLAB

Finite‐Difference Equation for Ez 30 EMLAB

Summary of Finite‐Difference Equations 31 Each equation is enforced separately for each cell in the grid. This is repeated for each time step until the simulation is finished. These equations get repeated a lot!! EMLAB

32 Governing Equations For One‐Dimensional FDTD EMLAB

Reduction to One Dimension 33 We saw in Lecture 3 that some problems composed of dielectric slabs can be described in just one dimension. In this case, the materials and the fields are uniform in two directions. Derivatives in these uniform directions will be zero. We will define the uniform directions to be the x and y axes. EMLAB

x and y Derivatives are Zero 34 EMLAB

Maxwell’s Equations Decouple Into Two Independent Modes 35 We see that the longitudinal field components Ez and Hz are always zero. We also see that Maxwell’s equations have decoupled into two sets of two equations. EMLAB

Two Remaining Modes are the Same 36 We see that the longitudinal field components Ez and Hz are always zero. We also see that Maxwell’s equations have decoupled into two sets of two equations. While these modes are physical and would propagate independently, the are numerically the same and will exhibit the same electromagnetic behavior. Therefore, it is only necessary to solve one. We will proceed with the Ey/Hx mode. EMLAB

No Longer Need i and j Array Indices 37 Ex/Hy Mode Ey/Hx Mode EMLAB

38 Derivation of the Basic Update Equations EMLAB

Update Equation for Ex 39 Start with the finite‐difference equation which has Ex in the time‐derivative: Solve this for Ex at the future time value. EMLAB

Update Equation for Ey 40 Start with the finite‐difference equation which has Ey in the time‐derivative: Solve this for Ey at the future time value. EMLAB

Update Equation for Hx 41 Start with the finite‐difference equation which has Hx in the time‐derivative: Solve this for Hx at the future time value. EMLAB

Update Equation for Hy 42 Start with the finite‐difference equation which has Hy in the time‐derivative: Solve this for Hy at the future time value. EMLAB

43 Implementation of the Basic Update Equations for the Ey/Hx Mode EMLAB

Efficient Implementation of the Update Equations 44 The update coefficients do not change their value during the simulation. They should be computed only once before the main FDTD loop and not at each iteration inside the loop. The finite‐difference equations in terms of the update coefficients are: Hx, Ey, �εyy, �μxx, m. Hx, and m. Ey are all stored in 1 D arrays of length Nz. c 0, Δt, and Δ z are single scalar numbers, not arrays. EMLAB

The Basic 1 D‐FDTD Algorithm Calculate the update coefficients before the main loop to save computations. 45 Includes: • Basic FDTD engine Excludes: • Dirichlet BC’s • Calculate source parameters • Simple soft source • Perfectly absorbing BC’s • TF/SF source • Fourier transforms • Reflectance/Transmittance • Calculate grid parameters • Incorporate device EMLAB

Equations → MATLAB Code 46 Update Coefficients Update Equations You will need to update the fields at every point in the grid so these equations are placed inside a loop from 1 to Nz. EMLAB

Each Cell Has Its Own Set of Update Equations 47 Each cell has its own update equation and its own update coefficients. They are implemented separately for each cell. All of these equations have the same general form so it is more efficient to implement them using a loop. For a 1 D grid with 10 cells, think of it this way… EMLAB

48 Numerical Boundary Conditions EMLAB

A Problem at the Boundary of the Grid 49 We must implement the update equations for every point in the grid. What do we do at k = Nz ? Ey(k=Nz+1) does not exist. What do we do at k = 1? Hx(k=0) does not exist. EMLAB

Dirichlet Boundary Condition 50 One easy thing to do is assume the fields outside the grid are zero. This is called a Dirichlet boundary condition. We modify the update equations as follows. EMLAB

Equations �MATLAB Code Update Equations 51 DO NOT use ‘if’ statements to implement boundary conditions. EMLAB

Periodic Boundary Condition 52 Some devices are periodic along a particular direction. When this is the case, the field is also periodic. We modify the update equations as follows. EMLAB

53 Grid Resolution EMLAB

Consideration #1: Wavelength 54 The grid resolution must be sufficient to resolve the shortest wavelength. First, you must determine the smallest wavelength: nmax is the largest refractive index found anywhere in the grid. fc is the maximum frequency in your simulation. Second, you resolve the wave with at least 10 cells. Nλ Comments 10~20 Low contrast dielectrics 20~30 High contrast dielectrics 40~60 Most metallic structures 100~200 Plasmonic devices EMLAB

Consideration #2: Mechanical Features 55 The grid resolution must be sufficient to resolve the smallest mechanical features of the device. First, you must determine the smallest feature: Second, you must divide this by 1 to 4. EMLAB

Calculating the Initial Grid Resolution 56 Must resolve the minimum wavelength Must resolve the minimum structural dimension Set the initial grid resolution to the smallest number computed above EMLAB

“Snap” grid to critical dimensions 57 Decide what dimensions along each axis are critical → Typically this is a lattice constant or grating period along x → Typically this is a layer thickness along y Compute how many grid cells comprise dx, y, and round UP Adjust grid resolution to fit this dimension in grid EXACTLY EMLAB

58 Courant Stability Condition EMLAB

Numerical Propagation Through Grid 59 During a single iteration, a disturbance in the electric field at one point can only be felt by the immediately adjacent magnetic fields. It takes at least two time steps before that disturbance is felt by an adjacent electric field. This is simply due to how the update equations are implemented during a single iteration. EMLAB

Physical Propagation Through Space 60 Electromagnetic waves propagate at the speed of light in a vacuum. Inside a material, they propagate at a reduced speed. EMLAB

A Limit on Δt 61 Over a time duration of one time step �t, an electromagnetic disturbance will travel: Numerical distance covered in one time step: Δz Physical distance covered in one time step: c 0Δt/n • Because of the numerical algorithm, it is not possible for a disturbance to travel farther than one unit cell in a single time step. • We need to make sure that a physical wave would not propagate farther than a single unit cell in one time step. This places an upper limit on the time step. n should be set to the smallest refractive index found anywhere in the grid. Usually this is just made to be 1. 0 and dropped from the equation. EMLAB

The Courant Stability Condition 62 Refractive index is greater than or equal to one, so our condition on Δt can be written more simply as: Sort of a worst case. n=1 produces the fastest possible physical wave. For 2 D or 3 D grids, the condition can be generalized as The Courant stability condition EMLAB

Practical Implementation of the Stability Condition 63 The stability condition will be most restrictive along the shortest dimension of the grid unit cell. A good equation to ensure stability and accuracy on any grid is then Note: A factor of 0. 5 was included here as a safety margin. This can be generalized to account for special cases. 1. Your grid is filled with dielectric and travels slower everywhere. 2. Your model includes dispersive materials with refractive index less than one. EMLAB

Time Step for Our 1 D Grid 64 nbc = refractive index at the grid boundaries. This means a wave will travel the distance of one grid cell in exactly two time steps. It is a necessary condition for the perfect boundary condition we will soon implement. This implies that we cannot have different materials at the two boundaries using this boundary condition. EMLAB

65 Perfect 1 D Boundary Condition EMLAB

The Problem 66 The finite‐difference equation here requires knowledge of the electric field outside of the grid. EMLAB

Implementing the Perfect Boundary Condition 67 If… • the field is only travelling outward at the boundaries • the materials at the boundaries are linear, homogeneous, isotropic and non‐dispersive • The refractive index at both boundaries is nbc. • Δt = nbc Δ z/(2 c 0) exactly. then EMLAB

Visualizing the Perfect Boundary Condition 68 EMLAB

Summary of the 1 D Perfect Boundary Condition 69 Conditions • Waves at the boundaries are only travelling outward. • Materials at the boundaries are linear, homogeneous, isotropic and non‐dispersive • The refractive index is the same at both boundaries and is nbc. • Time step is chosen so physical waves travel 1 cell in exactly two time steps • Δ t = nbc Δ z/(2 c 0) Implementation at z‐Low Boundary At the z‐low boundary, we need only modify the E‐field update equation. Implementation at z‐High Boundary At the z‐high boundary, we need only modify the H‐field update equation. EMLAB

70 Sources EMLAB

The Gaussian Pulse 71 In FDTD, we typically use short pulses as sources. This approximates an impulse function so we can excite a single simulation with a broad range of frequencies at the same time. EMLAB

Frequency Content of Gaussian Pulse 72 The Fourier transform of a Gaussian pulse is another Gaussian function. The frequency content of the Gaussian pulse extends from DC up to around B. EMLAB

73 Designing the Pulse Source (1 of 2) Step 1: Determine the bandwidth B of the source to cover all frequencies of interest. Step 2: Compute the pulse width to encompass the needed bandwidth. Step 3: You may need to reduce your time step. Your Gaussian pulse should be resolved by at least 10 to 20 time steps. Typically, you determine a first Δt based on the Courant stability condition, then determine a second Δ t based on resolving the maximum frequency, and finally go with the smallest value of the two Δ t’s. All of this should be satisfied automatically if Δ t = n Δ z/(2 c 0). EMLAB

Visualizing the Arrays 74 EMLAB

Designing the Pulse Source (2 of 2) 75 Step 4: Compute the delay time t 0 The pulse source must start at zero and gradually increase. NO STEP FUNCTIONS!! WRONG!! The step function at the beginning will induce very large field gradients. STILL WRONG!! While better, this source still starts with a step function that will produce large field gradients. CORRECT!! This source “eases” into the Gaussian source. EMLAB

Two Ways to Incorporate a Source 76 Source #1: Simple Hard Source The simple hard source is the easiest to implement. After updating the field across the entire grid, one field component at one point on the grid is overwritten with the source. This approach injects power into the model, but the source point behaves like a perfect electric conductor or perfect magnetic conductor and will scatter waves. OVERWRITE Not usually a practical source. We won’t use it. Source #2: Simple Soft Source The simple soft source is better than the hard source because it is transparent to scattered waves passing through it. After updating the field across the entire grid, the source function is added to one field component at one point on the grid. This approach injects power into the model in both directions. It is great for testing boundary conditions. ADD TO Rarely used. Use this until we learn TF/SF. EMLAB

77 Total Number of Iterations EMLAB

Considerations for Estimating the Total Number of Iterations 78 The total number of iterations depends heavily on the device being modeled and what properties of it are being calculated. Device Considerations 1. Highly resonant devices typically require more iterations. 2. Purely scattering devices typically require very few iterations. 3. More iterations are needed the more times the waves bounce around in the grid. Information Considerations 1. Calculating abrupt spectral shapes requires many iterations. 2. Calculating fine frequency resolution requires many iterations. 3. Calculating the approximate position of resonances often requires fewer iterations. Great for hunting resonances! EMLAB

A Rule of Thumb 79 How long does it take a wave to propagate across the grid (worst case)? Simulation time T must include the entire pulse of duration τ. Simulation time should allow for � 5 bounces. A rule‐of‐thumb for total simulation time is then Note: For highly resonant devices, this will NOT be enough time! Given the time step Δt, the total number of iterations is then This must be an integer quantity. EMLAB

80 Revised FDTD Algorithm EMLAB

Revised FDTD Algorithm 81 EMLAB

Equations → MATLAB Code 82 EMLAB