$\overrightarrow{\par}$

# Numerical Methods in Poisson Geometry and their Application to Mechanics

Oscar Cosserat  O. Cosserat: LaSIE – CNRS & La Rochelle University, Av. Michel
Crépeau, 17042 La Rochelle Cedex 1, France<EMAIL_ADDRESS>,
Camille Laurent-Gengoux  C. Laurent-Gengoux: Institut Elie Cartan de Lorraine
(IECL), UMR 7502 – 3 rue Augustin Fresnel, 57000 Technopôle Metz, France
<EMAIL_ADDRESS>and  Vladimir Salnikov  V. Salnikov:
LaSIE – CNRS & La Rochelle University, Av. Michel Crépeau, 17042 La Rochelle
Cedex 1, France<EMAIL_ADDRESS>
###### Abstract.

We recall the question of geometric integrators in the context of Poisson
geometry, and explain their construction. These Poisson integrators are tested
in some mechanical examples. Their properties are illustrated numerically and
they are compared to traditional methods.

###### Contents

  1. 1 An introductory example
  2. 2 Geometric integrators for Poisson Hamiltonian systems.
    1. 2.1 Poisson structures: definition
    2. 2.2 The underlying geometry of a Poisson structure
    3. 2.3 Examples of Hamiltonian equations – first candidates for integrators.
    4. 2.4 Poisson Hamiltonian integrators.
  3. 3 Bi-realisations for Poisson manifolds
    1. 3.1 Bi-realisation I: Definition and existence.
    2. 3.2 Bi-realisations II: Explicit constructions.
    3. 3.3 Examples
  4. 4 Explicit construction of Poisson Hamiltonian integrators
  5. 5 Numerical tests
    1. 5.1 The Rigid Body
    2. 5.2 The Lotka-Volterra System

## Introduction

Poisson geometry allows to describe a large class of conservative systems in
mechanics, both for discrete and continuous media. Those may be obtained as a
result of a reduction procedure or as an ad hoc model for evolution of some
natural systems. Just to cite a few non-standard examples: chemistry (polymer
dynamics, [8]), plasticity (elastoplasticity, [18]), population dynamics
([16]), liquid crystals theory ([9]), thermodynamics (GENERIC formalism,
[12]), control theory (active and kinematic constraints, [21]). Moreover, the
geometry of the Poisson structure matters to express symmetries, conservation
laws and qualitative behaviour of dynamical systems.

On top of the purely mathematical significance of Poisson structures, for
Hamiltonian differential equations, they provide the appropriate geometry to
be preserved in numerical simulation thus potentially resulting in a very
broad class of geometric integrators. The state of the art in this context
concerns mostly _symplectic integrators_ , known for decades now ([26]), they
correspond to symplectic structures, i.e. non-degenerate Poisson for systems
defined on a canonical phase space, geometrically meaning on a cotangent
bundle. To the best of our knowledge there are very few works treating more
general Poisson structures. Most of them rely on Weinstein’s splitting theorem
(see for example [30]), which states that a Poisson manifold can be foliated
by symplectic leaves, the natural idea that emerges is: “restrict to the leaf
containing the trajectory, apply symplectic integration on it”.

There are two conceptual difficulties with this “naive” approach to construct
_Poisson integrators_ : the mentioned foliation is very rearly explicitly
known, and in the generic situation it is singular, meaning that the
successfull applications of the above strategy are rather exceptional
particular cases. A way out was proposed in [1], it is based on integration of
Poisson manifolds to (local) symplectic groupoids, geometrically this
procedure can be viewed as a kind of desingularisation.

The goal of the current paper is to explain the subtleties occurring in this
integration procedure while making it explicit and also present the features
of resulting Poisson integrators. We will not go into mathematical proofs (a
motivated reader is invited to consult [1] for details), but rather focus on
the precise behaviour of the constructed numerical methods. We will however
provide all the necessary building blocks to make the paper self-consistent.

We can address one important detail already here: the mere definition of what
is an appropriate notion of structure preserving numerical method in the
context of Poisson geometry. For the symplectic case the usual strategy is to
say that the the preservation of the symplectic form guarantees the energy
conservation. In fact it is a bit subtler than that: a discrete symplectic
flow indeed preserves the level surface of some Hamiltonian (suppose there is
no topological obstruction for its global existence), but not necessarily the
same that corresponds to the energy of the system and was defining the
evolution of it. The difference however can be estimated: it will be
oscillating around zero with an amplitude that can be made small for
exponentially long time. For the Poisson case this subtlety is even more
pronounced: the preservation of the Poisson structure alone does not guarantee
the energy conservation, roughly speaking this is due to the existence of the
“degenerate” directions.

We thus introduce a stronger notion of a _Poisson Hamiltonian integrator_
(PHI) – a numerical method for which discretising a trajectory of a
Hamiltonian system preserves the Poisson structure and coincides exactly with
the trajectory of some time-dependent Hamiltonian, which equals to the initial
one at the order of the method. This latter condition seems to be very strong,
but in fact it is not: the Poisson Hamiltonian integrator requires the
existence of such a time-dependent Hamiltonian but not constructively, in
fact, only a theoretical existence is enough: we do not need, in general, to
compute it. This unspecified existence guarantees (see again [1]) that the
discretisation preserves the geometric invariants of the phase space
(symplectic leaves, leafwise symplectic structures) together with the
qualitative physical properties of the system – we will illustrate this in
different situations.

The paper is organized as follows: We start with an illustration of our PHI
technique for a simple Lotka-Volterra system – we observe a better behaviour
of a numerical solution in comparison with the standard Runge-Kutta method. In
section 2 we recall the minimal working knowledge from Poisson geometry to
formulate precisely the definition and give the strategy of construction of
the Poisson Hamiltonian integrators. Then in section 3 we introduce more
advanced geometric notions that are needed for the construction which is made
explicit in section 4. The last section 5 is devoted to numerical results of
comparison of the constructed Poisson Hamiltonian integrators with standard
methods for various typical situations.

##  1\. An introductory example

Let us look at a particular case of Lotka-Volterra type equations

$\begin{array}[]{lll}\dot{x}_{1}&=x_{1}(x_{2}+x_{3})\\\
\dot{x}_{2}&=x_{2}(-x_{1}+x_{3})\\\
\dot{x}_{3}&=-x_{3}(x_{1}+x_{2}).\end{array}$

This system of differential equations appears in [27] (page 97, equation (16))
as a model in population dynamics, and similar systems have been extensively
studied since then. For this particular one, an explicit solution was computed
in [16]: it allows to compare, for any initial point and desired time,
numerical simulations with the exact solution
${x}^{\hbox{\tiny{exact}}}(t)=(x_{1}(t),x_{2}(t),x_{3}(t))$.
Some integral curves ${x}^{\hbox{\tiny{exact}}}(t)$ go to infinity, exploding
exponentially fast while approaching some specific time. For instance, using
the exact formulas given in [16], one observes that for the initial value

$\begin{pmatrix}x_{1}(0)\\\ x_{2}(0)\\\
x_{3}(0)\end{pmatrix}=\begin{pmatrix}-3\\\ 5\\\ 10^{-3}\end{pmatrix},$

the trajectory ${x}^{\hbox{\tiny{exact}}}(t)$ starts exploding around
$T_{sing}\simeq 0.23$. We use this singularity to test two numerical methods
of different nature.

The first one is the standard explicit 2nd order Runge–Kutta method (RK-2);
the second one is a numerical scheme (PHI-1) where
${{x}}_{n+1}^{\hbox{\tiny{PHI}}}$ is implicitely computed from
${{x}}_{n}^{\hbox{\tiny{PHI}}}$ by the following two step procedure:

