Numerical study of condensation in a Fermi-like model of conterflowing particles via Gini coefficient

TL;DR

This paper investigates a Fermi-like model of counterflowing particles, analyzing how asynchronous updates influence clogging transitions and demonstrating the Gini coefficient's effectiveness in studying mobility phenomena.

Contribution

It introduces a detailed comparison of synchronous and asynchronous dynamics in a Fermi-like particle model and applies the Gini coefficient to analyze clogging transitions.

Findings

Asynchronous updates support the mobile-clogging transition.

The transition depends on initial conditions.

Gini coefficient effectively characterizes clogging phenomena.

Abstract

The collective motion of self-driven particles shows interesting novel phenomena such as swarming and the emergence of patterns. We have recently proposed a model for counterflowing particles that captures this idea and exhibits clogging transitions. This model is based on a generalization of the Fermi-Dirac statistics wherein the maximal occupation of a cell is used. Here we present a detailed study comparing synchronous and asynchronous stochastic dynamics within this model. We show that an asynchronous updating scheme supports the mobile-clogging transition and eliminates some mobility anomalies that are present in synchronous Monte Carlo simulations. Moreover, we show that this transition is dependent upon its initial conditions. Although the Gini coefficient was originally used to model wealth inequalities, we show that it is also efficient for studying the mobile-clogging transition. Finally, we compare our stochastic simulation with direct numerical integration of partial differential equations used to describe this model.

Equations

$$\begin{array}[]{ccc}\frac{\partial\rho_{A}}{\partial t}&=&K_{A}\frac{\partial}{\partial x}\left[f_{A}(\rho_{A},\rho_{B})\rho_{A}\right]+D_{A}\frac{\partial^{2}}{\partial x^{2}}\rho_{A}\\ &&\\ \frac{\partial\rho_{B}}{\partial t}&=&K_{B}\frac{\partial}{\partial x}\left[f_{B}(\rho_{A},\rho_{B})\rho_{B}\right]+D_{B}\frac{\partial^{2}}{\partial x^{2}}\rho_{B}\end{array}$$

$$\frac{\partial {\rho }_{A,B}}{\partial t}=K_{A,B}\frac{\partial {\rho }_{A,B}}{\partial x}+D_{A,B}\frac{{\partial }^2{\rho }_{A,B}}{\partial x^2}.$$

$$\begin{array}[]{ccc}\frac{\partial\rho_{A}}{\partial t}&=&-k_{1}\frac{\partial\rho_{A}}{\partial x}+k_{2}\frac{\partial}{\partial x}\left(\frac{\rho_{A}\rho_{B}}{\rho_{A}+\rho_{B}}\right)\\ &&\\ \frac{\partial\rho_{B}}{\partial t}&=&k_{1}\frac{\partial\rho_{B}}{\partial x}-k_{2}\frac{\partial}{\partial x}\left(\frac{\rho_{A}\rho_{B}}{\rho_{A}+\rho_{B}}\right)\end{array}$$

$$f_{A,B}({\rho }_A,{\rho }_B)=\frac{{\rho }_{A,B}}{{\rho }_A+{\rho }_B}.$$

$$f_{A,B}({\rho }_A,{\rho }_B)=\frac{1}{1+e^{\alpha ({\rho }_A+{\rho }_B-{\sigma }_{\max })}},$$

$${\rho }_A(m,n)=p_{m-1,m}^{(n-1)}{\rho }_A(m-1,n-1)+p_{m,m}^{(n-1)}{\rho }_A(m,n-1)\text{,}$$

$$\begin{array}[]{lll}\rho_{A}(m,n)-\rho_{A}(m,n-1)&=&p_{m-1,m}^{(n-1)}\rho_{A}(m-1,n-1)\\ &&\\ &&-\ p_{m,m+1}^{(n-1)}\rho_{A}(m,n-1)\end{array}$$

$$p_{i,j}^{(n)}=\frac{1}{1+e^{\alpha ({\sigma }_{j,n}-{\sigma }_{\max })}}.$$

