On the Diffusion of Sticky Particles in 1D
Joshua D M Hellier Graeme J Ackland





 J.D.M.Hellier@sms.ed.ac.uk

 G.J.Ackland@ed.ac.uk

 SUPA, School of Physics and Astronomy, University of Edinburgh, Mayfield Road, Edinburgh EH9 3JZ, United Kingdom



             


The 1D Ising model is the simplest Hamiltonian-based model in statistical mechanics. The sim-
plest interacting particle process is the Symmetric Exclusion Process (SEP), a 1D lattice gas of
particles that hop symmetrically and cannot overlap. Combining the two gives a model for sticky
particle diffusion, SPM, which is described here. SPM dynamics are based on SEP with short-range
interaction, allowing flow due to non-equilibrium boundary conditions. We prove that SPM is also
a detailed-balance respecting, particle-conserving, Monte Carlo (MC) description of the Ising model. Neither the Ising model nor SEP have a phase transition in 1D, but the SPM exhibits a non-equilibrium
transition from a diffusing to a blocked state as stickiness increases.  This transition manifests in peaks in the MC density fluctuation, a change in the dependency of flow-rate on stickiness, and odd structure in the eigenspectrum of the transition rate matrix.
We derive and solve a fully non-linear, analytic, mean-field solution, which has a crossover from a positive to a negative diffusion constant where the MC shows the transition.
The negative diffusion constant in fact indicates a breakdown of the mean-field approximation, with close to zero flow and breaking into a two-phase mixture, and thus the mean field theory successfully predicts its own demise.
We also present an analytic solution for the flow via direct analysis of the transition rate matrix.  The
simplicity of the model suggest a wide range of possible applications.





                              
On The Diffusion of Sticky Particles in 1-D
    Graeme J Ackland
    December 30, 2023
===========================================



 



§ INTRODUCTION


Lattice gases are a ubiquitous tool for modeling complex systems from
biology to traffic <cit.>.  Analytically solvable cases involve
non-interacting or excluding particles <cit.>, but in any real system of interest the moving
objects interact. Many models tackle the situation where the diffusing
objects interact with the substrate, but despite clear
application-relevance there is surprisingly little work considering
interactions between the moving particles themselves.  One reason for
this is that the interactions introduce nonlinearities in analytical
models, which makes them challenging to solve, at least outside of
limits in which they can be linearized. This is unfortunate because it
is precisely these nonlinearities which introduce interesting
behaviors such as discontinuities at interfaces or
diffusion instability <cit.>.

Another feature of previous models is that the flow is driven by
either asymmetry in the dynamics (e.g. the Asymmetric Simple Exclusion Process, ASEP) or an external field
which permeates the system (e.g. the Katz-Lebowitz-Spohn model, KLS).  In either case, the particles
always see a local asymmetry.  However in many systems the flow is
driven by a pressure or chemical potential difference applied at the
boundaries, so that any local asymmetry arises from self-organization.
This is the situation that is addressed here.

In this paper we will investigate what we believe to be the simplest interacting one-dimensional exclusion model, the
“Sticky Particle Model” or SPM, which is specified in Fig. <ref>. We
will explore the impact this has on particle behavior, in particular
when observed in the large-scale limit. The model will be fully defined in Section <ref>, and we will show that it obeys detailed balance away from boundaries. We will then discuss a method to solve the system exactly for small systems in Section <ref>, and Monte-Carlo techniques for larger systems in Section <ref>. In Section <ref> we will attempt to understand our results by means of a simple MFT approximation.  One might contrast this
approach (making a simple microscopic model and trying to learn about large-scale properties) with approaches such as the KPZ
equation <cit.>
(where one analyses the extreme large-scale dynamics using
universality classes).
Throughout the simulations and analytics, we aim to address three kinds of questions: how does stickiness affect the flow, how does the
particle current depend on the driving force, and what determines the density of particles in the system?












§ THE STICKY PARTICLE MODEL
 

The SPM is an excluded-particle model in which adjacent particles
separate with rate λ and single particles move at rate 1.


It differs from the symmetric exclusion
process <cit.> in that particles “unstick”
with rate λ instead of their normal hopping rate, 1.  Low
λ corresponds to sticky particles, high λ to repelling
particles.
It could be regarded as a version of the KLS model <cit.> in 1-dimension without an applied field, which
is itself similar to the dynamics used to analyze the Ising model by
Kawasaki <cit.>; it is a member of the class of models considered in
<cit.>.  The KLS model has a field that
introduces an asymmetry which in turn drives a flow.  It is tacitly assumed
that this field is required for the nontrivial flow behavior seen in
the model, so the simpler symmetric model has received less attention.




 §.§ Detailed Balance Proof
 
The SPM is intended to study flow, so it is defined by the
hopping rates.  In Fig. <ref> we
show all the possible transitions which may occur between local
configurations. 

Assume that the system is now on a ring, with L
lattice sites and N particles.  Let us label possible system
configurations by ξ and let the number of adjacencies (or
“bonds”) between particles be b(ξ). Now for our ansatz, assume
that the probability of the system being in state ξ is
λ^-b(ξ).  In the top and bottom diagrams of
Fig. <ref> we can see that the number of bonds on
both sides is the same, as are the transition rates back and forth;
thus our ansatz holds for these states, as it predicts the
probabilities of the left and right configurations are the same. The
middle two diagrams are equivalent; in the upper diagram a
bond is formed going left to right and then broken going right to
left, so the probability of being in the left state is λ times
that of being in the right state. This is again in agreement with the
detailed balance criterion. As these are the only types of transition
that may occur on a ring, we have proven that the closed system obeys
detailed balance, with an energy proportional to b(ξ) lnλ.



 §.§ Equivalence to Ising Model in Equilibrium


Since the model obeys detailed balance, there must be an associated
Hamiltonian.  This is simply proportional to the number of particles stuck together:

    H = q ∑_i p_ip_i+1 = q b(ξ)

where p_i is 1 for a particle or zero for an empty site.  A simple
transformation s_i = 2p_i-1 shows this to be equivalent to the Ising
Hamiltonian

    H = q ∑_i [ 1/4 s_i s_i+1 + ( 1/2 s_i + 1/4) ].

The contribution from the final term is constant, and because only hopping moves are allowed
∑ s_i is constant<cit.>. 
Therefore in equilibrium the SPM samples the 1D
Ising model with fixed magnetization, a fact which was used to
validate the codes. Compatibility between the detailed balance condition, the standard Boltzmann distribution at temperature T and our defined rates requires that q = k_B Tlnλ.



 §.§ Nonequilibrium Calculations


