Finite-Difference Time-Domain (FDTD) simulations are one of the most common methods for numerically simulating the wave equation. The canonical example is the motion of a string when you wiggle one end and hold the other fixed, but in this case we are working with light (electromagnetic waves). In particular, the question "If we have a small particle in water, about 100 nm across, what color will it be?" arises with some regularity. It can be surprisngly difficult to find a reference that has your exact geometry, material, and dielectric environment, so this code lets you answer that question by running a simple 2-D simulation right in your browser. While we make no claims that it produces publication-quality results, it should be fairly accurate for a wide range of parameters (but we make no guarantees) and should definitely be enough to get you started on your research again (or pique your curiousity).
The simulation here is transverse-electric and polarized along the vertical direction (perpendicular to the direction the initial pulse travels). That means that the non-zero field components are: \( E_{x}, E_{y}, \) and \( H_{z} \). The timestepping method is a leapfrog alternating-direction implicit (LADI) FDTD, as described in this paper by Prokopidis and Zografopoulos. The grid is terminated using a Mur absorbing boundary condition (ABC) implemented in the LADI FDTD using the formulation of Gan and Tan. Note that while an ABC is not normally sufficient above 1D, it has been demonstrated that in 2D using the LADI methodology, the ABC has little to no additional reflection as compared with a traditional PML. The major advantage of the (L)ADI method is that we can increase the effective speed of light without introducing numberical instabilities and so this simulation uses a Courant Number of 1.5.
The initial pulse is a hard-source Ricker wavelet of the \(E_{y}\) field centered at 600 nm. The grid step is 10 nm (uniform in X and Y). The simulation runs for 1300 time steps (by default, see Adaptive Timing) and "running" discrete Fourier-Transforms (DFTs) are kept to track the reflected and transmitted electric field strengths. Note that while we only run for a thousand timesteps, we assume that this is sufficient time for the high-frequency response (that we care about) to decay completely. Under this assumption, we then 0-pad our total "running time" to be much longer (16,384 steps). That length was picked to be a power of 2 in the event that I needed to run a post-processing FFT (which are frequently most efficient for powers of 2), but for the implementation of a running DFT any number of steps could be used (provided that the fields have sufficient time to propogate). The transmission/reflection results are normalized to a pre-computed run without a structure, which is why you are not allowed to change the total area or number of steps in the simulation (and also why only discrete background indices are allowed).
One of the major advantages of FDTD methods is that they are able to simulate a wide band of frequencies with a single run of the simulation. However, the results from such a simulation are, generally speaking, meaningless without some way to compare them against what the source would have done in the absence of any scattering media in the simulation volume. So, we need to normalize our results to account for the behavior of the source. For transmission spectra, this is straightforward -- one just divides the value of the transmitted flux for a "real" run by that of an "empty" run. For reflected spectra, the situation is substantially more complicated as interference effects can occur. To account for this, one has to subtract off the "empty" fields (not fluxes) from those for a "real" run when computing the flux. For a much more thorough description of the mathematics invloved, see this part of the meep documentation.
Adaptive timing means that the simulation won't complete until the field has sufficiently decayed. Specifically, it is checking for when the absolute value of the field has decayed to one millionth of what the maximum input value was. This has the potential to make some simulations (e.g. mirrors) complete faster and enables the simulation to keep running so that long-lived resonances have a chance to radiate away all of their energy. Note that this breaks the progress bar at the bottom, as there is no longer any way to know ahead of time how many steps will be required (though it does its best to estimate).
This code uses a critical points model for the dielectric function of materials. The general form of the dielectric function is: $$ \epsilon(\omega) = \epsilon_{\infty} - \frac{\sigma}{i \omega} + \sum_{p=1}^{N_{p}} \left[ \frac{c_{p}}{-i \omega - a_{p}} + \frac{c_{p}^{*}}{-i \omega - a_{p}^{*}} \right] , $$ where \(\epsilon_{\infty}\) is the infinite-frequency permittivity of the material (for metals should be \( 1 \)); \(\sigma\) is the generalized static conductivity (in reality it gets divided by \( \epsilon_{0} \), but the fitting works better with reasonable-sized numbers); and \( c_{p} \) and \( a_{p} \) are the residues and poles, respectively, for each critical point in the approximation (with \(N_{p}\) terms). Note that for lossless materials or materials described only by a static conductivity, \( N_{p} = 0 \) and we ignore the critical points entirely.
Physically, each critical point term corresponds to things like interband transitions and allows highly complex dielectric functions to be represented. Each additional critical point term noticeably slows down the simulation (as is evidenced by the difference in speed between Silica and any of the metals).
This formulation is among the most general that exist and can accurately represent materials that are well-defined by the Drude, Lorentz, or Debye models (as well as any other model I've come across).
For an example of how to fit your own material data, there is Mathematica code for generating the fits here. It is important that you fit the real and imaginary parts of the data simultaneously as well as employ the heuristic constraints \( \epsilon_{\infty} \ge 1 \) and \( \sigma \ge 0 \). Finally, note that poor material fits can make the simulation unstable (it makes for very nice images, but very poor science). Generally I've found the above restrictions to produce stable fits, but there is not a mathematical proof of that fact.
If you're particularly interested (or frustrated that I haven't fixed a bug you found), please stop by the GitHub page to see how it all fits together. The simulation itself is coded in C and has been ported to a browser-compatible format using emscripten to convert to WebAssembly.
If you'd like to use this code to embed a simulation on your site, please do! So long as you respect the (fairly permissive) license and link people back to the GitHub, feel free to use the code as you like.
This simulation leans heavily on previous work. In particular, I appreciate the work by John B. Schneider and Doug Neubauer who provided sufficient online resources for an experimental physicist to cobble together a functioning FDTD simulation. If you're interested in learning more about FDTD simulations, I highly recommend reading through John's online book and trying to implement some of it yourself.