# Conservative Integrators for Many–body Problems

Andy T. S. Wan<EMAIL_ADDRESS>Alexander Bihlo<EMAIL_ADDRESS>Jean-
Christophe Nave<EMAIL_ADDRESS>Department of Mathematics and
Statistics, University of Northern British Columbia,
Prince George, BC, V2N 4Z9, Canada  Department of Mathematics and Statistics,
Memorial University of Newfoundland,
St. John’s, NL, A1C 5S7, Canada  Department of Mathematics and Statistics,
McGill University, Montréal, QC, H3A 0B9, Canada

###### Abstract

Conservative symmetric second–order one–step schemes are derived for dynamical
systems describing various many–body systems using the Discrete Multiplier
Method. This includes conservative schemes for the $n$-species Lotka–Volterra
system, the $n$-body problem with radially symmetric potential and the
$n$-point vortex models in the plane and on the sphere. In particular, we
recover Greenspan–Labudde’s conservative schemes for the $n$-body problem.
Numerical experiments are shown verifying the conservative property of the
schemes and second–order accuracy.

###### keywords:

dynamical systems , conserved quantity , first integral , conservative methods
, Discrete Multiplier Method , long-term stability , divided difference ,
many–body system , Lotka–Volterra equations , point vortex equations

##  1 Introduction

The general many–body problem is of great importance in the mathematical
sciences. Except for particular configurations such as in [1], the governing
equations of $n\geqslant 3$ bodies in the classical $n$-body problem cannot be
integrated analytically and in general one has to resort to numerical
simulations. Similar statements are also true for other nonlinear many–body
problems, which also require numerical integrations for large number of
bodies. Furthermore, the underlying equations of many–body problems typically
have rich geometric structures, such as simplified real-world phenomena
including the $n$-species Lotka–Volterra systems and the point vortex system
described in this paper. Examples of geometric structures include variational
formulations, the existence of first integrals and the invariance under
certain coordinate transformations. Note that throughout this paper, we use
the terms first integral, invariants of motion, and conserved quantities
interchangeably. In order to preserve these structures numerically, special
classes of numerical methods, called geometric numerical integrators [2, 3, 4]
are often employed for this purpose.

Within the field of geometric numerical integration, finding numerical schemes
that preserve an underlying Hamiltonian structure of a system of ordinary
differential equations (ODEs), so-called symplectic integrators, has been
historically of prime interest in recent decades. Indeed, early examples of
such numerical schemes date back to the 1950s [5], with other important early
contributions found in [6] and [7], and a first monograph written in the 1990s
by [8]. For more recent expositions, see for example in [2, 3, 4, 9].

While Hamiltonian systems are important in the mathematical sciences, there
are some important restrictions that limit the general applicability of the
Hamiltonian framework. Amongst the most important restrictions are systems
exhibiting dissipation and systems that cannot be brought easily into a
canonical Hamiltonian form, such as dynamical systems with an odd number of
degrees of freedom. While for the latter, the more general Poisson geometry
and associated Poisson integrators are available [3], these numerical schemes
are not as universally applicable as the symplectic schemes for canonical
Hamiltonian systems.

Symplectic integrators also do not preserve arbitrary first integrals of
differential equations111In fact, it is known from [10] that for Hamiltonian
systems without additional first integrals, if a symplectic integrator using
an uniform time step size is also energy–preserving, then it is an exact
integrator; that is, the discrete flow reproduces the exact flow up to a time
reparameterization.. While the Hamiltonian is nearly conserved over
exponentially long time periods [11] and linear and quadratic first integrals
can be preserved by some symplectic Runge–Kutta methods [3, 4], higher–order
polynomial conserved quantities or first integrals of arbitrary form are
generally not preserved with symplectic methods [7]. If one is interested in
the exact preservation of first integrals of arbitrary forms, then
conservative methods would need to be utilized, i.e. methods that numerically
preserved conserved quantities exactly up to machine precision. What sets
conservative methods apart from other geometric integrators is that they may
possess long-term stability over arbitrarily long time periods [12].

Conservative numerical schemes have also been extensively investigated in the
literature. In [13], the average vector field method was proposed which allows
the exact preservation of the energy of arbitrary form in a Hamiltonian
system. A generalization of this idea to higher–order schemes using
collocation is found in [14], with applications to energy preserving schemes
for Poisson systems discussed in [15]. A further class of general conservative
schemes is the discrete gradient method, originally proposed in [16]. This
method relies on expressing a first–order system of ODEs in a skew-symmetric
gradient form. While some dynamical systems, in particular Hamiltonian systems
and so-called Nambu systems [17], admit natural skew-symmetric gradient
representations, many other systems have to be first brought into a skew-
gradient form before the discrete gradient method can be applied. Moreover,
one main drawback when applying the discrete gradient method to large
dimensional systems with multiple invariants is that the order of the
associated skew-symmetric tensor increases with the number of conserved
quantities.

One further class of exactly conservative methods is given by projection
methods [3]. Here one applies a standard (usually explicit) integrator over
one or more time steps and subsequently projects the resulting numerical
approximation onto the manifold spanned by the conserved quantities. As
discussed in [12], the projection step can become problematic if the manifold
spanned by the invariants consists of several connected components, since the
projection may then bring the numerical solution onto the wrong connected
component.

In [18], we have introduced the Discrete Multiplier Method (DMM), which is a
general purpose method for finding conservative schemes for dynamical systems
with arbitrary forms of conserved quantities. The proposed method rests on
discretizing the characteristic [19], also called conservation law multiplier
[20], of conservation laws so that the discrete conserved quantities hold.
This idea was originally proposed in [21] for both PDE and ODE systems, with a
systematic framework for constructing conservative finite difference schemes
for general ODE systems derived in [18]. There it was also shown that the
average vector field method corresponds to a special choice of the
discretization of the conservation law multiplier of Hamiltonian systems.
Several examples of conservative schemes for classical dynamical systems were
presented in [18], but many–body systems were not considered there. As several
important dynamical systems are in fact many–body problems, the purpose of
this paper is to demonstrate that the constructive framework proposed in [18]
is also suitable for large dynamical systems.

We note here that for the purpose of the present paper, ‘many–body systems’
refers to dynamical systems with at most a few thousand degrees of freedom.
While the DMM for finding conservative integrators is not dependent on the
number of degrees of freedom of the underlying dynamical system, the resulting
conservative schemes are typically implicit. As such, a practical
implementation of these schemes generally relies on an (efficient) implicit
solver which renders the case of very many bodies (i.e. millions and more)
computationally challenging. We do not aim to address this computational
challenge in the present work, where we exclusively use a standard fixed point
iteration for solving these implicit conservative schemes, and reserve a more
in-depth study of this computational problem for future work.

The further organization of this paper is as follows. In Section 2, we give a
brief review of DMM for conservative discretizations as originally proposed in
[21]. We then propose several second–order conservative schemes derived using
DMM for many–body systems in the following sections. Specifically, Section 3.1
is devoted to a conservative schemes for the general $n$-species
Lotka–Volterra system of population dynamics. In Section 3.2, we present a
conservative scheme for the $n$-body problem with general radially symmetric
potential and recover Greenspan–Labudde’s conservative scheme [22, 23]. The
celestial $n$-body problem and the Lennard–Jones potential from molecular
dynamics are considered as special cases. Section 3.3 details conservative
schemes for the $n$-point vortex model on the plane and on the sphere. Section
4 features numerical results of the various schemes derived in this paper.
Finally, we make some concluding remarks in Section 5. A includes theoretical
verifications of conservative properties of all the schemes presented and
shows that they are all second–order accurate.

##  2 Construction of exactly conservative integrators

Before discussing the theory of DMM presented in [18], we first fix some
notations which will be used throughout this article.

###  2.1 Notations and conventions

Let $U\subset\mathbb{R}^{n}$ and $V\subset\mathbb{R}^{m}$ be open subsets
where here and in the following $n,m,p\in\mathbb{N}$. $f\in C^{p}(U\rightarrow
V)$ means $f$ is a $p$-times continuously differentiable function with domain
in $U$ and range in $V$. We often use boldface to indicate a vector quantity
$\boldsymbol{f}$. If $\boldsymbol{f}\in C^{1}(U\rightarrow V)$,
$\partial_{\boldsymbol{x}}\boldsymbol{f}:=\left[\frac{\partial f_{i}}{\partial
x_{j}}\right]$ denotes the Jacobian matrix. Let $I\subset\mathbb{R}$ be an
open interval and let $\boldsymbol{x}\in C^{1}(I\rightarrow U)$,
$\dot{\boldsymbol{x}}$ denotes the derivative with respect to time $t\in I$.
Also if $\boldsymbol{x}\in C^{p}(I\rightarrow U)$, $\boldsymbol{x}^{(q)}$
denotes the $q$-th time derivative of $\boldsymbol{x}$ for $1\leq q\leq p$.
For brevity, the explicit dependence of $\boldsymbol{x}$ on $t$ is often
omitted with the understanding that $\boldsymbol{x}$ is to be evaluated at
$t$. If $\boldsymbol{\psi}\in C^{1}(I\times U\rightarrow V)$,
$D_{t}\boldsymbol{\psi}$ denotes the total derivative with respect to $t$, and
$\partial_{t}\boldsymbol{\psi}$ denotes the partial derivative with respect to
$t$. $M_{m\times n}(\mathbb{R})$ denotes the set of all $m\times n$ matrices
with real entries.

###  2.2 Conserved quantities of quasilinear first–order ODEs

Consider a quasilinear first–order system of ODEs,

$\displaystyle\dot{\boldsymbol{x}}(t)$
$\displaystyle=\boldsymbol{f}(t,\boldsymbol{x}),$ (1)
$\displaystyle\boldsymbol{x}(t_{0})$ $\displaystyle=\boldsymbol{x}_{0}.$

where $t\in I$, $\boldsymbol{x}=(x_{1}(t),\dots,x_{n}(t))\in U$. For $1\leq
p\in\mathbb{N}$, if $\boldsymbol{f}\in C^{p-1}(I\times
U\rightarrow\mathbb{R}^{n})$ and is Lipschitz continuous in $U$, then standard
ODE theory implies there exists an unique solution $\boldsymbol{x}\in
C^{p}(I\rightarrow U)$ to the first–order system (1) in a neighborhood of
$(t_{0},\boldsymbol{x}_{0})\in I\times U$.

######  Definition 1.

Let $m\in\mathbb{N}$ with $1\leq m\leq n$. A vector-valued function
$\boldsymbol{\psi}\in C^{1}(I\times U\rightarrow\mathbb{R}^{m})$ is a vector
of conserved quantities222By quasilinearity of (1), it suffices to consider
conserved quantities depending only on $t,\boldsymbol{x}$, see [18]. (or
equivalently first integrals) if

$\displaystyle D_{t}\boldsymbol{\psi}(t,\boldsymbol{x})=\boldsymbol{0},\text{
for any $t\in I$ and $C^{1}(I\rightarrow U)$ solution }\boldsymbol{x}\text{ of
}\eqref{odeEqn}.$ (2)

In other words, $\boldsymbol{\psi}(t,\boldsymbol{x})$ is constant on any
$C^{1}(I\rightarrow U)$ solution $\boldsymbol{x}$ of (1).

A generalization of integrating factors is known as _characteristics_ by [19]
or equivalently, _conservation law multipliers_ by [20]. We will adopt the
terminology of conversation law multiplier or just multiplier when the context
is clear.

######  Definition 2.

Let $m\in\mathbb{N}$ with $1\leq m\leq n$ and $U^{(1)}$ be an open subset of
$\mathbb{R}^{n}$. A conservation law multiplier of $\boldsymbol{F}$ is a
matrix-valued function $\Lambda\in C(I\times U\times U^{(1)}\rightarrow
M_{m\times n}(\mathbb{R}))$ such that there exists a function
$\boldsymbol{\psi}\in C^{1}(I\times U\rightarrow\mathbb{R})$ satisfying,

$\displaystyle\Lambda(t,\boldsymbol{x},\dot{\boldsymbol{x}})(\dot{\boldsymbol{x}}(t)-\boldsymbol{f}(t,\boldsymbol{x}))=D_{t}\boldsymbol{\psi}(t,\boldsymbol{x}),\text{
for $t\in I$, }\boldsymbol{x}\in C^{1}(I\rightarrow U).$ (3)

Here, we emphasize that condition (3) is satisfied as an identity for
arbitrary $C^{1}$ functions $\boldsymbol{x}$; in particular $\boldsymbol{x}$
need not be a solution of (1). It follows from the definition of conservation
law multiplier that existence of multipliers implies existence of conservation
laws. Conversely, given a known vector of conserved quantities
$\boldsymbol{\psi}$, there can be many conservation law multipliers which
correspond to $\boldsymbol{\psi}$. It was shown in [18] that it suffices to
consider multipliers of the form $\Lambda(t,\boldsymbol{x})$ where a one-to-
one correspondence exists between conservation law multipliers and conserved
quantities of (1).

######  Theorem 1 (Theorem 4 of [18]).

Let $\boldsymbol{\psi}\in C^{1}(I\times U\rightarrow\mathbb{R}^{m})$. Then
there exists a unique conservation law multiplier of (1) of the form
$\Lambda\in C(I\times U\rightarrow M_{m\times n}(\mathbb{R}))$ associated with
the function $\boldsymbol{\psi}$ if and only if $\boldsymbol{\psi}$ is a
conserved quantity of (1). And if so, $\Lambda$ is unique and satisfies for
any $t\in I$ and $\boldsymbol{x}\in C^{1}(I\rightarrow U)$,

$\displaystyle\Lambda(t,\boldsymbol{x})=\partial_{\boldsymbol{x}}\boldsymbol{\psi}(t,\boldsymbol{x}),$
(4a)
$\displaystyle\Lambda(t,\boldsymbol{x})\boldsymbol{f}(t,\boldsymbol{x})=-\partial_{t}\boldsymbol{\psi}(t,\boldsymbol{x}).$
(4b)

To construct conservative methods for (1) with conserved quantities (2), we
shall discretize the time interval $I$ by a uniform time size
$\tau\in\mathbb{R}$, i.e. $t^{k+1}=t^{k}+\tau$ for $k\in\mathbb{N}$, and focus
on one–step conservative methods333Analogous results hold for variable time
step sizes and multi-step methods, see [18] for more details.. First, we
recall some definitions from [18].

######  Definition 3.

Let $W$ be a normed vector space, such as $\mathbb{R}^{m}$ with the Euclidean
norm or $M_{m\times n}(\mathbb{R})$ with the operator norm. A function
$g^{\tau}:I\times U\times U\rightarrow W$ is called a one–step function if
$g^{\tau}$ depends only on $t^{k}\in I$ and the discrete approximations
$\boldsymbol{x}^{k+1},\boldsymbol{x}^{k}\in U$.

######  Definition 4.

