Snowflakes

by
Richard J. Gaylord and Kazume Nishidate

Mathematica in Education
Vol. 3 No.4
Copyright 1994
TELOS/Springer-Verlag

Introduction

Crystallization is a liquid-solid phase transition in which the solid grows by adding material from the liquid  environment adjacent to it. This process, which can be classified as a type of spreading phenomena [Gaylord and Wellin,1994], is complicated by the fact that there is a growth inhibition effect whereby the crystallization of a site inhibits the crystallization of nearby sites. The molecular causes of crystal growth inhibition, including interfacial surface tension and latent heat diffusion effects, can be incorporated in a detailed lattice model [Gaylord and Nishidate, 1994] but it is also possible to use a very simple cellular automaton model [Packard, 1986] incorporating growth inhibition in a  naive fashion which produces patterns that bear a remarkable resemblance to the observed structures of snowflakes.  We will develop this cellular automaton model for the formation of snowflakes on a hexagonal lattice.

Working with a Hexagonal Lattice

We can display the centers of a small lattice of hexagons, using:

Show[Graphics[
 {PointSize[0.06],  
  Map[Point, 
     Flatten[Table[{i Sqrt[3], j},
                  {i, 1, 10}, {j, Mod[i, 2], 7, 2}],
            1]]
 }], PlotRange -> All, AspectRatio -> Automatic]

The figure illustrates that a hexagonal lattice is essentially a distorted rectangular lattice. Thus it should be possible to model a hexagonal lattice CA using a rectangular lattice CA and converting the results onto a hexagonal lattice. To do this, the neighborhood of a site in a hexagonal lattice needs to be mapped onto the neighborhood of a site in a rectangular lattice. 

Consider the following small hexagonal lattice:

11       12         13       14        15        16

      21       22         BB        24        25        26
          
31        AA        33        34         35       36

     41       42        43         44         45        46

The nearest neighbor sites of element AA in the third row are 21, 22, 31, 33, 41, 42. The nearest neighbor sites of element BB in the second row are 13, 14, 22, 24, 33, 34.

When this hexagonal lattice is converted into a rectangular lattice, these nearest neighbors sites become elements in a Moore neighborhood, which consists of the sites N, NE, E, SE, S, SW, W, NW of a site. The positioning of the nearest neighbor sites differs for sites in odd- and even-numbered rows. We can demonstrate this using our simple hexagonal lattice.

For site AA, the six nearest neighbor sites (21, 22, 31, 33, 41, 42) are positioned as follows:

11   12   13   14   15   16

21  22   BB   24   25   26

31   AA   33   34   35   36

41   42   43   44   45   46

In general, the  shape of the nearest neighbor sites for a site in an odd-numbered row looks like

*  *
*      *
*  *

Hence, the six nearest neighbors of a site in an odd-numbered row on  a hexagonal lattice become the nearest neighbors of the site in the N, E, S, SW, W, SW directions on a rectangular lattice.

For site BB, the six nearest neighbor sites (13, 14, 22, 24, 33, 34) are positioned as follows:

11   12   13   14   15   16

21  22   BB   24   25   26

31   AA   33   34   35   36

41   42   43   44   45   46

In general, the  shape of the nearest neighbor sites for a site in an even-numbered row looks like

*  *
*      *
    *  *

Hence, the six nearest neighbors of a site in an even-numbered row on a hexagonal lattice become the nearest neighbors of the site in the N, NE, E, SE, S, W directions on a rectangular lattice.

We can now develop the algorithm for snowflake formation using a rectangular lattice.

The Snowflake Algorithm

