A thermodynamic description of the hysteresis in specific-heat curves in glass transitions

By refining the definition of thermodynamic equilibrium and state variables (thermodynamic coordinate, TC) for solids, it is determined that the state of a glass substance transforms into an equilibrium state after it is solidified. In contrast, the state of a glass substance during the glass transition is a nonequilibrium state. The specific-heat (C) versus temperature (T) curve exhibits hysteresis, which is traditionally believed to invalidate thermodynamic methods. However, the glass transition slowly occurs in a manner such that structural change is decoupled with the fast process of thermal relaxation of phonons, which enables us to describe the hysteresis by thermodynamic methods. The hysteresis is caused by the structural relaxation and the time of relaxation is determined by the energy barrier, which depends solely on the current value of TCs. Therefore, the state in hysteresis can be described by the information of the current structure alone: history-dependent response functions are unnecessary. On the basis of these conclusions, the behavior of the C-T curve with changing heating/cooling rate γ is simulated. The main features of the hysteresis, the shift of C to higher temperatures with increasing γ, the hump structure, and the memory effect are well reproduced from a structure-dependent energy barrier. In view of the structural dependence of the energy barrier, it is not surprising to observe deviations from the Arrhenius law. However, only the terms that are higher than linear in T are observed in Arrhenius plot as the deviation. An important finding of this study is that the apparent energy barrier obtained using the Arrhenius plot significantly overestimates the real value. The extraordinarily large values of the pre-exponential factor of the relaxation time can also be understood on this basis.


Introduction
The mechanism of glass transition has been a central issue in glass physics for over a century. A difficult problem regarding the transition is the thermal-history dependence on its properties. A prominent example for this is the hysteresis appearing in the specific-heat (C) versus temperature (T) curve. Although hysteresis itself does not necessarily imply irreversibility, the hysteresis of the present case is directly linked to irreversibility, because the transition accompanies a spontaneous relaxation in the structure. Although early studies on the timedependence of enthalpy H were successful in capturing the nature of relaxation [1,2], the nonlinear [3] and nonexponential [4] behaviors made the analysis very difficult. The most successful treatment may be the Tool-Narayanaswamy-Moynihan method [5][6][7][8][9]. By refining the formula with a few parameters, the method becomes a standard for analyzing C-T curves. However, the nonequilibrium nature of the glass transition hides the physics behind these phenomenological parameters. A difficulty in interpreting very large activation energies leaves these parameters as no more than fitting parameters [10]. The current trend of study has shifted more towards microscopic dynamics [11][12][13][14][15][16][17]. Although many developments, such as dynamic heterogeneity, modecoupling theory, and random first-order transition theory have been achieved, these microscopic theories have yet to provide a good explanation for the overestimation of activation energy.
Currently, the view that glass states are nonequilibrium states is standard for many researchers. However, it is strange to speak of equilibrium while lacking an understanding of what it should imply. The author approaches the issue of glass states by framing it via the fundamental problem of equilibrium in thermodynamics. In a previous paper [18], an accurate definition of equilibrium for solids and its relationship to state variables (thermodynamic coordinates, TC), were established on rigorous grounds. There are as many equilibrium states for a given solid as the number of the constituent atoms, and accordingly the thermodynamic states of a solid are expressed by the time-averaged positions R j {¯} of all the constituent atoms. This foundation of equilibrium and TCs is substantiated in glass in [19] (hereafter called paper (I)). In this study, it was shown that (A1) an enormously large number of glass states are indeed equilibrium states once they are solidified: different states can be thermodynamically identified by the time-averaged positions R j {¯} of the atoms. On the other hand, the states of a glass substance during the glass transition are nonequilibrium states, as is widely accepted, because R j {¯} are changing with time t. Despite this, during the transition, the changes in t R j {¯( )}are decoupled with the change in their vibrational part t u j {¯( )}, so that temperature T(t) is well defined at each time t; namely, temporal equilibrium holds. It follows from this that (A2) thermodynamic quantities for the transition state are well defined as a function of t, and can be expressed solely by the current values of t R j {¯( )}, irrespective of the past history. In addition, it has been shown, from the positive-definite property of specific heat, that (A3) the equilibrium state in the transition region is not a supercooled-liquid state, but is an inhomogeneous mixture of liquid and solid phases.
In this paper, the hysteresis in the C-T curve during the glass transition is studied, on the basis of the previous conclusions (A1) to (A3). The analysis contrasts quite strongly to traditional methods, which assume historydependent response functions. The significant consequence of the present approach is that the activation energy and the pre-exponential factor for the relaxation time can be reasonably interpreted on physical grounds, which are standard in solid state physics for crystalline solids. The present paper assumes the background knowledge of the previous papers [18,19]. However, for those readers who are reluctant to read further references, a glossary of the keywords in Supplemental material may be useful.

