FDTD simulation algorithm of magnetized plasma media

TL;DR

This paper introduces a new FDTD simulation algorithm for magnetized plasma media, integrating gyrotropic properties into MEEP to model wave propagation in such materials.

Contribution

A novel FDTD simulation method for magnetized plasma media with gyrotropic properties implemented in MEEP, enabling accurate modeling of polarized wave propagation.

Findings

Simulation results match analytical solutions

Demonstrates gyrotropic effects in wave polarization

Validates the new FDTD algorithm

Abstract

A technique of magnetized plasma simulation have been implemented into a well known FDTD solver MEEP as a child class of "susceptibility" class. Magnetized plasma posses gyrotrophic properties, polarization vector is being rotated while the wave propagates through the material. Solid state plasma can be found in semiconductor materials which rises a question of using gyrotrophic properties of plasma in modern photonic development. To test the software, a circular polarized wave propagation through 9mm layer of magnetized plasma has been simulated and compared to the analytical solution.

Equations

$$\frac{B_z^{n+1/2}(i,j+\frac{1}{2},k+\frac{1}{2})-B_z^{n-1/2}(i,j+\frac{1}{2},k+\frac{1}{2})}{\Delta t}$$

$$\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}$$

$$\epsilon (\omega )=\left(\begin{array}{ccc}{\epsilon }_{xx}& j{\epsilon }_{xy}& 0\\ j{\epsilon }_{yx}& {\epsilon }_{yy}& 0\\ 0& 0& {\epsilon }_{zz}\end{array}\right)$$

$${\epsilon }_{xx}={\epsilon }_{yy}=1+{\chi }_{xx}=1+{\chi }_{yy}$$

$${\epsilon }_{zz}=1+{\chi }_{zz}$$

$${\epsilon }_{xy}=-{\epsilon }_{yx}={\chi }_{xy}$$

$${\hat{\chi }}_{xx}(\tau )=\frac{{\omega }_p^2}{{\nu }_c^2+{\omega }_b^2}({\nu }_c+j{\omega }_b)\left[1-e^{-({\nu }_c-j{\omega }_b)\tau }\right]U(\tau )$$

$${\hat{\chi }}_{xy}(\tau )=\frac{{\omega }_p^2}{{\nu }_c^2+{\omega }_b^2}({\omega }_b-j{\nu }_c)\left[1-e^{-({\nu }_c-j{\omega }_b)\tau }\right]U(\tau )$$

$${\hat{\chi }}_{zz}(\tau )=\frac{{\omega }_p^2}{{\nu }_c}\left[1-e^{-{\nu }_c\tau }\right]U(\tau )$$

$$P_i^{(n)}={\vec{P}}_{i,diag}^{(n)}+{\vec{P}}_{i,offdiag}^{(n)}$$

$$\begin{array}{r}{P}_{i,offdiag}=P_{i,offdiag,1}+P_{i,offdiag,2}\end{array}$$

$$P_{i,offdiag,1}^{(n+1)}={\vec{E}}_j^{(n+1)}\frac{{\omega }_p^2}{{\nu }_c^2+{\omega }_b^2}({\omega }_b-j{\nu }_c)\Delta t+P_{i,offdiag,1}^{(n)}$$

$$\begin{array}{r}{P}_{i,offdiag,2}^{(n+1)}=\frac{{\omega }_p^2}{{\nu }_c^2+{\omega }_b^2}\frac{{\omega }_b-j{\nu }_c}{j{\omega }_b-{\nu }_c}\left[E_j^{(n+1)}\left(e^{-({\nu }_c-j{\omega }_b)\Delta t}-1\right)\right]+\\ +P_{i,offdiag,2}^{(n)}\cdot e^{-({\nu }_c-j{\omega }_b)\Delta t}\end{array}$$

$$\begin{array}{r}{P}_{i,diag,2}^{(n+1)}=E_i^{(n+1)}\frac{{\omega }_p^2}{{\nu }_c^2+{\omega }_b^2}\frac{{\nu }_c+j{\omega }_b}{j{\omega }_b-{\nu }_c}\cdot \\ \cdot \left(1-e^{-({\nu }_c-j{\omega }_b)\Delta t}\right)+P_{i,diag,2}\cdot e^{-({\nu }_c-j{\omega }_b)\Delta t}\end{array}$$

$$P_z=P_{z,1}+P_{z,2}$$

$$P_{z,1}^{(n+1)}=E_z^{(n+1)}\frac{{\omega }_p^2\Delta t}{{\nu }_c}+P_{z,1}^{(n)}$$

$$P_{z,2}^{(n+1)}=E_z^{(n+1)}\frac{{\omega }_p^2}{{\nu }_c^2}(e^{-{\nu }_c\Delta t}-1)+P_{z,2}^{(n)}\cdot e^{-{\nu }_c\Delta t}$$

$$P_{diag/offdiag/z,1/2}=Con{v}_{diag/offdiag/z,1/2}$$