$\left\\{\begin{array}[]{lll}x_{1,n+1}&=e^{(\frac{\Delta
t}{2}({y_{2,n}}+{y_{3,n}}))}{y_{1,n}}\\\ x_{2,n+1}&=e^{(\frac{\Delta
t}{2}(-{y_{1,n}}+{y_{3,n}}))}{y_{2,n}}\\\ x_{3,n+1}&=e^{(\frac{-\Delta
t}{2}({y_{1,n}}+{y_{2,n}}))}{y_{3,n}}\end{array}\right.$

with the intermediate point $\mathbf{y}_{n}=(y_{1,n},y_{2,n},y_{3,n})$ subject
to the relation:

$\left\\{\begin{array}[]{lll}e^{(\frac{-\Delta
t}{2}(y_{2,n}+{y_{3,n}}))}{y_{1,n}}&=x_{1,n}\\\ e^{(\frac{\Delta
t}{2}({y_{1,n}}-{y_{3,n}}))}{y_{2,n}}&=x_{2,n}\\\ e^{(\frac{\Delta
t}{2}({y_{1,n}}+{y_{2,n}}))}{y_{3,n}}&=x_{3,n}\end{array}\right.$

This numerical scheme is of order 1 only. Nevertheless, it performs much
better than the order two Runge-Kutta method near the singularity $T_{sing}$
described above (see Fig. 1):

Figure 1. Comparison of the error for RK-2 and PHI-1

We observe that the PHI-1 method approximates the solution much better than
the RK-2. In fact comparing the values of the variables, we see that the RK-2
misses the singularity completely, in a sense that it goes off the exact
solution much earlier than it may tend to infinity, so using it alone one
would not even notice that the solution is singular; while the PHI-1 method
pushes the solution up to the last step before hitting the singularity, where
it goes to what one can call “numerical infinity”.

It is actually quite unexpected that an order $1$ integrator behaves better
than an order $2$ integrator. The same situation happens when one compares
PHI-1 one with a 4th order Runge-Kutta scheme for bigger values of $T.$ The
explanation is that as often in such situations when such a behavior happens,
we are comparing an integrator that does not preserve the underlying geometric
structure of the differential equations with an integrator that does. We will
see below that the second method is a _Poisson Hamiltonian integrator_ , so we
will understand what is preserved exactly and why this guarantees an
appropriate trajectory. The aim of the following sections is to explain what
are such integrators, which properties they satisfy, and how one can construct
and implement them in a wide class of examples.

##  2\. Geometric integrators for Poisson Hamiltonian systems.

###  2.1. Poisson structures: definition

Let us briefly recall111A reader familiar with Poisson geometry may safely
skip this section up to 2.3. what are Poisson bivector fields and Hamiltonian
differential equations, and why they matter. In mechanics, quite some
differential equations governing a motion $x(t)$ valued in an open subset
$\mathcal{U}\subset\mathbb{R}^{n}$ take the form:

$\left\\{\begin{array}[]{rcl}\dot{x}_{1}(t)&=&\sum\limits_{j=1}^{n}\pi_{1j}(x(t))\,\frac{\displaystyle\partial
H(x(t))}{\displaystyle\partial x_{j}}\\\ &\vdots&\\\
\dot{x}_{n}(t)&=&\sum\limits_{j=1}^{n}\pi_{nj}(x(t))\,\frac{\displaystyle\partial
H(x(t))}{\displaystyle\partial x_{j}}\end{array}\right.$ (1)

where $H\colon\mathcal{U}\longrightarrow\mathbb{R}$ is a smooth function, that
it is customary to call _Hamiltonian_ in this context, and
$(\pi_{ij}),\;i,j\in\\{1,\dots,n\\}$ is a family of smooth functions which
satisfy

$\pi_{ji}=-\pi_{ij}\;\text{ and
}\;\sum_{a=1}^{n}\frac{\partial\pi_{ij}}{\partial
x_{a}}\,\pi_{ak}+\circlearrowleft=0,$ (2)

where all indices $i,j,k\in\\{1,\dots,n\\}$ and $\circlearrowleft$ stands for
their cyclic permutations..

Let us first explain these conditions: there is a theorem (see chapter 1 of
[23]) claiming that a family $(\pi_{ij}),i,j\in\\{1,\dots,n\\}$ of smooth
functions on $\mathcal{U}\subset\mathbb{R}^{n}$ satisfy (2) if and only if the
bilinear map:

$\begin{array}[]{rcll}\\{\cdot,\cdot\\}\colon&\mathcal{C}^{\tiny\infty}(\mathcal{U})\times\mathcal{C}^{\tiny\infty}(\mathcal{U})&\longrightarrow&\mathcal{C}^{\tiny\infty}(\mathcal{U})\\\
&f,g&\mapsto&\sum\limits_{i,j=1}^{n}\frac{\displaystyle 1}{\displaystyle
2}\pi_{ij}(x)\left(\frac{\displaystyle\partial f}{\displaystyle\partial
x_{i}}\frac{\displaystyle\partial g}{\displaystyle\partial
x_{j}}-\frac{\displaystyle\partial g}{\displaystyle\partial
x_{i}}\frac{\displaystyle\partial f}{\displaystyle\partial x_{j}}\right)\\\
\end{array}$

is a Lie bracket, i.e. is anti-symmetric and satisfies the Jacobi identity:

$\\{f,g\\}=-\\{g,f\\}\hbox{ and
}\\{\\{f,g\\},h\\}+\\{\\{g,h\\},f\\}+\\{\\{h,f\\},g\\}=0$ (3)

for all $f,g,h\in\mathcal{C}^{\tiny\infty}(M)$. Equivalently, the functions
$\pi_{ij}(x)$ above can be considered as being a _tensor_ (i.e. a map from
$\mathcal{U}\subset\mathbb{R}^{n}$ to ${\mathfrak{gl}}_{n}(\mathbb{R})$,
viewed as $n\times n$ matrices depending on $x$:

$\begin{array}[]{rcll}\pi\colon&\mathcal{U}&\longrightarrow&{\mathfrak{gl}}_{n}(\mathbb{R})\\\
&x&\mapsto&(\pi_{ij}(x)_{i,j=1}^{n})\end{array}$

######  Definition 2.1.

A _Poisson structure_ on an open subset $\mathcal{U}\subset\mathbb{R}^{n}$ is
a tensor $\pi$ on $\mathcal{U}$ such that the bilinear map
$\mathcal{C}^{\tiny\infty}(\mathcal{U})\times\mathcal{C}^{\tiny\infty}(\mathcal{U})\to\mathcal{C}^{\tiny\infty}(\mathcal{U})$
defined by

$\\{f,g\\}(x)={}^{t}\nabla_{x}f.\pi(x).\nabla_{x}g$

satisfies the skew-symmetry and Jacobi identities (3). Equivalently, Poisson
structures can be seen as tensors whose coefficients
$(\pi_{ij}(x)),i,j\in\\{1,\dots,n\\}$ satisfy Equations (2).

The long history behind this notion comes with a vocabulary, which is
sometimes confusing: it is customary to call the bilinear map
$\\{\cdot,\cdot\\}$ the _Poisson bracket_. Also, functions in
$\mathcal{C}^{\tiny\infty}(\mathcal{U})$ are – depending on the context –
sometimes all called _Hamiltonian functions_ or simply _Hamiltonians_. Last,
the easy to check relation:

$\\{f,gh\\}=\\{f,g\\}h+g\\{f,h\\}$

is called _Leibniz identity_. For examples of Poisson structures that appear
in mechanics, see Section 3.3.

We saw that a Poisson structure associates to _two_ Hamiltonian function
$f,g\in\mathcal{C}^{\tiny\infty}(\mathcal{U})$ another Hamiltonian function
$\\{f,g\\}$. But it also allows to associate to _one_ Hamiltonian function
$H\in\mathcal{C}^{\tiny\infty}(\mathcal{U})$ a first order autonomous
differential equation, as in (1). More abstractly, (1) means that to a
Hamiltonian function $H$ we associate the (vector) differential equation:

$\dot{x}(t)=\pi(x(t))\cdot\nabla_{x(t)}H.$ (4)

We say that a first order autonomous differential equation of the form (4)
above is a _Hamiltonian differential equation for $(\pi,H)$_.

######  Remark 2.2.

In differential geometry, a first order differential equation on an open
subset $U$ of $\mathbb{R}^{n}$ with an independent parameter having the
meaning of time

$\dot{x}=F(x)$ (5)

is generally referred to as a vector field on $U$. Moreover, rather than
considering only open subsets of $\mathbb{R}^{n}$, the phase space is often
assumed to be a differential manifold.

###  2.2. The underlying geometry of a Poisson structure

It is natural to ask why it matters that behind an autonomous first order
differential equation, there is a Poisson structure and a Hamiltonian
function. What does one gain by knowing that a given differential equation is
Hamiltonian of $(\pi,H)$? The classical answer is that many “quantities”
related to the $\pi$ or $H$ are preserved under the flow of the differential
equation.

More precisely, any solution $x(t)$ (often called “integral curve” in
mathematics) of a differential equation Hamiltonian for $(\pi,H)$ preserves
several functions:

  1. (1)

$H$ is a constant of motion, i.e. $H(x(t))=H(x(0))$. In words, it means that
the flow of a Hamiltonian differential equation for $(\pi,H)$ preserves the
level sets of the Hamiltonian function $H$.

  2. (2)

More generally, any function $G$ such that $\\{G,H\\}=0$ is a constant of
motion.

  3. (3)

Even more generally, for any function
$G\in\mathcal{C}^{\tiny\infty}(\mathcal{U})$, one hav

$\frac{\mathrm{d}G(x(t))}{\mathrm{d}t}=\\{G,H\\}(x(t)).$

Now, let us recall some properties of the time $t$ flow of a differential
equation which is Hamiltonian for $(\pi,H)$, i.e. the map $\phi_{t}\colon
x(0)\mapsto x(t)$, which is well-defined in a neighborhood of any
$m\in\mathcal{U}$ for $t$ small enough:

  1. (1)

$\phi_{t}$ preserves $\pi$:
$\pi(\phi_{t}(x))={}^{t}\nabla_{x}\phi_{t}(x).\pi(x).\nabla_{x}\phi_{t}(x)$.
In words, it means that the flow of a Hamiltonian differential equation for
$(\pi,H)$ preserves the Poisson structure $\pi$.

  2. (2)

The previous condition can also be stated as meaning that the pull-back map
$f\mapsto\phi_{t}^{*}f$, i.e. the map assigning to a smooth Hamiltonian
function $f$ the Hamiltonian function $f\circ\phi_{t}$, is a Lie algebra
morphism, i.e.

$\phi_{t}^{*}\\{f,g\\}=\\{\phi_{t}^{*}f,\phi_{t}^{*}g\\}$

for all functions $f,g$.

  3. (3)

Item 1 above means in particular that the geometry of $\pi$ is preserved. For
instance if at the initial point $x(0)$, the matrix $\pi(x(0))$ has some given
rank, it has this same rank at every point along the integral curve $x(t)$.
Below, we will give a much stronger statement, using the notion of symplectic
leaves.

A _symplectic singular foliation_ on $\mathcal{U}\subset\mathbb{R}^{n}$ is a
partition $\mathcal{U}=\sqcup_{c\in I}\mathcal{F}_{c}$ by submanifolds, such
that each $\mathcal{F}_{c}$ is equipped with a symplectic structure
$\omega_{c})$. The pair $(\mathcal{F}_{c},\omega_{c})$ is called a _symplectic
leaf_.

######  Theorem 2.3 ([23, 30]).

Any Poisson structure on $\mathcal{U}\subset\mathbb{R}^{n}$ induces a natural
foliation by symplectic leaves characterized by the following two properties:

  1. (1)

two points belong to the same symplectic leaf if and only if they can be
connected by a sequence of Hamiltonian trajectories, i.e. by integrating the
equation 4 for some choice of Hamiltonian functions.

  2. (2)

for every $c\in I$, the inclusion $i\colon\mathcal{F}_{c}\xhookrightarrow{}M$
is a Poisson morphism.

In addition, the tangent space of the symplectic leaf $\mathcal{F}_{c}$ at a
point $m$ coincides with the image of the matrix $\pi(m)$.

Figure 2. A foliation of $\mathbb{R}^{3}$ by concentric spheres and the
origin.

As observed in some of the following examples, the foliation is generically
singular. Two neighbouring leaves do not necessarily have the same dimension
and can differ from a topological point of view. Therefore, its study is an
active field of research and motivates one of the long term applications of
the numerical tools we present here.

The last reason explaining the importance of knowing that a differential
equation is Hamiltonian for $(\pi,H)$ is the following: a solution $x(t)$ of a
differential equation Hamiltonian for $(\pi,H)$ can not “jump” from one
symplectic leaf to another. That is if $x(0)$ belongs to a leaf
$\mathcal{F}_{c}$, then the solution $x(t)$ belongs to the same leaf for every
$t$.

######  Remark 2.4.

This last “constraint” is maybe less studied for numerical methods than the
previous ones, because when the Poisson is symplectic, it is not a constraint
at all: the foliation contains only one leaf being the whole space. But for
non-symplectic Poisson structures this is a very important feature to take
under account.

In conclusion, for any Poisson structure $\pi$:

  1. $\diamondsuit$

for any first order autonomous differential equation which is Hamiltonian for
$(\pi,H)$, the Hamiltonian $H$ is constant along any integral curve;

  2. $\heartsuit$

each integral curve stays on the same symplectic leaf of the foliation defined
by $\pi$;

  3. $\spadesuit$

the flow of this differential equation preserves $\pi$, i.e. is a Poisson
morphism;

  4. $\clubsuit$

the converse is not necessarily true: preservation of $\pi$ does not guarantee
that the flow is Hamiltonian.

###  2.3. Examples of Hamiltonian equations – first candidates for
integrators.

The goal of what follows is to explain the logic behind the construction of
numerical schemes that take into account the geometric features described
above. We illustrate it on simple cases yet instructive examples.

Important examples of Poisson structures are the symplectic ones in their
canonical form, e.g. where $\pi=\begin{pmatrix}0&-I\\\ I&0\end{pmatrix}$. For
those, a wide example of symplectic integrators are already available in the
literature. One construction of such integrators uses the principle of
symplectic Runge–Kutta schemes ([29]):

$\begin{array}[]{ll}x_{n+1}&=x_{n}+\Delta t\sum\limits_{i=1}^{s}b_{i}k_{i}\\\
k_{i}&=\pi(x_{n}+\Delta t\sum\limits_{j=1}^{n}a_{ij}k_{j}).\nabla
H(x_{n}+\Delta t\sum\limits_{j=1}^{n}a_{ij}k_{j})\end{array}$ (6)

where slopes $k_{i}$ are implicitly defined and coefficients $b_{i}$ and
$a_{ij}$ are chosen such that the discrete flow preserves $\pi$. For this
precise $\pi$, any trajectory preserving it is necessary a time-dependent
Hamiltonian one, at least locally.

When $\pi$ is degenerate, the same principle can be applied ([14]) and leads
to a discrete flow that preserves the Poisson tensor. However, it would lead
to non-physical simulations, e.g. non-Hamiltonian ones. It does not guarantee
the Hamiltonian property of the discrete trajectory anymore because of the
existence of outer Poisson automorphisms, as illustrated in the following
example.

######  Example 2.5.

Consider $\mathcal{U}=\mathbb{R}^{3}$. The system of differential equations

$\left\\{\begin{array}[]{rcl}\dot{x}&=&-\frac{(x+y-z)}{8}\left((-x-y+z)^{2}+(x+y-z)^{2}\right)\\\
\dot{y}&=&\frac{(-y+z)}{4}\left((x-y+z)^{2}+(x+y-z)^{2}\right)\\\
\dot{z}&=&\frac{(x-y+z)}{8}\left((x-y+z)^{2}+(x+y-z)^{2}\right)\end{array}\right.$
(7)

is Hamiltonian with respect to the Poisson structure

$\pi(x,y,z)=\frac{(x-y+z)^{2}+(x+y-z)^{2}}{4}\begin{pmatrix}0&-1&-1\\\
1&0&-1\\\ 1&1&0\end{pmatrix}$

and the Hamiltonian $H\colon(x,y,z)\mapsto\frac{(x-y+z)^{2}+(x+y-z)^{2}}{8}.$
For any $\Delta t>0$, the system of equations:

$\left\\{\begin{array}[]{lllll}x_{n+1}=&x_{n}\cos{\left(\Delta
t\frac{x_{n}^{2}+y_{n}^{2}-2y_{n}z_{n}+z_{n}^{2}}{2}\right)}+y_{n}\sin{\left(\Delta
t\frac{x_{n}^{2}+y_{n}^{2}-2y_{n}z_{n}+z_{n}^{2}}{2}\right)}\\\
&-z_{n}\sin{\left(\Delta
t\frac{x_{n}^{2}+y_{n}^{2}-2y_{n}z_{n}+z_{n}^{2}}{2}\right)}\\\
y_{n+1}=&\frac{-x_{n}+y_{n}-z_{n}}{2}\sin{\left(\Delta
t\frac{(x_{n}-y_{n}+z_{n})^{2}+(x_{n}+y_{n}-z_{n})^{2}}{4}\right)}+\frac{-x_{n}+y_{n}+z_{n}}{2}\exp{(\Delta{t}^{k})}\\\
&+\frac{x_{n}+y_{n}-z_{n}}{2}\cos{\left(\Delta{t}\frac{(x_{n}-y_{n}+z_{n})^{2}+(x_{n}+y_{n}-z_{n})^{2}}{4}\right)}\\\
z_{n+1}=&\frac{-x_{n}+y_{n}+z_{n}}{2}\exp{(\Delta{t}^{k})}+\frac{x_{n}-y_{n}+z_{n}}{2}\cos{\left(\Delta{t}\frac{(x_{n}-y_{n}+z_{n})^{2}+(x_{n}+y_{n}-z_{n})^{2}}{4}\right)}\\\
&+\frac{x_{n}+y_{n}-z_{n}}{2}\sin{\left(\Delta{t}\frac{(x_{n}-y_{n}+z_{n})^{2}+(x_{n}+y_{n}-z_{n})^{2}}{4}\right)}\end{array}\right.$
(8)

is a discretisation of order $k$ of the differential equation (7). It is
routine to check that it is a Poisson integrator, i.e.
$(x_{n},y_{n})\mapsto(x_{n+1},y_{n+1})$ is a Poisson isomorphism. However, it
is _not_ a Poisson Hamiltonian integrator. This can be proven as follows: for
any vector field on $\mathbb{R}^{3}$ that vanishes at least quadratically at
$(0,0,0)$, the differential of its flow at $(0,0,0)$ is the identity map. In
particular, since the coefficients of the Poisson structure vanish at least
quadratically at $(0,0,0)$, so does any Hamiltonian vector field, so that any
Poisson Hamiltonian integrator should be made of a local diffeomorphism whose
differential at $(0,0,0)$ is the identity map. Since the differential at
$(0,0,0)$ of the map $(x_{n},y_{n},z_{n})\mapsto(x_{n+1},y_{n+1},z_{n+1})$ is
not equal to identity map, the latter Poisson integrator is not Hamiltonian.

(a) Flow of (7)

(b) Poisson integrator (8)

Figure 3. Discrepancy created by the Poisson integrator (8)

Figures 3.A 3.B show the difference between the actual flow of (7) and the
first iterations of the Poisson integrator (8) at order $k=2$ with initial
points $\begin{pmatrix}1\\\ 1\\\ 2\end{pmatrix}$ and timestep $\Delta
t=10^{-4}.$ The flow should be $4\pi$-periodic while an approximation of it at
order 2 destroys the topology of the curve, even while it preserves the
Poisson tensor. The geometric reason is that the Poisson integrator (8) does
not stay on a symplectic leaf of the Poisson structure, i.e. a hyperplane of
equation $\\{x-y+z=\text{constant}\\}.$

###  2.4. Poisson Hamiltonian integrators.

As in classical numerical analysis, we call _integrator_ of order $k$ for a
differential equation

$\dot{x}(t)=F(x(t))$

a family of diffeomorphisms222There is a subtle point here: we cannot assume
$\phi_{h}$ to be a well-defined diffeomorphism from $\mathcal{U}$ to
$\mathcal{U}$ for all $h$ small enough, but we can assume that for all $h$
small enough, there is an open subset $\mathcal{U}_{h}\subset\mathcal{U}$ on
which $\phi_{h}$ is a diffeomorphism onto its image.

$\phi_{h}\colon\mathcal{U}\longrightarrow\mathcal{U},$

depending smoothly on a real parameter $h$ such that the exact solution of the
differential equation coincides with $\phi_{h}(x)$ up to order $k$ in $h$,
i.e.
$\|\phi_{h}(x)-x(h)\|=\mathop{}\mathopen{}o\mathopen{}\left(h^{k}\right)$. The
_numerical scheme of timestep $\Delta t$_ associated to an integrator consists
in the recursive sequence

$x_{0}=x\hbox{ and }x_{n+1}=\phi_{\Delta t}(x_{n}).$

######  Remark 2.6.

Since a numerical scheme is defined by its iterations, the words integrator,
numerical method and numerical scheme can be understood without ambiguity as
synonyms all along this article.

Consider now a differential equation (4) on $\mathcal{U}\subset\mathbb{R}^{n}$
which is Hamiltonian for a Poisson structure $\pi$ and a Hamiltonian $H\in
C^{\infty}(\mathcal{U})$.

######  Definition 2.7.

An integrator of order $k$ for the Hamiltonian differential equation (4)

$\phi_{h}\colon\mathcal{U}\longrightarrow\mathcal{U},$

is said to be a _Poisson integrator_ if $\phi_{h}$ is a Poisson diffeomorphism
of $(\mathcal{U},\pi)$ for all $h$ for which it is defined333Again, it is more
rigorous to say $\phi_{h}\colon\mathcal{U}_{h}\to\phi_{h}(\mathcal{U}_{h})$ is
a Poisson diffeomorphism.

As explained in the Example 2.5, Poisson integrators can have a flow: the
trajectories may jump from one symplectic leaf to another, and thus have non-
physical behaviour. Hence, we formulate the following definition.

######  Definition 2.8.

A Poisson integrator of order $k$ for (4)

$\phi_{h}\colon\mathcal{U}\longrightarrow\mathcal{U},$

is said to be a _Poisson Hamiltonian integrator of order $k$ for (1)_ if there
exists a time-dependent Hamiltonian function $(H_{t})_{t}$, depending smoothly
on $t$, that coincides with $H$ up to order $k$, i.e.
$H_{t}-H=\mathop{}\mathopen{}o\mathopen{}\left(h^{k}\right)$ whose integral
curve coincides with the curve $h\longrightarrow\phi_{h}(x)$.

The following proposition claims that this second definition is strictly
stronger.

######  Proposition 2.9.

A Poisson Hamiltonian integrator is a Poisson integrator.

###### Proof.

The flow of a time dependent Hamiltonian differential equation is a Poisson
diffeomorphism, as long as it is well-defined. Also, if $H_{t}$ and $H$
coincide up to order $k$, then so do their Hamiltonian flows. ∎

######  Remark 2.10.

As mentioned in the introduction, in general, we will not need in Definition
2.8 to describe explicitely the family $(H_{t})$. All what matters at this
point is that it exists.

##  3\. Bi-realisations for Poisson manifolds

Having set the preliminaries and the framework in the previous sections, we
are now ready to address the core of the paper: construct Poisson Hamiltonian
integrators for a wide class of Poisson structures and any Hamiltonian
differential equation on them. One of the important "tool" for the procedure
is the notion of local symplectic groupoid associated to a Poisson structure
that arose in [3]. See [4] for a modern introduction to the matter. Notice
that we will mostly not need the whole groupoid structure but a neighborhood
of the identity of the latter, which can be considered to be closer from the
version of Karasev [15]. Several authors [10], [7] (to cite a few) have
already used studied symplectic groupoids444A Lie groupoid over a manifold,
roughly speaking, is a “higher” analog of a Lie group, where to each element
one associate two mappings to this base manifold: source $\alpha$ and target
$\beta$. Then two elements are composable when the source of one matches the
target of the other, and for those the standard group axioms are satisfied.
The symplectic form is compatible with this composition. Details are for
example in [30]. to construct numerical integrators: the relation is explained
in [1].

Symplectic groupoids of Poisson manifolds are neither easy to understand as a
notion, nor easy to construct as an object. Although we also somehow follow
the same path, our method does not use the Lie groupoid structures (product,
inverse) but only the source and target, so that all we need is what we call a
_bi-realisation_ , so we will not have to define the notion in full
generality.

In the first subsections of what follows, we explain under which circumstances
this bi-realisation, whose existence is guaranteed by theoretical arguments,
is explicitly constructable. Then, we explain why the graph of $\epsilon dH$
gives a decent Poisson integrator for the differential equation (4). It is
moreover possible to get a better Poisson integrator at an arbitrary order by
replacing $H$ by a polynomial in $\epsilon$ whose terms are computed by an
easy recursion, solving Hamilton-Jacobi equation at the desired order. Details
are developed at the beginning of section 4. The modified Hamiltonian is also
computed by recursion.

While the justification of existence and estimates are guaranteed by
complicated mathematical theorems, for implementation this section is
sufficient and it does not assume advanced differential geometry knowledge; we
again orient an interested reader to [1] for details. So throughout the
presentation in this section we will systematically make remarks on what is
computable and with what precision.

###  3.1. Bi-realisation I: Definition and existence.

Assume we are given a _bi-surjection_ , i.e. the following data:

  1. (1)

an open subset $\mathcal{U}\subset\mathbb{R}^{n}$ – the _phase space_

  2. (2)

a subset $\mathcal{W}\subset\mathcal{U}\times\mathbb{R}^{n}$ containing
$\mathcal{U}\times\\{0\\}$,

  3. (3)

two surjective submersions, called _source_ and _target_ ,
$\alpha,\beta\colon\mathcal{W}\to\mathcal{U}$ such that for all
$x\in\mathcal{U}$

$\alpha(x,0)=\beta(x,0)=x.$

We denote a bi-surjection by $\mathcal{W}\rightrightarrows\mathcal{U}$.

Bi-surjections allow to recover a diffeomorphism of $\mathcal{U}$ out of any
_bi-section_ , i.e. any submanifold $\Sigma$ of dimension $n$ in $\mathcal{W}$
to which the restrictions of the source $\alpha$ and the target $\beta$ are
diffeomorphisms onto $\mathcal{U}$. (Sometimes, we only assume that $\alpha$
and $\beta$ are diffeomorphisms onto their images, which are open subsets of
$\mathcal{U}$.) A bi-section $\Sigma$ of a bi-surjection
$\mathcal{W}\rightrightarrows\mathcal{U}$ induces a diffeomorphism
$\underline{\Sigma}\colon\mathcal{U}\longrightarrow\mathcal{U}$ defined as
$\beta\circ(\alpha_{|_{\Sigma}})^{-1}$:

$\textstyle{\mathcal{W}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\alpha}$$\scriptstyle{\beta}$$\textstyle{\ignorespaces\ignorespaces\ignorespaces\ignorespaces\Sigma\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\simeq}$$\scriptstyle{\simeq}$$\textstyle{\mathcal{U}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\underline{\Sigma}}$$\textstyle{\mathcal{U}}$

The crucial remark is that, if $\Sigma$ and
$\mathcal{W}\rightrightarrows\mathcal{U}$ are explicitly known, then the
computation $\underline{\Sigma}$ only requires to invert a diffeomorphism.
This operation, in general, can be done numerically with machine precision and
with reasonable cost, so that the diffeomorphism $\underline{\Sigma}$ can be
easily computed. The discretisations that we are going to construct are
families $(\underline{\Sigma}_{h})$ of diffeomorphisms, depending on a “small”
real parameter $h$, associated to a family $\Sigma_{h}$ of bi-sections such
that $\Sigma_{0}=\mathcal{U}\times\\{0\\}$, so that $\underline{\Sigma}_{0}$
is the identity map.

Then $\mathcal{W}\subset\mathbb{R}^{n}\times\mathbb{R}^{n}$ comes equipped
with a symplectic structure:

$\omega_{can}:=\sum_{i=1}^{n}dp_{i}\wedge dx_{i}$ (9)

with $x_{1},\dots,x_{n},p_{1},\dots,p_{n}$ being the natural variables on
$\mathcal{U}\times\mathbb{R}^{n}$, labeled in that order. The corresponding
Poisson structure satisfies:

$\\{x_{i},x_{j}\\}_{\omega}=\\{p_{i},p_{j}\\}_{\omega}=0\hbox{ and
}\\{x_{i},p_{j}\\}_{\omega}=\delta_{i}^{j}.$

We can now state the main definition:

######  Definition 3.1 (Bi-realisation).

Let $\pi$ be a Poisson structure on an open subset
$\mathcal{U}\subset\mathbb{R}^{n}$. A bi-realisation of $(\mathcal{U},\pi)$ is
a bi-surjection $\mathcal{W}\rightrightarrows\mathcal{U}$, with source
$\alpha$ and target $\beta$, for which the auxiliary dimension coincides with
the dimension of $\mathcal{U}$, satisfying the following:

  1. (1)

$\alpha$ is a Poisson map,

  2. (2)

$\beta$ is an anti-Poisson map,

  3. (3)

the fibers of $\alpha$ and $\beta$ are symplectically orthogonal to each
other.

In all three items above, $\mathcal{W}$ comes equipped with the Poisson
bracket $\\{\cdot,\cdot\\}_{\omega}$ associated to the symplectic structure
(9).

######  Remark 3.2.

Conditions 1 – 3 in Definition 3.1 mean that for all functions $F,G\in
C^{\infty}(\mathcal{U})$:

$\alpha^{*}\\{F,G\\}_{\pi}=\\{\alpha^{*}F,\alpha^{*}G\\}_{\omega},\quad\beta^{*}\\{F,G\\}_{\pi}=-\\{\beta^{*}F,\beta^{*}G\\}_{\omega},$
$\\{\alpha^{*}F,\beta^{*}G\\}_{\omega}=0.$

We will quote the following two results (without proof), for completness and
future references.

######  Proposition 3.3.

Let $\mathcal{W}\rightrightarrows\mathcal{U}$ be a bi-realisation. For any bi-
section $\Sigma\subset\mathcal{W}$ which is Lagrangian with respect to (9),
the induced diffeomorphism
$\underline{\Sigma}\colon\mathcal{U}\longrightarrow\mathcal{U}$ is a Poisson
diffeomorphism.
Moreover, for any Lagrangian submanifold of the form

$\\{(x_{1},\ldots,x_{n},\partial_{x_{1}}F,\ldots,\partial_{x_{n}}F),(x_{1},\ldots,x_{n})\in\mathcal{U}\\}\subset\mathcal{W}$

for $F\in C^{\infty}(\mathcal{U})$ a small enough555More precisely, the
$\infty$-norm of the derivative at a point must be smaller than some local
bound. smooth function, this Poisson diffeomorphism is the value at time $1$
of the flow of a time-dependent Hamiltonian vector field.

######  Theorem 3.4 (Existence and uniqueness).

Any Poisson structure on an open subset $\mathcal{U}\subset\mathbb{R}^{n}$
admits a bi-realisation. Furthermore, it is canonical in the following sense:
two different bi-realisations above a Poisson structure are symplectomorphic
through a unique symplectomorphism fixing $\mathcal{U}\times\\{0\\}.$

###  3.2. Bi-realisations II: Explicit constructions.

Theorem 3.4 states that bi-realisations do exist and are unique. Below we
explain how one can construct them.

#### Cotangent paths

Let $(M,\pi)$ be a Poisson manifold666By this we mean a collection of open
sets $\mathcal{U}$ equipped with Poisson bivectors in a consistent way. The
result being essentially local, one may think of just one open set as before,
and in the applications we will work in one coordinate chart anyway.. Out of a
path $\alpha$ valued in $T^{*}M$, two paths valued in $TM$ can be constructed:

  1. (1)

consider $\dot{\gamma}(t)$ with $\gamma=\tau\circ\alpha$ the base path of
$\alpha$, $\tau$ being the projection defining $T^{*}M$.

  2. (2)

$\alpha(t)\in T^{*}_{\gamma(t)}M$ with
$\pi_{\gamma(t)}\in\wedge^{2}T_{\gamma(t)}M$: the contraction
$\pi^{\\#}(\alpha)$ is a path valued in $TM$.

We call a path _cotangent_ when both $TM$-valued paths above coincide.

When an affine connection $\nabla$ is given on $T^{*}M$, every $\xi\in T^{*}M$
is a starting point of a parallel cotangent path $t\mapsto\xi(t)$ required to
satisfy the additional condition:

$\nabla_{\dot{\gamma}(t)}\xi=0.$

There exists a neighborhood777For purpose of notation, we denote it
$\mathcal{W}$ again, even though it is a collection of such open sets
described before. $\mathcal{W}$ of $M$ in $T^{*}M$ for which the parallel
cotangent path above is defined for all $t\in[-1,1]$. We call _geodesic flow_
the map:

${\Xi}_{\nabla,\pi}:\xi\in\mathcal{W}\to\xi(t)\in\\{\hbox{Cotangent
Paths}\\}.$

#### Karasev’s construction

Let $(M,\pi)$ be a Poisson structure with $M$ a subset of $\mathbb{R}^{n}$, so
that $T^{*}M$ can be identified with pairs $m\in M$ and
$\xi\in\mathbb{R}^{n}$. Let $\nabla$ be the trivial affine connection. For
every $\xi\in\mathbb{R}^{n}$, ${\Xi}_{\nabla,pi}(m,\xi)$ defines a path $m(t)$
solution of the differential equation:

$\dot{m_{\xi}}(t)=\pi_{m_{\xi}(t)}^{\\#}\xi.$

For every $\xi$, consider the open subsets $M^{\alpha}_{\xi}\subset M$ (resp.
$M^{\beta}_{\xi}\subset M$) on which this path $m_{\xi}$ is well defined for
all $t\in[0,1]$ (resp. all $t\in[-1,0]$).

The idea of Karasev consists in looking at the following two equations whose
unknown $\alpha,\beta$ are in $M$, for a given $(m,\xi)\in T^{*}M$:

$\int_{0}^{1}\Xi(\alpha,\xi)(t)dt=m\hbox{ and
}\int_{-1}^{0}\Xi(\beta,\xi)(t)dt=m.$

Since, for $\xi=0$, the unique solutions are $\alpha=\beta=m$, there exists a
neighborhood $\mathcal{W}$ of $M$ in $T^{*}M$ on which the two previous
equation have a unique solution, defining therefore two maps $\mathcal{W}\to
M$ that we denote $\alpha$ and $\beta$.

######  Proposition 3.5 (Karasev).

The triple $(\mathcal{W},\alpha,\beta)$ is a bi-realisation for a Poisson
structure $(M,\pi)$.

######  Remark 3.6.

This bi-realisation is explicit provided that the geodesic flow $\Xi$ and its
integral can be computed. It is computable by quadrature if so is the geodesic
flow, which is the case for a large class of Poisson structures.

#### Poisson Spray and Moser’s trick

A sligthly more academic (axiomatic) and thus conceptual approach to the above
construction may be presented using the notion of Poisson spray.

######  Definition 3.7.

Let $(M,\pi)$ be a Poisson manifold, $\tau\colon T^{*}M\to M$ the cotangent
projection and for $\lambda\in\mathbb{R}^{*},$ $m_{\lambda}\colon\xi\in
T^{*}M\mapsto\lambda\xi\in T^{*}M$ the fiberwise multiplication by $\lambda.$
$X\in\mathfrak{X}(T^{*}M)$ is said to be a Poisson spray if it verifies the
following two conditions:

  1. (1)

$\forall\xi\in T^{*}M,$ $\text{d}_{\xi}\tau.X(\xi)=\pi^{\\#}(\xi),$

  2. (2)

$X$ is homogeneous of degree 1: $\forall\lambda\in\mathbb{R}^{*},$
$\text{m}_{\lambda}^{*}X=\lambda X,$
i.e. $\text{d}_{\alpha\xi}m_{\lambda^{-1}}.X(\alpha\xi)=\lambda X(\xi).$

######  Example 3.8.

For some choice of coordinates $x_{i}$, the Poisson tensor has the form

$\pi(x)=\sum_{1\leq i<j\leq
n}\pi_{ij}(x)\partial_{x_{i}}\wedge\partial_{x_{j}}.$

Denoting $(x,\xi)$ the induced cotangent coordinates,

$X=\sum_{1\leq i<j\leq n}\pi_{ij}(x)\xi_{j}\partial_{x_{i}}$

is a Poisson spray.

The second point of its definition implies that $X$ vanishes on the zero
section $0_{T^{*}M}$. Consequently, there exists a neighborhood $\mathcal{W}$
of $0_{T^{*}M}$ such that the time-1 flow of $X$ $\Phi^{1}_{X}\colon
U\to\Phi^{1}_{X}(U)$ is a well-defined global diffeomorphism onto its image.

######  Remark 3.9.

For a given Poisson structure, Poisson sprays always exist (see [5]). However,
Poisson sprays are far from being unique. For instance, one can add a term of
the form "$f(x)\xi_{i}\xi_{j}\partial_{\xi_{j}}$" to it – this is an important
freedom that allows to construct explicit integration of the flow above in a
lot of important cases.

######  Theorem 3.10.

Any Poisson spray induces target, source and multiplicative form of the local
symplectic groupoid near $0_{T^{*}M}$ in the following way :

  1. (1)

$\bar{\alpha}=\tau\colon T^{*}M\to M,$

  2. (2)

$\bar{\beta}=\tau\circ\Phi^{1}_{X},$

  3. (3)

$\Omega=\int_{0}^{1}\Phi^{\,s}_{X}{}^{*}\omega\,\text{d}s,$ where $\omega$ is
the canonical symplectic form.

Note that $\Omega$ is symplectic up to shrinking of $\mathcal{W}.$

######  Theorem 3.11.

Any Poisson spray induces a bi-realisation.

###### Proof.

By Moser’s trick, $\omega$ and $\Omega$ are symplectomorphic in a neighborhood
of the zero section in $T^{*}M$: $\omega=\psi^{*}\Omega.$ Moreover, $\psi$ is
the identity map on $M\subset T^{*}M$, and so is its differential at any point
of $M$. Then, a bi-realisation on this neighborhood is given by :
$\alpha=\psi^{*}\bar{\alpha},\quad\beta=\psi^{*}\bar{\beta}.$ ∎

However, for a generic Poisson structure, the Poisson spray, and its flow
$\psi$, may not be explicitly computable.

###  3.3. Examples

In what follows we construct a symplectic bi-realisation for several classes
of Poisson structures, using various techniques, including Poisson sprays. We
start with the simplest Poisson structure given by a symplectic form written
in canonical (Darboux) coordinates, to recover some symplectic integrators.
Then we continue with a couple of constructions that will later be used in the
numerical tests.

#### Symplectic case.

Let $M=\mathbb{R}^{2n}=\\{(q,p)\\}$, then
$\pi=\partial_{p}\wedge\partial_{q}.$ Let $(q,p,\xi_{q},\xi_{p})$ be cotangent
coordinates on $R^{4n}.$ A Poisson spray is
$X(q,p,\xi_{q},\xi_{p})=\xi_{p}\partial_{q}-\xi_{q}\partial_{p}.$
The objects of theorem 3.10 are:

  1. (1)

$\bar{\alpha}\colon T^{*}M\to M\colon q,p,\xi_{q},\xi_{p}\mapsto(q,p),$

  2. (2)

$\bar{\beta}\colon T^{*}M\to M\colon
q,p,\xi_{q},\xi_{p}\mapsto(q+\xi_{p},p-\xi_{q}),$

  3. (3)

$\Omega=\omega+\frac{1}{2}\text{d}p\wedge\text{d}\xi_{p}-\frac{1}{2}\text{d}\xi_{q}\wedge\text{d}q-\frac{1}{3}\text{d}\xi_{q}\wedge\text{d}\xi_{p}$.

The symplectomorphism between $\Omega$ and $\omega$ is given by

$\Psi\colon\begin{pmatrix}q\\\ p\\\ \xi_{q}\\\
\xi_{p}\end{pmatrix}\mapsto\begin{pmatrix}q-\frac{\xi_{p}}{2}\\\
p+\frac{\xi_{q}}{2}\\\ \xi_{q}\\\ \xi_{p}\end{pmatrix},$

and the resulting bi-realisation is

$\left\\{\begin{array}[]{ll}\alpha:(q,p,\xi_{q},\xi_{p})\mapsto(q-\frac{1}{2}\xi_{p},p+\frac{1}{2}\xi_{q})\\\
\beta:(q,p,\xi_{q},\xi_{p})\mapsto(q+\frac{1}{2}\xi_{p},p-\frac{1}{2}\xi_{q})\end{array}\right..$

#### Quadratic Poisson structures

The following example will be important for Lotka-Volterra systems. Consider
$M=\mathbb{R}^{n}$ and a quadratic Poisson structure:

$\pi=\sum_{1\leq i,j\leq
n}a_{ij}x_{i}x_{j}\partial_{x_{i}}\wedge\partial_{x_{j}}$ (10)

Using the (natural) Poisson spray of [17]:

$X=\sum_{1\leq i,j\leq n}a_{ij}x_{i}x_{j}\xi_{i}\partial_{x_{j}}-\sum_{1\leq
i,j\leq n}a_{ij}x_{i}\xi_{i}\xi_{j}\partial_{\xi_{j}},$

and the symplectomorphism:

$\Psi\colon\begin{pmatrix}x_{j}\\\
p_{j}\end{pmatrix}\mapsto\begin{pmatrix}e^{-\frac{1}{2}\sum_{i}a_{ij}x_{i}p_{i}}x_{j}\\\
e^{\frac{1}{2}\sum_{i}a_{ij}x_{i}p_{i}}p_{j}\end{pmatrix},$

one constructs the following global bi-realisation:

$\left\\{\begin{array}[]{ll}\alpha:(x,p)\mapsto\left(e^{-\frac{1}{2}\sum_{i}a_{ij}x_{i}p_{i}}.x_{j}\right)_{j=1,\dots,n}\\\
\beta:(x,p)\mapsto\left(e^{\frac{1}{2}\sum_{i}a_{ij}x_{i}p_{i}}.x_{j}\right)_{j=1,\dots,n}\end{array}\right..$

#### Dual of a Lie algebra: cotangent lifts

In the case of the linear Poisson structure on the dual of a Lie algebra,
there is another way of constructing bi-realisations.

######  Proposition 3.12.

Let $\mathfrak{g}$ be a Lie algebra of a Lie group $G$, and $\varphi\colon
G\to\mathfrak{g}$ a local diffeomorphism in $1_{G}$, bijective on an open
subset $\mathcal{V}$ containing the unit $1_{G}$, such that :

  1. (1)

$\varphi(1_{G})=0$

  2. (2)

$T_{1_{G}}\varphi=Id.$

Let us denote by $\psi$ the inverse of $\varphi$ and
$\bar{\mathcal{V}}=\varphi(\mathcal{V}).$
Then a bi-realisation of the Lie-Poisson structure on $\mathfrak{g}^{*}$ is
given by:

$\left\\{\begin{array}[]{ll}\alpha\colon\bar{\mathcal{V}}\times\mathfrak{g}^{*}\to\mathfrak{g}^{*}:(\eta,\xi)\mapsto
R^{*}_{{\psi(\eta)}}\left(T_{\psi(\eta)}^{*}\varphi\right).\xi\\\
\beta\colon\bar{\mathcal{V}}\times\mathfrak{g}^{*}\to\mathfrak{g}^{*}:(\eta,\xi)\mapsto
L^{*}_{{\psi(\eta)}}\left(T_{\psi(\eta)}^{*}\varphi\right).\xi=\text{Ad}^{*}_{\psi(\eta)}.\alpha(\eta,\xi)\end{array}\right..$

Let us describe more precisely these source and target maps. They are the dual
of the inverse of the so-called right and left logarithmic derivatives of
$\psi$. Since $\varphi$ maps $\mathcal{V}\subset G$ to
$\bar{\mathcal{V}}\subset\mathfrak{g}$, its differential $T\varphi$ maps
$T\mathcal{V}$ to
$T\bar{\mathcal{V}}\simeq\bar{\mathcal{V}}\times\mathfrak{g}$. Composing this
map with the right and left identifications of $T\mathcal{V}$ with
$\mathcal{V}\times\mathfrak{g}$ and using the diffeomorphism $\psi,$ one gets
two families indexed by $\eta\in\bar{\mathcal{V}}$ of linear invertible
endomorphisms of $\mathfrak{g}.$ The source and targets above are the dual of
these maps.888The reader familiar with the notion of logarithmic derivative
will notice that those maps are the inverse of the dual of the logarithmic
derivative of $\psi$ after right and left trivialisations of $TG.$

######  Remark 3.13.

Notice that we no not assume $\varphi$ to be the logarithm, i.e. the inverse
of the exponential map. It may be any local diffeomorphism. In fact, the
logarithm map may not be a good choice since its differential may be too
complicated to compute.

######  Remark 3.14.

If $\varphi$ and its inverse are explicitly computable, then so are $\alpha$
and $\beta$.

######  Example 3.15.

Let us spell-out this construction in the case of the algebra $so(n)$ of anti-
symmetric matrices. The scalar product $<.,.>\colon(X,Y)\in
so(n)\mapsto\text{Tr}(X^{T}.Y)\in\mathbb{R}$ induces an isomorphism between
$so(n)$ and its dual. The local diffeomorphism we use is

$\varphi\colon SO(n)_{+}\to so(n):Q\mapsto 4\frac{Q-I}{Q+I}$

with the inverse

$\psi\colon so(n)\to SO(n)_{+}:A\mapsto\frac{4+A}{4-A}.$

Its derivative is

$T_{Q}\varphi\colon so(n)\to so(n):H\mapsto 4(I+Q^{-1})^{-1}.H.(I+Q)^{-1}$

and the transpose of it by $<.,.>$ is its cotangent lift:

$\begin{array}[]{cccc}T^{*}\varphi\colon&so(n)\times so(n)^{*}&\to&SO(n)\times
so(n)\\\
&(H,x)&\mapsto&(\psi(H),(I+\frac{A}{4}).x.(I-\frac{A}{4}).\end{array}$

Since the metric is Ad-invariant, $\text{Ad}^{*}_{Q}x=Q^{-1}.x.Q$. And the
source and target are:

$\left\\{\begin{array}[]{ll}\alpha\colon so(n)\times so(n)\to
so(n):(A,x)\mapsto(1+\frac{A}{4}).x.(1-\frac{A}{4})\\\ \beta\colon so(n)\times
so(n)\to
so(n)\colon(A,x)\mapsto(1-\frac{A}{4}).x.(1+\frac{A}{4})\end{array}\right..$

##  4\. Explicit construction of Poisson Hamiltonian integrators

We are now ready to put together all what has been discussed in the context of
Poisson geometry in the two previous sections, and make the final step to
construction of the appropriate structure preserving integrators. To sum it
up, we start with a Poisson structure $\pi$ defined on an open subset
$\mathcal{U}\subset\mathbb{R}^{n}$. The only assumption that we need it that
it admits an _explicit_ bi-realisation
$\mathcal{W}\rightrightarrows\mathcal{U}$.

We recall that $\mathcal{W}$ is an open subset of
$\mathcal{U}\times\mathbb{R}^{n}$ containing $\mathcal{U}\times\\{0\\}$. We
denote its source by $\alpha$, its target by $\beta$, and its base map by
$\tau$. We denote by ${\mathbf{0}}:\mathcal{U}\to\mathcal{W}$ the map
${\mathbf{0}}(x)=(x,0)$.

$\textstyle{\mathcal{W}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\tau}$$\scriptstyle{\alpha}$$\scriptstyle{\beta}$$\textstyle{\mathcal{U}}$$\textstyle{\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\mathcal{U}}$$\scriptstyle{{\mathbf{0}}}$$\textstyle{\mathcal{U}}$

######  Remark 4.1.

We recall that for any $(x,p)\in\mathcal{W}$, $\alpha(x,p)$ and $\beta(x,p)$
are in the same symplectic leaf of $\pi$. This leaf is not the same symplectic
leaf at the one containing $\tau(x,p)=x$.

Consider again the Hamiltonian differential equation

$\dot{x}(t)=\pi(x(t))(\nabla H(x(t)))$ (11)

for some Hamiltonian $H\in C^{\infty}(\mathcal{U})$.

We claim that we can construct an explicit Poisson Hamiltonian integrator of
order $k$ for (11). There are several steps that we now present.

  1. Step 1.

To start with, one needs to compute the first $k$ terms of the _Hamilton-
Jacobi transform_ of $H$. The latter is a formal series of smooth functions on
$\mathcal{U}$ of the form

$\mathcal{S}_{h}(H)(x)=hS_{1}(x)+\frac{h^{2}}{2}S_{2}(x)+\frac{h^{3}}{6}S_{3}(x)+\dots,$

whose coefficients are computed by recursion as follows:

    1. (a)

Set $S_{1}(x)=H(x)$.

    2. (In particular, for $k=1$, the truncation of the generating transform of $H(x)$ is simply $hH(x)$.)

    3. (b)

The smooth function $S_{k+1}(x)$ is then given by the recursive formula:

$S_{i+1}(m)=\left.\frac{d^{i}}{dt^{i}}\right|_{t=0}H\left(\alpha\left(d_{m}S_{t}^{(i)}\right)\right)$

where we write $S_{t}^{(i)}=\sum_{j=1}^{i}\frac{t^{j}}{j!}S_{j}.$

    4. Since the bi-realisation is supposed to be explicitly known, the construction of these terms can be done explicitly as well.

  2. Step 2.

Now starts the construction of the Poisson Hamiltonian integrator itself.
Choose a timestep, i.e. fix a small positive real number $h$. We define a
numerical scheme approximating the integral curve of (11) with initial value
$x_{0}$ by constructing the sequence $(x_{n})$ according to the following
recursion:

    1. (a)

Assume that for every $n\in\mathbb{N}$, the equation

$\alpha\left(y_{n},\sum_{i=1}^{k}h^{i}\,\nabla S_{i}(y_{n}))\right)=x_{n}$

admits a unique solution $y_{n}$ (otherwise, it means that the time step is
too large).

    2. (b)

Set

$x_{n+1}:=\beta\left(y_{n},\sum_{i=1}^{k}h^{i}\,\nabla S_{i}(y_{n})\right).$

######  Remark 4.2.

The computations related to formal power series in Step 1. can be done
efficiently with computer algebra tools. The resolution of the implicit
relation in Step 2 is done approximately (for example by fixed point
techniques), but can eventually be done with machine precision.

######  Example 4.3.

For $k=1$, this numerical scheme consists in mapping $x_{n}$ to
$\beta(y_{n},h\nabla H(y_{n}))$ where $y_{n}$ is the unique solution of
$\alpha(y_{n},\nabla H(y_{n}))=x_{n}$.

######  Example 4.4.

In the case of a Lie-Poisson structure on a Lie algebra $\mathfrak{g}$
equipped with a local diffeomorphism $\varphi\colon G\to\mathfrak{g}$ with
inverse $\psi\colon\mathfrak{g}\to G$, for $k=1$, our Hamiltonian Poisson
integrator consists in

  1. (1)

Compute $a\mapsto\left(\psi(\eta)\right)^{-1}T_{\eta}\psi(a)$999where
$\left(\psi(\eta)\right)^{-1}$ is the inverse of $\psi(\eta)$ for the group
law of $G$, which is a family depending on $\eta\in\bar{\mathcal{V}}$ of
diffeomorphisms $\mathfrak{g}\simeq\mathfrak{g}$. Then consider the dual of
its inverse, which is now a family of maps
$D_{\eta}\colon\mathfrak{g}^{*}\simeq\mathfrak{g}^{*}$ being the identity map
for $\eta=0$. Then solve $D_{\Delta t\nabla H(y_{n})}(y_{n})=x_{n}$.

  2. (2)

Consider

$x_{n+1}={\mathrm{Ad^{*}}}_{\psi(\Delta t\nabla y_{n}H)}x_{n}$

######  Remark 4.5.

By construction, $x_{n}$ and $x_{n+1}$ belong to the same symplectic leaf. But
to go from $x_{n}$ to $x_{n+1}$ one uses a point $y_{n}$ which is not, in
general, on that common symplectic leaf. This is extremely counter-intuitive.

######  Remark 4.6.

Let us give the first terms of the Hamilton-Jacobi transform:

  1. (1)

$S_{1}=H$

  2. (2)

$S_{2}={\mathbf{0}}^{*}\left(\frac{1}{2}\\{\alpha^{*}H,\tau^{*}S_{1}\\}\right)$

  3. (3)

$S_{3}={\mathbf{0}}^{*}\left(\frac{1}{3}\\{\alpha^{*}H,\tau^{*}S_{2}\\}+\frac{1}{6}\\{\\{\alpha^{*}H,\tau^{*}S_{1}\\},\tau^{*}S_{1}\\}\right)$

  4. (4)

$S_{4}={\mathbf{0}}^{*}\left(\frac{1}{4}\\{\alpha^{*}H,\tau^{*}S_{3}\\}+\frac{1}{3!}\\{\\{\alpha^{*}H,\tau^{*}S_{2}\\},\tau^{*}S_{1}\\}+\right.$
$\left.+\frac{1}{4!}\\{\\{\\{\alpha^{*}H,\tau^{*}S_{1}\\},\tau^{*}S_{1}\\},\tau^{*}S_{1}\\}\right)$

Above, ${\mathbf{0}}^{*}$ means that the function on $\mathcal{W}$ is
restricted to $\mathcal{U}\times\\{0\\}$, and therefore considered as a
function on $\mathcal{U}$.

######  Theorem 4.7.

The above numerical scheme defines a Hamiltonian Poisson integrator at order
$k$ for the Hamiltonian differential equation (11).

###### Proof.

(See [1]) Let us give a brief outline of the proof of this theorem, which will
also explain how the time-dependent Hamiltonian $H_{t}$ whose flow at time $h$
matches exactly $x_{n}$ to $x_{n+1}$ is constructed.

To start with, recall the two points that we saw in Section 3 about the set
$Gr(dS):=\left\\{(x,d_{x}S),x\in\mathcal{U}\right\\}$ for $S\in
C^{\infty}(\mathcal{U})$

  1. (1)

It is a Lagrangian subset of $\mathcal{U}\times\mathbb{R}^{n}$.

  2. (2)

Provided that the differential of the function $S$ is small enough, it is a
bisection101010Meaning that restrictions of $\alpha$ and $\beta$ to $Gr(dS)$
are invertible. of $\mathcal{W}\rightrightarrows\mathcal{U}$.

As a consequence, as we saw in Proposition 3.3, the map
$\beta\circ\alpha^{-1}|_{Gr(dS)}$ is a Poisson map.

It is a more subtle result that if the function $S$ has a small enough
differential, then the Poisson morphism is the time $1$-flow for a time
dependent Hamiltonian function. For instance, if $S$ depends on a parameter
$h$, i.e. $S=S_{h}$ where $S_{t}$ is a time dependent function with $S_{0}=0$,
then this time dependent Hamiltonian function is given by:

$\tilde{S}_{t}(x):=\frac{\partial S_{t}}{\partial t}(y)$

where $y$ is chosen such that $(y,\nabla S_{t}(y))\in Gr(dS_{t})$. Afterwards,
the question reduces to finding $S_{t}$ such that the flow of $\tilde{S}_{t}$
at time $h$ matches the flow of $H$ at time $h$ up to order $k$ in the
variable $h$, as in Step 2. ∎

######  Remark 4.8.

At order $k=1$ in the symplectic case (i.e. non-degenerate constant Poisson
structure), it is easy to check that for the harmonic oscillator
$H=\frac{p^{2}+q^{2}}{2}$, one recovers the symplectic mid-point scheme. For a
general Hamiltonian $H,$ the present construction gives the fact that an
implicit Euler scheme of timestep $\frac{\Delta t}{2}$ composed with an
explicit Euler scheme of timestep $\frac{\Delta t}{2}$ is a symplectic
integrator of order $1$ and timestep $\Delta t.$ More generally, for higher
orders the constructed Poisson Hamiltonian integrators for symplectic
structures will be symplectic integrators, but a priori different from the
standard symplectic Runge-Kutta methods.

######  Remark 4.9.

We have mentioned in the introduction that the naive idea “restrict to a leaf,
be symplectic there” to recover Poisson globally, does not work because is
almost never constructive. But the other way around it is actually fruitful:
now having constructed a Poisson Hamiltonian integrator forcing the trajectory
to stay on the correct leaf, one can apply the backward analysis techniques
(restricted to leaves) for error estimates.

######  Remark 4.10.

Recall that in the case of linear Poisson structures of Proposition 3.12, the
construction of the bi-realisation amounts to computation of the coadjoint
action of $G$ on $\mathfrak{g}$, and construction of a local diffeomorphism:
$\phi\colon G\to\mathfrak{g}$ with its differential at $1$ being the identity.

The obtained Hamiltonian Poisson integrator of order $1$ is of the form:

$x_{n+1}={\mathrm{Ad^{*}}}_{\phi^{-1}(\Delta t\nabla y_{n}H)}x_{n}$

which is certainly not surprising: any such numerical scheme stays in the
symplectic leaf where one starts from. The same remark about the point $y_{n}$
outside this leaf holds.

An obvious choice for $\phi$ is the inverse of the exponential map, but there
is some freedom in it: any such a local diffeomorphism can be used to compute
an Hamiltonian Poisson integrator up to order $k$. It is important, however,
to be able to compute easily its differential and its inverse.

##  5\. Numerical tests

In this last section we illustrate the advantages of Poisson Hamiltonian
integrators on a couple of examples.

###  5.1. The Rigid Body

First turn to the linear Poisson structures – a good example of those can be
provided by the dynamics of a rigid body about a periodic orbit.

The equations governing the system read:

$\dot{x}=-x\wedge J.x,$

where: $\wedge$ denotes the vector product in $\mathbb{R}^{3},$ and the
symmetric positive matrix $J$ is the inertia tensor of the body.

It is a Hamiltonian differential equation for
$\pi(x)=\begin{pmatrix}0&-x_{3}&x_{2}\\\ x_{3}&0&-x_{1}\\\
-x_{2}&x_{1}&0\end{pmatrix}$ and
$H(x)=\frac{1}{2}\text{Tr}\Big{(}j(x)^{T}.J.j(x)\Big{)}$ where
$j\colon\mathbb{R}^{3}\widetilde{\longrightarrow}so(3)$ given by:

$x\longmapsto\begin{pmatrix}0&-x_{3}&x_{2}\\\ x_{3}&0&x_{1}\\\
-x_{2}&-x_{1}&0\end{pmatrix}.$

We consider the inertia tensor $J=\begin{pmatrix}1&0&0\\\ 0&\pi&0\\\
0&0&100\end{pmatrix}$ and $x_{0}=\begin{pmatrix}1\\\ 1\\\ 1\end{pmatrix}$ so
that the trajectory is given by Figure 4.

Figure 4. The trajectory of the angular velocity of a rigid body in
$\mathbb{R}^{3}.$

Numerical simulations are for timestep $\Delta t=10^{-4}.$ The Poisson
Hamiltonian integrator of order 2 behaves much better than the Runge-Kutta
method of order 4 in the preservation of both Casimir and Hamiltonian levels
(Figure 5).

Details are a bit more involved. The error of the traditional method depends
linearly on the number of iterations and so diverges from the continuous
(closed) trajectory. A Poisson Hamiltonian integrator preserves Casimir level
at machine precision and oscillates around a Hamiltonian value with the
amplitude depending on $\Delta t^{k},$ $k$ being the order of the method. We
emphasize that this distance does not depend on the amount of iteration. One
recovers a typical stability phenomenon of symplectic integrators, already
noticed and explained in [11]. A zoom on Hamiltonian errors is made in Figure
5(a). For Poisson Hamiltonian integrators, a theoretical explanation relies on
the Magnus formula for Poisson structures introduced in [1]. Those phenomena
are illustrated on the schematic section of the trajectory – Figure 6.

(a) Errors on Hamiltonian values for PHI-2

(b) Error on Casimir values

(c) Errors on Hamiltonian values for RK-4

Figure 5. Comparison between Runge-Kutta 4 and our Poisson Hamiltonian
integrator at order 2 for the Rigid Body dynamics

######  Remark 5.1.

The Casimir is the square of the norm. Hence Figure 5(b) indicates that RK-4
iterations will converge to $0$ in $\mathbb{R}^{3},$ which is a fixed point of
the dynamics as well as a singular leaf of the foliation of the total space.
This lead on long run simulations to pathological behaviours. In turn, it
stresses the importance of numerical methods preserving leaves of a singular
foliation such as the ones appearing in Poisson structures.

(a) A Poincaré section of the trajectory (green plane)

(b) Intersection between (interpolated) discrete trajectories and the Poincaré
section

Figure 6. Illustration of the difference of behaviour of two numerical methods
with respect to Hamiltonian and Casimir levels

###  5.2. The Lotka-Volterra System

Recall the behaviour of the PHI and RK-2 from section 1 – we now have the
correct language to explain it, studying in particular the Casimirs of the
system.

The Poisson structure of the generic Lotka-Volterra system coincides with the
quadratic one of Equation (10), fully encoded in an $n\times n$ matrix $A$.
The dynamics is governed by a linear Hamiltonian $H=\sum_{i=1}^{n}x_{i}$.

######  Proposition 5.2.

Let $u\in\text{Ker}A.$

$\displaystyle f\colon\mathbb{R}^{n}$ $\displaystyle\to\mathbb{R}$
$\displaystyle x$ $\displaystyle\mapsto\prod_{1\leq i\leq n}x_{i}^{u_{i}}$

is a local Casimir of the quadratic Poisson structure given by the matrix $A.$

In the numerical test, we considered $n=3$ and, $A=\begin{pmatrix}0&1&1\\\
-1&0&1\\\ -1&-1&0\end{pmatrix}$ so that the local Casimir is $C\colon
x\mapsto\frac{x_{1}x_{3}}{x_{2}}.$ Generic symplectic leaves are hyperbolas,
and intersecting them with the level surfaces of $H$ one obtains the geometry
of the real trajectory (up to time parametrisation). A Runge-Kutta method does
not preserve $C$, while the constructed Poisson Hamiltonian integrator can
preserve the Casimir value with machine precision and the Hamiltonian up to
any given order in timestep. Clearly the singular behaviour of the trajectory
can be observed only provided these conservation laws are respected. Figure 7
enlightens the stability of a Poisson Hamiltonian integrator in the
neighborhood of a singularity observed in Section 1: on top of preserving the
Poisson structure, it stays on a symplectic leaf along iterations.

Figure 7. Comparison on Casimir values between PHI-1 and RK-2

## Conclusion / perspectives

In this paper, we have explained how the idea of the groupoid construction
from [1] can be implemented for design of Poisson integrators. Let us stress
again that the term _Poisson Hamiltonian integrators_ we have introduced is
important – it explains the conceptual difference to straightforward
constructions present in literature.

We have seen that even for _simple academic examples_ in _generic situations_
constructed Poisson Hamiltonian integrators proved to be more accurate than
even higher order classical methods, especially on long run simulations. But a
similar strategy can be implemented with no changes for more complicated
systems of ordinary differential equations – we are working on a symbolic
package for automatic generation of the simulation source codes for that
([2]). Moreover, similar methods can be designed even for Poisson Hamiltonian
partial differential equations, which often appear in fluid dynamics and waves
simulations. The key idea there is to use the locality of discretisation in
space to spell-out the groupoid structure maps – we intend to explore this
direction in further works.

Acknowledgments. We are thankful to participants of the Geometry and Mechanics
working group (La Rochelle, M2N team) and the Seminar on Geometry, Mechanics
and Control (ICMAT – IMAULL) for their valuable feedback. The last section
benefits from fruitful remarks of the CNRS Research Regroupement “Differential
Geometry and Mechanics”. We appreciate enlightening discussion with Pol
Vanhaecke, Dina Razafindralandy and Aziz Hamdouni at various stages of this
work. We are also thankful to Antoine Falaize for his help in implementation
of the symbolic computations mentioned in remark 4.2.

## References

  * [1] O. Cosserat, _Symplectic groupoids for Poisson integrators_ , Journal of Geometry and Physics, 2022\.
  * [2] O. Cosserat, A. Falaize, V. Salnikov, _On automatic generation of Poisson numerical methods of higher order_ , in preparation.
  * [3] A. Coste, P. Dazord, A. Weinstein, _Groupoïdes symplectiques_ , Publications du Département de Mathématiques de Lyon, 1987\.
  * [4] M. Crainic, R. Fernandes, I. Mărcuţ, _Lectures on Poisson Geometry_ , American Mathematical Society.
  * [5] M. Crainic, I. Mărcuţ, _On the existence of symplectic realizations_ , Journal of symplectic geometry, 2010\.
  * [6] J.-P. Dufour, N. Tien Zung, _Poisson Structures and their Normal Forms_ , Birkhäuser Verlag, 2005\.
  * [7] S. Ferraro, M. de Leon, J. C. Marrero, D. Martın de Diego, M. Vaquero, _On the Geometry of the Hamilton–Jacobi Equation and Generating Functions_ , Archive for Rational Mechanics and Analysis, 2017\.
  * [8] F. Gay-Balmaz, D. D. Holm, V. Putkaradze, T. S. Ratiu, _Exact geometric theory of dendronized polymer dynamics_ , Advances in Applied Mathematics, 2011\.
  * [9] F. Gay-Balmaz, T. S. Ratiu, C. Tronci, _Equivalent Theories of Liquid Crystal Dynamics_ , Archive for Rational Mechanics and Analysis, 2013\.
  * [10] Z. Ge, _Generating Functions, Hamilton-Jacobi Equations and Symplectic Groupoids on Poisson Manifolds_ , Indiana University Mathematics Journal, 1990\.
  * [11] G. Benettin, A. Giorgilli, _On the Hamiltonian Interpolation of Near-to-the-Identity Symplectic Mappings with Application to Symplectic Integration Algorithms_ , Journal of Statistical Physics, 1994\.
  * [12] M. Grmela, _GENERIC guide to the multiscale dynamics and thermodynamics_ , Journal of Physics Communications, 2018\.
  * [13] E. Hairer, C. Lubich, G. Wanner, _Geometric Numerical Integration_ , Springer Series in Computational Mathematics, 2002\.
  * [14] L.O. Jay, _Preserving Poisson structure and orthogonality in numerical integration of differential equations_ Computers & Mathematics with Applications, 2004\.
  * [15] M. V. Karasev, _Analogues of the objects of Lie group theory for nonlinear Poisson brackets_ Mathematics of the USSR-Izvestiya, 1987\.
  * [16] T. E. Kouloukas, G. R. W. Quispel and P. Vanhaecke, _Liouville integrability and superintegrability of a generalized Lotka-Volterra system and its Kahan discretization_ , Journal of Physics A: Mathematical and Theoretical, Volume 49, Number 22, 2016\.
  * [17] S. Li and D. Rupel, _Symplectic groupoids for cluster manifolds_ , Journal of Geometry and Physics, 2018\.
  * [18] C. Liu, _A Lie–Poisson bracket formulation of plasticity and the computations based on the Lie-group $SO(n)$_, International Journal of Solids and Structures, 2013\.
  * [19] P. Libermann, C.-M. Marle, _Symplectic Geometry and Analytical Mechanics_ , Kluwer Academic Publishers, 1987\.
  * [20] K. Mackenzie, _Lie Groupoids and Lie Algebroids in Differential Geometry_ , Cambridge University Press, 1987\.
  * [21] C.-M. Marle, _Géométrie des systèmes mécaniques à liaisons actives_ , Compte-Rendu de l’Académie des Sciences, 1990\.
  * [22] R. I. McLachlan, _Explicit Lie-Poisson Integration and the Euler Equations_ , Physical Review Letters, 1993\.
  * [23] A. Pichereau, C. Laurent-Gengoux, P. Vanhaecke, _Poisson structures_ , Springer-Verlag, 2012\.
  * [24] D. Razafindralandy, A. Hamdouni, _A review of some geometric integrators_ , Advanced Modeling and Simulation in Engineering Sciences, 5:16, 2018.
  * [25] V. Salnikov, A. Hamdouni, D. Loziienko, _Generalized and graded geometry for mechanics: a comprehensive introduction_ , Mathematics and Mechanics of Complex Systems, Vol. 9, No. 1, 2021\.
  * [26] Loup Verlet, _Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules_ , Phys. Rev. 159, 98, 1967\.
  * [27] V. Volterra, _Leçons sur la Théorie Mathématique de la Lutte pour la Vie_ , Gauthier-Villars et Cie, 1931\.
  * [28] P. Xu, _Morita Equivalence of Poisson Manifolds_ , Communications in Mathematical Physics, 1991\.
  * [29] H. Yoshida, _Construction of higher order symplectic integrators_ , Phys. Lett. A 150 (5–7): 262, 1990\.
  * [30] A. Cannas da Silva, A. Weinstein, Geometric models for noncommutative algebras, Berkeley Mathematics Lecture Notes, 10, AMS, 1999


# Chemical kinetics and stochastic differential equations

Chiara Pezzotti and Massimiliano Giona<EMAIL_ADDRESS>Dipartimento di Ingegneria Chimica, Materiali, Ambiente La Sapienza Università
di Roma
Via Eudossiana 18, 00184 Roma, Italy

###### Abstract

We propose a general stochastic formalism for describing the evolution of
chemical reactions involving a finite number of molecules. This approach is
consistent with the statistical analysis based on the Chemical Master
Equation, and provides the formal setting for the existing algorithmic
approaches (Gillespie algorithm). Some practical advantages of this
formulation are addressed, and several examples are discussed pointing out the
connection with quantum transitions (radiative interactions).

All the chemical physical processes involve, in an atomistic perspective, a
stochastic description of the events, be them reactive or associated with a
change of phase (for instance adsorption) chemfluct . Nonetheless, in the
overwhelming majority of the cases of practical and laboratory interest, the
number of molecules involved is so large to justify a mean field approach,
essentially based on the Boltzmannian hypothesis of molecular chaos (the
“stosszahlansatz”) boltzmann . The mean field formulation represents the
backbone of the classical theory of chemical reaction kinetics chemkin1 ;
chemkin2 .

It is well known that, in all the cases where the number of molecule is small
(and this occurs in subcellular biochemical reactions, in nanoscale systems,
or in the growth kinetics of microorganisms bio1 ; bio2 ; bio3 ), the effects
of fluctuations become significant, motivating a stochastic description of
chemical kinetic processes, involving the number of molecules present in the
system, thus explicitly accounting for due to their finite number mcquarrie ;
stoca ; stoca1 ; stoca2 . The statistical theory of chemical kinetics in these
conditions is grounded on the Chemical Master Equation (CME) gillespierigorous
; cme3 , expressing the evolution equation for the probabilities $p({\bf
N},t)$ of all the possible number-configurations ${\bf
N}(t)=(N_{1}(t),\dots,N_{s}(t))$, where $N_{h}(t)$ is the number of molecules
of the $h$-th reacting species at time $t$, $h=1,\dots,s$. However, apart from
a handful of simple cases, for which the CME can be solved analytically cmecpl
, numerical methods should be applied to it in order to compute mean values
and higher-order moments. But also this choice reveals itself to be unfeasible
in most of the situations of practical and theoretical interests, due to the
extremely large number of configurations involved, making the multi-index
matrix $p({\bf N},t)$ so huge to exceed reasonable computational facilities.

In order to solve this problem, Gillespie proposed an algorithmic solution to
the numerical simulation of stochastic reacting systems, based on the
Markovian nature of the reactive events gillespiegeneral ; gillespieexamples .
The original Gillespie algorithm has been extended and improved over time,
providing a variety of slightly different computational alternatives. A common
denominator of the first family of the Gillespie algorithms (namely those
based on the direct method, the first reaction method or their derivates ff1 ;
ff2 ; ff3 ) is to associate to every time step the occurrence of just one
reaction. This formulation comes directly from the assumption that, if the
time step is small enough, the probability that more than one reaction will
occur is negligible. While correct, this choice brings to significant
computational costs for complex reaction schemes. This problem has been
highlighted several times, from the Gillespie group itself, as _stiffness_ in
stochastic chemical reacting systems stiffness . A brilliant way to overcome
this limit originates the famous tau-leaping method, which, unfortunately,
requires to check that the propensity functions remain _almost constant_ at
each iteration and can be applied just if this condition is verified tau1 ;
tau2 . The algorithmic solution associated with the formalism here introduced
combines the accuracy of the first SSA with the computational advantages of
the $\tau$-leaping method.

There is, moreover, a missing link between the CME theory and the Gillespie
algorithm, consisting in the straight mathematical formulation of the
stochastic differential equations associated with a chemical reacting system,
the statistical description of which would correspond to the CME. To clarify
this issue, consider the conceptually analogous problem of particle diffusion
over the real line, the statistical description of which is expressed by the
parabolic equation $\partial p(x,t)/\partial t=D\,\partial^{2}p(x,t)/\partial
x^{2}$, for the probability density $p(x,t)$ of finding a particle at position
$x$ at time $t$. Setting $x_{n}=x(n\,\Delta t)$, an algorithm describing this
process can be simply expressed by the discrete evolution equation
$x_{n+1}=x_{n}+\sqrt{2\,D\,\Delta t}\,r_{n+1}$, where $r_{h}$, $h=1,2,\dots$
represent independent random variables sampled from a normal distribution
(with zero mean, and unit variance) diffusion . This represents an efficient
algorithmic solution of the problem, whenever the time resolution $\Delta t$
is small enough. Nevertheless, the mere algorithmic approach cannot be
considered physically satisfactory, in a comprehensive formulation of
transport theory embedded in a continuous space-time (in which both position
$x$ and time $t$ are real valued). In point of fact, only with the
mathematical formulation due to K. Ito of stochastic differential equations
driven by the increments $dw(t)$ of a Wiener process (Langevin equations) ito
, namely $dx(t)=\sqrt{2\,D}\,dw(t)$ the theory of diffusive motion has found a
proper mathematical physical setting.

A similar situation applies to the case of stochastic models of chemical
reaction kinetics, and the present Letter is aimed at filling this gap. The
basic idea is that any reactive process corresponds to a system of elementary
events (the single reaction) possessing a Markovian transitional structure,
and, consequently, amenable to a description by means of the increments of
counting processes (Poisson processes, in the Markovian case). This topic has
been also pointed out in closed1  in terms of Poisson measures, although the
latter formulation is much less simple and physically intuitive than the
approach proposed in the present Letter.

To begin with, consider the simple case of a first-order chemical reaction
$A\underset{k_{-1}}{\stackrel{{\scriptstyle k_{1}}}{{\rightleftharpoons}}}B$
(for instance, an isomerization). This model is perfectly analogous to the
radiative transition of a molecule possessing two energy states, due to
emission and adsorption of an energy quantum (figure 1). Let
$N_{A}(0)+N_{B}(0)=N_{g}$ the total number of molecules at time $t=0$. The
state of the system is characterized by the state functions $\sigma_{h}(t)$,
$h=1,\dots,N_{g}$ for each molecule, attaining values $\\{0,1\\}$, and such
that $\sigma_{h}(t)=0$ if the energy state at time $t$ is $E_{0}$ (or
equivalently if the molecule finds itself in the state $A$), and
$\sigma_{h}(t)=1$ in the opposite case (energy state $E_{1}$, or isomeric
state $B$).

Figure 1: Schematic representation of the analogy between a two-level quantum
system and a first-order chemical kinetics, such as an isomerization.

Let $\\{\chi_{h}^{(1)}(t,k_{1}),\chi_{h}^{(2)}(t,k_{-1})\\}_{h=1}^{N_{g}}$ be
two systems of independent Poisson processes, characterized by the transition
rates $k_{1}$, and $k_{-1}$, respectively. The evolution of $\sigma_{h}(t)$
can be expressed via the stochastic differential equation

$\frac{d\sigma_{h}(t)}{dt}=\left(1-\sigma_{h}(t)\right)\,\frac{d\chi_{h}^{(1)}(t,k_{1})}{dt}-\sigma_{h}(t)\,\frac{d\chi_{h}^{(2)}(t,k_{-1})}{dt}$
(1)

$h=1,\dots,N_{g}$, where $d\chi(t,\lambda)/dt$ is the distributional
derivative of the Poisson process $\chi(t,\lambda)$, corresponding to a
sequence of unit impulsive functions at the transition instants $t_{i}^{*}$,
$i=1,2,\dots$, $0<t_{i}^{*}<t_{i+1}^{*}$, where for $\varepsilon>0$,
$\lim_{\varepsilon\rightarrow
0}\int_{t_{i}^{*}-\varepsilon}^{t_{i}^{*}+\varepsilon}d\chi(t,\lambda)=1$.
Summing over $h=1,\dots N_{g}$, and observing that
$N_{A}(t)=\sum_{h=1}^{N_{g}}\left(1-\sigma_{h}(t)\right)$,
$N_{B}(t)=\sum_{h=1}^{N_{g}}\sigma_{h}(t)$, we have

$\frac{dN_{B}(t)}{dt}=\sum_{h=1}^{N_{A}(t)}\frac{d\chi_{h}^{(1)}(t,k_{1})}{dt}-\sum_{h=1}^{N_{B}(t)}\frac{d\chi_{h}^{(2)}(t,k_{-1})}{dt}$
(2)

and $dN_{A}(t)/dt=-dN_{B}(t)/dt$, representing the evolution equation for
$N_{A}(t)$ and $N_{B}(t)$, attaining integer values. The stochastic evolution
of the number of molecules $N_{A}(t)$, $N_{B}(t)$ is thus expressed as a
differential equation with respect to the continuous physical time
$t\in{\mathbb{R}}^{+}$, over the increments of a Poisson process. Intepreted
in a mean-field way, if $c_{\rm tot}$ is the overall concentration of the
reactants at time $t=0$, then the concentrations $c_{\alpha}(t)$ at time $t$
can be recovered from eq. (2) as

$c_{\alpha}(t)=c_{\rm tot}\,\frac{N_{\alpha}(t)}{N_{g}}\,,\qquad\alpha=A,B$
(3)

representing the calibration relation connecting the stochastic description in
terms of number of molecules $N_{\alpha}(t)$ and the concentrations
$c_{\alpha}(t)$, $\alpha=A,\,B$ entering the mean-field description.

The analytical formulation of a stochastic differential equation for chemical
kinetics, expressed in terms of the number of molecules of the chemical
species involved, rather than an algorithm defined for discretized times,
permits to develop a variety of different numerical strategies, that naturally
perform a modified tau-leaping procedure, as the occurrence of several
distinct reactive events in any elementary time step $\Delta t$ is
intrinsically accounted for. This can be easily seen by considering the simple
reaction defined by the evolution equation (2). In terms of increments, eq.
(2) can be written as
$dN_{B}(t)=\sum_{h=1}^{N_{A}(t)}d\chi^{(1)}(t,k_{1})-\sum_{h=1}^{N_{B}(t)}d\chi^{(2)}(t,k_{-1})$.
If $\Delta t$ is the chosen time step, it follows from this formulation, a
simple numerical approximation for eq. (2), namely,

$\Delta N_{B}(t)=N_{B}(t+\Delta
t)-N_{B}(t)=\sum_{h=1}^{N_{A}(t)}\xi_{h}^{(1)}(k_{1}\,\Delta
t)-\sum_{h=1}^{N_{B}(t)}\xi_{h}^{(2)}(k_{-1}\,\Delta t)$ (4)

where $\xi^{(1)}(k_{1}\,\Delta t)$, $\xi_{h}^{(2)}(k_{-1}\,\Delta t)$
$h=1,2,\dots$, are two families of independent binary random variables, where

$\xi^{(\alpha)}_{h}(p)=\left\\{\begin{array}[]{lll}1&&\mbox{with
probability}\;p\\\ 0&&\mbox{otherwise}\end{array}\right.$ (5)

$\alpha=1,2$, $h=1,2,\dots$. The time step $\Delta t$, can be chosen in eq.
(4) from the condition

$K\,\Delta t<1\,,\qquad K=\max\\{k_{1},k_{-1}\\}$ (6)

In practice, we choose $\Delta t=0.1/K$. As can be observed, the choice of
$\Delta t$ is limited by the intrinsic rates of the process. The advantage of
deriving different algorithmic schemes for solving numerically the stochastic
kinetic equations becomes more evident in dealing with bimolecular reactions
(addressed below). Due to the intrinsic limitations of this communication, a
thorouh discussion of this issue is postponed to a future more extensive
article PG_next .

The same approach can be extended to include amongst the elementary events not
only the reactive steps, but also feeding conditions, thus representing the
evolution of chemically reacting systems with a finite number of molecules in
a perfectly stirred open reactor. This is the case of the tank-loading
problem, in which a tracer is injected in an open vessel assumed perfectly
mixed, for which, in the absence of chemical reactions, the mean field
equation for the concentration of the tracer reads

$\frac{dc(t)}{dt}=D\,\left(c_{0}-c(t)\right)$ (7)

where $c_{0}$ is the inlet concentration and $D$ the dilution rate (reciprocal
of the mean retention time), and $c(0)=0$. Fixing $N_{g}$ so that
$c(t)=c_{0}\,N(t)/N_{g}$, the corresponding stochastic differential equation
for the integer $N(t)$ involves, also in this case, two families of counting
processes, one for the loading at constant concentration $c_{0}$, and the
other for tracer discharge in the outlet stream, characterized by the same
transition rate $D$,

$\frac{dN(t)}{dt}=\sum_{h=1}^{N_{g}}\frac{d\chi_{h}^{(1)}(t,D)}{dt}-\sum_{k=1}^{N(t)}\frac{d\chi_{h}^{(2)}(t,D)}{dt}$
(8)

starting from $N(0)=0$. Figure 2 depicts several realizations of the tank-
loading process, obtained by discretizing eq. (8) with a time step $\Delta
t=10^{-3}$.

Figure 2: $c(t)=N(t)/N_{g}$ vs $t$ from a single realization of the tank-
loading process eq. (8) with $D=1$, $c_{0}=1$ a.u.. Panel (a): $N_{g}=30$,
panel (b) $N_{g}=100$, panel (c) $N_{g}=1000$. The solid horizontal lines
represent the steady-state value $c^{*}=1$.

Despite the simplicity of the process, this example permits to highlight the
role of $N_{g}$, that can be referred to as the granularity number, and the
way stochastic models of chemical reactions can be fruitfully applied. Indeed,
there is a two-fold use of the stochastic formulation of chemical kinetic
schemes. The first refers to a chemical reacting system involving a small
number of molecules, and in this case $N_{g}$ represents the effective number
of molecules present in the system. The other is to use stochastic algorithms
for simulating reacting systems in an alternative (and sometimes more
efficient way) with respect to the solution of the corresponding mean-field
equations. In the latter case, the granularity number $N_{g}$ represents
essentially a computational parameter, tuning the intensity of the
fluctuations. Two choices are then possible: (i) it can be chosen large
enough, in order to obtain from a single realization of the process an
accurate approximation of the mean-field behavior, or (ii) it can be chosen
small enough in order, to deal with extremely fast simulations of a single
realization of the process, that could be averaged over a statistically
significant number of realizations in due time. These two choices are depicted
in figure 2 (panel c), choosing $N_{g}=10^{3}$, and in figure 3 panel (a)
obtained for $N_{g}=30$. Of course, the latter approach is valid as long as
the low-granularity (low values of $N_{g}$) does not influence the qualitative
properties of the kinetics.

Figure 3: Panel (a): $\langle c\rangle(t)$ vs $t$ at $N_{g}=30$ (symbols)
averaged over $[10^{6}/N_{g}]$ realizations of the tank-loading process with
$D=1$, $c_{0}=1$ a.u. Here, $[\cdot]$ indicates the integer part of its
argument. The solid line represents the mean-field result $\langle
c\rangle(t)=1-e^{-t}$. Panel (b): Variance $\sigma_{c}(t)$ vs $t$ for the
tank-loading process. Symbols are the results of stochastic simulations of eq.
(8) averaged over $[10^{6}/N_{g}]$ realizations, lines the solutions of eq.
(10). Line (a) refers to $N_{g}=30$, line (b) to $N_{g}=100$.

The second (computational) use of stochastic simulations of chemical kinetics
requires a further discussion. At a first sight, it may appear that any
stochastic simulation would be computationally less efficient than the
solution of the corresponding mean-field equations. This is certainly true for
classical chemical reaction schemes in a perfectly mixed system, for which the
mean-field model reduces to a system of ordinary differential equations for
the concentrations of the reactants. But there are kinetic problems e.g.,
associated with the growth of microorganisms and eukaryotic cell lines in
bioreactors (these growth phenomena, are indeed amenable to a description in
terms of equivalent chemical reactions), the mean-field model of which is
expressed in the form of higher-dimensional nonlinear integro-differential
equations . For this class of problems, the use of stochastic simulations is
the most efficient, if not the only way to achieve a quantitative description
of the process, in those cases where the number $n_{p}$ of internal parameters
describing the physiological state of an eukaryotic cell becomes large enough,
$n_{p}\geq 3$. This issue is addressed in detail in PBGbio . This case is
altogether similar to some transport problems, such as Taylor-Aris dispersion
for high Péclet numbers or the analysis of microfluidic separation processes
(DLD devices) for which the stochastic simulation of particle motion is far
more efficient that the corresponding solution of the corresponding mean-field
model expressed in the form of advection-diffusion equations taris ; dld .

To complete the analysis of the tank-loading problem, the associated CME reads

$\frac{dp(n,t)}{dt}=D\,N_{g}\,\left[p(n-1,t)\,\eta_{n-1}-p(n,t)\right]+D\left[(n+1)\,p(n+1,t)-n\,p(n,t)\right]$
(9)

where $\eta_{h}=1$ for $h\geq 0$ and $\eta_{h}=0$ otherwise. It follows that
$\langle c\rangle(t)=c_{0}\sum_{n=1}^{\infty}n\,p(n,t)/N_{g}$ satisfies
identically the mean-field equation (due to the linearity of the problem),
while the variance $\sigma_{c}(t)$, with
$\sigma_{c}^{2}(t)=c_{0}^{2}\sum_{n=1}^{\infty}n^{2}\,p(n,t)/N_{g}^{2}-\left(c_{0}\sum_{n=1}^{\infty}n\,p(n,t)/N_{g}\right)^{2}$,
satisfies the equation

$\frac{d\sigma_{c}^{2}}{dt}=-2\,D\,\sigma_{c}^{2}+D\,\left(\frac{1}{N_{g}}+\frac{\langle
c\rangle}{N_{g}}\right)$ (10)

Figure 3 panel (b) compares the results of stochastic simulations against the
solutions of eq. (10) for two values of $N_{g}$.

The above approach can be extended to any system of nonlinear reaction schemes
involving unimolecular and bimolecular reaction, and in the presence of
slow/fast kinetics. The structure of the reaction mechanism can be arbitrarily
complicated without adding any further complexity (other than purely
notational) in the formulation of the stochastic evolution expressed in terms
of number of molecules. The only practical issue, is that the number of
different families of stochastic processes grows with the number of elementary
reactive processes considered. For instance, in the case of the subtrate-
inhibited Michaelin-Menten kinetics

$\displaystyle E+S\underset{k_{-1}}{\stackrel{{\scriptstyle
k_{1}}}{{\rightleftharpoons}}}ES$ $\displaystyle
ES\overset{k_{2}}{\rightarrow}E+P$ (11) $\displaystyle
ES+S\underset{k_{-3}}{\stackrel{{\scriptstyle
k_{3}}}{{\rightleftharpoons}}}ESS$

there are five reactive processes (five channels in the language of the
Gillespie algorithm) and consequently five families of counting processes
$\\{\chi_{i_{h}}^{(h)}(t,\cdot)\\}$, $h=1,\dots,5$, should be introduced, so
that the formulation of the discrete stochastic dynamics reads

$\displaystyle\frac{dN_{S}(t)}{dt}$ $\displaystyle=$
$\displaystyle-\sum_{i=1}^{N_{S}(t)}\frac{d\chi_{i}^{(1)}(t,\widetilde{k}_{1}\,N_{E}(t))}{dt}+\sum_{j=1}^{N_{ES}(t)}\frac{d\chi_{j}^{(2)}(t,k_{-1})}{dt}$
$\displaystyle\frac{dN_{E}(t)}{dt}$ $\displaystyle=$
$\displaystyle-\sum_{i=1}^{N_{S}(t)}\frac{d\chi_{i}^{(1)}(t,\widetilde{k}_{1}\,N_{E}(t))}{dt}+\sum_{j=1}^{N_{ES}(t)}\frac{d\chi_{j}^{(2)}(t,k_{-1})}{dt}+\sum_{h=1}^{N_{ES}(t)}\frac{d\chi_{h}^{(3)}(t,k_{2})}{dt}$
$\displaystyle\frac{dN_{ES}(t)}{dt}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N_{S}(t)}\frac{d\chi_{i}^{(1)}(t,\widetilde{k}_{1}\,N_{E}(t))}{dt}-\sum_{j=1}^{N_{ES}(t)}\frac{d\chi_{j}^{(2)}(t,k_{-1})}{dt}-\sum_{h=1}^{N_{ES}(t)}\frac{d\chi_{h}^{(3)}(t,k_{2})}{dt}-\sum_{k=1}^{N_{S}(t)}\frac{d\chi_{k}^{(4)}(t,\widetilde{k}_{3}\,N_{ES}(t))}{dt}$
(12) $\displaystyle+$
$\displaystyle\sum_{l=1}^{N_{ESS}(t)}\frac{d\chi_{l}^{(5)}(t,k_{-3})}{dt}$
$\displaystyle\frac{dN_{ESS}(t)}{dt}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{N_{S}(t)}\frac{d\chi_{k}^{(4)}(t,\widetilde{k}_{3}\,N_{ES}(t))}{dt}-\sum_{l=1}^{N_{ESS}(t)}\frac{d\chi_{l}^{(5)}(t,k_{-3})}{dt}$
$\displaystyle\frac{dN_{P}(t)}{dt}$ $\displaystyle=$
$\displaystyle\sum_{h=1}^{N_{ES}(t)}\frac{d\chi_{h}^{(3)}(t,k_{2})}{dt}$

equipped with the initial conditions $c_{S}(0)=c_{S,0}$, $c_{E}(0)=c_{E,0}$,
$c_{ES}(0)=c_{ESS}(0)=c_{P}(0)=0$. Observe that for the bimolecular steps we
have used a number-dependent rate coefficient. This is just one possibility,
out of other fully equivalent alternatives, of defining bimolecular reacting
processes, and out of tem a numerical algorithm for solving them. This issue,
and its computational implications will be addressed elsewhere PG_next . The
granularity number $N_{g}$ can be fixed, so that

$N_{S}(0)=[c_{S,0}\,N_{g}]\,,\qquad N_{E,0}=[c_{E,0}\,N_{g}]$ (13)

where $[\xi]$ indicates the integer part of $\xi$, thus defining the relation
betwen $N_{\alpha}(t)$ and $c_{\alpha}(t)$, namely
$c_{\alpha}(t)=N_{\alpha}(t)/N_{g}$, $\alpha=S,\,E,\,ES,\,ESS,\,P$. This
implies also that the effective rate parameters entering the discrete
stochastic evolution equation (12), and associated with the two bimolecular
reactive steps, are given by $\widetilde{k}_{1}=k_{1}/N_{g}$, and
$\widetilde{k}_{3}=k_{3}/N_{g}$.

Consider the case $k_{-1}=k_{2}=k_{3}=k_{-3}=1$, $c_{S,0}=4$, $c_{E,0}=0.1$.
In this case the quasi steady-state approximation of the $c_{ES}$-$c_{S}$
curve (representing the slow manifold of the kinetics takes the expression

$c_{ES}=\frac{c_{E,0}\,c_{S}}{K_{M}+c_{S}+\beta\,c_{S}^{2}}\,,\qquad
K_{M}=\frac{k_{-1}+k_{2}}{k_{1}}\,,\quad\beta=\frac{k_{-3}}{k_{3}}$ (14)

Figure 4 depicts the $c_{ES}$-$c_{S}$ graph obtained from a single realization
of the stochastic process eq. (11) at several values of $k_{1}$ so as to
modify the Michaelis-Menten constant $K_{M}$ for a value $N_{g}=10^{6}$ of the
granularity number.

Apart from the initial transient giving rise to an overshot in the values of
$c_{ES}$ near $c_{S}\simeq c_{S,0}$, the dynamics rapidly collapses towards
the slow manifold and the stochastic simulations at high $N_{g}$-value provide
a reliable description of the mean-field behavior starting from a single
stochastic realization.

Figure 4: $c_{ES}$ vs $c_{S}$ plot of the substrate-inhibited enzymatic
kinetics discussed in the main text. Symbols (in color) are the results of
stochastic simulations of a single realization of the process eq. (11),
(black) solid lines the graph of the quasi steady-state approximation. The
arrow indicates increasing values of $K_{M}$, i.e. decreasing values of
$k_{1}=20,\,6,\,2$.

To conclude, we want to point out some advantages and extensions of the
present approach:

  * •

it shows a direct analogy between chemical reaction kinetics, radiative
processes and stochastic formulation of open quantum systems, thus, paving the
way for a unified treatment of the interpaly between these phenomena, that is
particularly important in the field of photochemistry, and in the foundation
of statistical physics PGradia ; petruccione ;

  * •

it can be easily extended to semi-Markov transition. This is indeed the case
of the growth kinetics of eukaryotic microorganisms, the physiological state
of which can be parametrized with respect to internal (hidden) parameters such
as the age, the cytoplasmatic content, etc.;

  * •

it can be easily extended to include transport phenomena. In point of fact,
the occurrence of Markovian or semi-Markovian transitions in modeling chemical
kinetics is analogous to the transitions occurring in the direction of motion
(Poisson-Kac processes, Lévy flights, Extended Poisson-Kac processes) or in
the velocity (linearized Boltzmannian schemes) GPK ; EPK ; Levy .

  * •

it is closely related to the formulation of stochastic differential equations
for the thermalization of athermal system athermo1 , in which the classical
mesoscopic description of thermal fluctuations, using the increments of a
Wiener process, is replaced by a dynamic model involving the increments of a
counting process.

Due to the limitations of a Letter, all these issues will be addressed in
forthcoming works. But apart for these extensions and improvements, the
proposed formulation indicates that the stochastic theory of chemical
reactions can be built upon a simple and consistent mathematical formalism
describing the elementary reactive events as Markovian or semi-Markovian
counting processes semimarkovcounting , that perfectly fits with the
description of molecular non reactive events (molecular collisions), providing
an unifying stochastic formalism of elementary (classical and quantum)
molecular events.

## References

  * (1) P. L. Krapivsky, S. Redner, E. Ben-Naim, A Kinetic View to Statistical Physics, Cambridge University Press, Cambridge (2010).
  * (2) L. Boltzmann, Weitere Studien $\ddot{u}$ber das W$\ddot{a}$rmeglichgenicht unter Gas-molek$\ddot{u}$len, Sitzungsberichte Akademie der Wissenschaften 66 (1872) 275-370.
  * (3) G. B. Marin, G. S. Yablonsky, D. Constales, Kinetics of chemical reactions: decoding complexity, John Wiley & Sons, New York, (2019).
  * (4) O. Levenspiel, Chemical Reaction Engineering, J. Wiley & Sons (1998).
  * (5) Z. Wang, Z. Hou, H. Xin, Internal noise stochastic resonance of synthetic gene network, Chemical Physics Letters, 401 (1-3) (2005) 307-311.
  * (6) M. Perc, M. Gosak, and M. Marhl, From stochasticity to determinism in the collective dynamics of diffusively coupled cells, Chemical Physics Letters, 421 (1-3) (2006) 106–110.
  * (7) G. Lente, A binomial stochastic kinetic approach to the michaelis–menten mechanism, Chemical Physics Letters, 568 (2013) 167–169.
  * (8) D. A. McQuarrie, Stochastic approach to chemical kinetics, Journal of Applied Probability, 4 (3) (1967) 413-478.
  * (9) D. T. Gillespie, Stochastic simulation of chemical kinetics, Annual Review of Physical Chemistry 58 (1) (2007) 35-55.
  * (10) M. Delbrück, Statistical fluctuations in autocatalytic reactions, The Journal of Chemical Physics 8 (1) (1940) 120-124.
  * (11) A. F. Bartholomay, A stochastic approach to statistical kinetics with application to enzyme kinetics, Biochemistry 1 (2) (1962) 223-230.
  * (12) D. T. Gillespie, A rigorous derivation of the chemical master equation, Physica A 188 (1-3) (1992) 404-425.
  * (13) J. Keizer, On the necessity of using the master equation to describe the chemical reaction $X+A\rightleftharpoons B+X$, Chemical Physics Letters, 10 (4) (1971) 371–374.
  * (14) B. J. Gaynor, R. G. Gilbert, K. D. King, Solution of the master equation for unimolecular reactions, Chemical Physics Letters, 55 (1) (1978) 40-43.
  * (15) D. T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, Journal of Computational Physics 22 (4) (1976) 403-434.
  * (16) D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, The Journal of Physical Chemistry 81 (25) (1977) 2340-2361.
  * (17) M. A. Gibson, J. Bruck, Efficient exact stochastic simulation of chemical systems with many species and many channels, The Journal of Physical Chemistry A 104 (9) (2000) 1876-1889.
  * (18) L. Lok , R. Brent, Automatic generation of cellular reaction networks with Moleculizer, Nature Biotechnology 23 (2005) 131–36
  * (19) Y. Cao, H. Li,L. R. Petzold, Efficient formulation of the stochastic simulation algorithm for chemically reacting systems, The Journal of Chemical Physics 121 (2004) 4059–67.
  * (20) M. Rathinam, L. R. Petzold, Y. Cao, D. T. Gillespie, Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method, The Journal of Chemical Physics, 119 (24) (2003) 12784-12794.
  * (21) C. Yang, D. T. Gillespie, L. R. Petzold, Efficient step size selection for the tau-leaping simulation method, The Journal of Chemical Physics 124 (4) (2006) 044109.
  * (22) C. Yang, D. T. Gillespie, L. R. Petzold, Adaptive explicit-implicit tau-leaping method with automatic tau selection, The Journal of Chemical Physics 126 (22) (2007) 224101.
  * (23) D. C. Venerus and H. C. Öttinger, A modern Course in Transport Phenomena, Cambridge University Press, Cambridge (2018).
  * (24) K. Ito and H. P. McKean Jr., Diffusion Processes and their Sample Paths, Springer, Berlin (1974).
  * (25) F. Campillo, M. Chebbi, S. Toumi, Stochastic modeling for biotechnologies Anaerobic model AM2b, Revue Africaine de la de la Recherche en Informatique et Mathématiques Appliqués, INRIA 28 (2018 - 2019), Mathematics for Biology and the Environment 13-23.
  * (26) C. Pezzotti, M. Giona, Stochastic chemical reactions: from algorithmic approaches to stochastic differential models, in preparation (2022).
  * (27) C. Pezzotti, G. Procopio, A. Brasiello, M. Giona, Stochastic simulations of bioreactors in the presence of biomass heterogeneity and structured eukaryotic populations, in preparation (2022).
  * (28) R. Aris, ”On the dispersion of a solute in a fluid flowing through a tube, Proceedings of the Royal Society of London A (235) (1956) 67-77.
  * (29) S. Cerbelli, M. Giona, F. Garofalo, Quantifying dispersion of finite-sized particles in deterministic lateral displacement microflow separators through Brenner’s macrotransport paradigm, Microfluidics and nanofluidics 15 (2013) 431-449.
  * (30) C. Pezzotti and M. Giona, Particle-photon radiative interactions and thermalization, in preperation (2022).
  * (31) H.-P. Breuer, F. Petruccione, The Theory of Open Quantum Systems, Clarendon Press, Oxford (2002).
  * (32) M. Giona, A. Brasiello, S. Crescitelli, Stochastic foundations of undulatory transport phenomena: Generalized Poisson–Kac processes—Part I basic theory, Journal of Physics A (50) (2017) 335002.
  * (33) M. Giona, A. Cairoli, R. Klages, Extended Poisson-Kac theory: A unifying framework for stochastic processes with finite propagation velocity, Physical Review X (12) (2022) 021004.
  * (34) K.-I. Sato, Lévy processes and infinitely divisible distributions, Cambridge University Press, Cambridge (1999).
  * (35) K. Kanazawa, Statistical Mechanics for Athermal Fluctuation, Springer Nature, Singapore (2017).
  * (36) D. Cocco, M. Giona, Generalized Counting Processes in a Stochastic Environment. Mathematics 9 (2021) 25-73.