# Port-Hamiltonian structure of interacting particle systems and its mean-
field limit

Birgit Jacob111Research Group Functional Analysis<EMAIL_ADDRESS>IMACM, School of Mathematics and Natural Sciences,
University of Wuppertal, Germany  Claudia Totzeck222Research Group
Optimization<EMAIL_ADDRESS>IMACM, School of Mathematics and
Natural Sciences,
University of Wuppertal, Germany

(January 2023)

###### Abstract

We derive a minimal port-Hamiltonian formulation of a general class of
interacting particle systems driven by alignment and potential-based force
dynamics which include the Cucker-Smale model with potential interaction and
the second order Kuramoto model. The port-Hamiltonian structure allows to
characterize conserved quantities such as Casimir functions as well as the
long-time behaviour using a LaSalle-type argument on the particle level. It is
then shown that the port-Hamiltonian structure is preserved in the mean-field
limit and an analogue of the LaSalle invariance principle is studied in the
space of probability measures equipped with the 2-Wasserstein-metric. The
results on the particle and mean-field limit yield a new perspective on
uniform stability of general interacting particle systems. Moreover, as the
minimal port-Hamiltonian formulation is closed we identify the ports of the
subsystems which admit generalized mass-spring-damper structure modelling the
binary interaction of two particles. Using the information of ports we discuss
the coupling of difference species in a port-Hamiltonian preserving manner.

AMS classification: 37K45, 82C22, 93A16.

Keywords: Port-Hamiltonian systems, interacting particle systems, mean-field
limit, long-time behaviour

##  1 Introduction

Since the seminal works by Reynolds [33] in 1987, Vicsek, Czirók, Ben-Jacob,
Cohen, Shochet [40] in 1995, and Cucker, Smale [19] in 2007, mathematical
modelling of interacting particle systems and the structural analysis of these
models attracts the attention of researchers from applied mathematics. One of
the fascinating aspects is that simple interaction rules imposed for the
binary interaction of two particles lead to collective behaviour of the whole
crowd. In fact, often slight parameter changes can turn dynamics of ordinary
differential equations (ODEs), where the particles form rings that resemble
the milling of birds into clumps [20].

As the analysis of these pattern formation is difficult on the particle level,
where the position and velocity information of each member of the crowds is
explicitly captured by the equations, more abstract formulations of the
dynamics were proposed. Sending the number of particles to infinity leads to
the so-called mean-field formulation of the crowd [24]. Here, the exact
position and velocity information is averaged and only the probability of
finding a particle at a certain time in a certain position with a certain
velocity is described. The binary interaction structure on the particle level
turns into an convolution that yields a nonlinear, nonlocal partial
differential equation (PDE) as evolution equation on the mean-field level.
Instead of a huge system of ODEs the mean-field equation requires the solution
of a high-dimensional PDE. To further reduce the dimensionality of the
problem, also hydro-dynamic descriptions were proposed, there the velocity
information is averaged over space, leading to a coupled PDE-system that is
only space-dependent [9]. In contrast to the passage from mean-field level to
hydro-dynamic formulations, which in general require formal closure relations,
the rigorous relationship of particle and mean-field limit is well-understood,
see for example [24]. For an detailed overview of different models on the
various scales we refer to [34].

Many interacting particle systems are driven by two mechanisms: alignment in
velocity and attraction/repulsion in space [12]. Alignment is the main
component of (generalized) Cucker-Smale dynamics and leads to bird-like
behaviour. Many different interaction kernels were proposed to fit the model
closer to reality, see for example [1, 11, 18, 28]. On the other hand, there
are models that incorporate only attraction/repulsion forces leading to a
three-phase interaction behaviour of long-range attraction, short-range
repulsion and a mid-range comfort zone, where no interaction forces are
present [3, 20, 8]. These models often admit a gradient-structure, that means
the binary interaction forces are gradients of a prescribed interaction
potential like the Morse potential [9] and the force acting on one particle is
the average of all these binary interaction forces.

The gradient-structure opens the toolbox of gradient flows and large deviation
theory to study stability and long-time behaviour of the systems [4]. In the
absence of a gradient structure, the analysis of the alignment part of the
dynamics is generally more complex. However, starting with [13] stability
results for different setups in the Cucker-Smale context were established [29,
31, 25, 17, 16, 23, 5, 32, 2]. In the literature stable states of alignment
dynamics are often called flocking, clustering or consensus solutions.

In [8] flocking behaviour of a three-zone model was first investigated with
the help of the energy of the system. Although we address a similar question,
in this contribution we propose a novel viewpoint on interacting particle
dynamics, namely a port-Hamiltonian one. Building on the discussion of [27],
where the port-Hamiltonian system (PHS) formulation of a Cucker-Smale dynamics
with repulsion and attraction is interpreted as generalized mass-spring-damper
system, we propose a minimal PHS representation of interacting particle
systems that does not require an increase of the phase space. Indeed, the
formulation in [27] introduces the relative positions of all particles as new
variables, thereby increasing the state space dimension from $2Nd$ to
$N(N-1)d/2+Nd$ where $N$ denotes the number of particles and $d$ the space
dimension. In the following, the port-Hamiltonian reformulation preserves the
state space dimension which allows in particular to pass to the mean-field
limit.

The port-Hamiltonian formulation of interacting particles opens the door to
the well-established theory of PHS for finite-dimensional systems, see [39,
22, 21]. Well-known facts from PHS theory include the property that PHS are
closed under network interconnection, that is, coupling of port-Hamiltonian
systems again leads to a port-Hamiltonian system. Furthermore, the port-
Hamiltonian approach is suitable for the investigation of the qualitative
solution behavior such as asymptotic stability and control questions, as it
provides an energy balance. Moreover, the port-Hamiltonian structure allows to
identify several conserved quantities namely the Hamiltonian and the so-called
Casimir functions. The mean-field limit yields a connection to nonlinear,
nonlocal infinite-dimensional PHS systems which are yet less explored.

Our main contributions are: the port-Hamiltonian formulation of an interacting
particle systems with same state space dimension; structure preserving mean-
field limit; characterization of long-time behavior on the particle and mean-
field level and the characterization of conserved quantities such as the
Hamiltonian and Casimir functions. This yields a new perspective on uniform
stability and the dissipativity of interacting particle systems using the PH
dissipativity inequality. Moreover, the idetification of ports allows for PH
structure preserving coupling of interacting particle systems of (different)
species.

The article is organized as follows: we recall some background information on
interacting particle systems in Section 2. Then we motivate derive the port-
Hamiltonian reformulation in Section 3, that is used to discuss the Casimir
function and stability properties of the systems based on LaSalle theory. In
Section 4 we dicuss the mean-field limit sending the number of particles
$N\rightarrow\infty$ and show that the PH structure is conserved in the limit.
Section 5 discusses the coupling of (different) species in a PH structure
preserving manner. The article concludes with a summary of the main ideas and
an outlook to future work.

##  2 Background on interacting particle systems

We recall the classical formulation of interacting particle systems in
position and velocity coordinates and the corresponding mean-field limit. Let
us consider $N\in\mathbb{N},N\geq 2$ interacting particles in space dimension
$d.$ We denote their positions by $x_{i}\colon[0,T]\rightarrow\mathbb{R}^{d}$
and their velocities by $v_{i}\colon[0,T]\rightarrow\mathbb{R}^{d}$ for
$i=1,\dots,N,$ respectively. We collect the position and velocity information
of all particles in the vectors $x=(x_{i})_{i=1}^{N}$ and
$v=(v_{i})_{i=1}^{N},$ respectively. The dynamics of the $i$-th particle is
given by

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}x_{i}$ $\displaystyle=v_{i},$ (1a)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}v_{i}$
$\displaystyle=\frac{1}{N}\sum_{j=1}^{N}\psi(|x_{j}-x_{i}|)(v_{j}-v_{i})-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}(x_{i}-x_{j})=:F_{i}(t,x,v),$
(1b) $\displaystyle x_{i}(0)$ $\displaystyle=\hat{x}_{i},\qquad
v_{i}(0)=\hat{v}_{i}.$ (1c)

Here, $\psi\colon\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}$ models the strength
of the velocity alignment and
$\mathcal{V}\colon\mathbb{R}^{d}\rightarrow\mathbb{R}$ denotes the potential
modelling the binary interactions among the particles. For the forces
resulting from the interactions we require that $\mathcal{V}$ is continuously
partial differentiable satisfying

$\nabla\mathcal{V}(x)=-\nabla\mathcal{V}(-x),\qquad x\in\mathbb{R}^{d}.$ (2)

This general class of interaction models contains well-known examples: [19],
[27], Morse-interactions (sheep flocks, double and single milling birds) as
proposed in [20], or herding dynamics [7].

To obtain the existence and uniqueness of a global solution to the particle
system by standard results from ODE theory, we make the following

######  Assumption 2.1.

  * (1)

$F_{i}$ is continuous on $[0,T]\times\mathbb{R}^{dN}\times R^{dN}$ for all
$i\in\\{1,\dots,N\\}$

  * (2)

For some $C>0$ it holds

$|F_{i}(t,x,v)|\leq C(1+|x|+|v|)\quad\text{ for all }\quad
i\in\\{1,\dots,N\\},t\in[0,T]\text{ and }x,v\in\mathbb{R}^{dN}$

  * (3)

For all $i\in\\{1,\dots,N\\}$, $F_{i}(t,x,v)$ is locally Lipschitz continuous
w.r.t. $x$ and $v$. In particular, for every compact set
$K\subset\mathbb{R}^{dN}\times\mathbb{R}^{dN}$ there exists some $L_{K}>0$
such that $F_{i}$ is on $K$ Lipschitz continuous w.r.t. $x$ and $v$ uniformly
for all $i\in\\{1,\dots,N\\}$.

The following proposition ensures that (1) admits unique solutions in
$\mathcal{C}^{1}([0,\infty),\mathbb{R}^{dN}\times\mathbb{R}^{dN})$ the set of
all continuously differentiable functions
$z:[0,\infty)\rightarrow\mathbb{R}^{dN}\times\mathbb{R}^{dN}$.

######  Proposition 2.2.

