Condensation in stochastic mass transport models: beyond the zero-range process
We consider an extension of the zero-range process to the case where the hop rate depends on the state of both departure and arrival sites. We recover the misanthrope and the target process as special cases for which the probability of the steady state factorizes over sites. We discuss conditions which lead to the condensation of particles and show that although two different hop rates can lead to the same steady state, they do so with sharply contrasting dynamics. The first case resembles the dynamics of the zero-range process, whereas the second case, in which the hop rate increases with the occupation number of both sites, is similar to instantaneous gelation models. This new “explosive” condensation reveals surprisingly rich behaviour, in which the process of condensate’s formation goes through a series of collisions between clusters of particles moving through the system at increasing speed. We perform a detailed numerical and analytical study of the dynamics of condensation: we find the speed of the moving clusters, their scattering amplitude, and their growth time. We finally show that the time to reach steady state decreases with the size of the system.
pacs:89.75.Fb, 05.40.-a, 64.60.Ak
Date July 5, 2019
Condensation is one of the most ubiquitous forms of phase transition experienced in everyday life. Water vapour condenses into a liquid on cold surfaces such as windows or a cold beverage taken out of the fridge. Fog is an aerosol that occurs when temperature drops below the dew point and water vapour condenses into tiny droplets suspended in air. Condensation can also be seen as small “clouds” in low-pressure zones above aircraft wings when flying in damp weather conditions. Condensation is also extensively used in industry for making liquid oxygen and nitrogen, production of liquefied natural gas, or in petroleum distillation.
All these examples describe a transition from a gas to a liquid state of matter. However, in statistical physics, the notion of condensation is more widely used to describe any process in which a finite fraction of some conserved quantity becomes localized in the phase space. So understood, condensation can occur either in real space or in momentum space. In fact, the most prominent example of statistical-physics condensation - Bose-Einstein condensation - takes place in the momentum space where a finite fraction of all particles present in the system assumes the lowest energy state. Such generalized condensation occurs in many different systems, both in- and out-of-equilibrium: granular clustering, wealth condensation , quantum gravity , polydisperse hard spheres , hub formation in complex networks , or even jamming in traffic flow [5, 6].
Perhaps surprisingly, many features of the generalized condensation transitions observed in the systems listed in the previous paragraph, can be captured by a very simple model called the zero-range process (ZRP) . ZRP is a well-studied, non-equilibrium model of particles hopping between sites of a lattice. In the simplest one-dimensional asymmetric version of this model, a particle moves from site to of a one-dimensional periodic lattice with rate , where is the occupancy of the departure site . As the rates are totally asymmetric a current always flows and detailed balance cannot be satisfied, thus the steady state is non-equilibrium in nature. When the hop rate decays slowly enough with and the density of particles is larger than a certain critical density (the exact value of which depends on the function ), a finite fraction of all particles condenses onto a single, randomly chosen site.
The advantage of ZRP is that its non-equilibrium steady state has a simple factorised form which is amenable to analysis and has furnished our understanding of the condensation transition. This is why the model is often used as an effective description of more complicated systems, for example systems with more microscopic degrees of freedom, where such a complete analytical description is often impossible.
The characteristic feature of the ZRP is that the hop rate depends only on the departure site and not on the destination site. A natural generalisation is to consider a hop rate where is the occupancy of the departure site and is the occupancy of the destination site . Such processes are sometimes referred to as migration processes in the mathematical literature . An example of such process is the ‘misanthrope process’  whose name originally derived from considering a hop rate as an increasing function of and decreasing function ; an effective repulsion between individuals results. The name ‘misanthrope process’ is now generally used for any hop rate of the form . Such a process has been used, for example, to model link rewiring in complex network models . Another example is the target process introduced in Ref. . Here, the rate , with some function , depends on the occupation of the arrival site.
Our aim in this paper is to analyse a more general case of . It has been shown that the misanthrope process has a factorised steady state when certain conditions on (expressed below in Eqs. (14) and (15)) are satisfied. Here we provide examples of hop rates when these conditions are satisfied. A particularly interesting case that we discuss is when the hop rate factorizes:
We also analyse the conditions for condensation to occur. As the steady state has a factorised form the conditions on the weights will be the same as for the zero-range process. However, there is more freedom in the misanthrope dynamical rules with which to realise this asymptotic behaviour. Moreover, the existence of condensation depends not only on the asymptotic behaviour of as it does for the ZRP, but also on which gives the hop rate to an empty site. In the case of factorized hop rate (1), the condensation criterion turns out to depend on .
We also expand on our previous work  and study the dynamics of condensation for two different choices of the factorized hop rate (1). Both forms of lead to the same steady-state probability distribution but contrasting dynamical behaviour. A particularly interesting case, in which increases with both and , leads to “explosive condensation” where the time to condensation decreases with increasing system size. This is very similar to “instantaneous gelation” known from the theory of coagulation processes, and we discuss similarities and differences between the two models in the last section.
Beyond factorised steady states it is known that pair-factorised states (also referred to as Markov measures) may occur in some stochastic models. A simple example is the 1+1d solid-on-solid model [13, 14, 15] with a “pinning potential”, which describes the height of the interface between two phases and whose microstate probability reads
with two parameters . The -term describes the tendency of the system to form straight interface (“surface tension”) whereas gives the energy cost for the interface having zero height. The model can be used for example to describe the formation of a new layer of atoms on a solid substrate. The probability (2) can be recast into
where is the pairwise weight. Although this model is an equilibrium model, more recently a class of nonequilibrium models which exhibit such steady states have also been proposed . One could hope that if conditions (14), (15) are not satisfied there might be some choices of rates which yield pair-factorised states. We shall show that actually this is not the case and that either (14), (15) is satisfied or the steady state has an unknown form.
2 Misanthrope process
We now define the model that we consider. particles reside on sites of a one-dimensional closed, periodic chain of length . Each site carries particles, so that . A particle hops from site to site with rate which depends on the occupancies of both the departure and the arrival site, see Fig. 1. Periodicity implies .
The misanthrope process can be mapped to the total asymmetric simple exclusion process (TASEP) on a closed chain. In Fig. 2 we show how this can be done. In TASEP, each site of a periodic 1d lattice is occupied by at most one particle. We can identify the number of particles in the misanthrope process with the number of vacancies between the particles in TASEP. To obtain the same dynamics as in the misanthrope process, we assume that the particle jumps to the left with rate depending on the number of empty sites in front () and behind () the particle.
It is known that under certain conditions the steady state probability of the misanthrope process is given by a factorised form
where is a certain non-negative weight function and is the normalization
In (4) the Kronecker delta imposes the constraint that the total mass (number of particles) in the system is conserved. The conditions for factorisation can be derived as follows. In the steady state the equation balancing probability currents from and to a given configuration is
where is assumed to be given by Eq. (4). Here we have taken the conventions that (the hop rate is zero if there are no particles at the departure site), , and . Dividing both sides by and shifting indices in terms on the right-hand, we obtain
In Appendix A we show that the most general solution to an equation of the form
reads for . Therefore, from Eq. (7) we obtain that
with some (yet unknown) function . Let us consider now the case . Because is identically zero, the second term in the above equation vanishes and we obtain
Similarly, for and arbitrary , the first term of Eq. (9) vanishes and we find
The last relation (11) may be rearranged to give a recursion relation for as follows
It can be proved (see Appendix B) that this actually reduces to two conditions:
To obtain insight into the physical meaning of conditions (14,15), let us examine the matrix . The first row is made up of zeros since . Let us define two vectors and of infinite length which specify the second row and the first column, i.e.,
The matrix assumes now the following form:
By iterating these equations one obtains unambiguous expressions for all , for example, given for , (20) implies and (19) gives . Thus the hop rate is fully determined by fixing the rate of hopping from a site with particles to an empty site, and the rate of hopping from a site with only one particle to a site with particles. There is, however, an additional condition that must be non-negative for all . This imposes some constraints on which we were not able to express in a closed form and hence we do not have an expression for the most general form of that implies non-negative . In what follows we will consider some special cases of for which one can prove non-negativeness of .
It is important to note that any exponential factor in does not change the steady-state properties since it appears in as a constant prefactor due to the fixed number of sites, , and particles . Thus we may neglect the factor in Eq. (21) and write
2.1 Lack of Pair-factorized states
One may wonder whether conditions (14,15) could be relaxed if, for example, instead of insisting on the factorization of the steady state over sites we assumed a much less-restrictive factorization over pairs of sites (Eq. (3)) as in Ref. [16, 19]. There, the hop rate depended on the occupation number of the departure site as well as two neighbouring sites . Surprisingly, it turns out that for the misanthrope process, whose hop rate depends only on the departure and arrival site, pair-factorization is not possible except for a trivial case when the steady state fully factorizes as in Eq. (4). We give a proof of this result in Appendix C.
3 Physical examples of the hop rate
3.1 Factorized hop rate
One can check that for a factorized hop rate of the form
with some arbitrary, non-zero constant . Thus
is the form of a factorised hop rate that yields a factorised steady state. For non-negativeness of , it suffices that one of the following conditions is fulfilled:
, for ,
In both these cases is always positive. One could additionally consider cases where is always negative but as is invariant under these are trivially equivalent to the two cases above.
Our first simple example of the factorised hop rate (25) is and with some integer . This leads to
If this reduces to the asymmetric simple exclusion process where the occupancy of each site is limited to 1 and . Similarly the case of general integer corresponds to ‘partial exclusion’  where each site of a lattice contains at most particles and each particle attempts hops forward to the next site with rate one which succeed with probability where is the occupancy of the destination site.
The case for has been studied and is referred to as the Inclusion Process. In a prescribed limit , the limit of zero hopping rate onto empty sites, a form of condensation occurs .
3.2 Non-factorized choices of
The factorized form (25) is not the only allowed form for that solves (19,20) and gives rise to a factorized steady state. To see this, let us consider an example in which , with . In this case, one can check that
The corresponding weight function (22) reads for
where in the last expression we have suppressed a term exponential in which as, previously noted, is unimportant. This expression assumes a form similar to the weight function in the zero-range process. In fact, this model is closely related to a recently considered facilitated exclusion process . In that work an exclusion process is considered wherein, for a particle to hop, the neighbouring site behind has to be occupied. Using the standard mapping of section 1 (see Fig. 2) the facilitated exclusion process corresponds to the limit and of the misanthrope process.
4 Condensation in the misanthrope process
An important feature of models such as ZRP  or balls-in-boxes model , in which the steady-state probability takes the factorized form (4), is that for some choices of the weight function the system exhibits a phase transition in the thermodynamic limit . If the density of particles is above some critical value , the surplus of particles accumulates at a single site and is called the condensate. The sufficient condition for the condensation above some finite is the appropriate asymptotic behaviour of . As any exponential dependence in is irrelevant (see previous section), we just need to consider the asymptotic behaviour of modulo any exponential factor. There are two generic cases (see e.g. ):
with . The critical density is finite but its numerical value depends on the particular form of and not only on its asymptotic behaviour. The fraction of all particles goes to the condensate. We will refer to this as standard condensation.
increases with more quickly than exponentially, e.g., as . This leads to so called strong (or complete) condensation - the critical density and a fraction of particles tending to one in the thermodynamic limit is located at one site.
There are also some other specific examples, such as the Backgammon model [25, 26], which corresponds to taking and exhibits condensation in the limit , and the inclusion process in the limit mentioned earlier .
To see why the condensation happens in the two generic cases highlighted above, we shall follow a standard approach . Treating the steady-state probability as the statistical weight of a given configuration, and defining the grand-canonical partition function
we see that the phase transition, signalled by a singularity of at some , is possible only if the series has a finite radius of convergence . Moreover, the density calculated as a function of fugacity from the grand-canonical partition function,
must be finite as . This is only possible if decays as a power law in which case we may have a finite (case I) or grows very fast with in which case and (case II). Any exponential factor in only shifts the position of .
4.1 Factorized hop rate
We will now discuss which choices of the hop rate of the misanthrope process lead to the generic cases I, II described above. We begin by considering factorized hop rates , as in Eq. (25). Using
equation (22) yields
where we have removed a constant factor that generates an exponential factor in . Let us first consider strong condensation (case II), that is for large . From Eq. (34) we see that for this to happen, has to tend to in the limit of . For example, for
we obtain and strong condensation occurs.
In standard condensation (case I) we require a power law decaying as
Expression (34) fulfils this condition for two distinct large , asymptotic behaviours of :
where is a constant. Therefore the asymptotic behaviour of the single-site weight is where the exponent reads
Condensation is possible if is larger than 2. This leads to the following conditions on , and for condensation to occur:
Interestingly, the existence of condensation depends not only on the asymptotic behaviour of but also on . This should be contrasted with the ZRP for which it is only the asymptotic decay of the hop rate, given by in , that determines whether condensation is possible or not and the exact form of affects only the value of the critical density, above which condensation occurs. Here, the influence of the particular form of is much more significant.
4.2 Illustrative example of a factorized hop rate
As an illustrative example consider the hop rate
which differs from Eq. (35) only in the value of . In the previous case, we have strong condensation, whereas in Eq. (44) corresponds to from Eq. (39). Thus, by (43), for there is no condensation, but for there is standard condensation and for there is strong condensation, even though is the same in all cases for .
It is possible to compute exactly the critical density in this case by using (22), which yields
after removing exponential factors in . Here
and we use the Pochhammer symbol . Then the generating function (31) is a hypergeometric function
which converges for . Using standard identities for the hypergeometric function (which hold when the series converge) ,
one may determine the critical density as the maximum allowed density (32), given by :
We see indeed that as and as .
4.3 General case with arbitrary
In order to have strong condensation (case II) we require for large . For standard condensation with (36) (case I) we have again two distinct possible asymptotic behaviours for :
which gives , and as the necessary condition for condensation, and
5 Dynamics of the condensate
We return now to the different microscopic dynamics that yield the same steady state and exhibit condensation above the same . To gain insight into the dynamics, we simulated this system using a modified continuous time Monte Carlo in which the time interval to the next event is stochastically generated. We prepare the system in a state in which particles are randomly distributed among the sites. Then, in each time step we choose a random site and move a particle from this site to its right neighbour. The site is selected at random with probability . To do this effectively, we sort all possible “moves” with respect to their probabilities of occurrence. Lastly, we update the time by adding an exponentially distributed random variable with mean . To minimize computational time, after each move we recalculate only the rates , , and which have been affected by the move. We also use tabulated values of to avoid time-consuming algebraic operations such as calculating for non-integer .
5.1 ZRP-like rates given by Equation (39)
and . Then taking in (25), yields the asymptotic behaviour
for sufficiently large. In this figure, the positions of the five most-occupied sites are plotted against time, together with snapshots of the system at three different times. We see that at first, small condensates are formed (multiple horizontal lines in the upper diagram). These condensates coalesce quickly until only two are left. Finally, the two condensates merge into a single one; this last process is the slowest. The size of the condensate initially grows quickly and then slows down, see Fig. 4, left. This is very similar to the condensation in the zero-range process with , which also leads to the power-law . The dynamics of such ZRP has been extensively studied in the past [28, 29]. It has been argued that the typical time to reach the steady state in an asymmetric 1d system scales as with the number of sites if the density is kept constant and are sufficiently large.
One can argue that the same scaling should hold in our system, that is the -dependence of in Eq. (55) is not important. The argument goes as follows. In the last, slowest stage of condensation when only two condensates are left, occupation numbers on sites surrounding these condensates fluctuate very rapidly. Thus, from each condensate’s point of view, the hop rates to and from the condensate at site are effectively averaged over . According to formula (55), the average inflow rate is while the outflow rate reads . The net flow between both condensates having and particles, respectively, is . Since both condensates have particles, the flow is at most . Thus, statistical fluctuations which are of order dominate and perform effectively unbiased random walks, at least as long as they are of order . The typical time scale is defined by the time it takes to reach either or , which scales diffusively as . The finite distance between the condensates does not matter, because due to the finite average hop rate in the bulk, the particles traverse it in steps which is smaller than the time scale related to fluctuations.
To confirm this prediction, we carried out simulations for and various sizes , and measured the time it took before the maximal occupation number reached the mean steady-state value . One clearly sees from Fig. 6, left, that scales approximately as , independently of and other parameters.
5.2 Explosive-condensation rates given by Equation (38)
We now turn to the case from Eq. (38), for which for large . To be more specific, we choose so that
with (typically) . In Fig. 5 we show results of computer simulations. If we compare it with Fig. 3, we see some striking differences. First, the site carrying the maximal number of particles moves unidirectionally, in the same direction as hopping particles. The motion of the condensate is similar to the “slinky”-like motion of a non-Markovian model [30, 31] or a generalization of ZRP to non-factorizing probabilities . Second, the evolution speeds up in time; the condensate moves faster and faster as it gains particles, see also Fig. 4, right. This is why in the previous work  we termed it “explosive condensation”. Third, the speed at which the condensate travels through the system stabilizes after the system reaches the steady state. Finally, other (smaller) peaks move in the opposite direction to the main condensate each time they interact with it.
The dynamics thus differs significantly from the ZRP-like case. Figure 6, right, shows another striking feature of this new condensation. If we fix the density of particles and increase the system size , the time to reach steady state decreases as a function of . This is in contrast to the ZRP-like behaviour from Fig. 6, left, and to the ordinary ZRP, where the time to steady state increases with . In Ref.  we focused on deriving the relation between and , assuming that the time scale was dominated by the slow process of initial coalescence, in contrast to the ZRP where this process is the fastest one. Our reasoning was based on several assumptions on how clusters of particles interact, which we inferred from numerical simulations. In what follows we shall provide more evidence in support of these arguments.
First, we show that a cluster of particles of size moves through the system at speed . As it moves, the cluster collides with other clusters and exchanges particles with them. On average, particles flow from smaller to larger clusters. Next, we show that the rate of accumulation increases with increasing speed of a cluster and, as a result, it takes a finite time for the cluster to reach macroscopically large occupation (condensate). We also show that the distribution of this time is quite broad and arbitrarily short times are possible (although not very likely). Finally, we show that, since there are initially clusters and each of them can potentially become the condensate in a finite time, the time to reach steady state (condensation) decreases with as a result of a simple extreme value statistics argument.
We shall now discuss our arguments in detail, assuming that . Let us first calculate the speed of an isolated cluster of particles travelling in an otherwise empty system. We assume that the cluster is initially at site , so that , and . The only allowed transition for the cluster is to lose one particle to site , which happens with rate . The average time for this to happen is . In the second step, the condensate may lose another particle with rate , which will take time . Alternatively, the particle at site may jump to site , but this is much less likely because and this is much smaller than . The process in which site loses particles in favour of site will thus continue with rate (with other processes having negligible rates) until and . The average times for each step will be , where . Therefore, the total time of the process will be the sum of times :
For which is the case of condensation and for , we can approximate the above sum by the last term only, and we obtain that . This shows that a cluster with particles moves sites per unit time. In particular, the speed of the condensate of size is and increases as a power of the system size . Although in our analysis we have neglected a small probability of losing particles by the cluster as it moves through the system, this is justified for both small clusters and the macroscopic cluster (the condensate). For small clusters of size (initial stage of condensation), the distance between the clusters is of order and the probability of losing a particle between collisions is negligible. For the condensate (), it is possible to lose particles but as it moves through the system, the condensate will also gain particles from the background and its size will fluctuate only a little. The above formula is in excellent agreement with simulations; for example, for parameters from Fig. 5, the measured speed is sites per unit time, whereas the formula gives .
We shall now analyse what happens in a collision between two clusters of size and . The process can be conveniently studied in numerical simulations by performing “scattering experiments” in which a larger condensate collides with a slower, smaller condensate , see Fig. 7. We prepare the system in a state such that and are initially located at sites and , respectively, and we simulate the system using the same kinetic Monte-Carlo algorithm as before until one of the condensates reaches site . We then measure its mass and calculate the mass transferred in any collision111In principle, no collision may occur but the probability of this is vanishingly small for large enough from the smaller to the larger condensate: . In Fig. 8 we show that although the distribution of looks almost symmetrical, it is biased toward positive values. The average and depends only weakly on the size of the clusters involved in the collision and on the value of . The average mass transfer does not also depend on the value of , although is necessary to have a non-zero hop rate for a single particle. Thus, in our further analysis below, we shall assume whenever possible as this greatly simplifies calculations.
5.2.1 Deterministic description of scattering
To understand the numerical results, we assume (as before) that the system is completely empty, save for the two clusters. As the larger cluster approaches the smaller cluster , the collision is initiated when and become separated by just a single site with one particle that has hopped out from , see Fig. 9. The initial configuration is therefore: particles at the left site, particle at the middle site, and particles at the right site. In the collision process, the left site loses particles to the middle site, which in turn loses particles in favour of the right site, with the transition rates given by and , respectively, where is the total number of particles in the system (this is a conserved quantity). As a first step, it is useful to describe the dynamics of this process using a deterministic approach in which stochasticity is neglected. Our expectation would be that this description should be valid for large clusters, for which random fluctuations of do not matter. The equations for the average read
with the initial condition . We note that
which means that the quantity
is a conserved quantity determined by the initial condition. We also note that upon time reversal , the variables and exchange their roles:
and thus if we change and , we will recover Eqs. (58)-(59). This means that the solution for must be equal to the initial condition at , and hence , (where is small). In Fig. 10 we plot numerical solutions of Eqs. (58)-(59) which confirm our prediction. The deterministic, continuous approximations predicts zero mass transfer because the size of the larger cluster is the same as before the collision. This is at variance with the results of stochastic simulations, even for very large , hence we see that stochasticity is important.
5.2.2 Stochastic description of scattering
It is convenient to express the evolution of the fully stochastic problem not in real time , but as a function the number of particles that have hopped since the beginning of the process. Assume that corresponds to the state . From the deterministic problem and Fig. 10 we know that the occupation of the middle site first rises and then decreases again. The collision ends when the middle site reaches zero occupation, so that the final state is with a random variable. The mass transferred from the smaller to the larger condensate in the collision is given by . Let us denote by the probability that the state of the three sites is after hops. To calculate the mass transferred in the process, we need to calculate
where is the distribution of the mass transferred. The probability can be found recursively:
with the initial condition . In Fig. 11 we plot calculated numerically using the above formula and Eq. (65), for different sizes of colliding condensates, and . We can see that for , the mass transferred in the collision tends to a constant value of about , depending on . In Fig. 12 we show that stabilises at for and . We also see in Fig. 11 that for , the average mass calculated from the recursion relation is actually negative, as if particles flew from the larger to the smaller condensate. This does not agree with Monte-Carlo simulations, see Fig. 8, which indicate that , even for a small difference between . Since our numerical calculation is exact, this discrepancy implies that in the simulation, clusters collide more than once, until the right cluster becomes sufficiently large compared to the left cluster, and manages to escape.
5.2.3 Coalescence process of clusters
Our simulations and numerical calculations thus indicate that we can describe the relaxation towards steady state as a process of coalescence of small clusters into the condensate in a series of collisions. In each collision, only a few particles are exchanged on average, and the net flow of particles is from smaller to larger clusters. Let us analyse what happens if the system is initially prepared in a state with clusters separated by runs of empty sites. Numerically, this can be realised by placing clusters at sites , where is some fixed number larger than 1 and independent of . This initial state is not the same as the state with randomly distributed particles that we typically use in simulations, but the existence of well-defined clusters at all times makes the effective description from the previous subsections (in terms of collisions between isolated clusters) applicable to the whole process, from the initial to the final state. Moreover, on a large scale the system still looks homogeneous, as if the particles were distributed uniformly.
Let us first look at the evolution of one cluster. It collides with other clusters, gains or loses particles, and finally it either dissolves into the background or it reaches the steady-state occupation and becomes the condensate. Let be the time it takes to reach this steady-state size for this particular condensate ( if the cluster dissolves into the background). differs among different clusters and different realisations of the system, and it can be thought of as a random variable with some distribution . In Figure 13 we show a hypothetical form of such a distribution. Since all clusters are initially identical, is the probability density function of a randomly chosen cluster having evolved into a condensate after time . We have, however, such clusters and each of them should have the same chance to become the final condensate, although only one of them (the fastest evolving one) will make it. Therefore, the time to reach steady state will be the minimal time out of times for individual clusters:
and we can estimate it as