Coalescing particle systems and applications to nonlinear Fokker-Planck equations

TL;DR

This paper introduces a stochastic particle system with singular interactions, relating it to nonlinear Fokker-Planck equations, and develops an efficient numerical method to simulate collision dynamics and singularity formation in models like Keller-Segel.

Contribution

It presents a novel particle system model with inelastic collisions and a collision-detection numerical method that avoids pairwise computations, linking microscopic dynamics to nonlinear PDEs.

Findings

The numerical method effectively detects collisions without pairwise calculations.

Hydrodynamic limit of the particle system solves nonlinear Fokker-Planck equations.

Numerical simulations capture finite-time singularities and blow-up behavior in Keller-Segel models.

Abstract

We study a stochastic particle system with a logarithmically-singular inter-particle interaction potential which allows for inelastic particle collisions. We relate the squared Bessel process to the evolution of localized clusters of particles, and develop a numerical method capable of detecting collisions of many point particles without the use of pairwise computations, or very refined adaptive timestepping. We show that when the system is in an appropriate parameter regime, the hydrodynamic limit of the empirical mass density of the system is a solution to a nonlinear Fokker-Planck equation, such as the Patlak-Keller-Segel (PKS) model, or its multispecies variant. We then show that the presented numerical method is well-suited for the simulation of the formation of finite-time singularities in the PKS, as well as PKS pre- and post-blow-up dynamics. Additionally, we present numerical evidence that blow-up with an increasing total second moment in the two species Keller-Segel system occurs with a linearly increasing second moment in one component, and a linearly decreasing second moment in the other component.

Equations