[36] Let Assumption 2.1 hold, then for every initial condition
$(\hat{x},\hat{v})\in\mathbb{R}^{dN}\times\mathbb{R}^{dN}$ there exists a
unique global solution
$(x,v)\in\mathcal{C}^{1}([0,\infty),\mathbb{R}^{dN}\times\mathbb{R}^{dN})$ to
(1).

As we are interested in the limiting behaviour as $N\rightarrow\infty$ we
introduce the Wasserstein or Monge-Kantorovich-Rubinstein metric as distance
measure between the particle and the density perspective. Let us denote by
$\mathcal{P}_{2}(\mathbb{R}^{d})$ the space of Borel probability measures on
$\mathbb{R}^{d}$ with finite $2$-nd moment. Equipping with the 2-Wasserstein
distance, makes $\mathcal{P}_{2}(\mathbb{R}^{d})$ a complete metric space. We
further denote $\mathcal{P}_{2}^{ac}(\mathbb{R}^{d})$ the subset of
$\mathcal{P}_{2}(\mathbb{R}^{d})$ containing probability measures with
Lebesgue density. For the sake of completeness we recall the 2-Wasserstein
distance:

$W_{2}^{2}(\mu,\nu):=\inf\limits_{\pi\in\Pi(\mu,\nu)}\biggl{\\{}\int_{\mathbb{R}^{d}}|x-y|^{2}\mathrm{d}\pi(x,y)\biggr{\\}},\qquad\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d}),$

where $\Pi(\mu,\nu)$ denotes the set of all Borel probability measures on
$\mathbb{R}^{2d}$ that have $\mu$ and $\nu$ as first and second marginals
respectively, i.e.

$\pi(B\times\mathbb{R}^{d})=\mu(B),\qquad\pi(\mathbb{R}^{d}\times
B)=\nu(B)\quad\text{for }B\in\mathcal{B}(\mathbb{R}^{d}).$

We emphasize that throughout the article we denote the integral of a function
$\varphi\in\mathcal{C}(\mathbb{R}^{d})$ with respect to a probability measure
$\mu\colon\mathbb{R}^{d}\rightarrow[0,1]$ by
$\int\varphi(x)\mu(x)\mathrm{d}x$, even if the probability measure is not
absolutely continuous with respect to the Lebesgue measure, and hence does not
have an associated density. Moreover, for evaluations of
$f\in\mathcal{C}([0,T],\mathcal{P}(\mathbb{R}^{d}\times\mathbb{R}^{d}))$ at
time $t$, position $x$ and velocity $v$ we write $f_{t}(x,v):=f(t,x,v)$. In
particular, we use this for the empirical measure of the dynamics (1), which
is given by

$f^{N}_{t}(x,v):=\sum_{i=1}^{N}\delta(x-x_{i}(t))\otimes\delta(v-v_{i}(t)).$

It is well-known that $f^{N}$ satisfies the so-called mean-field equation

$\displaystyle\partial_{t}f^{N}+\nabla_{x}\cdot(vf^{N})=\nabla_{v}\cdot\Big{(}\int\psi(|y-x|)(w-v)-\nabla\mathcal{V}(y-x)\mathrm{d}f^{N}(t,y,w)f^{N}\Big{)}$
(3)

in the weak sense. To be more precise, we consider the following notion of
solution, where $\mathcal{C}_{c}^{\infty}$ denotes the set of all infinitely
differentiable function with compact support.

######  Definition 2.3.

We call
$\mu\in\mathcal{C}([0,T],\mathcal{P}_{2}(\mathbb{R}^{d}\times\mathbb{R}^{d}))$
a weak measure solution of (3) with initial condition
$\hat{\mu}\in\mathcal{P}_{2}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ if and only
if for any test function
$h\in\mathcal{C}_{c}^{\infty}([0,T)\times\mathbb{R}^{d}\times\mathbb{R}^{d})$
we have

$\displaystyle\int_{0}^{T}\int_{\mathbb{R}^{d}}\big{(}\partial_{t}h_{t}+v\cdot\nabla_{x}h_{t}\big{)}\mathrm{d}\mu_{t}\mathrm{d}t+\int_{\mathbb{R}^{d}}h_{0}\mathrm{d}\hat{\mu}$
$\displaystyle=\int_{0}^{T}\int_{\mathbb{R}^{d}}\Big{(}\int\psi(|y-x|)(w-v)-\nabla\mathcal{V}(y-x)\mathrm{d}\mu_{t}(y,w)\cdot\nabla_{v}h\,\mathrm{d}\mu_{t}\mathrm{d}t.$

For notational convenience we introduce a short hand notation for the force
term, let

$F\ast\mu\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d},\quad\big{(}F\ast\mu\big{)}(x,v)=-\int\psi(|y-x|)(w-v)-\nabla\mathcal{V}(y-x)\mathrm{d}\mu(t,y,w).$

For the empirical measure $f^{N}$ it holds $\big{(}F\ast
f^{N})(x_{i},v_{i})=F_{i}(x,v),$ where $F_{i}$ was already defined in (1). In
particular, the assumption on $F_{i}$ in Assumption 2.1 imply the same
properties for $F\ast\mu$ for compactly supported measures
$\mu\in\mathcal{P}_{c}(\mathbb{R}^{d}\times\mathbb{R}^{d}).$ Moreover, we have
the following well-posedness result on the mean-field level.

######  Proposition 2.4.

[15] Let Assumption 2.1 hold and
$f_{0}\in\mathcal{P}_{c}(\mathbb{R}^{d}\times\mathbb{R}^{d}).$ Then there
exists a unique solution
$f\in\mathcal{C}([0,T],\mathcal{P}_{2}(\mathbb{R}^{d}\times\mathbb{R}^{d}))$
to (3) in the sense of Definition 2.3 with initial condition $f_{0}$.
Moreover, the solution remains compactly supported for all $t\in[0,T]$.

Moreover, in the limit $N\rightarrow\infty$ we have the well-known convergence
of the empirical measure to the solution of the PDE in Wasserstein sense, see
for example [24].

######  Proposition 2.5 (Dobrushin).

Let Assumption 2.1 hold, $\mu$ and $\mu^{\prime}$ be solutions to the
continuity equation (3) for compactly supported initial data $\hat{\mu}$,
$\hat{\mu}^{\prime}\in\mathcal{P}_{c}(\mathbb{R}^{d})$, respectively. Then,
there exists a positive constant $c$ such that

$W_{2}^{2}(\mu_{t},\mu_{t}^{\prime})\leq
W_{2}^{2}(\hat{\mu},\hat{\mu}^{\prime})e^{ct}\quad\text{for all\,
$t\in[0,T]$}.$

In the following we reformulate the ODE system as port-Hamiltonian system. The
port-Hamiltonian structure opens the door for an alternative investigation of
asymptotic flocking and uniform stability of general alignment-interaction
models, which is discussed in the seminal article [13] for the Cucker-Smale
model. Moreover, we show that the port-Hamiltonian structure is preserved
while passing to the mean-field limit and characterize Casimir functions of
the dynamics.

##  3 Port-Hamiltonian formulation

In this section, we derive two port-Hamiltonian formulations of the
interacting particle system. We emphasize that both port-Hamiltonian
formulations preserve the state space dimension. This is in contrast to [27]
where all relative positions between the particles are considered. First, we
discuss the reformulation of the system given in $(x,v)$ coordinates. Then we
present a variant that exploits the translational invariance of the systems.
This enables us to benefit from LaSalle’s invariance principle to analyse the
flocking behaviour.

Let us define $z=(x_{1},\dots,x_{N},v_{1},\dots,v_{N})$. The Hamiltonian of
the system (1) is given by the sum of kinetic and potential energy

$\mathcal{H}^{N}(z)=\frac{1}{2}\sum_{i=1}^{N}\Big{(}v_{i}^{\top}v_{i}+\frac{1}{N}\sum_{j=1}^{N}\mathcal{V}(x_{i}-x_{j})\Big{)},\quad
z=(x_{1},\dots,x_{N},v_{1},\dots,v_{N}).$

Then we calculate

$\displaystyle\frac{\partial\mathcal{H}^{N}(z)}{\partial
z}=\left(-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}(x_{j}-x_{1}),\ldots,-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}(x_{j}-x_{N}),v\right)^{\top}$

and obtain

$\frac{\mathrm{d}}{\mathrm{d}t}z=\left(\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}-\begin{pmatrix}0&0\\\
0&\Psi(z)\end{pmatrix}\right)\frac{\partial\mathcal{H}^{N}(z)}{\partial
z},\quad z(0)=z_{0},$ (4)

where

$\Psi(z)=\frac{1}{N}\begin{pmatrix}\displaystyle\sum_{j=1}^{N}\psi(|x_{j}-x_{1}|)&-\psi(|x_{2}-x_{1}|)&\ldots&-\psi(|x_{N}-x_{1}|)\\\
-\psi(|x_{1}-x_{2}|)&\displaystyle\sum_{j=1}^{N}\psi(|x_{j}-x_{2}|)&\ddots&\vdots\\\
\vdots&\ddots&\ddots&&\\\ &&&-\psi(|x_{N}-x_{N-1}|)\\\
-\psi(|x_{1}-x_{N}|)&\ldots&-\psi(|x_{N-1}-x_{N}|)&\displaystyle\sum_{j=1}^{N}\psi(|x_{j}-x_{N}|)\end{pmatrix}.$

The state space of the port-Hamiltonian system (4) is given by
$X:=\mathbb{R}^{Nd+Nd}$.

######  Remark.

For later use, we remark that $\Psi(z)$ is diagonally dominant with non
negative diagonal elements, hence positive semi-definite. Further,
$\\{v\in\mathbb{R}^{Nd}:v=(\tilde{v},\dots,\tilde{v})\text{ for some
}\tilde{v}\in\mathbb{R}^{d}\\}$ is contained in the kernel of $\Psi(z)$.
Moreover, if $\psi(x)>0$ for $x\in\mathbb{R}$, then the eigenspace
corresponding to the zero eigenvector is spanned by
$\\{v\in\mathbb{R}^{Nd}:v=(\tilde{v},\dots,\tilde{v})\text{ for some
}\tilde{v}\in\mathbb{R}^{d}\\}$.

