Diffusion of particles with short-range interactions
Maria Bruna S. Jonathan Chapman Martin Robinson


Definition of the relativistic geoid in terms of isochronometric surfaces
    Dennis Philipp, Volker Perlick, Dirk Puetzfeld, Eva Hackmann, and Claus Lämmerzahl
    
======================================================================================





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.




diffusion, soft spheres, closure approximations, particle systems, matched asymptotic expansions




35Q84, 60J70, 82C31





§ 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

    X_i(t) =  √(2D)  W_i(t) +  f  (  X_i (t) )  t  - ∑_j  i∇_ x_i u( X_i(t) -  X_j(t))  t,

for i = 1, …, N,
where X_i (t) ∈Ω⊂ℝ^d is the position of the ith particle at time t,  D is the diffusion constant, W_i(t) denotes a d-dimensional Brownian motion, and f: Ω→ℝ^d is an external force. 

Despite the conceptual simplicity of the stochastic model <ref>, 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 <ref>. 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( x,t) of finding a given particle at position x at time t. The standard mean-field approximation (MFA)  procedure applied to the microscopic model <ref> gives

    ∂ p/∂ t  =  ∇_ x [ D  ∇_ x p -  f ( x) p  + (N-1) p ∇_ x (u∗ p) ],
 
where  u ∗ p = ∫ u( x -  y) p( y) 
  y. The main assumption in writing down <ref> 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 (ρ = N p) 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→∞, N/V = ρ_0. This is known as the thermodynamic
limit <cit.>, and results in equation
(<ref>) for ρ without the factor N-1.

While the mean-field approximation <ref> is convenient and leads to an accurate description for long-range interactions, such as with Coulomb interactions <cit.>, it fails when considering relatively strong repulsive short-range potentials. Sometimes <ref> 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 <ref> does not make sense with hard-core repulsive interactions, which are commonly used to model excluded-volume effects in biological and social contexts <cit.>. 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 <cit.>. One can derive an analogous continuous limit to <ref> using Taylor expansions <cit.>. However, it turns out that with identical particles the interaction terms (analogous to the convolution term in <ref>) cancel out unless the external force 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 <cit.> 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” <cit.> 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 nth and (n+1)th density functions in the BBGKY hierarchy (see <ref>). 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 <cit.>. 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 <cit.>. 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 <cit.>. 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 <ref> 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 ϵ arising from the short-range potential and the typical separation between particles. The result is a nonlinear advection-diffusion equation of the form 

    ∂ p/∂ t = ∇_  x· [ D ∇_  x  p -   f(  x) p + α_u ϵ^d (N-1)  p ∇_  x  p    ],

where the coefficient α_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 <ref> 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 <ref>. In <ref> we discuss three common closure approximations which reduce the Fokker–Planck equation to  a population-level PDE.  In <ref> we present our alternative approach to closure based on MAE, and derive equation <ref>. In <ref> 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 <ref>, we present our conclusions. 





§ INDIVIDUAL-BASED MODEL
 


We consider a set of N identical particles evolving according to the Langevin stochastic differential equation <ref> in a domain Ω⊂ℝ^d, with d ≤ 3. We nondimensionalise time and space such that the diffusion coefficient D = 1, and the volume of the domain  | Ω | = 1. We suppose the interaction potential u(r) is repulsive and short range, with range ϵ≪ 1. 
The interaction potential of a system of N particles is, assuming pairwise additivity, the sum of isolated pair interactions

    U(x⃗) =  ∑_1≤ i<j≤ N u ( x_i -  x_j),

where x⃗ = ( x_1, …,  x_N) is the N-particle position vector. The interaction force acting on the ith particle due to the other N-1 particles is given by

    g_i (x⃗) = - ∇_ x_i U (x⃗) = - ∑_j  i∇_ x_i u( x_i -  x_j).

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 X_i(0) are random and identically distributed. 


The counterpart of <ref> in probability space is the Fokker–Planck equation



    ∂ P/∂ t (x⃗, t) = ∇⃗_x⃗· [  ∇⃗_x⃗ P - F⃗ (x⃗) P +  ∇_x⃗ U (x⃗) P  ]     in    Ω^N,

where P(x⃗, t) is the joint probability density function of the N particles being at positions x⃗ = ( x_1, …,  x_N)∈Ω^N at time t and F⃗ (x⃗) = ( f_1( x_1), …,  f_N( x_N) ). Since we want to conserve the number of particles, on the domain boundaries ∂Ω ^N we require either no-flux or periodic boundary conditions. Throughout this work we use the latter. Accordingly, the potential u will be a periodic function in Ω. 
The initial condition is 

    P(x⃗, 0) = P_0(x⃗),

with P_0 invariant to permutations of the particle labels. 


We proceed to reduce the dimensionality of the problem <ref> by looking at the marginal density function of one particle (the first particle, say) given by  

    p( x_1,t) = ∫_Ω^N-1 P(x⃗, t)   x_2 ⋯ x_N.

