Dynamic fracture of icosahedral model quasicrystals:
A molecular dynamics study
Frohmut Rösch, Christoph Rudhart, Johannes Roth, and Hans-Rainer Trebin
Institut für Theoretische und Angewandte Physik,
Universität Stuttgart, Pfaffenwaldring 57, 70550 Stuttgart, Germany
Peter Gumbsch
Institut für Zuverlässigkeit von Bauteilen und Systemen,
Universität Karlsruhe, Kaiserstr. 12, 76131 Karlsruhe, Germany and
Fraunhofer Institut für Werkstoffmechanik, Wöhlerstr. 11, 79108 Freiburg, Germany
Ebert et al. [Phys. Rev. Lett. 77, 3827 (1996)] have fractured icosahedral Al-Mn-Pd single crystals
in ultrahigh vacuum and have investigated the cleavage planes in-situ by scanning tunneling mi-
croscopy (STM). Globular patterns in the STM-images were interpreted as clusters of atoms. These
are significant structural units of quasicrystals. The experiments of Ebert et al. imply that they are
also stable physical entities, a property controversially discussed currently. For a clarification we
performed the first large scale fracture simulations on three-dimensional complex binary systems.
We studied the propagation of mode I cracks in an icosahedral model quasicrystal by molecular
dynamics techniques at low temperature. In particular we examined how the shape of the cleavage
plane is influenced by the clusters inherent in the model and how it depends on the plane structure.
Brittle fracture with no indication of dislocation activity is observed. The crack surfaces are rough
on the scale of the clusters, but exhibit constant average heights for orientations perpendicular to
high symmetry axes. From detailed analyses of the fractured samples we conclude that both, the
plane structure and the clusters, strongly influence dynamic fracture in quasicrystals and that the
clusters therefore have to be regarded as physical entities.
PACS numbers: 62.20.Mk, 61.44.Br, 02.70.Ns
Keywords: fracture, quasicrystals, molecular dynamics simulations
I. INTRODUCTION
Quasicrystals are intermetallic compounds with long-
range quasi-periodic translational order. They possess
well-defined atomic planes and hence diffract electromag-
netic and matter waves into sharp Bragg spots. But they
also display atomic clusters as basic building blocks1,2,
whose arrangement in space is compatible with the pla-
nar structure. These clusters consist for example of sev-
eral shells of icosahedral symmetry (Bergman-, Mackay-,
pseudo-Mackay-clusters). Or they form polytopes, e.g.
decagonal prisms, which like the unit cells of periodic
crystals fill space, although with large overlaps (“quasi-
unit-cell picture”)3. Janot and others4,5,6 have postu-
lated that a self-similar hierarchical assembly of the clus-
ters is responsible for the stability of quasicrystals and
for many physical properties, like the low electric conduc-
tivity. However, it is a controversial and persistent dis-
cussion, whether the clusters are merely structural units
or whether they represent physical entities. The discus-
sion was fueled by an experiment of Ebert et al.7, where
icosahedral Al-Mn-Pd was fractured under ultrahigh vac-
uum conditions at room temperature. Scanning tunnel-
ing microscopy images of the cleavage planes revealed
elements of 0.6 to 1 nm in diameter8. The authors ar-
gue that these are the clusters which were circumvented
by the crack and hence form highly stable aggregates of
matter. Others point out that flat terraces evolve on
fivefold surfaces of i-Al-Mn-Pd when annealed at high
temperatures9,10,11,12. As this requires truncated rows of
Bergman and Mackay clusters it is stated13 that these
therefore could not represent firm entities.
In the present article we report on molecular dynamics
simulations of crack propagation in a three-dimensional
icosahedral model quasicrystal at low temperature. Seed
cracks are inserted along different planes and therein
along different directions. The fracture planes are care-
fully analysed to answer the role of clusters in dynamic
fracture. In Sec. II we provide some requirements on the
theoretical description of fracture. In Sec. III the model
quasicrystal, the molecular dynamics technique, and the
methods to visualise the results of the simulations are
outlined. Subsequently, in Sec. IV the simulation results
are presented and then discussed in Sec. V.
II. FRACTURE
The stress concentration and the strength of the load-
ing at a crack tip are determined by the macroscopic
geometry and dimensions of a sample. In linear elastic
continuum mechanics, a sharp mode I (opening mode)
crack is characterised by a singular stress field and a cor-
responding displacement field, which both are propor-
tional to the stress intensity factor K. This factor is
proportional to the applied external load and contains
the geometry of the sample. A simple energy based con-
dition for crack propagation is the Griffith criterion14.
It states that a crack is in equilibrium when the change
in mechanical energy per unit area of crack advance –
the energy release rate G – equals the change in surface
energy of the two fracture surfaces, 2γ. In continuum
mechanics the energy release rate is proportional to the
square of the stress intensity factor for a given mode. A
crack then should start moving when the stress intensity
factor exceeds the critical Griffith value.
A continuum mechanical description of fracture, how-
ever, has a few drawbacks. First, the requirements for
linear elasticity are no longer valid near the crack tip
where atomic bonds clearly become non-linear and even-
tually break. Second, a continuum theory neglects the
discrete nature of the lattice. Thus, it is fully ignored
that fracture of materials is ultimately caused by bond
breaking processes on the atomic scale.
A way to understand the processes is to perform nu-
merical experiments, since experimental information on
this length scale is difficult to obtain. Molecular dynam-
ics studies have provided useful insight into crack prop-
agation in pure metals and simple intermetallic alloys,
whereas in complex metallic alloys the mechanisms are
not yet so clear. Atomistic studies show for example
that cracks remain stable in a region around the critical
stress intensity factor due to the discrete nature of the
lattice. This effect is called lattice trapping15. A further
consequence of the discrete lattice is that the fracture
behaviour in one and the same plane can depend on the
crack propagation direction16. Such observations cannot
be explained by a simple continuum mechanical descrip-
tion.
III. MODEL AND METHOD
A. Icosahedral binary model
In the numerical experiments we use a three-
dimensional model quasicrystal which has been proposed
by Henley and Elser17 as a structure model for the icosa-
hedral phase of (Al,Zn)Mg. This is the simplest possi-
ble model quasicrystal that is stabilised by pair poten-
tials. Furthermore it allows Burgers circuit analysis and
is a prototype of Bergman-type quasicrystals. As we
do not distinguish between Al and Zn atoms, we term
this decoration icosahedral binary model. It can be ob-
tained by decorating the structure elements of the three-
dimensional Penrose tiling, the oblate and the prolate
rhombohedra (see Fig. 1, top). Al and Zn atoms (A
atoms) are placed on the vertices and the midpoints of
the edges of the rhombohedra. Two Mg atoms (B atoms)
divide the long body diagonal of each prolate rhombohe-
dron in ratios τ :1:τ , where τ is the golden mean. Two
prolate and two oblate rhombohedra with a common ver-
tex form a rhombic dodecahedron18,19. To obtain the
icosahedral binary model, in these dodecahedra the atom
at the common vertex is removed and the four neigh-
bouring A atoms are transformed into B atoms. Finally
these atoms are shifted to the common vertex to divide
the edges of the corresponding rhombohedra in a ratio
of 1:τ . Fig. 1 (bottom left) shows the final decoration of
these dodecahedra, in which the B atoms form hexagonal
bipyramides. This modification increases the number of
Bergman-type clusters (see Fig. 1, bottom right) inher-
ent in the structure, leads to a higher stability with the
potentials used, and takes better into account the experi-
mentally observed stoichiometry of the quasicrystal. The
Bergman-type clusters may also be interpreted as build-
ing units of the quasicrystal and are the main feature of
the structure apart from the plane structure.
B. Molecular dynamics technique
As is often done in fracture simulations (see e.g. Abra-
ham20) we use simple Lennard-Jones potentials to model
the interactions. The following facts led to this choice:
First, the very few potentials available for quasicrys-
tals are unsuitable for fracture simulations: These pair
potentials stabilise the quasicrystal only for a fixed vol-
ume. With free surfaces atoms sometimes simply evap-
orate (e.g. with potentials based on those of Hafner et
al.21). Very long-ranged potentials with Friedel oscilla-
tions (e.g. those of Al-Lehyani et al.22) frequently display
such large cohesive energies that nearly no elastic defor-
mation is possible and instead intrinsic cracks develop.
Second, many simulations have proven that model po-
tentials are helping to understand the elementary pro-
FIG. 1: Tiles of the icosahedral binary model decorated with
two types of atoms and the Bergman-type cluster. Top: pro-
late rhombohedron (left) and oblate rhombohedron (right).
Bottom: rhombic dodecahedron (left) and the 45 atoms build-
ing the Bergman-type cluster (right) inherent in the model
quasicrystal. Small grey and large black spheres denote A
and B atoms respectively.
cesses in fracture (see e.g. Abraham20) and so are rea-
sonable when qualitative mechanisms are the centre of
interest. For quantitative results, a known limitation is
the neglect of non-local and many-body interactions.
Third, the Lennard-Jones potentials used23,24 keep the
model stable even under strongest mechanical deforma-
tions or irradiation (introduction of point defects) and
have been used in our group in many simulations of dis-
location motion25 or even shock waves26. The structure
is robust under a wide variation of the potential depths.
Very similar potentials have shown to stabilise the icosa-
hedral atomic structure in a simpler model27. It is also
known since the early fifties that Lennard-Jones poten-
tials favour icosahedral clusters28, indicating that these
potentials are useful for structures like icosahedral qua-
sicrystals.
The minima of the potentials for interactions between
atoms of the same type are set to �0, whereas the mini-
mum of the potential for the interactions between atoms
of different kind is set to 2�0. The conclusions drawn
from our simulations, however, remained essentially un-
changed if all binding energies are set equal. The shortest
distance between two A atoms is denoted r0. All masses
are set to m0. The time is then measured in units of
t0 = r0
m0/�0.
The molecular dynamics simulations were carried out
using the microcanonical ensemble with the program
code IMD29,30. It performs well on a large variety of
hardware, including single and dual processor worksta-
tions and massively parallel supercomputers.
First, we searched for the potential cleavage planes.
According to the Griffith criterion they should be planes
of low surface energy31. To identify these surfaces we re-
lax a specimen and split it into two parts. Subsequently,
the two regions are shifted apart rigidly. The surface en-
ergy is then calculated from the energy difference of the
artificially cleaved and the undisturbed specimen.
In contrast to simple periodic crystals, the atomistic
structure of the planes and therefore also the surface en-
ergy in quasicrystals varies strongly within the family
of planes perpendicular to a fixed axis. In the present
model, surfaces with low interface energy are located
perpendicular to two- and fivefold directions at certain
heights (Fig. 2). Perpendicular to other directions the
plane structure is less pronounced and the minimal sur-
face energies are higher.
Since we are interested in the morphology of fracture
surfaces we apply a sample form that allows us to fol-
low the dynamics of the running crack for a long time.
For this purpose, a strip geometry is used to model crack
propagation with constant energy release rate32. The
samples consist of about 4 to 5 million atoms, with di-
mensions of approximately 450r0×150r0×70r0. Periodic
boundary conditions33 are applied in the direction paral-
lel to the crack front. For the other directions, all atoms
in the outermost boundary layers of width 2.5r0 are held
fixed. An atomically sharp seed crack is inserted at a
plane of lowest surface energy, from one side to about
FIG. 2: Surface energy of cleavage planes perpendicular to
two- and fivefold axes.
one quarter of the strip length. The system is uniaxially
strained perpendicular to the crack plane up to the Grif-
fith load and is relaxed to obtain the displacement field of
a stable crack at zero temperature. Then a temperature
of about 10−4 of the melting temperature is applied34 to
the configurations with and without the relaxed crack.
From the resulting configurations we obtain an averaged
displacement field for this temperature. The crack now is
further loaded by linear scaling of this displacement field.
The response of the system then is monitored by molec-
ular dynamics techniques. The sound waves emitted by
the propagating crack (see Fig. 3 and online movie35)
are damped away outside of an elliptical stadium32 to
prevent reflections.
C. Visualisation
To study crack propagation on an atomic level the se-
lection and reduction of data is of crucial importance.
Due to the large number of atoms required for the study
of crack propagation in three-dimensional systems, it is
not feasible to write out the positions of all atoms during
the simulation, and even less to display all of them. To
obtain a first insight into crack propagation only atoms
FIG. 3: Kinetic energy density. Sound waves emitted by the
propagating crack and the elliptical region without damping
are clearly visible. See online movie35.
FIG. 4: Snapshot of a simulation with some 4 million atoms.
Only atoms with low coordination number are displayed. See
online movie36.
near the fracture surfaces are of interest. Whereas they
can be visualised in periodic crystals by plotting only
those atoms whose potential energy exceeds a certain
threshold, this technique is not applicable for quasicrys-
tals. Because of the largely varying environments the
potential energy varies significantly from atom to atom,
even for atoms of the same type in a defect-free sample.
A solution to this problem is to display only those
atoms whose coordination number is smaller than a cer-
tain threshold. This number is evaluated by counting
atoms within the nearest neighbour distance. Like the
potential energy, in quasicrystals the coordination num-
ber varies from atom to atom, but to a much smaller de-
gree. As atoms near defects have a significantly lower co-
ordination number, it becomes possible to visualise frac-
ture surfaces and dislocation cores. For the A atoms
the threshold for the coordination number is set to 12,
whereas for the B atoms it is set to 14. With this method,
the number of atoms to write out or to display can be
reduced by three orders of magnitude. Fig. 4 shows a
snapshot of a simulation with some 4 million atoms fil-
tered by this technique, which was also used in a movie
that is available online36.
As can be seen from Fig. 4 the fracture surfaces are
rough. Therefore to decide to which fracture surface an
atom belongs we apply the displacement vectors between
the initial configuration with the built-in seed crack and
the fractured sample. The morphology of the fracture
planes is then analysed via geometrical scanning of the
atoms forming the surfaces. The roughness can be visu-
alised by colour coding the height of the surface in a view
normal to the fracture surface.
To investigate the influence of the Bergman-type clus-
ters on cleavage they have to be displayed together with
the fracture surfaces intersecting them. This is done by
restoring the initial sample without crack at zero tem-
perature. The atoms forming the two sample parts are
taken back to their positions in this initial sample and
FIG. 5: Bergman-type cluster cut by a flat surface.
then scanned geometrically. In addition all atoms form-
ing clusters in the vicinity of this surfaces are known. By
displaying only these atoms and the scanned surface one
directly can see where and to which amount clusters are
cut. A problem of this kind of visualisation is shown in
Fig. 5, where a cluster is cut by a flat surface. When
looked-at from above it is obvious that only four atoms
are separated from the rest of the cluster. When looked-
at from below one could get the impression that the clus-
ter is heavily intersected. On a real fracture surface clus-
ters with centres above and below the crack surface are
present. Therefore both views of Fig. 5 are appearing
at the same time in a two-dimensional projection. Thus
such pictures are not very intuitive.
A way out of this dilemma is presented in Fig. 6. Clus-
ters with midpoints above the crack surface are displayed
together with the upper geometrically scanned fracture
surface only, the other clusters are shown together with
the lower fracture surface. As a result we get two pictures
with clusters cut by surfaces. Note that for a qualitative
and quantitative analysis always both pictures or sets of
data are needed.
FIG. 6: Visualisation of the clusters cut by the dynamic crack.
Midpoints of clusters are indicated as black dots, the clusters
are idealised as spheres. The crack propagated from the left to
the right. The upper and lower fracture surfaces are projected
onto an x− z plane.
IV. RESULTS
In this section the results of the numerical simulations
are presented. For practical purposes we define k as the
stress intensity factor K relative to the stress intensity
factor KG at the Griffith load:
k = K/KG.
The orientations of the samples are characterised by the
notation yx, where y is an axis perpendicular to the cleav-
age plane and x is an axis in the crack propagation di-
rection (see Fig. 6). An axis perpendicular to a fivefold
(5) and a twofold (2) axis is denoted pseudo-twofold (p2)
axis.
A. Crack velocities
Simulations were performed for different orientations
with loads in a range from k = 1.1 to k = 2 (see nota-
tions and loads in Fig. 7). Brittle fracture without any
crack tip plasticity is observed irrespective of the orien-
tation of the fracture plane. For loads below k = 1.2,
the crack propagates only a few atomic distances, and
then stops. Hence the energy needed for initiating crack
propagation is about 1.4 times the value predicted by the
Griffith criterion. Therefore, a simple global thermody-
namic description of fracture is not applicable. The min-
imal velocity for brittle crack propagation is about 10%
of the shear wave velocity37 vs. For loads k ≥ 1.2 the
velocity increases monotonically with the applied load.
The crack velocities are in a range of 10-45% of vs (see
Fig. 7). At high loads the crack velocities on fivefold
cleavage planes show higher average velocities than on
the other planes. Velocities for the two different crack
propagation directions on the fivefold planes differ sig-
nificantly at intermediate loads (k = 1.3). This coincides
with ledges that are produced in the fracture surface (see
section IV B and Fig. 8, bottom).
B. Fracture surfaces
To analyse the morphology of the fracture surfaces,
they are displayed as described in section III C. In Fig. 8
the average height is shown in grey, heights above +2r0
are shown in white and heights below −2r0 are shown
in black. The crack propagation direction is from the
left to the right. The initial fracture surface is flat, as
can be seen from the homogeneous regions on the left.
The surfaces resulting from the propagation of the crack,
however, show pronounced patterns of regions with dif-
ferent heights. From the observation of the fracture sur-
faces it is already evident that they are rough and that
the peak-to-valley roughness is of the order of the di-
ameter of the clusters. The peak-to-valley roughness
and the root-mean-square roughness of the height pro-
files both increase38 for higher loads for surfaces without
FIG. 7: Average crack velocities for different loads and orien-
tations.
ledge formation. For fracture surfaces perpendicular to
twofold and fivefold axes the crack fluctuates about a
constant height (in the areas without ledges). In con-
trast, a crack inserted perpendicular to a pseudo-twofold
direction seems to deviate from this plane39.
C. Anisotropy
As can be seen already from the fracture surfaces per-
pendicular to a fivefold axis in Fig. 8 crack surfaces for the
same cleavage plane differ significantly for different in-
plane crack propagation directions. For the orientation
5p2 ledges are produced, while no ledges form for the ori-
entation 52. By visual inspection of the fracture surfaces
(see Fig. 8) and evaluation of height-height-correlation
functions23 it becomes evident that for every orientation
there are distinct angles in the height profile, which show
pronounced height correlations that correspond to mark-
ings in the fracture surfaces or to ledges. These angles
are given in Table I.
D. Clusters
In Fig. 9 the clusters cut by the fracture surfaces are
presented as described in section III C and Fig. 6. It is
obvious from Fig. 9 that the dynamic crack does not per-
fectly circumvent the clusters, but intersects them par-
tially (right side of Fig. 9). These intersections, however,
TABLE I: Angles observed in the height profiles of the frac-
ture surfaces. Angles measured clockwise from the crack prop-
agation direction get a negative sign.
orientation 22 25 52 5p2 p25
angles 0◦,±32◦ 0◦, +32◦,−58◦ 0◦,±36◦ ±18◦ 0◦,±90◦
FIG. 8: Height profiles of sections of the simulated fracture surfaces. Load: k = 1.3, orientation 22 (top), 52 (middle), 5p2
(bottom). The height increases from black (≤ −2r0) to white (≥ +2r0). The scanning sphere has the same size as an atom of
type B.
FIG. 9: Clusters cut by the fracture surfaces. The visualisation technique is described in Sec. III C and Fig. 6. Load: k = 1.3,
orientation: 52.
FIG. 10: Surface energy and density of cluster centres for the
orientation 22. The corresponding fracture surface is shown
in Fig. 8, top.
are much less frequent than for the flat seed cracks (left
side of Fig. 9). More detailed analyses for different ori-
entations validate this statement. For the orientations
perpendicular to twofold and fivefold axes at k = 1.3 the
ratio of clusters cut by the dynamic crack to clusters in-
tersected by flat cuts is approximately 0.6. Additionally
the absolute value of clusters cut by the crack for the
fivefold surfaces is lower than for the twofold surfaces.
Fig. 10 and Fig. 11 display bottom up: The density of
the cluster centres, the surface energy, a cluster in the
proper length scale, the grey coding of the heights in
Fig. 8 (top, middle), and the position of the seed crack
(dashed vertical line). For the twofold fracture surface
the low energy seed crack is located between two peaks
in the cluster density, whereas for the fivefold surface
this seed crack is situated close to a peak of this density.
It is evident from the figures that it is not possible to
circumvent all clusters by a planar cut.
The grey coding is adjusted to the average height of the
fracture surfaces. It is therefore evident from Fig. 11 that
the crack deviates for the orientation 52 from the low en-
ergy cleavage plane of the seed crack to a parallel plane,
reducing the number of cluster intersections dramatically
(see also Fig. 9 and Fig. 8, middle). Samples cut flat at
the new height show slightly higher surface energy (see
also Fig. 11). However, for low loads and low roughness
the actual fracture surfaces of the dynamic cracks have
about 5-15% higher energies than those of the low energy
seed cracks. To assure that the dynamic crack is depart-
ing from the initial plane not in a random manner the
seed crack was built in at the position colour coded as
medium grey in Fig. 11. The resulting fracture surface
had a similar roughness but the crack did not change to
a parallel plane.
FIG. 11: Surface energy and density of cluster centres for the
orientation 52. The corresponding fracture surface is shown
in Fig. 8, middle.
V. DISCUSSION
Taken together the results of our simulations presented
above indicate that the distribution of the clusters is
crucial for the fracture behaviour: First, circumventing
the clusters or intersecting them disturbs the propagat-
ing crack and leads to additional radiation. This man-
ifests itself in the crack speed. Dynamic cracks prop-
agating in fivefold planes with few cluster intersections
are faster than those in twofold planes, where the abso-
lute number of cluster intersections is higher. Also the
generation of ledges for low loads in the orientation 5p2
costs energy and therefore slows down the crack even
further. Second, circumvention of the clusters leads to
characteristic height-variations. Therefore the peak-to-
valley roughness of the fracture surfaces is of the order
of the diameter of the clusters. Third, the observed pat-
terns in the fracture surfaces correspond to lines along
which the clusters are located. The associated angles are
given by the icosahedral symmetry of the sample, namely
0◦, 18◦, 31.72◦, 36◦, 58.28◦, and 90◦ (see Table I). Ledges
seem to be produced only for the smallest angles mea-
sured from the crack propagation direction. Fourth, less
clusters are intersected by the fracture surfaces than by
the flat seed cracks. Fifth, a seed crack at a low energy
cleavage plane deviates to a parallel plane to reduce the
number of cluster intersections in spite of the higher en-
ergy required to form the fracture surfaces. In contrast,
a crack built-in at this new position does not show such
a deviation.
Another observation of the simulations is that the
plane structure of the quasicrystal also influences frac-
ture. The fracture surfaces that are located perpendicu-
lar to the twofold and fivefold symmetry axes show con-
stant average heights.
The three-dimensional quasicrystals give perfect cleav-
age fracture with no indication of any dislocation activ-
ity. This is in contrast to results on two-dimensional
decagonal quasicrystals, where a dislocation enhanced
fracture mechanism has been observed40. However a
corresponding three-dimensional decagonal quasicrystal
would have a periodic direction with a straight dislo-
cation line. In the simulations presented here this di-
rection is also quasiperiodic. As the clusters have a
strong influence on fracture they also may bend and
hinder dislocation lines. So dislocation emission in the
three-dimensional icosahedral quasicrystal modelled here
should be less likely than in the two-dimensional decago-
nal model quasicrystal. Very high stresses are indeed
needed to move dislocations in our model quasicrystal in
molecular dynamics simulations25.
There are also indications for the stability of the
clusters from the electronic structure of quasicrystals:
First, experiments and ab-initio calculations show that
directional bonding may be present within clusters of
quasicrystals41,42. Second, the electrons may addition-
ally stabilise the clusters because of their hierarchical
structure4. Therefore they should be even more stable
than we have modelled them with simple pair potentials.
With this evidence the results concerning the clusters
seem reasonable and should even underestimate their sta-
bility.
So far fracture experiments in ultrahigh vacuum have
only been performed on icosahedral Al-Mn-Pd, which
has a more complicated atomic structure than the icosa-
hedral binary model. Additionally the clusters are not
Bergman-type. Therefore we cannot expect to represent
this structure on an atomic level, when comparing exper-
iments in Al-Mn-Pd with our simulations. Nevertheless,
the size of the clusters, the icosahedral symmetry, and a
distinct plane structure are common features and qual-
itative aspects should be reproduced well, namely the
size and shape of the patterns and the appearance of dis-
tinct angles on the fracture surfaces. This is indeed the
case, as can be seen in Fig. 12. There a geometrically
scanned fracture surface generated in our simulations is
confronted to an STM-image of Ebert et al.7,43 at the
same length scale. As we were able to correlate the ob-
served structures to the clusters in our model, the simi-
larities corroborate the assumption that the clusters are
responsible for the globular structures observed in ex-
periment. More detailed comparisons to fracture experi-
ments in icosahedral Al-Zn-Mg-type quasicrystals would
be desirable, but such comparative data is currently not
available.
We also performed fracture simulations in a C15 Laves
phase with the same Lennard-Jones potentials to give fur-
ther evidence, that the roughness is correlated to the clus-
ters. This C15 structure is built up from deformed pro-
late rhombohedra, one of the major tiles of the quasicrys-
tal (see Fig. 1 and Sec. III A). However, no Bergman-
type clusters are present. When colour coded like the
quasicrystal, fracture surfaces of the C15 cracks lack any
roughness at low loads and cleave smoothly on high sym-
metry planes44.
FIG. 12: Fracture surfaces perpendicular to twofold axes. In
the left picture (simulation) the atomically sharp seed crack
can be seen on the top, whereas below this area the simulated
fracture surface appears. The orientation of the sample is 22,
the load was k = 1.3. The surface was geometrically scanned
with an atom of type B. Thus atomic resolution is achieved.
The right picture (experiment, adopted from Ebert et al.43)
shows an STM-image of icosahedral Al-Pd-Mn cleaved in ul-
trahigh vacuum. As 20r0 is approximately 5 nm the surfaces
are displayed at the same length scale.
When a crack traverses a solid, it leaves a typical non-
equlibrium surface. Thus, up to this point, we have not
dealt with any equilibrium surface. Even the flat Griffith
cuts we introduced to perform the numerical experiments
are no explicit equilibrium surfaces31 by definition. It is
known that i-Al-Pd-Mn surfaces sputtered and annealed
up to about 900 K are rough with cluster-like protru-
sions. Fivefold surfaces annealed at higher temperatures
exhibit flat terraces9,10,11, which are believed to be bulk
terminated45. As fivefold surfaces are not as rough as
twofold surfaces, they are often studied in experiment.
Adsorption on these surfaces, nevertheless, is often very
site-sensitive. A recent review on quasicrystal surfaces
was given by McGrath et al.12. All of these observations
are consistent with our simulations of non-equilibrium
surfaces. Fivefold surfaces experience less roughness than
twofold surfaces (also in agreement with the experiments
of Ebert et al.7,43; see Fig. 8), where more clusters are
intersected. Additionally, since clusters are cut in the
simulations, their binding energy is not so large as to
term them supermolecules. Thus annealing at very high
temperatures can favour flat surfaces11,12. Nevertheless,
selective adsorption may then be related to completion
of clusters.
VI. CONCLUSIONS
We have simulated crack propagation in an icosahe-
dral model quasicrystal. Brittle fracture without any
signature of dislocation emission is observed. The frac-
ture surfaces are rough on the scale of the clusters and
show constant average heights for orientations perpen-
dicular to twofold and fivefold axes. Thus both the plane
structure and the clusters play an important role in frac-
ture. The influence of the clusters is also seen in the
average crack velocities for different orientations, the ob-
served patterns in the fracture surfaces, the anisotropy
with respect to the in-plane propagation direction, and
the smaller amount of clusters cut by the propagating
crack than by planar cuts. The clusters, too, are a rea-
son why the positions of the cleavage planes cannot be
predicted by a simple energy criterion. Since partial clus-
ter intersections occur, the binding energy of the clusters
is not so large as to term them supermolecules. Never-
theless our observations clearly show that they are not
only structural units but physical entities.
Acknowledgments
Financial support from the Deutsche Forschungsge-
meinschaft under contract numbers TR 154/13 and TR
154/20-1 is gratefully acknowledged.
1 C. L. Henley, Phys. Rev. B 43 (1), 993 (1991).
2 V. Elser, Phil. Mag. B 73 (4), 641 (1996).
3 P. J. Steinhardt, H.-C. Jeong, K. Saitoh, M. Tanaka, E.
Abe, and A. P. Tsai, Nature 396 (6706), 55 (1998).
4 C. Janot and M. de Boissieu, Phys. Rev. Lett. 72 (11),
1674 (1994).
5 C. Janot, J. Phys.: Condens. Matter 9, 1493 (1997).
6 H.-C. Jeong and P. J. Steinhardt, Phys. Rev. Lett. 73 (14),
1943 (1994).
7 Ph. Ebert, M. Feuerbacher, N. Tamura, M. Wollgarten,
and K. Urban, Phys. Rev. Lett. 77 (18), 3827 (1996).
8 Ph. Ebert, M. Yurechko, F. Kluge, K. Horn, and K. Ur-
ban, Quasicrystals – Structure and Physical Properties,
WILEY-VCH, ed. H.-R. Trebin, 6.1, 572 (2003).
9 T. M. Schaub, D. E. Bürgler, H.-J. Güntherodt, and J. B.
Suck, Phys. Rev. Lett. 73 (9), 1255 (1994).
10 M. Gierer, M. A. Van Hove, A. I. Goldman, Z. Shen, S.-L.
Chang, C. J. Jenks, C.-M. Zhang, and P. A. Thiel, Phys.
Rev. Lett. 78 (3), 467 (1997).
11 G. Cappello, J. Chevrier, F. Schmithüsen, A. Stierle, V.
Formoso, F. Comin, M. de Boissieu, M. Boudard, T. Lo-
grasso, C. Jenks, and D. Delaney, Phys. Rev. B 65, 245405
(2002).
12 R. McGrath, J. Ledieu, E. J. Cox, and R. D. Diehl, J.
Phys.: Condens. Matter 14 (4), R119 (2002).
13 G. Kasner, Z. Papadopolos, P. Kramer, and D. E.
Bürgler, Quasicrystals – Structure and Physical Properties,
WILEY-VCH, ed. H.-R. Trebin, 2.4, 123 (2003).
14 A. A. Griffith, Philos. Trans. R. Soc. Lond. Ser. A 221, 163
(1921).
15 R. Thomson, C. Hsieh, and V. Rana, J. Appl. Phys. 42 (8),
3154 (1971).
16 P. Gumbsch and R. M. Cannon, Mat. Res. Soc. Bull.
25 (5), 15 (2000).
17 C. L. Henley and V. Elser, Phil. Mag. B 53 (3), L59 (1986).
18 The icosahedral binary model is described and named BIB
in J. Roth, Eur. Phys. J. B 15 (1), 7 (2000).
19 In rings of oblate rhombohedra the number of clusters is
maximised. No overlapping rhombic dodecahedra are gen-
erated. Remaining oblate rhombohedra stay unchanged.
20 F. F. Abraham, Advances in Physics 52 (8), 727 (2003).
21 J. Hafner and M. Krajč́ı, Europhys. Lett. 13 (4), 335
(1990).
22 I. Al-Lehyani, M. Widom, Y. Wang, N. Moghadam, G.
M. Stocks, and J. A. Moriarty, Phys. Rev. B 64, 075109
(2001).
23 F. Rösch, diploma thesis, Universität Stuttgart (2003),
http://elib.uni-stuttgart.de/opus/volltexte/2004/1899/
24 F. Rösch, Ch. Rudhart, P. Gumbsch, and H.-R. Trebin,
Mat. Res. Soc. Symp. Proc. 805, LL9.3.1 (2004).
25 G. D. Schaaf, J. Roth, and H.-R. Trebin, Phil. Mag.
83 (21), 2449 (2003).
26 J. Roth, Phys. Rev. B 71, 064102 (2005).
27 J. W. Roth, R. Schilling, and H.-R. Trebin, Phys. Rev. B
51 (22), 15833 (1995).
28 F. C. Frank, Proc. Roy. Soc. Lond. Ser. A 215 (1120), 43
(1952).
29 J. Stadler, R. Mikulla, and H.-R. Trebin, Int. J. Mod.
Phys. C 8 (5), 1131 (1997).
30 IMD, the ITAP Molecular Dynamics Program.
http://www.itap.physik.uni-stuttgart.de/˜imd
31 To simulate fracture with molecular dynamics it is nec-
essary to determine the optimal place for a low energy
seed crack and to find out the critical strain for setting
the crack into motion. This is always done by determining
the energies of flat cuts. Note, that these do not necessar-
ily represent equilibrium surfaces which may be bent and
curved.
32 P. Gumbsch, S. J. Zhou, and B. L. Holian, Phys. Rev. B
55 (6), 3445 (1997).
33 As quasicrystals show non-periodic translational order we
generate periodic approximants to apply periodic bound-
ary conditions. The 4 to 5 million atoms of our sample
then form the unit cell. Because of this high number of
atoms the configuration should mechanically behave like
the perfect quasicrystal.
34 Here we use equation (2) of Gumbsch et al.32 with a “filled
stadion damping” f ≡ 1.
35 Cracks propagate by breaking bonds between atoms. In
this process, energy is dissipated in the form of acoustic
waves, as clearly visible in the online movie. The kinetic
energy density is colour coded. Regions of high kinetic en-
ergy density are shown in red whereas blue indicates areas
of lower kinetic energy density.
36 Movie (available online) of a cracking model quasicrystal
at a high load (5p2, k = 1.6). By displaying only atoms
with a low coordination number fracture surfaces are visu-
alised. Blue (black in Fig. 1) and red (grey in Fig. 1) balls
correspond to the atoms in the icosahedral binary model.
Atoms on the border, which are not allowed to move, are
coloured in yellow and green. On the left the atomically
sharp seed crack can be seen. Due to the high deformation
near the crack tip, blue atoms loose their neighbours and
show up for short instances.
37 Concerning only the pure phonon term of the elastic en-
ergy in linear elastic continuum theory, the icosahedral
http://elib.uni-stuttgart.de/opus/volltexte/2004/1899/
http://www.itap.physik.uni-stuttgart.de/~imd
quasicrystal behaves like an isotropic solid.
38 The height profiles were determined via geometrical scan-
ning of the fracture surfaces with an atom of type B at
equidistant points separated by 0.2r0. For these meshes
and e.g. for the orientation 52 the root-mean-square rough-
ness increases from 0.39r0 at a load of k = 1.3 (see Fig. 8,
middle) to 0.54r0 for k = 1.6. The peak-to-valley rough-
ness increases likewise from 3.7r0 to 5.1r0. Although the
exact values may depend on the stepwidth (resolution),
scanning sphere size, and scanned area a general tendency
for an increased roughness for increased loads is evident.
39 The question whether this crack finally chooses a plane
perpendicular to a twofold or fivefold axis cannot be an-
swered, as simulations with the required sample-size would
exceed present computer capabilities.
40 R. Mikulla, J. Stadler, F. Krul, H.-R. Trebin, and P.
Gumbsch, Phys. Rev. Lett. 81 (15), 3163 (1998).
41 Only clusters in the bulk give reliable information on their
relative stability. In ab-initio studies, where the number of
atoms is limited, one therefore is restricted to small ap-
proximants.
42 K. Kirihara, T. Nagata, K. Kimura, K. Kato, M. Takata,
E. Nishibori, and M. Sakata, Phys. Rev. B 68, 014205
(2003).
43 Ph. Ebert, F. Yue, and K. Urban, Phys. Rev. B 57 (5),
2821 (1998).
44 F. Rösch, P. Gumbsch, and H.-R. Trebin, unpublished.
45 Z. Papadopolos, G. Kasner, J. Ledieu, E. J. Cox, N. V.
Richardson, Q. Chen, R. D. Diehl, T. A. Lograsso, A. R.
Ross, and R. McGrath, Phys. Rev. B 66, 184207 (2002).
	INTRODUCTION
	FRACTURE
	MODEL AND METHOD
	Icosahedral binary model
	Molecular dynamics technique
	Visualisation
	RESULTS
	Crack velocities
	Fracture surfaces
	Anisotropy
	Clusters
	DISCUSSION
	CONCLUSIONS
	Acknowledgments
	References

  Ebert et al. [Phys. Rev. Lett. 77, 3827 (1996)] have fractured icosahedral
