# The aromatic bicomplex for the description of divergence-free aromatic forms and volume-preserving integrators

Adrien Laurent, Robert I. Mclachlan, Hans Z. Munthe-Kaas, Olivier Verdier

► To cite this version:

Adrien Laurent, Robert I. Mclachlan, Hans Z. Munthe-Kaas, Olivier Verdier. The aromatic bicomplex for the description of divergence-free aromatic forms and volume-preserving integrators. Forum of Mathematics, Sigma, 2023, 11, pp.E69. 10.1017/fms.2023.63 . hal-04347717

HAL Id: hal-04347717

https://hal.science/hal-04347717v1

Submitted on 15 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# The aromatic bicomplex for the description of divergence-free aromatic forms and volume-preserving integrators

Adrien Laurent<fnref n="1" />, Robert I. McLachlan<fnref n="2" />, Hans Z. Munthe-Kaas<fnref n="1" />, and Olivier Verdier<fnref n="3" />

February 6, 2023

## Abstract

Aromatic B-series were introduced as an extension of standard Butcher-series for the study of volume-preserving integrators. It was proven with their help that the only volume-preserving B-series method is the exact flow of the differential equation. The question was raised whether there exists a volume-preserving integrator that can be expanded as an aromatic B-series. In this work, we introduce a new algebraic tool, called the aromatic bicomplex, similar to the variational bicomplex in variational calculus. We prove the exactness of this bicomplex and use it to describe explicitly the key object in the study of volume-preserving integrators: the aromatic forms of vanishing divergence. The analysis provides us with a handful of new tools to study aromatic B-series, gives insights on the process of integration by parts of trees, and allows to describe explicitly the aromatic B-series of a volume-preserving integrator. In particular, we conclude that an aromatic Runge-Kutta method cannot preserve volume.

**Keywords:** geometric numerical integration, volume-preservation, aromatic B-series, divergence, aromatic bicomplex, aromatic forms, solenoidal forms, Euler operator, homotopy operator, integration by parts, Euler-Lagrange complex.

AMS subject classification (2020): 65L06, 41A58, 58J10, 58A12, 37M15, 05C05.

## 1 Introduction

Let $f: \mathbb{R}^d \to \mathbb{R}^d$ be a smooth Lipschitz vector field and let $y: [0, T] \to \mathbb{R}^d$ be the solution to the ordinary differential equation

$$y'(t) = f(y(t)), \quad t \in [0, T[, \quad y(0) = y_0. \qquad (1.1)$$

If the vector field $f$ is divergence-free, that is, if $\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)$. Divergence-free vector fields appear in a variety of concrete dynamical systems, for instance in fluid dynamics, meteorology or molecular dynamics. To integrate such

<footnote n="1">¹Department of Mathematics, University of Bergen, Bergen, Norway.
adrien.laurent@uib.no. hans.munthe-kaas@uib.no.</footnote>

<footnote n="2">²Institute of Fundamental Sciences, Massey University, Palmerston North, New Zealand.
r.mclachlan@massey.ac.nz.</footnote>

<footnote n="3">³Department of Computing, Electrical Engineering and Mathematical Sciences, Western Norway University of Applied Sciences, Bergen, Norway.
olivier.verdier@hvl.no.</footnote>
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 time stepsize, and we choose a one-step integrator

$$y_{n+1} = \Phi(y_n, h). \qquad (1.2)$$

For large classes of integrators, the standard methodology to study volume-preservation uses backward error analysis (see, for instance, the textbook [23]). The integrator (1.2) can be interpreted formally as the exact solution of a modified ODE

$$\tilde{y}'(t) = \tilde{f}(\tilde{y}(t)), \qquad (1.3)$$

where the modified vector field $\tilde{f}$ typically depends on $f$ and its derivatives. Then, the integrator is volume-preserving if and only if $\text{div}(\tilde{f}) = 0$.

There is a considerable literature on volume-preserving methods, with many applications for solving a variety of dynamical systems. The existing volume-preserving integrators rely either on a specific form of the vector field $f$ [23, 18], on splitting methods [33, 43, 42, 61], or on generating functions [53, 55, 51]. For quadratic differential equations, we mention the works [49, 16, 15, 14, 7] that study the Kahan-Hirota-Kimura discretization [32, 28, 29] 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 non-reversible problems [5]. 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 volume-preserving 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 volume-preserving aromatic B-series method.

Introduced in [10, 24] (see also the textbooks [23, 11, 12] and the review [41]), the Butcher-series formalism is an important tool in numerical analysis. Originally used for the calculation of order conditions for Runge-Kutta methods, the use of B-series was quickly extended to a variety of applications, in particular in geometric numerical integration [23], and more recently in the approximation of stochastic evolutionary problems [9, 34, 52, 19, 36] or in the theory of rough paths [22, 25]. For a large class of integrators, such as Runge-Kutta methods, the modified vector field $\tilde{f}$ in (1.3) is a B-series in $f$. An extension of B-series, called aromatic B-series, was introduced in [18, 31] for the study of volume-preserving integrators. They allow to compute the divergence of a B-series, and were later studied in [40, 45, 6, 21] for their algebraic and geometric properties. In [18, 31], it is showed that no non-trivial B-series is divergence-free, so that the only volume-preserving B-series method is the exact flow. In particular, no Runge-Kutta method can preserve volume exactly (see also [33]). However, the space of divergence-free aromatic B-series is infinite-dimensional, the simplest non-trivial example being:

$$\sum_{i,j,k=1}^{d} (f_k^j f_j^k f^i + f_{jk}^j f^k f^i - f_j^j f_k^i f^k - f_{jk}^i f^j f^k) \partial_i \in \text{Ker}(\text{div}). \qquad (1.4)$$

In [45], the question whether there exists a volume-preserving integrator that has an expansion as an aromatic B-series is raised. The present article gives a handful of tools to answer this question and describes explicitly the aromatic B-series of a volume-preserving integrator. In particular, we prove that no aromatic Runge-Kutta method can preserve-volume exactly.
When using B-series, 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 [38]. For instance, in dimension 3, denote $\Omega_n^{(\text{diff})}(\mathbb{R}^3)$ the space of smooth differential $n$-forms on $\mathbb{R}^3$, then the De Rham complex writes

