Confined Sandpile in two Dimensions: Percolation and Singular Diffusion

TL;DR

This paper studies a two-dimensional confined sandpile model, revealing two scale-invariant regimes related to percolation and giant cluster formation, with implications for understanding critical phenomena in confined granular systems.

Contribution

It derives a diffusion equation from microdynamics and identifies two distinct scale-invariant regimes in a confined sandpile model, linking them to percolation and cluster coalescence phenomena.

Findings

Power-law cluster size distribution near percolation threshold

Evidence of gradient percolation characteristics

Power-law jump size distribution in giant cluster regime

Abstract

We investigate the properties of a two-state sandpile model subjected to a confining potential in two dimensions. From the microdynamical description, we derive a diffusion equation, and find a stationary solution for the case of a parabolic confining potential. By studying the systems at different confining conditions, we observe two scale-invariant regimes. At a given confining potential strength, the cluster size distribution takes the form of a power law. This regime corresponds to the situation in which the density at the center of the system approaches the critical percolation threshold. The analysis of the fractal dimension of the largest cluster frontier provides evidence that this regime is reminiscent of gradient percolation. By increasing further the confining potential, most of the particles coalesce in a giant cluster, and we observe a regime where the jump size distribution takes the form of a power law. The onset of this second regime is signaled by a maximum in the fluctuation of energy.

Equations

$${\rho }_{i,j|n+1}$$

$${\Omega }_{i,j|k}^{\pm }=\prod _{k^{\prime }=1}^k{\rho }_{i\pm k^{\prime },j}\text{and}{\Psi }_{i,j|k}^{\pm }=\prod _{k^{\prime }=1}^k{\rho }_{i,j\pm k^{\prime }},$$

$${\Xi }_{i,j}^{\pm }=\sum _{i^{\prime }=1}^{\infty }{\Theta }_{i,i\pm i^{\prime }|j}{\Omega }_{i,j|i^{\prime }}^{\pm }$$

$${\Lambda }_{i,j}^{\pm }=\sum _{i^{\prime }=1}^{\infty }(1-{\rho }_{i\pm i^{\prime },j}){\Theta }_{i\pm i^{\prime },i|j}{\Omega }_{i,j|i^{\prime }-1}^{\pm }$$

$$\frac{{\rho }_{i,j|n+1}-{\rho }_{i,j|n}}{\tau }$$

$$\frac{\partial \rho }{\partial t}=D\nabla \cdot \left[\frac{(1+\rho )}{(1-\rho {)}^3}\nabla \rho +\beta \frac{(1+\rho )}{(1-\rho {)}^2}\rho \nabla \Phi \right],$$

$$D={\underset{\delta \to 0}{\lim }}_{\tau \to 0}\frac{{\delta }^2}{4\tau },$$

$$\mathbf{J}=-D\left[\frac{(1+\rho )}{(1-\rho {)}^3}\nabla \rho +\beta \frac{(1+\rho )}{(1-\rho {)}^2}\rho \nabla \Phi \right].$$

$$r\left[\frac{(1+\rho )}{(1-\rho {)}^3}\frac{d\rho }{dr}+\beta \frac{(1+\rho )}{(1-\rho {)}^2}\rho \frac{d\Phi }{dr}\right]=0.$$

$$\frac{1}{\beta (1-\rho )\rho }\frac{d{\rho }_{\text{st}}}{dr}=-\frac{d\Phi }{dr}=\frac{d}{dr}\left[{\int }_{u_0}^{{\rho }_{\text{st}}}\frac{1}{\beta (1-u)u}du\right].$$

$${\rho }_{\text{st}}(r)=\frac{1}{1+e^{\beta [\Phi (r)-\mu ]}},$$

$${\int }_0^{\infty }\frac{1}{1+e^{\beta [\Phi (r)-\mu ]}}rdr=\frac{N_p{\delta }^2}{2\pi }.$$

$${\rho }_{\text{st}}(r)=\frac{(e^{\beta \frac{\kappa N_p{\delta }^2}{\pi }}-1)}{e^{\beta \frac{\kappa N_p{\delta }^2}{\pi }}+(e^{\beta \kappa r^2}-1)}.$$