Of course, the real point of our research into the SPM is to see how it behaves out of equilibrium. We drive a system of length L out of equilibrium by attaching it to
reservoirs with densities ρ_0 and ρ_L at the x=0 and x=L ends respectively,
with a single fixed value for λ throughout the system.
We have computed the main properties of this system (principally the current and the average density of particles) by using three methods: numerical analysis of the transition rate matrix (TRM) for a small system, Monte Carlo (MC) methods for somewhat larger systems and a mean-field
approximation for large systems in the continuum limit. These methods and their results are discussed in Sections <ref>, <ref> and <ref>. respectively.



§ TRANSITION RATE MATRIX CALCULATIONS


The SPM, just like the KLS model and ASEP, can exhibit flow. In this paper we will be principally concerned with the computation of the properties of the steady state of an SPM system of L sites
connected to (unequal) particle reservoirs at either end.

It is possible to “analytically” solve the SPM on a finite domain with fixed boundary densities by analyzing the transition rate matrix which represents that system. For a system of length L+4 this is a sparse matrix with dimension 
2^L+4× 2^L+4.
The non-zero matrix elements relating to bulk motion are λ or 1 while those relating to the boundary conditions are given by Eqs. <ref> and <ref>. 
The matrix represents a finite irreducible continuous time Markov process, so it is guaranteed to have a unique zero eigenvalue,
whose eigenvector represents the long-term steady state flow. All other eigenvalues have negative real part and represent fluctuations.
It is actually possible
to “analytically” solve the SPM on a finite domain with fixed boundary densities by analyzing the transition rate matrix which represents that system.
To do this, we wrote a Python code, in  which specific configurations of a system of length L+4 are represented by strings of binary numbers.
For example, the configuration “” corresponds to a system
with L=4 where there is a particle in sites with a “” and a vacancy in sites with a “”. As there is a clear correspondence between binary strings and positive integers, in our code we loop over
the configurations; in each configuration, we can loop over the particles, see which moves they can make and with what rates, and which configurations these moves would take us to. Doing this allows us to add
entries to a sparse matrix representation of the SPM with boundary conditions.
There are a pair of additional sites at each boundary
representing the particle reservoirs, hence the overall system size being recorded as L + 4. In these boundary layers particles rapidly pop in and out of existence,
with rates such that the mean occupation should be the desired boundary density. The choice of how to achieve this is not unique, but in our computations the rate at which particles appear in empty spaces
at the boundary is 

    C_0 (1+λ)√(λρ_0/1-ρ_0),

and the rate at which particles disappear is

    C_0 (1+λ)√(λ1-ρ_0/ρ_0).


In isolation, the above rates would ensure that the average occupation was always ρ_0. We have an equivalent set of creation and annihilation rates at the other boundary attempting to maintain the density at ρ_L. In our overall SPM system, we attempt to ensure that all of these rates are larger than, or at the very least not smaller than,
the other rates in the system (1 and λ); this should ensure that the occupation of the boundary sites is always approximately ρ_0, regardless of what is happening internally.
Thus, C_0 is chosen to be quite large, and the λ-dependence exists primarily in order to ensure that λ does not dominate the boundary rates when it becomes large.
It is worth noting that with the addition of these boundary conditions, the model as a whole no longer obeys detailed balance, even though the bulk still does in terms of
local behavior.

Once we have built the correct transition rate matrix, we then use the
Python routine  to find the
eigenvector and associated eigenvalue with the most positive real
part; this eigenvalue
should be zero, up to numerical error, and the corresponding
eigenvector describes the system's
stationary state.  The stationary state contains most of the key
information we desire about the system's behavior during steady flow,
including the current and particle density profile, which may be
accessed by the application of occupation and current operators
(constructed along with the transition rate matrix). In this steady
state, we are guaranteed to have a uniform current flowing through the
interior of the system; this is implied by the uniqueness of the
steady state and the homogeneity of the internal dynamics.  The next
least-negative eigenpairs dominate the approach to equilibrium, and
determine the timescale over which it is achieved.



The near-zero eigenspectrum of the transition rate operator, and its
dependence upon λ for fixed boundary conditions, is extremely
interesting in its own right, as shown in Fig. <ref>.
The structure is rich, with eigenvalues frequently crossing over each
other, as well the appearance and disappearance of bands.

Repeat calculations show that large changes in C_0 cause very small
changes in the low-lying eigenspectrum of the transition rate matrix,
so in this sense the slow dynamics are insensitive to C_0. It does
affect the group of eigenvalues of larger magnitude,
which form a u-shape towards the top of Fig. <ref>,
suggesting that this group is associated with rapid changes in
configuration at the boundaries. The
extremely large gap between the low-lying and upper sequences of
eigenvalues reinforces that notion.  


One of the most notable features of the eigenspectrum is the presence
of a relatively small number of low-lying eigenvalues which seem to
split off from the main sequence as λ becomes much smaller
than 1; they have 𝒪(λ) scaling, whilst the main
sequence has 𝒪(1) scaling. This split coincides with, and
may be related to, a suspected transition, which we will discuss in
more detail in Section <ref>.  System size variations
suggest that the number of eigenvalues which belong to this grouping
scale as 𝒪(2^∼ 0.8L), slower than the total number of
eigenvalues, 2^L+4, suggesting that they describe long-range
phenomena.  The split begins as the lowest non-zero eigenvalue takes
its lowest values: this may correspond to a change in the character of
the flow (i.e. eigenvector of the zero-eigenvector) to include more
microstates corresponding to high particle density.




Unfortunately the space complexity of the sparse transition rate
matrix scales as 𝒪(2^L L), whilst the matrix density
scales as 𝒪(2^-L L).  This makes computations rapidly
become unfeasible for large L. The overall time complexity for this
process is 𝒪(2^2LL).  In spite of this, we have managed
to compute properties such as the steady-state flow rate for systems
with sizes as large as 10 internal lattice sites, which are large
enough to produce consistent outputs comparable to our other methods
(as shown in Figures <ref> and <ref>),
whilst being computationally cheap enough to allow us to perform
calculations with many different parameter configurations.




§ SIMULATIONS


