$U(\tau)$ if $\tau$ is sufficiently small. Most importantly, if all the $H_i$ are real and skew-symmetric, $U_1(\tau)$ is orthogonal by construction. Therefore, by construction, a numerical scheme based on Eq. (\[tsapprox\]) will be unconditionally stable. Using the fact that both $U(\tau)$ and $U_1(\tau)$ are orthogonal matrices, it can be shown that [@DeRaedt87] $$\|U(\tau)-U_1(\tau)\|\leq\frac{\tau^2}{2}\sum_{i<j}\|[{H}_i,{H}_j]\|\,,
\label{tserror}$$ where $[{H}_i,{H}_j]=H_i H_j - H_j H_i$. From Eq. (\[tserror\]) it follows that, in general, the Taylor series of $U(\tau)$ and $U_1(\tau)$ are identical up to first order in $\tau$. We will call $U_1(\tau)$ a first-order approximation to $U(\tau)$.

The product-formula approach provides simple, systematic procedures to improve the accuracy of the approximation to $U(\tau)$ without changing its fundamental symmetries. For example the orthogonal matrix $$U_2(\tau)={U_1(-\tau/2)}^TU_1(\tau/2)=
e^{\tau{H}_p/2}\ldots e^{\tau{H}_1/2}e^{\tau{H}_1/2}\ldots e^{\tau{H}_p/2},
\label{secordapp}$$ is a second-order approximation to $U(\tau)$ [@Suzuki8591; @DeRaedt83]. Suzuki’s fractal decomposition approach [@Suzuki8591] gives a general method to construct higher-order approximations based on $U_1(\tau)$ or $U_2(\tau)$. A particularly useful fourth-order approximation is given by [@Suzuki8591] $$U_4(\tau)=U_2(a\tau)U_2(a\tau)U_2((1-4a)\tau)U_2(a\tau)U_2(a\tau),
\label{fouordapp}$$ where $a=1/(4-4^{1/3})$. The approximations Eqs.(\[tsapprox\]) and (\[secordapp\]), and (\[fouordapp\]) have proven to be very useful in many applications  [@Suzuki77; @Chorin78; @DeRaedt83; @DeRaedt87; @Koboyashi94; @DeRaedt94; @Rouhi95; @Shadwick97; @Krech98; @Tran98; @Michielsen98; @DeRaedt00] and, as we show below, turn out to be equally 