$$ \Omega_0^{(\text{diff})}(\mathbb{R}^3) \xrightarrow{d} \Omega_1^{(\text{diff})}(\mathbb{R}^3) \xrightarrow{d} \Omega_2^{(\text{diff})}(\mathbb{R}^3) \xrightarrow{d} \Omega_3^{(\text{diff})}(\mathbb{R}^3), \quad (1.5) $$

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 divergence-free 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 [58, 57, 56, 59, 60] as a natural and general development of the variational chain. It has a variety of applications in the areas of differential geometry and topology, differential equations, mathematical physics and PDEs (see, for instance, the textbooks [1, 47], the introductory article [2], 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 $\text{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 $\text{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 [36, 37]. The main application of this work is the explicit description of the aromatic B-series of a volume-preserving 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 volume-preserving aromatic B-series methods.

The paper is organised as follows. Section 2 introduces the aromatic forests and forms, 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 divergence-free 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 volume-preserving integrators to obtain an explicit description of the B-series of an aromatic volume-preserving integrator.
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 [2] for further details (see also the textbooks [47, 35, 38]). The infinite jet bundle $J^\infty(\mathbb{R}^d)$ is used to describe coordinate-free Taylor expansions. It is the infinite-dimensional vector bundle over $\mathbb{R}^d$ where elements of $J_x^\infty(\mathbb{R}^d)$ have coordinates

$$ (f^i(x), f_j^i(x), f_{jk}^i(x), \dots), \quad f^i: \mathbb{R}^d \to \mathbb{R}, \quad f_{j_1, \dots, j_m}^i = \frac{\partial f^i}{\partial x^{j_1} \dots \partial x^{j_m}}. \qquad (2.2) $$

We consider differential forms on $J^\infty(\mathbb{R}^d)$ and define the following contact forms

$$ \theta_{j_1, \dots, j_p}^i = df_{j_1, \dots, j_p}^i - f_{j_1, \dots, j_p, k}^i dx^k. $$

The set $\{dx^i, \theta^i, \theta_j^i, \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

$$ g(f)dx^{i_1} \wedge \dots \wedge dx^{i_n} \wedge \theta_{I_1}^{j_1} \wedge \dots \wedge \theta_{I_p}^{j_p}, $$

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_{n,p}^{(\text{diff})}(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: \mathbb{R}^d \to \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:

$$ F(\gamma)(f) = \sum_{i_w \in \{1, \dots, d\}} \prod_{w \in V \setminus R} f_{I_{\Pi(v)}}^{i_w} \partial_{i_{r_1}} \otimes \dots \otimes \partial_{i_{r_n}} \otimes \theta_{I_{\Pi(\circled{1})}}^{i_1} \otimes \dots \otimes \theta_{I_{\Pi(\circled{p})}}^{i_p}, $$

where $\Pi(v)$ is the set of predecessors of $v \in V$. We extend $F$ on $\text{Span}(\mathcal{F}_{n,p})$ by linearity.

For instance, we find

$$ F(\bullet)(f) = \text{div}(f) = \sum_{i=1}^{d} f_i^i, \quad F(\bullet)(f) = \sum_{i,j=1}^{d} f_j^i f^j \partial_i, \quad F(\circled 1)(f) = \sum_{i=1}^{d} \partial_i \otimes \theta^i. $$

The forest $\gamma$ in (2.1) represents the elementary differential

$$ F(\gamma)(f) = \left( \sum_{j,k,l=1}^{d} f_j^k f_k^l f_l^j \right) \left( \sum_{j,k=1}^{d} f_{jk}^j f^k \right) \sum_{i_{r_1}, i_{r_2}, i_{r_3}, j, k=1}^{d} f_j^{i_{r_1}} f^k f^{i_{r_2}} \partial_{i_{r_1}} \otimes \partial_{i_{r_2}} \otimes \partial_{i_{r_3}} \otimes \theta^{i_{r_3}} \otimes \theta_k^j. $$

In the following, we work with specific linear combinations of aromatic forests in $\text{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 $S_n^\bullet$ (resp. $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

$$ \wedge^\bullet \gamma = \frac{1}{n!} \sum_{\sigma \in S_n^\bullet} \varepsilon(\sigma) \sigma \gamma, $$
where $\varepsilon(\sigma)$ is the signature of the permutation $\sigma$. Similarly, the covertices wedge is

$$ \wedge^{\circ} \gamma = \frac{1}{p!} \sum_{\sigma \in S_p^{\circ}} \varepsilon(\sigma) \sigma \gamma, $$

and the total wedge is $\wedge = \wedge^{\bullet} \wedge^{\circ} = \wedge^{\circ} \wedge^{\bullet}$.

We extend the wedge operations to $\mathrm{Span}(\mathcal{F}_{n,p})$ by linearity and we denote the set of aromatic forms $\Omega_{n,p} = \wedge \mathrm{Span}(\mathcal{F}_{n,p})$, respectively $\Omega_{n,p}^N = \wedge \mathrm{Span}(\mathcal{F}_{n,p}^N)$, and $\Omega_n = \Omega_{n,0}$. As $\wedge^2 = \wedge$, the operator $\wedge: \mathrm{Span}(\mathcal{F}_{n,p}) \to \Omega_{n,p}$ is a projection on $\Omega_{n,p}$.

**Example.** Let $\gamma = \bullet \dot{\bullet} \in \mathcal{F}_2$, then

$$ \wedge \gamma = \frac{1}{2} (\bullet \dot{\bullet} - \dot{\bullet} \bullet) \in \Omega_2. $$

For $\gamma$ given by (2.1), the associated aromatic form is

$$ \begin{aligned} \wedge \gamma = \frac{1}{12} & (\bullet \circ \bullet \circ \bullet \circ \bullet \circ \bullet \circ \bullet + \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ + \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ - \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ) \\ & - \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ - \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ - \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ - \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ \\ & - \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ + \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ + \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ + \bullet \circ \circ \circ \circ \circ \circ \circ \circ \circ \circ). \end{aligned} $$

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 \to 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 \to u}\gamma$. Let $\gamma \in \mathcal{F}_{n,p}$ and $v \in V^*$, we define $\gamma_{v \to \mathcal{k}}$ as the forest obtained by replacing the node $v$ by a new covertex $\mathcal{k}$. Similarly, $\gamma_{\mathcal{k} \to \bullet}$ is the forest obtained by replacing the covertex $\mathcal{k}$ by a vertex.

**Example.** Let $\gamma = \bullet^{\circ}$ and $r$ its root, then

$$ D^r \gamma = \bullet^{\circ} + \bullet^{\circ}, \quad \gamma_{r \to \bullet} = \bullet^{\circ}, \quad \gamma_{\bullet \to \bullet} = \bullet^{\circ}. $$

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

$$ d_H \gamma = D^{r_n} \gamma, \quad d_V \gamma = \wedge \sum_{v \in V^*} \gamma_{v \to \overline{(p+1)}}. $$

We extend $d_H$ and $d_V$ on $\Omega_{n,p}$ by linearity into $d_H: \Omega_{n,p} \to \Omega_{n-1,p}$ and $d_V: \Omega_{n,p} \to \Omega_{n,p+1}$. The aromatic forms in $\mathrm{Ker}(d_H)$ are called solenoidal forms, $\Psi = \mathrm{Ker}(d_H|_{\Omega_1})$ are the solenoidal combinations of trees, and $\Psi^N = \mathrm{Ker}(d_H|_{\Omega_1^N})$ are the solenoidal combinations of trees of order $N$.

The operator $d_H: \Omega_1 \to \Omega_0$ is often called the divergence of an aromatic tree, as for $\gamma \in \Omega_1$, $d_H$ satisfies $\mathrm{div}(g) = F(d_H(\gamma))(f)$, where $F(\gamma)(f) = \sum_i g^i \partial_i$. The operator $d_H: \Omega_1 \to \Omega_0$ is often called the divergence of an aromatic tree.
**Example.** Consider $\gamma_1 = \bullet \in \Omega_1$, $\gamma_2 = \wedge \bullet \dot{\bullet} \in \Omega_2$, and $\gamma_3 = \wedge \bullet \textcircled{1} \in \Omega_{2,1}$, then

$$d_H \gamma_1 = \circled{0}, \quad d_H \gamma_2 = \frac{1}{2} (\circled{0} \bullet + \circled{0} \dot{\bullet} \bullet - \circled{0} \circled{0} - \circled{0} \circled{0} \dot{\circ}), \quad d_H \gamma_3 = \frac{1}{2} (\circled{0} \circled{0} \bullet + \textcircled{1} \circled{0} \bullet - \circled{0} \textcircled{1} \circled{0} - \textcircled{1} \circled{0}), \\ d_V \gamma_1 = \textcircled{1}, \quad d_V \gamma_2 = \wedge \textcircled{1} \dot{\circ} + \wedge \textcircled{1} \text{\textcircled{1}} + \wedge \textcircled{1} \text{\textcircled{1}}, \quad d_V \gamma_3 = \wedge \textcircled{2} \text{\textcircled{1}} = \frac{1}{2} (\textcircled{2} \text{\textcircled{1}} - \textcircled{1} \text{\textcircled{2}}).$$

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

$$d_H^2 = 0, \quad d_V^2 = 0, \quad d_V d_H = d_H d_V.$$

**Remark 2.7.** The horizontal and vertical derivatives were already defined on the set of aromatic trees $\Omega_1 = \text{Span}(\mathcal{F}_1)$ respectively in [18, 31] for $d_H$ and in [21] for $d_V$. In this last work, the trace operator $\text{Tr}: \Omega_{1,1} \to \Omega_0$ on aromatic trees is studied. For $\gamma \in \Omega_{1,1}$, it is given with our notations by

$$\text{Tr} \gamma = (D^{r \to \textcircled{1}} \gamma)_{\textcircled{1} \to \bullet},$$

and it makes the following diagram commute.

<figure><img src="image_2.png" /></figure>

## 2.2 The aromatic bicomplex: exactness and description of solenoidal forms

The variational bicomplex is a powerful tool of variational calculus [2]. 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 divergence-free 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 $\text{Im}(d_H|_{\Omega_{n,p}})$, as we will see in Subsection 4.1. We refer the reader to the appendix B 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 [2] in the following way. Let $\sigma$ be the standard duality map between the type $n$-alternating contravariant tensors and the type $(d-n)$-differential forms with respect to the standard Euclidean coordinates and volume form $dx^1 \wedge \dots dx^d$ (see, for instance, [38]). Let $f: \mathbb{R}^d \to \mathbb{R}^d$, the map $\sigma F(\cdot)(f)$ sends an aromatic form in $\Omega_{n,p}$ to a differential form in $\Omega_{d-n,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

$$d_H^{(\text{diff})} \sigma F(\gamma)(f) = \sigma F(d_H \gamma)(f), \quad d_V^{(\text{diff})} F(\gamma)(f) = (-1)^{n+p} F(d_V \gamma)(f).$$
<figure><img src="image_3.png" /><figcaption>Figure 1: The aromatic bicomplex (left) and its subcomplex of order N (right).</figcaption></figure>

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(\cdot)(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 \ge 0$,

$$ \mathrm{Im}(d_H|_{\Omega_{n+1,p}^N}) = \mathrm{Ker}(d_H|_{\Omega_{n,p}^N}), \quad \mathrm{Im}(d_V|_{\Omega_{n,p}^N}) = \mathrm{Ker}(d_V|_{\Omega_{n,p+1}^N}). $$

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

$$ 2d_H \wedge \bullet \cdot \bullet = \bullet \cdot \bullet + \circ \cdot \bullet - \circ \circ \bullet - \circ \bullet \circ, \qquad (2.3) $$

and a basis of $\Psi^4$ is given by the forests

$$ 2d_H \wedge \bullet \bullet \bullet = \circ \circ \bullet + \circ \circ \circ + \circ \circ \bullet - \circ \circ \circ - \circ \circ \bullet - \circ \circ \circ - \circ \circ \bullet, $$

$$ 2d_H \wedge \bullet \circ \circ = \circ \circ \bullet + 2 \circ \circ \circ + \circ \circ \circ - 2 \circ \circ \circ - \circ \circ \circ - \circ \circ \circ, $$

$$ 2d_H \wedge \circ \circ \bullet \circ = \circ \circ \bullet + \circ \circ \circ \circ + \circ \circ \circ \circ - \circ \circ \circ \circ - \circ \circ \circ \circ - \circ \circ \circ \circ - \circ \circ \circ \circ. $$

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 volume-preserving integrators for solving (1.1), it is fundamental to assume that the vector field $f$ satisfies $\mathrm{div}(f) = 0$. With aromatic forms, it amounts to sending all forests
containing an aroma with a 1-loop to 0. We consider the vector space $\mathcal{A}$ spanned by aromatic forests containing at least one 1-loop, $\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 $\tilde{\mathcal{F}}_{n,p}$ the set of aromatic forests in $\mathcal{F}_{n,p}$ without 1-loops, $\tilde{\Omega}_{n,p} = \Omega_{n,p}/\mathcal{A}_{n,p}$, $\tilde{\Omega}_{n,p}^N = \Omega_{n,p}^N/\mathcal{A}_{n,p}^N$, $\tilde{\Omega}_n = \tilde{\Omega}_{n,0}$, $\tilde{\Psi} = \text{Ker}(d_H|_{\tilde{\Omega}_1})$, and $\tilde{\Psi}^N = \text{Ker}(d_H|_{\tilde{\Omega}_1^N})$. The divergence-free aromatic bicomplex with $N$ nodes is drawn in Figure 2. In the simplest case $N=1$ (see Figure 2), the aromatic bicomplex

<figure><img src="image_4.png" /><figcaption>Figure 2: The divergence-free aromatic bicomplex of order $N$ (left) and of order $N=1$ (right).</figcaption></figure>

is not exact. Indeed, we have $\bullet \in \text{Ker}(d_H)$, but $\bullet \notin \text{Im}(d_H)$. One of the main results of this paper is that the case $N=1$ is the only case where the divergence-free aromatic bicomplex is not exact.

**Theorem 2.10.** The divergence-free 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 $\text{div}(f) = 0$ does not create new non-trivial solenoidal forms.

**Theorem 2.11.** For $N \neq 1$, and all $n \ge 1$, $p \ge 0$, the kernel of the divergence operator $d_H$ satisfies

$$ \text{Ker}(d_H|_{\tilde{\Omega}_{n,p}^N}) = \text{Ker}(d_H|_{\Omega_{n,p}^N})/\mathcal{A}_{n,p}^N, $$

that is, the solenoidal aromatic forms in the divergence-free context exactly correspond to the solenoidal aromatic forms in the standard context, except for the forms spanned by $\bullet$ and $\textcircled{1}$.

*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

$$ \text{Im}(d_H|_{\tilde{\Omega}_{n+1,p}^N}) = \text{Im}(d_H|_{\Omega_{n+1,p}^N})/\mathcal{A}_{n,p}^N. $$

Theorem 2.9 and Theorem 2.10 yield the desired identity. $\square$

A generating set of the solenoidal forms $\tilde{\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 divergence-free context are given by

$$ \tilde{\Psi}^3 = \text{Span}(\bullet \cdot - \bullet \bullet). $$
We give generators of the solenoidal forms $\tilde{\Psi}^N$ for the first orders in Appendix A. This explicit description of solenoidal forms is especially useful in the numerical study of volume-preserving integrators, as discussed in Subsection 4.5.

We enumerate in Subsection 4.2 the dimensions of the $\Omega_{n,p}^N$ and $\tilde{\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 $\tilde{\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 1-loops $\tilde{\Omega}_0^N$ (see Table 1 for examples).

<table><caption>Table 1: Dimensions of the space of solenoidal forms in both contexts for the first orders N (see Theorems 4.2 and 4.4). Note how $|\Psi^N| = |\Omega_1^N| - |\tilde{\Omega}_0^N|$.</caption><thead><tr><th>N</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th></tr></thead><tbody><tr><td>|Ω<sub>1</sub><sup>N</sup>|</td><td>1</td><td>2</td><td>6</td><td>16</td><td>45</td><td>121</td><td>338</td><td>929</td><td>2598</td><td>7261</td><td>20453</td><td>57738</td><td>163799</td><td>465778</td></tr><tr><td>|Ω̂<sub>0</sub><sup>N</sup>|</td><td>1</td><td>2</td><td>5</td><td>13</td><td>34</td><td>90</td><td>243</td><td>660</td><td>1818</td><td>5045</td><td>14102</td><td>39639</td><td>111982</td><td>317533</td></tr><tr><td>|Ψ<sup>N</sup>|</td><td>0</td><td>0</td><td>1</td><td>3</td><td>11</td><td>31</td><td>95</td><td>269</td><td>780</td><td>2216</td><td>6351</td><td>18099</td><td>51817</td><td>148245</td></tr><tr><td>|Ψ̃<sup>N</sup>|</td><td>1</td><td>0</td><td>1</td><td>2</td><td>7</td><td>16</td><td>48</td><td>123</td><td>346</td><td>937</td><td>2626</td><td>7284</td><td>20533</td><td>57804</td></tr></tbody></table>

# 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 divergence-free 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 \to u}\gamma \in \mathcal{F}_{n-1,p}^{\diamond}$ is the marked aromatic forest obtained by adding the edge linking r to u to $\gamma$. The set of detached roots of $D^{r \to u}\gamma$ is $R_0$ if $r \notin R_0$ and $R_0 \setminus \{r\}$ else. We define

$$D^r \gamma = \sum_{u \in V \setminus \{v_{\diamond}\}} D^{r \to u} \gamma.$$

Let $q$ be a non-negative 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 \to u}$ by

$$D^I \gamma = \sum_{\varphi: I \to V} \prod_{w \in I} D^{w \to \varphi(w)} \gamma, \quad D^{I \to u} \gamma = \prod_{w \in I} D^{w \to u} \gamma,$$
and $D^q$ and $D^{q \to u}$ by

$$D^q \gamma = \sum_{I \subset R_0, |I|=q} D^I \gamma, \quad D^{q \to u} \gamma = \sum_{I \subset R_0, |I|=q} D^{I \to u} \gamma.$$

Let $\gamma \in \mathcal{F}_{n,p}$ and $v \in V$, then $\gamma_{v^\circ} \in \mathcal{F}_{n+|\Pi(v)|,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_\diamond$ and put them in a set $R_0$. As long as the symbol $\diamond$ is present, the grafting operators $D^I$ and $D^q$ graft the detached nodes on every node except $v_\diamond$. In particular, the following operation is trivial:

$$D^{|\Pi(v)| \to v}(\gamma_{v^\circ})^\diamond = \gamma.$$

**Example.** Consider the forest $\gamma = \bigcirc \bullet^2$ that we label for the sake of clarity. The marked forest $\gamma_{1^\circ}$ has the set of detached nodes $R_0 = \{1\}$ and satisfies

$$\gamma_{1^\circ} = \bullet^2 \bullet^1, \quad (D\gamma_{1^\circ})^\diamond = \bullet^1 \bullet^2, \quad D(\gamma_{1^\circ})^\diamond = \bullet^1 \bullet^2 + \bullet^1 \bullet^2.$$

For simplicity, we write $\gamma_{v \to \hat{k}}^\diamond$ for $(\gamma_{v \to \hat{k}})^\diamond$, 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

$$\binom{n}{p} D^n = D^p D^{n-p}. \qquad (3.1)$$

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)^{|\Pi(v)|}\gamma_{v^\circ}$, $\Pi(v)$ denotes the number of predecessors of $v$ in $\gamma$, not in $\gamma_{v^\circ}$.

In the context of differential geometry, the Euler operator describes the differential forms that are divergences [47]. It has a variety of applications such as the study of conservation laws and Lagrangians for PDEs [48, 54, 26, 27, 50]. 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: \mathcal{F}_{n,p} \to \Omega_{n,p}$ and $\mathcal{E}_v^\circ: \mathcal{F}_0 \to \Omega_{0,1}$ are given by

$$\mathcal{E}_v \gamma = (-1)^{|\Pi(v)|} (D^{\Pi(v)} \gamma_{v^\circ})^\diamond, \quad \mathcal{E}_v^\circ \gamma = (-1)^{|\Pi(v)|} (D^{\Pi(v)} \gamma_{v \to \hat{1}})^\diamond$$

and are extended by linearity on $\Omega_{n,p}$. The Euler operator $\mathcal{E}: \Omega_{n,p} \to \Omega_{n,p}$ and its variant $\mathcal{E}^\circ: \Omega_0 \to \Omega_{0,1}$ are given by

$$\mathcal{E} \gamma = \sum_{v \in V} \mathcal{E}_v \gamma, \quad \mathcal{E}^\circ \gamma = \sum_{v \in V} \mathcal{E}_v^\circ \gamma.$$
The output of $\mathcal{E}^\circ$ on the forms of first orders is

$$ \mathcal{E}^\circ \odot = 0, \quad \mathcal{E}^\circ \bullet = 2 \odot, \quad \mathcal{E}^\circ \bullet = -2 \odot, \quad \mathcal{E}^\circ \odot = -2 \odot. $$

From these examples, we observe that the Euler operator vanishes on aromatic forms that are divergences:

$$ \mathcal{E}^\circ d_H \bullet = 0, \quad \mathcal{E}^\circ d_H \odot = 0, \quad \mathcal{E}^\circ d_H \odot \bullet = 0. $$

**Proposition 3.4.** Let $\gamma \in \mathcal{F}_{n,p}$ and $v \in V$, the Euler operators satisfy

$$ \mathcal{E}_v d_H \gamma = 0, \quad \mathcal{E}_v d_H \gamma = 0 $$

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

$$ \begin{aligned} \mathcal{E}_v d_H \gamma &= \sum_{u \in V} \mathcal{E}_v D^{r_n \to u} \gamma \\ &= (-1)^{|\Pi(v)|+1} (D^{\Pi(v)} D^{r_n} \gamma_{v^\circ})^\diamond + \sum_{u \neq v} (-1)^{|\Pi(v)|} (D^{\Pi(v)} D^{r_n \to u} \gamma_{v^\circ})^\diamond = 0, \end{aligned} $$

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. $\square$

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.

$$ \Omega_1 \xrightarrow{d_H} \Omega_0 \xrightarrow{\mathcal{E}^\circ} \Omega_{0,1} \qquad (3.2) $$

The complex (3.2) is exact if $\mathrm{Im}(d_H) = \mathrm{Ker}(\mathcal{E}^\circ)$; that is, if $\gamma \in \Omega_0$ is a divergence if and only if $\mathcal{E}^\circ \gamma = 0$. A fundamental question in variational calculus [47] 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}_v^q \gamma$ of order $q \ge 0$ on $v \in V$ is

$$ \mathcal{E}_v^q \gamma = (-1)^{|\Pi(v)|-q} (D^{|\Pi(v)|-q} \gamma_{v^\circ})^\diamond, $$

that we extend on $\Omega_{n,p}$ by linearity. The higher Euler operators are

$$ \mathcal{E}^q \gamma = \sum_{v \in V} \mathcal{E}_v^q \gamma. $$

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}_v^q \gamma = 0$ if $q > |\gamma|$. Note also that $\mathcal{E}_v^0 = \mathcal{E}_v$.
**Proposition 3.6.** The Euler operators satisfy for all $\gamma \in \mathcal{F}_{n,p}$ and all vertices $v \in V$,

$$ \gamma = \sum_{q=0}^{\infty} D^q \mathcal{E}_v^q \gamma. \qquad (3.3) $$

In particular, we have

$$ |\gamma| \gamma = \sum_{q=0}^{\infty} D^q \mathcal{E}^q \gamma. \qquad (3.4) $$

In addition, if $\gamma$ satisfies $\mathcal{E}^p \gamma = 0$ for $p = 0, \dots, n-1$, then $\gamma = D^n(h^{(n)}\gamma)$, where

$$ h^{(n)}\gamma = \frac{1}{|\gamma|} \sum_{p=n}^{\infty} \binom{p}{n}^{-1} D^{p-n} \mathcal{E}^p \gamma. $$

The proof shares similarities with the approach of [27] 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 = |\Pi(v)| - 1$. Then, the following holds:

$$ (q+1)D^{q+1}D^{k \to v}(D^p \gamma_{v^\diamond})^\diamond = (k+1)D^q D^{k+1 \to v}(D^p \gamma_{v^\diamond})^\diamond + (p+1)D^q D^{k \to v}(D^{p+1} \gamma_{v^\diamond})^\diamond. $$

More precisely, we have

$$ D^q D^{k \to v} (D^p \gamma_{v^\diamond})^\diamond = \sum_{n=0}^{k} (-1)^{k-n} \binom{q+n}{q} \binom{p+k-n}{p} D^{q+n} (D^{p+k-n} \gamma_{v^\diamond})^\diamond. \qquad (3.5) $$

**Proof of Proposition 3.6.** Let $\gamma \in \mathcal{F}_{n,p}$. With the help of Lemma 3.7, we get

$$ \begin{aligned} \mathcal{E}_v \gamma &= (-1)^{|\Pi(v)|} (D^{|\Pi(v)|} \gamma_{v^\diamond})^\diamond \\ &= (\gamma_{v^\diamond})^\diamond \mathbf{1}_{|\Pi(v)|=0} - (-1)^{|\Pi(v)|-1} (D^{|\Pi(v)|} \gamma_{v^\diamond})^\diamond \mathbf{1}_{|\Pi(v)| \ge 1} \\ &= \gamma \mathbf{1}_{|\Pi(v)|=0} - \frac{(-1)^{|\Pi(v)|-1}}{|\Pi(v)|} D(D^{|\Pi(v)|-1} \gamma_{v^\diamond})^\diamond + \frac{(-1)^{|\Pi(v)|-1}}{|\Pi(v)|} D^{1 \to v}(D^{|\Pi(v)|-1} \gamma_{v^\diamond})^\diamond \end{aligned} $$

We apply the same reasoning to the last term

$$ \frac{(-1)^{|\Pi(v)|-1}}{|\Pi(v)|} D^{1 \to v} (D^{|\Pi(v)|-1} \gamma_{v^\diamond})^\diamond = \gamma \mathbf{1}_{|\Pi(v)|=1} - \frac{(-1)^{|\Pi(v)|-2}}{|\Pi(v)|(|\Pi(v)|-1)} D D^{1 \to v} (D^{|\Pi(v)|-2} \gamma_{v^\diamond})^\diamond \\ + \frac{2(-1)^{|\Pi(v)|-2}}{|\Pi(v)|(|\Pi(v)|-1)} D^{2 \to v} (D^{|\Pi(v)|-2} \gamma_{v^\diamond})^\diamond. $$

Iterating this reasoning, we find

$$ \mathcal{E}_v \gamma = \gamma - \sum_{p=1}^{\infty} \frac{(-1)^{|\Pi(v)|-p}}{p^{|\Pi(v)|}} D D^{p-1 \to v} (D^{|\Pi(v)|-p} \gamma_{v^\diamond})^\diamond. \qquad (3.6) $$

Using the formula (3.5) in (3.6), standard combinatorics [44] yield

$$ \gamma = \mathcal{E}_v \gamma + \sum_{p=1}^{\infty} \frac{(-1)^{|\Pi(v)|-p}}{p^{|\Pi(v)|}} \sum_{n=0}^{p-1} (-1)^{p-n-1} (n+1) \binom{|\Pi(v)|-n-1}{|\Pi(v)|-p} D^{n+1} (D^{|\Pi(v)|-n-1} \gamma_{v^\diamond})^\diamond $$
$$ = \mathcal{E}_v \gamma + \sum_{n=0}^{\infty} \sum_{p=n+1}^{\infty} \frac{n+1}{p^{|\Pi(v)|}} \binom{|\Pi(v)|-n-1}{|\Pi(v)|-p} D^{n+1} \mathcal{E}_v^{n+1} \gamma $$

$$ = \mathcal{E}_v \gamma + \sum_{n=1}^{|\Pi(v)|} \sum_{p=n}^{|\Pi(v)|} \frac{n(|\Pi(v)|-n)!(p-1)!}{|\Pi(v)|!(p-n)!} D^n \mathcal{E}_v^n \gamma $$

$$ = \mathcal{E}_v \gamma + \sum_{n=1}^{|\Pi(v)|} \binom{|\Pi(v)|}{n}^{-1} \sum_{p=0}^{|\Pi(v)|-n} \binom{p+n-1}{p} D^n \mathcal{E}_v^n \gamma = \sum_{n=0}^{\infty} D^n \mathcal{E}_v^n \gamma. $$

We sum (3.3) on all the nodes $v \in V$ to obtain (3.4). The last claim of Proposition 3.6 is obtained by using the formula (3.1) in (3.4). □

**Proof of Lemma 3.7.** Using the identity (3.1), we distribute the derivatives on $v$ and the other nodes.

$$ \begin{aligned} D^{q+1} D^{k \to v} (D^p \gamma_{v^\diamond})^\diamond &= \frac{1}{q+1} D D^q D^{k \to v} (D^p \gamma_{v^\diamond})^\diamond \\ &= \frac{1}{q+1} \sum_{\substack{S_1 \sqcup S_2 \sqcup S_3 \cup \{u\} = \Pi(v) \\ |S_1| = q, |S_2| = k, |S_3| = p}} D^u D^{S_1} D^{S_2 \to v} (D^{S_3} \gamma_{v^\diamond})^\diamond \\ &= \frac{1}{q+1} \sum_{\substack{S_1 \sqcup S_2 \sqcup S_3 \cup \{u\} = \Pi(v) \\ |S_1| = q, |S_2| = k, |S_3| = p}} D^{S_1} \left[ D^{S_2 \cup u \to v} (D^{S_3} \gamma_{v^\diamond})^\diamond + D^{S_2 \to v} (D^{S_3 \cup u} \gamma_{v^\diamond})^\diamond \right] \\ &= \frac{k+1}{q+1} \sum_{\substack{S_1 \sqcup S_2 \sqcup S_3 = \Pi(v) \\ |S_1| = q, |S_2| = k+1, |S_3| = p}} D^{S_1} D^{S_2 \to v} (D^{S_3} \gamma_{v^\diamond})^\diamond \\ &\quad + \frac{p+1}{q+1} \sum_{\substack{S_1 \sqcup S_2 \sqcup S_3 = \Pi(v) \\ |S_1| = q, |S_2| = k, |S_3| = p+1}} D^{S_1} D^{S_2 \to v} (D^{S_3} \gamma_{v^\diamond})^\diamond \\ &= \frac{k+1}{q+1} D^q D^{k+1 \to v} (D^p \gamma_{v^\diamond})^\diamond + \frac{p+1}{q+1} D^q D^{k \to v} (D^{p+1} \gamma_{v^\diamond})^\diamond. \end{aligned} $$

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

$$ h_V \gamma = \frac{1}{|\gamma|} \gamma^{(1) \to \bullet}, $$

and for $\gamma \in \mathcal{F}_0$, let the horizontal homotopy operator $h_H$ be

$$ h_H \gamma = \frac{1}{|\gamma|} \sum_{q=1}^{\infty} \frac{1}{q} D^{q-1} \mathcal{E}^q \gamma. $$

We extend the definition on $\Omega_{0,1}$ and $\Omega_0$ by linearity. These operators satisfy for all $\gamma \in \Omega_0$,

$$ (d_H h_H + h_V \mathcal{E}^\circ) \gamma = \gamma. \qquad (3.7) $$

In particular, the variational complex (3.2) is exact.
Proof. Let $\gamma \in \Omega_0$. From Proposition 3.6, we obtain

$$ \gamma = \frac{1}{|\gamma|} \sum_{q=0}^{\infty} D^q \mathcal{E}^q \gamma = \frac{1}{|\gamma|} \mathcal{E} \gamma + \frac{1}{|\gamma|} \sum_{q=1}^{\infty} \frac{1}{q} D D^{q-1} \mathcal{E}^q \gamma = h_V \mathcal{E}^{\circ} \gamma + d_H h_H \gamma. $$

Proposition 3.4 yields that $\mathrm{Im}(d_H) \subset \mathrm{Ker}(\mathcal{E})$. If $\gamma \in \mathrm{Ker}(\mathcal{E})$, then the identity (3.7) becomes $\gamma = d_H(h_H \gamma) \in \mathrm{Im}(d_H)$. The exactness of the chain follows straightforwardly. $\square$

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

$$ \begin{aligned} d_H^2 \wedge \gamma &= \frac{1}{n!} \wedge \sum_{\sigma \in S_n^\bullet} \varepsilon(\sigma) D^{r_{n-1}} D^{r_n} \sigma \gamma \\ &= \frac{1}{n!} \wedge \sum_{\sigma \in S_n^\bullet} \varepsilon((n(n-1))\sigma) D^{r_{n-1}} D^{r_n} (n(n-1))\sigma \gamma \\ &= -\frac{1}{n!} \wedge \sum_{\sigma \in S_n^\bullet} \varepsilon(\sigma) D^{r_n} D^{r_{n-1}} \sigma \gamma \\ &= -d_H^2 \wedge \gamma, \end{aligned} $$

where $(n(n-1)) \in S_n^\bullet$ is a transposition and where we used that $D^{r_{n-1}}$ and $D^{r_n}$ commute. For the vertical derivative, a similar approach on $\gamma \in \mathcal{F}_{n,p-1}$ gives

$$ \begin{aligned} d_V^2 \wedge \gamma &= \sum_{\substack{v,w \in V \\ v \neq w}} \wedge (\wedge_{\gamma v \to \mathring{p}})_{w \to \mathring{p+1}} \\ &= \frac{1}{p!(p+1)!} \wedge \cdot \sum_{v,w \in V} \sum_{\sigma \in S_p^{\circ}} \sum_{\tilde{\sigma} \in S_{p+1}^{\circ}} \varepsilon(\tilde{\sigma}\sigma)\tilde{\sigma}(\sigma(\gamma_{v \to \mathring{p}})_{w \to \mathring{p+1}}) \\ &= \frac{1}{p!(p+1)!} \wedge \cdot \sum_{v,w \in V} \sum_{\sigma \in S_p^{\circ}} \sum_{\tilde{\sigma} \in S_{p+1}^{\circ}} \varepsilon(\tilde{\sigma}\sigma)\tilde{\sigma}\sigma(\gamma_{v \to \mathring{p}}, w \to \mathring{p+1}) \\ &= \frac{1}{p!(p+1)!} \wedge \cdot \sum_{v,w \in V} \sum_{\sigma \in S_p^{\circ}} \sum_{\tilde{\sigma} \in S_{p+1}^{\circ}} \varepsilon(\tilde{\sigma}\sigma)\tilde{\sigma}(\sigma(\gamma_{v \to \mathring{p}})_{w \to \mathring{p+1}}) \\ &= -\frac{1}{p!(p+1)!} \wedge \cdot \sum_{v,w \in V} \sum_{\sigma \in S_p^{\circ}} \sum_{\tilde{\sigma} \in S_{p+1}^{\circ}} \varepsilon(\tilde{\sigma}\sigma)\tilde{\sigma}(\sigma(\gamma_{v \to \mathring{p}})_{w \to \mathring{p+1}}) \\ &= -d_V^2 \wedge \gamma. \end{aligned} $$

As $\Omega_{n,p} = \wedge(\mathcal{F}_{n,p})$, we get the desired identities. The commutativity of $d_H$ and $d_V$ is straightforward. $\square$
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: \Omega_{n,p} \to \Omega_{n,p-1}$ given by

$$h_V \gamma = \frac{p}{|\gamma|} \gamma(\mathring{\gamma}) \to \bullet,$$

satisfies for $p \ge 1$ and $\gamma \in \Omega_{n,p}$,

$$(d_V h_V + h_V d_V) \gamma = \gamma.$$

*Proof.* Let $\gamma \in \mathcal{F}_{n,p}$. On the first hand, we have

$$d_V h_V \gamma = \frac{p}{|\gamma|} \wedge \gamma + \frac{p}{|\gamma|} \sum_{v \in V \cdot} \wedge \gamma(\mathring{\gamma}) \to \bullet, v \to \mathring{\gamma} \\ = \frac{p}{|\gamma|} \wedge \gamma + \frac{1}{|\gamma| (p-1)!} \wedge \bullet \sum_{v \in V \cdot} \sum_{\sigma \in S_p^{\circ}} \varepsilon(\sigma) \sigma(\gamma(\mathring{\gamma}) \to \bullet, v \to \mathring{\gamma}).$$

On the other hand, we have

$$h_V d_V \gamma = \frac{1}{|\gamma|} \sum_{v \in V \cdot} h_V(\wedge \gamma_{v \to (p+1)}) \\ = \frac{p+1}{|\gamma| (p+1)!} \wedge \bullet \sum_{v \in V \cdot} \sum_{\sigma \in S_{p+1}^{\circ}} \varepsilon(\sigma) \sigma(\gamma_{v \to (p+1)})(p+1) \to \bullet \\ = \frac{1}{|\gamma| p!} \wedge \bullet \sum_{v \in V \cdot} \sum_{\substack{\sigma \in S_{p+1}^{\circ} \\ \sigma((p+1)) = (p+1)}} \varepsilon(\sigma) \sigma(\gamma_{v \to (p+1)})(p+1) \to \bullet \\ + \frac{1}{|\gamma| p!} \wedge \bullet \sum_{v \in V \cdot} \sum_{\substack{\sigma \in S_{p+1}^{\circ} \\ \sigma((p+1)) \neq (p+1)}} \varepsilon(\sigma) \sigma(\gamma_{v \to (p+1)})(p+1) \to \bullet \\ = \frac{1}{|\gamma| p!} \wedge \bullet \sum_{v \in V \cdot} \sum_{\sigma \in S_p^{\circ}} \varepsilon(\sigma) \sigma \gamma \\ + \frac{1}{|\gamma| p!} \wedge \bullet \sum_{i=1}^{p} \sum_{v \in V \cdot} \sum_{\substack{\sigma \in S_p^{\circ} \\ \sigma(\mathring{i}) = \mathring{p}}} \varepsilon(\sigma) \sigma(\gamma(\mathring{i}) \to \bullet, v \to \mathring{p}) \\ = \frac{|\gamma| - p}{|\gamma|} \wedge \gamma + \frac{1}{|\gamma| p!} \wedge \bullet \sum_{i=1}^{p} \sum_{v \in V \cdot} \sum_{\sigma \in S_p^{\circ}} \varepsilon(\sigma(\mathring{i}) \mathring{p}) \sigma(\gamma(\mathring{p}) \to \bullet, v \to \mathring{p}) \\ = \frac{|\gamma| - p}{|\gamma|} \wedge \gamma - \frac{1}{|\gamma| (p-1)!} \wedge \bullet \sum_{v \in V \cdot} \sum_{\sigma \in S_p^{\circ}} \varepsilon(\sigma) \sigma(\gamma(\mathring{p}) \to \bullet, v \to \mathring{p}).$$

We deduce that

$$(d_V h_V + h_V d_V) \gamma = \wedge \gamma,$$

and we obtain the desired equality by linearity as $\Omega_{n,p} = \wedge \operatorname{Span}(\mathcal{F}_{n,p})$ and $\wedge^2 = \wedge$. $\square$
The expression of the horizontal homotopy operator is much more technical. We refer the reader to [3, 57, 47] 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: \Omega_{n,p} \to \Omega_{n+1,p}$ by

$$h_H \gamma = \frac{1}{|\gamma|} \sum_{q=0}^{\infty} \frac{n+1}{q+n+1} \wedge D^q \mathcal{E}^{q+1} \gamma.$$

It satisfies for $n \ge 1$ and $\gamma \in \Omega_{n,p}$,

$$(d_H h_H + h_H d_H) \gamma = \gamma.$$

*Proof.* For the sake of simplicity, we consider $\gamma \in \mathcal{F}_n$. We define for $v \in V$

$$J_v \gamma = \sum_{q=0}^{\infty} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{n+1}{q+n+1} D^I \mathcal{E}_v^{I, \hat{r}} \gamma,$$

where we denote $\mathcal{E}_v^{I, \hat{r}} \gamma = (-1)^{|\Pi(v)|-(q+1)} (D^J \gamma_{v^\circ})^\diamond$ and where $\hat{r}$ becomes the new root of $J_v \gamma$. On the first hand, we have

$$d_H \wedge \gamma = d_H \sum_{\sigma \in S_n} \frac{\varepsilon(\sigma)}{n!} \sigma \gamma = \sum_{\sigma \in S_n} \frac{\varepsilon(\sigma)}{n!} D^{r_{\sigma^{-1}(n)}} \sigma \gamma,$$

and $J_v d_H \wedge \gamma$ is given by

$$J_v d_H \wedge \gamma = \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{u \neq v} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q+n)(n-1)!} D^I \mathcal{E}_v^{I, \hat{r}} D^{r_{\sigma^{-1}(n)}} \to u \sigma \gamma \\ + \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\Pi(v) \cup \{r_k\} = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q+n)(n-1)!} D^I \mathcal{E}_v^{I, \hat{r}} D^{r_{\sigma^{-1}(n)}} \to v \sigma \gamma.$$

The expression of $D^I \mathcal{E}_v^{I, \hat{r}} D^{r \to u} \gamma$ satisfies for $u \in V$:

$$D^I \mathcal{E}_v^{I, \hat{r}} D^{r \to u} \gamma = \begin{cases} D^I \mathcal{E}_v^I \gamma & \text{if } \hat{r} = r, \ u = v, \\ \sum_{w \in V} D^{r \to w} D^{I \setminus \{r\}} \mathcal{E}_v^{I \setminus \{r\}, \hat{r}} \gamma & \text{if } r \in I, \ u = v, \\ - \sum_{w \neq v} D^{r \to w} D^I \mathcal{E}_v^{I, \hat{r}} \gamma & \text{if } r \in J, \ u = v, \\ D^{r \to u} D^I \mathcal{E}_v^{I, \hat{r}} \gamma & \text{if } u \neq v. \end{cases}$$

Thus, for $\gamma \in \mathcal{F}_{n,p}$, $J_v d_H \wedge \gamma$ is given by

$$J_v d_H \wedge \gamma = \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\Pi(v) = I \cup J} \frac{\varepsilon(\sigma)}{(q+n)(n-1)!} D^I \mathcal{E}_v^I \sigma \gamma \\ + \sum_{q=1}^{\infty} \sum_{\sigma \in S_n} \sum_{w \in V} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q+n)(n-1)!} D^{r_{\sigma^{-1}(n)}} \to w D^I \mathcal{E}_v^{I, \hat{r}} \sigma \gamma$$
$$
\begin{align*}
& - \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{w \neq v} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q+n)(n-1)!} D^{r_{\sigma^{-1}(n)}} \mathcal{E}_v^{I, \hat{r}} \sigma \gamma \\
& + \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{u \neq v} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q+n)(n-1)!} D^{r_{\sigma^{-1}(n)}} \mathcal{E}_v^{I, \hat{r}} \sigma \gamma \\
&= \sum_{q=0}^{\infty} \frac{n}{q+n} D^q \mathcal{E}_v^q \wedge \gamma \\
&+ \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q+n+1)(n-1)!} D^{r_{\sigma^{-1}(n)}} \mathcal{E}_v^{I, \hat{r}} \sigma \gamma.
\end{align*}
$$