As the positions of the particles appear only relatively in the dynamics,
interacting particle systems are translational invariant w.r.t. $x$. This
motivates us to consider the dynamics with center of mass shifted to zero. Let
us therefore define the center of mass by
$\bar{x}=\frac{1}{N}\sum_{j=1}^{N}x_{j}$ and consider
$z=(r,v)=(r_{1},\dots,r_{N},v_{1},\dots,v_{N})$ with $r_{i}=x_{i}-\bar{x}$
being the position relative to the center of mass. Note that the velocity of
the center of mass is conserved. Indeed, using
$\nabla\mathcal{V}(-r)=-\nabla\mathcal{V}(r)$, we find

$\displaystyle\sum_{i=1}^{N}\frac{\mathrm{d}}{\mathrm{d}t}v_{i}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}\psi(|r_{j}-r_{i}|)(v_{j}-v_{i})-\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}\nabla\mathcal{V}(r_{i}-r_{j})=0,$
(5)

Let us denote $\bar{v}:=\frac{1}{N}\sum_{j=1}^{N}v_{j}(0)$. The shifted
dynamics is given by

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}r_{i}$
$\displaystyle=v_{i}-\bar{v},$ $\displaystyle r_{i}(0)=(r_{0})_{i},$ (6a)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}v_{i}$
$\displaystyle=\frac{1}{N}\sum_{j=1}^{N}\psi(|r_{i}-r_{j}|)(v_{j}-v_{i})-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}(r_{i}-r_{j}),$
$\displaystyle v_{i}(0)=(v_{0})_{i}.$ (6b)

We want to emphasize that system (1) contains explicit position information of
the particles, which we loose in the PHS formulation as we shift by the center
of mass. However, pattern formation and uniform stability which are of
interest in the context of interacting particle systems are translational
invariant, hence the explicit position information plays a minor role.

The mean velocity $\bar{v}$ can be incorporated in the Hamiltonian leading to

$\mathcal{H}^{N}(z)=\frac{1}{2}\sum_{i=1}^{N}\Big{(}(v_{i}-\bar{v})^{\top}(v_{i}-\bar{v})+\frac{1}{N}\sum_{j=1}^{N}\mathcal{V}(r_{i}-r_{j})\Big{)}$

with partial derivative

$\displaystyle\frac{\partial\mathcal{H}^{N}(z)}{\partial
z}=\left(-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}(x_{j}-x_{1}),\ldots,-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}(x_{j}-x_{N}),v_{1}-\bar{v},\dots,v_{N}-\bar{v}\right)^{\top}.$

For notational convenience, we define
$v-\textbf{1}\bar{v}:=(v_{1}-\bar{v},\dots,v_{N}-\bar{v})^{\top}$.

To find the port-Hamiltonian structure we note that the upper-left part of the
skew-symmetric matrix is predefined by the differential equation for $r$.
Indeed, we obtain

$\frac{\mathrm{d}}{\mathrm{d}t}r=I(v-\textbf{1}\bar{v}).$

Thanks to the structure of $\Psi$, see Remark Remark, it holds

$\Psi(z)\textbf{1}\bar{v}=0,$ (7)

and hence we obtain

$\frac{\mathrm{d}}{\mathrm{d}t}v=-I\frac{\partial\mathcal{H}^{N}(z)}{\partial
r}-\Psi(z)\frac{\partial\mathcal{H}^{N}(z)}{\partial v}$

with

$\Psi(z)=\frac{1}{N}\begin{pmatrix}\displaystyle\sum_{j=1}^{N}\\!\psi(|r_{j}-r_{1}|)&-\psi(|r_{2}-r_{1}|)&\ldots&-\psi(|r_{N}-r_{1}|)\\\
-\psi(|r_{1}-r_{2}|)&\displaystyle\sum_{j=1}^{N}\\!\psi(|r_{j}-r_{2}|)&\ddots&\vdots\\\
\vdots&\ddots&\ddots&&\\\ &&&-\psi(|r_{N}\\!-\\!r_{N-1}|)\\\
-\psi(|r_{1}-r_{N}|)&\ldots&-\psi(|r_{N-1}\\!-\\!r_{N}|)&\displaystyle\sum_{j=1}^{N}\\!\psi(|r_{j}-r_{N}|)\end{pmatrix}.$

Combining the two equations yields the port-Hamiltonian structure

$\frac{\mathrm{d}}{\mathrm{d}t}z=\left[\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}-\begin{pmatrix}0&0\\\
0&\Psi(z)\end{pmatrix}\right]\frac{\partial\mathcal{H}^{N}(z)}{\partial z}.$
(8)

The state space of the port-Hamiltonian system (8) is given by
$X:=\mathbb{R}^{Nd+Nd}$.

The main results of the article are based on this formulation. The following
theorem discusses the well-posedness and some characteristics of solutions to
the port-Hamiltonian system.

######  Theorem 3.1.

Let Assumption 2.1 hold, then for every initial condition
$z_{0}=(r_{0},v_{0})\in\mathbb{R}^{Nd+Nd}$ the port-Hamiltonian system (8)
possesses a unique global solution $z=(r,v)$ satisfying the dissipativity
inequality

$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{H}^{N}(z)=-\langle(v-\textbf{1}\bar{v}),\Psi(z)(v-\textbf{1}\bar{v})\rangle\leq
0,$ (9)

where

$\bar{v}:=\sum_{j=1}^{N}v_{j}(0).$

Further, the velocity of the center of mass is conserved, that is,
$\bar{v}=\sum_{j=1}^{N}v_{j}(t)$, $t\geq 0$.

###### Proof.

The existence of a unique global solution $z=(r,v)$ follows from Proposition
2.2. The non-negativity of $\psi$ yields $\Psi(z)$ is real, symmetric and
diagonally dominant, hence positive semi-definite. Thus we compute for the
solution $z$

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{H}^{N}(z)$
$\displaystyle=\left\langle\frac{\partial\mathcal{H}^{N}(z)}{\partial
z},\frac{\mathrm{d}}{\mathrm{d}t}z\right\rangle=\left\langle\frac{\partial\mathcal{H}^{N}(z)}{\partial
z},\left[\begin{pmatrix}0&I\\\ -I^{\top}&0\end{pmatrix}-\begin{pmatrix}0&0\\\
0&\Psi(z)\end{pmatrix}\right]\frac{\partial\mathcal{H}^{N}(z)}{\partial
z}\right\rangle$ $\displaystyle=-\langle
v-\textbf{1}\bar{v},\Psi(z)v-\textbf{1}\bar{v}\rangle\leq 0.$

This shows the dissipativity inequality. Finally, let
$z_{0}=(r_{0},v_{0})\in\mathbb{R}^{Nd}\times\mathbb{R}^{Nd}$. Then (5) implies
that $\sum_{i=1}^{N}v_{i}$ is constant. This completes the proof. ∎

In the following we assume that Assumption 2.1 is satisfied, and thus the
port-Hamiltonian system (8) has for every initial condition a unique global
solution. Next we investigate conserved quantities of our port-Hamiltonian
formulation of interacting particle systems.

######  Definition 3.2 (Definition 6.4.1 in [38]).

A function $C:X\rightarrow\mathbb{R}$ that is partially differentiable is
called a _Casimir function_ for the port-Hamiltonian system (8) if

$\left(\frac{\partial C}{\partial
z}(z)\right)^{\top}\left[\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}-\begin{pmatrix}0&0\\\ 0&\Psi(z)\end{pmatrix}\right]=0,\qquad
z\in X.$

A Casimir function is a conserved quantity as for solutions $z$ we obtain

$\frac{\mathrm{d}}{\mathrm{d}t}C(z)=\left(\frac{\partial C}{\partial
z}(z)\right)^{\top}\left[\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}-\begin{pmatrix}0&0\\\
0&\Psi(z)\end{pmatrix}\right]\frac{\partial\mathcal{H}^{N}(z)}{\partial z}=0$

independently of the Hamiltonian $\mathcal{H}^{N}$.

######  Proposition 3.3.

A function $C:X\rightarrow\mathbb{R}$ that is partially differentiable is a
Casimir function for the port-Hamiltonian system (8) if and only if
$C(r,v)=\gamma$ for some $\gamma\in\mathbb{R}$.

###### Proof.

In [38] it is shown that a function $C:X\rightarrow\mathbb{R}$ that is
partially differentiable is a Casimir function for the port-Hamiltonian system
(8) if and only if

$\left(\frac{\partial C}{\partial z}(z)\right)^{\top}\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}=0\quad\text{and}\quad\left(\frac{\partial C}{\partial
z}(z)\right)^{\top}\begin{pmatrix}0&0\\\ 0&\Psi(z)\end{pmatrix}=0,\qquad z\in
X.$

This is equivalent to

$\frac{\partial C}{\partial r}(z)=0\quad\text{and}\quad\frac{\partial
C}{\partial v}(z)=0,\qquad z\in X,$

and thus the statement of the proposition follows. ∎

######  Remark.

Note that the system matrices of the formulations in $(x,v)$ and $(r,v)$
coincide. This allows to conclude that the the Casimir functions of the
different formulations coincide as well.

######  Remark.

We want to emphasize that, in contrast to [27], we do not require a null space
condition in order to define the Casimir function in neither of the two
formulations. This is due to the different choice of the port-Hamiltonian
formulation.

We conclude this section with a sufficient condition for asymptotic stability.
Here the function $\mathcal{H}^{N}$ serves as candidate for a suitable
Lyapunov function. The following lemma will be useful.

######  Lemma 3.4.

Let Assumption 2.1 hold, $\psi(x)>0$ for $x\in\mathbb{R}$ and we define

$\mathcal{N}:=\\{z=(r,v)\in X\mid-\langle
v-\mathbf{1}\bar{v},\Psi(z)(v-\mathbf{1}\bar{v})\rangle=0\\}.$

Let $z:[0,\infty)\rightarrow X$, $z=(r,v)$, be a solution of (6) with
$z(t)\in\mathcal{N}$ for $t\geq 0$. Then the function $v$ is constant,

$v_{i}(t)=\frac{1}{N}\sum_{j=1}^{N}v_{j}(0)\quad\text{and}\quad\sum_{j=1}^{N}\nabla\mathcal{V}(r_{i}(t)-r_{j}(t))=0,\qquad
i=1,\ldots,N,\,t\geq 0.$

###### Proof.

