Project 291 - John Benjamin

# Two dimensional simulation of electric field magnitude for dielectric device

In this project the propegation behaviour of electric fields were investigated in a system which has a dielectric and absorbing boundary. The investigated system was an infinite rectangular rod, using such a construct, the dimensionality of the governing equations (Maxwell's equations) can be reduced. The infinite rod is shown in the figure below:

## 1. The procedure and approach to implementation of Maxwell's equation

The approach in the finite difference time domain (FDTD) is determined by the style of materials used. Since perfectly matched layers (PML) is included in the model to absorb and remove energy from the system, the Maxwell equations has to be considered in the frequency domain. PML is strictly a frequency domain property, and has to be moved to the time domain using inverse Laplace transforms. Using inverse Laplace transforms the frequency domain is then moved to the time domain in the following manner:

$$j \omega [\,\,\,] = \frac{\partial}{\partial t}[\,\,\,]$$

and,

$$\frac{1}{j \omega} [\,\,\,] = \int_{0}^{t} [\,\,\,]dt'$$

The PML functions as an absorbing boundary placed at the edge of the problem space to simulate waves leaving the system. This is done by matching electromagnetic impedance with the problem space while gradually increasing the conductivity as the wave propegate deeper into the PML. The interal structure is shown below:

The dielectric has the following properties: Diagonally anisotropic permittivity $$[ \epsilon_{r} ]$$, anisotropic permeability $$[ \mu_{r} ]$$, and anisotropic conductivity $$[ \sigma ]$$. The brackets represent that the quantities are forming a matrix when interacting with the field vectors of Maxwells equation.

### 1.1 Maxwell's equation from 3D to 2D

The Maxwell's curl equations (and constitutive relations) in frequency space with PML matrix $$[ S ]$$ is given by:

$$\nabla \times E(\omega) = -j \omega \mu_{0} [\mu_{r}][S] H(\omega)$$ $$\nabla \times H(\omega) = [\sigma]E(\omega) + j \omega [S] D(\omega)$$ $$D(\omega) = \epsilon_{0} [\epsilon_{r}] E(\omega)$$

where the Electric field is given by E, the magnetic field by H, the displacement field by D, the vacuum permittivity $$\epsilon_{0}$$, and the vacuum permeability $$\mu_{0}$$. Before reducing the dimensions, the electric field is normalized: $$\tilde{E}=\sqrt{ \frac{\epsilon_{0}}{\mu_{0}} } E=\frac{1}{\eta_{0}} E,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \tilde{D}=\sqrt{ \frac{1}{ \epsilon_{0} \mu_{0}} }=c_{0}D$$ Using a Cartesian coordinate system we note that the rod is extended infinitely far in the z-direction. This implies that all the functions are constant when displaced in the z-direction, i.e. $$\partial / \partial z = 0$$. Expanding the Maxwell's equation in terms of its cartesian coordinate system using the vectors x,y,z:

$$C_{x}^{E}=\frac{\partial \tilde{E}_{z}}{\partial y}=-j\omega \big( \frac{\mu_{rx}}{c_{0}} S_{y} H_{x} \big)$$ $$C_{y}^{E}=-\frac{\partial \tilde{E}_{z}}{\partial x}=-j\omega \big( \frac{\mu_{ry}}{c_{0}} S_{y}^{-1} H_{y} \big)$$ $$C_{z}^{H}=\frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial x}=(\eta_{0} \sigma_{z})\tilde{E}_{z} + j\omega \big( \frac{1}{c_{0}} S_{y} \tilde{D}_{z} \big)$$

and,

$$C_{x}^{H}=\frac{\partial H_{z}}{\partial y}=(\eta_{0}\sigma_{x}) \tilde{E}_{x} + j\omega \big( \frac{1}{c_{0}} S_{y} \tilde{D}_{x} \big)$$ $$C_{y}^{H}=-\frac{\partial H_{z}}{\partial x}=(\eta_{0}\sigma_{y}) \tilde{E}_{y} + j\omega \big( \frac{1}{c_{0}} S_{y}^{-1} \tilde{D}_{y} \big)$$ $$C_{z}^{E}=\frac{\partial \tilde{E}_{y}}{\partial x}-\frac{\partial \tilde{E}_{x}}{\partial x}=-j\omega \big( \frac{\mu_{rz}}{c_{0}} S_{y} H_{z} \big)$$

