construct orthogonal approximations to the time-evolution operator of the Maxwell equations. However, this framework does not specify how to implement the algorithm.

Recently, a spectral-domain split-operator technique has been proposed to solve the TDME [@Harsh00]. The split-operator approach is based on one of the many forms of the Lie-Trotter-Suzuki product formulae. The spectral-domain method makes use of Fast Fourier Transforms to compute the matrix exponentials of the displacement operators. The choice made in Ref. [@Harsh00] yields an approximation to the time-evolution operator that is no longer orthogonal and hence unconditional stability is not automatically guarantueed [@COMMENT1]. In contrast, the methodology that we propose yields efficient, explicit, unconditionally stable schemes that operate on the EM fields defined on the real space grid only. This renders the algorithms rather flexible, avoids wrap-around effects [@Harsh00], and naturally allows for the spatial variations in both the permittivity and the permeability. On the other hand, the implementation described in this paper is by no means unique and leaves a lot of room for further improvements.

For EM fields in a homogeneous medium, Zheng et al. showed that there is an alternating-direction-implicit (ADI) time-stepping algorithm that is unconditionally stable [@Zheng1; @Zheng2]. Conceptually this approach is different 