Periodic layers of a dodecagonal quasicrystal and a floating hexagonal crystal in

sedimentation-diffusion equilibria of colloids

Harini Pattabhiraman, and Marjolein Dijkstra

Citation: The Journal of Chemical Physics 147, 104902 (2017);

View online: https://doi.org/10.1063/1.4993521

View Table of Contents: http://aip.scitation.org/toc/jcp/147/10

Published by the American Institute of Physics

Articles you may be interested in

Phase behaviour of quasicrystal forming systems of core-corona particles

The Journal of Chemical Physics 146, 114901 (2017); 10.1063/1.4977934

Phase and vacancy behaviour of hard “slanted” cubes

The Journal of Chemical Physics 147, 124501 (2017); 10.1063/1.5001483

The effect of finite pore length on ion structure and charging

The Journal of Chemical Physics 147, 104708 (2017); 10.1063/1.4986346

Diffusion and interactions of point defects in hard-sphere crystals

The Journal of Chemical Physics 146, 244905 (2017); 10.1063/1.4990416

Liquid bridging of cylindrical colloids in near-critical solvents

The Journal of Chemical Physics 147, 104701 (2017); 10.1063/1.4986149

Variation of ionic conductivity in a plastic-crystalline mixture

The Journal of Chemical Physics 147, 104502 (2017); 10.1063/1.5001946

THE JOURNAL OF CHEMICAL PHYSICS 147, 104902 (2017)

Periodic layers of a dodecagonal quasicrystal and a ﬂoating hexagonal

crystal in sedimentation-diﬀusion equilibria of colloids

Harini Pattabhiraman

a)

and Marjolein Dijkstra

b)

Department of Physics, Soft Condensed Matter, Debye Institute for Nanomaterials Science, Utrecht University,

Princetonplein 5, 3584 CC Utrecht, The Netherlands

(Received 29 June 2017; accepted 30 August 2017; published online 14 September 2017)

We investigate the behaviour of a system of colloidal particles interacting with a hard-core and a

repulsive square shoulder potential under the inﬂuence of a gravitational ﬁeld using event-driven

Brownian dynamics simulations. We use a ﬁxed square shoulder diameter equal to 1.4 times the hard-

core diameter of the colloids, for which we have previously calculated the equilibrium phase diagram

considering two-dimensional disks [H. Pattabhiraman et al., J. Chem. Phys. 143, 164905 (2015) and

H. Pattabhiraman and M. Dijkstra, J. Phys.: Condens. Matter 20, 094003 (2017)]. The parameters

in the simulations are chosen such that the pressure at the bottom of the sediment facilitates the

formation of phases in accordance with the calculated phase diagram of the two-dimensional system.

It is surprising that we observe the formation of layers with dodecagonal, square, and hexagonal

symmetries at the relevant pressures in the three-dimensional sedimentation column. In addition,

we also observe a re-entrant behaviour exhibited by the colloidal ﬂuid phase, engulﬁng a hexagonal

crystal phase, in the sedimentation column. In other words, a ﬂoating crystal is formed between the

colloidal ﬂuid regions. Published by AIP Publishing. [http://dx.doi.org/10.1063/1.4993521]

I. INTRODUCTION

In the case of colloidal suspensions consisting of particles

with sizes on the order of micro-meters, the effect of the grav-

itational force is not negligible. Under these conditions, the

gravitational energy and the thermal energy of the colloids are

comparable. This leads to the formation of a spatially inho-

mogeneous suspension in which the density of the particles

varies along the height of the suspension. The inhomogeneous

distribution of the colloidal particles along the height under

the inﬂuence of gravity is termed as sedimentation.

As a result of this inhomogeneous density distribution in

colloidal suspensions, the particles at the bottom of the sedi-

ment can crystallise when they reach a certain density. Exper-

imentally, sedimentation is regarded as a prevalent tool to

extract information regarding the equilibrium phase behaviour

of the system. For example, the measured concentration pro-

ﬁles can be inverted to obtain the osmotic equation of state.

Sedimentation processes can also be used the other way

around, i.e., they can be used to validate the theoretically calcu-

lated equilibrium phase behaviour of a system and thereby its

bulk phase behaviour. For example, a system of hard spheres

has been a model system for sedimentation studies. This sys-

tem exhibits a phase behaviour characterised by a ﬂuid phase at

low densities and a face-centered cubic phase at high densities.

This behaviour, which was earlier theoretically predicted,

1

has

later been corroborated by sedimentation studies.

2

Sedimentation behaviour of various charged particles,

3

mixtures of hard particles,

4

and particles of different shapes

5

which result in the formation of various colloidal crystals has

a)

