matrix hydrodynamics, backward error analysis, Butcher series, biplanar forests, Zeitlin’s model, 2-D Euler, Hamiltonian PDEs, Lie–Poisson reduction, isospectral flows, symplectic Runge–Kutta methods

1 Introduction↩︎

The success of symplectic numerical integration schemes for long-term simulations of Hamiltonian systems can conceptually be explained as follows. Given a symplectic phase space \(M\) and a Hamiltonian system on it \[\label{eq:ham95sys} \dot{y} = X_H(y) \qquad \big(\text{if M=\mathbb{R}^{2d} then X_H =J^{-1}\nabla H for J = \begin{bmatrix}0 & -I \\ I & 0 \end{bmatrix}}\big),\tag{1}\] let \(\Phi_h\colon M\to M\) denote a numerical integration method for the system 1 . If we try to reinterpret the discrete integrator map \(\Phi_h\) as the exact flow of a modified vector field \(\tilde{X}_h\) (this is the basic idea of backward error analysis), then, if \(\Phi_h\) is a symplectic map, the vector field \(\tilde{X}_h\) must be symplectic. In particular, it (locally) corresponds to a modified Hamiltonian function \(\tilde{H}_h\) on \(M\), namely \(\tilde{X} = X_{\tilde{H}_h}\). This means that the discrete trajectory traced out by \(y_{k+1}= \Phi_h(y_k)\) corresponds to the exact flow of a true Hamiltonian system. All properties shared by solutions to Hamiltonian systems are thereby shared by the numerical trajectory. More specifically, if the method is consistent and the step size is small enough, all properties shared by Hamiltonian systems nearby the system 1 are shared by the numerical trajectory. For example, if the system 1 is integrable in the Arnold–Liouville sense, KAM-theory tells us that nearby Hamiltonian system are “almost” integrable (cf. [1]), leading to near conservation of the first integrals. As another and more generic example, as long as the numerical trajectory remains on a compact subset of \(M\), we have that \(H(y_k) - \tilde{H}_h(y_k) = \mathcal{O}(h^p)\) where \(p\) is the order of the method. Consequently, since \(\tilde{H}_h\) is exactly preserved by the modified vector field, the original Hamiltonian \(H\) is nearly conserved.

This conceptual story of backward error analysis for symplectic integrators is pleasing and the consequences drawn from it are essentially correct for a finite-dimensional phase space \(M\). But the story is not quite right; the devil is in the details. The core issue is a deep result in the analysis of infinite-dimensional Lie groups. The set of diffeomorphisms on \(M\) form a Lie group (for example in the Fréchet topology of smooth functions, cf. [2]). Its Lie algebra consists of vector fields on \(M\). The corresponding exponential map is obtained by integrating the vector field to a diffeomorphism at time one. And herein lies the core issue: this exponential map is never locally surjective (contrary to finite-dimensional Lie groups where it is always locally surjective). Thus, there are diffeomorphisms arbitrarily close to the identity map that cannot be generated by time one integration of a vector field. The same is true for the subgroup of symplectic diffeomorphisms and its Lie algebra of symplectic vector fields. Thus, the very outset of backward error analysis, that there exists a modified vector field, is invalid in general. As a result, modified Hamiltonians do not exist in general. Instead, the rigorous theory of backward error analysis looks for an “almost” modified Hamiltonian, obtained as an optimally truncated Taylor expansion \[\tilde{H}_h = H + h H_1 + \cdots + h^{N-1}H_N.\] This theory, given by Benettin and Giorgilli [3], is technical. In particular, it requires a real analytic setting and, because of the truncation, the results are valid only for exponentially long time-intervals \(t \leq e^{c/h}\), where the constant \(c > 0\) depends on the vector field \(X_H\), the numerical method, and the compact domain of the numerical trajectory, but it is independent of the step size \(h\). In practise, the exponentially long time-interval means that other numerical artifacts, such as drift due to round-off errors, will dominate over the truncation error in the modified Hamiltonian. The conceptual backward error analysis story, as told in the first paragraph, is thereby valid, in essence, for finite-dimensional Hamiltonian systems.

At first sight, the situation for infinite-dimensional Hamiltonian PDEs seems similar: if the spatial discretization gives a finite-dimensional Hamiltonian system we can use a symplectic integration scheme to obtain near conservation of energy. But this approach fails, because all error constants in the estimates depend on the spatial discretization parameter \(n\) in an uncontrolled way as \(n\to\infty\). Furthermore, the framework depends on qualitative arguments that automatically fail as \(n\to\infty\), for example that “the numerical trajectory remains on a compact subset of \(M\)”. These hurdles for rigorous backward error analysis of Hamiltonian PDEs have persisted since the birth of symplectic numerical integration. So far, the only rigorous results are for semi-linear PDEs, such as the nonlinear Schrödinger equation, and are based on splitting [4]–[6].

In this paper, we give rigorous results on backward error analysis and near conservation of energy for full discretizations of the two-dimensional (2-D) Euler equations—a fully nonlinear Hamiltonian PDE. The spatial discretization is based on matrix hydrodynamics (cf. [7]), leading to an isospectral flow of matrices which we discretize in time using isospectral symplectic Runge–Kutta (ISOSYRK) methods [8]. The core theoretical tool behind our results is a new theory of Butcher series over biplanar forests, applicable to isospectral flows. The theory is obtained by lifting discrete Lie–Poisson reduction theory (cf. [8]) from individual vector fields and Hamiltonian functions to the level of Butcher series. Interestingly, the resulting series for ISOSYRK methods, via biplanar forests, are distinct from Butcher series for Runge–Kutta–Munthe-Kaas (RKMK) methods. Although the focus in this paper is on the 2-D Euler equations, the approach works also for other Hamiltonian PDEs admitting matrix discretizations as isospectral flows.

We continue the introduction with a brief recollection of the Hamiltonian (Lie–Poisson) structure of the 2-D Euler equations, spatial discretizations via matrix hydrodynamics, the ISOSYRK methods, and a statement of the main result. Thereafter, we continue in 2 with a review of Butcher series and modified vector fields for symplectic Runge–Kutta methods for canonical Hamiltonian systems. This theory is the basis for Butcher series and backward error analysis over biplanar forests developed in 3 and 4, providing a proof of the main theorem. Finally, in 5, we give numerical verifications of the theoretical results.

1.1 The 2-D Euler equations and matrix hydrodynamics↩︎

For the incompressible Euler equations on the sphere \(\mathbb{S}^2\), if \(v = v(x,t)\) is the velocity field for the fluid and \(\omega = \operatorname{curl}v\) is the corresponding vorticity function, the equations are \[\label{eq:euler2D} \dot{\omega} + \{\omega, \psi \} = 0, \quad -\Delta\psi = \omega,\quad \omega(\cdot,0) = \omega_0\in C^\infty(\mathbb{S}^2),\tag{2}\] where \(\dot{\omega} = \frac{\partial\omega}{\partial t}\) is the time derivative, \(\{\cdot,\cdot\}\) is the Poisson bracket on \(\mathbb{S}^2\), \(\Delta\) is the Laplace-Beltrami operator, and \(\psi\) is the stream function. Euler’s equations 2 constitute a Lie–Poisson system for the infinite-dimensional Lie algebra of divergence free vector fields on \(\mathbb{S}^2\) (see e.g. [7]). The Hamiltonian for the system is \[H(\omega) = \frac{1}{2}\int_{S^2} \omega (-\Delta)^{-1}\omega \, \mathrm{d}x .\]

The core idea of matrix hydrodynamics, initiated by Zeitlin [9], [10], is to obtain spatial discretizations for Lie–Poisson PDEs via quantization theory: replace the infinite-dimensional Poisson algebra of smooth functions \((C^\infty(\mathbb{S}^2), \{ \cdot,\cdot\})\) by the finite-dimensional Lie algebra of skew-Hermitian matrices \(\mathfrak{u}(n)\). This can be achieved via the Berezin–Toeplitz quantization operator \(\mathcal{T}_n\colon C^\infty(\mathbb{S}^2) \to \mathfrak{u}(n)\) (cf. [11]). For the Euler equations 2 , it leads to the Euler–Zeitlin equations, namely the isospectral matrix flow \[\label{eq:euler-zeitlin} \dot{W} + \frac{1}{\hbar_n}[W, S] = 0, \quad -\Delta_n S = W, \quad W(0) = W_0\in \mathfrak{u}(n),\tag{3}\] where \(W=W(t)\) and \(S = S(t)\) are skew-Hermitian matrices, \(\hbar_n = 2/\sqrt{n^2-1}\) is a scaling constant, \([\cdot,\cdot]\) is the commutator, and \(\Delta_n\) is the Hoppe–Yau operator (see 6). The system 3 constitute itself a Lie–Poisson system for the Lie algebra \(\mathfrak{u}(n)\) with the scaled bracket \(\frac{1}{\hbar_n}[\cdot,\cdot]\) and the Hamiltonian \[\label{eq:zeitlin95hamiltonian} H_n(W) = \frac{2\pi}{n}\mathop{\mathrm{Tr}}(W^\dagger (-\Delta_n)^{-1} W).\tag{4}\] For \(\omega\in C^\infty(S^2)\) we have \(H_n\big(\mathcal{T}_n(\omega)\big) \to H(\omega)\) as \(n\to \infty\). That the flow 3 is isospectral means geometrically that it preserves the coadjoint orbits (or symplectic leaves) of the Lie–Poisson structure. For a full account of matrix hydrodynamics, including convergence results as \(n\to\infty\), we refer to the work of Modin and Viviani [7] and references therein.