As $\psi(x)>0$ for $x\in\mathbb{R}$, we obtain $v_{i}(t)=v_{j}(t)$ for
$i,j=1,\ldots,N$ and $t\geq 0$. This together with Theorem 3.1 implies that
$v$ is constant and satisfies $v_{i}(t)=\frac{1}{N}\sum_{j=1}^{N}v_{j}(0)$ for
$i=1,\ldots,N$, and $t\geq 0$. The remaining statement follows directly from
equation (6). ∎

######  Theorem 3.5.

Let Assumption 2.1 hold, $\psi(x)>0$ for $x\in\mathbb{R}$ and
$\nabla\mathcal{V}$ bounded. Then for every initial condition
$z_{0}=(r_{0},v_{0})\in\mathbb{R}^{Nd}\times\mathbb{R}^{Nd}$ the corresponding
solution $z$ of (8) satisfies

$\lim_{t\rightarrow\infty}\operatorname{dist}(z(t),L)=0,$

where

$L:=\Big{\\{}(\tilde{r},\tilde{v})\in\mathbb{R}^{Nd}\times\mathbb{R}^{Nd}\mid\tilde{v}_{i}=\frac{1}{N}\sum_{j=1}^{N}v_{j}(0),\,\,\sum_{j=1}^{N}\nabla\mathcal{V}(\tilde{r}_{i}-\tilde{r}_{j})=0,\,i=1,\ldots,N\Big{\\}}.$

###### Proof.

Our goal is to apply LaSalle’s stability theorem [26, Theorem 3.2.11], hence
we have to show that the trajectories are contained in a compact subset of
$\mathbb{R}^{Nd}\times\mathbb{R}^{Nd}$.

Let $\bar{v}$ the mean velocity as above. We estimate

$\displaystyle\|v(t)-\textbf{1}\bar{v}\|^{2}$
$\displaystyle=\|v(0)-\textbf{1}\bar{v}\|^{2}+\int_{0}^{t}2(v-\textbf{1}\bar{v})\cdot\frac{d}{ds}(v-\textbf{1}\bar{v})\mathrm{d}s$
$\displaystyle=\|v(0)-\textbf{1}\bar{v}\|^{2}+\int_{0}^{t}2(v-\textbf{1}\bar{v})(-\frac{\partial}{\partial
r}\mathcal{H}^{N}(z)-\Psi(z)(v-\textbf{1}\bar{v}))\mathrm{d}s,$

where we used $\frac{d}{dt}\bar{v}=0$ and $\Psi(z)\textbf{1}\bar{v}=0$. The
boundedness of $\nabla\mathcal{V}$ and the Peter-Paul inequality allow us to
estimate for $\varepsilon>0$

$\displaystyle\|v(t)-\textbf{1}\bar{v}\|^{2}$
$\displaystyle\leq\|v(0)-\textbf{1}\bar{v}\|^{2}+\int_{0}^{t}\frac{1}{\varepsilon}\|\nabla\mathcal{V}\|_{\infty}^{2}+\varepsilon\|v-\textbf{1}\bar{v}\|^{2}-2(v-\textbf{1}\bar{v})\Psi(z)(v-\textbf{1}\bar{v})\mathrm{d}s$
$\displaystyle\leq\|v(0)-\textbf{1}\bar{v}\|^{2}+\frac{t}{\varepsilon}\|\nabla\mathcal{V}\|_{\infty}^{2}-(2\lambda_{2}-\varepsilon)\int_{0}^{t}\|v-\textbf{1}\bar{v}\|^{2}\mathrm{d}s,$

where we denote by $\lambda_{2}$ the second smallest eigenvalue of $\Psi(z)$.
Note that $\lambda_{2}>0$ since the $(n-1)\times(n-1)$ submatrix
$\big{(}\Psi(z)_{ij}\big{)}_{i=1,\dots,n-1;j=1,\dots,n-1}$ is strictly
diagonally dominant and therefore positive definite.

We define

$\alpha(t):=\|v(0)-\textbf{1}\bar{v}\|^{2}+\frac{t}{\varepsilon}\|\nabla\mathcal{V}\|_{\infty}^{2},$

and then Gronwall’s inequality yields

$\|v(t)-\textbf{1}\bar{v}\|^{2}\leq\alpha(t)e^{-(2\lambda_{2}-\varepsilon)t}.$

Choosing $\varepsilon$ such that $2\lambda_{2}>\varepsilon>0$, we find that
$v(t)$ relaxed towards $\textbf{1}\bar{v}$. In particular,

$\sup\\{\|v(t)-\textbf{1}\bar{v}\|^{2}\mid t\geq 0\\}<\infty,$

which shows the boundedness of the solution trajectory w.r.t. the velocity
variable.

We are left to show the boundedness of the trajectory w.r.t. the positions. We
estimate

$\displaystyle\|r(t)\|\leq\|r(0)\|+\int_{0}^{t}\|v(s)-\textbf{1}\bar{v}\|\mathrm{d}s\leq\|r(0)\|+\int_{0}^{\infty}\alpha(s)e^{-(2\lambda_{2}-\varepsilon)s}\mathrm{d}s.$

As the integral on the right hand side converges, we conclude that also the
position variables of the trajectories are bounded uniformly for all $t\geq
0$. Altogether, for any initial data the trajectories are contained in a
compact set for all times. The statement follows by LaSalle’s stability
theorem. ∎

######  Example 3.6.

Let us consider the example of the Cucker-Smale dynamics with potential, where
the alignment function [14]

$\psi(|r_{j}-r_{i}|)=\frac{K}{(\delta^{2}+|r_{j}-r_{i}|^{2})^{\beta}},\qquad
K,\delta>0\text{ and }\beta\geq 0$

and (regularized) Morse potential [20]

$\mathcal{V}(d)=Re^{-|d|^{2}/r}-Ae^{-|d|^{2}/a},\qquad R,A\geq 0\text{ and
}r,a>0.$

are explicitly given. Note that $\psi\colon\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}_{>0}$ is strictly positive and suppose that $r>a$ and
$A>0$. Then $\mathcal{H}^{N}$ satisfies the assumption of Theorem 3.5. Thus
the Hamiltonian decreases as long as the velocities of the swarm members are
not aligned. In particular, this yields unconditional flocking. Comparing our
result to [13] we find that the PHS structure allows us to boil the argument
for flocking down to an application of LaSalle’s stability theorem. However,
we note that we require $\mathcal{V}\neq 0$ in the proof. Instead the argument
in [13] exploits structure of the support of the solution to the particle
system and holds for $\mathcal{V}\equiv 0$.

######  Remark.

Also the Kuramoto model with inertia and fully connected incidence matrix [35]
fits into the setting with $\mathcal{V}(r_{i}-r_{j})=-\cos(r_{i}-r_{j})$
leading to

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}r_{i}$
$\displaystyle=v_{i}-\bar{v},\qquad i=1,\dots,N$ (10a)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}v_{i}$ $\displaystyle=-\gamma
v_{i}+\frac{1}{N}\sum_{j=1}^{N}\sin(r_{i}-r_{j})$ (10b) $\displaystyle
r_{i}(0)$ $\displaystyle=\hat{r}_{i},\qquad v_{i}(0)=\hat{v}_{i}.$ (10c)

with friction parameter $\gamma>0$, can be formulated as port-Hamiltonian
system. We obtain the Hamiltonian

$\mathcal{H}^{N}(z)=\frac{1}{2}\sum_{i=1}^{N}\Big{(}(v_{i}-\bar{v})^{\top}(v_{i}-\bar{v})-\frac{1}{N}\sum_{j=1}^{N}\cos(r_{i}-r_{j})\Big{)},\quad
z=(r_{1},\dots,r_{N},v_{1},\dots,v_{N})$

and the dynamics

$\frac{\mathrm{d}}{\mathrm{d}t}z=(J-R)\frac{\partial\mathcal{H}^{N}(z)}{\partial
z},\;z(0)=z_{0},\quad J=\begin{pmatrix}0&I\\\
-I&0\end{pmatrix},\;R=\begin{pmatrix}0&0\\\ 0&\gamma I\end{pmatrix}.$

As $R$ has higher rank compared to the one discussed above, the stability
result with equilibrium point $v_{i}=0$ is easier to show for this particular
case of the Kuramoto model.

##  4 Mean-field limit

We obtain a candidate for the mean-field equation of the shifted dynamics
following the standard derivation and moreover the mean-field Hamiltonian by
rescaling with $\frac{1}{N}$. In fact, using the empirical measure

$f^{N}(t,r,v)=\frac{1}{N}\sum_{i=1}^{N}\delta(r-r_{i}(t))\otimes\delta(v-v_{i}(t)),$

were $(r_{i},v_{i})_{i=1}^{N}$ denotes the solution of (6), we obtain the PDE
describing the mean-field dynamics given by

$\partial_{t}f^{N}+\nabla_{r}\cdot\Big{(}\big{(}v-\bar{v}\big{)}f^{N}\Big{)}=\nabla_{v}\cdot\Big{(}\big{(}\int\psi(r-\hat{r})(v-\hat{v})+\nabla\mathcal{V}(r-\hat{r})\mathrm{d}f^{N}(\hat{r},\hat{v})\big{)}f^{N}\Big{)},$
(11)

where $\bar{v}=\int\hat{v}\mathrm{d}f^{N}(0,\hat{r},\hat{v})$. Moreover, for
the mean-field Hamiltonian we obtain

$\lim\limits_{N\rightarrow\infty}\frac{1}{N}\mathcal{H}^{N}(z)=\lim\limits_{N\rightarrow\infty}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\frac{1}{2}\left((v-\bar{v})^{\top}(v-\bar{v})+\mathcal{V}(r-\hat{r})\right)\mathrm{d}f^{N}(t,\hat{r},\hat{v})\mathrm{d}f^{N}(t,r,v).$

This motivates to define the Hamiltonian for the mean-field equation as

$\displaystyle\mathcal{H}(f)$
$\displaystyle=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\frac{1}{2}\left((v-\bar{v})^{\top}(v-\bar{v})+\mathcal{V}(r-\hat{r})\right)\mathrm{d}f(t,\hat{r},\hat{v})\mathrm{d}f(t,r,v)$
$\displaystyle=\frac{1}{2}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}f\ast\big{(}v^{\top}v+\mathcal{V}(r)\big{)}\mathrm{d}f(t,r,v),$

