# Heterogeneity-induced lane and band formation in self-driven particle
systems

Basma Khelfa  School for Mechanical Engineering and Safety Engineering,
University of Wuppertal, Wuppertal, Germany  Raphael Korbmacher  School for
Mechanical Engineering and Safety Engineering, University of Wuppertal,
Wuppertal, Germany  Andreas Schadschneider  Institute for Theoretical Physics,
University of Cologne, Cologne, Germany  Antoine Tordeux  School for
Mechanical Engineering and Safety Engineering, University of Wuppertal,
Wuppertal, Germany  e-mail<EMAIL_ADDRESS>
###### Abstract

The collective motion of interacting self-driven particles describes many
types of coordinated dynamics and self-organisation. Prominent examples are
alignment or lane formation which can be observed alongside other ordered
structures and nonuniform patterns. In this article, we investigate the
effects of different types of heterogeneity in a two-species self-driven
particle system. We show that heterogeneity can generically initiate
segregation in the motion and identify two heterogeneity mechanisms.
Longitudinal lanes parallel to the direction of motion emerge when the
heterogeneity statically lies in the agent characteristics (quenched
disorder). While transverse bands orthogonal to the motion direction arise
from dynamic heterogeneity in the interactions (annealed disorder). In both
cases, non-linear transitions occur as the heterogeneity increases, from
disorder to ordered states with lane or band patterns. These generic features
are observed for a first and a second order motion model and different
characteristic parameters related to particle speed and size. Simulation
results show that the collective dynamics occur in relatively short time
intervals, persist stationary, and are partly robust against random
perturbations.

##  1 Introduction

The emergence of coordinated movement and self-organisation in the collective
motion of large groups of interacting autonomous individuals is ubiquitous in
nature[1, 2, 3, 4]. Flock and swarm behaviour arises e.g. in bacterial and
cell migration[5, 6] (see also the pioneering works on morphogenesis by A.M.
Turing[7]), animal aggregation [8, 9, 10, 11], rods and non-living systems
[12, 13, 14], microswimmers [15, 16, 17], but also in robotic systems [18,
19]. This includes a large class of non-equilibrium systems of self-driven
particles often called _active matter_ in the literature. The microscopic
interaction rules initiate _motility-induced phase separation_ from disordered
states to non-uniform macroscopic patterns displaying order (at least
partially)[20, 21]. Pedestrian and vehicle flows are also systems of
interacting self-driven agents that can exhibit phase transitions and self-
organisation[22, 23, 24, 25, 26]. Indeed, pedestrians and drivers interact
locally with their environment and the neighbourhood and have by the way
density-dependent motility parameters. Similar forms of collective dynamics
can arise in more abstract conceptions of self-driven particles like social
systems, social networks, and opinion formation[27, 28, 29, 30]. Regardless of
the important differences in the specific type and interaction of the
individuals, collective motions exhibit similarities suggesting the existence
of generic underlying mechanisms and systemic phase transitions[20, 31, 32].

The terminology active, self-driven, or self-propelled particles dates back to
the 1970s when it was introduced for the description of swimming micro-
organisms[33]. It became popular during the 1990s with the Vicsek model[8] and
is nowadays used for animal aggregation [10], bacterial migration[6],
suspensions of microswimmers[16], or vehicle and pedestrian dynamics[34].
Self-driven particle models are mainly designed to describe collective
dynamics. In contrast to Brownian colloids dominated by physical (exclusion)
rules, the motion of self-driven particles is mainly governed by local
interaction rules with the environment. The interactions may result from
external chemical potentials, for instance in cell behaviour (chemotaxis),
while they come from social and proxemics interaction rules for pedestrians
(sociotaxis)[35]. The social rules of pedestrians can initiate collective
phenomena improving the system performance. The collective motions have
positive effects on transport properties and may be referred to as
_intelligent collective dynamics_. Examples of intelligent collective dynamics
of pedestrians are lane formation[36, 37, 38, 32], collective motion and
alignment[39, 21], or intermittent flows at bottlenecks and diagonal, chevron
or circular patterns at intersections[40, 41]. Intelligent collective dynamics
in transportation are also of practical relevance, e.g. for the control of
autonomous driving systems[42, 43]. Yet, self-organisations may also induce
negative effects on traffic safety and performance, for instance through stop-
and-go waves[44, 45, 46], herding[35] or clogging effects[47, 48] among other
segregation phenomena[49]. Besides scientific interests, extracting the
essential features of collective motions from individual behaviors is
fundamental to authorities for the control of crowd and traffic dynamics and
the development of intelligent transportation strategies.

In self-driven particle systems, collective dynamics can result from
heterogeneity effects in the microscopic behaviour of the particles, upon
other inertia or delay mechanisms. Pedestrian dynamics describe for instance
lane formation for counter-flow or for pedestrians walking in the same
direction but with different speeds[32, 50]. Other examples are stripe,
diagonal travelling band or chevron patterns for crossing flows[40, 41]. In
this article, we show by simulation that heterogeneity effects can generically
initiate segregation and spontaneous formation of lane or band patterns in
two-species flows of polarised agents. Two heterogeneity mechanisms are
identified: static heterogeneity in the agent characteristics and dynamic
heterogeneity in the interactions. Static heterogeneity refers to _quenched
disorder_ in solid state physics and the terminology of random walks, when
dynamic heterogeneity relies on _annealed disorder_ (see[51, 52] and
references therein). Interestingly, lanes spontaneously occur when the
heterogeneity relies statically on the agent features (quenched disorder),
while bands emerge if the heterogeneity operates dynamically in the
interactions (annealed disorder). The features are generically observed with
different microscopic motion models, namely the first order collision-free
speed model[53] and the inertial second order social force model[54], and
different types of parameters related to agent speed or agent size.

###  1.1 Models

We consider in the following two types of agents evolving on a torus. We
denote $n=1,\ldots,N$ the agent’s ID while $k_{n}={1,2}$ is the agent’s type.
The agent’s motion is given by a dynamic model
$F_{\mathbf{p}}(\mathbf{X}_{n})$ that defines the agent speed as in the
collision-free model[53] or the agent acceleration as in the social force
model[54] according to local spatio-temporal variables $\mathbf{X}_{n}$ (e.g.
the position and speed differences with the neighbours) and a set of
parameters $\mathbf{p}$ (namely, desired speed, desired time gap, repulsion
rate, agent size, and so on). We assume two different settings
$\mathbf{p}_{1}$ and $\mathbf{p}_{2}$ for the parameters. Two types of
heterogeneity are then considered.

  1. 1.

_Heterogeneity in the agent characteristics_ – We attribute statically the two
parameter settings $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$ to the two types of
agents:

$M_{1}(n,k_{n})=F_{\mathbf{p}_{k_{n}}}(\mathbf{X}_{n}).$ (1)

We aim here to model different types of agents (for instance pedestrians and
bicycles) with specific characteristics in term of desired speed, agent size,
etc. This kind of heterogeneity is usually called _quenched disorder_ in solid
state physics. It refers to static heterogeneity features remaining constant
(i.e. quenched) over the time.

  2. 2.

_Heterogeneity in the interactions_ – We attribute dynamically the two
parameter settings $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$ according to the type
of the closest agent in front. The parameter setting is $\mathbf{p}_{1}$ if
the agent in front is of the same type, while it is $\mathbf{p}_{2}$ in case
of interaction with an another agent type:

