













\section{Introduction}
\input{S1_Introduction.tex}

\section{Fourth-order Magnus--Zassenhaus splitting}
\label{sec:MZh4}
\input{S2_FirstMethod.tex}




\section{Sixth-order Magnus--Zassenhaus splitting}
\label{sec:CFModd}
\input{S7_TimeSymmetry_HigherOrder.tex}





\section{Numerical examples}
\label{sec:numerics}
\input{S4_Experiments.tex}




\section{Conclusions}
\input{S8_Conclusions.tex}



\section*{Acknowledgments}
The work of Karolina Kropielnicka in this project was financed by The National Center for
Science, based on grant no. 2016/22/M/ST1/00257.

\bibliographystyle{agsm}


\subsection{Intro from Numerische -- kill or keep whichever parts are useful, if ANY}



\subsection{The  approach of this paper}
In this paper we develop an effective computational approach for Schr\"odinger  equations in
semiclassical regime featuring time-dependent potentials that may be highly
oscillatory. We choose optimal tools for this purpose, introducing new ones
where required for further  efficiency. In particular:

\begin{enumerate}
\item We start the derivation of our schemes from the integral-preserving
    simplified-commutator Magnus expansions of \cite{IKS18sinum}.
\item Instead of directly resorting to Lanczos iterations for
    exponentiating the Magnus expansion (which proves to be a costly
    approach in the semiclassical regime), we perform a symmetric
    Zassenhaus splitting on this Magnus expansion. The use of this asymptotic
    splitting algorithm separates scales. Consequently, only very small
    exponents (in terms of spectral radius) need  be exponentiated via
    Lanczos iterations. This enables a far cheaper computation of
    exponentials.
   
\item
    Lastly, we introduce a technique to lower the number of exponentials by
    carefully exploiting the time-symmetry of Magnus expansions. This
    involves starting from a slightly modified version of the Magnus
     expansions of \cite{IKS18sinum}, which will be expressed solely in odd
    powers of the time step. %
\end{enumerate}

\subsection{Magnus--Zassenhaus splittings of order two and four}
The simplest method in the family of methods introduced in this paper is the
splitting
\begin{equation}
 \label{eq:expmidpoint}
 \MM{u}^1 = \exp\!\left(\Frac12 {\mathrm i} \Hc{h} \varepsilon \mathcal{K}_2\right) \exp\!\left(- {\mathrm i} \Hc{h} \varepsilon^{-1} \mathcal{D}_{\Hc{\mu_{0,0}(\dt)}} \right) \exp\!\left(\Frac12 {\mathrm i} \Hc{h} \varepsilon \mathcal{K}_2\right) \MM{u}^0.
\end{equation}
where, under spectral collocation as a choice of spatial discretisation,
$\mathcal{K}_k$ is a discretisation of $\partial_x^k$, $\mathcal{D}_f$ is the diagonal matrix with
values of $f$ along the diagonal, and $\Hc{\mu_{0,0}(\dt)} = \int_0^\Hc{h} V(\xi) d\xi$.
The approximation of $\Hc{\mu_{0,0}(\dt)}$ at the midpoint by $\Hc{h} V(\Hc{h}/2)$, reminiscent of the familiar {\em Strang splitting,\/} results in the
so-called  {\em exponential midpoint rule\/} with $\O{h^3}$ local
error\footnote{In the context of our Magnus expansions, which are {\em
power-truncated}, time-symmetric and consequently feature only odd powers in
the Taylor expansion, any $\O{h^{2k}}$ quadrature method automatically
becomes $\O{h^{2k+1}}$. In particular, evaluation at middle of the interval
incurs $\O{h^3}$ error. See \cite{IKS18sinum} for a more detailed
narrative.}, which is
well known and has been in use for a long while \cite{lubich08fqc}.\\

The first non-trivial splitting in the family of Magnus--Zassenhaus
splittings introduced in this paper is the order four splitting which will be
derived in Section~\ref{sec:MZh4},
\begin{equation}
\label{eq:Zass2}
	 \MM{u}^1 = {\mathrm e}^{\frac12 {\mathrm i}  \Hc{h} \varepsilon \mathcal{K}_2} {\mathrm e}^{-\frac12 {\mathrm i}  \varepsilon^{-1}  \mathcal{D}_{\Hc{\mu_{0,0}(\dt)}}} {\mathrm e}^{\mathcal{W}}  {\mathrm e}^{-\frac12 {\mathrm i}  \varepsilon^{-1}  \mathcal{D}_{\Hc{\mu_{0,0}(\dt)}}} {\mathrm e}^{\frac12 {\mathrm i}  \Hc{h} \varepsilon \mathcal{K}_2} \MM{u}^0,
\end{equation}
where $\mom{1}{1} = \int_0^\Hc{h} (\xi - h/2) V(\xi) d\xi$ and the central
exponent is
\begin{Eqnarray*}
\mathcal{W} & = & \Frac{1}{6}{\mathrm i}  \Hc{h} \varepsilon^{-1} \mathcal{D}_{\left(\partial_x \Hc{\mu_{0,0}(\dt)} \right)^2}
+\left(\mathcal{K}_1 \mathcal{D}_{\partial_x \mom{1}{1}}+ \mathcal{D}_{\partial_x
\mom{1}{1}}\mathcal{K}_1\right)\\
&&-\Frac{1}{3}{\mathrm i} \Hc{h}^2 \varepsilon \left(\mathcal{K}_2 \mathcal{D}_{\partial_x^{2} \Hc{\mu_{0,0}(\dt)} } + \mathcal{D}_{\partial_x^{2}
\Hc{\mu_{0,0}(\dt)} } \mathcal{K}_2\right).
\end{Eqnarray*}



\subsection{Organisation of the paper}
Section~2 is devoted to the derivation of our fourth-order scheme
\R{eq:Zass2}. The derivation commences in Subsection~2.1 with a
straightforward rewriting of the fourth-order integral-preserving
simplified-commutator Magnus expansion of \cite{IKS18sinum} in the
semiclassical regime. An analysis of the size of various terms is carried out
in Subsection~2.2, while keeping discretisation considerations in mind,
followed by a symmetric Zassenhaus splitting in Subsection~2.3.
Implementation details for this scheme are discussed briefly in
Subsection~2.4.

The derivation of the sixth-order Magnus--Zassenhaus scheme of Section~3 is
more involved than the approach in Section~2. In contrast to the fourth-order
Magnus expansion, the sixth and higher order Magnus expansions of
\cite{IKS18sinum} feature even powers of the time step. In Subsection~3.2 we
discuss how these even powers result in a large number of exponential stages.
An approach for remedying this issue, which involves expressing the Magnus
expansion solely in odd powers of the time step, is developed in
Subsections~3.3--3.6. We conclude the derivation by performing a Zassenhaus
splitting on such a modified Magnus expansion in Subsection~3.7.

Numerical examples confirming the efficiency of our schemes are provided in
Section 4, while in the last section we briefly summarise the main
conclusions.

The rules for simplifying commutators from \cite{psalgebra} and the symmetric
Zassenhaus algorithm of \cite{bader14eaf}, which are both utilised in the
derivation of the Zassenhaus splitting in Subsections~2.3 and 3.7, are confined
to Appendices A and B, respectively.




















   
   
   
   
   
   
    %
    %
   
   
   
   
   
   
    %
   
   
   
   
   