where $\bar{v}=\int\hat{v}\mathrm{d}f(t,\hat{r},\hat{v})$.

To check if the Hamiltonian structure is preserved in the mean-field limit we
compute the variation of $\mathcal{H}(f)$ in the space of probability
measures, see also [6]. In order to preserve the normalization of the
measures, we consider the push forward of $f$ w.r.t.
$(\text{id}+\varepsilon\zeta)$ for
$\zeta\in\mathcal{C}_{c}(\mathbb{R}^{2d},\mathbb{R}^{2d})$,
$\zeta(r,v)=(\zeta_{r}(r,v),\zeta_{v}(r,v))$ and $\varepsilon>0$ to find

$\displaystyle\mathcal{H}((\text{id}+\varepsilon\zeta)_{\\#}f)-\mathcal{H}(f)$
$\displaystyle=\frac{1}{2}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\left(v+\varepsilon\zeta_{v}(r,v)-\bar{v}\right)^{\top}\left(v+\varepsilon\zeta_{v}(r,v)-\bar{v}\right)\mathrm{d}f(t,\hat{r},\hat{v})\mathrm{d}f(t,r,v)$
$\displaystyle\quad+\frac{1}{2}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\mathcal{V}((r+\varepsilon\zeta_{r}(r,v))-\hat{r}-\varepsilon\zeta_{r}(\hat{r},\hat{v})))\mathrm{d}f(t,\hat{r},\hat{v})\mathrm{d}f(t,r,v)$
$\displaystyle\quad-\frac{1}{2}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\left((v-\bar{v})^{\top}(v-\bar{v})+\mathcal{V}(r-\hat{r})\right)\mathrm{d}f(t,\hat{r},\hat{v})\mathrm{d}f(t,r,v)$
$\displaystyle=\varepsilon\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}(v-\bar{v})^{\top}\zeta_{v}(r,v)+\nabla\mathcal{V}(r-\hat{r})\zeta_{r}(r,v)\,\mathrm{d}f(t,\hat{r},\hat{v})\mathrm{d}f(t,r,v)+o(\varepsilon),$

where we used (2) and Fubini’s Theorem. In the limit $\varepsilon\rightarrow
0$ we obtain

$\lim_{\varepsilon\rightarrow
0}\frac{1}{\varepsilon}\Big{(}\mathcal{H}((\text{id}+\varepsilon\zeta)_{\\#}f)-\mathcal{H}(f)\Big{)}=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\begin{pmatrix}\nabla\mathcal{V}(r-\hat{r})\\\
v-\bar{v}\end{pmatrix}\cdot\zeta(r,v)\mathrm{d}f(t,\hat{r},\hat{v})\mathrm{d}f(t,r,v).$

Following [4] we can identify

$\nabla_{f}\mathcal{H}(f)=\begin{pmatrix}\left(\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\nabla\mathcal{V}(r-\hat{r})\,\mathrm{d}f(t,\hat{r},\hat{v})\right)\\\
v-\bar{v}\end{pmatrix}=\begin{pmatrix}\left(\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\nabla\mathcal{V}(r-\hat{r})\,\mathrm{d}f(t,\hat{r},\hat{v})\right)\\\
\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}(v-\hat{v})\;\mathrm{d}f(t,\hat{r},\hat{v})\end{pmatrix}=f\ast\begin{pmatrix}\nabla\mathcal{V}(r)\\\
v\end{pmatrix}.$ (12)

We can rewrite the mean-field equation as

$\partial_{t}f=-\nabla_{(r,v)}\cdot\left(f\left(f\ast\left(\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}-\begin{pmatrix}0&0\\\
0&\psi(r)\end{pmatrix}\right)\begin{pmatrix}\nabla\mathcal{V}(r)\\\
v\end{pmatrix}\right)\right).$ (13)

As expected we obtain the dissipativity inequality also on the mean-field
level.

######  Theorem 4.1.

Let $\psi(x)>0$ for all $x\in\mathbb{R}$ and Assumption 2.1 hold. Then the
dynamics (13) admits for every initial condition
$f_{0}\in\mathcal{P}_{c}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ a unique
solution $f$ in the sense of Definition 2.3. Moreover, the solution is
dissipative, i.e. it holds

$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{H}(f)\leq 0.$

Further, for every initial condition
$f_{0}\in\mathcal{P}_{c}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ there exists a
vector $\bar{v}\in\mathbb{R}^{d}$ such that

$\int v\,\mathrm{d}f(t,r,v)=\bar{v}\quad\text{for }t\geq 0.$

###### Proof.

We prove the last statement first. Indeed, using integration by parts we
calculate

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\int v\,\mathrm{d}f(t,r,v)$
$\displaystyle=\int
v\,\nabla_{v}\cdot\Big{(}\int\psi(|r-\hat{r}|)(v-\hat{v})-\nabla\mathcal{V}(r-\hat{r})\mathrm{d}f(t,\hat{r},\hat{v})f\Big{)}\mathrm{d}r\mathrm{d}v$
$\displaystyle=\iint\psi(|r-\hat{r}|)(v-\hat{v})-\nabla\mathcal{V}(r-\hat{r})\mathrm{d}f(t,\hat{r},\hat{v})\mathrm{d}f(t,r,v)$
$\displaystyle=-\frac{1}{2}\iint\nabla\mathcal{V}(r-\hat{r})\mathrm{d}f(t,\hat{r},\hat{v})\mathrm{d}f(t,r,v)+\frac{1}{2}\iint\nabla\mathcal{V}(r-\hat{r})\mathrm{d}f(t,\hat{r},\hat{v})\mathrm{d}f(t,r,v)$
$\displaystyle=0,$

which proves that the first moment with respect to the velocity is preserved.
The well-posedness of (13) is obtained with the same arguments as (3) in
Proposition 2.4.

For the dissipativity we use the product rule, the symmetry of $\psi$ and the
anti-symmetry of $\nabla\mathcal{V}$ to obtain

$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{H}(f)=-\frac{1}{2}\iint(v-\hat{v})\cdot(\psi(|r-\hat{r}|)(v-\hat{v}))\mathrm{d}f(t,\hat{r},\hat{v})\mathrm{d}f(t,r,v)\leq
0.\qed$

Next we investigate the characteristics and their port-Hamiltonian
formulation. The characteristics $Z_{t}(z)=(R_{t}(r),V_{t}(v))$ read

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}R_{t}(r)$
$\displaystyle=V_{t}(v)-\int
V_{t}(v)\mathrm{d}f_{0}(r,v),\qquad\text{law}(z)=f(0,r,v),$ (14a)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}V_{t}(v)$
$\displaystyle=\int\psi(|R_{t}(r)-R_{t}(\hat{r})|)(V_{t}(\hat{v})-V_{t}(v))-\nabla\mathcal{V}(R_{t}(r)-R_{t}(\hat{r}))\,\mathrm{d}f(0,\hat{r},\hat{v}),$
(14b)

and we can rewrite a solution to the mean-field equation as
$f(t)=Z_{t}{{}_{\\#}}f(0),$ where $\\#$ denotes the push-forward operator
given by

$\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}f(z)\,\mathrm{d}f(t,z)=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}f(Z_{t}(z))\,\mathrm{d}f(0,z).$

A natural question is if the Dobrushin inequality in combination with the
stability result for the ODE dynamics yields a stability result on the mean-
field level. This would require the interchangebility of the limits
$N\rightarrow\infty$ and $t\rightarrow\infty$ which is beyond the scope of
this article. However, we prove the stability on the mean-field limit with the
help of LaSalle’s theorem in metric spaces [26, Theorem 3.2.11] in the
following. Clearly, the precompactness argument is more involved in this
setting. We begin with a the mean-field analogue of Lemma 3.4. For notational
convenience we write in the following $f_{t}$ for $f(t,\cdot)$ as we already
did in Section 2.

######  Lemma 4.2.

Let Assumption 2.1 hold, $\psi(x)>0$ for $x\in\mathbb{R}$ and

${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathcal{N}:=\\{h\in\mathcal{P}_{2}(\mathbb{R}^{d}\times\mathbb{R}^{d})\mid\iint(v-\hat{v})^{\top}\psi(|r-\hat{r}|)(v-\hat{v})\mathrm{d}h(\hat{r},\hat{v})\mathrm{d}h(r,v)=0\\}.}$

Let
$f\colon[0,\infty)\rightarrow\mathcal{P}_{2}(\mathbb{R}^{d}\times\mathbb{R}^{d})$
be a solution of (13) with $f_{t}\in\mathcal{N}$ for $t\geq 0$. Then the
support of $f_{t}$ in the velocity space is concentrated at
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bar{v}:=\int
v\mathrm{d}f_{0}(r,v)}$, i.e. there exists
$g_{t}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ such that it holds

$f_{t}(r,v)=g_{t}(r)\delta(v-{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bar{v}})\quad\text{and}\quad\int\nabla\mathcal{V}(r-\hat{r})\mathrm{d}g_{t}(\hat{r})=0,\quad\text{for
all }r\in\mathrm{supp}(g_{t}),t\geq 0.$

###### Proof.

As $\psi(x)>0$ for $x\in\mathbb{R}$, the support of $f_{t}$ in the velocity
variable has to be concentrated at some velocity
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bar{v}}$
for all $t\geq 0$ in order to satisfy the integration condition in the
definition of $\mathcal{N}$. Together with Theorem 4.1 this implies that
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bar{v}=\int
v\mathrm{d}f_{0}(r,v)}$. The remaining statement follows directly from
equation (13). ∎

######  Theorem 4.3.

Let Assumption 2.1 hold and there exists $\underline{\psi}$ with
$\psi(x)\geq\underline{\psi}>0$ for $x\in\mathbb{R}$. Assume that
$\mathcal{H}$ is bounded from below and that $\nabla\mathcal{V}$ is bounded.
Then for every initial condition
$f_{0}\in\mathcal{P}_{c}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ the
corresponding solution $f$ of (13) satisfies

$\lim\limits_{t\rightarrow\infty}\mathrm{dist}(f_{t},\mathcal{L})=0,$

where