$$p_{i,j}^{(n)}=\left\{\begin{array}[]{ll}1&\text{if }\sigma_{j,n}=\rho_{A}(j,n)+\rho_{B}(j,n)<\sigma_{\max}\\ &\\ 1/2&\text{if }\sigma_{j,n}=\sigma_{\max}\\ &\\ 0&\text{otherwise}\end{array}\right.$$

$$\begin{array}[]{lll}a_{m,n}&=&\rho_{A}(m,n)/\left[1+e^{\alpha(\rho_{A}(m+1,n)+\rho_{B}(m+1,n)-\sigma_{\max})}\right],\text{and}\\ &&\\ b_{m,n}&=&\rho_{B}(m,n)/[1+e^{\alpha(\rho_{A}(m-1,n)+\rho_{B}(m-1,n)-\sigma_{\max})}]\end{array}$$

$$\begin{array}[]{lll}\rho_{A}(m,n)&=&\rho_{A}(m,n-1)+a_{m-1,n-1}-a_{m,n-1}\\ &&\\ \rho_{B}(m,n)&=&\rho_{B}(m,n-1)+b_{m+1,n-1}-b_{m,n-1}\end{array}$$

$$\frac{\partial {\rho }_{A(B)}(x,t)}{\partial t}=-(+)C\frac{\partial }{\partial x}\left[\frac{{\rho }_{A(B)}(x,t)}{1+e^{\alpha ({\rho }_A(x,t)+{\rho }_B(x,t)-{\sigma }_{\max })}}\right],$$

$$G(n)=\frac{1}{(L-1)}\left[L+1-2\left(\frac{{\sum }_{j=1}^L(L+1-j){\rho }_{A,B}(j,n)}{{\sum }_{j=1}^N{\rho }_{A,B}(j,n)}\right)\right].$$

$$M(n)=\frac{1}{N_{part}}\sum _{i=1}^{N_{part}}{\xi }_i(n),$$

$$\delta_{x,y}=\left\{\begin{array}[]{ccc}1&\text{se}&x=y\\ &&\\ 0&\text{se}&x\neq y\end{array}\right.$$

$${\alpha }_c\sim {\sigma }_{\max }^{-\Delta },$$

Full text

††journal: Journal of Statistical Mechanics

1. 1 Introduction

2. 2 The model

3. 2.1 Monte Carlo simulations

4. 3 Results

5. 3.1 DDP Initial Conditions

6. 3.2 UD Initial Conditions

7. 3.3 Gini Coefficient

8. 3.4 Occupation Effects and Dilution

9. 3.5 Finite Size Scaling and Effects of the Updating Schemes on the Mobility

10. 4 Summary and conclusions

1 Introduction

The motion of counterflowing streams of particles [1, 2, 3] is a current important topic in the physics of complex systems. Examples range from the motion of oppositely charged colloidal particles[4, 5] to the flow of pedestrians[2, 3, 6]. The formation of patterns is a general property of such systems and different models have been proposed to explain the emergence of structures such as lanes, clogging and jamming. [3, 7, 8].

Three dimensional flows often can be studied by one dimensional transport models [9, 10] or decomposed in separate flows through narrow pipes or independent cells. For instance, we have proposed a model for the flow or counterflow of particles that is governed by a two-species partial differential equation system (TSPDES):[11, 7, 12]

[TABLE]

where $\rho_{A,B}$ is the mass concentration, $K_{A,B}$ is an intra/interspecies interaction constant, $D$ is a diffusion constant, and $f_{A,B}(\rho_{A},\rho_{B})$ is a damping function dependent upon the type of interaction between species $A$, and $B$. It is possible to use these equations to model either the full counterflowing motion of particles or the motion of a single species. For instance, if the damping function $f=1$, $K_{A,B}<0$, and $D_{A,B}>0$, one obtains a standard diffusion equation for a single species without interaction:

[TABLE]

The solution for this equation is given by setting the initial concentration to a known value $\rho_{A,B}(x,t=0)=g_{A,B}(x)$. On the other hand, the modeling of two counterflows is given by working with both equations and setting $K_{A}\cdot K_{B}\leq 0$.

As found in some alternatives derived for the interaction between particles[11, 13], we consider the counterflow of two species in a directed random walk. In this model, the particles do never return, thus $D_{A}=D_{B}=0$. Furthermore, the probability of a particle hopping to the next cell (or remaining in the same cell) is $1/2$ if this cell is empty. Moreover, this probability decreases in a rate $\alpha$ multiplied by the ratio between the occupation of one specie and the total occupation of a cell. Eq. 1 is then rewritten as:

[TABLE]

where $k_{1}\geq 0$, $k_{2}\geq 0$. These coupled differential equations describe interacting counterflowing streams of two oppositely charged species under a field with magnitude proportional to $\alpha$ in a one dimensional discretized ring. The damping factor is given in this case by:

[TABLE]

This problem has been extended for two-dimensional circular crowns[7]. Lanes have been observed in the steady state regime of these systems exactly as observed in pedestrian dynamics [6] or in the motion of charged colloids[4, 5]. Notwithstanding, such models exclude very important effects such as the self exclusion of particles. Thus, we have recently proposed a model[12] where the occupation of the next cells is given by a modified Fermi-Dirac distribution. This leads to the following choice for the damping function:

[TABLE]

where $\sigma_{\max}$ works similarly to the Fermi level as the indicator for the maximum occupation of the cell. This is weighted by the randomness parameter $\alpha$.

Here we explore the impact of the initial conditions on the TSPDES considering the choice of $f_{A,B}(\rho_{A},\rho_{B})$ made in Eq. 5. In order to conduct this study we use the Gini coefficient to investigate the phase transition between a mobile and a clogging phase (or condensate phase) as partially explored earlier[12]. We show that the Gini coefficient is capable of indicating important details of this transition.

We complete this study by performing MC simulations for asynchronous and synchronous dynamics. We compare both alternatives and also compare them with numerical solutions of the TSPDES according to the $f_{A,B}(\rho_{A},\rho_{B})$ given in Eq. 5.

In section 2 we present the Fermi-Dirac directed random walk (FDDEW) model as a derivation of the TSPDES to a problem previously defined [12]. We obtain recurrence relations for this model and extend it to the continuous limit.

Section 3 covers the solutions of TSPDES where we obtain the steady state Gini coefficient as function of $\alpha$. We determine a critical $\alpha_{c}$ coefficient and analyse its dependence upon the initial distribution of particles $A$ and $B$. We also determine the Gini coefficient via Monte Carlo (MC) simulations with both synchronous and asynchronous updating schemes. We finish this work comparing these two strategies with regards to a mobility parameter previously defined[12].

Section 4 concentrates the most important conclusions of this work.

2 The model

This work focuses on a two-specie model of particles that drift in counterflows on an annular system composed of $L$ cells. Considering that the number of particles (of whichever species) in the following cell affects the movement of the particles in the present cell, the concentration of the target objects can be written by the recurrence relation:

[TABLE]

where $p_{j,k}$ is the number (or density) of particles of species $A$ in cell $j$ at time $k$. Moreover, $p_{m,m}^{(n-1)}+p_{m,m+1}^{(n-1)}=1$, since $p_{i,j}^{(n)}$ denotes the probability that a particle in cell $i$ (position $x=i\varepsilon$) hops to cell $j$ (position $x=j\varepsilon$) at $t=$ $n\tau$. In these equations $\tau$ is the time necessary to perform such transition and $\varepsilon$ is the step length.

Combining these equations, one obtains:

[TABLE]

We use a generalized Fermi-Dirac occupation function to model a behavior similar to that found in solids. The peculiarity of our model is that a Fermi level-like parameter indicates the maximum occupation of a cell. Moreover, the inverse temperature relates to a parameter $\alpha$. Thus, the probability of finding a total of $\sigma_{j,n}=\rho_{A}(j,n)+\rho_{B}(j,n)$ particles in cell $j$ at time $n$ is given by:

[TABLE]

A process leading to a total number of particles $\sigma_{j,n}>\sigma_{max}$ in a cell should not be as likely as a process that leads to the opposite result. Also, the likelihood of any of these processes is regulated by $\alpha$, which is not necessarily the inverse temperature, but a control parameter.

Two regimes of this model are of special attention. The first one happens when $\alpha\rightarrow 0$. In this case the occupation probability $p_{i,j}^{(n)}\rightarrow 1/2$, and this is equivalent to a situation where the particles are small enough so that the interparticle collision is completely random.

Another limiting situation happens when $\alpha\rightarrow\infty$. This produces an occupation probability given by:

[TABLE]

This implies that no more than $\sigma_{\max}+1$ objects per cell are allowed. This corresponds to a situation where particles move in narrow tube, guided by a field of constant intensity, but due to their large sizes with respect to the environment, the penetration in the next cells is more restricted and deterministic. As depicted in Fig. 1 such dense systems may lead to the formation of condensates.

There is a critical parameter $\alpha_{C}$ in the intermediate regime $\left(\alpha\in(0,\infty)\right)$ where the system loses mobility causing an abrupt change in the distribution of particles. Setting:

[TABLE]

in Eq. 7, one finds:

[TABLE]

Combining the terms, and taking the continuous limit one then obtains:

[TABLE]

where $C=\lim_{\tau,\varepsilon\rightarrow 0}\frac{\varepsilon}{\tau}$ corresponds to a particular case of Eq. 1 for the damping factor shown in Eq. 5.

2.1 Monte Carlo simulations

It is possible to either solve Eq. 11 directly or via Monte Carlo simulations. For the latter, we consider that a particle of species $A$ at instant $i$ in cell $j$ will occupy the cell $j+1$ with probability $p_{i,j+1}^{(n)}=\frac{1}{1+e^{\alpha(\sigma_{j,n}-\sigma_{\max})}}$.

This simulation can be performed with two different updating schemes:

Synchronous update: All particles are simultaneously updated. This choice was performed in [12]; 2. 2.

Asynchronous update: We select $n_{part}$ particles (the total number of the particles) distributed in $L$ cells. The update is performed for each particle leaving the remaining ones unchanged.

In order to quantify the fluctuations and find the transition between a mobile and a clogged phase we propose the use of the Gini coefficient $G(t)$[15]. Since the number $\rho_{A,B}(j,n)$ of particles from species $A$ and $B$ at cell $j=1,...,L$ at an instant $n$ is known, the Gini coefficient for the concentration of particles is given by:

[TABLE]

Unlike other order parameters such as the mobility[12], the Gini coefficient can be calculated during both Monte Carlo simulations and the direct solution of the TSPDES. Nonetheless, the mobility is given by:

[TABLE]

where $\xi_{i}(n)$ is a binary variable associated to particle $i$ that assumes 0 if the particle stays still at time $n$ and 1 if this same particle hops to the next cell at that period. This quantity cannot be calculated from the solution of the recurrent relations. Nonetheless, it is easily obtained from the MC simulations since they are performed directly on the particles.

It is important to notice that unlike the direct solutions of TSPDES, the steady state values of $G_{\infty}$ and $M_{\infty}$ for MC simulations depend only slightly on the run showing little ensemble variability. Furthermore, the time $t_{steady}$ necessary to reach the steady state regime depends on the value of $\alpha$.

This stochastic behavior was not observed when we directly solved the TSPDES. In this case $10^{5}$ iterations are enough to obtain $G_{\infty}$ and $M_{\infty}$. The slopes of linear fittings for the last $\Delta=10^{3}$ values of $G$ and $M$ were calculated for MC simulations until both of them were smaller than $\eta=10^{-7}$.

In the next section we present the main results obtained from both the direct solution of the TSPDES and MC simulations.

3 Results

We studied the time evolution of particles of species $A$ counterflowing particles of species $B$ under three distinct initial conditions:

Dirac Delta Pulses (DDP): All particles of species $A$ are placed in the cell $i$ and all particles of species $B$ are placed in cell $i+L-1$, where $L$ is the number of cells in the simulation; 2. 2.

Uniform Distribution (UD): The same number of particles of species $A$ and $B$ are uniformly distributed over the $L$ cells; 3. 3.

Constant Occupation (CO): Each cell has two particles: one of species $A$ and another of species $B$.

3.1 DDP Initial Conditions

We simulated a ring with $L=128$ cells by integrating the TSPDES. The DDP initial conditions in this case are explicitly given by $\rho_{A}(x,t=0)=L\delta_{x,0}$ and $\rho_{B}(x,t=0)=L\delta_{x,L}$, with:

[TABLE]

Figures 2, 3, and 4 show the distributions of concentration for $\alpha=0.4$, $3.0$, and $20.0$ respectively for different periods ($t=0,10^{3},10^{4}$ and $10^{5}$).

For $\alpha=0.4$, the steady state regime is reached when $t>10^{4}$. In this situation the particles flow without any clogging effects and are uniformly distributed along the ring, i.e., $\rho_{A}=\rho_{B}=1$ for all cells. This behavior is shown in Fig. 2.

The onset of condensation happens at $\alpha=3$ as depicted in Fig. 3. At this level of disorder it is possible to observe peaks of approximately 6 particles in a cell. Complete clogging is obtained at $\alpha=20$ as shown in Fig. 3. In this case, clogging happens when the two species meet for the first time.

3.2 UD Initial Conditions

The steady state regime shows a peculiar behavior for UD initial conditions. Here we analyzed the distributions for three distinct values of $\alpha$ ($0.45,0.55,$ and $10$). These results are shown in Fig. 5 a, b, and c.

It is possible to observe in Fig. 5a that the system evolves to a completely mobile phase without condensation. As depicted by a pronounced peak in Fig. 5b, condensation appears when $\alpha=0.55$. This characterizes a clogging phase.

The system then moves to a state composed of small condensates for values of $\alpha$ as high as 10. This can be observed in Fig. 5c. Unlike to pronounced peak found for $\alpha=0.55$, here the particles are uniformly distributed.

3.3 Gini Coefficient

The collective phenomena can be further accessed by the Gini coefficient. Since it measures the heterogeneity of the distribution, we use it to quantify the condensation of the system. For instance, a single peak corresponds to a maximum Gini coefficient ($G\rightarrow 1$). A distribution of condensates yields intermediate values ($G\sim 1/2$). Finally, a homogeneous distribution corresponding to the absence of condensation produces small Gini values ($G\rightarrow 0$).

Simulations based on the integration of the TSPDES were performed for $t_{\infty}=10^{4}$ iterations. This is sufficient to reach a steady state value for $G_{\infty}$. On the other hand, it was not possible to achieve convergence for a single $G_{\infty}$ value in MC simulations. We used $G_{\infty}\ \text{or}\ M_{\infty}<10^{-7}$ as a stop condition.

In order to verify whether $G_{\infty}$ is a good estimator for the clogging-mobile transitions, it was measured as a function of $\alpha$ for i) asynchronous MC simulations, ii) synchronous MC simulations, and iii) numerical TSPDES integration.