h.pattabhiraman@uu.nl

b)

m.dijkstra@uu.nl

been extensively studied. However, similar studies involving

the formation of quasicrystals (QCs) by sedimentation has

received less attention. Quasicrystals are solids with long-

range orientational order and no periodicity and may exhibit

intriguing properties including the formation of photonic band

gaps.

6–11

Hence, the experimental realisation of quasicrystals

is a topic of intense research.

12–18

In our previous studies,

19,20

we have theoretically cal-

culated the phase diagram of a two-dimensional system of

particles interacting with a hard-core and repulsive square

shoulder (HCSS) potential. At a shoulder width of 1.4 times

the hard core diameter of the particle, we ﬁnd a stable random-

tiling dodecagonal quasicrystal. In this work, we investigate

how the phase behaviour in a two-dimensional HCSS sys-

tem extends to three-dimensional systems and whether or not

multiple layers of this quasicrystal can be self-assembled by

sedimentation in three dimensions.

This paper is organised as follows. In Sec. II, we present

the simulation and analysis methods used. In Sec. III, we indi-

vidually discuss the formation of layers with various symme-

tries, and ﬁnally we present the conclusions and the direction

of future studies in Sec. IV.

II. METHODS

We ﬁrst explain the simulation model and computational

methods used for this study in Sec. II A and then give an

account of the analysis methods in Sec. II B.

A. Computational methodology

We perform Event-Driven Brownian Dynamics (EDBD)

simulations of N spherical particles of diameter σ

HS

and buoy-

ant mass m interacting with the HCSS potential in the NVT

0021-9606/2017/147(10)/104902/11/$30.00 147, 104902-1 Published by AIP Publishing.

104902-2 H. Pattabhiraman and M. Dijkstra J. Chem. Phys. 147, 104902 (2017)

ensemble. The HCSS potential can be written as a sum of a

hard-sphere potential V

HS

(r) and a square-shoulder potential

V

SS

(r), i.e.,

V

HCSS

(r) = V

HS

(r) + V

SS

(r), (1)

where

V

HS

(r) =

∞, r ≤ σ

HS

0, r > σ

HS

(2)

and

V

SS

(r) =

, r ≤ δ

0, r > δ

, (3)

where r is the interparticle center-of-mass distance and > 0

is the height of the square shoulder.

In the EDBD method, a sequence of collision events

involving only two particles at any given instant is com-

puted. During the simulation, the velocities of the particles

are randomly adjusted at regular intervals ∆t as

v(t + ∆t) = α

t

v(t) + β

t

v

R

(t), (4)

where v(t) and v(t + ∆t) are, respectively, the velocities of the

particles before and after the stochastic velocity adjustment,

v

R

(t) is a 3-D Gaussian variable with a mean of 0 and vari-

ance of k

B

T/m, with k

B

as the Boltzmann constant and T as

the temperature. Further, α

t

has a value 1/

√

2 with a proba-

bility ν∆t and 1 otherwise. The temperature is kept constant

by setting β

t

=

q

1 − α

2

t

. In accordance to previous EDBD

simulations,

21,22

we set ν to 10τ

−1

MD

and ∆t to 0.01τ

MD

, where

τ

MD

is the unit of time of an event-driven molecular dynamics

simulation given as τ

MD

=

√

m/k

B

T σ

HS

.

The present simulation method is an event driven molecu-

lar dynamics (EDMD) method which is extended with Brow-

nian motion. Previous experiments and simulations show that

the structure of the crystalline sediment depends strongly on

the Brownian time, the sedimentation time, and the initial vol-

ume fraction.

21,23

The faster the sedimentation, the higher the

initial volume fraction, and the less time the system has avail-

able for equilibration during sedimentation and the crystalline

sediment will become more defective. Brownian motion is,

thus, important for equilibration during sedimentation. In this

case, the Brownian time is set by α

t

and ν. It is convenient to

deﬁne a dimensionless Peclet number g

*

that is equal to the

inverse gravitational length and a dimensionless particle ﬂux

