APS/123-QED
Simulations of Aging and Plastic Deformation in Polymer Glasses
Mya Warren∗ and Jörg Rottler
Department of Physics and Astronomy, The University of British Columbia,
6224 Agricultural Road, Vancouver, BC, V6T 1Z1, Canada
(Dated: August 14, 2018)
We study the effect of physical aging on the mechanical properties of a model polymer glass
using molecular dynamics simulations. The creep compliance is determined simultaneously with
the structural relaxation under a constant uniaxial load below yield at constant temperature. The
model successfully captures universal features found experimentally in polymer glasses, including
signatures of mechanical rejuvenation. We analyze microscopic relaxation timescales and show
that they exhibit the same aging characteristics as the macroscopic creep compliance. In addition,
our model indicates that the entire distribution of relaxation times scales identically with age.
Despite large changes in mobility, we observe comparatively little structural change except for a
weak logarithmic increase in the degree of short-range order that may be correlated to an observed
decrease in aging with increasing load.
PACS numbers: 81.40.Lm, 81.40.Lg, 83.10.Rs
I. INTRODUCTION
Glassy materials are unable to reach equilibrium over
typical experimental timescales [1, 2, 3]. Instead, the
presence of disorder at temperatures below the glass tran-
sition permits only a slow exploration of the configura-
tional degrees of freedom. The resulting structural re-
laxation, also known as physical aging [4], is one of the
hallmarks of glassy dynamics and leads to material prop-
erties that depend on the wait time tw since the glass was
formed. While thermodynamic variables such as energy
and density typically evolve only logarithmically, the re-
laxation times grow much more rapidly with wait time
[3, 4, 5].
Aging is a process observed in many different glassy
systems, including colloidal glasses [6], microgel pastes
[7], and spin glasses [8], but is most frequently studied
in polymers due to their good glass-forming ability and
ubiquitous use in structural applications. Of particular
interest is therefore to understand the effect of aging on
their mechanical response during plastic deformation [5].
In a classic series of experiments, Struik [4] studied many
different polymer glasses and determined that their stiff-
ness universally increases with wait time. However, it has
also been found that large mechanical stimuli can alter
the intrinsic aging dynamics of a glass. Cases of both
decreased aging (rejuvenation) [4] and increased aging
(overaging) [9, 10] have been observed, but the interpre-
tation of these findings in terms of the structural evolu-
tion remains controversial [11, 12].
The formulation of a comprehensive molecular model
of the non-equilibrium dynamics of glasses has been im-
peded by the fact that minimal structural change oc-
curs during aging. Traditional interpretations of aging
presume that structural relaxation is accompanied by a
∗Electronic address: mya@phas.ubc.ca
decrease in free volume available to molecules and an as-
sociated reduction in molecular mobility [4]. While this
idea is intuitive, it suffers from several limitations. First,
the free volume has been notoriously difficult to define ex-
perimentally. Also, this model does not seem compatible
with the observed aging in glassy solids under constant
volume conditions [13], and cannot predict the aging be-
havior under complex thermo-mechanical histories. Mod-
ern energy landscape theories describe the aging process
as a series of hops between progressively deeper traps in
configuration space [14, 15]. These models have had some
success in capturing experimental trends, but have yet to
directly establish a connection between macroscopic ma-
terial response and the underlying molecular level pro-
cesses. Recent efforts to formulate a molecular theory of
aging are promising but require knowledge of how local
density fluctuations control the relaxation times in the
glass [16].
Molecular simulations using relatively simple models
of glass forming solids, such as the binary Lennard-Jones
glass [17] or the bead spring model [18] for polymers,
have shown rich aging phenomenology. For instance,
calculations of particle correlation functions have shown
explicitly that the characteristic time scale for particle
relaxations increases with wait time [19]. Recent work
[13, 20] has focused on the effect of aging on the mechan-
ical properties; results showed that the shear yield stress
(defined as the overshoot or maximum of the stress-strain
curve) in deformation at constant strain rate generally in-
creases logarithmically with tw. Based on a large number
of simulations at different strain rates and temperatures,
a phenomenological rate-state model was developed that
describes the combined effect of rate and age on the shear
yield stress for many temperatures below the glass tran-
sition [21].
In contrast to the strain-controlled studies described
above, experiments on aging typically impose a small,
constant stress and measure the resulting creep as a func-
tion of time and tw [4]. In this study, we perform molecu-
http://arxiv.org/abs/0704.1840v1
mailto:mya@phas.ubc.ca
lar dynamics simulations on a coarse grained, glass form-
ing polymer model in order to investigate the relation-
ship between macroscopic creep response and microscopic
structure and dynamics. In Section IIIA, we determine
creep compliance curves for different temperatures and
applied loads (in the sub-yield regime) and find that, as
in experiments, curves for different ages can be super-
imposed by rescaling time. The associated shift factors
exhibit a power-law dependence on the wait time, and
the effect of aging can be captured by an effective time
as originally envisioned by Struik [4]. In Section III B, we
compute microscopic mobilities and the full spectrum of
relaxation times and show their relationship to the creep
response. Additionally, we study several parameters that
are sensitive to the degree of short-range order in Section
III C. We detect very little evolution toward increased
local order in our polymer model, indicating that short
range order is not a sensitive measure of the mechanical
relaxation times responsible for the creep compliance of
glassy polymers.
II. SIMULATIONS
We perform molecular dynamics (MD) simulations
with a well-known model polymer glass on the bead-
spring level. The beads interact via a non-specific van
der Waals interaction given by a 6-12 Lennard-Jones po-
tential, and the covalent bonds are modeled with a stiff
spring that prevents chain crossing [22]. This level of
modeling does not include chemical specificity, but al-
lows us to study longer aging times than a fully atomistic
model and seems appropriate to examine a universal phe-
nomenon found in all glassy polymers. All results will be
given in units of the diameter a of the bead, the mass m,
and the Lennard-Jones energy scale, u0. The characteris-
tic timescale is therefore τLJ =
ma2/u0, and the pres-
sure and stress are in units of u0/a
3. The Lennard-Jones
interaction is truncated at 1.5a and adjusted vertically
for continuity. All polymers have a length of 100 beads,
and unless otherwise noted, we analyze 870 polymers in a
periodic simulation box. Results are obtained either with
one large simulation containing the full number of poly-
mers, or with several smaller simulations, each starting
from a unique configuration, whose results are averaged.
The large simulations and the averaged small simulations
provide identical results. The small simulations are used
to estimate uncertainties caused by the finite size of the
simulation volume.
To create the glass, we begin with a random dis-
tribution of chains and relax in an ensemble at con-
stant volume and at a melt temperature of 1.2u0/kB.
Once the system is fully equilibrated, it is cooled over
750τLJ to a temperature below the glass transition at
Tg ≈ 0.35u0/kB [18]. The density of the melt is cho-
sen such that after cooling the pressure is zero. We
then switch to an NPT ensemble - the pressure and
temperature are controlled via a Nosé-Hoover thermo-
  750
 2250
 7500
22500
75000
  500
 1500
 5000