$M_{2}(n,k_{n})=\left\\{\begin{array}[]{ll}F_{\mathbf{p}_{1}}(\mathbf{X}_{n}),&\text{if
$\tilde{k}(\mathbf{X}_{n})=k_{n}$,}\\\\[2.84526pt]
F_{\mathbf{p}_{2}}(\mathbf{X}_{n}),&\text{otherwise,}\end{array}\right.$ (2)

with $\tilde{k}(\mathbf{X}_{n})$ the type of the closest agent in front (see
Sec. 4 for details). Such a mechanism may correspond for instance to a
pedestrian adapting his/her behaviour when interacting with a bicycle and
inversely. The heterogeneity features are here time-dependent. They are
usually called _annealed disorder_ in the literature of solid state
physics[51, 52].

In contrast to the model Eq. (1) for which the heterogeneity statically lies
in agent characteristics, the model Eq. (2) induces a dynamic heterogeneity
mechanism taking place in the interactions. See Fig. 1 for an illustrative
example in one dimension.

Model $M_{1}$Model $M_{2}$Agent type $1$Agent type $2$Parameter
$\mathbf{p}_{1}$Parameter $\mathbf{p}_{2}$ Figure 1: Illustrative scheme in
one dimension for the two heterogeneity models. The parameter setting (orange
or blue) depends on the type of agents (represented as disc and triangle) for
the model Eq. (1) while the setting depends on the type of the agent in front
for the model Eq. (2): it is orange when the agent in front is of the same
type and blue if it is of another type.

###  1.2 Analysis

We qualitatively observe by simulation that the static heterogeneity model
$M_{1}$ Eq. (1) initiates the formation of lanes in the system, while the
dynamic heterogeneity model $M_{2}$ Eq. (2) allows the formation of bands (see
Fig. 2 below). To classify the state of the system, we measure the agent’s
mean speed and also order parameters for the lane and band formation. The
order parameter has been introduced to detect lanes in a colloidal
suspension[55] and used in pedestrian dynamics[56]. We denote in the following
$(x_{n},y_{n})$ the positions of the agents $n=1,\ldots,N$. The order
parameter for lane formation is

$\Phi_{L}=\frac{1}{N}\sum_{n}\phi_{n}^{L}\qquad\text{with}\qquad\phi^{L}_{n}=\Big{[}\frac{L_{n}-\bar{L}_{n}}{L_{n}+\bar{L}_{n}}\Big{]}^{2}\qquad\left|\begin{array}[]{l}L_{n}=\text{card}\big{(}m,\,|y_{n}-y_{m}|<\Delta,\leavevmode\nobreak\
k_{n}=k_{m}\big{)}\\\\[2.84526pt]
\bar{L}_{n}=\text{card}\big{(}m,\,|y_{n}-y_{m}|<\Delta,\leavevmode\nobreak\
k_{n}\neq k_{m}\big{)}\end{array}\right..$ (3)

Here $L_{n}$ is the number of agents on a lane of width $\Delta>0$ in front of
the agent $n$ with the same type, $\text{card}(A)$ being the operator counting
the elements of an ensemble $A$, while $\bar{L}_{n}$ is the number of agents
with different types. The order parameter $\Phi_{L}$ tends by construction to
be close to one when the system describes lanes. Assuming a disordered state
for which the agents are uniformly randomly distributed on a $w\times h$
rectangle with $h>\Delta>0$ the system’s height and $w>0$ the system’s width,
the number $L_{n}$ of agents with the same type is distributed according to
the binomial model $\mathcal{B}(m,p)$, with $m=N_{k_{n}}$ and $p=\Delta/h$.
Here $N_{k_{n}}$ is the total number of agents with type $k_{n}$. The
distribution of the number $\bar{L}_{n}$ of agents with different types can be
deduced similarly.

For band formation, the order parameter is

$\Phi_{B}=\frac{1}{N}\sum_{n}\phi_{n}^{B}\qquad\text{with}\qquad\phi^{B}_{n}=\Big{[}\frac{B_{n}-\bar{B}_{n}}{B_{n}+\bar{B}_{n}}\Big{]}^{2}\qquad\left|\begin{array}[]{l}B_{n}=\text{card}\big{(}m,\,|x_{n}-x_{m}|<\Delta
w/h,\leavevmode\nobreak\ k_{n}=k_{m}\big{)}\\\\[2.84526pt]
\bar{B}_{n}=\text{card}\big{(}m,\,|x_{n}-x_{m}|<\Delta
w/h,\leavevmode\nobreak\ k_{n}\not=k_{m}\big{)}\end{array}\right..$ (4)

The band order parameter includes a term $w/h$, $w$ and $h$ being the width
and height of the system. The distribution of the order parameters for lanes
and bands is by construction the same in cases of random positions of the
agents. Indeed for disordered states, the number $B_{n}$ of agents on the
sides with the same type has a binomial distribution $\mathcal{B}(m,p)$ with
$m=N_{k_{n}}$ and $p=\Delta w/(hw)=\Delta/h$ as well. This makes the lane and
band order parameters directly comparable. In particular,
$\mathbb{E}\Phi_{L}=\mathbb{E}\Phi_{B}\geq(1-p)/(1-p+pm)$ systematically holds
for disordered states.

##  2 Simulation results

We carry out simulations of two-species flows on a $9\times 5$ m rectangular
with top-down and right-left periodic boundary conditions (torus). We simulate
the evolution of $N=45$ agents (density of $1$ agent$/$m2) from random initial
conditions using the first order collision-free (CF) pedestrian model[53] and
the inertial social force (SF) model[54] in the Supplementary Materials. All
agent desired direction motions are polarised to the right. The heterogeneity
in the two parameter settings $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$ is
introduced by varying model parameters related to the speed (i.e. desired
speed or time gap parameters) or to the size of the agents. We quantify the
heterogeneity level in the two-species system using the index $\delta_{s}$
when we vary parameters related to agent speed, and the index $\delta_{l}$
when we vary parameters related to agent size. The definitions of the
microscopic motion model and details on the setting of the model’s parameters
and heterogeneity indexes are provided in Sec. 4.

###  2.1 Preliminary experiment

We first present single simulation histories of the two-species system with
the two heterogeneity models Eq. (1) and Eq. (2). We simulate the evolution of
the agents using the collision-free model[53] for given heterogeneity indexes
$h_{s}$ on the parameters related to the agent speed, namely the desired speed
and the time gap parameters (see Sec. 4 for details on the setting of the
model parameters). Successive snapshots of the system are presented in Fig. 2.
The evolution of the system with the static heterogeneity model Eq. (1) is
shown in the left panels while the evolution with the dynamic heterogeneity
model Eq. (2) is displayed in the right panels. The bottom panels provide the
evolution of the order parameters for lane and band formation. We observe fast
formation of lanes within the first heterogeneity model, while bands emerge
with the second model. The order parameters converge after a transient phase
to stationary performances with lanes or bands where they are polarised to one
or zero. The duration of the transient states is approximately 40 seconds of
simulation. Note that the duration of the transient states varies from a
simulation to another but the system systematically converges to a stationary
state with lanes or bands. Similar performances are observed when using the
social force model instead of the collision-free model (see Fig. S1 in the
Supplementary Materials). Here, the heterogeneity of the two parameter
settings $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$ and corresponding index $h_{s}$
are relatively high. Reducing the heterogeneity index can result in a longer
transient phase or even no formation of lanes and bands. We may expect that
lanes and bands progressively emerge as the heterogeneity index increases.
This is however not the case. As described in the next section, we observe in
stationary states an abrupt phase transition from disorder states to order
states with lanes or bands as the heterogeneity index increases.

Model $M_{1}$ – Static heterogeneity

$t=0,\qquad\bar{v}=0.63,\qquad\Phi_{L}=0.23,\qquad\Phi_{B}=0.14$

$t=10,\qquad\bar{v}=1.36,\qquad\Phi_{L}=0.51,\qquad\Phi_{B}=0.24$

$t=50,\qquad\bar{v}=1.35,\qquad\Phi_{L}=0.97,\qquad\Phi_{B}=0.23$

Model $M_{2}$ – Dynamic heterogeneity

$t=0,\qquad\bar{v}=0.51,\qquad\Phi_{L}=0.19,\qquad\Phi_{B}=0.29$

$t=10,\qquad\bar{v}=1.22,\qquad\Phi_{L}=0.13,\qquad\Phi_{B}=0.5$

$t=50,\qquad\bar{v}=1.28,\qquad\Phi_{L}=0.12,\qquad\Phi_{B}=1$

020406080Time, sOrder parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$020406080Time,
sOrder parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$

Figure 2: Typical screenshots for the model Eq. (1) with heterogeneity in the
agent characteristics for which lanes emerge: the lane order parameter tends
to one while the band order parameter is close to zero (left panels,
$\delta_{s}=18$), and for the model Eq. (2) with heterogeneity in the
interactions where bands emerge with opposite characteristics for the order
parameter (right panels, $\delta_{s}=8$). Flow direction from left to right,
periodic boundary conditions, random initial conditions.

###  2.2 Stationary performances

The preliminary experiment shows that lanes tend to emerge in the dynamics
when the heterogeneity relies on agent characteristics (quenched disorder
model $M_{1}$ Eq. (1)) while bands occur when the heterogeneity takes place in
the interactions (annealed disorder model $M_{2}$ Eq. (2)). The results
presented in Fig. 2 are obtained for given values of the heterogeneity index
$\delta_{s}$ between the two parameter settings $\mathbf{p}_{1}$ and
$\mathbf{p}_{2}$. The index is sufficiently high to rapidly observe the
formation of lanes or bands. In this section, we analyse the performances by
progressively increasing the heterogeneity indexes $\delta_{s}$ and
$\delta_{l}$. We repeated one thousand Monte-Carlo simulations from
independent random initial configurations for the two heterogeneity models
$M_{1}$ Eq. (1) and $M_{2}$ Eq. (2) by varying the heterogeneity indexes
$\delta_{s}$ and $\delta_{l}$ over twenty levels. The differences between the
two parameter settings $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$ are zero at the
lower heterogeneity level, while they are important at the higher level. We
measure the system during $60$ s after a simulation time $t_{0}=600$ s to
observe stationary performances. The remaining Figs. 3–4 and Figs. S2–S5 in
the Supplementary Materials show the median estimates of the Monte-Carlo
simulations with inter-quartile ranges of the agent mean speed and order
parameters for lane and band formation. Details on the setting of the model’s
parameters and heterogeneity indexes are provided in Sec. 4.

Model $M_{1}$Order parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$Model $M_{2}$Order
parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$Heterogeneity index $\delta_{s}$ for the
desiredspeed and time gap parametersMean speed, m/s01201019Type 1Type
2AllHeterogeneity index $\delta_{s}$ for the desiredspeed and time gap
parametersMean speed, m/s01201019 Figure 3: Lane and band order parameters
(top panels) and mean speed (bottom panels) according to the heterogeneity
index $\delta_{s}$ for agent speed features. Longitudinal lanes emerge with
the model $M_{1}$ Eq. (1) for which the heterogeneity relies on agent
characteristics: $\Phi_{L}$ tends to $1$ while $\Phi_{B}$ is low (top left
panel). Orthogonal bands appear with the model $M_{2}$ Eq. (2) for which the
heterogeneity operates in the interaction: $\Phi_{B}$ tends to $1$ while
$\Phi_{L}$ is low (top right panel). The lane pattern allows the speed of the
agents with faster characteristics to be higher than the speed of agents with
slower features (bottom left panel). In contrast, the band pattern acts as a
gridlock for which all the agents have slower speed features (bottom right
panel).

We first present the Monte-Carlo simulation results obtained by varying the
heterogeneity index $\delta_{s}$ relying on agent speed features (desired
speed and time gap parameters). We observe the emergence of lanes when the
heterogeneity relies on agent characteristics (quenched disorder model $M_{1}$
Eq. (1), see Fig. 3, left panels) while band occurs when the heterogeneity
operates in the interactions (annealed disorder model $M_{2}$ Eq. (2), see
Fig. 3, right panels). An abrupt phase transition occurs as the heterogeneity
index $\delta_{s}$ increases from a disordered state for which the order
parameters are close to 0.2 (dotted line in Fig. 3, top panels) to an ordered
dynamics with lanes or bands for which the order parameters are polarised on
zero or one. A critical heterogeneity index can be identified. The lane
patterns allow the speed of the agents with faster characteristics to be
higher than the speed of agents with slower features (Fig. 3, bottom left
panel). This makes the agent speed on average close to the mean speed of a
homogeneous flow (dotted line). In contrast, the band patterns in the model
$M_{2}$ Eq. (2) correspond to gridlocks for which the speed of all the agents
have slower features (Figs. 3, bottom right panel). Similar performances occur
by varying parameters relying on agent size (see Figs. 4 below), or by using
the social force inertial model instead of the first order collision-free
model (see Figs. S2 and S3 in the Supplementary Materials).

The lane formation for counter-flow is an universal collective mechanism in
pedestrian dynamics [36, 37, 38, 32]. It relies on heterogeneity of agent
velocity: some of them have a desired speed $v$ while the others have a
desired speed $-v$. No counter-flow occurs in the presented simulation
experiments since all agents are polarised to the right. The agents of type 1
have a speed $v_{1}>0$ and agents of type 2 have a speed $v_{2}>0$ with
$v_{1}\not=v_{2}$. However, using a Galilean transformation with
$\tilde{v}=(v_{1}+v_{2})/2$ and neglecting anisotropic effects the situation
is the same as in the counter-flow with $v=|v_{2}-v_{1}|/2$. Therefore, it is
not surprising to observe longitudinal lane formation in uni-directional flows
with the quenched disorder model $M_{1}$ Eq. (1) even if no counter-flow
arises (see also recent experiments[50]). The dynamic heterogeneity features
of the annealed disorder model $M_{2}$ Eq. (2) result in transversal
segregation in space and the formation of bands vertical to the direction of
motion. Band structures, as stop-and-go dynamics[57, 44, 45], are self-
organised density waves. However here the agents have a constant speed and the
bands propagate downstream, while the bands propagate upstream with a
characteristic speed and enforce agents to deceleration and acceleration
phases in case of stop-and-go dynamics.

Heterogeneity in agent speed featuring results in transitions from disordered
states to ordered patterns with lanes or bands according to, respectively,
model $M_{1}$ Eq. (1) and model $M_{2}$ Eq. (2). Similar although smother
transitions occur when varying the agent size characteristic index
$\delta_{l}$ instead of the speed index $\delta_{s}$ (compare Figs. 3 and 4).
Lanes emerge for the static heterogeneity model Eq. (1), see Fig. 4, left
panels, while bands arise for the dynamic heterogeneity model Eq. (2), see
Fig. 4, right panels. In contrast to heterogeneous models relying on agent
speed, varying the agent size induces bi-dimensional steric effects making the
average speed in the presence of lanes less than the mean speed of a
homogeneous flow (see Fig. 4, bottom left panel). On contrast, the mean speed
can be higher than the homogeneous one in the presence of bands (see Fig. 4,
bottom right panel). Indeed, varying the agent size acts in two dimensions,
reducing or increasing the available space in case of presence of lanes or
bands. Similar performances occur when using the social force model instead of
the collision-free model (compare Fig. 4 and Fig. S3 in the Supplementary
Materials).

Model $M_{1}$Order parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$Model $M_{2}$Order
parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$Heterogeneity index $\delta_{l}$ for the
agentsize parameterMean speed, m/s01201019Type 1Type 2AllHeterogeneity index
$\delta_{l}$ for the agentsize parameterMean speed, m/s01201019 Figure 4: Lane
and band order parameters (top panels) and mean speed (bottom panels)
according to the heterogeneity index $\delta_{l}$ for agent size. Comparable
although smother transitions to patterns with lanes or bands arise when
varying parameters rela- ted to agent size (compare to Fig. 3) suggesting
universal characteristics of the heterogeneity mechanisms Eq. (1) and Eq. (2).

###  2.3 Transient states and perturbed systems

The simulations above describe stationary situations. Yet, it is interesting
to observe the transient states of the system and the time required for the
emergence of lanes or bands. In Fig. 5, we carry out different simulation
times $t_{0}=0$, $60$, $600$, $1200$ and $3000$ s before starting the
measurements. The initial conditions are random. The lanes and bands
spontaneously emerge during the first minute of simulation when the
heterogeneity index $\delta_{s}$ is sufficiently high. Similar phase
transition to lane and band patterns occur for $t_{0}=600$, $t_{0}=1200$ and
$t_{0}=3000$ s suggesting that the dynamics can be considered stationary as
soon as $t\geq 600$ s. The simulation times required to obtain stationary
performances fluctuate from one simulation to another. They also depend on the
size of the system and the density level. Generally speaking, bigger or more
dense systems require on average longer simulation times to reach a stationary
state than smaller or least dense systems.

Model $M_{1}$Heterogeneity index $\delta_{s}$ for the desiredspeed and time
gap parametersLane order parameter $\Phi_{L}$0.10.50.901019$t_{0}$
(s)06060012003000 Model $M_{2}$Heterogeneity index $\delta_{s}$ for the
desiredspeed and time gap parametersBand order parameter
$\Phi_{B}$0.10.50.901019$t_{0}$ (s)06060012003000  Figure 5: Lane order
parameter for the model $M_{1}$ Eq. (1) (left panel) and band order parameter
for the model $M_{2}$ Eq. (2) (right panel) according to the heterogeneity
index $\delta_{s}$ of agent speed features. The different curves correspond to
different simulation times $t_{0}=0$, $60$, $600$, $1200$ and $3000$ s before
starting the measurements (random initial conditions). The phase transition to
lane and band patterns relatively fast emerges. It can be observed during the
first minutes of simulation. More precisely, the order parameters are similar
for $t_{0}=600$, $1200$ and $3000$ s suggesting the system stationary as soon
as $t\geq 600$ s.

So far, the modelling approach is deterministic. Analysing whether the
collective motion is robust against random noising may reveal unexpected
behaviours. In Fig. 6, we present the order parameter for stochastic systems
for which the agent speeds are subject to independent Brownian noises.
Simulations are carried out for a noise amplitude $\sigma=0.1$, $0.2$ and
$0.5$ m/s. The noise monotonically perturbs the lane formation in the static
heterogeneity model $M_{1}$ Eq. (1) (Fig. 6, left panel). No phase transition
occurs for $\sigma=0.5$ m/s. This phenomenon is well known in the literature
as the freezing-by-heating effect [58]. The concept is borrowed from the plant
growth stimulation process. On the contrary, introducing a low noise in the
dynamics allows improving the band formation in the dynamic heterogeneity
model $M_{2}$ Eq. (2). Indeed, the critical heterogeneity indexes for the
phase transition are smaller with $\sigma=0.1$ or $\sigma=0.2$ m/s than for
the deterministic model with $\sigma=0$, see Fig. 6, right panel. The band
pattern is not only robust to the noise, it exists noise ranges improving the
band formation. In contrast to the freezing-by-heating effect, this can be
interpreted as an example of _noise-induced ordering effect_ [59, 60].

Model $M_{1}$Heterogeneity index $\delta_{s}$ for the desiredspeed and time
gap parametersLane order parameter $\Phi_{L}$0.10.50.901019$\sigma$
(m/s)00.10.20.5 Model $M_{2}$Heterogeneity index $\delta_{s}$ for the
desiredspeed and time gap parametersBand order parameter
$\Phi_{B}$0.10.50.901019$\sigma$ (m/s)00.10.20.5  Figure 6: Lane order
parameter for the model $M_{1}$ Eq. (1) (left panel) and band order parameter
for the model $M_{2}$ Eq. (2) (right panel) according to the heterogeneity
index $\delta_{s}$ of agent speed features. The different curves correspond to
different noise amplitudes $\sigma=0$, $0.1$, $0.2$ and $0.5$ m/s in the
dynamics. The noise clearly perturbs the lane formation (left panel). Such a
phenomenon is well-known in the literature as freezing-by-heating-effect[58].
Oppositely, the noise initially improves the band formation for $\sigma=0.1$
or $0.2$ m/s (noise-induced-ordering-effect), before partly altering it for
$\sigma=0.5$ m/s (right panel).

##  3 Discussion

Two microscopic heterogeneity mechanisms are identified for the formation of
collective segregation in two-species systems of polarised agents. Quenched
disorder model $M_{1}$ Eq. (1) with static heterogeneity in the agent
characteristics initiates the formation of longitudinal lanes. While annealed
disorder model $M_{2}$ Eq. (2) for which the heterogeneity dynamically lies in
the interactions initiates transversal bands perpendicular to the motion
direction. These features are observed when varying parameters related to
agent speed (e.g. desired speed or time gap parameters, see Figs. 2 and 3) or
to agent size (see Fig. 3), and by using the first order collision-free
model[53] or the second order social force model[58] (see Fig. S1–S3 in the
Supplementary Materials). The lanes and bands already emerge rather early in
the simulations (see Figs. 5 and S4 in the Supplementary Materials) and are
partly robust against random perturbations (see Figs. 6 and S5 in the
Supplementary Materials). The band formation is especially robust and may even
be improved by perturbing the agent dynamics (noise-induced ordering effect),
when the lane formation does not (freezing-by-heating effect). Further
simulation results show similar behaviours when the noise is introduced in the
agent polarity (i.e. the desired direction) instead of the speed.

Regardless of the modelling order of the motion models (speed-based or
acceleration-based models) and related parameters, the static quenched
disorder and dynamic annealed disorder heterogeneity mechanisms results in
collective formation of lanes or bands. Relying the heterogeneity on agent
characteristics or on the interactions initiates generically segregation and
the formation of lanes and bands in the dynamics. Such results corroborate the
universality of the lane formation observed in counter-flows [32] and opens
new explanation perspectives for the formation of stripes or bands. Further
theoretical investigations remain necessary to rigorously characterise the
phase transitions. A possibility is to analyse mean-field instability
phenomena of discrete lattice representations of the model [41]. The presence
of walls and obstacles and the role of the geometry of given facilities may
also be of interest. Preliminary simulation results show segregation effects
of slower or bigger agents in case of bottleneck within the static
heterogeneity model. These are expelled at the edges of the system and are
obstructed by the presence of walls. These simulation results require more
attention, notably for elderly people, people with motor disabilities, or in
the current context of social distanciation.

##  4 Methods

The two agent motion models used in the simulations are the collision-free
(CF) model[53] and the social force (SF) model[54] (see the Supplementary
Materials). The CF model is a speed-based model of first order while the SF
model is an inertial acceleration-based model. The agent dynamics are
polarised in both models, i.e. all the agents have identical desired
direction.

#### Collision-free model

In the collision-free model, the dynamics of an agent $n$ with position
$\mathbf{x}_{n}$ is given by the first order differential equation

$\dot{\mathbf{x}}_{n}=F^{\textsc{\tiny
cf}}_{\mathbf{p}_{j}}(\mathbf{X}_{n})=V(\mathbf{X}_{n},p_{j})\;\mathbf{e}(\mathbf{X}_{n},p_{j})+\sigma\xi_{n},$
(5)

composed of the scalar speed model

$V(\mathbf{X}_{n},p_{j})=\max\big{\\{}0,\min\big{\\{}V_{j},\big{(}s(\mathbf{X}_{n})-\ell_{j}\big{)}/T_{j}\big{\\}}\big{\\}},$
(6)

here $V_{j}\geq 0$ is the desired speed, $T_{j}>0$ denotes the desired time
gap, and $\ell_{j}\geq 0$ the agent size, the index $j=1,2$ representing the
two parameter setting $\mathbf{p}_{1}$ and $\mathbf{p}_{2}$, and the direction
model

$\mathbf{e}(\mathbf{X}_{n},p_{j})=\mathbf{e}_{0}+a\sum_{m\not=n}\exp\big{(}(\ell_{j}-\|\mathbf{x}_{n}-\mathbf{x}_{m}\|)B^{-1}\big{)}\mathbf{e}_{nm},$
(7)

with $\mathbf{e}_{0}=0$ the desired direction (polarity), $a=5$ and $B=0.1$ m
parameters for the repulsion with the neighbours, $\mathbf{e}_{nm}$ the
direction from agent $m$ to agent $n$ and $C>0$ a normalisation constant. A
bi-dimensional white noise $\xi_{n}$ (i.e. the time derivatives of two
independent Wiener processes) with amplitude $\sigma>0$ is used for the
stochastic model in Sec. 2.3.

The function
$s(\mathbf{X}_{n})=\|\mathbf{x}_{n}-\mathbf{x}_{m_{0}(\mathbf{X}_{n})}\|$ in
the scalar speed model determines the minimal distance in front,

$m_{0}(\mathbf{X}_{n})=\text{arg}\min_{m\in
N_{n}}\|\mathbf{x}_{n}-\mathbf{x}_{m}\|,\qquad
N_{n}=\big{\\{}m,\leavevmode\nobreak\ \mathbf{e}_{n}\cdot\mathbf{e}_{nm}\leq
0\ \mbox{ and }\
|\mathbf{e}_{n}^{\perp}\cdot\mathbf{e}_{nm}|\leq\ell/\|\mathbf{x}_{n}-\mathbf{x}_{m}\|\big{\\}},$
(8)

being the closest agent in front of the agent $n$. Note that in the definition
of the dynamic heterogeneity type Eq. (2), the type of the closest agent in
front is $\tilde{k}(\mathbf{X}_{n})=k_{m_{0}(\mathbf{X}_{n})}$.

The simulations are carried out using an explicit Euler numerical scheme in
deterministic cases, and using an Euler-Maruyama scheme for the simulations
including the stochastic noise. The time step is $\delta t=0.01$ s in both
cases.

#### Setting of the parameters

The default values for the parameters $\mathbf{p}=(\ell,V,T)$ of the CF model
are based on the setting proposed in the literature[53]

$\ell=0.3\leavevmode\nobreak\ \text{m},\quad V=1.5\leavevmode\nobreak\
\text{m/s},\quad T=1\leavevmode\nobreak\ \text{s}.$ (9)

Note that $\Delta=0.6$ m in the order parameters corresponds approximately to
two times the size of a pedestrian. Starting from the default values, we vary
using heterogeneity indexes the parameter settings
$\mathbf{p}_{1}=(\ell_{1},V_{1},T_{1})$ and
$\mathbf{p}_{2}=(\ell_{2},V_{2},T_{2})$.

  * •

In the analysis of the speed heterogeneity (heterogeneity index $\delta_{s}$,
cf. Figs. 2, 3, 5 and 6) the time gap parameter $T$ ranges into $[0.05,1.95]$
s by step of $0.05$ s,

$T_{1}=T+0.05\delta_{s},\qquad
T_{2}=T-0.05\delta_{s},\qquad\qquad\delta_{s}=0,\ldots,19,$ (10)

while the desired (maximal) speed $V$ ranges into $[1,2]$ m/s by step of
$0.025$ m/s,

$V_{1}=V-0.025\delta_{s},\qquad
V_{2}=V+0.025\delta_{s},\qquad\qquad\delta_{s}=0,\ldots,19.$ (11)

The parameters $\ell_{1}=\ell_{2}=\ell=0.3$ m for the agent size remain
constant.

  * •

In the analysis of the size heterogeneity (heterogeneity index $\delta_{l}$,
cf. Fig. 4) the parameter $\ell$ ranges into $[0,0.9]$ m by step of $0.015$ m
decreasing and $0.03$ m increasing,

$\ell_{1}=\ell-0.015\delta_{l},\qquad\ell_{2}=\ell+0.03\delta_{l},\qquad\qquad\delta_{l}=0,\ldots,19.$
(12)

The remaining parameters for the agent speed $V_{1}=V_{2}=V=1.5$ m/s and
$T_{1}=T_{2}=T=1$ s are constant.

## References

  * [1] Castellano, C., Fortunato, S. & Loreto, V.  Statistical physics of social dynamics.  _Rev. Mod. Phys._ 81, 591–646 (2009).
  * [2] Ramaswamy, S.  The mechanics and statistics of active matter.  _Annual Review of Condensed Matter Physics_ 1, 323–345 (2010).
  * [3] Vicsek, T. & Zafeiris, A.  Collective motion.  _Physics Reports_ 517, 71–140 (2012).
  * [4] Elgeti, J., Winkler, R. G. & Gompper, G.  Physics of microswimmers—single particle motion and collective behavior: a review.  _Reports on Progress in Physics_ 78, 056601 (2015).
  * [5] Friedl, P. & Gilmour, D.  Collective cell migration in morphogenesis, regeneration and cancer.  _Nature reviews Molecular cell biology_ 10, 445–457 (2009).
  * [6] Chen, C., Liu, S., Shi, X.-q., Chaté, H. & Wu, Y.  Weak synchronization and large-scale collective oscillation in dense bacterial suspensions.  _Nature_ 542, 210–214 (2017).
  * [7] Turing, A. M.  The chemical basis of morphogenesis.  _Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences_ 237, 37–72 (1952).
  * [8] Vicsek, T., Czirók, A., Ben-Jacob, E., Cohen, I. & Shochet, O.  Novel type of phase transition in a system of self-driven particles.  _Phys. Rev. Lett._ 75, 1226–1229 (1995).
  * [9] Parrish, J. K. & Edelstein-Keshet, L.  Complexity, pattern, and evolutionary trade-offs in animal aggregation.  _Science_ 284, 99–101 (1999).
  * [10] Buhl, J. _et al._ From disorder to order in marching locusts.  _Science_ 312, 1402–1406 (2006).
  * [11] Needleman, D. & Dogic, Z.  Active matter at the interface between materials science and cell biology.  _Nature reviews materials_ 2, 1–14 (2017).
  * [12] Narayan, V., Ramaswamy, S. & Menon, N.  Long-lived giant number fluctuations in a swarming granular nematic.  _Science_ 317, 105–108 (2007).
  * [13] Giomi, L.  Geometry and topology of turbulence in active nematics.  _Phys. Rev. X_ 5, 031003 (2015).
  * [14] Bär, M., Großmann, R., Heidenreich, S. & Peruani, F.  Self-propelled rods: Insights and perspectives for active matter.  _Annual Review of Condensed Matter Physics_ 11, 441–466 (2020).
  * [15] Janssen, L. M., Kaiser, A. & Löwen, H.  Aging and rejuvenation of active matter under topological constraints.  _Scientific reports_ 7, 1–13 (2017).
  * [16] Jahanshahi, S., Lozano, C., ten Hagen, B., Bechinger, C. & Löwen, H.  Colloidal brazil nut effect in microswimmer mixtures induced by motility contrast.  _The Journal of Chemical Physics_ 150, 114902 (2019).
  * [17] Grauer, J., Löwen, H., Be’er, A. & Liebchen, B.  Swarm hunting and cluster ejections in chemically communicating active mixtures.  _Scientific reports_ 10, 1–11 (2020).
  * [18] Ibele, M., Mallouk, T. & Sen, A.  Schooling behavior of light-powered autonomous micromotors in water.  _Angewandte Chemie International Edition_ 48, 3308–3312 (2009).
  * [19] Ozkan-Aydin, Y., Goldman, D. I. & Bhamla, M. S.  Collective dynamics in entangled worm and robot blobs.  _Proceedings of the National Academy of Sciences_ 118 (2021).
  * [20] Cates, M. E. & Tailleur, J.  Motility-induced phase separation.  _Annual Review of Condensed Matter Physics_ 6, 219–244 (2015).
  * [21] Caprini, L., Marini Bettolo Marconi, U. & Puglisi, A.  Spontaneous velocity alignment in motility-induced phase separation.  _Phys. Rev. Lett._ 124, 078001 (2020).
  * [22] Chowdhury, D., Santen, L. & Schadschneider, A.  Statistical physics of vehicular traffic and some related systems.  _Phys. Rep._ 329, 199–329 (2000).
  * [23] Bellomo, N., Piccoli, B. & Tosin, A.  Modeling crowd dynamics from a complex system viewpoint.  _Mathematical Models and Methods in Applied Sciences_ 22, 1230004 (2012).
  * [24] Moussaïd, M. _et al._ Traffic instabilities in self-organized pedestrian crowds.  _PLOS Comput. Biol._ 8, 1–10 (2012).
  * [25] Shahhoseini, Z. & Sarvi, M.  Collective movements of pedestrians: How we can learn from simple experiments with non-human (ant) crowds.  _PLOS ONE_ 12, 1–20 (2017).
  * [26] Boltes, M., Zhang, J., Tordeux, A., Schadschneider, A. & Seyfried, A.  _Empirical Results of Pedestrian and Evacuation Dynamics_ , 1–29 (Springer, Berlin, Heidelberg, 2018).
  * [27] Hermann, G. & Touboul, J.  Heterogeneous connections induce oscillations in large-scale networks.  _Phys. Rev. Lett._ 109, 018702 (2012).
  * [28] Muchnik, L., Aral, S. & Taylor, S. J.  Social influence bias: A randomized experiment.  _Science_ 341, 647–651 (2013).
  * [29] Moussaïd, M., Kämmer, J. E., Analytis, P. P. & Neth, H.  Social influence and the collective dynamics of opinion formation.  _PLOS ONE_ 8, 1–8 (2013).
  * [30] Touboul, J. D.  The hipster effect: When anti-conformists all look the same.  _Discrete & Continuous Dynamical Systems - B_ 24, 4379 (2019).
  * [31] Fu, X. _et al._ Spatial self-organization resolves conflicts between individuality and collective migration.  _Nature communications_ 9, 1–12 (2018).
  * [32] Cristìn, J., Méndez, V. & Campos, D.  General scaling in bidirectional flows of self-avoiding agents.  _Sci. Rep._ 9, 18488 (2019).
  * [33] Childress, S., Levandowsky, M. & Spiegel, E. A.  Pattern formation in a suspension of swimming microorganisms: equations and stability theory.  _Journal of Fluid Mechanics_ 69, 591–613 (1975).
  * [34] Helbing, D.  Traffic and related self-driven many-particle systems.  _Phys. Mod. Phys._ 73, 1067–1141 (2001).
  * [35] Kirchner, A. & Schadschneider, A.  Simulation of evacuation processes using a bionics-inspired cellular automaton model for pedestrian dynamics.  _Physica A: Statistical Mechanics and its Applications_ 312, 260–276 (2002).
  * [36] Burstedde, C., Klauck, K., Schadschneider, A. & Zittartz, J.  Simulation of pedestrian dynamics using a two-dimensional cellular automaton.  _Physica A_ 295, 507–525 (2001).
  * [37] Nakayama, A., Hasebe, K. & Sugiyama, Y.  Instability of pedestrian flow and phase structure in a two-dimensional optimal velocity model.  _Phys. Rev. E_ 71, 036121 (2005).
  * [38] Feliciani, C., Murakami, H. & Nishinari, K.  A universal function for capacity of bidirectional pedestrian streams: Filling the gaps in the literature.  _PLOS ONE_ 13, 1–31 (2018).
  * [39] Garcimartín, A., Pastor, J. M., Martín-Gómez, C., Parisi, D. & Zuriguel, I.  Pedestrian collective motion in competitive room evacuation.  _Sci. Rep._ 7, 1–9 (2017).
  * [40] Helbing, D., Buzna, L., Johansson, A. & Werner, T.  Self-organized pedestrian crowd dynamics: Experiments, simulations, and design solutions.  _Transportation Science_ 39, 1–24 (2005).
  * [41] Cividini, J., Appert-Rolland, C. & Hilhorst, H.-J.  Diagonal patterns and chevron effect in intersecting traffic flows.  _EPL_ 102, 20002 (2013).
  * [42] Schwarting, W., Pierson, A., Alonso-Mora, J., Karaman, S. & Rus, D.  Social behavior for autonomous vehicles.  _Proceedings of the National Academy of Sciences_ 116, 24972–24978 (2019).
  * [43] Cerotti, D., Distefano, S., Merlino, G. & Puliafito, A.  A crowd-cooperative approach for intelligent transportation systems.  _IEEE Transactions on Intelligent Transportation Systems_ 18, 1529–1539 (2017).
  * [44] Stern, R. E. _et al._ Dissipation of stop-and-go waves via control of autonomous vehicles: Field experiments.  _Transportation Research Part C: Emerging Technologies_ 89, 205–221 (2018).
  * [45] Bain, N. & Bartolo, D.  Dynamic response and hydrodynamics of polarized crowds.  _Science_ 363, 46–49 (2019).
  * [46] Friesen, M., Gottschalk, H., Rüdiger, B. & Tordeux, A.  Spontaneous wave formation in stochastic self-driven particle systems.  _SIAM J. Appl. Math._ 81, 853–870 (2021).
  * [47] Zuriguel, I. _et al._ Clogging transition of many-particle systems flowing through bottlenecks.  _Sci. Rep._ 4, 7324 (2014).
  * [48] Nicolas, A., Bouzat, S. & Kuperman, M. N.  Pedestrian flows through a narrow doorway: Effect of individual behaviours on the global flow and microscopic dynamics.  _Transportation Research Part B: Methodological_ 99, 30–43 (2017).
  * [49] Foulaadvand, M. E. & Neek-Amal, M.  Asymmetric simple exclusion process describing conflicting traffic flows.  _Europhysics Letters (EPL)_ 80, 60002 (2007).
  * [50] Fujita, A., Feliciani, C., Yanagisawa, D. & Nishinari, K.  Traffic flow in a crowd of pedestrians walking at different speeds.  _Phys. Rev. E_ 99, 062307 (2019).
  * [51] Krüsemann, H., Godec, A. & Metzler, R.  First-passage statistics for aging diffusion in systems with annealed and quenched disorder.  _Physical Review E_ 89, 040101 (2014).
  * [52] Tateishi, A., Ribeiro, H., Sandev, T., Petreska, I. & Lenzi, E.  Quenched and annealed disorder mechanisms in comb models with fractional operators.  _Physical Review E_ 101, 022135 (2020).
  * [53] Tordeux, A., Chraibi, M. & Seyfried, A.  Collision-free speed model for pedestrian dynamics.  In _Traffic and Granular Flow ’15_ , 225–232 (Springer, Cham, 2016).
  * [54] Helbing, D. & Molnár, P.  Social force model for pedestrian dynamics.  _Phys. Rev. E_ 51, 4282–4286 (1995).
  * [55] Rex, M. & Löwen, H.  Lane formation in oppositely charged colloids driven by an electric field: Chaining and two-dimensional crystallization.  _Phys. Rev. E_ 75, 051402 (2007).
  * [56] Nowak, S. & Schadschneider, A.  Quantitative analysis of pedestrian counterflow in a cellular automaton model.  _Phys. Rev. E_ 85, 066128 (2012).
  * [57] Nakayama, A., Hasebe, K. & Sugiyama, Y.  Effect of attractive interaction on instability of pedestrian flow in a two-dimensional optimal velocity model.  _Physical Review E_ 77, 016105 (2008).
  * [58] Helbing, D., Farkas, I. & Vicsek, T.  Freezing by heating in a driven mesoscopic system.  _Phys. Rev. Lett._ 84, 1240–1243 (2000).
  * [59] Wiesenfeld, K., Pierson, D., Pantazelou, E., Dames, C. & Moss, F.  Stochastic resonance on a circle.  _Physical Review Letters_ 72, 2125 (1994).
  * [60] D’Huys, O., Veltz, R., Dolcemascolo, A., Marino, F. & Barland, S.  Canard resonance: on noise-induced ordering of trajectories in heterogeneous networks of slow-fast systems.  _Journal of Physics: Photonics_ 3, 024010 (2021).

## Acknowledgements

RK and AT acknowledge the Franco-German research project MADRAS funded in
France by the Agence Nationale de la Recherche (ANR, French National Research
Agency), project number ANR-20-CE92-0033, and in Germany by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation), project number
446168800.

