Diffusion of particles with short-range interactions

TL;DR

This paper derives a nonlinear diffusion equation for interacting Brownian particles with short-range repulsive potentials, comparing its accuracy with mean-field, Kirkwood, and Monte Carlo methods.

Contribution

It systematically derives a continuum model for short-range interactions using matched asymptotic expansions and evaluates its performance against other approximation methods.

Findings

The continuum model performs best for very repulsive short-range potentials.

Mean-field approximation is suitable for long-range interactions.

Kirkwood superposition approximation is accurate for both short- and long-range potentials, but computationally intensive.

Abstract

A system of interacting Brownian particles subject to short-range repulsive potentials is considered. A continuum description in the form of a nonlinear diffusion equation is derived systematically in the dilute limit using the method of matched asymptotic expansions. Numerical simulations are performed to compare the results of the model with those of the commonly used mean-field and Kirkwood-superposition approximations, as well as with Monte Carlo simulation of the stochastic particle system, for various interaction potentials. Our approach works best for very repulsive short-range potentials, while the mean-field approximation is suitable for long-range interactions. The Kirkwood superposition approximation provides an accurate description for both short- and long-range potentials, but is considerably more computationally intensive.

Tables

Table 1. Table 1: Values of the nonlinear diffusion coefficient obtained via matched asymptotics [ α u subscript 𝛼 𝑢 \alpha_{u} in Eq.   49b ] or via closure [ α ¯ u subscript ¯ 𝛼 𝑢 \overline{\alpha}_{u} in Eq.   13b ] for various types of interaction potential in 2 or 3 dimensions.

	$𝒅=\mathrm{𝟏}$		$𝒅=\mathrm{𝟐}$		$𝒅=\mathrm{𝟑}$
	${\alpha }_u$	${\overline{\alpha }}_u$	${\alpha }_u$	${\overline{\alpha }}_u$	${\alpha }_u$	${\overline{\alpha }}_u$
$u_{\text{HS}}={\left(ϵ/r\right)}^{\infty }$	2	$\nexists$	$\pi$	$\nexists$	$1.33\pi$	$\nexists$
$u_{\text{SS}}={\left(ϵ/r\right)}^4$	$2.45$	$\nexists$	$1.77\pi$	$\nexists$	$4.83\pi$	$\nexists$
$u_{\text{LJ}}={\left(\frac{ϵ}{r}\right)}^{12}-{\left(\frac{ϵ}{r}\right)}^6$	1.63	$\nexists$	$0.555\pi$	$\nexists$	$0.154\pi$	$\nexists$
$u_{\text{EX}}=e^{-r/ϵ}$	$1.59$	$2$	$1.78\pi$	$2\pi$	$7.54\pi$	$8\pi$
$u_{\text{YU}}=\frac{ϵ}{r}{e}^{-r/ϵ}$	$1.80$	$\nexists$	$1.25\pi$	$2\pi$	$3.47\pi$	$4\pi$

Equations

$$\mathrm{d}{\mathbf{X}}_i(t)=\sqrt{2D}\mathrm{d}{\mathbf{W}}_i(t)+\mathbf{f}({\mathbf{X}}_i(t))\mathrm{d}t-\sum _{j\ne i}{\nabla }_{{\mathbf{x}}_i}u(\Vert {\mathbf{X}}_i(t)-{\mathbf{X}}_j(t)\Vert )\mathrm{d}t,$$

$$\frac{\partial p}{\partial t}={\nabla }_{\mathbf{x}}\left[D{\nabla }_{\mathbf{x}}p-\mathbf{f}(\mathbf{x})p+(N-1)p{\nabla }_{\mathbf{x}}(u\ast p)\right],$$

$$\frac{\partial p}{\partial t}={\nabla }_{\mathbf{x}}\cdot \left[D{\nabla }_{\mathbf{x}}p-\mathbf{f}(\mathbf{x})p+{\alpha }_u{ϵ}^d(N-1)p{\nabla }_{\mathbf{x}}p\right],$$

$$U(\vec{x})=\sum _{1\le i<j\le N}u(\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert ),$$

$${\mathbf{g}}_i(\vec{x})=-{\nabla }_{{\mathbf{x}}_i}U(\vec{x})=-\sum _{j\ne i}{\nabla }_{{\mathbf{x}}_i}u(\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert ).$$

$$\frac{\partial P}{\partial t}(\vec{x},t)={\vec{\nabla }}_{\vec{x}}\cdot \left[{\vec{\nabla }}_{\vec{x}}P-\vec{F}(\vec{x})P+{\nabla }_{\vec{x}}U(\vec{x})P\right]\text{in}{\Omega }^N,$$

$$P(\vec{x},0)=P_0(\vec{x}),$$

$$p({\mathbf{x}}_1,t)={\int }_{{\Omega }^{N-1}}P(\vec{x},t)\mathrm{d}{\mathbf{x}}_2\cdots {\mathbf{x}}_N.$$

$$\frac{\partial p}{\partial t}({\mathbf{x}}_1,t)={\nabla }_{{\mathbf{x}}_1}\cdot \left[{\nabla }_{{\mathbf{x}}_1}p-\mathbf{f}({\mathbf{x}}_1)p+\mathbf{G}({\mathbf{x}}_1,t)\right],$$

$$P_2({\mathbf{x}}_1,{\mathbf{x}}_2,t)={\int }_{{\Omega }^{N-2}}P(\vec{x},t)\mathrm{d}{\mathbf{x}}_3\cdots \mathrm{d}{\mathbf{x}}_N$$

$$P_2({\mathbf{x}}_1,{\mathbf{x}}_2,t)=p({\mathbf{x}}_1,t)p({\mathbf{x}}_2,t).$$

$$\mathbf{G}({\mathbf{x}}_1,t)=-{ϵ}^{d-1}(N-1)p({\mathbf{x}}_1,t){\int }_{{\mathbb{R}}^d}\left[p({\mathbf{x}}_1,t)+ϵ\tilde{\mathbf{x}}\cdot {\nabla }_{{\mathbf{x}}_1}p({\mathbf{x}}_1,t)\right]{\nabla }_{\tilde{\mathbf{x}}}\tilde{u}(\Vert \tilde{\mathbf{x}}\Vert )\mathrm{d}\tilde{\mathbf{x}}+\cdots ,$$

$$\mathbf{G}({\mathbf{x}}_1,t)\sim {ϵ}^d(N-1)p({\mathbf{x}}_1,t){\nabla }_{{\tilde{\mathbf{x}}}_1}p({\mathbf{x}}_1,t){\int }_{{\mathbb{R}}^d}\tilde{u}(\Vert \tilde{\mathbf{x}}\Vert )\mathrm{d}\tilde{\mathbf{x}}+O({ϵ}^{d+\delta }).$$

$$\displaystyle\frac{\partial p}{\partial t}=\nabla_{{\bf x}_{1}}\cdot\left\{\big{[}1+\overline{\alpha}_{u}(N-1)\epsilon^{d}p\big{]}{\nabla}_{{\bf x}_{1}}p-{\bf f}({\bf x}_{1})p\right\},$$

$${\overline{\alpha }}_u={\int }_{{\mathbb{R}}^d}u(ϵ\Vert \mathbf{x}\Vert )\mathrm{d}\mathbf{x}.$$

$$\frac{\partial \rho }{\partial t}={\nabla }_{{\mathbf{x}}_1}\cdot \left({\nabla }_{{\mathbf{x}}_1}\rho +\int Q({\mathbf{x}}_1,{\mathbf{x}}_2,t){\nabla }_{{\mathbf{x}}_1}u(r)\mathrm{d}{\mathbf{x}}_2\right),$$