1.2 Isospectral symplectic Runge–Kutta (ISOSYRK) methods↩︎

Consider isospectral flows of the form \[\label{eq:isospectral} \dot{W} = [f(W), W] \,, \quad W(0) = W_0\in \mathfrak{u}(n),\tag{5}\] for a matrix \(W = W(t)\) and a linear function \(f\colon \mathfrak{u}(n)\to\mathfrak{u}(n)\).³ When \(f(W) = \frac{\delta H}{\delta W}\) for some function \(H\), the flow 5 constitutes a Hamiltonian system on \(\mathfrak{u}(n)^*\simeq \mathfrak{u}(n)\) for the Lie–Poisson structure \[\prec F,G\succ(W) = \Big\langle W, \Big[\frac{\delta F}{\delta W}, \frac{\delta G}{\delta W}\Big]\Big\rangle, \quad F,G\in C^\infty(\mathfrak{u}(n)^*).\] The Euler–Zeitlin equations 3 corresponds to the case \(f(W) = \frac{\delta H}{\delta W} = -\frac{1}{\hbar}\Delta_n^{-1}W\) for the quadratic Hamiltonian 4 scaled by \(n/(2\pi \hbar)\).

Given a Butcher tableau \(A = (a_{ij})_{i,j=1}^s\) and \(b = (b_i)_{i=1}^s\) of a symplectic Runge–Kutta method, the corresponding ISOSYRK method \(W_k \mapsto W_{k+1}\) for the flow 5 is defined by \[\label{eq:isosyrk} \begin{align} W_{k+1} &= W_k + h \sum_{i=1}^s b_i [f(\tilde{W}_i), \tilde{W}_i] \,, &\tilde{W}_i &= W_k + h \sum_{j=1}^s a_{ij} (X_j + Y_j + K_{ij}) \,, \\ X_i &= - (W_n + h \sum_{j=1}^s a_{ij} X_j) f(\tilde{W}_i) \,, &Y_i &= f(\tilde{W}_i) (W_k + h \sum_{j=1}^s a_{ij} Y_j) \,, \\ K_{ij} &= h \sum_{k=1}^s f(\tilde{W}_i) (a_{ik} X_k + a_{jk} K_{ik}) \,, & & \end{align}\tag{6}\] where \(h>0\) is the timestep and \(\tilde{W}_{i}, X_i, Y_i, K_{ij}\) are intermediate variables.

The proof of this results uses Lie–Poisson reduction theory. A brief summary goes as follows. Via the momentum map \[\mu\colon T^*\mathrm{U}(n) \to \mathfrak{u}(n)^*, \quad \mu(Q,P) = Q^\dagger P ,\] the Hamiltonian \(H\) on \(\mathfrak{u}(n)^*\) is lifted to a Hamiltonian \(\mathcal{H} = H\circ\mu\) in the variables \((Q,P)\). In these variables, the Lie–Poisson system 5 becomes a canonical Hamiltonian system on \(T^*\mathrm{U}(n)\) \[\label{eq:canonical95equations} \dot{Q} = \frac{\delta \mathcal{H}}{\delta P} = Q f(Q^\dagger P)^\dagger ,\quad \dot{P} = -\frac{\delta \mathcal{H}}{\delta Q} = -P f(Q^\dagger P), \quad Q(0) = Q_0, \quad P(0) = P_0.\tag{7}\] Its solutions \((Q(t),P(t))\) thus map to solutions \(W(t) = \mu(Q(t),P(t))\) of equation 5 . Similarly, via the momentum map the ISOSYRK method 6 corresponds to the underlying symplectic Runge–Kutta method applied to the canonical Hamiltonian system 7 .

1.3 Main result↩︎

Given that ISOSYRK methods are symplectic, it is natural to ask if the energy is nearly conserved for long times \(t_k\gg 1\) while at the same time \(n\to\infty\). As explained in the introduction, this is a question about the modified Hamiltonian and backward error analysis. Our main result is that arbitrary truncations of the modified Hamiltonian for the Euler–Zeitlin equations are well-behaved as \(n\to\infty\) as long as the time step is scaled as \(h = \mathcal{O}(\epsilon/n)\).

2 Butcher series and backward error analysis↩︎

The theory of Butcher series was introduced in the 1960’s by Butcher [12] and later developed by Hairer and Wanner [13] in the 1970’s. Consider an ordinary differential equation in \(\mathbb{R}^d\) of the following form \[y^\prime = f(y), \quad f : \mathbb{R}^d \to \mathbb{R}^d, \quad y(0) = y_0.\] The solution \(y(h)\) as well as the approximation \(y_1\) computed using a Runge–Kutta method can be Taylor expanded around \(h=0\). It was noted by Cayley [14] in 1857 that the terms (also called elementary differntials) appearing in such Taylor expansions can be represented by non-planar trees. Let \(T\) be the set of non-planar rooted trees as defined in 1.

All non-planar trees up to size \(4\) can be found below with all edges being directed downwards, \[\tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,.\] We note that since the trees are non-planar, the order of branches does not matter and we have \[\tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }= \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }.\]

Let \(F\) be the set of non-planar forests defined in 2.

Note that the empty graph \(\mathbf{1}\) is a forest. Let \(B^+ : F \to T\) be a surjective map that adds a new vertex to a forest and connects all roots of the forest to the new vertex thus creating a tree whose root is the new vertex. For example, \[B^+(\mathbf{1}) = \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; }\,, \quad B^+ ( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad B^+( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; }) = \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,.\]

The correspondence between non-planar trees and elementary differentials is denoted by a map \(\mathbb{F}_f\) defined in 3

Some examples of elementary differentials are, \[\mathbb{F}_f( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; }) = f \,, \quad \mathbb{F}_f( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = \sum_{i,j,k=1}^d f^i f^k \frac{\partial f^j}{\partial x_i} \frac{\partial^2 f}{\partial x_j \partial x_k} = f^{(2)} (f, f^\prime f).\]

Let \(\sigma : T \to \mathbb{R}\) denote the symmetry of a tree defined as the number of automorphisms of the tree, for example, \[\sigma( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; }) = 1 \,, \quad \sigma( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = 2 \,, \quad \sigma( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = 4 \,, \quad \sigma ( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = 3! \,.\] Given a tree \(\tau = B^+(\tau_1^{r_1} \cdots \tau_n^{r_n})\) where all \(\tau_i\) are distinct with multiplicities \(r_i\) for \(i=1,\dots,n\), we have the following recursive formula, \[\sigma(\tau) = r_1! \cdots r_n! \sigma(\tau_1)^{r_1} \cdots \sigma(\tau_n)^{r_n} \,.\] Let \(|\tau|\) denote the number of vertices of a tree \(\tau\), for example \(| \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }| = 5\).

The Taylor expansions of the exact solution and of Runge–Kutta methods \(\Phi_h\) are written using Butcher series as \[y(h) = y_0 + B_f(1/\gamma)(y_0), \quad \Phi_h (y_0) = y_0 + B_f(a) (y_0),\] with the coefficient map \(a : T \to \mathbb{R}\) being defined in [prop:coeff95a95np] and the map \(\gamma\) being the factorial of a tree \(\tau = B^+(\tau_1 \cdots \tau_n)\) defined as \[\gamma(\tau) = |\tau| \cdot \prod_{i=1}^n \gamma(\tau_i) \,, \quad \gamma( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; }) = 1 \,.\]

2.1 Backward error analysis of a symplectic integrator↩︎

