1 Introduction
Let $f\colon \mathbb {R} ^d \rightarrow \mathbb {R} ^d$ be a smooth Lipschitz vector field, and let $y\colon [0,T]\rightarrow \mathbb {R} ^d$ be the solution to the ordinary differential equation
If the vector field f is divergence free, that is, if $\operatorname{\mathrm{div}}(f)=0$ , it is known that the solution of the ODE (1.1) is volume preserving, that is, for any measurable set of initial conditions D with respect to the Lebesgue measure $\lambda $ , for any $t>0$ , the flow $\varphi _t$ of the ODE (1.1) satisfies $\lambda (\varphi _t(D))=\lambda (D)$ . Divergencefree vector fields appear in a variety of concrete dynamical systems, for instance in fluid dynamics, meteorology or molecular dynamics. To integrate such systems, it is fundamental to use numerical integrators that also preserve the volume. We discretize the time interval $[0,T]$ into $N+1$ equidistant steps $t_n=nh$ , with h the timestep size, and we choose a onestep integrator
For large classes of integrators, the standard methodology to study volume preservation uses backward error analysis (see, for instance, the textbook [Reference Hairer, Lubich and Wanner23]). The integrator (1.2) can be interpreted formally as the exact solution of a modified ODE
where the modified vector field $\widetilde {f}$ typically depends on f and its derivatives. Then, the integrator is volume preserving if and only if $\operatorname{\mathrm{div}}(\widetilde {f})=0$ .
There is a considerable literature on volumepreserving methods, with many applications for solving a variety of dynamical systems. The existing volumepreserving integrators rely either on a specific form of the vector field f [Reference Hairer, Lubich and Wanner23, Reference Chartier and Murua18], on splitting methods [Reference Kang and Shang34, Reference McLachlan and Quispel45, Reference McLachlan, MuntheKaas, Quispel and Zanna44, Reference Xue and Zanna63], or on generating functions [Reference Scovel55, Reference Shang57, Reference Quispel53]. For quadratic differential equations, we mention the works [Reference Petrera, Pfadler and Suris51, Reference Celledoni, McLachlan, Owren and Quispel16, Reference Celledoni, McLachlan, McLaren, Owren and Quispel15, Reference Celledoni, Evripidou, McLaren, Owren, Quispel and Tapley14, Reference Bogfjellmo, Celledoni, McLachlan, Owren and Quispel7] that study the Kahan–Hirota–Kimura discretization [Reference Kahan33, Reference Hirota and Kimura29, Reference Hirota and Kimura30] for the preservation of measures. The splitting approach relies on the knowledge of the exact flows involved or depends on the dimension of the problem and is limited to order two in the case of nonreversible problems [Reference Blanes and Casas5]. The methodology with generating functions works with any vector field, but its complexity increases with the dimension of the problem and the approach requires the evaluation of multiple integrals per step. An important open question in geometric numerical integration is the creation of a volumepreserving method for solving a general ODE of the form (1.1) with a complexity independent of the dimension of the problem. We investigate in this work the volumepreserving aromatic Bseries method.
Introduced in [Reference Butcher10, Reference Hairer and Wanner24] (see also the textbooks [Reference Hairer, Lubich and Wanner23, Reference Butcher11, Reference Butcher12] and the review [Reference McLachlan, Modin, MuntheKaas and Verdier43]), the Butcherseries formalism is an important tool in numerical analysis. Originally used for the calculation of order conditions for Runge–Kutta methods, the use of Bseries was quickly extended to a variety of applications, in particular in geometric numerical integration [Reference Hairer, Lubich and Wanner23], and more recently in the approximation of stochastic evolutionary problems [Reference Burrage and Burrage9, Reference Komori, Mitsui and Sugiura35, Reference Rößler54, Reference Debrabant and Kværnø19, Reference Laurent and Vilmart38, Reference Laurent and MuntheKaas37] or in the theory of rough paths [Reference Gubinelli22, Reference Hairer and Kelly26]. For a large class of integrators, such as Runge–Kutta methods, the modified vector field $\widetilde {f}$ in equation (1.3) is a Bseries in f. An extension of Bseries, called aromatic Bseries, was introduced in [Reference Chartier and Murua18, Reference Iserles, Quispel and Tse32] for the study of volumepreserving integrators. They allow to compute the divergence of a Bseries and were later studied in [Reference McLachlan, Modin, MuntheKaas and Verdier42, Reference MuntheKaas and Verdier47, Reference Bogfjellmo6, Reference Fløystad, Manchon and MuntheKaas21] for their algebraic and geometric properties. In [Reference Chartier and Murua18, Reference Iserles, Quispel and Tse32], it is showed that no nontrivial Bseries is divergence free so that the only volumepreserving Bseries method is the exact flow. In particular, no Runge–Kutta method can preserve volume exactly (see also [Reference Kang and Shang34]). However, the space of divergencefree aromatic Bseries is infinite dimensional, the simplest nontrivial example being:
In [Reference MuntheKaas and Verdier47], the question whether there exists a volumepreserving integrator that has an expansion as an aromatic Bseries is raised. The present article gives a handful of tools to answer this question and describes explicitly the aromatic Bseries of a volumepreserving integrator. In particular, we prove that no aromatic Runge–Kutta method can preserve volume exactly.
When using Bseries, the study of volume preservation translates into the study of the linear combinations of aromatic forests of vanishing divergence that we call solenoidal forms. To the best of our knowledge, there does not exist any tool in the literature to describe these forests. In differential geometry, the study of differential forms of vanishing divergence is done with the De Rham complex [Reference Lee40]. For instance, in dimension 3, denote $\Omega _n^{\text {(diff)}}(\mathbb {R} ^3)$ the space of smooth differential nforms on $\mathbb {R} ^3$ , then the De Rham complex writes
where d is the exterior derivative. In this particular example, the arrows correspond in order to the gradient, the curl and the divergence operators. The chain (1.5) is called a complex as the composition of two successive maps vanishes. Moreover, the De Rham complex (1.5) is exact; that is, the image of a map is exactly the kernel of its successor. For instance, a divergencefree form $\omega \in \Omega _3(\mathbb {R} ^3)~$ is a curl $\omega =d\eta $ . This exactness property is typically proven via the use of homotopy operators. The analysis presented in this paper relies heavily on a generalisation of the De Rham complex, called the variational bicomplex, that we extend in the context of aromatic forests.
The variational bicomplex was originally introduced in the context of differential geometry [Reference Vinogradov60, Reference Tulczyjew59, Reference Tsujishita58, Reference Vinogradov61, Reference Vinogradov62] as a natural and general development of the variational chain. It has a variety of applications in the areas of differential geometry and topology, (partial) differential equations and mathematical physics (see, for instance, the textbooks [Reference Anderson1, Reference Olver49], the introductory article [Reference Anderson2] and references therein). In this work, we introduce an algebraic tool on aromatic forests that we call the aromatic bicomplex, in the spirit of the variational bicomplex. If one considers the elementary differentials associated to the aromatic forests in the bicomplex, it yields a subcomplex of the standard variational bicomplex. An originality of the approach is that the analysis of the aromatic bicomplex uses simple combinatorics and graph theory and avoids the technical details of differential geometry. We emphasize that, unlike the analysis of the variational bicomplex, the analysis of the aromatic bicomplex is completely independent of the dimension of the problem. We will also add the extra assumption $\operatorname{\mathrm{div}}(f)=0$ , whose effect has never been studied, to the best of our knowledge, in the context of the variational bicomplex. Thanks to the exactness of the aromatic bicomplex, we describe completely the solenoidal forms, both in the standard context and under the assumption $\operatorname{\mathrm{div}}(f)=0$ . We provide new operations on aromatic forests, such as the Euler operators and homotopy operators, and we draw links with the process of integration by parts of trees described in [Reference Laurent and Vilmart38, Reference Laurent and Vilmart39]. The main application of this work is the explicit description of the aromatic Bseries of a volumepreserving integrator. In particular, we show that no aromatic Runge–Kutta method can preserve volume, and we propose a possible new ansatz for the creation of volumepreserving aromatic Bseries methods.
The paper is organised as follows. Section 2 introduces the aromatic forests and forms and the aromatic bicomplex and presents the main theoretical results of this paper. In Section 3, we introduce the Euler operators and study the exactness of the aromatic bicomplex, in the standard context and in the case of a divergencefree vector field. In Section 4, we present different extensions and applications of the aromatic bicomplex. More precisely, we introduce the augmented aromatic bicomplex and the Euler–Lagrange complex, we compute exactly the number of solenoidal forms, we derive bases and properties on divergences and solenoidal forms, and we draw links with the different integration by parts process of trees existing in the literature. Finally, we apply our results to the study of volumepreserving integrators to obtain an explicit description of the Bseries of an aromatic volumepreserving integrator.
2 Preliminaries and main results
This section is devoted to the presentation of the new objects and to the associated main algebraic results. We first recall the definition of aromatic forests and present an extension, called aromatic forms, well suited for studying the divergencefree combination of forests. Using grafting operations, we define the aromatic bicomplex, a similar tool to the variational bicomplex in the context of differential geometry, and we use its exactness to describe aromatic forms of vanishing divergence. The proofs of the results of this section, as well as further tools, such as the Euler operators, are presented in Section 3, while we give more details and concrete applications in Section 4.
2.1 Aromatic forms: definition and operations
Bseries were introduced by Hairer and Wanner in [Reference Hairer and Wanner24], based on the work of Butcher [Reference Butcher10]. Their applications in the numerical analysis of deterministic differential equations are numerous (see, for instance, the textbooks [Reference Butcher12, Reference Hairer, Lubich and Wanner23]). The aromatic extension of Butcherseries was introduced independently in the works [Reference Chartier and Murua18, Reference Iserles, Quispel and Tse32] to study volumepreserving integrators. In particular, this extension allows us to represent the divergence of standard Bseries, which is a key tool in the study of volumepreserving methods (see [Reference Hairer, Lubich and Wanner23, Sect. VI.9]). In the spirit of differential geometry, we work in this paper with an extension of aromatic Bseries that is analogous to differential forms, and we follow the graph definition of aromatic Bseries of [Reference Bogfjellmo6].
Definition 2.1. Let V be a finite set of nodes and $E\subset V\times V$ a set of edges. If $a=(v,w)\in E$ , the edge a is going from v to w, and v is a predecessor of w. We split the set V into vertices $V^{\bullet }$ and covertices $V^{\circ }$ . The covertices are numbered from $1$ to p, while the vertices are not numbered. Each node in V is the source of exactly one edge, except the roots that have no outgoing edges, that we order and number from $1$ to n. Any connected component of such a graph either has a root, and is called a tree, or does not have a root, and is called an aroma. We call aromatic forests such graphs, up to equivalence of graphs that preserve the numbering of the covertices and the roots. We write $\mathcal {F} _{n,p}$ the set of aromatic forests with n roots and p covertices, $\mathcal {F} _{n,p}^N$ its subset with forests of exactly N nodes, and $\mathcal {F} _{n}=\mathcal {F} _{n,0}$ . The number of nodes $\left \vert \gamma \right \vert $ is called the order of the forest $\gamma $ . The elements of $\mathcal {F} _1$ are called aromatic trees and the subset $\mathcal {T} $ of $\mathcal {F} _1$ that contains the trees without aromas is the set of Butcher trees.
We draw the aromatic forests as follows. The vertices are represented as black nodes and the covertices as circles of the form , where i is the associated number. The trees are drawn in the ascending order of their roots, from left to right. The aromas are placed in front and their order does not matter. The orientation of the edges goes from top to bottom and in clockwise order for loops. A loop with K nodes is called a Kloop, in the spirit of [Reference Iserles, Quispel and Tse32]. For instance, the following forest $\gamma \in \mathcal {F} _{3,2}^{10}$ has two aromas, one $1$ loop and one $3$ loop,
The aromatic forests with $N=2$ nodes are
An aromatic forest in $\mathcal {F} _{n,p}$ represents a tensor on the infinite jet bundle through the application of the elementary differential map. The definition of this map relies on contact forms on the infinite jet bundle, of which we give a brief definition below. These details can be ignored as we shall work in the following with aromatic forests only, but they draw the link with the variational bicomplex (see Remark 2.8). We refer the reader to [Reference Anderson2] for further details (see also the textbooks [Reference Olver49, Reference Kushner, Lychagin and Rubtsov36, Reference Lee40]). The infinite jet bundle $J^\infty (\mathbb {R} ^d)$ is used to describe coordinatefree Taylor expansions. It is the infinitedimensional vector bundle over $\mathbb {R} ^d$ where elements of $J^\infty _x(\mathbb {R} ^d)$ have coordinates
We consider differential forms on $J^\infty (\mathbb {R} ^d)$ and define the following contact forms
The set $\{dx^i,\theta ^i,\theta ^i_j,\dots \}$ locally generates all differential forms on $J^\infty (\mathbb {R} ^d)$ . A differential form of type $(n,p)$ is a sum of terms of the form
where $g(f)$ is a functional of finitely many coordinates (2.2). The differential forms of type $(n,p)$ are gathered in the set $\Omega ^{\text {(diff)}}_{n,p}(J^\infty (\mathbb {R} ^d))$ . In agreement with the numerical analysis literature, the elementary differential map is defined with the help of the vector basis $\partial _i$ .
Definition 2.2. Let $\gamma \in \mathcal {F} _{n,p}$ , $f\colon \mathbb {R} ^d\rightarrow \mathbb {R} ^d$ a smooth vector field and $R=\{r_1, \dots , r_n\}\subset V$ the n roots of $\gamma $ , then the elementary differential $F(\gamma )(f)$ is the following tensor:
where $\Pi (v)$ is the set of predecessors of $v\in V$ . We extend F on $\operatorname{\mathrm{Span}}(\mathcal {F} _{n,p})$ by linearity.
For instance, we find
The forest $\gamma $ in equation (2.1) represents the elementary differential
In the following, we work with specific linear combinations of aromatic forests in $\operatorname{\mathrm{Span}}(\mathcal {F} _{n,p})$ . Given an elementary tensor, one can alternatize it to obtain a differential form. The same idea gives rise to the aromatic forms.
Definition 2.3. For $\gamma \in \mathcal {F} _{n,p}$ , let $\mathcal {S} _n^{\bullet }$ (resp. $\mathcal {S} _p^{\circ }$ ) be the set of permutations of the roots of $\gamma $ (resp.the covertices of $\gamma $ ). We define the roots wedge of $\gamma $ as
where $\varepsilon (\sigma )$ is the signature of the permutation $\sigma $ . Similarly, the covertices wedge is
and the total wedge is $\wedge =\wedge ^{\bullet }\wedge ^{\circ }=\wedge ^{\circ }\wedge ^{\bullet }$ .
We extend the wedge operations to $\operatorname{\mathrm{Span}}(\mathcal {F} _{n,p})$ by linearity and we denote the set of aromatic forms $\Omega _{n,p}=\wedge \operatorname{\mathrm{Span}}(\mathcal {F} _{n,p})$ , respectively $\Omega _{n,p}^N=\wedge \operatorname{\mathrm{Span}}(\mathcal {F} _{n,p}^N)$ , and $\Omega _n=\Omega _{n,0}$ . As $\wedge ^2=\wedge $ , the operator $\wedge \colon \operatorname{\mathrm{Span}}(\mathcal {F} _{n,p})\rightarrow \Omega _{n,p}$ is a projection on $\Omega _{n,p}$ .
We now define the grafting and replacing operations on aromatic forests.
Definition 2.4. Let $\gamma \in \mathcal {F} _{n,p}$ be a forest, r a root of $\gamma $ and $u\in V$ (possibly equal to r), then $D^{r\rightarrow u} \gamma $ returns a copy of $\gamma $ , where the node r is now a predecessor of u. We define the operator $D^r \gamma =\sum _{u\in V} D^{r\rightarrow u} \gamma $ . Let $\gamma \in \mathcal {F} _{n,p}$ and $v\in V^\bullet $ , we define as the forest obtained by replacing the node v by a new covertex . Similarly, is the forest obtained by replacing the covertex by a vertex.
Example. Let
and r its root, then
The horizontal derivative is defined using grafting operations, while the vertical derivative uses the replacing operation.
Definition 2.5. Let $\gamma \in \mathcal {F} _{n,p}$ , the horizontal and vertical derivatives are
We extend $d_H$ and $d_V$ on $\Omega _{n,p}$ by linearity into $d_H\colon \Omega _{n,p} \rightarrow \Omega _{n1,p}$ and $d_V\colon \Omega _{n,p} \rightarrow \Omega _{n,p+1}$ . The aromatic forms in $\operatorname{\mathrm{Ker}}(d_H)$ are called solenoidal forms, $\Psi =\operatorname{\mathrm{Ker}}(d_H_{\Omega _{1}})$ are the solenoidal combinations of trees and $\Psi ^N=\operatorname{\mathrm{Ker}}(d_H_{\Omega _{1}^N})$ are the solenoidal combinations of trees of order N.
The operator $d_H\colon \Omega _1\rightarrow \Omega _0$ is often called the divergence of an aromatic tree, as for $\gamma \in \Omega _1$ , $d_H$ satisfies $\operatorname{\mathrm{div}}(g)=F(d_H(\gamma ))(f)$ , where $F(\gamma )(f)=\sum _i g^i \partial _i$ .
Example. Consider
,
, and
, then
A calculation yields $d_H^2\gamma _2=0$ , so that $d_H\gamma _2\in \Psi $ is a solenoidal form.
The following result, proven in Subsection 3.2, allows us to define the bicomplex.
Proposition 2.6. The horizontal and vertical derivatives satisfy
Remark 2.7. The horizontal and vertical derivatives were already defined on the set of aromatic trees $\Omega _1=\operatorname{\mathrm{Span}}(\mathcal {F} _1)$ , respectively, in [Reference Chartier and Murua18, Reference Iserles, Quispel and Tse32] for $d_H$ and in [Reference Fløystad, Manchon and MuntheKaas21] for $d_V$ . In this last work, the trace operator $\operatorname{\mathrm{Tr}}\colon \Omega _{1,1}\rightarrow \Omega _0$ on aromatic trees is studied. For $\gamma \in \Omega _{1,1}$ , it is given with our notations by
and it makes the following diagram commute.
2.2 The aromatic bicomplex: exactness and description of solenoidal forms
The variational bicomplex is a powerful tool of variational calculus [Reference Anderson2]. We introduce in the context of aromatic forms a tool in the spirit of the variational bicomplex. This new complex, that we call the aromatic bicomplex, allows us in particular to describe explicitly the solenoidal forms in the standard context and in the divergencefree case.
The aromatic bicomplex is the diagram drawn in Figure 1. We also introduce its variant with forms of fixed order N. The aromatic bicomplex can be completed by an extra column on the right in order, for instance, to describe the aromatic forms in $\operatorname{\mathrm{Im}}(d_H_{\Omega _{1,p}})$ , as we will see in Subsection 4.1. We refer the reader to the appendix 2 for examples of the (augmented) aromatic bicomplex for the first values of N.
Remark 2.8. For a fixed dimension d, the elementary differential map F (Definition 2.2) makes a link between the aromatic bicomplex and a subcomplex of the variational bicomplex [Reference Anderson2] in the following way. Let $\sigma $ be the standard duality map between the type nalternating contravariant tensors and the type $(dn)$ differential forms with respect to the standard Euclidean coordinates and volume form $dx^1\wedge \dots dx^d$ (see, for instance, [Reference Lee40]). Let $f\colon \mathbb {R} ^d\rightarrow \mathbb {R} ^d$ , the map $\sigma F(.)(f)$ sends an aromatic form in $\Omega _{n,p}$ to a differential form in $\Omega _{dn,p}^{\text {(diff)}}(J^\infty (\mathbb {R} ^d))$ . The derivatives $d_H$ , $d_V$ on the aromatic bicomplex and $d_H^{\text {(diff)}}$ , $d_V^{\text {(diff)}}$ on the variational bicomplex satisfy
Note that the dimension d of the problem plays a role in the context of differential geometry but not in the context of aromatic forms. The exactness results presented in this paper translate directly to new results on a subcomplex of the variational bicomplex through the application of $\sigma F(.)(f)$ . Note also that $d_H$ and $d_V$ commute (see Proposition 2.6), while their counterparts from differential geometry anticommute.
The main property of the aromatic bicomplex is its exactness.
Theorem 2.9. The horizontal and vertical sequences of the aromatic bicomplex are exact, that is, for all n, $p\geq 0$ ,
With Theorem 2.9, it is straightforward to generate all the solenoidal aromatic forms by considering the $d_H \wedge \gamma $ for $\gamma \in \mathcal {F} _2$ . For example, the only basis element of $\Psi ^3$ (which corresponds to the vector field (1.4)) is
and a basis of $\Psi ^4$ is given by the forests
However, the set $\{d_H \wedge \gamma , \gamma \in \mathcal {F} _2\}$ does not form a basis of the solenoidal forms in general. We give a basis of $\Psi $ in Subsection 4.3 alongside bases of the image and kernel of $d_H$ and its dual $d_H^*$ .
For the study of volumepreserving integrators for solving equation (1.1), it is fundamental to assume that the vector field f satisfies $\operatorname{\mathrm{div}}(f)=0$ . With aromatic forms, it amounts to sending all forests containing an aroma with a 1loop to $0$ . We consider the vector space $\mathcal {A} $ spanned by aromatic forests containing at least one 1loop, $\mathcal {A} _{n,p}=\mathcal {A} \cap \Omega _{n,p}$ , and $\mathcal {A} _{n,p}^N=\mathcal {A} \cap \Omega _{n,p}^N$ . We write $\widetilde {\mathcal {F} }_{n,p}$ the set of aromatic forests in $\mathcal {F} _{n,p}$ without 1loops, $\widetilde {\Omega }_{n,p}=\Omega _{n,p}/\mathcal {A} _{n,p}$ , $\widetilde {\Omega }_{n,p}^N=\Omega _{n,p}^N/\mathcal {A} _{n,p}^N$ , $\widetilde {\Omega }_n=\widetilde {\Omega }_{n,0}$ , $\widetilde {\Psi }=\operatorname{\mathrm{Ker}}(d_H_{\widetilde {\Omega }_{1}})$ , and $\widetilde {\Psi }^N=\operatorname{\mathrm{Ker}}(d_H_{\widetilde {\Omega }_{1}^N})$ . The divergencefree aromatic bicomplex with N nodes is drawn in Figure 2. In the simplest case $N=1$ (see Figure 2), the aromatic bicomplex is not exact. Indeed, we observe that , but .
One of the main results of this paper is that the case $N=1$ is the only case where the divergencefree aromatic bicomplex is not exact.
Theorem 2.10. The divergencefree aromatic bicomplex with N nodes is exact if and only if $N\neq 1$ .
The main consequence of the exactness of the aromatic bicomplex in both contexts is that the assumption $\operatorname{\mathrm{div}}(f)=0$ does not create new nontrivial solenoidal forms.
Theorem 2.11. For $N\neq 1$ , and all $n\geq 1$ , $p\geq 0$ , the kernel of the divergence operator $d_H$ satisfies
that is, the solenoidal aromatic forms in the divergencefree context exactly correspond to the solenoidal aromatic forms in the standard context, except for the forms spanned by and .
Proof of Theorem 2.11
As $d_H$ does not decrease the number of $1$ loops in a forest, the image of $d_H$ satisfies
A generating set of the solenoidal forms $\widetilde {\Psi }^N$ of order $N>1$ is obtained by deleting the $1$ loops in the generating set $\{d_H \wedge \gamma , \gamma \in \mathcal {F} _2^N\}$ of $\Psi ^N$ . For instance, for order $N=3$ , we derive from the element (2.3) that the solenoidal forms in the divergencefree context are given by
We give generators of the solenoidal forms $\widetilde {\Psi }^N$ for the first orders in Appendix A. This explicit description of solenoidal forms is especially useful in the numerical study of volumepreserving integrators, as discussed in Subsection 4.5.
We enumerate in Subsection 4.2 the dimensions of the $\Omega _{n,p}^N$ and $\widetilde {\Omega }_{n,p}^N$ in the first two rows of the bicomplex and deduce the dimension of $\Psi ^N$ in Theorem 4.2 and of $\widetilde {\Psi }^N$ in Theorem 4.4. A surprising fundamental result is that the solenoidal forms in $\Psi ^N$ are enumerated by the difference between the number of aromatic trees in $\Omega _{1}^N$ and the number of aromas with 1loops $\mathring {\Omega }_{0}^N$ (see Table 1 for examples).
3 The aromatic bicomplex
This section is devoted to the study of the aromatic bicomplex. First, we introduce the Euler operators and prove the exactness of the variational chain in Subsection 3.1. We define the horizontal and vertical homotopy operators and prove the exactness of the aromatic bicomplex in Subsection 3.2. We study the aromatic bicomplex in the divergencefree context in Subsection 3.3.
3.1 The Euler operators and the variational complex
To define the Euler operators, we extend Definition 2.4 to allow one to detach edges and graft them back.
Definition 3.1. A graph $\gamma $ is a marked aromatic forest if it is an aromatic forest with exactly one of its node that holds the symbol $\diamond $ . We write $\mathcal {F} _{n,p}^\diamond $ the set of marked aromatic forests with n roots and p covertices. Given $\gamma \in \mathcal {F} _{n,p}^\diamond $ , $\gamma ^\diamond \in \mathcal {F} _{n,p}$ is the same forest $\gamma $ where the symbol $\diamond $ is removed.
Let $\gamma \in \mathcal {F} _{n,p}$ or $\gamma \in \mathcal {F} _{n,p}^\diamond $ , where $v_\diamond \in V$ is the node with the symbol $\diamond $ . Let $R_0\subset R$ be a given subset of roots that we call the set of detached nodes. For $\gamma \in \mathcal {F} _{n,p}^\diamond $ , $r\in R$ and $u\in V$ with $u\neq v_\diamond $ , $D^{r\rightarrow u} \gamma \in \mathcal {F} _{n1,p}^\diamond $ returns a copy of $\gamma $ with an additional edge going from r to u. The set of detached roots of $D^{r\rightarrow u} \gamma $ is $R_0$ if $r\notin R_0$ and $R_0\setminus \{r\}$ else. We define
Let q be a nonnegative integer, $u\in V$ with $u\neq v_\diamond $ if $\gamma \in \mathcal {F} _{n,p}^\diamond $ , and $I\subset R_0$ . We define $D^I$ and $D^{I\rightarrow u}$ by
and $D^q$ and $D^{q\rightarrow u}$ by
Let $\gamma \in \mathcal {F} _{n,p}$ and $v\in V$ , then $\gamma _{v^{\diamond }}\in \mathcal {F} _{n+\left \vert \Pi (v)\right \vert ,p}$ is the graph obtained by cutting off the edges of $\gamma $ pointing to v and placing the newly obtained roots after the roots in R (in an arbitrary order). In addition, we add a symbol $\diamond $ on the node v and we fix $R_0=\Pi (v)$ as the set of detached nodes.
The tools from Definition 3.1 are extended by linearity. In simple words, the $\diamond $ detaches the predecessors of a given node v and puts them in a set $R_0$ . As long as the symbol $\diamond $ is present on v, the grafting operators $D^I$ and $D^q$ graft the detached nodes on every node except v. In particular, the following operation is trivial:
Example. Consider the forest
that we label for the sake of clarity. The marked forest $\gamma _{1^{\diamond }}$ has the set of detached nodes $R_0=\{1\}$ and satisfies
For simplicity, we write for , and $D=D^1$ in the rest of the paper. Note that, in general, $D^2$ and $D\circ D$ are different operators. Instead, straightforward combinatorics yield
However, for u, $v\in R$ , the operations $D^u$ and $D^v$ commute.
Remark 3.2. The notations introduced in Definition 3.1 closely depend on the forest $\gamma $ they are applied to. For the sake of clarity, in the rest of the paper, the notations shall always relate to the original forest denoted $\gamma $ . For instance, in $(1)^{\left \vert \Pi (v)\right \vert } \gamma _{v^{\diamond }}$ , $\Pi (v)$ denotes the number of predecessors of v in $\gamma $ , not in $\gamma _{v^{\diamond }}$ .
In the context of differential geometry, the Euler operator describes the differential forms that are divergences [Reference Olver49]. It has a variety of applications such as the study of conservation laws and Lagrangians for partial differential equations [Reference Olver and Shakiban50, Reference Shakiban56, Reference Hereman, Colagrosso, Sayers, Ringler, Deconinck, Nivala and Hickman27, Reference Hereman, Deconinck and Poole28, Reference Poole and Hereman52]. We define a similar operator for aromatic forms.
Definition 3.3. For $\gamma \in \mathcal {F} _{n,p}$ and $v\in V$ , the Euler operators $\mathcal {E} _v \colon \mathcal {F} _{n,p}\rightarrow \Omega _{n,p}$ and $\mathcal {E} _v^{\circ }\colon \mathcal {F} _{0}\rightarrow \Omega _{0,1}$ are given by
and are extended by linearity on $\Omega _{n,p}$ . The Euler operators $\mathcal {E} \colon \Omega _{n,p}\rightarrow \Omega _{n,p}$ and its variant $\mathcal {E} ^{\circ }\colon \Omega _{0}\rightarrow \Omega _{0,1}$ are given by
The output of $\mathcal {E} ^{\circ }$ on the forms of first orders is
From these examples, we observe that the Euler operator vanishes on aromatic forms that are divergences:
Proposition 3.4. Let $\gamma \in \mathcal {F} _{n,p}$ and $v\in V$ , the Euler operators satisfy
and for $p=0$ , $\mathcal {E} ^{\circ }$ satisfies $\mathcal {E} ^{\circ } d_H\gamma =0.$
Proof. Let $\gamma \in \mathcal {F} _{n,p}$ and $r_n$ be its last root, then
where the two terms correspond to the cases $u=v$ and $u\neq v$ and where we used the convention of Remark 3.2. The identities with $\mathcal {E} $ and $\mathcal {E} ^{\circ }$ are direct consequences.
We know that the composition of the two maps $\mathcal {E} ^{\circ }$ and $d_H$ vanishes. One is then interested in a necessary and sufficient condition for an aromatic form $\gamma \in \Omega _0$ to be a divergence. The following chain is called the variational complex.
The complex (3.2) is exact when $\operatorname{\mathrm{Im}}(d_H)=\operatorname{\mathrm{Ker}}(\mathcal {E} ^{\circ })$ , that is, when $\gamma \in \Omega _0$ is a divergence if and only if $\mathcal {E} ^{\circ } \gamma =0$ . A fundamental question in variational calculus [Reference Olver49] is the exactness of this chain (in the context of differential forms). We prove in the rest of this subsection the exactness of the variational complex (3.2) in the context of aromatic forms. We rely on the use of the Euler operators of higher orders and homotopy operators. To the best of our knowledge, the use of such operators on aromatic forests is completely new.
Definition 3.5. For $\gamma \in \mathcal {F} _{n,p}$ , the Euler operator $\mathcal {E} ^q_v \gamma $ of order $q\geq 0$ on $v\in V$ is
that we extend on $\Omega _{n,p}$ by linearity. The higher Euler operators are
A fundamental result for our analysis is that any aromatic forest can be rewritten with the Euler operators. Note that all appearing series are finite, as $\mathcal {E} ^q_v \gamma =0$ if $q> \left \vert \gamma \right \vert $ . Note also that $\mathcal {E} ^0_v=\mathcal {E} _v$ .
Proposition 3.6. The Euler operators satisfy for all $\gamma \in \mathcal {F} _{n,p}$ and all vertices $v\in V$ ,
In particular, we have
In addition, if $\gamma $ satisfies $\mathcal {E} ^{p}\gamma =0$ for $p=0,\dots ,n1$ , then $\gamma =D^n(h^{(n)}\gamma )$ , where
The proof shares similarities with the approach of [Reference Hereman, Deconinck and Poole28] and uses the following intermediate result in the spirit of the Leibniz rule.
Lemma 3.7. Let $\gamma \in \mathcal {F} _{n,p}$ and $v\in V$ ; let three integers p, q, k such that $p+q+k=\left \vert \Pi (v)\right \vert 1$ . Then, the following holds:
More precisely, we have
Proof of Proposition 3.6
Let $\gamma \in \mathcal {F} _{n,p}$ . With the help of Lemma 3.7, we get
We apply the same reasoning to the last term
Iterating this reasoning, we find
Using the formula (3.5) in equation (3.6), standard combinatorics [Reference Merris46, Ex. 4.2.7] yield
We sum equation (3.3) on all the nodes $v\in V$ to obtain equation (3.4). The last claim of Proposition 3.6 is obtained by using the formula (3.1) in equation (3.4).
Proof of Lemma 3.7
Using the identity (3.1), we distribute the derivatives on v and the other nodes.
We obtain the formula (3.5) by induction on k.
A direct consequence of Proposition 3.6 is the exactness of the variational complex (3.2).
Theorem 3.8. For $\gamma \in \mathcal {F} _{0,1}$ , let the variational homotopy operator $h_V$ be
and for $\gamma \in \mathcal {F} _0$ , let the horizontal homotopy operator $h_H$ be
We extend the definition on $\Omega _{0,1}$ and $\Omega _0$ by linearity. These operators satisfy for all $\gamma \in \Omega _0$ ,
In particular, the variational complex (3.2) is exact.
Proof. Let $\gamma \in \Omega _0$ . From Proposition 3.6, we obtain
Proposition 3.4 yields that $\operatorname{\mathrm{Im}}(d_H)\subset \operatorname{\mathrm{Ker}}(\mathcal {E} )$ . If $\gamma \in \operatorname{\mathrm{Ker}}(\mathcal {E} )$ , then the identity (3.7) becomes $\gamma =d_H (h_H \gamma )\in \operatorname{\mathrm{Im}}(d_H)$ . The exactness of the chain follows straightforwardly.
We present an alternative horizontal homotopy operator on $\Omega _0$ in Subsection 4.4, with some examples of the outputs of both operators in Table 3.
3.2 Exactness of the aromatic bicomplex
This section is devoted to the proof of Theorem 2.9. We start by showing the aromatic bicomplex is indeed a bicomplex.
Proof of Proposition 2.6
Let $\gamma \in \mathcal {F} _{n,p}$ , and $(r_1,\dots ,r_n)$ its roots. The horizontal derivative satisfies
where $(n(n1))\in \mathcal {S} _{n}^{\bullet }$ is a transposition and where we used that $D^{r_{n1}}$ and $D^{r_n}$ commute. For the vertical derivative, a similar approach on $\gamma \in \mathcal {F} _{n,p1}$ gives
As $\Omega _{n,p}=\wedge (\mathcal {F} _{n,p})$ , we get the desired identities. The commutativity of $d_H$ and $d_V$ is straight forward.
In order to prove the exactness of the sequences appearing in the variational bicomplex, we rely on the use of two homotopy operators. We begin by the vertical homotopy as it is the easiest. Note that the variational homotopy operator and the vertical homotopy operator introduced in Theorem 3.8 coincide for $p=1$ .
Proposition 3.9. The vertical homotopy operator $h_V\colon \Omega _{n,p}\rightarrow \Omega _{n,p1}$ given by
satisfies for $p\geq 1$ and $\gamma \in \Omega _{n,p}$ ,
Proof. Let $\gamma \in \mathcal {F} _{n,p}$ . On the first hand, we have
On the other hand, we have
We deduce that
and we obtain the desired equality by linearity as $\Omega _{n,p}=\wedge \operatorname{\mathrm{Span}}(\mathcal {F} _{n,p})$ and $\wedge ^2=\wedge $ .
The expression of the horizontal homotopy operator is much more technical. We refer the reader to [Reference Anderson and Duchamp3, Reference Tulczyjew59, Reference Olver49] for its derivation in the context of differential geometry. Note that for $n=0$ , the horizontal homotopy operator coincides with the one introduced for the variational complex in Theorem 3.8.
Proposition 3.10. We define the horizontal homotopy operator $h_H\colon \Omega _{n,p}\rightarrow \Omega _{n+1,p}$ by
It satisfies for $n\geq 1$ and $\gamma \in \Omega _{n,p}$ ,
Proof. For the sake of simplicity, we consider $\gamma \in \mathcal {F} _n$ . We define for $v\in V$
where we denote $\mathcal {E} ^{I,\hat {r}}_{v} \gamma =(1)^{\left \vert \Pi (v)\right \vert (q+1)} (D^J \gamma _{v^{\diamond }})^{\diamond }$ and where $\hat {r}$ becomes the new root of $J_v \gamma $ . On the first hand, we have
and $J_v d_H \wedge \gamma ~$ is given by
The expression of $D^I \mathcal {E} ^{I,\hat {r}}_{v} D^{r\rightarrow u} \gamma $ satisfies for $u\in V$ :
Thus, for $\gamma \in \mathcal {F} _{n,p}$ , $J_v d_H \wedge \gamma ~$ is given by
As $D^q \mathcal {E} ^{q}_{v} \wedge \gamma $ is unchanged by the application of the wedge operator, we find
On the other hand, we have
where we split the cases $\sigma \nu (n+1)=n+1$ and $\sigma \nu (n+1)\neq n+1$ and where we substitute $\nu $ by $\nu (n(n+1))\eta $ in the latter case. The last equality is obtained by substituting $\sigma $ with $\nu \sigma $ . Using Proposition 3.6, we get
Summing on $v\in V$ gives the desired homotopy identity.
Remark 3.11. The identity (3.8) suggests simpler homotopy operators, obtained by fixing a node v and considering the operator $\wedge J_v$ . We emphasize that this approach does not work. Indeed, given an aromatic form $\gamma \in \Omega _{n,p}$ , there is no canonical choice of the node v so that $J_v \gamma $ is illdefined, where $h_V \gamma $ is well defined. In the proof of Proposition 3.10, $J_v \gamma $ makes sense since we consider a single aromatic forest $\gamma \in \mathcal {F} _{n,p}$ .
3.3 The aromatic bicomplex with a divergencefree vector field
This section is devoted to the proof of Theorem 2.10. Proposition 2.6 and Proposition 3.9 extend naturally to the divergencefree context. Proposition 3.10 is not valid anymore, and we replace it by the following result.
Proposition 3.12. We define
For $n>1$ and $\gamma \in \widetilde {\Omega }_{n,p}$ , the horizontal homotopy identity is
while for $n=1$ , the identity is
The remainder in equation (3.9) is the linear map $R\gamma =\frac {1}{\left \vert \gamma \right \vert }\mathcal {E} _r \gamma $ for $\gamma \in \widetilde {\mathcal {F} }_{1,p}$ , with r the root of $\gamma $ . Moreover, if $\gamma \in \widetilde {\Omega }_{1,p}^N$ satisfies $d_H \gamma =0$ , then $R\gamma =0$ if and only if $N>1$ . In particular, the horizontal sequences in the divergencefree aromatic bicomplex are exact.
The proof of Proposition 3.12 follows the structure and notations of the proof of Proposition 3.10. We recall that the notations in the proof follow the convention of Remark 3.2.
Proof. We consider for simplicity $\gamma \in \widetilde {\mathcal {F} }_n$ . For $v\in V$ , we define
We have if $r=v$ ,
and if $r\neq v$ ,
Thus, $\widetilde {J}_v d_H \wedge \gamma $ is given by
Then, we find for $n=1$ and $v\in R$ ,
and otherwise
where for $v\in R$ , the wedge adds the missing term $v= r_{\sigma ^{1}(n)}$ , and rescales the expression with a coefficient $\frac {n1}{n}$ . On the other hand, we have
Using Proposition 3.6, we get for $n>1$ ,
and for $n=1$ ,
Summing on $v\in V$ gives the desired homotopy identities.
Following the proof of Proposition 3.4, we observe that $\mathcal {E} _r d_H=0$ on $\Omega _{2,p}$ . Thus, we deduce from the identity (3.9) that $\gamma \in \operatorname{\mathrm{Im}}(d_H)$ if and only if $d_H \gamma =0$ and $\mathcal {E} _r \gamma =0$ in the case $n=1$ . Assume that $d_H\gamma =0$ for $\gamma \in \Omega _{1,p}^N$ , and apply $\mathcal {E} _r$ to the homotopy identity. As $\mathcal {E} _r^2=\mathcal {E} _r$ , it yields
We deduce that $\mathcal {E} _r \gamma =0$ or $N=1$ . Thus, the bicomplex of order N is exact if and only if $N\neq 1$ .
Remark 3.13. Let $\gamma \in \widetilde {\Omega }_1^N$ with $N>1$ , define the modified homotopy operators by
Then, the homotopy identity (3.9) on $\widetilde {\Omega }_1^N$ with $N>1$ is replaced by the simpler identity
The proof of (3.10) relies on the identities $\mathcal {E} _r^2=\mathcal {E} _r$ and $\mathcal {E} d_H=d_H \mathcal {E} _r$ on $\widetilde {\Omega }_1$ .
4 Applications and extensions
The study of the aromatic bicomplex brings a variety of new theoretical results, as presented in Section 2, but also direct applications in numerical analysis. Subsection 4.1 is devoted to the study of the generalised aromatic bicomplex, a natural extension of the aromatic bicomplex that includes the Euler–Lagrange complex. In Subsection 4.2, we deduce from the exactness of the aromatic bicomplex the dimensions of the spaces in the first two rows of the bicomplex, as well, as the number of solenoidal forms. We describe further the divergences and the solenoidal forms in Subsection 4.3. In Subsection 4.4, we draw a bridge between the existing notions of integration by parts of Butcher trees by defining an alternative horizontal homotopy operator. In Subsection 4.5, we give an explicit description of the Bseries of an aromatic volumepreserving integrator and we prove that an aromatic Runge–Kutta method cannot be volume preserving.
4.1 The augmented aromatic bicomplex
Following [Reference Anderson1, Reference Anderson2] and the work on the variational complex (3.2) of Subsection 3.1, we augment the aromatic bicomplex with the aromatic equivalent of the Euler–Lagrange complex.
For $p\geq 1$ , define the interior Euler operator $I\colon \Omega _{0,p} \rightarrow \Omega _{0,p}$ by
write $\mathcal {I} _p=I(\Omega _{0,p})$ , $\mathcal {I} _p^N=I(\Omega _{0,p}^N)$ and the variational derivative $\delta _V=I\circ d_V$ . The augmented aromatic bicomplex is drawn in Figure 3. The edge complex (4.1) is called the Euler–Lagrange complex, and it is the object of ultimate interest here. Note that the variational complex (3.2) is a subcomplex of the Euler–Lagrange complex as $\delta _V=\mathcal {E} ^{\circ }$ on $\Omega _0$ .
Theorem 4.1. Define the augmented homotopy operators as
Then, the maps I and $\delta _V$ satisfy
and the following identities hold
In particular, the horizontal and vertical sequences of the augmented aromatic bicomplex are exact, and the Euler–Lagrange complex (4.1) is exact.
We follow the approach of [Reference Anderson1, Chap. 4] for the proof of Theorem 4.1.
Proof. We first observe that for all $\gamma \in \mathcal {F} _{1,p}$ ,
Using equation (3.3) with
gives the augmented horizontal homotopy identity
Applying I to this last equality yields $I\circ I=I$ . Let us now look at the EulerLagrange complex. Let $\gamma \in \Omega _{0,p}$ , then $d_V \gamma \in \Omega _{0,p+1}$ . We apply the augmented horizontal homotopy identity to $d_V \gamma $ ,
We apply $d_V$ and use that $d_V$ and $d_H$ commute:
Since $I d_H=0$ , applying I yields
Let $\gamma \in \mathcal {I} _p$ , the augmented horizontal homotopy applied to $d_V\gamma $ gives
If $p=1$ , then we use the identity (4.2) in the vertical homotopy identity to get
where we used that $h_V$ and $d_H$ commute and where $\widetilde {\gamma }\in \Omega _{1,1}$ . Applying I yields
If $p> 1$ , we apply the augmented horizontal homotopy identity to $h_V\gamma $ ,
Applying the vertical derivative $d_V$ gives
where we used that $d_H$ and $d_V$ commute. We use the identities (4.2) and (4.3) in the vertical homotopy identity to get
where $\widetilde {\gamma }\in \Omega _{1,p}$ . As $I\circ d_H=0$ and $I \gamma =\gamma $ , we find
The exactness of the augmented aromatic bicomplex is a straighforward consequence of Theorem 2.9 and of the augmented homotopy identities.
4.2 Combinatorics on the aromatic bicomplex
In this subsection, we determine the dimensions of the bottom two rows of the augmented aromatic bicomplex in the standard and in the divergencefree case. The primary motivation was to compute the dimension of the space of solenoidal (i.e., divergencefree) aromatic trees of each order. However, the result revealed a surprisingly simple connection with another combinatorial object, the selflooped scalar aromas, which allowed the construction of the fundamental spaces associated with the divergence operator.
We recall the fundamental generating functions associated with graphical enumeration. The number $\mathcal {T} ^N$ of rooted trees of order N (sequence A000081 in the OEIS [48]) has generating function
and satisfies the functional equation
Considering a rooted tree as a directed graph, each node except the root has a single outgoing edge. Thus, rooted trees are equivalent to the class of ‘mapping patterns’ of functions $\{1,\dots ,N1\}\to \{0,\dots ,N1\}$ modulo the symmetric group $\mathcal {S} _{N1}$ (node 0 is the root, and we ‘forget the labels’).
The number $\Omega ^N_{0}$ of scalar aromas in addition to the empty aroma (sequence A001372) has generating function
and is related to $t(z)$ by the equation
The scalar aromas are mapping patterns of functions $\{1,\dots ,N\}\to \{1,\dots ,N\}$ mod $\mathcal {S} _{N}$ .
We introduce two new spaces, the aromatic forms with a $1$ loop, called the selflooped aromatic forms $\mathring {\Omega }_{n,k}$ , and their complements, the nonselflooped aromatic forms $\overline {\Omega }_{n,k}$ , and generating functions $\mathring {a}(z)$ of $\mathring {\Omega }_{0}$ and $\overline {a}(z)$ of $\overline {\Omega }_{0}$ . As for mapping patterns on $\{1,\dots ,N\}$ , the mappings in $\mathring {\Omega }_{0}$ (enumerated by sequence A217896) have at least one fixed point and those in $\overline {\Omega }_{0}$ (also known as the ‘functional digraphs’, enumerated by sequence A001373) have no fixed points.
Each pair consisting of one nonselflooped scalar aroma and one rooted tree generates a scalar aroma of one lower degree: Cut off the root and replace each new root by a selfloop. The process is invertible (redirect all selfloops in an arbitrary scalar aroma to a new root). Expressed in terms of generating functions, it writes
Therefore,
The following result enumerates the first two rows of the augmented bicomplex and the solenoidal forms. We refer to Table 1 and Table 2 for the dimensions for the first orders N.
Theorem 4.2. Let
be the bivariate generating function for row p of the aromatic bicomplex, let
be the generating function of the typep functional forms and let
be the generating function of the solenoidal aromatic trees. Then
Proof. We first note that the aromatic trees $\Omega _{1}$ are given by the product of a scalar aroma (enumerated by $a(z)$ ) with a rooted tree (enumerated by $t(z)$ ). This is the Cartesian product construction of enumerative combinatorics. Therefore, the aromatic trees are enumeratedFootnote ^{1} by $a(z)t(z)$ .
An element of $\Omega ^N_{k}$ is given by the product of a scalar (enumerated by $a(z)$ ) with a wedge product of k rooted trees, enumerated by $t(z)$ . Since the rooted trees are unordered and distinct, this is the power set construction; the expression for $b_0(u,z)$ then follows from [Reference Flajolet and Sedgewick20, Proposition III.5].
The bottom row of the bicomplex is exact (Theorem 2.9), so the dimension of $\Psi $ is given by the alternating sum
An element of $\Omega _{n,1}$ is obtained from an element of $\Omega _{n}$ by marking a node with the symbol . The marked node can be in one of the scalar aroma components or in one of the rooted tree components. Therefore, an element of $\Omega _{n,1}^N$ is either