A sufficiently smooth one–step function $g^{\tau}:I\times U\times U\rightarrow
W$ is consistent to a sufficiently smooth $g:I\times U\times
U^{(1)}\rightarrow W$ if for any $\boldsymbol{x}\in C^{2}(I\rightarrow U)$,
there is a constant $C>0$ independent of $\tau$ so that $\left\lVert
g(t^{k},\boldsymbol{x}(t^{k}),\dot{\boldsymbol{x}}(t^{k}))-g^{\tau}(t^{k},\boldsymbol{x}(t^{k+1}),\boldsymbol{x}(t^{k}))\right\rVert_{W}\leq
C\left\lVert\boldsymbol{x}\right\rVert_{C^{2}([t^{k},t^{k+1}])}\tau,$ where
$\left\lVert\boldsymbol{x}\right\rVert_{C^{2}([t^{k},t^{k+1}])}:=\displaystyle\max_{0\leq
i\leq
2}\left\lVert\boldsymbol{x}^{(i)}\right\rVert_{L^{\infty}([t^{k},t^{k+1}])}$.
If so, we write $g^{\tau}=g+\mathcal{O}(\tau)$.

We shall be considering the following consistent one–step functions for
$\dot{\boldsymbol{x}},D_{t}\boldsymbol{\psi},\partial_{t}\boldsymbol{\psi}$:

$\displaystyle
D_{t}^{\tau}\boldsymbol{x}(t^{k},\boldsymbol{x}^{k+1},\boldsymbol{x}^{k})$
$\displaystyle:=\frac{\boldsymbol{x}^{k+1}-\boldsymbol{x}^{k}}{\tau}=\dot{x}+\mathcal{O}(\tau),$
(5) $\displaystyle
D_{t}^{\tau}\boldsymbol{\psi}(t^{k},\boldsymbol{x}^{k+1},\boldsymbol{x}^{k})$
$\displaystyle:=\frac{\boldsymbol{\psi}(t^{k+1},\boldsymbol{x}^{k+1})-\boldsymbol{\psi}(t^{k},\boldsymbol{x}^{k})}{\tau}=D_{t}\boldsymbol{\psi}+\mathcal{O}(\tau),$
(6)
$\displaystyle\partial_{t}^{\tau}\boldsymbol{\psi}(t^{k},\boldsymbol{x}^{k+1},\boldsymbol{x}^{k})$
$\displaystyle:=\frac{\boldsymbol{\psi}(t^{k+1},\boldsymbol{x}^{k})-\boldsymbol{\psi}(t^{k},\boldsymbol{x}^{k})}{\tau}=\partial_{t}\boldsymbol{\psi}+\mathcal{O}(\tau).$
(7)

######  Definition 5.

Let $\boldsymbol{f}^{\tau}$ be a consistent one–step function to
$\boldsymbol{f}$. We say that the one–step method,

$D_{t}^{\tau}\boldsymbol{x}(t^{k},\boldsymbol{x}^{k+1},\boldsymbol{x}^{k})=\boldsymbol{f}^{\tau}(t^{k},\boldsymbol{x}^{k+1},\boldsymbol{x}^{k})$
(8)

is conservative in $\boldsymbol{\psi}$, if
$\boldsymbol{\psi}(t^{k+1},\boldsymbol{x}^{k+1})=\boldsymbol{\psi}(t^{k},\boldsymbol{x}^{k})$
on any solution $\boldsymbol{x}^{k+1}$ of (8) and $k\in\mathbb{N}$.

We now state two key conditions from [18] for constructing conservative
one–step methods, which can be seen as a discrete analogue of (4a) and (4b).

######  Theorem 2 (Theorem 17 of [18]).

Let
$D_{t}^{\tau}\boldsymbol{x},D_{t}^{\tau}\boldsymbol{\psi},\partial_{t}^{\tau}\boldsymbol{\psi}$
be as defined in (5)–(7), and let $\Lambda$ be the conservation law multiplier
of (1) associated with a conserved quantity $\boldsymbol{\psi}$. If
$\boldsymbol{f}^{\tau}$ and $\Lambda^{\tau}$ are consistent one–step functions
to $\boldsymbol{f},\Lambda$ satisfying

$\displaystyle\Lambda^{\tau}D^{\tau}_{t}\boldsymbol{x}$
$\displaystyle=D^{\tau}_{t}\boldsymbol{\psi}-\partial^{\tau}_{t}\boldsymbol{\psi},$
(9a) $\displaystyle\Lambda^{\tau}\boldsymbol{f}^{\tau}$
$\displaystyle=-\partial^{\tau}_{t}\boldsymbol{\psi},$ (9b)

then the one–step method defined by (8) is conservative in
$\boldsymbol{\psi}$.

In [18], condition (9a) was solved by the use of divided difference calculus
and (9b) was solved using a local matrix inversion formula. In the following
many–body problems, we shall directly verify (9a) and (9b) for specific
choices of $\boldsymbol{f}^{\tau}$ and $\Lambda^{\tau}$.

Lastly, we recall a well-known result for even order of accuracy of symmetric
schemes. For more details, see Chapter II.3 of [3].

######  Definition 6 (Symmetric Schemes [3]).

Let $\Phi^{\tau}$ be the discrete flow of a one–step numerical method for
system (1) with time step $\tau$. The associated adjoint method
$(\Phi^{\tau})^{*}$ of the one–step method $\Phi^{\tau}$ is the inverse of the
original method with reversed time step $-\tau$, i.e.
$(\Phi^{\tau})^{*}=(\Phi^{-\tau})^{-1}$. A method is symmetric if
$(\Phi^{\tau})^{*}=\Phi^{\tau}$.

######  Theorem 3 (Theorem II-3.2 of [3]).

A symmetric method is of even order.

In order words, combining with Theorem 2, the one–step conservative schemes
are at least second–order accurate if they are symmetric.

##  3 Examples of DMM for many–body systems

In the section, we review some examples of many–body systems from population
dynamics, classical mechanics, molecular dynamics and fluid dynamics.
Specifically, we present conservative schemes derived using DMM for the
$n$-species Lotka–Volterra systems, the $n$-body problem involving the
gravitational potential and the Lennard–Jones potential, and the $n$-point
vortex problem on the plane and on the unit sphere. For brevity and clarity,
we have included all the details of calculations for derivations and
verifications in A.

###  3.1 Conservative schemes for Lotka–Volterra systems

The $n$-species Lotka–Volterra equations describe a simplified dynamics among
$n$ competing species interacting in an environment [24]. Specifically, we
consider the $n$-species Lotka–Volterra system given in the form of,

$\displaystyle\boldsymbol{F}(\boldsymbol{x},\dot{\boldsymbol{x}}):=\begin{bmatrix}\dot{x}_{i}-x_{i}\sum_{j=1}^{n}a_{ij}(x_{j}-\xi_{j})\end{bmatrix}_{1\leq
i\leq n}={\boldsymbol{0}},$ (10)

where $\boldsymbol{x}=(x_{1},\dots,x_{n})^{T}$ is the population of each
species with $x_{i}>0$, $A=[a_{ij}]$ is an $n\times n$ interaction matrix with
real entries and $\boldsymbol{\xi}=(\xi_{1},\dots,\xi_{n})^{T}$ is a fixed
point of the system. It is known from [25] that (10) has the conserved
quantity

$V(\boldsymbol{x}):=\sum_{i=1}^{n}d_{i}(\xi_{i}\log x_{i}-x_{i}),$ (11)

if there exists an $n\times n$ real diagonal matrix
$D=\text{diag}(d_{1},\dots,d_{n})$ such that $DA$ is skew-symmetric.

A conservative scheme for (10) which preserves $V$ numerically was derived
using DMM and is given by,

$\displaystyle\boxed{\boldsymbol{F}^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{x}^{k}):=\begin{bmatrix}\dfrac{x_{i}^{k+1}-x_{i}^{k}}{\tau}-x_{i}^{\tau}\sum\limits_{j=1}^{n}a_{ij}x_{j}^{\tau}\left(1-\dfrac{\xi_{j}}{x_{j}^{k}}g\left(\dfrac{x_{j}^{k+1}}{x_{j}^{k}}\right)\right)\end{bmatrix}_{1\leq
i\leq n}=\boldsymbol{0},}$ (12)

where $g(z):=\dfrac{\log z}{z-1}$ and $x_{i}^{\tau}$ is any consistent
discretization of $x_{i}$ (i.e. $x_{i}^{\tau}\rightarrow x_{i}$, as
$\tau\rightarrow 0$). It is interesting to note that (12) is consistent to
(10) since $g(z)=1-\frac{1}{2}(z-1)+\mathcal{O}((z-1)^{2})\rightarrow 1$ as
$\dfrac{x_{j}^{k+1}}{x_{j}^{k}}\rightarrow 1$ with $\tau\rightarrow 0$.

The simplest consistent choices of $x_{i}^{\tau}$ would be
$x_{i}^{\tau}:=x_{i}^{k}$ or $x_{i}^{k+1}$, which will lead to first–order
conservative schemes. Thus, for improved accuracy with similar computational
costs, we choose $x_{i}^{\tau}$ so that the resulting scheme (12) is
symmetric, which will lead to a second–order scheme according to Theorem 3. In
particular, we show in A.1 that (12) is conservative and is symmetric if
$x_{i}^{\tau}$ itself is symmetric. Specifically, choosing
$x_{i}^{\tau}=\overline{x_{i}}:=\dfrac{1}{2}\left(x_{i}^{k}+x_{i}^{k+1}\right)$
leads to an “Arithmetic mean DMM” scheme for (10). Since the phase variables
$x_{i}$ for (10) are nonnegative, another choice is
$x_{i}^{\tau}:=\sqrt{x_{i}^{k}x_{i}^{k+1}}$ or a “Geometric mean DMM” scheme
for (10).

###  3.2 Many–body problem with pairwise radial potentials

The many–body problem with pairwise radial potential, i.e. with conservative
forces that only depend on the radial difference of each two point masses, is
one of the most fundamental models in classical mechanics. It describes, in an
idealized fashion, numerous physical phenomena, including the motion of
planets in a solar system and atoms of molecules.

We consider the many–body problem as the Hamiltonian system with $n$ particles
in $\mathbb{R}^{3}$ with radial interaction potentials444There are in general
$6n$ unknowns for (13). In particular, the two body problem is integrable
since there are 10 constants of motion with an additional two conserved
quantities provided by the Laplace–Runge–Lenz vector.,

$\boldsymbol{F}(\boldsymbol{q},\boldsymbol{p},\dot{\boldsymbol{q}},\dot{\boldsymbol{p}}):=\begin{pmatrix}\begin{bmatrix}\dot{\boldsymbol{q}}_{i}-\dfrac{\boldsymbol{p}_{i}}{m_{i}}\end{bmatrix}_{1\leq
i\leq n}\\\ \begin{bmatrix}\dot{\boldsymbol{p}}_{i}+\sum\limits_{j=1,j\neq
i}^{n}\dfrac{\partial V_{ij}}{\partial
q_{ij}}(q_{ij})\dfrac{\boldsymbol{q}_{i}-\boldsymbol{q}_{j}}{q_{ij}}\end{bmatrix}_{1\leq
i\leq n}\end{pmatrix}=\boldsymbol{0},$ (13)

where
$\boldsymbol{q}=(\boldsymbol{q}_{1},\dots,\boldsymbol{q}_{n})^{T},\boldsymbol{p}=(\boldsymbol{p}_{1},\dots,\boldsymbol{p}_{n})^{T}$
with
$\boldsymbol{q}_{i}\in\mathbb{R}^{3},\boldsymbol{p}_{i}\in\mathbb{R}^{3},m_{i}\in\mathbb{R}_{+}$
as the position, momenta and mass of the $i$-th particle. For each distinct
pair $(i,j)$ of particles,
$q_{ij}:=|\boldsymbol{q}_{i}-\boldsymbol{q}_{j}|=q_{ji}$ denotes their
Euclidean distance and $V_{ij}$ is their radial pairwise potential energy such
that $V_{ij}=V_{ji}$. From classical mechanics, it is well-known that there
are ten constants of motion for (13). Specifically, they are the Hamiltonian
$H$, total linear momentum $\boldsymbol{P}$, total angular momentum
$\boldsymbol{L}$ and initial center of mass $\boldsymbol{C}$ – given by,

$\displaystyle\begin{split}H(\boldsymbol{q},\boldsymbol{p})&:=\sum_{i=1}^{n}\dfrac{\boldsymbol{p}_{i}^{T}\boldsymbol{p}_{i}}{2m_{i}}+\sum_{1\leq
i<j\leq n}V_{ij}(|\boldsymbol{q}_{i}-\boldsymbol{q}_{j}|)\\\
\boldsymbol{P}(\boldsymbol{q},\boldsymbol{p})&:=\sum_{i=1}^{n}\boldsymbol{p}_{i}\\\
\end{split}\begin{split}\boldsymbol{L}(\boldsymbol{q},\boldsymbol{p})&:=\sum_{i=1}^{n}\boldsymbol{q}_{i}\times\boldsymbol{p}_{i}\\\
\boldsymbol{C}(t,\boldsymbol{q},\boldsymbol{p})&:=\dfrac{1}{M}\left(\sum_{i=1}^{n}m_{i}\boldsymbol{q}_{i}\right)-\dfrac{\boldsymbol{P}}{M}t,\end{split}$
(14)

where $M=\sum_{i=1}^{n}m_{i}$ is the total mass of the system.

A conservative scheme for (13) which preserves all ten first integrals was
derived using DMM and is given by,

$\displaystyle\boxed{\boldsymbol{F}^{\tau}(\boldsymbol{q}^{k+1},\boldsymbol{p}^{k+1},\boldsymbol{q}^{k},\boldsymbol{p}^{k}):=\begin{pmatrix}\begin{bmatrix}\dfrac{\boldsymbol{q}^{k+1}-\boldsymbol{q}^{k}}{\tau}-\dfrac{\overline{\boldsymbol{p}_{i}}}{m_{i}}\end{bmatrix}_{1\leq
i\leq n}\\\
\begin{bmatrix}\dfrac{\boldsymbol{p}^{k+1}-\boldsymbol{p}^{k}}{\tau}+\sum\limits_{j=1,j\neq
i}^{n}\dfrac{\Delta V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})\end{bmatrix}_{1\leq
i\leq n}\end{pmatrix}=\boldsymbol{0}.}$ (15)

We note that the discretization (15) was previously reported in [23, Section
2.2], although there no constructive derivation was given. It is thus the
added benefit of DMM that conservative schemes can be constructed
systematically as detailed in the Appendix.

Here, scalar or vector quantities with a line above denote the arithmetic mean
of the quantity at $t_{k}$ and $t_{k+1}$. Similar to the previous
Lotka–Volterra example, the reason for this symmetric choice is due to Theorem
3. As shown in the Appendix, (15) is conservative and symmetric provided
$\frac{\Delta V_{ij}}{\Delta q_{ij}}$ is symmetric.

We now present some particular cases of the potential $V$ that specify the
discretizations (15) for physically relevant problems in celestial mechanics
and molecular dynamics.

####  3.2.1 Gravitational Potential

The classical Newtonian form of the many–body problem as applies to the solar
system is given by the following specification of the pairwise radial
potential,

$\displaystyle V_{ij}(q_{ij})=-\dfrac{Gm_{i}m_{j}}{q_{ij}},$

