{'fraction_non_alphanumeric': 0.05236980690462258, 'fraction_numerical': 0.0155061439438268, 'mean_word_length': 3.958575112830432, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anharmonic Phonon Quasiparticle Theory of Zero-point and Thermal Shifts in Insulators: Heat Capacity, Bulk Modulus, and Thermal Expansion

The Quasi-harmonic (QH) approximation uses harmonic vibrational frequencies omega(H,Q,V), computed at volumes V near the volume where the Born-Oppenheimer (BO) energy is minimum. When this is used in the harmonic free energy, QH approximation gives a good zeroth order theory of thermal expansion, and first order theory of bulk modulus. Here, n-th order means smaller than the leading term by n powers of epsilon, where epsilon is the ratio hbar omega(Q)/E(el) or kT/E(el), and E(el) is an electronic energy scale, typically 2 to 10 eV. Experiment often shows evidence for next order corrections. When such corrections are needed, anharmonic interactions must be included. The most accessible measure of anhamonicity is the quasiparticle (QP) energy, omega(Q,V,T), seen experimentally by vibrational spectroscopy. However, this cannot just be inserted into the harmonic free energy F(H). In this paper, a free energy formula is found which corrects the double-counting of anharmonic interactions that is made when F is approximated by F(H,omega(Q,V,T)). The term"QP thermodynamics"is used for this way of treating anharmonicity. It enables (n+1)-order corrections, if QH theory is accurate to order n. This procedure is used to give corrections to specific heat and volume thermal expansion. The QH formulas for isothermal and adiabatic bulk moduli are clarified, and the route to higher order corrections is indicated.


I. INTRODUCTION
In non-magnetic, insulating materials, thermodynamic behavior is controlled by vibrational excitations, which are often close to harmonic. There is a unique and correct version of harmonic theory, based on Taylor-expanding the Born-Oppenheimer (BO) energy E el ({ R }) around the atomic coordinates { R ,0 } of a crystal with volume V . The BO energy is a "ground state" property, and the main target of density functional theory (DFT). This gives normal mode eigenvectors and frequencies ω Q,H (V ), where the index Q = { Q, j} labels the states; Q runs over the N wave-vectors of the Brillouin zone, and j runs over the 3n branches. These states can be called "non-interacting quasiparticles." However, in this paper, the term quasiparticle (QP) will be reserved for the vibrational resonances seen experimentally. These differ from harmonic eigenfrequencies because higher order ("anharmonic") terms in the Taylor expansion are not negligible. There is a curious contrast with electrons, where quasiparticles are also seen experimentally, but there is not a single unique definition of a non-interacting quasiparticle limit.
In this paper, the term "harmonic" will refer to the unique correct harmonic limit, further specialized to the case when the volume is chosen to be V 0 , where the BO energy is minimum. It is useful also to know the harmonic normal modes (and their frequencies ω Q,H (V )) at other volumes; this is the "quasi-harmonic" (QH) theory. The harmonic approximation is a good starting point, successfully implemented by "ab initio" DFT calculations 1 , and useful, often to good approximation, for things like specific heat, C P (T ). Vibrational spectroscopy 2 of reasonably pure crystals most often sees reasonably sharp Lorentzian resonances. They can be assigned a wave-vector Q, and are expected to show a one-to-one correspondence with the harmonic normal modes. They are the QP's of this paper. The central frequency ω Q (the QP frequency, or energy) is T -dependent. There is good evidence from theory-experiment comparisons 3 that the QH energy ω Q (V ), evaluated at the correct thermally-expanded volume V (T ) at higher T , does not reproduce well the measured thermal shifts ∆ω Q of QP energies ω Q (V, T ). There is an anharmonic thermal shift, additional to and different from, the pure volume-related shift of QH theory, and it is often significant at higher T .
The QP relaxation rate 1/τ Q is the full width at half maximum of the spectroscopic Lorentzian line. In pure crystals it lies outside harmonic theory, and is also Tdependent. This paper is about extracting additional thermodynamic information from the temperature dependence of ω Q , ignoring 1/τ Q . This suffices for most low-order thermal corrections. I will call this "quasiparticle thermodynamics". It differs from QH theory. Deviations from harmonic vibrations are responsible for thermal shifts of many properties. The ones of prime concern in this paper are C P and C V (constant pressure and constant volume), bulk modulus B T and B S (isothermal and adiabatic), and volume expansion V (T ) and α = (1/V )dV /dT (constant pressure.) Good general discussions are in older literature. [4][5][6][7] There are two main ideas in QP theory: (1) lowlying excitations correspond 1-to-1 with a non-interacting single-particle picture; they have QP energies ω Q (V, T ) and mode occupancies n Q ; and (2) low energy dynamics can be described as the dynamics, in space and time, of the mode occupancy. QP theory can fail in at least two ways. (a) The resonance can be very non-Lorentzian so that ω Q is poorly defined. (b) Correlated occupancy n QnQ − n Q n Q may become important. In this paper, n Q denotes the equilibrium (Bose-Einstein) mean occupancy [exp( ω Q (V, T )/k B T ) − 1] −1 , n Q,H its harmonic or QH version, and n Q denotes the actual occupancy in a general ensemble, not necessarily equilibrium. Entropy plays a special role, since it is just the logarithm of the number of ways of distributing mode occupancies 8 , The equilibrium occupancy n Q is the one which maximizes Eq.(1) at fixed energy, and the thermodynamic entropy S(T ) is given by Eq.(1) with n Q → n Q . When harmonic frequencies ω Q,H are used in n Q , the result is the harmonic entropy S H . When the T -dependent QP energies are used in n Q , and inserted in S(T ), an accurate improvement of the thermodynamics is achieved. This will be denoted S QP . The same replacement does not work for the harmonic free energy. If ω Q (T ) is inserted into F H , the resulting F does not obey −dF/dT = S QP . A "correct" QP free energy that does agree with S QP is found as follows. The thermodynamic energy U (T ) = F + T S in harmonic theory is When normal mode frequencies acquire an anharmonic correction, ω Q,H → ω Q (T ), the energy acquires a correction, 9 This corrects for double-counting of the interaction, but only in leading anharmonic approximation where 7 (4) Further details are in the Appendix. Here ∂ω Q /∂n Q is a T -independent anharmonic interaction function. Then an accurate and consistent thermodynamic free energy is F QP = U QP − T S QP . Notice that the anharmonic shift, Eq.(4), does not vanish at T = 0, but has a zero-point component where n Q +1/2 → 1/2. The quasiparticle frequencies are shifted from the harmonic frequencies even at T = 0. QH thermodynamics is a limiting case. It uses only volume-dependence of ω Q,H (V ). The correction factor in U QP , Eq.(3), vanishes, and the quasi-harmonic free energy is just the harmonic free energy with a volumedependent harmonic frequency. It improves pure harmonic theory and gives correct lowest-order thermal shifts for properties such as the bulk modulus which are volume derivatives of F . The reason it works to lowest order is because dω Q /dV differs from dω Q,H /dV only in next order. QH theory is computationally accessible, but QP theory much less so. QP theory suffers from the dif-ficulty that the anharmonic shift is not usually measured except at a few temperatures. It can be numerically computed using DFT for the anharmonic forces. It is not yet computed routinely, but this is possibly changing. 10,11 Exact theory associates vibrational resonances with poles of a phonon Green's function, a correlation function describing dynamics on the BO energy surface. Exact thermodynamics should be computed from a corresponding theory for the free energy. This can be computed perturbatively 6,12 . At high T , classical molecular dynamics (MD) describes dynamics non-perturbatively, if the BO forces are known. This is called ab initio MD, or AIMD. Thermodynamics generally then requires a tricky "thermodynamic integration" [13][14][15][16] . To do a correct nonperturbative computation at lower T requires quantum corrections, as in path-integral MD. 17 Zero-point and related thermal nuclear motions cause shifts and isotope-dependences in many measured physical properties [4][5][6][7][18][19][20][21][22]44 . Explicit formulas are given here for the first-order shifts of C P , B T , B S , and α. If no specification (adiabatic versus isothermal) is made, isothermal is implied. The adiabatic shift can be found by thermodynamic rules, as shown in Sec. VI. Many of these results can be found in some form in the literature. There is a lot of correct, 4-7 plus much partially correct, as well as incorrect or confusing literature on this subject. There are semi-empirical formulas that have a long history of enabling useful fitting, even though the formulas do not seem to be justifiable in detail. 23,24 The aim of this paper is a simplified, possibly less confusing, derivation of correct formulas.
The paper is organized as follows. In Sec. II, examples of thermal shifts from experiments are given. In Sec. III, extra complexities of non-cubic crystals, and crystals with internal coordinates, are discussed. In Sec. IV, the QH approximation and the QP approximation are discussed. Specific heat formulas are presented, showing that QP theory provides a thermal correction inaccessible in QH approximation. In Sec. V, two orders of thermal correction to the volume are discussed. This gives Grüneisen theory of thermal expansion α 0 plus a first order correction. In Sec. VI, the leading correction to the bulk modulus is derived (from QH theory). The Appendix reviews the microscopic theory of Eqs.(1-4).

II. EXPERIMENTAL EXAMPLES
Figures 1,2,3 illustrate the thermal shifts under discussion. Fig. 1 shows that the bulk modulus has surprisingly large vibrational corrections. These have serious significance for geoscience, for example. 18,19 The leading-order bulk modulus, B 0 , comes from electronic stiffness. The product B 0 V a (V a = V 0 /N n is the volume per atom in leading order theory) has order of magnitude 10 eV, a characteristic electron energy. Vibrational energies are two orders of magnitude smaller. I will define as the dimensionless ratio of phonon to electron energies. This will appear explicitly in the form ω/BV a or k B T /BV a in various results. A parameter like ≈ 0.01 controls the size of the vibrational corrections under discussion. Fig. 1 shows ≈ 10% shifts, indicating that there can be a significant prefactor multiplying . The small parameter can be estimated as The "O-H stretch" vibrations at the opposite end of the spectrum have ω/k B larger by 11. These modes also contribute to the zero-point shifts in ice. If optic modes are used to define , the value of this "small parameter" increases to 0.5.
It is perhaps worth mentioning that the representation of a physical property P as an expansion in (P = P 0 +P 1 +P 2 2 +...) is not forced to have universallydefined coefficients. Especially because the expansion is truncated after the P 1 or possibly P 2 term, it is normal that the last coefficient may (or may not, depending on the source) contain some higher effects (P 1 = P 10 + P 11 , for example.) The only rule is that P n n should contain nothing of lower order than n . This will be mentioned again in Secs. V and VI. The volume shifts shown in Fig. 2 are smaller in relative size. Fig. 3 shows that α, the temperature derivative of the volume, roughly follows a harmonic specific heat (C H ) type of T -dependence. This is the result of Grüneisen theory 43 . But at higher T , there is a very significant thermal shift of both C P (T ) and α(T ) away from the C H (T ) form.
Figs. 1 and 2 also illustrate zero-point shifts. Mean square thermal lattice displacements of the i th atom, in harmonic theory, are 4 where i|Q is the component of the phonon eigenvector |Q on the i th atom, and where n Q is the equilibrium occupation number. The zero-temperature value u 2 ∼ /2M ω is the quantum zero-point motion, which depends on nuclear mass, whereas the high-T value k B T /M ω 2 is classical and depends only on the force constant M ω 2 , not on the nuclear mass M . The low T u 2 causes zero-point shifts of atomic volume V (T = 0) and bulk modulus, which differ for different isotopes. Therefore, the DFT ("frozen-lattice") value V 0 or B 0 should differ from the actual value V (0) or B(0). It is interesting that the frozen-lattice value can sometimes be deduced from experiment. 22 This is because the thermal factor n + 1/2 at high T becomes x(1−1/12x 2 +· · ·) where x = k B T / ω. An asymptotic linear-in-T fit to n + 1/2 at high T passes through 0 at T = 0. It is of course difficult to find the "correct" experimental asymptote, since thermal corrections enter to alter it. However, the bulk modulus simplifies the fit if both isothermal and adiabatic versions are available, because each should extrapolate to the same T = 0 value B 0 , as is shown on Fig. 1 We have already seen that the anharmonic shift ∆ω Q = ω Q − ω Q,H of phonon frequencies has a similar form, Eq.(4). The T -independent coupling parameter ∂ω Q /∂n Q has contributions like V 4 and |V 3 | 2 / ω; see the Appendix, Eq.(A3). They have order-of-magnitude ω Q . For example, the term of the type |V 3 | 2 / ω involves a third order anharmonic coupling coefficient V 3 , of structure u 3 ∂ 3 E el /∂u 3 , and of order (u/a) 3 E el . The ratio (u/a) 2 of lattice displacement to interatomic distance is of order ω/M ω 2 a 2 ≈ . Putting all these factors together, we see that The smallness of the anharmonic shift is only a crude estimate which sometimes may fail. A failure is likely to cause anharmonic broadening Γ Q of vibrations to be bigger than the spacing of vibrational levels ω Q . In such cases, phonon quasiparticles are poorly defined, perturbative treatments may fail, and thermal shifts may not be described well by quasiparticle theory.
The validity of perturbative computation beyond harmonic approximation for thermodynamic properties is not a closed issue. Wallace 44 summarizes evidence for failure of anharmonic perturbation theory to reproduce apparently reliable MD. But Boltzmann equa-tion treatments of thermal conductivity are now very successful [45][46][47] , and are based on the third-order term in the same perturbation theory. Computations based on DFT anharmonic forces are becoming more common, and generally claim decent agreement with experiment. A nice example is theory and experiment for aluminum by Tang et al. 48 . A thermal conductivity κ > 5 W/mK is a good hint that phonon quasiparticles are mostly well behaved. This crude estimate comes from κ ∼ Cv /3 where the specific heat is C = 3k B /Ω a , Ω a is the volume per atom, v ∼ πω max /a is the sound velocity, a the lattice constant, and the phonon mean free path. Quasiparticle theory requires large compared to a, or κ large compared to κ min = Cva/3. If we choose ω max /k B ∼ Θ D to be 300K, and a to be 5Å, then κ min ∼ 1W/mK, and κ > 5κ min should be sufficient to trust quasiparticle theory for most of the phonons of the material. However, Θ D may be significantly bigger or smaller than 300K, and the criterion could be scaled to κ > (Θ D /300K)× 5 W/mK.

III. NON-CUBIC CRYSTALS AND INTERNAL COORDINATES
Pressure, volume, and temperature are not the only thermodynamic variables in crystals. One can also have anisotropic stress σ αβ and anisotropic strains αβ . Pressure and volume change are the traces of these tensors. This paper looks only at pressure and volume. The generalization to tensor properties complicates notations and results, but the principles are not changed. Consider hexagonal structures as a simple example of non-cubic. The separate a(T ) and c(T ) lattice parameters are relevant thermodynamic variables. They are also not considered in this paper, only V (T ) = √ 3a 2 c/2 is considered. When T changes, not only does V change, but also c/a. This can not be ignored, but is kept hidden in this paper. The volume-dependent phonon frequency ω Q0 (V ) is treated as a well-defined variable. There is a hidden assumption that this has been computed at various volumes, and for each volume, the correct c/a ratio has been found and used in the phonon calculation. Finally, consider wurtzite crystal structure (hexagonal symmetry and 4 atoms per cell). There is one "internal coordinate" u(T ), which determines the c-axis offset between the cation and anion sublattices. This also cannot be ignored. But it is hidden by the assumption that for a particular choice of V , the correct u(T ) as well as c/a have been computed, and used to find ω Q,H = ω Q (V 0 ), and ∂ω Q,H /∂V , etc. Cubic crystals can also have internal coordinates not fixed by symmetry, which need to be treated the same way.

IV. QUASIPARTICLE THERMODYNAMICS
Even when harmonic approximation is seriously perturbed by anharmonic effects, there may still be phonon quasiparticles, with effective interactions not too strong. Thermodynamics is then approximated by using the QP frequencies ω Q (V, T ) in the non-interacting entropy formula 8 , Because of the V and T -dependence of the QP energy, S QP has altered V and T -derivatives which give corrections in thermodynamic calculations 7,20,49-51 . The corresponding free energy is where ∆ Q = ω Q (V, T ) − ω Q,H was defined in Eq.(4). This version of F QP is the same as U QP − T S QP and Eq.(3) for U QP . The quasi-harmonic formulas S QH and F QH are the same except that ω Q (V, T ) is replaced by ω Q (V ), usually calculated by DFT. In that case, the last term of Eq. (7) vanishes. In a metal or a magnetic material, one should include additional terms in F QP for thermal excitation of electrons or magnons. Such effects are omitted here. The QH procedure of using just a volumedependent QP energy in the harmonic free-energy formula, does give correct first-order V -derivatives, but fails to give thermal shifts which depend on T -derivatives. In this sense, it can be regarded as an incomplete, rather than an incorrect, theory, and as a partially correct simplification of QP theory. It correctly contains the information available from DFT calculations of the frequencyspectra at different volumes. Ramirez et al. 64 made a careful study of the accuracy of the QH approximation by comparison with well-converged path-integral MD for three phases of ice. They find generally very good agreement between QH and PIMD As an example of QP thermodynamics, consider the specific heat, C = T dS/dT . The free energy is not needed; the correct QP entropy is Eq.(6) with QP frequencies in the equilibrium occupation functions.
Here X is pressure P or volume V . This gives where the subscript "H" means (dn Q /dT ) H = ( ω Q /k B T 2 )n Q (n Q + 1), obtained by differentiating n Q by the explicit T in the Bose function, but not by the implicit T contained in ω Q (V, T ). The first term of Eq. (9) is a harmonic specific heat C H , but not the purely harmonic C 0 , because the frequencies ω Q appearing in the formula are the renormalized T -dependent quasiparticle frequencies.
The difference between C H and C 0 is a gentle T -dependent stretching of the harmonic C 0 (T ) curve along the T axis. This does not affect the high T classical limit, 3N nk B . A serious high-T deviation from harmonic theory (see the measurement for Ar, Fig. 3) must be caused by the second term of Eq.(9). In QH theory, (∂ω Q /∂T ) V = 0, so C V,QH = C H . Also in QH theory, (∂ω Q /∂T ) P = (∂ω Q /dV ) T (∂V /∂T ) P , so there is a significant QH correction to C P . QH theory gives the correct difference, C P − C V , but it misses the anharmonic correction which appears in the correct QP theory for both C V and C P . Computational evidence 50 shows that QH theory shifts C P (T ) away from the harmonic value C 0 , but that experiment exhibits different shifts. 52

V. VOLUME EXPANSION
Grüneisen's papers of 1912 and 1918 43 were a remarkable advance, simultaneous with the first true understanding of crystals that came from Rutherford and Bohr, von Laue and the Braggs, Born and von Karman, Eucken, and Debye. In modern language, the leading-order (Grüneisen) relation for the thermal shift of the volume The Grüneisen parameter is Since γ Q is ∼ O(1), the volume shift ζ is ∼ O( ), and Eq.(10) is, in fact, correct to this order. The parameter ω Q /B 0 V 0 appears explicitly. It was mentioned in Sec. II that expressions like Eq.(10), nominally first-order in the small parameter , may contain higher order parts. This is a good example. Available Grüneisen parameters are usually γ Q,H , from DFT quasi-harmonic frequencies.
The anharmonic shift ∆ Q = ω Q − ω Q,H is of order ω Q , and the volume derivative d∆ Q /dV is ∼ O( ω Q /V a ). The relation between the QP γ Q and the QH γ Q,H is Evidently γ Q − γ Q,H is of order , and Eq.(10) contains, in addition to the correct order answer, order 2 corrections if experimental ω Q and γ Q are used. But there are additional order 2 corrections needed for a complete theory to order 2 . To derive Eq.(10) and in addition the full order 2 corrections, start from the free energy, Eq.(7). Thermodynamics (P = −(∂F/∂V ) T ) gives Here P el is the frozen lattice pressure −∂E el /∂V . Geophysicists and others like to define an average Grüneisen parameter γ and to write Eq.(13) as P = P el + γU vib /V where U vib is the harmonic vibrational energy, Eq.(2). This is called the "Grüneisen equation of state." It ignores the last term of Eq. (14), and is correct to order if a careful T -dependent average γ is defined.
The first term of P vib in Eq. (14) is ∼ O( ) relative to P el , and the second term, which is missing in previous theory, is ∼ O( 2 ). In fact, it is appropriate to drop the part of the second term involving ∂n Q /∂V , because this brings in an additional power of . There is no point in keeping terms of order 3 , because at that level, anharmonic theory has higher-order parts; ∆ Q in higher order contains pieces like (ω Q − ∆ Q /2) (n Q + 1/2). (16) Here ∂ω Q /∂V is a simpler but less conventional way to write −ω Q γ Q /V . The combination (ω Q − ∆ Q /2) can equally well be written (ω Q,H +∆ Q /2) or (ω Q,H +ω Q )/2. For results to order 2 , it is necessary to expand E el to third order in ζ, is the order 0 electronic contribution to the bulk modulus, andḂ 0 is the third derivative, V 2 0 (d 3 E el /dV 3 ) V0 . The notationḂ is used because the notation B means dB/dP = −1 −Ḃ/B in equation-of-state theory. 18,19 Then we get which can be combined with Eq.(16) and inserted in Eq. (13) to give an equation-of-state to one higher order of accuracy in powers of than the usual theory. Now set the pressure to zero. To order 1 we get This gives the usual Grüneisen theory of thermal expansion, (20) Evidently α 0 behaves likeγC(T )/B 0 V 0 where C(T ) is the specific heat, as illustrated in Fig. 3. It is of order /T , and saturates to a constant (T -independent) value To get ζ correct to order 2 , set the pressure to zero, but keep both terms in Eq. (18) for P el and in Eq. (16) for P vib . The resulting quadratic equation for the small quantity (V − V 0 )/V 0 is then solved iteratively to order 2 , giving Experiment at high T , such as the results for argon shown in Fig 3, prove that higher order thermal corrections to the volume expansion α are significant. The QH answer for these higher shifts has been computed in the past, 33,52 but is not enough. There is a term missing in the nextorder QH correction that is now available from Eq. (21). To order 2 /T , where ζ 0 is Eq. (19), and its temperature derivative contains only one term, proportional to dn Q /dT . The vibrational part is A QH version of Eqs. (22,23) uses ω Q,H in place of ω Q in all factors, and uses the Grüneisen (QH) ζ 0 , Eq. (19) in ω Q,H (V (T )). This contains some, but not all, of the order 2 terms which increase linearly with T at high T (as does experiment

VI. THERMAL CORRECTION TO BULK MODULUS
The literature about B(T ) is large because of its importance in geoscience. The bulk modulus is the simplest and most accessible part of the elastic constant matrix C ij , all components of which show related zeropoint and thermal alteration. This paper focuses on B for simplicity, but generalization to the full elasticity tensor is not hard. 5,7,32,33,52,[54][55][56][57][58][59][60] The leading-order shift of the isothermal bulk modulus is found correctly from QH theory. The simulations by Ramirez et al. 64 (using a fluctuation formula for B(T ) 8 ) confirm this for phases of ice.
Here the perturbative answer in terms of the volume derivatives of E el and ω Q are derived. The answer involves a higher-order Grüneisen parameter, Volume-dependence of γ Q (V ) has often been neglected. If the mode Grüneisen parameter were independent of volume, one could integrate to find As observed previously 61,62 , there is no justification for this except unwarranted optimism. The notationγ Q used here (Eq.24) is not the same as V (∂γ Q /∂V ) = γ Q + γ 2 Q +γ Q . For cases like ice, where volume shifts are relatively large, it is insufficient to compute only low-order derivatives of energy and vibrational frequency. But QH theory computed separately for different volumes along the QH V (T )-curve, has been shown to work. 64 The quasi-harmonic bulk modulus is derived from the usual definition, B T = −V (∂P/∂V ) T for the case where quasi-harmonic theory, the first term of Eq. (14), is used for P vib .
The last term uses the identity ∂(ω Q γ Q /V )/∂V = ω QγQ /V 2 . The next to last term uses the fact that V ∂n Q /∂V equals T (∂n Q /∂T ) H γ Q . From Eq.(18) we find In Eq. (25) it is implicit that all derivatives are evaluated at the volume V that minimizes F QH (V, T ) at finite T , rather than at the frozen lattice volume V 0 (which minimizes E el (V )). The first term is ∼ O(B 0 ) and all other terms are ∼ O(B 0 ) or smaller. The aim to use param-eters evaluated at V 0 while retaining accuracy ∼ O( 1 ). All terms except the first can be evaluated at V 0 without altering the answer except to order 2 . The term B el shown explicitly in Eq. (26), to first order in , is B 0 + (B 0 +Ḃ 0 )ζ 0 , where ζ 0 is given in Eq. (19). Then Eq.(25) becomes The parameter is not truly small for ice Ih. For this reason, Eq. 27 does not work particularly well 63,64 . Direct computation and minimization of the QH free energy (Eq. 7 without the last term) may work. This has done used for many years, even in cases where the shifts are small enough that Eq.(27) should be adequate. 32,33,55,57,65,66 In cases, like ice Ih, where is too large to use Eq. (27), there is no guarantee that truly anharmonic terms of order 2 and higher are not as important as QH terms found by direct minimization.
It is important to distinguish between adiabatic (B S ) and isothermal (B T ) conditions 54 . The definitions are where U and F are the internal energy and Helmholz free energy respectively. Thermodynamics gives exact relations 5,7,18,19,32,54,67,68 , where C V /V is the heat capacity per volume. The product αT is of order , and C V T /B T V is also of order , so the shift (B S − B T )/B is positive and order . The full tensor version also is available 7,69,70 . The vibrational corrections δ S (adiabatic) and δ T (isothermal) are both first order in , and they differ from each other in the same order. The leading order value of T α 2 B T V /C V is sufficient for correcting isothermal to adiabatic. Using Eq. (20) and the harmonic specific heat, the result is Figure 1 shows approximate high T slopes (dB/dT ) of both B S and B T for NaCl. In the high T limit where ω Q (dn Q /dT ) → k B , Eq.(31) reduces to d(B S − B T )/dT = 3k Bγ 2 /V a , whereγ = Q γ Q /3N . The slopes shown in Fig. 1 then requireγ ≈ 1.5, in good agreement with other estimates ofγ for NaCl.

Appendix A
This appendix tries to illuminate the quasiparticle thermodynamics of Eqs.(1-4) by using anharmonic thermal perturbation theory. According to Cowley 6 , the vibrational thermal Green's function matrix, in the basis of harmonic eigenstates |λ with frequencies ω λ,H , is (G −1 ) λλ = (ω 2 λ,H − ω 2 )δ λλ + 2(ω λ ω λ ) 1/2 [∆ λλ − iΓ λλ ] (A1) The eigenvalues ω 2 of the matrixĜ −1 + ω 21 , with the imaginary partΓ omitted, are denoted ω 2 λ . They are the renormalized squared normal mode frequencies. If anharmonicity is weak, then in leading approximation these eigenvalues are ω 2 λ = ω 2 λ,H + 2ω λ ∆ λλ . At the same level of approximation, ω λ = ω λ,H + ∆ λλ . Cowley gives an explicit formula from lowest-order perturbation theory, for the anharmonic shift ∆ λλ = ∆ Q = ω Q − ω Q,H . The normal mode index λ is now replaced by Q = ( Q, j). Cowley's formula can be written in the form This is an explicit form for Eq.(4). Here V (3) and V (4) are third and fourth derivatives of the BO potential taken around the periodic sites of the crystal of volume V 0 , and the frequencies ω Q and the occupation number n Q use anharmonic renormalization (computed selfconsistently.) Cowley also derives the anharmonic free energy at the same level of perturbation theory. His answer can be written as where F H,0 is the free energy of non-interacting (harmonic) quasiparticles. Now find the corresponding entropy, S = −dF/dT . The non-interacting part gives the harmonic entropy, Consider what happens if the "quasiparticle entropy" is constructed by replacing the harmonic frequencies in Eq.(A6) by the anharmonic frequencies ω Q,H +∆ Q . Taylor expanding to first order in ∆ Q , the answer is This is the same as the entropy dF/dT from Eq.(A4). The factor 1/2 in A4 disappears when using Eq.(A2) while differentiating Eq.(A4) for ∆F , because ∂ω Q /∂n Q is symmetric in Q and Q . An alternate derivation using a variational principle is given in ref. 50. This suggests that the use of QP energies in the harmonic entropy formula may be valid somewhat beyond low-order perturbation theory. Consider then what happens if the same substitution ω Q,H → ω Q is done in the harmonic free energy F H = Q ω Q,H (n Q + 1/2) − T S H . The answer is ∆F H = Q ∆ Q (n Q +1/2)−T ∆S. This differs from the correct anharmonic free energy, Eq.(A4), by not having the correct factor of 1/2. This is a proof of the doublecounting correction that was added to the internal energy in Eq.(3). The correct formula Eq.(A4) does differ from the QP theory of Eqs.(1-4) by a small T -independent term F A0 .