Stable lattice Boltzmann model for Maxwell equations in media

The present work shows a method for stable simulations via the lattice Boltzmann (LB) model for electromagnetic waves (EM) transiting homogeneous media. LB models for such media were already presented in the literature, but they suffer from numerical instability when the media transitions are sharp. We use one of these models in the limit of pure vacuum derived from Liu and Yan [Appl. Math. Model. 38, 1710 (2014)] and apply an extension that treats the effects of polarization and magnetization separately. We show simulations of simple examples in which EM waves travel into media to quantify error scaling, stability, accuracy, and time scaling. For conductive media, we use the Strang splitting and check the simulations accuracy at the example of the skin effect. Like pure EM propagation, the error for the static limits, which are constructed with a current density added in a first-order scheme, can be less than $1\%$. The presented method is an easily implemented alternative for the stabilization of simulation for EM waves propagating in spatially complex structured media properties and arbitrary transitions.

Figure 1

A D3Q7 lattice and the directions of the velocity vectors ${\mathit{v}}_i$ accompanying to their $f_i$ within one cell.

Figure 2

Error of total energy $\mathrm{\Delta }W$ at iteration step $N_{\mathrm{it}}$ of the model derived from Ref. [23] (unseparated) and our model (separated) of an EM wave propagating into a dielectric medium. Both simulations were defined via (19) and (20) and $N_x=100$.

Figure 3

Average absolute error $\mathrm{\Delta }W$ at different lattice sizes defined by $N_{\mathrm{it}}$ and its regression curve. The error is obtained by averaging $\mathrm{\Delta }W$ over ${10}^4$ iterations for each lattice size.

Figure 4

$200\times 200$ lattice with dielectric media ${\varepsilon }_r=10$ (gray area) in vacuum. EM fields are generated by the current density $j_z$ represented by the encircled black dot in the middle of the system (100,100).

Figure 5

Comparison of different techniques to stabilize the total energy $W$ of an EM wave generated by an electric current $j_z=sin\left(\pi N_{\mathrm{it}}/50\right)$ radiating in the middle of the system the first 50 iterations. The system is defined with periodic boundary conditions as shown in Fig. 4. The simulation with sharp media transition using the unseparated model is compared to the separated one, the smoothed transition, and the entropic filtering method.

Figure 6

Snapshots of distributed energy density of the electromagnetic wave simulations with unseparated (left column) and separated (right column) model at different iteration steps $N_{\mathrm{it}}$ The fields are generated with a current that is placed in the middle of the $200\times 200$ cells system in which the EM waves cross dielectric media (not visible here; see Fig. 4).

Figure 7

Comparison between simulations on the basis of Eq. (12) (a, left column) with those based on Eq. (10) (b, right column). An one-dimensional Gaussian-shaped EM pulse travels through a 600-cells-wide lattice with periodic boundaries initialized in vacuum ($\varepsilon =1$) and crossing a dielectric (${\varepsilon }^{\prime }=9$, transparent gray area) interface starting at cell number 300 and ending at 600. Simulations at $N_{\mathrm{it}}=0$ and at $N_{\mathrm{it}}=1000$ are almost identical between these two approaches. At $N_{\mathrm{it}}=3000$, one observes in (a) the first results of inaccuracies in the area of cells 200 to 300 (outside dielectric) and inside noised curves around cell 500. This unphysical behavior grows in strength over time as one can see at $N_{\mathrm{it}}=6000$ and clearly at $N_{\mathrm{it}}=7000$ as these inaccuracies tend to grow exponentially with iteration step.

Figure 8

According to case (a) of Fig. 7, the position of the maximum and minimum of $E_z$ in time is shown. Here $E_{z,\min }$ corresponds to the reflected part of the EM wave. $E_{z,\max }$ of $N_{\mathrm{it}}<600$ belongs the initial EM wave and the transmitted one to $N_{\mathrm{it}}\ge 600$.

Figure 9

An EM wave with $\omega =\frac{2\pi }{1000}$ penetrating a conductor. Theoretical prediction of the skin depth $d$, impedance $\left|Z\right|$, speed of propagation $c_{\sigma }$, wavelength $\lambda$, and phase shift $\mathrm{\Delta }\varphi$ depending on the conductivity $\sigma$ are compared with simulation results. The error bars refer to the minimum and maximum that were obtained during the monitoring.

Figure 10

The amplitudes of an EM wave (initialized with $\omega =2\pi /1000, E_y=1, H_z=1$) propagating in a (good) conductor with $\sigma =0.1$ in the $x$ direction is compared with the theoretical prediction. The amplitudes $E_y$ and $H_z$ are amplified by $exp\left(dx\right)$ with $d=0.017176716$, the real part of (27). The numerical noise being amplified is visible at $x=650$ and dominating the further curve.

Figure 11

Comparison of the simulated electrostatic field component $E_x$ of a point charge and the magnetostatic field component $H_y$ of an electrical current with the theoretical prediction.