where $G$ is Newton’s gravitational constant. For the numerical scheme (15),
the respective divided difference for the potential thus simplifies to

$\displaystyle\boxed{\frac{\Delta V_{ij}}{\Delta
q_{ij}}=\frac{Gm_{i}m_{j}}{q_{ij}^{k+1}q_{ij}^{k}}},$

which is symmetric under the permutation $k\leftrightarrow k+1$.

####  3.2.2 Lennard–Jones Potential

In classical molecular dynamics, forces are modeled to be attractive when
particles are far from each other and repulsive when they are close. A
classical example for a potential in molecular dynamics is the Lennard–Jones
potential, given by

$\displaystyle
V_{ij}(q_{ij})=4\epsilon\left(\dfrac{\sigma^{12}}{q_{ij}^{12}}-\dfrac{\sigma^{6}}{q_{ij}^{6}}\right),$

where $\epsilon$ and $\sigma$ are the potential well depth and the distance
where the potential becomes zero, respectively. For this particular form of
the potential, the divided difference for the potential in (15) becomes

$\displaystyle\boxed{\frac{\Delta V_{ij}}{\Delta
q_{ij}}=4\epsilon\left(-\frac{\sigma^{12}}{q_{ij}^{k+1}q_{ij}^{k}}\sum_{l=0}^{11}(q_{ij}^{k+1})^{l-11}(q_{ij}^{k})^{-l}+\frac{\sigma^{6}}{q_{ij}^{k+1}q_{ij}^{k}}\sum_{l=0}^{5}(q_{ij}^{k+1})^{l-5}(q_{ij}^{k})^{-l}\right)},$
(16)

which is again symmetric under the permutation $k\leftrightarrow k+1$.

###  3.3 Point vortex problem

The simplified modeling of the continuous description of fluid mechanics on
the plane by an ensemble of point vortices dates back more than 150 years, and
was first introduced by Helmholtz [26], with other important early
contributions due to Kirchhoff who in particular first established the
Hamiltonian representation of the equations of point vortex dynamics [27,
Lecture 20], see also [28, 29] for more detailed reviews.

We first consider the classical $n$-point vortex problem in the plane before
discussing the $n$-point vortex problem on the unit sphere. While in the
planar case the point vortex equations are a canonical Hamiltonian system, in
the spherical case the equations constitute a non-canonical Hamiltonian
system. From the point of view of geometric numerical integration, standard
symplectic integrators can be used in the planar case, see results in [30, 31,
32], while the spherical case requires the use of Poisson integrators, with
some recent examples found in [33, 34]. Here, we show the DMM framework
provides exactly conservative integrators for both cases.

####  3.3.1 Planar case

Consider the $n$-point vortex problem on the plane,

$\displaystyle\boldsymbol{F}(\boldsymbol{x},\boldsymbol{y},\dot{\boldsymbol{x}},\dot{\boldsymbol{y}}):=\begin{pmatrix}\begin{bmatrix}\dot{x}_{i}+\dfrac{1}{2\pi}\sum\limits_{j=1,j\neq
i}^{n}\Gamma_{j}\dfrac{y_{ij}}{r_{ij}^{2}}\end{bmatrix}_{1\leq i\leq n}\\\
\begin{bmatrix}\dot{y}_{i}-\dfrac{1}{2\pi}\sum\limits_{j=1,j\neq
i}^{n}\Gamma_{j}\dfrac{x_{ij}}{r_{ij}^{2}}\end{bmatrix}_{1\leq i\leq
n}\end{pmatrix}=\boldsymbol{0},$ (17)

where $\boldsymbol{x}=(x_{1},\dots,x_{n})^{T}$ and
$\boldsymbol{y}=(y_{1},\dots,y_{n})^{T}$ with $(x_{i},y_{i})$ being the
position of the $i$-th point vortex on the plane, and $\Gamma_{i}$ is the
vorticity strength of the $i$-th vortex. We abbreviate $x_{ij}:=x_{i}-x_{j}$
and $y_{ij}:=y_{i}-y_{j}$ and $r_{ij}=\sqrt{x_{ij}^{2}+y_{ij}^{2}}$. It is
well-known that (17) possesses four conserved quantities – linear momentum
$\boldsymbol{P}$, angular momentum $L$ and the Hamiltonian $H$ – given by,

$\displaystyle\begin{split}\boldsymbol{P}(\boldsymbol{x},\boldsymbol{y})&:=\begin{pmatrix}\sum\limits_{i=1}^{n}\Gamma_{i}x_{i}\\\
\sum\limits_{i=1}^{n}\Gamma_{i}y_{i}\end{pmatrix},\end{split}\begin{split}L(\boldsymbol{x},\boldsymbol{y})&:=\sum\limits_{i=1}^{n}\Gamma_{i}(x_{i}^{2}+y_{i}^{2}),\\\
H(\boldsymbol{x},\boldsymbol{y})&:=-\dfrac{1}{2\pi}\sum\limits_{1\leq i<j\leq
n}\Gamma_{i}\Gamma_{j}\log r_{ij}.\end{split}$ (18)

A conservative scheme that preserves these four conserved quantities was found
using DMM and is given by,

$\displaystyle\boxed{\boldsymbol{F}^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{y}^{k+1},\boldsymbol{x}^{k},\boldsymbol{y}^{k}):=\begin{pmatrix}\begin{bmatrix}\dfrac{x_{i}^{k+1}-x_{i}^{k}}{\tau}+\dfrac{1}{2\pi}\sum\limits_{j=1,j\neq
i}^{n}\dfrac{\Gamma_{j}\overline{y_{ij}}}{(r_{ij}^{k})^{2}}g\left(\left(\frac{r_{ij}^{k+1}}{r_{ij}^{k}}\right)^{2}\right)\end{bmatrix}_{1\leq
i\leq n}\\\
\begin{bmatrix}\dfrac{y_{i}^{k+1}-y_{i}^{k}}{\tau}-\dfrac{1}{2\pi}\sum\limits_{j=1,j\neq
i}^{n}\dfrac{\Gamma_{j}\overline{x_{ij}}}{(r_{ij}^{k})^{2}}g\left(\left(\frac{r_{ij}^{k+1}}{r_{ij}^{k}}\right)^{2}\right)\end{bmatrix}_{1\leq
i\leq n}\end{pmatrix}=\boldsymbol{0},}$ (19)

where $g(z)=\dfrac{\log z}{z-1}$ as in the Lotka–Volterra example. We show in
the Appendix that (19) is conservative and symmetric.

####  3.3.2 Spherical case

The $n$-point vortex problem on the unit sphere is governed by the equations

$\displaystyle\boldsymbol{F}(\boldsymbol{x},\dot{\boldsymbol{x}}):=\dot{\boldsymbol{x}}_{i}-\dfrac{1}{4\pi}\sum\limits_{j=1,j\neq
i}^{n}\Gamma_{j}\dfrac{{\boldsymbol{x}}_{j}\times{\boldsymbol{x}}_{i}}{1-{\boldsymbol{x}}_{i}\cdot{\boldsymbol{x}_{j}}}=\boldsymbol{0},$
(20)

where $\boldsymbol{x}=({\boldsymbol{x}}_{1},\dots,{\boldsymbol{x}}_{n})^{T}$
with ${\boldsymbol{x}}_{i}$ being the position of the $i$-th point vortex on
the sphere and $\Gamma_{i}$ being the vortex strength of the $i$-th vortex.
The point vortex equations on the unit sphere (20) possess four conserved
quantities, given by the Noether momentum $\boldsymbol{P}$ and the Hamiltonian
$H$, which are

$\displaystyle\begin{split}\boldsymbol{P}(\boldsymbol{x})&:=\sum\limits_{i=1}^{n}\Gamma_{i}{\boldsymbol{x}}_{i},\end{split}\begin{split}H(\boldsymbol{x})&:=-\dfrac{1}{4\pi}\sum\limits_{1\leq
i<j\leq
n}\Gamma_{i}\Gamma_{j}\log(2-2{\boldsymbol{x}}_{i}\cdot{\boldsymbol{x}}_{j}).\end{split}$
(21)

The conservative discretization for (20) was found using DMM and is given by

$\displaystyle\boxed{\boldsymbol{F}^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{x}^{k}):=\dfrac{{\boldsymbol{x}}_{i}^{k+1}-{\boldsymbol{x}}_{i}^{k}}{\tau}-\dfrac{1}{4\pi}\sum\limits_{1\leq
j\leq n,j\neq
i}\Gamma_{j}\dfrac{\overline{{\boldsymbol{x}}_{j}}\times\overline{{\boldsymbol{x}}_{i}}}{1-{{\boldsymbol{x}}_{i}^{k}}\cdot{{\boldsymbol{x}}_{j}^{k}}}g\left(\frac{1-{\boldsymbol{x}}_{i}^{k+1}\cdot{\boldsymbol{x}}_{j}^{k+1}}{1-{\boldsymbol{x}}_{i}^{k}\cdot{\boldsymbol{x}}_{j}^{k}}\right)=\boldsymbol{0},}$
(22)

where $g(z)=\dfrac{\log z}{z-1}$, as in the planar case. In the Appendix, (22)
is also shown to be conservative and symmetric.

##  4 Numerical results

In this section, we present numerical results for the conservative schemes
derived using DMM. We compare numerically the conservative property and
second–order accuracy with some classical and symplectic methods.
Specifically, as the derived DMM schemes are second–order implicit schemes, we
compare all examples with the implicit Midpoint method, which has similar
computational cost and is a second–order symplectic method with favourable
long-term properties [3]. Moreover, we also compare with some explicit methods
when applicable, such as the standard fourth-order explicit Runge–Kutta method
and the Störmer–Verlet method, which is a second–order symplectic method.

In the following, all implicit methods were solved by a fixed point iteration
using the right hand side of the respective schemes with an absolute tolerance
of $10^{-14}$, unless otherwise stated. For a final time $T$, we have used a
uniform time step of size $\tau$ with a total number of $N$ time steps. For a
conserved quantity $\psi$, we denote the $\ell^{\infty}$ norm in discrete time
of the error by

$\displaystyle\text{Error}[\psi(t,\boldsymbol{x})]:=\max_{k=1,\dots,N}|\psi(t^{k},\boldsymbol{x}^{k})-\psi(0,\boldsymbol{x}^{0})|.$

###  4.1 Lotka–Volterra systems

We begin with the Lotka–Volterra systems of Section 3.1. For reference, we
compare the arithmetic and geometric mean DMM schemes of (12) with the
Midpoint method and the standard explicit fourth-order Runge–Kutta method.

Specifically, we consider the three–species non-degenerate Lotka–Volterra
system represented by the $3\times 3$ interaction matrix

$\displaystyle A=\left(\begin{array}[]{ccc}1&1&1\\\ 0&0&-2\\\
0&1&0\end{array}\right),$

and a fixed point of
${\bf\xi}=\left(\frac{1}{2},\frac{1}{2},\frac{1}{2}\right)^{\rm T}$. This
system is non-degenerate since $\det\left(A\right)\neq 0.$ Additionally, since
$DA+A^{\rm T}D=0$ for $D=\text{diag}\left(0,1,2\right)$, it possesses _one_
conserved quantity $V\left(\boldsymbol{x}\right)=\frac{1}{2}\log
y-y+2\left(\frac{1}{2}\log z-z\right)$, as defined in Section 3.1. This is
precisely the conserved quantity that DMM was constructed to preserve exactly.

We choose the initial conditions of
${\boldsymbol{x}_{0}}=\left(\frac{1}{10},\frac{1}{10},\frac{1}{10}\right)^{\rm
T}$ with a final time of $T=50$ and $\tau=5\cdot 10^{-2}$. For the implicit
methods considered, we have applied one step of the forward Euler method as
the initial guess for the fixed point iterations with an absolute tolerance of
$10^{-15}$.

The convergence plot in Figure 1(a) confirms that all methods has indeed the
expected accuracy orders. In particular, the Arithmetic mean and Geometric
mean DMM are confirmed to be second order. Figure 1(b) shows that the error in
time of the conserved quantity $V\left(\boldsymbol{x}\right)$ for all four
methods considered. We see the exact conservation (up to machine precision)
for the two DMM schemes, in contrast to the Midpoint method and the fourth-
order Runge–Kutta method. As expected, we only observe a slow growth of error
in $V(\boldsymbol{x})$ for the DMM schemes due to accumulation of round-off
errors and the nonzero tolerance imposed by the fixed point iterations. Table
1 summarizes the $\ell^{\infty}$ error in $V(\boldsymbol{x})$ for all four
methods.

Method | Midpoint | RK4 | Arith. Mean DMM | Geo. Mean DMM
---|---|---|---|---
$\text{Error}[V(\boldsymbol{x})]$ | $1.76\cdot 10^{-3}$ | $3.68\cdot 10^{-6}$ | $1.78\cdot 10^{-15}$ | $6.22\cdot 10^{-15}$
Table 1: Error in conserved quantity $V\left(\boldsymbol{x}\right)$ for the
non-degenerate three–species Lotka–Volterra model.

(a) Convergence plot for all methods.

(b) Error in conserved quantity $V$ versus time for all methods.

Figure 2 shows the trajectories along the isosurface of
$\mathcal{V}=\left\\{\boldsymbol{x}\in\mathbb{R}^{3}\left|V\left(\boldsymbol{x}\right)=V\left(\boldsymbol{x}_{0}\right)\right.\right\\}.$
Qualitatively, we see that throughout a transient phase, all trajectories
remain close to the isosurface $\mathcal{V}$. The dynamics then settles to a
limit-cycle on the $y$–$z$ plane. While all trajectories are visually close to
$\mathcal{V}$, only the trajectories computed by the DMM schemes are within
machine-precision from $\mathcal{V}$. In contrast, trajectories of the
Midpoint and fourth-order Runge–Kutta method are respectively $\approx
10^{-3}$ and $\approx 10^{-6}$ away from $\mathcal{V}$, even though it is
difficult to discern visually from Figure 2.

Figure 2: Trajectories of all methods for the three–species Lotka–Volterra
system along $\mathcal{V}.$

###  4.2 Many–body problem: Gravitational potential

Next we consider the standard ten–body solar system model. We compare the
Greenspan–Labudde scheme, or equivalently the derived DMM scheme, with the
Störmer–Verlet and Midpoint schemes, which are two popular symplectic methods
for many–body problems. The initial conditions were obtained from planetary
data in [35]. We simulate the system to a final time $T=2\cdot 10^{7}$ days
($\sim 55,000$ years) using $N=4\cdot 10^{6}$ time steps, corresponding to a
uniform time step size of $\tau=5$ days. For the initial guess of the fixed
point iterations, we have used a perturbation of the discrete solution from
the previous time step in order to avoid singularities that can arise in
evaluating divided differences.