We present backward error analysis by following [15] and use the formalism of rooted trees to compute the modified equation and the modified Hamiltonian. Assume an integrator \(\Phi_h\) can be written using Butcher series as \[\Phi_h (y) = y + B_f (a)(y) \,.\] Let \(\tau \to\gamma\) denote the Butcher product of trees \(\tau, \gamma \in T\) which attaches the root of \(\tau\) to the root of \(\gamma\), for example, \[\tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\to \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }= \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,.\] Integrator \(\Phi_h\) is symplectic if the coefficient map \(a\) satisfies the following property \[a(\tau \to\gamma) + a(\gamma \to\tau) = a(\tau) a(\gamma) \quad \text{for } \tau, \gamma \in T \,.\] Integrator \(\Phi_h\) can be viewed as an exact flow of a modified vector field \(\tilde{f}_h\), that is, \[\Phi_h(y) = y + B_f(a)(y) = y + B_{\tilde{f}_h} (1/\gamma)(y) \,.\]

Property ?? is a consequence of the symplecticity of the integrator and is equivalent to the fact that the modified vector field \(\tilde{f}_h\) is Hamiltonian [15]. It is used to introduce an equivalence relation on \(T\), \[\tau \to\gamma \; \sim \; \gamma \to\tau \,, \quad \text{for all } \tau, \gamma \in T \,,\] which, together with the property \(b(\tau \to\tau) = 0\), leads to the introduction of non-superfluous free trees which are equivalence classes of rooted trees with respect to the equivalence relation \(\sim\) excluding equivalence classes which contain \(\gamma = \tau \to\tau\) for some \(\tau \in T\). The set of non-superfluous free trees is denoted by \(FT^\prime\).

We set the canonical representative of every equivalence class to be the maximal element of the class with respect to a total order of rooted trees defined in [18], [19]. Let the canonical representative of \(\hat{\tau} \in FT^\prime\) be denoted by \(\hat{\tau}_* \in T\).

The values of the coefficient map \(b\) are completely determined by the values of \(b\) on the canonical representatives of non-superfluous free trees which greatly reduces the number of computations that need to be performed to compute the values of \(b\).

Canonical representatives in \(FT^\prime\) up to size \(5\) are \[\tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \dots.\]

The values of the coefficient map \(b\) given in 2 for canonical representatives of order up to \(5\) are given in 1.

We recall that the modified equation of a symplectic integrator applied to a Hamiltonian system is a modified Hamiltonian system. That is, there exists a modified Hamiltonian \(\tilde{H}_h\) such that \(\tilde{f}_h = J^{-1} \nabla \tilde{H}_h\). Let us now give an explicit formula for \(\tilde{H}_h\) as it is derived in [15].

Elementary Hamiltonians are used to represent the modified Hamiltonian as a formal sum over non-superfluous non-rooted trees as follows.

2.1.1 Analytic bounds↩︎

We note that the series defining the modified Hamiltonian \(\tilde{H}_h\) is a formal power series in \(h\) and its convergence is not guaranteed. As shown in [15] (see also [3]), there exists an optimally truncated modified Hamiltonian \(\tilde{H}_h^N\) such that the numerical trajectory generated by the integrator \(\Phi_h\) preserves \(\tilde{H}^N_h\) over exponentially long times. More precisely, there exists a constant \(c > 0\) such that \[\tilde{H}_h^N (y_n) = \tilde{H}_h^N (y_0) + \mathcal{O}(e^{-c/h}) \,, \quad \text{for } nh \leq e^{c/h} \,.\] This result is obtained by assuming that the trajectory stays on a compact domain. This assumption does not hold in our setting as is discussed in the introduction and, therefore, we obtain this result with a different set of assumptions in 4.2. The key ingredients for our derivation are the bounds of the number of rooted trees of size \(j \in \mathbb{N}\) as well as the bounds of \(|b(\tau)|\) for \(\tau \in T\). We derive these bounds below.

Let us consider the product \(a \curvearrowright b\) of coefficient maps \(a, b : T \to \mathbb{R}\) defined as, \[(a \curvearrowright b) (\tau) = \sum_{e \in E(\tau)} a(\tau \setminus_1 e) b(\tau \setminus_2 e) \,,\] where \(E(\tau)\) is the set of edges of \(\tau\), \(\tau \setminus e\) is the pair of trees obtained by cutting the edge \(e\) of \(\tau\), \(\tau \setminus_2 e\) is the tree which contains the root of \(\tau\), and \(\tau \setminus_1 e\) is the tree which does not contain the root of \(\tau\). Coefficient map \(a^{\curvearrowright m}\) is defined as \[\begin{align} a^{\curvearrowright m} (\tau) &= \big( a \curvearrowright(a \curvearrowright\cdots (a \curvearrowright a) \cdots ) \big) (\tau) \\ &= \sum_{\substack{e_i \in E(\tau)\\ i = 1, \dots, m}} a(\tau \setminus_1 e_1) \cdots a(\tau \setminus_1 e_m) a\big(\tau \setminus_2 (e_1, \dots, e_m)\big) \,, \end{align}\] where \(\tau \setminus_2 (e_1, \dots, e_m) = ( \cdots (\tau \setminus_2 e_1) \cdots ) \setminus_2 e_m\) is the tree obtained by cutting the edges \(e_1, \dots, e_m\) of \(\tau\) and keeping only the part which contains the root of \(\tau\). We note that \(a^{\curvearrowright m} (\tau) = 0\) if \(m \geq |\tau|\) since there are only \(|\tau| - 1\) edges in \(\tau\). We denote \(a^{\curvearrowright 0} = 0\).

Let \(b : T \to \mathbb{R}\) denote the coefficient map for the modified equation of the symplectic Runge–Kutta method with coefficient map \(a : T \to \mathbb{R}\). Following the discussion in [15] (see also [17]), we have, \[a = \exp^\curvearrowright(b) = \sum_{m=1}^\infty \frac{1}{m!} b^{\curvearrowright m} \,,\] analogously to \(1/\gamma = \exp^\curvearrowright(\delta_\bullet)\) for the exact solution where \(\delta_\bullet ( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; }) = 1\) and \(\delta_\bullet (\tau) = 0\) for \(\tau \neq \bullet\).

3 Butcher series over biplanar forests↩︎

We introduce the algebraic tools necessary for the description of Lie–Poisson reduction of Butcher series which leads to the introduction of biplanar forests and forest momentum map. Lie–Poisson reduction of Butcher series is a key ingredient in the backward error analysis of ISOSRYK methods.

Let \(\mathcal{U}\) be the universal enveloping algebra of the Lie algebra \(\mathfrak{u}(n)\). Let us consider the commutative tensor product \(\mathcal{U}^2 = \mathcal{U}\cdot \mathcal{U}\) where \((A,B) = (B,A)\) for \((A,B) \in \mathcal{U}^2\). We endow \(\mathcal{U}^2\) with the product, \[(A, B) \cdot (C, D) = (AC, BD) + (AD, BC) \,, \quad \text{for } A, B, C, D \in \mathcal{U}\,.\] Let us introduce the following action \(\odot : \mathcal{U}^2 \otimes \mathfrak{u}(n) \to \mathfrak{u}(n)\). Let \(A, B \in \mathcal{U}\) and \(W \in \mathfrak{u}(n)\), then, \[(A, B) \odot W := AWB^\dagger + BWA^\dagger \,.\] We define \(A \odot W\) as \((A,I) \odot W\). This action arises naturally when considering directional derivatives of the momentum map \(\mu : T^*\text{U}(n) \to \mathfrak{u}(n)^* \cong \mathfrak{u}(n)\), \(\mu(Q,P) = Q^\dagger P\) as is done in 3.3.

In the rest of this section, we introduce Butcher series over biplanar forests which are constructed using the action of \(\mathcal{U}^2\) on \(\mathfrak{u}(n)\). We introduce the forest momentum map and use it to obtain the Butcher series of isospectral symplectic Runge–Kutta (ISOSYRK) methods from the corresponding Butcher series of the symplectic Runge–Kutta methods through the Lie–Poisson reduction of Butcher series. We use the obtained biplanar Butcher series to compare ISOSYRK methods to the well-known Runge–Kutta–Munthe-Kaas methods for Lie–Poisson systems.

3.1 Biplanar forests↩︎

Let us introduce biplanar trees and forests which are used to denote the terms appearing in the Taylor expansion of ISOSYRK methods as well as the exact solution of an isospectral differential equation 5 . In particular, biplanar forests represent elements of \(\mathcal{U}^2\) through the map \(\mathbb{F}_f\) defined in 7.

