Cellular Automata Explorations

Driven Diffusion of Two Species

by Richard J. Gaylord and Kazume Nishidate

The random motion of the particles in a two-species system in the presense of an external field, originally proposed as a model of certain kinds of materials known as superionic conductors, fast ion conductors, super ionic solids, or solid electrolytes, shows rather interesting steady-state, non-equilibrium behavior. In this "driven diffusive" system, the external field biases the movement of each species in opposite directions, and an order-disorder transition occurs as a function of density. We develop here  a CA model of  the driven diffusion of two species.

The driven diffusive, two-species CA

The system

The CA employs an 'n by n' square lattice with periodic boundary conditions. Lattice sites are either empty or occupied by a single particle that has either a positive or negative charge. The total density of occupied sites is p, and there are (approximately) the same number of positive and negative particles, so that total charge of the system is (approximately) zero.  Initially, the particles are randomly distributed on the lattice and each particle faces a randomly chosen direction. The system is created using the following code:

initPhase = 
   Table[{{1, 2, 3, 4}[[Random[Integer, {1, 4}] ]], {A, B}[[Random[Integer, {1, 2}]]] }
             * Floor[Random[] + p], 
             {n}, {n}];

The sites in initPhase have values which are ordered pairs, {x, y}. The first component of each ordered pair is an integer ranging from 0 to 4 where the values have the following meanings:

0  - an empty site
1  - a site occupied by a north-facing particle 
2  - a site occupied by a east-facing particle
3  - a site occupied by a south-facing particle
4  - a site occupied by a west-facing particle

The second component of each ordered pair is either  0, A or B where 0 indicates an empty site, A indicates a site occupied by an A type particle, and B indicates a site occupied by an B type particle.

The CA rules

The rules for the behavior of a particle during a time step, take thirteen arguments

driven[site,  N,  E,  S,  W,  NE,  SE,  SW,  NW,  Nn,  Ee,  
      Ss,  Ww]

The arguments represent the value of the site and the values of 12 neighbor sites : the nearest neighbors in the N, E, S, W directions, the nearest neighbors in the NE, SE, SW, NW directions, and the next nearest neighbors in the N, E, S, W directions (these values are needed to implement the 'excluded volume' constraint).

The presense of an external electrical field causes an A (B) particle to hop preferentially along (against) the field direction. The field points diagonally in the north-east direction which drives the A (B) particles upward (downward). The probability of a step is a in the preferred direction (north and east for an A particle, south and west for a B particle) and (1 - a) in the opposite direction.  

In the CA, a particle will always moves in the direction it is facing, and therefore we will implement the bias in taking a step in various directions by biasing the probabilities that a particle chooses to face those directions. An A particle will choose a direction to face using rndA, and a B particle will choose a direction to face using rndB where rndA and rndB are defined as follows:

rndA  := 
       {Random[Integer,  {1,  2}],  Random[Integer,  {3,  4}]}[
        [2  -  Floor[Random[]  +  a]]];

rndA chooses Random[Integer, {1, 2}] a% of the time and Random[Integer, {3, 4}]
(1-a)% of the time, so that  1, 2, 3, and 4 are chosen (a/2)%, (a/2)%, (1 - a)/2%, and 
(1 - a)/2% of the time respectively:

rndB := {Random[Integer,  {1,  2}],  Random[Integer,  {3,  4}]}[[1 + Floor[Random[]  +  a]]];

rndB chooses Random[Integer, {1, 2}] (1 - a)% of the time and Random[Integer, {3, 4}] a% of the time so that  1, 2, 3, and 4 are chosen (1 - a)/2%, (1 - a)/2%, (a/2)%, and (a/2)% of the time respectively.

At each time step, all of the particles simultaneously attempt to jump to adjacent sites. When a particle attempts to move to the nearest neighbor site it faces, there are three possible situations that will determine the outcome:  