## Data availability

Simulations of lane or band formations with the two heterogeneity models
$M_{1}$ Eq. (1) and $M_{2}$ Eq. (2) can be implemented online and in real time
on the dedicated web-page https://www.vzu.uni-
wuppertal.de/fileadmin/site/vzu/Lane-and- band-formation. Simulation module
and programming code (in Logo language) can be downloaded on the same link.

## Competing interests

The authors declare no competing financial interests.

## Author contributions

AS and AT conceived and designed the study. BK and RK performed the numerical
simulations. AT wrote the manuscript. All authors reviewed the article.

Supplementary Materials

We use the collision-free (CF) model in the manuscript to analyse by
simulation the two heterogeneity types $M_{1}$ Eq. (1) and $M_{2}$ Eq. (2). In
these supplementary materials, we present identical simulation experiments as
those carried out in the manuscript by using the Social Force (SF) model[54].
Generally speaking, similar phase transitions to lane and band patterns occur
when using, respectively, the static heterogeneity models $M_{1}$ or the
dynamic model $M_{2}$. In these Supplementary Materials, we first present the
simulation results before providing technical details on the social force
model.

#### Simulation results

We simulate using the social force model[54] the evolution of 45 agents on a
$9\times 5$ m rectangle with periodic boundaries identically to the
simulations presented in the manuscript. The preliminary experiment presents
two system histories obtained with the two heterogeneity models for given
heterogeneity indexes. As with CF model, we observe using the SF model rapid
formation of lanes with the static heterogeneity type $M_{1}$, while band
patterns emerge with the dynamic heterogeneity type $M_{2}$ (compare Fig. 2 in
the manuscript to Fig. 7).

