Single-parameter aging in the weakly nonlinear limit

Physical aging deals with slow property changes over time caused by molecular rearrangements. This is relevant for non-crystalline materials like polymers and inorganic glasses, both in production and during subsequent use. The Narayanaswamy theory from 1971 describes physical aging - an inherently nonlinear phenomenon - in terms of a linear convolution integral over the so-called material time $\xi$. The resulting"Tool-Narayanaswamy (TN) formalism"is generally recognized to provide an excellent description of physical aging for small, but still highly nonlinear temperature variations. The simplest version of the TN formalism is single-parameter aging according to which the clock rate $d\xi/dt$ is an exponential function of the property monitored [T. Hecksher et al., J. Chem. Phys. 142, 241103 (2015)]. For temperature jumps starting from thermal equilibrium, this leads to a first-order differential equation for property monitored, involving a system-specific function. The present paper shows analytically that the solution to this equation to first order in the temperature variation has a universal expression in terms of the zeroth-order solution, $R_0(t)$. Numerical data for a binary Lennard-Jones glass former probing the potential energy confirm that, in the weakly nonlinear limit, the theory predicts aging correctly from $R_0(t)$ (which by the fluctuation-dissipation theorem is the normalized equilibrium potential-energy time-autocorrelation function).


I. INTRODUCTION
The properties of non-crystalline materials like polymers and inorganic glasses change slightly over time.In many cases this aging is so slow that it cannot be observed, but sometimes it results in unwanted degradation of material properties.When aging is exclusively due to molecular rearrangements with no chemical reactions involved, one speaks about physical aging [1][2][3][4][5][6][7][8][9][10][11][12][13][14].The present-day understanding of physical aging is based on the century-old observation [15] that any glass is in an out-of-equilibrium state and, as a consequence, relaxes continuously toward the equilibrium state.
During physical aging the system's volume decreases slightly.This reflects the fact that the equilibrium metastable liquid is denser than the glass at the same temperature.Likewise, the enthalpy decreases during aging.Both effects are extremely difficult to observe because they are minute and take place over a very long time.Defining a glass as any non-equilibrium state of a liquid resulting from a thermodynamic perturbation, a good way of studying physical aging is the following: First, equilibrate the glass-forming liquid at some "annealing" temperature just below the calorimetric glass-transition temperature.Depending on the viscosity of the liquid, this may take long time -experiments often study liquids at temperatures at which the equilibrium relaxation time is hours or days [16][17][18][19][20][21].If the equilibrium relaxation time is one day, annealing the sample for a week ensures virtually complete thermal equilibrium.Once the sample has been equilibrated, temperature is changed rapidly to a new value and kept there for long enough time to allow for monitoring virtually the entire equilibration process.This defines an "ideal aging experiment" [17,22].Such an experiment requires the ability to monitor some quantity accurately and continuously as a function of time, excellent temperature control, and the ability to change temperature rapidly [17].Aging may be probed by measuring any property that can be monitored precisely, e.g., the electrical capacitance at a particular frequency [21,[23][24][25][26]; in conjunction with a Peltier-element-based fast and accurate temperature control this is our favorite method in Roskilde [21].Other quantities that have been monitored during physical aging include, e.g., density [27,28], enthalpy [29,30], Young's modulus [31], gas permeability [32], high-frequency mechanical moduli [16,33], dc conductivity [2], X-ray photon correlation spectroscopy [34], and nonlinear dielectric susceptibility [35].
The present paper develops the theory of aging by studying the so-called single-parameter aging framework, which is the simplest realization of the concept of a material time controlling aging in the Tool-Narayanaswamy (TN) formalism [3].An important prediction of the TN formalism is that if the aging rate is known as a function of the property monitored, knowledge of the linear limit of physical aging, e.g., following an infinitesimal temperature jump, is enough to quantitatively determine the aging resulting from any time-dependent temperature variation.According to the fluctuation-dissipation (FD) theorem any linear-response property is determined by thermal-equilibrium fluctuations quantified in terms of a time-autocorrelation function.The prospect for future experimental investigations is that one can make quantitative predictions of aging from a knowledge of equilibrium fluctuations.
Single-parameter aging results in a first-order differential equation for the normalized relaxation function following a temperature jump [18].This equation involves an a priori unknown, system-specific function that determines the linear limit of aging.In order to predict aging, however, it is enough to know the linear-limit relaxation function that, by the FD theorem, is the relevant equilibrium time-autocorrelation function.This paper derives an explicit expression for the weakly nonlinear limit of aging based on the relevant equilibrium time-autocorrelation function.After developing the theory in Sec.II, we give an example of how to calculate a specific quantity similar to the fragility of glass science in Sec.III; this section can be skipped in a first reading of the paper.The general first-order solution to single-parameter aging following a temperature jump is derived in Sec.IV.The validity of the formalism is illustrated in Sec.V by results from computer simulations of a binary Lennard-Jones system; the final section gives a brief discussion.

