## 2 Formulation
 \tag{9}\]  
We discretize (7a) at time \((n+1/2)\Delta t\) and obtain the updated equation of the \(Q_{pxx}\):  
\[Q_{pxx}^{n+1}=\frac{2+a_{p,xx}\Delta t}{2-a_{p,xx}\Delta t}Q_{pxx}^{n}+\frac{ \varepsilon_{0}\varepsilon_{p,xx}\Delta t}{2-a_{p,xx}\Delta t}\Big{(}E_{x}^{n+ 1}+E_{x}^{n}\Big{)}. \tag{10}\]  
Similar update expressions can be derived for \(Q_{pxy}\) and \(Q_{pxz}\). We also discretize Equation (9) at time \((n+1/2)\Delta t\), and using the update expressions of \(Q_{pxx}\), \(Q_{pxy}\) and \(Q_{pxz}\), we obtain  
\[\alpha_{xx}E_{x}^{n+1}+\alpha_{xy}E_{y}^{n+1}+\alpha_{xz}E_{z}^{n+1}=\gamma_{x }^{n}, \tag{11}\]  
where \(\alpha_{xx}\), \(\alpha_{xy}\) and \(\alpha_{xz}\) are constants defined as  
\[\alpha_{xx} =\frac{\varepsilon_{0}\varepsilon_{\infty,xx}}{\Delta t}+\frac{ \sigma_{xx}^{\varepsilon}}{2}+\frac{2}{\Delta t}\sum_{p=1}^{M_{xx}}\Re\bigg{(} \frac{\varepsilon_{0}\varepsilon_{p,xx}\Delta t}{2-a_{p,xx}\Delta t}\bigg{)} \tag{12a}\] \[\alpha_{xy} =\frac{\varepsilon_{0}\varepsilon_{\infty,xy}}{\Delta t}+\frac{ \sigma_{xy}^{\varepsilon}}{2}+\frac{2}{\Delta t}\sum_{p=1}^{M_{xy}}\Re\bigg{(} \frac{\varepsilon_{0}\varepsilon_{p,xy}\Delta t}{2-a_{p,xy}\Delta t}\bigg{)}\] (12b) \[\alpha_{xz} =\frac{\varepsilon_{0}\varepsilon_{\infty,xz}}{\Delta t}+\frac{ \sigma_{xz}^{\varepsilon}}{2}+\frac{2}{\Delta t}\sum_{p=1}^{M_{xz}}\Re\bigg{(} \frac{\varepsilon_{0}\varepsilon_{p,xz}\Delta t}{2-a_{p,xz}\Delta t}\bigg{)}, \tag{12c}\]  
and \(\gamma_{x}^{n}\) changes at every time step according to  
\[ \begin{split}&\gamma_{x}^{n}=(\
abla\times\mathbf{H})_{x}^{n+1/2} \\ &+\bigg{(}\frac{\varepsilon_{0}\varepsilon_{\infty,xx}}{\Delta t}- \frac{\sigma_{xx}^{\varepsilon}}{2}-\frac{2}{\Delta t}\sum_{p=1}^{M_{xx}}\Re \bigg{(}\frac{\varepsilon_{0}\varepsilon_{p,xx}\Delta t}{2-a_{p,xx}\Delta t} \bigg{)}\bigg{)}E_{x}^{n}-\frac{2}{\Delta t}\sum_{p=1}^{M_{xx}}\Re\bigg{(} \frac{2a_{p,xx}\Delta t}{2-a_{p,xx}\Delta t}Q_{pxx}^{n}\bigg{)}\\ &+\bigg{(}\frac{\varepsilon_{0}\varepsilon_{\infty,xy}}{\Delta t}- \frac{\sigma_{xy}^{\varepsilon}}{2}-\frac{2}{\Delta t}\sum_{p=1}^{M_{xy}}\Re \bigg{(}\frac{\varepsilon_{0}\varepsilon_{p,xy}\Delta t}{2-a_{p,xy}\Delta t} \bigg{)}\bigg{)}E_{y}^{n}-\frac{2}{\Delta t}\sum_{p=1}^{M_{xy}}\Re\bigg{(} \frac{2a_{p,xy}\Delta t}{2-a_{p,xy}\Delta t}Q_{pxy}^{n}\bigg{)}\\ &+\bigg{(}\frac{\varepsilon_{0}\varepsilon_{\infty,xz}}{\Delta t}- \frac{\sigma_{xz}^{\varepsilon}}{2}-\frac{2}{\Delta t}\sum_{p=1}^{M_{xz}}\Re \bigg{(}\frac{\varepsilon_{0}\varepsilon_{p,xz}\Delta t}{2-a_{p,xz}\Delta t} \bigg{)}\bigg{)}E_{z}^{n}-\frac{2}{\Delta t}\sum_{p=1}^{M_{xz}}\Re\bigg{(} \frac{2a_{p,xz}\Delta t}{2-a_{p,xz}\Delta t}Q_{pxx}^{n}\bigg{)}.\end{split}  \tag{13}\]Working in a similar way, two more equations are obtained and a \(3\times 3\) linear system is formed:  
\[ \begin{split}\alpha_{xx}E_{x}^{n+1}+\alpha_{xy}E_{y}^{n+1}+\alpha_{ xz}E_{z}^{n+1}&=\gamma_{x}^{n}\\ \alpha_{xy}E_{x}^{n+1}+\alpha_{yy}E_{y}^{n+1}+\alpha_{yz}E_{z}^{n +1}&=\gamma_{y}^{n}\\ \alpha_{zx}E_{x}^{n+1}+\alpha_{zy}E_{y}^{n+1}+\alpha_{zz}E_{z}^{n +1}&=\gamma_{z}^{n}.\end{split}  \tag{14}\]  
The update equation of \(E_{x}\) is obtained by the solution of Equation (14) using [PERSON]'s rule:  
\[E_{x}^{n+1}=\frac{1}{D}\big{[}\gamma_{x}^{n}(\alpha_{yy}\alpha_{zz}-\alpha_{ yz}\alpha_{zy})-\gamma_{y}^{n}(\alpha_{xy}\alpha_{zz}-\alpha_{zy}\alpha_{xz})+ \gamma_{z}^{n}(\alpha_{xx}\alpha_{yz}-\alpha_{yy}\alpha_{xz})\big{]}, \tag{15}\]  
where \(D\) is the determinant of the matrix with elements \(\alpha_{ij}\) as in the l.h.s. of Equation (14). The update equations of \(E_{y}\) and \(E_{z}\) are obtained in a similar manner. Finally, the field updating procedure in each iteration of the algorithm is summarized in the following steps:  
1. Update \(H_{x}\), \(H_{y}\) and \(H_{z}\);
2. Store the current values of \(E_{x}\), \(E_{y}\) and \(E_{z}\);
3. Update \(E_{x}\), \(E_{y}\) and \(E_{z}\);
4. Update \(Q_{pxx}\), \(Q_{pxy}\) and \(Q_{pxz}\);
5. Update \(Q_{pyx}\), \(Q_{pyy}\) and \(Q_{pyz}\);
6. Update \(Q_{pzx}\), \(Q_{pzy}\) and \(Q_{pzz}\).