# Nonequilibrium Phase Transition in Non-Local and Nonlinear Diffusion Model

\[

## Abstract

We present the results of analytical and numerical studies of a one-dimensional nonlocal and nonlinear diffusion equation describing non-equilibrium processes ranging from aggregation phenomena to cooperation of individuals. We study a dynamical phase transition that is obtained on tuning the initial conditions and demonstrate universality and characterize the critical behavior. The critical state is shown to be reached in a self-organized manner on dynamically evolving the diffusion equation subjected to a mirror symmmetry transformation.

\]

The study of equilibrium phase transitions is now a mature field  and, increasingly, attention is being paid to outstanding problems in nonequilibrium statistical physics . Such problems are often challenging because of inherent non-linearities. There are many examples of non-equilibrium phenomena which are intrinsically non-local such as the growth of thin films in the presence of shadowing  and the sculpting of the drainage basin of river networks due to erosional processes . A striking development in the field of non-equilibrium statistical physics is the development of the paradigm of self-organized criticality , entailing the competition between two dynamical processes leading to a critical state without any fine tuning of parameters.

The principal theme of this letter is the study of a simple one dimensional nonlinear and nonlocal diffusion equation to elucidate the nature of a non-equilibrium phase transition. The equation is essentially one describing biased diffusion  with the magnitude of the bias determined by the instantaneous configuration of the random walkers undergoing diffusion. Physically, the equation describes aggregation, population dynamics, and represents a simple model for cooperative behavior in game theory . Strikingly, the behavior changes from that of a conventional critical point (which requires tuning) to that of self-organized criticality  on considering the time evolution of a transformed equation obtained by a mirror symmetry transformation with $`x`$ replaced by $`x`$.

Our basic equation for $`P(x,t)`$, the probability that the diffusing particle is at position $`x`$ at time $`t`$, is

$$\frac{P}{t}=v(t)\frac{P}{x}+\frac{1}{4}\frac{^2P}{x^2}$$
(1)   
where both the nonlinearity as well as the nonlocality are introduced in the bias velocity $`v`$ defined by

$$v(t)=_0^{\mathrm{}}𝑑xP(x,t)1/2.$$
(2)   
On setting $`v`$ equal to zero in Eq. (1), one recovers the standard unbiased diffusion equation , whereas one obtains simple biased diffusion, when $`v`$ is a constant. In our equation, $`v`$ is a measure of the imbalance between the population of walkers in the right and left and the drift bias promotes further aggregation.

Eq. (1) describes the temporal evolution of the distribution function $`P(x,t)`$ and leads to one of two outcomes in the large time limit. Depending on the initial distribution, one ends up with the bias to the right or to the left winning so that $`P`$ becomes 1 either at $`+\mathrm{}`$ or at $`\mathrm{}`$. Our focus is on the non-equilibrium transition between these limiting behaviors. Note that there is a set of initial conditions (of measure zero) that correspond to the dynamical critical point – we will demonstrate that the critical behavior is universal.

Let us define a new variable $`w(t)=_0^t𝑑\tau v(\tau )`$, and introduce $`y(t)=xw(t)`$ so that Eq. (1) is cast in the form of a standard diffusion equation,

$$\frac{P}{t}=\frac{1}{4}\frac{^2P}{y^2}.$$
(3)   
For simplicity, let us first consider an initial Gaussian distribution of $`P(x,t)`$ centered around $`x_0`$ and with variance $`\sigma _0`$. The solution of Eq. (3) is then given by

$$P(x,t)=N\mathrm{exp}\left\{\frac{[xx_0w(t)]^2}{t+2\sigma _0^2}\right\}.$$
(4)   
where $`N=1/\sqrt{\pi (t+2\sigma _0^2)}`$ is the normalization constant.

With $`\sigma _0=0`$, $`P(x,t)`$ is also the fundamental solution , which will be used to obtain the distribution at time $`t`$, starting from more general initial conditions. Expression (4) is only a formal solution because $`w(t)`$ is itself a function of $`P(x,t)`$. The transition is between two phases corresponding to aggregation on the right or on the left and therefore one would expect that the critical point would correspond to a situation with no bias (i.e. $`x_0=0`$). In this case, the distribution is symmetric with respect to the origin at all times, and there is nothing to choose between left and right, thus an unbiased behaviour ensues. In order to probe the behaviour in the vicinity of this critical point, one could start with an initial distribution with a tiny value of $`x_0`$ (small bias) and watch how the system evolves. Combining Eqs. (1) and (2) and noting that $`\dot{w}=v`$, we find