$$\frac{\partial Q}{\partial t}={\nabla }_{{\mathbf{x}}_1}\cdot \left[{\nabla }_{{\mathbf{x}}_1}Q+Q{\nabla }_{{\mathbf{x}}_1}u(r)\right]+{\nabla }_{{\mathbf{x}}_2}\cdot \left[{\nabla }_{{\mathbf{x}}_2}Q+Q{\nabla }_{{\mathbf{x}}_2}u(r)\right].$$

$${\rho }_s({\mathbf{x}}_1)={\rho }_0,Q_s({\mathbf{x}}_1,{\mathbf{x}}_2)={\rho }_0^2{g}_0(\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert ),$$

$$g_0(r)=e^{-u(r)}.$$

$$Q({\mathbf{x}}_1,{\mathbf{x}}_2,t)=\rho ({\mathbf{x}}_1,t)\rho ({\mathbf{x}}_2,t)g({\mathbf{x}}_1,{\mathbf{x}}_2,t),$$

$$Q({\mathbf{x}}_1,{\mathbf{x}}_2,t)\approx {\rho }_0^2{g}_0(r)+{\rho }_1({\mathbf{x}}_1,t){\rho }_0{g}_0(r)+{\rho }_0{\rho }_1({\mathbf{x}}_2,t)g_0(r)+{\rho }_0^2{g}_1({\mathbf{x}}_1,{\mathbf{x}}_2,t).$$

$$\frac{\partial\rho_{1}}{\partial t}=\nabla_{{\bf x}_{1}}\cdot\Big{[}\nabla_{{\bf x}_{1}}\rho_{1}+\rho_{0}\int g_{0}(r)\rho_{1}({\bf x}_{2},t)\nabla_{{\bf x}_{1}}u(r)\,\mathrm{d}{\bf x}_{2}\\ +\rho_{0}^{2}\int g_{1}({\bf x}_{1},{\bf x}_{2},t)\nabla_{{\bf x}_{1}}u(r)\,\mathrm{d}{\bf x}_{2}\Big{]}.$$

$$\frac{\partial {\rho }_1}{\partial t}={\nabla }_{{\mathbf{x}}_1}\cdot \left[{\nabla }_{{\mathbf{x}}_1}{\rho }_1+{\rho }_0\int g_0(r){\rho }_1({\mathbf{x}}_2,t){\nabla }_{{\mathbf{x}}_1}u(r)\mathrm{d}{\mathbf{x}}_2\right].$$

$$\frac{\partial {\rho }_1}{\partial t}={\nabla }_{{\mathbf{x}}_1}\cdot \left[(1+{\alpha }_u{ϵ}^d{\rho }_0){\nabla }_{{\mathbf{x}}_1}{\rho }_1\right],$$

$$\begin{gathered} \frac{\partial P_2}{\partial t}={\nabla }_{{\mathbf{x}}_1}\cdot \left[{\nabla }_{{\mathbf{x}}_1}{P}_2-\mathbf{f}({\mathbf{x}}_1)P_2+{\nabla }_{{\mathbf{x}}_1}u(\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert )P_2+{\mathbf{G}}_2({\mathbf{x}}_1,{\mathbf{x}}_2,t)\right] \\ +{\nabla }_{{\mathbf{x}}_2}\cdot \left[{\nabla }_{{\mathbf{x}}_2}{P}_2-\mathbf{f}({\mathbf{x}}_2)P_2+{\nabla }_{{\mathbf{x}}_2}u(\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert )P_2+{\mathbf{G}}_2({\mathbf{x}}_2,{\mathbf{x}}_1,t)\right], \end{gathered}$$