15000
50000
FIG. 1: Simulated creep compliance J(t, tw) at a glassy tem-
perature of T = 0.2u0/kB for various wait times tw (indi-
cated in the legend in units of τLJ ). A uniaxial load of (a)
σ = 0.4u0/a
3 and (b) σ = 0.5u0/a
3 is applied to the aged
glasses. The strain during creep is monitored over time to
give the creep compliance.
stat/barostat - with zero pressure and age for various
wait times (tw) between 500 to 75,000τLJ. The aged
samples undergo a computer creep experiment where a
uniaxial tensile stress (in the z-direction) is ramped up
quickly over 75τLJ , and then held constant at a value
of σ, while the strain ǫ = ∆Lz/Lz is monitored. After
an initial elastic deformation, the glass slowly elongates
in the direction of applied stress due to structural relax-
ations. In the two directions perpendicular to the applied
stress, the pressure is maintained at zero.
III. RESULTS
A. Macroscopic Mechanical Deformation
Historically, measurements of the creep compliance
have been instrumental in probing the relaxation dynam-
ics of glasses, and continue to be the preferred tool in
investigating the aging of glassy polymers [15, 23, 24]. In
his seminal work on aging in polymer glasses, Struik [4]
performed an exhaustive set of creep experiments on dif-
ferent materials, varying the temperature and the applied
load. In this section, we perform a similar set of exper-
iments with our model polymer glass. The macroscopic
creep compliance is defined as
J(t, tw) =
ǫ(t, tw)
. (1)
short time long time
FIG. 2: The same data as Fig. 1 is shown with the curves
shifted by aJ (tw) to form a master curve. The dashed lines are
fits to the master curves using the effective time formulation,
and the dotted line is a short-time fit for comparison (see
text).
Compliance curves J(t, tw) for several temperatures and
stresses were obtained as a function of wait time since
the quench; representative data is shown in Figure 1.
The curves for different wait times appear similar and
agree qualitatively with experiment. An initially rapid
rise in compliance crosses over into a slower, logarith-
mic increase at long times. The crossover between the
two regimes increases with increasing wait time. Struik
showed that experimental creep compliance curves for
different ages can be superimposed by rescaling the time
variable by a shift factor, aJ ,
J(t, tw) = J(taJ , t
w). (2)
This result is called time-aging time superposition [4, 5].
Simulated creep compliance curves from Fig. 1 can simi-
larly be superimposed, and the resulting master curve is
shown in Fig. 2.
Shift factors required for this data collapse are plot-
ted versus the wait time in Fig. 3. All data fall along a
straight line in the double-logarithmic plot, clearly indi-
cating power law behavior:
w . (3)
This power law in the shift factor is characteristic of ag-
ing. µ is called the aging exponent, and has been found
experimentally to be close to unity for a wide variety of
glasses in a temperature range near Tg [4].
Figure 4 shows the effect of stress and temperature on
the aging exponent, as determined from linear fits to the
data in Fig. 3. At T = 0.2u0/kB, µ is close to one for
small stresses, but decreases strongly with stress. This
apparent erasure of aging by large mechanical deforma-
tions has been called “mechanical rejuvenation” [25]. Ex-
periments have frequently found a stress dependence of
3 4 5
3 4 5
σ = 0.2
σ = 0.4
σ = 0.5
σ = 0.1
σ = 0.2
σ = 0.4
T = 0.2 T = 0.3
FIG. 3: Plot of the shift factors found by superimposing the
creep compliance curves, aJ (circles), the mean-squared dis-
placement curves, aMSD (triangles), and the incoherent scat-
tering function curves, aC (×) at different wait times (see
text). The solid lines are linear fits to the data.
the aging exponent [4], although it is not always the case
that the aging process slows down with applied stress;
stress has been known to increase the rate of aging in
some circumstances as well [9, 10]. The structural ori-
gins of this effect are not well understood [11, 12].
At T = 0.3u0/kB, we find that the aging exponent is
somewhat smaller than at T = 0.2u0/kB and varies much
less with applied stress. This behavior is most likely
due to the fact that the temperature is approaching Tg.
Indeed, experiments show that µ rapidly drops to zero
above Tg. compliance is an order of magnitude larger at
T = 0.3u0/kB than at T = 0.2u0/kB and the data does
not fully superimpose in a master curve for long times
where J > 0.2u0/a
3. Shift factors were obtained from
the small creep portion of the curves.
The relatively simple relationship between shift factors
and wait time permits construction of an expression that
describes the entire master curve in Fig. 2. For creep
times that are short compared to the wait time - such
that minimal physical aging occurs over the timescale of
the experiment - experimental creep compliance curves
can be fit to a stretched exponential (typical of processes
with a spectrum of relaxation times),
J(t) = J0 exp[(t/t0)
m] (4)
where t0 is the retardation factor, and the exponent, m,
has been found to be close to 1/3 for most glasses [4].
A fit of this expression to our simulated creep compli-
ance curves is shown in Fig. 2 (dotted line). This ex-
0 0.1 0.2 0.3 0.4 0.5 0.6
T = 0.2
T = 0.3
FIG. 4: The aging exponent, µ, determined from the slopes
of log(aJ ) versus log(tw) (from Fig. 3) plotted versus stress
(open symbols). The solid symbols at zero stress refer to
shift factors determined from aMSD (eq. 7) and aC (eq. 6)
data only. The dashed lines are guides to the eye.
pression is clearly only consistent with the data at times
t < tw. At times much longer than the wait time, the
creep compliance varies more slowly due to the stiffen-
ing caused by aging during the course of the experiment.
Struik suggested that eq. (4) could be extended to the
long-time creep regime, where the experimental timescale
may be longer than the wait time, by introducing an ef-
fective time to account for the slowdown in the relaxation
timescales:
teff =
tw + t′
dt′ (5)
Upon replacing t with teff , eq. (4) may be used to de-
scribe the entire creep curve. Creep compliance curves
from Fig. 2 can indeed be fit to this form (dashed lines)
for a known wait time, tw, and aging exponent, µ, as
obtained from the master curve. We find m ≈ 0.5 ± 0.1
for all stresses at T = 0.2u0/kB, and a relatively broad
range of values for J0 and t0 are consistent with the data.
For the simple thermo-mechanical history prescribed by
the creep experiment, Struik’s effective time formulation
appears to work quite well.
The present results parallel those of a recent simula-
tion study of the shear yield stress in glassy solids [21].
In this work, the glassy solid was deformed at constant
strain rate, and two different regimes of strong and weak
rate dependence emerged depending on the time to reach
the yield point relative to the wait time. In order to ra-
tionalize these results, a rate-state model was developed
that accounted for the internal evolution of the material
with age through a single state variable Θ(t). This formu-
lation successfully collapses yield stress data for different
ages and strain rates in a universal curve by adapting the
evolution of the state variable during the strain interval.
We note here that this state variable is closely related to
Struik’s effective time, as it tries to subsume the modified
aging dynamics during deformation in a single variable
and in particular easily includes the effects of overaging
or rejuvenation.
B. Microscopic Dynamics
The aging behavior of the simulated mechanical re-
sponse functions agrees remarkably well with experiment.
Additional microscopic information from simulations al-
lows us to obtain more directly the relevant timescales
of the system, and the relevant microscopic processes re-
sponsible for aging. One parameter which has been use-
ful in studying glassy dynamics is the “self” part of the
incoherent scattering factor [19],
Cq(t, tw) =
exp(i~q · [~rj(tw + t)− ~rj(tw)]) (6)
where ~rj is the position of the j
th atom, and ~q is the
wave-vector. Cq curves as a function of age are shown
in Fig. 5 and exhibit three distinct regions. Initially,
Cq decreases as particles make very small unconstrained
excursions about their positions. There follows a long
plateau, where the correlation function does not change
considerably. In this regime, atoms are not free to dif-
fuse, but are trapped in local cages formed by their near-
est neighbours. For this reason, the time spent in the
plateau regime is often associated with a “cage time”.
The plateau region ends when particles finally escape
from local cages (α-relaxation), and larger atomic rear-
rangements begin to take place. The cage time corre-
sponds closely to the transition from short-time to long-
time regime observed in the creep compliance. Structural
rearrangements taking place in the α-relaxation regime
are clearly associated with the continued aging observed
in the creep compliance, as well as plastic deformations
occurring in that region.
The correlation functions for different ages are similar
in form, but the time spent in the plateau region in-
creases with age. Just as creep compliance curves can be
shifted in time to form a master curve, we may overlap
the long-time, cage-escape regions of Cq by rescaling the
time variable of the correlation data at different ages (see
inset of Fig. 5). The corresponding shift factors aC(tw)
are also shown in Fig. 3, where we see that the increase
in cage time with age follows the same power law as the
shift factors of the creep compliance. These results are
qualitatively similar to the scaling of the relaxation times
with age found in [19] with no load.
The real space quantity corresponding to Cq is the
mean squared displacement,
〈~r(t, tw)
∆~rj(t, tw)
2 (7)
FIG. 5: Incoherent scattering factor (eq. 6) for different
wait times measured under the same loading conditions as
in Fig. 1(a) for q = (0, 0, 2π). The inset shows the master
curve created by rescaling the time variable by aC . Symbols
as in Fig. 1(b).
FIG. 6: Mean-squared displacement (eq. 7) for different wait
times measured under the same loading conditions as in
Fig. 1(a). The inset shows the master curve created by rescal-
ing the time variable by aMSD. Symbols as in Fig. 1(b).
where ∆~rj(t, tw) = ~rj(tw + t) − ~rj(tw). This function
is shown in Fig. 6. Again we see three characteristic re-
gions of unconstrained (ballistic), caged, and cage-escape
behavior. The departure from the cage plateau likewise
increases with age, and a master curve can be constructed
by shifting the curves with a factor aMSD (see inset of
Fig. 6). Shift factors aMSD are plotted in Fig. 3, along
with shifts for creep compliance and incoherent scatter-
ing function. As anticipated, the shifts versus wait time
fully agree with those obtained from Cq and J .
This clearly demonstrates that for our model, the cage
escape time is indeed the controlling factor in the aging
dynamics of the mechanical response functions.
0 0.5 1 1.5 2 2.5 3 3.5
FIG. 7: The displacement probability distribution versus time
measured under the same loading conditions as in Fig. 1(a),
with a wait time of 500τLJ . The solid lines from left to right
are obtained at times t of 75, 750, 7500, and 75000τLJ . The
dashed lines show fits to the double Gaussian distribution (see
text, eq. 8).
Additional information about microscopic processes
can be obtained by studying not only the mean of the dis-
placements, but also the full spectrum of relaxation dy-
namics as a function of time and wait time. To this end,
we measure the probability distribution P (∆r(t, tw)
2) of
atomic displacements during time intervals, t, for glasses
at various ages, tw. This quantity is complementary to
the measurements of dynamical heterogeneities detailed
in [26], where the spatial variations of the vibrational am-
plitudes were measured at a snapshot in time to show the
correlations of mobile particles in space. In our study, we
omit the spatial information, but retain all of the dynam-
ical information.
Representative distribution functions are shown in
Fig. 7 for a constant wait time of tw = 500τLJ and var-
ious time intervals t. The distributions were obtained
from a smaller system of only 271 polymer chains due to
memory constraints. The data does not reflect a simple
Gaussian distribution, but clearly exhibits the presence
of two distinct scales: there is a narrow distribution of
caged particles and a wider distribution of particles that
have escaped from their cages. This behavior can be de-
scribed by the sum of two Gaussian peaks,
P (∆r2) = N1 exp
+N2 exp
. (8)
Deviations from purely Gaussian behavior are common in
glassy systems and are a signature of dynamical hetero-
geneities [26, 27]. Experiments on colloidal glasses [28]
show a similar separation of displacement distributions
into fast and slow particles.
A fit of the normalized distributions to eq. (8) (dashed
lines in Fig. 7) requires adjustment of three parameters:
the variance of caged and mobile particle distributions,
1 2 3 4 5
FIG. 8: The Gaussian fit parameters for the distribution of
displacements (see text, eq. 8), (a) N1/N , (b) σ
1 , and (c) σ
measured under the same loading conditions as in Fig. 1(a).
The curves are for wait times increasing from left to right
from 500τLJ to 15000τLJ .
and σ2
, as well as their relative contributions N1/N ,
where N = N1 +N2. These parameters are sufficient to
describe the full evolution of the displacement distribu-
tion during aging. In Fig. 8, we show the fit parameters
as a function of time and wait time. Again two distinct
time scales are evident. At short times, most of the parti-
cles are caged (N1/N ≈ 1), and the variance of the cage
peak is also changing very little. There are few rear-
rangements in this regime, however Fig. 8(c) shows that
a small fraction of particles are mobile at even the short-
est times. At a time corresponding to the onset of cage
escape, the number of particles in the cage peak begins
to rapidly decay, and the variance of the cage peak in-
creases. This indicates that the cage has become more
malleable - small, persistent rearrangements occur lead-
ing to eventual cage escape. In this regime, the variance
of the mobile peak increases very little. Note that the
typical length scale of rearrangements is less than a par-
ticle diameter even in the cage escape regime, but the
number of particles undergoing rearrangements changes
by more than 50%.
Similar to the compliance and mean-squared displace-
ment curves, the data in Fig. 8(a) and (b) can also be
0 0.5 1 1.5 2
  500
 1500
 5000