Let \(BT\) and \(BF\) denote the sets of biplanar trees and forests defined in 6. We draw biplanar forest \((\pi, \eta)\) by separating \(\pi\) and \(\eta\) by \(\times\). If \(\pi = \mathbf{1}\) or \(\eta = \mathbf{1}\), then \(\times\) is omitted. Some examples of biplanar forests can be found below, \[\mathbf{1}, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {\times} edge from parent[ -, solid, solid, draw=none ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.5, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.5, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.5, 0.0) {} ; }\,, \quad \text{where } \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.5, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.5, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.5, 0.0) {} ; }\neq \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.5, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.5, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.5, 0.0) {} ; }\,, \quad \text{but } \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.5, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.5, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.5, 0.0) {} ; }= \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.5, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,.\] The vector spaces spanned by biplanar forests and trees are denoted by \(\mathcal{BF}\) and \(\mathcal{BT}\), respectively. Let \(T(\mathcal{BT})\) denote the tensor algebra of biplanar trees, then, \(\mathcal{BF}\) is the commutative tensor product of \(T(\mathcal{BT})\) with itself, that is, \(\mathcal{BF}= T(\mathcal{BT})^2 = T(\mathcal{BT}) \cdot T(\mathcal{BT})\) and \(\mathcal{BT}= B^+ (\mathcal{BF})\) as vector spaces. Let us abuse the notation and denote the set of ordered monomials of biplanar trees by \(T(BT)\).

For example, \[\begin{align} \mathbb{F}_f( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; }) \odot W = &\mathbb{F}_f( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) (W) \cdot W \cdot \mathbb{F}_f( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; })^\dagger + \mathbb{F}_f( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; }) (W) \cdot W \cdot \mathbb{F}_f( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; })^\dagger \\ = &f(f(W) \cdot W) \cdot W \cdot f(W)^\dagger + f(W \cdot f(W)^\dagger) \cdot W \cdot f(W)^\dagger \\ &+ f(W) \cdot W \cdot f(f(W) \cdot W)^\dagger + f(W) \cdot W \cdot f(W \cdot f(W)^\dagger)^\dagger \,. \end{align}\]