The Gini coefficient changes abruptly as a function of $\alpha$ as shown in Figs. 6 (a), (b), and (c). Although this transition does not occurr exactly at the same value of $\alpha$, there is a qualitative agreement between the MC simulations and the numerical integration of the TSPDES. Moreover, the transition always occurr for $\alpha<1$.

The Gini coefficient tends to 0 for the numerical integration of the TSPDES and $1/2$ for MC simulations in the mobile phase (low $\alpha$). The latter case corresponds to a situation where there are two types of occupation of the cells, whereas the former corresponds to a homogeneous occupation. Thus, it is possible to state that MC simulations produce statistically homogeneous steady state regimes, whereas the numerical integration of the TSPDES produces completely homogeneous steady state regimes.

Initial uniform distributions do not change in time if $\alpha$ is sufficiently high. Thus, $G_{\infty}\approx 1/2$ for both MC simulations and the numerical integration as depicted in Fig. 6a. Fig. 6b makes explicit the wave behavior of the model. The system departs from two concentrated groups and produces a transition simultaneously for both synchronous MC simulation and TSPDES integration. This happens because of the synchronicity of both schemes. After the two groups collide they produce a permanent condensation that does not dissipate. On the other hand, asynchronous MC simulations have a diffusive character. Hence, particles are dispersed and the condensation is less strong than found in the other cases. Therefore, the Gini coefficient is smaller. Furthermore, both the synchronous and asynchronous MC simulations show a perfect agreement in the limit of low interaction for $\alpha<3$.