Given access to large-scale computational hardware, it is relatively straightforward to perform Monte Carlo simulations of the SPM, with a little adjustment to the precise implementation of the model. 
We chose to calculate primarily by using the  <cit.>
package, which implements the Kinetic Monte Carlo algorithm
(essentially the same as the Gillespie algorithm <cit.>) on lattice systems. The codes used are kept
here <cit.>.

In the bulk, the transition rates are simply
those described in Fig. <ref>. At the boundaries we used a similar two-layer boundary trick to that described in Sec.<ref>, only this time the particles appear and disappear with rates

    √(λρ_0/1-ρ_0),

and

    √(λ1-ρ_0/ρ_0),

and likewise for ρ_L, respectively. Whilst it would be nice for these rates to be large in order to force the boundaries to maintain the correct density more accurately,
we choose to simply keep them proportional to the geometric mean of 1 and λ,
otherwise they tend to either happen too rarely or far too frequently, causing the KMC algorithm to be inefficient.


Independent of the KMCLib code, we wrote a simple Metropolis-Hastings
algorithm <cit.> which randomly selects single particle hops; this is more efficient than using KMC for values of
λ which are relatively close to 1 as the algorithm is simple and the acceptance rate is relatively high, so we can generate better statistics using this method in that regime.
We calculate the flow from the number of particles entering and
leaving the system at the boundaries.  Since the model is defined in
terms of rates, the flow in “particles per unit time” is a well
defined quantity.



Using KMClib we studied systems of length 64 (lengths 32, 128 and 256 give similar
results as shown in Fig. <ref>.), running them for 400000 Gillespie steps for equilibration
followed by 10000 measurement runs of 1000 steps interspersed with
relaxation runs of 16000 steps. This way we could gather statistics
about flow rates and densities in a well-equilibrated
system. Specifically, we generate a pool of 10000 samples of flow
rate and density, from which we can calculate estimates of the
descriptive statistics of both quantities.
These calculations, in which we held the boundary densities constant whilst varying λ, were repeated using Metropolis-Hastings for a length 100 system, and we have transition rate matrix currents
for a system of length 10. These are plotted on both logarithmic and linear axes in Fig. <ref>.
respectively.





 §.§ Diffusion Coefficient

We compare an MFT prediction (see <ref>) and the KMC numerical results for the
diffusion constant in Fig. <ref>. We see that MFT and
simulation agree well for low stickiness, and both show the symmetry
about ρ_M = 2/3. For high stickiness, where the MFT
prediction gives a negative diffusion constant, we actually see low
positive values for the current.  It should be noted that the MFT
assumes that ⟨ρ⟩ = ρ_M throughout, whereas in the simulations ⟨ρ⟩ tends to be much higher.



 §.§ Particle Motion During Flow

It is instructive to get an overview of how the particles move during
flow. Fig. <ref> show a plot of the flow structure in
the slow-flow regime.  In very short time averages the “striped”
pattern indicates separation into dense and sparse regions, with an
overall concentration gradient arising from the relative size of such
regions.  Particle/vacancy diffusion through the empty/full regions
can be seen.  As the averages are taken over longer times the blocks themselves appear to diffuse.  Averaged over the entire simulation (not shown) the averaging simply gives a smooth density gradient.



Interesting structure is visible. The dynamics look like a random walk with
some tendency for particles to clump; over longer timescales the
diffusive behavior is more evident, with a textured structure
suggesting characteristic velocity of particles or vacancies through
emergent correlated clumps.  Additional plots can be found in the
supplementary materials.






 §.§ Transition in flow character


Fig. <ref> shows how the current varies with stickiness
at fixed driving for a very broad range of λ.  At high λ the rate is simply proportional
to λ for any forcing.  This result is far from trivial - it
means that the overall flow is determined by the faster rate
(λ), not the slower rate (1).  That the MFT averages over the
two is unsurprising (Eq. <ref>), but for the
simulation to avoid having a “rate limiting step” requires the
system to fill to sufficient density that there are always particles
in contact to repel one another, which suggests that the system is performing some kind of self-organization.

When λ = λ_c∼0.2 the simulations show a transition to a different
behavior, with higher or lower λ_c depending on the boundary conditions and the strength of the driving force.
The current and its fluctuation remains finite, but there a distinctive peak in the particle density fluctuation
⟨ (ρ-ρ̅)^2 ⟩ (Fig. <ref>).  The width of the peak is independent of the
system size 
which suggests a continuous
phase transition from the free flowing to the “stuck” regime.



Figs. <ref> and <ref> reveal that in the (0.3, 0.1) and (0.75, 0.25) boundary
configurations the current switches from being 𝒪(λ) at high λ to being
𝒪(λ^k), k∼4 for low λ. In the high-density case, with boundary densities (0.9, 0.7),
the current actually starts flowing backwards for low λ in the transition rate matrix calculation. Therefore, we propose that the
peak in the density fluctuation and the odd behavior in the current suggest the existence of a nonequilibrium
phase transition in the SPM.





 §.§ Effective diffusion constant



Using boundary conditions

    (ρ_0,ρ_L)=(ρ_M+1/2δρ, ρ_M-1/2δρ)

we demonstrate the
dependence of current on boundary density difference and λ (Fig. <ref>).  This
shows that the transition to the stuck phase is suppressed by stronger
driving forces (large |δρ|).

We use the limit of small δρ to calculate the effective
diffusion constant (Fig. <ref>).  This is normalized to 1 in the case of
λ=1, which is just SEP.  The ρ_M=0 limit corresponds to
free flow of particles, so the diffusion constant here does to
1. Similarly the ρ_M=1 limit is flow of vacancies, so
D→λ.  One might expect a monotonic variation between
these limits, but surprisingly the simulations show that there is
always an extremal value for D close to ρ_M=2/3: this is a
minimum for λ<1 and a maximum for λ>1.





 §.§ Self-organized density



We calculated the time-averaged total number of particles in
the system by updating a histogram of particle numbers
as the simulation progresses. In each of our calculations, we make the
initial configuration by randomly filling the system with particles
and vacancies in such a way that the initial density should be
1/2(ρ_0 + ρ_L), and then run the system for a
sufficient number of equilibration steps to destroy any initial
transients.

In SEP, (λ=1) the density varies linearly across the system
from ρ_0 to ρ_L, as one expects for a diffusion process.
However, for sticky particles, λ<1, the density rises sharply
near to the boundary, then has a linear profile anchored about a value higher
than ρ_M  (Fig. <ref>). 
One might view this as particles being sucked into the
system to lower their Hamiltonian Energy, but such a notion can be
dismissed since the internal density also rises for repelling
  particles, λ>1.  The asymptotic values for the bulk density at