$\mathcal{L}=\left\\{f\in\mathcal{P}_{2}(\mathbb{R}^{d}\times\mathbb{R}^{d})\mid
f=g(x)\delta(v-\hat{v}),\quad\int\nabla\mathcal{V}(x-\bar{x})\mathrm{d}g(x)=0,\quad\hat{v}=\int
v\,\mathrm{d}f_{0}(x,v)\right\\}.$

###### Proof.

The proof follows the lines of the finite-dimensional result. The aim is again
to employ LaSalle’s invariance theorem. Hence we have to show that the
probability measure corresponding to the solution stays compactly supported
for all times $t\geq 0$, since then the precompactness w.r.t.
$\mathcal{W}_{2}$ follows by [30, Proposition 2.2.3].

We first consider the 2-Wasserstein distance of $f_{t}$ and
$\delta_{\bar{v}}$, where $\bar{v}=\int v\mathrm{d}f_{t}(r,v)$, given by

$\displaystyle\mathcal{W}_{2}$ $\displaystyle(f_{t},\delta_{\bar{v}})^{2}$
$\displaystyle=\int|v-\bar{v}|^{2}\mathrm{d}f_{t}=\int|v-\bar{v}|^{2}\mathrm{d}f_{0}+\int_{0}^{t}\int|v-\bar{v}|^{2}\partial_{s}f_{s}\mathrm{d}z\mathrm{d}s$
$\displaystyle=\int|v-\bar{v}|^{2}\mathrm{d}f_{0}-\int_{0}^{t}\int|v-\bar{v}|^{2}\cdot\nabla_{(r,v)}\cdot\left(f\ast\left(\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}-\begin{pmatrix}0&0\\\
0&\psi(r)\end{pmatrix}\right)\begin{pmatrix}\nabla\mathcal{V}(r)\\\
v\end{pmatrix}\right)\mathrm{d}f_{s}\mathrm{d}s$
$\displaystyle=\int|v-\bar{v}|^{2}\mathrm{d}f_{0}-2\int_{0}^{t}\int(f\ast
v)(f\ast(\psi(r)v+f\ast\nabla\mathcal{V}(r)))\mathrm{d}f_{s}\mathrm{d}s.$

A simple computation shows that $(f\ast
v)(f\ast\psi(r)v)=\frac{1}{2}f\ast(v\psi(r)v)$ which will be helpful in the
following.

This and Peter and Paul inequality applied to the term with the potential
interactions allows us to further estimate

$\displaystyle\int|v-\bar{v}|^{2}\mathrm{d}f_{t}$
$\displaystyle\leq\int|v-\bar{v}|^{2}\mathrm{d}f_{0}-\int_{0}^{t}\int\int(\underline{\psi}-\varepsilon)|v-\hat{v}|^{2}\mathrm{d}\hat{f}_{s}\mathrm{d}f_{s}\mathrm{d}s+\frac{\|\nabla\mathcal{V}\|_{\infty}^{2}}{\varepsilon}$
$\displaystyle\leq\int|v-\bar{v}|^{2}\mathrm{d}f_{0}-\int_{0}^{t}\int(\underline{\psi}-\varepsilon)|v-\bar{v}|^{2}\mathrm{d}f_{s}\mathrm{d}s+\frac{\|\nabla\mathcal{V}\|_{\infty}^{2}}{\varepsilon}.$

An application of Gronwall’s inequality yields

$\mathcal{W}_{2}(f_{t},\delta_{\bar{v}})^{2}=\int|v-\bar{v}|^{2}\mathrm{d}f_{t}\leq\left(\int|v-\bar{v}|^{2}\mathrm{d}f_{0}+\frac{\|\nabla\mathcal{V}\|_{\infty}^{2}}{\varepsilon}\right)e^{-(\underline{\psi}-\varepsilon)t}.$

This shows that the support of $f_{t}$ w.r.t. $v$ strictly decays over time.

We use this result to show the boundedness of the support w.r.t. $r$. Indeed,
we obtain

$\displaystyle\int|r|^{2}\mathrm{d}f_{t}$
$\displaystyle=\int|r|^{2}\mathrm{d}f_{0}-2\int_{0}^{t}\int
r\cdot(v-\bar{v})\mathrm{d}f_{s}\mathrm{d}s$ $\displaystyle\leq
1+\int|r|^{2}\mathrm{d}f_{0}+\int_{0}^{t}\left(\int|r|^{2}\mathrm{d}f_{s}\right)^{1/2}\;\left(\int|v-\bar{v}|^{2}\mathrm{d}f_{s}\right)^{1/2}\mathrm{d}s,$

and thus

$\displaystyle\left(\int|r|^{2}\mathrm{d}f_{t}\right)^{1/2}$
$\displaystyle\leq
1+\int|r|^{2}\mathrm{d}f_{0}+\int_{0}^{t}\left(\int|r|^{2}\mathrm{d}f_{s}\right)^{1/2}\;\left(\int|v-\bar{v}|^{2}\mathrm{d}f_{s}\right)^{1/2}\mathrm{d}s,$

where we used $x>\sqrt{x}$ for $x>1$.

We define $\alpha:=(1+\int|r|^{2}\mathrm{d}f_{0})$ and
$\beta:=\left(\int|v-\bar{v}|^{2}\mathrm{d}f_{0}+\frac{\|\nabla\mathcal{V}\|_{\infty}^{2}}{\varepsilon}\right)$
Then we obtain with Gronwall inequality

$\left(\int|r|^{2}\mathrm{d}f_{t}\right)^{1/2}\leq\alpha\exp(\int_{0}^{t}\beta
e^{-(\underline{\psi}-\varepsilon)s}\mathrm{d}s).$

Since the integral converges, also the support of $f_{t}$ w.r.t. $r$ is
compactly supported for all $t\geq 0$.

As the support of $f_{t}$ is compactly supported for all times, we are allowed
to apply LaSalle’s invariance principle to obtain the result. ∎

##  5 PHS preserving coupling of subsystems

To discuss strategies for the coupling of subsystems, we begin with the
identification of the ports of the generalized mass, spring and damper
components which model the interacting particle system.

###  5.1 Identification of ports

The interaction dynamics can be interpreted as a generalized mass-spring-
damper system. In order to identify this we decouple the system into its
smallest parts. This allows us then to discuss PHS preserving coupling of
different subsystems.

The $i$-th mass is described by its position $x_{i}$ and its velocity $v_{i}$.
Its evolution is driven by its kinetic energy
$H^{m}(x_{i},v_{i})=v_{i}^{\top}v_{i}$ and deviations from this velocity are
due to external forces $f_{i}$ which are called flows in the PHS framework.
Together with the effort $e_{i}=\frac{\partial}{\partial
v_{i}}H^{m}(x_{i},v_{i})=v_{i}$ this is leading to the dynamics

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}x_{i}=v_{i},\quad\frac{\mathrm{d}}{\mathrm{d}t}v_{i}=f_{i},\quad\Rightarrow\quad\frac{\mathrm{d}}{\mathrm{d}t}z_{i}=J\nabla
H^{m}(z_{i})+B_{i}u_{i},\quad y_{i}=B_{i}^{\top}\nabla H^{m}(z)=v_{i},$

where $z_{i}=(x_{i},v_{i})$, $J=\left(\begin{smallmatrix}0&I\\\
-I&0\end{smallmatrix}\right)$, $B_{i}=\left(\begin{smallmatrix}0\\\
I\end{smallmatrix}\right)$ and $u_{i}=f_{i}$ is the input and
$y_{i}=v_{i}=e_{i}$ the output.

The spring and damper connecting mass $i$ and mass $j$ are described by the
relative position $q_{ij}:=x_{i}-x_{j}$ and the relative velocity
$v_{ij}:=v_{i}-v_{j}$. Note that the damper is a purely dissipative element,
hence it admits no Hamiltonian but the force
$f_{ij}^{\mathrm{damper}}=-\psi(|q_{ij}|)v_{ij}$. The Hamiltonian of the
spring is given by $H^{sd}(q_{ij},v_{ij})=\mathcal{V}(q_{ij})$. Altogether
this leads to the dynamics

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}q_{ij}=v_{ij},\qquad\frac{\mathrm{d}}{\mathrm{d}t}v_{ij}=-\frac{\partial}{\partial
q_{ij}}H^{sd}(q_{ij},v_{ij})-\psi(|q_{ij}|)v_{ij}$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}z_{ij}=(J-R(z_{ij}))\nabla
H^{sd}(z_{ij})+B_{ij}u_{ij},$

where $z_{ij}=(q_{ij},v_{ij})$, $J=\left(\begin{smallmatrix}0&I\\\
-I&0\end{smallmatrix}\right)$, $R(z_{ij})=\left(\begin{smallmatrix}0&0\\\
0&\psi(|q_{ij}|)\end{smallmatrix}\right),B_{ij}=\left(\begin{smallmatrix}I&-I\\\
0&0\end{smallmatrix}\right)$ and $u_{ij}=\left(\begin{smallmatrix}v_{i}\\\
v_{j}\end{smallmatrix}\right)$. The corresponding output is given by

$y_{ij}=B_{ij}^{\top}\nabla H^{sd}(z_{ij})=\begin{pmatrix}I&0\\\
-I&0\end{pmatrix}\nabla H^{sd}(z_{ij})=\begin{pmatrix}I&0\\\
-I&0\end{pmatrix}\begin{pmatrix}\nabla\mathcal{V}(q_{ij})\\\
0\end{pmatrix}=\begin{pmatrix}\;\;\nabla\mathcal{V}(q_{ij})\\\
-\nabla\mathcal{V}(q_{ij})\end{pmatrix}.$

The flow and effort variable of the spring-damper system connecting mass $i$
and mass $j$ are given by

$f_{ij}=v_{i}-v_{j}\quad\textrm{ and }\quad e_{ij}=\frac{\partial}{\partial
q_{ij}}H^{sd}(q_{ij},v_{ij})+\psi(|q_{ij}|)v_{ij}=\nabla\mathcal{V}(q_{ij})+\psi(|q_{ij}|)v_{ij}.$

To couple the $i$-th mass with spring-damper system connecting mass $i$ and
mass $j$ we define $f_{i}=-e_{ij}$ and $e_{i}=f_{ij}$ leading to

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}x_{i}=v_{i},\quad\frac{\mathrm{d}}{\mathrm{d}t}v_{i}=-\nabla\mathcal{V}(q_{ij})-\psi(|q_{ij}|)v_{ij}.$