The Gini coefficient for a homogeneous initial distribution is shown in Fig. 6. Unlike the behavior described for uniform and concentrated initial conditions, $G_{\infty}\rightarrow 0$ for high values of $\alpha$ when the system departs from a homogeneous occupation. This happens because the system does not change its state when $\alpha$ is too high. A numerical strategy was used to perform the numerical integration in this situation. The center site was initially set to zero particles, whereas the remaining ones were set to one particle: $\rho_{A,B}(x,t=0)=1-\delta_{x,L/2}$.

The integration period $t_{\infty}$ is particularly important for the numerical integration of the TSPDES. Fig. 6d shows the Gini coefficient calculated for $t_{\infty}=1.0\times 10^{3},5.0\times 10^{3},\ \text{and }1.0\times 10^{4}$. Since there is no considerable differences between the Gini coefficients calculated for $t_{\infty}=5.0\times 10^{3},\ \text{and }1.0\times 10^{4}$, one can assume that the latter is a suitable value for the simulations.

3.4 Occupation Effects and Dilution

In order to study dilution and occupation effects, we integrated the TSPDES for a set of $N_{part}=128$ particles of each species uniformly distributed over $L=128$ cells. $G_{\infty}$ as a function of $\alpha$ was monitored while $\sigma_{\max}$ was varied from 1 to 20.