\subsection{Simplified-commutator Magnus expansions}
In the semiclassical regime, the order two and order four
simplified-commutator Magnus expansions of \cite{IKS18sinum} that preserve
the integral become
\begin{Eqnarray*}
\label{eq:Omt1} \Theta_0(h) &=& {\mathrm i} \Hc{h} \varepsilon \partial_x^2 - {\mathrm i}  \varepsilon^{-1} \Hc{\mu_{0,0}(\dt)},\\
\label{eq:Omt3} \Theta_2(h) &=& {\mathrm i} \Hc{h} \varepsilon \partial_x^2 - {\mathrm i}  \varepsilon^{-1} \Hc{\mu_{0,0}(\dt)} -
2\Ang{1}{\partial_x \mom{1}{1}},
\end{Eqnarray*}
respectively.  This follows directly from the results of \cite{IKS18sinum}
by accounting for the additional scalar factors  $\varepsilon$ and $\varepsilon^{-1}$. The quantities $\Hc{\mu_{0,0}(\dt)}$ and $\mom{1}{1}$ have been already  introduced
earlier and
\begin{equation}\label{eq:angbra}
  \ang{f}{k} := \Frac12 \left( f \circ \partial_x^k + \partial_x^k \circ f \right), \quad  k \geq 0,\ f \in
  \CC{C}_p^{\infty}(I;\BB{R})
\end{equation}
are the symmetrised differential operators which first appeared in
\cite{bader14eaf}, and have been studied in detail in
\cite{psalgebra,Singh17}. These operators are used here for the same reasons as in \cite{bader14eaf,BIKS,IKS18sinum} -- namely, for
preservation of stability under discretisation after simplification of
commutators.




\subsection{Analysis of size}
As a consequence of \R{eq:dxve}, the symmetrised differential operators scale
as
\begin{equation}\label{eq:angbradisc}
 \ang{f}{k} \leadsto \Frac12 \left( \mathcal{D}_f \mathcal{K}_k + \mathcal{K}_k \mathcal{D}_f \right) = \O{\varepsilon^{-k}}
\end{equation}
upon discretisation, assuming that $f$ is independent of $\varepsilon$. For
convenience, we abuse notation\footnote{This notion can be made more rigorous
by noting that in the semiclassical regime the derivatives of the solution
grow as $\norm{2}{\partial_x^k u}\leq C_k \varepsilon^{-k}$, and consequently
$\norm{2}{\ang{f}{k}(u)} \leq C_k \norm{}{f} \varepsilon^{-k}$. The analysis of the
methods presented here can be made rigorous in accordance with these
observations. See \cite{psalgebra,Singh17} for details. } somewhat and use
the shorthand $\ang{f}{k} = \O{\varepsilon^{-k}}$ instead of \R{eq:angbradisc}.



\begin{theorem}\label{thm:ve}
All nested commutators of $\varepsilon \partial_x^2$, $\varepsilon^{-1} V$ and, in general, terms of
the form $\varepsilon^{k} \ang{f}{k+1}$ are $\O{\varepsilon^{-1}}$ or smaller
\cite{psalgebra,Singh17}.
\end{theorem}
Thus a grade $p$ (i.e. $p-1$ times nested) commutator  of ${\mathrm i} h \varepsilon \partial_x^2$
and $-{\mathrm i} h \varepsilon^{-1} V$, which is typically considered to be $\O{h^p}$,
without accounting for $\varepsilon$ or growth corresponding to finer spatial
discretisation, becomes $\O{h^p \varepsilon^{-1}}$ when analysed taking all factors
into account.

To unify the analysis of magnitude of terms, we express the time step in
terms of the inherent parameter $\varepsilon$,
\begin{equation}\label{eq:hve} h = \O{\varepsilon^{\sigma}},
\quad \sigma>0,
\end{equation}
where a smaller $\sigma$ corresponds to larger time steps. We will be aiming
to develop methods that work for small values of $\sigma$ and at the very
least we may assume that $\sigma \leq 1$.

Taylor expansion for an analytic potential, $V(\xi) = \sum_{k=0}^\infty
[V^{(k)}(0)/k!] \xi^k$, shows that $\mom{1}{1}=\int_0^\Hc{h} (\xi - h/2) V(\xi)
d\xi$ is $\O{h^3}$ for smooth potentials, since its $\O{h^2}$ term,
$V^{(0)}(0)\int_0^\Hc{h} (\xi - h/2) d\xi$, vanishes.

Combining all these ingredients, we conclude that
\begin{equation}
\label{eq:Omt3size}
\Theta_2(h) = \overbrace{{\mathrm i} \Hc{h} \varepsilon \partial_x^2}^{\O{\varepsilon^{\sigma-1}}} - \overbrace{{\mathrm i}  \varepsilon^{-1} \Hc{\mu_{0,0}(\dt)}}^{\O{\varepsilon^{\sigma-1}}} -
\overbrace{2\Ang{1}{\partial_x \mom{1}{1}}}^{\O{\varepsilon^{3\sigma-1}}}.
\end{equation}
Here it becomes evident that the two leading terms in the Magnus expansion
are significantly larger than the trailing terms. As mentioned in the
introduction, the direct usage of Lanczos iterations for $\exp(\Theta_2)$
 is costly due to the $\O{\varepsilon^{\sigma-1}}$ growth of the
leading terms.

On the other hand, the large term ${\mathrm i} \Hc{h} \varepsilon \partial_x^2$ can be exponentiated
via FFTs at a $\O{N \log N}$ cost since it discretises as a circulant matrix
and the diagonal term ${\mathrm i} \varepsilon^{-1} \Hc{\mu_{0,0}(\dt)}$ can be exponentiated directly at the cost of
$\O{N}$ operations. Thus, we resort to Zassenhaus splittings to separate these
vastly disparate scales, whereby only the small $\O{\varepsilon^{3\sigma -1}}$ terms
need  be exponentiated via Lanczos iterations. The small size of the
remaining exponents means we can work with large time steps (a small
$\sigma$, say $\sigma=1/2$ for instance) and still require very few Lanczos
iterations.

\subsection{Symmetric Zassenhaus splitting of order four}

In order to split the exponential of the fourth-order Magnus expansion,
$\Theta_2$, using Zassenhaus splitting algorithm, we follow the procedure
outlined in Appendix~\ref{sec:Zassenhaus}. In particular, the exponential we
wish to approximate here is ${\mathrm e}^{\Theta_2}$. Thus, we take
$\mathcal{W}^{[0]}=\Theta_{2}$, which needs to be split in two components, $X
= W^{[0]}$ and $Y = \mathcal{W}^{[0]}- W^{[0]}$. We usually prefer to extract
the largest terms first yet we may choose between $W^{[0]} = {\mathrm i} \Hc{h} \varepsilon
\partial_x^2$ and $W^{[0]} = -{\mathrm i} \varepsilon^{-1} \Hc{\mu_{0,0}(\dt)}$, which are both
$\O{h\varepsilon^{-1}}=\O{\varepsilon^{\sigma-1}}$.