$$\vec{D}({\vec{r}}^{\prime })={\epsilon }^{\prime }\vec{E}({\vec{r}}^{\prime })$$

$$T\vec{D}(\vec{r})={\epsilon }^{\prime }T\vec{E}(\vec{r})$$

$${\epsilon }^{\prime }=T\epsilon T^{-1}$$

$$T=I+[v{]}_{\times }+[v{]}_{\times }^2\frac{1-c}{s^2}$$

$$[v{]}_{\times }=\left(\begin{array}{ccc}0& -v_3& v_2\\ v_3& 0& -v_1\\ -v_2& v_1& 0\end{array}\right)$$

$$[v{]}_{\times }=\left(\begin{array}{ccc}0& H_z^0& -H_y^0\\ -H_z^0& 0& H_x^0\\ H_y^0& -H_x^0& 0\end{array}\right)$$

$$\begin{array}{r}c=H_z^0\\ s=\sqrt{(H_x^0{)}^2+(H_y^0{)}^2}\end{array}$$

$$T=\left(\begin{array}{ccc}1& 0& 0\\ 0& -1& 0\\ 0& 0& -1\end{array}\right)$$

$$T=\left(\begin{array}{ccc}-1& 0& 0\\ 0& 1& 0\\ 0& 0& -1\end{array}\right)$$

$$\chi$$

$$+{\chi }_{xy}\left(\begin{array}{ccc}0& 1& 0\\ -1& 0& 0\\ 0& 0& 0\end{array}\right)+{\chi }_{zz}\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right)=$$

$$={\chi }_{xx}{M}_{xx}+{\chi }_{xy}{M}_{xy}+{\chi }_{zz}{M}_{zz}$$

$$\begin{array}{r}P{r}_{diag}=T\cdot M_{xx}\cdot T^{-1}\\ P{r}_{offdiag}=T\cdot M_{xy}\cdot T^{-1}\\ P{r}_z=T\cdot M_{zz}\cdot T^{-1}\end{array}$$

$$P_i(\vec{r})$$

$$+P{r}_{ij,offdiag}\cdot Con{v}_{ij,offdiag}+P{r}_{ij,z}\cdot Con{v}_{ij,z})$$

$$t=\frac{{\tau }_{12}{\tau }_{21}{e}^{-jh\cdot {\hat{n}}_2\cdot \omega /c}}{1+{\rho }_{21}\cdot {\rho }_{12}{e}^{-2jh\cdot {\hat{n}}_2\cdot \omega /c}}$$

Full text

I Introduction

Plasma media are widespread in solid state physics. These media can be magnetized and behave as magnetoactive materials in electromagnetic systems. One of important properties of magnetized plasmas is Parity Time symmetry (PT-symmetry) breaking. That leads to many interesting phenomena such as a new class of Talbot effects [1], unidirectional invisibility [2][3] and absorption [4] control. For semiconductor industry especially important kind of plasmas for semiconductor industry is solid state plasmas. Using magnetized solid state plasma medias would allow to exploit PT symmetry breaking effect in integrated circuits. Apparently we have found out that there was no publicly available FDTD software capable to simulate electromagnetic effects in magnetized plasmas. To overcome this problem we initially modified MEEP to simulate magnetized plasmas with cyclotron resonance caused by Z-axle directed static magnetic field [5]. In this work we are presenting a next version of the algorithm allowing to set any direction of the static magnetic field.

II FDTD basics

Initially FDTD algorithm has been designed by Kane Yee [6] as a technique to solve Maxwell equations in isotropic media. The fields are spatially distributed across a specific grid, so called Yee lattice, shown in Fig. 1. Time distribution is also different for electric and magnetic fields. The electric field stored at $k\cdot\Delta t$ time moments where $k$ is integer, while the magnetic fields keep $(k+\frac{1}{2})\cdot\Delta t$ time positions. This pattern of fields allocation might seems strange, but shows robust stability and reduces required memory [7].

As it has been shown at [6] in isotropic non-dispersive case fields can be calculated as:

[TABLE]

The same way we obtain equations for $\vec{D}$ and $\vec{H}$ accounting also current $\vec{J}$.

When electric and magnetic inductions are known, we have to calculate the intensity fields ($\vec{H}$ and $\vec{E}$) for the next time step. The way it is implemented in MEEP [8] can be assumed as a two step procedure. First we subtract material polarizations from the induction fields ($\vec{D}$ and $\vec{B}$), after we divide this value by $\varepsilon$ or $\mu$. The used approach allows to account material responses of different physical nature at the same time. Each non-trivial material response corresponds to a polarization.

III Recursive convolutions

The simplest case of magnetized plasma is a situation when cyclotron resonance is caused by Z-directed static magnetic field. The corresponding permittivity tensor can be found in [9]:

[TABLE]

where

[TABLE]

As it has been shown in [10] the corresponding electric susceptibilities can be considered as real parts of complex variables:

[TABLE]

Corresponding polarization equations have been derived in the previous work [5]:

[TABLE]