Taking all binary interactions into account and summation over all binary
interactions yields (1), which leads to the PHS formulation in closed form as
given in (4). Alternatively, instead of considering the actual positions
$x_{i}$ of the particles, we can use the relative positions $q_{ij}$ which
then yields the PHS formulation studied in [27].

###  5.2 Coupling of identical subsystems

On the mean-field level we can easily couple two interacting particle system
of same type by adding and rescaling their probability measures. Indeed, let
$f^{1},f^{2}$ be two interacting particle systems with identical interaction
behaviour. Then

$\displaystyle\partial_{t}f^{1}=-\nabla_{(x,v)}\cdot\Big{(}\big{(}f^{1}\ast(A\nabla_{f}\mathcal{H}(f^{1}))\big{)}f^{1}\Big{)},\qquad\partial_{t}f^{2}=-\nabla_{(x,v)}\cdot\Big{(}\big{(}f^{2}\ast(A\nabla_{f}\mathcal{H}(f^{2}))\big{)}f^{2}\Big{)}$

Let us now define $f^{\mathrm{sum}}=(f^{1}+f^{2})/(\int df^{1}+\int df^{2})$,
where we rescale to obtain a probability measure. Note that if $f^{1}$ and
$f^{2}$ are probability measures, it holds

$f^{\mathrm{sum}}=\frac{1}{2}(f^{1}+f^{2}).$

The Hamiltonian of the coupled system is given by

$H(f_{t}^{\mathrm{sum}})=\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\left(\frac{1}{2}v^{T}v+\mathcal{V}(x-\bar{x})\right)\mathrm{d}f^{\mathrm{sum}}(t,\bar{x},\bar{v})\mathrm{d}f^{\mathrm{sum}}(t,x,v)$

leading to the dynamics

$\displaystyle\partial_{t}f^{\mathrm{sum}}=-\nabla_{(x,v)}\cdot\Big{(}\big{(}f^{\mathrm{sum}}\ast(A\nabla_{f}\mathcal{H}(f^{\mathrm{sum}}))\big{)}f^{\mathrm{sum}}\Big{)}.$

If we sample now $N$ particles from $f^{1}(0,x,v)$ and $N$ particles from
$f^{2}(0,x,v)$ we obtain the systems

$\displaystyle\frac{d}{dt}x_{i}^{1}$ $\displaystyle=v_{i}^{1},$
$\displaystyle\frac{d}{dt}v_{i}^{1}$
$\displaystyle=-\frac{1}{N}\sum_{j=1}^{N}\psi(|x_{i}^{1}-x_{j}^{1}|)(v_{i}^{1}-v_{j}^{1})-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}(x_{i}^{1}-x_{j}^{1})$

and

$\displaystyle\frac{d}{dt}x_{i}^{2}$ $\displaystyle=v_{i}^{2},$
$\displaystyle\frac{d}{dt}v_{i}^{2}$
$\displaystyle=-\frac{1}{N}\sum_{j=1}^{N}\psi(|x_{i}^{2}-x_{j}^{2}|)(v_{i}^{2}-v_{j}^{2})-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}(x_{i}^{2}-x_{j}^{2}).$

Sampling $2N$ particles from $f^{\mathrm{sum}}$ yields

$\displaystyle\frac{d}{dt}x_{i}$ $\displaystyle=v_{i},$
$\displaystyle\frac{d}{dt}v_{i}$
$\displaystyle=-\frac{1}{2N}\sum_{j=1}^{2N}\psi(|x_{i}-x_{j}|)(v_{i}-v_{j})-\frac{1}{2N}\sum_{j=1}^{2N}\nabla\mathcal{V}(x_{i}-x_{j})$

a fully coupled system. Here, we stack the vectors $x=(x^{1},x^{2})$ and
$v=(v^{1},v^{2})$. The Hamiltonian as well as the system matrices $J$ and $R$
admit the same structure in the dimension $2N$. Note that the generalization
to $K\in\mathbb{N}$ coupled interacting particle systems of same type is
straightforward.

######  Remark.

We want to stress the fact that in general the value of the Hamiltonian of the
coupled systems is greater than the value of the sum of the Hamiltonians of
the subsystems. This is due to the fact that additional generalized springs
are needed to to define the interaction behaviour of the individuals of the
different swarms. However, in case of two identical swarms as above the
rescaling of the Hamiltonian yields that the values of the Hamiltonian
coincide. In the following we describe one approach that allows the coupling
different subsystems in a Hamiltonian preserving way. In fact, the interacting
across subsystems influences only the alignment.

###  5.3 Coupling of different species

Let us consider interacting particle of two different species, which are
modelled with the help of different interaction potentials $\mathcal{V}_{k}$
and $\psi_{k}$ for $k=1,2.$ For simplicity we assume that both subsystem
consist of $N\in\mathbb{N}$ particles. The dynamics of the subsystems read

$\displaystyle\frac{d}{dt}x_{i}^{k}$ $\displaystyle=v_{i}^{k},$
$\displaystyle\frac{d}{dt}v_{i}^{k}$
$\displaystyle=-\frac{1}{N}\sum_{j=1}^{N}\psi_{k}(|x_{i}^{k}-x_{j}^{k}|)(v_{i}^{k}-v_{j}^{k})-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}_{k}(x_{i}^{k}-x_{j}^{k})$

with Hamiltonian

$\mathcal{H}^{N}_{k}(z^{k}(t))=\sum_{i=1}^{N}(v_{i}^{k}(t))^{\top}v_{i}^{k}(t)+\frac{1}{N}\sum_{i,j=1}^{N}\mathcal{V}_{k}(x_{i}^{k}(t)-x_{j}^{k}(t)).$

In order to interconnect the systems in a power conserving way, we have to
define how particles of different species interact with each other. By
$\psi_{c}:\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}$ we denote this interaction
term. Then we obtain

$\displaystyle\frac{d}{dt}x_{i}^{k}$ $\displaystyle=v_{i}^{k},$
$\displaystyle\frac{d}{dt}v_{i}^{k}$
$\displaystyle=-\frac{1}{N}\sum_{j=1}^{N}\psi_{k}(|x_{i}^{k}-x_{j}^{k}|)(v_{i}^{k}-v_{j}^{k})-\frac{1}{N}\sum_{j=1}^{N}\psi_{c}(|x_{i}^{k}-x_{j}^{\ell}|)(v_{i}^{k}-v_{j}^{\ell})$
$\displaystyle\qquad-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}_{k}(x_{i}^{k}-x_{j}^{k}),$

where $i=1,\ldots,N$, $k,\ell\in\\{1,2\\}$ and $k+\ell=3$. The Hamiltonian is
then given by

$\mathcal{H}^{N}(z(t))=\mathcal{H}_{1}^{N}(z^{1}(t))+\mathcal{H}_{2}^{N}(z^{2}(t)).$

Let us define $z=(x^{1},v^{1},x^{2},v^{2}).$ Then

$\frac{\mathrm{d}}{\mathrm{d}t}z=\left[\begin{pmatrix}0&I&0&0\\\ -I&0&0&0\\\
0&0&0&I\\\ 0&0&-I&0\end{pmatrix}-\begin{pmatrix}0&0&0&0\\\
0&\Psi_{1}(z)&0&\Psi_{c}(z)\\\ 0&0&0&0\\\
0&\Psi_{c}^{\top}(z)&0&\Psi_{2}(z)\end{pmatrix}\right]\frac{\partial\mathcal{H}^{N}(z)}{\partial
z}.$ (15)

with

$\displaystyle\Psi_{k}(z)=$
$\displaystyle\frac{1}{N}\\!\begin{pmatrix}\displaystyle\sum_{j=1}^{N}\\!\psi_{k}(|x^{k}_{j}-x^{k}_{1}|)&-\psi_{k}(|x^{k}_{2}-x^{k}_{1}|)&\ldots&-\psi_{k}(|x^{k}_{N}-x^{k}_{1}|)\\\
-\psi_{k}(|x^{k}_{1}-x^{k}_{2}|)&\displaystyle\sum_{j=1}^{N}\\!\psi_{k}(|x^{k}_{j}-x^{k}_{2}|)&\ddots&\vdots\\\
\vdots&\ddots&\ddots&-\psi_{k}(|x^{k}_{N}\\!-\\!x^{k}_{N-1}|)\\\
-\psi_{k}(|x^{k}_{1}-x^{k}_{N}|)&\ldots&-\psi_{k}(|x^{k}_{N-1}\\!-\\!x^{k}_{N}|)&\displaystyle\sum_{j=1}^{N}\\!\psi_{k}(|x^{k}_{j}-x^{k}_{N}|)\end{pmatrix}$
$\displaystyle+\frac{1}{N}\\!\begin{pmatrix}\displaystyle\sum_{j=1}^{N}\\!\psi_{c}(|x^{\ell}_{j}-x^{k}_{1}|)&0&\ldots&0\\\
0&\displaystyle\sum_{j=1}^{N}\\!\psi_{c}(|x^{\ell}_{j}-x^{k}_{2}|)&\ddots&\vdots\\\
\vdots&\ddots&\ddots&0\\\
0&\ldots&0&\displaystyle\sum_{j=1}^{N}\\!\psi_{c}(|x^{\ell}_{j}-x^{k}_{N}|)\end{pmatrix},$

and

$\Psi_{c}(z)=-\frac{1}{N}\\!\begin{pmatrix}\psi_{c}(|x_{1}^{2}-x_{1}^{1}|)&\psi_{c}(|x_{2}^{2}-x_{1}^{1}|)&\ldots&\psi_{c}(|x_{N}^{2}-x_{1}^{1}|)\\\
\psi_{c}(|x_{1}^{2}-x_{2}^{1}|)&\psi_{c}(|x_{2}^{2}-x_{2}^{1}|)&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\psi_{c}(|x_{N}^{2}\\!-\\!x_{N-1}^{1}|)\\\
\psi_{c}(|x_{1}^{2}-x_{N}^{1}|)&\ldots&\psi_{c}(|x_{N-1}^{2}\\!-\\!x_{N}^{1}|)&\psi_{c}(|x_{N}^{2}-x_{N}^{1}|)\end{pmatrix}.$

