"Mathematica in Education and Research" COLUMN

by

Richard J. Gaylord and Kazume Nishidate

Phase Changes in Materials

I. Spinodal Decomposition in Binary Mixtures

Introduction

Binary alloy systems undergo non-equilibrium phase separation in which domains of two stable phases grow from a thermodynamically unstable homogeneous phase as a result of quenching below a critical temperature. 

Traditionally, phase ordering phenomena has been theoretically modeled by numerically solving nonlinear partial differential equations (viz the Ginzburg-Landau and Cahn-Hilliard equations), which amounts to space-time discretizing a continuum model. An alternative approach is to use a space-time discrete lattice from the start. This type of model is referred to as a cell dynamic scheme or CDS (coupled maps and cellular automata are examples of CDS's).

The following three features need to be incorporated into an CDS model of spinoidal decomposition: 

The local tendency to segregate (ie., to increase the local concentration difference between the species). 

The stability of the segregated bulk phases.

The conservation of matter. 

The order parameter, which is the concentration difference between the two components in a binary mixture, normalized between -1 and 1, can be used to describe the configuration of the binary system. The spinoidal decomposition CDS computes the evolution of the order parameter over time using the constraint that the order parameter be conserved.

In addition to using the CDS to model domain formation in spinoidal decomposition, it can be used to explain the periodic patterns formed in the microphase separation in block copolymers, which are long chain molecules (resembling an unclasped pearl necklace) consisting n consecutive monomeric units of type A and m consecutive monomeric units of type B joined together by covalent chemical bonds (this system differs from the spinoidal system in that the covalent bonds between A and B limit the size of the domains).
 
We will present two programs for the deterministic CDS model of the evolution of the order parameter field for phase ordering phenomenon: one program in which the order parameter is conserved and one program in which the order parameter is not conserved.

The Phase Ordering Algorithm

The non-conserved order parameter case:

The model employs a square lattice of size, m, with periodic boundary conditions. The lattice sites have values between -1 and 1 where the value represents the order parameter.

Initially, the order parameter is randomly distributed throughout the lattice in a uniform manner between -0.1 and 0.1.

Table[Random[Real, {-0.1, 0.1}], {m}, {m}]

During each time step, the order parameter values on all of the lattice sites change. In the absense of the constraint that the order parameter be conserved, the change in the order parameter value of a site arises from (a) a diffusive effect due to the difference of the order parameter from the average order parameter of neighboring sites, and (b) a local driving force due to the chemical potential. 

The order parameter values on the lattice sites are updated using the following function whose arguments are the values of a site and the eight nearest neighbors in its Moore neighborhood.

separate[x_, n_, e_, s_, w_, ne_, 
                           se_, sw_, nw_] := 
  1.3 Tanh[x] + d ((n + e + s + w)/6 + 
    (ne + se + sw + nw)/12 - x )

The quantity 1.3 Tanh[x] is a symmetric hyperbolic map. It has one hyperbolic source (which is an unstable fixed point corresponding to the disordered state before quenching) located at the origin and two hyperbolic sinks (which are stable fixed points corresponding to the two ordered states after quenching) symmetrically placed on either side of the unstable fixed point .

Note: The precise form of the hyperbolic map is not important for the CDS modeling.

The quantity d is a positive real number (proportional to the diffusion constant) representing the stability of bulk uniformity.

The term (n + e + s + w)/6 + (ne + se + sw + nw)/12 calculates an isotropic average over all of the sites in a Moore neighborhood except the center site, summing the total value of the neighbors, with sites in the N, S, E, and W direction given twice the weight of sites in the NE, SE, SW, and NW directions.

The separate rule is applied to the lattice over t time steps using

Nest[Moore[separate, #]&, initconfig, t] ]

where

Moore[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]

Finally, a graphics display (greyscale image) of the order parameter field is created using the ListDensityPlot function.

ListDensityPlot[newlat, Mesh->False, 
  FrameTicks->None]

These pieces of code can be combined into a program.

The non-conserved order parameter program:

PhaseOrderingNonConserved[d_, m_, t_]:= 
 Module[{initconfig, Moore, separate},
      
  initconfig = 
   Table[Random[Real, 
    {-0.1, 0.1}], {m}, {m}] ;
       
  separate[x_, n_, e_, s_, w_, ne_, 
                      se_, sw_, nw_] := 
   1.3 Tanh[x] + d ((n + e + s + w)/6 + 
    (ne + se + sw + nw)/12 - x );

Moore[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];
                     
newlat = Nest[Moore[separate, #]&, 
          initconfig, t];
         
ListDensityPlot[newlat, Mesh -> False, 
                FrameTicks -> None]    
]

The conserved order parameter case:

The constraint that the order parameter be conserved essentially demands that when there is an exchange of order parameter values between a site and its neighbors, there is no net change of the order parameter inside the neighborhood of the site. In this case, the lattice sites are updated at each time step using

[Function[y, (# + y - Moore[nnave, y])][Itn[#]]&

where

nnave[x_, n_, e_, s_, w_, ne_, se_, sw_, nw_] := (n + e + s + w)/6 + (ne + se + sw + nw)/12

and

Itn[mat_] := 
  d (Moore[nnave, mat] - mat) + 
     1.3 Tanh[mat] - mat

The conserved order parameter program:

PhaseOrderingConserved[d_, m_,  t_] :=
 Module[{nnave, Itn, newlat},   

  initconfig = 
    Table[Random[Real, {-0.1, 0.1}], {m}, {m}];
      
  nnave[x_, n_, e_, s_, w_, ne_, se_, 
                          sw_, nw_] := 
    (n + e + s + w)/6 + (ne + se + sw + nw)/12;
     
  Moore[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];
       
  Itn[mat_] := 
    d (Moore[nnave, mat] - mat) + 
       1.3 Tanh[mat] - mat;
        
  newlat = 
   Nest[Function[y, 
    (# + y - Moore[nnave, y])][Itn[#]]&,
     initconfig, t];
      
  ListDensityPlot[newlat, Mesh -> False,
    FrameTicks -> None]                      
                              ]

Domain Formation

We show below some typical output from running the phase ordering programs for the conserved order parameter constraint  case and the non-conserved order parameter constraint case using a 100 by 100 lattice with d equal to 0.5, after 0, 20, and 100 time steps.

SeedRandom[8];
PhaseOrderingConserved[0.5, 100, 0];

SeedRandom[8];
PhaseOrderingConserved[0.5, 100, 20];

SeedRandom[8];
PhaseOrderingConserved[0.5, 100, 100];

SeedRandom[8];
PhaseOrderingConserved[0.5, 100, 500];

SeedRandom[8]
PhaseOrderingNonConserved[0.5, 100, 0];

SeedRandom[8]
PhaseOrderingNonConserved[0.5, 100, 20];

SeedRandom[8]
PhaseOrderingNonConserved[0.5, 100, 100];

Computer Simulation Projects

1. The ListDensityPlot indicates variations in the local concentration difference between species  in a montonic fashion (ie., as the order parameter varies from -1 to 1, the greyscale of the ListDensityPlot changes monotonically from 0 to 1). As such, it doesn't allow us to clearly see where one or the other species dominates in an interfacial region. To do this, we can assign different colors to positive and negative values of the order parameter. Use the ColorFunction option with ListDensityPlot to do this.

Hints:
When the ListDensityPlot is applied to the order parameter lattice, the values which vary from -1 to 1 are automatically rescaled to vary from 0 to 1. Therefore in using the ColorFunction, the value 0.5 denotes the point at which there is a change in which species exists in greater local concentration.

ColorFunction -> 
 (If[# < 0.5, RGBColor[2 #, 0, 0],  
   RGBColor[0, 2 (# - 0.5), 0]]&)

where the second argument in the anonymous function results in order parameter values between -1 and 0 (which have been rescaled by the ListDensityPlot function to run from 0 to 0.5) being colorized between black (since 2 0.0 = 0) and red (since 2 0.5 = 1) and the third argument in the anonymous function results in order parameter values between 0 and 1 (which have been rescaled by the ListDensityPlot function to run from 0.5 to 1) being colorized between black (since 2 (0.5 - 0.5) = 0) and green (since 2 (1.0 - 0.5) = 1 ).  

One problem with using the above ColorFunction option is that  both colors shade to black at one limit. We can make the coloring range from dark red to red and from dark green to green for the two species by using

ColorFunction -> 
 (If[# < 0.5, 
    RGBColor[(0.2 + 0.8 (2 #)), 0, 0], 
    RGBColor[0, (0.2 + 0.8 (2 (# - 0.5))), 
   0]]&)

Note:  We can use CMYKColor or any of Mathematica's other colour functions in a similar way.

2. A very simple cellular automaton, known as the Vote CA, shows the same 'phase ordering'  phenomenon produced by the detailed CDS models given here.

The Vote CA can be described as follows:

The sites of a square lattice have values of 0 and 1. At each update, a site's value is replaced by the value possessed by the majority of the nine sites in its Moore neighborhood, with the following exceptions: if there is one more neighborhood site with a value of 1 than with a value of 0 (ie., 5 sites with value 1 and 4 sites with value 0 ), the site value is updated to 0 and if  there is one less neighborhood site with a value of 1 than with a value of 0  (ie., 4 sites with value 1 and 5 sites with value 0), the site value is updated to 1. The rules of this CA are given in the table below:

Print["                    The CA Vote Rule Table"];
TableForm[{Range[0, 9], 
 {0, 0, 0, 0, 1, 0, 1, 1, 1, 1}}, 
TableHeadings->
 {{"Sum over neighborhood", 
   "New cell value"},
  None}]

The Vote CA program can be written as:

Vote[s_, t_]:=
 Module[{rule, init, Moore},
  init = Table[Random[Integer], {s}, {s}];
       
  Moore[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];
       
 rule[4] := 1;
 rule[5] := 0;
 rule[x_] := Floor[x/5]; 
 Attributes[rule] = Listable;                  

 newlat = Nest[rule[Moore[Plus, #]]&, init, t];
       
 ListDensityPlot[newlat, Mesh -> False, 
    FrameTicks -> None]  
               ]

Entering and running the Vote CA program yields:

SeedRandom[3]
Vote[100, 0];

SeedRandom[3]
Vote[100, 50];

Compare the graphics produced by the PhaseOrdering and Vote programs.

Hint: To make a proper comparision of the PhaseOrdering programs with the Vote CA program, the quantity newlat in the ListDensityPlot function should be replaced with Sign[newlat].

General References

Oono, Yoshitsugu and Puri, S. 1987. "Computationally Efficient Modeling of Ordering of Quenched Phases", Phys. Rev. Lett., 58, 836 - 839.

Oono, Y. and Puri S. 1988. "Study of phase-separation dynamics by use of cell dynamical systems. I. Modeling", Phys. Rev. A, 38, 434 - 453.

Note: This article contains the statement " Nature gives physicists phenomena, not equations". This could serve as a rallying cry of the new field of algorithmic physics in which phenomena are directly modeled by algorithms or computer programs rather than by equations.

Puri S. and Oono, Y. 1988. "Study of phase-separation dynamics by use of cell dynamical systems. II. Two-dimensional demonstrations", Phys. Rev. A, 38, 1542 - 1565.

Oono, Yoshitsugu and Shiwa, Y. 1987. "Computationally Efficient Modeling of Block Copolymer and Benard Pattern Formations", Mod. Phys. Lett. B, 1, 49 - 55.