Let the symmetry \(\sigma(\pi, \eta)\) of a biplanar forest \((\pi, \eta) \in BF\) be defined as, \[\sigma(\pi, \eta) = \begin{cases} 2 \sigma(\pi)^2 \,, \quad \text{if } \pi = \eta \,, \\ \sigma(\pi) \sigma(\eta) \,, \quad \text{otherwise,} \end{cases}\] and \(\sigma(\pi \cdot \eta) = \sigma(\pi) \sigma(\eta) \,, \sigma \big(B^+(\pi,\eta)\big) = \sigma(\pi, \eta) \,, \sigma(\mathbf{1}) = 1\). Some examples of the values of the symmetry \(\sigma(\pi, \eta)\) for \((\pi, \eta) \in BF\) can be found below, \[\sigma( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {\times} edge from parent[ -, solid, solid, draw=none ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; }) = 2 \,, \quad \sigma( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {\times} edge from parent[ -, solid, solid, draw=none ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.5, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {\times} edge from parent[ -, solid, solid, draw=none ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = 8 \,, \quad \sigma( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = 2 \,, \quad \sigma( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {\times} edge from parent[ -, solid, solid, draw=none ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; }) = 1 \,.\]

We note that we denote the coefficient maps over classical trees and forests by latin letters, for example, \(a : F \to \mathbb{R}\,, \; b : T \to \mathbb{R}\), and the coefficient maps over biplanar trees and forests by greek letters, for example, \(\alpha : BF \to \mathbb{R}\,, \; \beta : T(BT) \to \mathbb{R}\). This also allows us to distinguish between classical Butcher series and biplanar Butcher series.

3.1.1 Coadjoint coefficient maps↩︎

In this section, we introduce coadjoint coeficient maps and infinitisimal coadjoint coefficient maps which correspond to biplanar Butcher series with nice geometric properties.

Recall that the tensor algebra \(T(\mathcal{BT})\) is a Hopf algebra with the deconcatenation coproduct \(\Delta_\cdot : T(\mathcal{BT}) \to T(\mathcal{BT}) \otimes T(\mathcal{BT})\) defined, for \(\tau_i \in BT\), as, \[\Delta_\cdot (\tau_1 \cdots \tau_n) = \mathbf{1}\otimes \tau_1 \cdots \tau_n + \sum_{i=1}^{n-1} (\tau_1 \cdots \tau_i) \otimes (\tau_{i+1} \cdots \tau_n) + \tau_1 \cdots \tau_n \otimes \mathbf{1}\,,\] and antipode \(S : T(\mathcal{BT}) \to T(\mathcal{BT})\) defined as, \[S(\tau_1 \cdots \tau_n) = (-1)^{n} \tau_n \cdots \tau_1 \,.\] For example, \(\Delta_\cdot (\tau_1 \tau_2 \tau_3) = \mathbf{1}\otimes \tau_1 \tau_2 \tau_3 + \tau_1 \otimes \tau_2 \tau_3 + \tau_1 \tau_2 \otimes \tau_3 + \tau_1 \tau_2 \tau_3 \otimes \mathbf{1}\) and \(S(\tau_1 \tau_2 \tau_3) = - \tau_3 \tau_2 \tau_1\).

Let \(\pi \in T(\mathcal{BT})\), we note the following property of the elementary differential \(\mathbb{F}_f(\pi)\), \[\label{eq:F95antipode} \mathbb{F}_f (\pi)^\dagger = \mathbb{F}_f \big(S(\pi)\big) \,.\tag{8}\]

Let \(\alpha : BF \to \mathbb{R}\) be a coefficient map over biplanar forests, then, we denote by \(\alpha_T\) the restriction of \(\alpha\) to \(T(BT) \cong T(BT) \cdot \mathbf{1}\subset BF\), that is, \(\alpha_T := \alpha \raisebox{-.5ex}{|}_{T(BT)}\). Given two coefficient maps \(\alpha, \tilde{\alpha} : BF \to \mathbb{R}\), we define \(\alpha_T \cdot \tilde{\alpha}_T : T(BT) \to \mathbb{R}\) as, \[(\alpha_T \cdot \tilde{\alpha}_T)(\pi) = \sum_{(\pi)} \alpha_T (\pi_{(1)}) \tilde{\alpha}_T(\pi_{(2)}) \,,\] where \(\Delta_\cdot (\pi) = \sum_{(\pi)} \pi_{(1)} \otimes \pi_{(2)}\) is the deconcatenation coproduct of \(\pi\) written using Sweedler’s notation. Denote by \(\delta_\mathbf{1}: T(BT) \to \mathbb{R}\) the coefficient map \(\delta_\mathbf{1}(\mathbf{1}) = 1\) and \(\delta_\mathbf{1}(\pi) = 0\) for \(\pi \neq \mathbf{1}\).

3.2 Forest momentum map↩︎

Let us introduce the forest momentum map which plays a central role in the Lie–Poisson reduction of Butcher series. Let \(T^2\) denote the set of non-planar forests with at most two trees.

For example, \[\psi( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; }) = \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; }\,, \quad \psi( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; }) = \psi( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \psi( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = \psi( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {\times} edge from parent[ -, solid, solid, draw=none ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,,\] \[\psi( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; }) = \psi( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; }) = \psi( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = \psi( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) = \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,.\]

The map \(\psi\) is related to the Butcher product \(\to: \mathcal{T}\otimes \mathcal{T}\to \mathcal{T}\) in the following way, \[\psi \big(\pi \cdot \tau \big) = \psi(\pi) \to\psi(\tau) \,,\] where \(\pi \in T(BT), \tau \in BT\). Let the Butcher product be extended to \(\to: \mathcal{T}\cup \{\mathbf{1}\} \otimes \mathcal{T}\to \mathcal{T}\) as, \[\mathbf{1}\to\tau = \tau \,, \quad \text{for } \tau \in \mathcal{T}\,.\]

Let us abuse the notation and denote by \(\langle {-}, {-}\rangle_\sigma\) inner products over non-planar and biplanar forests \(\mathcal{F}\) and \(\mathcal{BF}\) defined as \[\langle \pi, \eta \rangle_\sigma = \begin{cases} \sigma(\pi) \,, \quad \text{if } \pi = \eta \,, \\ 0 \,, \quad \text{otherwise.} \end{cases}\]

Let \(\psi^* : \mathcal{T}^2 \to \mathcal{BF}\) denote the adjoint of the momentum forest map \(\psi\) with respect to the inner products \(\langle {-}, {-}\rangle_\sigma\), that is, \[\langle \psi(\pi, \eta) , \nu \rangle_\sigma = \langle (\pi, \eta), \psi^* (\nu) \rangle_\sigma \,, \quad \text{for } \pi, \eta \in T(BT) \,, \nu \in \mathcal{T}^2 \,.\]

The rest of this section is devoted to the derivation of an explicit formula for \(\psi^*\). Let us start by introducing the necessary algebraic ingredients. Let \(\Delta_\shuffle : \mathcal{F}\to \mathcal{F}\otimes \mathcal{F}\) denote the deshuffle coproduct defined as \[\Delta_\shuffle (\tau) = \mathbf{1}\otimes \tau + \tau \otimes \mathbf{1}\,, \quad \Delta_\shuffle (\pi \cdot \eta) = \Delta_\shuffle (\pi) \cdot \Delta_\shuffle (\eta) \,,\] where \((a \otimes b) \cdot (c \otimes d) = (a \cdot c) \otimes (b \cdot d)\). For example, \[\begin{align} \Delta_\shuffle (\tau_1 \tau_2 \tau_3) &= \tau_1 \tau_2 \tau_3 \otimes \mathbf{1}+ \tau_1 \tau_2 \otimes \tau_3 + \tau_1 \tau_3 \otimes \tau_2 + \tau_2 \tau_3 \otimes \tau_1 \\ &\quad + \tau_1 \otimes \tau_2 \tau_3 + \tau_2 \otimes \tau_1 \tau_3 + \tau_3 \otimes \tau_1 \tau_2 + \mathbf{1}\otimes \tau_1 \tau_2 \tau_3 \,. \end{align}\] Let \(\Delta_\to: \mathcal{T}\to \mathcal{T}\cup \{\mathbf{1}\} \otimes \mathcal{T}\) denote the map which, for \(\tau = B^+(\tau_1 \cdots \tau_n)\), is defined as \[\Delta_\to(\tau) = \mathbf{1}\otimes \tau + \sum_{i=1}^n \tau_i \otimes B^+(\tau_1 \cdots \tau_{i-1} \tau_{i+1} \cdots \tau_n) \,.\]

We prove some properties of \(\psi^*\) in 3 that are used in [prop:psi4295formula] to derive a formula for \(\psi^*\).

3.3 Lie–Poisson reduction of Butcher series↩︎

In this section, we obtain the biplanar Butcher series of ISOSYRK methods from the Butcher series of the corresponding symplectic Runge–Kutta methods by performing Lie–Poisson reduction of Butcher series.

Let \(V \in \mathfrak{u}(n)\) be an element of the Lie algebra, then \(V\) acts on \((Q,P) \in T^* \text{U}(n)\) as \[V \cdot (Q, P) = (Q V^\dagger, -P V) \,.\] The lift of a vector field \(f : \mathfrak{u}(n) \to \mathfrak{u}(n)\) to \(T^* \text{U}(n)\) is then given by the vector field, \[(\mu_\sharp f) (Q, P) := f(\mu(Q,P)) \cdot (Q,P) = \Big( Q f(Q^\dagger P)^\dagger, -P f(Q^\dagger P) \Big) \,.\] We extend the action of \(\mathfrak{u}(n)\) on \(T^* \text{U}(n)\) to an action of the universal enveloping algebra \(\mathcal{U}\) on \(T^* \text{U}(n)\) in the following way, \[\begin{align} \mathbf{1}\cdot (Q, P) &= (Q, P) \,, \\ (V A) \cdot (Q, P) &= V \cdot \big(A \cdot (Q,P) \big) \,. \end{align}\] where \(V \in \mathfrak{u}(n), A \in \mathcal{U}, (Q,P) \in T^* \text{U}(n)\) which can be written explicitly as \[A \cdot (Q, P) = \big(Q A^\dagger, P S(A)\big) \,,\] where \(S\) is the antipode of the universal enveloping algebra \(\mathcal{U}\), that is, \[S(V_1 \cdots V_n) = (-1)^n V_n \cdots V_1 \,,\] for \(V_i \in \mathfrak{u}(n)\). We use this extension to define, for \(\pi \in T(BT)\), \[\mu_\sharp\mathbb{F}_f (\pi) (Q, P) = \Big( Q \cdot \mathbb{F}_f(\pi)^\dagger (W), P \cdot \mathbb{F}_f\big(S(\pi)\big) (W) \Big) \,,\] where \(W = Q^\dagger P\). We note that \(\mathbb{F}_{\mu_\sharp f} (\tau) : T^*\text{U}(n) \to TT^*\text{U}(n)\) for \(\tau \in T\) is an elementary differential defined on \(T^* \text{U}(n)\), while \(\mathbb{F}_f (\pi) : \mathfrak{u}(n) \to \mathcal{U}\) for \(\pi \in T(BT)\) is an elementary differential defined on \(\mathfrak{u}(n)\).

Let \(k \in \{0, 1, 2\}\), then the directional derivative of the momentum map \(\mu : T^* \text{U}(n) \to \mathfrak{u}(n)\) is denoted by \[\mathbb{F}_{\mu_\sharp f} (\pi) [\mu] := \mu^{(k)} \Big( \mathbb{F}_{\mu_\sharp f} (\tau_1), \dots, \mathbb{F}_{\mu_\sharp f} (\tau_k) \Big) \,, \quad \text{for } \pi = \tau_1 \cdots \tau_k \,, \tau_i \in T \,.\] Let \(T(BT)_* := T(BT) \setminus \{\mathbf{1}\}\) be the set of non-empty monomials \(\tau_1 \dots \tau_n\) with \(\tau_i \in BT\).

Given a coefficient map \(a : T \to \mathbb{R}\), the integrator with the Taylor expansion given by \(y_0 + B_{\mu_\sharp f} (a)(y_0)\) is denoted by \(\Phi_{\mu_\sharp f}^{(a)}\). Similarly, given a coefficient map \(\alpha : BF \to \mathbb{R}\), the integrator with the Taylor expansion given by \(B_f (\alpha) \odot W_0\) is denoted by \(\Phi_f^{(\alpha)}\).

4 shows that the biplanar Butcher series with coefficient map \(a \circ \psi\) is related to the Butcher series with coefficient map \(a\) via the momentum map \(\mu : T^* \text{U}(n) \to \mathfrak{u}(n)\), and the composition and substitution laws of the biplanar Butcher series follow from the composition and substitution laws of the classical Butcher series.

We note that the Butcher product on \(\mathcal{T}\) satisfies the relation, \[\label{eq:butcher95product95relation} \tau \to(\gamma \to\eta) = \gamma \to(\tau \to\eta) \,.\tag{17}\]

We prove that the coefficient map \(a \circ \psi\) of biplanar Butcher series is a coadjoint coefficient map as is defined in [prop:coadjoint95coeff].

3.4 Runge–Kutta–Munthe-Kaas methods for Lie–Poisson systems↩︎

In [21], the authors proposed a class of numerical integrators for Lie–Poisson systems which are obtained by using the Runge–Kutta–Munthe-Kaas (RKMK) framework [22], [23]. In this section, we show that RKMK methods applied to the isospectral equation 5 can be expanded using biplanar Butcher series by introducing a subclass of coadjoint coefficient maps called character coefficient maps.

The solution of the isospectral equation 5 evolves on the coadjoint orbit \(\mathcal{O}_{W_0}\) defined as \[\mathcal{O}_{W_0} = \{ \mathop{\mathrm{Ad}}^*_Q W_0 \; : \; Q \in \text{U}(n) \} \,,\] in particular, the solution can be written as \(W(t) = \mathop{\mathrm{Ad}}^*_{\exp(\omega(t))} W_0\) for some curve \(\omega : \mathbb{R}\to \mathfrak{u}(n)\) with \(\omega(0) = 0\) which solves, \[\dot{\omega}(t) = \mathop{\mathrm{dexp}}^{-1}_{\omega(t)} \Big( f\big(W(t)\big) \Big) \,, \quad \text{where } \mathop{\mathrm{dexp}}^{-1}_u (v) = \sum_{k=0}^\infty \frac{B_k}{k!} \mathop{\mathrm{ad}}^k_u (v) \,,\] with \(\mathop{\mathrm{ad}}^k_u (v) = [u, \mathop{\mathrm{ad}}^{k-1}_u (v)]\), \(\mathop{\mathrm{ad}}_u^0 (v) = v\), and \(B_k\) are the Bernoulli numbers. Therefore, the solution of the isospectral equation can be obtained by solving a differential equation on \(\mathfrak{u}(n)\) and applying the coadjoint action of \(\exp\big(\omega(t)\big)\) to \(W_0\).

This perspective is used in [21] to propose numerical integrators of the form, \[\begin{align} W_{n+1} &= \mathop{\mathrm{Ad}}^*_{\exp(\omega_{n+1})} W_n \,, \\ \omega_{n+1} &= h \sum_{i=1}^s b_i \mathop{\mathrm{dexp}}^{-1}_{\omega_n^i} \big( f(W_n^i) \big) \,, \\ W^i_n &= \mathop{\mathrm{Ad}}^*_{\exp(\omega_n^i)} W_n \,, \\ \omega_n^i &= h \sum_{j=1}^s a_{ij} \mathop{\mathrm{dexp}}^{-1}_{\omega_n^j} \big( f(W_n^j) \big) \,, \end{align}\] where \(b = (b_i)_{i=1}^s\) and \(A = (a_{ij})_{i,j=1}^s\) are coefficients of a Runge–Kutta method of order \(p\). The resulting integrator \(W_n \mapsto W_{n+1}\) is of order \(p\).

Let us consider the shuffle product \(\shuffle : T(\mathcal{BT}) \otimes T(\mathcal{BT}) \to T(\mathcal{BT})\) defined as, \[\tau\pi \shuffle \gamma\eta := \tau (\pi \shuffle \gamma\eta) + \gamma (\tau\pi \shuffle \eta) \,, \quad \text{for } \tau, \gamma \in BT, \;\; \pi, \eta \in T(BT) \,,\] for example, \(\tau_1 \tau_2 \shuffle \tau_3 \tau_4 = \tau_1 \tau_2 \tau_3 \tau_4 + \tau_1 \tau_3 \tau_2 \tau_4 + \tau_1 \tau_3 \tau_4 \tau_2 + \tau_3 \tau_1 \tau_2 \tau_4 + \tau_3 \tau_1 \tau_4 \tau_2 + \tau_3 \tau_4 \tau_1 \tau_2\).

We introduce a subclass of coadjoint coefficient maps, see [prop:coadjoint95coeff], that correspond to RKMK methods.

Note that the converse of [prop:character95coeff] does not hold. Indeed, given a character coefficient map \(\alpha : BF \to \mathbb{R}\), the corresponding method is not necessarily an RKMK method. The reason is that the space \(\mathbb{F}_f \Big(\mathop{\mathrm{Prim}}\big(T(\mathcal{BT})\big)\Big)\) is much larger than the subspace accessible by \(\omega_{n+1}\) and \(\omega^i_n\). In other words, the stages of an RKMK method probe only a restricted family of primitive elements.

4 Backward error analysis over biplanar forests↩︎

In this section, we study the geometric properties of ISOSYRK methods, in particular, we perform backward error analysis and obtain the modified equation as a Butcher series over biplanar forests. We also give an explicit formula for the modified Hamiltonian by introducing fair biplanar trees.

Let us use the substitution law described in 4 to perform the backward error analysis. We recall that the backward error analysis of a Runge–Kutta method consists in finding a modified vector field \(\tilde{f}_h\) such that the exact flow of \(\tilde{f}_h\) coincides with the numerical flow of the Runge–Kutta method. The modified vector field \(\tilde{f}_h\) is expressed as a Butcher series with coefficients given by the substitution law. In particular, we find the coefficient map \(b: T \to \mathbb{R}\) such that, \[B_f (a) = B_{\frac{1}{h} B_f (b)} (1/\gamma) = B_f \big(b \star (1/\gamma)\big) \,.\] Following [26], we can endow the set of coefficient maps together with the substitution law as a product with a group structure, therefore, \(b = a \star (1/\gamma)^{\star-1}\). Using 4 and given a coefficient map \(b\) computed for a symplectic Runge–Kutta method, the coefficient map of the modified vector field of the corresponding ISOSYRK method is given by \(b \circ \psi\), that is, \[B_f (a \circ \psi) = B_{\frac{1}{h} B_f (b \circ \psi)} \big( (1/\gamma) \circ \psi \big) \,.\]

We note that the coefficient map \(b\) satisfies the identity, \[\label{eq:Hamiltonian95vf} b(\tau \to\gamma) + b(\gamma \to\tau) = 0 \,,\tag{20}\] [prop:modified95isospectrality] uses this identity to show that the modified vector field of an ISOSYRK method is isospectral.

We recall 2.1 in which the identity 20 is used to introduce non-rooted trees and the set of non-superfluous non-rooted trees \(FT^\prime\). Non-superfluous non-rooted trees are then used to represent the modified vector field \(\tilde{f}_h\). This approach is used to dramatically decrease the amount of computation necessary to compute the values of the coefficient map \(b : T \to \mathbb{R}\) of the modified vector field as the set of non-superfluous non-rooted trees is much smaller than the set of rooted trees.

We do not represent the biplanar Butcher series of the modified equation using a subset of biplanar forests corresponding to the non-superfluous non-rooted trees since the values of the coefficient map \(b \circ \psi\) are already given by the values of the coefficient map \(b: T \to \mathbb{R}\). However, we introduce such a set in the next section to obtain an explicit biplanar Butcher series representation of the modified Hamiltonian of an ISOSYRK method.

4.1 Modified Hamiltonian↩︎

Recall 3 in which the modified Hamiltonian of a symplectic Runge–Kutta method is expressed as \[\label{eq:modified95Hamiltonian95series} \tilde{\mathcal{H}}_h (Q, P) = \sum_{\hat{\tau} \in FT^\prime} h^{|\hat{\tau}|-1} \frac{b(\hat{\tau}_*)}{\sigma(\hat{\tau}_*)} \mathcal{H}(\hat{\tau}_*)(Q, P) \,,\tag{21}\] with appropriately defined coefficient map \(b: T \to \mathbb{R}\) where elementary Hamiltonian \(\mathcal{H}(\hat{\tau}_*)\) is given by \[\mathcal{H}(\hat{\tau}_*) (Q, P) = \mathbb{F}_{\mu_\sharp f} (\pi) [\mathcal{H}] (Q, P) \,, \quad \text{for } \hat{\tau}_* = B^+(\pi) \,.\]

Similarly to how the directional derivatives of the momentum map \(\mu\) give rise to biplanar forests, the directional derivatives of the Hamiltonian \(\mathcal{H}\) give rise to the space \(\mathcal{BF}^{\leq 2} := \mathcal{BF}\oplus \mathcal{BF}^2\) which is a subspace of the symmetric algebra \(S(\mathcal{BF})\). Let the basis of \(\mathcal{BF}^2\) be denoted by \((\pi_1, \eta_1) \cdot (\pi_2, \eta_2)\) with \((\pi_i, \eta_i) \in BF\). We have, \[(\pi_1, \eta_1) \cdot (\pi_2, \eta_2) = (\pi_2, \eta_2) \cdot (\pi_1, \eta_1) \,, \quad \text{for } (\pi_i, \eta_i) \in BF \,.\] We extend the forest momentum map \(\psi : BF \to T^2\) to \(\psi : BF^2 \to T^4\) by defining for \((\pi_i, \eta_i) \in BF\), \[\psi\big( (\pi_1, \eta_1) \cdot (\pi_2, \eta_2) \big) = \psi(\pi_1, \eta_1) \psi(\pi_2, \eta_2) \,,\] and \(\sigma : BF \to \mathbb{N}\) to \(\sigma : BF^2 \to \mathbb{N}\) by, \[\sigma \big( (\pi_1, \eta_1) \cdot (\pi_2, \eta_2) \big) = \begin{cases} 2 \sigma(\pi_1, \eta_1)^2 \,, \quad \text{if } (\pi_1, \eta_1) = (\pi_2, \eta_2) \,, \\ \sigma(\pi_1, \eta_1) \sigma(\pi_2, \eta_2) \,, \quad \text{otherwise} \,. \end{cases}\] Let us also define \(B^+\) on \(BF^2\) with \(\psi\) commuting with \(B^+\), that is, \[\psi \circ B^+ = B^+ \circ \psi \,, \quad \text{and, moreover,} \quad \sigma \circ B^+ = \sigma \,.\]

The set of fair biplanar trees is denoted by \(FBT\). The name fair is motivated by the fact that the canonical representatives \(\hat{\tau}_*\) of non-rooted trees in \(FT^\prime\) are maximal with respect to a total order [18], [19] according to which a tree is bigger if it has more branches with each branch having a similar number of vertices. This corresponds to the forest \((\pi_1, \eta_1) \cdot (\pi_2, \eta_2) \in BF^{\leq2}\) having \(\pi_1, \pi_2, \eta_1, \eta_2\) with the same number of vertices if possible. A systematic study of fair biplanar trees is left for future work. All non-empty fair biplanar trees up to order \(4\) are listed below, where we separate \((\pi_1, \eta_1)\) and \((\pi_2, \eta_2)\) in \((\pi_1, \eta_1) \cdot (\pi_2, \eta_2)\) by \(\otimes\), \[\tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; }\,,\] \[\tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; }\,,\] \[\tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (5.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (6.0, 0.0) {} ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (5.0, 0.0) {} ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (5.0, 0.0) {} ; }\,.\] \[\tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; }\,, \quad \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; }\,.\]