(a) the adjacent site is occupied by another particle, in which case the particle remains in place and randomly chooses a direction to face.
(b) the adjacent site is empty and is faced by one or more other particles, in which case the particle remains in place and randomly chooses a direction to face.
(c) the adjacent site is empty and no other particles face it, in which case the particle abandons its current site and moves into the site and randomly chooses a direction to face.

Regardless of whether or not a particle is able to move, it picks randomly picks a direction to face using rndA if it is an A particle and using rndB if it is a B particle.

The following rules describe the behavior of lattice sites during a time step:

- A particle facing a site occupied by another particle (of either type) remains in place and randomly chooses a direction to face.

driven[{1,  A},  {x_?Positive,  _},  _,  _,  _,  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rndA,  A}driven[{2,  A},  _,  {x_?Positive,  _},  _,  _,  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rndA,  A}driven[{3,  A},  _,  _,  {x_?Positive,  _},  _,  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rndA,  A}driven[{4,  A},  _,  _,  _,  {x_?Positive,  _},  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rndA,  A}driven[{1,  B},  {x_?Positive,  _},  _,  _,  _,  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rndB,  B}driven[{2,  B},  _,  {x_?Positive,  _},  _,  _,  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rndB,  B}driven[{3,  B},  _,  _,  {x_?Positive,  _},  _,  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rndB,  B}driven[{4,  B},  _,  _,  _,  {x_?Positive,  _},  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rndB,  B}

- A particle facing an empty adjacent site that is also faced by another particle of either type (eg., an right-facing particle that lies to the left of an empty site that is faced by an left-facing particle lying to the right of the empty site) stays put and randomly chooses a direction to face.

driven[{1,  A},  {0,  0},  _,  _,  _,  {4,  _},  _,  _,  _, 
         _,  _,  _,  _]  :=  {rndA,  A}driven[{1,  A},  {0,  0},  _,  _,  _,  _,  _,  _,  {2,  _}, 
         _,  _,  _,  _]  :=  {rndA,  A}driven[{1,  A},  {0,  0},  _,  _, _ , _ , _ , _ ,  _, 
         {3,  _},  _,  _,  _]  :=  {rndA,  A}driven[{2,  A},  _,  {0,  0},  _,  _,  {3,  _},  _,  _,  _, 
         _,  _,  _,  _]  :=  {rndA,  A}driven[{2,  A},  _,  {0,  0},  _,  _,  _,  {1,  _},  _,  _, 
         _,  _,  _,  _]  :=  {rndA,  A}driven[{2,  A},  _,  {0,  0},  _,  _,  _,  _,  _,  _,  _, 
         {4,  _},  _,  _]  :=  {rndA,  A}driven[{3,  A},  _,  _,  {0,  0},  _,  _,  {4,  _},  _,  _, 
         _,  _,  _,  _]  :=  {rndA,  A}driven[{3,  A},  _,  _,  {0,  0},  _,  _,  _,  {2,  _},  _, 
         _,  _,  _,  _]  :=  {rndA,  A}driven[{3,  A},  _,  _,  {0,  0},  _,  _,  _,  _,  _,  _, 
         _,  {1,  _},  _]  :=  {rndA,  A}driven[{4,  A},  _,  _,  _,  {0,  0},  _,  _,  {1,  _},  _, 
         _,  _,  _,  _]  :=  {rndA,  A}driven[{4,  A},  _,  _,  _,  {0,  0},  _,  _,  _,  {3,  _}, 
         _,  _,  _,  _]  :=  {rndA,  A}driven[{4,  A},  _,  _,  _,  {0,  0},  _,  _,  _,  _,  _, 
         _,  _,  {2,  _}]  :=  {rndA,  A}driven[{1,  B},  {0,  0},  _,  _,  _,  {4,  _},  _,  _,  _, 
         _,  _,  _,  _]  :=  {rndB,  B}driven[{1,  B},  {0,  0},  _,  _,  _,  _,  _,  _,  {2,  _}, 
         _,  _,  _,  _]  :=  {rndB,  B}driven[{1,  B},  {0,  0},  _,  _, _ , _ , _ , _ ,  _, 
         {3,  _},  _,  _,  _]  :=  {rndB,  B}driven[{2,  B},  _,  {0,  0},  _,  _,  {3,  _},  _,  _,  _, 
         _,  _,  _,  _]  :=  {rndB,  B}driven[{2,  B},  _,  {0,  0},  _,  _,  _,  {1,  _},  _,  _, 
         _,  _,  _,  _]  :=  {rndB,  B}driven[{2,  B},  _,  {0,  0},  _,  _,  _,  _,  _,  _,  _, 
         {4,  _},  _,  _]  :=  {rndB,  B}driven[{3,  B},  _,  _,  {0,  0},  _,  _,  {4,  _},  _,  _, 
         _,  _,  _,  _]  :=  {rndB,  B}driven[{3,  B},  _,  _,  {0,  0},  _,  _,  _,  {2,  _},  _, 
         _,  _,  _,  _]  :=  {rndB,  B}driven[{3,  B},  _,  _,  {0,  0},  _,  _,  _,  _,  _,  _, 
         _,  {1,  _},  _]  :=  {rndB,  B}driven[{4,  B},  _,  _,  _,  {0,  0},  _,  _,  {1,  _},  _, 
         _,  _,  _,  _]  :=  {rndB,  B}driven[{4,  B},  _,  _,  _,  {0,  0},  _,  _,  _,  {3,  _}, 
         _,  _,  _,  _]  :=  {rndB,  B}driven[{4,  B},  _,  _,  _,  {0,  0},  _,  _,  _,  _,  _, 
         _,  _,  {2,  _}]  :=  {rndB,  B}