1. The crystallization lattice begins as a Boolean lattice with m rows and n columns, where m and n are even numbers (while not specifically necessary for the snowflake CA, using even values of m and n makes it possible to employ periodic boundary conditions for other hexagonal lattice CA's).  All the  sites have value 0 except the center site which has value 1 (this site serves as the seed for the crystallization process)

seed = 
ReplacePart[Table[0, {m}, {n}], 
           1, {Ceiling[m/2], Ceiling[n/2]}]

2. The rules for snowflake evolution  are:

(1) the value of a sites changes to 1 if it has exactly one nearest neighbor site whose value is 1. 

(2) the values of other sites remain unchanged.

These rules can be expressed in the following set of thirteen rules:

snowflake[a_?OddQ, 0, 1, 0, 0, 0, 0, 0, _, _] := 1;
snowflake[a_?OddQ, 0, 0, 1, 0, 0, 0 ,0, _, _] := 1;
snowflake[a_?OddQ, 0, 0, 0, 1, 0, 0 ,0, _, _] := 1;
snowflake[a_?OddQ, 0, 0, 0, 0, 1, 0 ,0, _, _] := 1;
snowflake[a_?OddQ, 0, 0, 0, 0, 0, 1 ,0, _, _] := 1;
snowflake[a_?OddQ, 0, 0, 0, 0, 0, 0 ,1, _, _] := 1;

snowflake[a_?EvenQ, 0, 1, 0, 0, 0, _, _, 0, 0] := 1;
snowflake[a_?EvenQ, 0, 0, 1, 0, 0, _, _, 0, 0] := 1;
snowflake[a_?EvenQ, 0, 0, 0, 1, 0, _, _, 0, 0] := 1;
snowflake[a_?EvenQ, 0, 0, 0, 0, 1, _, _, 0, 0] := 1;
snowflake[a_?EvenQ, 0, 0, 0, 0, 0, _, _, 1, 0] := 1;
snowflake[a_?EvenQ, 0, 0, 0, 0, 0, _, _, 0, 1] := 1;

snowflake[_, b_, _, _, _, _, _, _, _, _] := b;

Attributes[snowflake] = Listable

The ten arguments to snowflake are, in order:

argument 1: The value of a site in the reference lattice, OddEven, which indicates the number of the row that a site in the crystallization lattice lies in (this is used to determine if a site is in an odd- or even-numbered row).

OddEven = Table[i, {i, m}, {n}]

argument 2: The value of a site in the crystallization lattice.

arguments 3 - 10: The values of the 8 nearest neighbor sites in the Moore neighborhood of a site in the crystallization lattice. The matrices of these nearest neighbors can be written as follows,  where lat represents the crystallization lattice.

RotateRight[lat, {-1, 0}]   - the matrix of nearest neighbor sites in the West direction
RotateRight[lat, {1, 0}]    - the matrix of nearest neighbor sites in the East direction
RotateRight[lat, {0, 1}]    - the matrix of nearest neighbor sites in the North direction
RotateRight[lat, {0, -1}]   - the matrix of nearest neighbor sites in the South direction
RotateRight[lat, {1, 1}]    - the matrix of nearest neighbor sites in the Northeast direction
RotateRight[lat, {-1, 1}]   - the matrix of nearest neighbor sites in the Northwest direction
RotateRight[lat, {1, -1}]   - the matrix of nearest neighbor sites in the Southeast direction
RotateRight[lat, {-1, -1}]  - the matrix of nearest neighbor sites in the Southwest direction

Note: For a given site, the values of only six of the last eight argumnets of snowflake are relevant for updating the site: if the site is in an odd row, the 3rd to 8th arguments of snowflake are the relevant sites and if the site is in an even row, the 3rd to 6th and 9th to 10th arguments of snowflake are the relevant sites.

We can now write the program for snowflake formation.

The Snowflake Program

Snow[m_?EvenQ, n_?EvenQ, t_Integer] :=
        Module[{seed, OddEven, snowflake, crystallize},

snowflake[a_?OddQ, 0, 1, 0, 0, 0, 0, 0, _, _] := 1;
snowflake[a_?OddQ, 0, 0, 1, 0, 0, 0 ,0, _, _] := 1;
snowflake[a_?OddQ, 0, 0, 0, 1, 0, 0 ,0, _, _] := 1;
snowflake[a_?OddQ, 0, 0, 0, 0, 1, 0 ,0, _, _] := 1;
snowflake[a_?OddQ, 0, 0, 0, 0, 0, 1 ,0, _, _] := 1;
snowflake[a_?OddQ, 0, 0, 0, 0, 0, 0 ,1, _, _] := 1;

snowflake[a_?EvenQ, 0, 1, 0, 0, 0, _, _, 0, 0] := 1;
snowflake[a_?EvenQ, 0, 0, 1, 0, 0, _, _, 0, 0] := 1;
snowflake[a_?EvenQ, 0, 0, 0, 1, 0, _, _, 0, 0] := 1;
snowflake[a_?EvenQ, 0, 0, 0, 0, 1, _, _, 0, 0] := 1;
snowflake[a_?EvenQ, 0, 0, 0, 0, 0, _, _, 1, 0] := 1;
snowflake[a_?EvenQ, 0, 0, 0, 0, 0, _, _, 0, 1] := 1;

snowflake[_, b_, _, _, _, _, _, _, _, _] := b;

Attributes[snowflake] = Listable;

seed = ReplacePart[Table[0, {m}, {n}], 1,
                           {Ceiling[m/2], Ceiling[n/2]}];

OddEven = Table[i,{i, m}, {n}];

crystallize = snowflake[OddEven,
                        #,
                        RotateRight[#, {-1, 0}],
                        RotateRight[#, {1, 0}],
                        RotateRight[#, {0, 1}],
                        RotateRight[#, {0, -1}],
                        RotateRight[#, {1, 1}],
                        RotateRight[#, {-1, 1}],
                        RotateRight[#, {1, -1}],
                        RotateRight[#, {-1, -1}]]&;

Apply[Plus, NestList[crystallize, seed, t]]
]

Note: While we could create a GraphicsArray showing the snowflake structure after each step of the crytstallization process using the result of the NestList operation in the program, this would require an immense amount of memory for as few as 30 times steps. Instead, we will create a single graphics in which the sites added to the crystal at the various time steps are distinguished by their coloring. This is done by applying Plus to the NestList  result to produce a single list of lattice site values in which the sites that crystallize in the sth time step have the value (t + 1 - s) and these values are taken as arguments to the Hue function in the graphics program below.

Creating A Hexagonal Lattice Graphic

SquareToHex maps from the rectangular lattice to the hexagonal lattice.

SquareToHex[test_] :=
  Module[{},
   Show[Graphics[
    Table[
      {Hue[#, #, #]&[test[[j]][[i]]/Max[Flatten[test]]],
         
       Polygon[
               Map[(# + {i 2 Sin[Pi/3] - Mod[j, 2] Sin[Pi/3],
                   3(j - 1) Cos[Pi/3]})&,
                    Table[{Sin[n Pi/3], Cos[n Pi/3]},{n,0,5}]
                  ]
              ]
       },
         {i, Length[test]}, {j, Length[test]}]],
       Frame -> False,
       AspectRatio -> Automatic, DisplayFunction -> Identity
        ]
         ]

Here, we show the crystal formed in five steps.

Show[SquareToHex[Snow[26, 26, 5]], 
                DisplayFunction -> $DisplayFunction];

Computer Simulation Projects

1. Create graphics of snowflake growth over 35 time steps, showing that as the crystal grows, its structure cycles through the sequence, plate -> dendritic arms -> sidebranches -> plate ->... .

General References

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

Gaylord, Richard J. and Nishidate, Kazume. Solidification. Mathematica World, November (1994).

Packard, Norman H. Lattice models for solidification and aggregation, Proceedings of the First International Symposium for Science on Form, Tsukuba University   (1986)  (reprinted in Wolfram, Stephen. Theory and Applications of Cellular Automata, World Scientific Pub. Co. (1986) pp. 305-310).