15000
FIG. 9: The displacement probability distribution measured
under the same loading conditions and wait times as in
Fig. 1(a) plotted at times corresponding to < r2(t, tw) >=
0.7, shown in the inset as a dashed line. The legend indicates
the wait time.
superimposed by shifting time. Shift factors for N1/N
and σ2
coincide exactly with shifts for the mean; how-
ever, data for σ2
(Fig. 8(c)) seems to be much less af-
fected by the wait time. The aging dynamics appears to
be entirely determined by the cage dynamics, and not by
larger rearrangements within the glass.
Since the fit parameters exhibit the same scaling with
wait time as the mean, one might expect that the en-
tire distribution of displacements under finite load scales
with the evolution of the mean. In Fig. 9, we plot dis-
placement distributions for several wait times at time in-
tervals chosen such that the mean squared displacements
are identical (see inset). Indeed, we find that all curves
overlap, indicating that the entire relaxation spectrum
ages in the same way. A similar observation was recently
made in simulations of a model for a metallic glass aging
at zero stress [29], although in this study the tails of the
distribution were better described by stretched exponen-
tials.
In order to study the effect of load on the relaxation
dynamics, we compare in Figure 10 the fit parameters for
a sample undergoing creep (replotted from Fig. 8) and a
reference sample without load. It is clear that the dy-
namics are strongly affected by the applied stress only in
the region characterized by α-relaxations. For the stress
applied here, the onset of cage-escape does not appear to
be greatly modified by the stress, however the decay in
N1/N and the widening of the cage peak are accelerated.
The stress does not modify the variance of the mobile
peak, confirming again the importance of local rearrange-
ments as compared to large-scale motion in the dynamics
of the system. The accelerated structural rearrangements
caused by the stress result in creep on the macroscopic
scale, but may also be responsible for the modification of
the aging dynamics with stress as observed in Fig. 4.
1 2 3 4 5
σ = 0
σ = 0.4 (a)
FIG. 10: The Gaussian fit parameters to the displacement
distributions (see text, eq. 8) (a) N1/N , (b) σ
1 , and (c) σ
a sample aged at T = 0.2u0/kB for tw = 500τLJ , and then
either undergoing a creep experiment at σ = 0.4u0/a
3 (black),
or simply aging further at zero stress (red).
C. Structural evolution
The connection between the dynamics and the struc-
ture of a glass during aging remains uncertain, mostly
because no structural parameter has been found that
strongly depends on wait time. Recent simulation stud-
ies of metallic glasses have shown the existence of sev-
eral short range order parameters that can distinguish
between glassy states created through different quench-
ing protocols [30, 31, 32]. A strong correlation has been
found between “ordered” regions of the glass and strain
localization. Many metallic glasses are known to form
quasi-crystalline structures that optimize local packing.
It remains to be seen whether the short-range order
evolves in the context of aging and in other glass for-
mers such as polymers and colloids. A recent experimen-
tal study of aging in colloidal glasses found no change in
the distribution function of a tetrahedral order parame-
ter [33]. Below, we investigate several measures of local
order in our model as they evolve with age and under
load.
Since Lennard-Jones liquids are known to condense
into a crystal with fcc symmetry at low temperatures,
it is reasonable to look for the degree of local fcc order in
our polymer model. The level of fcc order can be quan-
tified via the bond orientational parameter [34],
. (9)
This parameter has been successfully used to character-
ize the degree of order in systems of hard sphere glasses.
Q6 is determined for each atom by projecting the bond
angles of the nearest neighbours onto the spherical har-
monics, Y6m(θ, φ). The overbar denotes an average over
all bonds. Nearest neighbours are defined as all atoms
within a cutoff radius, rc, of the central atom. For all of
the order parameters discussed here, the cutoff radius is
defined by the first minimum in the pair correlation func-
tion, in this case 1.45a. Q6 is approximately 0.575 for a
perfect fcc crystal; for jammed structures, it can exhibit
a large range of values less than about 0.37 [34]. The full
distribution of Q6 for our model glass is shown for several
ages as well as an initial melt state in Fig. 11(a). We see
that there is very little difference even between melt and
glassy states in our model, and no discernible difference
at all with increasing age.
Locally, close-packing is achieved by tetrahedral order-
ing and not fcc ordering, however, tetrahedral orderings
cannot span the system. The glass formation process has
been described in terms of frustration between optimal
local and global close-packing structures. To investigate
the type of local ordering in this model, we investigate
a 3-body angular correlation function, P (θ). θ is de-
fined as the internal angle created by a central atom and
individual pairs of nearest-neighbours, and P (θ) is the
probability of occurrence of θ. Results for this corre-
lation are shown in Fig. 11(b). Two peaks at approxi-
mately 60◦ and 110◦ indicate tetrahedral ordering. The
peaks sharpen under quenching from the melt, but the
distribution does not evolve significantly during aging.
In contrast, simulations of metallic glass formers showed
a stronger sensitivity of this parameter to the quench
protocol [31], but most of those changes may be due to
rearrangements in the supercooled liquid state and not
in the glassy state.
Another parameter that has been successful in classi-
fying glasses is the triangulated surface order parameter
[32],
(6 − q)νq (10)
which measures the degree of quasi-crystalline order. The
surface coordination number, q, is defined for each atom
of the coordination shell as the number of neighbouring
atoms also residing in the coordination shell; νq is the
number of atoms in the coordination shell with surface
coordination q. Ordered systems have been identified
with S equal to 12 (icosahedron), 14, 15 and 16. Figure
11(c) shows the probability distribution for P (S) for the
melt and for glassy states with short and long wait times.
0.1 0.2 0.3 0.4 0.5
0 50 100 150
θ (degrees)
5 10 15 20 25 30 35
FIG. 11: Short-range order parameters: (a) the bond-
orientational parameter, (b) the three-body angular corre-
lations, (c) the surface triangulated order (see text for dis-
cussion). x’s show the melt state, circles show the sample
aged for tw = 500τLJ , and triangles show the sample aged for
500, 000τLJ .
The peak of the distribution moves toward lower S (more
ordered) upon cooling, and continues to evolve slowly in
the glass. The mean of S relative to the as-quenched
state, 〈S〉, is shown in Fig. 12 as a function of wait time
at two temperatures. We see that 〈S〉 is a logarithmi-
cally decreasing function of wait time. Even though this
is not a strong dependence, this order parameter is sig-
nificantly more sensitive to age than the others that have
been investigated.
Figure 12 also shows the order parameter 〈S〉 after the
ramped-up stress has been applied to the aged samples.
We can see that at T = 0.2u0/kB, some of the order
that developed during age is erased, while no appreciable
change occurs at the higher temperature T = 0.3u0/kB.
These results correlate well with the behavior of the ag-
ing exponent found in Fig. 4, where mechanical rejuve-
nation was found at lower temperatures but was much
less pronounced at higher T . The load is applied very
quickly, and most of the deformation in this regime is
affine, however, the strain during this time was similar
for both temperatures, therefore the effect is not simply
due to a change in density. More work is needed to clarify
the nature of the structural changes during rejuvenation.
T = 0.2 T = 0.3
FIG. 12: Precent change in the triangulated surface order
parameter with wait time as compared to the just-quenched
sample. Circles are for samples aged at zero pressure for the
time tw. Triangles are for the same samples immediately after
ramping up to the creep stress. For T = 0.2u0/kB this stress
is 0.4u0/a
3, and for T = 0.3u0/kB the stress is 0.1u0/a
IV. CONCLUSIONS
We investigate the effects of aging on both macroscopic
creep response and underlying microscopic structure and
dynamics through simulations on a simple model polymer
glass. The model qualitatively reproduces key experi-
mental trends in the mechanical behavior of glasses under
sustained stress. We observe a power-law dependence of
the relaxation time on the wait time with an aging expo-
nent of approximately unity, and a decrease in the aging
exponent with increasing load that indicates the presence
of mechanical rejuvenation. The model creep compliance
curves can be fit in the short and long-time regimes using
Struik’s effective time formulation. Additionally, inves-
tigation of the microscopic dynamics through two-time
correlation functions has shown that, for our model glass,
the aging dynamics of the creep compliance exactly corre-
sponds to the increase in the cage escape or α-relaxation
time.
A detailed study of the entire distribution of parti-
cle displacements yields an interesting picture of the mi-
croscopic dynamics during aging. The distribution can
be described by the sum of two Gaussians, reflecting
the presence of caged and mobile particles. The frac-
tion of mobile particles and the amplitude of rearrange-
ments in the cage strongly increase at the cage escape
time. However, in analogy with results in colloidal glasses
[35], structural rearrangements occur even for times well
within the “caged” regime, and fairly independent of wait
time and stress. For our model glass, we find that the
entire distribution of displacements scales with age in
the same way as the mean. At times when the long-
time portion of the mean squared displacement overlaps,
the distribution of displacements at different wait times
completely superimpose, confirming that all of the me-
chanical relaxation times scale in the same way with age.
To characterize the evolution of the structure during
aging, we investigate several measures of short-range or-
der in our model glass. We find that the short-range order
does not evolve strongly during aging. The triangulated
surface order [32], however, shows a weak logarithmic de-
pendence on age. Results also show a change in structure
when a load is rapidly applied, and this seems to be cor-
related with an observed decrease in the aging exponent
under stress.
This study has characterized the dynamics of a model
glass prepared by a rapid quench below Tg, followed by
aging at constant T and subsequent application of a con-
stant load. For such simple thermo-mechanical histo-
ries, existing phenomenological models work well, how-
ever, the dynamics of glasses are in general much more
complex. For instance, large stresses in the non-linear
regime modify the aging dynamics and cause nontrivial
effects such as mechanical rejuvenation and over-aging
[10, 11]. Also, experiments have shown that the time-age
time superposition no longer holds when polymer glasses
undergo more complex thermal histories such as a tem-
perature jump [23]. The success of our study in analyzing
simple aging situations indicates that the present simula-
tion methodology will be able to shed more light on these
topics in the near future.
Acknowledgments
We thank the Natural Sciences and Engineering Coun-
cil of Canada (NSERC) for financial support. Computing
time was provided by WestGrid. Simulations were per-
formed with the LAMMPS molecular dynamics package
[36].
[1] J.-L. Barrat, M. V. Feigelman, J. Kurchan, and J. D.
(Eds.), Slow Relaxations and Nonequilibrium Dynam-
ics in Condensed Matter: Les Houches Session LXXVII
(Springer, 2002).
[2] C. A. Angell, Science 267, 1924 (1995).
[3] C. A. Angell, K. L. Ngai, G. B. McKenna, P. F. McMil-
lan, and S. W. Martin, J. Applied Physics 88, 3113
(2000).
[4] L. C. E. Struik, Physical Aging in Amorphous Polymers
and Other Materials (Elsevier, Amsterdam, 1978).
[5] J. M. Hutchinson, Prog. Polym. Sci. 20, 703 (1978).
[6] B. Abou, D. Bonn, and J. Meunier, Phys. Rev. E 64,
021510 (2001).
[7] M. Cloitre, R. Borrega, and L. Leibler, Phys. Rev. Lett.
85, 4819 (2000).
[8] T. Jonsson, J. Mattsson, C. Djurberg, F. A. Khan,
P. Nordblad, and P. Svedlindh, Phys. Rev. Lett. 75, 4138
(1995).
[9] V. Viasnoff and F. Lequeux, Phys. Rev. Lett. 89, 065701
(2002).
[10] D. J. Lacks and M. J. Osborne, Phys. Rev. Lett. 93,
255501 (2004).
[11] G. B. McKenna, J. Phys.: Condens. Matter 15, S737
(2003).
[12] L. C. E. Struik, Polymer 38, 4053 (1997).
[13] M. Utz, P. G. Debenedetti, and F. H. Stillinger, Phys.
Rev. Lett. 84, 1471 (2000).
[14] C. Monthus and J. P. Bouchaud, J. Phys. A: Math. Gen.
29, 3847 (1996).
[15] P. Sollich, F. Lequeux, P. Hebraud, and M. E. Cates,
Phys. Rev. Lett. 78, 2020 (1997).
[16] K. Chen and K. S. Schweizer (2007), preprint, personal
communication.
[17] W. Kob, J. Phys.: Condens. Mat. 11, R85 (1999).
[18] C. Bennemann, W. Paul, K. Binder, and B. Duenweg,
Phys. Rev. E 57, 843 (1998).
[19] W. Kob and J.-L. Barrat, Phys. Rev. Lett. 78, 4581
(1997).
[20] F. Varnik, L. Bocquet, and J.-L. Barrat, J. Chem. Phys
120, 2788 (2004).
[21] J. Rottler and M. O. Robbins, Phys. Rev. Lett. 95,
225504 (2005).
[22] K. Kremer and G. S. Grest, J. Chem. Phys 92, 5057
(1990).
[23] H. Montes, V. Viasnoff, S. Juring, and F. Lequeux, J.
Stat. Mech.: Theory and Experiment p. P03003 (2006).
[24] L. C. Brinson and T. S. Gates, Int. J. Solids and Struc-
tures 32, 827 (1995).
[25] G. B. McKenna and A. J. Kovacs, Polym. Eng. Sci. 24,
1131 (1984).
[26] K. Vollmayr-Lee and A. Zippelius, Phys. Rev. E 72,
041507 (2005).
[27] W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and
S. C. Glotzer, Phys. Rev. Lett. 79, 2827 (1997).
[28] E. R. Weeks and D. A. Weitz, Chem. Phys. 284, 361
(2002).
[29] H. E. Castillo and A. Parsaeian, Nature Physics 3, 26
(2007).
[30] Y. Shi and M. L. Falk, Phys. Rev. Lett. 95, 095502
(2005).
[31] F. Albano and M. L. Falk, J. Chem. Phys 122, 154508
(2005).
[32] Y. Shi and M. L. Falk, Phys. Rev. B 73, 214201 (2006).
[33] G. C. Cianci, R. E. Courtland, and E. R. Weeks, in
Flow Dynamics: The Second International Conference
on Flow Dynamics (AIP Conf. Proc., 2006), vol. 832,
pp. 21–25.
[34] S. Torquato, T. M. Truskett, and P. G. Debenedetti,
Phys. Rev. Lett. 84, 2064 (2000).
[35] R. E. Courtland and E. R. Weeks, J. Phys.:Condens.
Matter 15, S359 (2003).
[36] http://lammps.sandia.gov.

  We study the effect of physical aging on the mechanical properties of a model