$$\dot{w}=\frac{1}{2}\text{Erf}\left\{\frac{w(t)+x_0}{\sqrt{t+2\sigma _0^2}}\right\}$$
(5)   
At criticality ($`x_0=0`$), the velocity is zero at all times and $`w`$ vanishes too. In the critical regime, one expects that $`w(t)1`$, for all $`t`$ less than a crossover time $`T_r`$, so that one may linearize Eq. (5) to find that

$$\dot{w}=\frac{1}{\sqrt{\pi }}\frac{w+x_0}{\sqrt{t+2\sigma _0^2}}.$$
(6)   
This equation can be easily solved and yields

$$w(t)=x_0\left\{\text{e}^{2(\sqrt{t+2\sigma _0^2}\sqrt{2\sigma _0^2})/\sqrt{\pi }}1\right\}.$$
(7)   
One may define a characteristic transient time, $`T_r`$, spent in the critical region during which the linearization approximation holds. From (7), one finds that $`T_r`$ diverges as $`x_00`$ as

$$T_r\mathrm{ln}^2|x_0|.$$
(8)   
In the critical region, the characteristic length-scale $`\xi `$ is expected to follow the diffusion law, and therefore

$$\xi \sqrt{T_r}|\mathrm{ln}|x_0||,$$
(9)   
as we will verify numerically later on.

Alternatively, in the initial condition, one may introduce a bias by fixing $`x_0=0`$ and instead letting $`P(\mathrm{},0)=\mathrm{\Phi }_0`$ as an effective asymmetric boundary condition. With the same procedure used to derive Eq. (5) we have

$$\dot{w}=\frac{\mathrm{\Phi }_0}{2}+\frac{1\mathrm{\Phi }_0}{2}\text{Erf}\left\{\frac{w(t)}{\sqrt{t+2\sigma _0^2}}\right\}.$$
(10)   
$`\mathrm{\Phi }_0=0`$ corresponds to the critical point, and by linearizing Eq. (10) for small $`\mathrm{\Phi }_0`$, it is straightforward to show that one again obtains $`T_r\mathrm{ln}^2(\mathrm{\Phi }_0)`$, with $`\xi \sqrt{T_r}`$.

We have verified that the same critical behaviour holds for other families of initial conditions. For example, when

$$P(x,0)=(1/2\epsilon )\delta (x1)+(1/2+\epsilon )\delta (x+1),$$
(11)   
we find that

$$T_r\mathrm{ln}^2|\epsilon |$$
(12)   
and $`\xi |\mathrm{ln}|\epsilon ||`$. Note that $`\epsilon `$ is now a measure of the deviation from the critical point and $`\epsilon =0`$ correponds to the unbiased, zero drift situation.

We now turn to the results of numerical experiments on a 1-dimensional lattice which are useful for validating our analytic predictions and for probing the nonlinear regime. The discrete version of Eq. (1), used in our simulations, reads

$$P_x(t+1)=\frac{P_x(t)}{2}+\frac{\mathrm{\Phi }(t)}{2}P_{x1}(t)+\frac{1\mathrm{\Phi }(t)}{2}P_{x+1}(t)$$
(13)   
with $`\mathrm{\Phi }(t)=_{x0}P_x(t)`$. The velocity is thus given by $`v(t)=\mathrm{\Phi }(t)1/2`$. This equation was proposed by Nowak and Sigmund  as a simplified model of the evolution of indirect reciprocity by image scoring. Indirect reciprocity is determined by reputation and status and is characterized by each individual having an image score. A potential donor and recipient of an altruistic act have an interaction in which the donor helps the recipient provided the recipient’s image score is positive. Such an altruistic act increases the image score of the donor by $`1`$ (the selfish act would have decreased it by $`1`$) and the image score of the recipients is unchanged. Eq. (13) is the governing equation for the time evolution of $`P_x`$, the fraction of players with image score $`x`$. The two phases that we have considered correspond to cooperation and defection and our finding is that not much time is spent agonizing over which phase to select even in the vicinity of the critical point. Indeed, the time scale to decide on one of the two different phases only diverges logarithmically as one approaches the critical point.

