Structure factors of harmonic and anharmonic Fibonacci chains by
  molecular dynamics simulations
Structure factors of harmonic and anharmonic Fibonacci chains
by molecular dynamics simulations
Michael Engel,∗ Steffen Sonntag, Hansjörg Lipp, and Hans-Rainer Trebin
Institut für Theoretische und Angewandte Physik,
Universität Stuttgart, Pfaffenwaldring 57, 70550 Stuttgart, Germany
The dynamics of quasicrystals is characterized by the existence of phason excitations in addition
to the usual phonon modes. In order to investigate their interplay on an elementary level we resort
to various one-dimensional model systems. The main observables are the static, the incoherent,
and the coherent structure factor, which are extracted from molecular dynamics simulations. For
the validation of the algorithms, results for the harmonic periodic chain are presented. We then
study the Fibonacci chain with harmonic and anharmonic interaction potentials. In the dynamic
Fibonacci chain neighboring atoms interact by double-well potentials allowing for phason flips. The
difference between the structure factors of the dynamic and the harmonic Fibonacci chain lies in
the temperature dependence of the phonon line width. If a bias is introduced in the well depth,
dispersionless optic phonon bands split off.
PACS numbers: 63.20.Ry, 61.44.Br, 02.70.Ns
Keywords: Phonons; Phason flip; Quasicrystal.
I. INTRODUCTION
A. Motivation
Quasicrystals are long-range ordered materials lacking
translational symmetry1. Their diffraction patterns ex-
hibit a dense set of sharp Bragg reflections, that can
be indexed by an integer linear combination of a finite
number of basis vectors which is larger than the dimen-
sion of space. As a consequence the atomic configuration
of quasicrystals is described with reference to a higher-
dimensional analog of a periodic lattice. Elementary dy-
namic excitations within this ‘hyperspace’ description are
phonons2 and phasons3.
Phasons involve rearrangements of the structure by
atomic jumps over short distances, denoted ‘phason
flips’. They are connected with many physical proper-
ties of quasicrystals as for example elastic deformations4,
dislocations5,6, diffusion7,8, and phase transformations9.
Recently, indications for phason flips10,11 have been ob-
served by in situ transmission electron microscopy. A
coherent set of phason flips may form a static phason
field, e.g. during a phase transformation or in the neigh-
borhood of a dislocation12.
By investigating the dynamics of quasicrystals one
can find out the influence of the quasiperiodicity on the
phonon spectrum2 and one may gain a deeper under-
standing of phason flips13. Both points can be studied
in x-ray or neutron diffraction experiments by measur-
ing the response of the system in frequency (ω) and
momentum (q) space. Depending on the experimental
setup, different functions can be obtained from the scat-
tering experiments: (1) The static structure factor S(q)
is the usual – not energy resolved – diffraction image,
measured with either x-rays or neutrons. It is used for
the determination of the atomic structure. (2) The co-
herent structure factor14 S(q, ω) is studied via coherent
inelastic neutron scattering15 or alternatively via inelas-
tic x-ray scattering16. It allows the determination of the
phonon dispersion relations. The experiments on icosa-
hedral quasicrystals show well defined acoustic phonon
modes at small wave-vectors17 and dispersionless broad
optic bands at larger wave vectors18. The cross-over be-
tween the two regions is very sharp. (3) The incoherent
structure factor14 Si(q, ω) can be measured in quasielas-
tic neutron scattering. Neutrons are exclusively used
here, due to the necessity of a high energy resolution.
The technique also allows the investigation of phason
flips. In a series of experiments Coddens et al.13,19,20,21
have found an anomalous q-dependence of the quasielas-
tic signal in icosahedral quasicrystals. They interpreted
it as correlated simultaneous jumps of several atoms.21
Up to now various calculations of the coherent struc-
ture factor of quasiperiodic model systems have been
published, see Ref. 2. Amongst them are the perfect one-
dimensional Fibonacci chain22, by static phason fields
disordered Fibonacci chains23, and three-dimensional
tilings24. In these studies the dynamical matrix is diago-
nalized, which is a purely analytic method and yields the
phonon dispersion relations only. The results are highly
structured excitation spectra with a hierarchical system
of gaps22. The influence of anharmonicities, however, es-
pecially the dynamics of phason flips has not been taken
into account.
This ‘missing link’ markes the starting point of our
study. Here we present calculations of the structure fac-
tors of special one-dimensional quasiperiodic model sys-
tems by use of molecular dynamics (MD) simulations
with either harmonic potentials or potentials that allow
for phason flips. Although structure factors play such a
central role in the dynamics of solids, not much seems
to be known about their exact forms for one-dimensional
chains. Even for the simple harmonic chain only few ar-
ticles exist25,26,27,28.
FIG. 1: Double-well potential V (x) = x4−2x2 of the dynamic
Fibonacci chain. The equilibrium distances are S and L.
B. Model systems
As a simple one-dimensional model for a quasiperiodic
system we consider the Fibonacci chain. It consists of
particles arranged at two different distances: large ones
(L) and small ones (S). The length ratio L/S equals
the number of the golden mean τ = 1
5 + 1). The
sequence of the distances is created recursively by the
mapping {L, S} 7→ {LS,L} with starting condition L.
For example, after four iterations the resulting sequence
is LSLLSLSL.
We want to study chains consisting of identical parti-
cles with nearest-neighbor interactions. The Hamiltonian
has the form
+ V (xj − xj+1 − aj), (1)
where xj and pj are position and momentum of the jth
particle. The dynamic Fibonacci chain (DFC) is defined
by the choices aj ≡ a0 = τ3 for the equilibrium distances
and V (x) = x4 − 2x2 for the interaction potential, re-
spectively. The latter forms a double-well potential with
minima at ±1 and a potential hill of height ∆E = 1
as shown in Fig. 1. Because neighboring particles sit in
either of the potential minima, two nearest-neighbor dis-
tances L = a0 + 1 and S = a0 − 1 are possible. They
fulfil the constraint L/S = τ of the Fibonacci chain.
The DFC shows two types of elementary excitations:
Phonon vibrations in the local minima and phason flips
that interchange the particle distances L and S. At low
temperatures only phonons are excited, phason flips have
to be activated thermally. With its neighbors at rest the
activation energy of a particle for a phason flip is 2∆E.
This value is a result of the perfect superposition of the
potential hills of both neigbors. The value is lowered
when the neighbors assist by stepping simultaneously to
the inside or outside during the phason flip thus creating
a non-perfect superposition of the potential hills. Since
the particle distances L and S are energetically degen-
erate, the total equilibrium potential energy is invariant
under a phason flip.
The occurrence of phason flips makes nonlinearity an
intrinsic feature of the DFC and an analytical treatment
of the dynamics impossible. To understand the influence
of the nonlinearity, we study four model systems with
increasing complexity concerning their dynamical behav-
• Harmonic periodic chain (HPC):
VHPC(x) = 4x2 and aj = a = 2
• Harmonic Fibonacci chain (HFC):
VHFC(x) = 4x2 and aj = L or S according to the
Fibonacci sequence.
• Dynamic Fibonacci chain (DFC):
VDFC(x) = x4 − 2x2 and aj = τ3.
• Asymmetric Fibonacci chain (AFC):
VAFC(x) = VDFC +χ(x2−1)2(�x+x2/2−1/2) and
aj = τ3 with χ ∈ [0, 1] and � = ±1.
The potentials of the HPC, HFC, and DFC are chosen to
be identical in the harmonic approximation around the
equilibrium separation. The average particle distance a is
the same for all four systems. In the case of the Fibonacci
chain the occurrence probabilities for L and S are given
by τ−1 and τ−2, hence a = 2
For the solution of the equations of motions we use
a special molecular dynamics (MD) code. The code is
introduced in Sec. II together with a short theoretical
background. The simplest system is, of course, the HPC.
Exact solutions for the equation of motion exist as a su-
perposition of plane waves. We study the dynamics of the
HPC in Sec. III as a reference system. The HFC con-
sists of particles arranged on the Fibonacci chain with
distances L and S interacting with the same harmonic
potentials as the HPC, see Sec. IV. The DFC will then
be studied in Sec. V. In the case of the AFC the particles
in the two potential minima have different eigenfrequen-
cies. The parameters χ and � determine the degree of
asymmetry. For more details we refer to Sec. VI. We
finish with a discussion and conclusion in Sec. VII.
II. STRUCTURE FACTORS FROM
MOLECULAR DYNAMICS
A. Definition of the structure factors
We write the particle number density of the chain with
N particles as a sum of delta functions positioned along
the particle trajectories xl(t), n(x, t) =
l=1 δ(x−xl(t)).
The time dependent density-density correlation function
and the density-density autocorrelation function are de-
fined as
G(x, t) =
〈n(x′, t)n(x+ x′, 0)〉 dx′
〈δ(x− xj(t) + xl(0))〉 , (2a)
Ga(x, t) =
〈δ(x− xl(t) + xl(0))〉 , (2b)
where the brackets denote the thermal average.29 The
coherent and incoherent structure factor are the space-
time Fourier transforms30,
S(q, ω) =
e−iωt
e−iqxj(t)eiqxl(0)
Si(q, ω) =
e−iωt
e−iqxl(t)eiqxl(0)
Both functions are symmetric about q = 0 and ω = 0.
The static structure factor is the integral of the coher-
ent structure factor, S(q) =
S(q, ω) dω, i.e. the Fourier
transform of G(x, 0),
S(q) =
e−iqxj(0)eiqxl(0)
. (4)
B. Molecular dynamics simulations
For further calculations the particle trajectories are re-
quired as solutions of the equations of motion. Since in
the case of the anharmonic chains only numerical solu-
tions exist, we use a simple MD code. Initially the par-
ticles are placed on the equilibrium positions of a finite
chain of length L with periodic boundary conditions. The
velocities are initialized according to a Gaussian distribu-
tion. Its width determines the total energy and thus the
temperature of the system. The equations of motion are
integrated by a Verlet algorithm running for a simulation
time Tsim. After starting the simulation, the dynamics is
not controlled by a thermostat or in any other way.
For the direct numerical calculation of the Eqs. (3) we
must compute a fourfold sum: two sums over the par-
ticle number N and two over the time Tsim, one sum
for the Fourier transform and one for the time average.
Note, that by assuming ergodicity the thermal average
〈 〉 can be replaced by a time average 1
dt and ad-
ditionally by an average over several independent MD
runs. For the sake of clarity the averaging over the MD
runs is suppressed in the following notation. We intro-
duce a more compact notation by defining the functions
fl(q, t) = eiqxl(t). Let us assume tentatively that these
functions are periodic in time with period Tsim and in
space with period L. Then the Eqs. (3) and (4) are
greatly simplified to
S(q, ω) =
2πNTsim
∥∥∥∥∥
e−iωt
fl(q, t) dt
∥∥∥∥∥
, (5a)
Si(q, ω) =
2πNTsim
e−iωtfl(q, t) dt
∥∥∥∥2 , (5b)
S(q) =
NTsim
∫ ∥∥∥∥∥∑
fl(q, t)
∥∥∥∥∥
dt. (6)
The equations differ in the order of the absolute square
and the particle sum. Since only two sums are left, an
efficient numerical computation of the structure factors
is possible. Furthermore a fast Fourier transform is used
for the time integrals in Eqs. (5).
It is left to discuss the periodicity conditions. The spa-
tial periodicity follows from the periodic boundaries used
in the simulation. Therefore the chain acts like a ring and
the excitations can go round during the simulation. To
avoid such a behavior that would lead to unwanted cor-
relations, we limit the maximum simulation time by the
quotient of the length of the chain and the sound velocity
cs to Tmax = L/cs. For all the model systems HPC, HFC,
and DFC the sound velocity is the same, cs = a
8, and
Tmax = N/
8. We use Tsim = Tmax.
The functions fl(q, t) are in general not periodic in
time. There is no reason, why the particles should be
at the same positions at the end of the simulation as
at the beginning. To avoid this problem, we multiply
fl(q, t) with a window function w(t) to enforce an arti-
ficial periodicity. The function w(t) has to decrease fast
enough – both in direct as in Fourier space – towards the
boundaries of its domains. We use a normalized broad
Gaussian function. Its width is chosen as large as pos-
sible with the constraint that the Gaussian has decayed
to a small enough value at the interval boundaries. The
effect of the Gaussian is a smoothing of the structure fac-
tors by convolution with a narrow Gaussian. The exact
value of the width has no influence on the results.
III. HARMONIC PERIODIC CHAIN AS
REFERENCE SYSTEM
A. Analytic calculations
The harmonic periodic chain (HPC) is used as a ref-
erence system to test our algorithms since its equa-
tions of motions can be solved analytically. If we put
xl(t) = ul(t)+ la, then the ul(t) are expressed by a linear
combination of normal modes. The wave vector q and
the frequency ω are related according to the dispersion
relation ω(q) = 2ω0| sin(qa/2)|. Here, w0 is the eigen-
frequency of a single particle. In the case of the model
systems HPC, HFC, and DFC we have ω0 =
For the HPC the thermal averages in the structure
factors, Eqs. (3) can be calculated to be
S(q, ω) =
e−iωt
e−iqal exp
q2σ2l (t)
Si(q, ω) =
e−iωt exp
q2σ20(t)
dt. (7b)
where we used from the literature25
σ2l (t) =
[ul(t)− u0(0)]2
∫ 2ω0t
J2l(s)(2ω0t− s) ds
. (8)
Jn(s) is the Bessel function of the first kind of order
n. The only external parameter in these functions is
the temperature T . The particle sum and the Fourier
transform have to be evaluated numerically. Due to the
translation invariance of the HPC the double sum of the
Eqs. (3) is reduced to a single sum. The static structure
factor of the HPC can also be calculated analytically26
S(q) =
sinh(q2σ2/2)
cosh(q2σ2/2)− cos(qa)
, (9)
where σ2 = kBT/ω20 .
Let us take a closer look at the incoherent struc-
ture factor Si(q, ω). In the limit of small T the term
exp(− 1
q2σ20(t)) decays slowly with t and we substi-
tute σ0(t) with its approximation for large t: σ0(t) =
|t|kBT/ω0 for |t| → ∞. This leads to a Lorentzian peak
Si(q, ω) =
Γ2 + ω2
, Γ =
q2kBT
In the limit of large T the term σ0(t) is approximated for
small t: σ0(t) = t2kBT for |t| � ω−10 . Hence there is a
Gaussian peak
Si(q, ω) =
, γ =
q2kBT . (11)
The transition temperature between these two limiting
cases is kBT = 4ω20/q
The Fourier transform of the Eqs. (7) yields the corre-
lation functions
G(x, t) =
σl(t)
(x+ la)2
2σ2l (t)
(12a)
Ga(x, t) =
σ0(t)
2σ20(t)
. (12b)
The function G(x, t) is shown in Fig. 2. It consists of a
sum of Gaussians centered at the equilibrium positions
of the particles. The width of the Gaussians increases
with temperature, as well as with time and in space:
limx,t→∞G(x, t) = ρ and ρ = 1/a. This means that
there is no longe-range order. Indeed, for the particle
number density we have 〈n(x)〉 = ρ, which is uniform as
in liquids25. The autocorrelation function Ga(x, t) corre-
sponds to the center peak at x = 0.
B. Simulation results
In the case of the incoherent structure factor the vari-
ables T and q only appear as combination q2kBT in
FIG. 2: Density-density correlation function G(x, t) of the
HPC for kBT = 0.5. The Gaussians are centered at integer
multiples of x = a ≈ 4.47.
FIG. 3: Incoherent structure factor Si(q, ω) of the HPC for
q = π/a, and different temperatures. The symbols mark the
data from MD simulations with N = 1000 particles, the lines
the result from the analytical fomula, Eq. (7b). The peaks
and edges are at integer multiples of ω = 2ω0 = 2
Eqs. (7) and (8). Therefore it suffices to examine Si at a
fixed wave vector for different temperatures. We choose
q = π/a arbitrarily. The results from MD simulation
and the numerical integration of the analytical formula,
Eq. (7b) are shown in Fig. 3 for temperatures ranging
from 0.01 to 100.0. There is a maximum at ω = 0, called
the quasielastic peak. In the low temperature regime the
maximum has a Lorentzian shape. Furthermore a one-
phonon peak, a two-phonon edge, a three-phonon edge,
etc. are found at 2ω0 ≈ 5.7, 4ω0 ≈ 11.3, 6ω0 ≈ 17.0.31
Note that the multi-phonon contributions rapidly decay
for larger ω (logarithmic scale). At higher temperatures
the curve smoothes and approaches a Gaussian profile.
The coherent structure factor S(q, ω) is shown in Fig.
4. For a low temperature of kBT = 0.01 a one-phonon
branch, a two-phonon branch, and very weakly a three-
phonon branch are observed. The one-phonon branch
has a Lorentzian line shape and follows the phonon-
dispersion relation. At higher q-values the branch broad-
ens with a width proportional to kBTq2. This is similar
to the temperature behavior for the incoherent structure
factor. The multi-phonon branches follow the modified
relations ω(q) = 2nω0| sin(qa/2n)| with n = 2, 3.
In Fig. 4(b) the comparison of the MD simulation (left
side) and the analytical formula Eq. (7a) (right side) is
shown. The temperature in this figure is kBT = 0.1,
which is higher than in Fig. 4(a). As a consequence
the one-phonon branch is broader. For both methods
of calculating S(q, ω) a high accuracy over 12 orders of
magnitude is possible. The accuracy is only limited by
the internal floating point precision.
The static structure factor from Eq. (9) is compared
to the results from the MD simulation in Fig. 5. S(q)
consists of a sequence of Lorentzian peaks at the recipro-
cal lattice points. The increasing width for larger wave
vectors shows again that no long-range order is present
in the one-dimensional model system.
The MD simulations and analytical formulas show a
perfect agreement. Therefore we conclude that MD sim-
ulations are a well-suited numerical tool for calculating
the structure factors of the one-dimensional model sys-
tems. Although we integrate the equations of motions
with a good precision only on a short time scale using the
simple Verlet-algorithm, the statistics extracted from the
trajectories are correct. This confirms our approach and
encourages us to proceed studying the phason dynamics
of Fibonacci chains in the next section.
IV. HARMONIC FIBONACCI CHAIN:
INFLUENCE OF THE QUASIPERIODICITY
By changing the interparticle equilibrium distances of
the HPC to those of a Fibonacci sequence with separa-
tions L and S we obtain the harmonic Fibonacci chain
(HFC). The interaction potential is left unchanged. Since
the incoherent structure factor is a function of the single
particle motion only, it does not depend on the equilib-
rium distances of the particles but only the interaction
potential. Hence, the incoherent structure factor of the
HPC and of the HFC are identical. For the coherent
structure factor the interparticle distances become im-
portant. Instead of Eq. (7a) we now have
S(q, ω) =
e−iωt
e−iqx
l exp
−q2σ2l (t)/2
Here x0l =
j=1 aj for l > 0, x
j=l aj for l < 0 and
x00 = 0 denote the equilibrium positions of the particles.
For kBT = 0 this gives the Fourier transform of the static
Fibonacci chain
S(q, ω) = δ(ω)
e−iqx
l (14)
which is well known32. It consists of a dense set of delta
peaks with varying intensity, positioned at the reciprocal
lattice points
q = q0 (h+ τh
′) , h, h′ ∈ N (15)
with q0 = 2πa ≈ 1.40.
As shown in Fig. 6 for the temperature kBT = 0.02,
the coherent structure factor S(q, ω) of the HFC consists
of many different branches all following the one-phonon
dispersion relation. Relative to each other the branches
are displaced in the q direction. They start at the recip-
rocal lattice points with the intensity of the respective
delta peak. Two-phonon branches are also found. The
broadening of the branches proportional to kBTq2 has
already been discussed for the HPC.
From these findings and Fig. 5 one can guess the
form of the static structure factor S(q) for the HFC:
Lorentzians are positioned at the reciprocal lattice points
with varying intensity. In Fig. 7 S(q) is shown for
q ∈ [0, 3π]. For small q a large number of Lorentzians
occur. With increasing q-value their width increases and
the stronger ones hide the weaker ones. Comparing the
position of the strong peaks and ignoring the change of
their widths and heights, there is a self-similarity in S(q).
The deflation factor is ∆q1/∆q2 = τ3 as indicated in
Fig. 7.
It is interesting to note that there are regions with very
few peaks. They are positioned around q = π, q = 2π,
etc. . The same regions are also special for the coherent
structure factor. As seen in Fig. 6 all the weak one-
phonon branches vanish towards q = π. Only the strong
one-phonon branch starting from the Bragg peak at q =
τq0 ≈ 2.27 remains.
V. DYNAMIC FIBONACCI CHAIN:
OCCURRENCE OF PHASON FLIPS
A. Phason flips
Let us now proceed to the anharmonic chains with pha-
son flips by investigating the dynamic Fibonacci chain
(DFC). It is built from identical particles that interact
with a symmetric double-well potential VDFC = x4−2x2.
First the notion of a phason flip has to be specified. To
do so we identify the position of the changes from L to
S and from S to L of the interparticle distances along
the particle trajectories. This is done in Fig. 8. Often
a L → S change and a S → L change lie next to each
other (particle distance 1) and the sequences LS and SL
are interchanged. But there are also many cases where
the positions of the two changes are separated by 0, 2,
3, or even more particle distances as marked by lines in
Fig. 8. Sometimes it is not possible to find a partner lo-
cally. Only in the long time average every L→ S change
will eventually cancel with a S → L change.
In the literature on the Fibonacci chain a phason flip
is understood as the exchange of a L and a neighboring
FIG. 4: (Color online) Coherent structure factor S(q, ω) of the HPC with N = 6500 particles from MD simulation. The
temperatures are for kBT = 0.01 (a) and for kBT = 0.1 (b). One-, two- and three-phonon branches are observed. They start
at the reciprocal lattice points 2πn/a. In (b) the output from the MD simulation (left side) is compared to the output of the
analytical formula Eq. (7a) (right side).
FIG. 5: Static structure factor S(q) of the HPC for kBT =
0.02. The symbols mark the data from a MD simulation with
N = 1000 particles, the line the result from the analytical
fomula, Eq. (9).
S particle distance. As we learned above there are also
other types of exchanges in the DFC. In the following we
denote by phason flip every pair of flip partners as those
connected by lines in Fig. 8. Note that the times and
positions of the phason flips are not well defined. Only
their number can be estimated by counting the changes
in the particle distances as we will do now.
The temperature dependence of the average phason flip
frequency ωflip is shown in Fig. 9. Phason flips start to
appear at about kBT = 0.1. At low kBT the average pha-
son flip frequency increases rapidly by thermal excitation
and ωflip � ω0. At higher temperatures kBT > 0.4 the
average phason flip frequency slowly saturates. In this
region the internal energy is comparable to the potential
hill.
B. Results for the structure factors
For the anharmonic chains no analytic results are avail-
able, in particular not at elevated temperatures when
phason flips occur. The incoherent structure factors of
the DFC and of the HFC/HPC differ remarkably and are
shown in Fig. 10 for q = π/a. The comparison leads to
the following conclusions:
1. At a fixed temperature, there are ω-ranges where
the curve for the DFC lies below the curve for
the HFC/HPC and vice versa. Since we have∫
Si(q, ω) dω = 1 from Eq. (3b), the integral area
between the two curves has to vanish.
2. At very low kBT and ω < 2ω0 ≈ 5.7 the curves
of the DFC and the HFC/HPC cannot be distin-
guished in logarithmic scale except for two small
bumps. They are a consequence of the anharmonic-
ity of the interaction potential of the DFC and not
related to the phason flips. At larger ω values the
multi-phonon edges have different heights.
3. Above kBT = 0.1 the one-phonon peak and the
multi-phonon edges in the curves for the DFC
broaden and weaken considerably faster than in the
curves for the HFC/HPC. They finally disappear at
kBT = 1.0.
4. No additional peaks or edges occur at any temper-
ature.
Further MD runs show that different q-values change
the temperature dependence, but the general features re-
main unchanged: The phonon peaks and edges broaden
and weaken much faster with increasing temperature for
the DFC than for the HFC/HPC.
Similar conclusions follow for the coherent structure
factor of the DFC. S(q, ω) for the DFC looks qualitatively
FIG. 6: (Color online) Coherent structure factor S(q, ω) of the HFC with N = 13000 particles for kBT = 0.02. It consists of a
dense set of phonon branches starting from the reciprocal lattice points.
FIG. 7: Static structure factor S(q) of the HFC from a MD simulation with 2000 particles at the temperature kBT = 0.02.
similar to S(q, ω) for the HFC except that the branches
broaden more quickly with temperature. To compare the
broadening let us look at plane cuts through S(q, ω) for
ω = 1.0 fixed. Three cuts for the temperature kBT =
0.05, 0.2, and 0.3 are shown in Fig. 11. A one-phonon
branch is located inside the cut resulting in a peak with
approximate Lorentzian line shape as indicated by the
fits in the figure. Only for lower temperatures the line
shape deviates from a Lorentzian, which is seen at the
base of the peak for kBT = 0.05.
The width of the Lorentzian in Fig. 11 is shown as
a function of temperature in Fig. 12 for both, the HFC
and the DFC. In the case of the HFC the width increases
linearly with temperature, as discussed in Sec. IV. There
is a deviation from linearity at low temperatures. This
is an artefact from the method of calculation from MD
simulation data. As explained in Sec. II B the structure
factor is convoluted with a Gaussian due to the finite
FIG. 8: Snapshot of the particle trajectories of the DFC at
the temperature kBT = 0.6. Changes in the particle distances
from L to S and S to L are marked with a cross (×) and a
plus (+). A L → S change and a S → L change combined
form a phason flip.
FIG. 9: Average flip frequency as a function of the temper-
ature kBT . In the temperature range of the figure: ωflip <
ω0 ≈ 2.83.
simulation time. The (narrow) Gaussian generates the
observed offset. The convolution with the Gaussian is
also responsible for the shape of the cuts in Fig. 11 at
low temperature, deviating from the Lorentzian shape.
There is one aspect of the DFC that has been ignored
up to now. Due to the phason flips the original per-
fect quasiperiodic long-range order is slowly decaying.
With progressing simulation time the chain becomes ran-
domized which, however, has no effect on the incoherent
structure factor. To test the influence of the randomiza-
tion on the coherent structure factor, a MD simulation
was started with a totally randomized Fibonacci chain.
The interaction potentials and the occurrence ratio of L
and S were not adapted. The result of the simulation is
shown in Fig. 13. Most of the details in S(q, ω) are lost by
the randomization process and the branches are strongly
FIG. 10: Incoherent structure factor Si(q, ω) of the DFC
(solid) and the HFC/HPC (dashed) for q = π/a at different
temperatures. The data was calculated using MD simulations
with N = 1000 particles.
FIG. 11: Cuts through S(q, ω) of the DFC for a fixed ω = 1.0
including a one-phonon peak. The solid curves are fits with
Lorentzians.
broadened, although the simulation has been carried out
at a low temperature of kBT = 0.02.
We summarize the results of this section: By the in-
troduction of the anharmonic double-well potential of
the DFC the phonon peaks, edges, and branches are
strongly broadened and weakened with increasing tem-
perature. There are two effects responsible for the broad-
ening: (1) The anharmonicity and single phason flips
affect the propagation of phonons. This is seen in the
incoherent structure factor. (2) The destruction of the
long-range order by a large number of phason flips. This
is the main effect of broadening for the coherent structure
factor. The phason flips appear more or less uniformly
distributed over the simulation time and the particles.
FIG. 12: Width of the Lorentzian peak of Fig. 11 as a function
of temperature for the HFC (inset) and the DFC.
FIG. 13: (Color online) Coherent structure factor S(q, ω) of
a randomized Fibonacci chain with 2000 particles at the tem-
perature kBT = 0.02. The same interaction potentials as in
Fig. 6 have been used.
VI. ASYMMETRIC FIBONACCI CHAIN:
COMPETING EIGENFREQUENCIES
A. Band gaps
In the next step we modify the double well potential of
the DFC. Like all the previous model systems the asym-
metric Fibonacci chain (AFC) is built of identical par-
ticles, but they interact with the more complicated po-
tential VAFC = VDFC + ∆V , see Fig. 14. The additional
term is given by
∆V (x) = χ(x2 − 1)2(�x+ x2/2− 1/2) (16)
with χ ∈ [0, 1] and � = ±1. This term has been cho-
sen in such a way that the positions of the minima at
x = ±1 are left invariant, but the curvatures of the poten-
tial at those points is changed to V ′′AFC(±1) = 8(1± �χ).
FIG. 14: Interaction potential VAFC(x) of the AFC for differ-
ent values of the parameter χ and � = 1. VDFC is shown with
dashed lines.
In harmonic approximation a particle will feel different
eigenfrequencies, depending on the nearest-neighbor con-
figuration SS, SL, LS, or LL. In the case of � = 1
it is ωLL > ωLS = ωSL > ωSS and for � = −1 it is
ωLL < ωLS = ωSL < ωSS. The sign change of � corre-
sponds to a mirror operation of VAFC about the y-axis.
The coherent structure factor of the AFC with χ = 0.1,
0.3, 0.5 and � = ±1 is shown in Fig. 15. Band gaps of
different widths appear and broaden with increasing val-
ues for χ. They are positioned at the frequencies where
one-phonon branches intersect each other. Their posi-
tions and widths are different for � = 1 and � = −1.
For � = 1 three large gaps and several smaller gaps ap-
pear, whereas for � = −1 only one very large gap, one
medium gap and several small gaps appear. The band
gaps are a consequence of the competing eigenfrequen-
cies due to the asymmetric potential in the same way as
band gaps appear in periodic systems with several parti-
cles and eigenfrequencies per unit cell.
B. Density of states
The gaps are seen more clearly in the density of states
(DOS) D(ω). The DOS per particle is calculated from
the velocity autocorrelation function by a Fourier trans-
D(ω) =
e−iωt
〈vl(t)vl(0)〉dt, ω ≥ 0. (17)
It is obtained from MD simulation data similar to the
incoherent structure factor: fl(q, t) has to be substited by
vl(t) in Eq. (5b). By interchanging the Fourier transform
and two time derivatives we can alternatively write
D(ω) =
e−iωt
〈xl(t)xl(0)〉dt. (18)
FIG. 15: (Color online) Coherent structure factor S(q, ω) of the AFC for different values of χ and � from MD simulations with
6500 particles at kBT = 0.01. Band gaps are observed. The band gaps also appear in the DOS D(ω).
After a taylor expansion of the exponentials in Eq. (3b)
a connection to the incoherent structure factor is found,
D(ω) = 2ω2 lim
Si(q, ω)
, ω 6= 0. (19)
For comparison we present the DOS of a harmonic
chain (HPC or HFC)
D(ω) =
4ω20 − ω2
for 0 < ω < 2ω0 (20)
and 0 elsewhere. The total number of states is normalized
D(ω) dω = kBT .
It can be checked, that the DOS of the harmonic chain
fits well to the DOS of the AFC for � = 0.0, except
small bumps that originate from the anharmonicity of the
potentials. In Fig. 15 the DOS of the AFC for different
values of χ and � are drawn. The band gaps appear at the
same frequencies and the same widths as in the coherent
structure factor.
VII. DISCUSSION AND CONCLUSION
At the end we would like to make some general re-
marks concerning phason flips and phason modes: In the
context of a hydrodynamic theory, phason flips can be as-
sociated with phason modes. It was noted quite early33,
that phason modes are diffusive, in contrary to the prop-
agating phonons. This means that phason flips are only
weakly coherent in space and in time, which of course we
have also observed here. As a result their influence on the
structure factor is small making it difficult though still
interesting to study them by scattering experiments. An-
other point concerns the connection of phason flips and
quasiperiodicity. There is no reason why phason flips
should only occur in quasicrystals. Since interaction po-
tentials are not sensitive on the long-range order, phason
flips in the form of atomic jumps can also occur in pe-
riodic complex intermetallic phases, which is supported
by recent experimental results34,35. In the case of our
model systems, it is equally possible to compare simula-
tions of a periodic LSLSLS . . . chain with harmonic and
double-well potentials respectively.
In conclusion, we have investigated the dynamics of
phonons and phason flips in one-dimensional model sys-
tems with molecular dynamics simulations. An efficient
algorithm made it possible to calculate the structure fac-
tors with high precision and in great detail. As a result
multi-phonon contributions, – although weak in compar-
ison to the one-phonon peaks and branches – have been
identified. By introducing phasons in the model systems
we were able to study their influence on the structure
factors, which is mainly a broadening of the character-
istic peaks, edges, and branches with temperature. The
broadening can be further split into a broadening due to
the disorder as a result of collective phason flips, i.e. a
static phason field, and a broadening due to the anhar-
monicity of the interaction potentials and single phason
flips. The work presented here is a first step. Further
studies on two-dimensional and three-dimensional model
systems with phason flips are under way.
Acknowledgments
We wish to thank T. Odagaki and M. Umezaki for
stimulating discussions.
∗ Electronic address: mengel@itap.uni-stuttgart.de
1 An introduction on quasicrystals can be found in: C.
Janot, Quasicrystals: A Primer (Oxford Science Publica-
tions, Oxford, 1992); M. de Boissieu, P, Guyot, and M.
Audier, in Lectures on Quasicrystals, edited by D. G. F.
Hippert (Les Éditions de Physique, Aussois, 1994), p. 1.
2 For a review on phonon excitations in quasicrystals see M.
Quilichini and T. Janssen, Rev. Mod. Phys. 69, 277 (1997)
and references therein.
3 P. Bak, Phys. Rev. Lett. 54, 1517 (1985); P. A. Kalugin,
A. Y. Kitayev, and L. S. Levitov, J. Physique Lett. 46,
L-601 (1985); D. Levine, T. C. Lubensky, S. Ostlund, S.
Ramaswamy, P. J. Steinhardt, J. Toner, Phys. Rev. Lett.
54, 1520 (1985).
4 M. de Boissieu, M. Boudard, B. Hennion, R. Bellissent, S.
Kycia, A. Goldman, C. Janot, and M. Audier, Phys. Rev.
Lett. 75, 89 (1995).
5 J. E. S. Socolar, T. C. Lubensky, and P. J. Steinhardt,
Phys. Rev. B 34, 3345 (1986).
6 R. Rosenfeld, M. Feuerbacher, B. Baufeld, M. Bartsch, M.
Wollgarten, G. Hanke, M. Beyss, U. Messerschmidt, K.
Urban, Phil. Mag. Lett. 72, 375 (1995).
7 P. A. Kalugin and A. Katz, Europhys. Lett. 21, 921 (1993).
8 R. Bluher, P. Scharwaechter, W. Frank, and H.
Kronmüller, Phys. Rev. Lett. 80, 1014 (1998).
9 W. Steurer, Acta. Cryst. A 61, 28 (2005).
10 K. Edagawa, K. Suzuki, S. Takeuchi, Phys. Rev. Lett 85,
1674 (2000).
11 E. Abe, S. J. Pennycook, and A. P. Tsai, Nature 421, 347
(2003).
12 M. Engel and H.-R. Trebin, Phil. Mag 85, 2227 (2005).
13 G. Coddens, S. Lyonnard, B. Hennion, Y. Calvayrac, Phys.
Rev. B 62, 6268 (2000).
14 The coherent structure factor S(q, ω) is also called ‘(total)
scattering function’, ‘response function’, or ‘dynamic(al)
structure factor’. The incoherent structure factor Si(q, ω)
is alternatively known as ‘incoherent scattering function’
and ‘self part of the structure factor’.
15 M. de Boissieu, K. Shibata, R. Currat, A. Q. R. Baron, S.
Tsutsui, A.-P. Tsai, J. Non-Cryst. Solids 334&335, 303
mailto:mengel@itap.uni-stuttgart.de
(2004) and references therein.
16 M. Krisch, R. A. Brand, M. Chernikov, H. R. Ott, Phys.
Rev. B 65, 134201 (2002).
17 M. de Boissieu, M. Boudard, R. Bellissent, M. Quilichini,
B. Hennion, R. Currat, A. I. Goldman, and C. Janot, J.
Phys.: Condens. Matter 5, 4945 (1993).
18 M. Boudard, M. de Boissieu, S. Kycia, A. I. Goldman, B.
Hennion, R. Bellissent, M. Quilichini, R. Currat, and C.
Janot, J. Phys.: Condens. Matter 7, 7299 (1995).
19 S. Lyonnard, G. Coddens, Y. Calvayrac, and D. Gratias,
Phys. Rev. B 53, 3150 (1996).
20 G. Coddens, S. Lyonnard, and Y. Calvayrac, Phys. Rev.
Lett. 78, 4209 (1997).
21 G. Coddens, S. Lyonnard, B. Hennion, and Y. Calvayrac,
Phys. Rev. Lett. 83, 3226 (1999).
22 J. A. Ashraff and R. B. Stinchcombe, Phys. Rev. B 39,
2670 (1989).
23 G. G. Naumis, C. Wang, M. F. Thorpe, and R. A. Barrio,
Phys. Rev. B 59, 14302 (1999).
24 J. Hafner, M. Krajĉ́ı, J. Phys.: Condens. Matter 5, 2489
(1993).
25 T. Yoshida, K. Shobu, and H. Mori, Prog. Theor. Phys.
66, 759 (1981).
26 V. J. Emery and J. D. Axe, Phys. Rev. Lett. 40, 1507
(1978).
27 J. Florencio and M. H. Lee, Phys. Rev. A 31, 3231 (1985).
28 G. Radons, J. Keller, and T. Geisel, Z. Phys. B 61, 339
(1985).
29 The correlation functions have the following interpretation:
G(x, t)∆x is the probability to find a particle in [x, x+∆x]
at time t, if any particle has been at the origin at t = 0.
Ga(x, t)∆x is the probability to find a particle in [x, x+∆x]
at time t, if the same particle has been at the origin at
t = 0.
30 M. Bée, Quasielastic neutron scattering, IOP Publishing,
Bristol (1988).
31 The notation ‘x-phonon’ is adopted freely from the lan-
guage of neutron scattering. However in the case of a one-
dimensional system there is no regular phonon-expansion
possible, as it is done for three-dimensional crystals.
32 D. Levine and P. J. Steinhardt, Phys. Rev. Lett. 53, 2477
(1984).
33 T. C. Lubensky, S. Ramaswamy, and J. Toner, Phys. Rev.
B 32, 7444 (1985).
34 U. Dahlborg, W. S. Howells, M. Calvo-Dahlborg, and J.
M. Dubois, J. Phys.: Condens. Matter 12, 4021 (2000).
35 J. Dolinšek, T. Apih, P. Jeglič, M. Klanǰsek, J. Non-Cryst.
Solids 334&335, 280 (2004).
	Introduction
	Motivation
	Model systems
	Structure factors from molecular dynamics
	Definition of the structure factors
	Molecular dynamics simulations
	Harmonic periodic chain as reference system
	Analytic calculations
	Simulation results
	Harmonic Fibonacci chain: Influence of the quasiperiodicity
	Dynamic Fibonacci chain: Occurrence of phason flips
	Phason flips
	Results for the structure factors
	Asymmetric Fibonacci chain: Competing eigenfrequencies
	Band gaps
	Density of states
	Discussion and conclusion
	Acknowledgments
	References