high and low λ appear to be 1 and 2/3 respectively, regardless of the 
mean boundary density ρ_M (see Fig. <ref>). 







 §.§ Review of the simulations


MC simulation and TRM analysis of the sticky particle model reveals a number of curious
and unexpected features which have no equivalent in either SEP or the
1D Ising model.

 
  * A transition in the dependence of the current on λ from ∝λ to ∝λ^4.  

  * The mean density of the system  has complicated variation with  λ.  

  * The effective diffusion constant has an extremal value at
  intermediate boundary density.


To obtain further understanding of these behaviors, we tackle the system
analytically, using mean field theory. Although MFT isn't a particularly
sophisticated technique, it does give us a benchmark which we can compare to our other methods.



§ MEAN FIELD THEORY FOR FLOW
 
  We will be working in the style described in Sec. 2.2 of <cit.> and Sec. 5.2.2 of <cit.>.
Let the spacing between lattice
sites be a, let τ_0 be the free-particle hopping timescale, and
the time-averaged (or ensemble-averaged, assuming ergodicity)
occupation probability of the i^th lattice site be
ρ_i.  We will introduce ζ = 1 - λ here for convenience.



 §.§ Lattice Mean Field Theory
 Let the ensemble-averaged occupation
probability of the i^th site at time t be ρ_i
(t). In the mean-field approximation this is assumed to be
independent of ρ_j(t) for j ≠ i at equal times. Therefore,
if the i^th site is occupied, then the rate at which it empties
is the sum of contributions from the
four permutations of the (i+1)^th and
(i-1)^th occupancy: 


    1/τ_0  (1-ρ_i-1)[ (1 - ρ_i+1) + λρ_i+1] 
    
     +   1/τ_0  (1-ρ_i+1)[ (1 - ρ_i-1) + λρ_i-1] .


Similarly, if the (i)^th site is unoccupied, it fills with
rate which depends on the occupation probability of the neighbours,
and whether they are stuck to the (i+2)^th and (i-2)^th sites
respectively. Therefore, an unoccupied i^th site fills with rate

    1/τ_0 {ρ_i+1[ λρ_i+2 + (1-ρ_i+2) ] + ρ_i-1[ λρ_i-2 + (1-ρ_i-2) ] }.


If we now multiply the filling/emptying rates of site i by the
probability of being empty/full respectively, we obtain the
final equation for the site occupation 

    τ_0 ρ_it   = ( 1-ρ_i ) [ (1-ζρ_i-2) ρ_i-1 + (1-ζρ_i+2) ρ_i+1] 
       - ρ_i [ 2 ζρ_i-1ρ_i+1  - (3-ζ)(ρ_i-1 + ρ_i+1) + 2 ].




 §.§ Diffusion Equation: MFT Continuum Limit

 To obtain the continuum limit of the MFT we substitute ρ_i(t)
 →ρ(x, t), ρ_i+m(t) →ρ(x + am, t)
 into Eq. <ref>.  Then, using a Taylor expansion around x for
 small a, neglecting terms of 𝒪(a^4), and collecting
 terms we find that

    τ_0 ρt =    a^2 [ 1-ζρ (4-3ρ)  ] ^2 ρx^2
    
      +    2 a^2 ζ (3ρ-2) (ρx)^2 + 𝒪(a^4) ,

which can be factorized into the more familiar form of a continuity equation

    ρt = a^2/τ_0x{[1 - ζρ(4-3ρ) ] ρx},

having dropped the higher-order terms.  From this we can identify the
flow J(x) as the term in the curly brackets, and define a mean field
diffusion constant (see Fig.<ref>):


    D(ζ,ρ) =  a^2/τ_0[ 1 - ζρ(4-3ρ)  ]


We note that the MFT diffusion constant can become negative. 
The density at
which this occurs is given by finding the roots ρ_c of the RHS of Eq. <ref>:

    ρ_c=2/3±1/3√(()4-3/ζ).

The possible values for ρ_c are only real (and therefore physically possible) for ζ<0.75.  So in the limit
of very sticky particles (λ<0.25) the MFT predicts a
transition between forward and backwards diffusion at some densities.  




 §.§ Limiting Cases


In order to understand the implications of the MFT, let us consider
some limits. As ζ→ 0 (i.e. as the model becomes a
simple exclusion model), D →a^2/τ_0. Likewise,
in the dilute limit ρ→ 0, D →
  a^2/τ_0, reflecting the fact that it becomes a dilute lattice
gas and therefore the interactions between particles become irrelevant
as they never meet.  Conversely, in the full limit ρ→
1, D →λ a^2/τ_0; this is because we now
have a dilute gas of vacancies, which hop with rate
λ/τ_0.  

One may observe that the 
MFT has a symmetry under ρ↦4/3 - ρ; thus, the
dynamics should be symmetric under a density profile reflection around
ρ = 2/3. This is where D always attains its extremal
value, a^2/τ_0[1 - 4/3ζ], hence
for ζ>3/4 the diffusion coefficient becomes negative in regions
with 2/3 - √(ζ(4ζ - 3))/3ζ
< ρ < 2/3 + √(ζ(4ζ -
    3))/3ζ. 

Finally, it is possible to show that
solutions to the continuum MFT containing domains with a negative
diffusion coefficient are linearly unstable; thus, if we try to have a
flow containing ρ for which D(ρ)<0, the density of the
medium should gravitate towards one of the two densities for which D(ρ)∼
0. Instead of observing “backwards diffusion” we would see an
extremely slow flow or no flow at all. 



 §.§ Steady State Flow


It is possible to solve the continuum MFT in a steady state on a
finite domain, say x∈(0, L). Steady state implies that there is no
build-up of particles ρt=0, so the flow is constant through the system.