In the specific method presented here, we opt for the former,
whereby $X = W^{[0]}
= {\mathrm i} \Hc{h} \varepsilon \partial_x^2={\mathrm i} \Hc{h} \varepsilon \ang{1}{2}$ and $Y = \mathcal{W}^{[0]}-
W^{[0]} =- {\mathrm i} \varepsilon^{-1} \ang{\Hc{\mu_{0,0}(\dt)}}{0} - 2 \Ang{1}{\partial_x \mom{1}{1}}$. (Observe that $\ang{1}{k} = \partial_x^k$ and $\ang{f}{0} = f$. We
will switch between these  equivalent notations at our convenience.) Using
the sBCH formula \R{eq:sBCH}, we write
\begin{equation}\label{eq:3:sbch_inverted}
  {\mathrm e}^{\Theta_{2}} = {\mathrm e}^{\mathcal{W}^{[0]}} ={\mathrm e}^{\frac12 {\mathrm i} \Hc{h} \varepsilon \partial_x^2} {\mathrm e}^{\mathcal{W}^{[1]} } {\mathrm e}^{\frac12  {\mathrm i} \varepsilon \Hc{h} \partial_x^2} + \O{h^5 \varepsilon^{-1}},
\end{equation}
where, using the order four expression for $\mathcal{W}^{[1]} =
\CC{sBCH}(-X,X+Y)$ from \R{eq:sBCHshifted},
\begin{Eqnarray}
\label{eq:W1solve}\mathcal{W}^{[1]} &=& - {\mathrm i} \varepsilon^{-1} \ang{\Hc{\mu_{0,0}(\dt)}}{0} - 2
\Ang{1}{\partial_x \mom{1}{1}} +  \Frac{1}{24}  {\mathrm i} \Hc{h}^2 \varepsilon  [[ \ang{\Hc{\mu_{0,0}(\dt)}}{0} ,\ang{1}{2}], \ang{1}{2}] \\
\nonumber &&+ \Frac{1}{12} \Hc{h}^2  \varepsilon [[ \Ang{1}{\partial_x \mom{1}{1}},\ang{1}{2}],
\ang{1}{2}]
- \Frac{1}{12} {\mathrm i} \Hc{h}   \varepsilon^{-1} [[ \ang{\Hc{\mu_{0,0}(\dt)}}{0}, \ang{1}{2}],  \ang{\Hc{\mu_{0,0}(\dt)}}{0}]\\
\nonumber &&- \Frac{1}{6} \Hc{h} [[\ang{\Hc{\mu_{0,0}(\dt)}}{0}, \ang{1}{2}],\Ang{1}{\partial_x
\mom{1}{1}}] - \Frac{1}{6}\Hc{h}  [[
\Ang{1}{\partial_x \mom{1}{1}}, \ang{1}{2}], \ang{\Hc{\mu_{0,0}(\dt)}}{0}]\\
\nonumber && + \Frac{1}{3}{\mathrm i}  \Hc{h} \varepsilon   [[ \Ang{1}{\partial_x \mom{1}{1}},
\ang{1}{2}], \Ang{1}{\partial_x \mom{1}{1}}] + \O{h^5 \varepsilon^{-1}},
\end{Eqnarray}
and we have ignored all $\O{h^5\varepsilon^{-1}}=\O{\varepsilon^{5\sigma-1}}$ terms in the
sBCH since we are applying the Zassenhaus splitting to the fourth-order
Magnus expansion.


Further, we observe that
\begin{Eqnarray}
\label{eq:W1solve2}\mathcal{W}^{[1]} &=& - {\mathrm i} \varepsilon^{-1} \ang{\Hc{\mu_{0,0}(\dt)}}{0} - 2
\Ang{1}{\partial_x \mom{1}{1}} +  \Frac{1}{24}  {\mathrm i} \Hc{h}^2 \varepsilon  [[ \ang{\Hc{\mu_{0,0}(\dt)}}{0} ,\ang{1}{2}], \ang{1}{2}] \\
\nonumber && - \Frac{1}{12} {\mathrm i} \Hc{h}   \varepsilon^{-1} [[ \ang{\Hc{\mu_{0,0}(\dt)}}{0},
\ang{1}{2}],  \ang{\Hc{\mu_{0,0}(\dt)}}{0}] + \O{h^5 \varepsilon^{-1}},
\end{Eqnarray}
since most commutators are $\O{h^5 \varepsilon^{-1}}$ or smaller due to
Theorem~\ref{thm:ve} and the observation that $\Hc{\mu_{0,0}(\dt)} = \O{h}$ and $\mom{1}{1}
= \O{h^3}$.
These nested commutators can be simplified using the identities \R{eq:id},
the detailed theory for which is presented in \cite{psalgebra} and which were
utilised in the derivation of the simplified-commutator Magnus expansion
$\Theta_2$ in \cite{IKS18sinum}. Applying these identities, we find
\begin{Eqnarray}
\label{eq:commutators} [\ang{\Hc{\mu_{0,0}(\dt)}}{0} ,\ang{1}{2}] & = & - 2 \ang{\partial_x \Hc{\mu_{0,0}(\dt)}}{1}, \\
\nonumber\left[ \left[\ang{\Hc{\mu_{0,0}(\dt)}}{0} ,\ang{1}{2}\right], \ang{1}{2}\right] & = &  4 \ang{\partial_x^2 \Hc{\mu_{0,0}(\dt)}}{2} -  \ang{\partial_x^4 \Hc{\mu_{0,0}(\dt)}}{0},\\
\nonumber \left[ \Ang{1}{\partial_x \mom{1}{1}},\ang{1}{2}\right] &= & - 2\ang{\partial_x^2  \mom{1}{1} }{2} + \Frac12 \ang{\partial_x^4 \mom{1}{1}}{0},\\
\nonumber \left[ \left[ \ang{\Hc{\mu_{0,0}(\dt)}}{0}, \ang{1}{2}\right],
\ang{\Hc{\mu_{0,0}(\dt)}}{0}\right] & = &  -2 \ang{(\partial_x \Hc{\mu_{0,0}(\dt)})^2}{0}.
\end{Eqnarray}
Thus, the central exponent in \R{eq:3:sbch_inverted} up to $\O{h^5
\varepsilon^{-1}}=\O{\varepsilon^{5\sigma-1}}$ is
\begin{Eqnarray}\label{eq:3:central_exponent}
\mathcal{W}^{[1]} &=&  - {\mathrm i} \varepsilon^{-1}  \ang{\Hc{\mu_{0,0}(\dt)}}{0} +  \Frac{1}{6} {\mathrm i}
\Hc{h}^2 \varepsilon \ang{\partial_x^2
\Hc{\mu_{0,0}(\dt)}}{2} - 2 \Ang{1}{\partial_x \mom{1}{1}}\\
\nonumber && -  \Frac{1}{24} {\mathrm i} \Hc{h}^2  \varepsilon \ang{\partial_x^4 \Hc{\mu_{0,0}(\dt)}}{0}+\Frac{1}{6}
{\mathrm i} \Hc{h} \varepsilon^{-1}
 \ang{(\partial_x \Hc{\mu_{0,0}(\dt)})^2}{0}.
\end{Eqnarray}
In order to approximate the exponential of $\mathcal{W}^{[1]}$ (which is
required in \R{eq:3:sbch_inverted}), we apply the sBCH again along the lines
of \R{eq:3:sbch_inverted}, writing
\begin{equation}\label{eq:4:sbch_inverted}
  {\mathrm e}^{\mathcal{W}^{[1]}} ={\mathrm e}^{-\frac12  {\mathrm i} \varepsilon^{-1} \Hc{\mu_{0,0}(\dt)}} {\mathrm e}^{\mathcal{W}^{[2]} } {\mathrm e}^{-\frac12  {\mathrm i} \varepsilon^{-1} \Hc{\mu_{0,0}(\dt)}} + \O{h^5\varepsilon^{-1}}.