Table 2 shows the error in $\ell^{\infty}$ norm of the ten conserved
quantities for the methods considered. While all three methods perform quite
similarly, the Greenspan–Labudde or DMM scheme is the only one which preserves
the Hamiltonian up to machine precision. We also note that due to the time-
dependent nature of the initial center of mass
$\boldsymbol{C}(t,\boldsymbol{q},\boldsymbol{p})$, the error for this
conserved quantity is not up to machine precision for the
Greenspan–Labudde/DMM scheme. Specifically, since $\boldsymbol{C}$ is growing
linearly in time, the discrete counterparts, which are zero up to machine
precision, are being multiplied by a linear factor in time555For more details,
see the calculations on the time-dependent terms of DMM in (24) of A.2.,
resulting in $\sim 10^{-11}$ error for $\boldsymbol{C}$.

Method | Midpoint | Störmer–Verlet | Greenspan–Labudde/DMM
---|---|---|---
$\text{Error}[H(\boldsymbol{q},\boldsymbol{p})]$ | $7.03\cdot 10^{-12}$ | $1.72\cdot 10^{-12}$ | $2.03\cdot 10^{-17}$
$\text{Error}[P_{x}(\boldsymbol{p})]$ | $1.04\cdot 10^{-18}$ | $1.44\cdot 10^{-18}$ | $1.26\cdot 10^{-18}$
$\text{Error}[P_{y}(\boldsymbol{p})]$ | $8.84\cdot 10^{-19}$ | $1.36\cdot 10^{-18}$ | $9.44\cdot 10^{-19}$
$\text{Error}[P_{z}(\boldsymbol{p})]$ | $2.45\cdot 10^{-19}$ | $9.41\cdot 10^{-19}$ | $2.12\cdot 10^{-19}$
$\text{Error}[L_{x}(\boldsymbol{q},\boldsymbol{p})]$ | $2.14\cdot 10^{-17}$ | $4.15\cdot 10^{-18}$ | $7.44\cdot 10^{-18}$
$\text{Error}[L_{y}(\boldsymbol{q},\boldsymbol{p})]$ | $3.03\cdot 10^{-16}$ | $4.94\cdot 10^{-18}$ | $1.23\cdot 10^{-16}$
$\text{Error}[L_{z}(\boldsymbol{q},\boldsymbol{p})]$ | $5.41\cdot 10^{-16}$ | $1.36\cdot 10^{-17}$ | $2.24\cdot 10^{-16}$
$\text{Error}[C_{x}(t,\boldsymbol{q},\boldsymbol{p})]$ | $1.07\cdot 10^{-11}$ | $1.24\cdot 10^{-11}$ | $1.50\cdot 10^{-11}$
$\text{Error}[C_{y}(t,\boldsymbol{q},\boldsymbol{p})]$ | $8.89\cdot 10^{-12}$ | $1.75\cdot 10^{-11}$ | $1.11\cdot 10^{-11}$
$\text{Error}[C_{z}(t,\boldsymbol{q},\boldsymbol{p})]$ | $3.51\cdot 10^{-12}$ | $1.06\cdot 10^{-11}$ | $4.56\cdot 10^{-12}$
Table 2:  Error in conserved quantities for the ten–body solar system model.

Figure 3 shows the error in time for the ten conserved quantities of the three
methods considered. Furthermore, Figure 4 shows that the trajectories of all
three methods are in good agreement qualitatively.

Figure 3:  Error in conserved quantities versus time for the ten–body solar
system using the Midpoint method, Störmer–Verlet method and
Greenspan–Labudde/DMM scheme, plotted every 500 days.

Figure 4:  Trajectories for all three methods, plotted every 2500 days.

###  4.3 Many–body problem: Lennard–Jones Potential

Next we consider an example from molecular dynamics using the Lennard–Jones
Potential. For simplicity, we have chosen to look at a frozen Argon crystal
model in 2D. This is to avoid complications which arises when restricting
molecules to a finite domain, such as to truncate the potential to a finite
radial length and to handle periodic boundary conditions in a conservative
manner. We mentioned that [36] recently resolved these difficulties in the
case of general pairwise and three–body interaction potentials.

As with the gravitational potential example, we compare the derived DMM
scheme, with the Störmer–Verlet and the Midpoint method. For this problem, we
used the initial conditions and parameters provided by [3, Chapter I.4]. The
units are chosen to be nanoseconds $\left[ns\right]$ for time, nanometers
$\left[nm\right]$ for length, and $\left[kg\right]$ for mass. The number of
atoms is set to $n=7$ and the mass of each atom is uniformly set to
$m_{i}=66.34\cdot 10^{-27}\>\left[kg\right].$ We further set
$\sigma=0.341\>\left[nm\right],$ and $\epsilon=119.8\,k_{B}~{}\left[J\right],$
where $k_{B}=1.380658\cdot 10^{-23}\>\left[J\cdot K^{-1}\right]$ is
Boltzmann’s constant. In the actual simulation, we have rescaled the equations
by $k_{B}$ in order to mitigate round-off errors. As a result, the conserved
quantities presented are in these rescaled units. We simulate the system to a
final time $T=0.2~{}[ns]$ and $N=4\cdot 10^{4}$ time steps, corresponding to a
time step size of $\tau=50$ femtoseconds. Similar to the gravitational
potential example, we used a perturbation of the discrete solution from the
previous time step for the initial guess in the fixed point iterations.

Since the problem is two-dimensional, only six conserved quantities are
relevant and their respective $\ell^{\infty}$ norm errors are shown in Table
3. Similarly to previous examples, we see that the derived DMM scheme is the
only one preserving all conserved quantities. We note that there are larger
accumulation of round-off errors for the Hamiltonian due to cancellation
errors from the more complex divided difference expression (16) arising from
the Lennard–Jones potential. Furthermore, we observed good qualitative
agreement of trajectories with other methods in Figure 5 and favorable
conservative properties of DMM in Figure 6.

| Midpoint | Störmer-Verlet | DMM
---|---|---|---
$\text{Error}[H(\boldsymbol{q},\boldsymbol{p})]$ | $3.17\cdot 10^{-2}$ | $4.60\cdot 10^{-2}$ | $7.84\cdot 10^{-11}$
$\text{Error}[P_{x}(\boldsymbol{q},\boldsymbol{p})]$ | $9.66\cdot 10^{-15}$ | $1.95\cdot 10^{-14}$ | $8.16\cdot 10^{-15}$
$\text{Error}[P_{y}(\boldsymbol{q},\boldsymbol{p})]$ | $5.50\cdot 10^{-15}$ | $3.08\cdot 10^{-14}$ | $5.47\cdot 10^{-15}$
$\text{Error}[L_{z}(\boldsymbol{q},\boldsymbol{p})]$ | $4.41\cdot 10^{-15}$ | $4.86\cdot 10^{-15}$ | $3.36\cdot 10^{-15}$
$\text{Error}[C_{x}(t,\boldsymbol{q},\boldsymbol{p})]$ | $3.53\cdot 10^{-14}$ | $9.37\cdot 10^{-14}$ | $2.10\cdot 10^{-14}$
$\text{Error}[C_{y}(t,\boldsymbol{q},\boldsymbol{p})]$ | $3.41\cdot 10^{-14}$ | $1.11\cdot 10^{-13}$ | $2.67\cdot 10^{-14}$
Table 3:  Error in conserved quantities for the frozen Argon crystal model.
Figure 5: Trajectories of the seven Argon atoms for all three methods. Figure
6: Error in conserved quantities versus time for the frozen Argon crystal
model.

###  4.4 Point vortex problem: Planar case

Next, we investigate and compare numerical results of DMM with the Midpoint
method and the standard explicit fourth-order Runge–Kutta method.

The test consists of evolving $n=1000$ randomly distributed vortices. The
initial locations were sampled from a uniform distribution on
$\left[-5,5\right]^{2}$ and post-processed to ensure that no two vortices were
closer than a minimum distance of $10/n=10^{-2}$. Furthermore, the vorticity
strengths $\Gamma_{i}$ were sampled uniformly from $[-1,1]/n$. We run the
simulation to a final time of $T=100$ for a total of $N=1000$ time steps,
corresponding to a time step size $\tau=10^{-1}$. As with the many–body
examples, we have used a perturbation of the discrete solution from the
previous time step for the initial guess of the fixed point iterations.

Figure 7:  Error in conserved quantities versus time for the planar point
vortex problem.

Figure 7 shows the error in conserved quantities over time and Figure 8 shows
the trajectories of all $n=1000$ vortices for all three methods considered. We
observe that most vortices have qualitatively similar trajectories for all
three methods. However, there are a few vortices, with more subtle
interactions, showing rather different trajectories. Specifically, we
highlight these differences in Figure 8 and by zooming in on their dynamics in
Figure 9.

Figure 8: Trajectories for all three methods for the planar point vortex problem with $n=1000$ vortices. |
---|---
Figure 9: Closeup of trajectory differences for all three methods for the
planar point vortex problem with $n=1000$ vortices.

While the final time here is relatively short, we emphasize that for long
integrations of these point vortex equations, these small differences in the
trajectories will likely amplify, with the DMM being the only method
preserving the energy up to machine precision.

Table 4 shows the error in $\ell^{\infty}$ norm for all four conserved
quantities.

| Midpoint | RK4 | DMM
---|---|---|---
$\text{Error}[P_{x}(\boldsymbol{x})]$ | $4.15\cdot 10^{-15}$ | $2.64\cdot 10^{-16}$ | $2.78\cdot 10^{-16}$
$\text{Error}[P_{y}(\boldsymbol{y})]$ | $3.16\cdot 10^{-15}$ | $1.22\cdot 10^{-16}$ | $1.20\cdot 10^{-16}$
$\text{Error}[L(\boldsymbol{x},\boldsymbol{y})]$ | $4.69\cdot 10^{-14}$ | $8.02\cdot 10^{-9}$ | $1.11\cdot 10^{-15}$
$\text{Error}[H(\boldsymbol{x},\boldsymbol{y})]$ | $3.49\cdot 10^{-9}$ | $7.96\cdot 10^{-8}$ | $2.14\cdot 10^{-17}$
Table 4: Error in conserved quantities for $n=1000$ point vortices on the
plane.

This table shows that all methods preserve the linear momentum up to machine
precision, with the Midpoint method additionally preserving the angular
momentum666This is expected as Midpoint method preserves all quadratic
invariants, see [3, Chapter IV.2]. The DMM is the only method that also
preserves the Hamiltonian up to machine precision. This is particularly
noteworthy as the DMM method is only a second–order method, in comparison to
the fourth-order Runge–Kutta method.

###  4.5 Point vortex problem: Spherical case

In this final example, similar to the previous section, we investigate and
compare numerical results of DMM with the Midpoint method and the standard
explicit fourth-order Runge–Kutta method applied to the point vortex equations
on the unit sphere.

As with the planar case, the test consists of evolving $n=1000$ randomly
distributed vortices. The initial locations were sampled from a uniform
distribution on the unit sphere according to the procedure proposed in [37].
As in the planar case, we filtered the initial sampled locations to ensure
that no two vortices are closer than a minimum distance of $4\pi/n\approx
1.3\cdot 10^{-2}$. Furthermore, the vorticity strengths $\Gamma_{i}$ were
sampled uniformly from $[-1,1]/n$. We run the simulation to a final time of
$T=100$, with $N=1000$ time steps, corresponding to a time step size of
$\tau=10^{-1}$. As before, we used a perturbation of the discrete solution
from the previous time step for the initial guess in the fixed point
iterations.

We show in Table 5 the error in $\ell^{\infty}$ norm for all four conserved
quantities.

| Midpoint | RK4 | DMM
---|---|---|---
$\text{Error}[P_{x}(\boldsymbol{x})]$ | $5.03\cdot 10^{-17}$ | $3.99\cdot 10^{-17}$ | $4.16\cdot 10^{-17}$
$\text{Error}[P_{y}(\boldsymbol{x})]$ | $3.95\cdot 10^{-17}$ | $4.12\cdot 10^{-17}$ | $4.47\cdot 10^{-17}$
$\text{Error}[P_{z}(\boldsymbol{x})]$ | $3.47\cdot 10^{-17}$ | $4.86\cdot 10^{-17}$ | $5.38\cdot 10^{-17}$
$\text{Error}[H(\boldsymbol{x})]$ | $5.44\cdot 10^{-10}$ | $1.55\cdot 10^{-11}$ | $3.73\cdot 10^{-18}$
Table 5: Error in conserved quantities versus time for $n=1000$ point vortices
on the sphere.

Figure 10 depicts the error in conserved quantities over time for all three
methods and Figure 11 shows the trajectories of all $1000$ vortices. As with
the planar case, we observe that most trajectories of the three methods are
qualitatively indistinguishable, although certain vortices do exhibit
increasingly diverging trajectories over the short integration interval. This
is highlighted by zooming in on specific areas presented in Figure 12.

Figure 10: Error in conserved quantities for the point vortex problem on the
unit sphere with $n=1000$ vortices.

Figure 11: Trajectories for all three methods for the point vortex problem on
the unit sphere with $n=1000$ vortices.

|
---|---
Figure 12: Closeup on some diverging trajectories for all three methods for
the point vortex problem on the unit sphere with $n=1000$ vortices.

##  5 Conclusions

In this paper, we have constructed several conservative numerical schemes for
important mathematical models arising from a wide variety of different fields.
The framework used to derive such numerical schemes is the Discrete Multiplier
Method, or DMM, originally developed in [18]. DMM is not only suitable for
low–dimensional dynamical systems with multiple conserved quantities, as shown
in [18], but also as we showed here, it is applicable for constructing
conservative schemes for many–body systems. As the DMM does not require
additional geometric structures from the underlying dynamical system other
than the presence of invariants themselves, this approach can potentially be
applied to a wide variety of other dynamical systems arising in the
mathematical sciences.

With the derived conservative schemes for many–body systems, there are still
practical drawbacks which we wish to improve in future work. The main drawback
of the derived DMM schemes so far is that they are implicit, as it is not
currently known if there are explicit DMM schemes for general dynamical
systems. For large many-body Hamiltonian systems, explicit methods such as
Störmer-Verlet method and higher–order symplectic splitting methods are often
preferred due to their lower computational costs when comparing accuracy
versus number of force evaluations. To reduce the computational costs of the
implicit DMM schemes, one can use high–order DMM schemes and specific quasi-
Newton methods aimed at improving the efficiency of solving nonlinear
equations arising from DMM.

Despite these current practical limitations, we emphasize that our numerical
results, specifically in the point vortex examples, indicate that a
higher–order method, while more accurate at approximating the solution, does
not necessarily imply it is more accurate at preserving conserved quantities
when compared to a lower–order conservative method.

## Acknowledgements

ATSW was partially supported by the CRM and the NSERC Discovery Grant program.
AB was supported by the Canada Research Chairs program, the InnovateNL
Leverage R&D program and the NSERC Discovery Grant program. JCN was supported
by the NSERC Discovery Grant program. Additionally, JCN would like to thank
Prof. Wenjun Ying of the Natural Science Institute of the Shanghai Jiao Tong
University for his kind hosting. The environment provided by the NSI was
invaluable during the early stages of this project.