As Fig. 7 shows, the critical transition initially moves to higher values of $\alpha$ as $\sigma$ increases. When $\sigma=4$ this trend ends and the transitions becomes smoother, often showing irregularities with an onset that reduces with an increasing $\sigma$. This is a consequence of the system having more degrees of freedom because of larger capacities of the cells. Thus, the system becomes more mobile. This new trend continues until $\sigma_{\max}=20$. At this stage, the transition becomes continuous, its order changes, and $G_{\infty}$ reduces.

Here we used a constant density of particles of each species: $\overline{\rho}=N_{part}/L=N_{0}$. Another way of carrying this simulation is by keeping the occupation constant. This is defined by $C=N_{part}/(\sigma_{\max}L)=C_{0}$. Thus, $L$ remains constant, but $\sigma_{\max}$ and $L_{part}$ can be varied during the simulation.

Fig. 8 shows the behavior of $G_{\infty}$ as a function of $\alpha$. As the insect indicates, the critical value of $\alpha$ has a power law dependence with $\sigma_{\max}$:

[TABLE]

where $\Delta\approx 0.86$.

3.5 Finite Size Scaling and Effects of the Updating Schemes on the

Mobility

$G_{\infty}$ was measured as a function of $\alpha$ for different values of $L$ using DDP initial conditions. The results shown in Fig. 9 indicate that no significant variations are observed for $L>128$. Therefore, $L=128$ is an appropriate value for our simulations.