The particle choice is unimportant since P is invariant with respect to permutations of particle labels. Integrating <ref> over x_2, …,  x_N and applying the divergence theorem gives

    ∂ p/∂ t ( x_1,t) = ∇_  x_1·[ ∇_  x_1   p  -  f( x_1) p  +   G( x_1,t)  ],

where G is the d-vectorial function

    G( x_1,t)    = ∫_Ω^N-1  P(  x_1,  x_2, …,   x_N, t) ∑_j=2^N ∇_ x_1 u ( x_1 -  x_j)     x_2 ⋯ x_N 
       = (N-1) ∫_Ω P_2(  x_1,  x_2, t) ∇_ x_1 u( x_1 -  x_2 )     x_2,

and

    P_2(  x_1,  x_2, t) = ∫_Ω^N-2 P(x⃗, t)   x_3 ⋯ x_N

is the two-particle density function, which gives the joint probability density of particle 1 being at position x_1 and particle 2 being at x_2. An equation for P_2 can be written from <ref>, 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, …, N), the last of which is <ref> 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 <ref> we present an alternative approach based on matched asymptotic expansions.




§ CLOSURE APPROXIMATIONS







 §.§ 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 G, that is,

    P_2(  x_1,  x_2, t) = p(  x_1, t) p(  x_2, t).

Substituting <ref> into <ref> gives

    G ( x_1,t)= (N-1 )   p(  x_1,t) ∫_Ω p(  x_2, t) ∇_ x_1 u( x_1 -  x_2)     x_2.

Combining this with the equation for p in <ref> gives equation <ref> presented in the introduction, the mean-field approximation (MFA). However, one should keep in mind that <ref> 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 <ref> is when x_1 is close to 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 <ref> written down directly rather than being derived from <ref> <cit.>. The reasoning goes as follows: if p( x,t) is the probability of finding a particle at x, the force on a particle at x_1 is given by multiplying the force due to another particle at x_2 by the density of particles at x_2 and integrating over all positions x_2. 

If we suppose the pair potential is short-ranged, we can approximate the integral in <ref> to remove of the convolution term. In particular, we suppose that u = O(r^-(d+δ)) for some δ>0 as r→∞ and rewrite the potential as u(r) = ũ(r/ϵ) with ϵ≪ 1.  Introducing the change of variable x_2 =  x_1 + ϵx̃ and expanding p( x_2,t) about x_1 gives

    G (  x_1,t) = - ϵ^d-1 (N-1 )   p(   x_1,t) ∫_ℝ^d[ p(   x_1, t) + ϵx̃·∇_ x_1 p(   x_1, t)  ]   ∇_x̃ ũ( x̃ )  x̃ + ⋯,

where we can extend the integral with respect to variable x̃ to the whole space since the potential ũ 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

    G (  x_1,t) ∼ϵ^d (N-1 )   p(   x_1,t) ∇_x̃_1 p(   x_1,t)   ∫_ℝ^dũ( x̃)  x̃ + O(ϵ^d+δ).

Inserting <ref> into <ref>, we find that the marginal density function satisfies the following nonlinear Fokker–Planck equation



    ∂ p/∂ t = ∇_  x_1·{ [1 + α_u (N-1) ϵ^d p ] ∇_  x_1  p -    f(  x_1) p },

where the nonlinear coefficient is given by

    α_u =  ∫_ℝ^d   u( ϵ x )   x.


We will refer to (<ref>), 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 <ref> 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 <ref> diverges because of the behavior at the origin then the MFA itself will diverge, that is, the integral in <ref> will not exist. Examples of the inappropriate use of MFA for such short range potentials exist in the literature <cit.>. 




 §.§ Closure at the pair correlation function
 

A more elaborated closure is suggested by Felderhof <cit.>. 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, f≡ 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 = ρ_0 fixed). Because of this, instead of working with probability densities, it is convenient to switch to number densities: ρ ( x_1,t) : = ρ_0 p( x_1,t) and Q ( x_1, x_2, t) := ρ_0^2 P_2( x_1, x_2, t). In what follows we outline his derivation for hard spheres ignoring hydrodynamic interactions. The equation for the one-particle number density ρ( x_1,t) is

    ∂ρ/∂ t  = ∇_  x_1·( ∇_  x_1 ρ  + ∫ Q( x_1,  x_2, t) ∇_  x_1 u(r)   x_2  ),

where r =  x_1-  x_2 is the interparticle distance and Q( x_1,  x_2, t) (the two-particle number density) satisfies, to lowest order in ρ_0, 

    ∂ Q/∂ t = ∇_  x_1·[ ∇_  x_1   Q + Q ∇_  x_1 u(r) ] + ∇_  x_2·[ ∇_  x_2   Q + Q ∇_  x_2 u(r) ].

Equations <ref> and <ref> have the following time-independent equilibrium solutions

    ρ_s( x_1) = ρ_0,      Q_s ( x_1, x_2) = ρ_0^2    g_0 (  x_1-  x_2 ),

where ρ_0 is constant  and g_0 is the pair correlation function

    g_0(r) = e^-u(r).