Now, there are different cases: passing both species to the mean-field limit,
passing just one to the limit and the other remains finite. For the mixed case
we obtain

$\displaystyle\partial_{t}f$ $\displaystyle+\nabla_{x}\cdot\left(vf\right)$
$\displaystyle=\nabla_{v}\cdot\left(\left(\int\psi_{1}(|x-\bar{x}|)(v-\bar{v})+\nabla\mathcal{V}_{1}(x-\bar{x})\mathrm{d}f(t,\bar{x},\bar{v})\right.\right.$
$\displaystyle\qquad\qquad\quad+\left.\left.\frac{1}{N}\sum_{j=1}^{N}\psi_{c}(|x-x_{j}|)(v-v_{j})\right)f\right),$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}x_{i}$ $\displaystyle=v_{i},$
$\displaystyle\frac{d}{dt}v_{i}$
$\displaystyle=-\frac{1}{N}\sum_{j=1}^{N}\psi_{2}(|x_{i}-x_{j}|)(v_{i}-v_{j})-\int\psi_{c}(|x_{i}-\bar{x}|)(v_{i}-\bar{v})\mathrm{d}f(\bar{x},\bar{v})$
$\displaystyle\qquad-\frac{1}{N}\sum_{j=1}^{N}\nabla\mathcal{V}_{2}(x_{i}^{k}-x_{j}^{k}).$

The Hamiltonian is a combination as well
$\mathcal{H}(f_{t},z(t))=\mathcal{H}(f_{t})+\mathcal{H}^{N}(z(t))$ as the
coupling across subsystems affects only the alignment terms, see Remark Remark
for more details.

##  6 Conclusion and outlook

We derived a minimal port-Hamiltonian formulation of interacting particle
systems and showed that the port-Hamiltonian structure is preserved in the
mean-field limit. The Hamiltonian is used as Lyapunov function to characterize
the long-time behavior of the systems on the particle and the mean-field
level. Hence the PHS formulation opens a new perspective on the well-studied
particle and mean-field description of interacting particle system. Moreover,
the identification of Casimir functions prepares the ground to define port-
Hamiltonian preserving control strategies in future work. On the other hand,
the LaSalle-type argument for the long-term behavior may open the door for
convergence results for general consensus dynamics for optimization and
sampling tasks in the spirit of [10, 37].

## Acknowledgments

We thank the anonymous referee for their critical and constructive comments
that helped us to strengthen the results on the long-time behavior in both,
the ODE and the PDE setting.

## References

  * [1] S. M. Ahn, H. Choi, S.-Y. Ha, and H. Lee.  On collision-avoiding initial configurations to cucker-smale type flocking models.  Comm. Math. Sci., 10(2):625–643, 2012.
  * [2] G. Albi, D. Balagué, J. A. Carrillo, and J. von Brecht.  Stability analysis of flock and mill rings for second order models in swarming.  SIAM J Appl. Math., 74(3):794–818, 2014.
  * [3] G. Albi and L. Pareschi.  Modeling of self-organized systems interacting with a few individuals: from microscopic to macroscopic dynamics.  Appl. Math. Lett., 26(4):397–401, 2013.
  * [4] L. Ambrosio, N. Gigli, and G. Savaré.  Gradient flows: in metric spaces and in the space of probability measures.  Springer Science & Business Media, 2005.
  * [5] A. B. Barbaro, J. A. Canizo, J. A. Carrillo, and P. Degond.  Phase transitions in a kinetic flocking model of cucker–smale type.  Multiscale Modeling Simul., 14(3):1063–1088, 2016.
  * [6] M. Burger, R. Pinnau, C. Totzeck, and O. Tse.  Mean-field optimal control and optimality conditions in the space of probability measures.  SIAM J Control Opt., 59(2):977–1006, 2021.
  * [7] M. Burger, R. Pinnau, C. Totzeck, O. Tse, and A. Roth.  Instantaneous control of interacting particle systems in the mean-field limit.  J Comp. Phys., 405:109181, 2020.
  * [8] F. Cao, S. Motsch, A. Reamy, and R. Theisen.  Asymptotic flocking for the three-zone model.  Mathematical Biosciences and Engineering, 2020.
  * [9] J. Carrillo, M. Fornasier, G. Toscani, and F. Vecil.  Particle, kinetic, and hydrodynamic models of swarming.  In T. G. Naldi G., Pareschi L., editor, Mathematical Modeling of Collective Behavior in Socio-Economic and Life Sciences. Springer, 2010.
  * [10] J. Carrillo, F. Hoffmann, A. Stuart, and U. Vaes.  Consensus-based sampling.  Studies in Appl. Math., 148(3):1069–1140, 2022.
  * [11] J. A. Carrillo, Y.-P. Choi, P. B. Mucha, and J. Peszek.  Sharp conditions to avoid collisions in singular cucker–smale interactions.  Nonlinear Anal. Real World Appl., 37:317–328, 2017.
  * [12] J. A. Carrillo, Y.-P. Choi, and S. P. Perez.  A review on attractive–repulsive hydrodynamics for consensus in collective behavior.  Active Particles, Volume 1, pages 259–298, 2017.
  * [13] J. A. Carrillo, M. Fornasier, J. Rosado, and G. Toscani.  Asymptotic flocking dynamics for the kinetic cucker–smale model.  SIAM J Math. Anal., 42(1):218–236, 2010.
  * [14] J. A. Carrillo, M. Fornasier, G. Toscani, and F. Vecil.  Particle, kinetic, and hydrodynamic models of swarming.  In Mathematical modeling of collective behavior in socio-economic and life sciences, pages 297–336. Springer, 2010.
  * [15] J. A. Cañizo, J. A. Carrillo, and J. Rosado.  A well-posedness theory in measures for some kinetic models of collective motion.  Math. Models Methods Appl. Sci., 21(03):515–539, 2011.
  * [16] J. Cho, S.-Y. Ha, F. Huang, C. Jin, and D. Ko.  Emergence of bi-cluster flocking for the cucker–smale model.  Math. Models Methods Appl. Sci., 26(06):1191–1218, 2016.
  * [17] Y.-P. Choi and J. Haskovec.  Cucker-smale model with normalized communication weights and time delay.  arXiv preprint arXiv:1608.06747, 2016.
  * [18] F. Cucker and J.-G. Dong.  A general collision-avoiding flocking framework.  IEEE Trans. Autom. Control, 56(5):1124–1129, 2011.
  * [19] F. Cucker and S. Smale.  On the mathematics of emergence.  Japan J. Math., 2:197–227, 2007.
  * [20] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes.  Self-propelled particles with soft-core interactions: Patterns, stability, and collapse.  Phys. Rev. Lett., 96:104302, Mar 2006.
  * [21] V. Duindam, A. Macchelli, S. Stramigioli, and H. Bruyninckx, editors.  Modeling and Control of Complex Physical Systems.  Springer, Germany, 2009.
  * [22] D. Eberard, B. M. Maschke, and A. J. van der Schaft.  An extension of Hamiltonian systems to the thermodynamic phase space: towards a geometry of nonreversible processes.  Rep. Math. Phys., 60(2):175–198, 2007.
  * [23] R. Erban, J. Haskovec, and Y. Sun.  A cucker–smale model with noise and delay.  SIAM J Appl. Math., 76(4):1535–1557, 2016.
  * [24] F. Golse.  On the dynamics of large particle systems in the mean field limit.  In Macroscopic and large scale phenomena: coarse graining, mean field limits and ergodicity, pages 1–144. Springer, 2016.
  * [25] S.-Y. Ha, J. Kim, J. Park, and X. Zhang.  Complete cluster predictability of the cucker–smale flocking model on the real line.  Arch. Ration. Mech. Anal., 231(1):319–365, 2019.
  * [26] D. Hinrichsen and A. J. Pritchard.  Mathematical systems theory. I, volume 48 of Texts in Applied Mathematics.  Springer-Verlag, Berlin, 2005.
  * [27] I. Matei, C. Mavridis, J. Baras, and M. Zhenirovskyy.  Inferring particle interaction physical models and their dynamical properties.  IEEE 58th Conference on Decision and Control, 2019.
  * [28] S. Motsch and E. Tadmor.  A new model for self-organized dynamics and its flocking behavior.  Journal Stat. Phys., 144:923–947, 2011.
  * [29] S. Motsch and E. Tadmor.  Heterophilious dynamics enhances consensus.  SIAM Review, 56(4):577–621, 2014.
  * [30] V. M. Panaretos and Y. Zemel.  An invitation to statistics in Wasserstein space.  Springer Nature, 2020.
  * [31] J. Park, H. J. Kim, and S.-Y. Ha.  Cucker-smale flocking with inter-particle bonding forces.  IEEE Trans. Autom. Control, 55(11):2617–2623, 2010.
  * [32] C. Pignotti and E. Trélat.  Convergence to consensus of the general finite-dimensional cucker-smale model with time-varying delays.  arXiv preprint arXiv:1707.05020, 2017.
  * [33] C. W. Reynolds.  Flocks, herds and schools: A distributed behavioral model.  In Proceedings of the 14th annual conference on Computer graphics and interactive techniques, pages 25–34, 1987.
  * [34] R. Shvydkoy.  Dynamics and Analysis of Alignment Models of Collective Behavior.  Nečas Center Series. Birkhäuser, 2021.
  * [35] H.-A. Tanaka, A. J. Lichtenberg, and S. Oishi.  First order phase transition resulting from finite inertia in coupled oscillator systems.  Physical review letters, 78(11):2104, 1997.
  * [36] G. Teschl.  Ordinary differential equations and dynamical systems, volume 140\.  American Mathematical Soc., 2012.
  * [37] C. Totzeck.  Trends in consensus-based optimization.  In Active Particles, Volume 3, pages 201–226. Springer, 2022.
  * [38] A. Van der Schaft.  L2-gain and passivity techniques in nonlinear control.  Springer, 2000.
  * [39] A. van der Schaft.  Port-Hamiltonian systems: an introductory survey.  Proceedings on the International Congress of Mathematicians, Vol. 3, pags. 1339-1366, 01 2006.
  * [40] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet.  Novel type of phase transition in a system of self-driven particles.  Phy. Rev. Lett., 75(6):1226, 1995.