As $D^q \mathcal{E}_v^q \wedge \gamma$ is unchanged by the application of the wedge operator, we find

$$
\begin{align*}
\wedge J_v d_H \wedge \gamma &= \sum_{q=0}^{\infty} \frac{n}{q+n} D^q \mathcal{E}_v^q \wedge \gamma \\
&\quad + \sum_{q=0}^{\infty} \sum_{\nu \in S_n} \sum_{\sigma \in S_n} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma\nu)}{(q+n+1)(n-1)!n!} \nu D^{r_{\sigma^{-1}(n)}} D^I \mathcal{E}_v^{I,\hat{r}} \sigma\gamma.
\end{align*}
$$

On the other hand, we have

$$
\begin{align*}
d_H \wedge J_v \wedge \gamma &= \sum_{\sigma \in S_n} \frac{\varepsilon(\sigma)}{n!} d_H \wedge J_v \sigma \gamma \\
&= d_H \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\nu \in S_{n+1}} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma\nu)}{(q+n+1)(n!)^2} \nu D^I \mathcal{E}_v^{I,\hat{r}} \sigma \gamma \\
&= \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\nu \in S_{n+1}} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma\nu)}{(q+n+1)(n!)^2} D^{r_{(\sigma\nu)-1(n+1)}} \nu D^I \mathcal{E}_v^{I,\hat{r}} \sigma \gamma \\
&= \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q+n+1)n!} D^{I,\hat{r}} \mathcal{E}_v^{I,\hat{r}} \sigma \gamma \\
&\quad - \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\nu \in S_n} \sum_{\eta \in S_n} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma\nu\eta)}{(q+n+1)(n-1)!(n!)^2} D^{r_{(\sigma\nu)-1(n)}} \eta(n(n+1)) \nu D^I \mathcal{E}_v^{I,\hat{r}} \sigma \gamma \\
&= \sum_{q=0}^{\infty} \frac{q+1}{q+n+1} D^{q+1} \mathcal{E}_v^{q+1} \wedge \gamma \\
&\quad - \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\eta \in S_n} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma\eta)}{(q+n+1)(n-1)!n!} \eta D^{r_{\sigma^{-1}(n)}} D^I \mathcal{E}_v^{I,\hat{r}} \sigma \gamma \\
&= \sum_{q=1}^{\infty} \frac{q}{q+n} D^q \mathcal{E}_v^q \wedge \gamma
\end{align*}
$$
$$ -\sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\eta \in S_n} \sum_{\Pi(v)=\{\hat{r}\} \cup I \cup J} \sum_{|I|=q} \frac{\varepsilon(\sigma\eta)}{(q+n+1)(n-1)!n!} \eta D^{r_{\sigma^{-1}(n)}} D^I \mathcal{E}_v^{I, \hat{r}} \sigma\gamma, $$

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