Al-Mn-Pd single crystals in ultrahigh vacuum and have investigated the cleavage
planes in-situ by scanning tunneling microscopy (STM). Globular patterns in the
STM-images were interpreted as clusters of atoms. These are significant
structural units of quasicrystals. The experiments of Ebert et al. imply that
they are also stable physical entities, a property controversially discussed
currently. For a clarification we performed the first large scale fracture
simulations on three-dimensional complex binary systems. We studied the
propagation of mode I cracks in an icosahedral model quasicrystal by molecular
dynamics techniques at low temperature. In particular we examined how the shape
of the cleavage plane is influenced by the clusters inherent in the model and
how it depends on the plane structure. Brittle fracture with no indication of
dislocation activity is observed. The crack surfaces are rough on the scale of
the clusters, but exhibit constant average heights for orientations
perpendicular to high symmetry axes. From detailed analyses of the fractured
samples we conclude that both, the plane structure and the clusters, strongly
influence dynamic fracture in quasicrystals and that the clusters therefore
have to be regarded as physical entities.

Dynamic fracture of icosahedral model quasicrystals:
A molecular dynamics study
Frohmut Rösch, Christoph Rudhart, Johannes Roth, and Hans-Rainer Trebin
Institut für Theoretische und Angewandte Physik,
Universität Stuttgart, Pfaffenwaldring 57, 70550 Stuttgart, Germany
Peter Gumbsch
Institut für Zuverlässigkeit von Bauteilen und Systemen,
Universität Karlsruhe, Kaiserstr. 12, 76131 Karlsruhe, Germany and
Fraunhofer Institut für Werkstoffmechanik, Wöhlerstr. 11, 79108 Freiburg, Germany
Ebert et al. [Phys. Rev. Lett. 77, 3827 (1996)] have fractured icosahedral Al-Mn-Pd single crystals
in ultrahigh vacuum and have investigated the cleavage planes in-situ by scanning tunneling mi-
croscopy (STM). Globular patterns in the STM-images were interpreted as clusters of atoms. These
are significant structural units of quasicrystals. The experiments of Ebert et al. imply that they are
also stable physical entities, a property controversially discussed currently. For a clarification we
performed the first large scale fracture simulations on three-dimensional complex binary systems.
We studied the propagation of mode I cracks in an icosahedral model quasicrystal by molecular
dynamics techniques at low temperature. In particular we examined how the shape of the cleavage
plane is influenced by the clusters inherent in the model and how it depends on the plane structure.
Brittle fracture with no indication of dislocation activity is observed. The crack surfaces are rough
on the scale of the clusters, but exhibit constant average heights for orientations perpendicular to
high symmetry axes. From detailed analyses of the fractured samples we conclude that both, the
plane structure and the clusters, strongly influence dynamic fracture in quasicrystals and that the
clusters therefore have to be regarded as physical entities.
PACS numbers: 62.20.Mk, 61.44.Br, 02.70.Ns
Keywords: fracture, quasicrystals, molecular dynamics simulations
I. INTRODUCTION
Quasicrystals are intermetallic compounds with long-
range quasi-periodic translational order. They possess
well-defined atomic planes and hence diffract electromag-
netic and matter waves into sharp Bragg spots. But they
also display atomic clusters as basic building blocks1,2,
whose arrangement in space is compatible with the pla-
nar structure. These clusters consist for example of sev-
eral shells of icosahedral symmetry (Bergman-, Mackay-,
pseudo-Mackay-clusters). Or they form polytopes, e.g.
decagonal prisms, which like the unit cells of periodic
crystals fill space, although with large overlaps (“quasi-
unit-cell picture”)3. Janot and others4,5,6 have postu-
lated that a self-similar hierarchical assembly of the clus-
ters is responsible for the stability of quasicrystals and
for many physical properties, like the low electric conduc-
tivity. However, it is a controversial and persistent dis-
cussion, whether the clusters are merely structural units
or whether they represent physical entities. The discus-
sion was fueled by an experiment of Ebert et al.7, where
icosahedral Al-Mn-Pd was fractured under ultrahigh vac-
uum conditions at room temperature. Scanning tunnel-
ing microscopy images of the cleavage planes revealed
elements of 0.6 to 1 nm in diameter8. The authors ar-
gue that these are the clusters which were circumvented
by the crack and hence form highly stable aggregates of
matter. Others point out that flat terraces evolve on
fivefold surfaces of i-Al-Mn-Pd when annealed at high
temperatures9,10,11,12. As this requires truncated rows of
Bergman and Mackay clusters it is stated13 that these
therefore could not represent firm entities.
In the present article we report on molecular dynamics
simulations of crack propagation in a three-dimensional
icosahedral model quasicrystal at low temperature. Seed
cracks are inserted along different planes and therein
along different directions. The fracture planes are care-
fully analysed to answer the role of clusters in dynamic
fracture. In Sec. II we provide some requirements on the
theoretical description of fracture. In Sec. III the model
quasicrystal, the molecular dynamics technique, and the
methods to visualise the results of the simulations are
outlined. Subsequently, in Sec. IV the simulation results
are presented and then discussed in Sec. V.
II. FRACTURE
The stress concentration and the strength of the load-
ing at a crack tip are determined by the macroscopic
geometry and dimensions of a sample. In linear elastic
continuum mechanics, a sharp mode I (opening mode)
crack is characterised by a singular stress field and a cor-
responding displacement field, which both are propor-
tional to the stress intensity factor K. This factor is
proportional to the applied external load and contains
the geometry of the sample. A simple energy based con-
dition for crack propagation is the Griffith criterion14.
It states that a crack is in equilibrium when the change
in mechanical energy per unit area of crack advance –
the energy release rate G – equals the change in surface
energy of the two fracture surfaces, 2γ. In continuum
mechanics the energy release rate is proportional to the
square of the stress intensity factor for a given mode. A
crack then should start moving when the stress intensity
factor exceeds the critical Griffith value.
A continuum mechanical description of fracture, how-
ever, has a few drawbacks. First, the requirements for
linear elasticity are no longer valid near the crack tip
where atomic bonds clearly become non-linear and even-
tually break. Second, a continuum theory neglects the
discrete nature of the lattice. Thus, it is fully ignored
that fracture of materials is ultimately caused by bond
breaking processes on the atomic scale.
A way to understand the processes is to perform nu-
merical experiments, since experimental information on
this length scale is difficult to obtain. Molecular dynam-
ics studies have provided useful insight into crack prop-
agation in pure metals and simple intermetallic alloys,
whereas in complex metallic alloys the mechanisms are
not yet so clear. Atomistic studies show for example
that cracks remain stable in a region around the critical
stress intensity factor due to the discrete nature of the
lattice. This effect is called lattice trapping15. A further
consequence of the discrete lattice is that the fracture
behaviour in one and the same plane can depend on the
crack propagation direction16. Such observations cannot
be explained by a simple continuum mechanical descrip-
tion.
III. MODEL AND METHOD
A. Icosahedral binary model
In the numerical experiments we use a three-
dimensional model quasicrystal which has been proposed
by Henley and Elser17 as a structure model for the icosa-
hedral phase of (Al,Zn)Mg. This is the simplest possi-
ble model quasicrystal that is stabilised by pair poten-
tials. Furthermore it allows Burgers circuit analysis and
is a prototype of Bergman-type quasicrystals. As we
do not distinguish between Al and Zn atoms, we term
this decoration icosahedral binary model. It can be ob-
tained by decorating the structure elements of the three-
dimensional Penrose tiling, the oblate and the prolate
rhombohedra (see Fig. 1, top). Al and Zn atoms (A
atoms) are placed on the vertices and the midpoints of
the edges of the rhombohedra. Two Mg atoms (B atoms)
divide the long body diagonal of each prolate rhombohe-
dron in ratios τ :1:τ , where τ is the golden mean. Two
prolate and two oblate rhombohedra with a common ver-
tex form a rhombic dodecahedron18,19. To obtain the
icosahedral binary model, in these dodecahedra the atom
at the common vertex is removed and the four neigh-
bouring A atoms are transformed into B atoms. Finally
these atoms are shifted to the common vertex to divide
the edges of the corresponding rhombohedra in a ratio
of 1:τ . Fig. 1 (bottom left) shows the final decoration of
these dodecahedra, in which the B atoms form hexagonal
bipyramides. This modification increases the number of
Bergman-type clusters (see Fig. 1, bottom right) inher-
ent in the structure, leads to a higher stability with the
potentials used, and takes better into account the experi-
mentally observed stoichiometry of the quasicrystal. The
Bergman-type clusters may also be interpreted as build-
ing units of the quasicrystal and are the main feature of
the structure apart from the plane structure.
B. Molecular dynamics technique
As is often done in fracture simulations (see e.g. Abra-
ham20) we use simple Lennard-Jones potentials to model
the interactions. The following facts led to this choice:
First, the very few potentials available for quasicrys-
tals are unsuitable for fracture simulations: These pair
potentials stabilise the quasicrystal only for a fixed vol-
ume. With free surfaces atoms sometimes simply evap-
orate (e.g. with potentials based on those of Hafner et
al.21). Very long-ranged potentials with Friedel oscilla-
tions (e.g. those of Al-Lehyani et al.22) frequently display
such large cohesive energies that nearly no elastic defor-
mation is possible and instead intrinsic cracks develop.
Second, many simulations have proven that model po-
tentials are helping to understand the elementary pro-
FIG. 1: Tiles of the icosahedral binary model decorated with
two types of atoms and the Bergman-type cluster. Top: pro-
late rhombohedron (left) and oblate rhombohedron (right).
Bottom: rhombic dodecahedron (left) and the 45 atoms build-
ing the Bergman-type cluster (right) inherent in the model
quasicrystal. Small grey and large black spheres denote A
and B atoms respectively.
cesses in fracture (see e.g. Abraham20) and so are rea-
sonable when qualitative mechanisms are the centre of
interest. For quantitative results, a known limitation is
the neglect of non-local and many-body interactions.
Third, the Lennard-Jones potentials used23,24 keep the
model stable even under strongest mechanical deforma-
tions or irradiation (introduction of point defects) and
have been used in our group in many simulations of dis-
location motion25 or even shock waves26. The structure
is robust under a wide variation of the potential depths.
Very similar potentials have shown to stabilise the icosa-
hedral atomic structure in a simpler model27. It is also
known since the early fifties that Lennard-Jones poten-
tials favour icosahedral clusters28, indicating that these
potentials are useful for structures like icosahedral qua-
sicrystals.
The minima of the potentials for interactions between
atoms of the same type are set to �0, whereas the mini-
mum of the potential for the interactions between atoms
of different kind is set to 2�0. The conclusions drawn
from our simulations, however, remained essentially un-
changed if all binding energies are set equal. The shortest
distance between two A atoms is denoted r0. All masses
are set to m0. The time is then measured in units of
t0 = r0
m0/�0.
The molecular dynamics simulations were carried out
using the microcanonical ensemble with the program
code IMD29,30. It performs well on a large variety of
hardware, including single and dual processor worksta-
tions and massively parallel supercomputers.
First, we searched for the potential cleavage planes.
According to the Griffith criterion they should be planes
of low surface energy31. To identify these surfaces we re-
lax a specimen and split it into two parts. Subsequently,
the two regions are shifted apart rigidly. The surface en-
ergy is then calculated from the energy difference of the
artificially cleaved and the undisturbed specimen.
In contrast to simple periodic crystals, the atomistic
structure of the planes and therefore also the surface en-
ergy in quasicrystals varies strongly within the family
of planes perpendicular to a fixed axis. In the present
model, surfaces with low interface energy are located
perpendicular to two- and fivefold directions at certain
heights (Fig. 2). Perpendicular to other directions the
plane structure is less pronounced and the minimal sur-
face energies are higher.
Since we are interested in the morphology of fracture
surfaces we apply a sample form that allows us to fol-
low the dynamics of the running crack for a long time.
For this purpose, a strip geometry is used to model crack
propagation with constant energy release rate32. The
samples consist of about 4 to 5 million atoms, with di-
mensions of approximately 450r0×150r0×70r0. Periodic
boundary conditions33 are applied in the direction paral-
lel to the crack front. For the other directions, all atoms
in the outermost boundary layers of width 2.5r0 are held
fixed. An atomically sharp seed crack is inserted at a
plane of lowest surface energy, from one side to about
FIG. 2: Surface energy of cleavage planes perpendicular to
two- and fivefold axes.
one quarter of the strip length. The system is uniaxially
strained perpendicular to the crack plane up to the Grif-
fith load and is relaxed to obtain the displacement field of
a stable crack at zero temperature. Then a temperature
of about 10−4 of the melting temperature is applied34 to
the configurations with and without the relaxed crack.
From the resulting configurations we obtain an averaged
displacement field for this temperature. The crack now is
further loaded by linear scaling of this displacement field.
The response of the system then is monitored by molec-
ular dynamics techniques. The sound waves emitted by
the propagating crack (see Fig. 3 and online movie35)
are damped away outside of an elliptical stadium32 to
prevent reflections.
C. Visualisation
To study crack propagation on an atomic level the se-
lection and reduction of data is of crucial importance.
Due to the large number of atoms required for the study
of crack propagation in three-dimensional systems, it is
not feasible to write out the positions of all atoms during
the simulation, and even less to display all of them. To
obtain a first insight into crack propagation only atoms
FIG. 3: Kinetic energy density. Sound waves emitted by the
propagating crack and the elliptical region without damping
are clearly visible. See online movie35.
FIG. 4: Snapshot of a simulation with some 4 million atoms.
Only atoms with low coordination number are displayed. See
online movie36.
near the fracture surfaces are of interest. Whereas they
can be visualised in periodic crystals by plotting only
those atoms whose potential energy exceeds a certain
threshold, this technique is not applicable for quasicrys-
tals. Because of the largely varying environments the
potential energy varies significantly from atom to atom,
even for atoms of the same type in a defect-free sample.
A solution to this problem is to display only those
atoms whose coordination number is smaller than a cer-
tain threshold. This number is evaluated by counting
atoms within the nearest neighbour distance. Like the
potential energy, in quasicrystals the coordination num-
ber varies from atom to atom, but to a much smaller de-
gree. As atoms near defects have a significantly lower co-
ordination number, it becomes possible to visualise frac-
ture surfaces and dislocation cores. For the A atoms
the threshold for the coordination number is set to 12,
whereas for the B atoms it is set to 14. With this method,
the number of atoms to write out or to display can be
reduced by three orders of magnitude. Fig. 4 shows a
snapshot of a simulation with some 4 million atoms fil-
tered by this technique, which was also used in a movie
that is available online36.
As can be seen from Fig. 4 the fracture surfaces are
rough. Therefore to decide to which fracture surface an
atom belongs we apply the displacement vectors between
the initial configuration with the built-in seed crack and
the fractured sample. The morphology of the fracture
planes is then analysed via geometrical scanning of the
atoms forming the surfaces. The roughness can be visu-
alised by colour coding the height of the surface in a view
normal to the fracture surface.
To investigate the influence of the Bergman-type clus-
ters on cleavage they have to be displayed together with
the fracture surfaces intersecting them. This is done by
restoring the initial sample without crack at zero tem-
perature. The atoms forming the two sample parts are
taken back to their positions in this initial sample and
FIG. 5: Bergman-type cluster cut by a flat surface.
then scanned geometrically. In addition all atoms form-
ing clusters in the vicinity of this surfaces are known. By
displaying only these atoms and the scanned surface one
directly can see where and to which amount clusters are
cut. A problem of this kind of visualisation is shown in
Fig. 5, where a cluster is cut by a flat surface. When
looked-at from above it is obvious that only four atoms
are separated from the rest of the cluster. When looked-
at from below one could get the impression that the clus-
ter is heavily intersected. On a real fracture surface clus-
ters with centres above and below the crack surface are
present. Therefore both views of Fig. 5 are appearing
at the same time in a two-dimensional projection. Thus
such pictures are not very intuitive.
A way out of this dilemma is presented in Fig. 6. Clus-
ters with midpoints above the crack surface are displayed
together with the upper geometrically scanned fracture
surface only, the other clusters are shown together with
the lower fracture surface. As a result we get two pictures
with clusters cut by surfaces. Note that for a qualitative
and quantitative analysis always both pictures or sets of
data are needed.
FIG. 6: Visualisation of the clusters cut by the dynamic crack.
Midpoints of clusters are indicated as black dots, the clusters
are idealised as spheres. The crack propagated from the left to
the right. The upper and lower fracture surfaces are projected
onto an x− z plane.
IV. RESULTS
In this section the results of the numerical simulations
are presented. For practical purposes we define k as the
stress intensity factor K relative to the stress intensity
factor KG at the Griffith load:
k = K/KG.
The orientations of the samples are characterised by the
notation yx, where y is an axis perpendicular to the cleav-
age plane and x is an axis in the crack propagation di-
rection (see Fig. 6). An axis perpendicular to a fivefold
(5) and a twofold (2) axis is denoted pseudo-twofold (p2)
axis.
A. Crack velocities
Simulations were performed for different orientations
with loads in a range from k = 1.1 to k = 2 (see nota-
tions and loads in Fig. 7). Brittle fracture without any
crack tip plasticity is observed irrespective of the orien-
tation of the fracture plane. For loads below k = 1.2,
the crack propagates only a few atomic distances, and
then stops. Hence the energy needed for initiating crack
propagation is about 1.4 times the value predicted by the
Griffith criterion. Therefore, a simple global thermody-
namic description of fracture is not applicable. The min-
imal velocity for brittle crack propagation is about 10%
of the shear wave velocity37 vs. For loads k ≥ 1.2 the
velocity increases monotonically with the applied load.
The crack velocities are in a range of 10-45% of vs (see
Fig. 7). At high loads the crack velocities on fivefold
cleavage planes show higher average velocities than on
the other planes. Velocities for the two different crack
propagation directions on the fivefold planes differ sig-
nificantly at intermediate loads (k = 1.3). This coincides
with ledges that are produced in the fracture surface (see
section IV B and Fig. 8, bottom).
B. Fracture surfaces
To analyse the morphology of the fracture surfaces,
they are displayed as described in section III C. In Fig. 8
the average height is shown in grey, heights above +2r0
are shown in white and heights below −2r0 are shown
in black. The crack propagation direction is from the
left to the right. The initial fracture surface is flat, as
can be seen from the homogeneous regions on the left.
The surfaces resulting from the propagation of the crack,
however, show pronounced patterns of regions with dif-
ferent heights. From the observation of the fracture sur-
faces it is already evident that they are rough and that
the peak-to-valley roughness is of the order of the di-
ameter of the clusters. The peak-to-valley roughness
and the root-mean-square roughness of the height pro-
files both increase38 for higher loads for surfaces without
FIG. 7: Average crack velocities for different loads and orien-
tations.
ledge formation. For fracture surfaces perpendicular to
twofold and fivefold axes the crack fluctuates about a
constant height (in the areas without ledges). In con-
trast, a crack inserted perpendicular to a pseudo-twofold
direction seems to deviate from this plane39.
C. Anisotropy
As can be seen already from the fracture surfaces per-
pendicular to a fivefold axis in Fig. 8 crack surfaces for the
same cleavage plane differ significantly for different in-
plane crack propagation directions. For the orientation
5p2 ledges are produced, while no ledges form for the ori-
entation 52. By visual inspection of the fracture surfaces
(see Fig. 8) and evaluation of height-height-correlation
functions23 it becomes evident that for every orientation
there are distinct angles in the height profile, which show
pronounced height correlations that correspond to mark-
ings in the fracture surfaces or to ledges. These angles
are given in Table I.
D. Clusters
In Fig. 9 the clusters cut by the fracture surfaces are
presented as described in section III C and Fig. 6. It is
obvious from Fig. 9 that the dynamic crack does not per-
fectly circumvent the clusters, but intersects them par-
tially (right side of Fig. 9). These intersections, however,
TABLE I: Angles observed in the height profiles of the frac-
ture surfaces. Angles measured clockwise from the crack prop-
agation direction get a negative sign.
orientation 22 25 52 5p2 p25
angles 0◦,±32◦ 0◦, +32◦,−58◦ 0◦,±36◦ ±18◦ 0◦,±90◦
FIG. 8: Height profiles of sections of the simulated fracture surfaces. Load: k = 1.3, orientation 22 (top), 52 (middle), 5p2
(bottom). The height increases from black (≤ −2r0) to white (≥ +2r0). The scanning sphere has the same size as an atom of
type B.
FIG. 9: Clusters cut by the fracture surfaces. The visualisation technique is described in Sec. III C and Fig. 6. Load: k = 1.3,
orientation: 52.
FIG. 10: Surface energy and density of cluster centres for the
orientation 22. The corresponding fracture surface is shown
in Fig. 8, top.
are much less frequent than for the flat seed cracks (left
side of Fig. 9). More detailed analyses for different ori-
entations validate this statement. For the orientations
perpendicular to twofold and fivefold axes at k = 1.3 the
ratio of clusters cut by the dynamic crack to clusters in-
tersected by flat cuts is approximately 0.6. Additionally
the absolute value of clusters cut by the crack for the
fivefold surfaces is lower than for the twofold surfaces.
Fig. 10 and Fig. 11 display bottom up: The density of
the cluster centres, the surface energy, a cluster in the
proper length scale, the grey coding of the heights in
Fig. 8 (top, middle), and the position of the seed crack
(dashed vertical line). For the twofold fracture surface
the low energy seed crack is located between two peaks
in the cluster density, whereas for the fivefold surface
this seed crack is situated close to a peak of this density.
It is evident from the figures that it is not possible to
circumvent all clusters by a planar cut.
The grey coding is adjusted to the average height of the
fracture surfaces. It is therefore evident from Fig. 11 that
the crack deviates for the orientation 52 from the low en-
ergy cleavage plane of the seed crack to a parallel plane,
reducing the number of cluster intersections dramatically
(see also Fig. 9 and Fig. 8, middle). Samples cut flat at
the new height show slightly higher surface energy (see
also Fig. 11). However, for low loads and low roughness
the actual fracture surfaces of the dynamic cracks have
about 5-15% higher energies than those of the low energy
seed cracks. To assure that the dynamic crack is depart-
ing from the initial plane not in a random manner the
seed crack was built in at the position colour coded as
medium grey in Fig. 11. The resulting fracture surface
had a similar roughness but the crack did not change to
a parallel plane.
FIG. 11: Surface energy and density of cluster centres for the
orientation 52. The corresponding fracture surface is shown
in Fig. 8, middle.
V. DISCUSSION
Taken together the results of our simulations presented
above indicate that the distribution of the clusters is
crucial for the fracture behaviour: First, circumventing
the clusters or intersecting them disturbs the propagat-
ing crack and leads to additional radiation. This man-
ifests itself in the crack speed. Dynamic cracks prop-
agating in fivefold planes with few cluster intersections
are faster than those in twofold planes, where the abso-
lute number of cluster intersections is higher. Also the
generation of ledges for low loads in the orientation 5p2
costs energy and therefore slows down the crack even
further. Second, circumvention of the clusters leads to
characteristic height-variations. Therefore the peak-to-
valley roughness of the fracture surfaces is of the order
of the diameter of the clusters. Third, the observed pat-
terns in the fracture surfaces correspond to lines along
which the clusters are located. The associated angles are
given by the icosahedral symmetry of the sample, namely
0◦, 18◦, 31.72◦, 36◦, 58.28◦, and 90◦ (see Table I). Ledges
seem to be produced only for the smallest angles mea-
sured from the crack propagation direction. Fourth, less
clusters are intersected by the fracture surfaces than by
the flat seed cracks. Fifth, a seed crack at a low energy
cleavage plane deviates to a parallel plane to reduce the
number of cluster intersections in spite of the higher en-
ergy required to form the fracture surfaces. In contrast,
a crack built-in at this new position does not show such
a deviation.
Another observation of the simulations is that the
plane structure of the quasicrystal also influences frac-
ture. The fracture surfaces that are located perpendicu-
lar to the twofold and fivefold symmetry axes show con-
stant average heights.
The three-dimensional quasicrystals give perfect cleav-
age fracture with no indication of any dislocation activ-
ity. This is in contrast to results on two-dimensional
decagonal quasicrystals, where a dislocation enhanced
fracture mechanism has been observed40. However a
corresponding three-dimensional decagonal quasicrystal
would have a periodic direction with a straight dislo-
cation line. In the simulations presented here this di-
rection is also quasiperiodic. As the clusters have a
strong influence on fracture they also may bend and
hinder dislocation lines. So dislocation emission in the
three-dimensional icosahedral quasicrystal modelled here
should be less likely than in the two-dimensional decago-
nal model quasicrystal. Very high stresses are indeed
needed to move dislocations in our model quasicrystal in
molecular dynamics simulations25.
There are also indications for the stability of the
clusters from the electronic structure of quasicrystals:
First, experiments and ab-initio calculations show that
directional bonding may be present within clusters of
quasicrystals41,42. Second, the electrons may addition-
ally stabilise the clusters because of their hierarchical
structure4. Therefore they should be even more stable
than we have modelled them with simple pair potentials.
With this evidence the results concerning the clusters
seem reasonable and should even underestimate their sta-
bility.
So far fracture experiments in ultrahigh vacuum have
only been performed on icosahedral Al-Mn-Pd, which
has a more complicated atomic structure than the icosa-
hedral binary model. Additionally the clusters are not
Bergman-type. Therefore we cannot expect to represent
this structure on an atomic level, when comparing exper-
iments in Al-Mn-Pd with our simulations. Nevertheless,
the size of the clusters, the icosahedral symmetry, and a
distinct plane structure are common features and qual-
itative aspects should be reproduced well, namely the
size and shape of the patterns and the appearance of dis-
tinct angles on the fracture surfaces. This is indeed the
case, as can be seen in Fig. 12. There a geometrically
scanned fracture surface generated in our simulations is
confronted to an STM-image of Ebert et al.7,43 at the
same length scale. As we were able to correlate the ob-
served structures to the clusters in our model, the simi-
larities corroborate the assumption that the clusters are
responsible for the globular structures observed in ex-
periment. More detailed comparisons to fracture experi-
ments in icosahedral Al-Zn-Mg-type quasicrystals would
be desirable, but such comparative data is currently not
available.
We also performed fracture simulations in a C15 Laves
phase with the same Lennard-Jones potentials to give fur-
ther evidence, that the roughness is correlated to the clus-
ters. This C15 structure is built up from deformed pro-
late rhombohedra, one of the major tiles of the quasicrys-
tal (see Fig. 1 and Sec. III A). However, no Bergman-
type clusters are present. When colour coded like the
quasicrystal, fracture surfaces of the C15 cracks lack any
roughness at low loads and cleave smoothly on high sym-
metry planes44.
FIG. 12: Fracture surfaces perpendicular to twofold axes. In
the left picture (simulation) the atomically sharp seed crack
can be seen on the top, whereas below this area the simulated
fracture surface appears. The orientation of the sample is 22,
the load was k = 1.3. The surface was geometrically scanned
with an atom of type B. Thus atomic resolution is achieved.
The right picture (experiment, adopted from Ebert et al.43)
shows an STM-image of icosahedral Al-Pd-Mn cleaved in ul-
trahigh vacuum. As 20r0 is approximately 5 nm the surfaces
are displayed at the same length scale.
When a crack traverses a solid, it leaves a typical non-
equlibrium surface. Thus, up to this point, we have not
dealt with any equilibrium surface. Even the flat Griffith
cuts we introduced to perform the numerical experiments
are no explicit equilibrium surfaces31 by definition. It is
known that i-Al-Pd-Mn surfaces sputtered and annealed
up to about 900 K are rough with cluster-like protru-
sions. Fivefold surfaces annealed at higher temperatures
exhibit flat terraces9,10,11, which are believed to be bulk
terminated45. As fivefold surfaces are not as rough as
twofold surfaces, they are often studied in experiment.
Adsorption on these surfaces, nevertheless, is often very
site-sensitive. A recent review on quasicrystal surfaces
was given by McGrath et al.12. All of these observations
are consistent with our simulations of non-equilibrium
surfaces. Fivefold surfaces experience less roughness than
twofold surfaces (also in agreement with the experiments
of Ebert et al.7,43; see Fig. 8), where more clusters are
intersected. Additionally, since clusters are cut in the
simulations, their binding energy is not so large as to
term them supermolecules. Thus annealing at very high
temperatures can favour flat surfaces11,12. Nevertheless,
selective adsorption may then be related to completion
of clusters.
VI. CONCLUSIONS
We have simulated crack propagation in an icosahe-
dral model quasicrystal. Brittle fracture without any
signature of dislocation emission is observed. The frac-
ture surfaces are rough on the scale of the clusters and
show constant average heights for orientations perpen-
dicular to twofold and fivefold axes. Thus both the plane
structure and the clusters play an important role in frac-
ture. The influence of the clusters is also seen in the
average crack velocities for different orientations, the ob-
served patterns in the fracture surfaces, the anisotropy
with respect to the in-plane propagation direction, and
the smaller amount of clusters cut by the propagating
crack than by planar cuts. The clusters, too, are a rea-
son why the positions of the cleavage planes cannot be
predicted by a simple energy criterion. Since partial clus-
ter intersections occur, the binding energy of the clusters
is not so large as to term them supermolecules. Never-
theless our observations clearly show that they are not
only structural units but physical entities.
Acknowledgments
Financial support from the Deutsche Forschungsge-
meinschaft under contract numbers TR 154/13 and TR
154/20-1 is gratefully acknowledged.
1 C. L. Henley, Phys. Rev. B 43 (1), 993 (1991).
2 V. Elser, Phil. Mag. B 73 (4), 641 (1996).
3 P. J. Steinhardt, H.-C. Jeong, K. Saitoh, M. Tanaka, E.
Abe, and A. P. Tsai, Nature 396 (6706), 55 (1998).
4 C. Janot and M. de Boissieu, Phys. Rev. Lett. 72 (11),
1674 (1994).
5 C. Janot, J. Phys.: Condens. Matter 9, 1493 (1997).
6 H.-C. Jeong and P. J. Steinhardt, Phys. Rev. Lett. 73 (14),
1943 (1994).
7 Ph. Ebert, M. Feuerbacher, N. Tamura, M. Wollgarten,
and K. Urban, Phys. Rev. Lett. 77 (18), 3827 (1996).
8 Ph. Ebert, M. Yurechko, F. Kluge, K. Horn, and K. Ur-
ban, Quasicrystals – Structure and Physical Properties,
WILEY-VCH, ed. H.-R. Trebin, 6.1, 572 (2003).
9 T. M. Schaub, D. E. Bürgler, H.-J. Güntherodt, and J. B.
Suck, Phys. Rev. Lett. 73 (9), 1255 (1994).
10 M. Gierer, M. A. Van Hove, A. I. Goldman, Z. Shen, S.-L.
Chang, C. J. Jenks, C.-M. Zhang, and P. A. Thiel, Phys.
Rev. Lett. 78 (3), 467 (1997).
11 G. Cappello, J. Chevrier, F. Schmithüsen, A. Stierle, V.
Formoso, F. Comin, M. de Boissieu, M. Boudard, T. Lo-
grasso, C. Jenks, and D. Delaney, Phys. Rev. B 65, 245405
(2002).
12 R. McGrath, J. Ledieu, E. J. Cox, and R. D. Diehl, J.
Phys.: Condens. Matter 14 (4), R119 (2002).
13 G. Kasner, Z. Papadopolos, P. Kramer, and D. E.
Bürgler, Quasicrystals – Structure and Physical Properties,
WILEY-VCH, ed. H.-R. Trebin, 2.4, 123 (2003).
14 A. A. Griffith, Philos. Trans. R. Soc. Lond. Ser. A 221, 163
(1921).
15 R. Thomson, C. Hsieh, and V. Rana, J. Appl. Phys. 42 (8),
3154 (1971).
16 P. Gumbsch and R. M. Cannon, Mat. Res. Soc. Bull.
25 (5), 15 (2000).
17 C. L. Henley and V. Elser, Phil. Mag. B 53 (3), L59 (1986).
18 The icosahedral binary model is described and named BIB
in J. Roth, Eur. Phys. J. B 15 (1), 7 (2000).
19 In rings of oblate rhombohedra the number of clusters is
maximised. No overlapping rhombic dodecahedra are gen-
erated. Remaining oblate rhombohedra stay unchanged.
20 F. F. Abraham, Advances in Physics 52 (8), 727 (2003).
21 J. Hafner and M. Krajč́ı, Europhys. Lett. 13 (4), 335
(1990).
22 I. Al-Lehyani, M. Widom, Y. Wang, N. Moghadam, G.
M. Stocks, and J. A. Moriarty, Phys. Rev. B 64, 075109
(2001).
23 F. Rösch, diploma thesis, Universität Stuttgart (2003),
http://elib.uni-stuttgart.de/opus/volltexte/2004/1899/
24 F. Rösch, Ch. Rudhart, P. Gumbsch, and H.-R. Trebin,
Mat. Res. Soc. Symp. Proc. 805, LL9.3.1 (2004).
25 G. D. Schaaf, J. Roth, and H.-R. Trebin, Phil. Mag.
83 (21), 2449 (2003).
26 J. Roth, Phys. Rev. B 71, 064102 (2005).
27 J. W. Roth, R. Schilling, and H.-R. Trebin, Phys. Rev. B
51 (22), 15833 (1995).
28 F. C. Frank, Proc. Roy. Soc. Lond. Ser. A 215 (1120), 43
(1952).
29 J. Stadler, R. Mikulla, and H.-R. Trebin, Int. J. Mod.
Phys. C 8 (5), 1131 (1997).
30 IMD, the ITAP Molecular Dynamics Program.
http://www.itap.physik.uni-stuttgart.de/˜imd
31 To simulate fracture with molecular dynamics it is nec-
essary to determine the optimal place for a low energy
seed crack and to find out the critical strain for setting
the crack into motion. This is always done by determining
the energies of flat cuts. Note, that these do not necessar-
ily represent equilibrium surfaces which may be bent and
curved.
32 P. Gumbsch, S. J. Zhou, and B. L. Holian, Phys. Rev. B
55 (6), 3445 (1997).
33 As quasicrystals show non-periodic translational order we
generate periodic approximants to apply periodic bound-
ary conditions. The 4 to 5 million atoms of our sample
then form the unit cell. Because of this high number of
atoms the configuration should mechanically behave like
the perfect quasicrystal.
34 Here we use equation (2) of Gumbsch et al.32 with a “filled
stadion damping” f ≡ 1.
35 Cracks propagate by breaking bonds between atoms. In
this process, energy is dissipated in the form of acoustic
waves, as clearly visible in the online movie. The kinetic
energy density is colour coded. Regions of high kinetic en-
ergy density are shown in red whereas blue indicates areas
of lower kinetic energy density.
36 Movie (available online) of a cracking model quasicrystal
at a high load (5p2, k = 1.6). By displaying only atoms
with a low coordination number fracture surfaces are visu-
alised. Blue (black in Fig. 1) and red (grey in Fig. 1) balls
correspond to the atoms in the icosahedral binary model.
Atoms on the border, which are not allowed to move, are
coloured in yellow and green. On the left the atomically
sharp seed crack can be seen. Due to the high deformation
near the crack tip, blue atoms loose their neighbours and
show up for short instances.
37 Concerning only the pure phonon term of the elastic en-
ergy in linear elastic continuum theory, the icosahedral
http://elib.uni-stuttgart.de/opus/volltexte/2004/1899/
http://www.itap.physik.uni-stuttgart.de/~imd
quasicrystal behaves like an isotropic solid.
38 The height profiles were determined via geometrical scan-
ning of the fracture surfaces with an atom of type B at
equidistant points separated by 0.2r0. For these meshes
and e.g. for the orientation 52 the root-mean-square rough-
ness increases from 0.39r0 at a load of k = 1.3 (see Fig. 8,
middle) to 0.54r0 for k = 1.6. The peak-to-valley rough-
ness increases likewise from 3.7r0 to 5.1r0. Although the
exact values may depend on the stepwidth (resolution),
scanning sphere size, and scanned area a general tendency
for an increased roughness for increased loads is evident.
39 The question whether this crack finally chooses a plane
perpendicular to a twofold or fivefold axis cannot be an-
swered, as simulations with the required sample-size would
exceed present computer capabilities.
40 R. Mikulla, J. Stadler, F. Krul, H.-R. Trebin, and P.
Gumbsch, Phys. Rev. Lett. 81 (15), 3163 (1998).
41 Only clusters in the bulk give reliable information on their
relative stability. In ab-initio studies, where the number of
atoms is limited, one therefore is restricted to small ap-
proximants.
42 K. Kirihara, T. Nagata, K. Kimura, K. Kato, M. Takata,
E. Nishibori, and M. Sakata, Phys. Rev. B 68, 014205
(2003).
43 Ph. Ebert, F. Yue, and K. Urban, Phys. Rev. B 57 (5),
2821 (1998).
44 F. Rösch, P. Gumbsch, and H.-R. Trebin, unpublished.
45 Z. Papadopolos, G. Kasner, J. Ledieu, E. J. Cox, N. V.
Richardson, Q. Chen, R. D. Diehl, T. A. Lograsso, A. R.
Ross, and R. McGrath, Phys. Rev. B 66, 184207 (2002).
	INTRODUCTION
	FRACTURE
	MODEL AND METHOD
	Icosahedral binary model
	Molecular dynamics technique
	Visualisation
	RESULTS
	Crack velocities
	Fracture surfaces
	Anisotropy
	Clusters
	DISCUSSION
	CONCLUSIONS
	Acknowledgments
	References