\end{equation}
Here we have chosen to extract the largest term $X = W^{[1]} = - {\mathrm i} \varepsilon^{-1}
\Hc{\mu_{0,0}(\dt)} = \O{h\varepsilon^{-1}}$ and the remainder is $Y = \mathcal{W}^{[1]}- W^{[1]} =
\O{h^3\varepsilon^{-1}}$. The choice of X is more obvious in this case since there is
only one term of size $\O{h\varepsilon^{-1}}$. Keeping the sizes of $X$ and $Y$ in
mind, we find that
\begin{Eqnarray} \label{eq:4:central_exponent}
\nonumber \mathcal{W}^{[2]} &:=& \CC{sBCH}(-X,X+Y)+ \O{h^5 \varepsilon^{-1}} = (-X + (X+Y)) = Y + \O{h^5 \varepsilon^{-1}}\\
&=& \Frac{1}{6} {\mathrm i} \Hc{h}^2 \varepsilon \ang{\partial_x^2 \Hc{\mu_{0,0}(\dt)}}{2} - 2 \Ang{1}{\partial_x
\mom{1}{1}}\\
\nonumber && \mbox{}-  \Frac{1}{24}{\mathrm i} \Hc{h}^2 \varepsilon \ang{\partial_x^4 \Hc{\mu_{0,0}(\dt)}}{0}+\Frac{1}{6}
{\mathrm i} \Hc{h} \varepsilon^{-1}
 \ang{(\partial_x \Hc{\mu_{0,0}(\dt)})^2}{0} + \O{h^5 \varepsilon^{-1}}
\end{Eqnarray}
suffices up to $\O{h^5 \varepsilon^{-1}}$ since $[[Y,X],X]$ and $[[Y,X],Y]$ are
$\O{h^5\varepsilon^{-1}}$ and $\O{h^7\varepsilon^{-1}}$, respectively. Further, we discard
the $\O{h^3 \varepsilon}$ term $-  \Frac{1}{24}{\mathrm i} \Hc{h}^2 \varepsilon \ang{\partial_x^4 \Hc{\mu_{0,0}(\dt)}}{0}$
from \R{eq:4:central_exponent} since, under the assumption $\sigma \leq 1$,
it is effectively $\O{\varepsilon^{5\sigma-1}}$ or smaller.

Finally, we complete the derivation of our fourth-order Magnus--Zassenhaus
splitting by combining \R{eq:3:sbch_inverted} and \R{eq:4:sbch_inverted},
\begin{equation}
\label{eq:MZ2}
	 \exp(\Theta_{2}) = {\mathrm e}^{\frac12 W^{[0]}} {\mathrm e}^{\frac12 W^{[1]}}
									 {\mathrm e}^{\mathcal{W}^{[2]}} {\mathrm e}^{\frac12 W^{[1]}} {\mathrm e}^{\frac12 W^{[0]}} + \O{\varepsilon^{5\sigma-1}},
\end{equation}
which has the same form as the splittings described in
Appendix~\ref{sec:Zassenhaus}. Here the exponents (chosen and/or computed
during the derivation) are
\begin{Eqnarray}
W^{[0]} & = & {\mathrm i} \Hc{h} \varepsilon \partial_x^2  = \O{\varepsilon^{\sigma-1}},\label{eq:MZ2W} \\
W^{[1]} & = & - {\mathrm i}  \varepsilon^{-1}  \Hc{\mu_{0,0}(\dt)} = \O{\varepsilon^{\sigma-1}}, \nonumber\\
\mathcal{W}^{[2]} & = & \Frac{1}{6}{\mathrm i}  \Hc{h} \varepsilon^{-1} \left(\partial_x \Hc{\mu_{0,0}(\dt)} \right)^2
-2\Ang{1}{\partial_x \mom{1}{1}}+\Frac{1}{6}{\mathrm i}  \Hc{h}^2 \varepsilon \Ang{2}{\partial_x^{2} \Hc{\mu_{0,0}(\dt)} } =
\O{\varepsilon^{3\sigma-1}}.\nonumber
\end{Eqnarray}

It is only by this stage -- having arrived at an asymptotic splitting
expressed in operatorial terms -- that we start considering discretisation
issues. Apart from analysis of size, considerations of spatial and temporal
discretisation (in the form of approximation through quadrature) are entirely
independent of the procedure leading to \R{eq:MZ2W} and of each other: one may proceed to address them in any order.


\begin{remark}
Note that, in contrast to the Magnus--Zassenhaus splitting of \cite{BIKS},
our exponents retain the integrals intact (in $\Hc{\mu_{0,0}(\dt)}$ and $\mom{1}{1}$). This
provides us with the flexibility to choose the most suitable approximation method
 including the possibility of analytic approximations. This
flexibility proves particularly helpful in the case of potentials with high
temporal oscillations.
\end{remark}

\begin{remark}
As a quick sanity check, we do arrive at the standard symmetric Zassenhaus
splitting \cite{bader14eaf} for time-independent potentials, $V(x,t) :=
V(x)$.
\end{remark}

\subsection{Implementation}

\label{sec:MZ2implementation} The operatorial splitting given by \R{eq:MZ2}
and \R{eq:MZ2W} can be converted into a concrete numerical scheme by
resorting to spatial discretisation. As typical, we use spectral collocation,
whereby $W^{[0]}$ is discretised as the circulant matrix $\widetilde{W^{[0]}}
= {\mathrm i} \Hc{h} \varepsilon \mathcal{K}_2$. Note that $\mathcal{K}_2$ is diagonalisable via Fourier transform
since
\begin{equation}
\label{eq:KkFFT}
\mathcal{K}_k = \mathcal{F}^{-1} \mathcal{D}_{c_k} \mathcal{F},
\end{equation}
where $c_k$ is the
symbol of $\mathcal{K}_k$, $\mathcal{F}$ is the Fast Fourier Transform (FFT) and $\mathcal{F}^{-1}$
is its inverse (IFFT). Thus, we can evaluate
\[
{\mathrm e}^{\frac12 \widetilde{W^{[0]}}} \MM{u}  = \mathcal{F}^{-1} \mathcal{D}_{\exp(\frac12 {\mathrm i} \Hc{h} \varepsilon c_2)} \mathcal{F} \MM{u},
\]
in $\O{N \log N}$ operations. The next exponent, $W^{[1]}$, is discretised as
the diagonal matrix $\widetilde{W^{[1]}} = - {\mathrm i} \varepsilon^{-1} \mathcal{D}_{\Hc{\mu_{0,0}(\dt)}}$ and
exponentiated pointwise in $\O{N}$ operations,
\[
{\mathrm e}^{\frac12 \widetilde{W^{[1]}}} \MM{u} = \mathcal{D}_{\exp(\frac12 {\mathrm i} \varepsilon^{-1} \Hc{\mu_{0,0}(\dt)})} \MM{u}.
\]
Lastly, the innermost exponent,
\begin{Eqnarray*}
\mathcal{W}^{[2]} \leadsto  \widetilde{\mathcal{W}^{[2]}} & = & \Frac{1}{6}{\mathrm i}  \Hc{h} \varepsilon^{-1}
\mathcal{D}_{\left(\partial_x \Hc{\mu_{0,0}(\dt)} \right)^2} +\left(\mathcal{K}_1 \mathcal{D}_{\partial_x \mom{1}{1}}+ \mathcal{D}_{\partial_x
\mom{1}{1}}\mathcal{K}_1\right)\\
&&-\Frac{1}{3}{\mathrm i} \Hc{h}^2 \varepsilon \left(\mathcal{K}_2 \mathcal{D}_{\partial_x^{2} \Hc{\mu_{0,0}(\dt)} } + \mathcal{D}_{\partial_x^{2}
\Hc{\mu_{0,0}(\dt)} } \mathcal{K}_2\right),
\end{Eqnarray*}
which is neither circulant nor diagonal, can be very effectively
exponentiated via Lanczos iterations because (unlike the first two exponents)
it is very small in size, even when combined with large time steps.

