matrix exponential $U(\tau)$ that can be written as $$\tilde{U}_{\hbox{Yee}}(\tau)=I+\tau H_1+\tau^2 H_2\,,$$ where $H_1$ and $H_2$ are matrices, the structure of which depend on the lattice dimensionality. The presence of the second order contribution can render the algorithm stable under certain conditions. It seems difficult to determine these conditions for arbitrary (skew-symmetric) $H_1$ and (symmetric) $H_2$, i.e. without making use of very specific knowledge about the elements of $H_1$ and $H_2$. For EM fields moving in free space, a Fourier-space stability analysis of the Yee algorithm yields

$$\tau\leq\frac{\Delta}{c\sqrt{d}}\,,$$

as the condition for stability [@Taflove]. Here $c$ is the light velocity in vacuum and $\Delta$ denotes the spatial mesh size [@Taflove].

A systematic approach to construct orthogonal approximations to matrix exponentials, i.e. to construct unconditionally stable algorithms, is to make use of the Lie-Trotter-Suzuki formula [@Trotter59; @Suzuki77] $$e^{t(H_1+\ldots+H_p)}=
\lim_{m\rightarrow\infty}
\left(\prod_{i=1}^p e^{t{H}_i/m}\right)^m,
\label{TROT}$$ and generalizations thereof [@Suzuki8591; @DeRaedt83]. Applied to the case of interest here, the success of this approach relies on the basic but rather trivial premise that the matrix $H$ can be written as $${H}=\sum_{i=1}^{p}{H}_i,
\label{decompos}$$ where each of the matrices ${H}_i$ is real and skew-symmetric.

The expression Eq. (\[TROT\]) suggests that $$U_1(\tau)=e^{\tau{H}_1}\ldots e^{\tau{H}_p},
\label{tsapprox}$$ might be a good approximation to 