& 0
 \end{array}\right)
 \left(\begin{array}{c} \bX(t) \\ \bY(t)  \end{array} \right)
 \equiv {\mathcal H}
 \left(\begin{array}{c} \bX(t) \\ \bY(t)  \end{array} \right).
\label{eqn:mtxeqn}$$ Writing $\Psi(t)=(\bX(t),\bY(t))^T$, Eq. (\[eqn:mtxeqn\]) becomes $${\frac{\partial}{\partial t}}\Psi(t)={\mathcal H}\Psi(t). \label{TDME}$$ It is easy to show that ${\mathcal H}$ is skew-symmetric, i.e. ${\mathcal H}^T=-{\mathcal H}$, with respect to the inner product $$\langle\Psi|\Psi^\prime\rangle\equiv\int_V\Psi^T\cdot\Psi^\prime\, d\br,
\label{INNER}$$ where $V$ denotes the volume of the enclosing box.

The formal solution of Eq. (\[TDME\]) is given by $$\Psi(t) = e^{t{\mathcal H}}\Psi(0), \label{eqn:formal}$$ where $\Psi(0)$ represents the initial state of the EM fields and the operator $$U(t)= e^{t{\mathcal H}}, \label{Ut}$$ determines their time evolution. By construction $$\|\Psi(t)\|^2 = \langle \Psi(t)|\Psi(t)\rangle =\int_V
\left[\varepsilon{\bE}^2(t) +\mu{\bH}^2(t)
\right] \, d{\br},
\label{normpsi}$$ relating the length of $\Psi(t)$ to the energy density $$w(t)\equiv
\varepsilon{\bE}^2(t) +\mu{\bH}^2(t),
\label{wt}$$ of the EM fields [@BornWolf]. From $U(t)^T=U(-t)=U^{-1}(t)=e^{-t{\mathcal H}}$ it follows that $\langle U(t)\Psi(0)|U(t)\Psi(0)\rangle=\langle\Psi(t)|\Psi(t)\rangle=
\langle\Psi(0)|\Psi(0)\rangle$. Hence the time-evolution operator $U(t)$ is an orthogonal transformation, rotating the vector $\Psi(t)$ without changing its length $\|\Psi\|$. In physical terms this means that the energy density of the EM fields does not change with time, as expected on physical grounds [@BornWolf].

The fact that $U(t)$ is an orthogonal transformation is essential for the development of an unconditionally stable algorithm to solve the Maxwell equations. 