Fig. 1 shows the divergence of $`T_r`$ as $`\epsilon 0`$ for the initial distributions given in Eq. (11) and for a lattice of $`6`$x$`10^5`$ sites. Numerical results are in excellent accord with our theoretical predictions (see Eq. (12))

Let us consider now, an interesting generalization of the initial conditions that we studied analytically: $`P(\mathrm{},0)=\mathrm{\Phi }_0`$ and $`P(z,0)=1\mathrm{\Phi }_0`$ and $`P=0`$ at all other locations, initially. The critical value of $`\mathrm{\Phi }_0`$ increases monotonically to a non-zero value $`\mathrm{\Phi }_c(z)`$ for positive $`z`$ with $`\mathrm{\Phi }_c(z\mathrm{})=1/2`$. The smallest value of $`\mathrm{\Phi }_c(z)`$ on a discrete lattice occurs when $`z=1`$, which is the case we focus on. $`\mathrm{\Phi }_c\mathrm{\Phi }_c(z=1)`$ is found to be 0.261970531164…, a result that was noted earlier by Nowak and Sigmund .

Fig. 2 shows the behaviour of the bias velocity as a function of time as $`\mathrm{\Phi }_0`$ approaches $`\mathrm{\Phi }_c`$ from above and from below.

$`T_r`$ can be identified as the time after which $`v(t)`$ becomes equal to its asymptotic values of either $`1/2`$ and $`\mathrm{\Phi }_01/2`$ and its scaling behaviour is shown in Fig. 3.

The average location of the random walkers (excluding the number fixed at $`x=\mathrm{}`$) behaves with time (in the vicinity of the fixed point) as

$$x(t)\sqrt{t+2\sigma _0^2}.$$
(14)   
In order to derive this result, we note that the velocity increases very slowly in the linear critical regime and can be approximately considered constant (this is consistent with the assumption that, for time $`tT_r`$, $`w`$ is much smaller than $`1`$). The derivative of the velocity given by Eq. (2) is then vanishing, and using Eq. (1) to eliminate $`\dot{P}`$ we obtain, after integrating over $`x`$, the expression

$$4vP(0,t)_xP(0,t)=0.$$
(15)   
Using the formal solution (4), we obtain

$$w(t)\sqrt{t+2\sigma _0^2}$$

The result (14) then follows on noting that the average position of the walkers is given by $`w(t)+x_0`$.

This critical regime behaviour of $`x(t)`$ crosses over to a linear temporal behavior when the bias reaches a sufficient strength. There is indeed a sharp onset of the linear behaviour at a value of $`x`$, which one may identify with $`\xi `$. The scaling behaviour of $`\xi `$ is shown in Fig. 4.

We now turn to a simple mechanism for obtaining self-organized critical behavior in our model. In Eq. (1) the transformation $`xx`$ is equivalent to a change of the sign of the bias velocity. In this situation, the system spontaneously organizes in such a way that the aggregation of walkers is disfavoured. As a consequence, the asymptotic distribution becomes symmetric (characterized by $`v=0`$) and this corresponds to a self-tuning to the critical state, a behaviour typical of self-organized criticality .

In this asymptotic regime, the scaling $`x(t)\sqrt{t}`$ derived in (14) still ought to hold and is confirmed by simulations performed on the discretized diffusion equation (see Fig. 5).

In this letter, we have introduced and studied a diffusion equation with nonlinear and nonlocal features. In this model, spontaneous fluctuations in the population of walkers are able to drive the entire population towards one of the two boundaries located at $`(x=\pm \mathrm{})`$. This mechanism is of interest as the basis for the development of more realistic models of self-aggregation and self-organization in cooperative states of populations of interacting individuals. A mirror symmetry transformation applied to the equation reveals a dynamical evolution corresponding to generic behavior associated with self-organized criticality.

This work was supported by INFN, NASA, NATO and The Donors of the Petroleum Research Fund administered by the American Chemical Society.