Comparable phase transitions occur in stationary states as the heterogeneity
indexes increase, from disorder states to ordered states with lane or band
patterns and polarised order parameters, see Fig. 8. This holds for
heterogeneity relying on agent speed features (heterogeneity index
$\delta_{s}$, see Fig. 8) or on the agent size (heterogeneity index
$\delta_{l}$, see Fig. 9). However the transition to lane pattern is specially
laborious when dealing with the size parameter with SF model (see Fig. 9, top
left panel). Note that the size parameter $\ell$ with SF model does not
describe a hard-core exclusion between the agents as CF model does. On the
other hand the phase transitions to band patterns (heterogeneity type $M_{2}$)
occur for a lower heterogeneity feature again with SF model (compare Figs. 3
and 4 in the manuscript to Figs. 8 and 9).

The phase transitions to lane and band patterns can be observed during the
first minutes of simulations with SF model (see Fig. 10). However, in contrast
to the results obtained with CF model, the system is not completely stabilised
after 600 s, especially for the heterogeneity type $M_{1}$ (compare Fig. 5 in
the manuscript to Fig. 10). SF model, being an inertial model of the second
order, requires longer simulation times to describe stationary performances.

Similarly to CF model, noising the dynamic clearly perturbs the lane formation
(freezing-by-heating-effect, compare Fig. 6 in the manuscript to Fig. 11, left
panels). While oppositely, the band formation is robust against the noise (see
Fig. 11, right panel).