$$S(s)\approx \underset{t\to \infty }{\lim }\left[\frac{{\mathcal{N}}_{[t]}(s)}{{\mathcal{N}}_{[t]T}}\right],$$

$$S(s)\equiv S_o\sum _i\sum _j\frac{{\rho }_{i,j}}{4}\left[{\lambda }_{i,j|s}^{+}+{\lambda }_{i,j|s}^{-}+{\upsilon }_{i,j|s}^{+}+{\upsilon }_{i,j|s}^{-}\right],$$

$${\lambda }_{i,j|s}^{\pm }=(1-{\rho }_{i\pm s,j}){\Theta }_{i\pm s,i|j}{\Omega }_{i,j|s-1}^{\pm }.$$

$$⟨\Delta {\varphi }^2⟩$$

$$\mathcal{A}(ϵ)\propto {ϵ}^{2-d_f}.$$

Full text

I Introduction

Anomalous diffusion is observed in many physical scenarios such as fluid transport in porous media Kuntz2001 ; Lukyanov2012 , diffusion in crowded fluids Szymanski2009 , diffusion in fractal-like substrates Stephenson1995 ; Andrade1997 ; Buldyrev2001 ; Costa2003 ; Havlin2002 , turbulent diffusion in the atmosphere Richardson1926 ; Hentschel1984 , diffusion of proteins due to molecular crowding Banks2005 , systems including ultra-cold atoms Sagi2012 , analysis of heartbeat histograms Grigolini2001 , diffusion in “living polymers” Ott1990 and study of financial transactions Plerou2000 . Anomalous diffusion can also manifest its non-Gaussian behavior in terms of nonlinear Fokker-Plank equations Lenzi2001 ; Malacarne2001 ; Malacarne2002 ; DaSilva2004 ; Lenzi2005 , which is the case, for example, of the dynamics of interacting vortices in disordered superconductors Zapperi2001 ; Moreira2002 ; Miguel2003 ; Andrade2010 , diffusion in dusty plasma Liu2008 ; Barrozo2009 , and pedestrian motion Barrozo2009 . A very interesting case of anomalous diffusion is surely the singular diffusion which is identified as having a divergent diffusion coefficient Carlson1990 ; Carlson1993 ; Carlson1995 ; Kadanoff1992 ; Barbu2010 . This kind of behavior happens in nature in some physical situations, for instance when adsorbates diffuse on a adsorbent surface, its diffusion can be very nonlinear with a diffusion coefficient which depends on the local coverage $\theta$ as, $D\propto{\left|\theta-\theta_{c}\right|^{-\alpha}}$ Myshlyavtsev1995 . Therefore, the study of the basic mechanisms behind surface diffusion is of large importance for understanding technologically important processes like physical adsorption Vidali1991 and catalytic surface reactions Manandhar2003 ; Hofmann2005 ; Schmidtbauer2012 .

A special class of singular diffusion models were intensively studied by Carlson et al. Carlson1990 ; Carlson1993 in two state 1D sandpile models for which they derive diffusion equations with singularities in the diffusion coefficient of the form $D\propto(1+\rho)/(1-\rho)^{3}$, where $\rho$ is the local density. They suggest that some open driven systems present self-organized criticality Bak1987 because in their continuum limit singularities appear in the diffusion constant at a critical point Carlson1990 . Some characteristics of this model change drastically when a confining potential is applied Pires2015 . The jump-size distribution, for instance, starts to exhibit a power-law behavior which suggests a scale-invariant behavior of the system Pires2015 . Scale-invariant behavior in diffusive systems were also observed in gradient percolation diffusion fronts in 2D Sapoval1985 , that have been shown to display fractal diffusion fronts with characteristic dimension similar to the boundary of critical percolation clusters Sapoval1985 .