for the PML oriented such that it absorbes radiation in the y-direction:

$$[ S ] =\begin{bmatrix} S_{y} & 0 & 0\\0 & S_{y}^{-1} & 0\\0 &0 &S_{y}\end{bmatrix},\,\,\,\,\,\,\,\, S_{y}=1 + \frac{1}{j \omega} \frac{\sigma'(y)}{\epsilon_{0}},\,\,\,\,\,\,\,\, \sigma'(y)= \frac{\epsilon_{0}}{2 \Delta t} \big( \frac{y'}{L_{y}} \big)^{3}$$

where $$L_{y}$$ is the length of the PML, and $$y'$$ is the current depth within the PML. Outside of the PML, $$[S] =1$$. The choice of the max value of $$[S]$$ is usually the one shown above, but will depend on the exact parameters (source intensity, geometry of the chamber, material distribution, etc...), and is usually tweaked untill it produces a stable absorption. Selecting the wrong value of the max value can make the system numerically unstable. From the expanded curl equations it is noticed system becomes seperated into two independent sub-system of differential equations. To investigate the wave properties one of these modes can then be selected (given that the source is polarized). We select the mode corresponding to $$\{E_{x} (x,y),E_{y}(x,y),H_{z}(x,y)\}$$. Moving the above equations to the time domain using the inverse Laplace transform:

$$C_{z}^{E}|_{t}=-\big( \frac{\mu_{ry}}{c_{0}} \big)\frac{\partial H_{z}}{\partial t}|_{t} - \big( \frac{\mu_{ry}\sigma'}{c_{0}\epsilon_{0}} \big) H_{z}|_{t}$$ $$C_{x}^{H}|_{t}=(\eta_{0}\sigma_{x})\tilde{E}_{x}|_{t} + \big( \frac{1}{c_{0}} \big)\frac{\partial \tilde{D}_{x}}{\partial t}|_{t} + \big( \frac{\sigma'}{c_{0}\epsilon_{0}} \big)\tilde{D}_{x}|_{t}$$ $$C_{y}^{H}|_{t} + \big( \frac{\sigma'}{\epsilon_{0}} \big) \int_{0}^{t} C_{y}^{H} dt = (\eta_{0}\sigma_{y})\tilde{E}_{y}|_{t} + \big( \frac{\eta_{0}\sigma_{y}\sigma'}{\epsilon_0} \big) \int_{0}^{t}\tilde{E}_{y}|_{t} dt + \big( \frac{1}{c_{0}} \big) \frac{\partial \tilde{D}_{y}}{\partial t}|_{t}$$

### 1.2a From a continuous model to a discrete model: Temporal discretization

The electric field and the magnetic field is temporally displaced relative to eachother. Using the electric field at $$t_{n}$$, the magnetic field defined at the next half step in time can then be calculated $$t_{n}+\frac{\Delta t}{2}$$. Using this magnetic field the electric field can be at the time step $$t_{n}+\Delta t =t_{n+1}$$. This is how the temporal displacements work and how the fields evolve with time when in discrete form. Using finite differences the differentials and the integrals can be substituted with differences and sums, and from the differentials the next time-steps can be evaluated. The update equations are then given by:

$$H_{z}|_{t+\frac{\Delta t}{2}}=-h_{z1}C_{z}^{E}|_{t} + h_{z2}H_{z}|_{t-\frac{\Delta t}{2}}$$ $$\tilde{D}_{x}|_{t+\Delta t}= d_{x1}C_{x}^{H}|_{t+\frac{\Delta t}{2}} + d_{x1}C_{x}^{H}|_{t-\frac{\Delta t}{2}} - d_{x2}\tilde{E}_{x}|_{t} + \tilde{D}_{x}|_{t-\Delta t} - d_{x3}\tilde{D}_{x}|_{t}$$ $$\tilde{D}_{y}|_{t+\Delta t}= d_{y1}C_{y}^{H}|_{t+\frac{\Delta t}{2}} + d_{y1}C_{y}^{H}|_{t-\frac{\Delta t}{2}} + d_{y2} \sum_{t_{i}=\frac{\Delta t}{2}}^{t-\frac{\Delta t}{2}}C_{y}^{H}|_{t_{i}} - d_{x3}\tilde{E}_{y}|_{t} -d_{y4} \sum_{t_{i}=0}^{t}\tilde{E}_{y}|_{t_{i}} + \tilde{D}_{y}|_{t-\Delta t}$$

Note that linear interpolation has to be used on the time steps that are not defined. The full derivation of these update equations are given here. The constant parameters shown in the update equation above are given by:

$$h_{z1}=\frac{2c_{0}\epsilon_{0}\Delta t}{2\mu_{ry}\epsilon_{0}+\mu_{ry}\sigma'\Delta t}, \,\,\,\, h_{z2}=\frac{2\mu_{ry}\epsilon_{0}-\mu_{ry}\sigma'\Delta t}{2\mu_{ry}\epsilon_{0}+\mu_{ry}\sigma'\Delta t}$$ $$d_{x1}=c_{0}\Delta t, \,\,\,\,\,\,\,\,\,\,\, d_{x2}=2c_{0}\eta_{0}\sigma_{x}\Delta t, \,\,\,\,\,\,\,\,\,\,\, d_{3}=\frac{2\sigma'\Delta t}{\epsilon_{0}}$$ $$d_{y1}=2c_{0}\Delta t\big(\frac{1}{2} + \frac{\sigma' \Delta t}{4 \epsilon_{0}} \big), \,\,\,\,\,\,\,\,\,\,\,\,\,\, d_{y2}=2\sigma' \eta_{0}(c_{0}\Delta t)^{2}$$ $$d_{y3}=2\sigma_{y} \eta_{0}c_{0}\Delta t, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, d_{y4}=2\sigma_{y}\sigma' (\eta_{0}c_{0}\Delta t)^{2}$$

### 1.2b From a continuous model to a discrete model: Spatially discretization

To discreteze space a Yee grid was used. A Yee grid is a spacial discretization where the electric fields and the magnetic fields are not defined in the same position, and thus spatially displaced from eachother. The motivation for using a Yee grid is that (1) the divergence equations are intrinsically fulfilled, (2) the material boundary conditions are intrinsically fulfilled. The Yee grid for the 2D system is show below:

The curl terms, $$\{C_{z}^{E},C_{x}^{H},C_{y}^{H}\}$$ can then be written as:

$$C_{z}^{E}(i_{x},i_{y})=\frac{\tilde{E}_{y}(i_{x}+1,i_{y})-\tilde{E}_{y}(i_{x},i_{y})}{\Delta x}-\frac{\tilde{E}_{x}(i_{x},i_{y}+1)-\tilde{E}_{x}(i_{x},i_{y})}{\Delta y}$$ $$=c_{1}\tilde{E}_{y}(i_{x}+1,i_{y})-c_{1}\tilde{E}_{y}(i_{x},i_{y})-c_{2}\tilde{E}_{x}(i_{x},i_{y}+1)+c_{2}\tilde{E}_{x}(i_{x},i_{y})$$ $$C_{x}^{H}(i_{x},i_{y})=\frac{H_{z}(i_{x},i_{y})-H_{z}(i_{x},i_{y}-1)}{\Delta y}=c_{2}H_{z}(i_{x},i_{y})-c_{2}H_{z}(i_{x},i_{y}-1)$$ $$C_{y}^{H}(i_{x},i_{y})=-\frac{H_{z}(i_{x},i_{y})-H_{z}(i_{x}-1,i_{y})}{\Delta x}=-c_{1}H_{z}(i_{x},i_{y})+c_{1}H_{z}(i_{x}-1,i_{y})$$

where, $$c_{1}=\frac{1}{\Delta x}$$ and $$c_{2}=\frac{1}{\Delta y}$$.

### 1.3 The spatial and the temporal resolution

The general spatial resolution $$ds$$ is given by the following:

$$ds:=\frac{\lambda_{0}}{10}=dx=dy,$$ where $$\lambda_{0}$$ is the wavelengh produced by the source while in the vacuum. The selection of the value for the divisor is the main component for a high resolution result, however, as seen later also influces some visual errors. The temporal resolution $$dt$$ is then given by the time light uses to propegate (since the propagation is restricted to 2D the fastest propegation is diagonally, which means it has to be scaled by $$(\sqrt{2})^{-1}$$): $$dt=\frac{ds}{\sqrt{2}c_{0}}.$$

But from certain books is also written as,

$$dt=\frac{1}{c_{0}\sqrt{\big(\frac{1}{\Delta x}\big)^{2}+\big(\frac{1}{\Delta y}\big)^{2}}}$$

### 1.4 Summary of the iteration procedure

The iteration process is shown in the figure below, and uses the formulas shown above.

For the full derivation of the theory see the following PDF form: [Theory.pdf]

## 2. The code constructed by C++ and ROOT

The simulation code was constructed using C++ and Cern ROOT. This code is shown below:

### 2.1 The results

Here are some examples with different input parameters showing some of the output that works well, as well as some graphical bugs and other odd behaviors. The roots of these odd behaviors were not found. The simplest case is an isotropic system with no dielectric and no PML with $$\frac{\lambda}{10}$$ defining the grid lengths dx and dy:

As seen above, there is one huge graphical bug present which is causing the weird patterns. This bug seem to be associated with the resolution of the system as the following sorted out the graphical bug: by adjusting the resolution governed by $$ds=\frac{\lambda}{10}$$ into the following, $$ds=\frac{\lambda}{19}$$, causes the pattern to almost completely go away:

To resolve this bug, different definitions of the grid resolution has been invoked and tested as seen above, however all of them produced some graphical issues (to a large or small extent), but most importantly they depend on the other parameters as well such as the dielectric constant, as seen later. However, we can notice further issues: (1) Larger than expected field amplitudes, (2) no divergence around the source. The latter actually arrising from the fact that the fields does not spread propperly from one cell to the next in the x-direction, and will cause issues in the continuity of the wave equation to be broken when we add the dielectric. However, wavebehavior is attained as predicted: (1) reflections from the boundary condition, (2) production of standing waves due these reflections. Since the Yee grid should set the governing divergence equations, the mistake had to be with the Yee Grid. Taking a second look at the derivation and other books on FDTD electrodynamics no mistake was found by the Author. It is still interesting to see how the system performs when adding dielectric elements as this might give some further clues. For the next test we add a dielectric element: $$\epsilon_{r}=2, \mu_{r}=1, \sigma=0, ds=\frac{\lambda}{19}$$:

This as in the first picture, we attain a lot of junk. Just as previously, we change the grid resolution, $$ds=\frac{\lambda}{25}$$ as a temporary solution:

As mentioned earlier, here we see that there is some discontinuity that we should not be attaining from the Maxwell's equation on the boundary of the dielectric, which seems to stem from the inability to propegate along the x-direction. Otherwise we also have no standing waves beneath the dielectric. There thus seems to be two mains bugs, one bug related to the divergence, and one bug related to the divergence and wave behavior in the x-direction that makes waves exhibit odd behavior. Since the PML was also included in the derivation, it would be a shame to not introduce it. Including the PML:

The PML does not seem to absorb the waves as it should be doing, even when manually increasing the strength of the PML to 20 times the original strength it does not abosorb the waves properly. However, the code that simulates free space has issues the field amplitudes, and these large field values might mask the effect of the PML. What is left is two definite bugs that has been characterized, but unsucessfully resolved: The error with the divergence, and the resolution error.

## Regarding the author

Created by John Benjamin in association with the course PHYS291 at the University of Bergen.
Email: jlo016@uib.no