where index $i$ can be $x$ or $y$.To treat polarizations as recursive convolutions and avoid expensive calculations, diagonal and off-diagonal parts are splitted into two parts each:

[TABLE]

Both of the terms of 10 can be calculated as a recursive convolution:

[TABLE]

[TABLE]

Diagonal parts of the polarization calculated similar way and their recursive convolution form can be written as:

[TABLE]

Material response to $E_{z}$ field is z-component of the polarization vector can be accounted a recursive convolution too:

[TABLE]

[TABLE]

[TABLE]

To distinguish physically valuable polarization from the parts, we will rename the parts and will call them ”convolutions”:

[TABLE]

IV Coordinate frame change

Consider two coordinate frames. The first $\vec{r}$ is our old coordinates where bias static magnetic field $\vec{H}^{0}$ directed along Z-axis. For simplicity sake we will consider $|\vec{H}^{0}|=1$. The second one $\vec{r}^{\prime}$ is a new frame where the bias field is arbitrary directed. There is a transformation matrix $T$ changing coordinates of vectors from $\vec{r}$ to $\vec{r}^{\prime}$. As tensor $\varepsilon$ projects vector $\vec{E}(\vec{r})$ into $\vec{D}(\vec{r})$, in the new coordinate frame $\vec{r}^{\prime}$ tensor $\varepsilon^{\prime}$ should project $\vec{E}(\vec{r}^{\prime})$ into $\vec{D}(\vec{r}^{\prime})$:

[TABLE]

If we rewrite $\vec{E}(\vec{r}^{\prime})$ and $\vec{D}(\vec{r}^{\prime})$ using $T$ matrix:

[TABLE]

or

[TABLE]

Now we have to find $T$ matrix, to do this we can use Rodrigues’ rotation formula [11]:

[TABLE]

where $c$ is a cosine between the Z axis and the bias $\vec{H}$ field, s – sinus os the same angle, $v=\vec{H}^{0}\times\vec{e}_{z}$ and $[v]_{\times}$ is defined as:

[TABLE]

or in our case this matrix can be written as:

[TABLE]

while

[TABLE]

The function (21) is not defined in a case when the sinus between $\vec{H}^{0}$ and the Z-axis is 0. To overcome this we manually define two cases. First, when $\vec{H}^{0}$ directed along the Z-axis, is trivial and the $T$ matrix is identity matrix. In the second case, when $\vec{H}^{0}$ is directed opposite the Z axis, described one of two possible matrices:

[TABLE]

[TABLE]

these matrices give the same result. For sake of certainty we will use (25).

V Implementation

As elements of tensor $\chi$ are formulas one can not directly use (20). Instead of that we can decompose $\chi$ as:

[TABLE]

Applying (20) on $M$ matrices from decomposition (27) we obtain projection matrices:

[TABLE]

The projection matrices (28) are generated in the constructor of ”gyroelectric-susceptibility” class.

Since a high precision 3D simulation can require billions of grid cells, memory allocation is an important task. Memory allocation algorithm is determined in ”new_internal_data” function of the class. The function allocates 2 complex double float numbers for each grid cell for every non-zero projection matrix element. Those recursive convolutions are calculated in ”update_P” function of ”gyroelectric-susceptibility” class. New values of the polarizations are calculated as sum of products of convolutions and according projections:

[TABLE]

Function ”subtract_P” and the next steps are not affected by the change of coordinate system.

VI Test case

The analytical model is 9mm layer of magnetized plasma with the bias field $\vec{H}^{0}$ directed along the Z-axis. Parameters for the analytical model are chosen to be: $\omega_{p}=2\pi\cdot 50\cdot 10^{9}$ rad/s, $\omega_{b}=2\pi\cdot 3\cdot 10^{11}$ rad/s, and $\nu_{c}=2\cdot 10^{10}$ Hz. As it has been shown in [12] transmission of a slab can be written as:

[TABLE]

where $\rho$ and $\tau$ are the complex Fresnel reflection and transmission amplitude coefficients for vacuum (index 1) and the plasma (index 2).

There are very strong dependence frequency of transmission $T=|t|^{2}$ for left and right circular polarized waves when we put according $n_{2}$ [9] into (30):

[TABLE]

[TABLE]

The analytical solution is compared to the FDTD calculations result in Fig.2 and Fig.3

To test the coordinate frame change algorithm the bias field $\vec{H}^{0}$ was directed opposite the Z-axis. In this case LCP and RCP transmission swapped as it should be according to the T matrix defined in (25). The reader might substitute it into equations (28) keeping in mind that in this case T matrix is involutory.

VII Summary

A new approach of FDTD simulation of magnetized plasma media allowing to set arbitrary direction of the bias $\vec{H}^{0}$ field is proposed. The implementation is based on MEEP, a widely used FDTD solver. The simulation results for four cases compared to the analytical solutions and matched.

The source code is available on https://github.com/alexfriman/offidiagonal-meep

Acknowledgment

The author would like to thank Nikolay Shubin for the discussion on linear algebra.