Model $M_{1}$ – Static heterogeneity

$t=0,\qquad\bar{v}=0.24,\qquad\Phi_{L}=0.31,\qquad\Phi_{B}=0.33$

$t=10,\qquad\bar{v}=1.39,\qquad\Phi_{L}=0.26,\qquad\Phi_{B}=0.19$

$t=50,\qquad\bar{v}=1.27,\qquad\Phi_{L}=1,\qquad\Phi_{B}=0.12$

Model $M_{2}$ – Dynamic heterogeneity

$t=0,\qquad\bar{v}=0.2,\qquad\Phi_{L}=0.22,\qquad\Phi_{B}=0.1$

$t=10,\qquad\bar{v}=1.27,\qquad\Phi_{L}=0.17,\qquad\Phi_{B}=0.21$

$t=50,\qquad\bar{v}=1.34,\qquad\Phi_{L}=0.1,\qquad\Phi_{B}=1$

020406080Time, sOrder parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$020406080Time,
sOrder parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$

Figure 7: Typical screenshots for the model $M_{1}$ with heterogeneity in the
agent characteristics for which lanes emerge: the lane order parameter
$\Phi_{L}$ tends to one while the band order parameter $\Phi_{B}$ is close to
zero (left panels, $\delta_{s}=18$), and for the model $M_{2}$ with
heterogeneity in the interactions where bands emerge, $\Phi_{L}$ is close to
zero while $\Phi_{B}$ tends to one (right panels, $\delta_{s}=10$). SF motion
model, flow direction from left to right, periodic boundary conditions, random
initial conditions. Model $M_{1}$Order
parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$Model $M_{2}$Order
parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$Heterogeneity index $\delta_{s}$ for the
desiredspeed and repulsion rate parametersMean speed, m/s01201019Type 1Type
2AllHeterogeneity index $\delta_{s}$ for the desiredspeed and repulsion rate
parametersMean speed, m/s01201019 Figure 8: Lane and band order parameters
(top panels) and mean speed (bottom panels) according to the heterogeneity
index of agent speed features for the static heterogeneity model $M_{1}$ (left
panels) and for the dynamic heterogeneity model $M_{2}$ (right panels) with
the SF model. We observe qualitatively similar phase transitions as those
obtained with the CF model, compare with Fig. 3 in the manuscript. Model
$M_{1}$Order parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$Model $M_{2}$Order
parameter0.00.51.0$\Phi_{L}$$\Phi_{B}$Heterogeneity index $\delta_{l}$ for the
agentsize parameterMean speed, m/s01201019Type 1Type 2AllHeterogeneity index
$\delta_{l}$ for the agentsize parameterMean speed, m/s01201019 Figure 9: Lane
and band order parameters (top panels) and mean speed (bottom panels)
according to the heterogeneity index of agent size for the SF model. In
contrast to CF model (see Fig. 4 in the manuscript), the transition to lane
patterns is slower in absence of strict exclusion rules between the agents.
Model $M_{1}$Heterogeneity index $\delta_{s}$ for the desiredspeed and time
gap parametersLane order parameter $\Phi_{L}$0.10.50.901019$t_{0}$
(s)06060012003000 Model $M_{2}$Heterogeneity index $\delta_{s}$ for the
desiredspeed and time gap parametersBand order parameter
$\Phi_{B}$0.10.50.901019$t_{0}$ (s)06060012003000  Figure 10: Lane order
parameter for the static heterogeneity model $M_{1}$ (left panel), and band
order parameter for the dynamic heterogeneity model $M_{2}$ (right panel)
according to the heterogeneity index $\delta_{s}$ of agent speed features. The
different curves correspond to different simulation times $t_{0}=0$, $60$,
$600$, $1200$ and $3000$ s before starting the measurements (random initial
conditions). As for the CF model, the phase transition to lane and band
patterns relatively fast emerges with the SF model. It can be observed during
the first minutes of simulation. However, the system is not completely
stabilised after 600 s. Model $M_{1}$Heterogeneity index $\delta_{s}$ for the
desiredspeed and time gap parametersLane order parameter
$\Phi_{L}$0.10.50.901019$\sigma$ (m/s)00.050.100.20 Model $M_{2}$Heterogeneity
index $\delta_{s}$ for the desiredspeed and time gap parametersBand order
parameter $\Phi_{B}$0.10.50.901019$\sigma$ (m/s)00.050.100.20  Figure 11: Lane
order parameter for the static heterogeneity model $M_{1}$ (left panel), and
band order parameter for the dynamic heterogeneity model $M_{2}$ (right panel)
according to the heterogeneity index $\delta_{s}$ of agent speed features. The
different curves correspond to different noise amplitudes $\sigma=0$, $0.1$,
$0.2$ and $0.5$ m/s in the dynamics. The noise clearly perturbs the lane
formation (left panel, freezing-by-heating-effect). Oppositely, the band
formation is robust against the noise (right panel).

#### Definition of the social force model

In the social force model, the dynamics of an agent $n$ with position
$\mathbf{x}_{n}$ and neighborhood $\mathbf{X}_{n}$ is given by the second
order differential equation

$\ddot{\mathbf{x}}_{n}=F^{\textsc{\tiny
sf}}_{\mathbf{p}_{j}}(\mathbf{X}_{n})=\frac{1}{\tau}(V_{j}\,\mathbf{e}_{0}-\dot{\mathbf{x}}_{n})+A_{j}\sum_{m\not=n}\varphi(\mathbf{e}_{mn})\exp\big{(}(\ell_{j}-\|\mathbf{x}_{n}-\mathbf{x}_{m}\|)B^{-1}\big{)}\mathbf{e}_{nm}+\sigma\xi_{n}.$
(13)

Here $V_{j}\geq 0$ is the desired speed, $\ell_{j}\geq 0$ the agent size and
$\tau=0.5$ s a relaxation time while $\mathbf{e}_{mn}$ is the direction from
$m$ to $n$, $\varphi(\mathbf{e})=1-\cos(\pi-\hat{\mathbf{e}})$ is the vision
field factor and $A_{j}$ and $B=0.2$ m are parameters for the repulsion with
the neighbours. As for CF model, $\mathbf{e}_{0}=0$ is the desired direction
(polarity). A bi-dimensional white noise $\xi_{n}$ (i.e. the time derivatives
of two independent Wiener processes) with amplitude $\sigma>0$ is used in Fig.
11.

The model is simulated using an explicit Euler numerical scheme in
deterministic cases, and using an Euler-Maruyama scheme for the simulation
including a stochastic noise. The time step is $\delta t=0.01$ s in both
cases.

#### Setting of the parameters

The default values for the parameters $\mathbf{p}=(\ell,V,A)$ of the SF model
are based on the setting proposed in the literature[54]

$\ell=0.3\leavevmode\nobreak\ \text{m},\qquad V=1.5\leavevmode\nobreak\
\text{m/s},\qquad A=3\leavevmode\nobreak\ \text{m/s${}^{2}$}.$ (14)

Starting from the default values, we vary using heterogeneity indexes the
parameter settings $\mathbf{p}_{1}=(\ell_{1},V_{1},A_{1})$ and
$\mathbf{p}_{2}=(\ell_{2},V_{2},A_{2})$.

  * •