polymer glass using molecular dynamics simulations. The creep compliance is
determined simultaneously with the structural relaxation under a constant
uniaxial load below yield at constant temperature. The model successfully
captures universal features found experimentally in polymer glasses, including
signatures of mechanical rejuvenation. We analyze microscopic relaxation
timescales and show that they exhibit the same aging characteristics as the
macroscopic creep compliance. In addition, our model indicates that the entire
distribution of relaxation times scales identically with age. Despite large
changes in mobility, we observe comparatively little structural change except
for a weak logarithmic increase in the degree of short-range order that may be
correlated to an observed decrease in aging with increasing load.

APS/123-QED
Simulations of Aging and Plastic Deformation in Polymer Glasses
Mya Warren∗ and Jörg Rottler
Department of Physics and Astronomy, The University of British Columbia,
6224 Agricultural Road, Vancouver, BC, V6T 1Z1, Canada
(Dated: August 14, 2018)
We study the effect of physical aging on the mechanical properties of a model polymer glass
using molecular dynamics simulations. The creep compliance is determined simultaneously with
the structural relaxation under a constant uniaxial load below yield at constant temperature. The
model successfully captures universal features found experimentally in polymer glasses, including
signatures of mechanical rejuvenation. We analyze microscopic relaxation timescales and show
that they exhibit the same aging characteristics as the macroscopic creep compliance. In addition,
our model indicates that the entire distribution of relaxation times scales identically with age.
Despite large changes in mobility, we observe comparatively little structural change except for a
weak logarithmic increase in the degree of short-range order that may be correlated to an observed
decrease in aging with increasing load.
PACS numbers: 81.40.Lm, 81.40.Lg, 83.10.Rs
I. INTRODUCTION
Glassy materials are unable to reach equilibrium over
typical experimental timescales [1, 2, 3]. Instead, the
presence of disorder at temperatures below the glass tran-
sition permits only a slow exploration of the configura-
tional degrees of freedom. The resulting structural re-
laxation, also known as physical aging [4], is one of the
hallmarks of glassy dynamics and leads to material prop-
erties that depend on the wait time tw since the glass was
formed. While thermodynamic variables such as energy
and density typically evolve only logarithmically, the re-
laxation times grow much more rapidly with wait time
[3, 4, 5].
Aging is a process observed in many different glassy
systems, including colloidal glasses [6], microgel pastes
[7], and spin glasses [8], but is most frequently studied
in polymers due to their good glass-forming ability and
ubiquitous use in structural applications. Of particular
interest is therefore to understand the effect of aging on
their mechanical response during plastic deformation [5].
In a classic series of experiments, Struik [4] studied many
different polymer glasses and determined that their stiff-
ness universally increases with wait time. However, it has
also been found that large mechanical stimuli can alter
the intrinsic aging dynamics of a glass. Cases of both
decreased aging (rejuvenation) [4] and increased aging
(overaging) [9, 10] have been observed, but the interpre-
tation of these findings in terms of the structural evolu-
tion remains controversial [11, 12].
The formulation of a comprehensive molecular model
of the non-equilibrium dynamics of glasses has been im-
peded by the fact that minimal structural change oc-
curs during aging. Traditional interpretations of aging
presume that structural relaxation is accompanied by a
∗Electronic address: mya@phas.ubc.ca
decrease in free volume available to molecules and an as-
sociated reduction in molecular mobility [4]. While this
idea is intuitive, it suffers from several limitations. First,
the free volume has been notoriously difficult to define ex-
perimentally. Also, this model does not seem compatible
with the observed aging in glassy solids under constant
volume conditions [13], and cannot predict the aging be-
havior under complex thermo-mechanical histories. Mod-
ern energy landscape theories describe the aging process
as a series of hops between progressively deeper traps in
configuration space [14, 15]. These models have had some
success in capturing experimental trends, but have yet to
directly establish a connection between macroscopic ma-
terial response and the underlying molecular level pro-
cesses. Recent efforts to formulate a molecular theory of
aging are promising but require knowledge of how local
density fluctuations control the relaxation times in the
glass [16].
Molecular simulations using relatively simple models
of glass forming solids, such as the binary Lennard-Jones
glass [17] or the bead spring model [18] for polymers,
have shown rich aging phenomenology. For instance,
calculations of particle correlation functions have shown
explicitly that the characteristic time scale for particle
relaxations increases with wait time [19]. Recent work
[13, 20] has focused on the effect of aging on the mechan-
ical properties; results showed that the shear yield stress
(defined as the overshoot or maximum of the stress-strain
curve) in deformation at constant strain rate generally in-
creases logarithmically with tw. Based on a large number
of simulations at different strain rates and temperatures,
a phenomenological rate-state model was developed that
describes the combined effect of rate and age on the shear
yield stress for many temperatures below the glass tran-
sition [21].
In contrast to the strain-controlled studies described
above, experiments on aging typically impose a small,
constant stress and measure the resulting creep as a func-
tion of time and tw [4]. In this study, we perform molecu-
http://arxiv.org/abs/0704.1840v1
mailto:mya@phas.ubc.ca
lar dynamics simulations on a coarse grained, glass form-
ing polymer model in order to investigate the relation-
ship between macroscopic creep response and microscopic
structure and dynamics. In Section IIIA, we determine
creep compliance curves for different temperatures and
applied loads (in the sub-yield regime) and find that, as
in experiments, curves for different ages can be super-
imposed by rescaling time. The associated shift factors
exhibit a power-law dependence on the wait time, and
the effect of aging can be captured by an effective time
as originally envisioned by Struik [4]. In Section III B, we
compute microscopic mobilities and the full spectrum of
relaxation times and show their relationship to the creep
response. Additionally, we study several parameters that
are sensitive to the degree of short-range order in Section
III C. We detect very little evolution toward increased
local order in our polymer model, indicating that short
range order is not a sensitive measure of the mechanical
relaxation times responsible for the creep compliance of
glassy polymers.
II. SIMULATIONS
We perform molecular dynamics (MD) simulations
with a well-known model polymer glass on the bead-
spring level. The beads interact via a non-specific van
der Waals interaction given by a 6-12 Lennard-Jones po-
tential, and the covalent bonds are modeled with a stiff
spring that prevents chain crossing [22]. This level of
modeling does not include chemical specificity, but al-
lows us to study longer aging times than a fully atomistic
model and seems appropriate to examine a universal phe-
nomenon found in all glassy polymers. All results will be
given in units of the diameter a of the bead, the mass m,
and the Lennard-Jones energy scale, u0. The characteris-
tic timescale is therefore τLJ =
ma2/u0, and the pres-
sure and stress are in units of u0/a
3. The Lennard-Jones
interaction is truncated at 1.5a and adjusted vertically
for continuity. All polymers have a length of 100 beads,
and unless otherwise noted, we analyze 870 polymers in a
periodic simulation box. Results are obtained either with
one large simulation containing the full number of poly-
mers, or with several smaller simulations, each starting
from a unique configuration, whose results are averaged.
The large simulations and the averaged small simulations
provide identical results. The small simulations are used
to estimate uncertainties caused by the finite size of the
simulation volume.
To create the glass, we begin with a random dis-
tribution of chains and relax in an ensemble at con-
stant volume and at a melt temperature of 1.2u0/kB.
Once the system is fully equilibrated, it is cooled over
750τLJ to a temperature below the glass transition at
Tg ≈ 0.35u0/kB [18]. The density of the melt is cho-
sen such that after cooling the pressure is zero. We
then switch to an NPT ensemble - the pressure and
temperature are controlled via a Nosé-Hoover thermo-
  750
 2250
 7500