$$ (\wedge J_v d_H + d_H \wedge J_v) \wedge \gamma = \sum_{q=0}^{\infty} D^q \mathcal{E}_v^q \wedge \gamma = \wedge \gamma. \quad (3.8) $$

Summing on $v \in V$ gives the desired homotopy identity. $\square$

**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 ill-defined, 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 divergence-free vector field

This section is devoted to the proof of Theorem 2.10. Proposition 2.6 and Proposition 3.9 extend naturally to the divergence-free context. Proposition 3.10 is not valid anymore and we replace it by the following result.

**Proposition 3.12.** We define

$$ \tilde{h}_H \gamma = \frac{1}{|\gamma|} \sum_{v \in V} \sum_{q=0}^{\infty} \frac{n+1}{q+n+1} \wedge D^q \mathcal{E}_v^{q+1} \gamma. $$

For $n > 1$ and $\gamma \in \tilde{\Omega}_{n,p}$, the horizontal homotopy identity is

$$ (d_H \tilde{h}_H + \tilde{h}_H d_H)\gamma = \gamma, $$

while for $n = 1$, the identity is

$$ (d_H \tilde{h}_H + \tilde{h}_H d_H)\gamma = \gamma - R\gamma. \quad (3.9) $$

The remainder in (3.9) is the linear map $R\gamma = \frac{1}{|\gamma|} \mathcal{E}_r\gamma$ for $\gamma \in \tilde{\mathcal{F}}_{1,p}$, with $r$ the root of $\gamma$. Moreover, if $\gamma \in \tilde{\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 divergence-free 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 \tilde{\mathcal{F}}_n$. For $v \in V$, we define

$$ \tilde{J}_v(\gamma) = \sum_{q=0}^{\infty} \sum_{\Pi(v)=\{\hat{r}\} \cup I \cup J} \frac{n+1}{q+n+1} \wedge D^I \mathcal{E}_v^{I, \hat{r}} \gamma. $$
We have if $r = v$,

$$D^I \mathcal{E}_v^{I, \hat{r}} D^{r \to u} \gamma = \begin{cases} 0 & \text{if } u = v, \\ D^{r \to u} D^I \mathcal{E}_v^{I, \hat{r}} \gamma & \text{if } u \neq v, \end{cases}$$

and if $r \neq v$,

$$D^I \mathcal{E}_v^{I, \hat{r}} D^{r \to u} \gamma = \begin{cases} D^I \mathcal{E}_v^{I, \hat{r}} \gamma & \text{if } \hat{r} = r, \quad u = v, \\ \sum_{w \in V} D^{r \to w} D^{I \setminus \{r\}} \mathcal{E}_v^{I \setminus \{r\}, \hat{r}} \gamma & \text{if } r \in I, \quad u = v, \\ - \sum_{w \neq v} D^{r \to w} D^I \mathcal{E}_v^{I, \hat{r}} \gamma & \text{if } r \in J, \quad u = v, \\ D^{r \to u} D^I \mathcal{E}_v^{I, \hat{r}} \gamma & \text{if } u \neq v. \end{cases}$$

Thus, $\tilde{J}_v d_H \wedge \gamma$ is given by

$$\begin{aligned} \tilde{J}_v d_H \wedge \gamma = {}& 1_{v \neq r_{\sigma^{-1}(n)}} \left[ \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\Pi(v)=I \cup J} \frac{\varepsilon(\sigma)}{(q + n - 1_{v \in R})(n-1)!} D^I \mathcal{E}_v^I \sigma \gamma \right. \\ & \left. + \sum_{q=1}^{\infty} \sum_{\sigma \in S_n} \sum_{w \in V} \sum_{\Pi(v)=\{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q + n - 1_{v \in R})(n-1)!} D^{r_{\sigma^{-1}(n)} \to w} D^I \mathcal{E}_v^{I, \hat{r}} \sigma \gamma \right. \\ & \left. - \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{w \neq v} \sum_{\Pi(v)=\{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q + n - 1_{v \notin R})(n-1)!} D^{r_{\sigma^{-1}(n)} \to w} D^I \mathcal{E}_v^{I, \hat{r}} \sigma \gamma \right] \\ & + \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{u \neq v} \sum_{\Pi(v)=\{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q + n - 1_{v \in R})(n-1)!} D^{r_{\sigma^{-1}(n)} \to u} D^I \mathcal{E}_v^{I, \hat{r}} \sigma \gamma \\ = {}& 1_{v \neq r_{\sigma^{-1}(n)}} \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \frac{\varepsilon(\sigma)}{(q + n - 1_{v \in R})(n-1)!} D^q \mathcal{E}_v^q \sigma \gamma \\ & + \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\Pi(v)=\{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q + n + 1_{v \notin R})(n-1)!} D^{r_{\sigma^{-1}(n)}} D^I \mathcal{E}_v^{I, \hat{r}} \sigma \gamma \end{aligned}$$