## References

  * [1] F. Calogero, Classical many-body problems amenable to exact treatments, Vol. 66 of Lecture notes in physics, Springer, Berlin, 2003.
  * [2] S. Blanes, F. Casas, A Concise Introduction to Geometric Numerical Integration, Vol. 23, CRC Press, Boca Raton, 2016.
  * [3] E. Hairer, C. Lubich, G. Wanner, Geometric numerical integration: structure-preserving algorithms for ordinary differential equations, Springer, Berlin, 2006.
  * [4] B. Leimkuhler, S. Reich, Simulating Hamiltonian dynamics, Cambridge University Press, Cambridge, 2004.
  * [5] R. de Vogelaere, methods of integration which preserve the contact transformation property of the hamiltonian equations, Tech. rep.
  * [6] R. D. Ruth, A canonical integration technique, IEEE Trans. Nucl. Sci 30 (4) (1983) 2669–2671.
  * [7] K. Feng, Difference schemes for Hamiltonian formalism and symplectic geometry, J. Comput. Math. 4 (3) (1986) 279–289.
  * [8] J.-M. Sanz-Serna, M.-P. Calvo, Numerical Hamiltonian problems, Vol. 7 of Applied Mathematics and Mathematical Computation, Chapman & Hall, London, 1994\.
  * [9] J. E. Marsden, M. West, Discrete mechanics and variational integrators, Acta Numer. 10 (2001) 357–514.
  * [10] G. Zhong, J. E. Marsden, Lie–Poisson Hamilton–Jacobi theory and Lie–Poisson integrators, Phys. Lett. A 133 (3) (1988) 134–139.
  * [11] M.-P. Calvo, E. Hairer, Accurate long-term integration of dynamical systems, Appl. Numer. Math. 18 (1-3) (1995) 95–105.
  * [12] A. T. S. Wan, J.-C. Nave, On the arbitrarily long-term stability of conservative methods, SIAM J. Numer. Anal. 56 (5) (2018) 2751–2775.
  * [13] G. R. W. Quispel, D. I. McLaren, A new class of energy-preserving numerical integration methods, J. Phys. A: Math. Theor. 41 (4) (2008) 045206.
  * [14] E. Hairer, Energy-preserving variant of collocation methods, J. Numer. Anal. Ind. Appl. Math. 5 (2010) 73–84.
  * [15] D. Cohen, E. Hairer, Linear energy-preserving integrators for poisson systems, BIT 51 (1) (2011) 91–101.
  * [16] G. R. W. Quispel, G. S. Turner, Discrete gradient methods for solving odes numerically while preserving a first integral, J. Phys. A: Math. Gen. 29 (13) (1996) L341.
  * [17] Y. Nambu, Generalized Hamiltonian dynamics, Phys. Rev. D 7 (8) (1973) 2405–2412.
  * [18] A. T. S. Wan, A. Bihlo, J.-C. Nave, Conservative methods for dynamical systems, SIAM J. Numer. Anal. 55 (5) (2017) 2255–2285.
  * [19] P. J. Olver, Application of Lie groups to differential equations, Springer, New York, 2000.
  * [20] G. W. Bluman, A. F. Cheviakov, S. C. Anco, Application of symmetry methods to partial differential equations, Springer, New York, 2010.
  * [21] A. T. S. Wan, A. Bihlo, J.-C. Nave, The multiplier method to construct conservative finite difference schemes for ordinary and partial differential equations, SIAM J. Numer. Anal. 54 (1) (2016) 86–119.
  * [22] R. A. Labudde, D. Greenspan, Energy and momentum conserving methods of arbitrary order for the numerical integration of equations of motion, Numer. Math. 25 (4) (1975).
  * [23] D. Greenspan, $N$-body Problems and Models, World Scientific, Singapore, 2004.
  * [24] J. Hofbauer, K. Sigmund, Evolutionary games and population dynamics, Cambridge University Press, 1998.
  * [25] R. Schimming, Conservation laws for Lotka–Volterra models, Math. Methods Appl. Sci. 26 (17) (2003) 1517–1528.
  * [26] H. Helmholtz, Über Integrale der hydrodynamischen Gleichungen, welche den Wirbelbewegungen entsprechen, J. Reine Angew. Math. 55 (1858) 25–55.
  * [27] G. R. Kirchhoff, Vorlesungen über mathematische Physik. Mechanik, Vol. 1, B.G. Teubner, Leipzig, 1883.
  * [28] P. K. Newton, The N-Vortex Problem: Analytical Techniques, Springer, New York, 2001\.
  * [29] H. Aref, Point vortex dynamics: a classical mathematics playground, J. Math. Phys. 48 (6) (2007) 065401.
  * [30] P. J. Channell, C. Scovel, Symplectic integration of hamiltonian systems, Nonlinearity 3 (2) (1990) 231.
  * [31] D. I. Pullin, P. G. Saffman, Long-time symplectic integration: the example of four-vortex motion, Proc. Math. Phys. Eng. Sci. 432 (1886) (1991) 481–494.
  * [32] C. Scovel, Symplectic numerical integration of Hamiltonian systems, in: The geometry of Hamiltonian systems, Springer, New York, 1991, pp. 463–496.
  * [33] K. W. Myerscough, J. Frank, Explicit, parallel Poisson integration of point vortices on the sphere, J. Comput. Appl. Math. 304 (2016) 100–119.
  * [34] J. Vankerschaver, M. Leok, A novel formulation of point vortex dynamics on the sphere: geometrical and numerical aspects, J. Nonlinear Sci. 24 (1) (2014) 1–37.
  * [35] W. M. Folkner, J. G. Williams, D. H. Boggs, R. S. Park, P. Kuchynka, The planetary and lunar ephemerides de430 and de431, The Interplanetary Network Progress Report 42 (196) (2014).
  * [36] M. Schiebl, I. Romero, Energy-momentum conserving integration schemes for molecular dynamics, Comput. Mech. 67 (2021) 915–935.
  * [37] G. Marsaglia, Choosing a point from the surface of a sphere, Ann. Math. Stat. 43 (2) (1972) 645–646.

##  Appendix A Verification of conservative and symmetric properties

In this Appendix, we present some details of the derivations of the
conservative schemes and verification of symmetric property used throughout
this paper.

Due to the common appearance of the function $g(z)=\dfrac{\log z}{z-1}$, we
first show the following elementary lemma useful in proving the symmetric
property of our schemes.

######  Lemma 1.

Let $g(z)=\dfrac{\log z}{z-1}$ and $h$ be any scalar function which maps from
the phase space into real numbers excluding the zero, i.e.
$h:\overbrace{\mathbb{R}^{d}\times\cdots\times\mathbb{R}^{d}}^{n\text{
copies}}\rightarrow\mathbb{R}\backslash\\{0\\}$. Then,

$\displaystyle\frac{1}{h(\boldsymbol{x}^{k})}g\left(\frac{h(\boldsymbol{x}^{k+1})}{h(\boldsymbol{x}^{k})}\right)=\frac{1}{h(\boldsymbol{x}^{k+1})}g\left(\frac{h(\boldsymbol{x}^{k})}{h(\boldsymbol{x}^{k+1})}\right).$
(23)

###### Proof.

$\displaystyle\frac{1}{h(\boldsymbol{x}^{k})}g\left(\frac{\boldsymbol{x}^{k+1}}{\boldsymbol{x}^{k}}\right)$
$\displaystyle=\frac{1}{h(\boldsymbol{x}^{k})}\frac{\log\left(\dfrac{h(\boldsymbol{x}^{k+1})}{h(\boldsymbol{x}^{k})}\right)}{\dfrac{h(\boldsymbol{x}^{k+1})}{h(\boldsymbol{x}^{k})}-1}=\frac{\log\left(h(\boldsymbol{x}^{k+1})\right)-\log\left(h(\boldsymbol{x}^{k})\right)}{h(\boldsymbol{x}^{k+1})-h(\boldsymbol{x}^{k})}$
$\displaystyle=\frac{1}{h(\boldsymbol{x}^{k+1})}\frac{\log\left(\dfrac{h(\boldsymbol{x}^{k})}{h(\boldsymbol{x}^{k+1})}\right)}{\dfrac{h(\boldsymbol{x}^{k})}{h(\boldsymbol{x}^{k+1})}-1}=\frac{1}{h(\boldsymbol{x}^{k+1})}g\left(\frac{h(\boldsymbol{x}^{k})}{h(\boldsymbol{x}^{k+1})}\right).$

∎

In other words,
$\dfrac{1}{h(\boldsymbol{x}^{k})}g\left(\dfrac{h(\boldsymbol{x}^{k+1})}{h(\boldsymbol{x}^{k})}\right)$
is symmetric under the permutation $k\leftrightarrow k+1$.

###  A.1 Lotka–Volterra systems

We first show that $V(\boldsymbol{x})$ is a conserved quantity of (10) and the
derived scheme is conservative and symmetric using Lemma 1.

Using (4a), the associated multiplier $1\times n$ matrix is given by

$\displaystyle\Lambda(\boldsymbol{x}):=\begin{bmatrix}d_{i}\left(\dfrac{\xi_{i}}{x_{i}}-1\right)\end{bmatrix}_{1\leq
i\leq n}$

which verifies $V(\boldsymbol{x})$ is indeed a conserved quantity since,

$\displaystyle\Lambda(\boldsymbol{x})\boldsymbol{F}(\boldsymbol{x})$
$\displaystyle=\sum_{i=1}^{n}d_{i}\left(\dfrac{\xi_{i}}{x_{i}}-1\right)\dot{x}_{i}-\sum_{i,j=1}^{n}d_{i}\left(\xi_{i}-x_{i}\right)a_{ij}(x_{j}-\xi_{j})$
$\displaystyle=D_{t}V(\boldsymbol{x})+\underbrace{(\boldsymbol{x}-\boldsymbol{\xi})^{T}DA(\boldsymbol{x}-\boldsymbol{\xi})}_{=0\text{
since }DA=-A^{T}D}=D_{t}V(\boldsymbol{x}).$

We employ DMM to derive conservative schemes for (10). Since $V$ is a linear
combination of single variable functions, for any permutation $\sigma\in
S_{n}$ of $x_{1},\dots,x_{n}$, the discrete multiplier matrix is given by,

$\displaystyle\Lambda^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{x}^{k})$
$\displaystyle=\begin{bmatrix}d_{i}\dfrac{\Delta}{\Delta
x_{i}}\left(\xi_{i}\log x_{i}-x_{i}\right)\end{bmatrix}_{1\leq i\leq
n}^{T}=\begin{bmatrix}d_{i}\left(\dfrac{\xi_{i}}{x_{i}^{k}}g\left(\dfrac{x_{i}^{k+1}}{x_{i}^{k}}\right)\right)\end{bmatrix}_{1\leq
i\leq n}^{T},$

with $g(z)=\dfrac{\log z}{z-1}$. Next, to discretize the right hand side
$\boldsymbol{f}$, we first rewrite

$\displaystyle\boldsymbol{f}(\boldsymbol{x})$
$\displaystyle=\begin{bmatrix}x_{i}\sum_{j=1}^{n}a_{ij}x_{j}\left(1-\dfrac{\xi_{j}}{x_{j}}\right)\end{bmatrix}_{1\leq
i\leq n}.$

Since $g(z)=1+\mathcal{O}(z-1)$, let us propose by consistency of
$\boldsymbol{f}^{\tau}$ to $\boldsymbol{f}$ that

$\displaystyle\boldsymbol{f}^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{x}^{k})$
$\displaystyle=\begin{bmatrix}x_{i}^{\tau}\displaystyle\sum_{j=1}^{n}a_{ij}x_{j}^{\tau}\left(1-\dfrac{\xi_{j}}{x_{j}^{k}}g\left(\dfrac{x_{j}^{k+1}}{x_{j}^{k}}\right)\right)\end{bmatrix}_{1\leq
i\leq n},$

where $x_{i}^{\tau}$ is any consistent discretization of $x_{i}$. Indeed, (9b)
is satisfied for this choice of $\boldsymbol{f}^{\tau}$,

$\displaystyle\Lambda^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{x}^{k})\boldsymbol{f}^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{x}^{k})$
$\displaystyle=\sum_{i,j=1}^{n}d_{i}\left(\dfrac{\xi_{i}}{x_{i}^{k}}g\left(\dfrac{x_{i}^{k+1}}{x_{i}^{k}}\right)-1\right)x_{i}^{\tau}a_{ij}x_{j}^{\tau}\left(1-\dfrac{\xi_{j}}{x_{j}^{k}}g\left(\dfrac{x_{j}^{k+1}}{x_{j}^{k}}\right)\right)$
$\displaystyle=-\boldsymbol{y}^{T}DA\boldsymbol{y}=0,$

where the components of $\boldsymbol{y}$ are given by
$y_{i}=x_{i}^{\tau}\left(1-\dfrac{\xi_{i}}{x_{i}^{k}}g\left(\dfrac{x_{i}^{k+1}}{x_{i}^{k}}\right)\right)$.

Thus, a conservative scheme for (10) is given by (12). Now by Lemma 1, since
$\dfrac{1}{x_{j}^{k}}g\left(\dfrac{x_{j}^{k+1}}{x_{j}^{k}}\right)$ is
symmetric with $h(\boldsymbol{x}^{k})=x_{j}^{k}$ for all $j=1,\dots,n$, the
discretization (12) is symmetric provided $x_{i}^{\tau},x_{i}^{\tau}$ are
symmetric. Specifically, the symmetric choices
$x_{i}^{\tau}:=\frac{1}{2}(x_{i}^{k}+x_{i}^{k+1})$ or
$\sqrt{x_{i}^{k}x_{i}^{k+1}}$ were selected in (12).

###  A.2 Many–body problem with pairwise radial potentials

Next, we verify the conserved quantities of (13) and show the derived scheme
is conservative and symmetric.

One can verify (13) has the conserved quantities (14) using (9a). Indeed,
define the vector of conserved quantities $\boldsymbol{\psi}$ as

$\displaystyle\boldsymbol{\psi}(t,\boldsymbol{q},\boldsymbol{p}):=\begin{pmatrix}H(\boldsymbol{q},\boldsymbol{p})\\\
\boldsymbol{P}(\boldsymbol{q},\boldsymbol{p})\\\
\boldsymbol{L}(\boldsymbol{q},\boldsymbol{p})\\\
\boldsymbol{C}(t,\boldsymbol{q},\boldsymbol{p})\end{pmatrix}.$

Then by (4a), the associated $10\times(6n)$ multiplier matrix is given by,

$\displaystyle\Lambda(t,\boldsymbol{q},\boldsymbol{p}):=\begin{pmatrix}\begin{bmatrix}\sum\limits_{j=1,j\neq
i}^{n}\dfrac{\partial V_{ij}}{\partial
q_{ij}}\dfrac{1}{q_{ij}}(\boldsymbol{q}_{i}-\boldsymbol{q}_{j})\end{bmatrix}_{1\leq
i\leq
n}^{T}&\begin{bmatrix}\dfrac{1}{m_{i}}\boldsymbol{p}_{i}\end{bmatrix}_{1\leq
i\leq n}^{T}\\\ \begin{bmatrix}0_{3\times 3}\end{bmatrix}_{1\leq i\leq
n}^{T}&\begin{bmatrix}I_{3\times 3}\end{bmatrix}_{1\leq i\leq n}^{T}\\\
-[\Omega(\boldsymbol{p}_{i})]_{1\leq i\leq
n}^{T}&[\Omega(\boldsymbol{q}_{i})]_{1\leq i\leq n}^{T}\\\
\begin{bmatrix}\dfrac{m_{i}}{M}I_{3\times 3}\end{bmatrix}_{1\leq i\leq
n}^{T}&-\begin{bmatrix}\dfrac{t}{M}I_{3\times 3}\end{bmatrix}_{1\leq i\leq
n}^{T}\end{pmatrix},$

