Theory of Specific Heat in Glass Forming Systems
H.G.E. Hentschel, Valery Ilyin, Itamar Procaccia, Nurith Schupper
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug Theory of Specific Heat in Glass Forming Systems
H.G.E. Hentschel ∗ , Valery Ilyin, Itamar Procaccia and Nurith Schupper Department of Chemical Physics, The Weizmann Institute of Science, Rehovot 76100, Israel ∗ Dept of Physics, Emory University, Atlanta Ga 30322. (Dated: November 1, 2018)Experimental measurements of the specific heat in glass-forming systems are obtained from thelinear response to either slow cooling (or heating) or to oscillatory perturbations with a given fre-quency about a constant temperature. The latter method gives rise to a complex specific heat withthe constraint that the zero frequency limit of the real part should be identified with thermodynamicmeasurements. Such measurements reveal anomalies in the temperature dependence of the specificheat, including the so called “specific heat peak” in the vicinity of the glass transition. The aimof this paper is to provide theoretical explanations of these anomalies in general and a quantitativetheory in the case of a simple model of glass-formation. We first present new simulation results forthe specific heat in a classical model of a binary mixture glass-former. We show that in additionto the formerly observed specific heat peak there is a second peak at lower temperatures which wasnot observable in earlier simulations. Second, we present a general relation between the specificheat, a caloric quantity, and the bulk modulus of the material, a mechanical quantity, and thusoffer a smooth connection between the liquid and amorphous solid states. The central result of thispaper is a connection between the micro-melting of clusters in the system and the appearance ofspecific heat peaks; we explain the appearance of two peaks by the micro-melting of two types ofclusters. We relate the two peaks to changes in the bulk and shear moduli. We propose that thephenomenon of glass-formation is accompanied by a fast change in the bulk and the shear moduli,but these fast changes occur in different ranges of the temperature. Lastly, we demonstrate howto construct a theory of the frequency dependent complex specific heat, expected from heteroge-neous clustering in the liquid state of glass formers. A specific example is provided in the contextof our model for the dynamics of glycerol. We show that the frequency dependence is determinedby the same α -relaxation mechanism that operates when measuring the viscosity or the dielectricrelaxation spectrum. The theoretical frequency dependent specific heat agrees well with experimen-tal measurements on glycerol. We conclude the paper by stating that there is nothing universalabout the temperature dependence of the specific heat in glass formers - unfortunately one needs tounderstand each case by itself. I. INTRODUCTION
The traditional measurements of the specific heat C V at constant volume or C P at constant pressure involvecooling (or heating) the sample at a constant rate [1, 2].When applied to glass-forming systems, this approachhas an inherent difficulty. Since glass-forming systemtend to relax to equilibrium slower and slower as the tem-perature is lowered, at some point the ‘constant rate’ ofcooling becomes too high for the system to respond to,and then the system does not reach equilibrium. Typi-cally the specific heat then drops abruptly, giving rise tothe “specific heat peak” at some temperature which issometimes identified as the glass transition temperature T g . Needless to say, such a definition of transition tem-perature is less than compelling, since it clearly dependson the rate of cooling and is not inherent to the systemproperties.In an attempt to overcome this difficulty Birge andNagel [3] introduced ”specific heat spectroscopy”. In thistechnique one keeps the sample close to a temperature T at all times, but perturbs it periodically with a small-amplitude oscillation of frequency ω . Linear responsetheory then relates the amount of heat exchanged at thatfrequency, δQ ( ω ) to the oscillatory temperature pertur- bation δT ( ω ) via the relation δQ ( ω ) = C ( ω ) δT ( ω ) (1)where C ( ω ) is the frequency dependent specific heat thatcan be measured at either constant volume or constantpressure. In order to find the thermodynamic specificheat one needs to extrapolate data to the ω → r − repulsive potential) and one experimentalexample (glycerol). In the context of the first examplewe present results of new simulations that exhibit twodistinct peaks in the curve of the specific heat vs. thetemperature. We present for this example various the-oretical results, culminating with a new scenario to ex-plain the specific heat peaks, i.e. the micro-melting ofclusters. We believe that this is the central point of thepresent paper. To understand the nature of the specificheat anomalies one must understand the physics that isbehind the glassy behavior of this model in general and FIG. 1: (Color online). A snapshot of the system at T = 0 . the existence of the two specific heat peaks in particu-lar. When the temperature is lowered at a fixed pressurethis system [4] (as well as many other glass-formers [5–10]tends to form micro-clusters of local order. In the presentcase large particles form long-lived patches of hexagonalordering first (starting at about T = 0 .
5, and at lowertemperatures (around T = 0 .
1) also the small particlesform long-lived hexagonal clusters. The clusters are notthat huge, with at most O(100) particles, (cf. Fig. 1),depending on the temperature and the aging time. Butwe have shown that the long time properties of correla-tion functions are entirely carried by the micro-clusters[4]. Below we will refer to the micro-clusters as curdsand the liquid phase as whey. We will argue that thespecific heat responds to the micro-melting of the clus-ters - those of small particles at the lowest temperaturesand those of the larger particles at higher temperatures.The large increase in the number of degrees of freedomwhen a particle leaves a crystalline cluster and joins theliquid background is the basic reason for the increase inentropy that is seen as a specific heat peak.In the context of the second example we show thatthe calculation of the frequency dependent specific heatis easy when we have a reasonable model of the glassyrelaxation. Having such a model for glycerol [11], wedemonstrate in section IV that the information gainedfrom the frequency dependent specific heat is very simi-lar to that learned from other linear response functionslike broad-band dielectric spectroscopy. We will be abletherefore to present spectra of the frequency dependentspecific heat in close correspondence with experiments.The specific heat has interesting relations with the me- chanical moduli of the material, and we present relations(which pertain to any system with an r − n potential) tothe bulk and shear moduli. As a result of our think-ing we conclude that the bulk and shear moduli changerapidly in the temperature range of the two distinct spe-cific heat peaks mentioned above. The relation to thebulk modulus is explicit, and is shown rigorously in Sub-sect. II C. The relation to the shear modulus is less ex-plicit, simply because one does not have an equation ofstate with strain (a quantity that is ill defined in thecontext of glasses). Nevertheless we present a conjecturethat the bulk and shear moduli in generic glasses maychange rapidly at two different temperature ranges.The structure of the paper is as follows: In Sect. IIwe present the model glass consisting of a binary mix-ture of point particles interacting via soft potentials, anddiscuss its thermodynamic properties. We derive exactequations for its specific heat at constant volume, whichare correct at all temperatures through the glass transi-tion. We compare these results to molecular dynamicssimulations in which special care had been taken to equi-librate the system, summarizing a computational effortof about two years. The main conclusion of this sectionis that the details of the interaction potential are crucialin determining what the specific heat does in the vicin-ity of the glass transition, and there is nothing universalabout it. For a simple enough potential we can derivea theory that is in excellent agreement with simulationsup to the first specific heat peak. To explain both peakswe must present a theory that takes into account explic-itly the tendency of the system to form micro-clusters[4]. The state of the system then becomes like curds oflocal crystalline order embedded in a whey of disorderedfluid. It is the freezing or melting of these curds thataccount quantitatively for the specific heat peaks, as isshown in Subsect. III. Below we use interchangeably thewords ‘clusters’ and ‘curds’. In Sect. IV we turn to dis-cussing the frequency-dependent complex specific heat.To construct a theory of this object one needs a modelof the dynamics of the system under study, be it glyc-erol or any other material. We demonstrate, using ourdynamical model of glycerol [11], how this measurementis equivalent in terms of its dynamical contents to anyother linear response to an oscillatory perturbation. Wepresent theoretical spectra and show satisfactory agree-ment with the experiments. The paper ends in Sect. Vwhere we draw conclusions and summarize the resultsand the implications of our calculations. II. THE BINARY MODEL AND ITS SPECIFICHEAT
The model discussed here is the classical example [12,14] of a glass-forming binary mixture of N particles ina 2-dimensional domain of area V , interacting via a soft1 /r repulsion with a ‘diameter’ ratio of 1.4. We referthe reader to the extensive work done on this system[12, 14–17]. The sum-up of this work is that the modelis a bona fide glass-forming liquid meeting all the criteriaof a glass transition.In short, the system consists of an equimolar mixture oftwo types of particles, “large” with ‘diameter’ σ = 1 . σ = 1, respectively, butwith the same mass m . In general, the three pairwiseadditive interactions are given by the purely repulsivesoft-core potentials φ ab ( r ) = ǫ (cid:16) σ ab r (cid:17) n , a, b = 1 , , (2)where σ aa = σ a and σ ab = ( σ a + σ b ) /
2. The cutoff radiiof the interaction are set at 4 . σ ab . The units of mass,length, time and temperature are m , σ , τ = σ p m/ǫ and ǫ/k B , respectively, with k B being Boltzmann’s con-stant. In numerical calculations the stiffness parameterof the potential (2) was chosen to be n = 12.We turn now to the analysis of the specific heat of thismodel as a function of the temperature. A. Specific heat (simulations)
The specific heat capacity (specific heat) at constantvolume is defined by: C V N = d ∂∂T h U i N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) V , (3)where d is the space-dimension and the potential energyof a binary mixture is given by: U = 12 X i = j φ ab ( r ij ) . (4)Here r ij is the distance between particles i and j . Theaverage value of the potential energy is defined by aver-aging over configurational space Γ : h U i = R U exp( − UT )dΓ R exp( − UT )dΓ . (5)Substitution of (5) into (3) yields the following expressionfor the specific heat: C V N = d h U i − h U i N T . (6)The specific heat of our binary mixture model was mea-sured at constant volume in [13, 14, 18] and by us. Insimulations one can measure the specific heat directlyfrom its definition (3) or (6). We have used the last equa-tion which allows one to estimate the specific heat in asingle run of the canonical ensemble Monte Carlo simu-lations. At each temperature the density was chosen inaccordance with the simulation results in an NPT ensem-ble as described in [14] with the pressure value fixed at U / N P ( U / N ) T=0.4T=0.6T=1T=3
FIG. 2: Color online: the energy distribution functions for thebinary mixture model, computed at constant volume, suchthat the volume agrees with the pressure P=13.5 [14] at eachtemperature. P = 13 .
5. As the initial configuration in the Monte Carloprocess the last configuration of the molecular dynamicsrun for this model at given temperature after 1 . × time steps was used. After short equilibration the po-tential energy distribution functions were measured dur-ing 2 × Monte Carlo sweeps. The acceptance ratewas chosen to be 30%. Simulations were performed with N = 1024 particles in a square cell with periodic bound-ery conditions.Examples of the spline interpolation of the potentialenergy distribution for a few temperatures are shownin Fig.2. The first and second moments of these distri-butions define the average value of the potential energy(Fig.3) and the specific heat (Fig4). We stress that theseresults were computed at constant volume, such that thevolume corresponds to simulations in NPT ensemble withthe pressure P=13.5 [14] at each temperature.One can see from these figures that the behavior ofboth quantities, the first and second moments of the dis-tribution, change abruptly in the vicinity of T ∼ .
5. Thespecific heat displays a maximum in the temperature de-pendence. Our simulations appear to provide trustablevalues of C V down to lowest temperatures where thevalue of the specific heat coincides with that of two-dimensional solid, i.e. C V = 2. What could not be seenin earlier simulations is that there is a much smaller sec-ond peak of the specific heat at lower temperatures. Toresolve it to the naked eye we present in Fig. 4 a blow-up of the region of lowest temperatures where the secondpeak is more obvious. To understand the nature of thespecific heat anomalies we turn now to a theoretical anal-ysis of the caloric equation of state in order to study thespecific heat using the definition (3). The physical originof the two peaks will be explained in Subsect. III. The T < U > / N Virial expansionMC−NVT
FIG. 3: Color online: in dots: the temperature dependence ofthe average potential energy per in the binary mixture model,computed at constant volume, such that the volume agreeswith the pressure P=13.5 [14] at each temperature. The con-tinuous line represents the approximation furnished by thevirial expansion, which obviously fails for
T < . . T C V / N T C V / N FIG. 4: Color online: in dots: the temperature dependenceof the specific heat in the binary mixture model, computed atconstant volume, such that the volume agrees with the pres-sure P=13.5 [14] at each temperature. The data indicate theexistence of two specific heat peaks, one prominent at about T = 0 . T = 0 .
1, and see the insetfor finer detail. The continuous line represents the approx-imation furnished by the virial expansion, which obviouslyfails for
T < . . reader who is mostly interested in the physical insight isinvited to jump to that subsection. B. Specific heat (series expansions)
The general expression for the pressure (thermal equa-tion of state) obtained from the virial theorem is givenby: P = ρT − ρ d · N h X i = j r ij ∂φ ab ( r ij ) ∂r ij i , (7)where ρ is the particle number density. The potential isa homogenous function of degree − n (Eq. (2)), therefore: r ∂φ ab ( r ) ∂r = − nφ ab ( r ) . (8)Due to this property of the interaction potential we finda connection between the pressure and the temperatureand mean energy: P = ρT + nd ρ h U i N . (9)This equation is exact for one component and multicom-ponent systems and is valid at all temperatures, fromliquid to solid.The next simplification for systems with an inversepower inter-particle interaction consists in the depen-dence of all excess thermodynamic properties relative tothe ideal gas on a single density-temperature variable[19, 20]. To see why, recall that for a one componentsystem the canonical partition function Z N is defined by: Z N = 1 N !Λ dN Z exp − σ n ǫk B T X r nij ! d ~r N . (10)Here Λ = h/ (2 πmk B T ) / is the thermal de Brogliewavelength. The typical distance between particles isgiven by l = VN ! /d . (11)Thus one can use new variables in the integral (10), ~s = ~r/l , and the canonical partition function can be rewrittenas: Z N = V N N !Λ dN N N Z exp − ρ n/d ǫk B T X s nij ! d ~s N , (12)where the dimensionless particle number density is de-fined by ρ = NV σ d . This way of writing the partitionfunction underlines the existence of the ideal gas contri-bution before the configurational integral, and the de-pendence of the configurational integral, in the case ofone component, on a natural parameter, Γ = ρ ( ǫk B T ) d/n .In the case of a multicomponent system the propertiesof the mixture can be approximated by those of a onecomponent reference fluid [21] with an effective diameterdefined by: σ de = X a,b x a x b σ dab , (13)where x a = N a /N is the particle number concentration.Therefore, the properties of a mixture are defined by theeffective parameter: Γ e = Γ σ e σ ! d . (14) Nevertheless, in [14] it was shown, that for soft poten-tial a more suitable definition of the effective diameter isgiven by: σ e = x σ + x σ . (15)Such a definition leads to a more accurate virial expan-sion, as obtained for the present model by [14] usingmolecular dynamics simulations in the temperature range0 . ≤ T ≤ PρT = 1 + 1 . e + 2 . e + 2 . e + 7 . e − . e + 27 . e − . e + 5 . e . (16)Substitution of the equation (16) to (9) yields the caloricequation of state, the specific heat is calculated after thatby (3). The density corresponding to a given temperatureis defined as the solution of the equation (16) at P = 13 . T = 0 . C. Specific heat (mechanical approach)
In this subsection we connect the specific heat to thebulk modulus of the system. To this aim we begin withthe microscopic definition of the stress-tensor in an NVTensemble (see, e.g. [22]): σ αβ = 1 V X i p αi p βi − X j = i ∂φ ab ( r ij ) ∂r ij r αij r βij r ij , (17)where p αi is the α component of the dimensionless mo-mentum of particle i and r αij is the α component of thevector joining particles i and j . The first invariant of the stress tensor is its trace: V ( σ xx + σ yy ) = ( X i ( p xi ) + X i ( p yi ) ) + nU. (18)In order to average (18) one has to take into account that: h X i ( p xi ) i = N h ( p xi ) i = N · T. (19)Thus the average of the first invariant (18) is:( h σ xx i + h σ yy i ) = 2 ρT + nρ h U i N . (20)The pressure is defined as P = ( h σ xx i + h σ yy i ) / V ( σ xx σ xx + σ yy σ yy + 2 σ xx σ yy ) (21)= X i ( p xi ) + X i ( p yi ) ! +2 n X i ( p xi ) + X i ( p yi ) ! U + n U . To compute the average of this equation we need to usethe fact that h ( p xi ) i = T and h ( p xi ) i = 3 T . Afteraveraging (21) is written as: V ( h σ xx σ xx i + h σ yy σ yy i + 2 h σ xx σ yy i )= 4 N T + 4 N T + 4 N T n h U i + n h U i . (22)The square of (20) is given by: V ( h σ xx i + h σ yy i + 2 h σ xx ih σ yy i )= 4 N T + 4 N T n h U i + n h U i (23)After substraction of (23) from (22) we have: V " ( h σ xx σ xx i − h σ xx i ) + ( h σ yy σ yy i − h σ yy i )+ 2( h σ xx σ yy i − h σ xx ih σ yy i ) = 4 N T + n ( h U i − h U i ) . (24)This equation has a well defined thermodynamic limitsince the quantity in square brackets scales like V − .This is seen explicitly in Eq. (25).These results allow us now to find exact relationshipsbetween the specific heat, a caloric quantity, and the elas-tic constants which are mechanical quantities. To do sowe recall the definitions of the elastic constants throughthe stress fluctuations (see, e.g., [23]): VT ( h σ αβ σ γδ i − h σ αβ ih σ γδ i ) = 2 ρT ( δ αγ δ βδ + δ αδ δ βγ )+ ( C Bαβγδ − C αβγδ ) , (25)where C αβγδ are the elastic constants and C Bαβγδ are theso called Born terms which determine the instantaneouselastic constants for any given configuration.Substitution of (25) to (24) yields: n ( h U i − h U i ) = 4 N T + V " ( C Bxxxx + C Byyyy + 2 C Bxxyy ) − ( C xxxx + C yyyy + 2 C xxyy ) (26)The compression (bulk) modulus in two-dimensionalsystems is: K = 14 ( C xxxx + C yyyy + 2 C xxyy ) . (27)Recalling Eq. (6) and substituting it and Eq. (27) into(26) yields: C V N = 1 + 4 K ∞ − Kn ρT ., (28)where the bulk modulus in the infinite frequency limit K ∞ = ρT + K B and the Born term is defined by theinterpaticle interactions [24]: K B = 14 V X i = j r ij r ij ∂ φ ab ( r ij ) ∂r ij − ∂φ ab ( r ij ) ∂r ij ! . (29)For the present model the “infinite frequency” term (cf.[24]) is given by: K ∞ = ρT + n ( n + 2)4 ρ h U i N . (30) This is as much as one can go using exact identities.We reiterate that Eq. (28) is very interesting, allowing usto connect the bulk modulus to the specific heat. In fact,this connection implies that the specific heat measuresthe difference between the bulk modulus and its infinitefrequency limit. At low temperatures this difference inthe harmonic approximation as given by: K ∞ − K = n ρT , (31)independent of the solid structure in contrast to the shearmodulus (cf. [17]).The bulk modulus K cannot be computed exactly us-ing identities, and we need further information to eval-uate it. Fortunately we can estimate the bulk modulusfrom the virial expansion (7) at T > .
5, since : K = ρ ∂P∂ρ . (32)Having done so we can compare the measurements towhat we expect theoretically. The specific heat as pre-dicted from the virial expansion is shown in Fig. 4 asthe blue (continuous) line. We should stress that com-puting C V from either Eq. (3) or Eq. (28) (using thevirial expansion (16), yield essentially identical resultsthat cannot be distinguished in the blue line in Fig. 4for T > . III. SPECIFIC HEAT - THE PHYSICALEXPLANATION
In this section we propose the physical picture behindthe existence of the specific heat peaks. We argue thatthe specific heat responds to the micro-melting of theclusters - those of small particles at the lowest tempera-tures and those of the larger particles at higher temper-atures. The large increase in the number of degrees offreedom when a particle leaves a crystalline cluster andjoins the liquid background is the basic reason for the in-crease in entropy that is seen as a specific heat peak. Wecan specialize these observations for the model at hand(with inverse power potential) or present the discussion ingreater generality for any model. These two approachesare presented in the two following subsections.
A. Mechanical Equation of State
In this subsection we employ the mechanical equationof state derived above from which the specific heat willbe computed. To start we define v ℓw , v sw , v ℓc and v sc re-spectively as the volume of large particle in the whey,small particle in the whey, large particle in the solid andsmall particle in the solid. Similarly we denote by ǫ ℓw , ǫ sw , ǫ ℓc and ǫ sc the energy of a large and small particle inthe in the whey and in the crystalline phase respectively.Needless to say, all these quantities are temperature andpressure dependent; we will therefore explicitly use ourlow temperature knowledge concerning v ℓc and v sc in thecrystalline phase, but treat the difference v ℓw − v ℓc and v sw − v sc as constants that we estimate below from oursimulation knowledge. Similarly we estimate ǫ ℓc and ǫ sc from our knowledge of the hexagonal lattices at T = 0.We assume that ǫ ℓw ≈ ǫ ℓc and similarly ǫ sw ≈ ǫ sc sinceour simulations indicate a very small change in these pa-rameters, see Table I. It should be stressed that the en-thalpy change at these pressures are almost all due to the P V term. This will result in a semi-quantitative theoryascribing the important changes in specific heat to thechanges in the fraction of particles in curds and whey. Inother words the number of particles in the whey and thenumber of clusters are all explicit functions of tempera-ture and pressure. ǫ ℓc ǫ sc ǫ ℓw ǫ sw v ℓc v sc v ℓw v sw As the condensed phase consists of clusters of large andsmall partilces, we use the notation N ℓn for the numberof clusters of n large particles and N sm for the clusters of m small particles. In the next subsection we write theenergy of our system explicitly in terms of these clusternumbers. Here however we only need the intensive vari-ables p ℓc = 2 P n N ℓn /N , p sc = 2 P m N ℓm /N p ℓw = 2 N ℓw /N and p sw = 2 N sw /N which stand for the fraction of largeparticles and small particles in the curds, and large parti-cles and small particles in the whey, such that p ℓc + p ℓw = 1and p sc + p sw = 1. Using these variables we can write anexpression for the volume per particle v ≡ V /N : v = v ℓw + v sw v ℓc − v ℓw p ℓc + v sc − v sw p sc . (33)At this point we need to derive expressions for p ℓc and p sc . To do so we need to remember that in the rel-evant range of temperatures the large particles in thewhey can occupy either hexagonal or heptagonal Voronoicells, whereas small particles can occupy only pentago-nal or hexagonal cells [4, 15, 16]. Accordingly there are g ℓw ≈ (2 − / / g sw ≈ (2 − / / g ℓw and g sw . This is not our aim here; we aim at aphysical understanding of the specific heat peaks ratherthan an accurate theory. We thus end up with the simple estimates p ℓc ( P, T ) ≈
11 + g ℓw e [( ǫ ℓc − ǫ ℓc )+ P ( v ℓc − v ℓw )] /T , (34) p sc ( P, T ) ≈
11 + g sw e [( ǫ sc − ǫ sw )+ P ( v sc − v sw )] /T . (35)It is important to note that the combination of Eq.(33) together with Eqs. (34) and (35) provides a me-chanical equation of state that is alternative to the virialexpansion presented above. Whereas the latter is best attemperature higher than T ≈ . C v directly fromEq. (28). The peaks in the specific heat are determinedby the temperature dependence of p ℓc ( P, T ) and p sc ( P, T )each which has a temperature and pressure derivativesthat peaks at a different temperature, denoted as T ℓ ( P )and T s ( P ), and see below for details. As said above wetake ∆ v ℓ ≡ v ℓw − v ℓc and ∆ v s ≡ v sw − v sc as approximatelyconstants (as a function of temperature and pressure).The constants are estimated from the condition that thesecond temperature derivative of p ℓc ( P, T ) and p sc ( P, T )should vanish. Explicitly:∆ v ℓ ≈ T ℓ ( P ∗ ) ln g ℓw /P ∗ , ∆ v s ≈ T s ( P ∗ ) ln g sw /P ∗ , (36)where P ∗ is the pressure for which the peaks in thederivatives are observed (13.5 in our simulations). Thisis equivalent to a linear dependence of the specific heatpeaks as a function of pressure, T ℓ ( P ) /T ℓ ( P ∗ ) = P/P ∗ and similarly for the small particles.In terms of these objects we can rewrite v = v c ( P, T ) + ∆ v ℓ (1 − p ℓc ) + ∆ v s (1 − p sc ) , (37) (cid:18) ∂v∂P (cid:19) T = (cid:18) ∂v c ∂P (cid:19) T − ∆ v ℓ (cid:18) ∂p ℓc ∂P (cid:19) T − ∆ v s (cid:18) ∂p ℓc ∂P (cid:19) T . (38)To compute the temperature dependence of (cid:0) ∂v∂P (cid:1) T weneed first to determine its T → T →
0. Since we have already exact results for the bulkmodulus for the present model, we return to Eqs. (28)and (30). We know on the one hand that lim T → C v = 2and that h U i /N ≈ .
94 over the whole interesting tem-prature range, cf. Fig. 3. The compressibility κ is re-lated to the bulk modulus via κ = − (cid:0) ∂v∂P (cid:1) T /v = 1 /K and therefore easily estimated as T → (cid:0) ∂v c ∂P (cid:1) T ≈ − / (123 . − T ). We use this approximationup to T ≈ . C v /N . Theparameters used were estimated from the numerical sim-ulation and are summarized in Table I. Since the aim ofthis subsection is only semi-quantitative, we do not makeany attempt of parameter fitting, and show the result ofthe calculation in Fig. 5. C v / N FIG. 5: Specific heat at constant volume as predicted by thesimple theory which is based on the mechanical equation ofstate supplied by Eqs. (55),(33) and (34) and (35). Note thatthe theory predicts the two peaks which are associated withthe micro-melting or micro-freezing of the clusters of largeand small particles respectively. The magnitude of the peaksis too high, reflecting terms missing in the simple approach,like the effect of anharmonicity at the lowest temperatureswhich are negative, tending to decrease the height of the low-temperature peak.
Indeed, the theoretical calculation exhibits the exis-tence of two, rather than one, specific heat peaks. Wecan now explain the origin of the peaks as resulting fromthe derivatives (cid:16) ∂p sc ∂P (cid:17) T and (cid:16) ∂p ℓc ∂P (cid:17) T . These derivativeschange most abruptly when the micro-clusters form (ordissolve), each at a specific temperature determined by( h sw − h sc ) / ln g sw and ( h ℓw − h ℓc ) / ln g ℓw . Note that there canbe pressures (both upper and lower boundaries) wherethe the sign of ( h sw − h sc ) or ( h ℓw − h ℓc ) change sign andthe peak can be lost. B. Caloric Equation of State
In this subsection we present a more general approachwhich does not take direct input from results derived forthe inverse power potential. Thus although we use belowsome parameters read from the simulation, the derivationis very general any pertains to any distribution of clus-ters. To this aim we derive a second equation of state,a caloric one. It is quite standard to have two equationsof state, only for ideal gas and inverse potentials the twoequations of state degenerate into one.Denote the total energy of the system as a sum of E c ,the energy of the clusters (or curds), and E w , the energyof the liquid background (or whey), i.e : E = E c + E w . (39)These energies are sums over the degrees of freedom - translational, rotational, vibrational and configurational: E c = E tr,c + E rot,c + E vib,c + E conf,c , (40) E w = E tr,w + E vib,w + E conf,w . (41)To estimate E tr,c we consider the number N ℓn of clustersof n large particles and N sm of clusters of m small particlesand write E tr,c = T [ X n N ℓn + X m N sm ] . (42)On the other hand in the whey we follow Eyring [25] andGranato [26] and write E tr,w = T [ N ℓw + N sw ] f . (43)where f ≡ − V ( s ) w /V w is defined as the fraction of freevolume in the liquid phase compared to the equivalentsolid crystalline phase. In other words, V w = N ℓw v ℓw + N sw v sw . (44) V ( s ) w = N ℓw v ℓc + N sw v sc . (45)Similarly, we write E rot,c = 12 T [ X n N ℓn + X m N sm ] , (46) E vib,c = 2 T [ X n N ℓn n + X m N sm m ] , (47) E vib,w = 2 T [ N ℓn + N sm ](1 − f ) , (48) E conf,c = ǫ ℓc X n N ℓn n + ǫ sc X m N sm m , (49) E conf,w = ǫ ℓw N ℓw + ǫ sw N sw . (50)In terms of these variable we can rewrite E tr = N T [ p ℓc h n i + p sc h m i + ( p ℓw + p sw ) f ] , (51) E rot,c = N T p ℓc h n i + p sc h m i ] , (52) E vib = N T [2(1 − f ) + ( p ℓc + p sc ) f ] , (53) E conf = N ǫ ℓw + ǫ sw + ( ǫ ℓc − ǫ ℓw ) p ℓc + ( ǫ sc − ǫ sw ) p sc ] . (54)Summing up all the contributions we need to pay at-tention to the order of magnitude of the various terms.Since we expect the average size of clusters, at the tem-peratures of interest, to be of the order of 30 or so, we canneglect safely all the terms that have average size clustersin the denominator. With this in mind the expression forthe energy of the system takes the form EN = T [2 − − p ℓc − p sc f ]+ ǫ ℓw + ǫ sw ǫ ℓc − ǫ ℓw p ℓc + ǫ sc − ǫ sw p sc . (55)This is the caloric equation of state that we were after.Using it, we can compute C v directly from the thermo-dynamic identity C v N = C p N − T vα κ , (56)where the thermal expansion coefficient is α ≡ (cid:18) v ∂v∂T (cid:19) P = (57)= 1 v [ ( v ℓc − v ℓw )2 (cid:18) ∂p ℓc ∂T (cid:19) P + ( v sc − v sw )2 (cid:18) ∂p sc ∂T (cid:19) P ] , and the compressibility is κ ≡ − (cid:18) v ∂v∂P (cid:19) T (58)= − v [ ( v ℓc − v ℓw )2 (cid:18) ∂p ℓc ∂P (cid:19) T + ( v sc − v sw )2 (cid:18) ∂p sc ∂P (cid:19) T ] . The last object that we need to obtain for evaluating C v is C p : C p = (cid:18) ∂E∂T (cid:19) P + P (cid:18) ∂V∂T (cid:19) P . (59)Using our expressions (55) and (33) we find C p N = (cid:20) − (2 − p ℓc − p sc f (cid:21) + T (cid:18) ∂p ℓc ∂T (cid:19) P (cid:20) f h ℓc − h ℓw )2 T (cid:21) + T (cid:18) ∂p sc ∂T (cid:19) P (cid:20) f h sc − h sw )2 T (cid:21) − (2 − p ℓc − p sc )2 T (cid:18) ∂f∂T (cid:19) P . (60)Having all the ingredients we can sum up the terms inEq. (56). The parameters used were estimated fromthe numerical simulation and are summarized in Table I.As before, we did not make any attempt for parameterfitting, and show the result of the calculation in Fig. 6. C. Discussion
The bottom line of the simple theory described inthe previous subsection is that there are two importantranges of temperature, first around T ≈ . T ≈ . K ∞ counterparts.We cannot at this point asses how general is this split between bulk and shear moduli, and whether it will beseen in generic glasses. We thus leave this point for fur-ther research, stressing that we expect this phenomenonto appear whenever there exist micro-clusters of preferredoredering in the scenario of glass-forming.Having an effective equation of state, albeit approxi-mate, we can easily compute any thermodynamic deriva-tive of interest. We wrote above explicit expressions forthe compressibility and the thermal expansion coefficient.Others are as easily calculated. The point to stress how-ever is how non-universal the thermodynamics is. In thismodel we have two specific heat peaks, in others we mighthave one or several. We would also expect a strong pres-sure dependence for these peaks. IV. FREQUENCY-DEPENDENT SPECIFICHEAT
By applying time dependent heat fluxes δQ ( t ) to theliquid and measuring the resulting temperature fluctua-tions δT ( t ), the specific heat can be measured δQ ( t ) = CδT ( t ). As mentioned in the introduction, measure-ments of the specific heat of glassy fluids at low tempera-tures can in principle be made under conditions of eitherconstant volume (isochoric conditions) or constant pres-sure (isobaric conditions), but experimentally isobaricconditions are the norm.The first and best known measurement of the fre-0 C / N C p C v FIG. 6: Specific heat at constant volume and at constantpressure as predicted by the simple theory which is based onthe caloric equation of state supplied by Eqs. (56)-(59).Inagreeement with the simulations and the theory based on themechanical equation of state (Cf. Fig. 5) the theory predictsthe two peaks to both C v and C p which are associated withthe micro-melting or micro-freezing of the clusters of largeand small particles respectively. Note that C p ≥ C v as isexpected from thermodynamics. quency dependent complex specific heat was performedin glycerol, and we take these experimental results asour motivation for this section. We stress from the be-ginning that our approach is not particular for glycerol,and it can be applied to any other material where, as weassume for glycerol, there exists clusters of various sizesthat determine the dynamical response. In order to de-velop a model of the frequency dependent specific heat inglycerol we will employ our own model of the glassy phaseof glycerol. This model assumes that glassy glycerol isa heterogeneous fluid on macroscopic timescales. Thatis, that while on very long timescales the liquid phase ishomogeneous, there exist localized mesoscale domains inthe fluid that have macroscopic lifetimes. Indeed, inho-mogeneities that appear to survive for 10 seconds con-tribute to the dielectric response in the Fourier domain atfrequencies as low as 10 − Hz in some cases. Clearly suchimhomogeneities will also contribute to anomalies in thefrequency dependent specific heat C p ( T, ω ). We developthe theory by deriving expressions for the time-dependententhalpy fluctuations h ∆ H ( t )∆ H (0) i that are related tothe frequency dependent specific heat at constant pres-sure in terms of the distribution of these heterogeneities.The reader is referred to [11] for an introduction to thedynamical model of glassy glycerol in which the dielectricspectra are computed in great detail. A. Frequency Dependent Specific Heat
By considering a temperature field T ( t ) = T + δT ( t ) , t < T ( t ) = T, t > − βH ) /Z , with the en-thalpy given by H = E + P V , Nielsen and Dyre [27]find that the frequency dependent specific heat is givenby the form C p ( T, ω ) = h ∆ H i k B T + iωk B T Z ∞ h ∆ H ( t )∆ H (0) i e + iωt dt. (61)In Eq. (61) ∆ H ( t ) = H ( t ) − h H i is an enthalpy fluctua-tion away from equilibrium. Therefore we write H ( t ) = X s N s ( t ) H s + M l ( t ) h l (62)where N s ( t ) is the number of clusters consisting of s molecules in the glassy phase and M l ( t ) are the remainingmolecules in the mobile liquid phase. H s is the enthalpyof a cluster of s molecules at a pressure P and tempera-ture TH s ( P, T ) = E s ( P, T ) +
P V s ( P, T ) = ( ǫ c + pv c ) s + σs / . (63)In Eq. (63) ǫ c ( P, T ) is the energy/molecule in the con-densed phase; v c ( P, T ) is the volume per molecule in thecondensed phase; and σ ( P, T ) is the surface energy permolecule. Finally h l ( P, T ) = ǫ l + P v l ( P, T ) is the en-thalpy per molecule in the mobile phase.Now in equilibrium we can write h H i = X s h N s i H s + h M l i h l (64)and we also have the sum rule X s sN s ( t ) + M l ( t ) = M (65)where M is the total number of molecule in the system.We can write the enthalpy fluctuations away from equi-librium at time t as∆ H ( t ) = H ( t ) −h H i = X s ( N s ( t ) −h N s i ) h s = X s ∆ N s h s (66)where h s = H s − sh l = ( ǫ c − ǫ l ) s + P ( v c − v l ) s + σs / . (67)Let us first calculate the equilibrium fluctuations h (∆ H ) i . To this end we assume that there are no cor-relations between the dynamics of clusters of differentsizes, implying h ∆ N s ∆ N ′ s i = 0. Then, using the ex-pression Eq. (66) for the enthalpy fluctuations, we canimmediately write that h (∆ H ) i = X s h ∆ N s i h s . (68)Similarly for the time dependent enthalpy fluctuations h ∆ H ( t )∆ H (0) i = X s h ∆ N s ( t )∆ N s (0) i h s . (69)1We have thus reduced the correlation functions for theenthalpy fluctuations into expressions involving the fluc-tuations in cluster number for clusters of different sizes s . To estimate these correlation functions we proceed asfollows. First we note that the number of molecules N c ( t )in the clusters can be written as a sum over the clustersas N c ( t ) − P s N s ( t ) s , and consequently the fluctuationsin the total number of particles away from equilibriumare ∆ N c ( t ) = P s ∆ N s ( t ) s . Then, assuming Gaussianfluctuations we estimate the mean square fluctuations ofthe number of particles within some small volume h (∆ N c ) i ≈ h N c i (70)or rewriting X s h (∆ N s ) i s = X s h N s i s . (71)From this equation we therefore see that h (∆ N s ) i = h N s i /s . (72)For the time dependent fluctuations therefor, assumingan independent Debye relaxation for each cluster, h ∆ N s ( t )∆ N s (0) i = h N s i e − t/τ s /s (73)where τ s is the lifetime of a cluster of size s . Thus weget our final expression for the enthalpy fluctuations inequilibrium in terms of the cluster size distribution as h (∆ H ) i = X s h N s i ( h s /s ) . (74)and for their time dependent correlations h ∆ H ( t )∆ H (0) i = X s h N s i e − t/τ s ( h s /s ) . (75)We now substitute these expressions into Eq. 61 withthe result C p ( T, ω ) = 1 k B T X s h N s i ( h s /s )1 − iωτ s (76)or splitting the specific heat into its real and imaginaryparts ℜ C p ( T, ω ) = 1 k B T X s h N s i ( h s /s )1 + ( ωτ s ) ℑ C p ( T, ω ) = 1 k B T X s h N s i ( h s /s )( ωτ s )1 + ( ωτ s ) (77)We can now use our previous results for the clusterdistributions in the case of glycerol to find the real andimaginary specific heat anomalies in the case of glycerol.We do not re-fit any of the parameters used in the calcula-tion of the BDS spectra [ ? ], we simply use the previous −4 −3 −2 −1 0 1 2 3 401234567 x 10 −3 log ( ω ) ℜ C p κ Real part of the specific heat times thermal conductivity
T=210KT=205KT=200KT=195K
FIG. 7: The theoretical real part of the specific heat ℜ C p ( ω )multiplied by the thermal conductivity for glycerol. knowledge at the temperatures indicated, and plot theresults, fitting only the heat conductivity of glycerol. Weapproximate h s /s ≈ ( h c − h ℓ ) s , such that Eq (76) isrewritten as C p ( T, ω ) ≈ ( h c − h ℓ ) k B T X s h N s i s − iωτ s . (78)Splitting the specific heat into its real and imaginaryparts ℜ C p ( T, ω ) ≈ ( h c − h ℓ ) k B T X s h N s i s ωτ s ) , (79) ℑ C p ( T, ω ) ≈ ( h c − h ℓ ) k B T X s h N s i s ( ωτ s )1 + ( ωτ s ) . (80)The resulting curves multiplied by the thermal conduc-tivity are shown in Figs. 7 and 8. These should be com-pared to Fig.2 of [3]. The reader can convince himselfthat the theory captures the experimental results quan-titatively.We would like to stress at this point that the resultsobtained here are equivalent in dynamical contents tothe computation of the dielectric spectra in [11]. In thatcalculation one focused on the dielectric response ǫ ( ω )and as here decomposed it into its real and imaginaryparts ǫ ( ω ) = ℜ ǫ ( ω ) + i ℑ ǫ ( ω ). It was found in [11] thatwithout the dc contribution (which is absent in the caseof specific heat) we could write ℜ ǫ ( ω ) − ǫ ∞ ( ǫ − ǫ ∞ ) = P s h N s i s/ [1 + ( ωτ s ) ] P s h N s i s (81) ℑ ǫ ( ω )( ǫ − ǫ ∞ ) = P s h N s i s ( ωτ s ) / [1 + ( ωτ s ) ] P s h N s i s . Once normalized the specific heat spectra are identicalthese spectra . The reason for the identity is in the as-sumptions that h m s · m s i ∼ s and h h s i ∼ s , cf. Eq.(67). On the other hand the role of the relaxation times τ s and the distribution of cluster sizes are exactly thesame in the two expressions.2 −4 −3 −2 −1 0 1 2 3 400.20.40.60.811.21.41.61.82 x 10 −3 log ( ω ) ℑ C p κ Imaginary part of the specific heat times thermal conductivity
T=210KT=205KT=200KT=195K
FIG. 8: The theoretical imaginary part of the specific heat ℑ C p ( ω ) multiplied by the thermal conductivity for glycerol. V. SUMMARY AND DISCUSSION
Probably the most glaring consequence of the calcu-lations presented in this paper is that the specific heatis a valuable indicator of the interesting physics that oc-curs during the glass transition, but this transition is inno way universal. The temperature dependence of thespecific heat is determined by details like inter-particlepotentials and micro-melting or micro-formation of clus-ters. In this sense any hope for universality is untenable.Nevertheless we have shown that the specific heat peaksherald interesting new physics, leading to fast changes inthe mechanical moduli which are also associated with fastchanges in the inhomogeneities that are crucial for theglassy behavior, i.e. the formation of micro-clusters. Wepropose that the appearance of two specific heat peaksin the case of the binary mixture indicates two differentranges for the increase in moduli, the bulk modulus athigher temperatures when the first type of clusters form,and the shear modulus when the other type of clustersform, and the ’lubricating’ effect that allows the systemto shear disappears. All this interesting physics is in-dicated by the behavior of the thermodynamics specificheat. As for the complex specific heat we have shown, inthe context of the example of glycerol, that the physicsrevealed by the complex specific heat compared to othermethods of linear response like Broad Dielectric Spec-troscopy are identical. In fact, a straightforward conse-quence of our model for glycerol is the prediction thatthe spectra measured from specific heat can be dividedby the spectra computed, say, from BDS and the resultshould be a constant number. We do not have data forexactly the same temperature, but such an experimentwould be very useful for the near future.It is interesting to see in future research whether thetwo specific heat peaks discussed above may be seen inother systems, or may be an even richer scenario canappear, with more peaks, when more types of clustersintervene in the process of glass-formation.
Acknowledgments
We are grateful to Sid Nagel and Sasha Voronel fordiscussion, both in person and electronic, regarding theirmeasurements of specific heat. This work had been sup-ported in part by the German Israeli Foundation, theMinerva Foundation, Munich, Germany and the IsraelScience Foundation.
APPENDIX A: THE ZERO TEMPERATURELIMIT t follows from the simulation results that at lower tem-peratures the specific heat is at least close to the valueof the solid. Therefore, we can write the energy of thesystem in the harmonic approximation: U = U + dN X i,j a ij q i q j , (A1)where U is the potential energy of the system in thereference state, a ij are expansion coefficients and q i is aCartesian coordinate of the deviation of the current posi-tion of a particle from equilibrium. There are dN degreesof freedom (neglecting translations and rotations of thesystem) in (A1), therefore it follows from the equipar-tition theorem that the average potential energy of thesystem in the solid state is [28]: h U i N = U N + d T. (A2)Substituting(A2) in (3) immediately yields C V = 2 fortwo dimensional solids. At the equilibrium configurationthe potential energy (4) of the binary mixture model (2)is given by: U = 12 ǫ X i = j (cid:18) σ ab r ij (cid:19) n (A3)Due to the scaling properties of the inverse power poten-tial it is possible to normalize the interparticle distancesby the typical distance (11) r ij → s ij = r ij /l [29] : U N = c M ρ nd , (A4)where the constant c M = (1 / N ) P ij /s nij is indepen-dent of the density. Note that this constant is known asthe Madelung constant in solid-state physics. Taking intoaccount (A4) the average potential energy of a harmonicsolid (A2) can be rewritten as: h U i N = c M ρ nd + d T. (A5)The value of the constant c M can be calculated simula-tionally using Eq. (A5); the results are shown in Fig.9.3 T c M FIG. 9: Evaluated values of the constant c M from the simu-lation results. One can see that below T = 0 . c M = 14 . PρT = (1 + n nd ( c M /σ e )Γ n/de . (A6)The value of the renormalized constant c M /σ e = 1 . .
268 [31]. The fact that this constant is expected toincrease in an amorphous solid was anticipated in [29].Finally, we note from (A6) that in contrast to a crys-talline solid the thermal (caloric) equation of state hereremains ambiguous because the value of c M depends onthe preparation protocol. With this in mind it becomesfruitless to seek the anharmonic corrections to equation(A6) as in the case of a one component system with awell-defined reference state at low temperatures. Never-theless we stress that the specific heat at constant volumedoes not suffer from any ambiguity and therefore can betaken good as a good indicator of the solidification. [1] P.F. Sullivan and G. Seidel, Phys. Rev. , 679 (1968).[2] C.A. Angell and W. Sichina, Ann. N.Y Acad. Sci. ,53 (1976).[3] N.O. Birge and S.R. Nagel, Phys. Rev. Lett. , 2674(1985).[4] E. Lerner and I. Procaccia, Phys. Rev. E , 020501(2008)[5] H. Shintani and H. Tanaka, Nature Physics , 200 (2006).[6] E. Lerner, I. Procaccia and I. Regev, “Quantitative The-ory of a Time-Correlation Function in a One-ComponentGlass-Forming Liquid with Anisotropic Potential”, sub-mitted to Phys. Rev E. Also arXiv:0806.3685.[7] A. Geiger and H.E. Stanley, Phys. Rev. Lett. , 1749(1982).[8] R.V. Chamberlin, Phys. Rev. B , 15638 (1993).[9] D. Kivelson and G. Tarjus, Phil. mag. B , 245 (1998)[10] G. Tarjus, D. Kivelson and P. Viot, J. Phys. Condens.Matter , 6497 (2000).[11] H.G.E. Hentschel and I. Procaccia, Phys. Rev. E, ,031507 (2008).[12] D. Deng, A.S. Argon and S. Yip, Philos. Trans. R. Soc.London Ser A , 549, 575,595, 613 (1989).[13] G.S. Grest and S.R. Nagel, J. Phys. Chem. , 4196(1987).[14] D.N. Perera and P. Harrowell, Phys. Rev. E , 5721(1999) and references therein.[15] E. Aharonov, E. Bouchbinder, H.G.E. Hentschel,V. Ilyin, N. Makedonska, I. Procaccia and N. Schupper, Europhys.Lett. , 56002 (2007) .[16] H. G. E. Hentschel, V. Ilyin, N. Makedonska, I. Procacciaand N. Schupper, Phys. Rev. E , 0050404(R) (2007) .[17] V. Ilyin, N. Makedonska, I. Procaccia, N. Schupper,Phys. Rev. E , 052401 (2007).[18] F.G. Padilla. P. Harrowell, H. Fynewever, J. Non-Cryst.Solids bf 307-310, 436 (2002).[19] W.G. Hoover, M.Ross, K.W. Johnson, J. Chem. Phys. ,4931 (1970).[20] W.G. Hoover, M. Ross, Contemp. Phys.