Skip to main content

FDTD Solver

Solver Physics​

FDTD (Finite-Difference Time-Domain) is a numerical computational method used to solve Maxwell's equations, primarily for analyzing the propagation and scattering of electromagnetic waves in the time domain. The algorithm is based on finite difference approximations of Maxwell's equations, simulating the propagation of electromagnetic waves by directly solving the time-varying values of the electric and magnetic fields in the time domain. It discretizes continuous space and time, using difference equations to approximate partial differential equations.

In the medium, the electromagnetic processes satisfy Maxwell's curl equations as follows:

∇×H=∂D∂t+J\nabla \times H = \frac{\partial D}{\partial t}+J
∇×E=∂B∂t\nabla \times E = \frac{\partial B}{\partial t}

Combining with the constitutive relations for an isotropic medium:

D=εED=\varepsilon E
B=μHB=\mu H
J=σEJ= \sigma E

In Cartesian coordinates, the scalar Maxwell's equations are expressed as follows:

∂Hx∂t=1μ(∂Ey∂z−∂Ez∂y)\frac{\partial H_x}{\partial t} =\frac{1}{\mu } (\frac{\partial E_y}{\partial z}- \frac{\partial E_z}{\partial y})
∂Hy∂t=1μ(∂Ez∂x−∂Ex∂z)\frac{\partial H_y}{\partial t} =\frac{1}{\mu } (\frac{\partial E_z}{\partial x}- \frac{\partial E_x}{\partial z})
∂Hz∂t=1μ(∂Ex∂y−∂Ey∂x)\frac{\partial H_z}{\partial t} =\frac{1}{\mu } (\frac{\partial E_x}{\partial y}- \frac{\partial E_y}{\partial x})
∂Ex∂t=1ε(∂Hz∂y−∂Hy∂z−σEx)\frac{\partial E_x}{\partial t} =\frac{1}{\varepsilon } (\frac{\partial H_z}{\partial y}- \frac{\partial H_y}{\partial z}-\sigma E_x)
∂Ey∂t=1ε(∂Hx∂z−∂Hz∂x−σEy)\frac{\partial E_y}{\partial t} =\frac{1}{\varepsilon } (\frac{\partial H_x}{\partial z}- \frac{\partial H_z}{\partial x}-\sigma E_y)
∂Ez∂t=1ε(∂Hy∂x−∂Hx∂y−σEz)\frac{\partial E_z}{\partial t} =\frac{1}{\varepsilon } (\frac{\partial H_y}{\partial x}- \frac{\partial H_x}{\partial y}-\sigma E_z)

For the 2D problems, assuming all physical quantities are independent of the z coordinate, the components of the electromagnetic field can be divided into two independent sets: Ex, Ey, Hz form one set, which corresponds to TE (transverse electric) waves; Hx, Hy, Ez form another set, which corresponds to TM (transverse magnetic) waves.

TE：

∂Hz∂y=ε∂Ex∂t+σEx\frac{\partial H_z}{\partial y} =\varepsilon \frac{\partial E_x}{\partial t}+ \sigma E_x
−∂Hz∂x=ε∂Ey∂t+σEy-\frac{\partial H_z}{\partial x} =\varepsilon \frac{\partial E_y}{\partial t}+ \sigma E_y
∂Ey∂x−∂Ex∂y=−μ∂Hz∂t\frac{\partial E_y}{\partial x}- \frac{\partial E_x}{\partial y}=-\mu \frac{\partial H_z}{\partial t}

TM：

∂Ez∂y=−μ∂Hx∂t\frac{\partial E_z}{\partial y} =-\mu \frac{\partial H_x}{\partial t}
∂Ez∂x=μ∂Hy∂t\frac{\partial E_z}{\partial x} =\mu \frac{\partial H_y}{\partial t}
∂Hy∂x−∂Hx∂y=ε∂Ez∂t+σEz\frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y} =\varepsilon \frac{\partial E_z}{\partial t}+\sigma E_z