II. THE TN FORMALISM AND SINGLE-PARAMETER AGING
The quantity probed during aging is denoted by χ(t).Following a temperature jump at t = 0, χ(t) gradually approaches its equilibrium value χ eq at the new temperature T 0 .We define the normalized relaxation function R(t) by By definition, R(t) is unity at t = 0 and approaches zero as t → ∞, i.e., for both temperature up and down jumps R(t) is a decreasing function of time.In practice, in both experiments and simulations there is always a rapid initial change of χ(t) immediately after t = 0 deriving from χ's dependence on the fast, vibrational degrees of freedom and/or one or more fast relaxation processes decoupled from the main (α) relaxation.For this reason, workers in the field often normalize the relaxation function by defining R(t) to be unity after the initial rapid change of χ(t).That is different from what is done in Eq. ( 1), which is our preference because it avoids introducing the extra parameters that comes from estimating the value of the short-time "plateau" of χ(t).
The TN material time is denoted by ξ.This quantity, which may be thought of as the time measured on a clock with a clock rate γ(t) that changes as the material ages, is related to the clock rate as follows dξ = γ(t)dt . ( According to the TN formalism, the material time ξ = ξ(t) controls the physical aging in such a way that the variation of χ(ξ), denoted by is a linear convolution integral over the temperature variation history T (ξ) − T 0 [3,29].Single-parameter aging (SPA) is the simplest version of the TN formalism [18].SPA assumes that the clock rate γ(t) is an exponential function of the monitored property, i.e., Here γ eq is the equilibrium relaxation rate at T 0 and χ 0 is a constant with the same dimension as χ.In conjunction with the TN prediction that physical aging is linear in the temperature variation when formulated in terms of the material time, SPA may be applied to any relatively small (continuous or discontinuous) temperature variation around T 0 , not just to the discontinuous temperature jumps to which the below discussion is limited.Since ∆χ(t) = ∆χ(0)R(t) by the definition of R(t), Eq. ( 4) may be rewritten For temperature jumps the TN fundamental result is that [3,29] R(t) = Φ(ξ) (6) in which the function Φ(ξ) is the same for all temperature jumps of a given system.In view of the nonlinearity of physical aging, this is a highly nontrivial prediction.Keeping in mind the definition of γ(t) (Eq.( 2)), Eq. ( 6) implies Ṙ(t) = Φ (ξ)γ(t).Since according to Eq. ( 6) ξ is the same function of R for all jumps, defining F (R) ≡ −Φ (ξ(R)) leads to the "jump differential equation" in which F (R) is the same for all jumps of a given system.The negative sign in the definition of F (R) is introduced in order to make F (R) positive.Equation (7) has been confirmed in experiments on a silicone oil and several organic liquids [18,20,21] aged to equilibrium just below their calorimetric glass transition temperature.Even though the largest temperature jumps studied were just a few percent, this is enough to exhibit a strongly nonlinear response with more than one decade of relaxation-time variation.One experimental test of Eq. ( 7) involved rewriting it as [18,20] − Ṙ(t) and showing that the left-hand side is the same function of R for different jumps.A second test confirmed the consequence of Eq. ( 7) that R(t) for an arbitrary jump may be predicted from the data for a single jump [18,20,21].
This paper develops the SPA formalism based on Eq. ( 7) that for simplicity is rewritten by adopting the unit system in which γ eq = 1 at the reference temperature T 0 : with

III. CALCULATION OF A GENERALIZED FRAGILITY
Each value of Λ leads to a unique solution denoted by R(t, Λ) of the jump differential equation Eq. ( 9) with the initial condition R(0, Λ) = 1.As a first illustration of how perturbation theory may be applied when |Λ| 1, we determined the Λ dependence of the average relaxation time defined by From τ (Λ) a fragility-like [36] parameter m a (subscript "a" for aging) may be defined by The minus ensures that m a > 0 because Λ > 0 from Eq. ( 9) leads to a faster than equilibrium relaxation.We proceed to derive the following expression in which R 0 (t) ≡ R(t, Λ = 0) Note that whenever 0 < R 0 (t) < 1, which is usually the case [18], one has m a < 1.Note also that since γ eq (Λ) = exp(Λ) from Eq. ( 9), the equilibrium relaxation time τ eq (Λ) ≡ 1/γ eq (Λ) obeys d ln τ eq = −dΛ.Using this one can transform Eq. ( 13) into an expression for how the relative change of τ (Λ) from its value at T 0 depends on the relative change of the equilibrium relaxation time between the two temperatures involved in the jump, i.e., The fact that m a < 1 is now intuitively obvious since the graph of R(t) obviously falls between the equilibrium relaxation function graphs at the two temperatures.
To derive Eq. ( 13), note that Eq. ( 9) implies dt = − exp(−ΛR) dR/F (R).Thus From this we get and By substituting R = R 0 into both integrals and switching back to time as the integration variable one finds An alternative proof of Eq. ( 13) makes use of an integral criterion derived in Ref. 18 (Appendix).
For the calculation of m a from experimental or computer simulation data on R 0 (t) one proceeds as follows.Given a sequence of times (∆t, 2∆t, 3∆t, ..., n∆t) at which the equilibrium normalized relaxation function (R 0,1 , R 0,2 , R 0,3 , ..., R 0,n ) is known, we have As an illustration of the above we consider the case in which the linear-limit relaxation function is a stretched exponential with exponent β (for simplicity, a dimensionless time is used below), Defining the function we have we get If the normalized relaxation function in the experimental time window is described by a stretched exponential with short-time plateau C < 1, i.e., by the function C exp(−t β ), one finds

IV. SOLVING THE JUMP DIFFERENTIAL EQUATION TO FIRST ORDER IN THE TEMPERATURE CHANGE ∆T
To find the solution, R(t, Λ), of the jump differential equation in first-order perturbation theory we proceed as follows.A first-order expansion of R(t), is substituted into Eq.( 9) where R 0 (t) is the normalized relaxation function corresponding to an infinitesimal jump, i.e., to the linear limit of aging (this function is discussed below in Sec.V A).To first order in Λ one has which leads to the following zeroth-and first-order equations: Due to the zero-time normalization of both R(t) and R 0 (t), the initial condition of R 1 is R 1 (0) = 0.For t > 0 one has R 1 (t) < 0 because Λ > 0 as mentioned implies a faster relaxation toward equilibrium, i.e., R(t) < R 0 (t).Consequently, since R 1 (0) = R 1 (t → ∞) = 0, the function R 1 (t) is non-monotonous.
The solution to Eq. ( 28) obeying the initial condition R 1 (0) = 0 is To derive this we proceed as follows.First, we note that the inverse of R(t, Λ) is given by which follows by rewriting Eq. ( 9) as dt = − exp(−ΛR) dR/F (R) and integrating.Next we note that because in which it here and henceforth is understood that all functions are evaluated at Λ = 0, implying that one should put R = R 0 in the final evaluations.Noting that dR 0 /F (R 0 ) = −dt, using Eq. ( 30) the numerator is evaluated as follows For Λ = 0 the denominator of Eq. ( 31) is given by Combining these results one arrives at Eq. ( 29).To confirm that Eq. ( 29) indeed solves Eq. ( 28), one differentiates: Since R0 = −F (R 0 ) Ṙ0 by Eq. ( 27), we see that As an illustration, we show how Eq. ( 29) leads to Eq. ( 13).From Eq. ( 11), Eq. ( 12), and Eq. ( 25) one easily derives This is simplified by performing a partial integration: We thus arrive at Eq. ( 13).

V. NUMERICAL RESULTS FOR A BINARY LENNARD-JONES MODEL
A. The relevant fluctuation-dissipation theorem When the externally controlled variable is temperature itself, a slightly modified derivation of the FD theorem is required [37].In the end, however, the result looks much like in the standard FD case: If ∆β(t) is the variation of β ≡ 1/k B T from its equilibrium value at the reference temperature, δβ(t) ≡ β(t + dt) − β(t), and sharp brackets denote standard canonical averages, the variation of the potential energy is given [37] by Following a small inverse-temperature step of magnitude ∆β, ∆U (t) → − (∆U ) 2 ∆β for t → ∞.Since ∆β = −∆T 0 /k B T 2 0 , this leads to the well-known Einstein expression for the specific heat, c = (∆U ) 2 /k B T 2 0 .Equation (37) implies that the response to a jump at t = 0 for t > 0 is given by ∆U (t) = − (∆U ) 2 + ∆U (0)∆U (t) ∆β (note that right after the temperature step one has ∆U (t) ∼ = 0 because of continuity).Therefore, the linear-response normalized relaxation function, R 0 (t), is given by
The potential-energy time-autocorrelation function appearing in Eq. ( 38) was calculated at the reference temperature T 0 = 0.60 as follows.First, 10 7 time steps were taken for equilibration, which was confirmed from two consecutive runs comparing the self-part of the intermediate scattering function.After that a run of 5 × 10 6 time steps was carried out dumping the potential energy every 32 time steps.The potential-energy time-autocorrelation function was calculated using Fast Fourier Transform as implemented in RUMD [41].
In SPA the constant Λ of Eq. ( 9) is assumed to be proportional to the change in the monitored property, in casu ∆U .Λ was determined using the integral criterion of Ref. [18], which considers two jumps to the same temperature, an up and a down jump.For this we used the jumps from the temperatures 0.55 and from 0.65 to T 0 = 0.60 [21], leading to Here ∆U is the equilibrium potential energy at the starting temperature minus the corresponding quantity at the final temperature (following the tradition in the field).The Λ of Eq. ( 39) was used for all predictions.25) and Eq. ( 29)) defined from the potential energy U after a temperature jump at t = 0 starting from a state of thermal equilibrium (blue filled circles).(a) and (b) show results for magnitude 0.10 temperature up and down jumps to the reference temperature T0 = 0.60; (c) and (d) show results for magnitude 0.05 temperature up and down jumps to T0 = 0.60; (e) and (f) show results for magnitude 0.03 temperature up and down jumps to T0 = 0.60.The orange filled circles are the first-order predictions of the jump differential equation Eq. ( 9) (given in Eq. (25) in which R0(t) is the normalized equilibrium potential-energy time-autocorrelation function at T0 = 0.60, R1(t) is given by Eq. ( 29), and Λ is given by Eq. ( 39)).
For reference, R0(t) in all figures is plotted as small black filled circles.
Temperature-jump simulations were carried out as follows.First, 5 × 10 8 time steps were taken to ensure equilibration at the given starting temperature.A total of 1000 configurations were generated from a subsequent 5 × 10 8 simulation by dumping configurations every 2 19 time step.This ensures that the configurations are statistically in-dependent at the lowest temperature studied (T = 0.50).For each of the 1000 configurations an aging simulation of 10 6 time steps was performed and the potential energy was dumped every eighth time steps.The curves shown in Fig. 1 represent averages over these 1000 aging simulations.
The averages were smoothed using a Gaussian function.Each point represents an average calculated over all data points using R avg (t) = t R(t ) exp(−(t − t ) 2 /σ)/ t exp(−(t − t ) 2 /σ) in which t is the time-step number and σ = 15, 000.In order to reduce the number of points shown in Fig. 1, the data were divided into 24 bins per decade.
Figure 1 shows the simulation results (blue circles) for the normalized relaxation function of the potential energy for up and down jumps to T 0 = 0.60.The orange circles are the predictions of the first-order theory.In all figures the small filled black circles are the normalized equilibrium potential-energy time-autocorrelation function at T 0 = 0.60, which is the linear-limit normalized relaxation function R 0 (t) (Eq.( 38)).This function is faster than R(t) for up jumps and slower for down jumps.This is expected since relaxation is initially slow for the up jumps to T 0 = 0.60 because the fictive temperature [3] in this case is below 0.60, while the opposite happens for the down jumps to T 0 = 0.60.
We see that the theory generally fits data well, even for the fairly large temperature jumps of magnitude 0.05.Deviations between prediction and simulations is observed for larger up jumps, though.A similar pattern has been observed in experiments but there the observed relaxation function is faster than predicted, not slower [21].In both cases these deviations serve to emphasize that SPA is not accurate for large jumps.The TN formalism implies that the long-time decay of the normalized relaxation function for infinitesimal jumps to different temperatures are identical except for an overall scaling of the time, i.e., obeys time-temperature superposition (TTS) [42].In other words, TTS is a necessary condition for TN to apply and thus, in particular, for SPA to apply.We test TTS by plotting the normalized potential-energy time-autocorrelation functions R 0 (t) at temperatures ranging from 0.50 to 0.70 (Fig. 2(a)) and scaling these on the time axis (Fig. 2(b)).Except for the short-time signals that are not relevant to aging, we see that TTS indeed applies to a very good approximation.

