described by ${\mathcal H}$ onto a lattice problem defined by a matrix $H$. Ideally, this mapping should not change the basic symmetries of the original problem. The underlying symmetry of the TDME suggests to use matrices $H$ that are real and skew-symmetric. Formally the time evolution of the EM fields on the lattice is given by $$\Psi(t+\tau) = U(\tau)\Psi(t) = e^{\tau{H}}\Psi(t).
\label{TIMESTEP}$$ The second ingredient of the numerical procedure is to choose an approximation of the time-step operator $U(\tau)$. A standard procedure is to truncate the Taylor series of the matrix exponential [@BELLMAN; @GOLUB] $$U(\tau)=e^{\tau H}=\sum_{n=0}^{\infty}\frac{(\tau{H})^n}{n!}.
\label{taylorU}$$ Retaining terms up to first order in $\tau$ yields $$\tilde{U}(\tau)=I+\tau{H},
\label{approxU0}$$ where $I$ denotes the identity operator. As $\tilde{U}(\tau)\tilde{U}(\tau)^T=I-(\tau{H})^2\not=I$ for $\tau\not=0$, it is clear that the matrix Eq. (\[approxU0\]) is not orthogonal. Making use of the symmetry of $H$ and the positivity of the inner product, we find that $\langle\Psi(\tau)|\Psi(\tau)\rangle=
\langle\tilde{U}(\tau)\Psi(0)|\tilde{U}(\tau)\Psi(0)\rangle=
\langle\Psi(0)|\Psi(0)\rangle+\tau^2\langle H\Psi(0)|H \Psi(0)\rangle >
\langle\Psi(0)|\Psi(0)\rangle$, implying that Eq. (\[UNCOND\]) does not hold. Hence, according to the arguments given above, the (Euler) scheme Eq. (\[approxU0\]) is unstable, a fact which is of course well-known [@Smith85].

The Yee algorithm is based on a leapfrog arrangement [@Yee66] and formally corresponds to an approximation to the 