Next, it is necessary to discretize the FDTD numerical method about both the spatial domain and time domain. Yee sampled the components of the electromagnetic fields alternately in space, so that each electric/magnetic field component is surrounded by four magnetic/electric field components. As illustrated in the figure below , this is known as the Yee grid (or Yee Cell). Applying this to the scalar form of Maxwell's curl equations allows for the numerical solution of electromagnetic problems.

'fde1'

By using the Yee grid sampling method and applying second-order central differencing to the equations above, along with sampling the electric and magnetic fields at half-step time shifts, we can obtain the time-domain update equations for the electromagnetic fields.Take the first equation of the scalar Maxwell's curl equations as an example:

Hxn+12(i,j+12,k+12)=Hxn−12(i,j+12,k+12)H_x^{n+\frac{1}{2}}(i,j+\frac{1}{2},k+\frac{1}{2}) = H_x^{n-\frac{1}{2}}(i,j+\frac{1}{2},k+\frac{1}{2})
+Δtμ(Eyn(i,j+12,k+1)−Eyn(i,j+12,k)Δz−Ezn(i,j+1,k+12)−Ezn(i,j,k+12)Δy)+\frac{\Delta t}{\mu }(\frac{E_y^n(i,j+\frac{1}{2},k+1)-E_y^n(i,j+\frac{1}{2},k)}{\Delta z}-\frac{E_z^n(i,j+1,k+\frac{1}{2})-E_z^n(i,j,k+\frac{1}{2})}{\Delta y})

In the equations, Δx,Δy,Δz\Delta x, \Delta y,\Delta z represent the spatial step in the Cartesian coordinate, and Δt\Delta t is the time step. The index notation (i,j,k)(i,j,k)is the coordinate of field component at the Yee grid points, as shown in Figure above. Here, σ\sigma and ϵ\epsilon represent the electrical conductivity and permittivity respectively. And the superscript n denotes the time step index.

The spatial step and time step must satisfy the Courant stability condition.

Δt≤1c1(Δx)2+1(Δy)2+1(Δz)2 \Delta t \leq \frac{1}{c \sqrt{\frac{1}{(\Delta x)^2} + \frac{1}{(\Delta y)^2} + \frac{1}{(\Delta z)^2}}} \

where Δx,Δy,Δz\Delta x, \Delta y,\Delta z are the smallest grid size of the simulation in the X, Y, and Z direction, c is the speed of light in the material.

Features Description​

Adds or sets FDTD simulation region and boundary conditions.

Notes : When the FDTD solver is added to the project, user cannot add FDE solver and EME solver at the same time.

Settings Description​

1 General​

'fde1'

1) Dimension: Number of dimensions of the simulation region.The options include 2D and 3D. (Default: 3D)

2) Using Optical Path Estimate Time: It is the switch button that the estimation of simulation time based on optical path.

3) Simulation Time: Simulation time indicates the maximum duration of the simulation to be implemented. In reality, the simulation may end earlier when some of the auto-shutoff conditions are satisfied before running till this maximum simulation time. (Default: 1000 fs)

4) Background Material: The combo box allow user to set the background material from drop down menu. Project, Object Defined Dielectric, and Go to Material Library can be operated.

  • Project: Materials that have been imported into the project can be selected and used as background media in simulation. We can also select "Object Defined Dielectric" to customize the background material index in the refractive index.

    • Object Defined Dielectric: The object defined dielectric material, a default setting if user forgets to set background material, is defined for the current object background material setting. Once users chooses this option, they does not need to set any material from the standard, user, or project material database. And the object-defined dielectric will not be loaded into any material database.
      • Refractive Index: specify background index manually in stead of choosing form library or project. (Default :1).

  • Go to Material Library: If this option is selected, user can go to standard material database to set background material according to needs. The selected material will be imported into the project automatically.

2 Geometry​

1) X, Y, Z: The center position of the simulation region.

2) X Min, X Max: X min, X max position.

3) Y Min, Y Max: Y min, Y max position.

4) Z Min, Z Max: Z min, Z max position.

5) X Span, Y Span, Z Span: X, Y, Z span of the simulation region.

3 Mesh Settings​