as deﬁned by the ratio of the Brownian time and sedimentation

time. In our simulations, the Peclet number ranges from 1 to

5, whereas the particle ﬂux is in the range of 0.01-0.05, where

previous simulations and experiments show that there is ample

time for equilibration and least defective crystalline sediments

are obtained.

21,23

The simulation box of volume V has a square cross section

of area A and is elongated in the z-direction. Periodic boundary

conditions are applied along the cross section. In the elongated

direction, the particles are conﬁned between two smooth hard

walls at z = 0 and z = H with H as the height of the sedimen-

tation column. The height H is chosen such that the density

at z =

(

H − σ

HS

/2

)

is sufﬁciently small, which allows us to

consider the system to be inﬁnitely long in the z-direction. We

perform simulations starting with a non-overlapping isotropic

ﬂuid state ﬁlling the entire sedimentation column with packing

fraction η = 0.01. To mimic sedimentation experiments, these

particles are further subjected to a gravitational ﬁeld, which is

expressed as an external potential φ(z) written as

φ(z) =

mgz, σ

HS

/2 ≤ z ≤ (H − σ

HS

/2)

∞, otherwise

, (5)

where g is the acceleration due to gravity and z is the vertical

coordinate of the particle. The effect of the gravitational ﬁeld

on the particles is quantiﬁed in terms of the gravitational Peclet

number deﬁned as g

∗

= mgσ

HS

/k

B

T.

In this work, we scrutinise the kinetic formation of

the thermodynamic stable phases described for the two-

dimensional HCSS system with δ = 1.40σ

HD

, where σ

HD

is the hard-disk diameter, as used in our previous works.

19,24

To do so, we perform simulations such that the pressure mea-

sured at the bottom of the sedimentation column, i.e., at z = 0,

corresponds to the region of stability of a particular phase in the

2-D phase diagram. This pressure is calculated as P

∗

= βP(z

= 0)σ

3

HS

= g

∗

.ρ

∗

A

, where ρ

∗

A

is the mean areal density deﬁned

as the number of particles at the bottom of the sample per

area ρ

∗

A

= N σ

2

HS

/A, and we compare P

*

with P

∗

2D

= βPσ

2

HD

directly. However, we note that the direct comparison of pres-

sures between two- and three-dimensional systems may not

always hold true. For example, the ﬂuid (hexatic)-solid transi-

tion in the case of hard disks is at P

∗

2D

= 9.17,

25

whereas that

for hard spheres is at βPσ

3

HS

= 11.57.

26

It is good to remind

the readers that we use a three-dimensional system of spheres

here to validate the phase behaviour of a two-dimensional sys-

tem of disks. The underlying reason for this will be explained

in Secs. II B and III A–III C.

The phases considered in this study are dodecagonal qua-

sicrystal, square, low-density hexagonal, and ﬂuid phases. We

especially focus on the possibility of the formation of the qua-

sicrystal and, thus, consider the cross section to be squares with

sides of length 58σ

HS

, which can accommodate a random-

tiling dodecagonal quasicrystal at a density ρ

∗

= 1.07, similar

to that used in Ref. 19. The list of parameters used to simu-

late the different phases is given in Table I. The corresponding

values of pressures at the bottom P

∗

2D

= βP(z = 0)σ

3

HS

are

marked in the phase diagram given in Fig. 1. As a supple-

mentary study, we use two parameter sets having different

Peclet numbers to simulate the quasicrystal. As mentioned

above, a higher value of the Peclet number results in a condi-

tion of high settling rate of the particles, i.e., the particles do

not have enough time to rearrange and equilibrate, and vice

versa. Using systems with different Peclet numbers allows us

TABLE I. System parameters used in the EDBD simulations of a HCSS

system with δ = 1.40σ

HS

under gravity.

N k

B

T/ g

*

βP(z = 0)σ

3

HS

Stable phase at bottom

5 × 10

4

0.25 2.00 30.0 Quasicrystal

2 × 10

4

0.25 5.00 30.0 Quasicrystal

2 × 10

4

0.25 4.00 23.8 Square

5 × 10

4

0.15 0.67 10.0 Low-density hexagonal

1 × 10

5

0.15 0.67 20.0 Fluid

104902-3 H. Pattabhiraman and M. Dijkstra J. Chem. Phys. 147, 104902 (2017)