Using the constant-flow requirement for continuity
J(x)=J_0=const., and by integrating both sides of
Eq. <ref> with respect to x we find that

    ∫_x_0^x J(x')dx'    =    -a^2/τ_0∫dρ[1 - ζρ(4-3ρ) ]
     (x-x_0)J_0     =     -a^2/τ_0[ ρ + ζ (ρ -
     2) ρ^2  ].

The density profile across the system is given by solving for
ρ(x).  Since Eq.<ref> is cubic, the solution for density
ρ(x) is non-unique for cases of high ζ.  Thus in the limit of
high stickiness, the MFT is unable to make a unique prediction for the
density.  Furthermore, except in the SEP case ζ=0, the density
will not vary linearly across the system.



 §.§ Dirichlet Boundary Conditions


The constants x_0 and J_0 relate to the boundary conditions for the
flow, but they need not correspond to a physically-realizable
situation.  For driven systems it is more convenient to consider fixing
the density at each end. 

If we impose such Dirichlet boundary conditions on this system, say
ρ(0)=ρ_0 and ρ(L)=ρ_L, we find that

    J_0 = a^2/L τ_0[ ρ_0 - ρ_L + ζ( ρ_0[ρ_0^2-2] - ρ_L[ρ_L^2-2] ) ].


This equation can be used for direct comparison with the simulations
(Fig. <ref> and Fig. <ref>).  In general, the agreement is
good, except for the region where the MFT predicts negative flow.




 §.§ Interpretation of MFT


The mean field theory enables us to predict flow behavior of the SPM.
It recovers well-known limiting cases such as SEP (λ=0),
however at high stickiness (λ<0.25) it predicts its own
demise, with unphysical negative diffusion constants and by having
multiple solutions for the density at some positions in the system,
meaning that the MFT does not give a unique density
profile.  This breakdown of the MFT corresponds to the transition to
slow flow observed in the simulation.


In some
conditions MFT predicts densities greater than 1.  One might guess that when the
MFT offers two possible values for the density, it will correspond to
a phase separation in the actual system. Furthermore,
the MFT prediction of a maximum diffusion constant at a density of
ρ=2/3 suggests that this value of ρ might be favored for
strongly driven flows.

A curious feature of the MFT diffusion (Fig. <ref>) is
that for fixed stickiness there are two possible densities
giving the same diffusion constant.  Thus it is possible to have a
steady state flow with phase separation into regions of high and low
density. This echoes the situation seen in short time-averages of the
simulation (Fig.<ref>) where blocks of high and low
density are evident, as well as the oscillatory behavior near the boundaries at high λ (Fig. <ref>). 

Regarding the phase transition, we should note that we only have an
MFT prediction for the flow rate for λ>1/4, since
ρ(x) stops being unique when λ drops below 1/4. For higher λ, the MFT is broadly in good agreement with the
simulations, and this continues as λ stretches into the
thousands, where the mean flow varies as 𝒪(λ^1).


The MFT prediction for the mean flow continues to be a good fit until λ becomes sufficiently small,
when the simulations show no evidence of negative diffusion; rather the flow becomes critically slow for very sticky particles.
The higher moments of the simulated flow (e.g. variance) do not show peaks, indicating that hard transitions are not occurring.
Finally, the  simulated density is very close to the average of the boundary densities until λ drops below 1/4, at which point the system fills.





§ DISCUSSION


We have analysed the model using three methods: Direct simulation,
Transition rate matrices and mean field theory. Each has advantages
and drawbacks, but all show a nonequilibrium transition between a
free-flowing and a blocked phase.  In direct simulation this appears
as a peak in the energy and density fluctuations and a change from
λ to λ^4 dependence in the flow rate.  In MFT it
appears as a negative diffusion constant.  In the TRM it corresponds
to a minimum in the minimum eigenvalue, and a crossover in the
character of the associated eigenvector, as clearly illustrated in
Figure.<ref>

The MFT suggests a symmetry between vacancy-type and particle-type
flow at density of 2/3, which is observed in the MC
simulations and TRM calculations.  The density self-organizes to this
value in the limit of strongly repelling particles.

The continuum limit solution to
the MFT is a good predictor of the bulk flow behavior of the SPM for
high λ.  The negative diffusion constant found in MFT at high
stickiness indicates that the assumption of homogeneous density breaks
down: thus the MFT predicts its own demise at a point which agrees well with
the numerics.

Above a certain level of stickiness, the model exhibits a
nonequilibrium phase transition to a slow-flowing phase.  The required
stickiness for the transition is dependent on the strength of the
driving force - large differences between the boundary densities push
the transition towards higher stickiness.  Mean field analysis, together
with visualization of the flowing system, suggest that the transition
occurs when the density becomes inhomogeneous.

The flow exhibits a foamy pattern with time-and-space
correlations of intermediate size.  
The open boundary condition means the system self-organizes its
density.  At zero interaction (λ=1) this must be the mean of
the boundary conditions, and it is unsurprising that high stickiness
means the system fills.  Strong repulsion between the particles leads
to a density of 2/3, which appears to maximize the flow rate
at high λ; there is no known physical principle which requires this.



§ CONCLUSIONS


The sticky particle model is the combination of the Ising model with
the symmetric exclusion process.  It represents the simplest possible
flow model for interacting particles in a homogeneous medium, and as
such is a model for many physical systems where particles are driven
by a pressure or concentration difference across the boundaries.

We have shown that it exhibits complicated and discontinuous behaviour
absent in either 1D Ising or SEP model, specifically, a transformation
from fast to slow flow with stickiness. 

A number of questions remain open.  Is there a physical principle
which determines the density (Fig. <ref>)?  Why does
the strongly repelling system produce a density with maximum flow
(Fig. <ref>)?  Can one derive the λ^4
dependence of the slow-flow (Fig. <ref>)?  These clear,
emergent results from our detailed work strongly suggest that there is a higher
level description of complex flow which could predict
them, but to our knowledge no such unifying theory of nonequilibium
thermodynamics exists.






§ ACKNOWLEDGEMENTS

We would like to thank EPSRC (student grant 1527137) and Wolfson
Foundation and ERC for providing funding, Mikael Leetmaa for producing
, and the  team here at Edinburgh for
maintaining the hardware used.  We would also like to thank Martin
Evans, Bartek Waclaw and Richard Blythe for some very helpful
discussions.



§ LINEAR STABILITY OF MFT SOLUTION.
 Let ρ_0 (x) be a solution to the steady-state continuum MFT, and let us apply a small perturbation δρ(x, t). Let us also assume that ρ_0 is also slowly-varying in x,
in the sense that δρ = o(ρ_0) and ρ_0x = o(δρ), which should be approximately true for a large domain.
Then substituting ρ = ρ_0 + δρ into the continuum-limit MFT and only the highest-order terms, we find that

    δρt =  a^2/τ_0ρ_0 (4 - 3 ζρ_0) ^2 δρx^2 + o(δρ).

Performing a  Fourier transform  (or rather, a suitable local equivalent) with respect to x to yield δ̂ρ̂(k, t) = ℱ[ δρ(x, t) ], this becomes

    δ̂ρ̂t = -k^2 a^2/τ_0ρ_0 (4 - 3 ζρ_0)δ̂ρ̂.

This means that so long as ζ < 3/4, the RHS is always negative, and therefore small perturbations decay exponentially and all is well. However, if ζ>3/4, there may be regions of the solution where
ρ_0 (4 - 3 ζρ_0) becomes negative, causing small perturbations to grow exponentially with an emphasis on high-wavenumber (short-lengthscale) modes, which indicates linear instability.




§ ADDITIONAL FLOW DATA.
 Fig. <ref>, <ref> and <ref> display additional data and information from computations discussed in the main body of the paper. These were performed with KMC using the same run parameters as those used to create Fig. <ref>. 








§ MORE FLOW PLOTS.
 For the interested reader we have included for spacetime flow diagrams, show in Fig. <ref>. When λ=0.05, the medium consists of solid blocks surrounded by empty spaces containing a dilute
gas of particles; at λ=0.35 it is not so different from normal diffusion.
The most interesting images are those for the intermediate λ; here we see a “lumpy” or “foamy” structure, in which small blocks
of particles are being constantly created and destroyed whilst a rather minimal flow occurs across the system.
The simulations did not show any hard phase transition as we vary λ; rather, it seems that this “foamy”
behavior is part of a continuous range between the extremes, containing medium-range correlations between particles.
Unfortunately, computing equal-time correlation functions to the accuracy required
to draw conclusions about these correlations has proven to be extremely difficult, so we cannot find a quantitative description of the foam beyond the flow rate and density data we have discussed in this paper.
In all images in Fig. <ref>, long straight segments of white of black can be seen.  The represent coherent motion at a characteristic velocity given by their gradient. There is nothing in the MFT to suggest what this velocity
should be, and it is much smaller than the simulated system's length divided by the elapsed time,  L/T, thus it must be an emergent property arising from correlated motion of self-assembled regions of  high- or low-density material.
However, it has again proved difficult to analyze this numerically.






§ NON-UNIQUENESS OF THE MFT DENSITY PROFILE.

Recall that, according to the MFT,

    ∫dx  J(x) = (x-x_0)J_0 = -a^2/τ_0ρ[1+ζρ(ρ-2)],

which we would like to solve for ρ(x). We can do this uniquely so long as the right hand side is monotonic for ρ∈ (0, 1). Monotonicity requires that the sign of the derivative of the RHS wrt ρ does not change in that region.
This derivative is of course

    dRHS/dρ   = -a^2/τ_0ρ[1+ζρ(ρ-2)] 
       = -a^2/τ_0[ 3 ζ (ρ-2/3)^2 + (1-4/3ζ) ].

The term in square brackets is always 1 at ρ=0 and λ at ρ=1, which are both positive.
If ζ<0, the term is an n-shaped parabola with both boundaries positive; therefore it is always positive, and there is no sign change.
For ζ>0 (the case of interest for us), we can see from completing the square that the term is always positive unless ζ > 3/4, in which case there is a sign change and therefore a loss of monotonicity.
Thus for ζ > 3/4 (and correspondingly λ<1/4) the MFT does not have a unique steady state solution given a set of Dirichlet boundary conditions, and so we cannot expect the MFT to predict the density profile
in that regime.



§ PARTICLE DENSITY IN BOUNDED DOMAIN AT EXTREME Λ-VALUES.

We scanned across a wide range of λ with three sets of boundary conditions: (3/10, 1/10), (3/4, 1/4) and (9/10, 7/10). The resulting mean flow rate is shown in
Fig. <ref>, and the corresponding density results were already displayed. We once again see the transitions between power law behaviors, as discussed in the main body of the paper. Note that the MFT is never a particularly good
fit for the (3/4, 1/4) configuration; this may be because the difference between the boundaries is greater than in the other cases. The other MFTs are good fits in the high-λ regime until we start reaching λ ∼ 1000,
at which point they seem to start converging to the same flow regardless of boundary conditions. In each case for low-λ we have a power-law regime, each with with an exponent around 4, and then at extreme low-λ we lose
a consistent signal because of noise, which we attribute to lack of system convergence due to critically slow flow.

So far as the density is concerned, at λ=1 in each case it is about what we would expect (the average of the two boundary densities). It is minimized around λ = 0.6, and then as λ is reduced each density appears to 
converge to 1 regardless of boundary conditions (before we encounter the same convergence issue at extreme low-λ). As we go to high-λ, on the other hand, we see a similar convergence, but this time towards ρ∼ 0.7.
This is interesting, as a density of ρ=2/3 is predicted by the MFT to be where maximal flow should occur for λ>1; thus it would appear that the system self-organizes to facilitate maximal flow, which is something which
is hypothesized to occur in other nonequilibrium statmech systems. This also helps to explain the large deviations from the MFT flow predictions we encounter at high-λ, as the system in practice has a different density to the one assumed
by the MFT.













37
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Belitsky and Schütz(2011)]1742-5468-2011-07-P07007
  author author V Belitsky and author G M Schütz, title title Cellular
  automaton model for molecular traffic jams, http://stacks.iop.org/1742-5468/2011/i=07/a=P07007 journal
  journal Journal of Statistical Mechanics: Theory and
  Experiment volume 2011, pages P07007
  (year 2011)NoStop