where $\Omega(\boldsymbol{x})$ denotes the skew-symmetric $3\times 3$ matrix,

$\displaystyle\Omega(\boldsymbol{x}):=\begin{bmatrix}0&-z&y\\\ z&0&-x\\\
-y&x&0\end{bmatrix},$

associated with the cross product of
$\boldsymbol{x}=(x,y,z)^{T}\in\mathbb{R}^{3}$ such that
$\Omega(\boldsymbol{x})\boldsymbol{y}=\boldsymbol{x}\times\boldsymbol{y}$ for
$\boldsymbol{y}\in\mathbb{R}^{3}$. Thus, (4b) is satisfied since,

$\displaystyle\Lambda\boldsymbol{f}$
$\displaystyle=\begin{pmatrix}\sum\limits_{1\leq i\neq j\leq n}\dfrac{\partial
V_{ij}}{\partial
q_{ij}}\dfrac{1}{q_{ij}}\left[\dfrac{1}{m_{i}}\boldsymbol{p}_{i}^{T}(\boldsymbol{q}_{i}-\boldsymbol{q}_{j})-(\boldsymbol{q}_{i}-\boldsymbol{q}_{j})^{T}\dfrac{\boldsymbol{p}_{i}}{m_{i}}\right]\\\
\sum\limits_{1\leq i\neq j\leq n}\dfrac{\partial V_{ij}}{\partial
q_{ij}}\dfrac{1}{q_{ij}}(\boldsymbol{q}_{i}-\boldsymbol{q}_{j})\\\
-\sum\limits_{1\leq i\leq
n}\dfrac{1}{m_{i}}\Omega(\boldsymbol{p}_{i})\boldsymbol{p}_{i}+\sum\limits_{1\leq
i\neq j\leq n}\dfrac{\partial V_{ij}}{\partial
q_{ij}}\dfrac{1}{q_{ij}}\Omega(\boldsymbol{q}_{i})(\boldsymbol{q}_{i}-\boldsymbol{q}_{j})\\\
\dfrac{1}{M}\sum\limits_{1\leq i\leq n}\boldsymbol{p}_{i}-t\sum\limits_{1\leq
i\leq n}\dfrac{\partial V_{ij}}{\partial
q_{ij}}\dfrac{1}{q_{ij}}(\boldsymbol{q}_{i}-\boldsymbol{q}_{j})\end{pmatrix}$
$\displaystyle=\begin{pmatrix}0\\\ \sum\limits_{1\leq i<j\leq
n}\dfrac{\partial V_{ij}}{\partial
q_{ij}}\left(\dfrac{\boldsymbol{q}_{i}-\boldsymbol{q}_{j}}{q_{ij}}+\dfrac{\boldsymbol{q}_{j}-\boldsymbol{q}_{i}}{q_{ji}}\right)\\\
-\sum\limits_{1\leq i\leq
n}\dfrac{1}{m_{i}}\boldsymbol{p}_{i}\times\boldsymbol{p}_{i}+\sum\limits_{1\leq
i<j\leq n}\dfrac{\partial V_{ij}}{\partial
q_{ij}}\left(\dfrac{\boldsymbol{q}_{i}}{q_{ij}}\times(\boldsymbol{q}_{i}-\boldsymbol{q}_{j})+\dfrac{\boldsymbol{q}_{j}}{q_{ji}}\times(\boldsymbol{q}_{j}-\boldsymbol{q}_{i})\right)\\\
\dfrac{\boldsymbol{P}}{M}-t\sum\limits_{1\leq i\leq n}\dfrac{\partial
V_{ij}}{\partial
q_{ij}}\dfrac{1}{q_{ij}}\left(\dfrac{\boldsymbol{q}_{i}-\boldsymbol{q}_{j}}{q_{ij}}+\dfrac{\boldsymbol{q}_{j}-\boldsymbol{q}_{i}}{q_{ji}}\right)\end{pmatrix}=\begin{pmatrix}0\\\
\boldsymbol{0}\\\ \boldsymbol{0}\\\ \dfrac{\boldsymbol{P}}{M}\end{pmatrix}$
$\displaystyle=-\partial_{t}\boldsymbol{\psi},$

where the second last equality follows from
$q_{ji}=q_{ij},\boldsymbol{q}_{j}-\boldsymbol{q}_{i}=-(\boldsymbol{q}_{i}-\boldsymbol{q}_{j})$
and properties of the cross-product.

We now employ DMM to derive conservative schemes for (13). For simplicity, we
will propose consistent choices of $D_{t}^{\tau}\boldsymbol{x}$,
$D_{t}^{\tau}\boldsymbol{\psi}$, $\partial_{t}^{\tau}\boldsymbol{\psi}$,
$\Lambda^{\tau}$ and $\boldsymbol{f}^{\tau}$ such that conditions (9a) and
(9b) are satisfied. To accomplish this, we define
$D_{t}^{\tau}\boldsymbol{x}$, $D_{t}^{\tau}\boldsymbol{\psi}$ and
$\partial_{t}^{\tau}\boldsymbol{\psi}$ by,

$\displaystyle
D_{t}^{\tau}\boldsymbol{x}:=\dfrac{1}{\tau}\begin{pmatrix}\Delta\boldsymbol{q}\\\
\Delta\boldsymbol{p}\end{pmatrix},\hskip
14.22636ptD_{t}^{\tau}\boldsymbol{\psi}:=\dfrac{1}{\tau}\begin{pmatrix}\Delta
H\\\ \Delta\boldsymbol{P}\\\ \Delta\boldsymbol{L}\\\
\Delta\boldsymbol{C}\end{pmatrix},\hskip
14.22636pt\partial_{t}^{\tau}\boldsymbol{\psi}:=\begin{pmatrix}0\\\
\boldsymbol{0}\\\ \boldsymbol{0}\\\
-\dfrac{\overline{\boldsymbol{P}}}{M}\end{pmatrix}.$

For the discrete multiplier $\Lambda^{\tau}$ and the discrete right hand side
$\boldsymbol{f}^{\tau}$, we define

$\displaystyle\Lambda^{\tau}(t^{k+1},\boldsymbol{q}^{k+1},\boldsymbol{p}^{k+1},t^{k},\boldsymbol{q}^{k},\boldsymbol{p}^{k})$
$\displaystyle:=\begin{pmatrix}\begin{bmatrix}\sum\limits_{j=1,j\neq
i}^{n}\dfrac{\Delta V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})\end{bmatrix}_{1\leq
i\leq
n}^{T}&\begin{bmatrix}\dfrac{1}{m_{i}}\overline{\boldsymbol{p}_{i}}\end{bmatrix}_{1\leq
i\leq n}^{T}\\\ \begin{bmatrix}0_{3\times 3}\end{bmatrix}_{1\leq i\leq
n}^{T}&\begin{bmatrix}I_{3\times 3}\end{bmatrix}_{1\leq i\leq n}^{T}\\\
-[\Omega(\overline{\boldsymbol{p}_{i}})]_{1\leq i\leq
n}^{T}&[\Omega(\overline{\boldsymbol{q}_{i}})]_{1\leq i\leq n}^{T}\\\
\begin{bmatrix}\dfrac{m_{i}}{M}I_{3\times 3}\end{bmatrix}_{1\leq i\leq
n}^{T}&-\begin{bmatrix}\dfrac{t^{k+1}+t^{k}}{2M}I_{3\times
3}\end{bmatrix}_{1\leq i\leq n}^{T}\end{pmatrix},$
$\displaystyle\boldsymbol{f}^{\tau}(t^{k+1},\boldsymbol{q}^{k+1},\boldsymbol{p}^{k+1},t^{k},\boldsymbol{q}^{k},\boldsymbol{p}^{k})$
$\displaystyle:=\begin{pmatrix}\begin{bmatrix}\dfrac{\overline{\boldsymbol{p}_{i}}}{m_{i}}\end{bmatrix}_{1\leq
i\leq n}\\\ -\begin{bmatrix}\sum\limits_{j=1,j\neq i}^{n}\dfrac{\Delta
V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})\end{bmatrix}_{1\leq
i\leq n}\end{pmatrix}.$

It can be seen that $\Lambda^{\tau}$ and $\boldsymbol{f}^{\tau}$ are
consistent to $\Lambda$ and $\boldsymbol{f}$. We now verify the above choices
satisfy condition (9b).

$\displaystyle\Lambda^{\tau}\boldsymbol{f}^{\tau}$
$\displaystyle=\begin{pmatrix}\sum\limits_{1\leq i\neq j\leq n}\dfrac{\Delta
V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}\left[\dfrac{1}{m_{i}}\overline{\boldsymbol{p}_{i}}^{T}(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})-(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})^{T}\dfrac{\overline{\boldsymbol{p}_{i}}}{m_{i}}\right]\\\
\sum\limits_{1\leq i\neq j\leq n}\dfrac{\Delta V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})\\\
-\sum\limits_{1\leq i\leq
n}\dfrac{1}{m_{i}}\Omega(\overline{\boldsymbol{p}_{i}})\overline{\boldsymbol{p}_{i}}+\sum\limits_{1\leq
i\neq j\leq n}\dfrac{\Delta V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}\Omega(\overline{\boldsymbol{q}_{i}})(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})\\\
\dfrac{1}{M}\sum\limits_{1\leq i\leq
n}\overline{\boldsymbol{p}_{i}}-\dfrac{t^{k+1}+t^{k}}{2M}\sum\limits_{1\leq
i\neq j\leq n}\dfrac{\Delta V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})\end{pmatrix}$
$\displaystyle=\begin{pmatrix}0\\\ \sum\limits_{1\leq i<j\leq n}\dfrac{\Delta
V_{ij}}{\Delta
q_{ij}}\left(\dfrac{\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}}}{\overline{q_{ij}}}+\dfrac{\overline{\boldsymbol{q}_{j}}-\overline{\boldsymbol{q}_{i}}}{\overline{q_{ji}}}\right)\\\
-\sum\limits_{1\leq i\leq
n}\dfrac{1}{m_{i}}\overline{\boldsymbol{p}_{i}}\times\overline{\boldsymbol{p}_{i}}+\sum\limits_{1\leq
i<j\leq n}\dfrac{\Delta V_{ij}}{\Delta
q_{ij}}\left(\dfrac{\overline{\boldsymbol{q}_{i}}}{\overline{q_{ij}}}\times(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})+\dfrac{\overline{\boldsymbol{q}_{j}}}{\overline{q_{ji}}}\times(\overline{\boldsymbol{q}_{j}}-\overline{\boldsymbol{q}_{i}})\right)\\\
\dfrac{\overline{\boldsymbol{P}}}{M}-\dfrac{t^{k+1}+t^{k}}{2M}\sum\limits_{1\leq
i<j\leq n}\dfrac{\Delta V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}\left(\dfrac{\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}}}{\overline{q_{ij}}}+\dfrac{\overline{\boldsymbol{q}_{j}}-\overline{\boldsymbol{q}_{i}}}{\overline{q_{ji}}}\right)\end{pmatrix}$
(24) $\displaystyle=\begin{pmatrix}0\\\ \boldsymbol{0}\\\ \boldsymbol{0}\\\
\dfrac{\overline{\boldsymbol{P}}}{M}\end{pmatrix}=-\partial_{t}^{\tau}\boldsymbol{\psi}.$

By direct computation or using divided difference calculus, one can show that

$\displaystyle\Delta\left(\dfrac{\boldsymbol{p}_{i}^{T}\boldsymbol{p}_{i}}{2}\right)$
$\displaystyle=\overline{\boldsymbol{p}_{i}}^{T}\Delta\boldsymbol{p}_{i},$
$\displaystyle\Delta V_{ij}$ $\displaystyle=\dfrac{\Delta V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})^{T}\Delta(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}}),$
$\displaystyle\Delta\left(\boldsymbol{q}_{i}\times\boldsymbol{p}_{i}\right)$
$\displaystyle=\overline{\boldsymbol{q}_{i}}\times\Delta\boldsymbol{p}_{i}+\Delta\boldsymbol{q}_{i}\times\overline{\boldsymbol{p}_{i}},$
$\displaystyle\Delta(t\boldsymbol{p}_{i})$
$\displaystyle=\left(\dfrac{t^{k+1}+t^{k}}{2}\right)\Delta\boldsymbol{p}_{i}+\tau\overline{\boldsymbol{p}_{i}}.$

So by linearity of the forward difference $\Delta$, condition (9a) also holds
since,

$\displaystyle\Lambda^{\tau}D_{t}^{\tau}\boldsymbol{x}$
$\displaystyle=\dfrac{1}{\tau}\begin{pmatrix}\sum\limits_{1\leq i\leq
n}\dfrac{1}{m_{i}}\overline{\boldsymbol{p}_{i}}^{T}\Delta\boldsymbol{p}_{i}+\sum\limits_{1\leq
i\neq j\leq n}\dfrac{\Delta V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})^{T}\Delta\boldsymbol{q}_{i}\\\
\sum\limits_{1\leq i\leq n}\Delta\boldsymbol{p}_{i}\\\ \sum\limits_{1\leq
i\leq
n}\overline{\boldsymbol{q}_{i}}\times\Delta\boldsymbol{p}_{i}-\overline{\boldsymbol{p}_{i}}\times\Delta\boldsymbol{q}_{i}\\\
\dfrac{1}{M}\sum\limits_{1\leq i\leq
n}m_{i}\Delta\boldsymbol{q}_{i}-\dfrac{1}{M}\sum\limits_{1\leq i\leq
n}\left(\dfrac{t^{k+1}+t^{k}}{2}\right)\Delta\boldsymbol{p}_{i}\end{pmatrix}$
$\displaystyle=\dfrac{1}{\tau}\begin{pmatrix}\sum\limits_{1\leq i\leq
n}\Delta\left(\dfrac{\boldsymbol{p}_{i}^{T}\boldsymbol{p}_{i}}{2m_{i}}\right)+\sum\limits_{1\leq
i<j\leq n}\dfrac{\Delta V_{ij}}{\Delta
q_{ij}}\dfrac{1}{\overline{q_{ij}}}(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})^{T}\Delta(\overline{\boldsymbol{q}_{i}}-\overline{\boldsymbol{q}_{j}})\\\
\Delta\left(\sum\limits_{1\leq i\leq n}\boldsymbol{p}_{i}\right)\\\
\sum\limits_{1\leq i\leq
n}\Delta\left(\boldsymbol{q}_{i}\times\boldsymbol{p}_{i}\right)\\\
\Delta\left(\dfrac{1}{M}\sum\limits_{1\leq i\leq
n}m_{i}\Delta\boldsymbol{q}_{i}\right)-\dfrac{1}{M}\sum\limits_{1\leq i\leq
n}\Delta(t\boldsymbol{p}_{i})-\tau\overline{\boldsymbol{p}_{i}}\end{pmatrix}$
$\displaystyle=\dfrac{1}{\tau}\begin{pmatrix}\sum\limits_{1\leq i\leq
n}\Delta\left(\dfrac{\boldsymbol{p}_{i}^{T}\boldsymbol{p}_{i}}{2m_{i}}\right)+\sum\limits_{1\leq
i<j\leq n}\Delta V_{ij}\\\ \Delta\boldsymbol{P}\\\ \Delta\boldsymbol{L}\\\
\Delta\left(\dfrac{1}{M}\left(\sum\limits_{1\leq i\leq
n}m_{i}\boldsymbol{q}_{i}\right)-\dfrac{P}{M}t\right)+\dfrac{\tau}{M}\overline{\boldsymbol{P}}\end{pmatrix}=\dfrac{1}{\tau}\begin{pmatrix}\Delta
H\\\ \Delta\boldsymbol{P}\\\ \Delta\boldsymbol{L}\\\
\Delta\boldsymbol{C}\end{pmatrix}+\begin{pmatrix}0\\\ \boldsymbol{0}\\\
\boldsymbol{0}\\\
\dfrac{\overline{\boldsymbol{P}}}{M}\end{pmatrix}=D_{t}^{\tau}\boldsymbol{\psi}-\partial_{t}^{\tau}\boldsymbol{\psi}.$