FIG. 1. Phase diagram in the (reduced) pressure-temperature plane for a two-

dimensional HCSS system with shoulder width δ = 1.40σ

HD

. The reduced

quantities are deﬁned as P

∗

2D

= βPσ

2

HD

and T

∗

= k

B

T/. The phases marked

are ﬂuid (FL), low-density hexagonal (LDH), square (SQ), dodecagonal qua-

sicrystal (QC), and high-density hexagonal (HDH). The crosses denote the

state points corresponding to the pressures at the bottom of the sediment.

to study the effect of the settling rate on the formation of the

quasicrystal.

B. Structural analysis

In order to characterise the different phases, we employ

an analysis method that is two-dimensional in nature since the

phases observed in the sedimentation column have a layered

structure and resemble the phases observed for the HCSS sys-

tem in 2-D. Speciﬁcally, we ﬁrst identify different layers of the

sediment and then carry out the following analysis procedure

in these layers. We construct the polygonal tiling of the layer

and calculate the two-dimensional m-fold bond orientational

order parameter (BOO) of each particle j in layer l, χ

l

m

(j), and

the average BOO of each layer χ

l

m

as explained in Ref. 20.

The polygonal tiling of each layer is constructed by draw-

ing bonds between the neighbouring particles of each particle j,

which are at a center-of-mass distance smaller than the square

shoulder diameter δ from particle j.

We, then, calculate the m-fold BOO of each particle j in

layer l as

χ

l

m

(j) =

1

N

B

(j)

N

B

(j)

X

k=1

exp(imθ

r

jk

)

2

, (6)

where m is the symmetry of interest, r

jk

is the center-of-mass

distance vector between two neighbours j and k, θ

r

jk

is the angle

TABLE II. Method of classiﬁcation of particle j belonging to layer l

according to its bond orientational order (BOO) χ

l

m

(j).

Symmetry BOO conditions Colour

Fluid/other (OT) χ

l

4

(j), χ

l

6

(j), χ

l

12

(j) < 0.5 Orange

Crystal χ

l

4

(j), χ

l

6

(j), χ

l

12

(j) > 0.5

Square (SQ) χ

l

4

(j) > χ

l

6

(j), χ

l

12

(j) Purple

Hexagonal (HX) χ

l

6

(j) > χ

l

4

(j), χ

l

12

(j) Green

Dodecagonal (QC) χ

l

12

(j) > χ

l

4

(j), χ

l

6

(j) Red

FIG. 2. Colour scheme for classes of particles based on the BOO classiﬁcation

described in Table II.

between r

jk

and an arbitrary axis, and N

B

(j) is the number of

neighbours of particle j in the same layer. For each particle j, we

calculate χ

l

4

(j), χ

l

6

(j), and χ

l

12

(j), respectively, representing

square, hexagonal, and dodecagonal symmetries.

The particles are classiﬁed based on their BOO according

to the method given in Table II. We consider a particle to be

ﬂuid-like if each of the three χ

l

m

(j) is less than 0.5. On the other

hand, if each of χ

l

m

(j) is greater than 0.5, then a particle is said

to have symmetry m1 if χ

l

m1

(j) is greater than the other two,

namely, χ

l

m2

(j) and χ

l

m3

(j). Further, we identify and colour the

particles according to the following scheme: particles of square

symmetry in purple, those of hexagonal in green, dodecagonal

in red, and ﬂuid-like in orange as shown in Fig. 2.

After calculating the BOO of each particle, the average

BOO of each layer is then evaluated as

27

χ

l

m

=

1

N

l

N

l

X

j=1

χ

l

m

(j), (7)

where N

l

is the number of particles in each layer.

III. RESULTS AND DISCUSSION

In this section, we consider individually the different sed-

imentation simulations carried out to obtain the various stable

phases calculated for the two-dimensional HCSS system.

A. Formation of layers with dodecagonal symmetry

We start with the formation of layers with dodecagonal

symmetry. In this section, we present the sedimentation sim-

ulations using two different Peclet numbers in order to assess

the effect of the settling rate on the formation of the quasicrys-

tal. We ﬁrst compare the formation of the quasicrystal formed

in these simulations and then analyse the driving force behind

the formation of these layers. Finally, we review the valid-

ity of the phase diagram given in Fig. 1 by comparing the

phases formed in the sedimentation column. The dodecago-