For instance, when using time steps as large as $h = \O{\varepsilon^{1/2}}$, i.e.\
under the scaling $\sigma = 1/2$, the exponent $\widetilde{\mathcal{W}^{[2]}}$ is
$\O{\varepsilon^{1/2}}$. As the semiclassical parameter $\varepsilon \ll 1$ becomes small,
this spectral radius approaches $0$. Thus the superlinear accuracy of Lanczos
iterations is seen immediately and merely four Lanczos iterations are required
for an exponentiation to an accuracy of $\O{\varepsilon^{5\sigma-1}}=\O{\varepsilon^{3/2}}$,
which is the accuracy of the splitting. A more detailed analysis of the
number of Lanczos iterations required for such exponents is carried out in
\cite{bader14eaf,BIKS}.

Each Lanczos iteration requires the computation of a matrix--vector product
of the form $\widetilde{\mathcal{W}^{[2]}} \MM{v}$. In general Lanczos iterations for
the discretisation of $W = \sum_{k=0}^n {\mathrm i}^{k+1} \ang{f_k}{k}$ to
$\widetilde{W} = {\mathrm i} \mathcal{D}_{f_0} + \frac12 \sum_{k=1}^n {\mathrm i}^{k+1}
\left(\mathcal{D}_{f_k} \mathcal{K}_k +\mathcal{K}_k \mathcal{D}_{f_k} \right)$ can be efficiently computed
in $2n+2$ Fast Fourier Transforms (included inverses) using the expression
 \begin{Eqnarray}
\label{eq:Wv}  \widetilde{W}\MM{v} = {\mathrm i} \mathcal{D}_{f_0} \MM{v} + \Frac12 \left(
\sum_{k=1}^n {\mathrm i}^{k+1} \mathcal{D}_{f_k} \mathcal{F}^{-1} \mathcal{D}_{c_k} \right) \mathcal{F} \MM{v} +
\Frac12 \mathcal{F}^{-1} \left( \sum_{k=1}^n {\mathrm i}^{k+1}  \mathcal{D}_{c_k} \mathcal{F} \mathcal{D}_{f_k}
\MM{v}\right).
\end{Eqnarray}
Thus, the computation of $\widetilde{\mathcal{W}^{[2]}} \MM{v}$ in each Lanczos
iteration requires six FFTs and the total cost of its exponentiation is 24
FFTs. This completes the derivation of the splitting described in
\R{eq:Zass2}.

\begin{remark}
$\Theta_0$ is an order two Magnus expansion. An order two
Zassenhaus splitting is merely the Strang splitting. Thus an order two
Magnus--Zassenhaus splitting is, trivially, \R{eq:expmidpoint}, which has been
presented in the introduction.
\end{remark}

\begin{remark}
We assume that the integrals $\Hc{\mu_{0,0}(\dt)}$ and $\mom{1}{1}$ are either available
analytically or can be adequately approximated via quadrature formulae such
as Gauss--Legendre quadratures, along the lines of \cite{IKS18sinum}.

\end{remark}

\begin{remark} Along the lines of \cite{BIKS,IKS18sinum}, we also assume that spatial derivatives of the potential, as well as
its integrals, are either available or inexpensive to compute.
\end{remark}











\subsection{Size analysis}

Size of the terms is analysed with the aid of Lemmas~\ref{lem:powersofh_mu}
and \ref{lem:powersofh_lambda}.

\begin{lemma} \label{lem:powersofh_mu}
    For an analytic potential $V$,
   
    \begin{equation}\label{eq:powersofh_mu}
    \mom{j}{k} = \O{\Hc{h}^{jk+1}}, \quad k\neq 1, \qquad \mom{j}{1} = \O{h^{j+2}}.
    \end{equation}

    \begin{proof} \label{pf:powersofh_mu}
    Since the $j${\em th} rescaled Bernoulli polynomial scales as $\O{\Hc{h}^j}$, we
    expect $\mom{j}{k} = \O{\Hc{h}^{jk+1}}$ . The term $\mom{j}{1}$ gains an extra power of $\Hc{h}$ and is
    $\O{\Hc{h}^{j+2}}$ since  integrals of  Bernoulli   polynomials vanish,
    $\Int{\zeta}{0}{\Hc{h}}{\tilde{B}_j(\Hc{h},\zeta) }= 0$. To see this, consider the Taylor expansion,
    $ V(\zeta) = \sum_{n=0}^\infty [V^{(n)}/n!] \zeta^n$,
    and note that in the expansion for $\mom{j}{1}$,
    the leading order term $ \Int{\zeta}{0}{\Hc{h}}{V^{(0)}\tilde{B}_j(\Hc{h},\zeta) }= 0$ vanishes.

    \end{proof}
\end{lemma}


\begin{lemma} \label{lem:powersofh_lambda}
    In general, for a polynomial $r_n(\Hc{h},\zeta,\xi)$ featuring only degree-$n$
    terms in $\Hc{h}, \zeta$ and $\xi$, the linear (integral) functional is
    \begin{equation} \label{eq:powersofh_lambda1}
    \Lf{r_n}{a}{b}=\O{\Hc{h}^{n+2}}.
    \end{equation}
    For a polynomial $p_n$ whose integral over
    the triangle vanishes,
    $\Int{\zeta}{0}{\Hc{h}}{\Int{\xi}{0}{\zeta}{p_n(\Hc{h},\zeta,\xi)}}=0$,
    we gain an extra power of $\Hc{h}$ and
    \begin{equation} \label{eq:powersofh_lambda1}
    \Lf{p_n}{a}{b} = \O{\Hc{h}^{n+3}}.
    \end{equation}
    In particular, this happens for $\psi, \phi_1$ and $\phi_2$, but not for
    $\varphi_1$ and $\varphi_2$.
\end{lemma}










\subsection{Reduction of stages in Zassenhaus splitting}