Then, we find for $n = 1$ and $v \in R$,

$$\wedge \tilde{J}_v d_H \wedge \gamma = \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\nu \in S_n} \sum_{\Pi(v)=\{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma \nu)}{(q + n + 1_{v \notin R})(n-1)!n!} \nu D^{r_{\sigma^{-1}(n)}} D^I \mathcal{E}_v^{I, \hat{r}} \sigma \gamma,$$

and otherwise

$$\begin{aligned} \wedge \tilde{J}_v d_H \wedge \gamma = {}& \sum_{q=0}^{\infty} \frac{n - 1_{v \in R}}{q + n - 1_{v \in R}} D^q \mathcal{E}_v^q \wedge \gamma \\ & + \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\nu \in S_n} \sum_{\Pi(v)=\{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma \nu)}{(q + n + 1_{v \notin R})(n-1)!n!} \nu D^{r_{\sigma^{-1}(n)}} D^I \mathcal{E}_v^{I, \hat{r}} \sigma \gamma, \end{aligned}$$
where for $v \in R$, the wedge adds the missing term $v = r_{\sigma^{-1}(n)}$, and rescales the expression with a coefficient $\frac{n-1}{n}$. On the other hand, we have

$$
\begin{align*}
d_H \wedge \tilde{J}_v \wedge \gamma &= d_H \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\nu \in S_{n+1}} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma\nu)}{(q + n + 1_{v \notin R})(n!)^2} \nu D^I \mathcal{E}_v^{I, \hat{r}} \sigma\gamma \\
&= \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\nu \in S_{n+1}} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma\nu)}{(q + n + 1_{v \notin R})(n!)^2} D^{r(\sigma\nu)-1(n+1)} \nu D^I \mathcal{E}_v^{I, \hat{r}} \sigma\gamma \\
&= \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma)}{(q + n + 1_{v \notin R})n!} D^{I, \hat{r}} \mathcal{E}_v^{I, \hat{r}} \sigma\gamma \\
&\quad - \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\nu \in S_n} \sum_{\eta(n)=n} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma\nu\eta)}{(q + n + 1_{v \notin R})(n-1)!(n!)^2} \\
&\quad \cdot D^{r(\sigma\nu)-1(n)} \eta(n(n+1)) \nu D^I \mathcal{E}_v^{I, \hat{r}} \sigma\gamma \\
&= \sum_{q=0}^{\infty} \frac{q+1}{q+n+1_{v \notin R}} D^{q+1} \mathcal{E}_v^{q+1} \wedge \gamma \\
&\quad - \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\eta \in S_n} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma\eta)}{(q+n+1_{v \notin R})(n-1)!n!} \eta D^{r_{\sigma^{-1}(n)}} D^I \mathcal{E}_v^{I, \hat{r}} \sigma\gamma \\
&= \sum_{q=1}^{\infty} \frac{q}{q+n-1_{v \in R}} D^q \mathcal{E}_v^q \wedge \gamma \\
&\quad - \sum_{q=0}^{\infty} \sum_{\sigma \in S_n} \sum_{\eta \in S_n} \sum_{\Pi(v) = \{\hat{r}\} \cup I \cup J} \frac{\varepsilon(\sigma\eta)}{(q+n+1_{v \notin R})(n-1)!n!} \eta D^{r_{\sigma^{-1}(n)}} D^I \mathcal{E}_v^{I, \hat{r}} \sigma\gamma.
\end{align*}
$$

Using Proposition 3.6, we get for $n > 1$,

$$ (\wedge \tilde{J}_v d_H + d_H \wedge \tilde{J}_v) \wedge \gamma = \sum_{q=0}^{\infty} D^q \mathcal{E}_v^q \wedge \gamma = \wedge \gamma, $$

and for $n = 1$,