22500
75000
  500
 1500
 5000
15000
50000
FIG. 1: Simulated creep compliance J(t, tw) at a glassy tem-
perature of T = 0.2u0/kB for various wait times tw (indi-
cated in the legend in units of τLJ ). A uniaxial load of (a)
σ = 0.4u0/a
3 and (b) σ = 0.5u0/a
3 is applied to the aged
glasses. The strain during creep is monitored over time to
give the creep compliance.
stat/barostat - with zero pressure and age for various
wait times (tw) between 500 to 75,000τLJ. The aged
samples undergo a computer creep experiment where a
uniaxial tensile stress (in the z-direction) is ramped up
quickly over 75τLJ , and then held constant at a value
of σ, while the strain ǫ = ∆Lz/Lz is monitored. After
an initial elastic deformation, the glass slowly elongates
in the direction of applied stress due to structural relax-
ations. In the two directions perpendicular to the applied
stress, the pressure is maintained at zero.
III. RESULTS
A. Macroscopic Mechanical Deformation
Historically, measurements of the creep compliance
have been instrumental in probing the relaxation dynam-
ics of glasses, and continue to be the preferred tool in
investigating the aging of glassy polymers [15, 23, 24]. In
his seminal work on aging in polymer glasses, Struik [4]
performed an exhaustive set of creep experiments on dif-
ferent materials, varying the temperature and the applied
load. In this section, we perform a similar set of exper-
iments with our model polymer glass. The macroscopic
creep compliance is defined as
J(t, tw) =
ǫ(t, tw)
. (1)
short time long time
FIG. 2: The same data as Fig. 1 is shown with the curves
shifted by aJ (tw) to form a master curve. The dashed lines are
fits to the master curves using the effective time formulation,
and the dotted line is a short-time fit for comparison (see
text).
Compliance curves J(t, tw) for several temperatures and
stresses were obtained as a function of wait time since
the quench; representative data is shown in Figure 1.
The curves for different wait times appear similar and
agree qualitatively with experiment. An initially rapid
rise in compliance crosses over into a slower, logarith-
mic increase at long times. The crossover between the
two regimes increases with increasing wait time. Struik
showed that experimental creep compliance curves for
different ages can be superimposed by rescaling the time
variable by a shift factor, aJ ,
J(t, tw) = J(taJ , t
w). (2)
This result is called time-aging time superposition [4, 5].
Simulated creep compliance curves from Fig. 1 can simi-
larly be superimposed, and the resulting master curve is
shown in Fig. 2.
Shift factors required for this data collapse are plot-
ted versus the wait time in Fig. 3. All data fall along a
straight line in the double-logarithmic plot, clearly indi-
cating power law behavior:
w . (3)
This power law in the shift factor is characteristic of ag-
ing. µ is called the aging exponent, and has been found
experimentally to be close to unity for a wide variety of
glasses in a temperature range near Tg [4].
Figure 4 shows the effect of stress and temperature on
the aging exponent, as determined from linear fits to the
data in Fig. 3. At T = 0.2u0/kB, µ is close to one for
small stresses, but decreases strongly with stress. This
apparent erasure of aging by large mechanical deforma-
tions has been called “mechanical rejuvenation” [25]. Ex-
periments have frequently found a stress dependence of
3 4 5
3 4 5
σ = 0.2
σ = 0.4
σ = 0.5
σ = 0.1
σ = 0.2
σ = 0.4
T = 0.2 T = 0.3
FIG. 3: Plot of the shift factors found by superimposing the
creep compliance curves, aJ (circles), the mean-squared dis-
placement curves, aMSD (triangles), and the incoherent scat-
tering function curves, aC (×) at different wait times (see
text). The solid lines are linear fits to the data.
the aging exponent [4], although it is not always the case
that the aging process slows down with applied stress;
stress has been known to increase the rate of aging in
some circumstances as well [9, 10]. The structural ori-
gins of this effect are not well understood [11, 12].
At T = 0.3u0/kB, we find that the aging exponent is
somewhat smaller than at T = 0.2u0/kB and varies much
less with applied stress. This behavior is most likely
due to the fact that the temperature is approaching Tg.
Indeed, experiments show that µ rapidly drops to zero
above Tg. compliance is an order of magnitude larger at
T = 0.3u0/kB than at T = 0.2u0/kB and the data does
not fully superimpose in a master curve for long times
where J > 0.2u0/a
3. Shift factors were obtained from
the small creep portion of the curves.
The relatively simple relationship between shift factors
and wait time permits construction of an expression that
describes the entire master curve in Fig. 2. For creep
times that are short compared to the wait time - such
that minimal physical aging occurs over the timescale of
the experiment - experimental creep compliance curves
can be fit to a stretched exponential (typical of processes
with a spectrum of relaxation times),
J(t) = J0 exp[(t/t0)
m] (4)
where t0 is the retardation factor, and the exponent, m,
has been found to be close to 1/3 for most glasses [4].
A fit of this expression to our simulated creep compli-
ance curves is shown in Fig. 2 (dotted line). This ex-
0 0.1 0.2 0.3 0.4 0.5 0.6
T = 0.2
T = 0.3
FIG. 4: The aging exponent, µ, determined from the slopes
of log(aJ ) versus log(tw) (from Fig. 3) plotted versus stress
(open symbols). The solid symbols at zero stress refer to
shift factors determined from aMSD (eq. 7) and aC (eq. 6)
data only. The dashed lines are guides to the eye.
pression is clearly only consistent with the data at times
t < tw. At times much longer than the wait time, the
creep compliance varies more slowly due to the stiffen-
ing caused by aging during the course of the experiment.
Struik suggested that eq. (4) could be extended to the
long-time creep regime, where the experimental timescale
may be longer than the wait time, by introducing an ef-
fective time to account for the slowdown in the relaxation
timescales:
teff =
tw + t′
dt′ (5)
Upon replacing t with teff , eq. (4) may be used to de-
scribe the entire creep curve. Creep compliance curves
from Fig. 2 can indeed be fit to this form (dashed lines)
for a known wait time, tw, and aging exponent, µ, as
obtained from the master curve. We find m ≈ 0.5 ± 0.1
for all stresses at T = 0.2u0/kB, and a relatively broad
range of values for J0 and t0 are consistent with the data.
For the simple thermo-mechanical history prescribed by
the creep experiment, Struik’s effective time formulation
appears to work quite well.
The present results parallel those of a recent simula-
tion study of the shear yield stress in glassy solids [21].
In this work, the glassy solid was deformed at constant
strain rate, and two different regimes of strong and weak
rate dependence emerged depending on the time to reach
the yield point relative to the wait time. In order to ra-
tionalize these results, a rate-state model was developed
that accounted for the internal evolution of the material
with age through a single state variable Θ(t). This formu-
lation successfully collapses yield stress data for different
ages and strain rates in a universal curve by adapting the
evolution of the state variable during the strain interval.
We note here that this state variable is closely related to
Struik’s effective time, as it tries to subsume the modified
aging dynamics during deformation in a single variable
and in particular easily includes the effects of overaging
or rejuvenation.
B. Microscopic Dynamics
The aging behavior of the simulated mechanical re-
sponse functions agrees remarkably well with experiment.
Additional microscopic information from simulations al-
lows us to obtain more directly the relevant timescales
of the system, and the relevant microscopic processes re-
sponsible for aging. One parameter which has been use-
ful in studying glassy dynamics is the “self” part of the
incoherent scattering factor [19],
Cq(t, tw) =
exp(i~q · [~rj(tw + t)− ~rj(tw)]) (6)
where ~rj is the position of the j
th atom, and ~q is the
wave-vector. Cq curves as a function of age are shown
in Fig. 5 and exhibit three distinct regions. Initially,
Cq decreases as particles make very small unconstrained
excursions about their positions. There follows a long
plateau, where the correlation function does not change
considerably. In this regime, atoms are not free to dif-
fuse, but are trapped in local cages formed by their near-
est neighbours. For this reason, the time spent in the
plateau regime is often associated with a “cage time”.
The plateau region ends when particles finally escape
from local cages (α-relaxation), and larger atomic rear-
rangements begin to take place. The cage time corre-
sponds closely to the transition from short-time to long-
time regime observed in the creep compliance. Structural
rearrangements taking place in the α-relaxation regime
are clearly associated with the continued aging observed
in the creep compliance, as well as plastic deformations
occurring in that region.
The correlation functions for different ages are similar
in form, but the time spent in the plateau region in-
creases with age. Just as creep compliance curves can be
shifted in time to form a master curve, we may overlap
the long-time, cage-escape regions of Cq by rescaling the
time variable of the correlation data at different ages (see
inset of Fig. 5). The corresponding shift factors aC(tw)
are also shown in Fig. 3, where we see that the increase
in cage time with age follows the same power law as the
shift factors of the creep compliance. These results are
qualitatively similar to the scaling of the relaxation times
with age found in [19] with no load.
The real space quantity corresponding to Cq is the
mean squared displacement,
〈~r(t, tw)
∆~rj(t, tw)
2 (7)
FIG. 5: Incoherent scattering factor (eq. 6) for different
wait times measured under the same loading conditions as
in Fig. 1(a) for q = (0, 0, 2π). The inset shows the master
curve created by rescaling the time variable by aC . Symbols
as in Fig. 1(b).
FIG. 6: Mean-squared displacement (eq. 7) for different wait
times measured under the same loading conditions as in
Fig. 1(a). The inset shows the master curve created by rescal-
ing the time variable by aMSD. Symbols as in Fig. 1(b).
where ∆~rj(t, tw) = ~rj(tw + t) − ~rj(tw). This function
is shown in Fig. 6. Again we see three characteristic re-
gions of unconstrained (ballistic), caged, and cage-escape
behavior. The departure from the cage plateau likewise
increases with age, and a master curve can be constructed
by shifting the curves with a factor aMSD (see inset of
Fig. 6). Shift factors aMSD are plotted in Fig. 3, along
with shifts for creep compliance and incoherent scatter-
ing function. As anticipated, the shifts versus wait time
fully agree with those obtained from Cq and J .
This clearly demonstrates that for our model, the cage
escape time is indeed the controlling factor in the aging
dynamics of the mechanical response functions.
0 0.5 1 1.5 2 2.5 3 3.5
FIG. 7: The displacement probability distribution versus time
measured under the same loading conditions as in Fig. 1(a),
with a wait time of 500τLJ . The solid lines from left to right
are obtained at times t of 75, 750, 7500, and 75000τLJ . The
dashed lines show fits to the double Gaussian distribution (see
text, eq. 8).
Additional information about microscopic processes
can be obtained by studying not only the mean of the dis-
placements, but also the full spectrum of relaxation dy-
namics as a function of time and wait time. To this end,
we measure the probability distribution P (∆r(t, tw)
2) of
atomic displacements during time intervals, t, for glasses
at various ages, tw. This quantity is complementary to
the measurements of dynamical heterogeneities detailed
in [26], where the spatial variations of the vibrational am-
plitudes were measured at a snapshot in time to show the
correlations of mobile particles in space. In our study, we
omit the spatial information, but retain all of the dynam-
ical information.
Representative distribution functions are shown in
Fig. 7 for a constant wait time of tw = 500τLJ and var-
ious time intervals t. The distributions were obtained
from a smaller system of only 271 polymer chains due to
memory constraints. The data does not reflect a simple
Gaussian distribution, but clearly exhibits the presence
of two distinct scales: there is a narrow distribution of
caged particles and a wider distribution of particles that
have escaped from their cages. This behavior can be de-
scribed by the sum of two Gaussian peaks,
P (∆r2) = N1 exp
+N2 exp
. (8)
Deviations from purely Gaussian behavior are common in
glassy systems and are a signature of dynamical hetero-
geneities [26, 27]. Experiments on colloidal glasses [28]
show a similar separation of displacement distributions
into fast and slow particles.
A fit of the normalized distributions to eq. (8) (dashed
lines in Fig. 7) requires adjustment of three parameters:
the variance of caged and mobile particle distributions,
1 2 3 4 5
FIG. 8: The Gaussian fit parameters for the distribution of
displacements (see text, eq. 8), (a) N1/N , (b) σ
1 , and (c) σ
measured under the same loading conditions as in Fig. 1(a).
The curves are for wait times increasing from left to right
from 500τLJ to 15000τLJ .
and σ2
, as well as their relative contributions N1/N ,
where N = N1 +N2. These parameters are sufficient to
describe the full evolution of the displacement distribu-
tion during aging. In Fig. 8, we show the fit parameters
as a function of time and wait time. Again two distinct
time scales are evident. At short times, most of the parti-
cles are caged (N1/N ≈ 1), and the variance of the cage
peak is also changing very little. There are few rear-
rangements in this regime, however Fig. 8(c) shows that
a small fraction of particles are mobile at even the short-
est times. At a time corresponding to the onset of cage
escape, the number of particles in the cage peak begins
to rapidly decay, and the variance of the cage peak in-
creases. This indicates that the cage has become more
malleable - small, persistent rearrangements occur lead-
ing to eventual cage escape. In this regime, the variance
of the mobile peak increases very little. Note that the
typical length scale of rearrangements is less than a par-
ticle diameter even in the cage escape regime, but the
number of particles undergoing rearrangements changes
by more than 50%.
Similar to the compliance and mean-squared displace-
ment curves, the data in Fig. 8(a) and (b) can also be
0 0.5 1 1.5 2
  500
 1500
 5000
