Crystal growth

A parametric study of 3D crystal growth with Monte Carlo.

Monte Carlo is a powerful computational method allowing calculation of plethora of physical —but not only— processes, such as: calculus, market behavior, radioactive decay, diffusion, clusters and percolation, Ising spin model, thermodynamical ensembles, domain and crystal growth, and many more topics in biology, sociology, economy, and finances.

The modern version of Monte Carlo involving Markov chain was invented in the late 1940s by Stanislaw Ulam, a Polish mathematician, while he was working on nuclear weapons projects at the Los Alamos National Laboratory. It was then developed by a number of famous mathematicians and physicists, including John von Neumnann and Nicholas Metropolis. The latter proposed name of the method Monte Carlo referring to a well-known gambling place in Europe (stochastic processes in both cases).

3D Visualization of a crystal growth on a 100 by 100 surface. 107 steps rendered every 105.

The beauty of Monte Carlo lies in its simplicity. In the widely used, Metropolis version of Monte Carlo the randomly chosen steps or direction are accepted or rejected based on a well specified probability. Because it is a Markov process, the probability does not depend on system's history (previous steps). In physics, typically, a step is associated with some change of potential energy, ΔE and the corresponding probability p=-ΔE/kBT. As kBT, Boltzmann constant times temperature, has the meaning of thermal fluctuation and is always positive, a negative ΔE leads to p greater than 1. Conversely, a transition leading an increase in energy will always have p smaller than 1 — the greater the energy change, the smaller the probability. Thus a negative energy transition would be always accepted and a positive one only sometimes.

In this study we are going to simulate a crystal growth on a simple cubic lattice and later we will demonstrate how does bond strength and temperature influence the growth itself. Crystals growth from a melt or a vapor has been a topic of extensive study due to technological implications and more importantly because of a desire to understand the theoretical nature of the growth phenomenon (at the time when the resolution of experimental method was not yet at the atom level). The idea is as follows: the atoms can be adsorbed on or removed from (existing) crystal surface (lattice), or diffuse. The gas or melt phase is not simulated explicitly though so the 'adsorption' process means basically that a lattice site becomes occupied. The reverse argument works for removal / evaporation. No voids or overhangs are allowed and the resultant growth is ‘compact’ so the surface is being built block by block.

Crystal growth on 100 by 100 cubic lattice after 107 steps: α=12, β=1.5, and the maximum layer height is 4.

An adsorbed atom can form a maximum of four bonds on a cubic lattice and it is assumed that neither the 'diagonal' bonds or those with adjacent layers are not taken into account. When adsorbed on a flat surface it has therefore zero neighbors. So-called adatom has one with the terrace it joins. A corner has two neighbors and so-called kink site has three as it tucked into an indent of a terrace. Finally, a hole has the maximum of four neighbors.

We now define two parameters: α and β. α is the average bond strength normalized by the Boltzmann factor and is called thermal roughening. β is the ratio between the difference in chemical potential between the solid and liquid phase and the Boltzmann factor and is called kinetic roughening. It would depend on concentration and increase with saturation. Zero value means the system is at thermodynamical equilibrium. When the code runs, it first decides whether to do: 1) creation, 2) annihilation, or 3) diffusion step. In each case, all possible moves (in terms of energy) are precalculated:

p=exp(-α/4*(2-i)+β/2)

where i is the number of neighbors. The signs flip for annihilation.

An example simulation result is shown above, using a 100 by 100 cubic lattice after 107 steps and the following values: α=12, β=1.5. Within this number of step the maximum height or an increase in the number of layers is 4. We can see that the adsorbed atoms like to 'stick' to each other, forming rounded areas, looking like stains.

The configurations below are modifications of the above simulation. Now, if we modify, for example, β to 4, meaning the liquid above the surface is supersaturated, the highest layer is 32(!) and the overall look is quite chaotic.

Effect of supersaturation: α=12, β=4, and the maximum layer height is 32.

If we now make β=0, the equilibrium situation, we have only a couple etched holes and a few single deposed atoms — no major crystal growth.

Unsaturated liquid: α=12, β=0, and the maximum layer height is 1 but the original layer has been decreased by 5.

Leaving β aside, by making the bonds weaker (setting α to 4) we observe not only faster growth (78 layers!) but also the smaller size of the domains.

Effect of weak bonds: α=4, β=1.5, and the maximum layer height is 78.

In principle, the faster the growth, the more chaotically-looking our crystall. Thus, if we increase the temperature by a factor of two, both of our parameters decrease twice and the effect is shown below, resembling the weaker bonds case.

Effect of increasing temperature by a factor of two: α=6, β=0.75, and the maximum layer height is 24.

In Monte Carlo we can turn off the creation and annihilation steps to focus solely on the diffusion process. In the example below, the initial configuration was a square, single layer, slab in the center. The simulation has been carried out for two cases: one with bonds twice weaker than the other. (β does not matter here as we do not interface with liquid). The weaker bond case results obviously with a faster diffusion.

Only diffusion from a one-layer high central slab: relatively strong (α=12, left) and weak (α=6, right) bonds.

If instead we start from a regular grid of single atoms, we can observe how the atoms like to 'stick together', the effect that resembles bringing together foam spots on your drink. We can look at it as a way of the system to minimize its surface area that has some energetic penalty.

Tucking in together of atoms dispersed on the surface: effect of diffusion and strong bonds (α=20).

Finally, when we have a very high relative bond strength (here, α=20), the atoms are only deposited to the existing 'islands'. In the example below, we start again from a centrally located slab that grows its boundaries over time. The holed and single atoms appear here and there but are then immediatelly occupied (removed).

Slow and neat growth resembling spiral growth using high bond strength (α=20) and low kinetic roughening (β=1).

If you would like to know more details about this topic, feel free to contact us or check the provided links and references. If you enjoyed reading the article and find it useful, like us on Facebook and/or leave a comment.

author: Dr. Karol Kulasinski

References & External Links


G. H. Gilmer and P. Bennema, Computer Simulation of Crystal Surface Structure and Growth Kinetics. Journal of Crystal Growth 13/14 (1972) p. 148-153.

M. Rak, M Izebski, A. Brozi, Kinetic Monte Carlo study of crystal growth from solution. Computer Physics Communications 138 (2001) p. 250–263.

D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, Third Edition, Cambridge Univeristy Press, 2009.

Metropolis Monte Carlo