The expansion \R{eq:Th4} features both odd and even powers of $h$ -- in
particular, it features terms of size $\O{h^4 \varepsilon^{-1}}$. When we commence
the Zassenhaus procedure  directly from $\Theta_4$, the resulting exponential splitting
\begin{equation}
\label{eq:Zass4suboptimal}
 {\mathrm e}^{\frac12 W^{[0]}} {\mathrm e}^{\frac12 W^{[1]}} {\mathrm e}^{\frac12 W^{[2]}} {\mathrm e}^{\frac12 W^{[3]}}  {\mathrm e}^{\frac12 W^{[4]}}{\mathrm e}^{\mathcal{W}^{[5]}}{\mathrm e}^{\frac12 W^{[4]}} {\mathrm e}^{\frac12 W^{[3]}}  {\mathrm e}^{\frac12 W^{[2]}}{\mathrm e}^{\frac12 W^{[1]}}{\mathrm e}^{\frac12 W^{[0]}},
\end{equation}
features exponents $W^{[3]}=\O{h^4 \varepsilon^{-1}}$ and $\mathcal{W}^{[5]}=\O{h^6
\varepsilon^{-1}}$ and the power of $h$ is even. Similarly, symmetric Zassenhaus
splitting of a higher order Magnus expansion such as $\Theta_8$ also
features $\O{h^4 \varepsilon^{-1}}$, $\O{h^6 \varepsilon^{-1}}$, $\O{h^8 \varepsilon^{-1}}$ and
$\O{h^{10} \varepsilon^{-1}}$ exponents.

In contrast, the Zassenhaus splittings of \cite{bader14eaf} and \cite{BIKS},
only feature exponents in odd powers of $h$: even powers of
$h$ do not appear in \cite{BIKS} since the integrals in the Magnus expansion
are discretised at the outset by Gaussian quadrature. The appearance of
even powers of $h$ is suboptimal because it results in a larger number of
exponentials here.

To remedy this, we note that the results of \cite{IKS18sinum}, from which we
have commenced our analysis, are based on power truncated Magnus expansions.
These methods are known to be odd in $\Hc{h}$ around the midpoint of the
interval ($t_{1/2} = t + h/2$) due to time symmetry of the flow
\cite{iserles00lgm,iserles01tsa}.

It should, therefore, be possible to expand $\Theta_m$ solely in odd powers
of $\Hc{h}$ for any choice of the potential $V$. A Zassenhaus splitting
commencing from such an odd-powered expansion will never introduce even
powers of $\Hc{h}$ since underlying the procedure is a recursive application of
the symmetric BCH which features only odd-grade commutators. Consequently, in
contrast to the eleven exponentials appearing in \R{eq:Zass4suboptimal}, the
Magnus--Zassenhaus splitting \R{eq:Zass4} that we are about to develop later in this
section features just seven exponentials.

\begin{remark} The extra stages appearing in \R{eq:Zass4suboptimal} can also
be seen as a consequence of lack of symmetry of the scheme.
\end{remark}

\subsection{Time symmetry and shifting the origin}

In this section we derive an expression for $\Theta_4$
expanded solely in odd powers of $\Hc{h}$. Recall that, in order to advance from
time $t$ to $t+h$ we use the truncated Magnus expansion,
\begin{displaymath}
u(t+\Hc{h}) = {\mathrm e}^{\Theta_m(t+\Hc{h},t)}u(t),
\end{displaymath}
where $\Theta_m(t+\Hc{h},t)$ can be easily recovered from $\Theta_m(\Hc{h},0)$
(shortened to $\Theta_m(\Hc{h})$) by substituting all occurrences of $V(\zeta)$
by $V(t+\zeta)$. Let us define the midpoint of the interval $[t,t+\Hc{h}]$ as
$t_{1/2} = t + \Hc{h}/2$. Recalling that the expansion is odd about
$t_{1/2}$, we shift the origin to $t_{1/2}$ by defining
\[W(\zeta):= V(t_{1/2}+\zeta),\]
whereby we need to substitute $V(\zeta)$ with $W(\zeta-\Hc{h}/2)$ in the
Magnus expansion. Specifically, for methods described here, this only needs
to be done in the terms $\mom{j}{k}$ and $\Lf{f}{a}{b}$,
\begin{Eqnarray}
\label{eq:momN}
\mom{j}{k} &:=& \Int{\zeta}{0}{\Hc{h}}{\tilde{B}_{j}^k(\Hc{h},\zeta)\,W\!\left(\zeta-\frac{\Hc{h}}{2}\right)\! },\\
\label{eq:LfN}
\Lf{f}{a}{b} &:=& \Int{\zeta}{0}{\Hc{h}}{\Int{\xi}{0}{\zeta}{f(\Hc{h},\zeta,\xi) \left[ \partial_x^a W\!\left(\zeta-\frac{\Hc{h}}{2}\right)\ \partial_x^b W\!\left(\xi-\frac{\Hc{h}}{2}\right) \right] \!}}.
\end{Eqnarray}
As we demonstrate in the rest of this section, this substitution allows us to identify and discard even-powered terms in the Magnus expansion.

\begin{remark}
Since we have shifted the origin to $t_{1/2}$, all odd and even components are to be understood with respect to $0$ from this point onwards. This makes identification of the odd components of the Magnus expansion simpler, assuming that the odd and even components of $W$ can be found.
\end{remark}





\subsection{Odd and even components of $\mu$ and $\Lambda$}
\label{sec:oddeven_mu_lambda} Extending the new definitions of $\mu$ and
$\Lambda$ (\ref{eq:momN}, \ref{eq:LfN}), we define
\begin{equation}
\label{eq:momS}
\momS{\star}{j}{k} := \Int{\zeta}{0}{\Hc{h}}{\tilde{B}_{j}^k(\Hc{h},\zeta)\,W^\star\!\left(\zeta-\frac{\Hc{h}}{2}\right) },\quad \star \in \{e,o\},
\end{equation}
and
\begin{equation}
\label{eq:LfS}
\LfS{f}{\star}{a}{b} := \Int{\zeta}{0}{\Hc{h}}{\Int{\xi}{0}{\zeta}{f(\Hc{h},\zeta,\xi) \left[ \partial_x^a W\left(\zeta-\frac{\Hc{h}}{2}\right)\ \partial_x^b W\left(\xi-\frac{\Hc{h}}{2}\right) \right]^\star}},\quad \star \in \{e,o\}.
\end{equation}
The even and odd components of $\mom{j}{k}$ and $\Lf{f}{a}{b}$ can easily be
expressed in terms of these new definitions,
\[
\begin{tabular}{l l}
 
  $\left(\mom{j}{k}\right)^{o} = \begin{cases} \momS{e}{j}{k} &\mbox{if } jk \mbox{ is even},\\
\momS{o}{j}{k} &\mbox{if } jk \mbox{ is odd}, \end{cases}$
& \quad$\left(\mom{j}{k}\right)^{e}  =  \begin{cases} \momS{o}{j}{k} &\mbox{if } jk \mbox{ is even},\\
\momS{e}{j}{k} &\mbox{if } jk \mbox{ is odd}, \end{cases}$
\end{tabular} \]
where we exploit the fact that $\tilde{B}_j (\Hc{h}, \zeta)$ is odd in $(\Hc{h},
\zeta)$ for odd values of $j$ and even for even values. For instance,
$\Hc{\mu_{0,0}(\dt)}^o = \momOS{e} =
\Int{\zeta}{0}{\Hc{h}}{W^e\!\left(\zeta-\frac{\Hc{h}}{2}\right)}$. Note carefully
the difference of usage in subscript and superscript. The odd and even parts
of $\Lf{f}{a}{b}$ are
\begin{Eqnarray*}
\left(\Lf{f}{a}{b}\right)^{\!o} &=&  \LfS{f^o}{e}{a}{b}+ \LfS{f^e}{o}{a}{b},\\
\left(\Lf{f}{a}{b}\right)^{\!e} &=&  \LfS{f^o}{o}{a}{b}+ \LfS{f^e}{e}{a}{b}.
\end{Eqnarray*}
For an odd function such as $\psi(\Hc{h},\zeta,\xi): = \zeta - \xi -
\Frac{\Hc{h}}{3}$, for instance, this reduces to
\[\left(\Lf{\psi}{a}{b}\right)^{\!o}  = \LfS{\psi}{e}{a}{b}. \]