Thus, the discretization (15) is indeed a conservative scheme for (13).
Moreover, since $\overline{\boldsymbol{p}_{i}},\overline{\boldsymbol{q}_{i}}$
and $\overline{q_{ij}}$ are symmetric under the permutation of
$k\leftrightarrow k+1$, (15) is a symmetric scheme provided that
$\dfrac{\Delta V_{ij}}{\Delta q_{ij}}$ is symmetric under the permutation of
$k\leftrightarrow k+1$. Indeed, this is true by the definition of the divided
difference of $\dfrac{\Delta V_{ij}}{\Delta q_{ij}}$ and that $V_{ij}$ is a
function of $q_{ij}$.

###  A.3 Point vortex problem: Planar case

Here, we verify the conserved quantities of (17) and show the derived scheme
is conservative and symmetric using Lemma 1.

First, we verify that (18) are indeed conserved quantities of (17) by defining
the conserved vector $\boldsymbol{\psi}$ as,

$\displaystyle\boldsymbol{\psi}(\boldsymbol{x},\boldsymbol{y}):=\begin{pmatrix}\boldsymbol{P}(\boldsymbol{x},\boldsymbol{y})\\\
L(\boldsymbol{x},\boldsymbol{y})\\\
H(\boldsymbol{x},\boldsymbol{y})\end{pmatrix}.$

Then by (4a), the associated $4\times(2n)$ multiplier matrix is given by,

$\displaystyle\Lambda(\boldsymbol{x},\boldsymbol{y}):=\begin{pmatrix}\begin{bmatrix}\Gamma_{i}\end{bmatrix}_{1\leq
i\leq n}^{T}&\begin{bmatrix}0\end{bmatrix}_{1\leq i\leq n}^{T}\\\
\begin{bmatrix}0\end{bmatrix}_{1\leq i\leq
n}^{T}&\begin{bmatrix}\Gamma_{i}\end{bmatrix}_{1\leq i\leq n}^{T}\\\
\begin{bmatrix}2\Gamma_{i}x_{i}\end{bmatrix}_{1\leq i\leq
n}^{T}&\begin{bmatrix}2\Gamma_{i}y_{i}\end{bmatrix}_{1\leq i\leq n}^{T}\\\
\begin{bmatrix}-\dfrac{1}{2\pi}\Gamma_{i}\sum\limits_{j=1,j\neq
i}^{n}\Gamma_{j}\dfrac{x_{ij}}{r_{ij}^{2}}\end{bmatrix}_{1\leq i\leq
n}^{T}&\begin{bmatrix}-\dfrac{1}{2\pi}\Gamma_{i}\sum\limits_{j=1,j\neq
i}^{n}\Gamma_{j}\dfrac{y_{ij}}{r_{ij}^{2}}\end{bmatrix}_{1\leq i\leq
n}^{T}\end{pmatrix}$

So (4b) is satisfied since,

$\displaystyle\Lambda(\boldsymbol{x},\boldsymbol{y})\boldsymbol{f}(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\begin{pmatrix}\dfrac{1}{2\pi}\sum\limits_{1\leq i\neq j\leq
n}\Gamma_{i}\Gamma_{j}\dfrac{y_{ij}}{r_{ij}^{2}}\\\
-\dfrac{1}{2\pi}\sum\limits_{1\leq i\neq j\leq
n}\Gamma_{i}\Gamma_{j}\dfrac{x_{ij}}{r_{ij}^{2}}\\\
\dfrac{1}{\pi}\sum\limits_{1\leq i\neq j\leq
n}\Gamma_{i}\Gamma_{j}\dfrac{y_{ij}x_{i}-x_{ij}y_{i}}{r_{ij}^{2}}\\\
-\dfrac{1}{2\pi}\sum\limits_{1\leq i\leq
n}\Gamma_{i}\sum\limits_{j,l=1,j,l\neq
i}^{n}\Gamma_{j}\Gamma_{l}\left(\dfrac{y_{ij}}{r_{ij}^{2}}\dfrac{x_{il}}{r_{il}^{2}}-\dfrac{x_{ij}}{r_{ij}^{2}}\dfrac{y_{il}}{r_{il}^{2}}\right)\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\dfrac{1}{2\pi}\sum\limits_{1\leq i<j\leq
n}\Gamma_{i}\Gamma_{j}\left(\dfrac{y_{ij}}{r_{ij}^{2}}+\dfrac{y_{ji}}{r_{ji}^{2}}\right)\\\
-\dfrac{1}{2\pi}\sum\limits_{1\leq i<j\leq
n}\Gamma_{i}\Gamma_{j}\left(\dfrac{x_{ij}}{r_{ij}^{2}}+\dfrac{x_{ji}}{r_{ji}^{2}}\right)\\\
\dfrac{1}{\pi}\sum\limits_{1\leq i<j\leq
n}\Gamma_{i}\Gamma_{j}\left(\dfrac{y_{ij}x_{i}-x_{ij}y_{i}}{r_{ij}^{2}}+\dfrac{y_{ji}x_{j}-x_{ji}y_{j}}{r_{ji}^{2}}\right)\\\
-\dfrac{1}{2\pi}\sum\limits_{1\leq i\leq
n}\Gamma_{i}\left(\sum\limits_{j,l=1,j,l\neq
i}^{n}\Gamma_{j}\Gamma_{l}\dfrac{y_{ij}}{r_{ij}^{2}}\dfrac{x_{il}}{r_{il}^{2}}-\sum\limits_{j,l=1,j,l\neq
i}^{n}\Gamma_{j}\Gamma_{l}\dfrac{x_{ij}}{r_{ij}^{2}}\dfrac{y_{il}}{r_{il}^{2}}\right)\end{pmatrix}=\boldsymbol{0},$

where the last equality follows from

$\displaystyle r_{ji}$ $\displaystyle=r_{ij},$ $\displaystyle y_{ji}$
$\displaystyle=-y_{ij},$ $\displaystyle x_{ji}$ $\displaystyle=-x_{ij},$
$\displaystyle y_{ji}x_{j}-x_{ji}y_{j}$
$\displaystyle=y_{i}x_{j}-x_{i}y_{j}=-(y_{ij}x_{i}-x_{ij}y_{i}).$

Similar to the $n$-body problem, we will propose consistent choices of of
$D_{t}^{\tau}\boldsymbol{x}$, $D_{t}^{\tau}\boldsymbol{\psi}$,
$\partial_{t}^{\tau}\boldsymbol{\psi}$, $\Lambda^{\tau}$ and
$\boldsymbol{f}^{\tau}$ and verify that both conditions (9a) and (9b) are
satisfied. Analogous to the $n$-body problem, we define

$\displaystyle
D_{t}^{\tau}\boldsymbol{x}:=\dfrac{1}{\tau}\begin{pmatrix}\Delta\boldsymbol{x}\\\
\Delta\boldsymbol{y}\end{pmatrix},\hskip
14.22636ptD_{t}^{\tau}\boldsymbol{\psi}:=\dfrac{1}{\tau}\begin{pmatrix}\Delta\boldsymbol{P}\\\
\Delta L\\\ \Delta H\end{pmatrix},\hskip
14.22636pt\partial_{t}^{\tau}\boldsymbol{\psi}:=\boldsymbol{0}.$

Let us define the discrete multiplier $\Lambda^{\tau}$ and the discrete right
hand side $\boldsymbol{f}^{\tau}$ as,

$\displaystyle\Lambda^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{y}^{k+1},\boldsymbol{x}^{k},\boldsymbol{y}^{k}):=$
$\displaystyle\hskip
28.45274pt\begin{pmatrix}\begin{bmatrix}\Gamma_{i}\end{bmatrix}_{1\leq i\leq
n}^{T}&\begin{bmatrix}0\end{bmatrix}_{1\leq i\leq n}^{T}\\\
\begin{bmatrix}0\end{bmatrix}_{1\leq i\leq
n}^{T}&\begin{bmatrix}\Gamma_{i}\end{bmatrix}_{1\leq i\leq n}^{T}\\\
\begin{bmatrix}2\Gamma_{i}\overline{x_{i}}\end{bmatrix}_{1\leq i\leq
n}^{T}&\begin{bmatrix}2\Gamma_{i}\overline{y_{i}}\end{bmatrix}_{1\leq i\leq
n}^{T}\\\ \begin{bmatrix}-\dfrac{1}{2\pi}\Gamma_{i}\sum\limits_{j=1,j\neq
i}^{n}\Gamma_{j}\dfrac{\overline{x_{ij}}}{(r_{ij}^{k})^{2}}g\left(z_{ij}\right)\end{bmatrix}_{1\leq
i\leq n}^{T}&\begin{bmatrix}-\dfrac{1}{2\pi}\Gamma_{i}\sum\limits_{j=1,j\neq
i}^{n}\Gamma_{j}\dfrac{\overline{y_{ij}}}{(r_{ij}^{k})^{2}}g\left(z_{ij}\right)\end{bmatrix}_{1\leq
i\leq n}^{T}\end{pmatrix},$
$\displaystyle\boldsymbol{f}^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{y}^{k+1},\boldsymbol{x}^{k},\boldsymbol{y}^{k}):=\begin{pmatrix}-\begin{bmatrix}\dfrac{1}{2\pi}\sum\limits_{1\leq
j\leq n,j\neq
i}\dfrac{\Gamma_{j}\overline{y_{ij}}}{(r_{ij}^{k})^{2}}g\left(z_{ij}\right)\end{bmatrix}_{1\leq
i\leq n}\\\ \begin{bmatrix}\dfrac{1}{2\pi}\sum\limits_{1\leq j\leq n,j\neq
i}\dfrac{\Gamma_{j}\overline{x_{ij}}}{(r_{ij}^{k})^{2}}g\left(z_{ij}\right)\end{bmatrix}_{1\leq
i\leq n}\end{pmatrix},$

where for brevity we have denoted $\displaystyle
z_{ij}:=\left(\dfrac{r_{ij}^{k+1}}{r_{ij}^{k}}\right)^{2}$. Since
$\displaystyle g(z_{ij})\rightarrow 1$ as $z_{ij}\rightarrow 1$ when
$\tau\rightarrow 0$, $\Lambda^{\tau}$ and $\boldsymbol{f}^{\tau}$ are
consistent to $\Lambda$ and $\boldsymbol{f}$. Next, we verify condition (9b).

$\displaystyle\Lambda^{\tau}\boldsymbol{f}^{\tau}$
$\displaystyle=\begin{pmatrix}\dfrac{1}{2\pi}\sum\limits_{1\leq i\neq j\leq
n}\Gamma_{i}\Gamma_{j}\dfrac{\overline{y_{ij}}}{(r_{ij}^{k})^{2}}g\left(z_{ij}\right)\\\
-\dfrac{1}{2\pi}\sum\limits_{1\leq i\neq j\leq
n}\Gamma_{i}\Gamma_{j}\dfrac{\overline{x_{ij}}}{(r_{ij}^{k})^{2}}g\left(z_{ij}\right)\\\
\dfrac{1}{\pi}\sum\limits_{1\leq i\neq j\leq
n}\Gamma_{i}\Gamma_{j}\dfrac{\overline{y_{ij}}\hskip
2.84526pt\overline{x_{i}}-\overline{x_{ij}}\hskip
2.84526pt\overline{y_{i}}}{(r_{ij}^{k})^{2}}g\left(z_{ij}\right)\\\
-\dfrac{1}{2\pi}\sum\limits_{1\leq i\leq
n}\Gamma_{i}\sum\limits_{j,l=1,j,l\neq
i}^{n}\Gamma_{j}\Gamma_{l}\left(\dfrac{\overline{y_{ij}}}{(r_{ij}^{k})^{2}}\dfrac{\overline{x_{il}}}{(r_{il}^{k})^{2}}-\dfrac{\overline{x_{ij}}}{(r_{ij}^{k})^{2}}\dfrac{\overline{y_{il}}}{(r_{il}^{k})^{2}}\right)g\left(z_{ij}\right)g\left(z_{il}\right)\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\dfrac{1}{2\pi}\sum\limits_{1\leq i<j\leq
n}\Gamma_{i}\Gamma_{j}\left(\dfrac{\overline{y_{ij}}}{(r_{ij}^{k})^{2}}g\left(z_{ij}\right)+\dfrac{\overline{y_{ji}}}{(r_{ji}^{k})^{2}}g\left(z_{ji}\right)\right)\\\
-\dfrac{1}{2\pi}\sum\limits_{1\leq i<j\leq
n}\Gamma_{i}\Gamma_{j}\left(\dfrac{\overline{x_{ij}}}{(r_{ij}^{k})^{2}}g\left(z_{ij}\right)+\dfrac{\overline{x_{ji}}}{(r_{ji}^{k})^{2}}g\left(z_{ji}\right)\right)\\\
\dfrac{1}{\pi}\sum\limits_{1\leq i<j\leq
n}\Gamma_{i}\Gamma_{j}\left(\dfrac{\overline{y_{ij}}\hskip
2.84526pt\overline{x_{i}}-\overline{x_{ij}}\hskip
2.84526pt\overline{y_{i}}}{(r_{ij}^{k})^{2}}g\left(z_{ij}\right)+\dfrac{\overline{y_{ji}}\hskip
2.84526pt\overline{x_{j}}-\overline{x_{ji}}\hskip
2.84526pt\overline{y_{j}}}{(r_{ji}^{k})^{2}}g\left(z_{ji}\right)\right)\\\
-\dfrac{1}{2\pi}\sum\limits_{1\leq i\leq
n}\Gamma_{i}\left(\sum\limits_{j,l=1,j,l\neq
i}^{n}\Gamma_{j}\Gamma_{l}\dfrac{\overline{y_{ij}}}{(r_{ij}^{k})^{2}}\dfrac{\overline{x_{il}}}{(r_{il}^{k})^{2}}g\left(z_{ij}\right)g\left(z_{il}\right)-\sum\limits_{j,l=1,j,l\neq
i}^{n}\Gamma_{j}\Gamma_{l}\dfrac{\overline{x_{ij}}}{(r_{ij}^{k})^{2}}\dfrac{\overline{y_{il}}}{(r_{il}^{k})^{2}}g\left(z_{ij}\right)g\left(z_{il}\right)\right)\end{pmatrix}$
$\displaystyle=\boldsymbol{0},$