$${\mathbf{G}}_2({\mathbf{x}}_1,{\mathbf{x}}_2,t)=(N-2){\int }_{\Omega }{\nabla }_{{\mathbf{x}}_1}u(\Vert {\mathbf{x}}_1-{\mathbf{x}}_3\Vert )P_3({\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3,t)\mathrm{d}{\mathbf{x}}_3,$$

$$P_3({\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3,t)={\int }_{{\Omega }^{N-3}}P(\vec{x},t)\mathrm{d}{\mathbf{x}}_4\cdots \mathrm{d}{\mathbf{x}}_N$$

$$P_3({\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3,t)=\frac{P_2({\mathbf{x}}_1,{\mathbf{x}}_2,t)P_2({\mathbf{x}}_1,{\mathbf{x}}_3,t)P_2({\mathbf{x}}_2,{\mathbf{x}}_3,t)}{p({\mathbf{x}}_1,t)p({\mathbf{x}}_2,t)p({\mathbf{x}}_3,t)},$$

$$\begin{gathered} \frac{\partial P_2}{\partial t}={\nabla }_{{\mathbf{x}}_1}\cdot \left[{\nabla }_{{\mathbf{x}}_1}{P}_2-\mathbf{f}({\mathbf{x}}_1)P_2+{\nabla }_{{\mathbf{x}}_1}u(\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert )P_2\right] \\ +{\nabla }_{{\mathbf{x}}_2}\cdot \left[{\nabla }_{{\mathbf{x}}_2}{P}_2-\mathbf{f}({\mathbf{x}}_2)P_2+{\nabla }_{{\mathbf{x}}_2}u(\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert )P_2\right], \end{gathered}$$

$$P_{\text{o}}({\mathbf{x}}_1,{\mathbf{x}}_2,t)=p({\mathbf{x}}_1,t)p({\mathbf{x}}_2,t)+ϵP_{\text{o}}^{(1)}({\mathbf{x}}_1,{\mathbf{x}}_2,t)+\cdots ,$$

$$\begin{gathered} {ϵ}^2\frac{\partial \tilde{P}}{\partial t}=2{\nabla }_{\tilde{\mathbf{x}}}\cdot \left[{\nabla }_{\tilde{\mathbf{x}}}\tilde{P}+{\nabla }_{\tilde{\mathbf{x}}}\tilde{u}(\Vert \tilde{\mathbf{x}}\Vert )\tilde{P}\right]+ϵ{\nabla }_{\tilde{\mathbf{x}}}\cdot \left\{\left[\mathbf{f}({\tilde{\mathbf{x}}}_1)-\mathbf{f}({\tilde{\mathbf{x}}}_1+ϵ\tilde{\mathbf{x}})\right]\tilde{P}\right\} \\ -ϵ{\nabla }_{{\tilde{\mathbf{x}}}_1}\cdot \left[2{\nabla }_{\tilde{\mathbf{x}}}\tilde{P}+{\nabla }_{\tilde{\mathbf{x}}}\tilde{u}(\Vert \tilde{\mathbf{x}}\Vert )\tilde{P}\right]+{ϵ}^2{\nabla }_{{\tilde{\mathbf{x}}}_1}^2\tilde{P}-{ϵ}^2{\nabla }_{{\tilde{\mathbf{x}}}_1}\cdot \left[\mathbf{f}({\tilde{\mathbf{x}}}_1)\tilde{P}\right]. \end{gathered}$$

Full text

{AMS}

35Q84, 60J70, 82C31

1 Introduction

Nonlinear diffusion equations are often used to describe a system of interacting particles at the continuum level. These play a key role in various physical and biological applications, including colloidal systems and granular gases, ion transport, chemotaxis, neural networks, and animal swarms. These continuum models are important as tools to explain how individual-level mechanisms give rise to population-level or collective behavior. Closure approximations such as the mean-field closure are often used to obtain the continuum model, but depending on the type of interactions they can lead to substantial errors. In this paper we present a new approach that is suited for short-range repulsive interactions.

A typical model for a system of $N$ interacting particles is to assume each particle evolves according to the overdamped Langevin dynamics and interacts with the other particles via a pairwise interaction potential $u$, so that

[TABLE]

for $i=1,\dots,N$, where ${\bf X}_{i}(t)\in\Omega\subset\mathbb{R}^{d}$ is the position of the $i$th particle at time $t$, $D$ is the diffusion constant, ${\bf W}_{i}(t)$ denotes a $d$-dimensional Brownian motion, and ${\bf f}:\Omega\to\mathbb{R}^{d}$ is an external force.

Despite the conceptual simplicity of the stochastic model Eq. 1, it can be computationally intractable for systems with a large number of interacting particles $N$, since the interaction term has to be evaluated for all particle pairs. In such cases, a continuum description of the system, based on the evolution of the population-averaged spatial concentration instead of individual particles, becomes attractive. Depending on the nature of the interaction potential different averaging techniques may be suitable. Interactions can be broadly classified into local and nonlocal depending on the range of interaction between particles. Nonlocal interactions are associated with a long-range or ultra-soft interaction potential $u$ in Eq. 1. Then every particle can interact not only with its immediate neighbors but also with particles far away, and one can use a mean-field approximation to obtain a partial differential equation (PDE) for the one-particle probability density $p({\bf x},t)$ of finding a given particle at position $\bf x$ at time $t$. The standard mean-field approximation (MFA) procedure applied to the microscopic model Eq. 1 gives

[TABLE]

where $u\ast p=\int u(\|{\bf x}-{\bf y}\|)p({\bf y})\mathrm{d}{\bf y}$. The main assumption in writing down Eq. 2 is that, in deriving the interaction term, particles can be treated as though they were uncorrelated. Since one is often interested in system with a large number of particles $N$, it is common to consider the number density ($\rho=Np$) and take the limit in which the number of particles and the volume tend to infinity, while keeping the average number density constant, that is, $N,V\to\infty$, $N/V=\rho_{0}$. This is known as the thermodynamic limit [12, §2.3], and results in equation (2) for $\rho$ without the factor $N-1$.

While the mean-field approximation Eq. 2 is convenient and leads to an accurate description for long-range interactions, such as with Coulomb interactions [9], it fails when considering relatively strong repulsive short-range potentials. Sometimes Eq. 2 does not make sense since the convolution does not exist; in other cases, it results in a poor model of the system because the underlying assumptions of the method are not satisfied, as we will discuss later. In particular, the model Eq. 2 does not make sense with hard-core repulsive interactions, which are commonly used to model excluded-volume effects in biological and social contexts [4]. A common way to circumvent this is to assume particles are restricted to a lattice, giving rise to so-called on-lattice models. The most common of these is the simple exclusion model, in which a particle can only move to a site if it is presently unoccupied [19]. One can derive an analogous continuous limit to Eq. 2 using Taylor expansions [6]. However, it turns out that with identical particles the interaction terms (analogous to the convolution term in Eq. 2) cancel out unless the external force $\bf f$ is nonzero.

Because of the issues that the MFA has in particle systems with short-range interaction potentials, one must often resort to numerical regularisations [15] or alternative closure approximations. There are a large number of closure approximations, and choosing the right closure for a given pair potential is “an art in itself” [26] due to the phenomenological nature of the approach. Each choice of closure relation results in a different approximate equation for the density $p$, and it is not clear a priori whether it will be a good approximation. As discussed above, one may even obtain an equation that does not make sense or is ill-posed.

One class of closures, including the MFA and the Kirkwood superposition approximation (KSA), impose a relation between the $n$th and $(n+1)$th density functions in the BBGKY hierarchy (see Section 2). Whereas the MFA closes at $n=1$ (writing the two-particle density function in terms of the one-particle density), KSA closes at $n=2$, approximating the three-particle density function as a combination of the one- and two-particle density functions [18]. While the KSA originated in the field of statistical mechanics, where it has been the basis of a whole theory, it has recently also been used in biological applications to obtain closed equations for a system of biological cells [3, 21]. However, the resulting KSA model is quite complicated to solve and the MFA remains the most commonly used approximation.

Another class of closure relations, similar in spirit to KSA, is based on the Ornstein–Zernike (OZ) integral equation [12]. Here the pair correlation function is decomposed into a ‘direct’ part and an ‘indirect’ part. The latter is mediated through (and integrated over) a third particle. In addition, the OZ equation requires a further closure assumption providing an additional relation between the direct and indirect correlations. Commonly used closures, for hard spheres and soft spheres respectively, include the Percus–Yevick approximation and the hypernetted chain approximation.

In this paper we are interested in systems such as Eq. 1 with short-range repulsive interactions for which MFA fails. We will employ an alternative averaging method to obtain a continuum description of the system based on matched asymptotic expansions (MAE). Unlike the mean-field approach, this method is a systematic asymptotic expansion which does not rely on the system size being large. It is valid for low concentrations, exploiting a small parameter $\epsilon$ arising from the short-range potential and the typical separation between particles. The result is a nonlinear advection-diffusion equation of the form

[TABLE]

where the coefficient $\alpha_{u}$ depends on the interaction potential $u$ and $d$ is the dimension of the physical space.

The remainder of the paper is organized as follows. In Section 2 we introduce the Fokker–Planck PDE for the joint probability density of the particle system; this is another individual-based description equivalent to the Langevin stochastic differential equation Eq. 1. In Section 3 we discuss three common closure approximations which reduce the Fokker–Planck equation to a population-level PDE. In Section 4 we present our alternative approach to closure based on MAE, and derive equation Eq. 3. In Section 5 we test the models obtained from the different methods against each other and stochastic simulations of the stochastic particle system for various interaction potentials. Finally, in Section 6, we present our conclusions.

2 Individual-based model

We consider a set of $N$ identical particles evolving according to the Langevin stochastic differential equation Eq. 1 in a domain $\Omega\subset\mathbb{R}^{d}$, with $d\leq 3$. We nondimensionalise time and space such that the diffusion coefficient $D=1$, and the volume of the domain $|\Omega|=1$. We suppose the interaction potential $u(r)$ is repulsive and short range, with range $\epsilon\ll 1$. The interaction potential of a system of $N$ particles is, assuming pairwise additivity, the sum of isolated pair interactions

[TABLE]

where $\vec{x}=({\bf x}_{1},\dots,{\bf x}_{N})$ is the $N$-particle position vector. The interaction force acting on the $i$th particle due to the other $N-1$ particles is given by

[TABLE]

Here forces are non-dimensionalized with the mobility (the inverse of the drag coefficient) so that we can talk about a force acting on a Brownian particle. Finally, we suppose that the initial positions ${\bf X}_{i}(0)$ are random and identically distributed.

The counterpart of Eq. 1 in probability space is the Fokker–Planck equation

[TABLE]

We proceed to reduce the dimensionality of the problem Eq. 6 by looking at the marginal density function of one particle (the first particle, say) given by

[TABLE]

The particle choice is unimportant since $P$ is invariant with respect to permutations of particle labels. Integrating Eq. 6a over ${\bf x}_{2},\dots,{\bf x}_{N}$ and applying the divergence theorem gives

[TABLE]

where ${\bf G}$ is the $d$-vectorial function

[TABLE]

and

[TABLE]

is the two-particle density function, which gives the joint probability density of particle 1 being at position ${\bf x}_{1}$ and particle 2 being at ${\bf x}_{2}$. An equation for $P_{2}$ can be written from Eq. 6a, but this then depends on $P_{3}$, the three-particle density function. This results in a hierarchy (the BBGKY hierarchy) of $N$ equations for the set of $n$-particle density functions ($n=1,\dots,N$), the last of which is Eq. 6a itself (since $P$ is the $N$-particle density function). In order to obtain a practical model, a common approach is to truncate this hierarchy at a certain level to obtain a closed system. In particular, closure approximations in which the $n$-particle density function $P_{n}$ is replaced by an expression involving lower density functions $P_{s}$, $s<n$, are commonly used. However, because of their phenomenological nature, they can often lead to errors in the resulting model. In the next section we present three such closure approximations and highlight the issues they encounter when dealing with short-range repulsive potentials. In Section 4 we present an alternative approach based on matched asymptotic expansions.

3 Closure approximations

3.1 Mean-field closure

The simplest and most common closure approximation is to assume that particles are not correlated at all in evaluating the interaction term $\bf G$, that is,

[TABLE]

Substituting Eq. 10 into Section 2 gives

[TABLE]

Combining this with the equation for $p$ in Eq. 8 gives equation Eq. 2 presented in the introduction, the mean-field approximation (MFA). However, one should keep in mind that Eq. 10 might not always be valid when using such model. In particular, when $u(r)$ is a short-range interaction potential, the dominant contribution to the integral Section 2 is when ${\bf x}_{1}$ is close to ${\bf x}_{2}$, and this is exactly the region in which the positions of particles are correlated. We note that the mean-field closure is often used implicitly with Eq. 2 written down directly rather than being derived from Eq. 6a [15, 22]. The reasoning goes as follows: if $p({\bf x},t)$ is the probability of finding a particle at ${\bf x}$, the force on a particle at ${\bf x}_{1}$ is given by multiplying the force due to another particle at ${\bf x}_{2}$ by the density of particles at ${\bf x}_{2}$ and integrating over all positions ${\bf x}_{2}$.

If we suppose the pair potential is short-ranged, we can approximate the integral in Section 3.1 to remove of the convolution term. In particular, we suppose that $u=O(r^{-(d+\delta)})$ for some $\delta>0$ as $r\to\infty$ and rewrite the potential as $u(r)=\tilde{u}(r/\epsilon)$ with $\epsilon\ll 1$. Introducing the change of variable ${\bf x}_{2}={\bf x}_{1}+\epsilon\tilde{\bf x}$ and expanding $p({\bf x}_{2},t)$ about ${\bf x}_{1}$ gives

[TABLE]

where we can extend the integral with respect to variable $\tilde{\bf x}$ to the whole space since the potential $\tilde{u}$ is localized near the origin and decays at infinity. Noting that the potential is a radial function, the leading-order term in the integral vanishes, and, after integrating by parts in the next term, we obtain

[TABLE]

Inserting Eq. 12 into Eq. 8, we find that the marginal density function satisfies the following nonlinear Fokker–Planck equation

[TABLE]

We will refer to (13), in which the mean-field closure has been combined with the assumption of a short-range potential, as the localized MFA or LMFA. As we shall see later, the MFA or LMFA are not defined for many commonly used short-range repulsive potentials. If the integral in Eq. 13b does not exist because of the behavior at infinity then this is an indication that the potentials is too long range for the localization performed above to be valid and the full MFA integral needs to be retained. However, if the integral in Eq. 13b diverges because of the behavior at the origin then the MFA itself will diverge, that is, the integral in Section 3.1 will not exist. Examples of the inappropriate use of MFA for such short range potentials exist in the literature [15].

3.2 Closure at the pair correlation function

A more elaborated closure is suggested by Felderhof [11]. His derivation considers a more general context of interacting Brownian particles suspended in a fluid, including hydrodynamic interactions. His analysis is valid for zero external force, ${\bf f}\equiv 0$, and is based on the thermodynamic limit (in which the number of particles $N$ and the system volume $V$ tend to infinity, with the number density $N/V=\rho_{0}$ fixed). Because of this, instead of working with probability densities, it is convenient to switch to number densities: $\rho({\bf x}_{1},t):=\rho_{0}p({\bf x}_{1},t)$ and $Q({\bf x}_{1},{\bf x}_{2},t):=\rho_{0}^{2}P_{2}({\bf x}_{1},{\bf x}_{2},t)$. In what follows we outline his derivation for hard spheres ignoring hydrodynamic interactions. The equation for the one-particle number density $\rho({\bf x}_{1},t)$ is

[TABLE]

where $r=\|{\bf x}_{1}-{\bf x}_{2}\|$ is the interparticle distance and $Q({\bf x}_{1},{\bf x}_{2},t)$ (the two-particle number density) satisfies, to lowest order in $\rho_{0}$,

[TABLE]

Equations Eq. 14 and Eq. 15 have the following time-independent equilibrium solutions

[TABLE]

where $\rho_{0}$ is constant and $g_{0}$ is the pair correlation function

[TABLE]

Felderhof then looks for a linearized solution around the equilibrium values Eq. 16, by making the ansatz

[TABLE]

and considering the deviations $\rho_{1}$ and $g_{1}$,

)+g1(x1,x2,t). Then, to terms linear in $\rho_{1}$ and $g_{1}$, Eq. 18 becomes

[TABLE]

Substituting in Eq. 14 and linearizing gives

[TABLE]

At this stage, Felderhof supposes that perturbations from equilibrium Eq. 16 are small so that $\nabla_{{\bf x}_{1}}\rho_{1}({\bf x}_{1},t)\approx\nabla_{{\bf x}_{2}}\rho_{1}({\bf x}_{2},t)$ and the pair correlation function is at its equilibrium value, that is, $g({\bf x}_{1},{\bf x}_{2},t)\approx g_{0}(r)$. Then Eq. 23 simplifies to

[TABLE]

Now, using Eq. 17 and expanding $\rho_{1}({\bf x}_{2},t)$ about ${\bf x}_{1}$ and keeping only the first non-vanishing term gives [11]

[TABLE]

where $\alpha_{u}=\int_{\mathbb{R}^{d}}\left(1-e^{-u(\epsilon\|{\bf x}\|)}\right)\mathrm{d}{\bf x}.$ Note that this is the evolution equation for the perturbation $\rho_{1}$ from the uniform equilibrium $\rho_{0}$ (valid with ${\bf f}=0$).

3.3 Kirkwood closure

As we will see later in the results section, the MFA can only provide an adequate approximation for relatively soft interaction potentials and low densities. An alternative closure approximation is based on the Kirkwood superposition approximation (KSA) [18], and consists of truncating the hierarchy at the two-particle density function. To this end, we consider the equation satisfied by $P_{2}({\bf x}_{1},{\bf x}_{2},t)$ by integrating the $N$-particle Fokker–Planck equation Eq. 6a over ${\bf x}_{3},\dots,{\bf x}_{N}$, applying the divergence theorem and relabelling particles as before to obtain

[TABLE]

where

[TABLE]

and

[TABLE]

is the three-particle density function. We note that in writing Eq. 29 we are using that $P_{3}$ is invariant to particle relabelling. The KSA then approximates the three-particle density function as

[TABLE]

where $p$ and $P_{2}$ are the one- and two-particle density functions, respectively. Inserting Eq. 31 into Eq. 29 one can then solve the coupled system Eq. 8 and Eq. 28 for $p$ and $P_{2}$.

The KSA closure has been the basis of many subsequent closure approximations, and it can be derived as the maximum entropy closure in the thermodynamic limit [26]. Because it is thought to be superior to the MFA at high densities, it has been used in several biological applications such as on-lattice birth-death-movement processes with size-exclusion [1] and off-lattice cell motility processes with soft interactions [20, 21]. Middleton, Fleck and Grima [21] consider a system of Brownian particles evolving according to Eq. 1 in one dimension interacting via a Morse potential (see Eq. 55e) and compare the KSA closure to the MFA closure and simulations of the stochastic system. Markham et al. [20] use the KSA on a more general individual-based model, where the random jumps of particles are not Gaussian but depend on the positions of all particles, resulting in multiplicative noise. Berlyand, Jabin and Potomkin [2] use a variant of the KSA closure, approximating either $P_{2}({\bf x}_{1},{\bf x}_{3},t)$ or $P_{2}({\bf x}_{2},{\bf x}_{3},t)$ in Eq. 31 by their corresponding mean-field approximation $P_{2}({\bf x}_{i},{\bf x}_{j},t)=p({\bf x}_{i},t)p({\bf x}_{j},t)$, for a system of interacting deterministic particles.

It is worth noting that the KSA model is computationally expensive and complicated to solve, especially if the interaction potential $u$ is short ranged, requiring a fine discretization. For example, in 3 dimensions one must solve a 6-dimensional problem which, once discretized, involves a full discretization matrix because of the convolution terms ${\bf G}$ and ${\bf G}_{2}$. As we shall see in Section 5.3, even in one physical dimensional the KSA model is rather complicated to solve.

4 Matched asymptotic expansions

In this section we consider an approach based on matched asymptotic expansions (MAE) to obtain a closed equation for the one-particle density $p$ that is valid for short-range interaction potentials, and that is computationally practical to solve even in two or three dimensions.

We go back to the evolution equation Eq. 8 for the one-particle density $p$. Assuming that the pair potential $u$ is localized near ${\bf x}_{1}$, we can determine $\bf G$ in Section 2 using MAE. To do so, we first must obtain an expression for $P_{2}$.

For low-concentration solutions with short-range interactions, three-particle (and higher) interactions are negligible compared to two-particle interactions: when two particles are close to each other, the probability of a third particle being nearby is so small that it can be ignored. Mathematically, this means that the two-particle probability density $P_{2}({\bf x}_{1},{\bf x}_{2},t)$ is governed by the dynamics of particles 1 and 2 only, independently of the remaining $N-2$ particles. In other words, the terms ${\bf G}_{2}$ in Eq. 28 are negligible and the equation for $P_{2}$ reduces to

[TABLE]

for $({\bf x}_{1},{\bf x}_{2})\in\Omega^{2}$, complemented with periodic boundary conditions on $\partial\Omega^{2}$. Note that in approximating $P_{2}$ in this way it is no longer true that $p({\bf x}_{1},t)=\int P_{2}({\bf x}_{1},{\bf x}_{2},t)\mathrm{d}{\bf x}_{2}$, so that $p$ and $P_{2}$ need to be solved for as a coupled system. This equation is basically Eq. 15 with an added external force. Essentially our MAE approach aims to solve Eqs. 14 and 15 systematically asymptotically rather than through Felderhof’s linearization and approximations (see Section 3.2).

4.1 Inner and outer regions

By assumption, the pair interaction potential $u(r)$ is negligible everywhere except when the interparticle distance $r$ is of order $\epsilon$. Therefore, we suppose that when two particles are far apart ($\|{\bf x}_{1}-{\bf x}_{2}\|\gg\epsilon$) they are independent, whereas when they are close to each other ($\|{\bf x}_{1}-{\bf x}_{2}\|\sim\epsilon$) they are correlated. We designate these two regions of configuration space the outer region and inner region, respectively.

In the outer region we define $P_{\text{o}}({\bf x}_{1},{\bf x}_{2},t)=P_{2}({\bf x}_{1},{\bf x}_{2},t)$. By independence, we have that111Independence only tells us that $P_{\text{o}}({\bf x}_{1},{\bf x}_{2},t)\sim q({\bf x}_{1},t)q({\bf x}_{2},t)$ for some function $q$ but the normalization condition on $P$ implies $p=q+O(\epsilon)$.

[TABLE]

for some function $P_{\text{o}}^{(1)}$.

In the inner region, we set ${\bf x}_{1}=\tilde{\bf x}_{1}$ and ${\bf x}_{2}=\tilde{\bf x}_{1}+\epsilon\tilde{\bf x}$ and define $\tilde{P}(\tilde{\bf x}_{1},\tilde{\bf x},t)=P_{2}({\bf x}_{1},{\bf x}_{2},t)$ and $\tilde{u}(\|\tilde{\bf x}\|)=u(\|{\bf x}_{1}-{\bf x}_{2}\|)$. Rewriting Eq. 34 in terms of the inner coordinates gives

[TABLE]

The inner solution $\tilde{P}$ must match with the outer solution $P_{\text{o}}$ as $\|\tilde{\bf x}\|\rightarrow\infty$. Expanding $P_{\text{o}}$ in terms of the inner variables gives (omitting the time variable for ease of notation)

[TABLE]

We look for a solution of Eq. 38 matching with Eq. 39 as $\|\tilde{\bf x}\|\rightarrow\infty$ of the form $\tilde{P}\sim\tilde{P}^{(0)}+\epsilon\tilde{P}^{(1)}+\cdots$. The leading-order inner problem is

[TABLE]

with solution

[TABLE]

The $O(\epsilon)$ problem reads

[TABLE]

Using Eq. 41, we can rearrange Eq. 42a to give

[TABLE]

Solving Eq. 43 together with Eq. 42b gives

[TABLE]

Thus we find that the inner region solution is, to $O(\epsilon)$,

[TABLE]

4.2 Interaction integral

Now we go back to the interaction integral ${\bf G}({\bf x}_{1})$ in Section 2. Because of the short-range nature of the potential $u$, the main contribution to this integral is from the inner region. Therefore, we will use the inner solution Eq. 45 to evaluate it.

First we split the integration volume $\Omega$ for ${\bf x}_{2}$ into the inner and the outer regions defined in the previous section. Although there is no sharp boundary between the inner and outer regions, it is convenient to introduce an intermediate radius $\delta$, with $\epsilon\ll\delta\ll 1$, which divides the two regions. Then the inner region is $\Omega_{\text{i}}({\bf x}_{1})=\{{\bf x}_{2}\in\Omega\ :\ \|{\bf x}_{2}-{\bf x}_{1}\|<\delta\}$ and the outer region is the complimentary set $\Omega_{\text{o}}({\bf x}_{1})=\Omega\setminus\Omega_{\text{i}}({\bf x}_{1})$. The dominant contribution to Section 2 is then

[TABLE]

The first and third terms of the integral vanish using the divergence theorem and that $\tilde{u}$ is a radial function of $\|\tilde{\bf x}\|$. Integration by parts on the second component gives

[TABLE]

where we have used that $e^{-\tilde{u}(\|\tilde{\bf x}\|)}\approx 1$ at $\|\tilde{\bf x}\|=\delta/\epsilon$. Finally, rewriting the first term above as a volume integral,222The volume of a $d$-dimensional sphere of radius $\delta/\epsilon$ is equal to $2\pi(d-1)/d(\delta/\epsilon)^{d}$. Section 4.2 becomes

[TABLE]

Since $1-e^{-\tilde{u}(\|\tilde{\bf x}\|)}$ decays at infinity, we can extend the domain of integration to the entire $\mathbb{R}^{d}$ introducing only lower order errors. Therefore we can write

[TABLE]

with

[TABLE]

Note that we obtain the same coefficient $\alpha_{u}$ that appeared in equation Eq. 25 using the closure at the pair correlation function.

4.3 Reduced Fokker–Planck equation for soft spheres

Combining Eq. 48 with Eq. 8 we find that, to $O(\epsilon^{d})$,

[TABLE]

Therefore, we find that the MAE method yields the same type of equation as the LMFA (see Eq. 13), but with a different coefficient in the nonlinear diffusion term. Expanding the exponential in Eq. 49b for small $u$ gives

[TABLE]

that is, the LMFA closure is the leading contribution of the potential, provided that $u(\epsilon\|{\bf x}\|)$ is small. However, as we will see in Section 5, this is not always true. We also see that the equation obtained by Felderhof [11], Eq. 25 is the linearized version of Eq. 49a after taking $N\to\infty$ and setting the external force to zero.

The coefficient $\alpha_{u}$ can be related to basic concepts from statistical mechanics. Namely, the integrand in Eq. 49b is the negative of the total correlation function $h(r)=g(r)-1$, where $g(r)=\exp[-u(r)]$ is the low-density limit of the pair correlation function, the so-called Boltzmann factor of the pair potential [12, §2.6].

4.4 Two-particle density function via MAE

In contract to the MFA, using MAEs we have obtained an approximation of the two-particle density function $P_{2}({\bf x}_{1},{\bf x}_{2},t)$ in two regions of the configuration space, the outer and the inner regions, defined according to the separation between the two particles. It is convenient to have a uniformly valid expansion in the whole space–so that for example we can plot and compare it against simulations of the stochastic particle system. This can be done by constructing a so-called composite expansion, consisting of the inner expansion Eq. 45 plus the outer expansion Eq. 35 minus the common part [14, Chapter 5]. The result is

[TABLE]

5 Results

5.1 The hard-sphere potential

So far we have assumed that the set of Brownian particles interact via a soft pair potential $u(r)$. However, the resulting reduced Fokker–Planck equation Eq. 49 obtained via MAE can be used to model a system of hard-core interacting particles of diameter $\epsilon$. In particular, the model Eq. 49 for soft spheres has exactly the same structure as the counterpart model for hard spheres derived in [4]. Inserting the hard-sphere pair potential

[TABLE]

into the nonlinear diffusion coefficient in Eq. 49b gives

[TABLE]

that is, $\alpha_{\text{HS}}=2$ for $d=1$, $\alpha_{\text{HS}}=\pi$ for $d=2$, and $\alpha=4\pi/3$ for $d=3$. This is in agreement with our previous work where the reduced model was specifically derived for a system of hard spheres [4, 5]. For hard spheres the configuration space has holes due to the illegal configurations (associated with infinite energy) and the derivation is slightly different. We note that the MFA for hard spheres does not work and that, in particular, the coefficient $\overline{\alpha}_{\text{HS}}$ in Eq. 13b is not defined.

Because the MAE models for soft and hard spheres coincide, for every potential $u$ we can find an effective hard-sphere diameter $\epsilon_{\text{eff}}$ such that the continuum models of both systems – the soft-sphere system with pair potential $u(r)$ and the hard-sphere system with diameter $\epsilon_{\text{eff}}$ – are equivalent. In other words, a characterization of a given soft potential $u$ is to find $\epsilon_{\text{eff}}$ such that

[TABLE]

where $\alpha_{u}$ is given in Eq. 49b and $\alpha_{\text{HS}}$ given in Eq. 52. Rearranging, we find that

[TABLE]

This idea of finding the effective hard-sphere diameter associated to a soft-sphere system, known as the effective hard sphere diameter method, has been widely used to calculate both equilibrium and transport properties [23, §9.3.3]. The reason this method is appealing is that it allows us to “translate” a general system of interacting soft spheres (whose properties may have not been studied before) to the widely studied hard-sphere system, on which most theories are based. Moreover, the derived model Eq. 49 implies that, as far as the population-level dynamics are concerned, soft interactions may be incorporated into the effective hard particle model by adjusting the hard-sphere diameter with Eq. 54. This is provided $\alpha_{u}$ in Eq. 49b is well-defined and positive. In contrast, if $\alpha_{u}$ does not exist (or is negative), the MAE approach breaks down (or may become unstable) for the given pair potential $u$. Even this is instructive: it means that the potential is not decaying at infinity fast enough to be incorporated into a population-level equation of the form Eq. 49a and that the MFA is preferrable.

Another application of the effective hard-sphere diameter is to provide an effective volume fraction for the system of soft spheres: using $\epsilon_{\text{eff}}$, one can define the effective volume fraction of soft spheres and use it to check whether the “low-volume fraction” condition holds.

5.2 Comparison between MAE and LMFA

[TABLE]

To model some physical systems it is convenient to incorporate an attractive part to the repulsive pair potential. The most common situation is that particles repel each other in the short range and attract each other in a longer range. For example, the SS potential in Eq. 55a may be generalized to power-law repulsive–attractive potentials of the form $u(r)=(\epsilon/r)^{a}-C(\epsilon/r)^{b}$, with $a>b$ and $C>0$, known as the Mie potentials. The most famous example of this class of potentials is the Lennard–Jones potential,

[TABLE]

Another common repulsive–attractive potential is the Morse (MO) potential

[TABLE]

where $C$ and $l$ are, respectively, the relative strength and relative lengthscale of the repulsion to attraction. The most relevant situations for biological applications are given for $C>1$ and $l<1$, which correspond to short-range repulsion and weak long-range attraction [7].

Fig. 1 shows examples of all the interaction potentials above.

A simple way to compare between approaches for short-range interaction potentials is to consider the nonlinear diffusion term using either MAE or LMFA. Their respective coefficients $\alpha_{u}$ in Eq. 49b and $\overline{\alpha}_{u}$ in Eq. 13b for the potentials above are shown in Table 1. Since the behaviour of the integrals in $\alpha_{u}$ and $\overline{\alpha}_{u}$ at infinity is the same, the LMFA fails for any case in which the MAE fails. LMFA may fail also because of a singularity at the origin for which MAE is valid, as seen in some examples of Table 1. This is because the strongly repulsive short-range part of the potential results in correlations which violate the MFA that particles may be considered independent. Moreover, there are considerable discrepancies between $\alpha_{u}$ and $\overline{\alpha}_{u}$ for the cases for which the latter is defined.

The coefficient $\overline{\alpha}_{u}$ is undefined for inverse-power potentials such as SS and LJ, since the integral in Eq. 13b is either singular at the origin or at infinity for all possible powers. Therefore, the LMFA is not valid for these potentials. The MAE coefficient $\alpha_{u}$ exists for the SS potential in Eq. 55a for $\nu>d$, but is undefined for $\nu\leq d$. Therefore, we find that MAE is not valid for the Coulomb interaction [$\nu=1$ in Eq. 55a]. The interpretation is that the Coulomb interaction does not decay sufficiently quickly at infinity and hence the inner region spans the whole configuration space (so there is no outer region where the integral ${\bf G}({\bf x}_{1})$ is negligible as assumed in the MAE derivation).

Fig. 2 shows the variation of $\alpha_{\text{SS}}$ against the hardness parameter $\nu$. As expected, $\lim_{\nu\to\infty}\alpha_{\text{SS}}=\alpha_{\text{HS}}$ since the HS potential is the limiting case of the SS potential for $\nu\to\infty$. Also note that the strength of the nonlinear diffusion term, parametrized by $\alpha_{\text{SS}}$, decreases as the hardness $\nu$ increases. The effect of the softness parameter $1/\nu$ on the diffusion and other transport coefficients was studied in [13]; they found that for $\nu\geq 72$ the soft sphere system behaved essentially as a hard sphere system in molecular-dynamics simulations. We find that the relative error between $\alpha_{\text{SS}}$ and $\alpha_{\text{HS}}$ for $\nu=72$ in three dimensions is 2.6%.

The Morse potential Eq. 55e has well-defined coefficients $\alpha_{\text{MO}}$ and $\overline{\alpha}_{\text{MO}}$ for $C$ and $l$ positive. However, they differ substantially depending on the relative strength $C$ and relative lengthscale $l$ of the repulsion to the attraction, see Fig. 3. We note that both coefficients $\alpha_{\text{MO}}$ and $\overline{\alpha}_{\text{MO}}$ may be negative for some parameter values, which implies, from Eq. 49a, that the nonlinear component of the diffusion coefficient becomes negative. When this occurs, the system enters a so-called catastrophic regime and the particles collapse to a point as $N\to\infty$ [10].

Evaluation of $\alpha_{u}$ and $\overline{\alpha}_{u}$ provides a straightforward way to determine the type of dynamics the interacting particle system has. Initially it is not so easy to discern whether the potential is “short-range enough” and, as a result, which method is more suitable to obtain a reduced population-level model. The nonlinear diffusion coefficient Eq. 49a, $\alpha_{u}(N-1)\epsilon^{d}p$, gives an idea of the strength of the interaction and whether the MAE will be the appropriate method (in particular, one can compute the effective volume fraction as described in Section 5.1 and ascertain whether the low-volume fraction assumption holds). Regarding the MFA closure, one should keep in mind that $\overline{\alpha}_{u}$ not being defined does not necessarily mean that standard MFA fails, as mentioned in Section 4.3; it could be that the short-range assumption used to obtain the nonlinear diffusion equation Eq. 13 does not hold and that, instead, one should keep the original MFA integro-differential model Eq. 2. Conversely, if $\overline{\alpha}_{u}$ is not defined due to the singular behavior of the potential at the origin (which is the case for all those seen in Table 1), this is an indication that the convolution in the original MFA will not exist.

5.3 Comparison with the particle-level model

In this section we compare the macroscopic models obtained via the closures MFA and KSA, and via the MAE method we have introduced to each other as well as to numerical simulations of the stochastic particle system.

We use the open-source C++ library Aboria [24, 25] to perform the particle-level simulations. The overdamped Langevin equation Eq. 1 is integrated using the Euler-Maruyama method and a constant timestep $\Delta t$, leading to an explicit update step for each particle given by

[TABLE]

where $\Delta W^{m}$ is a $d$-dimensional normally distributed random variable with zero mean and unit variance. We choose the timestep $\Delta t$ so that the results are converged (that is, there is no change in the results for smaller timesteps). For all simulations in this paper a timestep of $\Delta t=(0.1\epsilon)^{2}/2D$ was sufficient for convergence, leading to an average diffusion step size of $0.1\epsilon$.

A naive implementation of the particle force interaction term over all particle pairs would lead to a large number ($N^{2}$) of potential pair interactions to perform. To improve the efficiency of the code, we take advantage of the compact nature of the potentials and restrict particle interaction to those pairs that are within a certain cutoff length $c<6\epsilon$. All particle pairs $i,j$ separated by a distance greater than this cutoff will implicitly be given a pair potential of $u(\|{\bf X}_{i}-{\bf X}_{j}\|)=0$.

In order to compare the particle-level models with the PDE models, we perform $R$ independent realizations and output the positions of all $NR$ particles at a set of output time points. A histogram of the positions is calculated and then scaled to produce a discretized density function that can be compared with the PDE models. To generate the two-particle density function, we create a two-dimensional histogram of the positions of each particle pair $(i,j)$ and scale it accordingly to produce a two-particle density.

First, we consider a one-dimensional problem in $\Omega=[0,1]$ with periodic boundary conditions. We compare estimates of the one-particle density $p(x_{1},t)$ as well as the two-particle density $P_{2}(x_{1},x_{2},t)$ obtained from simulating Eq. 63 to solutions of the same quantities using the KSA, MFA or MAE models. We start with a set of parameters in which the potential is sufficiently long range that the system is not really in the low-density limit with the number of particles we use, and then consider an example with a strongly repulsive short-ranged potential. In all the examples we set the external drift to zero.

We begin by presenting the models for the three approaches for $d=1$. The MFA reads, combining Eq. 8 and Section 3.1,

[TABLE]

where $p=p(x_{1},t)$ unless explicitly written and $f_{u}(x_{1},x_{2}):=\frac{\partial}{\partial x_{1}}u(|x_{1}-x_{2}|)$. To solve Eq. 64 we use a second-order accurate finite-difference approximation with $M$ grid points in space and the method of lines in time with the inbuilt Matlab ode solver ode15s. To evaluate the integral in Eq. 64, we use the periodic trapezoid rule (which converges exponentially fast for smooth integrands [27, Chapter 12]). To avoid evaluating the interaction potential at zero, we shift the grid for $x_{2}$ by $h/2$, where $h=1/M$ is the grid spacing. The density $p(x_{2},t)$ is approximated from $p(x_{1},t)$ using $p(x_{2},t)\approx(p(x_{2}-h/2,t)+p(x_{2}+h/2,t))/2$. Because of the convolution term, the discretization matrix of Eq. 64 is full, making the numerical solution of the MFA computationally expensive. An alternative would be to use a nonuniform mesh with more points near the diagonal or an adaptive grid scheme such as that presented by Carrillo and Moll [8], which uses a transport map between the uniform density and the unknown density $p$ such that more grid points are placed where the density is higher.

The KSA closure model is the coupled system for $p(x_{1},t)$ and $P_{2}(x_{1},x_{2},t)$, given by

[TABLE]

where $p=p(x_{1},t)$ and $P_{2}=P_{2}(x_{1},x_{2},t)$ unless explicitly written. We solve in the computational spatial domain $\Omega^{2}$ and therefore the KSA model requires a two-dimensional grid with full matrices, making it very computationally expensive. As with the MFA, we shift the grid for $x_{2}$ by $h/2$ to avoid any issues with $f_{u}(x_{1},x_{2})$ when the interaction potential $u$ is singular at the origin. In the first integral of Eq. 65b, the coordinate $x_{3}$ is evaluated on the same grid as $x_{2}$, and in the second integral it is evaluated on the same grid as $x_{1}$. An alternative implementation of the KSA system is to solve for $P_{2}$ only and evaluate the one-particle density as $p(x_{1},t)=\int_{\Omega}P_{2}(x_{1},x_{2},t)\,\mathrm{d}x_{2}$ (replacing Eq. 65a) [21]. While this is true in the infinite BBGKY hierarchy, we note that once the KSA closure Eq. 31 is adopted, the resulting model is not necessarily equal to the KSA system Eq. 65.

Finally, the reduced model from the MAE method reads

[TABLE]

where $\alpha_{u}=2\int_{0}^{\infty}(1-e^{-u(\epsilon r)})\,\mathrm{d}r$. Noting that the right-hand side in Eq. 66 can be written as $[p+\frac{1}{2}\alpha_{u}(N-1)\epsilon^{d}p^{2}]_{xx}$, the numerical implementation is straightforward and, since the discretization matrix is banded, very efficient.

In the first example, we use the exponential potential $u_{\text{EX}}$ Eq. 55b with $\epsilon=0.05$ in a system in $N=16$ particles. The particles are initially distributed according to $p(0)=0.5[\tanh(\beta(x-\theta))+\tanh(\beta(1-\theta-x))]$, with $\beta=30$ and $\theta=0.2$, and we let them diffuse until time $T_{f}=0.02$. The evolution in time of the one-particle density obtained with the various methods is shown in Fig. 4. We observe very good agreement between the stochastic simulations and the MFA and KSA models, whereas the MAE model slightly underestimates the diffusion strength. This can be explained since, for the chosen values of $\epsilon$ and $N$, the potential is not short-ranged and the assumptions of the MAE model break down. In particular, using Eq. 54 we find that the effective hard-sphere radius is $\epsilon_{\text{eff}}=0.04$ and hence the effective volume fraction is 0.64. Similarly, the solution obtained via the LMFA (shown at $T_{f}$ in Fig. 4(b)) differs noticeably to the MFA solution.

Figure 5 shows the corresponding two-particle density function at the final time $T_{f}$. For the KSA, $P_{2}$ is solved for as discussed above. For the MFA and MAE models, $P_{2}$ is given by Eq. 10 and Eq. 50 respectively. As before, $x_{2}$ is shifted by half the grid size and $p(x_{2},t)$ is approximated from $p(x_{1},t)$ using the centered average. The correlation between particles can be seen by the drop in probability at the diagonal $x_{1}=x_{2}$ in the simulations as well as in the KSA and MAE plots. Conversely, the MFA misses the correlation between particles, as expected from the ansatz $P_{2}(x_{1},x_{2},t)=p(x_{1},t)p(x_{2},t)$. The differences in $P_{2}$ between the models can be clearly seen in Fig. 8(a), which shows a plot of $P_{2}(x_{1},x_{2},T_{f})$ along $x_{2}=0.5$.

Next we consider an example using a more repulsive interaction potential. In particular, we consider a smoothed version of the Yukawa potential Eq. 55c, namely $u(r)=(\epsilon/\sqrt{r^{2}+\delta^{2}})\exp^{-r/\epsilon}$ with $\epsilon=0.01$ and $\delta=0.002$. We choose to smooth the potential so that particles can still swap positions and so that we can use the closure models KSA and MFA (in one dimension the singularity at zero poses problems in the convolution terms, see Table 1). We run a simulation with $N=20$ particles up to time $T_{f}=0.02$. For these parameters, the effective hard-sphere radius is $\epsilon_{\text{eff}}=0.009$ and the effective volume fraction is 0.18. Thus we expect the MAE model to perform better than in the previous example where the potential was not short range. The one-particle density plots are shown in Fig. 6, and the results for the two-particle density function are shown in Fig. 7. We see that the MFA closure solution diffuses substantially faster than the KSA closure solution, whereas the MAE model spreads slightly slower than the KSA one (Fig. 6). The particle-level simulations match the KSA solution closely over all times. The two-particle density functions for the particle simulation, the KSA and the MAE solution are indistinguishable at time $T_{f}=0.02$, whereas the MFA misses the drop in probability for closely spaced particles, as expected due to its assumption of no particle correlations (Figs. 7 and 8(b)).

Generally, we see that MAE is good for strongly repulsive interactions (such as in the example in Figs. 6 and 7) while MFA is better for softer or longer-range interactions (Figs. 4 and 5). The KSA provides a good approximation in all cases. When the two-particle density is required in a system with long-range interactions, our results may one lead to think that one should use the KSA, since the MFA does not capture the correlation in $P_{2}$. However, we suggest that a slight modification in the MFA can provide a reasonable approximation to $P_{2}$ also. Specifically, one could still use the standard MFA closure $P_{2}(x_{1},x_{2},t)=p(x_{1},t)p(x_{2},t)$ to compute the one-particle density, but then use $P_{2}(x_{1},x_{2},t)\approx p(x_{1},t)p(x_{2},t)\exp(-u(|x_{1}-x_{2}|))/C$ as an approximation of the two-particle density, where $C=\iint\exp(-u(|x_{1}-x_{2}|))\mathrm{d}x_{1}\mathrm{d}x_{2}$ is a normalization constant. This modification of $P_{2}$, using either the MFA or the LMFA to compute $p$, is shown as MFAe and LMFAe respectively in Figs. 5, 7 and 8. We see that MFAe and LMFAe provide a better approximation of $P_{2}$ than MFA and MAE in the first example with longer-range interactions (Fig. 8(a)), whereas in the second example (with strong repulsion) the three approximations MFA, MFAe, and LFMAe are poor (Fig. 8(b)).

Finally, we consider a two-dimensional example in $\Omega=[0,1]^{2}$. We choose a system of $N=400$ Yukawa-interacting particles, with interaction potential $u_{\text{YU}}$ Eq. 55c and $\epsilon=0.01$, initially distributed according to a normal distribution in $x$ with mean $0.5$ and standard deviation $\sigma=0.05$. The effective hard-sphere radius is $\epsilon_{\text{eff}}=0.0112$, giving an effective volume fraction of 0.04. Being in two dimensions, solving the KSA model would require solving a system of $M^{4}$ equations, where $M$ is the number of grid points in one direction. Because of the short-range nature of the pair potential $u_{\text{YU}}$, we require a large number of points $M$ to resolve the interaction near the origin. As a result, solving the KSA model for this system becomes computationally impractical, and we compare the MFA and the MAE methods only. We also solve the interaction-free model ($\epsilon=0$) for reference. We plot the comparison at times $T_{f}/2$ and $T_{f}$ in Fig. 9. We observe a very good agreement between the stochastic simulations of the particle system and the MAE model, whereas the MFA model overestimates the diffusion strength. Because of the short-range nature of the potential, we find no noticeable difference between the solution of the MFA (2) and the localized version LMFA Eq. 13, which is computationally much easier to solve (the two curves would lie on top of each other in Fig. 9, the norm of the relative error is of order $10^{-4}$). As expected, the interaction-free case spreads the slowest of all models since the nonlinear diffusion term is set to zero.

6 Discussion and conclusions

We have studied a system of Brownian particles interacting via a short-range repulsive potential $u$, and have discussed several ways to obtain a population-level model for the one-particle density. In particular, we have considered two common closure approximations and presented an alternative method based on matched asymptotic expansions (MAE). The MAE method has the advantage that it is systematic and works well for very short-ranged potentials, especially singular potentials for which common closure approximations can lead to ill-posed models.333These models are sometimes used regardless, with a numerical discretisation providing an ad hoc regularisation of the convolution integral [15]. The MAE result is a nonlinear diffusion equation similar to our previous work for hard-spheres [4], with the coefficient of the nonlinear term $\alpha_{u}$ depending on the potential through (49b).

We have performed Monte Carlo simulations of the stochastic particle system in one and two dimensions, and compared the results with the solution of the MAE model and two common closure approximations: the mean-field approximation (MFA) and the Kirkwood Superposition Approximation (KSA). While MFA closes the system at the level of the one-particle density $p$, the KSA closes it at the level of the two-particle density $P_{2}$. We have tested the models in examples with long- and short-range interactions.

We found that the KSA agreed well with the stochastic simulations in both scenarios, but we could only use it in the one-dimensional examples due to its high computational cost. This is because the discretisation of the convolution term yields a full matrix, making the method impractical to use in two or three dimensions, especially for strongly repulsive potentials that require a very fine mesh in the region where two particles are in close proximity. This is also true, but to a lesser extent, for the MFA model, which captured well the behavior of the system with a long-range potential but was outperformed by MAE in examples with a short-range potential. The MAE method results in a nonlinear diffusion model (which is thus local, with banded discretization matrix) for $p$, making it straightforward to solve.

Recently [21] argued in favour of the KSA because it gave the two-particle density as well as the one-particle density, and could therefore capture correlations in particle positions, which the MFA cannot. Our MAE method also gives an approximation for $P_{2}$, and it successfully captures the low likelihood of finding particles close to each other when there are strong short-range repulsions. We emphasize again that it does so at a fraction of the cost of KSA, so that the MAE method becomes particularly suited for problems in two or three dimensions where the KSA or similar higher-order closures are impractical. We noted also that the MFA can be extended simply (to what we called MFAe) to produce an approximation for $P_{2}$, so that it too becomes a viable alternative if the interactions are longer range.

Stochastic simulations of system of repulsive soft spheres are generally regarded as simpler than the simulation of hard spheres since one avoids the issue of how to approximate the collision between two particles. Instead, soft sphere simulations involve summing up the contribution to the interaction force of all the neighbors at each time-step and adding it as a drift term to the Brownian motion. However, for very repulsive potentials this needs to be done very carefully, since if the time-step is not small enough the repulsive part of the potential is not resolved correctly and easily missed. If this happened, the low-density diagonal in the two-particle density plots (see for example Fig. 7) would either not be well resolved or not be there at all. It is therefore important to do a proper convergence study for the stochastic simulations in order to decide on the appropriate time-step. We note in this respect that our coefficient $\alpha_{u}$ provides a natural way to determine the radius of the equivalent hard-sphere for any short range potential.

We have seen that the MAE method works well for repulsive short-range potentials, while the MFA provides a good approximation for long-range interactions. A natural question is to ask what to do for potentials with both characteristics, namely those which are very singular at the origin but with fat tails at infinity. This work provides a possible route to deal with such potentials: to combine the MAE and MFA methods. Specifically, one could break the potential into two parts and deal with each of them separately. The result would be an equation of the type Eq. 2 with the nonlocal convolution term due to the long-range component of the potential, and a nonlinear diffusion due to the short-range component.