\mbox{$i$ odd and $j$ even}
\\ X_x(i,j,t)=\sqrt{\mu_{i,j}}H_x(i,j,t), & \mbox{$i$ and $j$ odd}
\end{array} \right.\,.$$ Discretization of the differential operators that appear in Eq. (\[eqn:2DTM\]) yield expressions that have the same structure as Eq. (\[H\]), with extra subscripts to account for the second spatial dimension. It follows that on the lattice, $${\frac{\partial}{\partial t}}\Psi(t)={H} \Psi(t)= \mathop{{\sum}'}_{i=1}^{n_x}
\mathop{{\sum}'}_{j=1}^{n_y} \left[ H^{(x)}(i,j)+H^{(y)}(i,j)
\right]\Psi(t), \label{2Dsplit}$$ where $$\begin{aligned}
H^{(x)}(i,j)=&+& \frac{\be^{{\phantom{T}}}_{i,j+1}\be^T_{i+1,j+1}-
\be^{{\phantom{T}}}_{i+1,j+1}\be^T_{i,j+1}}
{\delta\sqrt{{\varepsilon}_{i+1,j+1}\mu_{i,j+1}}} +
\frac{\be^{{\phantom{T}}}_{i+1,j+1}\be^T_{i+2,j+1}-
\be^{{\phantom{T}}}_{i+2,j+1}\be^T_{i+1,j+1}}
{\delta\sqrt{{\varepsilon}_{i+2,j+1}\mu_{i+2,j+1}}}\,,
\label{eqn:2Dsplit_discr_x}
\\
H^{(y)}(i,j)=&-& \frac{\be^{{\phantom{T}}}_{i+1,j}\be^T_{i+1,j+1}-
\be^{{\phantom{T}}}_{i+1,j+1}\be^T_{i+1,j}}
{\delta\sqrt{{\varepsilon}_{i+1,j+1}\mu_{i+1,j}}} -
\frac{\be^{{\phantom{T}}}_{i+1,j+1}\be^T_{i+1,j+2}-
\be^{{\phantom{T}}}_{i+1,j+2}\be^T_{i+1,j+1}}
{\delta\sqrt{{\varepsilon}_{i+1,j+2}\mu_{i+1,j+1}}}\,,
\label{eqn:2Dsplit_discr_y}\end{aligned}$$ and the superscripts $(x)$ and $(y)$ refer to the derivative with respect to $x$ and $y$ respectively. Note that we use the pair $(i,j)$ to label the $n_xn_y$ unit vectors $\be_{i,j}$. In complete analogy with the 1D case we split Eqs.(\[eqn:2Dsplit\_discr\_x\]) and (\[eqn:2Dsplit\_discr\_y\]) into two parts and obtain for the first-order approximation to $U(\tau)$ $$\begin{aligned}
U_1(\tau) & = & e^{\tau{H_{1}^{(x)}}} e^{\tau{H_{2}^{(x)}}}
e^{\tau{H_{1}^{(y)}}} e^{\tau{H_{2}^{(y)}}}\,, \label{U1D2}\end{aligned}$$ where for instance, in formal analogy to Eq. (\[h2\]), we have $$H_2^{(x)}\,=\,\mathop{{\sum}'}_{i=1}^{n_x}
\mathop{{\sum}'}_{j=1}^{n_y}
\frac{\be^{{\phantom{T}}}_{i+1,j+1}\be^T_{i+2,j+1}-
\be^{{\phantom{T}}}_{i+2,j+1}\be^T_{i+1,j+1}}
{\delta\sqrt{{\varepsilon}_{i+2,j+1}\mu_{i+2,j+1}}}\,.$$ It is not difficult to convince oneself that approximation $U_1(\tau)$ and hence also $U_2(\tau)$ and $U_4(\tau)$ do not commute with the (lattice version of) the divergence. Therefore the divergence of the EM fields in 2D is not conserved. Although the 