Adrien Laurent
email: adrien.laurent@uib.no.hans.munthe-kaas@uib.no.

Robert I Mclachlan
email: r.mclachlan@massey.ac.nz.

Hans Z Munthe-Kaas

Olivier Verdier
email: olivier.verdier@hvl.no.

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

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

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

Introduction

Let f : R d Ñ R d be a smooth Lipschitz vector field and let y : r0, T s Ñ R d be the solution to the ordinary differential equation y I ptq f pyptqq, t s0, T r, yp0q y 0 .

(1.1) If the vector field f is divergence-free, that is, if divpf q 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 λ, for any t ¡ 0, the flow φ t of the ODE (1.1) satisfies λpφ t pDqq λpDq. Divergence-free vector fields appear in a variety of concrete dynamical systems, for instance in fluid dynamics, meteorology or molecular dynamics. To integrate such systems, it is fundamental to use numerical integrators that also preserve the volume. We discretize the time interval r0, T s into N 1 equidistant steps t n nh, with h the time stepsize, and we choose a one-step integrator y n 1 Φpy n , hq.

(1.2)

For large classes of integrators, the standard methodology to study volume-preservation uses backward error analysis (see, for instance, the textbook [START_REF] Hairer | Geometric numerical integration[END_REF]). The integrator (1.2) can be interpreted formally as the exact solution of a modified ODE r y I ptq r f pr yptqq, (1.3) where the modified vector field r f typically depends on f and its derivatives. Then, the integrator is volume-preserving if and only if divp r f q 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 [START_REF] Hairer | Geometric numerical integration[END_REF][START_REF] Chartier | Preserving first integrals and volume forms of additively split systems[END_REF], on splitting methods [START_REF] Kang | Volume-preserving algorithms for source-free dynamical systems[END_REF][START_REF] Mclachlan | Splitting methods[END_REF][START_REF] Mclachlan | Explicit volumepreserving splitting methods for linear and quadratic divergence-free vector fields[END_REF][START_REF] Xue | Explicit volume-preserving splitting methods for polynomial divergence-free vector fields[END_REF], or on generating functions [START_REF] Scovel | Symplectic numerical integration of Hamiltonian systems[END_REF][START_REF] Shang | Construction of volume-preserving difference schemes for source-free systems via generating functions[END_REF][START_REF] Quispel | Volume-preserving integrators[END_REF]. For quadratic differential equations, we mention the works [START_REF] Petrera | On integrability of Hirota-Kimura type discretizations[END_REF][START_REF] Celledoni | Geometric properties of Kahan's method[END_REF][START_REF] Celledoni | Integrability properties of Kahan's method[END_REF][START_REF] Celledoni | Detecting and determining preserved measures and integrals of birational maps[END_REF][START_REF] Bogfjellmo | Using aromas to search for preserved measures and integrals in Kahan's method[END_REF] that study the Kahan-Hirota-Kimura discretization [START_REF] Kahan | Unconventional numerical methods for trajectory calculations[END_REF][START_REF] Hirota | Discretization of the Euler top[END_REF][START_REF] Hirota | Discretization of the Lagrange top[END_REF] 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 [START_REF] Blanes | On the necessity of negative coefficients for operator splitting schemes of order higher than two[END_REF]. 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 [START_REF] Butcher | An algebraic theory of integration methods[END_REF][START_REF] Hairer | On the Butcher group and general multi-value methods[END_REF] (see also the textbooks [START_REF] Hairer | Geometric numerical integration[END_REF][START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF] Butcher | B-series: algebraic analysis of numerical methods[END_REF] and the review [START_REF] Mclachlan | Butcher series: a story of rooted trees and numerical methods for evolution equations[END_REF]), the Butcherseries formalism is an important tool in numerical analysis. Originally used for the calculation of order conditions for Runge-Kutta methods, the use of B-series was quickly extended to a variety of applications, in particular in geometric numerical integration [START_REF] Hairer | Geometric numerical integration[END_REF], and more recently in the approximation of stochastic evolutionary problems [START_REF] Burrage | High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations[END_REF][START_REF] Komori | Rooted tree analysis of the order conditions of ROW-type scheme for stochastic differential equations[END_REF][START_REF] Rößler | Stochastic Taylor expansions for the expectation of functionals of diffusion processes[END_REF][START_REF] Debrabant | B-series analysis of stochastic Runge-Kutta methods that use an iterative scheme to compute their internal stage values[END_REF][START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF] or in the theory of rough paths [START_REF] Gubinelli | Ramification of rough paths[END_REF][START_REF] Hairer | Geometric versus non-geometric rough paths[END_REF]. For a large class of integrators, such as Runge-Kutta methods, the modified vector field r f in (1.3) is a B-series in f . An extension of B-series, called aromatic B-series, was introduced in [START_REF] Chartier | Preserving first integrals and volume forms of additively split systems[END_REF][START_REF] Iserles | B-series methods cannot be volumepreserving[END_REF] for the study of volume-preserving integrators. They allow to compute the divergence of a B-series, and were later studied in [START_REF] Mclachlan | B-series methods are exactly the affine equivariant methods[END_REF][START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF][START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF][START_REF] Fløystad | The universal pre-Lie-Rinehart algebras of aromatic trees. In Geometric and harmonic analysis on homogeneous spaces and applications[END_REF] for their algebraic and geometric properties. In [START_REF] Chartier | Preserving first integrals and volume forms of additively split systems[END_REF][START_REF] Iserles | B-series methods cannot be volumepreserving[END_REF], 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 [START_REF] Kang | Volume-preserving algorithms for source-free dynamical systems[END_REF]). However, the space of divergence-free aromatic B-series is infinite-dimensional, the simplest non-trivial example being:

d i,j,k1 pf j k f k j f i f j jk f k f i ¡ f j j f i k f k ¡ f i jk f j f k qf i Kerpdivq.
(1.4)

In [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF], 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 [START_REF] Lee | Introduction to smooth manifolds[END_REF]. For instance, in dimension 3, denote Ω (diff) n pR 3 q the space of smooth differential n-forms on R 3 , then the De Rham complex writes Ω (diff) 0 pR 3 q

Ω (diff) 1 pR 3 q Ω (diff) 2 pR 3 q Ω (diff) 3 pR 3 q, d d d (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 ω Ω 3 pR 3 q is a curl ω dη. 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 [START_REF] Vinogradov | A spectral sequence that is connected with a nonlinear differential equation, and the algebraic-geometry foundations of the Lagrange field theory with constrain[END_REF][START_REF] Tulczyjew | The Euler-Lagrange resolution[END_REF][START_REF] Tsujishita | On variation bicomplexes associated to differential equations[END_REF][START_REF] Vinogradov | The C-spectral sequence, Lagrangian formalism, and conservation laws. I. The linear theory[END_REF][START_REF] Vinogradov | The C-spectral sequence, Lagrangian formalism, and conservation laws. II. The nonlinear theory[END_REF] 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 [START_REF] Anderson | The variational bicomplex[END_REF][START_REF] Olver | Applications of Lie groups to differential equations[END_REF], the introductory article [START_REF] Anderson | Introduction to the variational bicomplex[END_REF], 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 divpf q 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 divpf q 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 [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF][START_REF] Laurent | Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds[END_REF]. 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.

Preliminaries and main results

This section is devoted to the presentation of the new objects and to the associated main algebraic results. We first recall the definition of aromatic forests, and present an extension, called aromatic forms, well-suited for studying the divergence-free combination of forests. Using grafting operations, we define the aromatic bicomplex, a similar tool to the variational bicomplex in the context of differential geometry, and we use its exactness to describe aromatic forms of vanishing divergence. The proofs of the results of this section, as well as further tools, such as the Euler operators, are presented in Section 3, while we give more details and concrete applications in Section 4.

Aromatic forms: definition and operations

B-series were introduced by Hairer and Wanner in [START_REF] Hairer | On the Butcher group and general multi-value methods[END_REF], based on the work of Butcher [START_REF] Butcher | An algebraic theory of integration methods[END_REF]. Their applications in the numerical analysis of deterministic differential equations are numerous (see, for instance, the textbooks [START_REF] Butcher | B-series: algebraic analysis of numerical methods[END_REF][START_REF] Hairer | Geometric numerical integration[END_REF]). The aromatic extension of Butcher-series was introduced independently in the works [START_REF] Chartier | Preserving first integrals and volume forms of additively split systems[END_REF][START_REF] Iserles | B-series methods cannot be volumepreserving[END_REF] to study volume-preserving integrators. In particular, this extension allows us to represent the divergence of standard B-series, which is a key tool in the study of volume-preserving methods (see [START_REF] Hairer | Geometric numerical integration[END_REF]Sect. VI.9]). In the spirit of differential geometry, we work in this paper with an extension of aromatic B-series that is analogous to differential forms, and we follow the graph definition of aromatic B-series of [START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF]. Definition 2.1. Let V be a finite set of nodes, and E V ¢V a set of edges. If a pv, wq E, the edge a is going from v to w, and v is a predecessor of w. We split the set V into vertices V and covertices V ¥ . The covertices are numbered from 1 to p, while the vertices are not numbered.