[Mobilia et al.(2007)Mobilia, Georgiev, and Täuber]Mobilia2007
  author author Mauro Mobilia, author Ivan T. Georgiev,  and author Uwe C. Täuber, title title Phase
  transitions and spatio-temporal fluctuations in stochastic lattice
  lotka–volterra models, 10.1007/s10955-006-9146-3
  journal journal Journal of Statistical Physics  volume 128, pages 447–483 (year 2007)NoStop
[Tegner et al.(2015)Tegner,
  Zhu, Siemers, Saksl, and  Ackland]tegner2015high
  author author B. E. Tegner, author L. Zhu,
  author C. Siemers, author K. Saksl,  and author
  G. J. Ackland, title
  title High temperature oxidation resistance in
  titanium-niobium alloys, 10.1016/j.jallcom.2015.04.115
  journal journal Journal of alloys and compounds  volume 643, pages 100–105 (year 2015)NoStop
[Zhu et al.(2012)Zhu,
  Hu, Yang, and Ackland]zhu2012atomic
  author author Linggang  Zhu, author Qing-Miao  Hu, author Rui Yang,  and author Graeme J. Ackland, title title Atomic-scale
  modeling of the dynamics of titanium oxidation, 10.1021/jp309305n journal journal Journal of
  Physical Chemistry C volume 116, pages
  24201–24205 (year 2012)NoStop
