Anisotropic spatially heterogeneous dynamics on the $α$ and $β$ relaxation time scales studied via a four-point correlation function

4.1 MB·English

by  Elijah Flenner

Anisotropic spatially heterogeneous dynamics on the α and β relaxation time scales
studied via a four-point correlation function
Elijah Flenner and Grzegorz Szamel
Department of Chemistry, Colorado State University, Fort Collins, CO 80523
(Dated: January 5, 2009)
We examine the anisotropy of a four-point correlation function G4(~k,~r;t) and it’s associated
structure factor S4(~k,~q;t) calculated using Brownian Dynamics computer simulations of a model
glassforming system. Thesecorrelation functionsmeasure thespatialcorrelations oftherelaxation
of different particles, and we examine the time and temperature dependenceof theanisotropy. We
find that the anisotropy is strongest at nearest neighbor distances at time scales corresponding to
9 the peak of the non-Gaussian parameter α2(t)=3hδr4(t)i/[5hδr2(t)i2]−1, but is still pronounced
0
0 around the α relaxation time. We find that the structure factor S4(~k,~q;t) is anisotropic even for
thesmallestwavevectoraccessibleinoursimulationsuggestingthatoursystem(andothersystems
2
commonlyusedincomputersimulations)maybetoosmalltoextractthe~q→0limitofthestructure
n factor. Wefindthatthedetermination ofadynamiccorrelation lengthfromS4(~k,~q;t)isinfluenced
a
by the anisotropy. We extract an effective anisotropic dynamic correlation length from the small q
J
behaviorof S4(~k,~q;t).
5

