\magnification=1200
\baselineskip =16pt
\parindent=20pt
\parskip=3pt
\raggedbottom
\centerline {\bf A Walk in Phase Space: Solidification into
Crystalline and Amorphous States}
\centerline{Joan Adler}
\centerline{\sl Department of Physics, Technion, Haifa, Israel,
32000}
\medskip The formation of crystalline and amorphous solids is
described using a simple glass bead demonstration and a discussion
of annealing and rapid quenching on a computer.
\bigskip \noindent {\bf I. INTRODUCTION}
This article discusses the use of visualization as a tool for
exploring the solidification of a system into crystalline or
amorphous states. The demonstrations and the Monte Carlo simulation
can be used to help students gain an intuitive understanding of
these states without the need for a detailed discussion of
statistical mechanics. Alternatively, the material is suitable for a
statistical mechanics or a computational physics course after the
simulation of Ising models has been discussed.
A crystalline solid exhibits translational invariance so that the
atoms arrange themselves in a regular pattern, with a specific
lattice spacing between neighboring atoms, and the same number of
nearest neighbor atoms for each atom. In the amorphous state the
number of neighboring atoms varies, and there is no regular pattern
over long distances, although varying amounts of local order may be
present. The latter description is similar at
first glance to a liquid at a given time, but an amorphous state
differs from a liquid in detail. Atoms in the amorphous state
do not move far from their equilibrium sites, whereas in a liquid
such movement is common. If a liquid is cooled instantaneously
(quenched), it will usually go to an amorphous state. In contrast,
crystalline states are obtained by slow cooling, with frequent
small heatings during the cooling process to remove fault lines.
This process is known as annealing.
\medskip \noindent {\bf II. DEMONSTRATION}
We start with a square or rectangular
clear plastic container. Fishing tackle boxes
are perfect. We then obtain 20--100 glass beads
of the type used in chemical filters. A size similar to small ball
bearings is good. The beads can be obtained very cheaply at
chemical stock rooms in most universities. Spread the beads in a
single layer that covers about half the base of the box. The box
is held at 20--30 degrees to the horizontal and lowered towards the
horizontal.
A slow lowering of the box with occasional small shakings mimics the
annealing process and gives a perfect crystal (see Fig.~1(a)). Note
that the crystalline state in this case is one where each bead has
six nearest neighbors and is centered at the sites of a triangular
lattice. This lattice is the most efficient way to pack glass
beads into a two-dimensional space. Rapid lowering corresponds to a
quench and an amorphous structure is obtained (see Fig.~1(b)).
Defects such as vacancies and grain boundaries can be induced with
a little practice (a vacancy is shown in Fig.~1(c)).
These pictures were made by placing the box on a copy machine. In
class, multiple boxes can be used (closed boxes are recommended in
this case) or the demonstration can be made on an overhead
projector. If using an overhead or copy machine, a coin can be
placed to hold the box at a degree or two from horizontal. To use
an overhead, clear glass beads are essential; for other uses steel
ball bearings would be just as good, but are much
more expensive.
\medskip \noindent {\bf III. COMPUTER SIMULATION}
We start by formulating a model system. First, we choose
the interaction between two atoms. For simplicity, we
choose the two-body Lennard-Jones potential $V(r)$ between two atoms
a distance $r$ apart to be given by
$$V(r) = 4\epsilon \bigl [({r \over \sigma})^{12} - ({r \over
\sigma})^6 \bigr] . \eqno(1)$$
The parameter $\epsilon$ sets the depth
of the attractive well of the potential and the parameter $\sigma$
sets the range of the hard core repulsion. This potential is
accurate for noble gases (except Helium) with the standard example
being Argon, and often is used for generic simulations in other
systems owing to its relative simplicity. Simulations of real
materials, especially metals such as Aluminium, must be made with
more realistic potentials than Lennard-Jones.$^1$
Next we choose the number of atoms and their initial positions, and
the temperature and the geometry of the space in which they are
confined. We also need to think about what type of final state we
desire, and the positions and temperature needed in this state. A
common initial state is the completely random
state. If the desired state is the crystalline state, whose structure
may not be known in advance, then we must be very careful that we
actually reach equilibrium. The ground
state, the lowest energy state at zero temperature, is the
triangular lattice if the atoms are allowed to move only in a
two-dimensional space. We will consider the canonical ensemble for
which the temperature
$T$ and number of particles
$N$ is fixed, and the probability $P_s$ of a particular microstate
$s$ with energy $E_s$ is given by$^2$
$$P_s={1 \over Z}e^{-E_s/kT} . \eqno(2)$$
The function $Z=\sum^M_{s=1} e^{E_s/kT}$ is the partition function,
and
$M$ is the number of accessible states.
The Monte Carlo method is a way of generating states such that the
probability of a state is given by the appropriate distribution for
the chosen ensemble. The process of moving from one microstate to
another is a stochastic process. This process must sample the
states which are most important and obey the principle of
detailed balance whereby the ratio of the probability of going from
state
$j$ to state $i$ to that of moving from $i$ to $j$ is equal to
$\exp(-(E_j-E_i)/kT)$.$^2$ The algorithm commonly used to make these
moves in the canonical ensemble is the Metropolis algorithm.
The bead systems were constrained by hard walls, but free
boundaries with an effective neighbor correction at the edge are
used in the simulations. The effective neighbor correction is
implemented by using an effective temperature for each atom
depending on the atom's current coordination number (number of
neighbors within some distance), so that the probability of a change
of position is approximately independent of the actual number of
neighbors. This scheme facilitates the spatial rearrangement of the
particles and so is ideal for modeling the different
transformations within a solid. For simulations where spatial
reorganization is expected (for example when modeling stresses on
lattice mismatched adsorbed systems), this method is superior to
the commonly used periodic boundary conditions. However, the method
is unsuitable for simulations of a liquid because without walls or
an ``infinite'' sample to restrain it, a liquid would decompose
into droplets.
A common mistake is to allow the system to become stuck in a local
minimum or metastable state and to interpret this local minimum as
the ground state. One way to avoid this trap is to monitor a
variable such as the mean lattice spacing and to make sure that its
mean does not change if the simulation is run twice as long as was
originally thought to be sufficient. Another useful precaution,
especially for a random starting state, is to ensure that different
starting points give the same final result. We also can start at
different temperatures to reach the same final result. For long
simulations aimed at reaching the ground state, an easy check is to
take the final result, make random displacements in each atom's
position (unphysical of course), and check that the system returns
to the ground state. However, when the amorphous state is the
desired final state, we want to end in a metastable state and thus
only the route of repeated simulations from different starting
points provides a suitable check.
An earlier version of the program described below was given in
Ref.~3; the new version has some additional features and is
available at our Web site.$^5$ For the present examples we use
$N=100$ atoms. The initial states and temperatures will be varied
depending on the desired final state. At each temperature, we find
the equilibrium configuration by repeating the basic steps of the
Metropolis algorithm until no further changes in some averaged
quantity such as the energy are observed. At each step one particle
is randomly selected, and a random change is made to its
coordinates. The energy before ($E_0$) and after ($E_f$) the change
is compared, and if the energy is decreased, the change is
accepted. If the energy is increased, the step is accepted with
probability $p$ given by
$$p = e^{-6(E_f-E_0)/zkT}, \eqno(3)$$
where $T$ is the temperature,
$z$ is the coordination number of the particle (number of nearest
neighbors), and $k$ is Boltzmann's constant. The factor of six in
Eq.~(3) is included because the average coordination number in the
interior is six. We need to include the local value of $z$
explicitly to make the effective neighbor correction at the edges
of the sample. The possibility of accepting a change that gives a
local increase in energy enables the system to escape from local
minima. Local minima here correspond to fault lines, vacancies, or
other defects.
It is helpful to visualize this process while doing the simulations
by doing an animated presentation of points in a two-dimensional
space. We used PGPLOT,$^4$ which is a simple library of Fortran
routines that can be called from C or Fortran and is free for
non-commercial use. We have placed our code on the Technion
Computational Physics Web site,$^5$ where it may be downloaded for
non-commercial use. The code works on UNIX and Linux systems. Also
available are C and Java programs.
As in a laboratory experiment, we can grow the desired crystal or
amorphous sample by selecting a suitable initial state and
temperature, and varying the temperature until we reach the
desired final state. To model the solidification process into an
amorphous state, we drastically reduce the temperature. An example
of such a process is given in Fig.~2, where the initial state was
completely random (corresponding to infinite temperature), and the
temperature is immediately quenched. In this case we do not want to
achieve equilibration at each temperature, because our final state
is metastable.
An example of a crystalline state is shown in Fig.~3. All the
central atoms are neatly lined up in a triangular structure, but
the boundaries are not quite perfect. It is likely that the surface
also is not perfectly aligned in a crystal grown in the laboratory.
This particular state was begun from a random initial state, but
was very slowly cooled to the crystalline ground state. In another
attempt to reach the ground state a sample with two differently
oriented domains (Fig.~4) was obtained. This system was then
reheated and cooled again to obtain the perfect crystal of Fig.~3.
It also is possible to model the processes described above by
molecular dynamics. However, molecular dynamics requires more
sophisticated programming to solve the differential equations of
motion than is required by the present algorithm, and is less
directly connected to the principles of statistical mechanics. (The
introduction of temperature into molecular dynamics simulations is
far more subtle than in the present case, where the introduction of
temperature is very natural.) In realistic simulations of
solidification it usually is best to do molecular
dynamics for some time to model a process with a real dynamics
and then to find the local or global minimum with an anneal or a
quench using Monte Carlo.
The processes described above are the basis for much current
research in computational condensed matter/materials science.
Semiconductors, ceramics, Carbon, and polymer systems can be
modeled in this way with the proper choice of potentials. Of
course, an ab initio quantum mechanical potential that includes the
effects of all the electrons is most accurate, but severely limits
sample size. When the effects to be modeled are not strongly
dependent on the details of the electronic structure, it does not
matter which potential is used, if it has parameters which can be
adjusted to give the correct lattice spacing and ground state
structure. For example, we have seen that Lennard-Jones simulations
give results similar to the glass bead system, implying that the
hard core repulsion is the most important part of the potential.
Carbon must be modeled with a potential capable of describing the
insulating diamond structure and conducting graphite and the
transformations between them. A comparison between an ab initio and
the classical Tersoff potential for Carbon has recently
been presented.$^6$ On our Web site$^5$ there
are links to Web sites of other groups doing similar calculations.
In these calculations, the freedom to move the atoms away from
their ground state lattice structure is an essential feature.
\medskip \noindent {\bf IV. STUDENT ACTIVITIES}
\item{1.} One situation where crystalline order is not desired as a
result of the freezing process is the preparation of ice cream
(crystals would spoil the smooth creamy texture). A very fast
quench which would give the creamiest taste is sometimes
impractical. If we try to freeze ice cream in a freezer, ice
crystals start to form, so the ice cream needs to be repeatedly
removed and stirred to break up the crystals. This behavior
explains why an ice cream machine, which churns the ice cream while
it is freezing, gives such desirable ice cream. Suppose that the
ice cream at a supermarket has ice crystals in it. Suggest whether
it must have melted (perhaps en route to the supermarket), and if
so, whether it refroze rapidly or slowly.
\item{2.} You could prepare ice cream and justify the time
investment and calories with the idea that this preparation is part
of your physics homework. You could also try to find other examples
of the application of statistical mechanics in food preparation
processes. For example, preparing a jello mold is a polymerization
phase transition to an amorphous glassy state. Beautiful bubble
rafts, sometimes exhibiting crystalline order and illustrating
crystal defects can be made while doing the dishes if you are
lucky enough to find the right small opening to form regular
bubbles from the detergent. Think of other everyday examples.
\item{3.} Prepare a glass bead simulation as described above and
obtain the different states shown in Fig.~1.
\noindent The following problems require downloading the Monte
Carlo program$^5$ or writing your own.
\item{4.} Determine which annealing schedules lead to crystalline
ground states. Some examples are given in Ref. 3. Observe what
happens when you reheat from the crystalline phase. Do you see any
hysteresis?
\item{5.} What happens when you start from the square lattice which
is a saddle point structure. Are there differences between
annealing and quenching?
\item{6.} Think about the differences between crystals and amorphous
solids, and how you would determine whether you have a basically
crystalline sample with some grain boundaries, or a truly amorphous
system. Keep in mind that in nature there is a continuum between
small crystallites with grain boundaries and a truly amorphous
system with no preferred direction. How difficult is it to make
these distinctions which can be seen so easily by direct
visualization?
\item{7.} Are there observable differences between the structure of a
liquid and that of an amorphous solid?
\medskip \noindent {\bf ACKNOWLEDGEMENTS}
I thank Amihai Silverman for the collaboration leading to the
original version of the program described here, and David Saada,
Adham Hashibon, and Amit Kanigel for current collaborations in
computational condensed matter physics. I acknowledge
collaborations with Rafi Kalish and Wayne Kaplan for their
encouragement to do simulations to describe their experiments. I
acknowledge Alex Zunger for discussions about the importance of
moving simulations off-lattice. The support of the German-Israel
Foundation during the preparation of this manuscript has been
essential.
\medskip \noindent {\bf REFERENCES}
\item{1.}A detailed discussion of interatomic potentials and their
applicability to different materials can be found in N. W. Ashcroft
and N. D. Mermin, {\sl Solid State Physics} (Saunders College,
Philadelphia, 1976), p.\ 398 and H. Gould and J. Tobochnik, {\sl An
Introduction to Computer Simulation Methods}, 2nd ed.
(Addison-Wesley, Reading, MA, 1996), p.~214.
\item {2.} M. Plischke and B. Beregson, {\sl Equilibrium Statistical
Physics}, 2nd ed. (World Scientific, Singapore, 1994).
\item{3.} A. Silverman and J. Adler, ``Animated Simulated
Annealing,'' Computers in Physics {\bf 6}, 277--281 (1992).
\item{4.} Information on PGPLOT can be found at
http://astro.caltech.edu/$\sim$tjp/pgplot.
\item{5.} The Technion Computational Physics home page is at
http://phycomp.technion.ac.il and the page /$\sim$phr76ja/ajp.html
contains a summary of links relevant to this article.
\item{6.} D. Saada, J. Adler and R. Kalish, ``Computer simulation of damage in diamond
due to ion-impact and its annealing'', Phys. Rev. B {\bf 59},
6650 (1998).
\vfill\eject
\noindent {\bf FIGURE CAPTIONS}
\smallskip \noindent Figure 1. Structures found from the
glass beads demonstration. A perfect crystal (a), an
amorphous structure (b), and a crystal with a vacancy (c) is
shown. This figure can be viewed in color at
http://phycomp.technion.ac.il/$\sim$phr76ja/glassballs.html.
%on phjoan5 ~phr76ja/public.html/balls.gif
\smallskip \noindent Figure 2. The amorphous state obtained by a
rapid quench from a random start to a temperature just above
$T=0$. An animation of this figure can be
viewed in color by downloading the file,
ftp://phjoan12.technion.ac.il/pub/submissions/AJP/fig2.ps, and
viewing it with Ghostscript.
%on phjoan7 ~phr76ja/papers/theme/fig2.ps
\smallskip \noindent Figure 3. The crystalline state obtained
after a long slow cooling with occasional heatings.
%on phjoan7 ~phr76ja/papers/theme/fig3.ps
\smallskip \noindent Figure 4. An intermediate state in the
cooling process, which after reheating eventually resulted in a
picture similar to Fig.~3.
%on phjoan7 ~phr76ja/papers/theme/fig4.ps
\bye