Let \(\pi \in F\) be a forest and \(P(\pi)\) denote the set of partitions of trees in \(\pi\). Let \(P_{\leq2}(\pi)\) be the subset of \(P(\pi)\) consisting of partitions with at most two parts. The adjoint of \(\psi : BF^{\leq 2} \to T^4\) is denoted by \(\psi^* : T^4 \to BF^{\leq 2}\) and is given by, \[\psi^* (\pi) = \sum_{\substack{p \in P_{\leq2}(\pi) \\ p_i \in T^2}} \psi^*(p_1) \psi^*(p_2) \,,\] where \(\pi \in T^4\) and \(p_2 = \mathbf{1}\) if \(p\) has only one part. We also have \(B^+ \circ \psi^* = \psi^* \circ B^+\). We note that the set \(\{ p \in P_{\leq2} (\pi) \; | \; p_i \in T^2 \}\) has at most \(3\) elements, for example, given \(\pi = \tau_1 \cdots \tau_4\), \[\{ p \in P_{\leq2} (\pi) \; | \; p_i \in T^2 \} = \big\{ \{ \tau_1 \tau_2, \tau_3 \tau_4 \} \,, \{ \tau_1 \tau_3, \tau_2 \tau_4 \} \,, \{ \tau_1 \tau_4, \tau_2 \tau_3 \} \big\} \,.\]