In this paper we investigate a 2D confined sandpile model. Our model is the extension of the model introduced in Pires2015 to the case of two dimensions. We are able to deduce the continuous limit for the model, which culminate in a diffusion equation with a singular diffusion coefficient. We observe in the confined system the onset of two scale invariant regimes. The first one occurs when the concentration at the origin (the center of the confining region) reaches the critical percolation threshold. At this point we observe scale invariance in the cluster size distribution as well as fractality in the perimeter of the central cluster. At more confined regimes, when the concentration reaches the critical threshold of singular diffusion, there is another signature of scale invariance in the jump size distribution.

II Model Formulation

In the present model, $N_{p}$ particles are placed on a square lattice where we define the occupation $h_{i,j}$ for each site as [math] (if the site is empty) or $1$ (if the site is occupied). As shown in Fig. 1, at each iteration, a particle can move randomly to any of four directions of the square lattice. The particle moves from an occupied site (the source) and get past all occupied sites on the chosen direction until it reaches an empty site (the target) where it may stay with a given probability (see Fig. 1) in a movement that we label a jump. The probability that a jump is accepted is given by the Metropolis factor $\Theta=\min(1,\exp(-\beta[\Phi_{i_{t},j_{t}}-\Phi_{i_{s},j_{s}}]))$, where we define $\Phi_{i,j}$ as the external potential energy of a site. This probability introduces the effect of a confining potential on the particles. Here we present result using only four directions, since, for this case, a continuous limit of the model can be found analytically. However we also performed simulations with a model where particles can move in any arbitrary direction, leading to entirely similar results.

III Continuous limit of the model

It is possible to verify that the model we described is a Markov process, that is, the occupation probability of a site $(i,j)$ at a given step $n+1$, $\rho_{i,j|n+1}$, can be obtained from the occupation probability of all sites on the previous step,

[TABLE]

where $\Theta_{i|j,j^{\prime}}=\min(1,\exp(-\beta[\Phi_{i,j}-\Phi_{i,j^{\prime}}]))$ is the Metropolis factor, and $P_{i|j,j^{\prime}}$ is the probability of finding all sites between $j$ and $j^{\prime}$ in the $i$ column occupied. Similarly, $\Theta_{i,i^{\prime}|j}=\min(1,\exp(-\beta[\Phi_{i,j}-\Phi_{i^{\prime},j}]))$, where $P_{i,i^{\prime}|j}$ is the probability of finding all sites between column $i$ and $i^{\prime}$ on the $j$ line occupied. We define the probabilities of finding $k$ consecutive sites occupied in a given direction as

[TABLE]

where $\Omega^{\pm}_{i,j|k}$ is used for the horizontal direction and $\Psi^{\pm}_{i,j|k}$ is used for the vertical direction. We then define,

[TABLE]

as the contribution for the probability due to particles arriving from the horizontal direction. In a similar fashion, we define $\Gamma^{\pm}_{i,j}$ as the contribution due to particles arriving from vertical direction, and

[TABLE]

as the contribution due to particles leaving in the horizontal direction. Finally, $\Upsilon^{\pm}_{i,j}$ is the contribution due to particles leaving in vertical direction.

Using these definitions, Eq. (1) can be written as

[TABLE]

where $\tau$ is the time unit. The first factor on the right side of Eq. (2) accounts for particles arriving at a given site from each of four directions of the lattice, while the second factor on the left side of Eq. (2) accounts for particles leaving this site in each of the four directions of the lattice.

The continuous limit of Eq. (2) can be obtained similarly to what was done for 1D Pires2015 , resulting in the following non-linear diffusion equation

[TABLE]

where we define

[TABLE]

with $\delta$ been the space unit used for the space continuous limit and $\rho=\rho(x=i\delta,y=j\delta,t=n\tau)$. From Eq. 3, we see that our model obeys the continuity equation $\partial\rho/\partial t=-\nabla\mathbf{J}$, where we have

[TABLE]

In the case where the dependence of the potential is only radial, $\Phi=\Phi(r)$, the conditions for a stationary solution are $\partial\Phi/\partial r=0$, $\rho_{\text{st}}(r\rightarrow\infty)\rightarrow 0$ and $\mathbf{J}=\mathbf{0}$. Thus,

