Statistical Mechanics of Time Independent Non-Dissipative Nonequilibrium States
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l Statisti al Me hani s of Time Independent Non-Dissipative Nonequilibrium StatesStephen R. Williams ∗ and Denis J. Evans † Resear h S hool of Chemistry, Australian National University, Canberra, ACT 0200, Australia(Dated: 7th November 2018)We examine the question of whether the formal expressions of equilibrium statisti al me hani s an be applied to time independent non-dissipative systems that are not in true thermodynami equilibrium and are nonergodi . By assuming the phase spa e may be divided into time independent,lo ally ergodi domains, we argue that within su h domains the relative probabilities of mi rostatesare given by the standard Boltzmann weights. In ontrast to previous energy lands ape treatments,that have been developed spe i(cid:28) ally for the glass transition, we do not impose an a priori knowledgeof the inter-domain population distribution. Assuming that these domains are robust with respe t tosmall hanges in thermodynami state variables we derive a variety of (cid:29)u tuation formulae for thesesystems. We verify our theoreti al results using mole ular dynami s simulations on a model glassforming system. Non-equilibrium Transient Flu tuation Relations are derived for the (cid:29)u tuationsresulting from a sudden (cid:28)nite hange to the system's temperature or pressure and these are shown tobe onsistent with the simulation results. The ne essary and su(cid:30) ient onditions for these relationsto be valid are that the domains are internally populated by Boltzmann statisti s and that thedomains are robust. The Transient Flu tuation Relations thus provide an independent quantitativejusti(cid:28) ation for the assumptions used in our statisti al me hani al treatment of these systems.I. INTRODUCTIONThe formal expressions of equilibrium statisti al me hani s stri tly apply only to ergodi systems that are inthermodynami equilibrium. Thus these expressions only stri tly apply to systems whi h are at the global free energyminimum given the system Hamiltonian and the ma ros opi thermodynami state variables (number of parti les,temperature and pressure or density). For su h systems Gibbsian equilibrium statisti al me hani s provides an exa tpres ription for how to al ulate the various thermodynami quantities1. However, these pres riptions are routinelyapplied to systems that are not in true thermodynami equilibrium (for example to metastable liquids2, glasses3,polymorphs4 and allotropes). It is often observed empiri ally that within experimental un ertainties many expressionsfor thermodynami quantities yield onsistent results. In the present paper we provide arguments for why many ofthe results of equilibrium statisti al me hani s an be applied to su h time independent nondissipative nonequilibriumsystems. We also point out some of the limits inherent in the appli ation of the formulae of equilibrium statisti alme hani s to su h systems.We hoose to study the isothermal isobari ensemble5 (externally regulated pressure and temperature). The methodsand reasoning we use here an be dire tly transferred to other ensembles su h as the anoni al ((cid:28)xed volume andexternally regulated temperature). The Gibbs free energy G , whi h is the thermodynami potential for the isothermalisobari ensemble, is related to the partition fun tion ∆ by the equation G ( N, P , T ) = − k B T ln ∆( N, P , T ) , (1)and the partition fun tion is given by the integral ∆ = Z Z D dV d Γ exp[ − β ( H ( Γ ) + P V )] , (2)where Γ = ( q , p ) is the phase spa e ve tor des ribing the oordinates q and momenta p , of all the N parti les inthe system, P is the thermodynami pressure, and β = 1 /k B T where k B is Boltzmann's onstant and T is thetemperature. The integration domain D provides limits for both integrals and extends over all the available phasespa e ( Γ , V ) . This is ±∞ for every omponent of the generalized momentum, → ∞ for the volume V , and over thevolume for the Cartesian oordinates of the parti les. Sin e the system Hamiltonian H ( Γ , V ) is single valued, so toois the partition fun tion and in turn the free energy.If we require the distribution fun tion of a single thermodynami phase it is ne essary that other phases do not ontribute signi(cid:28) antly to the partition fun tion. The full integration domain D may in lude states that are hara -teristi of rystalline states or (cid:29)uids states. In the thermodynami limit this does not ause problems be ause, as weshall see, the partition fun tion will be ompletely dominated by those mi ros opi domains that have the lowest freeenergy. However the appli ation of these formulae to allotropes or metastable systems does present a problem. Thestandard equilibrium statisti al me hani al expressions for variables su h as the enthalpy I , the average volume h V i and se ond order quantities su h as the spe i(cid:28) heat at onstant pressure c P may all be omputed from a knowledgeof the partition fun tion Eq. 2 or equivalently the thermodynami potential Eq. 1. If other phases of lower freeenergy exist this omputation (from Eq. 2 as written) will stri tly speaking be in orre t.It is well known that the formulae for thermal properties su h as entropy, free energy, temperature and spe i(cid:28) heatdo not hold for dissipative nonequilibrium systems outside the linear response regime6,7. In this paper we examinethe question of whether they are orre t for any nondissipative nonequilibrium systems su h as allotropes, metastablesystems or history dependent glasses. We provide a statisti al me hani al theory of time independent, nondissipative,nonequilibrium systems. The theory is based on the fa t that these systems are nonergodi and individual samplesystems omprise ergodi domains that do not span all of phase spa e. We show that if these domains are robustwith respe t to small hanges in thermodynami state variables, a su essful statisti al me hani al treatment ofthese nonequilibrium systems an be given. We provide dire t eviden e, from mole ular dynami s simulations on amodel glass former, that the resulting statisti al me hani al formulae are satis(cid:28)ed within empiri al errors. Finallywe provide an independent test of the two key elements of our theory: Boltzmann weights within the phase spa edomains and the robustness of those domains. It happens that these two elements are the ne essary and su(cid:30) ient onditions for the appli ation of the transient (cid:29)u tuation relation to (cid:28)nite thermodynami quen hes (in temperatureor pressure) for su h systems8. While the appli ation of thermodynami s to a single time averaged system is usuallystraightforward the appli ation to a ensemble, whose members may be lo ked in di(cid:27)erent phase spa e domains, anrequire modi(cid:28) ation to the standard formulae.In the ase of glasses our treatment has some similarities with the energy lands ape approa h of Stillinger andWeber3,9,10. However there are signi(cid:28) ant di(cid:27)eren es; our treatment makes no referen e to the inherent stru tureand imposes no a priori knowledge of the inter-domain relative population levels. The energy lands ape approa h hasbeen extended to a ount for the phenomena of ageing or history dependen e by the addition of a (cid:28) tive parameter11.S iortino has onvin ingly shown that the addition of a single (cid:28) tive parameter is inadequate to deal with glasses,whi h may have di(cid:27)erent properties at the same temperature and pressure if they are prepared by a di(cid:27)erent proto ol(di(cid:27)erent history dependen e)11 and poses the hallenge to re over a thermodynami des ription (cid:16)by de omposing theageing system into a olle tion of substates(cid:17). The treatment we present here su eeds in doing just that by providinga rigorous development of equilibrium statisti al me hani s and thermodynami s for ensembles of systems where thephase spa e breaks up into ensembles of domains whose inter-domain dynami s is nonergodi and whose inter-domainpopulation levels may not be Boltzmann weighted.II. CONDITIONS FOR EQUILIBRIUMA dynami al system in equilibrium has the properties that it is nondissipative and that its ma ros opi propertiesare time independent. Thus the N-parti le phase spa e distribution fun tion, f ( Γ , V, t ) , must be a time independentsolution to the Liouville equation7, ∂∂t f ( Γ ′ , t ) = − ˙Γ ′ · ∇ f ( Γ ′ , t ) − f ( Γ ′ , t )Λ = 0 , (3)where Λ is the phase spa e ompression fa tor7 obtained by taking the divergen e of the equations of motion (seeEq. 5) and Γ ′ is the extended phase spa e ve tor whi h onsists of Γ and may in lude additional dynami al variablessu h as the volume V . Sin e the system is assumed to be nondissipative both the ensemble average h Λ i and thetime average Λ of the phase spa e ompression fa tor (whi h is dire tly proportional to the rate at whi h heat isex hanged with the (cid:28) titious thermostat) are zero. The time independent solution to Eq. 3 depends on the detailsof the equations of motion. Equilibrium solutions to Eq. 3 for the equations of motion, suitable for use in mole ulardynami s simulations, are ompatible with Gibbsian equilibrium statisti al me hani s7.Mi ros opi expressions for me hani al properties like the pressure, the internal energy, the enthalpy and the volume an be derived without referen e to Gibbsian statisti al me hani s and indeed an be proved to hold for nonequilibriumsystems in luding nonequilibrium dissipative systems.There are two ways in whi h the formulae derived from Gibbsian equilibrium statisti al me hani s an breakdown. The most obvious way is that the relative weights of mi rostates may be non-Boltzmann and the exponentialfa tor exp[ − βH ( Γ )] , may be repla ed by some other fun tion (either the exponential fun tion itself may be modi(cid:28)edas in Tsallis statisti s12 or the Hamiltonian may be modi(cid:28)ed to some new fun tion H ( Γ ) → B ( Γ , t ) H ( Γ ) ). Ineither ir umstan e the standard expressions for the thermal quantities derived from equilibrium Gibbsian statisti alme hani s will not be valid. This ertainly happens in dissipative nonequilibrium systems where the distributionfun tion is not a time independent solution to Eq. 3.In deterministi nonequilibrium steady states the phase spa e may break down into ergodi ally separated domains(Ea h of whi h will be fra tal and of lower dimension than the ostensible phase spa e dimension. This is a onse-quen e of dissipation.) However for these steady states, the domains are always exquisitely sensitive to ma ros opi thermodynami parameters sin e they are strange fra tal attra tors13. Often a deterministi nonequilibrium steadystate approa hes a unique fra tal attra tor. As time progresses the distribution fun tion ollapses ever loser to (butnever rea hing) the steady state attra tor.The se ond way that these expressions may fail is that the system may be ome nonergodi . In this ase three thingshappen. a) Most obviously time averages no longer equal full (domain D ) ensemble averages. b) If we take an initialmi rostate the subsequent phase spa e traje tory will span some phase spa e domain D α where the initial phase islabeled Γ ′ α (0) . In this ase for nondissipative nonequilibrium systems where the domains are robust (i.e. small hangesin thermodynami state parameters, to leading order do not hange the domain) the standard equations of equilibriumstatisti al me hani s may ontinue to be valid but in a slightly modi(cid:28)ed form. We will examine this in some detailbelow. ) Given robust domains the population densities between ea h domain may well depend on the history of thesystem. The ma ros opi history an be expe ted to ondition the ensemble's set of initial mi rostates { Γ ′ α (0) } fromwhi h the ma ros opi material is formed. This in turn an be expe ted to ondition the set of nonergodi domains { D α } that hara terize the ensemble. For a ma ros opi sample spanning a single ergodi domain D α , the free energy G α then satis(cid:28)es only a lo al extrema prin iple and thus looses mu h of its thermodynami meaning.III. THEORY AND METHODSA. Equations of MotionWe use the onstant pressure Nosé-Hoover equations of motion by ombining the Nosé-Hoover feedba k me hanismwith the so- alled SLLOD or DOLLS equations of motion, whi h are equivalent for dilation. It is known that theseequations of motion do not produ e artifa ts in the systems linear response to an external (cid:28)eld and that to leadingorder the e(cid:27)e t on the dynami al orrelation fun tions is at most O (1 /N ) , where N is the number of parti les7. Theequations of motion are, ˙ q i = p i m + α V q i ˙ p i = F i − α V p i − α T p i ˙ α V = (cid:18) V ( t ) N k B T ( P ( t ) − P ) + 1 N (cid:19) /τ V ˙ α T = P Ni =1 p i · p i mN k B T − N ! /τ T ˙ V = 3 α V V, (4)where q i is the position, p i is the momentum and F i is the for e on the i th parti le, m is the parti le mass, τ V is the barostat time onstant, τ T is the thermostat time onstant, T is the input temperature, P is the input(thermodynami ) pressure and the instantaneous (me hani al) pressure is given by P ( t ) = ( P Ni =1 p i · p i /m + P Ni =1 F i · q i ) / V . Be ause these equations of motion have additional dynami al variables the extended phase spa e ve tor is Γ ′ = ( Γ , V, α V , α T ) . In order to obtain the equilibrium distribution fun tion we (cid:28)rst de(cid:28)ne the Hamiltonian, inthe absen e of any external (cid:28)elds, dilation α V ( t ) = 0 or thermostats α T ( t ) = 0 , as H = Φ + P Ni =1 p i · p i /m ,where Φ is the total inter-parti le potential energy. To pro eed further we identify the extended Hamiltonian as H E = H + N α T τ T k B T + N α V τ V k B T and then obtain the phase spa e ompression fa tor Λ ≡ ∇ · ˙ Γ ′ = N X i =1 3 X γ =1 ∂ ˙ q i,γ ∂q i,γ + N X i =1 3 X γ =1 ∂ ˙ p i,γ ∂p i,γ + ∂ ˙ V∂V + ∂ ˙ α V ∂α V + ∂ ˙ α T ∂α T = β ( ˙ H E + P ˙ V ) , (5)where the index γ sums over the omponents of the Cartesian position and momentum ve tors. Using the Heisenbergstreaming representation (rather than the more usual S hrödinger representation Eq. 3) of the Liouville equation ddt ln (cid:2) f ( Γ ′ ( t ) , t ) (cid:3) = − Λ (cid:0) Γ ′ ( t ) (cid:1) , (6)we an obtain the parti ular time independent solution for the distribution fun tion f ( Γ ′ ) ∝ exp( − βI ) exp( − N ( α T τ T + α V τ V )) , (7)where I ( t ) = H ( t ) + P V ( t ) is the instantaneous enthalpy. The se ond exponential on the RHS of Eq. 7 with α V and α T in the argument, whi h has no dependen e on the input temperature T or the input pressure P , isstatisti ally independent from the rest of the distribution fun tion, whi h is the standard equilibrium isothermalisobari distribution. We an normalize Eq. 7 by integrating over all spa e to obtain the thermodynami equilibriumdistribution fun tion, f ( Γ ′ ) = 32 N τ V τ T π exp( − N ( α T τ T + α V τ V )) f ( Γ , V ) , (8)where the standard isothermal isobari distribution fun tion is f ( Γ , V ) = exp ( − β ( H + P V )) R ∞ dV R D d Γ exp ( − β ( H + P V )) . (9)It should be emphasized that the derivation of Eq. 7 says nothing about the existen e or otherwise of any domains.These must be onsidered when normalizing Eq. 7 and thus Eqs. 8 & 9 are only valid in thermodynami equilibrium.If we wish to use Eq. 7 outside thermodynami equilibrium we must onsider domains.We an also use the so alled SLLOD equations of motion7 to apply strain rate ontrolled Couette (cid:29)ow (planarshear) to our equations of motion. The ne essary modi(cid:28) ations to the (cid:28)rst two lines of Eq. 4 result in ˙ q i = p i m + α V q i + i ˙ γq yi ˙ p i = F i − α V p i − α T p i − i ˙ γp yi , (10)where ˙ γ is the strain rate and the last three lines of Eq. 4 remain un hanged.B. Equilibrium Statisti al Me hani s in a Single DomainAs we have stated in the introdu tion, the full phase spa e domain in ludes phase points from many di(cid:27)erentthermodynami phases (gases, liquids and rystals). In the thermodynami limit this does not ause problems. Tounderstand this suppose we an label mi rostates to be in either of two possible thermodynami phases 1 or 2 bound bytwo phase spa e domains D and D . By assumption we are not presently onsidering the possibility of o-existen e.The system is assumed to be ergodi : atoms in one thermodynami phase an in time, transform into the other phase.Assume the two thermodynami phases have di(cid:27)erent free energies: G is the Gibbs free energy of the (cid:28)rst phase and G is that of the se ond phase. For a su(cid:30) iently large N the free energy Eq. 1 is an extensive variable. We may thusexpress the partition fun tion as the sum of ontributions from the two phases ∆ = e − βG + e − βG = e − βNg + e − βNg , (11)where the lower ase g on the se ond line is used to represent the intensive free energies whi h do not hange withsystem size N in the thermodynami limit. If g is less than g then in the thermodynami limit, N → ∞ , theonly signi(cid:28) ant ontribution to the partition fun tion ∆ will be due to the (cid:16)equilibrium(cid:17) phase namely phase 1. Thusalthough the free energy de(cid:28)ned in Eq. 2, is given by an integral over all of phase spa e D , in the thermodynami limitthis integral an be approximated to arbitrary pre ision, as an integral over the domain that in ludes the most stablephase. Suppose D in ludes only rystalline phases and D in ludes only amorphous phases and further suppose aparti ular rystalline phase has a lower free energy that any amorphous phase. A ording to Eq. 2, we should al ulatethe free energy by integrating over all rystalline and all amorphous phases. In pra ti e in the thermodynami limitwe an ompute the free energy to arbitrary a ura y by integrating Eq. 2, only over that part of phase spa e withinwhi h the thermodynami ally stable state resides.If we onsider a nonergodi system that a ording to di(cid:27)erent preparative proto ols an be formed in either phase1 or phase 2. After preparation, be ause the system is nonergodi both phases are kineti ally stable inde(cid:28)nitely. Byrestri ting the phase spa e integrals for the free energy to those domains that ontain the kineti ally stable phasewe an ompute the free energy of that phase. However, although it may be possible to formally assign free energiesto nonergodi systems, these free energies learly fail to satisfy any global extremum prin iple. As we will showthese partition fun tions an be used formally to yield (cid:28)rst and se ond order thermodynami quantities by numeri aldi(cid:27)erentiation. The metastable domain is a subset of the thermodynami equilibrium domain whi h ontains allpossible atom positions in luding ones belonging to the metastable phase.Within a single domain the system is, by onstru tion, ergodi . Thus for almost all mi rostates Γ ′ α (0) ∈ D α ,ensemble averages, of some variable B , h B i , equal time averages B , for phase spa e traje tories that start at timezero, h B i α = B α ≡ lim t →∞ t Z t ds B ( Γ ′ ( s ); Γ ′ α (0)) . (12)Mi ros opi expressions for me hani al variables may be used as a test of ergodi ity in nondissipative systems whi hare out of equilibrium. (Note a nondissipative system does not on average ex hange heat with any thermal reservoirwith whi h it has been in onta t for a long time.) In the ase of metastable (cid:29)uids or allotropes we may introdu e asingle restri ted domain and by onstru tion the system remains ergodi within this domain.A gedanken experiment an be used to justify the Boltzmann weighting and the appli ability of the Zeroth Law ofThermodynami s for su h systems. Consider a double well potential with an inner and outer potential well. If thebarrier between the inner and outer wells is mu h greater than k B T , so that over the duration of observation (whi h ismu h greater than any relaxation time in the ergodi ally restri ted subsystem) no parti les ross the barrier, then thesystem onsidered as a double well system, will be, by onstru tion, nonergodi . For systems omposed of parti lesthat are solely found in the inner potential well, our hypotheses are that the distribution of states in the inner wellwill be given by a Boltzmann distribution taken over the inner domain only and that if su h a system is in thermal onta t with a body in true thermodynami equilibrium, then the temperature of the ergodi ally restri ted systemmust equal that of the system in true thermodynami equilibrium. We an justify these hypotheses by onsidering a(cid:28) titious system that only has the inner potential well and in whi h the potential fun tion is positive in(cid:28)nity for allseparations that are greater than the inner well (this in ludes the position of the outer well). In a ord with Gibbsianstatisti al me hani s the distribution of states is anoni al over this (single well) potential. Furthermore the ZerothLaw of thermodynami s will apply to this single well system. Now if we dynami ally generate the outer well, all theparti les lo ked inside the inner well annot (cid:16)know(cid:17) that the outer well has been formed so their dynami s will be ompletely un hanged by the time dependent generation of the new outer well. The generation of an ina essibleouter well will not alter the distribution of states in the inner well nor will it ause any (cid:29)ow of heat to the equilibriumheat bath surrounding the system. This provides a ompelling physi al justi(cid:28) ation for our domain hypotheses overa single ergodi sub domain of phase spa e.In order to re over many of the basi relationships of Gibbsian statisti al me hani s it is also ne essary that thesystem appears to be in dynami al equilibrium, i.e. | f ( Γ α , t ) − f ( Γ α , t + τ o ) | < ε, ∀ Γ α , ∈ D α , for some small ε , overthe longest observation time τ o . We use the de(cid:28)nition of the partition fun tion Eq. 2 as before but now the domain D α in the integral is over a single ontiguous hypervolume in the on(cid:28)guration spa e of the generalized position oordinate q and volume V . The domain over the generalized momentum p and multipliers α V and α T remainsun hanged. We then obtain the Gibbs free energy by use of Eqs. 1 & 2. Thus far all we have altered is our de(cid:28)nitionof the domain. In hanging the de(cid:28)nition of the domain we have opened a potential problem for Gibbsian statisti alme hani s. If we hange the temperature or the pressure of the system the domain may also hange. If the domain hanges this may make a ontribution to the derivatives of the partition fun tion, Eq. 2, and the dire t onne tionwith the standard out omes of ma ros opi thermodynami s will be lost. Thus the domains need to be robust withrespe t to hanges in thermodynami state variables.There are three means by whi h a system ould have robust domains. The (cid:28)rst and most obvious is that the domaindoes not hange when the pressure or the temperature hanges, ∂D α ( X, Y ) /∂X = 0 , where X is a thermodynami state variable and Y is the other thermodynami state variables. When we lower the temperature only the inversetemperature β in Eq. 7 hanges and when we hange the pressure only the parameter P hanges. If the domain'sboundary is determined by a surfa e on whi h I ( Γ ′ ) always has a very high value it will remain un hanged underin(cid:28)nitesimal hanges in P or β . We will refer to a surfa e domain that doesn't hange with the state variables as ompletely robust. The se ond way is that the distribution fun tion is always identi ally zero on the domain boundary, f ( Γ , V ) = 0 ∀ Γ ∈ S α , (13)where S α is the surfa e of the domain D α . Be ause the domain is ontiguous (required for it to be ergodi ) it musthave a single onne ted surfa e. Su h a domain will be robust. The third way the domain an be robust is lessrestri tive and allows for the possibility that the domain does hange when the thermodynami variables are hangedsubstantially. If δX is an in(cid:28)nitesimal hange in a thermodynami state variable then, δD α ( X + δX, Y ) = δD α ( X, Y ) + O ( δX ) n (14)where Y denotes the other thermodynami state variables, we require that n ≥ for (cid:28)rst order thermodynami property formulae to be orre t n ≥ for se ond order property formulae to be orre t et . Obviously this third waywill be satis(cid:28)ed in the (cid:28)rst two ases as well.Later in the paper we will introdu e an independent test of domain robustness. However, if a system was notrobust then we would expe t that small hanges in the state variables would hange the ma ros opi properties ofthe sample permanently - it would be as though the preparation history of the sample was ontinuing even for small hanges in the state variable. Quite obviously if we produ e huge hanges in the state variables we will of oursepermanently hange the properties of the system be ause we permanently deform the ergodi domain. Experien eshows however that very many nondissipative nonequilibrium systems are quite robust with respe t to small hanges instate variables. All that is required for (cid:29)u tuation formulae for (cid:28)rst se ond and third order thermodynami quantitiesto be valid, is that the domains be un hanged, to (cid:28)rst se ond or third order, by in(cid:28)nitesimal hanges in the statevariables. Obviously a robust domain is an ideal onstru t. However on the typi al time s ale of interest, whi h isusually orders of magnitude less than the time s ale on whi h the system will hange to a new phase of lower freeenergy, this an be a very good approximation.We are now able to re over most of the standard results of Gibbsian equilibrium statisti al me hani s. For examplewe may al ulate the enthalpy h I i α from the partition fun tion ∆ , Eq. 2, as h I i α = k B T (cid:18) ∂ln ∆ D α ∂T (cid:19) = R R D α dV d Γ I ( Γ , V ) exp ( − βI ( Γ , V )) R R D α dV d Γ exp ( − βI ( Γ , V )) = I α . (15)Here we are onsidering an ensemble of systems whi h o upy a single ergodi domain D α . Sin e this domain is selfergodi the ensemble average is equal to the orresponding time average.The term on the RHS of the se ond line is obviously the average value of the instantaneous enthalpy, I ( Γ , V ) = H ( Γ , V ) + P V , obtained by using the equilibrium distribution fun tion, Eq. 9 or equivalently Eq. 8 with theintegration limits restri ted to the domain D α , where P is the externally set thermodynami pressure. We an alsoobtain expressions for the average volume h V i and the onstant pressure spe i(cid:28) heat c P by taking the appropriatederivatives of the partition fun tion Eq. 2. In other ensembles we an use the same pro edure to (cid:28)nd other variablese.g. the internal energy, the average pressure and the onstant volume spe i(cid:28) heat in the anoni al (N,V,T ) ensemble.An important out ome is that this des ription remains ompatible with ma ros opi thermodynami s. Here theGibbs free energy is de(cid:28)ned as, G ≡ U − T S + P h V i , (16)where U = h H i is the internal energy and S is the entropy. If we take the derivative of Eq. 16 with respe t to oneof the isobari isothermal ensembles onjugate variables ( N, P , T ) while keeping the others (cid:28)xed we obtain (cid:18) dGdT (cid:19) N,P = − S (cid:18) dGdP (cid:19) N,T = h V i . (17)We now write down the mi ros opi equilibrium equation for the Gibbs entropy S = − k B Z Z D α dV d Γ f ( Γ , V ) ln f ( Γ , V ) . (18)It is an easy matter to show that Eqs. 1, 2 & 18 are onsistent with the two derivatives given in Eq. 17. Given our ondition of a robust boundary we thus have a form of Gibbsian statisti al me hani s for metastable states whi hremains in agreement with ma ros opi thermodynami s.C. Multiple Domains and Nonergodi ityWe now wish to onsider an ensemble of systems whi h is prepared from an initial ergodi (usually high temperature)equilibrium ensemble. There is some proto ol P , whi h involves a temperature quen h or a sharp pressure in reaseet , whi h breaks the ensemble into a set of sub-ensembles hara terized by di(cid:27)erent ma ros opi properties. Afterthe proto ol P , has been exe uted we allow all the ensemble members to relax to states whi h are ma ros opi allytime independent - to within experimental toleran es. We assume that the ensemble an be lassi(cid:28)ed into a setof sub-ensembles { α, α = 1 , N D } whose ma ros opi properties take on N D distin t sets of values. For the longestobservation times available a ma ros opi system lassi(cid:28)ed as an α system is not observed to transform into a β system, and vi e versa. The full ensemble of systems is thus non-ergodi . However, in ea h individual sub-ensemble,say sub-ensemble α , the onstituent members are ergodi (by onstru tion). Thus we an partition the full phasespa e into a set of domains, { D α } .From the arguments given above (in se tion B), after the relaxation of initial transients, we expe t to observe aBoltzmann distribution of states within an individual domain whi h is therefore independent of the quen h proto ol.However the distribution between domains annot be expe ted to be Boltzmann distributed and will instead bedependent on the quen h proto ol. Within a given ensemble the proportion of ensemble members ultimately foundin domain D α is given by a weight w α ( P ) whi h is subje t to the onstraint N D X α =1 w α = 1 . (19)We an al ulate the full ensemble average of some ma ros opi property B as, h B i = N D X α =1 w α R R D α dV d Γ B ( Γ , V ) exp( − βI ) R R D α dV d Γ exp( − βI ) . (20)Sin e the full ensemble of states is nonergodi the phase spa e breaks up into disjoint domains whi h in themselvesare ergodi . Thus ea h domain may be identi(cid:28)ed by any point in phase spa e, ( Γ , V ) , that is a member of it so thesubs ript α is a fun tion of the phase ve tor α ( Γ , V ) allowing the following expression for the distribution fun tion f ( Γ , V ) = N D X α =1 w α s ( Γ , D α ) f α ( Γ , V ) (21)where s ( Γ , D α ) = 1 if Γ ǫ D α and s ( Γ , D α ) = 0 otherwise, and f α ( Γ , V ) = exp( − βI ( Γ , V )) R ∞ dV R D α d Γ exp( − βI ( Γ , V )) . (22)The entropy is given by S = − k B R ∞ dV R D d Γ f ln( f ) and using Eq. 21 we have the following expressions for themultidomain entropy, S = − k B N D X α =1 w α (cid:20)Z ∞ dV Z D α d Γ f α ln( f α ) + ln( w α ) (cid:21) = N D X α =1 w α S α − k B N D X α =1 w α ln( w α ) . (23)The term − P N D α =1 k B w α ln( w α ) ≡ S D is the inter-domain entropy, whi h is maximized by an even distribution ofensemble members over all domains, while S α is the intra-entropy of domain α onsidered as a single N -parti lesystem.If we substitute Eq. 22 into Eqs. 23 we (cid:28)nd that, S = S D + T − N D X α =1 w α h I i α + k B N D X α =1 w α ln Z ∞ dV Z D α d Γ exp( − βI ( Γ , V ))= S D + T − h I i + k B N D X α =1 w α ln Z ∞ dV Z D α d Γ exp( − βI ( Γ , V )) , (24)where h B i α = R R D α dV d Γ B ( Γ ) f α ( Γ , V ) . Combining Eq. 16 with Eq. 24 we obtain the following expression for theGibbs free energy G = − k B T N D X α =1 w α (cid:20) ln Z Z D α dV d Γ exp( − βI ) − ln( w α ) (cid:21) = N D X α =1 w α G α − S D T. (25)It is easy to verify that if we hold the lo al domain weights (cid:28)xed and then vary the temperature or the pressurethat Eqs. 23 & 25 are ompatible with Eqs. 17. This means that if we have a (cid:28)xed number of robust domains,whose population levels or weights are non-Boltzmann distributed, Eqs. 25 & 23 provide a dire t mi ros opi link tostandard ma ros opi thermodynami s. On the extremely long time s ale the weighting fun tions w α may vary andthe system will tend towards the dire tion where the free energy Eq. 25 is redu ed. Without the inter-domain entropyterm, S D , Eq. 25 would be minimized when the domain with the lowest free energy has all the ensemble members init. It turns out that Eq. 25 is minimized when all the domain weights are Boltzmann distributed, ie when w α = R R D α dV d Γ exp ( − βI ) P N D β =1 R R D β dV d Γ exp ( − βI ) . (26)Here (i.e. upon obeying Eq. 26) the entropy and free energy given by Eqs. 23 & 25 oin ide with the standardequilibrium expressions so the free energy must be a minimum. To prove this we use Eq. 25 and we remove the (cid:28)rstweight w = 1 − P N D α =2 w α , so that the onstraint Eq. 19 is respe ted while the remaining weights are independent.This means that the free energy an be written as G (1 − N D X α =2 w α , w , w , ...w N D ) . The onstrained partial derivativesare then ∂G∂w α (cid:12)(cid:12)(cid:12)(cid:12) c = ∂G∂w ∂w ∂w α + ∂G∂w α = − ∂G∂w + ∂G∂w α , α ≥ − ∂G∂w + ∂G∂w α , α ≥ , = − G − kT ln ( w ) − k B T + G α + k B T ln ( w α ) + k B T, α ≥ . (27)Using the fa t that for Boltzmann weights, Eq. 26, w α = exp [ β ( G eq − G α )] (28)where G eq = − k B T ln P N D α =0 R R D α dV d Γ exp( − βI ) is the equilibrium free energy, we (cid:28)nd that at equilibrium, ∂G∂w α | c = 0 , α ≥ . (29)It remains to prove that this is indeed a minimum. Using the same approa h for treating the onstraint, we ontinuemaking the (cid:28)rst weight a fun tion of all the others, and obtain ∂ G∂w α ∂w γ (cid:12)(cid:12)(cid:12)(cid:12) c = ∂ G∂w − ∂ G∂w γ ∂w α − ∂ G∂w γ ∂w − ∂ G∂w ∂w α . (30)Using the Boltzmann weights it is easy to show that ∂ G∂w = kTw = k B T exp [ β ( G − G eq )] ∂ G∂w ∂w α = ∂ G∂w γ ∂w = 0 ∂ G∂w γ ∂w α = δ γα k B Tw α = δ γα k B T exp [ β ( G α − G eq )] . (31)From these results it is easy to see that the Hessian matrix ∂ G∂w α ∂w γ (cid:12)(cid:12)(cid:12) c is positive de(cid:28)nite, thus we have proved the freeenergy to be a minimum for the ase of equilibrium.We an use a knowledge of the multiple domain thermodynami potential, Eq. 25, to ompute averages. As anexample we onsider the average enthalpy again, h I i = − k B T ∂ βG∂T (cid:12)(cid:12)(cid:12)(cid:12) w α ,P ,N . (32)Eq. 32 an easily be derived from Eq. 16 and is in essen e the same as the (cid:28)rst line of Eq. 15. It is straightforwardto see that upon using Eq. 25 to al ulate the derivative in Eq. 32 one obtains the average enthalpy as given by Eq.20. One an do the same for other quantities su h as the spe i(cid:28) heat et .D. Appli ation to Mole ular Dynami s SimulationWe are now in a position to test various out omes, using omputer simulation, whi h may be derived by drawingon the previous material. If we start with an ensemble of systems whi h are initially in equilibrium at temperature T = T and then at time t = 0 we quen h them by setting T = T we an solve the Liouville equation Eq. 6 to obtain f ( Γ ( t ) , α V ( t ) , α T ( t ) , t ) = f ( Γ (0) , α V (0) , α T (0) , × exp [ β ( I (0) − I ( t ))] (33)where β = 1 /k B T . This nonequilibrium distribution fun tion, valid for t > , expli itly requires the solution of theequations of motion and is of limited utility. However it allows the identi(cid:28) ation of a formal ondition to identifythe amount of time, whi h must elapse after the quen h, before Eqs. 7, 9, 20 or 21 an be applied to the ensemble.That is the quantity I ( t ) = H ( t ) + P V ( t ) must be statisti ally independent of I (0) . Thus we are interested in the orrelation fun tion C ( I ( t ) , I (0)) = h I ( t ) I (0) i − h I ( t ) i h I (0) i C , (34)where C , = r(cid:16) h I ( t ) i − h I ( t ) i (cid:17) (cid:16) h I (0) i − h I (0) i (cid:17) . (35)This fun tion will equal 1 for a perfe tly orrelated system, -1 for a perfe tly anti orrelated system and 0 for anun orrelated system. When we onsider an ensemble of systems, o upying the various domains to di(cid:27)erent degrees,we see that Eq. 34 may not de ay to zero given the traje tories are unable to leave their domains. If the transients,due to the quen h, fully de ay Eq. 34 will de ay to zero and Eqs. 7 & 9 will be ome valid for the ensemble.If the orrelation fun tion, Eq. 34 de ays to a plateau then we may have a situation where Eqs. 20 and 21 arevalid. It may be that the material an still slowly age, due to pro esses, that o ur on a time s ale whi h is longerthan the one we are monitoring. For a glass we expe t that orrelation fun tion Eq. 34 will not fully de ay on areasonable time s ale. If we give the system time to age, su h that it appears to be time translationaly invariant, andthen ompute the following orrelation fun tion, C ( I ( τ ) , I (0)) = P N t α =0 h I ( t ) I (0) i α − h I ( t ) i α h I (0) i α C , (36)where C , = N t X α =0 r(cid:16) h I ( t ) i α − h I ( t ) i α (cid:17) (cid:16) h I (0) i α − h I (0) i α (cid:17) , (37)we may observe a full de ay. If this o urs Eqs. 20 and 21 will be valid. To ompute this orrelation fun tion N t traje tories are produ ed and for ea h of these the averages, (cid:10) ¯ B (cid:11) α , where B is an arbitrary dynami al variable,appearing in Eq. 36 are approximately obtained by time averaging. When the system is ergodi and time translationalyinvariant these two orrelation fun tions, Eqs. 34 & 36, will give the same result. However for a nonergodi system C will de ay to zero on a reasonable time s ale while C will not. Rather C may de ay to some plateau on a reasonabletime s ale and then de ay on a mu h slower time s ale. The pre eding se tions then rest on this lear separation oftime s ales in the orrelation fun tion Eq. 34. For metastable (cid:29)uids and allotropes this separation is so extreme thatwe probably annot observe, even the early stages of, the later slow de ay on any reasonable experimental time s ale.Further for these systems there will be only a single domain and thus they appear ergodi . For glasses some signs ofthe later de ay an often be observed, however it is still very mu h slower than the initial de ay. In the (cid:28)eld of glassydynami s the initial de ay is often alled the β de ay and the slower long time de ay is often alled the α de ay. Asthe glass is further aged this se ond stage de ay is observed to slow down dramati ally while the early de ay does not hange very mu h14.To allow the hypothesis of lo al phase spa e domains to be tested we will now introdu e several relations whosederivation draws upon the equilibrium distribution fun tion Eq. 9. We will also dis uss the e(cid:27)e t of the phase spa edomains on these relations.First we introdu e the on(cid:28)gurational temperature k B T = − h F · F ih∇ · F i (38)0where F is a N dimensional ve tor representing the inter-parti le for es on ea h atom F = −∇ Φ . This relation iseasily derived from Eq. 9 (set B to ∇ · F , integrate by parts and drop the boundary terms) and will remain valid fornonergodi systems where Eq. 26 is not obeyed. The relation involves the spatial derivative of the for e, whi h is notzero at the uto(cid:27) radius for the potential we use in our simulations (see below). This along with (cid:28)nite size e(cid:27)e ts anresult in a small disagreement between Eq. 38 and the kineti temperature for our system when in equilibrium.Equilibrium (cid:29)u tuation formulae may be easily derived from Eq. 20, see ref.1,5 for some examples of how this isdone. We onsider how the average enthalpy hanges with the temperature at onstant pressure, C P = h V i c P = ddT (cid:12)(cid:12)(cid:12)(cid:12) N,P h I i = βT h(cid:10) I (cid:11) − h I i i (39)where c P is the onstant pressure spe i(cid:28) heat. The al ulation of the RHS of Eq. 39 using the ensemble averagegiven by Eq. 20 is subtle. If we simply al ulate h I i and (cid:10) I (cid:11) with the use of Eq. 20 and then plug the results intoEq. 39 we obtain, what we will refer to as, a single domain average whi h does not give us the orre t hange inthe average enthalpy for a history dependent equilibrium ensemble. This is be ause the average (cid:10) I (cid:11) superimposesa ross the di(cid:27)erent domains while the quantity h I i ontains spurious ross terms. If we derive the heat apa ity bytaking the se ond derivative of the Gibbs free energy for the multiple domain distribution fun tion, Eqs. 25 & 32,with respe t to the temperature, C P = − ∂∂T (cid:12)(cid:12)(cid:12)(cid:12) N,P ,w α k B T ∂ βG∂T (cid:12)(cid:12)(cid:12)(cid:12) N,P ,w α , (40)we obtain the following C P = h V i c P = βT N D X α =1 w α h(cid:10) I (cid:11) α − h I i α i , (41)where the quantity (cid:10) I (cid:11) α −h I i α is obtained for ea h domain separately (here h . . . i α represents an average taken whereall ensemble members are in the α th domain). We will refer to this as a multidomain average whi h is onsistentwith the nonergodi statisti al me hani s and thermodynami s that we have introdu ed here. It is obvious thatboth a single and multiple domain average will give the same result in the ase of thermodynami equilibrium andmetastable equilibrium (single domain). The transition from the single domain average produ ing the orre t resultto an anomalous result is symptomati of an ergodi to a history dependent nonergodi transition. If we onsidera large ma ros opi system (the super system) to be made of N s independent subsystems the multidomain averageremains self- onsistent. To see this we apply Eq. 41 to (cid:29)u tuations in the super system and then we inquire howthis relates to (cid:29)u tuations in the subsystem. The enthalpy of one instan e of the super system will be given by I s = P N s α =1 I α . In prin iple an ensemble of super systems an be prepared by applying the same history dependentma ros opi proto ol to all members of this ensemble. Due to the statisti al independen e of the subsystems, upontaking an ensemble average of super systems, we have h I α I β i S = h I α i S h I β i S for all α = β . Here the average h . . . i S istaken over the ensemble of super systems and the α th subsystem in ea h super system is identi(cid:28)ed by its lo ation. It isthen easy to show that the spe i(cid:28) heat obtained from the ensemble average Eq. 41 of the super system is equivalentto that obtained from the subsystem due to the two quantities (cid:10) I s (cid:11) S = (cid:28)(cid:16)P N s α =1 I α (cid:17) (cid:29) S and h I s i S = DP N s α =1 I α E S possessing identi al ross terms whi h an el ea h other out (as a result of the independen e of the subsystems) uponapplying Eq 41.If we ignore (cid:28)nite size e(cid:27)e ts, due to assuming the equivalen e of ensembles, the onstant volume spe i(cid:28) heat isrelated to the onstant pressure spe i(cid:28) heat by the equation C V = h V i c V = C P − P dVdT (cid:12)(cid:12)(cid:12)(cid:12) T − (cid:18) ∂H∂P (cid:12)(cid:12)(cid:12)(cid:12) T (cid:30) ∂V∂P (cid:12)(cid:12)(cid:12)(cid:12) T (cid:19) ∂V∂T (cid:12)(cid:12)(cid:12)(cid:12) P . (42)We may also obtain an expression for the onstant volume spe i(cid:28) heat c V by deriving equilibrium (cid:29)u tuation formulafor ea h of the derivatives appearing in Eq. 42 in lieu of dire tly measuring them. We may then obtain a single domainexpression for c V whi h does not work for the history dependent glass and also a orre tly weighted ensemble average(multidomain average) whi h does. This is ompletely analogous to what has been shown in detail for c P . As theequations are unwieldy, and their derivation (given an understanding of the c P ase) is straightforward, we will notreprodu e them here.1E. Test of Domain Robustness: Transient Flu tuation TheoremThe appli ation of the Evans-Searles Transient Flu tuation Theorem8,15,16 to the systems treated in this paperprovides a sharp test of the assumptions used to develop the theory given in this paper. The Theorem des ribesa time reversal symmetry satis(cid:28)ed by a generalized entropy produ tion, namely the so- alled dissipation fun tion.The pre ise mathemati al de(cid:28)nition of this fun tion requires a knowledge of the dynami s and also of the initialdistribution fun tion. The three ne essary and su(cid:30) ient onditions for the Flu tuation theorem to be valid are thatthe initial distribution is known (here we assume the distribution is Boltzmann weighted over some initial domain ofphase spa e), that the dynami s is time reversible (all the equations of motion used here are time reversible) and lastlythat the system satis(cid:28)es the ondition known as ergodi onsisten y. When applied to the systems studied here thisrequires that the phase spa e domains should be robust with respe t to the sudden hanges imposed on the systemand that the number of inter-domain transitions remain negligible on the time over whi h the theory is applied. Ifany one of these three onditions fails then the Theorem annot be applied to the system and the orresponding(cid:29)u tuation relation will not be satis(cid:28)ed8.We an use the (cid:29)u tuation theorem to obtain relations for how the system responds upon suddenly hanging theinput temperature or pressure for a system, whi h is initially in equilibrium as spe i(cid:28)ed by Eq. 7 or 8. Firstly we onsider a hange in the pressure, while holding the temperature (cid:28)xed, by hanging the input variable P in Eq. 4(thermodynami pressure) to P = P at time t = 0 for a system initially in equilibrium with P = P . The probabilitydensity p (∆ V = A ) of observing a hange in volume of ∆ V ( t ) = V ( t ) − V (0) relative to a hange of equal magnitudebut opposite sign is then given by p (∆ V ( t ) = A ) p (∆ V ( t ) = − A ) = exp ( β ( P − P ) A ) . (43)To derive this expression we have had to assume that the intra-domain populations are Boltzmann distributed a ord-ing to Eq. 7. Ergodi onsisten y requires that for any initial phase spa e point Γ (0) that an be initially observedwith nonzero probability, there is a nonzero probability of initially observing the time reversal map M T of the endpoint Γ( t ) , (ie ∀ Γ (0) su h that f ( Γ (0) , = 0 , f ( M T ( Γ ( t )) , = 0 ). This ondition obviously requires that the phasespa e domains remain robust and the number of inter-domain transitions remain negligible for at least a time t , afterthe pressure (or temperature) quen h.If we sample all or our initial t = 0 , P = P on(cid:28)gurations from the one traje tory whi h remains lo ked in asingle domain even after the quen h, we expe t Eq. 43 to be valid. If we prepare an ensemble of initial on(cid:28)gurationsusing the same proto ol we still expe t Eq. 43 to remain valid even with di(cid:27)erent domain weightings w i , as de(cid:28)nedin Eqs. 20 & 21, provided the domains are robust over the time t appearing in Eq. 43. Note that suddenly redu ingthe pressure by a very large amount ould result in a breakdown of the robustness ondition. Eq. 43 may be partiallysummed to obtain what is referred to as the integrated (cid:29)u tuation theorem p (∆ V ( t ) > p (∆ V ( t ) <
0) = h exp ( β ( P − P )∆ V ) i ∆ V < . (44)For the ase where we hange the input temperature T in Eq. 4 while holding the input pressure P (cid:28)xed we obtaina relation for (cid:29)u tuations in the extended instantaneous enthalpy I E ( t ) = H E ( t ) + P V ( t ) . We start with a systeminitially in equilibrium at temperature T = 1 / ( k B β ) and we then subje t it to a temperature quen h by hanging theinput temperature in Eq. 4 to T = 1 / ( k B β ) , at time t = 0 , while holding the input pressure (cid:28)xed. The probabilitydensity p (∆ I E ( t ) = A ) of observing a hange in instantaneous enthalpy ∆ I E ( t ) = I E ( t ) − I E (0) relative to a hangeof equal magnitude but opposite sign is then given by p (∆ I E ( t ) = A ) p (∆ I E ( t ) = − A ) = exp (( β − β ) A ) . (45)Note that if we suddenly in rease the temperature by a very large amount we ould expe t to violate the robustnessor the negligible inter-domain transition ondition. In ommon with Eq. 43 we expe t that this expression will bevalid when all initial on(cid:28)gurations are sampled from a single ommon domain and also when sampled from multiplearbitrarily populated domains under the assumption that the domains are robust and the number of transitions arenegligible over time t . This equation may also be partially summed to obtain p (∆ I E ( t ) > p (∆ I E ( t ) <
0) = h exp (( β − β )∆ I E ) i ∆ I E < . (46)2Figure 1: A logarithmi plot of the self di(cid:27)usion oe(cid:30) ient for both spe ies A and spe ies B parti les as a fun tion of theseparation parameter T − T g with T g = 0 . .IV. SIMULATION DETAILSFor our simulations we use a variation on the Kob Andersen glass former17 featuring a purely repulsive potential18.The pairwise additive potential is u ij ( r ij ) = 4 ǫ αβ (cid:20)(cid:16) σ αβ r ij (cid:17) − (cid:16) σ αβ r ij (cid:17) + (cid:21) ∀ r ij < √ σ αβ u ij ( r ij ) = 0 ∀ r ij > √ σ αβ , (47)where the spe ies identities of parti les i and j , either A or B , are denoted by the subs ripts α and β . The energyparameters are set ǫ BB = 0 . ǫ AA , ǫ AB = 1 . ǫ AA and the parti le intera tion distan es σ BB = 0 . σ AA , σ AB =0 . σ AA . The energy unit is ǫ AA , the length unit is σ AA and the time unit is p m σ AA /ǫ AA with both spe ies havingthe same mass m . The omposition is set at X = N B /N A = 0 . , the number of parti les are N = N A + N B = 108 ,the pressure is set to P = 14 ǫ AA /σ and the temperature unit is ǫ AA /k B . The time onstants are set at τ V = 5 √ N and τ T = √ N . Note that the energy parameters are slightly di(cid:27)erent to the potential we used in ref.18. The equationsof motion were integrated using a fourth order Runge-Kutta method19. The time step used was dt = 0 . andsometimes dt = 0 . for very low temperatures.From previous work on binary mixtures we know the basi reason why this system is vary relu tant to rystallize20,21,22,23. The hosen nonadditivity of the spe ies A-B intera tion makes the mixture extremely mis i-ble; onsider the present value of σ AB = 0 . σ AA relative to the additive value of σ AB = 0 . σ AA . This e(cid:27)e tdominates over the hoi e of the energy parameters. Due to this extreme mis ibility the relatively large omposition(cid:29)u tuations ne essary, about the average omposition of X = 0 . , to form the rystal phases (either the pure spe iesA, X = 0 , FCC rystal or the binary, X = 0 . , CsCl rystal) are strongly suppressed and rystallization is stronglyfrustrated.We identify the nominal glass transition by al ulating the di(cid:27)usion oe(cid:30) ient as a fun tion of temperature. Thisis shown for both spe ies in Fig. 1 on a logarithmi plot demonstrating how the di(cid:27)usion oe(cid:30) ient approa heszero riti ally, D α ( T ) ∝ ( T − T g ) b where b is the riti al exponent, with a nominal glass transition temperature of T g = 0 . . It would be, perhaps, more ustomary to obtain a nominal glass transition temperature by analyzingthe riti al divergen e of the vis osity. Given that the Stokes Einstein relation is strongly violated upon approa hingthe glass transition one might be on erned that this would give a very di(cid:27)erent result. However the violation ofthe Stokes Einstein relation an largely be attributed to the exponent b being di(cid:27)erent between the vis osity and thedi(cid:27)usion oe(cid:30) ient rather than the nominal glass transition temperature T g T = 1 . The al ulation of the fun tion was started after the (cid:29)uid had been given time to equilibrate. Thestrong agreement between the two orrelation fun tions is indi ative of ergodi ity.b) The instantaneous enthalpy orrelation fun tions as de(cid:28)ned by Eqs. 34 & 36 as a fun tion of logarithmi time for thetemperature T = 0 . . The al ulation of the orrelation fun tion was started at various times after the quen h, all approximatelyat t = 8 × , in an attempt to age the system. The di(cid:27)eren e between the orrelation fun tions is indi ative of nonergodi ity.Noti e that C rea hes a near full de ay between t = 1 & 10 , while C rea hes a non-de aying plateau.V. RESULTS AND DISCUSSIONThe orrelation fun tions given in Eqs. 34 & 36 were al ulated from ensembles of independent simulations atthe two temperatures given in Fig. 2 ( T = 1 and T = 0 . ). In all ases the systems were subje t to an instantaneousquen h, from an initial equilibrium at T = 5 , by hanging the value of the input temperature in Eq. 4. The systemwas then run for a signi(cid:28) ant time, in the ase of the glass ensemble τ age > × , in an attempt to age it. Of oursethe longest time that an be a essed in a mole ular dynami s simulation is rather short, and so the system is notvery well aged, but we are still able to meaningfully treat it as a time invariant state. Ea h of the independentsimulations was interpreted as being stu k in its own domain D α and the orrelation fun tions were al ulated forea h of these domains using time averaging. The time averaging was approximately times longer than the longesttime t = 800 that the orrelation fun tions were al ulated out to. Obviously in the limit of an in(cid:28)nite number ofindependent simulations and the ase where the domains are robust we will obtain the exa t multidomain averagegiven by Eq. 20. We assume our limited ensemble of simulations is representative of this. The data from ea h domain(independent simulation) was then used to obtain the orrelation fun tions Eqs.34 & 36 as seen in Fig. 2. At thehigher temperature T = 1 it an be seen that the two fun tions are equivalent demonstrating how the system isergodi . It an also be seen that the orrelation fun tion has de ayed on a time s ale of t ∼ whi h is therefore (byEq. 33) the time s ale on whi h the ensemble be omes a urately represented by Eq. 9 with only one domain D whi hdoes not ne essarily extend over all phase spa e. The os illations, whi h an be seen in the orrelation fun tion, are4due to both ringing in the feedba k me hanisms of Eq. 4 and the frequen y dependent storage omponent of the bulkvis osity. The statisti al un ertainty in the orrelation fun tion be omes larger than these os illations somewherebetween a time of t = 1 & 10 . If we onstru ted an experiment where the pressure was regulated by a piston and aspring, the orrelation fun tions, Eqs.34 & 36, would depend on the details of the piston and spring parameters in asimilar way to the simulations dependen e on the details of the feed ba k me hanism.When the system is quen hed to the lower temperatures ( T = 0 . ) ergodi ity is lost and we obtain a glass. The omplete de ay of the (cid:28)rst orrelation fun tion, C , Eq. 34 may no longer o ur be ause the individual traje toriesremain stu k in lo al domains whi h have di(cid:27)erent average values for h I i α . This is similar to the mu h studieddensity orrelation fun tion25 (the intermediate s attering fun tion) whi h de ays to a (cid:28)nite plateau for a glass ormore generally a solid material. On the other hand there is nothing to stop the se ond orrelation fun tion, C , Eq.36, from fully de aying when the system is nonergodi . If the se ond orrelation fun tion fully de ays whilst the (cid:28)rstis only able to de ay to a plateau then we have a situation where Eqs. 21 & 20 are valid as an be seen from Eq. 33. InFig. 2 it an be seen that C does indeed fail to de ay while C omes very lose to fully de aying at the time where C rea hes the plateau. The reason C doesn't fully de ay here is due to the fa t that the inter-domain transition rates,while small, are not exa tly zero. If the system had been aged more extensively this problem would be signi(cid:28) antlyredu ed. This e(cid:27)e t is exa erbated by the time averaging, used to form the averages for ea h traje tory, being twoorders of magnitude longer than the longest time the orrelation fun tion was al ulated to. The e(cid:27)e t of the stateslowly evolving due to (cid:28)nite inter-domain transition rates is too small to seriously ompromise the modeling of thesystem as obeying Eqs. 20 & 21 and thus we have obtained dire t eviden e for the validity of these equations. Theheight of the plateau for C will depend on the history of the system, i.e. the proto ol used to prepare the ensemble.We move on to a omparison between the kineti temperature and the on(cid:28)gurational temperature, the results ofwhi h may be seen in Fig. 3a). The input temperature ranges from T = 3 to T = 0 . . Also shown are results for thesystem, undergoing onstant planar shear, Eq. 10, with a strain rate of ˙ γ = 0 . . At the higher temperatures we seea very small relative dis repan y between the two types of temperatures, whi h we attribute to the dis ontinuity inthe (cid:28)rst spatial derivative of the inter-parti le for e at the uto(cid:27) radius and to (cid:28)nite size e(cid:27)e ts. These e(cid:27)e ts appearto diminish a little at lower temperatures. For temperatures above T = 1.5 the hosen strain rate has no signi(cid:28) ante(cid:27)e t on the on(cid:28)gurational temperature indi ating that our system is in the linear response domain18. At the lowesttemperatures, well below the glass transition temperature, we observe good agreement between the on(cid:28)gurationaland input temperatures for the system without shear. This provides further eviden e of our assertion that the systemobeys Boltzmann statisti s in the glass Eqs. 20 & 21. However, at low temperatures, the system that is undergoingshear shows an in reasing relative dis repan y between the two temperatures. At low temperatures the system leavesthe linear response domain18 demonstrating the fundamental di(cid:27)eren e between the nonequilibrium distribution ofthe history dependent glassy state and that of a strongly driven steady state. If we drive the system hard enough,at any given temperature, we an always make a disagreement between the two types of temperature due to thesteady state no longer being a urately represented by a Boltzmann distribution i.e. due to a break down in lo althermodynami equilibrium. When the system is not driven by an external (cid:28)eld we have been unable to observe anydi(cid:27)eren e in the two temperatures by deeply super ooling a glass forming mixture apart from the initial transientde ay immediately following the quen h, whi h falls o(cid:27) surprisingly rapidly.In Fig 3b) results are presented for the heat apa ities (the spe i(cid:28) heat multiplied by the volume) at both onstantpressure C P and onstant temperature C V . Details of the proto ol used to obtain this data is given in the endnote26. The estimates from the multidomain average are ompared with those from the single domain average. Theresults from the multidomain averages exhibit the well-known peak, whi h is a signature of the onset of the glasstransition, and has been observed dire tly by alorimetry in many experiments on real glass forming materials2.The temperature, where the peak is observed, depends on the history of the system. No peak is observed for thesingle domain averages whi h ontinue to in rease as the temperature is lowered. While the two methods for formingaverages give the same results at temperatures above the peak, they diverge at temperatures below the peak. It isthe multidomain average that gives results onsistent with the a tual alorimetri behavior of the system. This maybe seen in the (cid:28)gure by omparing the data whi h has been omputed by numeri ally di(cid:27)erentiating the enthalpyusing entral di(cid:27)eren e. At the peak neither the entral di(cid:27)eren e (due to rapid rate of hange) or the multidomainaverage data (due to a la k of domain robustness) are reliable and they show signi(cid:28) ant di(cid:27)eren es. However belowthe peak they on e again show quantitative agreement providing strong eviden e that the domains are robust in thisregion. If we substantially in rease the duration the time averages (for ea h domain) are onstru ted on, the peak willbe shifted to lower temperatures as previously shown27. This requires the time average to be onstru ted over sometwo orders of magnitude more time than the de ay time for the orrelation fun tion in Eq. 34. This is ne essary inorder to obtain enough independent samples for a meaningful estimate of the varian e, of the instantaneous enthalpyappearing in Eq. 39. For a large ma ros opi system we would expe t that the spe i(cid:28) heat measured over the entireensemble would di(cid:27)er very little to that measured from any one of its members. We are now in a position to makean unambiguous interpretation of the peak in the spe i(cid:28) heat. The peak is observed at the temperature where the5Figure 3: a) The relative di(cid:27)eren e between the kineti temperature T ( ontrolled dire tly by the Nosé Hoover thermostat)and the on(cid:28)gurational temperature given by Eq. 38 as a fun tion of the kineti temperature. The results for the systemundergoing Couette (cid:29)ow (Shear) are for a onstant strain rate of ˙ γ = 0 . .b) The heat apa ity al ulated using the equilibrium (cid:29)u tuation formula by the single domain averaging method Eq. 39 andthe multidomain method Eq. 41. Also shown is data obtained by numeri ally di(cid:27)erentiating the enthalpy by entral di(cid:27)eren e.At the temperatures above the peak the three types of averages give very similar results. Also shown is equivalent equilibrium(cid:29)u tuation formula data for the onstant volume spe i(cid:28) heat Eq. 42.system leaves metastable equilibrium and enters a history dependent state that requires averages to be omputed byEq. 20 rather than by dire t use of Eq. 9. The al ulation of both h I i α and (cid:10) I (cid:11) α will be di(cid:27)erent for ea h domain.If we use time averaging to al ulate these variables on a time s ale that falls within the plateau region for Eq. 34,see Fig. 2b, then the amount of time hosen to form the average is not riti al. The peak o urs be ause the variousensemble members have be ome lo ked in lo al domains on the time s ale that we are able to a ess. Near the peakitself these domains are not expe ted to be robust.The multidomain average, Eq. 41, gives the heat apa ity for a glass with robust domains. At the lowest tem-peratures the heat apa ity rea hes the beginning of a plateau, Fig. 3b. For the onstant volume heat apa ity C V this plateau (within un ertainties due to (cid:28)nite size e(cid:27)e ts) has a value that is onsistent with the Dulong-Petit law28,as would be expe ted for an amorphous solid where the potential energy surfa e an be modeled as harmoni upontransformation to the orthogonal independent basis set. This is exa tly what we would expe t from our lo al domainmodel at low temperatures.Testing the integrated transient (cid:29)u tuation theorem (ITFT) for a sudden pressure hange Eq. 44 and a suddentemperature hange Eq. 46 provides further eviden e that the Boltzmann distribution may be used to a uratelydes ribe intra-domain statisti s in the glassy state and also that in a properly aged glass, the domains are robustwith respe t to the pressure and temperature hanges studied here, Fig. 4. These equations remain valid whether wesubje t an ensemble of simulations (multidomain) to a quen h or we sample from a single traje tory (single domain),whi h remains stu k in a single domain. The a ura y with whi h these relations are satis(cid:28)ed is powerful independenteviden e for the appli ability of our assumptions to the systems studied here. The fa t that over the times shown inFig. 4, the ITFT does indeed yield orre t results dire tly implies that, within experimental toleran e of the data,the phase spa e domains must be robust and the number of inter-domain transitions must be negligible. Unlike the6Figure 4: Results from applying the (cid:29)u tuation theorem to the glass for a) a sudden pressure hange, where the symbols are p (∆ V > /p (∆ V < and the solid line is h exp( β ( P − P )∆ V ) i ∆ V < and b) a sudden temperature hange, where the symbolsare p (∆ I E > /p (∆ I E < and the solid line is h exp(( β − β )∆ I E ) i ∆ I E < . Results from an ensemble of independent initialsystems and in addition from a single initial traje tory (with a time of t = 5 being omputed between ea h transient traje tory)are shown for a total of pressure or temperature hanges. The serial results have been shifted up, for larity, by adding 0.2to the dataspe i(cid:28) heat (cid:29)u tuation formula this requires that the domains are robust to (cid:28)nite hanges of the state rather thanin(cid:28)nitesimal hanges. Thus, given that we have aged the glass su(cid:30) iently that domains are robust the number oftransitions are negligible over the longest time the (cid:29)u tuation formulae are omputed, we expe t Eqs. 44 & 46 to be orre t. If we wished to apply the steady state (cid:29)u tuation theorem matters be ome more di(cid:30) ult18.VI. CONCLUSIONSWe have presented a rigorous development of statisti al me hani s and thermodynami s for nonergodi systemswhere the ma ros opi properties are sensibly time independent and the phase spa e for the ensemble is partitionedinto robust domains. Using omputer simulation we have arried out various tests on a glassy system and shown thatapart from the immediate vi inity of the glass transition, the omputed results are onsistent with our theory. Whilethe intra-domain populations are individually Boltzmann distributed, the inter-domain populations are not.A orrelation fun tion whose de ay to zero requires global Boltzmann weighting, has been derived and it has beenshown that it de ays on a reasonable time s ale for ergodi systems but not for nonergodi systems. A se ond orrelation fun tion whi h de ays to zero if the intra-domain populations are Boltzmann distributed but globally theinter-domain populations are not, has also been de(cid:28)ned. We have developed expressions for obtaining averages in amultiple domain ensemble and shown how single domain averages, whi h always give orre t results in metastableequilibrium, an give spurious results in a history dependent nonergodi ensemble. The statisti al me hani s andthermodynami s developed here allow the derivation of expressions for multidomain ensemble averages whi h give the orre t results for time nondissipative nonequilibrium ensembles. The fundamental origin of the peak in the spe i(cid:28) heat near the glass transition has been unambiguously shown to be a signature of a transition from metastable7equilibrium to a nonergodi multi-domain ensemble.We have shown that the transient (cid:29)u tuation relations for temperature and pressure quen hes provide independenttests of the fundamental hypotheses used in our theory: that intra-domain populations are individually Boltzmanndistributed, that ex ept in the immediate vi inity of the glass transition the domains are robust with respe t to smallbut (cid:28)nite variations in thermodynami state variables, and that the inter-domain transition rates are negligible. ∗ Ele troni address: swilliamsrs .anu.edu.au † Ele troni address: evansrs .anu.edu.au1 D. A. M Quarrie, Statisti al Me hani s (Harper & Row, New York, 1976).2 P. G. Debenedetti, Metastable Liquids Con epts and Prin iples (Prin eton University Press, Prin eton, 1996).3 P. G. Debenedetti and F. H. Stillinger, Nature 410, 259 (2001).4 H. W. Sheng, H. Z. Liu, Y. Q. Cheng, J. Wen, P. L. Lee, W. K. Luo, S. D. Shastri, and E. Ma, Nat. Mater. 6, 192 (2007).5 J. P. Hansen and I. R. M Donald, Theory of Simple Liquids (A ademi Press, London, 1986), 2nd ed.6 J. R. Dorfman, An Introdu tion to Chaos in Nonequilibrium Statisti al Me hani s (Cambridge University Press, Cambridge,1999).7 D. J. Evans and G. P. Morriss, Statisti al Me hani s of Nonequilibrium Liquids. (A ademi , London, 1990).8 D. J. Evans and D. J. Searles, Adv. Phys. 51, 1529 (2002).9 F. H. Stillinger and T. A. Weber, Phys. Rev. A 25, 978 (1982).10 F. H. Stillinger and T. A. Weber, S ien e 225, 983 (1984).11 F. S iortino, J. Stat. Me h. p. P05015 (2005).12 C. Tsallis, R. S. Mendes, and A. R. Plastino, Physi a A 261, 534 (1998).13 W. G. Hoover, Time Reversibility, Computer Simulation, and Chaos (World S ienti(cid:28) , Singapore, 1999).14 W. van Megen, T. C. Mortensen, S. R. Williams, and J. Muller, Phys. Rev. E 58, 6073 (1998).15 D. J. Evans, E. G. D. Cohen, and G. P. Morriss, Phys. Rev. Lett. 71, 2401 (1993).16 D. J. Searles and D. J. Evans, Aust. J. Chem. 57, 1119 (2004).17 W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 (1995).18 S. R. Williams and D. J. Evans, Phys. Rev. Lett. 96, 015701 (2006).19 J. C. But her and G. Wanner, Appl. Numer. Math. 22, 113 (1996).20 S. R. Williams, I. K. Snook, and W. van Megen, Phys. Rev. E 64, 021506 (2001).21 J. R. Fernandez and P. Harrowell, Phys. Rev. E 67, 011403 (2003).22 J. R. Fernandez and P. Harrowell, J. Chem. Phys. 120, 9222 (2004).23 J. R. Fernandez and P. Harrowell, J. Phys. Chem. B 108, 6850 (2004).24 A. M. Puertas, E. Za arelli, and F. S iortino, J. Phys.-Condes. Matter 17, L271 (2005).25 W. Kob and H. C. Andersen, Phys. Rev. E 52, 4134 (1995).26 An ensemble of 100 simulations was used to obtain the spe i(cid:28) heat data starting from an initial equilibrium at T = 5 .The ensemble was sequentially quen hed to the next lower temperature at the ompletion of ea h y le. For ea h y le theensemble was equilibrated for a duration of t = 2 × before omputing a produ tion run of duration t = 1 × . Thetemperatures were T ==