In the analysis of the speed heterogeneity (heterogeneity level $\delta_{s}$,
cf. Figs. 7, 8, 10 and 11) the expulsion rate $A$ ranges into $[2,4]$ m/s2 by
step of $0.05$ m/s2,

$A_{1}=A-0.05\delta_{s},\qquad
A_{2}=A+0.05\delta_{s},\qquad\qquad\delta_{s}=0,\ldots,19,$ (15)

while the desired (maximal) speed $V$, as with CF model, ranges into $[1,2]$
m/s by step of $0.025$ m/s,

$V_{1}=V-0.025\delta_{s},\qquad
V_{2}=V+0.025\delta_{s},\qquad\qquad\delta_{s}=0,\ldots,19.$ (16)

Note that $\tau_{j}=V_{j}/A_{j}$ systematically holds for all $j=1,2$. The
parameter $\ell_{1}=\ell_{2}=\ell$ for the agent size remains constant.

  * •

In the analysis of the size heterogeneity (heterogeneity level $\delta_{l}$,
cf. Fig. 9) the parameter $\ell$, as with CF model, ranges into $[0,0.9]$ m by
step of $0.015$ m decreasing and $0.03$ m increasing,

$\ell_{1}=\ell-0.015\delta_{l},\qquad\ell_{2}=\ell+0.03\delta_{l}\qquad\qquad\delta_{l}=0,\ldots,19.$
(17)

The remaining parameters for the agent speed $V_{1}=V_{2}=V$, and
$A_{1}=A_{2}=A$ are constant.


# Optimizing working parameters of the twin-range cutoff method in terms of
accuracy and efficiency

Han Wang  LMAM and School of Mathematical Sciences, Peking University, Beijing
<EMAIL_ADDRESS>Pingwen Zhang  LMAM and School of Mathematical Sciences,
Peking University, Beijing

###### Abstract

We construct a priori error estimation for the force error of the twin-range
cutoff method, which is widely used to treat the short-range non-bonded
interactions in molecular simulations. Based on the error and cost estimation,
we develop a work flow that can automatically determine the nearly most
efficient twin-range cutoff parameters (i.e. the cutoff radii and the neighbor
list updating frequency) prior to a simulation for a predetermined accuracy.
Both the error estimate and the parameter tuning method are demonstrated to be
effective by testing simulations of the standard Lennard-Jones 6-12 fluid in
gas, liquid as well as supercritical state. We recommend the tuned twin-range
cutoff method that can save precious user time and computational resources.

##  I Introduction

Non-bonded interactions are encountered in nearly every molecular simulation,
but their calculations are computationally expensive. It is thus important to
develop computational methods that boost both efficiency and accuracy at the
same time. Non-bonded interactions are mainly formed by a pairwise interaction
$u(r)$, where $r$ is the distance between two interacting particles. There are
two types of pairwise interactions, namely the long-range interaction and the
short-range interaction, depending on the rate at which $u(r)$ decay with
respect to the distance.

Most short-range interactions satisfy $|u(r)|\leq Cr^{-m}$, $m>3$, which
guarantees an absolute convergence of the energy. A naive idea to treat short-
range interactions is to explicitly calculate and sum all pairwise
interactions. This results in a computational cost scaling
$\mathcal{O}(N^{2})$ per time step, which rapidly becomes inefficient as the
number of particles $N$ grows large. A better way is to introduce a cutoff
radius, outside of which all pairwise interactions are simply neglected. In
combination with the cell list and the neighbor list algorithms frenkel02b ,
the total computational cost of the short-range interactions can be reduced to
an acceptable level of $\mathcal{O}(N)$.

It is well known that various physical properties show significant dependence
on the cutoff radius and the method to treat the discontinuity at the cutoff,
for example, the phase diagram of the Lennard-Jones fluid smit1992phase ;
baidakov2000effect ; ou2005molecular , the density profile of the liquid-vapor
interface panagiotopoulos1994molecular ; mecke1997molecular , the surface
tension trokhymchuk1999computer ; baidakov2000effect ; shen2007comparative
and the free energy calculation ou2005molecular . Ill-chosen cutoff can lead
to undesirable artifacts, for example, the phase diagrams of the Lennard-Jones
fluid have been demonstrated to be substantially different with different
choices of cutoff radii smit1992phase . A straightforward way to eliminate the
cutoff effects is to use an extremely large cutoff radius such that all
properties of interest satisfactorily converge. This has motivated the
adoption of a non-truncated potential (i.e. the cutoff radius is the same as
half of the simulation box) to study the critical point phenomena
lotfi1992vapour ; caillol1998critical ; potoff1998critical ; perez2006critical
. However, such a long cutoff radius drastically increases the computational
cost, thereby prohibiting long time and large size molecular simulations. An
alternative way is to use a small cutoff while applying long range correction
(LRC) allen87a  to eliminate the systematic error of the potential energy, the
pressure and the free energy. The idea is to integrate the thermodynamic
properties from the cutoff radius to infinity, assuming the radial
distribution function is equal to 1. Some more sophisticated methods have been
developed in recent years to improve LRC for the constant pressure simulation
lague2004pressure , the inhomogeneous systems janecek2006long  and the free
energy calculation of systems containing macromolecules shirts2007accurate .

A promising way to quantitatively analyze the undesirable cutoff effects is to
express these artifacts in terms of the difference between the cutoffed
interaction and the exact interaction. There are several error analysis works
on the long-range interaction algorithms kolafa1992cee ; hummer95a ;
petersen1995accuracy ; deserno1998mue2 ; wang2010optimizing , while the
analysis on the short-range interaction is scarce. In this paper, we develop
the error estimate of the short-range force introduced by the widely used
twin-range cutoff method (which is reduce the commonly used single-range
cutoff method by a special choice of parameters) in homogeneous systems.
Furthermore, to automatically determine the most efficient cutoff radii and
neighbor list updating frequency, a work flow optimizing these working
parameters with respect to speed and accuracy is proposed and tested.

##  II Theory

###  II.1 The twin range cutoff method

The twin-range cutoff method introduces two cutoffs, denoted by $r_{1}$ and
$r_{2}$, with $r_{1}\leq r_{2}$. Assuming the neighbor list updating frequency
is $M$, then at every $M$ steps, the neighbor list is built for all
neighboring particles that fall in the short cutoff radius $r_{1}$ and the
corresponding interactions are calculated and applied every next $M-1$ steps.
The neighbors fall in between the short cutoff and the long cutoff are assumed
to move slowly, so the interactions can be calculated less frequently. At
every $M$ steps, the interactions with these particles are calculated and
stored. In the subsequent $M-1$ steps, they are applied to the corresponding
particles without any change. When $r_{1}=r_{2}$, the twin-range cutoff method
reduces to the single-range cutoff method.

###  II.2 The error estimate of the twin-range cutoff method

For simplicity, we will only study the error estimate of single component
systems, because it is not difficult to derive the error estimate for
multicomponent systems in the same way. We consider a reference particle $i$
and denote the sets of neighbors fall in $r_{1}$, between $r_{1}$ and $r_{2}$,
and out of $r_{2}$ by $\Omega_{1}^{i}$, $\Omega_{2}^{i}$ and $\Omega_{3}^{i}$,
respectively. The position of any particle $j$ is denoted by
$\textbf{{r}}_{j}$ at the step when the neighbor list is built. In the
following several steps, the accumulated absolute displacement is denoted by
$\textbf{{d}}_{j}$. Then the exact force $\textbf{{F}}_{i}^{\ast}$ on
particles $i$ and the cutoffed force $\textbf{{F}}_{i}$ by the twin-range
cutoff method are

$\displaystyle\textbf{{F}}_{i}^{\ast}$
$\displaystyle=\sum_{j\in\Omega_{1}^{i}}\textbf{{F}}(\textbf{{r}}_{ij}+\textbf{{d}}_{ij})+\sum_{j\in\Omega_{2}^{i}}\textbf{{F}}(\textbf{{r}}_{ij}+\textbf{{d}}_{ij})+\sum_{j\in\Omega_{3}^{i}}\textbf{{F}}(\textbf{{r}}_{ij}+\textbf{{d}}_{ij}),$
(1) $\displaystyle\textbf{{F}}_{i}$
$\displaystyle=\sum_{j\in\Omega_{1}^{i}}\textbf{{F}}(\textbf{{r}}_{ij}+\textbf{{d}}_{ij})+\sum_{j\in\Omega_{2}^{i}}\textbf{{F}}(\textbf{{r}}_{ij}).$
(2)

Therein the relative position is
$\textbf{{r}}_{ij}=\textbf{{r}}_{i}-\textbf{{r}}_{j}$ and the relative
displacement is $\textbf{{d}}_{ij}=\textbf{{d}}_{i}-\textbf{{d}}_{j}$.

The difference between the exact force and the cutoffed force is therefore

$\displaystyle\Delta\textbf{{F}}$
$\displaystyle=\textbf{{F}}_{i}^{\ast}-\textbf{{F}}_{i}$
$\displaystyle=\sum_{j\in\Omega_{2}^{i}}[\,\textbf{{F}}(\textbf{{r}}_{ij}+\textbf{{d}}_{ij})-\textbf{{F}}(\textbf{{r}}_{ij})\,]+\sum_{j\in\Omega_{3}^{i}}\textbf{{F}}(\textbf{{r}}_{ij}+\textbf{{d}}_{ij})$
$\displaystyle\approx\sum_{j\in\Omega_{2}^{i}}\nabla\textbf{{F}}(\textbf{{r}}_{ij})\cdot\textbf{{d}}_{ij}+\sum_{j\in\Omega_{3}^{i}}[\textbf{{F}}(\textbf{{r}}_{ij})+\nabla\textbf{{F}}(\textbf{{r}}_{ij})\cdot\textbf{{d}}_{ij}].$
(3)

The last approximation holds by Taylor expansions, so it is important to keep
in mind that (3) is only valid when $|\textbf{{d}}_{ij}|$ is small. To
estimate the error of the force, it is crucial to provide a proper definition
of the word “error” first. We adopt the widely used root mean squared (RMS)
force error
$\mathcal{E}(r_{1},r_{2},M)=\sqrt{\langle|\Delta\textbf{{F}}|^{2}\rangle}$,
where the $\langle\cdot\rangle$ is the average over all positions
$\textbf{{r}}_{i},i=1,\cdots,N$ and all displacements
$\textbf{{d}}_{i},i=1,\cdots,N$.

To calculate this error we need some assumptions:

  1. 1.

$\textbf{{r}}_{i},\ i=1,\cdots,N$ are random variables with uniform
distributions over the space.

  2. 2.

$\textbf{{d}}_{i},\ i=1,\cdots,N$ are random variables having a normal
distribution with mean 0 and variance $d^{2}$, namely $\mathcal{N}(0,d^{2})$.

  3. 3.

$\textbf{{r}}_{i}$ and $\textbf{{r}}_{j}$ are independent, if $i\neq j$.

  4. 4.

$\textbf{{d}}_{i}$ and $\textbf{{d}}_{j}$ are independent, if $i\neq j$.

  5. 5.

$\textbf{{r}}_{i}$ and $\textbf{{d}}_{j}$ are independent, for any $i$ and
$j$.

The standard deviation $d$ of the random variable $\textbf{{d}}_{i}$ can be
approximately related to the neighbor list updating frequency $M$ by

$\displaystyle d\approx M\Delta t\sqrt{\frac{3k_{B}T}{m}},$ (4)