15000
FIG. 9: The displacement probability distribution measured
under the same loading conditions and wait times as in
Fig. 1(a) plotted at times corresponding to < r2(t, tw) >=
0.7, shown in the inset as a dashed line. The legend indicates
the wait time.
superimposed by shifting time. Shift factors for N1/N
and σ2
coincide exactly with shifts for the mean; how-
ever, data for σ2
(Fig. 8(c)) seems to be much less af-
fected by the wait time. The aging dynamics appears to
be entirely determined by the cage dynamics, and not by
larger rearrangements within the glass.
Since the fit parameters exhibit the same scaling with
wait time as the mean, one might expect that the en-
tire distribution of displacements under finite load scales
with the evolution of the mean. In Fig. 9, we plot dis-
placement distributions for several wait times at time in-
tervals chosen such that the mean squared displacements
are identical (see inset). Indeed, we find that all curves
overlap, indicating that the entire relaxation spectrum
ages in the same way. A similar observation was recently
made in simulations of a model for a metallic glass aging
at zero stress [29], although in this study the tails of the
distribution were better described by stretched exponen-
tials.
In order to study the effect of load on the relaxation
dynamics, we compare in Figure 10 the fit parameters for
a sample undergoing creep (replotted from Fig. 8) and a
reference sample without load. It is clear that the dy-
namics are strongly affected by the applied stress only in
the region characterized by α-relaxations. For the stress
applied here, the onset of cage-escape does not appear to
be greatly modified by the stress, however the decay in
N1/N and the widening of the cage peak are accelerated.
The stress does not modify the variance of the mobile
peak, confirming again the importance of local rearrange-
ments as compared to large-scale motion in the dynamics
of the system. The accelerated structural rearrangements
caused by the stress result in creep on the macroscopic
scale, but may also be responsible for the modification of
the aging dynamics with stress as observed in Fig. 4.
1 2 3 4 5
σ = 0
σ = 0.4 (a)
FIG. 10: The Gaussian fit parameters to the displacement
distributions (see text, eq. 8) (a) N1/N , (b) σ
1 , and (c) σ
a sample aged at T = 0.2u0/kB for tw = 500τLJ , and then
either undergoing a creep experiment at σ = 0.4u0/a
3 (black),
or simply aging further at zero stress (red).
C. Structural evolution
The connection between the dynamics and the struc-
ture of a glass during aging remains uncertain, mostly
because no structural parameter has been found that
strongly depends on wait time. Recent simulation stud-
ies of metallic glasses have shown the existence of sev-
eral short range order parameters that can distinguish
between glassy states created through different quench-
ing protocols [30, 31, 32]. A strong correlation has been
found between “ordered” regions of the glass and strain
localization. Many metallic glasses are known to form
quasi-crystalline structures that optimize local packing.
It remains to be seen whether the short-range order
evolves in the context of aging and in other glass for-
mers such as polymers and colloids. A recent experimen-
tal study of aging in colloidal glasses found no change in
the distribution function of a tetrahedral order parame-
ter [33]. Below, we investigate several measures of local
order in our model as they evolve with age and under
load.
Since Lennard-Jones liquids are known to condense
into a crystal with fcc symmetry at low temperatures,
it is reasonable to look for the degree of local fcc order in
our polymer model. The level of fcc order can be quan-
tified via the bond orientational parameter [34],
. (9)
This parameter has been successfully used to character-
ize the degree of order in systems of hard sphere glasses.
Q6 is determined for each atom by projecting the bond
angles of the nearest neighbours onto the spherical har-
monics, Y6m(θ, φ). The overbar denotes an average over
all bonds. Nearest neighbours are defined as all atoms
within a cutoff radius, rc, of the central atom. For all of
the order parameters discussed here, the cutoff radius is
defined by the first minimum in the pair correlation func-
tion, in this case 1.45a. Q6 is approximately 0.575 for a
perfect fcc crystal; for jammed structures, it can exhibit
a large range of values less than about 0.37 [34]. The full
distribution of Q6 for our model glass is shown for several
ages as well as an initial melt state in Fig. 11(a). We see
that there is very little difference even between melt and
glassy states in our model, and no discernible difference
at all with increasing age.
Locally, close-packing is achieved by tetrahedral order-
ing and not fcc ordering, however, tetrahedral orderings
cannot span the system. The glass formation process has
been described in terms of frustration between optimal
local and global close-packing structures. To investigate
the type of local ordering in this model, we investigate
a 3-body angular correlation function, P (θ). θ is de-
fined as the internal angle created by a central atom and
individual pairs of nearest-neighbours, and P (θ) is the
probability of occurrence of θ. Results for this corre-
lation are shown in Fig. 11(b). Two peaks at approxi-
mately 60◦ and 110◦ indicate tetrahedral ordering. The
peaks sharpen under quenching from the melt, but the
distribution does not evolve significantly during aging.
In contrast, simulations of metallic glass formers showed
a stronger sensitivity of this parameter to the quench
protocol [31], but most of those changes may be due to
rearrangements in the supercooled liquid state and not
in the glassy state.
Another parameter that has been successful in classi-
fying glasses is the triangulated surface order parameter
[32],
(6 − q)νq (10)
which measures the degree of quasi-crystalline order. The
surface coordination number, q, is defined for each atom
of the coordination shell as the number of neighbouring
atoms also residing in the coordination shell; νq is the
number of atoms in the coordination shell with surface
coordination q. Ordered systems have been identified
with S equal to 12 (icosahedron), 14, 15 and 16. Figure
11(c) shows the probability distribution for P (S) for the
melt and for glassy states with short and long wait times.
0.1 0.2 0.3 0.4 0.5
0 50 100 150
θ (degrees)
5 10 15 20 25 30 35
FIG. 11: Short-range order parameters: (a) the bond-
orientational parameter, (b) the three-body angular corre-
lations, (c) the surface triangulated order (see text for dis-
cussion). x’s show the melt state, circles show the sample
aged for tw = 500τLJ , and triangles show the sample aged for
500, 000τLJ .
The peak of the distribution moves toward lower S (more
ordered) upon cooling, and continues to evolve slowly in
the glass. The mean of S relative to the as-quenched
state, 〈S〉, is shown in Fig. 12 as a function of wait time
at two temperatures. We see that 〈S〉 is a logarithmi-
cally decreasing function of wait time. Even though this
is not a strong dependence, this order parameter is sig-
nificantly more sensitive to age than the others that have
been investigated.
Figure 12 also shows the order parameter 〈S〉 after the
ramped-up stress has been applied to the aged samples.
We can see that at T = 0.2u0/kB, some of the order
that developed during age is erased, while no appreciable
change occurs at the higher temperature T = 0.3u0/kB.
These results correlate well with the behavior of the ag-
ing exponent found in Fig. 4, where mechanical rejuve-
nation was found at lower temperatures but was much
less pronounced at higher T . The load is applied very
quickly, and most of the deformation in this regime is
affine, however, the strain during this time was similar
for both temperatures, therefore the effect is not simply
due to a change in density. More work is needed to clarify
the nature of the structural changes during rejuvenation.
T = 0.2 T = 0.3
FIG. 12: Precent change in the triangulated surface order
parameter with wait time as compared to the just-quenched
sample. Circles are for samples aged at zero pressure for the
time tw. Triangles are for the same samples immediately after
ramping up to the creep stress. For T = 0.2u0/kB this stress
is 0.4u0/a
3, and for T = 0.3u0/kB the stress is 0.1u0/a
IV. CONCLUSIONS
We investigate the effects of aging on both macroscopic
creep response and underlying microscopic structure and
dynamics through simulations on a simple model polymer
glass. The model qualitatively reproduces key experi-
mental trends in the mechanical behavior of glasses under
sustained stress. We observe a power-law dependence of
the relaxation time on the wait time with an aging expo-
nent of approximately unity, and a decrease in the aging
exponent with increasing load that indicates the presence
of mechanical rejuvenation. The model creep compliance
curves can be fit in the short and long-time regimes using
Struik’s effective time formulation. Additionally, inves-
tigation of the microscopic dynamics through two-time
correlation functions has shown that, for our model glass,
the aging dynamics of the creep compliance exactly corre-
sponds to the increase in the cage escape or α-relaxation
time.
A detailed study of the entire distribution of parti-
cle displacements yields an interesting picture of the mi-
croscopic dynamics during aging. The distribution can
be described by the sum of two Gaussians, reflecting
the presence of caged and mobile particles. The frac-
tion of mobile particles and the amplitude of rearrange-
ments in the cage strongly increase at the cage escape
time. However, in analogy with results in colloidal glasses
[35], structural rearrangements occur even for times well
within the “caged” regime, and fairly independent of wait
time and stress. For our model glass, we find that the
entire distribution of displacements scales with age in
the same way as the mean. At times when the long-
time portion of the mean squared displacement overlaps,
the distribution of displacements at different wait times
completely superimpose, confirming that all of the me-
chanical relaxation times scale in the same way with age.
To characterize the evolution of the structure during
aging, we investigate several measures of short-range or-
der in our model glass. We find that the short-range order
does not evolve strongly during aging. The triangulated
surface order [32], however, shows a weak logarithmic de-
pendence on age. Results also show a change in structure
when a load is rapidly applied, and this seems to be cor-
related with an observed decrease in the aging exponent
under stress.
This study has characterized the dynamics of a model
glass prepared by a rapid quench below Tg, followed by
aging at constant T and subsequent application of a con-
stant load. For such simple thermo-mechanical histo-
ries, existing phenomenological models work well, how-
ever, the dynamics of glasses are in general much more
complex. For instance, large stresses in the non-linear
regime modify the aging dynamics and cause nontrivial
effects such as mechanical rejuvenation and over-aging
[10, 11]. Also, experiments have shown that the time-age
time superposition no longer holds when polymer glasses
undergo more complex thermal histories such as a tem-
perature jump [23]. The success of our study in analyzing
simple aging situations indicates that the present simula-
tion methodology will be able to shed more light on these
topics in the near future.
Acknowledgments
We thank the Natural Sciences and Engineering Coun-
cil of Canada (NSERC) for financial support. Computing
time was provided by WestGrid. Simulations were per-
formed with the LAMMPS molecular dynamics package
[36].
[1] J.-L. Barrat, M. V. Feigelman, J. Kurchan, and J. D.
(Eds.), Slow Relaxations and Nonequilibrium Dynam-
ics in Condensed Matter: Les Houches Session LXXVII
(Springer, 2002).
[2] C. A. Angell, Science 267, 1924 (1995).
[3] C. A. Angell, K. L. Ngai, G. B. McKenna, P. F. McMil-
lan, and S. W. Martin, J. Applied Physics 88, 3113
(2000).
[4] L. C. E. Struik, Physical Aging in Amorphous Polymers
and Other Materials (Elsevier, Amsterdam, 1978).
[5] J. M. Hutchinson, Prog. Polym. Sci. 20, 703 (1978).
[6] B. Abou, D. Bonn, and J. Meunier, Phys. Rev. E 64,
021510 (2001).
[7] M. Cloitre, R. Borrega, and L. Leibler, Phys. Rev. Lett.
85, 4819 (2000).
[8] T. Jonsson, J. Mattsson, C. Djurberg, F. A. Khan,
P. Nordblad, and P. Svedlindh, Phys. Rev. Lett. 75, 4138
(1995).
[9] V. Viasnoff and F. Lequeux, Phys. Rev. Lett. 89, 065701
(2002).
[10] D. J. Lacks and M. J. Osborne, Phys. Rev. Lett. 93,
255501 (2004).
[11] G. B. McKenna, J. Phys.: Condens. Matter 15, S737
(2003).
[12] L. C. E. Struik, Polymer 38, 4053 (1997).
[13] M. Utz, P. G. Debenedetti, and F. H. Stillinger, Phys.
Rev. Lett. 84, 1471 (2000).
[14] C. Monthus and J. P. Bouchaud, J. Phys. A: Math. Gen.
29, 3847 (1996).
[15] P. Sollich, F. Lequeux, P. Hebraud, and M. E. Cates,
Phys. Rev. Lett. 78, 2020 (1997).
[16] K. Chen and K. S. Schweizer (2007), preprint, personal
communication.
[17] W. Kob, J. Phys.: Condens. Mat. 11, R85 (1999).
[18] C. Bennemann, W. Paul, K. Binder, and B. Duenweg,
Phys. Rev. E 57, 843 (1998).
[19] W. Kob and J.-L. Barrat, Phys. Rev. Lett. 78, 4581
(1997).
[20] F. Varnik, L. Bocquet, and J.-L. Barrat, J. Chem. Phys
120, 2788 (2004).
[21] J. Rottler and M. O. Robbins, Phys. Rev. Lett. 95,
225504 (2005).
[22] K. Kremer and G. S. Grest, J. Chem. Phys 92, 5057
(1990).
[23] H. Montes, V. Viasnoff, S. Juring, and F. Lequeux, J.
Stat. Mech.: Theory and Experiment p. P03003 (2006).
[24] L. C. Brinson and T. S. Gates, Int. J. Solids and Struc-
tures 32, 827 (1995).
[25] G. B. McKenna and A. J. Kovacs, Polym. Eng. Sci. 24,
1131 (1984).
[26] K. Vollmayr-Lee and A. Zippelius, Phys. Rev. E 72,
041507 (2005).
[27] W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and
S. C. Glotzer, Phys. Rev. Lett. 79, 2827 (1997).
[28] E. R. Weeks and D. A. Weitz, Chem. Phys. 284, 361
(2002).
[29] H. E. Castillo and A. Parsaeian, Nature Physics 3, 26
(2007).
[30] Y. Shi and M. L. Falk, Phys. Rev. Lett. 95, 095502
(2005).
[31] F. Albano and M. L. Falk, J. Chem. Phys 122, 154508
(2005).
[32] Y. Shi and M. L. Falk, Phys. Rev. B 73, 214201 (2006).
[33] G. C. Cianci, R. E. Courtland, and E. R. Weeks, in
Flow Dynamics: The Second International Conference
on Flow Dynamics (AIP Conf. Proc., 2006), vol. 832,
pp. 21–25.
[34] S. Torquato, T. M. Truskett, and P. G. Debenedetti,
Phys. Rev. Lett. 84, 2064 (2000).
[35] R. E. Courtland and E. R. Weeks, J. Phys.:Condens.
Matter 15, S359 (2003).
[36] http://lammps.sandia.gov.