Felderhof then looks for a linearized solution around the equilibrium values <ref>, by making the ansatz

    Q( x_1, x_2, t) = ρ( x_1,t) ρ ( x_2,t) g( x_1, x_2, t),

and considering the deviations ρ_1 and g_1,

    ρ( x_1,t)  = ρ_0 + ρ_1( x_1,t),      g( x_1, x_2, t) = g_0(r
    )  + g_1( x_1, x_2, t).

Then, to terms linear in ρ_1 and g_1,  <ref> becomes

    Q( x_1, x_2, t)  ≈ρ_0^2   g_0(r) +  ρ_1( x_1,t) ρ_0   g_0(r) + ρ_0 ρ_1 ( x_2,t) g_0(r) + ρ_0^2   g_1( x_1, x_2, t).

Substituting in <ref> and linearizing gives 

    ∂ρ_1/∂ t = ∇_  x_1·[   ∇_  x_1ρ_1  + ρ_0 ∫ g_0(r) ρ_1( x_2, t)   ∇_  x_1 u(r)   x_2 
    
    + ρ_0^2 ∫ g_1( x_1, x_2, t) ∇_  x_1 u(r)   x_2 ].

At this stage, Felderhof supposes that perturbations from equilibrium <ref> are small so that ∇_  x_1ρ_1( x_1,t) ≈∇_  x_2ρ_1( x_2,t) and the pair correlation function is at its equilibrium value, that is, g( x_1, x_2, t) ≈ g_0(r). Then <ref> simplifies to

    ∂ρ_1/∂ t  = ∇_  x_1· [ ∇_  x_1ρ_1  + ρ_0 ∫ g_0(r) ρ_1( x_2, t)   ∇_  x_1 u(r)   x_2  ].

Now, using <ref> and expanding ρ_1( x_2,t) about x_1 and keeping only the first non-vanishing term gives <cit.>

    ∂ρ_1/∂ t  = ∇_  x_1·[ (1+  α_u ϵ^d ρ_0) ∇_  x_1ρ_1   ],

where α_u = ∫_ℝ^d ( 1 - e^- u(ϵ x))   x. Note that this is the evolution equation for the perturbation ρ_1 from the uniform equilibrium ρ_0 (valid with f= 0). 





 §.§ 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) <cit.>, and consists of truncating the hierarchy at the two-particle density function. To this end, we consider the equation satisfied by P_2(  x_1,   x_2, t) by integrating the N-particle Fokker–Planck equation <ref> over x_3, …,  x_N, applying the divergence theorem and relabelling particles as before to obtain

    ∂ P_2/∂ t
    =   ∇_  x_1·[ ∇_ x_1   P_2 -    f(  x_1)  P_2 + ∇__1 u( _1 - _2  ) P_2 +   G_2(_1, _2, t) 
      ]  
    
     +   ∇_  x_2·[ ∇_ x_2   P_2 -    f(  x_2)  P_2 + ∇__2 u( _1 - _2  ) P_2 +   G_2(_2, _1, t)  ],

where

    G_2(_1, _2, t)  =  (N-2) ∫_Ω∇__1 u( _1 - _3  ) P_3( _1,_2,_3,t)   _3,

and 

    P_3(  x_1,  x_2, _3, t) = ∫_Ω^N-3 P(x⃗, t)   x_4 ⋯ x_N

is the three-particle density function. We note that in writing <ref> we are using that P_3 is invariant to particle relabelling. 
The KSA then approximates the three-particle density function as

    P_3( x_1, x_2, x_3,t) =  P_2 ( x_1, x_2,t) P_2 ( x_1, x_3,t) P_2 ( x_2, x_3,t) /p ( x_1,t) p ( x_2,t) p ( x_3,t) ,

