Particle dynamics subject to impenetrable boundaries: existence and uniqueness of mild solutions

TL;DR

This paper establishes the existence and uniqueness of mild solutions for particle systems confined by impenetrable boundaries, providing a new approach to analyze such systems with boundary interactions.

Contribution

It introduces the concept of mild solutions for particle dynamics with boundary collisions, simplifying analysis and enabling proofs of uniqueness and many-particle limits.

Findings

Existence of mild solutions for particle systems with boundary interactions.

Mild solutions facilitate proving uniqueness and limits.

Numerical simulations demonstrate the approach on non-convex domains.

Abstract

We consider the dynamics of particle systems where the particles are confined by impenetrable barriers to a bounded, possibly non-convex domain $\Omega$. When particles hit the boundary, we consider an instant change in velocity, which turns the systems describing the particle dynamics into an ODE with discontinuous right-hand side. Other than the typical approach to analyse such a system by using weak solutions to ODEs with multi-valued right-hand sides (i.e., applying the theory introduced by Filippov in 1988), we establish the existence of mild solutions instead. This solution concept is easier to work with than weak solutions; e.g., proving uniqueness of mild solutions is straight-forward, and mild solutions provide a solid structure for proving many-particle limits. We supplement our theory of mild solutions with an application to gradient flows of interacting particle energies with a singular interaction potential, and illustrate its features by means of numerical simulations on various choices for the (non-convex) domain $\Omega$.

Equations

$$m_0\ddot{x}(t)+\mu \dot{x}(t)=f(t)$$