where $k_{B}$ is the Boltzmann constant, and $m$ is the mass of the particles.
By the theorem of equipartition, ${3k_{B}T}/m$ is the mean squared velocity.
This relation is only valid when $d$ is smaller than the mean free path. In
all the test simulations studied in this paper, the largest deviation of $d$
from $M\Delta t\sqrt{3k_{B}T/m}$ is less than 10%. All the above assumptions
are ideal cases that facilitate the derivation of the error estimate. However,
in systems studied in real problems, these assumptions can be violated. In
section III, we will show when the real force error deviates from our
theoretical estimate and to what extent the theoretical estimate is reliable.

Based on the aforementioned assumptions, we reach the error estimate:

$\displaystyle\mathcal{E}^{2}(r_{1},r_{2},M)=$
$\displaystyle(4\pi\rho)^{2}d^{2}\bigg{\\{}\int_{r_{1}}^{\infty}\frac{1}{3}r^{2}u^{\prime\prime}(r)\textsf{d}r+\int_{r_{1}}^{\infty}\frac{2}{3}ru^{\prime}(r)\textsf{d}r\bigg{\\}}^{2}+$
$\displaystyle
8\pi\rho\,d^{2}\bigg{\\{}\int_{r_{1}}^{\infty}\frac{1}{3}\,[ru^{\prime\prime}(r)]^{2}\textsf{d}r+\int_{r_{1}}^{\infty}\frac{2}{3}\,[u^{\prime}(r)]^{2}\textsf{d}r\bigg{\\}}+$
$\displaystyle 4\pi\rho\int_{r_{2}}^{\infty}[ru^{\prime}(r)]^{2}\textsf{d}r.$
(5)

It is straightforward to develop the error estimate of the short-range
interaction that has the form $u(r)=4\varepsilon(\sigma/r)^{m}$ and $m>3$:

$\displaystyle\mathcal{E}^{2}(r_{1},r_{2},M)=\;$
$\displaystyle(\frac{16}{3}m\pi\rho\varepsilon\sigma^{m})^{2}d^{2}\bigg{(}\frac{1}{r_{1}^{2m-2}}\bigg{)}+$
$\displaystyle
2\frac{m^{2}+2m+3}{3(2m+1)}\pi\rho(8m\varepsilon\sigma^{m})^{2}d^{2}\bigg{(}\frac{1}{r_{1}^{2m+1}}\bigg{)}+$
$\displaystyle\frac{4}{2m-1}\pi\rho(4m\varepsilon\sigma^{6})^{2}\bigg{(}\frac{1}{r_{2}^{2m-1}}\bigg{)}.$
(6)

###  II.3 Tuning the working parameters for the twin-range cutoff method

In general, the calculation of the cutoffed short-range interaction with the
neighbor list method is performed in two steps: 1, generate the neighbor list
at every $M$ steps. 2, calculate the interactions by using the neighbor list
in the subsequent $M-1$ steps. We use the following formulas to estimate the
computational costs of the force calculation ($T_{F}$) and the neighbor list
generation ($T_{N}$):

$\displaystyle T_{F}$ $\displaystyle=c_{1}r_{1}^{3}+d_{1},$ (7) $\displaystyle
T_{N}$ $\displaystyle=c_{2}r_{2}^{3}+d_{2}.$ (8)

Where $c_{1}$, $d_{1}$, $c_{2}$ and $d_{2}$ are constants depending on the
system being studied, the software implementation, as well as the hardware
architecture. For a good estimate of these constants, we time a series of
short test runs with different cutoff radii and fit the computational costs in
the least square sense. The average computational cost of the short-range
interaction per step is

$\displaystyle T=T_{F}+\frac{1}{M}T_{N}.$ (9)

Provided with both the error and the computational cost estimates for the
twin-range cutoff method, it is possible to design a routine that determines
the most efficient combination of parameters $r_{1}$, $r_{2}$ and $M$. Here,
the phrase “most efficient” refers to a set of parameters that reaches a given
accuracy at minimal computational cost. From a mathematical point of view,
this is a constrained optimization problem that can be written in the
following way:

$\displaystyle\min\quad$ $\displaystyle T(r_{1},r_{2},M),$ (10) s.t.
$\displaystyle\mathcal{E}(r_{1},r_{2},M)=\mathcal{E}^{\ast}\quad\textrm{and}$
(11) $\displaystyle d\,(M)\leq d_{0}.$ (12)

Where $\mathcal{E}^{\ast}$ defines the required accuracy. Constraint (12) is
added because the Taylor expansions in (3) are good approximations only when
$d$ is small. In the multicomponent systems, the displacement should be
constrained for the lightest particles. We find $d_{0}=0.2\sigma$ will give
reasonable results. Since $M$ can be analytically solved by constraint (11),
saying $M=M(r_{1},r_{2},\mathcal{E}^{\ast})$, (10) – (12) is reduced to

$\displaystyle\min\quad$ $\displaystyle
T(r_{1},r_{2},M(r_{1},r_{2},\mathcal{E}^{\ast})),$ (13) s.t. $\displaystyle
d\,(M(r_{1},r_{2},\mathcal{E}^{\ast}))\leq d_{0}$ (14)

The constrained optimization problem (13) and (14) can be solved by standard
optimization algorithms.

##  III Results of Testing simulations

We ran molecular dynamics (MD) simulations on an Intel Xeon E5520 Processor
with Gromacs 4.0.7 compiled by GCC 4.5. The testing systems used in our
studies contained 16,000 particles interacting via the standard Lennard-Jones
6-12 interaction. The MD time step was $\Delta t=0.002\,\tau$, where
$\tau=\sigma\sqrt{m/\varepsilon}$. NVT simulations were performed by coupling
the systems to the Nosè-Hoover thermostat with a relaxation time $1\,\tau$.
Short test runs of 10,000 steps, $r_{1}=r_{2}-\sigma$ and $M=20$ were
performed at different $r_{2}$ to provide an estimate of the constants in the
computational cost expressions (7) and (8). This process lasted for about an
hour. It is worthwhile to spend this time because it is short comparing with
the time costs of real simulations. Moreover, the constants can be reused in
all simulations on the same machine with the same density. The error estimate
(6) ($m=6$ for Lennard-Jones interaction) was studied by the systems in the
gas ($T=1.20\varepsilon/k_{B}$, $\rho=0.05\sigma^{-3}$), liquid
($T=1.20\varepsilon/k_{B}$, $\rho=0.80\sigma^{-3}$) and supercritical state
($T=1.34\varepsilon/k_{B}$, $\rho=0.30\sigma^{-3}$). Target precisions
$10^{-2}\varepsilon/\sigma$, $10^{-3}\varepsilon/\sigma$ and
$10^{-4}\varepsilon/\sigma$ were tested. To measure the real accuracies, the
cutoff radii equal to half the simulation boxes were employed and the
resulting forces served as the exact forces. The tuned parameters, accuracies
and computational costs of the single- and twin-range cutoff methods are
listed in Tables 1 and 2.

$\rho$ [$\sigma^{-3}$] |  $\mathcal{E}^{\ast}$ [$\varepsilon/\sigma$] |  $r_{1}$ [$\sigma$] |  $d$ [$\sigma$] | $M$ |  $\mathcal{E}_{real}$ [$\varepsilon/\sigma$] |  $10^{4}\,T_{est}$ [s] |  $10^{4}\,T_{real}$ [s]
---|---|---|---|---|---|---|---
| $10^{-2}$ | 3.23 | 0.200 | 52 | $1.10\times 10^{-2}$ | 21.3 | 19.7
0.05 | $10^{-3}$ | 4.90 | 0.200 | 52 | $1.09\times 10^{-3}$ | 37.0 | 36.8
| $10^{-4}$ | 7.48 | 0.200 | 52 | $1.12\times 10^{-4}$ | 92.9 | 93.1
| $10^{-2}$ | 3.91 | 0.180 | 44 | $1.17\times 10^{-2}$ | 85.1 | 88.0
0.30 | $10^{-3}$ | 5.87 | 0.126 | 31 | $1.43\times 10^{-3}$ | 266 | 268
| $10^{-4}$ | 8.91 | 0.100 | 25 | $1.86\times 10^{-4}$ | 914 | 851
| $10^{-2}$ | 4.28 | 0.110 | 29 | $0.69\times 10^{-2}$ | 273 | 251
0.80 | $10^{-3}$ | 6.45 | 0.081 | 21 | $0.64\times 10^{-3}$ | 900 | 873
| $10^{-4}$ | 9.81 | 0.067 | 17 | $0.61\times 10^{-4}$ | 3168 | 3120
Table 1:  The tuned parameters of the single-range cutoff method. System in the gas ($\rho=0.05\sigma^{-3}$, $T=1.20\varepsilon/k_{B}$), supercritical ($\rho=0.30\sigma^{-3}$, $T=1.34\varepsilon/k_{B}$) as well as liquid state ($\rho=0.80\sigma^{-3}$, $T=1.20\varepsilon/k_{B}$) were tested. $\mathcal{E}^{\ast}$ is the target accuracy while the $\mathcal{E}_{real}$ is the real error calculated from the simulation. $T_{est}$ is the computational expense estimated by (9). $T_{real}$ is the real computational cost timed in the test simulations. The unit of the computational costs is second per step.  $\rho$ [$\sigma^{-3}$] |  $\mathcal{E}^{\ast}$ [$\varepsilon/\sigma$] |  $r_{1}$ [$\sigma$] |  $r_{2}$ [$\sigma$] |  $d$ [$\sigma$] | $M$ |  $\mathcal{E}_{real}$ [$\varepsilon/\sigma$] |  $10^{4}\,T_{est}$ [s] |  $10^{4}\,T_{real}$ [s]
---|---|---|---|---|---|---|---|---
| $10^{-2}$ | 2.83 | 3.69 | 0.200 | 52 | $1.08\times 10^{-2}$ | 19.7 | 19.3
0.05 | $10^{-3}$ | 4.24 | 5.61 | 0.200 | 52 | $1.03\times 10^{-3}$ | 30.9 | 33.0
| $10^{-4}$ | 6.05 | 8.24 | 0.126 | 33 | $1.02\times 10^{-4}$ | 69.7 | 77.4
| $10^{-2}$ | 3.41 | 4.31 | 0.127 | 31 | $1.06\times 10^{-2}$ | 72.9 | 82.0
0.30 | $10^{-3}$ | 4.96 | 6.36 | 0.082 | 20 | $1.20\times 10^{-3}$ | 215 | 235
| $10^{-4}$ | 7.52 | 9.54 | 0.065 | 16 | $1.49\times 10^{-4}$ | 738 | 738
| $10^{-2}$ | 3.76 | 4.64 | 0.080 | 21 | $0.85\times 10^{-2}$ | 238 | 228
0.80 | $10^{-3}$ | 5.60 | 6.90 | 0.057 | 14 | $0.76\times 10^{-3}$ | 767 | 793
| $10^{-4}$ | 8.59 | 10.41 | 0.047 | 12 | $0.75\times 10^{-4}$ | 2742 | 2882
Table 2:  The tuned parameters of the twin-range cutoff method. System in the
gas ($\rho=0.05\sigma^{-3}$, $T=1.20\varepsilon/k_{B}$), supercritical
($\rho=0.30\sigma^{-3}$, $T=1.34\varepsilon/k_{B}$) as well as liquid state
($\rho=0.80\sigma^{-3}$, $T=1.20\varepsilon/k_{B}$) were tested.
$\mathcal{E}^{\ast}$ is the target accuracy while the $\mathcal{E}_{real}$ is
the real error calculated from the simulation. $T_{est}$ is the computational
expense estimated by (9). $T_{real}$ is the real computational cost timed in
the test simulations. The unit of the computational costs is second per step.
The last column gives the acceleration ratios of the real computational costs
with respect to the real costs of the single-range cutoff method.