[Deal and Grove(1965)]DealGrove1965
  author author B. E. Deal and author A. S. Grove, title title General
  relationship for the thermal oxidation of silicon, 10.1063/1.1713945 journal journal Journal of
  Applied Physics volume 36, pages
  3770–3778 (year 1965), http://arxiv.org/abs/https://doi.org/10.1063/1.1713945
  https://doi.org/10.1063/1.1713945 NoStop
[Cabrera and Mott(1949)]MottCabrera1949
  author author N Cabrera and author N F Mott, title title Theory of the
  oxidation of metals, http://stacks.iop.org/0034-4885/12/i=1/a=308
  journal journal Reports on Progress in Physics  volume 12, pages 163 (year
  1949)NoStop
[Buzzaccaro et al.(2007)Buzzaccaro, Rusconi, and Piazza]Buzzaccaro2007
  author author Stefano Buzzaccaro, author Roberto Rusconi,  and author Roberto Piazza, title title “sticky” hard
  spheres: Equation of state, phase diagram, and metastable gels, 10.1103/PhysRevLett.99.098301 journal journal Phys. Rev. Lett. volume 99, pages 098301 (year 2007)NoStop
[Ladd et al.(1988)Ladd,
  Colvin, and Frenkel]ladd1988application
  author author Anthony J. C.  Ladd, author Michael E.  Colvin,  and author Daan  Frenkel, title title
  Application of lattice-gas cellular automata to the brownian motion of
  solids in suspension, 10.1103/PhysRevLett.60.975
  journal journal Phys. Rev. Lett. volume 60, pages 975–978 (year
  1988)NoStop
[Liggett(1985)]liggett1985interacting
  author author Thomas M  Liggett, @noop  title
  Interacting particle systems (publisher Springer-Verlag,
  Berlin, year 1985)NoStop
[Ben-Naim et al.(1999)Ben-Naim, Chen, Doolen, and Redner]BenNaim1999
  author author E. Ben-Naim, author S. Y. Chen,
  author G. D. Doolen,  and  author S. Redner, title title Shocklike dynamics of inelastic
  gases, 10.1103/PhysRevLett.83.4069 journal
  journal Phys. Rev. Lett. volume 83,  pages 4069–4072 (year 1999)NoStop
[Shandarin and Zeldovich(1989)]Shandarin1989
  author author S. F. Shandarin and author Ya. B. Zeldovich, title title The
  large-scale structure of the universe: Turbulence, intermittency, structures
  in a self-gravitating medium, 10.1103/RevModPhys.61.185
  journal journal Rev. Mod. Phys. volume 61, pages 185–220 (year
  1989)NoStop
[Frachebourg(1999)]Frachebourg1999
  author author L. Frachebourg, title title Exact
  solution of the one-dimensional ballistic aggregation, 10.1103/PhysRevLett.82.1502 journal journal
  Phys. Rev. Lett. volume 82, pages
  1502–1505 (year 1999)NoStop
[Frachebourg et al.(2000)Frachebourg, Martin, and Piasecki]Frachebourg2000
  author author L. Frachebourg, author P. A. Martin,  and author J. Piasecki, title title Ballistic
  aggregation: a solvable model of irreversible many particles dynamics,  10.1016/S0378-4371(99)00585-3 journal
  journal Physica A Statistical Mechanics and its Applications  volume 279, pages 69–99 (year 2000), http://arxiv.org/abs/cond-mat/9911346
  cond-mat/9911346 NoStop
[Obukhovsky et al.(2017)Obukhovsky, Kutsyk, Nikonova, and  Ilchenko]Obukhovsky2017
  author author Vyacheslav V.  Obukhovsky, author Andrii M.  Kutsyk, author Viktoria V.  Nikonova,  and author Oleksii O. Ilchenko, title title Nonlinear diffusion in multicomponent liquid solutions,  10.1103/PhysRevE.95.022133 journal journal Phys. Rev. E volume 95, pages 022133 (year 2017)NoStop
[Gorokhova and Melnik(2010)]Gorokhova2010
  author author N. V. Gorokhova and author O. E. Melnik, title title Modeling of the
  dynamics of diffusion crystal growth from a cooling magmatic melt, 10.1134/S0015462810050017 journal journal Fluid Dynamics volume 45, pages 679–690 (year 2010)NoStop
[Kardar et al.(1986)Kardar,
  Parisi, and Zhang]PhysRevLett.56.889
  author author Mehran Kardar, author Giorgio Parisi,  and author Yi-Cheng Zhang, title title Dynamic scaling
  of growing interfaces, 10.1103/PhysRevLett.56.889
  journal journal Phys. Rev. Lett. volume 56, pages 889–892 (year
  1986)NoStop
[Krug and Spohn(1988)]PhysRevA.38.4271
  author author J. Krug and author H. Spohn,  title title Universality classes for
  deterministic surface growth, 10.1103/PhysRevA.38.4271
  journal journal Phys. Rev. A volume 38, pages 4271–4283 (year
  1988)NoStop
[Sasamoto and Spohn(2010)]Sasamoto2010
  author author Tomohiro  Sasamoto and author Herbert  Spohn, title title
  One-dimensional kardar-parisi-zhang equation: An exact solution and its
  universality, 10.1103/PhysRevLett.104.230602 journal journal Phys. Rev. Lett. volume 104, pages 230602 (year
  2010)NoStop
[Sugden and Evans(2007)]sugden2007dynamically
  author author K E P Sugden and author M R Evans, title title A dynamically
  extending exclusion process, http://stacks.iop.org/1742-5468/2007/i=11/a=P11013 journal
  journal Journal of Statistical Mechanics: Theory and
  Experiment volume 2007, pages P11013
  (year 2007)NoStop
[Kollmann(2003)]Kollmann2003
  author author Markus Kollmann, title title Single-file
  diffusion of atomic and colloidal systems: Asymptotic laws, 10.1103/PhysRevLett.90.180602 journal journal Phys. Rev. Lett. volume 90, pages 180602 (year 2003)NoStop