\subsection{Expanding $\mu$ and $\Lambda$ in powers of $h$.}

\label{sec:powersofh_oddeven_mu_lambda} A natural consequence of separating
odd and even components is that the sizes of certain $\mu$ and $\Lambda$
integrals become smaller (i.e.\ they are associated with higher powers of $h$) than prescribed by
Lemmas~\ref{lem:powersofh_mu} and \ref{lem:powersofh_lambda}. Essentially, this
is due to certain terms disappearing in the relevant Taylor expansions.
\begin{lemma}\label{lem:powersofh_oddeven_mu}
    For an analytic potential $V$  (and thus $W$),
    \begin{equation} \label{eq:powersofh_even_mu}
        \momS{e}{j}{k} = \O{h^{jk+1}}, \quad k \neq 1, \qquad \momS{e}{j}{1} = \O{h^{j+3}},
    \end{equation}
    and
    \begin{equation} \label{eq:powersofh_odd_mu}
        \momS{o}{j}{k} = \O{h^{jk+2}}, \quad k \geq 0.
    \end{equation}
   


    \begin{proof}
    \label{pf:powersofh_oddeven_mu} For an analytic potential,
    $W\left(\zeta-\Hc{h}/2\right) = \sum_{n=0}^{\infty} [W^{(n)}(0)/n!] \left(\zeta-\Hc{h}/2\right)^n $,
    we expect a gain of a single power of $\Hc{h}$ in $\mom{j}{1}$ since
    $\Int{\zeta}{0}{\Hc{h}}{ \tilde{B}_j(\Hc{h},\zeta) W(0)}$ vanishes (see
    Lemma~\ref{lem:powersofh_mu}). Thus, instead of being $\O{\Hc{h}^{j+1}}$, this
    term turns out to be $\O{\Hc{h}^{j+2}}$.

    In the case of $\momS{e}{j}{1}$ where we take the even part of the potential
    $W$ at $\zeta-\frac{\Hc{h}}{2}$, we see a gain of two powers of $\Hc{h}$ since the
    smallest non-vanishing term is
    \[\frac12 \Int{\zeta}{0}{\Hc{h}}{ B_j(\Hc{h},\zeta) W^{(2)}(0) \left(\zeta-\frac{\Hc{h}}{2}\right)^2}.\]
    Thus, $\momS{e}{j}{1} = \O{h^{j+3}}$ instead of being $\O{h^{j+2}}$. This is
    the only exception to Lemma~\ref{lem:powersofh_mu} for $\momS{e}{j}{k}$.
    Similarly, for $\momS{o}{j}{k}$, observe that there is no $W^{(0)}$ term
    in the Taylor expansion.
   
   
   
   
    \end{proof}
\end{lemma}



\begin{lemma} \label{lem:powersofh_oddeven_lambda}
For a polynomial $r_n(\Hc{h},\zeta,\xi)$ featuring only degree-$n$ terms in
$\Hc{h}, \zeta$ and $\xi$,
\begin{equation} \label{eq:powersofh_even_lambda_r}
    \LfS{r_n}{e}{1}{1} = \O{h^{n+2}},
\end{equation}
while for a polynomial $p_n$ whose integral over the triangle vanishes,
\begin{equation*}
\label{eq:IntPhi}
 \Int{\zeta}{0}{\Hc{h}}{\Int{\xi}{0}{\zeta}{p_n(\Hc{h},\zeta,\xi)}}=0,
\end{equation*}
\begin{equation} \label{eq:powersofh_even_lambda_p}
    \LfS{p_n}{e}{1}{1} = \O{h^{n+4}}.
\end{equation}
In particular, the exception \R{eq:powersofh_even_lambda_p} holds for $\psi,
\phi_1$ and $\phi_2$, but not for $\varphi_1$ and $\varphi_2$. For the odd
part,

\begin{equation} \label{eq:powersofh_even_lambda_r}
\LfS{r_n}{o}{1}{1} = \O{h^{n+3}}, \end{equation}
whether the integral of $r_n$ vanishes or not.
\end{lemma}
The proof for Lemma~\ref{lem:powersofh_oddeven_lambda} follows directly
along the lines of Lemma~\ref{lem:powersofh_oddeven_mu}.

\subsection{Identifying the odd components of Magnus expansions}
\label{sec:extracting_odd_Magnus} Recall that the time symmetry of power
truncated Magnus expansions ensures that, having shifted the origin to $t_{1/2}$, $\Theta_m$ is odd about the origin.
Therefore the even part
$\Theta_{m}^{e}$ vanishes, and it suffices to work with $\Theta_{m} =
\Theta_{m}^{o}$.
For the sixth-order simplified-commutator Magnus expansion of
\cite{IKS18sinum}, $\Theta_4$, we can identify the odd components using the
results of Subsection~\ref{sec:oddeven_mu_lambda},
\begin{Eqnarray}
\label{eq:Theta_4o}
\Theta_4 = \Theta_4^o & = &  \overbrace{{\mathrm i} \Hc{h} \varepsilon \partial_x^2 - {\mathrm i} \varepsilon^{-1} \momOS{e}}^{\O{h \varepsilon^{-1}}} -  \overbrace{2 \Ang{1}{\partial_x \momS{o}{1}{1}}}^{\O{h^3 \varepsilon^{-1}}}\\
\nonumber &&\mbox{} + \overbrace{{\mathrm i} \varepsilon^{-1} \LfS{\psi}{e}{1}{1} + 2{\mathrm i} \varepsilon \Ang{2}{\partial_x^2 \momS{e}{2}{1}}}^{\O{h^5 \varepsilon^{-1}}} \\
\nonumber &&\mbox{} +\overbrace{\Frac{1}{6}\varepsilon^{-1} \Ang{1}{\LfS{\varphi_1 + \phi_1}{o}{1}{2} + \LfS{\varphi_2 + \phi_2}{o}{2}{1} }}^{\O{h^5 \varepsilon^{-1}}} \\
\nonumber &&\mbox{}+\overbrace{\Frac{4}{3} \varepsilon^2 \Ang{3}{\partial_x^3
\momS{o}{3}{1}}}^{\O{h^5 \varepsilon^{-1}}} +\overbrace{\Frac{1}{4} {\mathrm i} \varepsilon \partial_x^4
\momS{e}{2}{1}}^{\O{h^5 \varepsilon}} = \mTh{} + \O{h^7 \varepsilon^{-1}},
\end{Eqnarray}
where size analysis in powers of $h$ is done using
Lemmas~\ref{lem:powersofh_oddeven_mu} and \ref{lem:powersofh_oddeven_lambda},
and we have used the fact that $\left(\Lf{\varphi_j + \phi_j}{a}{b}\right)^o
= \LfS{\varphi_j + \phi_j}{o}{a}{b}$, since $\varphi_j + \phi_j$ is even. We
can now discard the $\O{h^5 \varepsilon}$ term $\Frac{1}{4} {\mathrm i} \varepsilon \partial_x^4
\momS{e}{2}{1}$ since it is $\O{\varepsilon^{7\sigma-1}}$ or smaller under the
assumption $\sigma\leq 1$.