[TABLE]

which can be written as

[TABLE]

The Eq. (5) can be easily solved and results in the stationary solution given by

[TABLE]

where $\mu$ can be obtained from the constraint $\int\frac{\rho(r)}{\delta^{2}}dS=N_{p}$, leading to

[TABLE]

In the particular case of a parabolic confining potential, $\Phi(r=R\delta)=\kappa r^{2}=\kappa(i^{2}+j^{2})\delta^{2}$, this integral can be solved, and it is possible to show that $\mu=(1/\beta)\ln(\exp(\beta\kappa N_{p}\delta^{2}/\pi)-1)$, so that Eq. (6) is now independent of $\mu$

[TABLE]

In Fig. 2, we can see the agreement between predictions from Eq. (7) and results from numerical simulations. Due to the intrinsic exclusion mechanism of the model, the stationary state given by Eq. (6) is a Fermi-Dirac distribution Fernando2010 , and as $\kappa$ increases, the occupation tends to saturate at $\rho_{\text{st}}=1$ near $r=0$. This behavior leads to the formation of a giant cluster of particles near the origin as $\kappa$ increases. The occupation at $r=0$ follows $\rho_{o}=1-\exp(-\beta\kappa N_{p}\delta^{2}/\pi)$, and, as the inset in Fig. 2 shows, this prediction follows closely the results from numerical simulations.

IV Jump Size Distribution and Mean Square Energy Fluctuation

At this point, it is important to define a quantity that can give information about the dynamics of the model, which is the jump size distribution

[TABLE]

where $s$ is the size of a jump, $S(s)$ is the probability of a jump of size $s$ to occur, $\mathcal{N}_{[t]}(s)$ is the number of jumps of size $s$ that appeared in a time interval $t$, and $\mathcal{N}_{[t]T}$ is the total number of jumps in a time interval. Assuming that it is possible to obtain the dynamics of the model from the occupation probability, we can use the following approximation:

[TABLE]

where $S_{o}$ is a normalization constant, and we obtain the occupation probability from the stationary condition, $\rho_{i,j}=\rho_{st}(x=i\delta,y=j\delta)$. The factors $\lambda^{\pm}_{i,j|s}$ account for the chance of a jump to be accepted in the horizontal direction and is defined as

[TABLE]

Similarly, $\upsilon^{\pm}_{i,j|s}$ accounts for the chance of a jump being accepted in the vertical direction and is defined in a similar fashion. We performed the summations of Eq.(8) numerically and compared it with the results from numerical simulations. The results are displayed in Fig. 3, showing very good agreement.

As $\kappa$ increases, the occupation at the center of the system saturates to one (see Fig. 2). Consequently, in very confined systems, jumps that pass through the whole system have larger probability, as indicated by the peak at a large value of $s$ in Fig. 3 when $\kappa=21.7$. The agreement between numerical simulations and the predictions from Eq. (8) shown in Fig. 3 allow us to confidently study larger systems and extract useful information without the need to run time expensive simulations.

Another useful information about the dynamics of the model is the mean square energy fluctuation of the system, which can be determined from the probabilities $\lambda^{+}_{i,j|s}$ and $\upsilon^{+}_{i,j|s}$ as,

[TABLE]

where the sum goes over all jump sizes $s$ and all sites $(i,j)$. Figure 4 shows the $\langle\Delta\phi^{2}\rangle$ as function of $\kappa$ for different number of particles. As can be seen, this function have a maximum in a specific value of $\kappa^{*}$ for each value of $N_{p}$. The shape of the jump size distribution at the condition of maximum energy fluctuation appears to be size invariant, as shown in Fig. 5. The distributions in Fig. 5 collapse when the jump size is scaled by the diameter $\mathcal{D}=2\sqrt{N_{p}/\pi}$ of a dense cluster with all particles. This result is suggestive of a scale invariant regime observed at $\kappa=\kappa^{*}$.

V Cluster Size Distribution and Diffusion Frontier