1) Mesh Type: Algorithms for generating the mesh are available, to be explained as follows:

  • Auto non-uniform (Default): By default, a non-uniform mesh is generated according to the mesh accuracy.

  • Uniform: A uniform mesh is used over the whole simulation region, disregarding any material properties. Note that when adding a Local Mesh with finer grid resolution than the global mesh, the uniform mesh will automatically become non-uniform.

2) Mesh Accuracy: Mesh Accuracy indicates the minimum number of cells per wavelength. (Default: 15).

3) Minimum Mesh Step Settings: Minimum Mesh Step indicates the minimum mesh step in the whole simulation region.(Default: 1e-4 μm)

4) Mesh Refinement: Select an approach to calculate refined mesh properties.

  • Curve Mesh: The effective permittivity can be derived using the contour path formula, which effectively takes into account the shape of the dielectric interface as well as the weighting of the material filling ratio within the mesh.

  • Staircase: It is possible to compute any point within the Yee grid to determine which material it is filled with, and the properties of that single material are used to describe the E-field at that point.As a result, the discrete structure can barely account for the structural variations within a single Yee grid, leading to a Staircase dielectric constant mesh that is fully consistent with the Cartesian grid. Moreover, all layers are effectively moved to the closest E-field position within the Yee grid, meaning that thicknesses cannot be resolved to finer than dx, and thus, cannot be resolved better than dx.

  • Grading: The Grading factor specifies the biggest ratio of the neighboring spatial grids. (Default: 1.2)

4 Boundary Conditions​

1) Apply in All Directions: Turning on this switch ensures that the boundary conditions are consistent in different directions.

2) PML(Perfect Matched Layer): Electromagnetic waves incident on the Perfect Matched Layer (PML) boundary will be completely absorbed, effectively simulating an ideal open (or non-reflecting) boundary. Compared to traditional boundary conditions, PML boundaries occupy a finite volume around the simulation region, thus having only a limited thickness where the light absorption process occurs.

  • PML Layers: For discretization purposes, the PML region is divided into multiple layers.

  • Kappa, Sigma: To control the absorption performance of the PML boundary according to simulation needs, Kappa and Sigma are used. As referenced, Kappa is dimensionless by definition, while Sigma needs to be normalized to be entered as a dimensionless value in the PML settings table. Specifically, both Kappa and Sigma are evaluated through polynomial variations based on their geometric positions within the PML region.

  • PML Polynomial: This specifies the order of the polynomial used for grading Kappa and Sigma.

2) PEC(Perfect Electric Conductor): The PEC boundary condition is introduced to simulate boundaries that behave as perfect electric conductor. Metal boundaries reflect all electromagnetic waves, thus preventing any energy from passing through the simulation volume enclosed by the metal.

3) Symmetry/Anti-symmetry: Symmetric/Anti-Symmetric boundary condition is used to reduce the simulation time with electromagnetic fields that are symmetric with respect to a plane. The choice between Symmetric and Anti-Symmetric conditions is based on the relationship between the source polarization and the symmetry plane. Select the Symmetric option if the normal to the symmetry plane is tangent to the source polarization. Otherwise, choose the Anti-Symmetric option.

5 Advanced Options​

1) Auto Shutoff: Stops the simulation when the energy in the simulation goes below the “Auto Shutoff Min” when the “Use Early Shutoff ” state is on, you can set the value of Auto Shutoff Min (Default: 1e-4) and Down Sample Time (Default: 200)

2) Down Sample Time: Set the multiple dt time steps down sampling, check whether satisfying the auto shutoff min conditions. (Default: 200)

3) Live Slice Field Display Settings

  • Show Field: Real-time filed slice display switch. (Default: On)

  • Select Field Section: To Choose the real-time field display with a 2D cross-section. (Default: 2D Z Normal)

  • Position：The coordinate position of 2D cross-section. (Default: 0 μm)

  • Select Component：To the components of the displaying field (Default: Ex)

  • Down Sample Time:Set the multiple dt time steps down sampling to renew the displaying field.(200 fs by default)