By this stage we have been able to eradicate all even powers of $\Hc{h}$ and are
now in a position to effectively apply the Zassenhaus splitting to the
sixth-order Magnus expansion $\Theta_4^o$.


\subsection{Sixth-order symmetric Zassenhaus splitting}
A Zassenhaus splitting commencing from the Magnus expansion $\Theta_4^o$
naturally yields an exponential splitting featuring only $\O{\varepsilon^{\sigma
-1}}$, $\O{\varepsilon^{3\sigma -1}}$ and $\O{\varepsilon^{5\sigma -1}}$ terms,
\begin{equation}
\label{eq:Zass4}
 \exp\left(\Theta_4\right) = {\mathrm e}^{\frac12 W^{[0]}} {\mathrm e}^{\frac12 W^{[1]}} {\mathrm e}^{\frac12 W^{[2]}} {\mathrm e}^{\mathcal{W}^{[3]}} {\mathrm e}^{\frac12 W^{[2]}}{\mathrm e}^{\frac12 W^{[1]}}{\mathrm e}^{\frac12 W^{[0]}} + \O{\varepsilon^{7\sigma -1}} ,
\end{equation}
where
\begin{Eqnarray*}
W^{[0]}  &=& {\mathrm i}  \Hc{h} \varepsilon \partial_x^2=\O{\varepsilon^{\sigma-1}},\\
W^{[1]} &=& -{\mathrm i}  \varepsilon^{-1}  \momOS{e} =\O{\varepsilon^{\sigma-1}},\\
W^{[2]} &=& \Frac16{\mathrm i}  \Hc{h} \varepsilon^{-1}  \left(\partial_x \momOS{e}\right)^2   - 2 \Ang{1}{\partial_x \momS{o}{1}{1}} + \Frac{1}{6}{\mathrm i}  \Hc{h}^{2} \varepsilon \Ang{2}{ \partial_x^2 \momOS{e} }=\O{\varepsilon^{3\sigma-1}},\\
\mathcal{W}^{[3]} &=& {\mathrm i}  \varepsilon^{-1}  \LfS{\psi}{e}{1}{1}  - \Frac{1}{24}{\mathrm i}  \Hc{h}^{2} \varepsilon \left(\partial_x^{4}  \momOS{e} \right)
-\Frac{1}{6}{\mathrm i}  \varepsilon^{-1}  \left(\partial_x \momOS{e}\right)^2 \left(\partial_x^2 \momS{e}{2}{1}\right)\\
&&\mbox{} - \Frac{2}{45} {\mathrm i} \Hc{h}^2 \varepsilon^{-1} \left(\partial_x \momOS{e} \right)^2 \left(\partial_x^2 \momOS{e} \right) \\
&& +\Frac16 \Ang{1}{\LfS{\varphi_1+\phi_1}{o}{1}{2} + \LfS{\varphi_2+\phi_2}{o}{2}{1} }\\
&&\mbox{}- \Hc{h} \Ang{1}{\left(\partial_x \momOS{e}\right) \left(\partial_x^{2} \momS{o}{1}{1} \right) -\Frac13 \left(\partial_x^2 \momOS{e}\right) \left(\partial_x \momS{o}{1}{1} \right)} \\
&&\mbox{}+\Frac{1}{30} {\mathrm i} \Hc{h}^3 \varepsilon \Ang{2}{\left(\partial_x^2 \momOS{e}\right)^2 - 2\left(\partial_x \momOS{e}\right) \left(\partial_x^3 \momOS{e}\right)} \\
&&\mbox{}+2 {\mathrm i}  \varepsilon\Ang{2}{\partial_x^{2} \momS{e}{2}{1}} + \Frac{4}{3} \varepsilon^{2}\Ang{3}{\partial_x^3 \momS{o}{3}{1}} \\
&&\mbox{}+\Frac13 \Hc{h}^2 \varepsilon^2 \Ang{3}{\partial_x^3 \momS{o}{1}{1}} - \Frac{1}{120}
{\mathrm i} \Hc{h}^4 \varepsilon^3 \Ang{4}{\partial_x^4 \momOS{e}} = \O{\varepsilon^{5\sigma-1}},
\end{Eqnarray*}

The term $\Frac{1}{24}{\mathrm i}  \Hc{h}^{2} \varepsilon \left(\partial_x^{4}  \momOS{e} \right)$ in
$\mathcal{W}^{[3]}$ is $\O{\varepsilon^{3\sigma+1}}$. We remind the reader that this
splitting is obtained subject to $\sigma\leq 1$. For $1/2 <\sigma \leq 1$,
this term can be combined with the $\O{\varepsilon^{5\sigma-1}}$ terms in
$\mathcal{W}^{[3]}$, whereas for $\sigma \leq 1/2$ it can be ignored since it is
smaller than $\O{\varepsilon^{7\sigma-1}}$.

\subsection{Remarks concerning implementation}
The approximation of the exponentials in \R{eq:Zass4} can be done along the
same lines as Subsection~\ref{sec:MZ2implementation}. For instance, for
$h=\O{\varepsilon^{1/2}}$, the exponents $W^{[2]}$ and $\mathcal{W}^{[3]}$ are
$\O{\varepsilon^{1/2}}$ and $\O{\varepsilon^{3/2}}$, respectively. Approximating their
exponentials to an accuracy of $\O{\varepsilon^{7\sigma-1}} = \O{\varepsilon^{5/2}}$ requires
six and two Lanczos iterations, respectively. Each Lanczos iteration requires
the computation of a matrix--vector product via \R{eq:Wv}.

In  each time step of the order-six scheme we need to approximate five line
integrals
($\momOS{e}$, $\momS{o}{1}{1}$, $\momS{e}{1}{2}$, $\momS{e}{2}{1}$, $\momS{o}{3}{1}$) %
and three integrals over the triangle ($\LfS{\phi}{e}{1}{1}$, $\LfS{\psi_1 +
\theta_1}{o}{1}{2}$ and $\LfS{\psi_2 + \theta_2}{o}{2}{1}$).
If analytic expressions are not available, it is possible to approximate
these through quadrature formulae such as Gauss--Legendre quadrature. For the
$\O{\varepsilon^{7 \sigma -1}}$ splitting here, for instance, we require merely three
Gauss--Legendre knots $t_k = h (1 + k \sqrt{3/5})/2 , k=-1,0,1$, with weights
$w_k = \Frac{5}{18}\Hc{h}, \Frac49 \Hc{h}, \Frac{5}{18} \Hc{h}$ \cite{davis84quad},
respectively, for the $\O{\Hc{h}^7}$ accuracy which is required of the
quadrature.

Evaluation for the integrals over the line, $\momS{\star}{j}{k}$, follows
directly using these quadrature knots. Integrals over the triangle are also
easily evaluated using the same three knots by substituting $W$ with its
interpolating polynomial. A more detailed narrative on the approximation of
such integrals can be found in \cite{IKS18sinum}.