In this section we investigate some of the geometric aspects of our model. Figure 6 shows the formation of a giant cluster centered at the origin of the confining potential, as the strength $\kappa$ increases. As $\kappa$ grows, there is a clear percolation-like behavior induced where the large cluster starts to grow at the origin. This can be observed quantitatively by investigating the mean largest cluster size $\mathcal{M}$ which can be obtained at each time step. Figure 7 shows how the mean largest cluster size changes when $\kappa$ increases. For convenience, we plot $\mathcal{M}$ as a function of the density at the origin $\rho_{o}$, which increases monotonically with $\kappa$, (see the inset of Fig. 2). As depicted, there is a sudden increase in the size of the largest cluster for $\rho_{o}\approx p_{c}$, where $p_{c}=0.59274621(13)$ is the percolation critical point Newman2000 . Figure 8 suggests that this increasing in the size of the largest cluster takes place due to a percolation-like transition that changes the cluster size distribution of this model. This conjecture may be supported by investigating the cluster size distribution at the condition where the confining potential starts to induce the formation of a larger cluster in the center of the system, $\kappa^{\prime}=\kappa(\rho_{o}=p_{c})\approx 2.82$. Figure 9 shows that the cluster size distribution at this value of $\kappa$ follows a power law.

Another relevant quantity that can be considered is the cluster external perimeter (the so called hull) of the largest cluster. As shown in Fig. 10, this structure is what we call the diffusion frontier and its characteristics are closely related to the gradient percolation diffusion front Sapoval1985 . We define the frontier radius $\mathcal{R}_{f}$ as the mean distance of the diffusion frontier sites to the origin. To study systems with a lager number of particles, we used random samples Knuth1981 to generate a sample system from Eq. 7. Figure 11 shows the probability of finding a diffusion frontier site at a given distance to the origin. The insets in the Fig. 11 shows $\sigma$, the standard deviation of the distance of the diffusion frontier sites to origin, and $\mathcal{R}_{f}$, the diffusion frontier radius as a function of the strength $\kappa$ of the confining potential. One important difference between our model and usual gradient percolation models is the relation between the linear size $\mathcal{R}_{f}$ and $\sigma$, since both are function of $\kappa$, in such way that we cannot fix one and vary the other. This condition limits the values of $\mathcal{R}_{f}$ and $\sigma$ that are close to the condition where finite size effects affect the diffusion frontier. In Fig. 12, we use multiple resolution analysis Peleg1984 ; Peli1990 to estimate the fractal dimension $d_{f}$ of the diffusion frontier. This method consists in considering all points at a distance less than a given $\epsilon$ from the set for which the fractal dimension is being estimated. These points form a new set with the area given by the Richardson law Peli1990

[TABLE]

One can see, at higher scales, that the frontier appears as a one dimensional line, crossing over to a higher fractal dimension $d_{f}$ at smaller scales. As $\kappa$ decreases, approaching the value $\kappa^{\prime}$ where $\rho_{0}\approx p_{c}$, the values of the higher dimension $d_{f}$ grows, but the scaling region decreases. In the inset of Fig. 12, we can see how $d_{f}$ changes with $\sigma/\mathcal{R}_{f}$.

VI Conclusions

We studied a 2D two-state confined sandpile model. From a microscopic dynamics we derived a singular diffusion equation and were able to obtain an analytical stationary solution for the particular case of parabolic potential. This system appears to display two scale-invariant regimes. The first is observed when the concentration at the origin approaches the critical value for percolation. This regime is similar to what is observed for gradient percolation, that is, power laws in the cluster size distribution are observed, as well as a fractal shape for the singular diffusion frontier. The second regime is associated with more intense confinements, when the concentration in the center approaches the maximum value, and a scale-invariant behavior is observed for the jump size distribution. We derived an analytical expression for the jump size distribution and find it to be in good agreement with our numerical solutions. We could, also, find a natural way to define the onset of the scale-invariant regime as the situation where the energy fluctuations are maximum.

VII Acknowledgments

We thank the Brazilian agencies CNPq, CAPES, FUNCAP, and the National Institute of Science and Technology for Complex Systems (INCT-SC) in Brazil for financial support.