2.
Analysis of a C-T curve in the heating process 2.1. Basic assumptions Let us apply the theory of paper (I) to the glass transition in a C-T curve. In this paper, C refers to isobaric specific heat, although the subscript p is omitted. The transition in a heating process begins at the lower bound T g,1 and ends at the upper bound T g,2 (figure 1 in paper (I)). The C-T curve in the transition region depends on the rate of change in temperature, γ. This dependence has been treated by many authors in a traditional manner, where history-dependent functions were phenomenologically introduced [7,8,[20][21][22][23].
Assumption (A1) states that the fundamental relation of equilibrium (FRE) for the glass state is expressed as a function of the internal energy U and the averaged position of the atoms; = S S U R , j ( {¯}). For the transition state in the glass transition, from assumption (A3), FRE is given by the sum of the liquid and glass (solid) parts, as S = S l + S g , where the subscripts l and g indicate the liquid and glass components, respectively. Similar notations are adapted for an extensive quantity X: X = X l + X g , where X l and X g are the liquid and solid components, respectively, of X. From assumption (A2), in the transition region, all the variables T and R j {¯} are well-defined TCs as functions of time t, leading to a relationship with constraints on the internal energy U l + U g = U, number of atoms N l + N g = N, and volume V l + V g = V.
Here, R j Treating a large number of variables is too great a task for the theory in its infant status. Fortunately, it is not necessary to do so, because many TCs are insensitive to properties of glass: see the discussion on the degeneration of the state space in [18]. As explained in paper (I) (section IIB2), any structure-dependent property = Y Y R j ({¯}) can be used as a TC. Here, we chose enthalpy H. This choice of TC may be the best choice because the dependence of H on the structure R j {¯} has rigorous grounds based on density functional theory. In addition, it is particularly suitable for comparing with previous studies on the relaxation of the glass transition, many of which have focused on enthalpy relaxation [7,8,10,20,21,24]. The state of the system is thus expressed in a two-dimensional state space (T, H), as shown in figure 1. Only when the system is in thermodynamic equilibrium is there a unique relationship between them: H(T) = H (eq) (T). In order to describe relaxation processes, the time-dependent enthalpy is separated by a structural H s and phonon part H t as in the sprit of the adiabatic approximation of the second kind (equation (10) in paper (I)). Here, R is a collective notation for R j {¯}. The time dependence of H s (R(t)) appears through the time dependence of R(t). Bearing this in mind, in the following H s (R(t)) is simply written as H s (t). On the other hand, H t (T(t)) is a unique function of T only, because the phonon subsystem is always in equilibrium, that is, temporal equilibrium.

Heating process of glass
In this subsection, the glass substance is referred to simply as a system. Consider a heating process with a constant heating rate γ > 0. The heating is accomplished by consecutive injection of heat pulses with period t p . The system is adiabatically isolated from the surroundings, except when heat pulses are injected. Most previous studies considered isothermal processes [7,8,20,21]. But, adiabatic processes are suitable for the present purpose, because it clarifies which part, H t or H s , is responsible to the relaxation. The relaxation occurs by adjusting the portions between H t and H s , while the total H is maintained constant. The state of a system changes to achieve the maximum entropy while the internal energy remains fixed [25]. Because of the adiabatic approximation of the second kind, the temperature T of the system is understood to be the temperature of the phonon subsystem, which is always in equilibrium.
A pulse of heat ΔQ i is injected into the system at t = t i . Just before the heat injection, the temperature of the system is T(t i − ) = T i . The duration t w of the pulse is very short, but much longer than τ t , so that the phonon subsystem can immediately follow the heat input and equilibrates. At this moment, the enthalpy of the system changes with the phonon part H t only; namely, ΔH t (t = t i + ) = ΔQ i , while ΔH s (t = t i + ) = 0 because of the slow response of the structure. Using the adiabatic approximation of the second kind, the increase in T is given by ΔT(t = t i + ) = ΔQ i /C t , so the temperature of the system becomes . Just after the heat injection, the system undergoes an adiabatic change toward the thermodynamic equilibrium state. Because the process is adiabatic, the total enthalpy of the system is preserved, so H = H i+1 ≡ H i + ΔQ i . The total enthalpy H(T) moves along the blue curve in figure 1. Only the portions between H t and H s vary with time t. The second law of thermodynamics says that the state of an isolated system changes so as to maximize the entropy of the system under the constraint of a fixed H. It is a common practice in many phenomena in physics to assume that the rate of change of H s (t) is proportional to the deviation ΔH s from the equilibrium value DH where the proportionality constant τ s is the structural relaxation time. Thus, H t (t) is given as where H T t i eq ( ) ( ) is the enthalpy of the phonon part when the thermodynamic equilibrium is reached for step i. The critical point that differs from previous studies is that the equilibrium state into which the system transitions is considered to be the mixed state of solid and liquid that is obtained by cooling as slow as possible; see section IIIC in paper (I). The adiabatic change causes a change in T given by For a given H i+1 , there is one and only one thermodynamic equilibrium state that satisfies . Before reaching the equilibrium state, the next heat pulse ΔQ i+1 is injected at t = t i+1 = t i + t p . The final temperature after heat pulse i is then T i+1 = T(t = t i + t p ) in equation (5). The specific heat at this step is This process continues until T becomes greater than T g,2 .
For cooling processes, the same procedure is applied, except the relaxation process proceeds toward restoration of the liquid state.

Activation energy
Another important feature of the present model lies in the determination of the relaxation time τ s . In usual relaxation phenomena, such as diffusion in solids, the relaxation time τ obeys the thermal activation law where E a is the activation energy and A is an appropriate constant. Normally, E a is considered to be independent of T. In the glass transition, the assumption of constant E a leads to poor agreement with experiment. The response of glass to external perturbations has been known since the early days to be nonlinear [3] and nonexponential [4]. To treat these properties, Narayanaswamy et al replaced the factor e t s in equation (4) with a time-dependent response function f ¢ t t , ( ) [6,20,26]. They took the thermal history into account through the form τ = τ(T, T * ), and the structural dependence is accounted for by introducing the fictive temperature T * , giving By using the reduced time ò z t = ¢ dt , the expression of equation (8) is rewritten as In accordance with the custom in the study of dielectric relaxation, the response function takes the form where β is a constant, with 0 < β 1 [27,28]. There remains the task of determining τ(T, T * ). Moynihan et al did this by replacing E a in equation (7) by where x and Δh å are fitting parameters [5-8, 26, 29]. This equation is known as the Tool-Narayanaswamy-Moynihan (TNM) formula. In the formula, Δh å is regarded as the activation energy instead of E a . In this manner, the calculation of the relaxation process became tractable, and this method using the four parameters A, Δh å , x, and β became a standard approach for today's glass research. However, by examining use of this fourparameter model for a wide range of glasses, Hodge concludes that no single set of these parameters provides satisfactory fitting for glasses with various histories [10].
In contrast with the previous studies, we retain the original exponential decay in equation (4), and accordingly E a in equation (7) is maintained as the activation energy. Instead of this recurrence to the standard formula equation (7), we assume that E a and therefore τ s depend on the structure of the solid at the present time, according to theorem 2 in paper (I). The activation energy for the glass transition is in a wider sense the energy barriers for atom migration and trapping. The behavior of atoms of glass in the transition region is expressed by an energy landscape, which is shown schematically in figure 2. As the glass substance cools, the energy barrier E b develops around atoms and impedes their migration. The calculation of its structural dependence, ,0 eq ( ) ( ) , and H g,0 eq ( ) is the enthalpy of the glass at T = T g,1 . Setting E b (H) = 0 for > H H l,0 eq ( ) is, of course, an approximation for a small value E b in that region. As usual in chemical reaction, the energy barrier varies depending on the direction of reaction, and hence E b is different for cooling and heating processes. Here, this difference is ignored for the sake of simplicity. In this manner, the complicated task of tracing the thermal history is tremendously simplified. In principle, there are as many τ j as N at . As usual, we proceed as far as possible with a simple model that uses a single relaxation time.

Simulation result
3.1. Dependence on cooling or heating rate We now present the results of the simulation of the glass transition. The specific heat C (eq) of the equilibrium state is prescribed for a hypothetical glass. Given the energy barrier E b for the structural relaxation in the form of equation (11), how the specific heat C depends on the cooling or heating rate γ was simulated. Figure 3 shows the assumed C (eq) -T curve for the hypothetical glass in the equilibrium state, which is supposed to be obtained by cooling at the slowest possible rate: γ → 0. The specific heat C (eq) (T) is assumed to be linear in T in the transition region. Around T g,2 , the curve for C (eq) connects smoothly to C l . The curve for C (eq) as a function of   obtaining reasonable results. E b0 = 0.2 eV was assumed in the following simulations. Using twice or half of this value would force us to use unrealistic values of γ to obtain a reasonable dependence of the C-T curve on γ. In this study, we assume a range of γ from 0.1 to 10 K/s. Figure 4(a) shows how the C-T curve depends on the cooling rate γ ( < 0). As |γ| increases, the specific heat C drops faster, and the enthalpy of the glass H g,0 at T g,1 becomes larger. In this manner, the final state of the glass differs from sample to sample when different cooling rates are used. Let us call the deviation from the equilibrium value D = - ,0 ,0 eq ( ) the residual enthalpy of the glass. Despite nonzero ΔH g,0 , the specific heat of the glass is the same once the atom motions are frozen in, because only H t contributes to specific heat C at T < T g,1 . Figure 4(b) shows how the C-T curve depends on the heating rate γ ( > 0) when starting from a sample in the lowest-energy state; that is, ΔH g,0 = 0. As the heating rate increases, the rise in specific heat C shifts toward high temperatures. A hump appears for large γ, which is often found experimentally. This is because the structural relaxation cannot follow high rates of change in temperature T. The requirement of the closed relation equation (20) of paper (I) raises the C near T g,2 to compensate for the delay in the onset T g,1 . This creates a hump near T g,2 when γ is large. In most experiments, the hump appears near T g,2 . However, the hump is observed even far from T g,2 , where the heating rate is very small [21]. The present simulation indicates that this situation in fact occurs. Figure 5 shows how the C-T curve depends on the residual enthalpy H g,0 . The sample is heated from the glass state at γ = 1 K/s. As the residual enthalpy ΔH g,0 increases, the raising of C is delayed. This relationship between the residual enthalpy and the delay in raising of C is known as a memory effect in glasses, because the current C − T curve reflects the past history [6,7]. When ΔH g,0 is large, the specific heat C initially decreases as T increases. This behavior is observed experimentally for samples prepared at a very fast cooling rate [20], which implies a large residual enthalpy.
In this manner, the main behaviors in the C-T curve as γ is changed have been reproduced. One concern is the shape of the hump feature in the C-T curve. When the sample is heating, experimental observations indicate that, subsequent to reaching the maximum in C, the C converges monotonically to C l as T increases further.  However, the present simulations produce somewhat oscillatory behavior in C-T in that region. This issue is examined in the next section.

Barrier spectroscopy
We now do the calculation in the reverse direction to deduce the barrier height E b (H) from an observed C-T curve. This approach may furnish a new method of spectroscopy to resolve how the barrier height E b depends on the enthalpy H. Figure 6(a) shows the input data for the C-T curve for a hypothetical glass. The specific heat C is obtained by first cooling the sample as slow as possible from the liquid state to ensure that the specific heat C is the equilibrium value C (eq) . After reaching the glass state, the sample is reheated at a finite rate γ. Although this curve does not represent any real data, it was prepared to meet the following conditions: (i) it exhibits a hump in the heating process, and (ii) it satisfies the closed condition equation (20) of paper (I) over a single heating cycle, which ensures that the same sample is used.
Given the heating rate γ = 1 K/s, E b is calculated and the result is plotted in figure 6(b). An unexpected feature of this spectrum is the presence of a plateau in E b before the hump in C; otherwise E b decreases monotonically with T, as expected. Whether this is due to unrealistic input or to the inadequacy of the model to describe the relaxation process is unclear-remember that only a single relaxation time is assumed. Presently, the author supposes that the problem is caused by the present treatment of the single TC: that is, the mean value of H is used for the two-component system; the liquid part H l and the solid part H g should be separated. The observed specific heat C is the sum of the glass part and the liquid part. When x mole fraction of the substance is in the liquid state, the observed specific heat C is given by the sum = - As T is raised, the mole fraction x of the liquid increases, and the total C is dominated by the liquid part, which has no energy barrier. Thus, applying equation (11) to the total C leads to an overestimate of E b near T g,2 . Equation (12) must be taken into account to improve the simulation, which is left for a future study.

Discussion of activation energy
Although a discussion of the details of the properties of glasses is not the main purpose of this paper, a few comments on the present results are warranted. A significant question in the current study of glasses is the deviation of the activation energy from the Arrhenius law. Depending on the magnitude of the deviation, glasses are classified as either strong or fragile [31][32][33]. From the present view, the change of the energy barrier E b with T is evident from the outset because E b is a function of atom position, E b = E b (R). Despite this, it is rather surprising, in the glass literature, that this is not reflected in the analysis of the Arrhenius plot: many authors have attempted to interpret the activation energy as a constant quantity for a given glass [10,23,[34][35][36][37][38][39]. The temperature dependence of the activation energy becomes clearer in molecular-dynamic simulations; for example, Debenedetti and Stillinger showed that the basins become deeper as T decreases [30]; Han et. al. more clearly show this dependence [40]. For chemical reactions, the structures of the reactants and products do not change, and there is no reason that E b should change. The question, therefore, is why E b appears not to deviate much from the Arrhenius law, as is expected: for strong glasses the Arrhenius law in fact holds well. The answer lies in the manner in which the Arrhenius plot is displayed. Equation (11) can be approximated as where b = E b0 /ΔT g is a constant. This approximation shows that E b varies significantly from E b0 to zero in a narrow range of T. However, in the Arrhenius plot, the term linear in T is canceled by the denominator in the exponent of equation (7). Thus, equation (7) becomes * . This form shows why the linear term disappears from the Arrhenius analysis. Importantly, equation (13) explains the large discrepancies in the related quantities. First, the apparent activation energy obtained from the Arrhenius plot is E b *, not E b . In the literature, the activation energies for the glass transitions of organic glasses are reported to be over 1 eV [10,21]; activation energies greater than 10eV are not rare. Although the definition of activation energy slightly differs from method to method, these values are too large to be reasonably understood, as noted by Hodge [41]. The migration enthalpy of a vacancy in silicon is about 0.1 eV [42]. The energy barrier for the diffusion in metals is at most a few eV [43]. Based on the empirical rule for the activation enthalpy of diffusion, E b is given by where K 1 = 1.5 × 10 −3 eV mol −1 K −1 and T m is in K ( [43], p.144). This gives a value E b = 0.45 eV/mol for a material having the melting temperature T m = 300 K. Although the relaxation in the glass transition is different from diffusion processes, both involve atom movements and hence the activation energies must be similar to each other. It is unrealistic that soft organic glasses, with T m often less than room temperature, have energy barriers larger than those of hard solids, even before complete solidification. In contrast, in section 3.1, E b = 0.2 eV was obtained, which indeed falls in a reasonable range in accordance with the empirical equation (15). The reason for this large discrepancy is that the apparent activation energy E b * obtained by the Arrhenius plot enlarges the true value by a factor b. The factor b = E b0 /ΔT g is very large because of the small value ΔT g : e.g. b ∼ 100 for the present case. Thus, the term bT g,1 becomes 2.8 eV for T g,1 = 280 K, which exceeds the original value E b0 by one order of magnitude.
Second, although the linear term is eliminated by the denominator, the effect of the linear term appears in the pre-exponential factor A * = Ae b . Although, even in the standard case for the Arrhenius analysis, it is not rare to find a quantitative discrepancy in the pre-exponential factor, a reasonable interpretation is possible. For example, for impurity diffusion, the pre-exponential factor D 0 in the diffusion constant )is related to the jumping frequency and jumping distance. This gives an order-ofmagnitude agreement. Conversely, for the glass transitions, analyses [10,21] show that A is on the order of 10 100 s −1 . No models can explain this extremely large value for A by combining the elemental frequencies; for example, the typical phonon frequency is only 10 14 s −1 . However, we have already seen that b is very large, of the order of 10 2 . This factor is multiplied by A, giving a total pre-exponential factor A * ≈ e 100 s −1 . This explains why quite large pre-exponential factors appear in the Arrhenius analyses.
Lastly, a comment is made on the relaxation time. In principle, there are relaxation times {τ j } as many as the number of TCs. In the present simulation, all the relaxation times are lumped together in a single parameter τ s . This may not be good approximation, on account of the frequency-dependent relaxation phenomena [11]. Recently, the relation to the low-frequency excitation, which is known as boson peak, is discussed [44]. Inclusion of characteristic modes for the relaxation times {τ j } is desirable for future study.

Conclusions
A new thermodynamic method developed in paper (I) has been applied to the hysteresis appearing in specificheat curves in the glass transition. The hysteresis arises from structural relaxation during the glass transition. This accompanies entropy production, rendering the transition state a nonequilibrium one. Despite this, because the glass transition slowly occurs in a manner such that the structural change is decoupled with the fast process of thermal relaxation of phonons (the adiabatic approximation of the second kind), the thermodynamic state during the transition can be specified solely by the current values of TCs. The current values of TCs determine the current energy barrier E b (R), which in turn determines the relaxation time τ(R). Thus, the hysteresis in specific heat during the transition can be described by the current value of the relaxation time and the heating/cooling rate.
Given a structure-dependent energy barrier E b (R), the present simulation of the C-T curve successfully reproduced its main features when the cooling/heating rate γ is changed. These include the shift C to higher temperatures with increasing γ, the hump structure, and the memory effect. Conversely, it is shown that by analyzing the C-T curve by changing γ, the structure dependence of the energy barrier can be obtained. In view of the structural dependence of E b (R), it is not surprising to observe deviations from the Arrhenius law for glasses. An important finding of this study is that the apparent energy barrier E b * obtained by the Arrhenius plot