- An site occupied by a particle facing an empty adjacent site that is not faced by another particle  is vacated.

driven[{1,  _},  {0,  0},  _,  _,  _,  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {0,  0}driven[{2,  _},  _,  {0,  0},  _,  _,  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {0,  0}driven[{3,  _},  _,  _,  {0,  0},  _,  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {0,  0}driven[{4,  _},  _,  _,  _,  {0,  0},  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {0,  0}

- An empty site faced by two or more particles remains empty.

driven[{0,  0},  {3,  _},  {4,  _},  _,  _,  _,  _,  _,  _, 
         _,  _,  _,  _]  :=  {0,  0}driven[{0,  0},  {3,  _},  _,  {1,  _},  _,  _,  _,  _,  _, 
         _,  _,  _,  _]  :=  {0,  0}driven[{0,  0},  {3,  _},  _,  _,  {2,  _},  _,  _,  _,  _, 
         _,  _,  _,  _]  :=  {0,  0}driven[{0,  0},  _,  {4,  _},  {1,  _},  _,  _,  _,  _,  _, 
         _,  _,  _,  _]  :=  {0,  0}driven[{0,  0},  _,  {4,  _},  _,  {2,  _},  _,  _,  _,  _, 
         _,  _,  _,  _]  :=  {0,  0}driven[{0,  0},  _,  _,  {1,  _},  {2,  _},  _,  _,  _,  _, 
         _,  _,  _,  _]  :=  {0,  0}

- An empty site faced by exactly one particle beomes occupied.

driven[{0,  0},  {3,  A},  _,  _,  _,  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {rndA,  A}driven[{0,  0},  _,  {4,  A},  _,  _,  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {rndA,  A}driven[{0,  0},  _,  _,  {1,  A},  _,  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {rndA,  A}driven[{0,  0},  _,  _,  _,  {2,  A},  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {rndA,  A}driven[{0,  0},  {3,  B},  _,  _,  _,  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {rndB,  B}driven[{0,  0},  _,  {4,  B},  _,  _,  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {rndB,  B}driven[{0,  0},  _,  _,  {1,  B},  _,  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {rndB,  B}driven[{0,  0},  _,  _,  _,  {2,  B},  _,  _,  _,  _,  _, 
         _,  _,  _]  :=  {rndB,  B}

- An empty site remains unchanged.

driven[{0,  0},  _,  _,  _,  _,  _,  _,  _,  _,  _,  _,  _, 
         _]  :=  {0,  0}

Note: The number of rules (51) given above is rather large and can be shortened to 31  (with some loss of readibilty and a small  decrease in efficiency) by making the following changes:

(a) using 1 and 2 as the second component of the ordered pair value of a lattice site, rather than A and B, to indicate the type of particle. 

(b) replacing rndA and rndB with

rnd := {{Random[Integer, {1, 2}], Random[Integer, {3, 4}]}[[2 - Floor[Random[] + a]]],
             {Random[Integer, {1, 2}], Random[Integer, {3, 4}]}[[1 + Floor[Random[] + a]]]}

(c) combining the rules that are the same except for the type of particle they apply to. For example, the two rules

driven[{1,  A},  {x_?Positive,  _},  _,  _,  _,  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rndA,  A}driven[{1,  B},  {x_?Positive,  _},  _,  _,  _,  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rndB,  B}

can be merged into the single rule

driven[{1,  y_},  {x_?Positive,  _},  _,  _,  _,  _,  _,  _, 
         _,  _,  _,  _,  _]  :=  {rnd[[y]],  y}

Applying the CA rules

The update rules are applied to the CA lattice using the anonymous function

MvonN[driven,  #]&

where

\(MvonN[func__, \ lat_]\  := \ \ \ 
\ \ 
      MapThread[func, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        Map[RotateRight[lat, \ #]&, 
          \ {{0, \ 0}, \ {1, \ 0}, \ {0, \ \(-1\)}, \ {\(-1\), \ 0}, 
            \ {0, \ 1}, \ 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {1, \ \ 
              \(-1\)}, \ {\(-1\), \ \(-1\)}, \ {\(-1\), \ 1}, \ {1, \ 1}, \ 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {2, \ 0}, 
            \ {0, \ \(-2\)}, \ {\(-2\), \ 0}, \ {0, \ 2}}

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ], \ 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2]\)

The positions of the random walkers evolve over t time steps, starting with the initial lattice configuration, using the following NestList operation.

NestList[MvonN[driven,  #]&, initPhase, t]

The program

DrivenDiffusion2Species[n_,  p_,  a_,  t_]:=
   Module[{initPhase,  driven,  rndA,  rndB},
   
initPhase = 
   Table[{{1,  2,  3,  4}[
                            [Random[Integer,  {1,  4}] ]], {A,  B}[[Random[Integer,  {1,  2}]]] }
             * Floor[Random[]  +  p], 
             {n}, {n}]; 

rndA := {Random[Integer,  {1,  2}],  Random[Integer,  {3,  4}]}
                  [[2 - Floor[Random[]  +  a]]];
      
rndB := {Random[Integer,  {1,  2}],  Random[Integer,  {3,  4}]}
                  [[1 + Floor[Random[]  +  a]]];
      
driven[{1,  A},  {x_?Positive,  _},  _,  _,  _,  _,  _, 
                 _,  _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{2,  A},  _,  {x_?Positive,  _},  _,  _,  _, 
                 _,  _,  _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{3,  A},  _,  _,  {x_?Positive,  _},  _, 
                 _,  _,  _,  _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{4,  A},  _,  _,  _,  {x_?Positive,  _},  _,  _, 
                 _,  _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{1,  B},  {x_?Positive,  _},  _,  _,  _,  _, 
                 _,  _,  _,  _,  _,  _,  _]  :=  {rndB,  B};
driven[{2,  B},  _,  {x_?Positive,  _},  _,  _, 
                 _,  _,  _,  _,  _,  _,  _,  _]  :=  {rndB,  B};
driven[{3,  B},  _,  _,  {x_?Positive,  _},  _,  _,  _, 
                 _,  _,  _,  _,  _,  _]  :=  {rndB,  B};
driven[{4,  B},  _,  _,  _,  {x_?Positive,  _},  _, 
                 _,  _,  _,  _,  _,  _,  _]  :=  {rndB,  B};

driven[{1,  A},  {0,  0},  _,  _,  _,  {4,  _},  _, 
                 _,  _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{1,  A},  {0,  0},  _,  _,  _,  _,  _,  _, 
                 {2,  _},  _,  _,  _,  _]  :=  {rndA,  A};
driven[{1,  A},  {0,  0},  _,  _, _ , _ , _ , _
                 ,  _,  {3,  _},  _,  _,  _]  :=  {rndA,  A};
driven[{2,  A},  _,  {0,  0},  _,  _,  {3,  _}, 
                 _,  _,  _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{2,  A},  _,  {0,  0},  _,  _,  _, 
                 {1,  _},  _,  _,  _,  _,  _,  _]  := 
               {rndA,  A};
driven[{2,  A},  _,  {0,  0},  _,  _,  _,  _,  _,  _, 
                 _,  {4,  _},  _,  _]  :=  {rndA,  A};
driven[{3,  A},  _,  _,  {0,  0},  _,  _,  {4,  _}, 
                 _,  _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{3,  A},  _,  _,  {0,  0},  _,  _,  _, 
                 {2,  _},  _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{3,  A},  _,  _,  {0,  0},  _,  _,  _,  _,  _, 
                 _,  _,  {1,  _},  _]  :=  {rndA,  A};
driven[{4,  A},  _,  _,  _,  {0,  0},  _,  _, 
                 {1,  _},  _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{4,  A},  _,  _,  _,  {0,  0},  _,  _,  _, 
                 {3,  _},  _,  _,  _,  _]  :=  {rndA,  A};
driven[{4,  A},  _,  _,  _,  {0,  0},  _,  _, 
                 _,  _,  _,  _,  _,  {2,  _}]  :=  {rndA,  A};
driven[{1,  B},  {0,  0},  _,  _,  _,  {4,  _},  _, 
                 _,  _,  _,  _,  _,  _]  :=  {rndB,  B};
driven[{1,  B},  {0,  0},  _,  _,  _,  _,  _,  _, 
                 {2,  _},  _,  _,  _,  _]  :=  {rndB,  B};
driven[{1,  B},  {0,  0},  _,  _, _ , _ , _ , _
                 ,  _,  {3,  _},  _,  _,  _]  :=  {rndB,  B};
driven[{2,  B},  _,  {0,  0},  _,  _,  {3,  _}, 
                 _,  _,  _,  _,  _,  _,  _]  :=  {rndB,  B};
driven[{2,  B},  _,  {0,  0},  _,  _,  _, 
                 {1,  _},  _,  _,  _,  _,  _,  _]  := 
               {rndB,  B};
driven[{2,  B},  _,  {0,  0},  _,  _,  _,  _,  _,  _, 
                 _,  {4,  _},  _,  _]  :=  {rndB,  B};
driven[{3,  B},  _,  _,  {0,  0},  _,  _,  {4,  _}, 
                 _,  _,  _,  _,  _,  _]  :=  {rndB,  B};
driven[{3,  B},  _,  _,  {0,  0},  _,  _,  _, 
                 {2,  _},  _,  _,  _,  _,  _]  :=  {rndB,  B};
driven[{3,  B},  _,  _,  {0,  0},  _,  _,  _,  _,  _, 
                 _,  _,  {1,  _},  _]  :=  {rndB,  B};
driven[{4,  B},  _,  _,  _,  {0,  0},  _,  _, 
                 {1,  _},  _,  _,  _,  _,  _]  :=  {rndB,  B};
driven[{4,  B},  _,  _,  _,  {0,  0},  _,  _,  _, 
                 {3,  _},  _,  _,  _,  _]  :=  {rndB,  B};
driven[{4,  B},  _,  _,  _,  {0,  0},  _,  _, 
                 _,  _,  _,  _,  _,  {2,  _}]  :=  {rndB,  B};

driven[{1,  _},  {0,  0},  _,  _,  _,  _,  _,  _,  _, 
                 _,  _,  _,  _]  :=  {0,  0};
driven[{2,  _},  _,  {0,  0},  _,  _,  _,  _,  _,  _, 
                 _,  _,  _,  _]  :=  {0,  0};
driven[{3,  _},  _,  _,  {0,  0},  _,  _,  _,  _,  _, 
                 _,  _,  _,  _]  :=  {0,  0};
driven[{4,  _},  _,  _,  _,  {0,  0},  _,  _,  _,  _, 
                 _,  _,  _,  _]  :=  {0,  0};

driven[{0,  0},  {3,  _},  {4,  _},  _,  _,  _,  _, 
                 _,  _,  _,  _,  _,  _]  :=  {0,  0};
driven[{0,  0},  {3,  _},  _,  {1,  _},  _,  _,  _, 
                 _,  _,  _,  _,  _,  _]  :=  {0,  0};
driven[{0,  0},  {3,  _},  _,  _,  {2,  _},  _,  _, 
                 _,  _,  _,  _,  _,  _]  :=  {0,  0};
driven[{0,  0},  _,  {4,  _},  {1,  _},  _,  _,  _, 
                 _,  _,  _,  _,  _,  _]  :=  {0,  0};
driven[{0,  0},  _,  {4,  _},  _,  {2,  _},  _,  _, 
                 _,  _,  _,  _,  _,  _]  :=  {0,  0};
driven[{0,  0},  _,  _,  {1,  _},  {2,  _},  _,  _, 
                 _,  _,  _,  _,  _,  _]  :=  {0,  0};

driven[{0,  0},  {3,  A},  _,  _,  _,  _,  _, 
                 _,  _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{0,  0},  _,  {4,  A},  _,  _,  _,  _,  _, 
                 _,  _,  _,  _,  _]  :=  {rndA,  A};
driven[{0,  0},  _,  _,  {1,  A},  _,  _,  _,  _,  _, 
                 _,  _,  _,  _]  :=  {rndA,  A};
driven[{0,  0},  _,  _,  _,  {2,  A},  _,  _,  _,  _, 
                 _,  _,  _,  _]  :=  {rndA,  A};
driven[{0,  0},  {3,  B},  _,  _,  _,  _,  _,  _,  _, 
                 _,  _,  _,  _]  :=  {rndB,  B};
driven[{0,  0},  _,  {4,  B},  _,  _,  _,  _,  _,  _, 
                 _,  _,  _,  _]  :=  {rndB,  B};
driven[{0,  0},  _,  _,  {1,  B},  _,  _,  _,  _,  _, 
                 _,  _,  _,  _]  :=  {rndB,  B};
driven[{0,  0},  _,  _,  _,  {2,  B},  _,  _,  _,  _, 
                 _,  _,  _,  _]  :=  {rndB,  B};

driven[{0,  0},  _,  _,  _,  _,  _,  _,  _,  _,  _, 
                 _,  _,  _]  :=  {0,  0};

\(MvonN[func__, \ lat_]\  := \ \ \ 
\ \ 
              MapThread[func, \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                Map[RotateRight[lat, \ #]&, 
                  \ {{0, \ 0}, \ {1, \ 0}, \ {0, \ \(-1\)}, \ {\(-1\), \ 0}, 
                    \ {0, \ 1}, \ 

                    \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {1, 
                      \ \ \(-1\)}, \ {\(-1\), \ \(-1\)}, \ {\(-1\), \ 1}, 
                    \ {1, \ 1}, \ 

                    \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {2, 
                      \ 0}, \ {0, \ \(-2\)}, \ {\(-2\), \ 0}, \ {0, \ 2}}], \ 
                2]\);

evolve = NestList[MvonN[driven,  #]&, initPhase, t]
]

Running the program

The lattice configuration at each time step, showing whether a site is empty or is occupied by an A particle or a B particle can be determined by extracting the second component of the ordered pair value of each lattice site, using

DrivenDiffusion2Species[n, p, a, t] /. {{_, 0} -> 0, {_, A} -> A, {_, B} -> B}

For example, we can look at the development over 100 time steps, of a 30 by 30 lattice system with particles having a 90% probability of moving in the preferred direction and a 10% probability of moving in the opposite direction. We can compare the low density (p = 0.1) and high density (p = 0.5) cases.

lowDensityCase  =  
      DrivenDiffusion2Species[30,  0.1,  0.9,  100]  /. 
         {{_,  0}  ->  0,  {_,  A}  ->  A,  {_,  B}  ->  B};

Show[GraphicsArray[Partition[Map[Show[Graphics[RasterArray[\(#\  /. \ 

                                  \ \ \ \ \ \ \ \ \ \ {
                                    0\  -> \ RGBColor[0.7, \ 0.7, \ 0.7], \ 
                                    A\  -> \ RGBColor[0, \ 1, \ 0], 

                                    \ \ \ \ \ \ \ \ \ \ \ 
                                    B\  -> \ RGBColor[1, \ 0, \ 1]}\)]], 
           AspectRatio  ->  Automatic, DisplayFunction  ->  Identity 
                 ]&, 
          lowDensityCase[[Range[1,  101,  20] ]] 
   ], 2] ]];

highDensityCase = DrivenDiffusion2Species[30, 0.5, 0.9, 100] /. {{_, 0} -> 0, {_, A} -> A, {_, B} -> B};

Show[GraphicsArray[Partition[Map[Show[Graphics[RasterArray[\(#\  /. \ 

                                \ \ \ \ \ \ \ \ \ \ {
                                  0\  -> \ RGBColor[0.7, \ 0.7, \ 0.7], \ 
                                  A\  -> \ RGBColor[0, \ 1, \ 0], 

                                  \ \ \ \ \ \ \ \ \ \ \ 
                                  B\  -> \ RGBColor[1, \ 0, \ 1]}\)]], 
           AspectRatio  ->  Automatic, DisplayFunction  ->  Identity 
                 ]&, 
          highDensityCase[[Range[1,  101,  20] ]]  
   ], 2] ]];

At low densities, the lattice configuration remains disordered over time with a homogeneous distribution of the holes and the A and B particles. The current in this system is relatively high. 

As the density of particles increases, the particles mutually block and impede one another, giving rise to inhomogeneities in the distribution of the holes and the A and B particles, and in the charge density and mass density  (particle density regardless of charge) distributions. Eventually, there is a transition resulting in phase segregation and the emergence of a striped, banded structure. The space between the strips consists primarily of holes. Within each strip, there are two substrips of approximately equal size; one comprised of A particles and one comprised of B particles. The  'locking up' the system by the 'mutual blocking' effect significantly reduces the current in the system.

Note: The transition in the distribution of particles (and hence of charge) as a function of particle density is analagous to the 'jamming' transition seen in another driven diffusive system, traffic flow [see Gaylord and Wellin 1995].

Computer Simulation Projects

1. We can locate the critical value at which the first order phase transition in the two-species driven diffusive system occurs by visually identifying the lowest density at which an inhomogeneous banded structure appears in the lattice. It should also be possible to use some quantitative measure whose value changes abruptly at the transition point (eg., the density profile of the holes is quite different in the homogeneous and inhomogeneous states). Think of some quantity that can be computed from the lattice configuration, and create a plot of that quantity as a function of density to locate the critical density at which the value of the quantity changes abruptly.

General References

Schmittmann, B. and Zia, R.K.P. Statistical Mechanics of Driven Diffusive Systems. Academic Press, 1995.

Schmittmann, B., Hwang, K., and Zia, R.K.P. "Onset of Spatial Structures in Biased Diffusion of Two Species". Europhysics Letters, vol. 19, 1992. pp. 19 - 25.

Kertesz, J., and Ramawamy, R. "Coarsening in a Driven Diffusive System with Two Species". Europhysics Letters, vol. 28, 1994. pp. 617 - 622.

Gaylord, Richard J. and Wellin, Paul R. Computer Simulations with Mathematica: Explorations in Complex Physical and Biological Systems, TELOS/Springer-Verlag (1995).

ABOUT THE EDITORS

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file DrivenDiffusion62.nb containing Mathematica code for the material described in this article.