Each node in V is the target of at most one node. The nodes that are not the target of any node are called roots. The roots are numbered from 1 to n. Any connected component of such a graph either has a root, and is called a tree, or does not have a root, and is called an aroma. We call aromatic forests such graphs, up to equivalence of graphs that preserve the numbering of the covertices and the roots. We write F n,p the set of aromatic forests with n roots and p covertices, F N n,p its subset with forests of exactly N nodes, and F n F n,0 . The number of nodes |γ| is called the order of the forest γ. The elements of F 1 are called aromatic trees, and the subset T of F 1 that contains the trees without aromas is the set of Butcher trees.

We draw the aromatic forests as follows. The vertices are represented as black nodes, and the covertices as circles of the form i , where i is the associated number. The trees are drawn in the ascending order of their roots, from left to right. The aromas are placed in front and their order does not matter. The orientation of the edges goes from top to bottom and in clockwise order for loops. A loop with K nodes is called a K-loop, in the spirit of [START_REF] Iserles | B-series methods cannot be volumepreserving[END_REF]. For instance, the following forest γ F 10 3,2 has two aromas, one 1-loop and one 3-loop, (2.1)

The aromatic forests with N 2 nodes are

F 2 2 t u, F 2 1 t , u, F 2 0 t , , u, F 2 2,1 t 1 , 1 u, F 2 1,1 t 1 , 1 , 1 , 1 u, F 2 0,1 t 1 , 1 , 1 , 1 u, F 2 2,2 t 1 2 , 2 1 u, F 2 1,2 t 1 2 , 2 1 , 1 2 , 2 1 u, F 2 0,2 t 1 2 , 2 1 , 1 2 , 1 2 u.
An aromatic forest in 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 [START_REF] Anderson | Introduction to the variational bicomplex[END_REF] for further details (see also the textbooks [START_REF] Olver | Applications of Lie groups to differential equations[END_REF][START_REF] Kushner | Contact geometry and non-linear differential equations[END_REF][START_REF] Lee | Introduction to smooth manifolds[END_REF]). The infinite jet bundle J V pR d q is used to describe coordinate-free Taylor expansions. It is the infinite-dimensional vector bundle over R d where elements of J V

x pR d q have coordinates pf i pxq, f i j pxq, f i jk pxq, . . . q, f i : R d Ñ R, f i j 1 ,...,jm ff i fx j 1 . . . fx jm . (2.2)
We consider differential forms on J V pR d q and define the following contact forms

θ i j 1 ,...,jp df i j 1 ,...,jp ¡ f i j 1 ,...,jp,k dx k .
The set tdx i , θ i , θ i j , . . . u locally generates all differential forms on J V pR d q. A differential form of type pn, pq is a sum of terms of the form

gpf qdx i 1 ¤ ¤ ¤ dx in θ j 1 I 1 ¤ ¤ ¤ θ jp Ip ,
where gpf q is a functional of finitely many coordinates (2.2). The differential forms of type pn, pq are gathered in the set Ω

(diff)
n,p pJ V pR d qq. In agreement with the numerical analysis literature, the elementary differential map is defined with the help of the vector basis f i . Definition 2.2. Let γ F n,p , f : R d Ñ R d a smooth vector field, and R tr 1 , . . . , r n u V the n roots of γ, then the elementary differential F pγqpfq is the following tensor:

F pγqpfq iwt1,...,du wV zR ¹ vV f iv I Πpvq f ir 1 ¤ ¤ ¤ f ir n θ i 1 I Πp 1 q ¤ ¤ ¤ θ i p I Πp p q ,
where Πpvq is the set of predecessors of v V . We extend F on SpanpF n,p q by linearity.

For instance, we find

F p qpfq divpf q d i1 f i i , F p qpfq d i,j1 f i j f j f i , F p 1 qpfq d i1 f i θ i .
The forest γ in (2.1) represents the elementary differential

F pγqpfq ¡ d j,k,l1 f k j f l k f j l ©¡ d j,k1 f j jk f k © d ir 1 ,ir 2 ,ir 3 ,j,k1 f ir 1 j f k f ir 2 f ir 1 f ir 2 f ir 3 θ ir 3 θ j k .
In the following, we work with specific linear combinations of aromatic forests in SpanpF n,p q.

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 γ F n,p , let S n (resp. S ¥ p ) be the set of permutations of the roots of γ (resp. the covertices of γ). We define the roots wedge of γ as

γ 1 n! σS n εpσqσγ,
where εpσq is the signature of the permutation σ. Similarly, the covertices wedge is

¥ γ 1 p! σS ¥ p εpσqσγ,
and the total wedge is ¥ ¥ .

We extend the wedge operations to SpanpF n,p q by linearity and we denote the set of aromatic forms Ω n,p SpanpF n,p q, respectively Ω N n,p SpanpF N n,p q, and Ω n Ω n,0 . As 2 , the operator :

SpanpF n,p q Ñ Ω n,p is a projection on Ω n,p . Example. Let γ F 2 , then γ 1 2 p ¡ q Ω 2 .
For γ given by (2.1), the associated aromatic form is

γ 1 12 ¡ 2 1 1 2 1 2 ¡ 2 1 ¡ 2 1 ¡ 1 2 ¡ 1 2 ¡ 2 1 ¡ 2 1 1 2 1 2 2 1 ©
.

We now define the grafting and replacing operations on aromatic forests. Example. Let γ 1 and r its root, then

D r γ 1 1 , γ rÑ 2 2 1 , γ 1 Ñ .
The horizontal derivative is defined using grafting operations, while the vertical derivative uses the replacing operation. Definition 2.5. Let γ F n,p , the horizontal and vertical derivatives are

d H γ D rn γ, d V γ v V γ vÑ p 1 .
We

d H γ 1 , d H γ 2 1 2 ¡ ¡ ¡ © , d H γ 3 1 2 ¡ 1 1 ¡ 1 ¡ 1 © , d V γ 1 1 , d V γ 2 1 1 1 , d V γ 3 2 1 1 2 p 2 1 ¡ 1 2 q. A calculation yields d 2 H γ 2 0, so that d H γ 2 Ψ 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 2 H 0, d 2 V 0, 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 Ω 1 SpanpF 1 q respectively in [START_REF] Chartier | Preserving first integrals and volume forms of additively split systems[END_REF][START_REF] Iserles | B-series methods cannot be volumepreserving[END_REF] for d H and in [START_REF] Fløystad | The universal pre-Lie-Rinehart algebras of aromatic trees. In Geometric and harmonic analysis on homogeneous spaces and applications[END_REF] for d V . In this last work, the trace operator Tr : Ω 1,1 Ñ Ω 0 on aromatic trees is studied. For γ Ω 1,1 , it is given with our notations by

Tr γ pD rÑ 1 γq 1 Ñ , and it makes the following diagram commute.

Ω 1,1 Ω 1 Ω 0 Tr d H d V

The aromatic bicomplex: exactness and description of solenoidal forms

The variational bicomplex is a powerful tool of variational calculus [START_REF] Anderson | Introduction to the variational bicomplex[END_REF]. 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 Impd H | Ω 1,p q, 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 . [START_REF] Anderson | Introduction to the variational bicomplex[END_REF] in the following way. Let σ be the standard duality map between the type n-alternating contravariant tensors and the type pd¡nq-differential forms with respect to the standard Euclidean coordinates and volume form dx 1 . . . dx d (see, for instance, [START_REF] Lee | Introduction to smooth manifolds[END_REF]). Let f : R d Ñ R d , the map σF p.qpfq sends an aromatic form in Ω n,p to a differential form in

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

Ω (diff) d¡n,p pJ V pR d qq. The derivatives d H , d V on the aromatic bicomplex and d (diff) H , d (diff) V on the variational bicomplex satisfy d (diff) H σF pγqpfq σF pd H γqpf q, d (diff) V F pγqpfq p¡1q n p F pd V γqpf q. . . . . . . . . . . . . Ω 2,2 Ω 1,2 Ω 0,2 . . . Ω 2,1 Ω 1,1 Ω 0,1 . . . Ω 2 Ω 1 Ω 0 0 0 0 d H d H d V d H d V d V d H d H d V d H d V d V d H d H d V d H d V d V 0 0 0 Ω N N,N . . . Ω N 0,N . . . . . . 0 Ω N N . . . Ω N 0 0 0 d H d H d V d V d H d V d H d V Figure 1:
The aromatic bicomplex (left) and its subcomplex of order N (right).

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

σF p.qpfq. 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 ¥ 0,

Impd H | Ω N n 1,p q Kerpd H | Ω N n,p q, Impd V | Ω N n,p q Kerpd V | Ω N n,p 1 q.
With Theorem 2.9, it is straightforward to generate all the solenoidal aromatic forms by considering the d H γ for γ F 2 . For example, the only basis element of Ψ 3 (which corresponds to the vector field (1.4)) is

2d H ¡ ¡ , (2.3)
and a basis of Ψ 4 is given by the forests

2d H ¡ ¡ ¡ , 2d H 2 ¡ 2 ¡ ¡ , 2d H ¡ ¡ ¡ .
However, the set td H γ, γ F 2 u does not form a basis of the solenoidal forms in general. We give a basis of Ψ 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 divpf q 0. With aromatic forms, it amounts to sending all forests containing an aroma with a 1-loop to 0. We consider the vector space A spanned by aromatic forests containing at least one 1-loop, A n,p A Ω n,p , and A N n,p A Ω N n,p . We write r F n,p the set of aromatic forests in F n,p without 1-loops, r