VI. SUMMARY AND OUTLOOK
We have solved the jump differential equation analytically to first order.The solution is Eq. ( 25) in which R 1 (t) is given by Eq. ( 29).The solution does not explicitly involve the function F (R); indeed R 1 (t) has a universal expression in terms of the zeroth-order solution, R 0 (t).Since the latter by the fluctuation-dissipation theorem is an equilibrium time-autocorrelation function, our results imply that within the single-parameter aging scheme knowledge of equilibrium fluctuations is enough to predict aging.The expression for R 1 (t) relevant for the weakly nonlinear limit was confirmed by computer simulations of the Kob-Andersen binary Lennard-Jones glass former monitoring the aging of the potential energy following temperature jumps of varying magnitude.
For future developments of the TN single-parameter aging formalism it would be most interesting to monitor the equilibrium fluctuations in experiments in order to check whether aging is predicted correctly from these fluctuations.This is experimentally very challenging, but should not be impossible.It would also be interesting to monitor other quantities in simulations than the potential energy used here, although it should be mentioned that many quantities relax in a very similar way for the Kob-Andersen system [43].APPENDIX Equation ( 13) is derived here from the integral criterion of Ref. [18] that considers two jumps to the same temperature: an up and a down jump.For two small such jumps of same magnitude to the temperature T 0 denoted by a and b, one has the two normalized relaxation functions The integral criterion [18] is Here Λ ab = −Λ ba is the difference in the value of Λ jumping from above and below, implying that Λ ab = 2Λ and Λ ba = −2Λ.When Eq. ( 40) is substituted into Eq.( 41), we get Expanding to second order in Λ leads to This reduces to which implies Eq. ( 36) and thereby Eq. ( 13).

FIG. 1 .
FIG. 1. Results from computer simulations of the Kob-Andersen binary Lennard-Jones model.The figures show the normalized relaxation function R(t) (Eq.(25) and Eq.(29)) defined from the potential energy U after a temperature jump at t = 0 starting from a state of thermal equilibrium (blue filled circles).(a) and (b) show results for magnitude 0.10 temperature up and down jumps to the reference temperature T0 = 0.60; (c) and (d) show results for magnitude 0.05 temperature up and down jumps to T0 = 0.60; (e) and (f) show results for magnitude 0.03 temperature up and down jumps to T0 = 0.60.The orange filled circles are the first-order predictions of the jump differential equation Eq. (9) (given in Eq.(25) in which R0(t) is the normalized equilibrium potential-energy time-autocorrelation function at T0 = 0.60, R1(t) is given by Eq. (29), and Λ is given by Eq. (39)).For reference, R0(t) in all figures is plotted as small black filled circles.

FIG. 2 .
FIG.2.Test of time-temperature superposition for the normalized potential-energy time-autocorrelation function R0(t) at the temperatures indicated in the legends.(a) shows the simulation data and (b) shows the same data empirically scaled on the time axis.We conclude that TTS applies except at the shortest times.