$$\{\begin{array}{ll}{\partial }_t{\rho }_1& =\nabla \cdot ({\mu }_1\nabla {\rho }_1-\chi {\rho }_1\nabla c),\\ & ⋮\\ {\partial }_t{\rho }_K& =\nabla \cdot ({\mu }_K\nabla {\rho }_K-\chi {\rho }_K\nabla c),\\ \mathcal{L}c& =-({\rho }_1+\cdots +{\rho }_K),\end{array}$$

$$\mathcal{L}c(x,t)=\nabla \cdot \left(G(x)\nabla c(x,t)\right)+F(x,c)$$

$$d{X}_t^{(n)}=-\chi \frac{\partial }{\partial X_t^{(n)}}\sum _{\begin{array}{c}i=1\\ i\ne n\end{array}}^N{m}_{i}V(X_t^{(n)},X_t^{(i)})dt+\sqrt{\frac{2\tilde{\mu }}{m_n}}d{W}_t^{(n)}.$$

$$V(x,y)$$

$$\{\begin{array}{ll}{\partial }_t\rho & =\nabla \cdot (\mu \nabla \rho -\chi \rho \nabla c),\\ \Delta c& =-\rho ,\end{array}$$

$$M>8\pi \mu /\chi ,$$

$$M<8\pi \mu /\chi ,$$

$$T<\frac{2\pi F(0)}{(\chi M-8\pi \mu )M},$$

$$F(t)={\int }_{{\mathbb{R}}^2}|x{|}^2\rho (x,t)dx$$

$$d{X}_t^{(n)}=-\frac{\chi }{2\pi }\sum _{i\ne n}{m}_i\frac{X_t^{(n)}-X_t^{(i)}}{{\left|X_t^{(n)}-X_t^{(i)}\right|}^2}dt+\sqrt{\frac{2\tilde{\mu }}{m_n}}d{W}_t^{(n)}.$$

$$Y_t=\frac{1}{2M^2}\sum _{i,j}{m}_i{m}_j{\left|X_t^{(i)}-X_t^{(j)}\right|}^2.$$

$$d{Y}_t=\alpha dt+2\beta \sqrt{Y_t}d{W}_t$$

$$\alpha$$

$$d{W}_t=\frac{1}{(M{)}^{3/2}\sqrt{Y_t}}\sum _{i,j=1}^N{m}_j\sqrt{m_i}\left(X_t^{(i)}-X_t^{(j)}\right)\cdot d{W}_t^{(i)}$$

$$d{\tilde{Y}}_t=2(\nu +1)dt+2\sqrt{{\tilde{Y}}_t}d{W}_t,$$

$$\nu (m_1,m_2,\cdots ,m_N)=(N-2)-\frac{\chi M^2}{8\pi \tilde{\mu }}\left(1-\sum _j{\left(\frac{m_j}{M}\right)}^2\right).$$

$$\nu (m_1,\dots ,m_N)<0.$$

$${\nu }_f-{\nu }_i=-\left(\overline{\nu }+1\right).$$

$$\tau \asymp \tilde{\mu }\frac{Y_0}{U},$$

$$d{x}_{\tau -ϵ}=-\chi \sum _{i=1}^{N^{\prime }}{m}_i\frac{\partial V}{\partial x}\left(x_{\tau -ϵ},X_{\tau -ϵ}^{(i)}\right)dt+\sigma (m_{N^{\prime }+1})d{W}_{\tau -ϵ}^{(N^{\prime }+1)}.$$

$$d{x}_{\tau }=-\chi M^{\prime }\frac{\partial V}{\partial x}\left(x_{\tau },Z_{\tau }\right)dt+\sigma (m_{N^{\prime }+1})d{W}_{\tau }^{(N^{\prime }+1)}.$$

$$Z_{\tau }=\underset{ϵ\to 0^{+}}{\lim }\frac{1}{M^{\prime }}\sum _{i=1}^{N^{\prime }}{m}_i{X}_{\tau -ϵ}^{(i)},$$

$$\sigma (M^{\prime })=\frac{1}{M^{\prime }}\sqrt{\sum _{i=1}^{N^{\prime }}{m}_i^2(\sigma (m_i){)}^2},$$

$$(M^{\prime }{)}^2\left(\sigma (M^{\prime }\right){)}^2=\sum _{i=1}^{N^{\prime }}{m}_i^2(\sigma (m_i){)}^2.$$

$$f(x)=f(x)+f(y).$$

$$\sigma (m)=\sqrt{\frac{2\tilde{\mu }}{m}},$$

$$W_t^{(cm)}=\frac{1}{\sqrt{M^{\prime }}}\sum _{i=1}^{N^{\prime }}\sqrt{m_i}{W}_t^{(i)}.$$

$$d{Z}_t=-\chi m\frac{\partial V}{\partial x}\left(Z_t,x_t\right)dt+\sqrt{\frac{2\mu }{M^{\prime }}}d{W}_t^{(cm)},$$

$$\widetilde{Y_t}=\frac{1}{2(M^{\prime }{)}^2}\sum _{i,j=1}^{N^{\prime }}{m}_i{m}_j{\left|X_t^{(i)}-X_t^{(j)}\right|}^2,$$

$$d\widetilde{Y_t}\approx d{Q}_t=\left(\frac{4\tilde{\mu }(N^{\prime }-1)}{M^{\prime }}-\frac{\chi M^{\prime }}{2\pi }\left(1-\sum _{j=1}^{N^{\prime }}{\left(\frac{m_j}{M^{\prime }}\right)}^2\right)\right)dt+2\sqrt{\widetilde{Y_t}}\sqrt{\frac{2\tilde{\mu }}{M^{\prime }}}d{\tilde{W}}_t,$$

Full text

{AMS}

35K58, 35Q83, 35Q92, 45G05, 60H30, 60H35, 65C35, 82C21, 82C22, 82C31, 82C80, 92C17

1 Introduction

1.1 Background

The connection between systems of interacting particles and kinetic-type PDEs was first investigated by Kac in his study of the motion of a tagged molecule in a bath of identical molecules [18], which arose as a simplified model of a Maxwellian gas [24]. This work introduced the property of “propagation of chaos”: as the number of molecules tends to infinity, the $N$-particle probability densities are well-approximated by the product of single particle marginals.

The connection between such processes and nonlinear parabolic equations, such as Boltzmann’s equation or Burgers’ equation, was then elaborated by McKean [25]. This line of research has continued since, and much more is now known about the duality between these processes and parabolic PDEs [28]. In particular, particle-based numerical methods have been developed for the solution of such PDEs [3] using the methods of “mean field Monte Carlo.” The solutions to these PDEs are approximated by the empirical density of $N$-particle systems. As the number of particles tends to infinity, such approximations become exact by the propagation of chaos property.

Rigorously proving propagation of chaos for particle systems with singular interaction coefficients is challenging, and has only been carried out in a few special cases, e.g. [17]. One PDE associated with a logarithmically-singular particle system is the Patlak-Keller-Segel chemotaxis model (PKS) [20, 26], which is reviewed extensively in [15, 16]. Despite the lack of a propagation of chaos result, the PDE has been numerically approximated using the associated particle system in several works, initially in [12, 13] and later in [10]. Various properties of the PKS, such as the formation of Dirac singularities in finite time [1], as well as interaction of singularities post-blow-up [30, 31, 7], can either be shown to be true in the particle system, or have considerable numerical evidence for their existence. Recent advacements in understanding this particle system include partial existence and uniqueness results for solutions to the subcritical (small mass) particle system [11, 4], and convergence of the empirical density of a similar particle system to the solution of a modified PKS system [2].

Singular interaction coefficients in the PKS particle system allow for particle collisions, and some type of regularization must be introduced in order to propagate the particle system past the first collision time. In [12], semi-deterministic heavy particles absorb light particles. In [10], collided particles are forced to move in unison due to a mean field. Broadly speaking, the two works take two different approaches to simulating the regularizations of the PKS derived in [7] and [30, 31]. The first work simulates the singular limit of the system, whereas the second work simulates the system with an effectively regularized Green’s function.

In [12], heavy particles corresponding to singularities in the PDE must be prescribed a priori and cannot arise as the result of a collision of many light particles. On the other hand, particles do not truly collide in [10], and the deterministic system approximated is closer to the one given in [30, 31], where singularities are replaced with regions of high density. In this work, we develop criteria for particle coalescence of particles of arbitrary masses, based on analytical estimates of exit times of the squared Bessel process. In this context, the particle system in [12] can be viewed as the limit of the particle system in [10] with collisions, as the number of particles tends to infinity.

1.2 Outline

We introduce a coalescing particle system with nonuniform particle masses and a logarithmic interaction kernel. Using estimates on the system’s second moment, we derive a criterion for a finite-time collision of the entire particle system. We then motivate the mass-dependence of the diffusion coefficient of a particle, and approximate the time evolution of a localized subsystem’s second moment. We then show that the hydrodynamic limit of such a system is the multispecies Patlak-Keller-Segel system, of which the PKS is a special case. Finally, we present a numerical method implementing many-body collisions and coalescence events, which is generally applicable to PDEs of the form

[TABLE]

where

[TABLE]

is an elliptic operator with a fundamental solution $V$ which has a logarithmic singularity. As an application, we apply it to the planar case with decaying (radiative) boundary conditions and $\mathcal{L}=\Delta$, though the method is equally applicable to bounded domains with Neumann boundary conditions. This special case is the planar PKS system, some properties of which we describe in Section 1.4, and whose measure-valued solutions we describe in Section 4.2. We also apply the numerical method to investigate blow-up in the components of the multispecies PKS.

1.3 The coalescing particle system

We study the $N$-particle systems described by the following equations

[TABLE]

Each particle has some mass $m_{n}$ and position $X^{(n)}_{t}\in\mathbb{R}^{2}$. The total mass is $M=\sum_{i}m_{i}$, and $\chi,\tilde{\mu}>0$ are parameters. The processes $W^{(n)}_{t}$ are independent Wiener processes.

The particle system in (1.3) is related to the PDE in (1.1) when $V$ is the fundamental solution of $\mathcal{L}$, e.g. if $\mathcal{L}=\Delta$ or $\mathcal{L}=\Delta-k^{2}$, we have

[TABLE]

where $K_{0}$ is the modified Bessel function of the second kind. When $m_{n}=M/N$ and $\tilde{\mu}=\mu M/N$, the empirical mass density of the particle system with (1.4a) approximates the PKS, and the particle system with (1.4b) is the one given in [10].

The dynamics prescribed in equation (1.3) allows for particle collisions provided that $V$ has logarithmic or stronger singularities. In this case, the SDE must be augmented with proper boundary conditions prescribing behavior when at least two particles’ coordinates are identical. Well-posedness and uniqueness results for these types of SDEs have not been rigorously established. We proceed formally, considering inelastic collisions: colliding particles merge into a single particle which absorbs their total mass.

1.4 Properties of the Patlak-Keller-Segel system

Since many of the applications of this work are related to the PKS, we give a short overview of its definition and properties here.

The PKS is prescribed by the following system of PDEs:

[TABLE]

and models a biological system consisting of amoeba, which spread across the plane with mass density $\rho(x,t)$ and produce a chemical (“chemoattractant”) of concentration $c(x,t)$. On average, amoeba diffuse in space with diffusivity $\mu$ and drift in the direction of $\nabla c$ with speed $\chi|\nabla c|$. The chemoattractant diffuses instantly. The boundary condition $\rho(x,t)\rightarrow 0$ as $|x|\rightarrow\infty$ is enforced, and mass is conserved: $\int\rho(x,t)dx=M$.

This system has been investigated extensively in the literature [15, 16], often in connection with the property that when

[TABLE]

solutions form singularities in finite time, and when

[TABLE]

solutions are global in time [1]. In the former case, an upper bound for the singularity formation time $T$ may be given as

[TABLE]

where

[TABLE]

and $F(0)$ is the system’s initial second moment [1].

2 Collisions and post-collision dynamics

2.1 Overview

Let us first carry out a moment-based computation for finding a criterion which predicts whether a particle system will coalesce into a single particle in finite time. Similar to the PKS mass criterion, this criterion only depends on the total mass of the system and the number of particles, and is otherwise independent of the distribution of particles in the plane. We then motivate the mass dependence of the diffusion coefficient of the newly created particle. Finally, we derive an approximate equation for the dynamics of the second moment of an isolated cluster of particles.

2.2 Collision criterion for the full system

Consider an $N$-particle system with masses and $V$ given as in (1.4a). The dynamics of the $n$th particle are then prescribed by

[TABLE]

To quantify the size of the system, consider its second moment

[TABLE]

By the positivity of $Y_{t}$, showing the total collision of the particles in finite time is equivalent to showing that $Y_{T}=0$ for some $T<\infty$.

It can be shown (by an application of Ito’s lemma) that

[TABLE]

where

[TABLE]

and

[TABLE]

is a Wiener process by the Lévy characterization. We stress that expression (2.12) is only valid between collision events, as $\alpha$ depends on the total number of particles and their masses, and must therefore be updated after each collision. Rescaling time as $t\rightarrow t/\beta^{2}$ and setting $\tilde{Y}_{t}=Y_{\beta^{2}t}$, we get

[TABLE]

where $\nu=\frac{\alpha}{2\beta^{2}}-1$. In terms of our original constants, $\nu$ is given by

[TABLE]

Equation (2.15) describes a squared Bessel process with index $\nu$. Its boundary behavior at $\tilde{Y}=0$ depends on its index [27, 19]:

When $\nu\in[0,+\infty)$, the origin is an entrance boundary, and $\tilde{Y}_{t}>0$ a.s. for all $t>0$ 2. 2.

When $\nu\in(-1,0)$, the origin is a regular boundary, and the behavior of the process at this point must be defined (e.g. absorbing boundary, reflective boundary) 3. 3.

When $\nu\in(-\infty,-1]$, the origin is an absorbing boundary which is hit in finite time

It then follows that a full, simultaneous collision of all the particles may occur if

[TABLE]

When $\nu\in(-1,0)$, we may choose the collision, which we call “soft,” to be fully inelastic, or fully elastic. Similarly, when $\nu~{}\in~{}(-\infty,-1]$, only an inelastic collision may occur.

The above is not sufficient for describing all collisions in the system. For instance, we expect the associated singular forces to force the subsystem inside the dashed line in Figure 1 to inelastically collide earlier than the full system. We will approximate the evolution of the second moment of such a colliding subsystem in Section 2.4, but already note here that a localized colliding subsystem’s second moment may be approximated as a separate squared Bessel process that’s independent of the particles not participating in the collision. As shown in Appendix A, the indices of the squared Bessel processes corresponding of the full system pre- and post-collision, and the index of the colliding subsystem, are related via a subtraction formula: if $\nu_{i}$ is the index of the full system described in Figure 1, $\nu_{f}$ is the index of the same system after the particles inside the dashed lines coalesce, and $\overline{\nu}$ is the index of the subsystem inside the dashed line, then

[TABLE]

From (2.18) we see that hard collisions, except in the critical $\overline{\nu}=-1$ case, always increase the system’s overall index, and soft collisions increase the system’s overall index.

To see the effect of this index change on the full system, let $\tau$ be the first hitting time of the origin for the SDE given in (2.12). This hitting time has the inverse gamma distribution [23],

[TABLE]

where $U\sim G(|\nu|,1)$ is distributed according to the gamma distribution with shape parameter $|\nu|$ and rate parameter $1$.

Intuitively, we see that increasing the index implies that a system contracts at a slower rate, and that a system with only hard inelastic collisions contracts at a slower rate after each collision (e.g. as in Figure 6). Furthermore, we expect many systems which can experience soft inelastic collisions to behave similarly, as a localized subsystem with an index $\overline{\nu}\in(-1,0)$ has a low probability of undergoing a collision in a time step (e.g. $\tau$ only has an expected value when $\nu<-1$), and may attract a sufficient number of additional particles into its aggregate to force the aggregate to experience a hard collision instead. Since in this work we will primarily focus on the large particle case, we prescribe that all collisions—soft and hard (i.e. $\nu<0$)—are inelastic.

We remark that the formula for the time derivative of the second moment of the PKS also only gives an upper bound for the formation of a singularity, since for a system of total mass $M$ greater than the system’s critical mass $M_{c}$, a second moment equal to zero implies the formation of a singularity of total mass $M>M_{c}$. However, singularities in the radially-symmetric PKS form with a mass of exactly $M_{c}$ [14, 29], after which the time derivative of the second moment changes [30].

2.3 Post-collision dynamics

The dynamics of the coalescing diffusion system, given by (1.3), are undefined at times when there exist two indices $i$ and $j$ such that $X^{(i)}_{t}=X^{(j)}_{t}$. If we prescribe that collisions only occur inelastically, we can propagate the system past collision times by coarsening the system: that is, by replacing each collided aggregate of particles with a single particle of the same mass as the aggregate. Let us now show the diffusion coefficient of the newly-created particle is inversely-dependent on the square root of the mass, as given in (1.3).

Consider an $N^{\prime}+1$ particle system, with the first $N^{\prime}$ particles positioned in a tight, pre-coalesced cluster at $X^{(n)}_{t}$ with masses $m_{n}$ totalling to $M^{\prime}$, and the last particle located far away at $x=X^{(N^{\prime}+1)}_{t}$ with mass $m=m_{N^{\prime}+1}$, as in Figure 2(a). In general, the diffusion coefficient of a particle may be given as a function of the particles mass, $\sigma_{n}=\sigma(m_{n})$. Let $\tau$ denote the time at which the first $N^{\prime}$ particles coalesce at $Z_{\tau}$, and fix $0<\epsilon\ll\tau$. Then

[TABLE]

At the moment the first $N^{\prime}$ particles coalesce, the system becomes a two-particle system, and so

[TABLE]

Let us assume the particle at $x_{t}$ should not experience an abrupt discontinuity in its drift at the moment of coalescence, i.e. we want $dx_{\tau-\epsilon}\rightarrow dx_{\tau}$ as $\epsilon\rightarrow 0^{+}$. Equating the right hand sides of (2.20) and (2.21) as $\epsilon\rightarrow 0^{+}$ and using the property that $X^{(n)}_{\tau-\epsilon}\rightarrow Z_{\tau}$ for all $n\leq N^{\prime}$, we get

[TABLE]

meaning the $N^{\prime}$ particles must coalesce at the center of mass of the subsystem. This suggests that the diffusion coefficient of the newly-created particle positioned at $Z_{\tau}$ should be the same as the diffusion coefficient of the center of mass process of the first $N^{\prime}$ particles for $t<\tau$. By the independence of the processes $W^{(i)}_{t}$ for $1\leq i\leq N^{\prime}$ and the definition of the center of mass inside the limit on the right hand side of (2.22), we get

[TABLE]

or equivalently,

[TABLE]

Since $M^{\prime}=\sum m_{i}$, it follows that $f(x)=x^{2}(\sigma(x))^{2}$ must be additive, i.e. satisfies Cauchy’s functional equation,

[TABLE]

Under the physically relevant assumption that $f$ is continuous, solutions to this functional equation must be linear [21]. We therefore get

[TABLE]

as in the dynamics given in the beginning of the work in (1.3).

By the same reasoning, we expect $Z_{t}$ to be driven by the weighted noise of the center of mass, $W^{(cm)}_{t}$, given by

[TABLE]

The dynamics of the coalesced particle of mass $M^{\prime}$ at $Z_{t}$ for $t\geq\tau$ are therefore

[TABLE]

which in the presence of additional particles generalize to (1.3).

2.4 Evolution of a subsystem’s second moment

Let us compute the local second moment of the highly localized subsystem of the first $N^{\prime}$ particles in Figure 1. First, we ignore all interactions with the outside particles not in the colliding cluster, and therefore approximate that the local second moment,

[TABLE]

evolves according to (2.12) with the summation being taken over the indices of the particle participating in the collision,

[TABLE]

where

[TABLE]

As shown in Figure 3, such an approach appears to be qualitatively correct, but introduces an error which appears to grow in time. Let us now find a higher order approximation.

As a model for the system in Figure 1, consider a system consisting of two nearby particles of masses $m_{1}$ and $m_{2}$, and a third, distant particle of mass $m_{3}$, i.e. $\left|X^{(1)}_{t}-X^{(2)}_{t}\right|\ll\left|X^{(1)}_{t}-X^{(3)}_{t}\right|\approx\left|X^{(2)}_{t}-X^{(3)}_{t}\right|$. We wish to investigate how the third particle affects the second moment of the subsystem consisting of the first two particles,

[TABLE]

Using (1.3) and an application of Ito’s lemma, we can get an exact correction to the deterministic part of the approximating process $Q_{t}$ given in (2.30):

[TABLE]

We introduce the small parameter

[TABLE]

through which (2.33) may be approximated as

[TABLE]

where we assume $X^{(1)}_{t}-X^{(3)}_{t}\approx X^{(2)}_{t}-X^{(3)}_{t}\approx X^{(cm)}_{t}-X^{(3)}_{t}$ and $\theta$ is the angle between $X^{(2)}_{t}-X^{(1)}_{t}$ and $X^{(cm)}_{t}-X^{(3)}_{t}$.

A similar monopole approximation may be used when there are $N-2$ particles affecting the evolution of the second moment of the first two particles. Then,

[TABLE]

where $\theta_{i}$ is the angle between $X^{(2)}_{t}-X^{(1)}_{t}$ and $X^{(cm)}_{t}-X^{(i)}_{t}$, and the shorthand $V(x,y)=V(|x-y|)$ is used to simplify the expression.

By a similar argument, for an $N$ particle system with a cluster consisting of the first $N^{\prime}$ particles, we have

[TABLE]

where $\theta_{ijk}$ is the angle between $X^{(j)}_{t}-X^{(k)}_{t}$ and $\tilde{X}^{(cm)}_{t}-X^{(i)}_{t}$, with

[TABLE]

Heuristically, we see that as $\tilde{Y_{t}}\rightarrow 0$, the corrections in (2.38) vanish, the subsystem essentially becomes decoupled from the rest of the system, and the subsystem’s second moment $\tilde{Y_{t}}$ becomes a squared Bessel process of negative index by (2.17). Since the collision process (before the collision time) does not involve the creation or annihilation of particles, it appears that a highly-localized aggregate which is not decoupled from the rest of the system, but is nontheless undergoing collision, should still satisfy (2.17), i.e.

[TABLE]

where $\nu$ is as in (2.16). This informal argument suggests that for a very tight cluster, this is a sufficient condition for an aggregate to undergo collision. For a less tight cluster (even if it is separated), the contributions of the higher order corrections may prevent a collision from occurring.

3 Simulation of particle coalescence and dynamics

3.1 Overview

We employ a grid-particle approach for computing interparticle interactions, which avoids pairwise computations in (1.3) by introducing a continuous global potential which varies in time. We remark that similar ideas have been developed in the particle-in-cell literature (e.g. [6], [32]), but without coalescing stochastic particles.

We sidestep the challenge of numerically detecting singular point collisions by introducing an adaptive grid which identifies highly localized aggregates, the second moment of which is computed and simulated using the appropriate Wiener process (given by (2.14)) in order to identify a collision inside a timestep.

3.2 Full numerical method

The numerical method for the simulation of the coalescing particle system (1.3) combines the upcoming sections at every timestep in the following order:

Detect highly isolated clusters of particles with negative indices, which may collide with high probability within the upcoming time step. For each such cluster, compute the local second moment, $\tilde{Y_{t}}$. 2. 2.

Simulate the particle dynamics, using adaptive timestepping when appropriate. For each particle in the above clusters, record the total increment of the driving Wiener process over the full time step. 3. 3.

For each cluster, simulate the second moment over a time step, using (2.14). If the second moment hits zero, coalesce the cluster’s particles at their center of mass.

3.3 Detection of isolated aggregates

To detect particle collisions, we first apply a density-based clustering algorithm for finding isolated particle aggregates. Such clusters are then checked for collisions, as described in Section 3.5.

To find clusters, we form a coarse mesh which covers all the particles (in practice, we use a $1\times 1$ mesh). For each cell, we compute the square of its diagonal, $s^{2}$, and the second moment of the particles inside the cell, $\tilde{Y}$. We call a cell “separated” if

[TABLE]

where $\eta$ is some fixed constant (in practice, the authors use $\eta=0.1$). If a cell is not separated, and has more than two particles, then we refine the cell into four equally-sized cells, and repeat this procedure with each subcell.

A separated cell is kept if it is “collidable,” otherwise it is refined as well. A cell is collidable if its index $\nu$ is negative, and the second moment satisfies

[TABLE]

where $\alpha$ and $\beta$ are given as in (2.13), $\Phi$ is the normal distribution function, and $0<p\ll 1$ is some small probability. The interpretation of this inequality is that it excludes cells which may collide within the time step with very low probabilities.

3.4 Particle dynamics

Since $V(x,y)=\frac{1}{2\pi}\ln|x-y|$ is the fundamental solution to the Laplace operator, we can get a global potential for the the dynamics given in (1.3):

[TABLE]

where $c$ is the interaction field and $P$ is the empirical mass density. The case of a different $V$ can be treated similarly.

For the simulation of these particle dynamics, we discretize a computational domain as in figure 4, and use the particles to interpolate a mass density field $P_{ij}$ onto the field. We then numerically solve for the mean field $C_{ij}$. To advance by $\Delta t$ forward in time, we introduce adaptive time steps $\Delta\tau_{1},\ldots,\Delta\tau_{K(n,t)}$ (this is needed for stability reasons—see below), and use a forward Euler-Maruyama scheme to simulate the dynamics of each particle:

[TABLE]

where $N_{i}^{(n)}(0,1)$ is a normal Gaussian random variable. This bookkeeping of the noise is helpful for numerically detecting collisions, where we need the quantity

[TABLE]

i.e. the increment of the $n$th Wiener process $W^{(n)}_{t}$ between $t$ and $t+\Delta t$ (see Section 3.5).

We approximate $\nabla c(x)$ in two steps. First, we construct the gradient field $\nabla C_{ij}~{}=~{}(CX_{ij},CY_{ij})$ using the second order approximation

[TABLE]

Then we approximate $\nabla c(x,t)$ using a bilinear interpolation of the values of $\nabla C$ at the four nearest grid points. In the case that $x$ is not inside the computational domain, we use a monopole approximation:

[TABLE]

where $x_{cm}$ is the center of mass of the system. Since the primary novelty of this numerical method is in its applicability to colliding systems, an appropriately-chosen computational domain (i.e. one which overlaps with most of the mass of the system) will make use of the monopole approximation rarely. Nonuniform meshes may be used as well, but have not been observed to make a significant improvement in systems with most of the mass sufficiently away from the boundaries .

As described in [10], an adaptive time step is dynamically chosen such that the expected length of a particle’s jump does not exceed the mesh size. This is necessary to prevent spurious mass oscillations around singularities in $c(x,t)$. Since $V$ has logarithmic singularities, we expect that time steps can get as small as $\Delta\tau\sim\Delta x^{2}$.

The mass density field $P_{ij}$ is computed by bilinearly interpolating the mass of each particle onto the four nearest grid points. The result is divided by $\Delta x^{2}$, to get a mass density. This first-moment-preserving approach prevents particles from “self-interacting,” a phenomenon which creates an artificial flux towards grid points, as described in [10].

The mean field $C_{ij}$ is solved for on the computational domain using a standard finite-differences scheme:

[TABLE]

The monopole approximation is used for the boundary conditions:

[TABLE]

for $X_{ij}$ on the boundary on the computational domain, and $x_{cm}$ the center of mass of the particle system.

3.5 Detection of collisions in isolated aggregates

After all the particles are propagated over one time step, we consider the terminal cells returned by the algorithm given in Section 3.3. We approximate the evolution of the second moment inside each cell which is both separated and collidable. To do this, for each cell, we compute the quantity

[TABLE]

and coalesce all the particles at their new center of mass if $\Delta\tilde{Y}\leq 0$.

The increment $\Delta\tilde{W}_{t}$ is given by (2.31), i.e.

[TABLE]

The cost of computing the above sum can be significantly reduced using the following identity:

[TABLE]

from which

[TABLE]

where $X^{(cm)}_{t}$ is the center of mass of the cell.

We note the dynamics of the second moment may be approximated more accurately by taking advantage of the first order correction presented in Section 2.4, but the necessity of such corrections may be avoided by simply choosing a very small localization parameter $\eta$, as in (3.41).

4 Finite-time blow-up in hydrodynamic limits

4.1 Overview

We first show how the PKS particle system described in the introduction fits in the context of the present work. We then formally derive the hydrodynamic limit of a particle system with masses approaching zero nonuniformly, which we call the multispecies Patlak-Keller-Segel system (MPKS), and derive a finite-time blow-up condition. Finally, we show how the hydrodynamic limit of the system may be taken in such a way that the limit is a regularized MPKS system after the time of blow-up.

4.2 The Patak-Keller-Segel particle system

As already described in Section 1.4, the PKS is given by the following system of PDEs:

[TABLE]

where the boundary condition $\rho(x,t)\rightarrow 0$ as $|x|\rightarrow\infty$ is enforced, and mass is conserved: $\int\rho(x,t)dx=M$.

The PKS may be rewritten more compactly as an integrodifferential equation:

[TABLE]

where $V(x)=\frac{1}{2\pi}\ln|x|$, as before, is the fundamental solution of the Laplace operator. Observe that if $c$ is predetermined, then the first equation in (4.55) is the Fokker-Planck equation for the process

[TABLE]

It follows that for an $N$-particle system with positions $X^{(n)}_{t}$, the empirical mass density

[TABLE]

approximates the solution to the PKS $\rho$.

Since $\nabla c$ is unknown, we approximate it by the mean field created by the particles themselves: this is readily done making the substitution $c\rightarrow-V\ast P_{N}$, as suggested by (4.56). We arrive at

[TABLE]

This is simply the particle system described in the bulk of this work, with $m_{n}=M/N$ and the diffusion coefficient

[TABLE]

Thus, the PKS with total mass $M$ and diffusion coefficient $\mu$ can be viewed as the hydrodynamic limit of the particle system with the above parameters.

The particle system described in this work collides only when the index of the system (2.16) is negative. Similarly, the PKS forms singularities when the total mass is above the critical mass $M_{c}=8\pi\mu/\chi$ [1]. Let us show that these two criteria coincide in the hydrodynamic limit.

Substituting the necessary diffusion coefficient (4.60) into the definition of the Bessel index (2.16), we get the PKS index:

[TABLE]

As per the classification of the origin for the second moment, listed in Section 2.2, we have that a finite-time collision will occur when $\nu\leq-1$. This criterion applied to (4.64) reduces exactly to $M>8\pi\mu/\chi$—the necessary and sufficient condition for finite-time blow-up in the PKS.

4.3 Post-blow-up PKS and particle coalescence

The PKS has been regularized and investigated post-blow-up in several works, including [30, 31] and [7]. Although the post-blow-up dynamics are slightly different in the two works, they share the common feature that the density becomes a measure, and splits into a regular, and an atomic component consisting of $K_{t}$ point masses:

[TABLE]

where the $n$th atomic component has a smoothly-evolving mass $M_{n}(t)\geq 8\pi\mu/\chi$, supported on a point moving along a smooth path. The point masses may emerge or collide, and thus their number $K_{t}$ varies in time. Mass is locally transferred from the regular component to each atomic component as

[TABLE]

With these dynamics, it can be shown [7] that the second moment of this system evolves as

[TABLE]

where $\bar{M}=M-\sum_{i=1}^{K_{t}}M_{i}(t)$ is the mass of the regular component (we note the quantity of interest in the PKS literature is typically the unnormalized second moment, which we choose to normalize, due to its geometric interpretation).

In the context of the PKS particle system, we expect light, uncoalesced particles to correspond to the regular component of the solution to the PKS, and each massive, coalesced particles to correspond to point mass in the atomic component of the solution to the PKS. By the previous section, such particles should only have mass above $8\pi\mu/\chi$, as in the PKS. Let us recover equation (4.67) using the particle system, assuming that this correspondence is true.

Consider a PKS system with smooth initial conditions, which blows up in finite time, and has an atomic component of mass $M_{1}$, consisting of one point mass, at time $t=T$. Now consider a PKS particle system, initialized with $N_{0}$ particles distributed according to the initial conditions given to the PKS PDE. The second moment $Y_{t}$ evolves according to

[TABLE]

where $\alpha$ and $\beta$ are given in (2.13), with $\tilde{\mu}=\mu M/N_{0}$, and $m_{n}=M/N_{0}$ initially. Near $t=T$, there should be one massive particle, consisting of $k$ coalesced light particles. Plugging this into (2.13), we get:

[TABLE]

As $N_{0}\rightarrow\infty$, we get $\beta\rightarrow 0$, and $Y_{t}$ becomes deterministic:

[TABLE]

consistent with (4.67) for a single point mass. A similar argument can be used to derive (4.67) fully.

4.4 Hydrodynamic limit to the multispecies PKS model

We remark that the sign of the PKS particle system’s index (4.63) becomes independent of $N$ as $N\rightarrow\infty$. This convenient property occurs only because $\tilde{\mu}\sim 1/N$, and is actually independent of the the particle masses, as long as the total sum of the particle masses is fixed and the mass of each individual particle approaches zero. Thus the question of the limiting system when individual particles approach [math] nonuniformly arises naturally.

As a first basic example, let us consider the system

[TABLE]

where $N=2N^{\prime}$, $M=M_{1}+M_{2}$, $m_{i}~{}=~{}M_{1}/N^{\prime}$ for $i\leq N^{\prime}$ and $m_{i}~{}=~{}M_{2}/N^{\prime}$ for $i>N^{\prime}$. That is, we break up the system into two families, the first family containing $N^{\prime}$ particles of uniform mass $m_{a}=M_{1}/N^{\prime}$, and the second family containing $N^{\prime}$ particles of uniform mass $m_{b}=M_{2}/N^{\prime}$. The particle dynamics are then given by

[TABLE]

where $P_{1}$ and $P_{2}$ are the empirical mass densities of the particles of the first and second mass:

[TABLE]

Appealing once more to the formal derivation of the hydrodynamic limit described earlier, we expect that $P_{i}$ approximates $\rho_{i}$ in the limit $N\rightarrow\infty$, where

[TABLE]

The above system can be seen as a “two species” PKS model, in which two species attract each other through the same mechanism, but have different average diffusion rates.

Similarly, we may break the system up into $K$ families, each family of total mass $M_{i}$ and containing $N_{i}$ particles of uniform mass $M_{i}/N_{i}$. We take the hydrodynamic limit by fixing $\eta_{i}>0$ for $1\leq i\leq K$ such that

[TABLE]

and letting $N\rightarrow\infty$ in such a way that

[TABLE]

Then $P_{i}\rightarrow\rho_{i}$, where

[TABLE]

with $\int\rho_{i}=M_{i}$ and

[TABLE]

which can be interpreted as $\mu$ scaled by the ratio of the overall system’s average particle mass, to the $i$th family’s particle mass. We will refer to (4.79) as the “multispecies Patlak-Keller-Segel system” (MPKS).

The excluded case $\eta_{i}=0$ corresponds to a mass of $M_{i}$ being supported entirely on a singular component of the solution post-blow-up.

4.5 Formation of singularities in the MPKS

As can be seen from (4.62), the sign of the index of a particle system that’s taken to its hydrodynamic limit becomes independent of the number of particles, and can therefore fully collide in finite time, if a specific mass condition is satisfied. In the PDE, this corresponds to a finite-time blow up. Let us verify that this is indeed the case.

Assume arbitrary diffusion coefficients $\mu_{i}$. Let $P(x,t)=\sum_{i=1}^{K}\rho_{i}(x,t)$ be the total mass density of an MPKS system. Then $\int_{\mathbb{R}^{2}}P(x,t)dx=M_{1}+\cdots+M_{K}=M$, and

[TABLE]

To show the existence of finite-time blow-up, define the second moment of the system,

[TABLE]

and compute its derivative:

[TABLE]

where the detailed computation is given in Appendix B. Thus for constants satisfying

[TABLE]

the second moment vanishes in finite time, but the total mass is conserved–thus implying the formation of a singularity.

As an aside, we remark that the the formula given by (4.83) remains valid when each component has a different chemosensitivity $\chi_{i}$. Furthermore, we note that the blow-up condition (4.84) is satifised when $M>\max(8\pi\mu_{i}/\chi)$, i.e. the MPKS forms a singularity when its total system mass is greater than the classic PKS critical mass for each separate components. Recalling the special structure of the diffusion coefficients in the hydrodynamic limit of the particle system (4.80), we see that the blow-up condition (4.84) coincides with the full particle system collision condition $\nu_{PKS}<0$, where $\nu_{PKS}$ is as in (4.62).

For two species, the system was investigated in [5], where initial data were classified in terms of having solutions which either blow up in finite time, or are global in time. Interestingly, that work showed that there exist initial data corresponding to finite time blow-up, for which the second moment is increasing, i.e. $F^{\prime}(t)>0$—in analogy with (2.18). An optimal classification was obtained for a disc domain in [8], though questions, such as if blow up occurs simultaneously in all components, remain (this question was affirmatively answered for the radial case in [9]). In Section 5.3, we investigate how the second moments of components of the two species MPKS evolve in the regime that a singularity forms in finite time with $F^{\prime}(t)>0$.

We expect that the MPKS can be regularized past blow-up times using a singular perturbation limit, as was done in [30, 31] for the PKS, and proposed in [22] for the MPKS. In this case, the presented method is well-suited for the investigation of this regularization.

4.6 More general $V$

As the particle system dynamics are equally valid for choices of $V$ which are not scaled logarithms, we left some formulas somewhat general, simply in terms of the derivatives of $V$. Particle coalescence, however, strongly depends on there being a logarithmic singularity in $V$. This is necessary to connect collisions to the Bessel process.

We note that, in the plane, the fundamental solution to a radially-symmetric, elliptic operator $\mathcal{L}$ with sufficiently regular coefficients (as in (1.2)) has logarithmic singularities. It therefore follows that the discussion above applies in the case when $V$ is such a fundamental solution. That is, suppose $V(x,y)\sim\gamma\ln|x-y|$ as $|x-y|\rightarrow 0$. Then the index formulas used in the previous sections should be replaced by the following index:

[TABLE]

Applying the same procedure as in Section 4.4 will result in a hydrodynamic limit which solves

[TABLE]

with post-blow-up dynamics similar to the ones given for the PKS in [30, 31] and [7].

5 Numerical simulations

5.1 Overview

One application of this work is in developing a numerical method for the PKS and PKS-like systems, which is able to handle the formation of singularities, as well as post-blow-up dynamics. Let us consider two example applications, for which we explicitly know the expected behavior: the evolution of the second moment for the PKS, pre- and post-blow-up, as given in (4.67), and blow-up with an increasing second moment in the two species MPKS, as described in Figure 5. In the first, we will show that the second moment of our particle approximation evolves as predicted by [7] both before and after blow-up, confirming that our numerical method correctly transitions from approximating smooth solutions to the PKS, to approximating measure-valued solutions. In the second, we will see how the second moment of the components of a two species MPKS system with masses inside the shaded region in Figure 5 evolve, thus giving numerical evidence to the idea that blow-up in this regime occurs via one contracting, and one expanding component.

We remark the presented numerical method is parallelizable, and scales approximately linearly with the number of particles. It can therefore be used to simulate a large number (on the order of millions) of particles very quickly. Averaging over such large ensembles reduces observed stochastic fluctuations to a minimum, as may be noted from the examples in this section.

5.2 Regularized PKS

For the first example, we reproduce the equation (4.67) for the PKS second moment:

[TABLE]

Thus, the graph of the second moment of a critical PKS system will initially appear linear, then decelerate, and then–depending on the mass distribution–will either become linear again (with a different slope), or continually change its slope due to nonstop mass transfer to the atomic component. Using the numerical method developed in this work, this second moment evolution can be observed. For a PKS system with mobility $\mu$ and chemosensitivity $\chi$, we associate an $N_{0}$-particle coalescing particle system with $\tilde{\mu}=\mu M/N_{0}$ and $m_{n}=M/N_{0}$, and approximate $\rho$ by the empirical mass density. As this particle approximation has been shown to be effective in approximating the PKS pre-blow-up [12, 10], we specifically concentrate on the formation and detection of singularities.

5.2.1 Mass transfer to singularity

In particular, we consider the case $\chi=\mu=1$, with total mass six times the critical mass, $M=6\cdot 8\pi$. We split the mass amongst a small bump function of mass $M_{1}=4\cdot 8\pi$ supported on a disc of unit radius, which is separated far away from a bump function of mass $M_{2}=2\cdot 8\pi$ that’s supported on an ellipse with axes $1$ and $7$. These initial initial conditions are chosen so as to make the solution initially exhibit a linear decay of the second moment, then a sudden change of slopes due to the rapid formation of a singularity caused by the first bump function, and finally–a continuous deceleration, due to continual mass transfer from the lighter bump function to the formed atomic component. With the chosen parameters, the first two rates of change of the second moment should be $-20$ and $-12$. This can be observed in Figure 6. Further in time, the gradual transfer of mass may be seen as well, as shown in Figure 7.

The underlying particle dynamics and collisions are illustrated in Figure 8, where each snapshot corresponds to qualitatively different rates of change of the second moment in Figure 6 and Figure 7: sudden mass coalescence of a tight aggregate (switch of slopes in Figure 7, and $t=0.050$ and $t=0.100$ in Figure 8), attraction of mass without coalescence (linear decay in Figure 6, and $t=0.100$ and $t=0.650$ in Figure 8), continuous slow and fast mass absorption (gradual deceleration in Figure 7, and $t=0.650$ and $t=0.950$ in Figure 8), and the transformation of the PKS system to being essentially singular (flat part of the figure in Figure 7, and $t=2.200$ in Figure 8).

5.2.2 Interaction of singularities

In another experiment, we initialize a system in which two singularities form and interact, as described in Figure 9. In this special case, the second moment is simply the square of the distance between the two singularities, the graph of which should be piecewise linear (as observed). We note that the numerical coalescence procedure avoids the “washing out” effect near the collision time in Figure 7 of [10].

5.3 Expanding MPKS with blow-up

For the second example, we simulate blow-up with an increasing total second moment for the two species Keller-Segel system:

[TABLE]

with $\int\rho_{1}=M_{1}$ and $\int\rho_{2}=M_{2}$. The interest in this phenomenon is described in Figure 5 and Section 4.4. In particular, we show that when a two species PKS system is in this regime, the second moment of one component increases linearly, while the other decreases. Such semi-decoupled behavior was suggested in [5]. We remark that the numerical method presented is well-suited for this investigation, as it can simulate the system in the entire plane.

We approximate this two system using $N_{0}$ particles, the first $N_{1}=\lfloor\eta_{1}N_{0}\rfloor$ of which have particle masses $M_{1}/N_{1}$, and distributed on the plane according to $\rho_{1}(\cdot,0)$. Similarly, the last $N_{2}=N_{0}-N_{1}$ particles have masses $M_{2}/N_{2}$, and are distributed on the plane according to $\rho_{2}(\cdot,0)$. Using (4.79) and (4.80), we see that

[TABLE]

where

[TABLE]

The particle system’s diffusion coefficient $\tilde{\mu}$ is then

[TABLE]

Thus, for a two species MPKS system with component masses $M_{1},M_{2}$ and diffusion coefficients $\mu_{1},\mu_{2}$, we associate an $N_{0}$ particle system with two different possible particle masses. The diffusion coefficient for (1.3) is given by (5.91). In this sense, the purpose of $\mu$ in (5.90) is auxiliary.

When $\mu_{1}>2\mu_{2}$, it is always possible to choose component masses which will force a radially-symmetric system to blow-up with increasing second moment. In this case, $M_{1}^{max}$ and $M_{2}^{max}$ in Figure 5 can be shown to be

[TABLE]

For the experiments in this section, we simulate the two species system as described above, and choose the convenient parameters

[TABLE]

which correspond to the auxiliary parameters

[TABLE]

For the above masses, we consider three different initial conditions. Each respective solution exhibits linear growth in the first component’s second moment, and decay in the second component’s second moments, but at rates which depend on the initial distribution of mass. In particular, we choose the following initial conditions:

Radially-symmetric component initial data. We initialize both components as bump functions supported on a disc of radius $a=0.35$ and centered at the origin. 2. 2.

Non-symmetric component initial data. We initialize the first component as a bump function of radius $a$ and centered at the origin, and the second component as a bump function supported on an ellipse centered at $(0.1,0)$ with axes $2a$ and $a/2$, with the major axis parallel to the $y$-axis. 3. 3.

Component initial data on disjoint support. We initialize each component on a bump function supported on a disc of radius $a$, where the first component is centered at $(a,-a)$, and the second at $(-a,a)$.

The results of these simulations can be seen in Figure 10. We note that although both components change linearly, their rates of change appear to depend on the initial conditions.

6 Conclusion

We investigated a planar particle system with nonuniform particle masses, in which particles interact via a logarithmically-singular kernel. As post-collision dynamics in such a system are undefined, we used the idea of particle coalescence in order to propagate the system further in time, and connected it to the theory of the squared Bessel process. We exploited this connection to develop an efficient numerical method for the simulation of the system, which has applications in the numerical approximation and regularization of a wide range of nonlinear Fokker-Planck equations, such as the multispecies Patlak-Keller-Segel model.

As mentioned before, properties of singularity formation in the MPKS are not fully understood, and have somewhat unexpected behavior, when compared to the PKS. For instance, singularities may form while the system’s second moment is increasing. It would be interesting to further connect existing results with predicting a nonuniform particle system’s behavior post-collision.

The question of coalescence in a system with memory arises naturally, as an analogue to the parabolic Keller-Segel model. In this case, the field $c(x,t)$ is replaced with the solution to the following equation,

[TABLE]

which has the more biologically-meaningful intepretation of a chemoattractant which thermalizes at a finite rate, diffuses, decays, and is produced by the particles. This system will be investigated in future works.

7 Acknowledgements

This material is based upon work supported by the National Science Foundation under Grant DMS-1056471. The support is gratefully acknowledged.

The authors would like to thank Hailiang Liu for references related to the two species Keller-Segel model.

Appendix A Subtraction formula for indices

If $\nu_{i}$ is the index of the full system described in Figure 1, $\nu_{f}$ is the index of the same system after the particles inside the dashed lines coalesce, and $\overline{\nu}$ is the index of the subsystem inside the dashed line, then using (2.16) we have

[TABLE]

from which it follows that

[TABLE]

Appendix B MPKS second moment

The evolution of the second moment of the MPKS can be computed as follows:

[TABLE]