where p and P_2 are the one- and two-particle density functions, respectively. Inserting <ref> into <ref> one can then solve the coupled system <ref> and <ref> 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 <cit.>. 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 <cit.> and off-lattice cell motility processes with soft interactions <cit.>.  Middleton, Fleck and Grima <cit.> consider a system of Brownian particles evolving according to <ref> in one dimension interacting via a Morse potential (see <ref>) and compare the KSA closure to the MFA closure and simulations of the stochastic system. Markham et al. <cit.> 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 <cit.> use a variant of the KSA closure, approximating either P_2(  x_1,  x_3, t) or P_2(  x_2,  x_3, t)  in <ref> by their corresponding mean-field approximation P_2(  x_i,  x_j, t) = p(  x_i, t) p(  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 G and G_2. As we shall see in <ref>, even in one physical dimensional the KSA model is rather complicated to solve. 
 
 
 


§ 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 <ref> for the one-particle density p. Assuming that the pair potential u is localized near x_1, we can determine G in <ref> 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(  x_1,  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 G_2 in <ref> are negligible and the equation for P_2 reduces to 

    ∂ P_2/∂ t
    =   ∇_  x_1·[ ∇_ x_1   P_2 -    f(  x_1)  P_2 + ∇__1 u( _1 - _2  ) P_2 
      ]  
    
     +   ∇_  x_2·[ ∇_ x_2   P_2 -    f(  x_2)  P_2 + ∇__2 u( _1 - _2  ) P_2   ],

for  (  x_1,   x_2) ∈Ω ^2, complemented with periodic boundary conditions on ∂Ω^2. Note that in approximating P_2 in this way it is no longer true that p(_1,t) = ∫ P_2(_1, _2,t) _2, so that p and P_2 need to be solved for as a coupled system. This equation is basically <ref> with an added external force. Essentially our MAE approach aims to solve <ref> systematically asymptotically rather than through Felderhof's linearization and approximations (see <ref>). 





 §.§ Inner and outer regions
 


By assumption, the pair interaction potential u(r) is negligible everywhere except when the interparticle distance r is of order ϵ. Therefore, we suppose that when two particles are far apart (x_1 -
  x_2≫ϵ) they are independent, whereas when
they are close to each other (x_1 -   x_2∼ϵ) 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_o(  x_1,   x_2, t) = P_2(   x_1,   x_2, t). By independence, we have that[Independence only tells us that P_o(  x_1,   x_2, t) ∼ q(  x_1, t) q(  x_2, t) for some function q but the normalization condition on P implies p = q + O(ϵ).]

    P_o(  x_1,   x_2, t) = p(  x_1, t) p(  x_2, t) + ϵ P_o^(1)(  x_1,   x_2, t) + ⋯,

for some function P_o^(1). 

In the inner region, we set  x_1 = x̃_1 and x_2 =
x̃_1 + ϵx̃ and define P̃ (x̃_1, x̃, t) = P_2(  x_1,   x_2 , t) and ũ (x̃ ) = u( x_1-   x_2). Rewriting <ref> in terms of the inner coordinates gives 

    ϵ^2 ∂P̃/∂ t
    =   2 ∇_x̃·[ ∇_x̃P̃ + ∇_x̃ũ( x̃ ) P̃] + ϵ∇ _x̃·{[   f(x̃_1) -     f(x̃_1 + ϵx̃) ] P̃}
     -  ϵ∇_x̃_1·[ 2 ∇_x̃P̃ + ∇_x̃ũ (  x̃ ) P̃] + ϵ^2 ∇_x̃_1^2  P̃ - ϵ ^2  ∇ _x̃_1·[   f(x̃_1)  P̃ ].

The inner solution P̃ must match with the outer solution P_o as x̃→∞.  Expanding P_o in terms of the inner variables gives (omitting the time variable for ease of notation)

    P_o (   x_1,  x_2)    ∼ p(x̃_1) p(x̃_1 + ϵx̃) + ϵ P_o^(1)(x̃_1, x̃_1 + ϵx̃) 
       ∼  p^2 (x̃_1) +  ϵ p(x̃_1)  x̃·∇ _x̃_1 p( x̃_1) + ϵ P_o^(1)(x̃_1, x̃_1)      as    x̃→∞.


We look for a solution of <ref> matching with <ref>  as x̃→∞ of the form
P̃∼P̃^(0) + ϵP̃^(1) +  ⋯.
The leading-order inner problem is 

    0 = 2 ∇_x̃·[ ∇_x̃P̃^(0) + ∇_x̃ũ( x̃ ) P̃^(0)],      P̃^(0)∼ p^2 (x̃_1)   as  x̃→∞,

with solution 

    P̃^(0) =  p^2 (x̃_1)   e^-ũ( x̃ ).
 

The O(ϵ) problem reads



    0   =2 ∇_x̃·[∇_x̃P̃^(1) + ∇_x̃ũ( x̃) P̃^(1)]  -    ∇_x̃_1·[ 2 ∇_x̃P̃^(0) + ∇_x̃ũ( x̃) P̃^(0)], 
    
    P̃^(1)   ∼  p(x̃_1)  x̃·∇ _x̃_1 p( x̃_1) + P_o^(1)(x̃_1, x̃_1)      as    x̃∼∞.


Using <ref>, we can rearrange <ref> to give

    ∇_x̃·[ ∇_x̃P̃^(1) + ∇_x̃ũ( x̃) P̃^(1)  - 12∇_x̃_1P̃^(0)]  = 0.

Solving <ref> together with <ref> gives

    P̃^(1) =  [ p (x̃_1)   x̃·∇ _x̃_1 p( x̃_1) + P_o^(1)(x̃_1, x̃_1) ] e^-ũ( x̃).
 
Thus we find that the inner region solution is, to O(ϵ),

    P̃∼[  p^2 (x̃_1,t) + ϵ  p (x̃_1,t)  x̃·∇ _x̃_1 p( x̃_1,t) + ϵ P_o^(1)(x̃_1, x̃_1,t) + ⋯] e^-ũ( x̃).






 §.§ Interaction integral
 

Now we go back to the interaction integral G( x_1) in <ref>. 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 <ref> to evaluate it. 

First we split the integration volume Ω for 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 δ, with ϵ≪δ≪ 1, which divides the two regions. Then the inner region is Ω_i( x_1) = { x_2 ∈Ω :   x_2-  x_1  <δ} and the outer region is the complimentary set Ω_o( x_1) = Ω∖Ω_i( x_1). The dominant contribution to <ref> is then

    G ( x_1, t)    = (N-1) ∫_Ω_i( x_1)   P( x_1,  x_2, t) ∇_ x_1 u(  x_1-  x_2 )    x_2 
       = -(N-1) ϵ^d-1∫_x̃ < δ/ϵP̃(  x_1, x̃, t)  ∇_x̃ũ( x̃)   x̃
       ∼  (N-1) ϵ^d-1  p (  x_1,t) ∫_x̃ < δ/ϵ{ p (  x_1,t) + ϵx̃·∇ _ x_1 p(   x_1,t) 
        +  ϵ P_o^(1)(   x_1,   x_1)   }∇_x̃   e^-ũ( x̃) x̃.

The first and third terms of the integral vanish using the divergence theorem and that ũ is a radial function of x̃. Integration by parts on the second component gives

    ∫_x̃ < δ/ϵ[ x̃·∇ _ x_1 p(  x_1,t) ]  
    ∇_x̃   e^-ũ( x̃) x̃
    ∼∇ _ x_1 p(  x_1,t)[ 2(d-1)πd ( δϵ)^d - ∫_x̃ <  δ/ϵ   e^-ũ( x̃) x̃ ],

where we have used that e^-ũ( x̃)≈ 1 at x̃ = δ/ϵ. Finally, rewriting the first term above as a volume integral,[The volume of a d-dimensional sphere of radius δ/ϵ  is equal to .]   <ref> becomes

    G( x_1, t)  ∼  (N-1) ϵ^d p(  x_1,t) ∇ _ x_1 p( x_1,t) ∫_x̃ <  δ/ϵ ( 1 - e^-ũ( x̃))  x̃.

Since 1 - e^-ũ( x̃) decays at infinity, we can extend the domain of integration to the entire ℝ^d introducing only lower order errors. Therefore we can write

    G(  x_1, t)  ∼α_u (N-1) ϵ^d p(  x_1,t) ∇ _ x_1 p(   x_1,t),

with

    α_u =  ∫_ℝ^d ( 1 - e^-ũ( x̃))  x̃ = ∫_ℝ^d ( 1 - e^- u( ϵ x))   x.

Note that we obtain the same coefficient α_u that appeared in equation <ref> using the closure at the pair correlation function. 




 §.§ Reduced Fokker–Planck equation for soft spheres
 


Combining <ref> with  <ref> we find that, to O(ϵ^d),



    ∂ p/∂ t=  ∇_  x_1·{ [1 +  α_u (N-1) ϵ^d p ] ∇_  x_1  p -    f(  x_1) p },

where 

    α_u = ∫_ℝ^d ( 1 - e^- u(ϵ x))   x.


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

    α_u = ∫_ℝ^d ( 1 - e^- u( ϵ x))   x = ∫_ℝ^d  u(ϵ x)  x - 1/2∫_ℝ^d u^2(ϵ x)   x + ⋯,

that is, the LMFA closure is the leading contribution of the potential, provided that u(ϵ x) is small. However, as we will see in <ref>, this is not always true. We also see that the equation obtained by Felderhof <cit.>, <ref> is the linearized version of <ref> after taking  N →∞ and setting the external force to zero. 

The coefficient α_u can be related to basic concepts from statistical mechanics. Namely, the integrand in <ref>  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 <cit.>. 




 §.§ 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(_1, _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 <ref> plus the outer expansion <ref> minus the common part <cit.>. The result is

    P_2( x_1,  x_2, t) ∼ p( x_1,t) p( x_2,t) e^-u(  x_1 -  x_2  ).





§ RESULTS
 


 




 §.§ 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 <ref> obtained via MAE can be used to model a system of hard-core interacting particles of diameter ϵ. 
In particular, the model <ref> for soft spheres has exactly the same structure as the counterpart model for hard spheres derived in <cit.>. Inserting the hard-sphere pair potential

    u_HS(r ) = {   ∞     r ≤ϵ, 
       0      r >ϵ,
     .

into the nonlinear diffusion coefficient in <ref> gives

    α_HS =  ∫_ x<1 x,

that is, α_HS = 2 for d=1, α_HS = π for d=2, and α = 4π /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 <cit.>. 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 α_HS in <ref> 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 ϵ_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 ϵ_eff – are equivalent. In other words, a characterization of a given soft potential u is to find ϵ_eff such that 

    α_u ϵ^d = α_HS ϵ_eff^d.

where α_u is given in <ref> and α_HS given in <ref>. Rearranging, we find that

    (ϵ_eff/ϵ)^d = d  ∫_0^∞ ( 1 - e^- u(ϵ r)) r^d-1  r.


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 <cit.>. 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 <ref> 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 <ref>. This is provided α_u in <ref> is well-defined and positive. In contrast, if α_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 <ref> 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 ϵ_eff, one can define the effective volume fraction of soft spheres and use it to check whether the “low-volume fraction” condition holds. 
 



 §.§ Comparison between MAE and LMFA
 






The simplest repulsive pair potential is the soft-sphere (SS) potential, which assumes the form

    u_SS (r) =  (ϵ/r)^ν,

where ϵ is a measure of the range of the interaction and ν is  the hardness parameter which characterizes the particles (the softness parameter is defined
as its inverse, 1/ν). We note that for ν =1, the SS potential corresponds to the Coulomb interaction (ν =1)  <cit.>. Other common purely repulsive potentials include the exponential (EX) potential <cit.>

    u_EX (r) = e^-r/ϵ,

and the repulsive Yukawa (YU) potential 

    u_YU (r) = ϵ/r e^-r/ϵ.

This potential, also known as screened Coulomb, is used to describe elementary particles, small charged “dust” grains observed in plasma environments, and suspensions of charge-stabilized colloids <cit.>. 

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 <ref> may be generalized to power-law repulsive–attractive potentials of the form u(r) = (ϵ/r)^a - C(ϵ/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,

    u_LJ (r) =  ( ϵ/r)^12 -  ( ϵ/r)^6.

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

    u_MO(r) = e^-r/ϵ - 1/C e^-lr/ϵ,

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 <cit.>. 

<ref> shows examples of all the interaction potentials above. 
0.75
1.0


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 α_u in <ref>  and α_u in <ref> for the potentials above are shown in <ref>. Since the behaviour of the  integrals in α_u and α_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 <ref>. 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 α_u and α_u for the cases for which the latter is defined.  



The coefficient  α_u is undefined for inverse-power potentials such as SS and LJ, since the integral in <ref> 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 α_u exists for the SS potential in <ref> for ν>d, but is undefined for ν≤ d. Therefore, we find that MAE is not valid for the Coulomb interaction [ν = 1 in <ref>]. 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 G ( x_1) is negligible as assumed in the MAE derivation). 

<ref> shows the variation of α_SS against the hardness parameter ν. As expected, lim_ν→∞α_SS = α_HS since the HS potential is the limiting case of the SS potential for ν→∞. Also note that the strength of the nonlinear diffusion term, parametrized by α_SS, decreases as the hardness ν increases. The effect of the softness parameter 1/ν on the diffusion and other transport coefficients was studied in <cit.>; they found that for ν≥ 72 the soft sphere system behaved essentially as a hard sphere system in molecular-dynamics simulations. We find that the relative error between α_SS and α_HS  for ν = 72 in three dimensions is 2.6%. 
0.7
1.0


The Morse potential <ref> has well-defined coefficients α_MO and α_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 <ref>.  We note that both coefficients  α_MO and α_MO may be negative for some parameter values, which implies, from <ref>, 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→∞ <cit.>.   
0.7
1.0


Evaluation of α_u and α_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 <ref>, α_u (N-1) ϵ^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 <ref> and ascertain whether the low-volume fraction assumption holds). Regarding the MFA closure, one should keep in mind that α_u not being defined does not necessarily mean that standard MFA fails, as mentioned in <ref>; it could be that the short-range assumption used to obtain the nonlinear diffusion equation <ref> does not hold and that, instead, one should keep the original MFA integro-differential model <ref>. 
Conversely, if α_u is not defined due to the singular behavior of the potential at the origin (which is the case for all those seen in <ref>), this is an indication that the convolution in the original MFA will not exist. 





 §.§ 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 <cit.> to perform the particle-level simulations. The overdamped Langevin equation <ref> is integrated using the Euler-Maruyama method and a constant timestep Δ t, leading to an explicit update step for each particle given by

    𝐗^m+1_i = 𝐗^m_i + √(2DΔ t)Δ W^m +  f  (  X_i^m ) Δ t  - ∑_j  i∇_ x_i u( X_i^m -  X_j^m) Δ t,

where Δ W^m is a d-dimensional normally distributed random variable with zero mean and unit variance. We choose the timestep Δ 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 Δ t = (0.1ϵ)^2/2D was sufficient for convergence, leading to an average diffusion step size of 0.1ϵ.

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ϵ. All particle pairs i,j separated by a distance greater than this cutoff will implicitly be given a pair potential of u( X_i -  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 Ω = [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 <ref> 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 <ref> and <ref>, 

    ∂ p/∂ t  = ∂/∂ x_1[ ∂ p /∂ x_1  + (N-1 )   p∫_Ω p( x_2, t) f_u(x_1, x_2)     x_2  ],

where p = p(x_1,t) unless explicitly written and  f_u(x_1, x_2) := ∂/∂ x_1  u(|x_1 - x_2|). To solve <ref> 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 . To evaluate the integral in <ref>, we use the periodic trapezoid rule (which converges exponentially fast for smooth integrands <cit.>). 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) ≈ (p(x_2-h/2,t) + p(x_2+h/2,t))/2. Because of the convolution term, the discretization matrix of <ref> 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 <cit.>, 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
 

    ∂ p/∂ t = ∂/∂ x_1[ ∂ p /∂ x_1 + (N-1 ) ∫_Ω P_2  f_u(x_1, x_2)    x_2  ],


    ∂ P_2/∂ t
    =       ∂/∂ x_1[ ∂ P_2/∂ x_1   +  f_u(x_1, x_2) P_2  
        +  (N-2)  P_2 (x_1,x_2,t) p (x_1,t) p (x_2,t)∫_Ω P_2 (x_1,x_3,t) P_2 (x_2,x_3,t) p (x_3,t)   f_u(x_1, x_3)    x_3
      ]   
        +  ∂/∂ x_2[ ∂ P_2/∂ x_2   +  f_u(x_2, x_1) P_2  
         +  (N-2)  P_2 (x_1,x_2,t) p (x_1,t) p (x_2,t)∫_Ω P_2 (x_1,x_3,t) P_2 (x_2,x_3,t) p (x_3,t)   f_u(x_2, x_3)    x_3
      ],


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 Ω^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 <ref>, 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) = ∫_Ω P_2(x_1,x_2,t)   x_2 (replacing <ref>) <cit.>. While this is true in the infinite BBGKY hierarchy, we note that once the KSA closure <ref> is adopted, the resulting model is not necessarily equal to the KSA system <ref>. 



Finally, the reduced model from the MAE method reads

    ∂ p/∂ t = ∂/∂ x_1{ [1+ α_u (N-1) ϵ^d p ] ∂ p /∂ x_1} ,

where α_u = 2 ∫_0^∞  (1-e^-u (ϵ r) )   r. Noting that the right-hand side in <ref> can be written as [p + 1/2α_u (N-1) ϵ^d p^2]_x x, the numerical implementation is straightforward and, since the discretization matrix is banded, very efficient.

In the first example, we use the exponential potential u_EX <ref> with ϵ = 0.05 in a system in N=16 particles. The particles are initially distributed according to p(0) = 0.5[tanh(β(x-θ))+tanh(β(1-θ-x))], with β = 30 and θ = 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 <ref>. 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 ϵ and N, the potential is not short-ranged and the assumptions of the MAE model break down. In particular, using <ref> we find that the effective hard-sphere radius is ϵ_eff = 0.04 and hence the effective volume fraction is 0.64. Similarly, the solution obtained via the LMFA (shown at T_f in <ref>(b)) differs noticeably to the MFA solution. 
0.7
1.0



<Ref> 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 <ref> and <ref> 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 <ref>(a), which shows a plot of P_2(x_1,x_2,T_f) along x_2 = 0.5. 
0.7
1.0


Next we consider an example using a more repulsive interaction potential. In particular, we consider a smoothed version of the Yukawa potential <ref>, namely u (r) =(ϵ /√(r^2 + δ^2)) exp^-r/ϵ with ϵ = 0.01 and δ = 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 <ref>). We run a simulation with N=20 particles up to time T_f = 0.02. For these parameters, the effective hard-sphere radius is ϵ_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 <ref>, and the results for the two-particle density function are shown in <ref>. 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 (<ref>). 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 (<ref>(b)).
0.7
1.0


0.7
1.0




Generally, we see that MAE is good for strongly repulsive interactions (such as in the example in <ref>) while MFA is better for softer or longer-range interactions (<ref>). 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) ≈ 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 =∬exp(-u(|x_1-x_2|))  x_1  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 <ref>.  
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 (<ref>(a)), whereas in the second example (with strong repulsion) the three approximations MFA, MFAe, and LFMAe are poor (<ref>(b)). 

Finally, we consider a two-dimensional example in Ω = [0,1] ^2. We choose a system of N =400 Yukawa-interacting particles, with interaction potential u_YU  <ref> and ϵ = 0.01, initially distributed according to a normal distribution in x with mean 0.5 and standard deviation σ = 0.05. The effective hard-sphere radius is ϵ_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_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 (ϵ=0) for reference. We plot the comparison at times T_f/2 and T_f in  <ref>. 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 (<ref>) and the localized version LMFA <ref>, which is computationally much easier to solve (the two curves would lie on top of each other in <ref>, 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.

0.7
1.0





§ 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.[These models are
   sometimes used regardless, with a numerical discretisation
   providing an ad hoc  regularisation of the
   convolution integral <cit.>.]
The MAE result is a nonlinear
 diffusion equation
similar to our previous work for hard-spheres <cit.>, with 
the  coefficient of the nonlinear term α_u depending on the
potential  through (<ref>).

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 <cit.> 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 <ref>) 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 α_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 <ref> with the nonlocal
convolution term due to the long-range component of the potential, and
a nonlinear diffusion due to the short-range component.  

10

Baker:2010is
R. E. Baker and M. J. Simpson, Correcting mean-field approximations
  for birth-death-movement processes, Phys. Rev. E, 82 (2010), p. 041905.

Berlyand:2016fu
L. Berlyand, P.-E. Jabin, and M. Potomkin, Complexity Reduction in
  Many Particle Systems with Random Initial Data, SIAM/ASA J. Uncertainty
  Quantification, 4 (2016), pp. 446–474.

Binny:2015ei
R. N. Binny, M. J. Plank, and A. James, Spatial moment dynamics for
  collective cell movement incorporating a neighbour-dependent directional
  bias, J. R. Soc. Interface, 12 (2015), pp. 20150228–20150228.

Bruna:2012cg
M. Bruna and S. J. Chapman, Excluded-volume effects in the
  diffusion of hard spheres, Phys. Rev. E, 85 (2012), p. 011103.

Bruna:2013ku
M. Bruna and S. J. Chapman, Diffusion of finite-size particles in
  confined geometries, Bull. Math. Biol., 76 (2014), pp. 947–982.

Burger:2010gb
M. Burger, M. Di Francesco, J.-F. Pietschmann, and B. Schlake, Nonlinear Cross-Diffusion with Size Exclusion, SIAM J. Math. Anal., 42
  (2010), p. 2842.

Carrillo:2010tx
J. A. Carrillo, M. Fornasier, G. Toscani, and F. Vecil, Particle,
  kinetic, and hydrodynamic models of swarming, in Mathematical Modeling of
  Collective Behavior in Socio-Economic and Life Sciences, G. Naldi,
  L. Pareschi, and G. Toscani, eds., Birkhäuser Basel, 2010, pp. 297–336.

Carrillo:2009dz
J. A. Carrillo and J. S. Moll, Numerical Simulation of Diffusive
  and Aggregation Phenomena in Nonlinear Continuity Equations by Evolving
  Diffeomorphisms, SIAM J. Sci. Comput., 31 (2009), pp. 4305–4329.

Dolbeault:2001jc
J. Dolbeault, P. A. Markowich, and A. Unterreiter, On Singular
  Limits of Mean-Field Equations, Arch. Ration. Mech. An., 158 (2001),
  pp. 319–351.

DOrsogna:2006ci
M. R. D’Orsogna, Y. Chuang, A. L. Bertozzi, and L. S.
  Chayes, Self-Propelled Particles with Soft-Core Interactions:
  Patterns, Stability, and Collapse, Phys. Rev. Lett., 96 (2006), p. 104302.

Felderhof:1978vn
B. U. Felderhof, Diffusion of interacting Brownian particles, J.
  Phys. A: Math. Gen., 11 (1978), p. 929.

Hansen:2006uv
J. P. Hansen and I. R. McDonald, Theory of Simple Liquids,
  Academic Press, London, 2006.

Heyes:2005ff
D. M. Heyes and A. C. Brańka, The influence of potential
  softness on the transport coefficients of simple fluids, J. Chem. Phys.,
  122 (2005), p. 234504.

Hinch:1991go
E. J. Hinch, Perturbation Methods, Cambridge University Press,
  1991.

Horng:2012io
T.-L. Horng, T.-C. Lin, C. Liu, and B. Eisenberg, PNP Equations
  with Steric Effects: A Model of Ion Flow through Channels, J. Phys. Chem.
  B, 116 (2012), pp. 11422–11441.

Hynninen:2003hr
A.-P. Hynninen and M. Dijkstra, Phase diagrams of hard-core
  repulsive Yukawa particles, Phys. Rev. E, 68 (2003), p. 021407.

Israelachvili:2011ug
J. N. Israelachvili, Intermolecular and Surface Forces, Revised
  Third Edition, Academic Press, 2011.

Kirkwood:1935is
J. G. Kirkwood, Statistical Mechanics of Fluid Mixtures, J. Chem.
  Phys., 3 (1935), pp. 300–313.

liggett1999stochastic
T. M. Liggett, Stochastic interacting systems: Contact, voter, and
  exclusion processes, vol. 324, Springer-Verlag, Berlin, 1999.

Markham:2013tv
D. C. Markham, M. J. Simpson, P. K. Maini, E. A. Gaffney, and R. E.
  Baker, Incorporating spatial correlations into multispecies mean-field
  models, Phys. Rev. E, 88 (2013), p. 052713.

Middleton:2014fa
A. M. Middleton, C. Fleck, and R. Grima, A continuum approximation
  to an off-lattice individual-cell based model of cell migration and
  adhesion, J. Theor. Biol., 359 (2014), pp. 220–232.

Mogilner:1999iy
A. Mogilner and L. Edelstein-Keshet, A non-local model for a
  swarm, J. Math. Biol., 38 (1999), pp. 534–570.

Mulero:2012vc
A. Mulero, Theory and simulation of hard-sphere fluids and related
  systems, vol. 753 of Lecture Notes in Physics, Springer-Verlag, Berlin,
  2008.

aboria
M. Robinson, Aboria library, 2016,
  <https://martinjrobins.github.io/Aboria/> (accessed 12/20/2015).
Version 0.1.

Robinson:2017vxa
M. Robinson, and M. Bruna, Particle-based and meshless Methods with Aboria, SoftwareX, In Press (2017), <http://dx.doi.org/10.1016/j.softx.2017.07.002>.

Singer:2004ec
A. Singer, Maximum entropy formulation of the Kirkwood
  superposition approximation, J. Chem. Phys., 121 (2004), p. 3657.

Trefethen:2000jp
L. N. Trefethen, Spectral Methods in MATLAB, Society for
  Industrial and Applied Mathematics, 2000.