In all cases presented, the constraint (11) is strictly satisfied, so the
deviation of the real errors from the target errors measures the quality of
the error estimates. In the gas state, the error estimates are sharp. In the
supercritical and liquid states, the error estimates do not deviate very far
from the real errors, and all of them fall in the range
$[\frac{1}{2}\mathcal{E}_{real},2\mathcal{E}_{real}]$. In the supercritical
cases, the errors tend to be underestimated, while in the liquid systems, the
errors are somehow overestimated. The estimated RMS errors are not exactly the
real values, because some of the assumptions in section II.2 are not preserved
in real simulations. The assumptions 3, 4 and 5 are obviously not satisfied
because any two particles cannot overlap due to the repulsive core of the
Lennard-Jones interaction. Moreover, they are also violated by the
correlations between particles due to the attractive dispersion term. It is
possible to include the pair distribution information (i.e. the radial
distribution function $g(r)$ in the homogeneous systems) to improve the
quality of the error estimates, but the estimates are then _posterior_ rather
than _prior_ , which is not convenient for the parameter tuning.

From Table 1 and 2, it is obvious that higher target accuracies and larger
system densities require more intensive parameters, i.e. larger cutoff radii
and higher neighbor list updating frequencies. In most gas cases, the
displacement $d$ hits the constraint (12). In all cases shown, the $M$ are
larger than the usually used value ($M=10$, default in Gromacs). In the gas
state, the $M$ can be as large as 52. The twin-range cutoff method is always
more efficient than the single-range cutoff method. However, in the liquid
state, the twin-range cutoff method is a little bit less accurate than the
single-range cutoff method. So the benefit of using twin-range cutoff should
be discounted in the sense of the same accuracy.

To test the benefit we gain by the parameter tuning, we used a set of
unoptimized single-range parameters ($r=4.00\sigma$ and $M=10$) in the gas
state ($\rho=0.05\sigma^{-3}$, $T=1.20\varepsilon/k_{B}$). The accuracy and
efficiency of these parameters are $\mathcal{E}=3.01\times
10^{-3}\varepsilon/\sigma$ and $T=4.10\times 10^{-3}\textsf{s}$, respectively.
To compare the efficiencies, we manually adjusted the target precisions
$\mathcal{E}^{\ast}$ so that the real accuracy of the tuned parameters is
nearly the same as the corresponding unoptimized parameters. The tuned single-
range parameters are $r=4.09\sigma$ and $M=52$ with computational cost
$T=2.71\times 10^{-3}\textsf{s}$. And the tuned twin-range parameters are
$r_{1}=3.54\sigma$, $r_{2}=4.68\sigma$ and $M=52$ with $T=2.56\times
10^{-3}\textsf{s}$. The tuned single- and twin-range parameters save 33% and
38% computational costs, respectively. The same unoptimized parameters result
in $\mathcal{E}=6.10\times 10^{-3}\varepsilon/\sigma$ and $T=2.74\times
10^{-2}\textsf{s}$ in the liquid state ($\rho=0.80\sigma^{-3}$,
$T=1.20\varepsilon/k_{B}$). At the same accuracy level, the tuned single-range
parameters are $r=4.36\sigma$ and $M=28$ with the computational cost
$T=2.76\times 10^{-2}\textsf{s}$. And the tuned twin-range parameters are
$r_{1}=3.94\sigma$, $r_{2}=4.86\sigma$ and $M=19$ with $T=2.76\times
10^{-2}\textsf{s}$. In this case, the tuned parameters do not improve the
efficiency at all, because the unoptimized parameters are good enough.
Moreover, since the error estimate (6) is not sharp in the liquid state and
the constants $c_{1}$, $d_{1}$, $c_{2}$ and $d_{2}$ in (7) and (8) are not
exactly measured, the tuned parameters are only nearly optimal rather than
really optimal. So the computational costs of the tuned parameters can be
trivially more expensive than the unoptimized ones. We also tested a set of
unoptimized twin-range parameters that are taken from one of our former
simulation settings: $r_{1}=2.5\sigma$, $r_{2}=4.0\sigma$ and $M=10$. In the
gas state, the error is $\mathcal{E}=4.48\times 10^{-3}\varepsilon/\sigma$
with $T=3.87\times 10^{-3}\textsf{s}$. The tuned parameters are
$r_{1}=3.28\sigma$, $r_{2}=4.31\sigma$ and $M=52$ with the computational cost
of $T=2.29\times 10^{-3}\textsf{s}$, 40% cheaper. While in the liquid state,
the accuracy and computational cost of the unoptimized parameters are
$\mathcal{E}=2.91\times 10^{-2}\varepsilon/\sigma$ and $T=1.85\times
10^{-2}\textsf{s}$, respectively. The tuned parameters are $r_{1}=3.12\sigma$,
$r_{2}=3.81\sigma$ and $M=27$ with the computational cost of $T=1.37\times
10^{-2}\textsf{s}$, 26% faster than the unoptimized parameters.

##  IV Conclusions and open questions

In this paper, an error estimate of the twin-range cutoff method was
developed. The error estimate of the single-range cutoff method was easily
derived from the twin-range cutoff method, because the former is only a
special case of the later. Equipped with both the error and the computational
cost estimates, we proposed a work flow that can automatically determine the
nearly optimal parameters demanding a certain target accuracy. We verified the
effectiveness of the error estimate and parameter tuning algorithm by
numerical simulations of Lennard-Jones 6-12 system in gas, liquid and
supercritical state. We also presented some examples to show, by parameter
tuning, to what extent the computational expense can be saved with respect to
the unoptimized parameters.

In the gas state, the error estimates are sharp. In the liquid and
supercritical state, the error estimates are reasonable: they are always
larger than half and smaller than twice of the real errors. The quality of the
error estimates depends on the extent to which the assumptions in section II.2
are satisfied. These estimated force errors provide a quantitative description
of the undesirable cutoff effects introduced by the single- and twin-range
cutoff methods. The parameter tuning algorithm enables an automatic searching
of the working parameters (i.e. the cutoff radii and the neighbor list
updating frequency) that reach a predetermined accuracy at nearly minimal
computational cost. Comparing with the unoptimized parameters, the tuned
parameters are always faster. The benefit ranged from 0% to 40%, depending on
the quality of the unoptimized parameters. We also demonstrated that the
optimized twin-range cutoff method is faster than the optimized single-range
cutoff method, so the former is recommended. Combining with the parameter
tuning algorithm developed for the long-range interactions wang2010optimizing
, all non-bonded interactions can be calculated at a nearly optimal
efficiency.

The error estimate and the parameter tuning for inhomogeneous systems were not
considered by the present paper, because the assumptions for the error
estimate are worse violated. However, how to handle the inhomogeneous systems
is a very important open question in this field, which requires further
studies.

## Acknowledgments

The authors gratefully acknowledge Florian Dommert for our valuable
discussions and Jia Deng for his help editing our English. P.Z. thanks the
financial support by the National Natural Science Foundation of China
(50930003).

## References

  * (1) Daan Frenkel and Berend Smit.  Understanding Molecular Simulation.  Academic Press, San Diego, second edition, 2002.
  * (2) B. Smit.  Phase diagrams of lennard-jones fluids.  The Journal of Chemical Physics, 96(11):8639, 1992.
  * (3) VG Baidakov, GG Chernykh, and SP Protsenko.  Effect of the cut-off radius of the intermolecular potential on phase equilibrium and surface tension in lennard-jones systems.  Chemical Physics Letters, 321(3-4):315–320, 2000.
  * (4) W.Z. Ou-Yang, Z.Y. Lu, T.F. Shi, Z.Y. Sun, and L.J. An.  A molecular-dynamics simulation study on the dependence of lennard-jones gas-liquid phase diagram on the long-range part of the interactions.  The Journal of chemical physics, 123:234502, 2005.
  * (5) AZ Panagiotopoulos.  Molecular simulation of phase coexistence: Finite-size effects and determination of critical parameters for two-and three-dimensional lennard-jones fluids.  International Journal of Thermophysics, 15(6):1057–1072, 1994.
  * (6) M. Mecke, J. Winkelmann, and J. Fischer.  Molecular dynamics simulation of the liquid–vapor interface: The lennard-jones fluid.  The Journal of Chemical Physics, 107:9264, 1997.
  * (7) A. Trokhymchuk and J. Alejandre.  Computer simulations of liquid/vapor interface in lennard-jones fluids: Some questions and answers.  The Journal of Chemical Physics, 111:8510, 1999.
  * (8) V.K. Shen, R.D. Mountain, and J.R. Errington.  Comparative study of the effect of tail corrections on surface tension determined by molecular simulation.  J. Phys. Chem. B, 111(22):6198–6207, 2007.
  * (9) A. Lotfi, J. Vrabec, and J. Fischer.  Vapour liquid equilibria of the lennard-jones fluid from the npt plus test particle method.  Molecular Physics, 76(6):1319–1333, 1992.
  * (10) JM Caillol.  Critical-point of the lennard-jones fluid: A finite-size scaling study.  The Journal of Chemical Physics, 109:4885, 1998.
  * (11) J.J. Potoff and A.Z. Panagiotopoulos.  Critical point and phase behavior of the pure fluid and a lennard-jones mixture.  The Journal of Chemical Physics, 109:10914, 1998.
  * (12) J. Pérez-Pellitero, P. Ungerer, G. Orkoulas, and A.D. Mackie.  Critical point estimation of the lennard-jones pure fluid and binary mixtures.  The Journal of chemical physics, 125:054515, 2006.
  * (13) Mike P. Allen and Dominik J. Tildesley.  Computer Simulation of Liquids.  Oxford Science Publications. Clarendon Press, Oxford, 1 edition, 1987\.
  * (14) P. Lagüe, R.W. Pastor, and B.R. Brooks.  Pressure-based long-range correction for Lennard-Jones interactions in molecular dynamics simulations: application to alkanes and interfaces.  J. Phys. Chem. B, 108(1):363–368, 2004.
  * (15) J. Janecek.  Long range corrections in inhomogeneous simulations.  J. Phys. Chem. B, 110(12):6264–6269, 2006.
  * (16) M.R. Shirts, D.L. Mobley, J.D. Chodera, and V.S. Pande.  Accurate and efficient corrections for missing dispersion interactions in molecular simulations.  J. Phys. Chem. B, 111(45):13052–13063, 2007.
  * (17) J. Kolafa and J.W. Perram.  Cutoff Errors in the Ewald Summation Formulae for Point Charge Systems.  Molecular Simulation, 9(5):351–368, 1992.
  * (18) Gerhard Hummer.  The numerical accuracy of truncated Ewald sums for periodic systems with long-range coulomb interactions.  Chem. Phys. Lett., 235:297–302, 1995.
  * (19) H.G. Petersen.  Accuracy and efficiency of the particle mesh Ewald method.  The Journal of Chemical Physics, 103:3668, 1995.
  * (20) M. Deserno and C. Holm.  How to mesh up Ewald sums. II. An accurate error estimate for the particle–particle–particle-mesh algorithm.  The Journal of Chemical Physics, 109:7694, 1998.
  * (21) H. Wang, F. Dommert, and C. Holm.  Optimizing working parameters of the smooth particle mesh Ewald algorithm in terms of accuracy and efficiency.  The Journal of chemical physics, 133:034117, 2010.