[Lin et al.(2005)Lin,
  Meron, Cui, Rice, and  Diamant]Lin2005
  author author Binhua Lin, author Mati Meron,
  author Bianxiao Cui, author Stuart A. Rice,  and author Haim Diamant, title
  title From random walk to single-file diffusion,  10.1103/PhysRevLett.94.216001 journal
  journal Phys. Rev. Lett. volume 94,  pages 216001 (year 2005)NoStop
[Hegde et al.(2014)Hegde,
  Sabhapandit, and Dhar]Hegde2014
  author author Chaitra Hegde, author Sanjib Sabhapandit,  and author Abhishek  Dhar, title title
  Universal large deviations for the tagged particle in single-file motion,  10.1103/PhysRevLett.113.120601 journal
  journal Phys. Rev. Lett. volume 113,  pages 120601 (year 2014)NoStop
[Krapivsky et al.(2014)Krapivsky, Mallick, and Sadhu]Krapivsky2014
  author author P. L. Krapivsky, author Kirone Mallick,  and author Tridib Sadhu, title title Large deviations
  in single-file diffusion, 10.1103/PhysRevLett.113.078101
  journal journal Phys. Rev. Lett. volume 113, pages 078101 (year
  2014)NoStop
[Imamura et al.(2017)Imamura, Mallick, and Sasamoto]Imamura2017
  author author Takashi Imamura, author Kirone Mallick,  and author Tomohiro Sasamoto, title title Large
  deviations of a tracer in the symmetric exclusion process, 10.1103/PhysRevLett.118.160601 journal journal Phys. Rev. Lett. volume 118, pages 160601 (year 2017)NoStop
[Katz et al.(1984)Katz,
  Lebowitz, and Spohn]Katz1984
  author author Sheldon Katz, author Joel L. Lebowitz,  and author Herbert Spohn, title title Nonequilibrium
  steady states of stochastic lattice gas models of fast ionic conductors,  10.1007/BF01018556 journal journal Journal of Statistical Physics volume
  34, pages 497–537 (year 1984)NoStop
[Zia(2010)]Zia2010
  author author R. K. P.  Zia, title title Twenty
  five years after kls: A celebration of non-equilibrium statistical
  mechanics, 10.1007/s10955-009-9884-0 journal journal Journal of Statistical Physics volume 138, pages 20–28 (year
  2010)NoStop
[Kafri et al.(2003)Kafri,
  Levine, Mukamel, Schütz,  and Willmann]Kafri2003
  author author Y. Kafri, author E. Levine,
  author D. Mukamel, author G. M. Schütz,  and author R. D. Willmann, title title Phase-separation transition in
  one-dimensional driven models, 10.1103/PhysRevE.68.035101 journal journal Phys.
  Rev. E volume 68, pages 035101
  (year 2003)NoStop
[Kawasaki(1966)]PhysRev.145.224
  author author Kyozi Kawasaki, title title Diffusion
  constants near the critical point for time-dependent ising models. i,  10.1103/PhysRev.145.224 journal journal Phys. Rev. volume 145, pages
  224–230 (year 1966)NoStop
[Spohn(1983)]spohn1983
  author author H Spohn, title title Long range
  correlations for stochastic lattice gases in a non-equilibrium steady
  state, http://stacks.iop.org/0305-4470/16/i=18/a=029 journal journal Journal of Physics A: Mathematical and
  General volume 16, pages 4275
  (year 1983)NoStop
[Leetmaa and Skorodumova(2014)]leetmaa2014kmclib
  author author M. Leetmaa and author N. V. Skorodumova, title title KMCLib:
  A general framework for lattice kinetic Monte Carlo (KMC) simulations,  10.1016/j.cpc.2014.04.017 journal journal Computer Physics Communications volume
  185, pages 2340–2349 (year 2014), http://arxiv.org/abs/1405.1221 arXiv:1405.1221 [physics.comp-ph]
  NoStop
[Gillespie(1977)]Gillespie1977
  author author Daniel T.  Gillespie, title title
  Exact stochastic simulation of coupled chemical reactions, 10.1021/j100540a008 journal journal
  The Journal of Physical Chemistry volume 81,  pages 2340–2361 (year 1977), http://arxiv.org/abs/http://dx.doi.org/10.1021/j100540a008
  http://dx.doi.org/10.1021/j100540a008 NoStop
[Bortz et al.(1975)Bortz,
  Kalos, and Lebowitz]Bortz1975
  author author A.B. Bortz, author M.H. Kalos,   and author J.L. Lebowitz,  title title A new algorithm for monte
  carlo simulation of ising spin systems, https://doi.org/10.1016/0021-9991(75)90060-1 journal journal Journal of Computational Physics volume
  17, pages 10 – 18 (year 1975)NoStop
[Prados et al.(1997)Prados,
  Brey, and Sánchez-Rey]Prados1997
  author author A. Prados, author J. J. Brey,   and author B. Sánchez-Rey,  title title A dynamical monte carlo
  algorithm for master equations with time-dependent transition rates,  10.1007/BF02765541 journal journal Journal of Statistical Physics volume
  89, pages 709–734 (year 1997)NoStop
[Hellier(2018)]jHellGitRepo
  author author Joshua Hellier, title title
  joshuahellier/phdstuff: Codes for sticky particles in 1d and 2d, 10.5281/zenodo.1162818   (year 2018),  10.5281/zenodo.1162818NoStop
[Robert(2015)]metHastAlg
  author author Christian P.  Robert, title The metropolis-hastings
  algorithm, in 10.1002/9781118445112.stat07834 booktitle Wiley StatsRef: Statistics Reference Online  (publisher American Cancer Society, year 2015)  pp. pages 1–15NoStop
[Blythe and Evans(2007)]blytheEvans2007
  author author R A Blythe and author M R Evans, title title Nonequilibrium
  steady states of matrix-product form: a solver's guide, http://stacks.iop.org/1751-8121/40/i=46/a=R01 journal
  journal Journal of Physics A: Mathematical and Theoretical  volume 40, pages R333 (year
  2007)NoStop
[Appert-Rolland et al.(2015)Appert-Rolland, Ebbinghaus, and Santen]AppertRolland2015
  author author Cécile  Appert-Rolland, author Maximilian Ebbinghaus,  and author
  Ludger Santen, title
  title Intracellular transport driven by cytoskeletal
  motors: General mechanisms and defects,  (year
  2015)NoStop