The elementary Hamiltonian \(H\big( B^+(\pi) \big)\) for \(\pi \in BF^{\leq 2}\) is defined as follows, \[\begin{align} H \big(B^+ ( \pi, \eta ) \big) (W) &:= \mathbb{F}_f (\pi, \eta) [H] (W) = H^\prime (W) \big(\mathbb{F}_f(\pi, \eta) \odot W \big) \,, \\ H \Big(B^+ \big( (\pi_1, \eta_1) \cdot (\pi_2, \eta_2) \big) \Big) (W) &:= \mathbb{F}_f \big((\pi_1, \eta_1) \cdot (\pi_2, \eta_2)\big) [H] (W) \\ &= H^{(2)} (W) \big(\mathbb{F}_f(\pi_1, \eta_1) \odot W \,, \;\mathbb{F}_f (\pi_2, \eta_2) \odot W\big) \,. \end{align}\]

This gives us an efficient way to compute the modified Hamiltonian of an ISOSYRK method. Using [prop:isospectral95modified95Hamiltonian], the modified Hamiltonian \(\tilde{H}_h\) can be written as, \[\begin{align} \tilde{H}_h = H &+ h^2 b( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) \mathbb{F}_f \Big( \frac{1}{2} \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; }+ \frac{1}{2} \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; }\Big)[H] + h^3 b( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) \mathbb{F}_f \Big( \frac{1}{2} \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; }\Big) [H] \\ &+ h^4 b( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) \mathbb{F}_f \Big( \frac{1}{2} \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }+ \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; }+ \frac{1}{2} \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; }\\ &\quad \quad \quad \quad \quad \;\; + \frac{1}{2} \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }+ \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; }+ \frac{1}{2} \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; }\Big) [H] \\ &+ h^4 b( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) \mathbb{F}_f \Big( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; }+ \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (5.0, 0.0) {} ; }+ \frac{1}{2} \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.0, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; }+ \frac{1}{2} \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (5.0, 0.0) {} ; }\Big) [H] \\ &+ h^4 b( \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-1.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (-0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (0.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } child { node [ b, label={[label distance=-1mm]0:{\scriptsize{}}} ] at (1.5, 1.0) {} edge from parent[ -, solid, solid, draw=black ] node [!l,right] {\scriptsize{}} } ; }) \mathbb{F}_f \Big(\frac{1}{8} \tikz[planar forest default, planar forest, ] { \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (0.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (1.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (2.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (3.0, 0.0) {\otimes} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (4.0, 0.0) {} ; \node [ .l, label={[label distance=-1mm]0:{\scriptsize{}}}] at (5.0, 0.0) {\times} ; \node [ b, label={[label distance=-1mm]0:{\scriptsize{}}}] at (6.0, 0.0) {} ; }\Big) [H] + \mathcal{O}(h^5) \,. \end{align}\]

See 5 for an explicit computation of the modified Hamiltonian of the ISOMP method applied to the Euler-Zeitlin equation.

4.2 Analytic bounds↩︎

Here we derive analytic bounds concerning the long-time behavior of the numerical flow of ISOSYRK methods. The procedure from [15] assumes that the trajectory stays on a compact domain and uses Cauchy estimates to bound the derivatives of the vector flow of the problem. In our setting, we drop the compactness assumption and replace it by the assumptions \[\label{eq:bound95f} \lVert f(W) \rVert_\infty \leq C_f \rVert W \rVert_\infty \,, \quad \lVert f^\prime(W) \rVert_\infty \leq C_f \,,\tag{23}\] for some constant \(C_f\) independent of \(W\). Additionally, we assume for \(a: T \to \mathbb{R}\) that \[\label{eq:bound95coeff} | a(\tau) | \leq C_a^{|\tau|} \,,\tag{24}\] for some constant \(C_a\) independent of \(\tau\).

First, we need a bound on the number of ordered monomials of biplanar trees \(T(BT)_j\) and biplanar forests \(BF_j\) of size \(j \in \mathbb{N}\).

Next, we expand the ISOSYRK methods \(\Phi_h\) as \[\Phi_h (W) = W + h f(W) + h^2 d_2 (W) + h^3 d_3 (W) + \cdots \,.\]

Let the modified vector field for an ISOSYRK method be written as \[\tilde{f}_h (W) = f(W) + h f_2 (W) + h^2 f_3(W) + \cdots \,.\]

Let us now consider the truncated modified vector field \(\tilde{f}_h^N\), \[\tilde{f}_h^N (W) = f(W) + h f_2 (W) + \cdots + h^{N-1} f_N(W) \,,\] and its flow \(\tilde{\varphi}_h^N\) for \(N \in \mathbb{N}\). We bound \(\|\tilde{f}^N_h(W)\|_\infty\) as follows, \[\begin{align} \|\tilde{f}^N_h(W)\|_\infty &\leq \|f(W)\|_\infty + \sum_{j=2}^N h^{j-1} \|f_j(W)\|_\infty \\ &\leq C_f \|W\|_\infty + \sum_{j=2}^N h^{j-1} \frac{e}{3} (j-1)! \big( 18 C_a C_f \|W\|_\infty \big)^j \,. \end{align}\] Similarly to the discussion in [15], we introduce \[h_0 := \frac{1}{18 e C_a C_f} \,,\] and assume \(N \leq \frac{e h_0}{h\|W\|_\infty}\). We obtain, \[\label{eq:bound95truncated95modeq} \|\tilde{f}^N_h(W)\|_\infty \leq C_f \|W\|_\infty + 6 e C_a C_f \|W\|_\infty \Big( \sum_{j=2}^N \frac{(j-1)!}{N^{j-1}} \Big) \leq C_f \|W\|_\infty \big( 1 + 6 e C_a \big) \,.\tag{25}\]

We now turn to the truncated modified Hamiltonian \[\tilde{H}^N_h (W) := H(W) + h^p H_{p+1} (W) + \cdots + h^{N-1} H_N (W) \,,\] where \(f_j(W) = \nabla H_j(W)^\dagger\). We use the bound 25 of \(\|f_h^N(W)\|_\infty\) as the Lipschitz constant of the truncated Hamiltonian \(\tilde{H}_h^N\) in the following proposition, which is used to prove the conservation of the modified Hamiltonian for exponential time in 6.

6 is an adaptation to our setting of the well-known result that symplectic integrators for finite-dimensional systems conserve the modified Hamiltonian over long times, see, e.g., [15].

5 Numerical experiments↩︎

In this section we carry out numerical simulations with the ISOMP method ?? applied to the Euler–Zeitlin equations 3 . We use the Python package QUFLOW.⁵ For initial data, we select random normally distributed spherical harmonic coefficients \(\omega_{\ell,m}\) for \(\ell \leq 16\). For each matrix size \(n\), this gives rise to an initial data matrix \(W_0 \in \mathfrak{su}(n)\). For a detailed example of how to solve the Euler–Zeitlin equation with QUFLOW, see Modin and Viviani [27]. Our simulation here follows the same setup.

As a first step, let us compute a few terms of the modified Hamiltonian. Let \(H_n(W)\) denote the Hamiltonian 4 and let \(\langle A,B\rangle_2 := \frac{4\pi}{n}\mathop{\mathrm{Tr}}(A^\dagger B)\) so that \(H_n(W) = \frac{1}{2}\langle W,(-\Delta)^{-1}W\rangle_2\). Following the discussion in 4.1 and using the notation, \[S = \mathcal{L}^{-1}_n W := (-\Delta_n)^{-1}W \,, \quad U = S^2 + \mathcal{L}^{-1}_n [S, W] \,, \quad \epsilon = h/\hbar_n\] the modified Hamiltonian \(\tilde{H}_{n,h}(W)\) corresponding to \(H_n(W)\) is given by, \[\begin{align} \tilde{H}_{n,h}(W) = H_n(W) &- \frac{1}{12} \epsilon^2 \langle S W S^\dagger, S \rangle_2 - \frac{1}{24} \epsilon^2 \langle [S, W], \mathcal{L}^{-1}_n [S, W] \rangle_2 \\ &+ \frac{1}{80} \epsilon^4 \langle U W U^\dagger, S \rangle_2 + \frac{1}{160} \epsilon^4 \langle U W + W U^\dagger, \mathcal{L}^{-1}_n (U W) + W U^\dagger \rangle_2 \\ &+ \frac{1}{240} \epsilon^4 \langle U W S^\dagger + S W U^\dagger, \mathcal{L}^{-1}_n[S, W] \rangle_2 \\ &+ \frac{1}{240} \epsilon^4 \langle U W + W U^\dagger, \mathcal{L}^{-1}_n (S W S^\dagger) \rangle_2 \\ &+ \frac{7}{480} \epsilon^4 \langle S W S^\dagger, \mathcal{L}^{-1}_n (S W S^\dagger) \rangle_2 + \mathcal{O}(\epsilon^6) \,. \end{align}\] Notice, as expected, that any truncation \(\tilde{H}_{n,h}^N(W)\) only depends on \(h\) via \(\epsilon\), and on \(n\) via \(\langle\cdot,\cdot\rangle_2\) and \(\mathcal{L}^{-1}_n\). Furthermore, under the quantization \(\mathcal{T}_n\), these operations converge to the corresponding infinite-dimensional operations as \(n\to \infty\) (see [7]).

Carrying out short simulations, with \(n=128\) and \(0 \leq t\leq 2\), for 5 different values of \(\epsilon\), we show in 2 the maximum error in the Hamiltonian \(\lvert H_{n}(W_k)-H_{n}(W_0)\rvert\) an the modified Hamiltonian \(\lvert \tilde{H}^N_{n}(W_k)-\tilde{H}^N_{n}(W_0)\rvert\) truncated to \(N=4\). As expected, the Hamiltonian \(H_n\) is preserved up to the order \(\mathcal{O}(\epsilon^2)\) of the ISOMP method, whereas the modified Hamiltonian \(\tilde{H}_{n,h}\) is preserved up to order \(\mathcal{O}(\epsilon^6)\) (there are no odd \(\epsilon^p\)-terms since the method is symmetric).

In much longer simulations, \(0\leq t\leq 2500\) and fixed \(\epsilon = 0.1\), we show in 3 how the relative error of the Hamiltonian \(H_n\) varies with time for 4 different choices of \(n\). Notice that it fluctuates in all simulations, but the magnitudes of the fluctuations do not grow nor decrease for growing \(n\).

6 Estimates for the Hoppe–Yau operator↩︎

The Hoppe–Yau operator [28] (or Hoppe–Yau Laplacian) is the matrix mapping \[\Delta_n\colon \mathfrak{u}(n) \to \mathfrak{u}(n)\] defined by \[\label{eq:hoppe95yau} \Delta_n P = \frac{1}{\hbar^2}\sum_{\alpha=1}^3 [[P, X_\alpha], X_\alpha],\tag{26}\] where the constant \(\hbar = 2/\sqrt{n^2-1}\) and the matrices \(X_1,X_2,X_3 \in \mathfrak{su}(n)\) are generators for an irreducible, unitary representation of \(\mathfrak{so}(3)\) on \(\mathbb{C}^n\) scaled so that \(\sum_{\alpha=1}^3 X_\alpha^\dagger X_\alpha = I\). In physics, these matrices give a spin-\(s\) representation for \(s = (n-1)/2\).

The Hoppe–Yau operator is an approximation of the Laplacian on the sphere: via Berezin–Toeplitz quantization, the Poisson algebra of smooth functions \[\big(C^\infty(S^2,\mathbb{R}), \{ \cdot,\cdot\}\big)\] is approximated by the finite-dimensional Lie algebra \[\big(\mathfrak{u}(n), \frac{1}{\hbar}[\cdot,\cdot]\big)\] and in this approximation the Hoppe–Yau operator \(\Delta_n\) corresponds to the Laplace–Beltrami operator \(\Delta\) (cf. [28]). From the point-of-view of representation theory, the Hoppe–Yau operator 26 is the Casimir element for the induced representation on \(\mathfrak{u}(n)\) generated by the operators \(\mathrm{ad}_{X_1},\mathrm{ad}_{X_2},\mathrm{ad}_{X_2}\). Thus, the spectrum of \(\Delta_n\) is \(\{-\ell(\ell+1) \}\) for \(\ell = 0,\ldots,n-1\) with eigenspaces \(V_\ell\) of dimension \(2\ell + 1\) (cf. [7]). It is precisely the truncation of the spectrum of \(\Delta\).

Just as the Laplace–Beltrami operator has a 1-dimensional kernel given by constant functions, the Hoppe–Yau operator has a 1-dimensional kernel spanned \(\mathrm{i}\mathbb{R}I\). To make it invertible on \(\mathfrak{u}(n)\) we make the following extension, where we also change the sign \[\label{eq:extededHoppeYau} \mathcal{L}_n\colon \mathfrak{u}(n) \to \mathfrak{u}(n), \qquad \mathcal{L}_nP = \frac{\operatorname{tr}(P)I}{n} - \Delta_n.\tag{27}\] This operator has positive spectrum \(1,2,6,\ldots,n(n-1)\). It is invertible and fulfills \(\operatorname{tr}(\mathcal{L}_n P) = \operatorname{tr}(P)\).

For \(p\in [1,2,\ldots,\infty]\), consider the scaled Schatten \(p\)-norms on \(\mathfrak{u}(n)\) \[\lVert P \rVert_p = \left\{ \begin{align} \left(\frac{1}{n}\sum_{i=1}^n \lvert \lambda_i \rvert^p \right)^{1/p} & \quad\text{if p< \infty} \\ \max_i(\lvert \lambda_1\rvert,\ldots,\lvert \lambda_n\rvert) & \quad\text{if p=\infty}\end{align} \right.\] where \(\lambda_1,\ldots,\lambda_n\) are the eigenvalues of \(P\). These norms correspond, in the limit \(n\to\infty\), to the \(L^p\)-norms on \(C^\infty(S^2,\mathbb{R})\). Since the smallest eigenvalue of \(\mathcal{L}_n\) is \(1\), it follows that \[\lVert \mathcal{L}_n^{-1}W \rVert_2 \leq \lVert W \rVert_2 .\] In other words, the operator norm of \(\mathcal{L}_n^{-1}\) is unitary relative to the \(2\)-norm on \(\mathfrak{u}(n)\). Our main result in this appendix is the analogous result for all \(p\)-norms.

6.1 Proof via positivity↩︎

Recall that a linear operator \(\Phi\colon \mathfrak{gl}(n,\mathbb{C}) \to \mathfrak{gl}(n,\mathbb{C})\) between complex matrices is positive if it maps positive Hermitian matrices to positive Hermitian matrices. For details on positive operators, we refer to the monograph by Bhatia [29].

The operator 27 naturally extends from \(\mathfrak{u}(n)\) to an operator on all complex matrices \(\mathfrak{gl}(n,\mathbb{C})\). Indeed, this extension corresponds to the complexification \(\mathfrak{u}(n)\otimes\mathbb{C} \simeq \mathfrak{gl}(n,\mathbb{C})\).

For the proof, we use the following result.

References↩︎

---

1. Department of Mathematical Sciences, Chalmers University of Technology and University of Gothenburg, Sweden (, )↩︎

2. Submitted to the editors 2026-07-14↩︎

3. For simplicity, we work here with the Lie algebra \(\mathfrak{u}(n)\), but the ISOSYRK methods work more generally for a subset \(V\subset \mathfrak{gl}(n,\mathbb{C})\) and \(f\colon V \to \{ \xi \in V \, | \, [\xi, V] \subset V \}\) (cf. [8]).↩︎

4. \(C_a = 1\) for the ISOMP method ?? .↩︎

5. Available at github.com/klasmodin/quflow.↩︎