---
abstract: 'In earlier work, Lomeli and Meiss [@lomeli09gff] used a generalization of the symplectic approach to study volume preserving generating differential forms. In particular, for the ${\mathbb{R}}^3$ case, the first to differ from the symplectic case, they derived thirty-six one-forms that generate exact volume preserving maps. In [@xue2014gf], Xue and Zanna studied these differential forms in connection with the numerical solution of divergence-free differential equations: can such forms be used to devise new volume preserving integrators or to further understand existing ones? As a partial answer to this question, Xue and Zanna showed how six of the generating volume form were naturally associated to consistent, first order, volume preserving numerical integrators. In this paper, we investigate and classify the remaining cases. The main result is the reduction of the thirty-six cases to five essentially different cases, up to variable relabeling and adjunction. We classify these five cases, identifying two novel classes and associating the other three to volume preserving vector fields under a Hamiltonian or Lagrangian representation. We demonstrate how these generating form lead to consistent volume preserving schemes for volume preserving vector fields in ${\mathbb{R}}^3$.'
author:
- Olivier Verdier
- Huiyan Xue
- Antonella Zanna
bibliography:
- 'volumepreserving.bib'
title: 'A classification of volume preserving generating forms in $\RR^3$'
---

Introduction and background {#laha}
===========================

The scope of this paper is the study volume preserving generating forms, with the ultimate goal of exploiting these differential forms to obtain consistent numerical methods that preserve volume for arbitrary volume preserving vector fields. This task is particularly hard: there exist *no-go* theorems [@chartier2007preserving; @iserles07] stating that it is not possible to construct volume preserving methods for generic $n$-dimensional volume preserving vector fields within the class of B-series methods, a class that includes classical integrators like Taylor-expansion based methods, Runge–Kutta methods and multistep methods. On the other hand, volume preserving methods can be constructed using the technique of *splitting* [@McLachlan2002]. Splitting methods correspond to P-series (“P” for partitioned systems), a generalization of B-series. Several splitting techniques can be adopted. The earliest and best known splitting consists in decomposing the vector field in 2D Hamiltonian sub-systems [@kang95vpa], which are then solved by a symplectic method. More recently, research has focussed on explicit splitting methods for classes of vector fields, like polynomial or trigonometric, which are wide enough to include most interesting cases, but not as large as the space of all possible vector fields [@quispel03evp; @mclachlan04egi; @mclachlan09evp; @Xue:2012mf; @zanna2014tensor].

Differently from the symplectic case, the generating form approach to generate volume preserving numerical methods is not well understood. Earlier work by [@shangXXgf1; @Shang1994] extends the Hamiltonian technique of [@F1986; @FW1989], using linear maps in the product space, to volume preserving forms, thus obtaining an equivalent of the Hamilton-Jacobi differential equation [@Shang1994]. To obtain a first and second order scheme, Shang had to impose simplifying conditions, requiring the transformation matrix to be a special case of Hadamard matrix. However, the numerical integrators by this approach are quite complicated, as they are defined via implicit maps, although the approach is valid for arbitrary vector fields. Another generating-functions related approach is due to [@GRW199526]: through a special combinations of explicit and implicit maps, Quispel shows that the resulting method is volume preserving. This is a “correction method”: starting from an arbitrary numerical integrator, one has to consider an extra term (the corrrection) for volume preservation. The above mentioned two approaches do not use differential forms directly, rather, they use the equivalent condition on the determinant of the Jacobian of the map. For this reason, they use the terminology of “generating functions” rather than “generating forms”.

More recently, L[ó]{}meli and Meiss [@lomeli09gff] have studied the problem of volume preserving maps using differential forms and generalization of the symplectic approach. They discussed in detail the ${\mathbb{R}}^3$ case, the first to differ from the symplectic case, and described how the generic volume preserving maps can be described by thirty-six one-forms. That paper paves the background for our investigations. In particular, we are interested in understanding how these differential forms are associated to numerical methods (if any) and whether some of these forms can lead to new techniques to obtain volume preserving maps. These questions were partially addressed in [@xue2014gf], where six of the thirty-six differential one forms were identified and associated to splitting methods. The scope of this paper is to discuss and classify the remaining cases.

The main result of the paper is the reduction of the thirty-six cases to five essentially different cases, using equivalence relations (global variable renaming and numerical adjoints). Thereafter, these five cases are classified and three of them associated to known techniques, based on Hamiltonian and/or Lagrangian formalism. In other words, the generating forms are associated to symplectic splitting methods like Symplectic Euler (SE) for Hamiltonian systems, or Discrete Lagrangian (DL) methods for appropriate Lagrangian functions, or a combination of both. We further identify two special classes, $S_1$ and $S_2$, that, to our knowledge, do not have a straightforward and direct mechanical interpretation. It is these two special classes that are of particular interest in the search of new volume preserving algorithms. We give explicit formulas of generating one-forms for the class $S_1$ and $S_2$ corresponding to volume preserving first order methods in the specific case of linear vector fields. A general approach is still unknown and will be the subject of future investigation.

Background and notation {#sec:background}
-----------------------

We consider a differentiable manifold $\mathcal{M}$. Let $\omega$ be a non degenerate $k$-differential form on $\mathcal{M}$, that is, a $k$-linear map, completely skew-symmetric with respect to its arguments. For each $p \in \mathcal{M}$, $\omega(p): (T_p\mathcal{M})^{\times k} \to {\mathbb{R}}$, namely the differential form takes as argument $k$ tangent vectors and returns a number. The coefficients of the differential form might depend on $p\in \mathcal{M}$. Let $d \omega$ be the $k+1$ form obtained with the usual rules of external derivation. Recall that $\omega$ is *closed* if $d\omega = 0$ and that $\omega$ is an exact differential, or simply *exact*, if $\omega = d\nu$, where $\nu$ is a $k-1$ form, called a *primitive*. By application of Stokes’ theorem, $\int_S d \omega = \int_{\partial S} \omega$, valid on any oriented manifold with oriented boundary $\partial S$, to the differential form $d\omega$, it follows that $d^2\omega=0$. Therefore, any exact form, $\omega = d \nu$, is closed, i.e. $d\omega = d^2\nu = 0$. The reverse statement is not true in general, but it holds on contractible manifolds, as explained from the following lemma.

\[th:poincare\] A closed form ($d\omega=0$) is locally exact ($\omega=d\nu$), that is, there is a neighborhood $U$ about each point on which $\omega=d\nu$. The statement is globally true on contractible manifolds.

In the sequel, we focus on volume forms and their primitives. We will assume, otherwise stated, that the differential forms are *non-degenerate*. This means that the coefficients of the form are never simultaneously zero.

A volume form $\Omega$ on a manifold $\mathcal{M}$ is preserved by a $C^{1}$-map $\mathbf{f}: \mathcal{M}\mapsto \mathcal{M}$ if $$\mathbf{f}^{*}\Omega=\Omega,
    \label{eq:vol_pres}$$ where $\mathbf{f}^{*}$ denotes the pull-back of $\mathbf{f}$. The map $\mathbf{f}$ is said to be *canonical* or *volume preserving*.

In what follows, we let $\mathcal{M}={\mathbb{R}}^n$. Let $\nu$ be any primitive form of $\Omega$, i.e. $d\nu=\Omega$. Then, condition becomes $\mathbf{f}^{*}d\nu-d\nu=0$ and implies $d(\mathbf{f}^{*}\nu-\nu)=0$, hence $\mathbf{f}^{*}\nu-\nu$ is the exact differential of a $n-2$ form, as a consequence of Lemma \[th:poincare\]. This motivates Definitions \[def:nu\_exact\] and \[df\] below, see [@lomeli09gff].

\[def:nu\_exact\] Let $\nu$ be a primitive of the volume form $\Omega$ and $\mathbf{f}: {\mathbb{R}}^n\mapsto {\mathbb{R}}^n$ an exact volume preserving diffeomorphism such that $$\mathbf{f}^{*}{\nu}-\nu=d\lambda
 \label{eq:nu_exact},$$ for a $n-2$ form $\lambda$. The differential form $\lambda$ is called a *generating form with respect to $\nu$.*

Primitives of forms are not uniquely determined: by choosing $\tilde \nu$ another primitive of $\Omega$, the volume preservation condition can be written as $\mathbf{f}^{*}d\nu-d\tilde\nu=0$. A procedure similar to the one just described above leads to:

\[df\] Let $\nu, \tilde \nu$ be two primitives of a volume form $\Omega$, i.e. $d\nu=d\tilde{\nu}=\Omega$ and $\mathbf{f}: {\mathbb{R}}^n\mapsto {\mathbb{R}}^n$ an exact volume preserving diffeomorphism such that $$\mathbf{f}^{*}\tilde{\nu}-\nu=d\lambda\label{eq:0806031},$$ for a $n-2$ form $\lambda$. The $n-2$ differential form $\lambda$ is called a *generating form with respect to $(\nu, \tilde \nu)$.*

We will consider the choice of canonical coordinates $x_1, \ldots, x_n$ in ${\mathbb{R}}^n$ and denote by $\mathbf{x} = (x_1, \ldots, x_n)^T$ the original (old) variables. Given a volume preserving map $\mathbf{f}$, we will denote the transformed (new) variables by uppercase letters, i.e. $\mathbf{X} = (X_1, \ldots, X_n)^T = \mathbf{f}(\mathbf{x})$.

Volume preservation in ${\mathbb{R}}^2$ is equivalent to preservation of area and volume forms are the same as symplectic forms. This case is well understood. The case $n=3$ is the first case for which volume forms and symplectic forms are different.

Generating forms in ${\mathbb{R}}^3$ {#volsec}
====================================

In [@lomeli09gff], Lomeli and Meiss studied in detail exact volume preserving mappings and generating forms in ${\mathbb{R}}^3$. Starting from and using canonical coordinates, they showed that for each choice of primitives $(\nu, \tilde{\nu})$ there are four different generating one-forms, $$\lambda= \phi dx_l + \Phi d X_m,  \qquad \phi \in \{A,B\}, \quad \Phi \in \{C,D\}\qquad l, m \in \{1,2,3\},
    \label{eq:oneform}$$ which they identified using four generating functions $A, B, C$ and $D$.

As $\nu$ and $\tilde{\nu}$ can be chosen in three different ways ($x_3dx_1\wedge dx_2$, $x_2dx_3\wedge dx_1$ and $x_1dx_2\wedge dx_3$), this approach gives a total of thirty-six generating one-forms. Four of them, corresponding to $\nu= \tilde \nu = x_3 dx_1 \wedge d x_2$, are shown in Table \[tab:123123\]. Each cell in the table is described by: a generating one-form $\lambda$; two *determining conditions*, determining a lowercase and an uppercase variable; a *compatibility condition* and two *twist conditions* to guarantee that the three equations are solvable and they give rise to a well-defined volume preserving map. Altogether, one there are nine such tables, obtained by even permutations of the $x_1,x_2,x_3$ and the $X_1,X_2,X_3$ variables.

  ---------------------------- ------------------------------------------------------------------------------------------ ------------------------------------------------------------------------------------------
  $X_3 dX_1 \wedge dX_2$ $-$                                                                                              

  $x_3 dx_1 \wedge dx_2$                                                $A dx_1$                                                                                   $B dx_2$

  $C dX_1$                                                        $\begin{array}{ccc}                                                                        $\begin{array}{ccc}
                                                           \lambda = A(x_1, x_2, X_1) dx_1 \\                                                         \lambda = B(x_1, x_2, X_1) dx_2 \\
                                                   \qquad \quad \mbox{}+ C(x_1, X_1, X_2) d X_1\\[5pt]                                       \qquad \quad \mbox{}+ C(x_2, X_1, X_2) d X_1\\[5pt]
                                                                 x_3=\partial_{x_2}A \\                                                                    x_3=-\partial_{x_1}B \\
                                                         \partial_{X_1} A = \partial_{x_1} C\\                                                      \partial_{X_1} B = \partial_{x_2} C\\
                                                               X_3 = -\partial_{X_2} C\\                                                                  X_3 = -\partial_{X_2} C\\
                                \frac{\partial X_1}{\partial x_3} \not=0, \quad \frac{\partial x_1}{\partial X_3} \not=0   \frac{\partial X_1}{\partial x_3} \not=0, \quad \frac{\partial x_2}{\partial X_3} \not=0
                                                                     \end{array} $                                                                              \end{array} $

  $D d X_2$                                                       $\begin{array}{ccc}                                                                        $\begin{array}{ccc}
                                                           \lambda = A(x_1, x_2, X_2) dx_1 \\                                                         \lambda = B(x_1, x_2, X_2) dx_2 \\
                                                  \qquad\quad  \mbox{}+ D(x_1, X_1, X_2) d X_2\\[5pt]                                         \qquad \quad \mbox{}+ D(x_2, X_1, X_2) d X_2\\[5pt]
                                                                 x_3=\partial_{x_2}A \\                                                                    x_3=-\partial_{x_1}B \\
                                                         \partial_{X_2} A = \partial_{x_1} D\\                                                      \partial_{X_2} B = \partial_{x_2} D\\
                                                                X_3 = \partial_{X_1} D\\                                                                   X_3 = \partial_{X_1} D\\
                                \frac{\partial X_2}{\partial x_3} \not=0, \quad \frac{\partial x_1}{\partial X_3} \not=0   \frac{\partial X_2}{\partial x_3} \not=0, \quad \frac{\partial x_2}{\partial X_3} \not=0
                                                                     \end{array} $                                                                              \end{array} $
  ---------------------------- ------------------------------------------------------------------------------------------ ------------------------------------------------------------------------------------------

  : The four basic types of generating 1-forms $\lambda$ for $\nu = \tilde \nu = x_3 dx_1 \wedge d x_2$, adapted from [@lomeli09gff]. All the other tables are obtained by applying cyclic even permutations to the variables $(x_1,x_2,x_3)$ in $\nu$ and $(X_1,X_2,X_3)$ in $\mathbf{f}^*\tilde{\nu}$.[]{data-label="tab:123123"}

In [@xue2014gf], Xue and Zanna studied these generating forms with the goal of associating them to volume preserving vector fields and numerical volume preserving integrators. They succeeded in identifying six cases and associating them to splitting methods in using two potential functions: each of the functions gave rise to a two-dimensional Hamiltonian approximated by a Symplectic Euler (SE) method. The six cases are not fundamentally different. First of all, having chosen one case, two of the other cases correspond to a global variable renaming (say, $(x_1,x_2,x_3) \to (x_2,x_3,x_1)$ and $(X_1,X_2,X_3) \to (X_2,X_3,X_1)$ simultaneously). Relabeling the variables in a numerical method does not give a new numerical method. Secondly, the remaining three cases were obtained by exchanging lower cases and reversing time, in other words, they corresponded to the *adjoint* numerical methods of the previous three cases. As it is known how to obtain the adjoint of a given method [@hairer87sod], these cases are not interesting per se either. All this indicates that there is redundancy in the thirty-six cases. To classify and understand which one forms are related to known methods and which forms can lead to genuinely new approaches, we need to establish equivalence classes, so that our search can be restricted to a single differential form for each class.

Problem statement
=================

We consider the ordinary differential equation $$\dot{\mathbf{x}}=\mathbf{a}(\mathbf{x}),  ~~\mathbf{x}(0)=\mathbf{x}_{0}, 
\label{eq:vf}$$ where $\mathbf{x} \in \mathcal{M}$ and $\mathbf{a} :\mathcal{M} \to
T_\mathbf{x} \mathcal{M}$, $\mathbf{a}(\mathbf{x}) = [a_1(\mathbf{x}),
\ldots, a_n(\mathbf{x})]^T$, is a smooth vector field. We denote by $\mathbf{a}^t$ the flow of and by $\omega$ a $k$-form.

Derivatives of differential forms along the flow are called *Lie derivative* and defined as $$L_{\mathbf{a}} \omega = \frac{d}{dt} (\mathbf{a}^t)^* \omega \Big|_{t=0}
\label{eq:lie_deriv}$$ Let $\Omega$ be a volume form on $\mathcal{M}$. We say that the vector field $\mathbf{a}$ is *volume preserving* if $$L_{\mathbf{a}} \Omega =0.
    \label{eq:volpres_a}$$ The flow $\mathbf{a}^t$ is volume preserving if its vector field $\mathbf{a}$ is volume preserving. By , this implies that $$(\mathbf{a}^t)^* \Omega = \Omega.$$ Let $i_\mathbf{a}$ a contraction, that is, for any $k$-form $\omega$, $i_\mathbf{a} \omega$ is the $k-1$ form $\omega (\mathbf{a}, \mbox{} \cdot \mbox{})$ obtained by inserting $\mathbf{a}$ in the first slot. By Cartan’s formula for Lie derivatives, $$L_{\mathbf{a}} \Omega = d (i_\mathbf{a} \Omega) + i_{\mathbf{a}} d \Omega,
\label{eq:Cartan}$$ it follows that a vector field $\mathbf{a}$ is volume preserving if $d (i_\mathbf{a} \Omega) =0$, that is, $ i_\mathbf{a} \Omega $ is closed (as $d \Omega =0$, being $\Omega$ a $n$-form).

Let $\lambda$ be a $n-2$ form and $\Omega$ a volume form on $\mathcal{M}$. A vector field $\mathbf{a}$ on $\mathcal{M}$ is *exact volume preserving* with respect to the potential form $\lambda$ if $$i_\mathbf{a} \Omega = d \lambda.
    \label{eq:exactvp}$$

In particular, it follows from Poincar[é]{}’s lemma \[th:poincare\] that, when $\mathcal{M}={\mathbb{R}}^n$, globally defined volume preserving vector fields are also exact.

Let us consider the case $n=3$ in more detail. For any vector $\mathbf{v}_i \in {\mathbb{R}}^3$, $i=1,2,3$, we have $\Omega(\mathbf{v}_1, \mathbf{v}_2,\mathbf{v}_3)= dx_1 \wedge dx_2 \wedge dx_3 (\mathbf{v}_1, \mathbf{v}_2,\mathbf{v}_3)= \det[\mathbf{v}_1, \mathbf{v}_2,\mathbf{v}_3]$. By direct computation, we see that $\Omega (\mathbf{a}, \mathbf{v}, \mathbf{w}) = [a_1 dx_2 \wedge dx_3 + a_2 dx_3 \wedge d x_1 + a_3 d x_1 \wedge d x_2] (\mathbf{v}, \mathbf{w})$ for any $\mathbf{v}, \mathbf{w}$. We deduce that $$i_\mathbf{a}\Omega = a_1 dx_2 \wedge dx_3 + a_2 dx_3 \wedge d x_1 + a_3 d x_1 \wedge d x_2.$$ There are three natural choices of the one-form $\lambda$, i.e. $\lambda_i = F^{i}(x_1,x_2,x_3) d x_i$, $i=1,2,3$, where the $F^i$s are arbitrary function. Then, $$d \lambda_i = \sum_{j=1}^3 \partial_{x_j} F^i d x_j \wedge d x_i, \qquad i=1,2,3,$$ and, by , we deduce that a volume preserving vector field, exact with respect to the form $\lambda_i$, must have the following form: $$\begin{aligned}
\label{eq:vf1}
 F^1 d x_1 &: \qquad
 \begin{array}{ccc}
a_1 &=& 0,\\
a_2 &=& \partial_{x_3} F^{1},\\
a_3 &=& -\partial_{x_2} F^1,
\end{array} \\
\label{eq:vf2}
 F^2 d x_2 &: \qquad
\begin{array}{ccc}
 a_1 &=& -\partial_{x_3} F^2,\\
a_2 &=& 0,\\
a_3 &=& \partial_{x_1} F^2,
\end{array}\\
\label{eq:vf3}
F^3 d x_3 &:\qquad
\begin{array}{ccc}
a_1 &=& \partial_{x_2} F^3,\\
a_2 &=& -\partial_{x_1} F^3,\\
a_3 &=& 0.\end{array}\end{aligned}$$ A generic three-dimensional volume preserving vector field will be a combination of - above. Note that only two of them are linearly independent, the same yields for the choices of $\lambda_i$: for instance, as long as $F^1, F^2$ depend on all variables, $d \lambda_1$ and $d \lambda_2$ generate all the $dx_i \wedge dx_j$, $i, j=1,2,3$ although $\lambda_1$ contains only $d x_1$ and $\lambda_2$ only $dx_2$.

This is a consequence of a well known result for volume preserving flows in ${\mathbb{R}}^n$: any $n$-dimensional volume preserving differential equation is described by $n-1$ independent potential functions (see for instance [@McLachlan2002; @kang95vpa; @FW1989; @xue2014gf] and discussion therein). One of the earliest normalization of the $n-1$ independent potential functions is due to Weyl [@weyl40tmo].

It is natural to draw a connection between the generating forms and the volume preserving vector field $\vf$. In particular, we address the following questions:

-   How do the generating functions $\phi, \Phi$ relate to the vector field $\vf$ and the $F^{i}$s, if there is any relation?

-   Can the generating forms be naturally associated to volume preserving numerical methods for whenever the vector field $\mathbf{a}$ in is volume preserving?

An important property of numerical methods is *consistency*. Numerical methods introduce a discrete time step $h$. A fundamental property required to the method is that, in the limit $h\to 0$, $\frac{\mathbf{X}-\mathbf{x}}{h} = \vf$, that is, the method solves the given differential equation.

We now define what is the main goal of this paper: to find suitable potential functions $\pot$ and $\Pot$ for a divergence free vector field $\vf$ which also are a solution of the *consistency problem*, below rephrased in the formalism of this paper for convenience.

The main goal of this section is to show that, subject to consistency, there are only five different classes of generating volume forms in ${\mathbb{R}}^3$.

Defining equations and compatibility conditions
-----------------------------------------------

Consider two arbitrary smooth functions $$\pot,\Pot \colon {\mathbb{R}}^3 \to {\mathbb{R}},$$ and an arbitrary sign, which will be denoted by $\sig$ in the equations.

Define the corresponding mapping $\volpot\colon (\yy1,\yy2,\yy3) \mapsto(\YY1,\YY2,\YY3)$ implicitly by the following equations.

\[eq:volpot\] $$\begin{aligned}
        \partial_{2} \pot(\YY3, \yy2, \yy3) &=\yy1 \label{eq:volpot1}  \\
        \partial_{3} \Pot(\YY3,\YY2,\yy3) &\sig \partial_{1} \pot(\YY3,\yy2, \yy3) = 0 \label{eq:volpot2}\\
        \YY1 &= \partial_{2} \Pot(\YY3, \YY2, \yy3). \label{eq:volpot3}\end{aligned}$$

The mapping $\volpot$ is well defined as soon as $\pd[\Pot]{32} \neq 0$ and $\pd[\pot]{21} \neq 0$.

We define the of a permutation $\perm\colon\set{1,2,3}\to\set{1,2,3}$ on an element $\mathbf{x}\in{\mathbb{R}}^3$ by $$\begin{aligned}
    (x \act \perm)_i \coloneqq x_{\perm(i)}
    .\end{aligned}$$ Consider two permutations $\perm$ and $\Perm$: $$    \perm,\Perm \colon \set{1,2,3} \to \set{1,2,3}  .$$ For any map $\volpota \colon {\mathbb{R}}^3 \to {\mathbb{R}}^3$, we define the map $\permact{\perm,\Perm}{\volpota}$ by $$\begin{aligned}
    \paren[\big]{\permact{\perm,\Perm}{\volpota}}(x) \coloneqq \volpota(x \act \perm) \act \Perm\inv
    .\end{aligned}$$

The consistency problem
-----------------------

\[def:conssol\] We say that a pair of maps $\pot, \Pot \colon \mathfrak{X}_0({\mathbb{R}}^3) \to \mathcal{C}^{\infty}({\mathbb{R}},{\mathbb{R}}^3)$, where $\mathfrak{X}_0({\mathbb{R}}^3)$ denotes the space of divergence-free vector fields on ${\mathbb{R}}^3$, is a for the pair of permutations $(\perm,\Perm)$ if the map $h \mapsto \permact{\perm,\Perm}{\volpot[\pot(h\vf),\Pot(h\vf)]}$ is *consistent* with $\vf$.

This means that $$\begin{aligned}
    \lim_{h\to 0} \frac{1}{h}\paren[\Big]{\paren[\big]{\permact{\perm,\Perm}{\volpot[\pot(h\vf),\Pot(h\vf)]}}(\mathbf{x}) - \mathbf{x}} = \vf
    .\end{aligned}$$

Note that for any fixed vector field $\vf$, the element $\pot(\vf)$ is a potential, i.e., $\pot(\vf)$ is itself a function from ${\mathbb{R}}^3$ to ${\mathbb{R}}$. Note also that a method obtained from will automatically have order one.

A classification and description of the five generating volume forms in ${\mathbb{R}}^3$
========================================================================================

Having proved that there are only five classes of generating volume forms in ${\mathbb{R}}^3$ (up to relabeling and adjunction) for a given vector field $\vf$ in ${\mathbb{R}}^3$, we proceed with a classification.

Class SE+SE {#sec:SE}
-----------

We commence with case SE+SE. Its generating differential form was already discussed and extended to the $n$-dimensional case in [@lomeli09gff]. This case was also discussed in detail in [@xue2014gf], where it was associated to numerical integrators consisting of combinations of Symplectic Euler methods.

With the notation of this paper, we have $\Sigma(+)=1, \Sigma(\circ)=2,\Sigma(-)=3$ and $\sigma(-)=1,\sigma(\circ)=2,\sigma(+)=3$, see Figure \[fig:triangle\]. From $\tau=\Sigma^{-1}\circ\sigma$, we see that $\mathrm{sign}(\tau)=-1$, therefore (\[eq:-o+general\]) becomes $$\begin{aligned}
  x_3&=\partial_{x_2}\phi(x_1,x_2,X_3),\\
 \partial_{X_3}\phi(x_1,x_2,X_3)&=\partial_{x_1}\Phi (x_1,X_2,X_3),\\
 X_1&=\partial_{X_2}\Phi (x_1,X_2,X_3),
 \label{eq:case5}
 \end{aligned}$$ with twist conditions $$\frac{\partial X_1}{\partial x_1}\neq 0,~~\frac{\partial x_3}{\partial X_3}\neq 0,$$ and generating form $$\lambda=\phi (x_1,x_2,X_3) dx_1+ \Phi (x_1,X_2,X_3)  dX_3.$$ As $\lambda$ is combination of differentials $dx_1$ and $d X_3$, it is natural to consider vector fields generated by and . Setting $\Phi =x_1X_2+\Delta t F^{3}(x_1, X_2,X_3)$ and $\phi=x_2X_3+\Delta tF^{1}(x_1,x_2,X_3)$, leads to the first order volume preserving scheme $$\begin{aligned}
X_1&=x_1+\Delta t \partial_{X_2}F^{3}(x_1,X_2,X_3),\\
X_2&=x_2-\Delta t \partial_{x_1}F^{3}(x_1,X_2,X_3)+\Delta t\partial_{X_3}F^{1}(x_1,x_2,X_3),\\
X_3&=x_3-\Delta t \partial_{x_2}F^{1}(x_1,x_2,X_3),\end{aligned}$$ which is equivalent to two steps of the Symplectic Euler (SE) method to solve each of the two 2D Hamiltonian systems and .

Class DL+SE {#t2}
-----------

We have $\Sigma(+)=1, \Sigma(-)=2,\Sigma(\circ)=3$, and, as for all cases under consideration in this paper, $\sigma(-)=1,\sigma(\circ)=2,\sigma(+)=3$. We have $\mathrm{sign}(\tau)=1$ and (\[eq:-o+general\]) becomes $$\begin{aligned}
     x_3&=\partial_{x_2}\phi(x_1,x_2,X_2),\\
    \partial_{X_2}\phi(x_1,x_2,X_2)&=\partial_{x_1}\Phi (x_1,X_2,X_3),\\
    X_1&=-\partial_{X_3}\Phi (x_1,X_2,X_3),
 \label{eq:case4}
 \end{aligned}$$ with twist conditions $$\frac{\partial x_1}{\partial X_1}\neq 0,~~\frac{\partial X_2}{\partial x_3}\neq 0,$$ and generating form $$\lambda=\phi(x_1,x_2,X_2) dx_1+  \Phi (x_1,X_2,X_3) dX_2.$$ As $\lambda$ is combination of differentials $dx_1$ and $d X_2$, it is natural to consider vector fields generated by and .

Note that $\phi$ is a function of both $x_2$ and $X_2$ and that $x_3$ is determined by $x_2, X_2$. This points to an interpretation of the Hamiltonian system defined by $H(x_2,x_3)=F^1(x_1, x_2, x_3) $, where $x_2 \equiv q$, $x_3\equiv p$, and $x_1$ is treated as a constant, $$\begin{aligned}
 \dot{x}_1&=0,\\
 \dot{x}_2&=\partial_{x_3}F^{1}(x_1,x_2,x_3),\\
 \dot{x}_3&=-\partial_{x_2}F^{1}(x_1,x_2,x_3).
 \end{aligned}$$ by a *Lagrangian* formulation, with Lagrangian function $$L^1(x_1, x_2, \dot x_2) = x_3 \dot x_2 - H = x_3 \dot x_2 - F^1(x_1, x_2,x_3)$$ [@goldstein01cm]. Consider a *discrete Lagrangian* $ L^1_d= \Delta t L^1(x_1, x_2, (X_2-x_2)/\Delta t)$ [@marsden01]. With the choice $$\phi(x_1, x_2, X_2) =L_d^1 (x_1, x_2, X_2),$$ we see that the first equation of is satisfied, and, moreover, an intermediate variable for $x_3$ is obtained, $$\tilde x_3 = - \partial_{X_2} \phi (x_1, x_2, X_2).$$

For the $\Phi$ function, choose $$\Phi(x_1, X_2, X_3) = -x_1X_3 + \Delta t F^2(x_1, X_2, X_3).$$ The third equation of gives $X_1 = x_1 - \Delta t \partial_{X_3} F^2(x_1,X_2,X_3)$ and $ \Phi_{x_1} = - X_3 + \Delta t \partial_{x_1}F^2(x_1,X_2,X_3)$.

Altogether, we obtain $$\begin{aligned}
    x_3 &= \partial_{x_2}  L_d^1 (x_1, x_2, X_2) \\
    \tilde x_3 &= - \partial_{X_2}  L_d^1 (x_1, x_2, X_2) \qquad (= - \partial_{X_2} \phi) \\
    X_3 &= \tilde x_3 + \Delta t \partial_{x_1}F^2(x_1,X_2,X_3)\qquad ( \hbox{compat.\ cond.\ } \partial_{X_2} \phi = \partial_{x_1} \Phi) \\
    X_1 &= x_1 - \Delta t \partial_{X_3} F^2(x_1,X_2,X_3)\end{aligned}$$ which is a combination of a Discrete Lagrangian method (DL) for and a Symplectic Euler (SE) for . As long as the discrete Lagrangian function $L_d^1$ is a consistent approximation to the continuous one, the composed method has at least order one.

Class DL+DL {#type4}
-----------

We have $\Sigma(-) = 2$, $\Sigma(\circ) = 1$, $\Sigma(+) = 3$ and $\mathrm{sign}(\tau) = -1$ and (\[eq:-o+general\]) becomes $$\begin{aligned}
    x_3 &= \partial_{x_2} \phi(x_1, x_2, X_2), \\
    \partial_{X_2} \phi(x_1, x_2, X_2) &= \partial_{x_1} \Phi(x_1, X_2, X_1), \\
    X_3 &= \partial_{X_1} \Phi(x_1, X_2, X_1),
\label{eq:case3}\end{aligned}$$ with twist conditions $$\frac{\partial x_1}{\partial X_3} \not=0, \quad \frac{\partial X_2}{\partial x_3} \not=0,$$ and generating form $$\lambda=\phi(x_1,x_2,X_2)  dx_1+\Phi (x_1,X_1,X_2) dX_2.$$ The generating form indicates that one should look for vector fields of the form and .

Similarly to the procedure described above, we choose $\phi = L_d^1$, generating the intermediate approximation $(x_1, X_2,\tilde x_3)$ (recall that $x_1$ is kept constant). Also the system is interpreted as a Lagrangian system, with Lagrangian function $L^2(x_1, \dot x_1, X_2)=  \tilde x_3 \dot x_1 - F^2(x_1,X_2, \tilde x_3)$. Now, $X_2$ is kept constant. We set $$\Phi(x_1, X_1, X_2) = -L_d^2 (x_1,X_1,X_2),$$ where $L_d^2 (x_1,X_1,X_2)= \Delta t L^2(x_1, (X_1-x_1)/\Delta t, X_2)$ is a discrete Lagrangian approximation to $L^2$.

Altogether, we obtain $$\begin{aligned}
    x_3 &= \partial_{x_2}  L_d^1 (x_1, x_2, X_2), \\
    \tilde x_3 &= - \partial_{X_2}  L_d^1 (x_1, x_2, X_2)  \qquad (= - \partial_{X_2} \phi) \\
    \tilde x_3 &=   \partial_{x_1} L_d^2 (x_1, X_1, X_2)  \qquad (\hbox{compat.\ cond.\ } \partial_{X_2} \phi = \partial_{x_1} \Phi)\\
    X_3 &=  -\partial_{X_1} L_d^2(x_1,X_1,X_2),\end{aligned}$$ which is a combination of a Discrete Lagrangian methods (DL) for and for . As long as the discrete Lagrangian functions $L_d^1, L_d^2$ are consistent approximations to the continuous ones, the composed method has at least order one.

Special classes: $S_1$ and $S_2$ {#type5}
--------------------------------

While all the cases discussed above can be interpreted as splitting in two two-dimensional Hamiltonian systems, either solved by a symplectic method or turned into Lagrangian systems solved by a discrete Lagrangian method, there is no obvious mechanical interpretation for cases $S_1$ and $S_2$. We are not aware of any numerical method for ordinary differential equations that is naturally related to these two generating forms in the same way as all the other cases discussed in this paper. In this respect, cases $S_1$ and case $S_2$ are novel cases.

Cases $S_1$ and $S_2$ are both associated to generating forms of type $$\lambda = \phi d x_1 + \Phi d X_1,$$ which would suggest the choice of two vector fields of the form , clearly a degenerate and not particularly interesting vector field. Are there non-degenerate vector fields for which such generating form gives consistent, non trivial maps? The answer is yes: it is possible to give explicit expressions for $\phi, \Phi$ so that the generating form $\lambda = \phi d x_1 + \Phi d X_1$ of cases $S_1$ and $S_2$ gives consistent, first order, volume preserving numerical methods at least in the case when the vector field $\mathbf{a}$ in is linear: for *linear vector fields*, $\phi, \Phi$ can be taken to be quadratic functions; thus the determining conditions in are linear in the unknown variables.

Hereafter, we restrict our attention to linear divergence-free vector fields $$\begin{aligned}
\label{linear}
\dot{x}_1&=a_1(x_1,x_2,x_3)=a_{11}x_1+a_{12}x_2+a_{13}x_3,\\
\dot{x}_2&=a_2(x_1,x_2,x_3)=a_{21}x_1+a_{22}x_2+a_{23}x_3, \qquad a_{11}+a_{22}+a_{33}=0. \\
\dot{x}_3&=a_3(x_1,x_2,x_3)=a_{31}x_1+a_{32}x_2+a_{33}x_3,\end{aligned}$$

### Class $S_1$

We have that $\Sigma(-)=1, \Sigma(\circ)=2,\Sigma(+)=3$, and $\tau$ is an even permutation. The generating one-form is $$\lambda=\phi(x_1, x_2, X_1) dx_1+ \Phi (x_1, X_1, X_2) d X_1,$$ where $\phi, \Phi $ satisfy $$\begin{aligned}
x_3&=\partial_{x_2}\phi(x_1, x_2, X_1),\\
\partial_{X_1} \phi(x_1, x_2, X_1) &= \partial_{x_1} \Phi (x_1,X_1,X_2),\\
X_3& = -\partial_{X_2} \Phi (x_1,X_1,X_2).
\label{eq:S1}\end{aligned}$$ In addition to assuming linearity of the vector field, we also assume $a_{13}\neq0$ (a twist condition, implying that $x_3$ can be determined from $X_1$, given $x_1, x_2$).

Two possible choices of $(\phi, \Phi)$ yielding first order numerical methods of linear vector fields are given below.

\[th:S1\] Consider the class $S_1$ generating one-form $\lambda = \phi dx_1 + \Phi dX_1$. Let $$\begin{aligned}
    \phi(x_1,X_1,x_2)&=\frac{X_1-x_1-\Delta ta_{11}x_1-\Delta t ^2a_{21}a_{12}x_1\frac{k_2}{k_1}}{\Delta t a_{13}+\Delta t^2a_{23}a_{12}\frac{k_2}{k_1}}x_2-\frac{a_{12}}{k_1a_{13}+\Delta ta_{23}a_{12}k_2}\frac{x_2^2}{2},\\
    \Phi (x_1,X_1,X_2)&=-\frac{X_1-x_1-\Delta ta_{11}x_1}{\Delta ta_{13}}(1+\Delta ta_{33})X_2+X_2^2\big(-  \frac{\Delta ta_{32}}{2}+\frac{a_{12}}{2a_{13}}(1+\Delta ta_{33})\big)\\
&-\Delta t a_{31}X_1X_2-\frac{2X_1x_1-x_1^2(1+\Delta ta_{11})\Delta ta_{23}k_2+\Delta t^2a_{13}a_{21}k_2x_1^2}{2k_3},
\label{eq:S1_quispel}\end{aligned}$$ where $k_1=1+\Delta t^2a_{11}a_{33}-\Delta ta_{22}$, $k_2=1+\Delta t^2\frac{a_{11}a_{33}}{1-\Delta ta_{22}}$, $k_3=\Delta ta_{13}(\Delta ta_{13}+\frac{k_2}{k_1}\Delta ta_{12}a_{23})$, and $$\begin{aligned}
\phi(x_1,X_1,x_2)&=\frac{X_1-(1+\Delta ta_{11})x_1}{\Delta a_{13}}x_2-\frac{1}{2}\frac{a_{12}}{a_{13}}x_2^2,\\
\Phi (x_1,X_1,X_2)&=-\frac{1}{l_1}\big((1+\Delta ta_{33})(\frac{X_1-(1+\Delta ta_{11})x_1}{\Delta ta_{33}}X_2-\frac{1}{2}\frac{a_{12}}{a_{13}}X_2^2)+\Delta ta_{31}x_1X_2\\
&+\frac{1}{2}\Delta ta_{32}X_2^2+\Delta t\frac{a_{12}}{a_{13}}(a_{21}x_1X_1+\frac{1}{2}a_{22}X_2^2)\big)\\
&+\Delta ta_{21}x_1^2/2+\Delta ta_{23}\frac{2X_1x_1-x_1^2(1+\Delta ta_{11})}{2\Delta ta_{13}},
\label{eq:S1_AZ}\end{aligned}$$ where $l_1=1-\Delta t\frac{a_{12}}{a_{13}}a_{23}$.

Both choices and yield first-order volume preserving integrators for the vector field (\[linear\]), provided that $a_{13}\not=0$.

The proof of the result can be found in Appendix \[app:S1\].

### Class $S_2$

This case corresponds to $\Sigma(-)=1, \Sigma(\circ)=3,\Sigma(+)=2$, with $\tau$ an odd permutation. The generating one-form is $$\lambda=\phi(x_1, x_2, X_1) dx_1+ \Phi (x_1, X_1, X_3) d X_1,$$ where $\phi, \Phi $ satisfy $$\begin{aligned}
x_3&=\partial_{x_2}\phi(x_1, x_2, X_1),\\
\partial_{X_1} \phi(x_1, x_2, X_1) &= \partial_{x_1} \Phi (x_1,X_1,X_3),\\
X_2& = \partial_{X_3} \Phi (x_1,X_1,X_3).
\label{eq:S2}\end{aligned}$$ As for the $S_1$ case, we also assume $a_{12}\neq0$. Below we give the explicit expression of a choice $(\phi, \Phi)$ yielding first order numerical methods of linear vector fields.

\[th:S2\] Consider the class $S_2$ generating one-form $\lambda = \phi dx_1 + \Phi dX_1$. Let$$\begin{aligned}
\label{eq:S2_quispel}
    \phi(x_1,X_1,x_2)&=\frac{\big(m_1(X_1-x_1-\Delta ta_{11}x_1)-\Delta t^2a_{31}a_{13}m_2 x_1\big)x_2}{\Delta ta_{13}}   - \frac{ m_1a_{12}x_2^2}{2a_{13}}-\frac{\Delta ta_{32}m_2 x_2^2}2   \\
    \Phi(x_1,X_1,X_3)&= \frac{(1+\Delta ta_{22})\big(X_1-(1+\Delta ta_{11})x_1\big)X_3}{\Delta ta_{12}}-\frac{a_{13}(1+\Delta ta_{22})X_3^2}{2a_{12}}\\
& \quad \mbox{} +\Delta ta_{21}X_1X_3+\frac{\Delta ta_{23}}2 X_3^2+ \frac{m_1(2X_1x_1-\big(1+\Delta ta_{11})x_1^2\big)}{2\Delta t^2 a_{12}a_{13}},\end{aligned}$$ where $m_1=1-\Delta t a_{33}+\Delta t^2a_{11}a_{22}$ and $m_2=1+\Delta t^2a_{11}a_{22}/(1-\Delta ta_{33})$. The choice yields first-order volume preserving integrators for the vector field (\[linear\]), provided that $a_{12}\not=0$.

The proof of this result is similar to that for . For completeness, it can be found in Appendix \[app:S2\].

Conclusions and remarks
=======================

In this paper, we have studied the thirty-six generating one-forms for volume preserving mappings in ${\mathbb{R}}^3$. This is the first $n$-dimensional case for which there is a difference between area preservation, well understood using the tools of symplectic geometry and symplectic forms, and volume preservation.

By imposing equivalence relations (equivalence under relabeling and equivalence under adjunction), we have shown that all cases can be generated by five classes of differential one forms.

We have classified these five cases in terms of known numerical methods that preserve volume and denoted them as SE+SE (already identified by Xue and Zanna in [@xue2014gf]), DL+SE, DL+DL, $S_1$ and $S_2$. Except for the special cases $S_1, S_2$, the classes can be naturally associated to the splitting of the vector field into two $2D$ Hamiltonian systems or into two Lagrangian system, solved by a symplectic method (symplectic Euler, SE) or a discrete Lagrangian approach (DL), or both.

Classes $S_1$ and $S_2$, both defined by a generating one-form of the type $\lambda = \phi dx_1 + \Phi d X_1$, are, to the best of our knowledge, novel cases. Whereas the other classes admit a natural mechanical interpretation (either Hamiltonian or Lagrangian mechanics), it is not clear whether there exists a natural mechanical interpretation for the classes $S_1$ and $S_2$. The corresponding generating forms can be used to generate well defined volume preserving maps, however, for general vector fields, these maps are hightly implicit and do not seem to lead to explicit or efficient numerical methods. For completeness, we have shown possible choices of functions $\phi, \Phi$, that yield consistent methods for linear vector fields. These two cases need a deeper understanding and will be the subject of future research.

Acknowledgement {#acknowledgement .unnumbered}
===============

The work has been supported by NFR grant no. 191178/V30, under the project Geometric Numeric Integration in Applications, by the SpadeACE Project and by the J.C. Kempe memorial fund (grant no. SMK-1238).

Appendix A {#app:S1}
==========

The two choices - correspond to two different techniques to determine $\phi, \Phi$.

The first one, , is based on the *correction method* by Quispel [@GRW199526]: we determine three maps $f_1, f_2, f_3$ that give a volume preserving transformation. Thereafter we invert $f_1$ and use integration, differentiation and some other algebraic manipulations to obtain suitable $\phi, \Phi$.

The second choice, , is based on the following idea: choose a consistent method for $X_1$, depending on $x_3$, say for instance Forward Euler. Because of the linearity of the vector field, this always determines $x_3 = \phi_{x_2}(x_1,X_1,x_2)$, hence $\phi$, up to a function depending only on $x_1,X_1$. Next, we think of $x_3$ as a function of $x_2$. Now, use a consistent map to obtain $X_3$, as function of $X_1, x_1, X_2$ and $x_3(x_1, X_1,s)$, where the occurrences of $x_2$ are replaced by $X_2$). Upon the replacement $x_2 \to X_2$ in $x_3$, the map is not necessarily consistent any longer, and some adjustments must be made, also to ensure consistency for $x_2$.

*\[Prop. \[th:S1\]\]* We commence with the case. Consider the implicit map $$\begin{aligned}
\label{ff}
X_1&=f_1(x_1,X_2,x_3),\\
x_2&=f_2(x_1,X_2,x_3),\\
X_3&=f_3(X_1,X_2,x_3),\end{aligned}$$ and Quispel’s correction method [@GRW199526], to obtain $$\begin{aligned}
\label{t5vp}
X_1&=x_1+\Delta t a_1(x_1,X_2,x_3),\\
X_2&=x_2+\Delta t a_2(x_1,X_2,x_3)- f_{correct}(x_1,X_2,x_3),\\
X_3&=x_3+\Delta t a_3(X_1,X_2,x_3),\end{aligned}$$ where $f_{correct}$ is determined to obtain a volume preserving scheme, $$\begin{aligned}
f_{correct}(x_1,X_2,x_3)&=\int^{X_2}_{const}\Delta t\frac{\partial a_3}{\partial x_3}(x_1+\Delta t a_1(x_1,X_2,x_3),X_2,x_3)-\Delta t\frac{\partial a_3}{\partial x_3}(x_1,X_2,x_3)\\
&+\Delta t^2\frac{\partial a_1}{\partial x_1}(x_1,X_2,x_3)\frac{\partial a_3}{\partial x_3}(x_1+\Delta t a_1(x_1,X_2,x_3),X_2,x_3)dX_2\\
&=\Delta t^2a_{11}a_{33}(X_2-const).\end{aligned}$$ The integration constant should satisfy $const=\Delta t a_2(x_1,const,x_3)$[^1], that is, $$const=\frac{\Delta t a_{21}x_1+\Delta ta_{23}x_3}{1-\Delta ta_{22}}.$$ First of all, we calculate $X_2$ from the second equation of (\[t5vp\]) which has the form, $$\label{x2}
X_2=\frac{x_2+(\Delta ta_{21}x_1+\Delta ta_{23}x_3)k_2}{k_1},$$ where $k_1=1+\Delta t^2a_{11}a_{33}-\Delta ta_{22}$ and $k_2=1+\Delta t^2\frac{a_{11}a_{33}}{1-\Delta ta_{22}}$. Substituting the above equation into the first equation of (\[t5vp\]), we can find $x_3$ in terms of variables $x_1, x_2, X_1$, denoted by $\tilde{x}_3(x_1,X_1,x_2)$. Substituting $\tilde{x}_3$ into the first equation of (\[eq:S1\]) and integrating both sides, we obtain $\phi$ $$\phi=\frac{X_1-x_1-\Delta ta_{11}x_1-\Delta t ^2a_{21}a_{12}x_1\frac{k_2}{k_1}}{\Delta t a_{13}+\Delta t^2a_{23}a_{12}\frac{k_2}{k_1}}x_2-\frac{a_{12}}{k_1a_{13}+\Delta ta_{23}a_{12}k_2}\frac{x_2^2}{2}+\tilde{A}(x_1,X_1).$$ From the first equation of (\[t5vp\]), we see that $x_3$ depends on the variables $x_1, X_1,X_2$. We then can solve $x_3$ in terms of $x_1, X_1,X_2$ from that equation since we assume $a_{13}\neq0$, and denote by $\hat{x}_3(x_1,x_2,X_1)$. We substitute $\hat{x}_3$ into the third equation of (\[t5vp\]). From the third equations of (\[eq:S1\]) and from (\[t5vp\]), we know that $$X_3(x_1,X_1,X_2)=\hat{x}_3(x_1,X_1,X_2)+\Delta t \mathbf{a}_3(X_1,X_2,\hat{x}_3(x_1,X_1,X_2))=-\partial_{X_2}\Phi .$$ From the above equation, we can integrate both sides with respect to $X_2$ and obtain $$\begin{aligned}
\Phi &=-\frac{X_1-x_1-\Delta ta_{11}x_1}{\Delta ta_{13}}(1+\Delta ta_{33})X_2+X_2^2(-\frac{\Delta ta_{32}}{2}+\frac{a_{12}}{2a_{13}}(1+\Delta ta_{33}))\\
&-\Delta t a_{31}X_1X_2+\tilde{C}(x_1,X_1).\end{aligned}$$ Without loss of generality, we set $\tilde{A}=0$. Using the second condition of (\[eq:S1\]) , we obtain $$\tilde{C}_{x_1}=\frac{x_2}{\Delta a_{13}+\Delta ta_{12}a_{23}\frac{k_2}{k_1}}-\frac{(1+\Delta ta_{11})(1+\Delta ta_{33})}{ha_{13}}X_2,$$ From the first equation in (\[t5vp\]), we have $$X_2=\frac{X_1-x_1-\Delta ta_{11}x_1-\Delta ta_{13}x_3}{\Delta ta_{12}},$$ By noticing the relation in (\[x2\]), we obtain $$\tilde{C}_{x_1}=-\frac{X_1-x_1(1+\Delta ta_{11})\Delta ta_{23}k_2+\Delta t^2a_{13}a_{21}k_2x_1}{k_3},$$ where $k_3=\Delta ta_{13}(\Delta ta_{13}+\frac{k_2}{k_1}\Delta ta_{12}a_{23})$. From it is not difficult to see that (\[eq:S1\_quispel\]) gives a first order volume preserving method.

Next, we consider the case. Using the forward Euler method to solve the first equation of (\[linear\]), that is $$X_1=(1+\Delta ta_{11})x_1+\Delta ta_{12}x_2+\Delta ta_{13}x_3.$$ Assuming that $a_{13}\neq0$ and solving for $x_3$, we have $$\label{x3tilde}
x_3=\frac{X_1-(1+\Delta ta_{11})x_1-\Delta ta_{12}x_2}{\Delta ta_{13}}.$$ Integrating both sides, we obtain $$\phi(x_1,X_1,x_2)=\frac{X_1-(1+\Delta ta_{11})x_1}{\Delta a_{13}}x_2-\frac{1}{2}\frac{a_{12}}{a_{13}}x_2^2+\tilde{A}(x_1,X_1),$$ where $\tilde{A}(x_1, X_1)$ is a function to be determined.

Using Euler method to solve the third equation of (\[linear\]), where $x_3$ in (\[x3tilde\]) is replaced by $x_3(x_1,X_1,X_2)$ and using $X_1=x_1+\Delta t a_1$ and $X_2=x_2+\Delta t a_2$, we obtain $$\begin{aligned}
    X_3&=(1+\Delta ta_{33})\frac{X_1-(1+\Delta ta_{11})x_1-\Delta ta_{12}X_2}{\Delta ta_{13}}+\Delta    ta_{31}x_1+a_{32}X_2\\
    &=(1+\Delta ta_{33})x_3+\Delta ta_{31}x_1+\Delta ta_{32}X_2-\Delta t\frac{a_{12}}{a_{13}}a_2+O(\Delta t^2).\end{aligned}$$ There is a problem with the term $a_2$ (as in (\[linear\])) in the above equation for consideration of consistency for $X_3$. So, we take $a_2=a_{21}x_1+a_{22}X_2+a_{23}X_3$ and substitute back into the above equation. Then, we obtain $$(1-\Delta t\frac{a_{12}}{a_{13}}a_{23})X_3=(1+\Delta ta_{33})x_3+\Delta ta_{31}x_1+\Delta ta_{32}X_2-\Delta t   \frac{a_{12}}{a_{13}}(a_{21}x_1+a_{22}X_2)$$ Solving $X_3$ from the above equation, substituting it back into the third equation of (\[eq:S1\]) and integrating both sides, we have $$\begin{aligned}
    \Phi &=-\frac{1}{l_1}\big((1+\Delta ta_{33})(\frac{X_1-(1+\Delta ta_{11})x_1}{\Delta ta_{33}}X_2-\frac{1}{2}\frac{a_{12}}{a_{13}}X_2^2)+\Delta ta_{31}x_1X_2\\
    &+\frac{1}{2}\Delta ta_{32}X_2^2+\Delta t\frac{a_{12}}{a_{13}}(a_{21}x_1X_1+\frac{1}{2}a_{22}X_2^2)\big)+\tilde{C}(x_1,X_1),\end{aligned}$$ where $l_1=1-\Delta t\frac{a_{12}}{a_{13}}a_{23}$ and $\tilde C$ is to be determined. We have two functions $\tilde A, \tilde C$ to be determined and one equation (compatibility condition). The two functions are not independent, hence we set $\tilde{A}=0$. Using the compatibility condition $\partial_{X_1} \phi= \partial_{x_1} \Phi $, we obtain $$\frac{x_2}{\Delta ta_{13}}=\frac{1}{l_1}(1+\Delta ta_{33})(1+\Delta ta_{11})X_2+\tilde{C}_{x_1}.$$ From the divergence-free condition $a_{22}=-(a_{11}+a_{33})$ we have $$X_2=x_2(1+\Delta ta_{22}+O(\Delta t^2))l_1+\tilde{C}_{x_1}.$$ Hence $X_2$ will be consistent provided that $$\tilde{C}=\Delta t a_{21}\frac{x_1^2}{2}+\Delta t a_{23}\frac{2X_1x_1-x_1^2(1+\Delta ta_{11})}{2\Delta ta_{13}}.$$ From the proof, we see that (\[eq:S1\_AZ\]) gives a first order volume preserving method.

Appendix B {#app:S2}
==========

*\[Prop. \[th:S2\]\]* The method generated by and is constructed in the same way as . The implicit map reads $$\begin{aligned}
\label{ffs2}
X_1&=f_1(x_1,x_2,X_3),\\
X_2&=f_2(X_1,X_2,x_3),\\
x_3&=f_3(x_1,x_2,X_3),\end{aligned}$$ and the corrections method is $$\begin{aligned}
\label{t5vps2}
X_1&=x_1+\Delta t a_1(x_1,x_2,X_3),\\
X_2&=x_2+\Delta t a_2(X_1,x_2,X_3),\\
X_3&=x_3+\Delta t a_3(x_1,x_2,X_3)- f_{correct}(x_1,x_2,X_3),\\\end{aligned}$$ where $$\begin{aligned}
    f_{correct}(x_1,x_2,X_3)&=\int^{X_3}_{const}\Delta t\frac{\partial a_2}{\partial x_2}(x_1+\Delta t  a_1(x_1,x_2,X_3),x_2,X_3)-\Delta t\frac{\partial a_2}{\partial x_2}(x_1,x_2,X_3)\\
    &+\Delta t^2\frac{\partial a_1}{\partial x_1}(x_1,x_2,X_3)\frac{\partial a_2}{\partial x_2}(x_1+\Delta t    a_1(x_1,x_2,X_3),x_2,X_3)dX_3\\
    &=\Delta t^2a_{11}a_{22}(X_3-const).\end{aligned}$$ The integration constant should satisfy $const=\Delta t a_3(x_1,x_2,const1)$, that is, $$const1=\frac{\Delta t a_{31}x_1+\Delta ta_{32}x_2}{1-\Delta ta_{33}}.$$ The rest of the proof is similar to that of Prop. \[th:S1\].

[^1]: Due to consideration of consistency for $X_2$, see more details in [@GRW199526].