$$ (\wedge \tilde{J}_v d_H + d_H \wedge \tilde{J}_v) \wedge \gamma = \sum_{q=0}^{\infty} D^q \mathcal{E}_v^q \wedge \gamma - \mathcal{E}_r \wedge \gamma 1_{v \in R} = \wedge \gamma - \mathcal{E}_r \wedge \gamma 1_{v \in R}. $$

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 \text{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

$$ \mathcal{E}_r \gamma = \frac{1}{N} \mathcal{E}_r \gamma. $$

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$. $\square$
**Remark 3.13.** Let $\gamma \in \tilde{\Omega}_1^N$ with $N > 1$, define the modified homotopy operators by

$$ \tilde{h}_H^1 = \tilde{h}_H \left( 1 + \frac{1}{N-1} \mathcal{E} \right), \quad \tilde{h}_H^2 = \tilde{h}_H \left( 1 + \frac{1}{N-1} \mathcal{E}_r \right). $$

Then, the homotopy identity (3.9) on $\tilde{\Omega}_1^N$ with $N > 1$ is replaced by the simpler identity

$$ (d_H \tilde{h}_H^2 + \tilde{h}_H^1 d_H)\gamma = \gamma. \qquad (3.10) $$

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 $\tilde{\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 B-series of an aromatic volume-preserving integrator and we prove that an aromatic Runge-Kutta method cannot be volume-preserving.

## 4.1 The augmented aromatic bicomplex

Following [1, 2] 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 \ge 1$, define the interior Euler operator $I: \Omega_{0,p} \to \Omega_{0,p}$ by

$$ I\gamma = \wedge \mathcal{E}_{\mathbb{P}}\gamma = (-1)^{|\Pi(\mathbb{P})|} \wedge (D^{|\Pi(\mathbb{P})|} \gamma_{(\mathbb{P})^\circ})^\diamond, $$

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$.

$$ \dots \xrightarrow{d_H} \Omega_2 \xrightarrow{d_H} \Omega_1 \xrightarrow{d_H} \Omega_0 \xrightarrow{\delta_V} \mathcal{I}_1 \xrightarrow{\delta_V} \mathcal{I}_2 \xrightarrow{\delta_V} \dots \qquad (4.1) $$

**Theorem 4.1.** Define the augmented homotopy operators as

$$ \mathfrak{h}_H\gamma = \wedge \sum_{q=1}^{\infty} \frac{1}{q} D^{q-1} \mathcal{E}_{\mathbb{P}}^q \gamma, \quad \mathfrak{h}_V\gamma = I \circ h_V. $$

Then, the maps $I$ and $\delta_V$ satisfy

$$ I^2 = I, \quad \mathrm{Id}_H = 0, \quad \delta_V^2 = 0, $$

and the following identities hold

$$ \gamma = (I + d_H \mathfrak{h}_H)\gamma, \quad \gamma \in \Omega_{0,p}, \quad p \ge 1, $$
<figure><img src="image_5.png" /><figcaption>Figure 3: The augmented aromatic bicomplex.</figcaption></figure>

$$ \begin{align*} \gamma &= (\delta_V h_V + h_V \delta_V)\gamma, \quad \gamma \in \mathcal{I}_1, \\ \gamma &= (\delta_V h_V + h_V \delta_V)\gamma, \quad \gamma \in \mathcal{I}_p, \quad p > 1. \end{align*} $$

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 [1, Chap. 4] for the proof of Theorem 4.1.

*Proof.* We first observe that for all $\gamma \in \mathcal{F}_{1,p}$,

$$ \mathrm{Id}_H \gamma = \wedge \sum_{v \neq \mathring{p}} (-1)^{|\Pi(\mathring{p})|} (D^{|\Pi(\mathring{p})|} D^{r_1 \to v} \gamma_{\mathring{p}}^{\circ})^{\diamond} + \wedge (-1)^{|\Pi(\mathring{p})|+1} (D^{\Pi(\mathring{p})} D^{r_1} \gamma_{\mathring{p}}^{\circ})^{\diamond} = 0. $$

Using Equation (3.3) with $v = \mathring{p}$ gives the augmented horizontal homotopy identity

$$ (I + d_H h_H) \wedge \gamma = \wedge \gamma. $$

Applying $I$ to this last equality yields $I \circ I = I$. Let us now look at the Euler-Lagrange 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$,

$$ (Id_V + d_H h_H d_V)\gamma = d_V \gamma. $$

We apply $d_V$ and use that $d_V$ and $d_H$ commute:

$$ (d_V Id_V + d_H d_V h_H d_V)\gamma = 0. $$

Since $Id_H = 0$, applying $I$ yields

$$ \delta_V^2 \gamma = 0. $$

Let $\gamma \in \mathcal{I}_p$, the augmented horizontal homotopy applied to $d_V \gamma$ gives

$$ d_V \gamma = \delta_V \gamma + d_H h_H d_V \gamma. \qquad (4.2) $$

If $p = 1$, then we use the identity (4.2) in the vertical homotopy identity to get

$$ \gamma = d_V h_V \gamma + h_V \delta_V \gamma + d_H \tilde{\gamma}, $$
where we used that $h_V$ and $d_H$ commute and where $\tilde{\gamma} \in \Omega_{1,1}$. Applying $I$ yields

$$ \gamma = (\delta_V h_V + h_V \delta_V)\gamma. $$

If $p > 1$, we apply the augmented horizontal homotopy identity to $h_V\gamma$,

$$ h_V\gamma = h_V\gamma + d_H h_H h_V\gamma. $$

Applying the vertical derivative $d_V$ gives

$$ d_V h_V \gamma = d_V h_V \gamma + d_H d_V h_H h_V \gamma, \qquad (4.3) $$

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

$$ \gamma = d_V h_V \gamma + h_V \delta_V \gamma + d_H \tilde{\gamma}, $$

where $\tilde{\gamma} \in \Omega_{1,p}$. As $I \circ d_H = 0$ and $I\gamma = \gamma$, we find

$$ \gamma = (\delta_V h_V + h_V \delta_V)\gamma. $$

The exactness of the augmented aromatic bicomplex is a straightforward consequence of Theorem 2.9 and of the augmented homotopy identities. $\square$

## 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 divergence-free case. The primary motivation was to compute the dimension of the space of solenoidal (i.e., divergence-free) aromatic trees of each order. However, the result revealed a surprisingly simple connection with another combinatorial object, the self-looped 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 [46]) has generating function

$$ t(z) = \sum_{N=1}^{\infty} \mathcal{T}^N z^n = z + z^2 + 2z^3 + 4z^4 + 9z^5 + 20z^6 + \dots $$

and satisfies the functional equation

$$ t(z) = z \exp \left( \sum_{k=1}^{\infty} \frac{1}{k} t(z^k) \right). \qquad (4.4) $$

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, N-1\} \to \{0, \dots, N-1\}$ modulo the symmetric group $S_{N-1}$ (node 0 is the root, and we “forget the labels”).

The number $|\Omega_0^N|$ of scalar aromas in addition to the empty aroma (sequence A001372) has generating function

$$ a(z) := 1 + \sum_{N=1}^{\infty} |\Omega_0^N| z^n = 1 + z + 3z^2 + 7z^3 + 19z^4 + 47z^5 + 130z^6 + \dots $$
and is related to $t(z)$ by the equation

$$a(z) = \prod_{k=1}^{\infty} (1 - t(z^k))^{-1}.$$

The scalar aromas are mapping patterns of functions $\{1, \dots, N\} \to \{1, \dots, N\} \bmod S_N$.

We introduce two new spaces, the aromatic forms with a 1-loop, called the self-looped aromatic forms $\tilde{\Omega}_{n,k}$, and their complements, the non-self-looped aromatic forms $\bar{\Omega}_{n,k}$, and generating functions $\dot{a}(z)$ of $\tilde{\Omega}_0$ and $\bar{a}(z)$ of $\bar{\Omega}_0$. As for mapping patterns on $\{1, \dots, N\}$, the mappings in $\tilde{\Omega}_0$ (enumerated by sequence A217896) have at least one fixed point and those in $\bar{\Omega}_0$ (also known as the “functional digraphs”, enumerated by sequence A001373) have no fixed points.

Each pair consisting of one non-self-looped 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 self-loop. The process is invertible (redirect all self-loops in an arbitrary scalar aroma to a new root). Expressed in terms of generating functions, it writes

$$za(z) = t(z)\bar{a}(z).$$

Therefore

$$\dot{a}(z) = a(z) - 1 - \bar{a}(z) = \frac{a(z)(t(z) - z)}{t(z)} - 1 = z + 2z^2 + 5z^3 + 13z^4 + 34z^5 + 90z^6 + \dots \quad (4.5)$$

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

$$b_p(u, z) = \sum_{k=0}^{\infty} \sum_{N=1}^{\infty} |\Omega_{k,p}^N| u^k z^N$$

be the bivariate generating function for row $p$ of the aromatic bicomplex, let

$$c_p(z) = \sum_{N=1}^{\infty} |\mathcal{I}_p^N| z^N$$

be the generating function of the type-$p$ functional forms, and let

$$s(z) = \sum_{N=1}^{\infty} |\Psi^N| z^N$$

be the generating function of the solenoidal aromatic trees. Then

$$b_0(u, z) = a(z) \exp \left( \sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{k} u^k t(z^k) \right),$$

$$b_1(u, z) = b_0(u, z) \frac{t(z)(1 + u - u t(z))}{(1 - t(z))^2},$$

$$c_1(z) = z b_1(0, z) = \frac{z a(z) t(z)}{(1 - t(z))^2},$$

$$s(z) = a(z) t(z) - \dot{a}(z).$$
*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 enumerated<fnref n="1" /> by $a(z)t(z)$.

An element of $\Omega_k^N$ 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 [20, 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

$$
\begin{align*}
s(z) &= \sum_{N=1}^{\infty} (|\Omega_2^N| - |\Omega_3^N| + |\Omega_4^N| - |\Omega_5^N| + \dots) z^N \\
&= \sum_{N=1}^{\infty} \left( |\Omega_1^N| - |\Omega_0^N| + \sum_{k=0}^{\infty} (-1)^k |\Omega_k^N| \right) z^N \\
&= a(z)t(z) - a(z) + 1 + b_0(-1, z) \\
&= a(z)t(z) - a(z) + 1 + a(z) \exp \left( - \sum_{k=1}^{\infty} \frac{1}{k} t(z^k) \right) \\
&= a(z)t(z) - a(z) + 1 + \frac{a(z)z}{t(z)} \tag*{using Eq. (4.4)} \\
&= a(z)t(z) - \hat{a}(z) \tag*{using Eq. (4.5)}
\end{align*}
$$

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; or

(ii) An unmarked scalar aroma times the wedge product of $n-1$ 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.

<footnote n="1">¹This is a right-shift of sequence A126285, the partial mapping patterns. From a partial mapping pattern on $N-1$ nodes, add a new node (the root of the tree) and point any nodes without outgoing edges to it to get an aromatic vector field; to invert, cut off the root.</footnote>
Putting this together, the rooted trees with one marked node are enumerated by the following (sequence A000107),

$$w(z) := t(z) + t(z)^2 + t(z)^3 + \cdots = t(z)/(1 - t(z)).$$

The connected scalar aromas with one marked node are enumerated by $w(z)+w(z)^2 = t(z)/(1-t(z))^2$ (sequence A038002). The scalar aromas with one marked node are enumerated by $a(z)t(z)/(1-t(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)/(1-t(z))^2$.

For type (ii), combining the components of an unmarked form with one fewer trees (i.e. an element of $\Omega_{n-1}$), enumerated by $ub_0(u, z)$, and a marked rooted tree, enumerated by $t(z)/(1-t(z))$, gives the contribution from type (ii) forms as $ub_0(u, z)t(z)/(1-t(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 $\textcircled{1}$). 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-1}| = |\Omega_{0,1}^{N-1}|$ and this gives the result for $c_1(z)$. $\square$

Note that

$$\Omega_{1,1}^{N-1} \cong \Omega_{0,1}^{N-1} \cong \mathcal{I}_1^N$$

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:

$$ \begin{align*} \left( \sum_{n=0}^{\infty} (-1)^n |\Omega_{n,1}| \right) - |\mathcal{I}_1| &= b_1(-1, z) - c_1(z) \\ &= b_0(-1, z) \frac{t(z)^2}{(1-t(z))^2} - \frac{a(z)t(z)z}{(1-t(z))^2} \\ &= \frac{a(z)z}{t(z)} \frac{t(z)^2}{(1-t(z))^2} - \frac{a(z)t(z)z}{(1-t(z))^2} \\ &= 0. \end{align*} $$

The following result extends Theorem 4.9 to the divergence-free case.

**Theorem 4.4.** Let $\tilde{b}_p(u, z)$, $\tilde{c}_p(z)$, and $\tilde{s}(z)$ be the divergence-free analogues of the generating functions $b_p(z)$, $c_p(z)$, and $s(z)$. Then

$$ \begin{align*} \tilde{b}_0(u, z) &= \frac{zb_0(u, z)}{t(z)}, \\ \tilde{b}_1(u, z) &= zb_0(u, z) \frac{t(z) + u - ut(z))}{(1 - t(z))^2}, \\ \tilde{c}_1(z) &= zc_1(z), \\ \tilde{s}(z) &= z + \frac{zs(z)}{t(z)} = z + \frac{z(a(z)t(z) - \dot{a}(z))}{t(z)}. \end{align*} $$
<table><caption>Table 2: Dimensions of the bottom two rows of the augmented aromatic bicomplex for orders 1 to 9.</caption><thead><tr><th rowspan="3">N</th><th colspan="5">First row: |Ω<sub>n</sub><sup>N</sup>|</th><th colspan="5">Second row: |Ω<sub>n,1</sub><sup>N</sup>| |I<sub>1</sub><sup>N</sup>|</th></tr><tr><th colspan="5">n</th><th colspan="5">n</th></tr><tr><th>4</th><th>3</th><th>2</th><th>1</th><th>0</th><th>4</th><th>3</th><th>2</th><th>1</th><th>0</th></tr></thead><tbody><tr><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td></tr><tr><td>2</td><td>0</td><td>0</td><td>0</td><td>2</td><td>3</td><td>0</td><td>0</td><td>1</td><td>4</td><td>4</td><td>1</td></tr><tr><td>3</td><td>0</td><td>0</td><td>1</td><td>6</td><td>7</td><td>0</td><td>0</td><td>4</td><td>15</td><td>15</td><td>4</td></tr><tr><td>4</td><td>0</td><td>0</td><td>3</td><td>16</td><td>19</td><td>0</td><td>1</td><td>16</td><td>52</td><td>52</td><td>15</td></tr><tr><td>5</td><td>0</td><td>0</td><td>11</td><td>45</td><td>47</td><td>0</td><td>5</td><td>57</td><td>175</td><td>175</td><td>52</td></tr><tr><td>6</td><td>0</td><td>2</td><td>33</td><td>121</td><td>130</td><td>0</td><td>22</td><td>197</td><td>571</td><td>571</td><td>175</td></tr><tr><td>7</td><td>0</td><td>7</td><td>102</td><td>338</td><td>343</td><td>2</td><td>85</td><td>654</td><td>1838</td><td>1838</td><td>571</td></tr><tr><td>8</td><td>0</td><td>29</td><td>298</td><td>929</td><td>951</td><td>11</td><td>310</td><td>2137</td><td>5834</td><td>5834</td><td>1838</td></tr><tr><td>9</td><td>1</td><td>99</td><td>878</td><td>2598</td><td>2615</td><td>53</td><td>1078</td><td>6859</td><td>18363</td><td>18373</td><td>5834</td></tr></tbody></table>

*Proof.* In each case we need to enumerate the non-self-looped elements of the aromatic bicomplex. The loops occur only in the scalar components. Recall that the non-self-looped scalars, $\tilde{\Omega}_{0,1}$, are enumerated by $\bar{a}(z) = za(z)/t(z)$. This gives the result for $\tilde{b}_0$ and, using Theorem 2.10, $\tilde{s}(z)$.

For the second row, we first consider $\tilde{\Omega}_{0,1}^N$, the scalar aromas with $N$ nodes, no self-loops, and one marked node indicated by ①. These are enumerated by the self-functions of $\{1, \dots, N-1\}$ with one marked node. That is, $\tilde{\Omega}_{0,1}^N \cong \Omega_{0,1}^{N-1}$. The construction is as follows. Consider an element of $\Omega_{0,1}^{N-1}$, that is, a directed graph with $N-1$ 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 self-loops to the new node; then, move the mark to the new node. This gives a marked non-self-looped scalar aroma of degree $N$. The process is invertible: given a marked scalar with no self-loops, 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 $\tilde{\Omega}_{0,1}^N$ has generating function $za(z)t(z)/(1-t(z))^2$.

For the rest of the second row, $\tilde{\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 divergence-free case is that the scalars must be non-self-looped, as just enumerated. For type (ii), we combine the three components of non-self-looped scalars (enumerated by $a(z)$), $n-1$ unmarked rooted trees (enumerated by $ub_0(u, z)/a(z)$), and one marked rooted tree (enumerated by $t(z)/(1-t(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 $\tilde{\mathcal{I}}_1^N = I(\tilde{\Omega}_{0,p}^N)$ are linear combinations of scalar aromas with $N$ nodes, one marked leaf, and no self-loops. 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 $\tilde{c}_1(z)$. $\square$

Note that we now have five isomorphic spaces,

$$ \Omega_{1,1}^{N-1} \cong \Omega_{0,1}^{N-1} \cong \mathcal{I}_1^N \cong \tilde{\mathcal{I}}_1^{N+1} \cong \tilde{\Omega}_{0,1}^N. $$

Each space is enumerated by the self-functions of $\{1, \dots, N-1\}$ with one marked node (and generating function $za(z)t(z)/(1-t(z))^2$), but in a different way in each case.

We remark that the second row of the divergence-free 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: \Omega_1 \to \Omega_0$ and we describe the image and the kernel of $d_H$ and $d_H^*$. We recall that $d_H^*: \Omega_0^* \to \Omega_1^*$ is the dual map of $d_H$. We described the dimension of $\Psi = \text{Ker } d_H$ in Subsection 4.2. The following result describes the dimensions of $\text{Im } d_H$, $\text{Ker } d_H^*$, and $\text{Im } d_H^*$.

**Theorem 4.5.** The scalar divergences $\text{Im } d_H$ have dimension $|\hat{\Omega}_0|$. The conditions $\text{Ker } d_H^*$ that a scalar must satisfy to be a divergence have dimension $|\bar{\Omega}_0|$. The conditions $\text{Im } d_H^*$ that a vector must satisfy to be divergence-free have dimension $|\hat{\Omega}_0|$.

*Proof.* Recall the fundamental theorem of linear algebra for a linear map $A: V \to W$:

$$|\text{Im } A| = |\text{Im } A^*| = \text{rank}A, \quad |\text{Im } A| + |\text{Ker } A| = |V|, \quad |\text{Im } A^*| + |\text{Ker } A^*| = |W|.$$

The space $\text{Im } A^* = \text{Ann}(\text{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 $\text{Ker } A^* = \text{Ann}(\text{Im } A)$ is the annihilator of the image of $A$. Choosing $A = d_H$ gives the result. $\square$

**Remark 4.6.** Consider the aromatic tree $\gamma \setminus 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 \setminus e)$, and $m_2(\gamma, e)$ the number of edges $\hat{e} \in E$ of $\gamma$ such that $\gamma \setminus \hat{e} = \gamma \setminus e$. Then, for $\gamma \in \Omega_0$, $d_H^* \gamma^*$ satisfies

$$d_H^* \gamma^* = \sum_{e \in E} \frac{m_1(\gamma, e)}{m_2(\gamma, e)} (\gamma \setminus e)^*.$$

We now construct bases of $\text{Ker } d_H$, $\text{Im } d_H$, $\text{Ker } d_H^*$, and $\text{Im } d_H^*$. We start with the basis of the solenoidal forms.

**Theorem 4.7.** Let $\varphi: \Omega_1 \to \hat{\Omega}_0$ be defined by attaching a self-loop to the root, extending by linearity. Define any total order on $\mathcal{T}$, then a basis of $\Psi = \text{Ker } d_H$ is

$$\mathcal{B}_{\Psi} = \{d_H(\phi\varphi(t_2) \dots \varphi(t_{l-1})\varphi(t_{l+1}) \dots \varphi(t_k)t_1 \wedge t_l), \quad \phi \in \mathcal{F}_0 \cup \{\emptyset\}, \quad t_1 < \dots < t_k \in \mathcal{T}\}.$$

*Proof.* The map $\varphi$ is surjective: given any self-looped scalar, removing one of the self-loops gives a pre-image. Therefore, $|\text{Ker } \varphi| = |\text{Ker } d_H|$. Let us first determine $\text{Ker } \varphi$. Consider a self-looped scalar whose distinct self-looped connected components are $\varphi(t_1), \dots, \varphi(t_k)$ for trees $t_1 < \dots < t_k$; it can be written $\varphi(t_1) \dots \varphi(t_k)$ where $\varphi$ is a scalar. Its distinct preimages under $\varphi$ are $\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 $k-1$ and we consider

$$\{\varphi(t_1\varphi(t_2) - \varphi(t_1)t_2)\varphi(t_3) \dots \varphi(t_k), \dots, \varphi(t_1\varphi(t_k) - \varphi(t_1)t_k)\varphi(t_2) \dots \varphi(t_{k-1})\}$$

as a basis of $\text{Ker } \varphi$. We now map $\text{Ker } \varphi$ to $\text{Ker } d_H$ by

$$\phi_l(t_1\varphi(t_l) - \varphi(t_1)t_l) \mapsto d_H(\phi_l t_1 \wedge t_l), \quad \text{where } \phi_l = \varphi(t_2) \dots \varphi(t_{l-1})\varphi(t_{l+1}) \dots \varphi(t_k),$$

extending by linearity. From exactness, the map is surjective. As $|\text{Ker } \varphi| = |\text{Ker } d_H|$, it is an isomorphism. $\square$
**Remark 4.8.** One could wonder whether the set {$d_H \wedge \gamma, \gamma \in \mathcal{F}_2$} is a basis of $\text{Ker}(d_H)$. This is not the case in general. For $N = 6$, one finds for instance the following identity

$$d_H(\bigcirc \cdots \bigcirc + \bigcirc \cdots \bigcirc + \bigcirc \cdots \bigcirc + \cdots \bigcirc \cdots \bigcirc + \cdots \bigcirc \cdots \bigcirc + \cdots \bigcirc \cdots \bigcirc + \cdots \bigcirc \cdots \bigcirc + \cdots \bigcirc \cdots \bigcirc + \cdots \bigcirc \cdots \bigcirc + \cdots \bigcirc \cdots \bigcirc + \cdots \bigcirc \cdots \bigcirc) = 0.$$

The following result shows that the divergences are a graph over the self-looped scalars.

**Theorem 4.9.** For $\alpha \in \check{\Omega}_0$, let $k(\alpha)$ be the number of self-loops in $\alpha$, and $\rho(\alpha)$ be the non-self-looped scalar obtained from $\alpha$ as the sum of the redirection of all 1-loops to other nodes in all possible ways. Then the map

$$\check{\Omega}_0 \to \text{Im } d_H, \quad \alpha \mapsto \alpha + (-1)^{k(\alpha)-1} \rho(\alpha)$$

is an isomorphism, and generates a basis of $\text{Im } d_H$.

*Proof.* Let $\check{V}$ be the set of nodes with self-loops of $\alpha$. From Proposition 3.6, we deduce that for $v \in \check{V}$, $\alpha - \mathcal{E}_v\alpha \in \text{Im}(d_H)$. If we have two nodes $v, w \in \check{V}$, then $\alpha - \mathcal{E}_v\alpha \in \text{Im}(d_H)$ and $\mathcal{E}_v\alpha - \mathcal{E}_w\mathcal{E}_v\alpha \in \text{Im}(d_H)$, so that $\alpha - \mathcal{E}_w\mathcal{E}_v\alpha \in \text{Im}(d_H)$. By applying iteratively this process, we find that

$$\alpha - \prod_{v \in \check{V}} \mathcal{E}_v \alpha = \alpha + (-1)^{k(\alpha)-1} \rho(\alpha) \in \text{Im}(d_H).$$

As each self-looped scalar appears once in the image, the map is injective. The map is an isomorphism as the domain and codomain have the same dimension. $\square$

**Remark 4.10.** The operation $\rho$ in Theorem 4.9 corresponds to removing all self-loops in $\alpha$ by repeated integration by parts, as illustrated in the following example on elementary differentials:

$$\begin{aligned} f_i^i f_j^j &= (f^i f_j^j)_i - f^i f_{ij}^j \\ &= (f^i f_j^j)_i - (f^i f_i^j)_j + f_j^i f_i^j \\ &= (f^i f_j^j - f^j f_i^j)_i + f_j^i f_i^j. \end{aligned}$$

That is, $f_i^i f_j^j - f_j^i f_i^j$ is a divergence. We describe this comparison with integration by parts further in Subsection 4.4.

**Corollary 4.11.** No non-trivial combination of non-self-looped scalars in $\Omega_0$ is a divergence.

**Corollary 4.12.** The conditions to be a divergence $\text{Ker } d_H^*$ are a graph over the dual of the non-self-looped forms in $\Omega_0$.

In the following theorem, this graph is realized explicitly.

**Theorem 4.13.** Let $\hat{E} \subset E$ be a set of edges of the scalar aroma $\beta \in \Omega_0$ and let $\beta \setminus \hat{E}$ be $\beta$ with edges $\hat{E}$ replaced by self-loops. Let $m(\beta, \hat{E})$ be the number of ways that redirecting self-loops of $\beta \setminus \hat{E}$ results in $\beta$. Let $\pi: \Omega_0 \to \Omega_0$ be defined by

$$\pi(\beta) = \sum_{\hat{E} \subset E} (-1)^{|\hat{E}|} m(\beta, \hat{E}) \beta \setminus \hat{E}$$

Then

$$\{\pi(\beta), \beta \in \bar{\Omega}_0^*\}$$

is a basis of $\text{Ann}(\text{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) : x \in \mathbb{R}^n\} \subset \mathbb{R}^n \times \mathbb{R}^m$ are $y - Ax = 0$, $A \in \mathbb{R}^{m \times n}$. A basis for these conditions is given by the rows of $y - Ax$, where Theorem 4.9 gives the columns of $A$. That is, for each non-self-looped scalar $\beta$ we need to determine the coefficient of $\beta$ in $\rho(\alpha)$ for each self-looped scalar $\alpha$. This is the expression for $\pi$: the term $\hat{E} = \emptyset$ gives $\beta$, and the terms from non-empty sets of edges $\hat{E}$ give the $\alpha$'s that can give rise to $\beta$. $\square$

Finally, we present a basis of 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

$$ \{d_H^* \phi^*, \phi \in \hat{\Omega}_0\} $$

is a basis of Im $d_H^*$, the conditions that a form in $\Omega_0$ must satisfy to be divergence free.

Proof. For any linear map $A: V \to W$, $\langle A^*w^*, v \rangle = \langle w^*, Av \rangle = 0$ for all $w^* \in W^*$ when $v \in$ Ker $A$. Thus, divergence-free 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^*|_{\hat{\Omega}_0^*}$ is injective.

Let $\text{pr}: \Omega_0 \to \hat{\Omega}_0$ be the natural projection to the self-looped scalars. From Theorem 4.9, the divergences form a graph over the self-looped scalars. That is, $\text{pr} \circ d_H$ is surjective. Therefore its dual $d_H^*|_{\hat{\Omega}_0^*}$ is injective. $\square$

## 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 [36, 37] on exotic aromatic B-series (see also [8]). 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 1-loop. 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 $\hat{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.

**Homotopy operator $\hat{h}_H$**

Given $\gamma \in \Omega_{0,p}^N$, initialize $\hat{h}_H\gamma = 0$ and $\hat{\gamma} = \gamma - |\gamma|^{-1}\varepsilon\gamma$.

while there is a 1-loop in a forest $\tau$ of $\hat{\gamma}$ on a vertex $v$ do

$\hat{h}_H\gamma \leftarrow \hat{h}_H\gamma + a(\tau)\theta_v(\tau)$,

$\hat{\gamma} \leftarrow \hat{\gamma} - a(\tau)d_H\theta_v(\tau)$.

end while

return $\hat{h}_H\gamma$

Note that each iteration in the algorithm reduces the number of 1-loops 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 [36], as there is an extra term involved in the integration by parts process.
<table><caption>Table 3: Comparison of the horizontal homotopy operators on Ω₀ for the first orders.</caption><thead><tr><th>γ ∈ Ω<sub>0</sub></th><th>h<sub>H</sub>γ</th><th>ĥ<sub>H</sub>γ</th></tr></thead><tbody><tr><td>○</td><td>•</td><td>•</td></tr><tr><td>○</td><td>0</td><td>0</td></tr><tr><td>○</td><td>•</td><td>•</td></tr><tr><td>○ ○</td><td>○ •</td><td>○ •</td></tr><tr><td>○</td><td>1/6 ○ • + 1/6 ○ • - 1/6 ○ • - 1/6 ○ •</td><td>1/3 ○ • - 1/3 ○ •</td></tr><tr><td>○</td><td>1/6 ○ • + 1/6 ○ • - 1/6 ○ • - 1/6 ○ •</td><td>1/3 ○ • - 1/3 ○ •</td></tr><tr><td>○</td><td>○</td><td>○</td></tr><tr><td>○</td><td>2/3 ○ • + 1/3 ○ •</td><td>○ • + 2/3 ○ • - 2/3 ○ •</td></tr><tr><td>○</td><td>0</td><td>0</td></tr><tr><td>○ ○</td><td>1/3 ○ • + 2/3 ○ •</td><td>○ • + 2/3 ○ • - 2/3 ○ •</td></tr><tr><td>○ ○</td><td>○</td><td>○</td></tr></tbody></table>

An aromatic B-series method can be seen as the exact solution of the modified ODE (1.3), and the modified flow is given by the aromatic B-series $B(b)$ satisfying $B(a) = B(b) \triangleright B(e)$. The operation $\triangleright$ is the substitution of B-series and $b: \tilde{\Omega}_1 \to \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 B-series coefficients [13, 17, 6]. The map $e$ is the coefficient of the exact flow of (1.1). Its expression is given for instance in [23, 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 [45] is the following: can we find an aromatic B-series method such that $d_H B(b) = 0$, that is, such that the modified B-series $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 B-series instead of standard B-series comes from the following negative result.

**Theorem 4.16** ([18, 31]). *The solenoidal combinations of rooted trees satisfy*

$$ \mathrm{Span}(\mathcal{T}) \cap \Psi = \emptyset, \quad \mathrm{Span}(\mathcal{T}) \cap \tilde{\Psi} = \mathrm{Span}(\bullet). $$

*In particular, the only volume-preserving consistent B-series 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 self-looped scalars (see Theorem 4.9), $d_H v = 0$ if the coefficient of each self-looped scalar in $d_H v$ is zero. But $\mathrm{pr}d_H v = \sum_i c_i \theta(t_i)$, where $\mathrm{pr}: \Omega_0 \to \tilde{\Omega}_0$ is the natural projection to the self-looped scalars, and the self-looped scalars $\theta(t_i)$ are linearly independent.
Following Theorem 4.16, one is interested in finding a class of non-trivial volume-preserving consistent aromatic B-series 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 volume-preserving integrator.

**Theorem 4.17.** If $B(a)$ is the aromatic B-series of a consistent volume-preserving integrator, then there exists $\eta \in \tilde{\Omega}_2$ such that the modified flow is a B-series of the form $B(\bullet + d_H\eta)$ and $B(a)$ is given by the substitution

$$B(a) = (\bullet + d\eta) \triangleright B(e).$$

More precisely, there exists a coefficient map $\alpha: \tilde{\Omega}_2 \to \mathbb{R}$ such that a is given by

$$a = (\bullet^* + \alpha d_H^* A_\sigma) \star e, \quad \text{where} \quad A_\sigma \gamma = \frac{1}{\sigma(\gamma)} \gamma. \qquad (4.6)$$

For the first orders, the B-series of a volume-preserving aromatic B-series method has the form:

$$B(a) = \bullet + \frac{1}{2}\bullet + \frac{1}{6}\bullet + \left(\frac{1}{6} - \frac{1}{2}\alpha(\wedge\bullet\cdot)\right)\bullet\cdot + \frac{1}{2}\alpha(\wedge\bullet\cdot)\bullet\cdot \\ + \frac{1}{24}\bullet\cdot + \left(\frac{1}{8} - \frac{1}{2}\alpha(\wedge\bullet\cdot\cdot) - \alpha(\wedge\bullet\cdot\cdot) - \frac{1}{2}\alpha(\wedge\bullet\cdot)\right)\bullet\cdot\cdot \\ + \left(\frac{1}{24} - \frac{1}{2}\alpha(\wedge\bullet\cdot\cdot\cdot) + \frac{1}{2}\alpha(\wedge\bullet\cdot\cdot\cdot) - \frac{1}{4}\alpha(\wedge\bullet\cdot)\right)\bullet\cdot\cdot\cdot + \frac{1}{2}\alpha(\wedge\bullet\cdot\cdot\cdot)\bullet\cdot\cdot + \frac{1}{2}\alpha(\wedge\bullet\cdot\cdot\cdot)\bullet\cdot\cdot\bullet \\ + \left(\frac{1}{24} - \frac{1}{2}\alpha(\wedge\bullet\cdot\cdot\cdot) - \frac{1}{4}\alpha(\wedge\bullet\cdot\cdot)\right)\bullet\cdot\cdot\cdot + \left(\frac{1}{2}\alpha(\wedge\bullet\cdot\cdot\cdot) + \alpha(\wedge\bullet\cdot\cdot\cdot) + \frac{1}{2}\alpha(\wedge\bullet\cdot\cdot)\right)\bullet\cdot\cdot\cdot + \dots$$

Note that the coefficients of the bamboo trees (or tall trees) $\mathcal{BT} = \{\bullet, \bullet, \bullet, \dots\}$ coincide with the ones of the exact flow. This fact has been noticed for standard B-series in particular in [33] (see also [23, 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 $\bullet$. Indeed $d_H^*$ vanishes on $\mathcal{BT}$ 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{BT}$, and no bamboo tree appears in these forms. According to the identity (4.6), the B-series of an aromatic volume-preserving method has to coincide with the exact flow $B(e)$ on the bamboo trees $\mathcal{BT}$. A Runge-Kutta integrator cannot be exact on all bamboo trees. As the aromatic forests represent different elementary differentials (see [31]), any aromatic integrator that reduces to a standard Runge-Kutta integrator when sending the aromas to zero cannot be volume-preserving. $\square$

The methodology proposed in [6, Sect. 7] to obtain volume-preservation of high order and the approach in [45, 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 volume-preserving 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, [4, 30, 39]) when sending the aromas to zero. This calls for future works that study the substitution law and the variational bi-complex directly on aromatic exponential B-series, in order to find a volume-preserving aromatic B-series 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 divergence-free case. To this end, we introduced the Euler operators, the homotopy operators, as well as an augmented bicomplex. The algebraic properties we proved have concrete consequences on the numerical analysis of volume-preserving integrators. They allow to describe completely the solenoidal forms and the B-series of an aromatic volume-preserving method. In particular, we proved that there are no volume-preserving 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 B-series, in order to find an aromatic exponential volume-preserving 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 volume-preserving 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 [36]. An exact formula would greatly benefit the creation of high-order methods for solving ergodic SDEs.

There is a considerable literature on the variational bicomplex and the De Rham cohomology (see [1] 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 divergence-free case (see Subsection 4.1). Two major applications of the variational bicomplex are the 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.

**Acknowledgements.** 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.

# References

[1] I. M. Anderson. The variational bicomplex. Unpublished, 1989.

[2] I. M. Anderson. Introduction to the variational bicomplex. In *Mathematical aspects of classical field theory (Seattle, WA, 1991)*, volume 132 of *Contemp. Math.*, pages 51–73. Amer. Math. Soc., Providence, RI, 1992.

[3] I. M. Anderson and T. Duchamp. On the existence of global variational principles. *Amer. J. Math.*, 102(5):781–868, 1980.

[4] H. Berland, B. Owren, and B. Skaflestad. B-series and order conditions for exponential integrators. *SIAM J. Numer. Anal.*, 43(4):1715–1727, 2005.
[5] S. Blanes and F. Casas. On the necessity of negative coefficients for operator splitting schemes of order higher than two. *Appl. Numer. Math.*, 54(1):23–37, 2005.

[6] G. Bogfjellmo. Algebraic structure of aromatic B-series. *J. Comput. Dyn.*, 6(2):199–222, 2019.

[7] G. Bogfjellmo, E. Celledoni, R. McLachlan, B. Owren, and R. Quispel. Using aromas to search for preserved measures and integrals in Kahan’s method. *Submitted*, 2022.

[8] E. Bronasco. Exotic B-series and S-series: algebraic structures and order conditions for invariant measure sampling. *Submitted*, 2022.

[9] K. Burrage and P. M. Burrage. High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations. *Appl. Numer. Math.*, 22(1-3):81–101, 1996. Special issue celebrating the centenary of Runge-Kutta methods.

[10] J. C. Butcher. An algebraic theory of integration methods. *Math. Comp.*, 26:79–106, 1972.

[11] J. C. Butcher. *Numerical methods for ordinary differential equations*. John Wiley & Sons, Ltd., Chichester, third edition, 2016.

[12] J. C. Butcher. *B-series: algebraic analysis of numerical methods*. Springer, 2021.

[13] D. Calaque, K. Ebrahimi-Fard, and D. Manchon. Two interacting Hopf algebras of trees: a Hopf-algebraic approach to composition and substitution of B-series. *Adv. in Appl. Math.*, 47(2):282–308, 2011.

[14] E. Celledoni, C. Evripidou, D. I. McLaren, B. Owren, G. R. W. Quispel, and B. K. Tapley. Detecting and determining preserved measures and integrals of birational maps. *Journal of Computational Dynamics*, 2022.

[15] E. Celledoni, R. I. McLachlan, D. I. McLaren, B. Owren, and G. R. W. Quispel. Integrability properties of Kahan’s method. *J. Phys. A*, 47(36):365202, 20, 2014.

[16] E. Celledoni, R. I. McLachlan, B. Owren, and G. R. W. Quispel. Geometric properties of Kahan’s method. *J. Phys. A*, 46(2):025201, 12, 2013.

[17] P. Chartier, E. Hairer, and G. Vilmart. Algebraic structures of B-series. *Found. Comput. Math.*, 10(4):407–427, 2010.

[18] P. Chartier and A. Murua. Preserving first integrals and volume forms of additively split systems. *IMA J. Numer. Anal.*, 27(2):381–405, 2007.

[19] K. Debrabant and A. Kværnø. B-series analysis of stochastic Runge-Kutta methods that use an iterative scheme to compute their internal stage values. *SIAM J. Numer. Anal.*, 47(1):181–203, 2008/09.

[20] P. Flajolet and R. Sedgewick. *Analytic combinatorics*. Cambridge University Press, Cambridge, 2009.

[21] G. Fløystad, D. Manchon, and H. Z. Munthe-Kaas. The universal pre-Lie-Rinehart algebras of aromatic trees. In *Geometric and harmonic analysis on homogeneous spaces and applications*, volume 366 of Springer Proc. Math. Stat., pages 137–159. Springer, Cham, 2021 ©2021.
[22] M. Gubinelli. Ramification of rough paths. *J. Differential Equations*, 248(4):693–721, 2010.

[23] E. Hairer, C. Lubich, and G. Wanner. *Geometric numerical integration*, volume 31 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2006. Structure-preserving algorithms for ordinary differential equations.

[24] E. Hairer and G. Wanner. On the Butcher group and general multi-value methods. *Computing (Arch. Elektron. Rechnen)*, 13(1):1–15, 1974.

[25] M. Hairer and D. Kelly. Geometric versus non-geometric rough paths. *Ann. Inst. Henri Poincaré Probab. Stat.*, 51(1):207–251, 2015.

[26] W. Hereman, M. Colagrosso, R. Sayers, A. Ringler, B. Deconinck, M. Nivala, and M. Hickman. Continuous and discrete homotopy operators and the computation of conservation laws. In *Differential equations with symbolic computation*, Trends Math., pages 255–290. Birkhäuser, Basel, 2005.

[27] W. Hereman, B. Deconinck, and L. D. Poole. Continuous and discrete homotopy operators: a theoretical approach made concrete. *Math. Comput. Simulation*, 74(4-5):352–360, 2007.

[28] R. Hirota and K. Kimura. Discretization of the Euler top. *J. Phys. Soc. Japan*, 69(3):627–630, 2000.

[29] R. Hirota and K. Kimura. Discretization of the Lagrange top. *J. Phys. Soc. Japan*, 69(10):3193–3199, 2000.

[30] M. Hochbruck and A. Ostermann. Explicit exponential Runge-Kutta methods for semilinear parabolic problems. *SIAM J. Numer. Anal.*, 43(3):1069–1090, 2005.

[31] A. Iserles, G. R. W. Quispel, and P. S. P. Tse. B-series methods cannot be volume-preserving. *BIT Numer. Math.*, 47(2):351–378, 2007.

[32] W. Kahan. Unconventional numerical methods for trajectory calculations. *Unpublished lecture notes*, 1:13, 1993.

[33] F. Kang and Z. J. Shang. Volume-preserving algorithms for source-free dynamical systems. *Numer. Math.*, 71(4):451–463, 1995.

[34] Y. Komori, T. Mitsui, and H. Sugiura. Rooted tree analysis of the order conditions of ROW-type scheme for stochastic differential equations. *BIT Numer. Math.*, 37(1):43–66, 1997.

[35] A. Kushner, V. Lychagin, and V. Rubtsov. *Contact geometry and non-linear differential equations*, volume 101 of *Encyclopedia of Mathematics and its Applications*. Cambridge University Press, Cambridge, 2007.

[36] A. Laurent and G. Vilmart. Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs. *Math. Comp.*, 89(321):169–202, 2020.

[37] A. Laurent and G. Vilmart. Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds. *Found. Comput. Math.*, 22(3):649–695, 2022.

[38] J. M. Lee. *Introduction to smooth manifolds*, volume 218 of *Graduate Texts in Mathematics*. Springer, New York, second edition, 2013.
[39] V. T. Luan and A. Ostermann. Exponential B-series: the stiff case. *SIAM J. Numer. Anal.*, 51(6):3431–3445, 2013.

[40] R. I. McLachlan, K. Modin, H. Munthe-Kaas, and O. Verdier. B-series methods are exactly the affine equivariant methods. *Numer. Math.*, 133(3):599–622, 2016.

[41] R. I. McLachlan, K. Modin, H. Munthe-Kaas, and O. Verdier. Butcher series: a story of rooted trees and numerical methods for evolution equations. *Asia Pac. Math. Newsl.*, 7(1):1–11, 2017.

[42] R. I. McLachlan, H. Z. Munthe-Kaas, G. R. W. Quispel, and A. Zanna. Explicit volume-preserving splitting methods for linear and quadratic divergence-free vector fields. *Found. Comput. Math.*, 8(3):335–355, 2008.

[43] R. I. McLachlan and G. R. W. Quispel. Splitting methods. *Acta Numer.*, 11:341–434, 2002.

[44] R. Merris. *Combinatorics*. Wiley-Interscience Series in Discrete Mathematics and Optimization. Wiley-Interscience [John Wiley & Sons], New York, second edition, 2003.

[45] H. Munthe-Kaas and O. Verdier. Aromatic Butcher series. *Found. Comput. Math.*, 16(1):183–215, 2016.

[46] OEIS Foundation Inc. The On-Line Encyclopedia of Integer Sequences, 2022.

[47] P. J. Olver. *Applications of Lie groups to differential equations*, volume 107 of *Graduate Texts in Mathematics*. Springer-Verlag, New York, second edition, 1993.

[48] P. J. Olver and C. Shakiban. A resolution of the Euler operator. I. *Proc. Amer. Math. Soc.*, 69(2):223–229, 1978.

[49] M. Petrera, A. Pfadler, and Y. B. Suris. On integrability of Hirota-Kimura type discretizations. *Regul. Chaotic Dyn.*, 16(3-4):245–289, 2011.

[50] D. Poole and W. Hereman. The homotopy operator method for symbolic integration by parts and inversion of divergences with applications. *Appl. Anal.*, 89(4):433–455, 2010.

[51] G. R. W. Quispel. Volume-preserving integrators. *Phys. Lett. A*, 206(1-2):26–30, 1995.

[52] A. Rößler. Stochastic Taylor expansions for the expectation of functionals of diffusion processes. *Stochastic Anal. Appl.*, 22(6):1553–1576, 2004.

[53] C. Scovel. Symplectic numerical integration of Hamiltonian systems. In *The geometry of Hamiltonian systems (Berkeley, CA, 1989)*, volume 22 of *Math. Sci. Res. Inst. Publ.*, pages 463–496. Springer, New York, 1991.

[54] C. Shakiban. A resolution of the Euler operator. II. *Math. Proc. Cambridge Philos. Soc.*, 89(3):501–510, 1981.

[55] Z. J. Shang. Construction of volume-preserving difference schemes for source-free systems via generating functions. *J. Comput. Math.*, 12(3):265–272, 1994.

[56] T. Tsujishita. On variation bicomplexes associated to differential equations. *Osaka Math. J.*, 19(2):311–363, 1982.
[57] W. M. Tulczyjew. The Euler-Lagrange resolution. In *Differential geometrical methods in mathematical physics (Proc. Conf., Aix-en-Provence/Salamanca, 1979)*, volume 836 of *Lecture Notes in Math.*, pages 22–48. Springer, Berlin, 1980.

[58] A. M. Vinogradov. A spectral sequence that is connected with a nonlinear differential equation, and the algebraic-geometry foundations of the Lagrange field theory with constrain. *Dokl. Akad. Nauk SSSR*, 238(5):1028–1031, 1978.

[59] A. M. Vinogradov. The *C*-spectral sequence, Lagrangian formalism, and conservation laws. I. The linear theory. *J. Math. Anal. Appl.*, 100(1):1–40, 1984.

[60] A. M. Vinogradov. The *C*-spectral sequence, Lagrangian formalism, and conservation laws. II. The nonlinear theory. *J. Math. Anal. Appl.*, 100(1):41–129, 1984.

[61] H. Xue and A. Zanna. Explicit volume-preserving splitting methods for polynomial divergence-free vector fields. *BIT Numer. Math.*, 53(1):265–281, 2013.

# Appendices

## A First solenoidal forms

We write the generators of the solenoidal forms $\tilde{\Psi}^N$ in the divergence-free 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 \tilde{\mathcal{F}}_2^N$, and by adding the trivial tree $\cdot$. For $N \le 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 divergence-free 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.
<table><caption>Table 4: Generators of the solenoidal forms $\tilde{\Psi}^N$ for the first orders N.</caption><thead><tr><th>N</th><th>γ ∈ F̃<sub>2</sub><sup>N</sup></th><th>Solenoidal form 2d<sub>H</sub> ∧ γ</th></tr></thead><tbody><tr><td>1</td><td>...</td><td>...</td></tr><tr><td>3</td><td>...</td><td>...</td></tr><tr><td rowspan="2">4</td><td>...</td><td>...</td></tr><tr><td>...</td><td>...</td></tr><tr><td rowspan="2">5</td><td>...</td><td>...</td></tr><tr><td>...</td><td>...</td></tr><tr><td rowspan="2">6</td><td>...</td><td>...</td></tr><tr><td>...</td><td>...</td></tr><tr><td rowspan="2">7</td><td>...</td><td>...</td></tr><tr><td>...</td><td>...</td></tr><tr><td rowspan="2">8</td><td>...</td><td>...</td></tr><tr><td>...</td><td>...</td></tr><tr><td rowspan="2">9</td><td>...</td><td>...</td></tr><tr><td>...</td><td>...</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 4: The augmented aromatic bicomplex for $N = 1$ and $N = 2$. The wedges are omitted for conciseness.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: The augmented aromatic bicomplex for $N = 3$. The wedges are omitted for conciseness.</figcaption></figure>