where the last equality follows from

$\displaystyle r_{ji}^{k}$ $\displaystyle=r_{ij}^{k},$ $\displaystyle z_{ji}$
$\displaystyle=z_{ij},$ $\displaystyle\overline{y_{ji}}$
$\displaystyle=-\overline{y_{ij}},$ $\displaystyle\overline{x_{ji}}$
$\displaystyle=-\overline{x_{ij}},$ $\displaystyle\overline{y_{ji}}\hskip
2.84526pt\overline{x_{j}}-\overline{x_{ji}}\hskip 2.84526pt\overline{y_{j}}$
$\displaystyle=\overline{y_{i}}\hskip
2.84526pt\overline{x_{j}}-\overline{x_{i}}\hskip
2.84526pt\overline{y_{j}}=-(\overline{y_{ij}}\hskip
2.84526pt\overline{x_{i}}-\overline{x_{ij}}\hskip 2.84526pt\overline{y_{i}}).$

To verify condition (9a), it follows from direct computation or divided
difference calculus that

$\displaystyle\Delta(x_{i}^{2}+y_{i}^{2})$
$\displaystyle=2\overline{x_{i}}\Delta x_{i}+2\overline{y_{i}}\Delta y_{i},$
$\displaystyle\Delta(r_{ij})^{2}$ $\displaystyle=2\overline{x_{i}}\Delta
x_{i}+2\overline{y_{i}}\Delta y_{i},$
$\displaystyle\dfrac{g\left(z_{ij}\right)}{(r_{ij}^{k})^{2}}$
$\displaystyle=\dfrac{2(\log r_{ij}^{k+1}-2\log
r_{ij}^{k})}{\Delta(r_{ij})^{2}},$ $\displaystyle\Delta\log(r_{ij})$
$\displaystyle=(\log r_{ij}^{k+1}-\log
r_{ij}^{k})\dfrac{2\overline{x_{ij}}\Delta x_{i}+2\overline{y_{ij}}\Delta
y_{i}}{\Delta(r_{ij})^{2}}=\left(\overline{x_{ij}}\Delta
x_{i}+\overline{y_{ij}}\Delta
y_{i}\right)\dfrac{g\left(z_{ij}\right)}{(r_{ij}^{k})^{2}}.$

Combining with linearity of $\Delta$, condition (9a) is satisfied since,

$\displaystyle\Lambda^{\tau}D_{t}^{\tau}\boldsymbol{x}$
$\displaystyle=\frac{1}{\tau}\begin{pmatrix}\sum\limits_{1\leq i\leq
n}\Gamma_{i}\Delta x_{i}\\\ \sum\limits_{1\leq i\leq n}\Gamma_{i}\Delta
y_{i}\\\ \sum\limits_{1\leq i\leq n}\Gamma_{i}(2\overline{x_{i}}\Delta
x_{i}+2\overline{y_{i}}\Delta y_{i})\\\ -\dfrac{1}{2\pi}\sum\limits_{1\leq
i\neq j\leq n}\Gamma_{i}\Gamma_{j}\left(\overline{x_{ij}}\Delta
x_{i}+\overline{y_{ij}}\Delta
y_{i}\right)\dfrac{g\left(z_{ij}\right)}{(r_{ij}^{k})^{2}}\end{pmatrix}$
$\displaystyle=\frac{1}{\tau}\begin{pmatrix}\Delta\left(\sum\limits_{1\leq
i\leq n}\Gamma_{i}x_{i}\right)\\\ \Delta\left(\sum\limits_{1\leq i\leq
n}\Gamma_{i}y_{i}\right)\\\ \Delta\left(\sum\limits_{1\leq i\leq
n}\Gamma_{i}(x_{i}^{2}+y_{i}^{2})\right)\\\ -\dfrac{1}{2\pi}\sum\limits_{1\leq
i<j\leq n}\Gamma_{i}\Gamma_{j}\left(\overline{x_{ij}}\Delta
x_{ij}+\overline{y_{ij}}\Delta
y_{ij}\right)\dfrac{g\left(z_{ij}\right)}{(r_{ij}^{k})^{2}}\end{pmatrix}=\frac{1}{\tau}\begin{pmatrix}\Delta\boldsymbol{P}\\\
\Delta L\\\ \Delta
H\end{pmatrix}=D_{t}^{\tau}\boldsymbol{\psi}-\partial_{t}^{\tau}\boldsymbol{\psi}.$

Thus, the scheme (19) is indeed conservative for (17). Moreover, (19) is
symmetric since $\overline{x}_{ij},\overline{y}_{ij}$ and
$\dfrac{1}{(r_{ij}^{k})^{2}}g(z_{ij})=\dfrac{1}{(r_{ij}^{k})^{2}}g\left(\left(\dfrac{r_{ij}^{k+1}}{r_{ij}^{k}}\right)^{2}\right)$
is symmetric under the permutation of $k\leftrightarrow k+1$, by Lemma 1 with
$h(\boldsymbol{x}^{k})=(r_{ij}^{k})^{2}$.

###  A.4 Point vortex problem: Spherical case

Finally, we verify the conserved quantities of (20) and show that the derived
scheme (22) is conservative and symmetric using Lemma 1.

The conserved vector $\boldsymbol{\psi}$ of (21) for point vortex problem on
the unit sphere (20) is

$\displaystyle\boldsymbol{\psi}(\boldsymbol{x}):=\begin{pmatrix}\boldsymbol{P}(\boldsymbol{x})\\\
H(\boldsymbol{x})\end{pmatrix}.$

Using (4a), the associated $4\times(3n)$ multiplier matrix is given by,

$\displaystyle\Lambda(\boldsymbol{x}):=\begin{pmatrix}\begin{bmatrix}\Gamma_{i}I\end{bmatrix}_{1\leq
i\leq n}^{T}\\\ \begin{bmatrix}\dfrac{1}{4\pi}\Gamma_{i}\sum\limits_{j=1,j\neq
i}^{n}\Gamma_{j}\dfrac{{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j}}{1-{\boldsymbol{x}}_{i}\cdot{\boldsymbol{x}}_{j}}\end{bmatrix}_{1\leq
i\leq n}^{T}\end{pmatrix}$

where $I$ is the $3\times 3$ identity matrix. So (4b) is satisfied since,

$\displaystyle\Lambda(\boldsymbol{x})\boldsymbol{f}(\boldsymbol{x})$
$\displaystyle=\begin{pmatrix}\dfrac{1}{4\pi}\sum\limits_{1\leq i\neq j\leq
n}\Gamma_{i}\Gamma_{j}\dfrac{{\boldsymbol{x}}_{j}\times{\boldsymbol{x}}_{i}}{1-{\boldsymbol{x}}_{j}\cdot{\boldsymbol{x}}_{i}}\\\
\dfrac{1}{4\pi}\Gamma_{i}\sum\limits_{j=1,j\neq
i}^{n}\Gamma_{j}\dfrac{{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j}}{1-{\boldsymbol{x}}_{i}\cdot{\boldsymbol{x}}_{j}}\cdot\dfrac{1}{4\pi}\sum\limits_{k=1,k\neq
i}^{n}\Gamma_{k}\dfrac{{\boldsymbol{x}}_{k}\times{\boldsymbol{x}}_{i}}{1-{\boldsymbol{x}}_{i}\cdot{\boldsymbol{x}_{k}}}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\dfrac{1}{4\pi}\sum\limits_{1\leq i<j\leq
n}\Gamma_{i}\Gamma_{j}\left(\dfrac{{\boldsymbol{x}}_{j}\times{\boldsymbol{x}}_{i}}{1-{\boldsymbol{x}}_{j}\cdot{\boldsymbol{x}}_{i}}+\dfrac{{\boldsymbol{x}}_{i}\times{\boldsymbol{x}}_{j}}{1-{\boldsymbol{x}}_{j}\cdot{\boldsymbol{x}}_{i}}\right)\\\
\dfrac{1}{16\pi^{2}}\sum\limits_{i=1}^{n}\Gamma_{i}\sum\limits_{1\leq j<k\leq
n}\Gamma_{j}\Gamma_{k}\left(\dfrac{{\boldsymbol{x}}_{j}\cdot({\boldsymbol{x}}_{k}\times{\boldsymbol{x}}_{i})+{\boldsymbol{x}}_{k}\cdot({\boldsymbol{x}}_{j}\times{\boldsymbol{x}}_{i})}{(1-{\boldsymbol{x}}_{i}\cdot{\boldsymbol{x}}_{j})(1-{\boldsymbol{x}}_{i}\cdot{\boldsymbol{x}}_{k})}\right)\end{pmatrix}=\boldsymbol{0},$

where the last equality follows from the antisymmetry of the cross product and
the antisymmetry of the scalar triple product.

As before, we choose $D_{t}^{\tau}\boldsymbol{x}$,
$D_{t}^{\tau}\boldsymbol{\psi}$, $\partial_{t}^{\tau}\boldsymbol{\psi}$,
$\Lambda^{\tau}$ and $\boldsymbol{f}^{\tau}$ such that both conditions (9a)
and (9b) are satisfied. We start by defining

$\displaystyle
D_{t}^{\tau}\boldsymbol{x}:=\dfrac{\Delta\boldsymbol{x}}{\tau},\hskip
14.22636ptD_{t}^{\tau}\boldsymbol{\psi}:=\dfrac{1}{\tau}\begin{pmatrix}\Delta\boldsymbol{P}\\\
\Delta H\end{pmatrix},\hskip
14.22636pt\partial_{t}^{\tau}\boldsymbol{\psi}:=\boldsymbol{0}.$

The discrete multiplier $\Lambda^{\tau}$ and the discrete right hand side
$\boldsymbol{f}^{\tau}$ are defined as,

$\displaystyle\Lambda^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{x}^{k}):=\begin{pmatrix}\begin{bmatrix}\Gamma_{i}I\end{bmatrix}_{1\leq
i\leq n}^{T}\\\ \begin{bmatrix}\dfrac{1}{4\pi}\Gamma_{i}\sum\limits_{j=1,j\neq
i}^{n}\Gamma_{j}\dfrac{\overline{{\boldsymbol{x}}_{ij}}}{1-{\boldsymbol{x}}_{i}^{k}\cdot{\boldsymbol{x}}_{j}^{k}}g(z_{ij})\end{bmatrix}_{1\leq
i\leq n}^{T}\end{pmatrix},$
$\displaystyle\boldsymbol{f}^{\tau}(\boldsymbol{x}^{k+1},\boldsymbol{x}^{k}):=\dfrac{1}{4\pi}\sum\limits_{1\leq
j\leq n,j\neq
i}\Gamma_{j}\dfrac{\overline{{\boldsymbol{x}}_{j}}\times\overline{{\boldsymbol{x}}_{i}}}{1-{{\boldsymbol{x}}_{i}^{k}}\cdot{{\boldsymbol{x}}_{j}^{k}}}g\left(z_{ij}\right),$

where for brevity we denoted $\displaystyle
z_{ij}:=\dfrac{1-{\boldsymbol{x}}_{i}^{k+1}\cdot{\boldsymbol{x}}_{j}^{k+1}}{1-{\boldsymbol{x}}_{i}^{k}\cdot{\boldsymbol{x}}_{j}^{k}}$.
As in the planar case, $\displaystyle g(z_{ij})\rightarrow 1$ as
$z_{ij}\rightarrow 1$ when $\tau\rightarrow 0$ and so both $\Lambda^{\tau}$
and $\boldsymbol{f}^{\tau}$ are consistent to $\Lambda$ and $\boldsymbol{f}$,
respectively.

Similar to the planar case, it is readily verified that we have
$\Lambda^{\tau}\boldsymbol{f}^{\tau}=\boldsymbol{0}$ and thus it remains to
check condition (9a). Indeed, this condition is satisfied as

$\displaystyle\Lambda^{\tau}D_{t}^{\tau}\boldsymbol{x}$
$\displaystyle=\frac{1}{\tau}\begin{pmatrix}\sum\limits_{1\leq i\leq
n}\Gamma_{i}\Delta{\boldsymbol{x}}_{i}\\\ \dfrac{1}{4\pi}\sum\limits_{1\leq
i\neq j\leq
n}\Gamma_{i}\Gamma_{j}\dfrac{\overline{{\boldsymbol{x}}_{ij}}\cdot\Delta{\boldsymbol{x}}_{i}}{1-{\boldsymbol{x}}_{i}^{k}\cdot{\boldsymbol{x}}_{j}^{k}}g\left(z_{ij}\right)\end{pmatrix}$
$\displaystyle=\frac{1}{\tau}\begin{pmatrix}\Delta\left(\sum\limits_{1\leq
i\leq n}\Gamma_{i}{\boldsymbol{x}}_{i}\right)\\\
\dfrac{1}{4\pi}\sum\limits_{1\leq i<j\leq
n}\Gamma_{i}\Gamma_{j}\dfrac{\overline{{\boldsymbol{x}}_{ij}}\cdot\Delta{\boldsymbol{x}}_{ij}}{1-{\boldsymbol{x}}_{i}^{k}\cdot{\boldsymbol{x}}_{j}^{k}}g\left(z_{ij}\right)\end{pmatrix}=\frac{1}{\tau}\begin{pmatrix}\Delta\boldsymbol{P}\\\
\Delta
H\end{pmatrix}=D_{t}^{\tau}\boldsymbol{\psi}-\partial_{t}^{\tau}\boldsymbol{\psi}.$

Hence, the discretization (22) is conservative for (20). Finally, (22) is
symmetric since
$\overline{\boldsymbol{x}}_{j}\times\overline{\boldsymbol{x}}_{i}$, and
$\dfrac{1}{1-{\boldsymbol{x}}_{i}^{k}\cdot{\boldsymbol{x}}_{j}^{k}}g(z_{ij})=\dfrac{1}{1-{\boldsymbol{x}}_{i}^{k}\cdot{\boldsymbol{x}}_{j}^{k}}g\left(\dfrac{1-{\boldsymbol{x}}_{i}^{k+1}\cdot{\boldsymbol{x}}_{j}^{k+1}}{1-{\boldsymbol{x}}_{i}^{k}\cdot{\boldsymbol{x}}_{j}^{k}}\right)$
is symmetric under the permutation of $k\leftrightarrow k+1$, by Lemma 1 with
$h(\boldsymbol{x}^{k})=1-{\boldsymbol{x}}_{i}^{k}\cdot{\boldsymbol{x}}_{j}^{k}$.