PACSnumbers:
]
t
f
o
I. INTRODUCTION corresponding the the first peak of the pair correlation
s
. function are most likely to move in the same direction,
t
a particles that start at a separation corresponding to the
m Itisnowgenerallyacceptedthatuponapproachingthe first minimum are more likely to initially move in oppo-
glass transition, the liquid’s dynamics are becoming in-
- site directions.
d creasingly heterogeneous [1, 2, 3]. However, the details
n of the spatial and temporal characteristics of dynamic
o heterogeneities are still being debated. In particular,
c the connection between heterogeneous dynamics and a In view of the experimental and simulationalevidence
[
growing dynamic correlation length has been the topic for anisotropic correlations of particle’s displacements,

1 of many simulations [4, 5, 6, 7, 8, 9, 10, 11] and a few it should not be a surprise that four-point correlation
v experimental studies [12, 13, 14]. Four-point correlation functions designed to study these dynamics can also be
3 functionshavebeenintroducedtofacilitatethequantita- anisotropic. However, this anisotropy is normally stud-
9 tive description of heterogeneous dynamics. The analy- ied for times less than the α relaxation time, thus it is
4
0 sis ofthe spatialdecayofthese correlationfunctions was uncertain if understanding this anisotropy is important
. used to extract a dynamic correlation length. Recently, forthestructuralrelaxationoftheliquid. Previously[19]
1 the mode-coupling theory has been extended and a the- we reported on a four-point correlation function that is
0
oretical treatment of four-point correlation functions is anisotropic on the α relaxation time scale as well as the
9
0 starting to emerge [6, 12, 15, 16, 17, 18]. However, in β relaxation time scale for a model glass forming liq-
: most simulation studies and in some theoretical treat- uid. Since the spatial decay of this correlation function
v ments these four-point correlation functions have been can be used to determine a dynamic length scale, the
i
X assumed to be isotropic or they are isotropic by design. anisotropyintroducesacomplicationindetermining this
length scale.
r
a Researchers have noticed anisotropy in the correlated
motionofparticlesontheβ relaxationtimescale,andre-
cently this anisotropic motion has also been reported on
the α relaxation time scale [19]. Doliwa and Heuer [20] In this paper we expand on previous work [19]. Af-
reported anisotropic correlated motion in a hard sphere ter describing the simulation in Sec. II, we explore the
system on the β relaxation time scale. Anisotropic mo- anisotropic correlated dynamics by examining a four-
tion has also been extensively studied by Donati et al. point correlation function G4(~k,~r;t), Sec. III, and the
andGebremichaelet al.[8,21]whodescribedthemotion associatedstructurefactorS (~k,~q;t),Sec.IV. Weexam-
4
of ”mobile” particles as ”string-like”, with mobile par- ine the anisotropy at around nearest neighbor distances,
ticles following each other in one dimensional ”strings”. which corresponds to local rearrangement of particles
Weeks et al. [22] have reportedanisotropicdynamics as- and its cage, and at large distances. We examine how
sociated with the break down of the ”cage” surrounding the anisotropyinfluences the determination of a growing
a particle. They found that the correlations of the par- length scale accompanying the glass transition, and de-
ticle’s displacements depends on the initial separationof termine an effective anisotropic correlation length. We
the particles. While particles that start at a separation finish with a discussion of the results in Sec. V.
2
II. SIMULATION relaxation of different particles. Consider the function
Fˆ (~k;t)=e−i~k·[~rn(t)−~rn(0)], (4)
We performed Brownian dynamics simulations of an n
80:20 binary mixture of 1000 particles introduced by where~r (t) is the position of particle n at a time t. The
n
Kob and Andersen [23, 24]. The interaction potential ensembleaverageofFˆ (~k;t)istheself-intermediatescat-
is Vαβ(r) = 4ǫαβ[(σαβ/r)12 − (σαβ/r)6] where α, β ∈ tering function F (k;nt), thus we will term Fˆ (~k;t) the
{A,B}, ǫ = 1.0, ǫ = 1.5, ǫ = 0.5, σ = 1.0, s n
AA AB BB AA microscopic self-intermediate scattering function. The
σ = 0.8, and σ = 0.88 and the interaction poten-
AB BB four-point correlationfunction
tialis cutat2.5σ . Periodicboundaryconditionswere
αβ
used with a box length of 9.4 σ . The equation of mo- V
AA G (~k,~r;t)= hFˆ (~k;t)Fˆ (−~k;t)δ[~r−~r (0)]i
tion for the position of particle i is 4 N2 n m nm
X
n6=m
1 (5)
~r˙ = F~ (t)+~η (t), (1)
i ξ i i measures the correlations between the microscopic self-
0
intermediate scattering function at time t, pertaining to
where ξ = 1.0 is the friction coefficient of an isolated
0 particles that are separated by a vector ~r at the initial
particle and the force acting on a particle i is
time. In Eq. (5) ~r = ~r −~r , V is the volume, and
nm n m
F~ =−∇ V (|~r −~r |) (2) N is the number of particles. Notice that G4(~k,~r;0) =
i i αβ i n
g(r) where g(r) is the pair correlation function. In this
X
n6=i
work we choose |~k| to have the same value as the one
with ∇i being the gradient operator with respect to ~ri. that determines the α relaxationtime, i.e. |~k| is located
Therandomforce~η(t)satisfiesthefluctuationdissipation
aroundthefirstpeakofthepartialstaticstructurefactor
relation for the A particles, |~k|=7.25.
h~ηi(t)~ηj(t′)i=2D0δij1, (3) It should be noted that G4(~k,~r;t) is, in general, com-
plex. Its real and imaginary parts can be written in the
where D0 =kBT/ξ0, kB is Boltzmann’s constant, and 1 following form
is the unit tensor. The results are presented in terms of
reduced units with σAA, ǫAA/kB, and σA2Aξ0/ǫAA being Re[G4(~k,~r;t)]= (6)
the units of length, energy, and time, respectively. Since V
cos{~k·[~r (t)−~r (0)]}δ[~r−~r (0)]
the equation of motion allows for diffusion of the center N2 X D nm nm nm E
ofmass,allresults arepresentedrelativeto the centerof n6=m
mass.
We present results for temperatures T = 0.45, 0.47, Im[G (~k,~r;t)]= (7)
4
0.5, 0.55, 0.6, 0.8, 0.9, and 1.0. The onset of supercool- V
ing is around T = 1.0 and we use Tc = 0.435 as the −N2 X Dsin{~k·[~rnm(t)−~rnm(0)]}δ[~r−~rnm(0)]E
mode coupling temperature. As a means to expand the n6=m
temperature scale, we will plot some quantities versus
Eqs. (7-8) show that particles which are getting closer
ǫ=(T −T )/T . The equation of motion was integrated
c c together or farther apart along the direction of vector
usingaHeunalgorithmwithasmalltimestepof5×10−5. ~k (i.e. are moving in the opposite direction or in the
We ran an equilibration run that was at least half as
same direction along ~k) make the same contribution to
long as a production run, and four production runs at
each temperature. The results are an average over the the real part of G4(~k,~r;t) but opposite contributions to
production runs. We present results only for the larger itsimaginarypart. Inparticular,particlesmovingfarther
and more abundant A particles. We define the α relax- apart along the direction of vector ~k make a negative
ation time τα as through relation Fs(~k;τα) = e−1 for a contribution to the imaginary part of G4(~k,~r;t).
wave vector around the first peak of the partial static In severalother simulationaland experimentalstudies
structurefactorfortheAparticles,whichcorrespondsto [9, 10, 25, 26] four-point correlation functions involving
|~k|=7.25. single-particle overlaps rather than the microscopic self-
intermediate scattering functions were investigated. For
example, Lacevic et. al [9] used the following function
III. FOUR-POINT CORRELATION FUNCTION [27]
G4(~k,~r;t)
V
gol(r;t)= hw (a;t)w (a;t)δ[~r−~r (0)]i, (8)
4 N2 n m nm
A. Definition and connection with overlap nX6=m
correlations
where w (a;t) is the overlap function pertaining to par-
n
ticle n,
We study a four-point correlation function that mea-
sures the spatial and temporal correlations between the w (a;t)=θ(a−|~r (t)−~r (0)|). (9)
n n n
3
We wouldliketo pointoutthat g (r;t) canbe expressed
4
in terms of functions which are generalizations of our
G (~k,~r;t),
4
RRee[[GG((kk,,rr;;tt))]]
44
d~k d~k 2.5
gol(r;t)= 1 2f(k ;a)f(k ;a)G (~k ,~k ,~r;t),
4 Z (2π)6 1 2 4 1 2 2
1.5
(10)
1
whereG (~k ,~k ,~r;t)isdefinedasthecorrelationfunction
4 1 2 0.5
ofthemicroscopicself-intermediatescatteringfunctionat 0 1
time t and calculated for different wave vectors,
-0.5
0.5 0 cos(θ)
1
G (~k ,~k ,~r;t)= V hFˆ (~k ;t)Fˆ (~k ;t)δ[~r−~r (0)]i, r 1.5 2 -1
4 1 2 N2 n 1 m 2 nm
X
n6=m
(11)
and f(k;a)=4πa2j (ka)/k with j denoting a spherical
1 1
Bessel function of the first kind.
The present work is mostly concerned with the
IImm[[GG((kk,,rr;;tt))]]
anisotropicnatureofdynamicheterogeneities,whichcan 44
0.8
be monitored using the four-point correlation function 0.6
given by Eq. (5). In this context we would like to em- 0.4
0.2
phasize that in principle the more general function (11)
0
isalsoanisotropic. However,anytraceofthisanisotropy -0.2
is lost after the integration over wave vectors ~k and ~k -0.4 1
1 2 -0.6
and thus the overlap correlation function (8) is, by con- -0.8
struction, isotropic. 0.5 1 0 cos(θ)
1.5
r 2 -1
B. Anisotropy of G4(~k,~r;t)
FIG. 1: The real part of thecorrelation function G4(~k,~r;τα)
Since the functions Fˆn(~k;t) are sensitive to displace- (upper figure) and the imaginary part of G4(~k,~r;τα) (lower
mentsofparticlesalongthedirectionof~k,thenG (~k,~r;t) figure) for T=0.45 calculated at theα relaxation time.
4
measures interparticle correlations weighted by the dis-
placements along the vector ~k. Particles which move in
~r/r, and
the direction perpendicular to ~k make a contribution to
G4(~k,~r;t) which is the same as their contribution to the L (k,r;t)= 2n+1 G (~k,~r;t)P (~kˆ·~rˆ)d~rˆ. (13)
pair correlation function g(r). We notice that for t > 0 n 4π Z 4 n
four-point function G (~k,~r;t) is not isotropic, but de-
pends on the angleθ b4etween~r and~k. Shownin the up- IfG4(~k,~r;t)doesnotdependontheanglebetween~kand
~r,thenL (k,r;t)iszeroforallnnotequaltozero. Since
perfigureinFig.1istherealpartG (~k,~r;t)forT =0.45 n
4 there are nonzero real and imaginary parts to G (~k,~r;t)
calculated at t = τ , and the lower figure shows the 4
α
fort>0,thentherearenonzerorealandimaginaryparts
imaginary part. The maximum value of the real part
to L (k,r;t). By symmetry, the imaginary part is zero
of G (~k,~r;τ ) occurs for values of cos(θ) corresponding n
4 α for even n, and the real part is zero for odd n.
to θ = 0◦ and θ = 180◦, which shows that the correla-
Shown in Fig. 2 is the real part L (k,r;τ ) for n = 0
tions are most pronounced for~r paralleland antiparallel n α
and 2, and the imaginary part for n = 1 at the alpha
to~k. Thus,thecorrelationsofthemicroscopicrelaxation
relaxation time τ for T = 0.45. There is a peak in
α
function isanisotropiconthe theαrelaxationtime scale
L (k,r;τ ) and L (k,r;τ ) around the first peak of the
2 α 0 α
and the correlations are the strongest when neighboring
pair correlation function g(r). The dashed lines in the
particles move in the same or in opposite directions. figure are g(r)e−2. Note that due to our definition of
To examine these anisotropic correlations at length the α relaxation time e−2 is the asymptotic limit of the
scales around nearest neighbor distances, we expand
isotropic component L at t=τ , lim L (k,r;τ )=
G4(~k,~r;t) into the Legendre polynomials Fs2(k,τα) = e−2. The0positiveαpeakri→n∞L2(0k,r;ταα) in-
dicates that particles that are initially separated by a
G4(~k,~r;t)= Ln(k,r;t)Pn(~kˆ·~rˆ), (12) distance corresponding to the first peak of g(r) have a
Xn tendency to move in the same direction or in opposite
directions, while the values close to zero around the first
where P is the nth Legendre Polynomial,~kˆ =~k/k,~rˆ= minimum of the static structure factor can result from
n
4
FIG. 3: The time dependence of the the first peak of
L2(k,r;τα) for T = 1.0, 0.9, 0.8, 0.6, 0.55, 0.5, 0.47, 0.45,
shown from left toright.
103
α 102
τ
FIG.2: TherealpartofLn(k,r;τα)forn=0and2,andthe , 2
L
timimaeg.inTarhyefdoarsnhe=d1linfoeriTn=th0e.4fi5gcuarlecsuliastegd(ra)te−th2ewαherreelagx(art)ioins τ, ng 101
τ
thepair correlation function.
100
motion which is perpendicular to the initial separation 10-1
10-2 10-1 100
vector. Thespatialvariationofthe correlatedmotionon (T - T )/T
these length scales has been reported previously in col- c c
loidal suspensions [22] and is related to the break up of
the cage surrounding a particle. FIG.4: ThetimeatwhichthefirstpeakofL2(k,r;t)reaches
The variation of the imaginary part of L (k,r;τ ) in- its maximum value, τL2, (squares) compared to the α relax-
1 α ation time, τα, (circles) and the peak time of the standard
dicatesthatparticlescloserthanthefirstpeakofg(r)are
non-gaussian parameter, τng, (triangles).
more likely to move apart, while particles at a distance
greater than this peak are more likely to move closer to-
gether. Ingeneral,negativevaluesofL (k,r;t) indicates
1
that particles move farther apart while positive values
indicate that particles move closer together. corresponding to the peak position of the standard non-
To look at the time dependence of the anisotropy, Gaussian parameter α (t) = 3hδr4(t)i/[5hδr2(t)i2]− 1,
2
we calculated the height of the first peak of L (k,r;t) τ (triangles in Fig. 4). Furthermore, the maximum
2 ng
as a function of time, which is shown in Fig. 3 for value does not monotonicallyincrease with a decreasein
T = 1.0, 0.9, 0.8, 0.6, 0.55, 0.5, 0.47 and 0.45. The the temperature, but rather reaches a maximum around
peak height starts at zero since the liquid is isotropic, T = 0.55, then begins to decrease with decreasing tem-
thenincreases,reachesamaximum,andfinallydecreases perature. Thus the anisotropy around nearest neighbor
to zero at long times. The height of the peak, τ , distancesinitiallyincreasesuponsupercoolingtheliquid,
L2
is around the α relaxation time for high temperatures, but reaches a maximum and begins to slowly decrease
Fig. 4, but its position increases slower with decreas- when the liquid is cooled further. It is not known if
ing temperature than the α relaxationtime and approx- the peakheightcontinuestodecreaseorsaturatesatlow
imately follows the temperature dependence of the time temperatures.
5
4 8
t = 10 τ T = 0.45 θ = 0ο
α
t = 1.0 τ θ = 90ο T = 0.45
3 t = 0.1 τα 6 θ = 180ο
α )
(k,q;t) 2 t = 0.0 τα τk,q;α4
S4 (4
S
1 2
0 0
0 5 10 0 5 10
q q
FIG.5: Thefour-pointcorrelationfunctionS4(~k,~q;t)forθ= FIG. 6: The four-point correlation function S4(~k,~q;τα) for
0, where θ is the angle between~k and ~q, calculated at t=0, θ =0◦, 90◦, and 180◦, where θ is the angle between ~q and~k
0.1 τα, τα and 10τα for a temperatureT =0.45. calculated for T =0.45.
IV. FOUR-POINT STRUCTURE FACTOR The usual interpretation of results shown in Fig. 5 is
S4(~k,~q;t) thattheincreaseofS4(~k,~q;t)atsmallqvaluessuggestsa
growing dynamic length scale ξ(t). To find the dynamic
A. Anisotropy of S4(~k,~q;t) lengthscale,itiscommontofitthesmallq behaviortoa
functional form and to examine the scaling of S (~k,~q;t)
4
for small q. In such a procedure it is implicitly assumed
To investigate the correlations between microscopic
that S (~k,~q;t) is isotropic.
self-intermediatescatteringfunctions atlargerdistances, 4
we examined the the structure factor corresponding to However, we find that S4(~k,~q;t) is not isotropic and
G (~k,~r;t), depends on the angle between~k and ~q. Shown in Fig. 6
4
is S (~k,~q,τ ) for T =0.45 and for θ =0◦, 90◦, and 180◦
4 α
S (~k,~q;t)=1+ NH (~k,~q;t) (14) where θ is the angle between~k and ~q. The anisotropyof
4 V 4 S (~k,~q;t) adds a complication in finding a unique ξ(t).
4
Since we do not expect any slowly-decaying with
whereH (~k,~q;t)istheFouriertransformofG (~k,~r;t)−
4 4 increasing distance spatial correlations between self-
F2(k;t). For ~q 6=0
s intermediate scattering functions pertaining to different
particles, we can safely assume that the ~q → 0 limit
1
S (~k,~q;t)= hFˆ (~k;t)Fˆ (−~k;t)e−iq~·~rnm(0)i. (15) of S (~k,~q;t) is well defined and it does not depend on
4 n m 4
N
Xn,m the angle between vectors ~q and ~k. However, the re-
sults shown in Fig. 6 suggest that the correlation length
Again, we fix |~k| to be around the position of the first may be anisotropic. We would like to emphasize that
peak of the static structure factor for the A particles, our results are consistent with such a possibility but do
|~k| =7.25. not prove it. To prove that the correlation length is
Functions similar to (15) have been used to examine anisotropic one would need to simulate bigger systems
a growing dynamic length scale in glass forming liquids in order to be able to examine the structure factor at
[4, 9, 10, 11, 16]. In Fig. 5 we show results similar to smaller wave vectors ~q.
those presented in, e.g. Ref. [9]. Specifically, we show in We examine the anisotropyof the four-pointstructure
S4(~k,~qk;t) for T = 0.45 at times t = 0, 0.1τα, τα, and factorbycalculatingtheprojectionofS4(~k,~q;t)ontothe
10τ . Notethatfort=0,S (~k,~q;0)=S(q)whereS(q)is Legendre polynomials,
α 4
the static structure factor for the A particles. We would
liketo emphasize thatresultsshowninFig.5 arefor one I (k,q;t)= 2n+1 S (~k,~q;t)P (~kˆ·~qˆ)d~qˆ. (16)
n 4π Z 4 n
specific angle between~k and ~q; the angle between ~q and
~k is zero. It should be noted that for this angle between ShowninFig.7(a)isI (k,q;τ )(i.e.,theangularaverage
0 α
vectors ~q and ~k, S4 does not depend on time for q = k. of S4) for T = 1.0, 0.8, 0.6, 0.55 and 0.45. In most sim-
Thisfollowsfromdefinition(15);S (~k,~k;t)=S(k)atall ulational studies of four-point correlation functions the
4
times. correlationfunctionsareshownasaveragesoverdifferent
6
6 8
(a) 0.45 T = 0.45
5 0.55
6
0.60
4 0.80 )
τ)α 1.0 q;t04
q; 3 k,
k, (4
( S
0 2
I 2
1 0
0 10-1 100 101 102 103 104 105
0 1 2 3 4 5 t
q
FIG. 8: Time dependence of S4(~k,~q0;t) for different angles
0 between ~k and ~q0 calculated for T = 0.45. The solid lines
◦ ◦ ◦ ◦ ◦
(b) correspond to θ = 0 , 30 , 45 , 60 , and 90 listed from
bottomtotop. Thedashedlinesare120◦,135◦,150◦and180◦
-0.2 listed from top to bottom. The vertical dotted line indicates
) 0 theα relaxation time.
α
τ
I(k,q;2-0.4 τ)I(k,q;α20---000...426
-0.6
-0.8
0.01 0.1 1
(T-T)/T
c c
-0.8
0 1 2 3 4 5
q
FIG. 7: The projections I0(k,q;t),(a), and I2(k,q;t), (b),as
describedinthetext. Theprojection I0(k,q;t)istheaverage
overanglesθ between~k and~q ofS4(~k,~q;t). IfS4(~k,~q;t)does
not depend on θ, then I2(k,q;t) would be zero. Shown in
the inset is I2(k,q0;τα) where q0 is the smallest wave vector
allowed due to periodic boundary conditions as a function of
temperature. The symbols in (a) and (b) correspond to the
same temperatures.
FIG.9: TimedependenceofI2(k,q0;t)whereq0 isthesmall-
est wave vector allowed due to periodic boundary conditions
for T =1.0, 0.9, 0.8, 0.6, 0.55, 0.5, 0.47, and 0.45 listed from
directions of wave vector ~q, thus the results are similar
left to right.
to what is shown in Fig. 7(a). Note, however, that an
averageoverdifferentdirectionsof~q maynotcorrespond
toanangularaverageifthesamenumberofwavevectors
T =0.5, then it remains approximately constant.
correspondingtoeachanglebetween~qand~karenotused
intheaverage. Therefore,differentroutinestodetermine
S (~k,~q;t) can lead to different conclusions, and our re-
su4lts demonstrate that the averagingprocedure needs to B. Time dependence of the anisotropy of S4(~k,~q0;t)
be performed with caution.
We now turn to the examination of the time depen-
Shown in Fig. 7(b) is I (k,q;τ ) for T = 1.0, 0.8, 0.6,
0.55and0.45. Thenon-ze2rovalueαsofI2 isaconsequence dence of the anisotropyof S4(~k,~q;t). To this end, we set
ofS (~k)beinganisotropicontheαrelaxationtimescale. |~q|tobeequaltothesmallestwavevectorallowedforour
4
finite size simulationbox,|~q|=q =2π/L,andcalculate
The anisotropy is largest for the smallest q values. The 0
temperaturedependenceofI2(k,q0;τα)isshownasanin- S4(~k,~q0;t) as a function of time for different angles be-
setto Fig.7(b). The anisotropyattheα relaxationtime tween~k and~q. Results for T =0.45 are shownin Fig. 8,
for q grows with decreasing temperature until around andtheverticallinemarkstheαrelaxationtime. Wesee
0
7
S (~k,~q;t)todifferentfunctionalforms[9, 11]. Lacevicet
4
al. [9] used an Ornstein-Zernicke form A/(1+(ξq)2) to
fit an overlapfunction Sol(q) that is isotropic by design,
103 while Toninelli et al. [114] used (A−C)/(1+(ξq)β)+C
to fit a function similar to the one studied in this work.
102 Lacevic et al. found a correlation length growing with
τα timeuntilthepeaktimeintheassociatedfour-pointsus-
, g ceptibility, and then decreasing. In contrast, Toninelli
, τ2n 101 et al. found a correlation length growing with time even
τI afterthepeakintheassociatedsusceptibility. Itispossi-
100 blethatthedifferencebetweenthesefindingswasrelated
th the presence of the new parameter β in the fit used
by Toninelli et al. More recently, Berthier et al. [6] used
10-1 10-2 10-1 100 A/(1 + (ξq)β) and found a value of β = 2.4 provided
(T-T )/T good fits to the same correlation function studied in the
c c
Ref. [6]. Here we focus on a possible anisotropy of the
correlation length at the time equal to the α relaxation
FIG. 10: The time corresponding to the maximum value of
timeandweleaveitstimedependenceforafuturestudy.
the magnitude of I2(k,q0;t), τI2, (circles) compared to the
We started with
peakposition ofthenon-Gaussianparameterτng (diamonds)
and the αrelaxation time τα (squares). S (~k,0;τ )−C
S (~k,~q;τ )= 4 α +C, (17)
4 α 1+(ξ q)2+(aq)4
θ
that S (~k,~q;t) grows with increasing time, then reaches
4
a maximum that depends on θ for a time around the α as a fitting function to extract the dynamic correlation
relaxation time and finally decays to one at long times. length ξθ. In Eq. (17) we added a constant C because of
Note that, while the position of the maximum is around thegrowingbaselinewhichcanbeseeninFig.5. Wenote
theαrelaxationtime,thespecifictimeatwhichthepeak that since we do not expect any slowly decaying spatial
is reached depends on the angle between~k and ~q. correlations,inthe limitq →0,S4(~k,~q;τα)shouldbein-
To determine the time dependence of the anisotropy, dependent on the angle between~k and ~q. In contrast, in
we examined I (k,q ;t) where q is the smallest wave Eq. (17) we allowed for the dependence of the dynamic
2 0 0
vector allowed due to periodic boundary condition, q0 = correlation length ξθ on the angle θ between ~k and ~q.
2π/L. As seen in Fig. 9, I2(k,q0;t) is zero at short and While fits to Eq. 17 were very good for q < 3, the re-
long times, but develops a peak at intermediate times. sults were not satisfactory. The values of S (~k,0;t) were
4
Note that the shape of I2(k,q0;t) is somewhat similar notconsistentfor different anglesθ between~k and~q and
to that of L (k,r ;t) shown in Fig. 3 except that
2 peak the length scales ξ obtained from the fits were greater
θ
I (k,q ;t) is negative (the last fact could be expected
2 0 than40atthe lowesttemperatures. To solvetheseprob-
fromtherelationbetweenL (k,r;t)andI (k,q;t)). The
2 2 lems we performed the procedure described below. We
peak height increases with decreasing temperature until
emphasize that simulations of larger systems need to be
T = 0.47, where it starts to decrease. However, as we
performed to test this procedure and its results.
show in the next subsection, the correlation length ob-
Initially, we attempted to set C and a to zero, thus
tainedfromthefitsatT =0.45areallclosetoorgreater
fitting functions to the Ornstein-Zernicke form. We set
than half the box length, and it is currently unknown if
a to zerosince itwas alwaysverysmallinthe previously
the decrease in the peak height is a finite size effect.
attempted fitting procedure. If this form is correct,then
To determine when the anisotropy is a maximum one could ideally find S (~k,0;t) by fitting S (~k,~q;t) for
at large distances, we found the time when I (k,q ;t) 4 4
2 0 differentanglesbetween~kand~qundertheconditionthat
reaches its maximum value, τ . Shown in Fig. 10 is
I2
one obtains consistentresults. We did notobtainconsis-
the temperature dependence ofτ (circles)comparedto
I2 tent results for S (~k,0;t) with this procedure and also
τ (squares) and the peak position of the standard non- 4
Gαaussian parameter τng (diamonds). We notice similar found that we needed to fix the value of S4(~k,0;t) to
trends as with the time corresponding to the maximum obtain values of ξθ less than 50. Therefore, to obtain
value of L2(k,rmax;t) in that the τI2 occurs around τng an estimate for S4(~k,0;t), we choose to fit I0(k,q;t) for
and has a similar temperature dependence. q < 1.5 to an Ornstein-Zernicke form and then set the
value of S (~k,0;t) = I (k,0;t) where I (k,0,t) is ob-
4 0 0
tained from the fits. Note that this is consistent with
C. Effective dynamic correlation length our assumption that the limit lim S (~k,~q;t) does not
q~→0 4
depend on the angle between~k and ~q.
There has been some effort to determine the dy- If the glass transition is governed by a growing dy-
namic correlation length by fitting functions similar to namic length scale, then it is expected that for small
8
100
100
(a)
ττI(k,q;)/I(k,0;)αα001100--21 ττI(k,q;)/I(k,0;)αα00111000--201 ττ)S(k,q;)/S(k,0;αα441100--21 ξΘ2468
10-3 10-130-1 100qξiso 101 10-3 0 10-(1T-Tc)/Tc 100
10-1 100 101 10-1 100 101
qξ qξ
iso Θ
FIG. 11: The four smallest wave vectors allowed
due to periodic boundary conditions of the projection 101
I0(k,q;τα)/I0(k,0;τα) versus qξiso for T = 0.45, 0.47, 0.5,
(b)
0.55, 0.6, and 0.8. The solid line is the scaling function
1/(1+x2). The inset shows all calculated wave vectors less
than 5.
o
s
ξi
enough ~q that S4(~k,~q;τα)/S4(~k,0;τα) versus ξq should ξ, Θ
be described by a universal function F(qξ) that is inde-
pendentoftemperature[7]. Tocheckifthisscalingholds 100
for I (k,q;τ ), we plotted I (k,q;τ )/I (k,0;τ ) versus
0 α 0 α 0 α
qξ where I (k,0;t) and ξ are obtained from the fits
iso 0 iso
describedaboveandthescalingfunction1/(1+x2),which 10-1 100 101 102 103
τ
is shownin Fig. 11. The subscriptiso inξ emphasizes α
iso
that this correlation length was obtained from the ori-
entationalaverageI0(k,q;τα) ofthe four-pointstructure FIG. 12: (a) S4(~k,~q;t)/S4(~k,0;t) versus qξθ for θ = 0◦, 45◦,
factorS4(~k,~q;τα). Itappearsthatthis scalingholds well and 90◦ where θ is the angle between~k and ~q calculated for
for the small q values, but we will again caution that the temperatures T =0.8, 0.6, 0.55, 0.5, 0.47, and 0.45. The
simulations of larger systems need to be performed to inset shows the length scales obtained from the different fits
verify this observation. Shown as the inset to the fig- (θ = 0◦, triangles; θ = 45◦, squares; θ = 90◦, circles; ξiso,
ure is I (k,q;τ )/I (k,0;τ ) versus qξ for wave vec- X’s). (b)Dynamiccorrelationlengthsversustheαrelaxation
0 α 0 α iso ◦ ◦ ◦
tors with a magnitude less than five, and the deviation time for θ = 0 (triangles), 45 (squares), and 90 (circles).
The solid lines are fits to ξ ∼ τγ. The X’s and the dashed
from the scaling behavior is obvious for the larger wave α
vectors. ThecorrelationlengthobtainedfromI (k,q;τ ) line corresponds to ξiso obtained from thefitsof I0(k,q;τα).
0 α
is on the order of a particle diameter at the larger tem-
peratures, but grows to about five particle diameters at
T = 0.45. This growth of the correlation length is con- periodic boundary conditions at each temperature and
sistentwithrecentresultsofBerthierandJack[5]. Note, angle. The overlap is very good for the 18 functions
however,thatatthe lowesttemperatureξiso iscompara- shown, and shown in the inset to Fig. 12 are the corre-
ble to the half the length of the simulation cell, which is lationlengths. They depend onthe anglebetween~k and
the largest length we expect to be able to extract from ~q,andthe correlationlengths arelargestforθ =90◦ and
the simulation without finite size effects. smallest for θ = 0◦. Again, we observe that for θ = 90◦,
WiththevaluesofS4(~k,0;τα)fixedusingthe fitsfrom the correlation lengths are larger than half the simula-
I0(~k,0;τα), we fit S4(~k,~q;τα) where the angle θ between tion cell for T = 0.5 (where ξ90 ≈ 4.5) and lower. This
~k and~q are0,45,and90degreestoanOrnstein-Zernicke strongly suggests that already at T = 0.5 simulations of
formwhereonlythecorrelationlengthisallowedtovary. larger systems are needed in order to verify the present
WeshowS (~k,~q;τ )/S (~k,0;τ )versusqξ ,whereξ de- results.
4 α 4 α θ θ
pends on the angle θ between ~k and ~q, for T = 0.8, 0.6, In previous studies it has been found that the correla-
0.55, 0.5, 0.47, and 0.45 in Fig. 12. Only wave vectors tion length is related to the α relaxation time according
with a magnitude less than 1.5 are shown, which corre- to a powerlaw, ξ ∼τγ [6, 9, 28]. Recently, this behavior
α
sponds to the four smallest wave vectors allowed due to was rationalized by the inhomogeneous mode-coupling
9
theory [7]. We fitted the the correlation lengths to a parable to the nearest neighbor distance the anisotropy
power law of the form aτγ and obtained values ranging initially increasesupon supercoolingthe liquid, but then
α
from γ = 0.22±0.01 for θ = 0◦ and γ = 0.18±0.01 for seems to saturate or evendecreaseat the lowesttemper-
θ = 90◦, Fig. 9. Also shown in Fig. 9 is ξ obtained atures. Furthermore, the time scale that this anisotropy
iso
from I (k,q;t); in this case we found γ = 0.21±0.01, is a maximum for nearest neighbor distances is around
0 0
which is very close to the previously reported value of theαrelaxationtimeathighertemperatures,butthenit
0.22, [28]. Using this analysis, we find that the dynamic increases slowerwith decreasing temperature than the α
correlationlengthisnotonlydifferentfordifferentangles relaxationtimeandroughlyfollowsthetimecorrespond-
between ~k and ~q, but they also grow at a different rate ing to the peak position of the non-Gaussian parameter
as the temperature is lowered and the α relaxation time α2(t), τng.
increases. The range of correlations for particles moving For larger distances, we also found anisotropy of the
in the same direction are longer than for particles mov- four-point correlation function. We studied the time
ing in different directions, but it increases slower with dependence of this longer ranged anisotropy and found
decreasing temperature. that the time at which it is the largest also approxi-
Another scaling prediction is that S4(~k,0;τα) ∼ τα∆. mately follows τng in the supercooled liquid. The longer
TotestthispredictionwefitI (k,0;τ )totheformaτ∆. range anisotropy introduces a challenge in determining
0 α α
Inthiswayweobtain∆≈0.37,whichisagainveryclose the growing dynamic length scale ξ in glass forming sys-
tothevalueof0.4reportedinRef.[28]. Thesevaluesare tems. This difficulty is compounded by the relatively
slightly smaller than the recent inhomogeneous mode- smallsystemsizesusuallyemployedinsimulationalstud-
coupling theory prediction of ∆=0.5 [7]. ies of the glass transition. We developed a procedure to
extract effective dynamic length scales, but larger sys-
tem sizes need to be simulated to verify our results. Our
V. CONCLUSIONS procedure suggests that the dynamic correlation length
is different depending on the relative direction of mo-
There have been many studies looking for a growing tion of two particles within the fluid. Furthermore, this
length scale that accompanies the drastic slowing down anisotropic length scale also increases at a different rate
ofthedynamicsinsupercooledandglassformingliquids. with decreasing temperature.
Recently, one suchpossibility wasexaminedby Biroliet. We hope that our present work will stimulate future
al [29]wheretheyassociatedagrowingcorrelationlength research in two different directions. First, we advocate
with a point-to-set correlation function in a model su- the need to study larger systems and to perform serious
percooled liquid. Ever since the observationof heteroge- finite-size analysis of the results [30, 31, 32]. In par-
neousdynamicsinsupercooledandglassysystems,ithas ticular, we expect that in the small ~q limit four-point
beensuggestedthatadynamiccorrelationlengthmaybe structurefactorS (~k,~q;t)isisotropicandwethusweex-
4
associatedwiththesizeofthedynamicallyheterogeneous pectitsanisotropiccomponentI (k,q;t)tovanishinthe
2
regions. Since two point correlation functions are inad- small~q limit. Theseexpectationsshouldbeconfirmedby
equate to describe the correlated motion of atoms and simulations of larger systems. Second, we hope that this
correlated relaxation of the fluid, four-point correlation workwillstimulateadevelopmentofatheoreticalmodel
functions have been developed to examine this cooper- that describes the anisotropy of four-point dynamic cor-
ative motion. Normally these correlation functions are relations.
assumedtobeisotropic,orareisotropicbydesign. How-
ever, it has been observed that correlated displacements
ofparticlesarenotisotropic,andthusitisnotsurprising
Acknowledgments
that the four-point correlation functions might also not
be isotropic.
Inthisworkweexaminedtheanisotropyofafour-point We gratefully acknowledge the support of NSF Grant
correlation function. We found that for distances com- No. CHE 0517709.
[1] M. Ediger, Annu.Rev.Phys. Chem. 51, 99 (2000). [7] G. Biroli, J. Bouchaud, K. Miyazaki, and D. Reichman,
[2] R.Richert, J. Phys.: Condens. Matter 14, R703 (2002). Phys. Rev.Lett. 97, 195701 (2006).
[3] H. Andersen, Proc. Natl. Acad. Sci. U.S.A. 102, 6686 [8] C. Donati, J. Douglas, W. Kob, S. Plimpton, P. Poole,
(2005). and S.Glotzer, Phys. Rev.Lett. 80, 2338 (1998).
[4] L. Berthier, Phys. Rev.E 69, 020201(R) (2004). [9] N. Lacevic, F. Starr, T. Schroder, and S. Glotzer, J.
[5] L.BerthierandR.Jack,Phys.Rev.E76,041509(2007). Chem. Phys.119, 7372 (2003).
[6] L. Berthier, G. Biroli, J. Bouchaud, W. Kob, [10] N. Lacevic, F. Starr, T. Schroder, V. Novikov, and
K. Miyazaki, and D. Reichman, J. Chem. Phys. 126, S. Glotzer, Phys. Rev.E 66, 030101(R) (2002).
184503 (2007). [11] C. Toninelli, M. Wyart, L. Berthier, G. Biroli, and
10
J. Bouchaud, Phys.Rev. E 71, 041505 (2005). [22] E. Weeks and D. Weitz, Phys. Rev. Lett. 89, 095704
[12] L. Berthier, G. Biroli, J. Bouchaud, L. Cipelletti, (2002).
D. Masri, D. L’Hote, F. Ladieu, and M. Pierno, Science [23] W.KobandH.Andersen,Phys.Rev.E51,4626(1995).
310, 1797 (2006). [24] W.KobandH.Andersen,Phys.Rev.E52,4134(1995).
[13] C. Dalle-Ferrier, C. Thibierge, C. Alba-Simionesco, [25] S. Abraham and B. Bagchi, Phys. Rev. E 78, 051501
L.Berthier,G.Biroli,J.Bouchaud,F.Ladieu,D.L’Hote, (2008).
and G. Tarjus, Phys. Rev.E 76, 041510 (2007). [26] T. Abete, A. de Candia, E. Del Gado, A. Fierro, and
[14] F. Lechenault, O. Dauchot, G. Biroli, and J. Bouchaud, A. Coniglio, Phys. Rev.E 78, 041404 (2008).
Europhys.Lett. 83, 46003 (2008). [27] InRef.[10]bothg4oldefinedinEq.(8)andamoregeneral,
[15] G.BiroliandJ.Bouchaud,Europhys.Lett.67,21(2004). collective version of thisfunction were discussed.
[16] L. Berthier, G. Biroli, J. Bouchaud, W. Kob, [28] S. Whitelam, L. Berthier, and J. Garrahan, Phys. Rev.
K. Miyazaki, and D. Reichman, J. Chem. Phys. 126, Lett. 92, 185705 (2004).
184503 (2006). [29] G. Biroli, J. Bouchaud, A. Cavagna, T. Grigera, and
[17] G. Szamel, Phys. Rev.Lett. 101, 205701 (2008). P. Verrochio, NaturePhys.4, 771 (2008).
[18] M.IwataandS.Sasa,Europhys.Lett.77,50008 (2007). [30] R. Stein, Ph.D.thesis, Stanford University(2007).
[19] E. Flenner and G. Szamel, J. Phys.: Condens. Matter [31] L. Berthier, Phys.Rev. Lett. 91, 055701 (2003).
19, 205125 (2007). [32] R. Stein and H. C. Andersen, Phys. Rev. Lett. 101,
[20] B. Doliwa and A. Heuer,Phys. Rev.E 61, 6898 (2000). 267802 (2008).
[21] Y. Gebremichael, M. Vogel, and S. Glotzer, J. Chem.
Phys.120, 4415 (2004).