(i) A scalar aroma with one marked node, times a wedge product of n distinct rooted trees;

(ii) An unmarked scalar aroma times the wedge product of $n1$ distinct rooted trees times a single rooted tree with a marked node. The marked rooted tree can coincide with one of the unmarked ones.
For type (i), we first enumerate the scalar aromas with one marked node in terms of $t(z)$ as follows. Consider the connected component containing the marked node. The marked node lies either on the cycle or on one of the trees attached to the cycle. Those in the first group are enumerated by the rooted trees with one marked node: Delete the edge of the cycle that points to the marked node (the construction is invertible). Those in the second group are enumerated by sequences of two rooted trees, each with a marked node: Add edges from the first root to the first marked node and from the second root to the first root, and remove the mark from the first marked node. The rooted trees with one marked node are enumerated by sequences (of any length) of rooted trees: Delete the outgoing edges from the nodes on the path from the root to the marked node.
Putting this together, the rooted trees with one marked node are enumerated by the following (sequence A000107),
The connected scalar aromas with one marked node are enumerated by $w(z) + w(z)^2 = t(z)/(1t(z))^2$ (sequence A038002). The scalar aromas with one marked node are enumerated by $a(z)t(z)/(1t(z))^2$ (sequence A027853). Including the unmarked scalar aroma component and the wedge product of rooted trees gives the contribution from type (i) forms as $b_0(u,z)t(z)/(1t(z))^2$ .
For type (ii), combining the components of an unmarked form with one fewer trees (i.e. an element of $\Omega _{n1}$ ), enumerated by $u b_0(u,z)$ , and a marked rooted tree, enumerated by $t(z)/(1t(z))$ , gives the contribution from type (ii) forms as $u b_0(u,z)t(z)/(1t(z))$ .
Summing the results for type (i) and type (ii) forms gives the expression for $b_1(u,z)$ .
The functional forms in $\mathcal {I}_1$ are associated with scalar aromas with one marked leaf (containing the symbol ). They are bijective to the scalar aromas with one marked node of degree one less: Delete the marked leaf and mark the node to which it points (the process is invertible). Therefore, we deduce $ \mathcal {I}_1^N = \Omega _{0,1}^{N1}$ and this gives the result for $c_1(z)$ .
Note that
and that all three spaces are enumerated by the mapping patterns on N (or $N+1$ ) points with one marked node.
Remark 4.3. As row 1 of the augmented bicomplex is exact, its alternating sum of dimensions is zero. This is verified directly:
The following result extends Theorem 4.9 to the divergencefree case.
Theorem 4.4. Let $\widetilde {b}_p(u,z)$ , $\widetilde {c}_p(z)$ , and $\widetilde {s}(z)$ be the divergencefree analogues of the generating functions $b_p(z)$ , $c_p(z)$ and $s(z)$ . Then
Proof. In each case, we need to enumerate the nonselflooped elements of the aromatic bicomplex. The loops occur only in the scalar components. Recall that the nonselflooped scalars, $\widetilde {\Omega }_{0}$ , are enumerated by $\overline {a}(z) = z a(z)/t(z)$ . This gives the result for $\widetilde {b}_0$ and, using Theorem 2.10, $\widetilde {s}(z)$ .
For the second row, we first consider $\widetilde {\Omega }_{0,1}^N$ , the scalar aromas with N nodes, no selfloops and one marked node indicated by . These are enumerated by the selffunctions of $\{1,\dots ,N1\}$ with one marked node. That is, $\widetilde {\Omega }_{0,1}^N\cong \Omega _{0,1}^{N1}$ . The construction is as follows. Consider an element of $\Omega _{0,1}^{N1}$ , that is, a directed graph with $N1$ nodes, each of which has exactly one outgoing edge, with one marked node. Add a new node with an outgoing edge going to the marked node, and redirect all selfloops to the new node; then, move the mark to the new node. This gives a marked nonselflooped scalar aroma of degree N. The process is invertible: Given a marked scalar with no selfloops, redirect the edges that point to the marked node to themselves, mark the node pointed to by the marked node and delete the marked node. Therefore, $\widetilde {\Omega }_{0,1}^N$ has generating function $z a(z)t(z)/(1t(z))^2$ .
For the rest of the second row, $\widetilde {\Omega }_{n,1}$ , recall the two types of forms, type (i) (scalar is marked) and type (ii) (tree is marked). For type (i), that the only change in the divergencefree case is that the scalars must be nonselflooped, as just enumerated. For type (ii), we combine the three components of nonselflooped scalars (enumerated by $a(z)$ ), $n1$ unmarked rooted trees (enumerated by $u b_0(u,z)/a(z)$ ) and one marked rooted tree (enumerated by $t(z)/(1t(z))$ ). The product of these three, plus the contribution from the forms of type (i), gives the result for $b_1(u,z)$ .
The elements of $\widetilde {\mathcal {I}}_1^N=I(\widetilde {\Omega }_{0,p}^N)$ are linear combinations of scalar aromas with N nodes, one marked leaf and no selfloops. As in the general case, these are bijective to the scalar aromas with one marked node of degree one less: Delete the marked leaf and mark the node to which it points. This gives the result for $\widetilde {c}_1(z)$ .
Note that we now have five isomorphic spaces,
Each space is enumerated by the selffunctions of $\{1,\dots ,N1\}$ with one marked node (and generating function $z a(z)t(z)/(1t(z))^2$ ) but in a different way in each case.
We remark that the second row of the divergencefree augmented aromatic bicomplex is not exact.
4.3 Bases of the kernel and image of $d_H$ and $d_H^*$
In this subsection, we work specifically with $d_H\colon \Omega _{1}\to \Omega _{0}$ and we describe the image and the kernel of $d_H$ and $d_H^*$ . We recall that $d_H^*\colon \Omega _{0}^*\to \Omega _{1}^*$ is the dual map of $d_H$ . We described the dimension of $\Psi =\operatorname{\mathrm{Ker}} d_H$ in Subsection 4.2. The following result describes the dimensions of $\mathrm{Im}\ d_H$ , ${\mathrm{Ker}}\ d_H^*$ , and $\mathrm{Im}\ d_H^*$ .
Theorem 4.5. The scalar divergences $\mathrm{Im} d_H$ have dimension $\mathring \Omega _{0}$ . The conditions $\operatorname{\mathrm{Ker}} d_H^*$ that a scalar must satisfy to be a divergence have dimension $\overline {\Omega }_{0}$ . The conditions $\mathrm{Im} d_H^*$ that a vector must satisfy to be divergence free have dimension $\mathring{\Omega}_{0}$ .
Proof. Recall the fundamental theorem of linear algebra for a linear map $A\colon V \to W$ :
The space $\mathrm{Im}\ A^*= \mathrm{Ann}({\mathrm{Ker}}\ A)$ is the annihilator of the kernel of A (i.e., the conditions that an element of A must satisfy in order to lie in the kernel) and ${\mathrm{Ker}}\ A^* = \mathrm{Ann}(\mathrm{Im}\ A)$ is the annihilator of the image of A. Choosing $A=d_H$ gives the result.
Remark 4.6. Consider the aromatic tree $\gamma \backslash e\in \mathcal {F} _{1}$ obtained by cutting one edge $e\in E$ of $\gamma \in \mathcal {F} _{0}$ , $m_1(\gamma ,e)$ the coefficient of $\gamma $ in $d_H(\gamma \backslash e)$ , and $m_2(\gamma ,e)$ the number of edges $\widehat {e}\in E$ of $\gamma $ such that $\gamma \backslash \widehat {e}=\gamma \backslash e$ . Then, for $\gamma \in \Omega _{0}$ , $d_H^*\gamma ^*$ satisfies
We now construct bases of $\operatorname{\mathrm{Ker}} d_H$ , $\mathrm{Im} d_H$ , ${\mathrm{Ker}} d_H^*$ , and $\mathrm{Im} d_H^*$ . We start with the basis of the solenoidal forms.
Theorem 4.7. Let $\varphi \colon \Omega _{1}\to \mathring \Omega _{0}$ be defined by attaching a selfloop to the root, extending by linearity. Define any total order on $\mathcal {T} $ , then a basis of $\Psi =\operatorname{\mathrm{Ker}} d_H$ is
Proof. The map $\varphi $ is surjective: given any selflooped scalar, removing one of the selfloops gives a preimage. Therefore, $\operatorname{\mathrm{Ker}}\varphi  = \operatorname{\mathrm{Ker}} d_H$ . Let us first determine $\operatorname{\mathrm{Ker}}\varphi $ . Consider a selflooped scalar whose distinct selflooped connected components are $\varphi (t_1),\dots , \varphi (t_k)$ for trees $t_1<\dots <t_k$ ; it can be written $\phi \varphi (t_1)\dots \varphi (t_k)$ where $\phi $ is a scalar. Its distinct preimages under $\varphi $ are $\phi \varphi (t_1) \dots t_l \dots \varphi (t_k)$ for $l=1,\dots ,k$ . Therefore, the kernel of $\varphi $ restricted to the span of these preimages has dimension $k1$ and we consider
as a basis of $\operatorname{\mathrm{Ker}}\varphi $ . We now map $\operatorname{\mathrm{Ker}}\varphi $ to $\operatorname{\mathrm{Ker}} d_H$ by
extending by linearity. From exactness, the map is surjective. As $\operatorname{\mathrm{Ker}}\varphi  = \operatorname{\mathrm{Ker}} d_H$ , it is an isomorphism.
Remark 4.8. One could wonder whether the set $\{d_H \wedge \gamma , \gamma \in \mathcal {F} _2\}$ is a basis of $\operatorname{\mathrm{Ker}}(d_H)$ . This is not the case in general. For $N=6$ , one finds for instance the following identity
The following result shows that the divergences are a graph over the selflooped scalars.
Theorem 4.9. For $\alpha \in \mathring {\Omega }_{0}$ , let $k(\alpha )$ be the number of selfloops in $\alpha $ , and $\rho (\alpha )$ be the nonselflooped scalar obtained from $\alpha $ as the sum of the redirection of all $1$ loops to other nodes in all possible ways. Then the map
is an isomorphism and generates a basis of $\mathrm{Im} d_H$ .
Proof. Let $\mathring {V}$ be the set of nodes with selfloops of $\alpha $ . From Proposition 3.6, we deduce that for $v\in \mathring {V}$ , $\alpha \mathcal {E} _v \alpha \in \operatorname{\mathrm{Im}}(d_H)$ . If we have two nodes $v,w\in \mathring {V}$ , then $\alpha \mathcal {E} _v \alpha \in \operatorname{\mathrm{Im}}(d_H)$ and $\mathcal {E} _v \alpha \mathcal {E} _w \mathcal {E} _v \alpha \in \operatorname{\mathrm{Im}}(d_H)$ so that $\alpha \mathcal {E} _w \mathcal {E} _v \alpha \in \operatorname{\mathrm{Im}}(d_H)$ . By applying this process iteratively, we find that
As each selflooped scalar appears once in the image, the map is injective. The map is an isomorphism as the domain and codomain have the same dimension.
Remark 4.10. The operation $\rho $ in Theorem 4.9 corresponds to removing all selfloops in $\alpha $ by repeated integration by parts, as illustrated in the following example on elementary differentials:
That is, $f^i_i f^j_jf^i_j f^j_i$ is a divergence. We describe this comparison with integration by parts further in Subsection 4.4.
Corollary 4.11. No nontrivial combination of nonselflooped scalars in $\Omega _0$ is a divergence.
Corollary 4.12. The conditions to be a divergence $\operatorname{\mathrm{Ker}} d_H^*$ are a graph over the dual of the nonselflooped forms in $\Omega _0$ .
In the following theorem, this graph is realized explicitly.
Theorem 4.13. Let $\widehat {E}\subset E$ be a set of edges of the scalar aroma $\beta \in \Omega _0$ , and let $\beta \backslash \widehat {E}$ be $\beta $ with edges $\widehat {E}$ replaced by selfloops. Let $m(\beta ,\widehat {E})$ be the number of ways that redirecting selfloops of $\beta \backslash \widehat {E}$ results in $\beta $ . Let $\pi \colon \Omega _{0}\to \Omega _{0}$ be defined by
Then
is a basis of $\mathrm{Ann}(\mathrm{Im}\ A)$ , the conditions that a scalar must satisfy to be a divergence.
Proof. The construction is directly related to that in Theorem 4.9. The conditions for $(x,y)$ to lie on the graph $\{x,Ax)\colon x\in \mathbb {R}^n\}\subset \mathbb {R} ^n\times \mathbb {R} ^m$ are $yAx=0$ , $A\in \mathbb {R} ^{m\times n}$ . A basis for these conditions is given by the rows of $yAx$ , where Theorem 4.9 gives the columns of A. That is, for each nonselflooped scalar $\beta $ we need to determine the coefficient of $\beta $ in $\rho (\alpha )$ for each selflooped scalar $\alpha $ . This is the expression for $\pi $ : The term $\widehat {E}=\emptyset $ gives $\beta $ , and the terms from nonempty sets of edges $\widehat {E}$ give the $\alpha $ ’s that can give rise to $\beta $ .
Finally, we present a basis of $\mathrm{Im} d_H^*$ . It is quite straightforward, as we can find a suitable subspace on which $d_H^*$ is injective.
Theorem 4.14. The set
is a basis of $\mathrm{Im} d_H^*$ , the conditions that a form in $\Omega _{0}$ must satisfy to be divergence free.
Proof. For any linear map $A\colon V\to W$ , $\langle A^*w^*,v\rangle = \langle w^*,A v\rangle = 0$ for all $w^*\in W^*$ when $v\in \operatorname{\mathrm{Ker}}\ A$ . Thus, divergencefree vectors do satisfy the given conditions. Furthermore, the dimension of the set is correct. It remains to show that the set is linearly independent. This is the same as showing that $d_H^*_{\mathring {\Omega }_{0}^*}$ is injective.
Let $\mathrm{pr}\colon \Omega _{0}\to \mathring \Omega _{0}$ be the natural projection to the selflooped scalars. From Theorem 4.9, the divergences form a graph over the selflooped scalars. That is, $\mathrm{pr}\circ d_H$ is surjective. Therefore, its dual $d_H^*_{\mathring {\Omega }_{0}^*}$ is injective.
4.4 Integration by parts of aromatic forests
The horizontal homotopy operator is often described in the differential geometry literature as an integration by parts operator. The concept of integration by parts of trees was also introduced in the context of stochastic numerical analysis in [Reference Laurent and Vilmart38, Reference Laurent and Vilmart39] on exotic aromatic Bseries (see also [Reference Bronasco8]). We show in this section that a similar integration by parts process can be adapted in the context of aromatic forms to define a different horizontal homotopy operator on $\Omega _{0,p}$ .
Let $\gamma \in \Omega _{0,p}^N$ a linear combination of forests; let $\tau \in \mathcal {F} _{0,p}^N$ one of these forests and v a vertex of $\tau $ on a 1loop. We denote $a(\tau )$ the coefficient of $\tau $ in $\gamma $ , and $\theta _v(\tau )$ the forest $\tau $ where we remove the edge linking v to itself and transform v into a root. The alternative horizontal homotopy operator $\widehat {h}_H$ on $\Omega _{0,p}$ is given by the following algorithm on $\Omega _{0,p}^N$ and extended to $\Omega _{0,p}$ by linearity.
Note that each iteration in the algorithm reduces the number of 1loops by one. Thus, the algorithm always terminates. We emphasize that the result of the algorithm is independent of the order in which we detach the $1$ loops. This is not the case in the similar algorithm proposed in [Reference Laurent and Vilmart38], as there is an extra term involved in the integration by parts process.
Theorem 4.15. For $\gamma \in \Omega _{0,p}$ , the output $\widehat {h}_H\gamma $ of the algorithm is the horizontal homotopy operator, up to a divergencefree term, that is,
Proof. After the algorithm terminates, $\widehat {\gamma }$ is given by
where $\widehat {\gamma }$ does not contain any 1loop and where we used Theorem 3.8. We deduce from Proposition 3.6 that $\widehat {\gamma }\in \operatorname{\mathrm{Im}}(d_H)$ . According to Corollary 4.11, we find $\widehat {\gamma }=0$ .
We now have two different ways to compute the horizontal homotopy operator on $\Omega _{0,p}$ . The first one, presented in Subsection 3.1, uses the Euler operators. The second one, in the spirit of [Reference Laurent and Vilmart38], is based on the repeated use of detaching and grafting operations on specific nodes. We emphasize that the expressions of the homotopy operator given by these two methods are different in general but are always equal up to a divergencefree term. The two homotopy operators can produce both concise and tedious outputs, and they outperform each other in this manner on different forests. We refer the reader to Table 3 for some examples. This difference in the number of terms increases rapidly with the order; for instance, for the form , $h_H \gamma $ has 6 terms, while $\widehat {h}_H \gamma $ has 26 terms. On the other hand, it is possible to find examples where $\widehat {h}_H~$ produces fewer terms than $h_H$ . It would be interesting to find a homotopy operator with a minimal number of terms in the output or a procedure to simplify the outputs of a homotopy operator in the spirit of [Reference Anderson1, Sect. IV.B]. Moreover, it is not known whether a similar approach in the divergencefree context could yield a different homotopy operator. This is matter for future work.
4.5 Explicit description of volumepreserving aromatic integrators
It is known that the only volumepreserving consistent Bseries method is the exact flow [Reference Chartier and Murua18, Reference Iserles, Quispel and Tse32]. In [Reference MuntheKaas and Verdier47], the question of the existence of a volumepreserving aromatic method is raised, where an aromatic method is a onestep integrator that has an expansion as an aromatic Bseries. In [Reference Bogfjellmo6], a methodology to create pseudovolumepreserving integrators is proposed by substituting in a standard Runge–Kutta method the vector field f with an aromatic Bseries. We give in this subsection an explicit expression of the form of a volumepreserving aromatic method, and we use it to prove that there does not exist any aromatic Runge–Kutta integrator and to discuss the form of a volumepreserving aromatic Bseries method.
Consider a consistent onestep integrator (1.2) for solving the differential equation (1.1) with the assumption $\operatorname{\mathrm{div}}(f)=0$ . We assume the integrator (1.2) has an expansion in aromatic Bseries given by the linear coefficient map $a\colon \widetilde {\Omega }_1 \rightarrow \mathbb {R} $ ; that is, its (formal) Taylor expansion has the form
where $\sigma (\tau )$ is the cardinal of the set of automorphisms on the set of nodes V of $\tau $ that leave $\tau $ unchanged (see [Reference Bogfjellmo6]). We call such a method an aromatic Bseries method. If in addition the integrator reduces to a standard Runge–Kutta method when choosing a vector field f that satisfies $F(\gamma )(f)=0$ for all $\gamma \in \mathcal {F} _1\setminus \mathcal {T} $ , we call the integrator an aromatic Runge–Kutta method.
An aromatic Bseries method can be seen as the exact solution of the modified ODE (1.3), and the modified flow is given by the aromatic Bseries $B(b)$ satisfying $B(a)=B(b)\triangleright B(e)$ . The operation $\triangleright $ is the substitution of Bseries and $b\colon \widetilde {\Omega }_1\rightarrow \mathbb {R} $ is the coefficient map of the modified flow. It is known that a and b satisfy $b \star e=a$ , where $\star $ is the substitution of Bseries coefficients [Reference Calaque, EbrahimiFard and Manchon13, Reference Chartier, Hairer and Vilmart17, Reference Bogfjellmo6]. The map e is the coefficient of the exact flow of equation (1.1). Its expression is given for instance in [Reference Hairer, Lubich and Wanner23, Chap. III] for Butcher trees. It is extended to the aromatic trees by $e(\gamma )=0$ if $\gamma $ is composed of at least an aroma. The question raised in [Reference MuntheKaas and Verdier47] is the following: Can we find an aromatic Bseries method such that $d_H B(b)=0$ , that is, such that the modified Bseries $B(b)$ is solenoidal? Note that choosing $a=e$ yields a simple solution to the problem, but there does not exist any reasonable numerical method whose coefficient map is given by the exact flow coefficient e.
The main motivation for considering aromatic Bseries instead of standard Bseries comes from the following negative result.
Theorem 4.16 ([Reference Chartier and Murua18, Reference Iserles, Quispel and Tse32])
The solenoidal combinations of rooted trees satisfy
In particular, the only volumepreserving consistent Bseries method is the exact flow.
Note that in the standard context, Theorem 4.16 is a consequence of Theorem 4.9. Indeed, let $v=\sum _i c_i t_i$ , $t_i\in \mathcal {T}$ , be a combination of rooted trees. As the divergences are graphs over the selflooped scalars (see Theorem 4.9), $d_H v=0$ if the coefficient of each selflooped scalar in $d_H v$ is zero. But $\mathrm{pr} d_H v=\sum _i c_i \theta (t_i)$ , where $\mathrm{pr}\colon \Omega _{0}\to \mathring \Omega _{0}$ is the natural projection to the selflooped scalars, and the selflooped scalars $\theta (t_i)$ are linearly independent.
Following Theorem 4.16, one is interested in finding a class of nontrivial volumepreserving consistent aromatic Bseries methods. We deduce from the previous discussion and Theorems 2.9, 2.10, and 2.11 the following explicit description of the coefficients of an aromatic volumepreserving integrator.
Theorem 4.17. If $B(a)$ is the aromatic Bseries of a consistent volumepreserving integrator, then there exists $\eta \in \widetilde {\Omega }_2$ such that the modified flow is a Bseries of the form
and $B(a)$ is given by the substitution
More precisely, there exists a coefficient map $\alpha \colon \widetilde {\Omega }_2\rightarrow \mathbb {R} $ such that a is given by
For the first orders, the Bseries of a volumepreserving aromatic Bseries method has the form:
Note that the coefficients of the bamboo trees (or tall trees)
coincide with the ones of the exact flow. This fact has been noticed for standard Bseries in particular in [Reference Kang and Shang34] (see also [Reference Hairer, Lubich and Wanner23, Lemma IV.3.2]). We deduce from this observation the following result.
Theorem 4.18. An aromatic Runge–Kutta method cannot be volume preserving.
Proof. The only bamboo tree that appears in solenoidal forms is . Indeed, $d_H^*$ vanishes on $\mathcal {B} \mathcal {T} $ as the only solenoidal forms where a bamboo tree can appear are of the form $d_H\wedge \tau _1 \tau _2$ , where $\tau _1$ , $\tau _2\in \mathcal {B} \mathcal {T} $ , and no bamboo tree appears in these forms. According to the identity (4.6), the Bseries of an aromatic volumepreserving method has to coincide with the exact flow $B(e)$ on the bamboo trees $\mathcal {B} \mathcal {T} $ . Such a method solves exactly linear problems and it is known (see [Reference Hairer and Wanner25, Chap. IV.3]) that Runge–Kutta integrators do not solve linear problems exactly in general. As the aromatic forests represent different elementary differentials (see [Reference Iserles, Quispel and Tse32]), any aromatic integrator that reduces to a standard Runge–Kutta integrator when sending the aromas to zero cannot be volume preserving.
The methodology proposed in [Reference Bogfjellmo6, Sect. 7] to obtain volumepreservation of high order and the approach in [Reference MuntheKaas and Verdier47, Sect. 9] give classes of aromatic integrators that can preserve volume up to a high order but that cannot be volume preserving. To build an aromatic volumepreserving method, it is fundamental to start with an ansatz that is exact on bamboo trees (that is, the method is exact for linear problems). A natural guess is to consider aromatic integrators that reduce to exponential Rosenbrock integrators (see, for instance, [Reference Berland, Owren and Skaflestad4, Reference Hochbruck and Ostermann31, Reference Luan and Ostermann41]) when sending the aromas to zero. This calls for future works that study the substitution law and the variational bicomplex directly on aromatic exponential Bseries, in order to find a volumepreserving aromatic Bseries method.
5 Conclusion and future work
In this work, we introduced a new algebraic object, called the aromatic bicomplex, for the study of aromatic forms. We studied the exactness of the bicomplex in the standard case and in the divergencefree case. To this end, we introduced the Euler operators and the homotopy operators, as well as an augmented bicomplex. The algebraic properties we proved have concrete consequences on the numerical analysis of volumepreserving integrators. They allow to describe completely the solenoidal forms and the Bseries of an aromatic volumepreserving method. In particular, we proved that there are no volumepreserving aromatic Runge–Kutta methods.
Many theoretical and applied questions arise from the present work. Following the results of Subsection 4.5, it would be interesting to rewrite the substitution and divergence operations in the context of exponential Bseries, in order to find an aromatic exponential volumepreserving method.
The integration by parts of (exotic) aromatic forests is a new operation that has applications in stochastic numerical analysis and in the study of volumepreserving integrators. To the best of our knowledge, few works study the structure of (exotic) aromatic forests equipped with the integration by parts process. In particular, there is no explicit expression for the output of the integration by parts process in the stochastic setting [Reference Laurent and Vilmart38]. An exact formula would greatly benefit the creation of highorder methods for solving ergodic stochastic differential equations.
There is a considerable literature on the variational bicomplex and the De Rham cohomology (see [Reference Anderson1] and references therein). It would be interesting to generalise some of the existing results in the context of aromatic forms. For instance, one could try to find simpler expressions for the homotopy operators (see Subsection 4.4) to find an augmented bicomplex in the divergencefree case (see Subsection 4.1). Two major applications of the variational bicomplex are Noether’s theorems and the study of the Laplace–De Rham operator $\Delta =d_H d_H^* +d_H^* d_H.$ It would be interesting to see how these results translate to aromatic forms.
A First solenoidal forms
We write the generators of the solenoidal forms $\widetilde {\Psi }^N$ in the divergencefree case for the first orders in Table 4. As a consequence of Theorem 2.11, we find all the generators by computing $d_H\gamma $ for $\gamma \in \widetilde {\mathcal {F} }_2^N$ and by adding the trivial tree . For $N\leq 5$ , we observe that they form a basis of the solenoidal forms (see Remark 4.8). Note how no bamboo trees appear in the solenoidal forms, as discussed in Subsection 4.5.
B The aromatic bicomplex for the first orders
We present in Figures 4 and 5 the augmented aromatic bicomplex for $N=1$ , $2$ , $3$ in the standard case, as defined in Subsection 2.2. The divergencefree aromatic bicomplex is deduced from it by deleting the $1$ loops and the extra column on the right. We give a basis of each space, and we omit for conciseness the trivial spaces surrounding the bicomplex and the wedge $\wedge $ when writing the aromatic forms in the diagrams. Note that the alternate sum of dimensions in each horizontal and vertical sequence, and in the Euler–Lagrange complex adds up to zero, as a consequence of Theorem 2.9 and Theorem 4.1.
Competing interest
The authors have no competing interest to declare.
Financial Support
The work of Adrien Laurent was supported by the Research Council of Norway through project 302831 ‘Computational Dynamics and Stochastics on Manifolds’ (CODYSMA). Robert McLachlan acknowledges the support of the Simons Foundation and the hospitality of the Isaac Newton Institute for the Mathematical Sciences through their program ‘Geometry, compatibility and structure preservation in computational differential equations’, where part of this work was conducted.