nal quasicrystal (QC) which, as seen in the phase diagram, is

sandwiched between two periodic crystal phases with square

and hexagonal symmetries. Thus, it is interesting to note if and

how the interfaces between the quasicrystal and the periodic

crystals are formed in the sedimentation column.

We ﬁrst present a typical conﬁguration of the sediment

forming quasicrystalline layers in Fig. 3. The panels on the

left correspond to simulations with a Peclet number g

*

= 5.0

and those on the right are obtained for g

*

= 2.0. The particles

here are coloured according to the convention explained in

Fig. 2. We notice the formation of about two quasicrystalline

layers for g

*

= 5.0 and about four quasicrystalline layers for

g

*

= 2.0. This difference in the number of layers is due to a

104902-4 H. Pattabhiraman and M. Dijkstra J. Chem. Phys. 147, 104902 (2017)

FIG. 3. Comparison of the quasicrys-

tal (QC) sediment formed for Peclet

numbers g

*

= 5.0 (left) and 2.0 (right).

Side view of a conﬁguration of the sed-

imentation column obtained at t/τ

MD

= 800 (top). The particles are coloured

according to their individual BOO: qua-

sicrystal (red), square (purple), hexag-

onal (green), and ﬂuid (orange). The

m-fold BOO of each layer with time cal-

culated for symmetries m = 4 (middle)

and 12 (bottom) showing, respectively,

the formation of layers with square and

dodecagonal symmetries.

difference in the height range that corresponds to the pres-

sure range of the stable quasicrystal phase. This height range

decreases with increasing g

*

. Additionally, we observe that

most of the particles seen in these layers are coloured either

in purple or red which, respectively, represent square or

dodecagonal symmetries. Therefore, we follow the dynam-

ics of the formation of these layers by calculating the BOO

χ

l

4

and χ

l

12

of each layer as a function of time. The time

evolution of χ

l

4

is given in the middle panel in Fig. 3 and

that of χ

l

12

is given at the bottom. In these time evolution

heat maps, the time scale t/τ

MD

is plotted on the horizon-

tal axis and the layer number is plotted along the vertical

axis.

We make the following observations from these plots: (1)

In both cases, the value of χ

l

12

is higher than that of χ

l

4

, which

conﬁrms the dodecagonal symmetry of these layers. (2) With

increasing time, we ﬁnd that the fraction of ﬂuid in the sed-

imentation column decreases, as seen by the receding blue

region in these plots. Alternatively, this means that more crys-

talline layers are formed with time. (3) The value of χ

l

12

at a

given time decreases as we go up in the sediment indicating that

the layers on the top are more ﬂuid-like than the bottom layers.

(4) Finally, we see that χ

l

12

obtained for the sediment at higher

Peclet numbers is larger than that at lower Peclet numbers. This

is counter-intuitive as this suggests that faster settling of the

particles result in the formation of a less-defective quasicrystal.

We investigate this further by plotting the polygonal tiling

constructed for the bottom two layers as a function of time for

both the sediments. The top view of these tilings is given in

Fig. 4 for g

*

= 5.0 and in Fig. 5 for g

*

= 2.0. Two striking fea-

tures are conspicuous from these polygonal tilings. First, there

are large portions of connected square tilings in the sediment

obtained for the lower Peclet number, while the square tilings

are more uniformly distributed in the case of the high Peclet

number sediment. Second, the position of the tiles in the ﬁrst

and second layers seems to be on top of each other. Let us now

evaluate each of these observations separately.

Firstly, we analyze the tilings and quantify the square tiles

by calculating the ratio of areas occupied by the square tiles

to those of the triangle tiles. This also relates to the square-

triangle tiling description of a dodecagonal quasicrystal, where

the maximum entropy of the tiling corresponds to equal areas

of squares and triangles.

28–30

In the current sediments, we ﬁnd

that the ratio of the areas of squares to triangles for g

*

= 2.0

is 1.40 ± 0.05 and for g

*

= 5.0 is 1.15 ± 0.03. In other words,

in both cases, there are more squares formed than in an ideal

dodecagonal tiling. This excess of squares is larger for the low

Peclet number sediment. This can be explained as follows.

A lower Peclet number refers to a lower rate of sedimentation

and, thus, a longer relaxation time for the particles to rearrange.

The density at the bottom of the sample increases slowly since

the beginning of the sedimentation simulation, as more and