Finally, the steady state mobility $M_{\infty}$ for synchronous [12] and asynchronous updating schemes were compared.

Fig. 10a shows that a transition between a clogged (condensation) and a mobile phase can be observed for large systems $L\geq 2^{4}$. As the system is made bigger, the transition becomes more abrupt. Small systems, though, show false mobile phases even for large values of $\alpha$.

This anomaly is more pronounced in the synchronous dynamics, where the system can recover its mobility for very large values of $\alpha$. This is explained by the symmetry of the Fermi-Dirac distribution. Nonetheless, this effect is severely reduced for $L>2^{4}$, which is still a small system. This is further explained elsewhere [12]. On the other hand, this is not observed for the asynchronous dynamics even for extremely small systems.

4 Summary and conclusions

This works extends a previous study on the synchronous dynamics of counterflowing particles.[12] Here we studied its asynchronous Monte Carlo dynamics and expanded the analysis for different initial conditions. Our results indicate that there is a transition from a mobile to a condensation phase. Furthermore we used the Gini coefficient as a non-conventional order parameter.

We began our discussion showing that the problem of counterflowing streams of particles is more general than previously stated. Moreover, we show that the Fermi-Dirac directed random walk used in this study is appropriate to model the clogging-mobile transition.

Furthermore, the level of randomness (or determinism) of the system is determined by the ratio between the size of the particles and the thickness of the tube wherein they flow. Thus, our study suggests that the steady state properties of the system are strongly dependent upon the initial conditions at which the system is prepared.

Acknowledgements

R. da Silva and E. V. Stock were financially supported by CNPq under grant numbers: 311236/2018-9, 424052/2018-0, and 154822/2016-7. This work was partly developed using the resources of Cluster Ada, IF-UFRGS.