$$\dot{x}(t)=\left\{\begin{array}[]{ll}F(t)&\mbox{if $x(t)\in\Omega$,}\\[5.0pt] F(t)-(F(t)\cdot\nu)_{+}\,\nu&\mbox{if $x(t)\in\partial\Omega$.}\end{array}\right.$$

$$\dot{x}_{i}(t)=\left\{\begin{aligned} &F_{i}\big{(}t,X(t)\big{)}&x_{i}(t)&\in\Omega,\\ &F_{i}\big{(}t,X(t)\big{)}-\big{(}F_{i}\big{(}t,X(t)\big{)}\cdot\nu(x_{i}(t))\big{)}_{+}\,\nu(x_{i}(t))&x_{i}(t)&\in\partial\Omega.\end{aligned}\right.$$

$$\Vert X\Vert :=\underset{1\le i\le n}{\max }|x_i|,$$

$$G:={\overline{\Omega }}^n$$

$$F_i\in C([0,T]\times G;\mathbb{R})(i=1,\dots ,n).$$

$$|F_i(t,X)-F_i(t,Y)|\le L\Vert X-Y\Vert (t\in [0,T],X,Y\in ({\mathbb{R}}^m{)}^n,i=1,\dots ,n).$$

$$P:\overline{\Omega}\times{\mathbb{R}}^{m}\to{\mathbb{R}}^{m},\quad P(x,f):=\left\{\begin{aligned} &f-\big{(}f\cdot\nu(x)\big{)}\nu(x)&&(x\in\partial\Omega,f\cdot\nu(x)>0),\\ &f&&\text{otherwise},\end{aligned}\right.$$

$$H_{i}:[0,T]\times G\to{\mathbb{R}}^{m},\quad H_{i}\big{(}t,X\big{)}:=P\big{(}x_{i},F_{i}(t,X)\big{)}\quad(t\in[0,T],X=\{x_{i}\}\in G).$$

$$\left\{\begin{aligned} &\dot{x}_{i}(t)=H_{i}\big{(}t,X(t)\big{)}\quad(0<t<T),\\ &x_{i}(0)=x_{i}^{\circ}.\end{aligned}\right.$$

$$x_{i}(t)=x_{i}^{\circ}+\int_{0}^{t}H_{i}\big{(}s,X(s)\big{)}ds\quad\big{(}t\in[0,T)\big{)},$$

$$\left\{\begin{aligned} &x_{i}\in\mathrm{Lip}([0,T];{\mathbb{R}}^{m})\\ &\dot{x}_{i}(t)=H_{i}(t,X(t))\quad\big{(}\text{a.e.~{}}t\in[0,T]\big{)}\\ &x_{i}(0)=x_{i}^{\circ}.\end{aligned}\right.$$

$$\text{dist}(x,K):=\underset{y\in K}{\inf }|x-y|,$$

$$N^{\epsilon }(K):=\{x\in {\mathbb{R}}^m;\text{dist}(x,K)<\epsilon \}.$$

$$d(x):=\mathrm{dist}(x,\overline{\Omega })(x\in {\mathbb{R}}^m),$$

$$d_s(x):=\{\begin{array}{rl}& \mathrm{dist}(x,\partial \Omega )(x\notin \overline{\Omega }),\\ & -\mathrm{dist}(x,\partial \Omega )(x\in \overline{\Omega }).\end{array}$$

$$\nabla d_{s}(x)=\nu(\bar{x}),~{}\bar{x}=x-d_{s}(x)\nabla d_{s}(x)\in\partial\Omega\quad\big{(}x\in\overline{N^{\varepsilon}(\partial\Omega)}\big{)},$$

$$d\nabla d(x):=\{\begin{array}{rlr}& d(x)\nabla d(x)& (x\in \overline{N^{\epsilon }(\partial \Omega )}\setminus \overline{\Omega }),\\ & 0& (x\in \overline{\Omega }).\end{array}$$

$$d_{s}^{\varepsilon}(x)\left\{\begin{aligned} &=d_{s}(x)&\big{(}x&\in N^{\varepsilon}(\partial\Omega)\big{)},\\ &\geq\varepsilon&\big{(}x&\in{\mathbb{R}}^{m}\setminus N^{\varepsilon}(\Omega)\big{)},\\ &\leq-\varepsilon&\big{(}x&\in\Omega\setminus N^{\varepsilon}(\partial\Omega)\big{)}.\end{aligned}\right.$$

$$d\nabla d\in \mathrm{Lip}(\overline{N^{\epsilon }(\Omega )};{\mathbb{R}}^m).$$

$$|d\nabla d(x)-d\nabla d(y)|\le C^{\prime }|x-y|$$

$$|d\nabla d(x)-d\nabla d(y)|=|d(y)\nabla d(y)|=d(y)=\underset{z\in \overline{\Omega }}{\min }|z-y|\le |x-y|.$$

$$\Vert F{\Vert }_{\infty }:=\underset{t\in [0,T]}{\max }\underset{X\in ({\mathbb{R}}^m{)}^n}{\max }\underset{1\le i\le n}{\max }|F_i(t,X)|.$$

$$(a{)}_{+}:=\max (0,a),(a{)}_{-}:=\max (0,-a).$$

$$\int_{A}|u^{\prime}(t)|dt=\int_{{\mathbb{R}}}\mathcal{H}^{0}\big{(}u^{-1}(y)\big{)}dy,$$

$$\{\begin{array}{rl}& {\dot{x}}_i^k(t)=F_i(t,X^k(t))-k(d\nabla d)(x_i^k(t))(0<t<T)\\ & x_i^k(0)=x_i^{\circ }.\end{array}$$

$$\frac{d}{dt}d(x_i^k(t))$$

$$=\nabla d(x_i^k(t))\cdot F_i(t,X^k(t))-kd(x_i^k(t)),$$

$$\displaystyle\frac{d}{dt}\big{(}e^{kt}d(x_{i}^{k}(t))\big{)}$$

$$\displaystyle=e^{kt}\big{(}\nabla d(x_{i}^{k}(t))\cdot F_{i}(t,X^{k}(t))\big{)}$$

$$\le |F_i(t,X^k(t))|e^{kt}$$

$$\le \Vert F{\Vert }_{\infty }{e}^{kt}.$$

$$e^{kt}d(x_i^k(t))\le e^{k{t}_0}d(x_i^k(t_0))+\Vert F{\Vert }_{\infty }{\int }_{t_0}^t{e}^{ks}ds=\Vert F{\Vert }_{\infty }\frac{e^{kt}-e^{k{t}_0}}{k},$$

Full text

We supplement our theory of mild solutions with an application to gradient flows of interacting particle energies with a singular interaction potential, and illustrate its features by means of numerical simulations on various choices for the (non-convex) domain $\Omega$.

1. 1 Introduction

2. 2 Discontinuous ODE system: existence and uniqueness

3. 2.1 Problem setting and main theorem

4. 2.2 Distance functions with respect to $\Omega$ and $\partial\Omega$

5. 2.3 Preliminaries and notation

6. 2.4 The approximate problem to Problem 2.1

7. 2.5 Existence of solutions

8. 2.6 Uniqueness of the solution

9. 2.7 Properties of the solutions

10. 3 Gradient flows of particle interaction energies

11. 3.1 Particle interaction energy and its minimization problem

12. 3.2 Gradient flow of the interaction energy

13. 3.3 Proof of Theorem 3.4: existence and uniqueness

14. 4 Numerical examples

15. 4.1 Circular domain $\Omega$

16. 4.2 A non-convex domain $\Omega$ (dumbbell)

17. 4.3 Channels with obstacles

18. 5 Conclusion and outlook

1 Introduction

In various case studies in the field of interacting particle systems, in particular SPH [GM77, LL04, LL10, Luc77, Mon92, AASS12], pedestrian dynamics [MV11, BBK13, DH10, HINT03, LSZ*+*14] and granular media [Ste98, Mau06], the particles are confined to a bounded, non-convex domain $\Omega\subset{\mathbb{R}}^{m}$ (see Figure 1). While there are various ways to implement such impenetrable boundaries in numerical schemes for computing particle dynamics, a satisfactory analytical framework seems to be largely missing. Such a missing structure makes, for example, any statement on many-particle limits out of reach. Our aim is therefore to set up an analytical framework for interacting particle systems which can handle impenetrable barriers, to define a satisfactory notion of a solution for the particle trajectories, and to prove existence and uniqueness of such solutions.

In this paper, we focus on first-order in time particle systems. To motivate formally and state the related system of ODEs, we first consider the motion of a point mass $x(t)$ with mass $m_{0}>0$ confined to $\Omega$. We subject the point mass to a given force $f(t)\in{\mathbb{R}}^{m}$ and impose a viscous resistance represented by a linear drag force $-\mu\dot{x}$ with a positive coefficient $\mu>0$. Then, using Newton’s second law, the equation of motion for the point mass is

[TABLE]

whenever $x(t)\in\Omega$. Next we assume that $m_{0}/\mu\ll 1$, and consider the overdamped limit. Setting $m_{0}=0$ and $F(t):=f(t)/\mu$ (we keep on calling it a “force” as in the above context), we obtain for the dynamics of the point mass that $\dot{x}(t)=F(t)$. Finally we motivate and include the effect of the confining domain $\Omega$. When $x(t)\in\partial\Omega$, from Newton’s third law, a counter force $-(F\cdot\nu)_{+}\,\nu$ acts on the particle from $\partial\Omega$, where $\nu$ is the outward normal vector on $\partial\Omega$ and $(a)_{+}:=\max(a,0)$. This leads to the following ODE with a “one-sided projection”:

[TABLE]

Using the overdamped dynamics of the point mass, we introduce the dynamics of the particle system. Denoting the positions of $n$ particles by $X(t)=\{x_{i}(t)\}_{t=1}^{n}\in\overline{\Omega}^{n}$, we consider a given force F_{i}\big{(}t,X(t)\big{)} which acts on particle $i$ for each $i=1,\ldots,n$. With the confinement described above, this yields the following set of ODEs:

[TABLE]

Because of the one-sided projection of the force at the boundary, the trajectories of particles typically have a discontinuity in their velocity at the time when they hit the boundary. Since classical $C^{1}$-solutions to (1.1) cannot describe such behavior, we cannot rely on classical ODE theory such as the Cauchy-Peano Theorem and the Picard-Lindel$\ddot{o}$f (Cauchy-Lipschitz) Theorem [Har82] for establishing the existence and uniqueness of solutions to (1.1).

Next we briefly review the literature on the existence and uniqueness of solutions to systems with a discontinuous right-hand side such as (1.1). When $\Omega$ is convex and $F_{i}$ is the gradient of a $\lambda$-convex energy with respect to $x_{i}$, then gradient flow theory [AGS08] provides a sufficiently strong analytical framework to describe (1.1). In the general case ($\Omega$ not convex or a more general discontinuous right-hand side than the one in (1.1)), the book by Filippov [Fil88] treats such discontinuities by making the right-hand side multi-valued. For the related differential inclusion, the existence and uniqueness of weak solutions is sought. Especially the conditions for the right-hand side to result in uniqueness of solutions can be difficult to check, and need to be considered differently on a case-to-case basis. Moreover, weak solutions to differential inclusions seem to lend themselves poorly for the purpose of studying many-particle limits.

As opposed to the general setting in [Fil88], (1.1) has a particular type of discontinuity in the right-hand side which, under appropriate regularity conditions of $F_{i}$ and $\Omega$, allows for a stronger sense of solutions. Indeed, our main result (Theorem 2.4) states that for enough regular $F_{i}$ and $\Omega$, (1.1) attains a unique, mild solution, without the need to transform (1.1) into a differential inclusion. In addition, we show that when particles detach from the boundary, their velocity is tangent to $\partial\Omega$ (Propsition 2.14).

Regarding the regularity condition on $F$, we require $F$ to be continuous for the existence of mild solutions. However, in many applications of interacting particle systems, $F$ blows up whenever two particles tend to the same position. To cover such scenarios, we extend Theorem 2.4 to the case in which $F$ is the gradient of an energy which blows up to $+\infty$ as any two particles tend to the same position (as is the case, for instance, for charged particles). This extension is our second main result, which is stated in Theorem 3.4.

The intended purpose of Theorem 3.4 is to be a stepping stone for proving the limit passage $n\to\infty$ (i.e., the many-particle limit) from (1.1) to a nonlocal and nonlinear equation for the particle density. The singularity of the particle interaction potential $V$ at [math] plays a crucial role here, since it has a regularising effect on the particle density. In recent studies, such limit passages $n\to\infty$ have been carried out for:

• •

convex $V$ (except at [math]) in one spatial dimension in [vMM14],

• •

logarithmic $V$ in two spatial dimensions in [Sch96], which was later extended to stronger singularities in [Due16],

• •

sub-Coulombic potentials (i.e., $V(x)\ll|x|^{2-m}$ as $x\to 0$) in arbitrary spatial dimension $m$ in [Hau09],

• •

multi-species particle systems with regular cross-interactions in [vM18, vMM18],

• •

multi-species particle systems with singular cross-interactions in [CXZ16, GvMPS18],

• •

a rate-independent velocity law in [MPS17].

Besides singular potentials, it is also possible to obtain a regular particle density by adding white noise to the particle dynamics. It has been proven in [Szn91, Phi07] that such systems satisfy the propagation of chaos property, which implies that the limit can be described by a deterministic equation for the particle density. However, in all these studies, the problem of particles attaching to or detaching from boundaries is side-stepped. With Theorem 3.4 in hand, we expect that it is possible to extend the techniques in the references cited above to establish the limit passage $n\to\infty$ in (1.1). We leave this challenge for future research.

Another interesting generalisation of our result is the extension to second order ODE systems with inertial term. This generalisation has been studied by [Sch78] for convex energies and convex domains. We believe that our current result provides a new viewpoint from which Schatzman’s result can be extended to non-convex energies (which allows, e.g., for repulsive particle interactions) and non-convex domains. An important application that we have in mind for this is a properly motivated boundary condition in particle simulations for fluids. Indeed, in many such simulations (in particular when SPH is used), there is no established theory on how to treat properly boundary conditions for the particles which effectively model slip and no-slip conditions in the equations for the fluid.

To gain insight in the solution trajectories of (1.1) (in particular, how particles attach to and detach from the boundary), we perform several numerical simulations of the solution to (1.1). We study the case where the force $F_{i}$ is the gradient of an unbounded interaction energy with an external potential. Our interest lies in the behavior of the system for several different, non-convex domains $\Omega$.

The paper is organised as follows. In Section 2 we state and prove Theorem 2.4, which provides the existence and uniqueness of mild solutions to (1.1). In Section 3 we extend Theorem 2.4 to Theorem 3.4, which allows for a class of discontinuous $F$. In Section 4 we perform and discuss our numerical simulations. We conclude in Section 5.

2 Discontinuous ODE system: existence and uniqueness

2.1 Problem setting and main theorem

Let $\Omega$ be a bounded domain in ${\mathbb{R}}^{m}$ with $C^{2}$ boundary $\partial\Omega$. We denote the outward unit normal vector on $\partial\Omega$ by $\nu$.

We consider $n$ particles $x_{1},\dots,x_{n}$ in $\overline{\Omega}$. We put $X:=\{x_{i}\}_{i=1}^{n}\in\overline{\Omega}^{n}\subset(R^{m})^{n}$ as the array of the particle’s positions. On $({\mathbb{R}}^{m})^{n}$ we define the norm

[TABLE]

where $|x_{i}|$ is the Euclidean distance in ${\mathbb{R}}^{m}$. We interpret

[TABLE]

as the set of admissible particle positions.

Next we introduce the dynamics of the particle system. We denote by $X(t)=\{x_{i}(t)\}_{i=1}^{n}$ the particle positions in time, and by $\dot{x_{i}}(t)=\frac{dx_{i}}{dt}(t)$ their velocity. We prescribe the velocity of any particle $x_{i}(t)$ depending on whether $x_{i}(t)$ is in the interior or at the boundary of $\Omega$. If $x_{i}(t)\in\Omega$, then its velocity is given by the force $F_{i}(t,X(t))$, where

[TABLE]

In some cases, we impose the following hypothesis on Lipschitz continuity of $F_{i}$: there exists an $L>0$ such that

[TABLE]

If $x_{i}(t)\in\partial\Omega$, then we prescribe a different rule. If the force $F_{i}(t,X(t))$ acting on $x_{i}(t)$ is directed to the interior of $\Omega$, then we prescribe its velocity by $F_{i}(t,X(t))$. Otherwise, we project $F_{i}(t,X(t))$ onto the tangent plane of $\partial\Omega$ at $x_{i}(t)$ (see Figure 2). To define this rule mathematically, we put

[TABLE]

and prescribe the velocity of $x_{i}(t)$ as $P(x_{i}(t),F_{i}(t,X))$. For convenience, we set

[TABLE]

To summarise the dynamics of the particles, we introduce the following problem:

Problem 2.1**.**

Given $T>0$ and an initial condition $X^{\circ}=\{x_{i}^{\circ}\}_{i=1}^{n}\in G$, for $i=1,2,\dots,n$, find a solution to

[TABLE]

Since $H_{i}$ is discontinuous, we do not expect solutions to Problem 2.1 to be $C^{1}$ in time. Instead, we consider mild solutions:

Definition 2.2**.**

We call $X(t)$ a (mild) solution to Problem 2.1 if it satisfies $x_{i}\in C^{0}([0,T];{\mathbb{R}}^{m})$, $x_{i}(t)\in\overline{\Omega}$ and

[TABLE]

for $i=1,2,\dots,n$.

The following proposition gives a convenient characterisation of mild solutions. The proof is obvious.

Proposition 2.3**.**

For Problem 2.1, $X(t)=\{x_{i}(t)\}_{i=1}^{n}$ is a mild solution if and only if

[TABLE]

Our main result is the following theorem, which grants the existence and uniqueness of mild solutions to Problem 2.1:

Theorem 2.4**.**

Let $T>0$ and $X^{\circ}\in G$. If $F_{i}$ satisfies (2.2), then there exists a solution to Problem 2.1. Moreover, if $F_{i}$ satisfies the Lipschitz condition (2.3), then the solution is unique.

This theorem is a natural extension of Cauchy-Peano theorem and the Picard-Lindel$\ddot{o}$f theorem to particles systems with impenetrable boundaries. The proof for the existence and the uniqueness is given in Sections 2.5 and 2.6, respectively.

2.2 Distance functions with respect to $\Omega$ and $\partial\Omega$

In the proof of Theorem 2.4, we define an approximation problem in which particles can leave $\overline{\Omega}$. To account for this, we define in this section distance functions on ${\mathbb{R}}^{m}$ with respect to $\Omega$ and $\partial\Omega$. For $K\subset{\mathbb{R}}^{m}$ and $\varepsilon>0$, we defne

[TABLE]

[TABLE]

For the domain $\Omega$, we define the distance function $d$ with respect to $\Omega$ and the signed distance function $d_{s}$ with respect to $\partial\Omega$ as

[TABLE]

Since $\partial\Omega\in C^{2}$, there exists $\varepsilon>0$ such that $d\in C^{2}(\overline{N^{\varepsilon}(\partial\Omega)}\setminus\Omega)$ and $d_{s}\in C^{2}(\overline{N^{\varepsilon}(\partial\Omega)})$. We fix this $\varepsilon$ throughout Section 2.5 and Section 2.6. It is known (see, e.g., [GT85, Kim08]) that the signed distance function $d_{s}$ satisfies

[TABLE]

where $\bar{x}\in\partial\Omega$ denotes the orthogonal projection from $x\in\overline{N^{\varepsilon}(\partial\Omega)}$ onto $\partial\Omega$. We set

[TABLE]

Let $\varepsilon>0$ as above. We put $d_{s}^{\varepsilon}\in C^{2}({\mathbb{R}}^{m})$ as a smooth extension of $d_{s}$ outside $N^{\varepsilon}(\partial\Omega)$:

[TABLE]

We note that, by construction, $\nabla d_{s}^{\varepsilon}$ is a $C^{1}$-extension of $\nu$ from $\partial\Omega$ to ${\mathbb{R}}^{m}$. Finally, we establish the following Lipschitz property of $d\nabla d$:

Proposition 2.5**.**

The function $d\nabla d(x)$ defined in (2.5) satisfies

[TABLE]

Proof.

We show that there exists a constant $C^{\prime}>0$ such that for all $x,~{}y\in\overline{N^{\varepsilon}(\Omega)}$,

[TABLE]

holds. We separate three cases.

Case 1: $x,~{}y\in\overline{\Omega}$. Since $d\nabla d\equiv 0$ on $\overline{\Omega}$, (2.6) is obvious.

Case 2: $x,~{}y\in\overline{N^{\varepsilon}(\Omega)}\setminus\overline{\Omega}$. Since $d\nabla d(x)\in C^{1}(N^{\varepsilon}(\partial\Omega)\setminus\overline{\Omega})$, (2.6) is obvious.

Case 3: $x\in\overline{\Omega}$ and $y\in\overline{N^{\varepsilon}(\Omega)}\setminus\overline{\Omega}$. We conclude (2.6) from

[TABLE]

∎

2.3 Preliminaries and notation

In this section, we introduce some notation for the proof of Theorem 2.4, and recall some well-known theorems.

Using the assumed regularity on $\Omega$ and $F_{i}$, we extend $F_{i}$ to $[0,\infty)\times({\mathbb{R}}^{m})^{n}$ such that $F_{i}$ is still continuous. In the remainder of the paper, we write $F_{i}$ for this extension.

We define the $L^{\infty}$-norm of $F:=(F_{1},\ldots,F_{n})$ as

[TABLE]

We set $L^{\infty}(0,T;[0,1]):=\{\chi\in L^{\infty}(0,T)~{};~{}\chi(t)\in[0,1]~{}\text{for a.e.}~{}t\}$, and denote the positive and negative part of a number or function by

[TABLE]

As preparation for the proof of Theorem 2.4, we recall two well-known results.

Theorem 2.6** (Ascoli-Arzela).**

Let $(\mathcal{X},\mathsf{m})$ be a compact metric space. A sequence $\{f_{n}\}\in C(\mathcal{X})$ has a uniformly convergent subsequence if it is bounded and equicontinuous, i.e., $\sup_{n}\|f_{n}\|<\infty$ and for every $\varepsilon>0$ there exists $\delta>0$ such that for all $n$ and all $x,~{}y\in\mathcal{X}$ with $\mathsf{m}(x,y)<\delta$, it holds that $|f_{n}(x)-f_{n}(y)|<\varepsilon$.

Theorem 2.7** (Co-area formula on ${\mathbb{R}}$ [EG92]).**

Let $A\subset{\mathbb{R}}$ be a Borel set, and $u:A\to{\mathbb{R}}$ be Lipschitz continuous. Then, y\mapsto\mathcal{H}^{0}\big{(}u^{-1}(y)\big{)} is measurable, and

[TABLE]

where $\mathcal{H}^{0}$ is the zero dimensional Hausdorff measure (counting measure) on ${\mathbb{R}}$.

2.4 The approximate problem to Problem 2.1

In preparation for proving the existence of solutions to Problem 2.1, we define Problem 2.8 to construct a sequence $\{X^{k}(t)\}_{k}$ of approximate solutions. We recall that we extended $F_{i}$ from $[0,T]\times G$ to $[0,\infty)\times({\mathbb{R}}^{m})^{n}$ in Section 2.3.

Problem 2.8**.**

Given $k\in{\mathbb{N}}$ and $X^{\circ}\in G$, find $X^{k}(t)=\{x_{i}^{k}(t)\}_{i=1}^{n}\in({\mathbb{R}}^{m})^{n}$ such that

[TABLE]

Figure 3 illustrates a typical situation for a solution to Problem 2.8. The following lemma describes properties of solutions to Problem 2.8.

Lemma 2.9**.**

Let $T>0$ and suppose that (2.2) holds and that $k>\frac{1}{\varepsilon}\|F\|_{\infty}$. Then for any $X^{\circ}\in G$, there exists a classical solution to Problem 2.8. Furthermore, any solution $X^{k}(t)=\{x_{i}^{k}(t)\}$ to Problem 2.8 satisfies $x_{i}^{k}(t)\in N^{\varepsilon}(\Omega)$ and $|\dot{x}_{i}^{k}(t)|\leq 2\|F\|_{\infty}$ for all $i=1,2,\dots,n$ and all $t\in[0,T]$.

Proof.

Since $F_{i}$ and $d\nabla d$ are continuous in $\overline{N^{\varepsilon}(\Omega)}$ and $X^{\circ}\in N^{\varepsilon}(\Omega)$, it follows from the Cauchy-Peano theorem that there exists $T_{0}>0$ for which Problem 2.8 has a classical solution $x_{i}^{k}(t)$ on $[0,T_{0}]$ such that $x_{i}^{k}(t)\in N^{\varepsilon}(\Omega)$ for all $0\leq t<T_{0}$ and such that either $x_{i}^{k}(T_{0})\in\partial N^{\varepsilon}(\Omega)$ for some $i\in\{1,2,\dots,n\}$ or $T_{0}>T$.

Next, we claim that $x_{i}^{k}(T_{0})\in\partial N^{\varepsilon}(\Omega)$ is impossible, from which it follows that $T_{0}>T$. To prove this claim, we suppose that $x_{i}^{k}(t)\in N^{\varepsilon}(\Omega)\setminus\overline{\Omega}~{}(0\leq t_{0}<t<t_{1}\leq T_{0})$ and $x_{i}^{k}(t_{0})\in\partial\Omega$ hold for some $i\in\{1,2,\dots,n\}$ and $[t_{0},t_{1}]\subset[0,T_{0}]$. We first compute, for any $t\in(t_{0},t_{1})$,

[TABLE]

and

[TABLE]

Integrating the above estimate from $t_{0}$ to $t\in(t_{0},t_{1}]$, we obtain

[TABLE]

and thus

[TABLE]

Hence, $x_{i}^{k}(t)\in N^{\varepsilon}(\Omega)$ for all $i\in\{1,2,\dots,n\}$ and $t\in[0,T_{0}]$.

Finally, we prove the bound on the velocity. By the computation above, $x_{i}^{k}(t)\in N^{\varepsilon}(\Omega)$ for all $t\in[0,T]$, and

[TABLE]

∎

2.5 Existence of solutions

In this section, we prove the existence statement of Theorem 2.4. Let $T>0$ and the initial condition $X^{\circ}=\{x_{i}^{\circ}\}_{i=1}^{n}\in G$ be given. Let $k_{0}\in{\mathbb{N}}$ with $k_{0}>\frac{1}{\varepsilon}\|F\|_{\infty}$, and $X^{k}(t)$ be the corresponding solution to Problem 2.8 for all $k\geq k_{0}$. We use limit points of the sequence $X^{k}$ as candidates for solutions of Problem 2.1. With this aim, we first specify the topology in which we seek such limit points (Lemma 2.10).

Lemma 2.10**.**

There exists a subsequence of $\{X^{k}\}_{k=k_{0}}^{\infty}\subset C^{1}([0,T];({\mathbb{R}}^{m})^{n})$ (not relabelled) and $X=\{x_{i}(\cdot)\}_{i=1}^{n}\in W^{1,\infty}(0,T;({\mathbb{R}}^{m})^{n})\cap C^{0}([0,T];G)$ such that

[TABLE]

for all $i=1,\dots,n$ as $k\to\infty$.

Proof.

From Lemma 2.9 we observe that $\{X^{k}\}_{k}$ is a bounded sequence in $C^{0}([0,T];{\mathbb{R}}^{m})$ and equicontinuous on $[0,T]$. From Ascoli-Arzela’s Theorem (Theorem 2.6), there exists a uniformly convergent subsequence $X^{k_{l}}$ of $X^{k}$. We define $X$ as its limit as $l\to\infty$, which completes the proof of the first of the two convergence statements in (2.9). Next we fix $i\in\{1,\ldots,n\}$. From the uniform bound on $\max_{0\leq t\leq T}|\dot{x}_{i}^{k}(t)|$, which is guaranteed by Lemma 2.9, there exists a subsequence of $\{X^{k_{l}}(t)\}_{l}$ (not relabeled) and $y_{i}\in L^{\infty}(0,T;{\mathbb{R}}^{m})$ such that $\dot{x}_{i}^{k_{l}}\rightharpoonup y_{i}$ in $L^{\infty}(0,T;{\mathbb{R}}^{m})$ as $l\to\infty$. We characterise $y_{i}$ as $\dot{x}_{i}$ by

[TABLE]

∎

Let $X$ be given by Lemma 2.10. In the remainder of the proof we pass to the limit $k\to\infty$ in the weak form of (2.7) in order to show that $X$ satisfies the second of the three conditions in Proposition 2.3. Let $\psi\in C_{0}^{\infty}((0,T);{\mathbb{R}}^{m})$ be a test function. To simplify notation, we define $F_{i}^{k}(t):=F_{i}(t,X^{k}(t))$ and $F_{i}(t):=F_{i}(t,X(t))$. We fix $i\in\{1,\dots,n\}$ and define the open set $I_{k}:=\{t\in[0,T]~{};~{}x_{i}^{k}(t)\notin\overline{\Omega}\}$. Starting from the weak form of (2.7), we use (2.8) to obtain

[TABLE]

where

[TABLE]

Using Lemma 2.10 it is easy to pass to the limit $k\to\infty$ in the left-hand side of (2.10). Passing to the limit in the right hand side requires more care. We start by proving that

[TABLE]

With this aim, we decompose $I_{k}=\bigcup_{\lambda=1}^{\infty}Q_{\lambda}$, where $Q_{\lambda}$ are open connected components of $I_{k}$. For any $\lambda$, we calculate

[TABLE]

Since $I_{k}=\bigcup_{\lambda=1}^{\infty}Q_{\lambda}\subset(0,T)$, we obtain

[TABLE]

which proves (2.11).

Next we prepare for passing to the limit in $B_{k}$ as $k\to\infty$. We define the indicator functions as

[TABLE]

We further set

[TABLE]

and observe that

[TABLE]

Lemma 2.11**.**

For all $t\in(0,T]$, there holds \chi_{\partial\Omega}(t)\big{(}g(t)\big{)}_{-}=0.

Proof.

We prove the assertion by the argument of contradiction. Suppose that the assertion is false. Then, there exists $t_{0}\in(0,T]$ such that $x_{i}(t_{0})\in\partial\Omega$ and $\alpha:=-g(t_{0})>0$ hold. By the continuity of $g$, there exists $\delta\in(0,t_{0})$ such that $g(t)\leq-\frac{\alpha}{2}$ for $t_{0}-\delta\leq t\leq t_{0}$. We choose $k_{1}\geq k_{0}$ such that $g_{k}(t)\leq-\frac{\alpha}{4}$ holds for all $k\geq k_{1}$ and all $t_{0}-\delta\leq t\leq t_{0}$. Since

[TABLE]

we get

[TABLE]

Passing to the limit $k\to\infty$, we obtain $d_{s}^{\varepsilon}(x_{i}(t_{0}))\leq d_{s}^{\varepsilon}(x_{i}(t_{0}-\delta))-\frac{\alpha\delta}{4}\leq-\frac{\alpha\delta}{4}$, which contradicts $x_{i}(t_{0})\in\partial\Omega$. ∎

Using the preparations above, we pass to the limit in $B_{k}$ as $k\to\infty$ in the next lemma.

Lemma 2.12**.**

There exists $\overline{\chi}\in L^{\infty}(0,T;[0,1])$ such that

[TABLE]

Proof.

We define

[TABLE]

such that we can split $B_{k}$ as

[TABLE]

It is easy to see that there exists $C>0$ such that $|\rho_{k}(t)|\leq C$ for $k\geq k_{0}$ and $t\in[0,T]$.

We claim that

[TABLE]

From this claim and the Dominated Convergence Theorem, we then obtain

[TABLE]

Next we prove the claim in (2.14). We fix an arbitrary $t_{0}\in(0,T)$. By Lemma 2.11 it is enough to consider the following three cases:

case 1: $x_{i}(t_{0})\in\Omega$. There exists $k_{1}\geq k_{0}$ such that for any $k\geq k_{1}$ there holds $x_{i}^{k}(t_{0})\in\Omega$. For any such $k$, we obtain $\chi_{k}(t_{0})=0$ and $\chi_{k}^{*}(t_{0})=0$, and thus $\rho_{k}(t_{0})=0$.

case 2: $x_{i}(t_{0})\in\partial\Omega$ and $g(t_{0})=0$. Since $\lim_{k\to\infty}g_{k}(t_{0})=g(t_{0})=0$, we obtain $\lim_{k\to\infty}\rho_{k}(t_{0})=0$.

case 3: $x_{i}(t_{0})\in\partial\Omega$ and $g(t_{0})>0$. There exists $k_{1}\geq k_{0}$ such that $g_{k}(t_{0})>0$ for all $k\geq k_{1}$. Let such a $k$ be fixed. If $x_{i}^{k}(t_{0})\notin\overline{\Omega}$, then since $\chi_{k}(t_{0})=\chi_{k}^{*}(t_{0})=\chi_{\partial\Omega}(t_{0})=1$, we obtain $\rho_{k}(t_{0})=0$. Otherwise, if $x_{i}^{k}(t_{0})\in\overline{\Omega}$, then $\chi_{k}(t_{0})=\chi_{k}^{*}(t_{0})=0$ and thus $\rho_{k}(t_{0})=0$.

Together, cases 1-3 imply (2.14).

In the final step we pass to the limit in $C_{k}$ as $k\to\infty$. Since $\chi_{k}^{*}(t)\in\{0,1\}$, there exists a subsequence of $\{\chi_{k}^{*}\}_{k}$ (not relabelled) and $\overline{\chi}\in L^{\infty}(0,T;[0,1])$ such that $\chi_{k}^{*}\rightharpoonup\overline{\chi}$ weakly-$*$ in $L^{\infty}(0,T)$ as $k\to\infty$. Then, from (2.12) and (2.13) we obtain

[TABLE]

From this and (2.15), Lemma 2.12 follows. ∎

Using (2.11) and Lemma 2.12, we pass to the limit $k\to\infty$ in (2.10). This yields

[TABLE]

Since $\psi\in C_{0}^{\infty}((0,T);{\mathbb{R}}^{m})$ is arbitrary, this implies that

[TABLE]

In the final step we show that (2.17) is equivalent to the second of the three conditions in Proposition 2.3. We put $A:=\{t\in[0,T]~{};~{}x_{i}(t)\in\partial\Omega\}$. Then, for a.e. $t\in A$, we obtain

[TABLE]

Applying the co-area formula (Theorem 2.7) to $A$ and $d_{s}^{\varepsilon}\circ x_{i}(t)$, we obtain

[TABLE]

In the last equality we have used that the integrand is [math] a.e. (in fact, the integrand is [math] for any $y\neq 0$). From (2.18), we obtain

[TABLE]

From (2.17), using (2.19) and Lemma 2.11, we obtain

[TABLE]

We conclude by applying Proposition 2.3.

2.6 Uniqueness of the solution

In this section, we prove the uniqueness statement of Theorem 2.4. We suppose that $F_{i}$ satisfies the Lipschitz condition (2.3). Let $T>0$ and $X^{\circ}\in G$ be given, and let $X(t)=\{x_{i}(t)\}_{i=1}^{n}$ and $Y(t)=\{y_{i}(t)\}_{i=1}^{n}$ be two mild solutions of Problem 2.1 subject to the initial condition $X^{\circ}$.

We prove uniqueness by applying Gronwall’s lemma. With this aim (recalling that $X,Y$ are Lipschitz continuous by Proposition 2.3) we set

[TABLE]

for $i\in\{1,2,\dots,n\}$ and a.e. $t\in(0,T)$. The following lemma provides a sufficient bound on $J_{i}(t)$.

Lemma 2.13**.**

There exists a constant $C>0$ independent of $X$, $Y$ such that

[TABLE]

Proof.

For a.e. $t\in(0,T)$, we obtain

[TABLE]

We prove the required estimate on $J_{i}(t)$ by separating three cases. Since the computations below do not depend on time or on $i$, we simplify notation by setting $x=x_{i}(t)$, $y=y_{i}(t)$, $X=X(t)$, $Y=Y(t)$, $f(X)=F_{i}(t,X(t))$ and $J=J_{i}$.

Case 1: $x\in\partial\Omega,~{}f(X)\cdot\nu(x)>0$ and $y\in\partial\Omega,~{}f(Y)\cdot\nu(y)>0$. Since $\nu$ and $f$ are Lipschitz continuous, we obtain

[TABLE]

Case 2: Both $x$ and $y$ do not satisfy the condition in Case 1. We obtain

[TABLE]

Case 3: $y\in\partial\Omega$, $f(Y)\cdot\nu(y)>0$ and $x$ does not satisfy the condition in Case 1. We compute

[TABLE]

To complete the estimate, we separate a further set of three cases.

Case 3.1: $|x-y|\geq\varepsilon$. We obtain

[TABLE]

where $C_{\varepsilon}=\frac{2\|F\|_{\infty}}{\varepsilon}$.

Case 3.2: $|x-y|<\varepsilon$ and $(x-y)\cdot\nu(y)\leq 0$ (see Figure 5 for a sketch). We obtain

[TABLE]

Case 3.3: $|x-y|<\varepsilon$ and $(x-y)\cdot\nu(y)>0$ (see Figure 5 for a sketch). We estimate

[TABLE]

Since $0<f(Y)\cdot\nu(y)<\|F\|_{\infty}$, it is enough to show that

[TABLE]

for some $C>0$ which is independent of $x$ and $y$. With this aim, we set

[TABLE]

where $\bar{x}$ is uniquely determined since $x\in N^{\varepsilon}(\partial\Omega)\cap\overline{\Omega}$. We notet that

[TABLE]

We expand the left-hand side of (2.20) as

[TABLE]

Next we show that each of the three terms in the right hand side of (2.22) is bounded from above by $C\|X-Y\|^{2}$. For the first term, we find $(x-\bar{x})\cdot\nu(\bar{x})=-|x-\bar{x}|\leq 0$. For the third term, we estimate

[TABLE]

For the second term in (2.22), we obtain by using a Taylor expansion of $d_{s}$ that

[TABLE]

where $|R(\bar{x},y)|\leq C|\bar{x}-y|^{2}$. Hence,

[TABLE]

This completes the proof of (2.20).

∎

To complete the proof of the uniqueness statement of Theorem 2.4, we continue from Lemma 2.13 by using Gronwall’s Lemma. Integrating the result from Lemma 2.13 from [math] to $t$, we obtain for any $1\leq i\leq n$ that

[TABLE]

Since the constant $C$ does not depend on $i$, we can take the maximum over $i$ on both sides to obtain

[TABLE]

Hence, by Gronwall’s Lemma,

[TABLE]

This completes the proof of the uniqueness statement of Theorem 2.4.

2.7 Properties of the solutions

We prove two properties of solutions to Problem 2.1. As we will demonstrate by simulations in Section 4, particles can attach to or detach from the boundary. The first property states that particles can detach only in tangential direction. Again, we will simply write $F_{i}(t)=F_{i}(t,X(t))$.

Proposition 2.14**.**

Suppose that (2.2) holds and $X(t)~{}(0\leq t\leq T)$ is a solution of Problem 2.1. If $x_{i}(t_{0})\in\partial\Omega$ at $t_{0}\in(0,T]$, then it holds that

[TABLE]

Moreover, if there exists a $\delta$ such that $x_{i}(t)\in\Omega$ for all $t_{0}<t<t_{0}+\delta\leq T$, then it holds that

[TABLE]

Proof.

As in (2.13), we set g(t):=\big{(}\nabla d_{s}^{\varepsilon}(x_{i}(t)))\cdot F_{i}(t)\big{)}. We prove $g(t_{0})=F_{i}(t_{0})\cdot\nu(x_{i}(t_{0}))\geq 0$ by contradiction. Suppose that $g(t_{0})<0$. Then there exists an $h>0$ such that $g(t)<0$ for all $t\in[t_{0}-h,t_{0}+h]$. From the definition of $H_{i}$, we then obtain that $\dot{x}_{i}(t)=F_{i}(t)$ for a.e. $t\in[t_{0}-h,t_{0}+h]$. Hence, $x_{i}$ is differentiable at $t_{0}$, and we obtain

[TABLE]

which contradicts with $x_{i}(t)\in\overline{\Omega}$. Therefore, $F_{i}(t_{0})\cdot\nu(x_{i}(t_{0}))\geq 0$.

Next we prove the second statement of Proposition 2.14. Let the asserted $\delta$ be given. As before, we obtain that $x_{i}\in C^{1}([t_{0},t_{0}+\delta];{\mathbb{R}}^{m})$, and thus

[TABLE]

and

[TABLE]

Combining this with (2.23), we also obtain $F_{i}(t_{0})\cdot\nu(x_{i}(t_{0}))=0$. ∎

The second property specifies the sense in which, under the Lipschitz condition (2.3) on $F$, how the solution to Problem 2.8 approximates the solution to Problem 2.1.

Theorem 2.15** (Approximation theorem).**

If $F_{i}$ satisfies (2.2) and the Lipschitz condition (2.3), the solution $X^{k}(t)$ of Problem 2.8 is unique for each $k>\frac{1}{\varepsilon}\|F\|_{\infty}$ and it uniformly converges to the unique solution $X(t)$ of Problem 2.1 on $[0,T]$ as $k\to\infty$.

Proof.

The uniqueness of $X^{k}(t)$ follows from the Lipschitz continuity of $F_{i}$ and $d\nabla d$ (Proposition 2.5). By the argument in the proof of the existence of $X$, any uniformly converging subsequence of $\{X^{k}\}_{k}$ converges to a solution $Y$ of Problem 2.1. Since the solution of Problem 2.1 is unique, we conclude that the whole sequence $X^{k}$ converges uniformly to $X$ on $[0,T]$ as $k\to\infty$. ∎

3 Gradient flows of particle interaction energies

3.1 Particle interaction energy and its minimization problem

Let $\Omega$ and $X=\{x_{i}\}_{i=1}^{n}\in G$ be as in Section 2.1. To each state $X\in G$ we assign an energy value which is the sum of an interaction potential $V(x_{i}-x_{j})$ (summed over each pair of two particles $x_{i}$ and $x_{j}$) and a potential energy $W(x_{i})$ (summed over $i$). We assume that the potentials $V$ and $W$ satisfy

[TABLE]

Figure 6 illustrates a typical example of $V$.

We set

[TABLE]

as the set of states where the interaction energy is finite. The, we define the particle energy $E:G_{*}\to{\mathbb{R}}$ as

[TABLE]

We consider the following minimization problem.

Problem 3.1**.**

Find $\overline{X}\in G_{*}$ such that

[TABLE]

Existence of solutions to Problem 3.1 follows from the lower semi-continuity of $E$ and $G$ being closed and bounded. In preparation for constructing a gradient flow of $E$, we prove that the particles within a particle configuration $X$ with finite energy have a minimal length of separation.

Lemma 3.2**.**

For all $X\in G_{*}$, there exists an $h=h(E(X))>0$ non-increasing with respect to the value of $E(X)$, such that

[TABLE]

Proof.

Since $V,W\geq 0$, we obtain

[TABLE]

Since $V(x)\to\infty$ as $|x|\to 0$, there exists an $h>0$ such that for all $|x|\leq h$, $V(x)>(n-1)nE(X)$. Together with (3.2), this also proves (3.2). ∎

3.2 Gradient flow of the interaction energy

We consider $n$ moving particles $X(t)\in G_{*}$. For a particle $x_{i}=(x_{i}^{1},\dots,x_{i}^{m})^{\top}$, we set

[TABLE]

and

[TABLE]

In this setting, the problem corresponding to Problem (2.1) is:

Problem 3.3**.**

Let $P$ be as in (2.4), $F_{i}$ as in (3.5), $X^{\circ}=\{x_{i}^{\circ}\}_{i=1}^{n}\in G_{*}$ and $T>0$. Find a solution to

[TABLE]

for $i=1,2,\dots,n$.

Similar to Problem 2.1, we do not expect solutions to Problem 3.3 to be of class $C^{1}$. Therefore, we seek again mild solutions (see Definition 2.2) to Problem 3.3. Our main result in this section (Theorem 3.4) is the existence and uniqueness of the solution to Problem 3.3.

Theorem 3.4**.**

Let $E$ be as in (3.1). Then for all $T>0$ and all $X^{\circ}\in G_{*}$, there exists a unique mild solution (see Definition 2.2) to Problem 3.3.

3.3 Proof of Theorem 3.4: existence and uniqueness

In this section we prove Theorem 3.4. Let $T>0$ and $X^{\circ}\in G^{*}$ be given. Since $V$ is singular at 0, $F_{i}$ does not satisfy the Lipschitz condition in Section 2.1. Hence, we cannot apply directly Theorem 2.4 to prove Theorem 3.4.

We construct a solution to Problem 3.3 by introducing an auxiliary problem which fits to the assumptions in Section 2.1 and whose solution will turn out to be a solution to Problem 3.3 too. With this aim, we first prove an energy-decay estimate under the assumption that Problem 3.3 attains a solution.

Proposition 3.5**.**

Let $X$ be a solution to Problem 3.3. Then for a.e. $t\in(0,T)$, $\frac{d}{dt}E(X(t))\leq 0$.

Proof.

Setting

[TABLE]

we rewrite (3.6) as

[TABLE]

for a.e. $t\in[0,T]$, $i=1,2,\dots,n$. Hence, we obtain

[TABLE]

for a.e. $t\in[0,T]$. ∎

Using Lemma 3.2, we obtain an $h=h(E(X^{\circ}))>0$ such that

[TABLE]

Using $h$, we choose a regular potential $\widetilde{V}\in C^{2}({\mathbb{R}}^{m};{\mathbb{R}})$ which satisfies

[TABLE]

Figure 7 illustrates a typical example of $\widetilde{V}$.

Using $\widetilde{V}$, we define the regularized particle energy $\widetilde{E}:G\to{\mathbb{R}}$ by

[TABLE]

By (3.7), we note that $\widetilde{E}=E$ on

[TABLE]

Using the preparations above, we finally introduce the auxiliary problem.

Problem 3.6**.**

Let $P$ be as in (2.4), $\widetilde{F}_{i}:=-\nabla_{i}\widetilde{E}$, $\widetilde{X}^{\circ}\in G$ and $\widetilde{T}>0$. Find a solution to

[TABLE]

for $i=1,2,\dots,n$.

If follows directly from Theorem 2.4 that Problem 3.6 attains a unique solution for any $\widetilde{X}^{\circ}\in G$ and $\widetilde{T}>0$. Choosing $\widetilde{X}^{\circ}=X^{\circ}$ and $\widetilde{T}=T$, we denote the corresponding solution to Problem 3.6 by $\widetilde{X}$. It is obvious from the proof of Proposition 3.5 that $\widetilde{E}$ and $\widetilde{X}$ also satisfy $\frac{d}{dt}\widetilde{E}(\widetilde{X}(t))\leq 0$ for all $t\in(0,T)$. Following a similar argument as in the proof of Lemma 3.2, we then obtain

[TABLE]

for all $t\in[0,T]$. From (3.3) and (3.7), we obtain for any $|x|\leq h$ that

[TABLE]

From the two inequalities above we conclude that $\widetilde{X}(t)\in G_{h}$ for all $t\in[0,T]$. Therefore, since $\widetilde{E}=E$ on $G_{h}$, $\widetilde{X}$ is also a solution to Problem 3.3.

Next we prove the uniqueness of solutions to Problem 3.3. If $X$ is a solution to Problem 3.3 with initial condition $X^{\circ}$, we obtain from Proposition 3.5 that $E(X(t))$ is non-increasing in time. Hence, Lemma 3.2 provides an $h=h(E(X^{\circ}))>0$ such that $X(t)\in G_{h}$ for all $t\in[0,T]$. Setting again $\widetilde{V}\in C^{2}({\mathbb{R}}^{m};{\mathbb{R}})$ and $\widetilde{E}$ as above, $X$ satisfies Problem 3.6. Since the solution to Problem 3.6 is unique, Problem 3.3 cannot have any other solution.

4 Numerical examples

In this section we compute and discuss numerical solutions to the gradient flow with confinement as given by Problem 3.3. We solve Problem 3.3 numerically by discretizing in time with a fixed time step $\Delta t$. We compute the solution for the next time step by first applying the classical fourth order Runge-Kutta method without considering the confinement of the particles to the domain $\Omega$, and then projecting all particles outside $\Omega$ along the normal $\nu$ to $\partial\Omega$.

We solve Problem 3.3 for several different scenarios. In all scenarios, we choose the interaction potential as

[TABLE]

We split the scenarios in three sections for different choices of the domain $\Omega$.

4.1 Circular domain $\Omega$

We set $\Omega:=B_{1}(0):=\{|x|<1\}$ and consider no external force (i.e., $W\equiv 0$). We further take $n=3000$ and $\Delta t=0.5$ as reference values. Regarding the initial positions $X^{\circ}$ for $n$ particles, we consider two cases:

• Case 1

Similar to $G_{*}$ in Section 3, we set

[TABLE]

where we interpret $0.025$ as the separation distance. We choose $X^{\circ}\in G_{1}$ randomly uniformly. Figure 8 illustrates the realisation of $X^{\circ}$ that we use in our simulation.

• Case 2

Similar to Case 1, we choose $X^{\circ}$ from the uniform distribution on

[TABLE]

Figure 9 illustrates an example of $X^{\circ}$.

Figure 8 illustrates the particle positions $X(t)$ at $t=3000$. Firstly, we observe that the density of the particles is higher near the boundary. This is in line with the potential $V$ being decreasing in radial direction. Also, from Proposition 2.14 we know that the particles cannot detach from the boundary after hitting it111To see this, suppose $x_{i}(t)\in\partial\Omega$. Then for all $j\neq i$, the force that particle $j$ exerts on particle $i$ has direction $x_{i}(t)-x_{j}(t)$. Together with the strict convexity of $\Omega$, we then conclude that $F_{i}(t,X(t))\cdot\nu(x_{i}(t))>0$.. Secondly, the color-coding suggests that no mixing of particles occurs.

Next we solve Problem 3.3 for the initial condition as in Case 2. Figure 9 illustrates the numerical solution $X$ at several time instances. As in Case 1, we observe no mixing between the red and the blue particles. Also, the particle positions at $t=3000$ are comparable. A new effect which we observe (especially in the second figure) is that several blue particles are separated from the bulk of the blue particles. A possible explanation for this effect is that we use an explicit time discretisation which is unstable if $\Delta t$ is not small enough.

To gain insight in the dependence of the numerical solution, we consider the energy value $E(X(t))$ along the solutions in Case 1 and Case 2 for different values of $\Delta t$ and $n$. The results are displayed in Figure 10. We first discuss the case where $\Delta t=3$. The graph displays peaks which violate the property of the analytical solution to Problem 3.3 on the decay of the energy (see Proposition 3.5). A possible explanation for these peaks is the (unstable) explicit time-discretization, which may cause the interaction energy to become large when particles come too close to each other. Interestingly, this instability occurs throughout the complete time interval.

Next we examine the figure where five different values of $n$ are considered. We observe for $t$ large enough that the energy increases as $n$ increases, but that the rate of the gap between these energy values decreases as $n$ increases. Not only does this observation suggest that solutions to Problem 3.3 (after rescaling time by a factor $n$) may converge as $n\to\infty$ to a dynamical particle density, it also hints to a convergence rate.

4.2 A non-convex domain $\Omega$ (dumbbell)

We choose $\Omega$ to be the non-convex domain as illustrated in Figure 11. $\Omega$ is chosen such that the initial positions $x_{i}^{\circ}$ of Case 2 of Section 4.1 fit inside $\Omega$ and such that $\Omega$ is included in the domain considered in Section 4.1. Hence, it is natural to compare this setting to Case 2 of Section 4.1. Other than $\Omega$, we will therefore make as few adjustments to the setting of Case 2 as possible. We only change $\Delta t$ to $\Delta t=0.25$ to prevent numerical instabilities.

Figure 11 illustrates the solution to Problem 3.3 at different time instances, and Figure 12 displays the corresponding decay of the energy. We observe that the particles, despite the non-convexity, spread out through the dumbbell in time such that the amount of particles on the left side is comparable to that on the right side. As in Section 4.1 we observe no mixing between the red part and blue particles. Regarding the energy decay, we observe that our confinement of the domain with respect to Case 2 in Section 4.1 yields a slower decrease of the energy.

4.3 Channels with obstacles

In this section we set $\Omega$ as a channel and consider 2 different obstacles in it. The channel is given by $(-\infty,\infty)\times(0,1.2)$ and the obstacles are displayed in Figures 14 and 15. We set $W(x)=-0.002x_{1}$ as a driving force to move the particles past the obstacle. We further put $\Delta t=0.5$ and consider $n=900$ particles, which are initially put as a regular square grid in the range of $-1.7\leq x_{1}\leq 3.1$. We are interested in the manner in which particles may detach from the boundary after passing the obstacle, and to which extend the obstacle slows down the flow of the particles.

As a reference, we first solve Problem 3.3 on the channel without obstacle. Figure 13 illustrates the dynamics.

Next we insert a small obstacle, given by a bump, in the channel. Figure 14 illustrates the channel and the dynamics. We observe that the particles starting at the lower part of the channel pass the barrier more slowly than those that start in the upper part. The plots also show how some particles detach from the boundary after the obstacle. We also observe that the flow of the particles is hardly slowed down when compared to the dynamics displayed in Figure 13.

Next we insert a large obstacle in the channel with a shape similar to a horse shoe. Figure 15 illustrates the dynamics. The detachment of the particles from the boundary of the horse shoe is clearly visible. We also observe that several particles get stuck in the left lower region around the obstacle. Also, the total flow of the particles through the channel per unit time is significantly less than the previous two cases. From the color coding we see that it is not so clear in which order the particles pass the obstacle. Also, after having passed the obstacle, the separation line between the red and blue particles becomes less clear, which indicates a small mixing effect.

Figure 16 compares the energy decay of the dynamics between all three channels. While the small obstacle has little influence on the decay of the energy, the large obstacle clearly slows down the decay until most particles have passed the obstacle.

5 Conclusion and outlook

In Theorem 2.4 we proved existence and uniqueness of mild solutions to interacting particle systems when the particles are confined to a given domain (Problem 2.1). We also showed that when particles detach from the boundary, they do so with velocity tangential to the boundary (Proposition 2.14). We applied the results to gradient flows with confinement (see Theorem 3.4). We have illustrated the dynamics by numerical simulations for various (non-convex) domains with a special focus on particles attaching to and detaching from the boundary.

The numerical illustrations spark several topics for future research. For instance, they suggest that the collective motion of the particles could be described by a continuum density. In future work we intend to obtain this density from Theorem 2.4 by passing to the many-particle limit $n\to\infty$ in Problem 2.1. Another example is that the coloring of the particles suggests that no mixing occurs (at least for convex domains); it would be interesting to formulate and prove a precise statement on the mixing of particles. Lastly, the group behavior of the particles shows a moving front in time for the support of the particle density; it would be interesting to describe this front with more precision.

Acknowledgments

This work is partially supported by JSPS KAKENHI Grant Numbers JP16H02155 and JP17H02857.