Ω n,p Ω n,p {A n,p , r Ω N n,p Ω N n,p {A N n,p , r Ω n r Ω n,0 , r Ψ Kerpd H | r Ω 1 q, and r Ψ N Kerpd H | r Ω N
1 q. 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 . . . . . . . . .

. . . r Ω N 2,2 r Ω N 1,2 r Ω N 0,2 . . . r Ω N 2,1 r Ω N 1,1 r Ω N 0,1 . . . r Ω N 2 r Ω N 1 r Ω N 0 0 0 0 d H d H d V d H d V d V d H d H d V d H d V d V d H d H d V d H d V d V 0 0 0 r Ω 1 1,1 Spanp 1 q r Ω 1 0,1 0 0 r Ω 1 1 Spanp q r Ω 1 0 0 0 0 d H d H d V d V Figure 2:
The divergence-free aromatic bicomplex of order N (left) and of order N 1 (right).

is not exact. Indeed, we have Kerpd H q, but Impd H q. 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 $ 1.

The main consequence of the exactness of the aromatic bicomplex in both contexts is that the assumption divpf q 0 does not create new non-trivial solenoidal forms. Theorem 2.11. For N $ 1, and all n ¥ 1, p ¥ 0, the kernel of the divergence operator d H satisfies

Kerpd H | r Ω N n,p q Kerpd H | Ω N n,p q{A N n,p
, 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 and 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

Impd H | r Ω N n 1,p q Impd H | Ω N n 1,p q{A N
n,p . Theorem 2.9 and Theorem 2.10 yield the desired identity.

A generating set of the solenoidal forms r Ψ N of order N ¡ 1 is obtained by deleting the 1loops in the generating set td H γ, γ F N 2 u of Ψ 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 r Ψ 3 Spanp ¡ q.

We give generators of the solenoidal forms r Ψ 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 Ω N n,p and r Ω N n,p in the first two rows of the bicomplex and deduce the dimension of Ψ N in Theorem 4.2 and of r Ψ N in Theorem 4.4. A surprising fundamental result is that the solenoidal forms in Ψ N are enumerated by the difference between the number of aromatic trees in Ω N 1 and the number of aromas with 1-loops ΩN 0 (see Table 1

|Ψ N | |Ω N 1 | ¡ | ΩN 0 |.

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.

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.

D rÑu γ is R 0 if r R 0 and R 0 ztru else. We define D r γ uV ztvu D rÑu γ.
Let q be a non-negative integer, u V with u $ v if γ F n,p , and I R 0 . We define D I and D IÑu by

D I γ φ: IÑV ¹ wI D wÑφpwq γ, D IÑu γ ¹ wI D wÑu γ,
and D q and D qÑu by

D q γ IR 0 ,|I|q D I γ, D qÑu γ IR 0 ,|I|q D IÑu γ.
Let γ F n,p and v V , then γ v F n |Πpvq|,p is the graph obtained by cutting off the edges of γ pointing to v and placing the newly obtained roots after the roots in R (in an arbitrary order). In addition, we add a symbol on the node v and we fix R 0 Πpvq as the set of detached nodes.

The tools from Definition 3.1 are extended by linearity. In simple words, the detaches the predecessors of a given node v and put them in a set R 0 . As long as the symbol is present, the grafting operators D I and D q graft the detached nodes on every node except v . In particular, the following operation is trivial:

D |Πpvq|Ñv pγ v q γ.
Example. Consider the forest γ 1 2 that we label for the sake of clarity. The marked forest γ 1 has the set of detached nodes R 0 t1u and satisfies

γ 1 2 1 , pDγ 1 q 2 1 , Dpγ 1 q 2 1 1 2 .
For simplicity, we write γ vÑ k for pγ vÑ k q k , and D D 1 in the rest of the paper. Note that, in general, D 2 In the context of differential geometry, the Euler operator describes the differential forms that are divergences [START_REF] Olver | Applications of Lie groups to differential equations[END_REF]. It has a variety of applications such as the study of conservation laws and Lagrangians for PDEs [START_REF] Olver | A resolution of the Euler operator[END_REF][START_REF] Shakiban | A resolution of the Euler operator[END_REF][START_REF] Hereman | Continuous and discrete homotopy operators and the computation of conservation laws[END_REF][START_REF] Hereman | Continuous and discrete homotopy operators: a theoretical approach made concrete[END_REF][START_REF] Poole | The homotopy operator method for symbolic integration by parts and inversion of divergences with applications[END_REF]. We define a similar operator for aromatic forms.

Definition 3.3. For γ F n,p and v V , the Euler operators E v : F n,p Ñ Ω n,p and E ¥ v : F 0 Ñ Ω 0,1
are given by E v γ p¡1q |Πpvq| pD Πpvq γ v q , E ¥ v γ p¡1q |Πpvq| pD Πpvq γ vÑ 1 q and are extended by linearity on Ω n,p . The Euler operator E : Ω n,p Ñ Ω n,p and its variant E ¥ : Ω 0 Ñ Ω 0,1 are given by

Eγ vV E v γ, E ¥ γ vV E ¥ v γ.
The output of E ¥ on the forms of first orders is

E ¥ 0, E ¥ 2 1 , E ¥ ¡2 1 , E ¥ ¡2 1 .
From these examples, we observe that the Euler operator vanishes on aromatic forms that are divergences:

E ¥ d H 0, E ¥ d H 0, E ¥ d H 0. Proposition 3.4. Let γ F n,p and v V , the Euler operators satisfy E v d H γ 0, Ed H γ 0 and for p 0, E ¥ satisfies E ¥ d H γ 0.
Proof. Let γ F n,p and r n be its last root, then

E v d H γ uV E v D rnÑu γ p¡1q |Πpvq| 1 pD Πpvq D rn γ v q u$v p¡1q |Πpvq| pD Πpvq D rnÑu γ v q 0,
where the two terms correspond to the cases u v and u $ v, and where we used the convention of Remark 3.2. The identities with E and E ¥ are direct consequences.

We know that the composition of the two maps E ¥ and d H vanishes. One is then interested in a necessary and sufficient condition for an aromatic form γ Ω 0 to be a divergence. The following chain is called the variational complex.

Ω 1 Ω 0 Ω 0,1 d H E ¥ (3.2)
The complex (3.2) is exact if Impd H q KerpE ¥ q; that is, if γ Ω 0 is a divergence if and only if E ¥ γ 0. A fundamental question in variational calculus [START_REF] Olver | Applications of Lie groups to differential equations[END_REF] 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 γ F n,p , the Euler operator E q v γ of order q ¥ 0 on v V is

E q v γ p¡1q |Πpvq|¡q pD |Πpvq|¡q γ v q ,
that we extend on Ω n,p by linearity. The higher Euler operators are

E q γ vV E q v γ.
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

E q v γ 0 if q ¡ |γ|. Note also that E 0 v E v .
Proposition 3.6. The Euler operators satisfy for all γ F n,p and all vertices v V ,

γ V q0 D q E q v γ. (3.3)
In particular, we have

|γ| γ V q0 D q E q γ. (3.4)
In addition, if γ satisfies E p γ 0 for p 0, . . . , n ¡ 1, then γ D n ph pnq γq, where

h pnq γ 1 |γ| V pn ¢ p n ¡1 D p¡n E p γ.
The proof shares similarities with the approach of [START_REF] Hereman | Continuous and discrete homotopy operators: a theoretical approach made concrete[END_REF] and uses the following intermediate result in the spirit of the Leibniz rule. Lemma 3.7. Let γ F n,p and v V , let three integers p, q, k such that p q k |Πpvq| ¡ 1.

Then, the following holds:

pq 1qD q 1 D kÑv pD p γ v q pk 1qD q D k 1Ñv pD p γ v q pp 1qD q D kÑv pD p 1 γ v q .
More precisely, we have

D q D kÑv pD p γ v q k ņ0 p¡1q k¡n ¢ q n q ¢ p k ¡ n p D q n pD p k¡n γ v q . (3.5)
Proof of Proposition 3.6. Let γ F n,p . With the help of Lemma 3.7, we get

E v γ p¡1q |Πpvq| pD |Πpvq| γ v q pγ v q 1 |Πpvq|0 ¡ p¡1q |Πpvq|¡1 pD |Πpvq| γ v q 1 |Πpvq|¥1 γ1 |Πpvq|0 ¡ p¡1q |Πpvq|¡1 |Πpvq| DpD |Πpvq|¡1 γ v q p¡1q |Πpvq|¡1 |Πpvq| D 1Ñv pD |Πpvq|¡1 γ v q
We apply the same reasoning to the last term

p¡1q |Πpvq|¡1 |Πpvq| D 1Ñv pD |Πpvq|¡1 γ v q γ1 |Πpvq|1 ¡ p¡1q |Πpvq|¡2 |Πpvq| p|Πpvq| ¡ 1q DD 1Ñv pD |Πpvq|¡2 γ v q 2p¡1q |Πpvq|¡2 |Πpvq| p|Πpvq| ¡ 1q D 2Ñv pD |Πpvq|¡2 γ v q .
Iterating this reasoning, we find

E v γ γ ¡ V p1 p¡1q |Πpvq|¡p p |Πpvq| p ¨DD p¡1Ñv pD |Πpvq|¡p γ v q . (3.6)
Using the formula (3.5) in (3.6), standard combinatorics [START_REF] Merris | Combinatorics. Wiley-Interscience Series in Discrete Mathematics and Optimization[END_REF] yield

γ E v γ V p1 p¡1q |Πpvq|¡p p |Πpvq| p ¨p¡1 ņ0 p¡1q p¡n¡1 pn 1q ¢ |Πpvq| ¡ n ¡ 1 |Πpvq| ¡ p D n 1 pD |Πpvq|¡n¡1 γ v q E v γ V ņ0 V pn 1 n 1 p |Πpvq| p ¨¢|Πpvq| ¡ n ¡ 1 |Πpvq| ¡ p D n 1 E n 1 v γ E v γ |Πpvq| ņ1 |Πpvq| pn np|Πpvq| ¡ nq!pp ¡ 1q! |Πpvq|!pp ¡ nq! D n E n v γ E v γ |Πpvq| ņ1 ¢ |Πpvq| n ¡1 |Πpvq|¡n p0 ¢ p n ¡ 1 p D n E n v γ V ņ0 D n E n v γ.
We sum (3.3) on all the nodes v 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.

D q 1 D kÑv pD p γ v q 1 q 1 DD q D kÑv pD p γ v q 1 q 1 Ş1 S 2 S 3 tuuΠpvq |S 1 |q,|S 2 |k,|S 3 |p D u D S 1 D S 2 Ñv pD S 3 γ v q 1 q 1 Ş1 S 2 S 3 tuuΠpvq |S 1 |q,|S 2 |k,|S 3 |p D S 1 D S 2 uÑv pD S 3 γ v q D S 2 Ñv pD S 3 u γ v q & k 1 q 1 Ş1 S 2 S 3 Πpvq |S 1 |q,|S 2 |k 1,|S 3 |p D S 1 D S 2 Ñv pD S 3 γ v q p 1 q 1 Ş1 S 2 S 3 Πpvq |S 1 |q,|S 2 |k,|S 3 |p 1 D S 1 D S 2 Ñv pD S 3 γ v q k 1 q 1 D q D k 1Ñv pD p γ v q p 1 q 1 D q D kÑv pD p 1 γ v q .
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 γ F 0,1 , let the variational homotopy operator h V be

h V γ 1
|γ| γ 1 Ñ , and for γ F 0 , let the horizontal homotopy operator h H be

h H γ 1 |γ| V q1 1 q D q¡1 E q γ.
We extend the definition on Ω 0,1 and Ω 0 by linearity. These operators satisfy for all γ Ω 0 , pd H h H h V E ¥ qγ γ.

(3.7)

In particular, the variational complex (3.2) is exact.

Proof. Let γ Ω 0 . From Proposition 3.6, we obtain

γ 1 |γ| V q0 D q E q γ 1 |γ| Eγ 1 |γ| V q1 1 q DD q¡1 E q γ h V E ¥ γ d H h H γ.
Proposition 3.4 yields that Impd H q KerpEq. If γ KerpEq, then the identity (3.7) becomes γ d H ph H γq Impd H q. The exactness of the chain follows straightforwardly.

We present an alternative horizontal homotopy operator on Ω 0 in Subsection 4.4, with some examples of the outputs of both operators in Table 3.

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 γ F n,p , and pr 1 , . . . , r n q its roots. The horizontal derivative satisfies

d 2 H γ 1 n! ¥ σS n εpσqD r n¡1 D rn σγ 1 n! ¥ σS n εppnpn ¡ 1qqσqD r n¡1 D rn pnpn ¡ 1qqσγ ¡ 1 n! ¥ σS n εpσqD rn D r n¡1 σγ
¡d 2 H γ, where pnpn ¡ 1qq S n is a transposition and where we used that D r n¡1 and D rn commute. For the vertical derivative, a similar approach on γ F n,p¡1 gives

d 2 V γ v,wV v$w pγ vÑ p q wÑ p 1 1 p!pp 1q! v,wV v$w σS ¥ p ŗ σS ¥ p 1
εpr σσqr σpσpγ vÑ p q wÑ p 1 q

1 p!pp 1q! v,wV v$w σS ¥ p ŗ σS ¥ p 1
εpr σσqr σσpγ vÑ p ,wÑ p 1 q

1 p!pp 1q! v,wV v$w σS ¥ p ŗ σS ¥ p 1 εpr σσp p p 1 qqr σσpγ vÑ p 1 ,wÑ p q ¡ 1 p!pp 1q! v,wV v$w σS ¥ p ŗ σS ¥ p 1
εpr σσqr σpσpγ wÑ p q vÑ p 1 q ¡d 2 V γ. As Ω n,p pF n,p q, we get the desired identities. The commutativity of d H and d V is straightforward.

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 : Ω n,p Ñ Ω n,p¡1 given by

h V γ p |γ| γ p Ñ , satisfies for p ¥ 1 and γ Ω n,p , pd V h V h V d V qγ γ.
Proof. Let γ F n,p . On the first hand, we have

d V h V γ p |γ| γ p |γ| v V γ p Ñ,vÑ p p |γ| γ 1 |γ| pp ¡ 1q! v V σS ¥ p εpσqσpγ p Ñ,vÑ p q.
On the other hand, we have

h V d V γ 1 |γ| v V h V pγ vÑ p 1 q p 1 |γ| pp 1q! v V σS ¥ p 1 εpσqσpγ vÑ p 1 q p 1 Ñ 1 |γ| p! v V σS ¥ p 1 σp p 1 q p 1 εpσqσpγ vÑ p 1 q p 1 Ñ 1 |γ| p! v V σS ¥ p 1 σp p 1 q$ p 1 εpσqσpγ vÑ p 1 q p 1 Ñ 1 |γ| p! v V σS ¥ p εpσqσγ 1 |γ| p! p i1 v V σS ¥ p σp i q p εpσqσpγ i Ñ,vÑ p q |γ| ¡ p |γ| γ 1 |γ| p! p i1 v V σS ¥ p εpσp i p qqσpγ p Ñ,vÑ p q |γ| ¡ p |γ| γ ¡ 1 |γ| pp ¡ 1q! v V σS ¥ p εpσqσpγ p Ñ,vÑ p q.
We deduce that pd V h V h V d V qγ γ, and we obtain the desired equality by linearity as Ω n,p SpanpF n,p q and 2 .

The expression of the horizontal homotopy operator is much more technical. We refer the reader to [START_REF] Anderson | On the existence of global variational principles[END_REF][START_REF] Tulczyjew | The Euler-Lagrange resolution[END_REF][START_REF] Olver | Applications of Lie groups to differential equations[END_REF] 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 : Ω n,p Ñ Ω n 1,p by h H γ 1 |γ| V q0 n 1 q n 1 D q E q 1 γ. It satisfies for n ¥ 1 and γ Ω n,p , pd H h H h H d H qγ γ.
Proof. For the sake of simplicity, we consider γ F n . We define for v V

J v γ V q0 ΠpvqtruIJ |I|q n 1 q n 1 D I E I,r v γ,
where we denote E I,r v γ p¡1q |Πpvq|¡pq 1q pD J γ v q and where r becomes the new root of J v γ.

On the first hand, we have

d H γ d H σSn εpσq n! σγ σSn εpσq n! D r σ ¡1 pnq σγ,
and J v d H γ is given by

J v d H γ V q0 σSn u$v ΠpvqtruIJ
D I E I v γ if r r, u v, °wV D rÑw D Iztru E Iztru,r v γ if r I, u v, ¡ °w$v D rÑw D I E I,r v γ if r J, u v, D rÑu D I E I,r v γ if u $ v. Thus, for γ F n,p , J v d H γ is given by J v d H γ V q0 σSn ΠpvqIJ |I|q εpσq pq nqpn ¡ 1q! D I E I v σγ V q1 σSn wV ΠpvqtruIJ |I|q¡1 εpσq pq nqpn ¡ 1q! D r σ ¡1 pnq Ñw D I E I,r v σγ ¡ V q0 σSn w$v ΠpvqtruIJ |I|q εpσq pq nqpn ¡ 1q! D r σ ¡1 pnq Ñw D I E I,r v σγ V q0 σSn u$v ΠpvqtruIJ |I|q εpσq pq nqpn ¡ 1q! D r σ ¡1 pnq Ñu D I E I,r v σγ V q0 n q n D q E q v γ V q0 σSn ΠpvqtruIJ |I|q εpσq pq n 1qpn ¡ 1q! D r σ ¡1 pnq D I E I,r v σγ.
As D q E q v γ is unchanged by the application of the wedge operator, we find

J v d H γ V q0 n q n D q E q v γ V q0 νSn σSn ΠpvqtruIJ |I|q εpσνq pq n 1qpn ¡ 1q!n! νD r σ ¡1 pnq D I E I,r v σγ.
On the other hand, we have

d H J v γ σSn εpσq n! d H J v σγ d H V q0
pq n 1qpn ¡ 1q!pn!q 2 D r pσνq ¡1 pnq ηpnpn 1qqνD I E I,r v σγ V q0 q 1 q n 1 D q 1 E q 1 v γ ¡ V q0 σSn η Sn ΠpvqtruIJ |I|q εpσηq pq n 1qpn ¡ 1q!n! ηD r σ ¡1 pnq D I E I,r v σγ V q1 q q n D q E q v γ ¡ V q0 σSn η Sn ΠpvqtruIJ |I|q εpσηq pq n 1qpn ¡ 1q!n! ηD r σ ¡1 pnq D I E I,r v σγ,
where we split the cases σνpn 1q n 1 and σνpn 1q $ n 1, and where we substitute ν by νpnpn 1qqη in the latter case. The last equality is obtained by substituting σ with νσ. Using Proposition 3.6, we get

pJ v d H d H J v q γ V q0 D q E q v γ γ. (3.8)
Summing on v V gives the desired homotopy identity.

Remark 3.11. The identity (3.8) suggests simpler homotopy operators, obtained by fixing a node v and considering the operator J v . We emphasize that this approach does not work. Indeed, given an aromatic form γ Ω n,p , there is no canonical choice of the node v, so that J v γ is ill-defined, where h V γ is well-defined. In the proof of Proposition 3.10, J v γ makes sense since we consider a single aromatic forest γ F n,p .

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

r h H γ 1 |γ| vV V q0 n 1 q n 1 vR D q E q 1 v γ.
For n ¡ 1 and γ r Ω n,p , the horizontal homotopy identity is

pd H r h H r h H d H qγ γ, while for n 1, the identity is pd H r h H r h H d H qγ γ ¡ Rγ.
(3.9)

The remainder in (3.9) is the linear map Rγ 1 |γ| E r γ for γ r F 1,p , with r the root of γ.

Moreover, if γ r Ω N 1,p satisfies d H γ 0, then Rγ 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 γ r F n . For v V , we define

r J v pγq V q0 ΠpvqtruIJ |I|q n 1 q n 1 vR D I E I,r v γ.
We have if r v,

D I E I,r v D rÑu γ 4 0 if u v, D rÑu D I E I,r v γ if u $ v, and if r $ v, D I E I,r v D rÑu γ
D I E I v γ if r r, u v, °wV D rÑw D Iztru E Iztru,r v γ if r I, u v, ¡ °w$v D rÑw D I E I,r v γ if r J, u v, D rÑu D I E I,r v γ if u $ v. Thus, r J v d H γ is given by r J v d H γ 1 v$r σ ¡1 pnq V q0 σSn ΠpvqIJ |I|q εpσq pq n ¡ 1 vR qpn ¡ 1q! D I E I v σγ V q1 σSn wV ΠpvqtruIJ |I|q¡1 εpσq pq n ¡ 1 vR qpn ¡ 1q! D r σ ¡1 pnq Ñw D I E I,r v σγ ¡ V q0 σSn w$v ΠpvqtruIJ |I|q εpσq pq n ¡ 1 vR qpn ¡ 1q! D r σ ¡1 pnq Ñw D I E I,r v σγ % V q0 σSn u$v ΠpvqtruIJ |I|q εpσq pq n ¡ 1 vR qpn ¡ 1q! D r σ ¡1 pnq Ñu D I E I,r v σγ 1 v$r σ ¡1 pnq V q0 σSn εpσq pq n ¡ 1 vR qpn ¡ 1q! D q E q v σγ V q0 σSn ΠpvqtruIJ |I|q εpσq pq n 1 vR qpn ¡ 1q! D r σ ¡1 pnq D I E I,r v σγ
Then, we find for n 1 and v R,

r J v d H γ V q0 σSn νSn ΠpvqtruIJ |I|q εpσνq pq n 1 vR qpn ¡ 1q!n! νD r σ ¡1 pnq D I E I,r v σγ,
and otherwise

r J v d H γ V q0 n ¡ 1 vR q n ¡ 1 vR D q E q v γ V q0 σSn νSn ΠpvqtruIJ |I|q εpσνq pq n 1 vR qpn ¡ 1q!n! νD r σ ¡1 pnq D I E I,r v σγ,

20

where for v R, the wedge adds the missing term v r σ ¡1 pnq , and rescales the expression with a coefficient n¡1 n . On the other hand, we have

d H r J v γ d H V q0 σSn νS n 1 ΠpvqtruIJ |I|q εpσνq pq n 1 vR qpn!q 2 νD I E I,r v σγ V q0 σSn νS n 1 ΠpvqtruIJ |I|q εpσνq pq n 1 vR qpn!q 2 D r pσνq ¡1 pn 1q νD I E I,r v σγ V q0 σSn ΠpvqtruIJ |I|q εpσq pq n 1 vR qn! D I,r E I,r v σγ ¡ V q0 σSn νSn η Sn ηpnqn ΠpvqtruIJ |I|q εpσνηq pq n 1 vR qpn ¡ 1q!pn!q 2 ¤ D r pσνq ¡1 pnq ηpnpn 1qqνD I E I,r v σγ V q0 q 1 q n 1 vR D q 1 E q 1 v γ ¡ V q0 σSn η Sn ΠpvqtruIJ |I|q εpσηq pq n 1 vR qpn ¡ 1q!n! ηD r σ ¡1 pnq D I E I,r v σγ V q1 q q n ¡ 1 vR D q E q v γ ¡ V q0 σSn η Sn ΠpvqtruIJ |I|q εpσηq pq n 1 vR qpn ¡ 1q!n! ηD r σ ¡1 pnq D I E I,r v σγ.
Using Proposition 3.6, we get for n ¡ 1, p r

J v d H d H r J v q γ V q0 D q E q v γ γ, and for n 1, p r J v d H d H r J v q γ V q0 D q E q v γ ¡ E r γ1 vR γ ¡ E r γ1 vR .
Summing on v V gives the desired homotopy identities. Following the proof of Proposition 3.4, we observe that E r d H 0 on Ω 2,p . Thus, we deduce from the identity (3.9) that γ Impd H q if and only if d H γ 0 and E r γ 0 in the case n 1. Assume that d H γ 0 for γ Ω N 1,p , and apply E r to the homotopy identity. As

E 2 r E r , it yields E r γ 1 N E r γ.
We deduce that E r γ 0 or N 1. Thus, the bicomplex of order N is exact if and only if N $ 1.

Remark 3.13. Let γ r

Ω N 1 with N ¡ 1, define the modified homotopy operators by

r h 1 H r h H ¡ 1 1 N ¡ 1 E © , r h 2 H r h H ¡ 1 1 N ¡ 1 E r © .
Then, the homotopy identity (3.9) on r Ω N

1 with N ¡ 1 is replaced by the simpler identity pd H r h 2 H r h 1 H d H qγ γ. (3.10)
The proof of (3.10) relies on the identities E 2 r E r and Ed H d H E r on r Ω 1 .

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.

The augmented aromatic bicomplex

Following [START_REF] Anderson | The variational bicomplex[END_REF][START_REF] Anderson | Introduction to the variational bicomplex[END_REF] 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 ¥ 1, define the interior Euler operator I : Ω 0,p Ñ Ω 0,p by Iγ E p γ p¡1q |Πp p q| pD |Πp p q| γ p q , write I p IpΩ 0,p q, I N p IpΩ N 0,p q, and the variational derivative δ V I ¥ 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 δ V E ¥ on Ω 0 .

. . .

Ω 2 Ω 1 Ω 0 I 1 I 2 . . . d H d H d H δ V δ V δ V (4.1)
Theorem 4.1. Define the augmented homotopy operators as

h H γ V q1 1 q D q¡1 E q p γ, h V γ I ¥ h V .
Then, the maps I and δ V satisfy

I 2 I, Id H 0, δ 2 V 0,
and the following identities hold

γ pI d H h H qγ, γ Ω 0,p , p ¥ 1, . . . . . . . . . . . . . . . Ω 2,2 Ω 1,2 Ω 0,2 I 2 0 . . . Ω 2,1 Ω 1,1 Ω 0,1 I 1 0 . . . Ω 2 Ω 1 Ω 0 0 0 0 d H d H d V d H d V d V I δ V d H d H d V d H d V d V I δ V d H d H d V d H d V d V δ V
γ pδ V h V h V δ V qγ, γ I 1 , γ pδ V h V h V δ V qγ, γ I p , p ¡ 1.
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 [START_REF] Anderson | The variational bicomplex[END_REF]Chap. 4] for the proof of Theorem 4.1.

Proof. We first observe that for all γ F 1,p ,

Id H γ v $ p
p¡1q |Πp p q| pD |Πp p q| D r 1 Ñv γ p q p¡1q |Πp p q| 1 pD Πp p q D r 1 γ p q 0.

Using Equation (3.3) with v p gives the augmented horizontal homotopy identity pI d H h H q γ γ.

Applying I to this last equality yields I ¥I I. Let us now look at the Euler-Lagrange complex. Let γ Ω 0,p , then d V γ Ω 0,p 1 . We apply the augmented horizontal homotopy identity to

d V γ, pId V d H h H d V qγ d V γ.
We apply d V and use that d V and d H commute:

pd V Id V d H d V h H d V qγ 0.
Since Id H 0, applying I yields

δ 2 V γ 0. Let γ I p , the augmented horizontal homotopy applied to d V γ gives d V γ δ V γ d H h H d V γ. (4.2)
If p 1, then we use the identity (4.2) in the vertical homotopy identity to get

γ d V h V γ h V δ V γ d H r γ,
where we used that h V and d H commute and where

r γ Ω 1,1 . Applying I yields γ pδ V h V h V δ V qγ. If p ¡ 1, we apply the augmented horizontal homotopy identity to h V γ, h V γ h V γ d H h H h V γ.
Applying the vertical derivative d V gives

d V h V γ d V h V γ d H d V h H h V γ, (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

γ d V h V γ h V δ V γ d H r γ,
where r

γ Ω 1,p . As I ¥ d H 0 and Iγ γ, we find

γ pδ V h V h V δ V qγ.
The exactness of the augmented aromatic bicomplex is a straighforward consequence of Theorem 2.9 and of the augmented homotopy identities.

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 T N of rooted trees of order N (sequence A000081 in the OEIS [START_REF]The On-Line Encyclopedia of Integer Sequences[END_REF]) has generating function tpzq V Ņ1 T N z n z z 2 2z 3 4z 4 9z 5 20z 6 . . . and satisfies the functional equation

tpzq z exp £ V ķ1 1 k tpz k q . ( 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 t1, . . . , N ¡ 1u Ñ t0, . . . , N ¡ 1u modulo the symmetric group S N ¡1 (node 0 is the root, and we "forget the labels").

The number |Ω N 0 | of scalar aromas in addition to the empty aroma (sequence A001372) has generating function

apzq : 1 V Ņ1 |Ω N 0 |z n 1 z 3z 2 7z 3 19z 4 47z 5 130z 6 . . .
and is related to tpzq by the equation

apzq V ¹ k1 ¡ 1 ¡ tpz k q © ¡1
.

The scalar aromas are mapping patterns of functions t1, . . . , N u Ñ t1, . . . , N u mod S N .

We introduce two new spaces, the aromatic forms with a 1-loop, called the self-looped aromatic forms Ωn,k , and their complements, the non-self-looped aromatic forms Ω n,k , and generating functions åpzq of Ω0 and apzq of Ω 0 . As for mapping patterns on t1, . . . , N u, the mappings in Ω0 (enumerated by sequence A217896) have at least one fixed point and those in Ω 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 zapzq tpzqapzq.

Therefore åpzq apzq ¡ 1 ¡ apzq apzqptpzq ¡ zq tpzq ¡ 1 z 2z 2 5z 3 13z 4 34z 5 90z 6 . . . . (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 pu, zq V ķ0 V Ņ1 § § Ω N k,p § § u k z N
be the bivariate generating function for row p of the aromatic bicomplex, let

c p pzq V Ņ1 § § I N p § § z N
be the generating function of the type-p functional forms, and let

spzq V Ņ1 § § Ψ N § § z N
be the generating function of the solenoidal aromatic trees. Then b 0 pu, zq apzq exp

£ V ķ1 p¡1q k¡1 k u k tpz k q , b 1 pu, zq b 0 pu, zq tpzqp1 u ¡ utpzqq p1 ¡ tpzqq 2 ,
c 1 pzq zb 1 p0, zq zapzqtpzq p1 ¡ tpzqq 2 , spzq apzqtpzq ¡åpzq.

Proof. We first note that the aromatic trees Ω 1 are given by the product of a scalar aroma (enumerated by apzq) with a rooted tree (enumerated by tpzq). This is the Cartesian product construction of enumerative combinatorics. Therefore the aromatic trees are enumerated1 by apzqtpzq.

An element of Ω N k is given by the product of a scalar (enumerated by apzq) with a wedge product of k rooted trees, enumerated by tpzq. Since the rooted trees are unordered and distinct, this is the power set construction; the expression for b 0 pu, zq then follows from [START_REF] Flajolet | Analytic combinatorics[END_REF]Proposition III.5].

The bottom row of the bicomplex is exact (Theorem 2.9), so the dimension of Ψ is given by the alternating sum

spzq V Ņ1 |Ω N 2 | ¡ |Ω N 3 | |Ω N 4 | ¡ |Ω N 5 | . . . ¨zN V Ņ1 £ |Ω N 1 | ¡ |Ω N 0 | V ķ0 p¡1q k |Ω N k | z N apzqtpzq ¡ apzq 1 b 0 p¡1, zq apzqtpzq ¡ apzq 1 apzq exp £ ¡ V ķ1 1 k tpz k qq apzqtpzq ¡ apzq 1 apzqz tpzq using Eq. (4.4)
apzqtpzq ¡åpzq using Eq. (

An element of Ω n,1 is obtained from an element of Ω n by marking a node with the symbol 1 . The marked node can be in one of the scalar aroma components or in one of the rooted tree components. Therefore, an element of Ω N n,1 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 tpzq as follows. Consider the connected component containing the marked node. The marked node lies either on the cycle or on one of the trees attached to the cycle. Those in the first group are enumerated by the rooted trees with one marked node: delete the edge of the cycle that points to the marked node (the construction is invertible). Those in the second group are enumerated by sequences of two rooted trees, each with a marked node: add edges from the first root to the first marked node and from the second root to the first root, and remove the mark from the first marked node. The rooted trees with one marked node are enumerated by sequences (of any length) of rooted trees: delete the outgoing edges from the nodes on the path from the root to the marked node.

Putting this together, the rooted trees with one marked node are enumerated by the following (sequence A000107), wpzq : tpzq tpzq 2 tpzq 3 ¤ ¤ ¤ tpzq{p1 ¡ tpzqq. The connected scalar aromas with one marked node are enumerated by wpzq wpzq 2 tpzq{p1 ¡ tpzqq 2 (sequence A038002). The scalar aromas with one marked node are enumerated by apzqtpzq{p1 ¡ tpzqq 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 pu, zqtpzq{p1 ¡ tpzqq 2 .

For type (ii), combining the components of an unmarked form with one fewer trees (i.e. an element of Ω n¡1 ), enumerated by ub 0 pu, zq, and a marked rooted tree, enumerated by tpzq{p1 ¡ tpzqq, gives the contribution from type (ii) forms as ub 0 pu, zqtpzq{p1 ¡ tpzqq.

Summing the results for type (i) and type (ii) forms gives the expression for b 1 pu, zq.

The functional forms in I 1 are associated with scalar aromas with one marked leaf (containing the symbol 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 |I N 1 | |Ω N ¡1 0,1 | and this gives the result for c 1 pzq.

Note that

Ω N ¡1 1,1 ! Ω N ¡1 0,1 ! I N 1
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:

V ņ0 p¡1q n |Ω n,1 | ¨¡ |I 1 | b 1 p¡1, zq ¡ c 1 pzq b 0 p¡1, zq tpzq 2 p1 ¡ tpzqq 2 ¡ apzqtpzqz p1 ¡ tpzqq 2 apzqz tpzq tpzq 2 p1 ¡ tpzqq 2 ¡ apzqtpzqz p1 ¡ tpzqq 2 0.
The following result extends Theorem 4.9 to the divergence-free case. 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, r Ω 0 , are enumerated by apzq zapzq{tpzq. This gives the result for r b 0 and, using Theorem 2.10, r spzq. For the second row, we first consider r Ω N 0,1 , the scalar aromas with N nodes, no self-loops, and one marked node indicated by 1 . These are enumerated by the self-functions of t1, . . . , N ¡ 1u with one marked node. That is, r

Ω N 0,1 ! Ω N ¡1 0,1 .
The construction is as follows. Consider an element of Ω N ¡1 0,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 r Ω N 0,1 has generating function zapzqtpzq{p1 ¡ tpzqq 2 .

For the rest of the second row, r Ω 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 apzq), n ¡ 1 unmarked rooted trees (enumerated by ub 0 pu, zq{apzq), and one marked rooted tree (enumerated by tpzq{p1 ¡ tpzqq).

The product of these three, plus the contribution from the forms of type (i), gives the result for b 1 pu, zq.

The elements of r

I N 1 Ip r Ω N 0,p q 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 r c 1 pzq.

Note that we now have five isomorphic spaces,

Ω N ¡1 1,1 ! Ω N ¡1 0,1 ! I N 1 ! r I N 1 1 ! r Ω N 0,1 .
Each space is enumerated by the self-functions of t1, . . . , N ¡ 1u with one marked node (and generating function zapzqtpzq{p1 ¡ tpzqq 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.

Bases of the kernel and image of

is B Ψ td H pϕφpt 2 q . . . φpt l¡1 qφpt l 1 q . . . φpt k qt 1 t l q, ϕ F 0 tru, t 1 ¤ ¤ ¤ t k T u.
Proof. The map φ is surjective: given any self-looped scalar, removing one of the self-loops gives a pre-image. Therefore, | Ker φ| | Ker d H |. Let us first determine Ker φ. Consider a self-looped scalar whose distinct self-looped connected components are φpt 1 q, . . . , φpt k q for trees t 1 ¤ ¤ ¤ t k ; it can be written ϕφpt 1 q . . . φpt k q where ϕ is a scalar. Its distinct preimages under φ are ϕφpt 1 q . . . t l . . . φpt k q for l 1, . . . , k. Therefore, the kernel of φ restricted to the span of these preimages has dimension k ¡ 1 and we consider tϕpt 1 φpt 2 q ¡ φpt 1 qt 2 qφpt 3 q . . . φpt k q, . . . , ϕpt 1 φpt k q ¡ φpt 1 qt k qφpt 2 q . . . φpt k¡1 qu as a basis of Ker φ. We now map Ker φ to Ker d H by ϕ l pt 1 φpt l q ¡ φpt 1 qt l q Þ Ñ d H pϕ l t 1 t l q, where ϕ l ϕφpt 2 q . . . φpt l¡1 qφpt l 1 q . . . φpt k q, extending by linearity. From exactness, the map is surjective. As | Ker φ| | Ker d H |, it is an isomorphism. Remark 4.8. One could wonder whether the set td H γ, γ F 2 u is a basis of Kerpd H q. This is not the case in general. For N 6, one finds for instance the following identity

d H p q 0.
The following result shows that the divergences are a graph over the self-looped scalars.

Theorem 4.9. For α Ω0 , let kpαq be the number of self-loops in α, and ρpαq be the non-selflooped scalar obtained from α as the sum of the redirection of all 1-loops to other nodes in all possible ways. Then the map

Ω0 Ñ Im d H , α Þ Ñ α p¡1q kpαq¡1 ρpαq
is an isomorphism, and generates a basis of Im d H .

Proof. Let V be the set of nodes with self-loops of α. From Proposition 3.6, we deduce that

for v V , α ¡ E v α Impd H q. If we have two nodes v, w V , then α ¡ E v α Impd H q and E v α ¡ E w E v α Impd H q, so that α ¡ E w E v α Impd H q.
By applying iteratively this process, we find that

α ¡ ¹ v V E v α α p¡1q kpαq¡1 ρpαq Impd H q.
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.

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

f i i f j j pf i f j j q i ¡ f i f j ij pf i f j j q i ¡ pf i f j i q j f i j f j i pf i f j j ¡ f j f i j q i f i j f j i . That is, f i i f j j ¡ f i j f j i is a divergence. We
γ γ ¡ 1 |γ| Eγ ¡ d H p h H γ d H ph H γ ¡ p h H γq,
where p γ does not contain any 1-loop and where we used Theorem 3.8. We deduce from Proposition 3.6 that p γ Impd H q. According to Corollary 4.11, we find p γ 0.

We now have two different ways to compute the horizontal homotopy operator on Ω 0,p . The first one, presented in Subsection 3.1, uses the Euler operators. The second one, in the spirit of [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF], is based on the repeated use of detaching and grafting operations on specific nodes. We emphasize that the expressions of the homotopy operator given by these two methods are different in general, but are always equal up to a divergence-free term. The two homotopy operators can produce both concise and tedious outputs, and they outperform each other in this manner on different forests. We refer the reader to Table 3 for some examples. This difference in the number of terms increases rapidly with the order; for instance, for the form γ , h H γ has 6 terms, while p h H γ has 26 terms. On the other hand, it is possible to find examples where p h H produces less terms than h H . It would be interesting to find a homotopy operator with a minimal number of terms in the output, or a procedure to simplify the outputs of a homotopy operator in the spirit of [1, Sect. IV.B]. Moreover, it is not known whether a similar approach in the divergence-free context could yield a different homotopy operator. This is matter for future work.

Explicit description of volume-preserving aromatic integrators

It is known that the only volume-preserving consistent B-series method is the exact flow [START_REF] Chartier | Preserving first integrals and volume forms of additively split systems[END_REF][START_REF] Iserles | B-series methods cannot be volumepreserving[END_REF]. In [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF], the question of the existence of a volume-preserving aromatic method is raised, where an aromatic method is a one-step integrator that has an expansion as an aromatic B-series. In [START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF], a methodology to create pseudo-volume-preserving interators is proposed, by substituting in a standard Runge-Kutta method the vector field f with an aromatic B-series. We give in this subsection an explicit expression of the form of a volume-preserving aromatic method, and we use it to prove that there does not exist any aromatic Runge-Kutta integrator and to discuss the form of a volume-preserving aromatic B-series method.

Consider a consistent one-step integrator (1.2) for solving the differential equation (1.1) with the assumption divpf q 0. We assume the integrator (1.2) has an expansion in aromatic Bseries given by the linear coefficient map a : r Ω 1 Ñ R; that is, its (formal) Taylor expansion has the form Φpy, hq y F pBpaqqphfq, Bpaq τF 1 apτ q σpτ q τ, where σpτ q is the cardinal of the set of automorphisms on the set of nodes V of τ that leave τ unchanged (see [START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF]). We call such a method an aromatic B-series method. If in addition the integrator reduces to a standard Runge-Kutta method when choosing a vector field f that satisfies F pγqpfq 0 for all γ F 1 zT , we call the integrator an aromatic Runge-Kutta method. 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 Bpbq satisfying Bpaq Bpbq Bpeq.

The operation is the substitution of B-series and b : r Ω 1 Ñ R is the coefficient map of the modified flow. It is known that a and b satisfy b e a, where is the substitution of B-series coefficients [START_REF] Calaque | Two interacting Hopf algebras of trees: a Hopf-algebraic approach to composition and substitution of B-series[END_REF][START_REF] Chartier | Algebraic structures of B-series[END_REF][START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF]. 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 epγq 0 if γ is composed of at least an aroma. The question raised in [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF] is the following: can we find an aromatic B-series method such that d H Bpbq 0, that is, such that the modified B-series Bpbq 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 SpanpT q Ψ r, SpanpT q r Ψ Spanp q.

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 °i c i t i , t i 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 prd H v °i c i θpt i q, where pr : Ω 0 Ñ Ω0 is the natural projection to the self-looped scalars, and the self-looped scalars θpt i q 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 volumepreserving integrator. For the first orders, the B-series of a volume-preserving aromatic B-series method has the form:

Bpaq 1 2 1 6 ¡ 1 6 ¡ 1 2 αp q © 1 2 αp q 1 24 ¡ 1 8 ¡ 1 2 αp q ¡ αp q ¡ 1 2 αp q © ¡ 1 24 ¡ 1 2 αp q 1 2 αp q ¡ 1 4 αp q © 1 2 αp q 1 2 αp q ¡ 1 24 ¡ 1 2 αp q ¡ 1 4 αp q © ¡ 1 2 αp q αp q 1 2 αp q © . . .
Note that the coefficients of the bamboo trees (or tall trees) BT t , , , . . . u coincide with the ones of the exact flow. This fact has been noticed for standard B-series in particular in [START_REF] Kang | Volume-preserving algorithms for source-free dynamical systems[END_REF] (see also [START_REF] Hairer | Geometric numerical integration[END_REF]Lemma IV.3.2]). We deduce from this observation the following result.

Theorem 4.18. An aromatic Runge-Kutta method cannot be volume-preserving.

Proof. The only bamboo tree that appears in solenoidal forms is . Indeed d ¦ H vanishes on BT as the only solenoidal forms where a bamboo tree can appear are of the form d H τ 1 τ 2 where τ 1 , τ 2 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 Bpeq on the bamboo trees BT . A Runge-Kutta integrator cannot be exact on all bamboo trees. As the aromatic forests represent different elementary differentials (see [START_REF] Iserles | B-series methods cannot be volumepreserving[END_REF]), any aromatic integrator that reduces to a standard Runge-Kutta integrator when sending the aromas to zero cannot be volume-preserving.

The methodology proposed in [START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF]Sect. 7] to obtain volume-preservation of high order and the approach in [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF]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, [START_REF] Berland | B-series and order conditions for exponential integrators[END_REF][START_REF] Hochbruck | Explicit exponential Runge-Kutta methods for semilinear parabolic problems[END_REF][START_REF] Luan | Exponential B-series: the stiff case[END_REF]) when sending the aromas to zero. This calls for future works that study the substitution law and the variational bicomplex directly on aromatic exponential B-series, in order to find a volume-preserving aromatic B-series method.

Definition 2 . 4 .

 24 Let γ F n,p be a forest, r a root of γ, and u V (possibly equal to r), then D rÑu γ returns a copy of γ where the node r is now a predecessor of u. We define the operator D r γ °uV D rÑu γ. Let γ F n,p and v V , we define γ vÑ k as the forest obtained by replacing the node v by a new covertex k . Similarly, γ k Ñ is the forest obtained by replacing the covertex k by a vertex.

Remark 3 . 2 .

 32 and D ¥ D are different operators. Instead, straightforward combinatorics yield ¢ n p D n D p D n¡p . (3.1) However, for u, v R, the operations D u and D v commute. The notations introduced in Definition 3.1 closely depend on the forest γ 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 γ. For instance, in p¡1q |Πpvq| γ v , Πpvq denotes the number of predecessors of v in γ, not in γ v .

  |I|q εpσq pq nqpn ¡ 1q! D I E I,r v D r σ ¡1 pnq Ñu σγ V q0 σSn Πpvqtr k utruIJ |I|q εpσq pq nqpn ¡ 1q! D I E I,r v D r σ ¡1 pnq Ñv σγ.The expression of D I E I,r v D rÑu γ satisfies for u V :

σSn νS n 1 V q0 σSn νS n 1

 11 ΠpvqtruIJ |I|q εpσνq pq n 1qpn!q 2 νD I E I,r v σγ ΠpvqtruIJ |I|q εpσνq pq n 1qpn!q 2 D r pσνq ¡1 pn 1q νD I E I,

Figure 3 :

 3 Figure 3: The augmented aromatic bicomplex.

Theorem 4 . 4 .

 44 Let r b p pu, zq, r c p pzq, and r spzq be the divergence-free analogues of the generating functions b p pzq, c p pzq, and spzq. Then r b 0 pu, zq zb 0 pu, zq tpzq , r b 1 pu,

Table 3 :

 3 Comparison of the horizontal homotopy operators on Ω 0 for the first orders.

Theorem 4 . 17 .

 417 If Bpaq is the aromatic B-series of a consistent volume-preserving integrator, then there exists η r Ω 2 such that the modified flow is a B-series of the form Bp d H ηq and Bpaq is given by the substitution Bpaq p dηq Bpeq.More precisely, there exists a coefficient map α : r Ω 2 Ñ R such that a is given by a p ¦ αd ¦ H A σ q e, where A σ γ 1

  extend d H and d V on Ω n,p by linearity into d H : Ω n,p Ñ Ω n¡1,p and dV : Ω n,p Ñ Ω n,p 1 .The aromatic forms in Kerpd H q are called solenoidal forms, Ψ Kerpd H | Ω 1 q are the solenoidal combinations of trees, and Ψ N Kerpd H | Ω N 1 q are the solenoidal combinations of trees of order N . The operator d H : Ω 1 Ñ Ω 0 is often called the divergence of an aromatic tree, as for γ Ω 1 , d H satisfies divpgq F pd H pγqqpfq, where F pγqpfq °i g i f i . The operator d H : Ω 1 Ñ Ω 0 is often called the divergence of an aromatic tree.Example. Consider γ 1 Ω 1 , γ 2 Ω 2 , and γ 3 1 Ω 2,1 , then

Table 1 :

 1 for examples). Dimensions of the space of solenoidal forms in both contexts for the first orders N (see Theorems 4.2 and 4.4). Note how

	N	1 2 3 4 5	6	7	8	9	10	11	12	13	14
	|Ω N 1 | 1 2 6 16 45 121 338 929 2598 7261 20453 57738 163799 465778
	| ΩN 0 | 1 2 5 13 34 90 243 660 1818 5045 14102 39639 111982 317533
	|Ψ N | 0 0 1 3 11 31 95 269 780 2216 6351 18099 51817 148245
	| r Ψ N | 1 0 1 2 7 16 48 123 346 937 2626 7284 20533 57804

  Let R 0 R be a given subset of roots that we call the set of detached nodes. For γ F n,p , r R and u V with u $ v , D rÑu γ F n¡1,p is the marked aromatic forest obtained by adding the edge linking r to u to γ. The set of detached roots of

Definition 3.1. A graph γ is a marked aromatic forest if it is an aromatic forest with exactly one of its node that holds the symbol . We write F n,p the set of marked aromatic forests with n roots and p covertices. Given γ F n,p , γ F n,p is the same forest γ where the symbol is removed. Let γ F n,p or γ F n,p , where v V is the node with the symbol .

  zq zb 0 pu, zq tpzq u ¡ utpzqq p1 ¡ tpzqq 2 , r c 1 pzq zc 1 pzq,

	First row: |Ω N n |			Second row: |Ω N n,1 |		|I N 1 |
	N	n					n
	4 3	2	1	0	4	3	2	1	0
	1 0 0	0	1	1	0	0	0	1	1	0
	2 0 0	0	2	3	0	0	1	4	4	1
	3 0 0	1	6	7	0	0	4	15	15	4
	4 0 0	3	16	19	0	1	16	52	52	15
	5 0 0 11	45	47	0	5	57	175	175	52
	6 0 2 33 121 130	0	22 197	571	571	175
	7 0 7 102 338 343	2	85 654 1838 1838	571
	8 0 29 298 929 951 11 310 2137 5834 5834 1838
	9 1 99 878 2598 2615 53 1078 6859 18363 18373 5834
		r spzq z zspzq tpzq z zpapzqtpzq ¡åpzqq tpzq	.

Table 2 :

 2 Dimensions of the bottom two rows of the augmented aromatic bicomplex for orders 1 to 9.

  d H and d ¦ In this subsection, we work specifically with d H : Ω 1 Ñ Ω 0 and we describe the image and the kernel of d H and d ¦ H . We recall that d ¦ is the dual map of d H . We described the dimension of Ψ Ker d H in Subsection 4.2. The following result describes the dimensions of Im d H , Ker d ¦ H , and Im d ¦ The scalar divergences Im d H have dimension | Ω0 |. The conditions Ker d ¦ H that a scalar must satisfy to be a divergence have dimension |Ω 0 |. The conditions Im d ¦ H that a vector must satisfy to be divergence-free have dimension | Ω0 |. Proof. Recall the fundamental theorem of linear algebra for a linear map A : V Ñ W :| Im A| | Im A ¦ | rankA, | Im A| | Ker A| |V |, | Im A ¦ | | Ker A ¦ | |W|.space Im A ¦ AnnpKer Aq 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 Ker A ¦ AnnpIm Aq is the annihilator of the image of A. Choosing A d H gives the result. Remark 4.6. Consider the aromatic tree γze F 1 obtained by cutting one edge e E of γ F 0 , m 1 pγ, eq the coefficient of γ in d H pγzeq, and m 2 pγ, eq the number of edges p e E of γ such that γzp e γze. Then, for γ Ω 0 , d ¦ eq m 2 pγ, eq pγzeq ¦ .We now construct bases of Ker d H , Im d H , Ker d ¦ Let φ : Ω 1 Ñ Ω0 be defined by attaching a self-loop to the root, extending by linearity. Define any total order on T , then a basis of Ψ Ker d H

	H .
	Theorem 4.5. The Theorem 4.7.

H H : Ω ¦ 0 Ñ Ω ¦ 1 H γ ¦ satisfies d ¦ H γ ¦ ȩE m 1 pγ, H ,

and Im d ¦

H . We start with the basis of the solenoidal forms.

Corollary 4.11. No

  describe this comparison with integration by parts further in Subsection 4.4. non-trivial combination of non-self-looped scalars in Ω 0 is a divergence. The conditions to be a divergence Ker d ¦ H are a graph over the dual of the non-self-looped forms in Ω 0 .In the following theorem, this graph is realized explicitly. Let p E E be a set of edges of the scalar aroma β Ω 0 and let βz p For γ Ω 0,p , the output p h H γ of the algorithm is the horizontal homotopy operator, up to a divergence-free term, that is,d H ph H γ ¡ p h H γq 0.

	Proof. After the algorithm terminates, p γ is given by
	p
	Corollary 4.12. Theorem 4.13. E be β
	with edges p E replaced by self-loops. Let mpβ, p Eq be the number of ways that redirecting self-loops of βz p E results in β. Let π : Ω 0 Ñ Ω 0 be defined by
	πpβq	p	p¡1q | p E| mpβ, p Eqβz p E
		EE
	Then	tπpβq, β Ω ¦ 0 u

is a basis of AnnpIm Aq, the conditions that a scalar must satisfy to be a divergence. Theorem 4.15.

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.

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.

Proof. The construction is directly related to that in Theorem 4.9. The conditions for px, yq to lie on the graph tx, Axq : x R n u R n ¢ R m are y ¡ Ax 0, A R m¢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 β we need to determine the coefficient of β in ρpαq for each self-looped scalar α. This is the expression for π: the term p E r gives β, and the terms from non-empty sets of edges p E give the α's that can give rise to β.

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

H , the conditions that a form in Ω 0 must satisfy to be divergence free.

Proof. For any linear map A : V Ñ W , xA ¦ w ¦ , vy xw ¦ , Avy 0 for all w ¦ W ¦ when v 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

Let pr : Ω 0 Ñ Ω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, pr

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 [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF][START_REF] Laurent | Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds[END_REF] on exotic aromatic B-series (see also [START_REF] Bronasco | Exotic B-series and S-series: algebraic structures and order conditions for invariant measure sampling[END_REF]). 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 Ω 0,p .

Let γ Ω N 0,p a linear combination of forests, let τ F N 0,p one of these forests and v a vertex of τ on a 1-loop. We denote apτ q the coefficient of τ in γ, and θ v pτq the forest τ where we remove the edge linking v to itself and transform v into a root. The alternative horizontal homotopy operator p h H on Ω 0,p is given by the following algorithm on Ω N 0,p , and extended to Ω 0,p by linearity.

Homotopy operator

while there is a 1-loop in a forest τ of p γ on a vertex v do

end while return p h H γ

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 [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF], as there is an extra term involved in the integration by parts process.

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 [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF]. 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 [START_REF] Anderson | The variational bicomplex[END_REF] and references therein). It would be interesting to generalise some of the existing results in the context of aromatic forms. For instance, one could try to find simpler expressions for the homotopy operators (see Subsection 4.4), to find an augmented bicomplex in the divergencefree case (see Subsection 4.1). Two major applications of the variational bicomplex are the

Appendices A First solenoidal forms

We write the generators of the solenoidal forms r Ψ 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 γ for γ r F N 2 , and by adding the trivial tree . For N ¤ 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 and5 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 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 4: Generators of the solenoidal forms r

Ψ N for the first orders N . ,