# A simple thermodynamic description of the combined Einstein and elastic models

###### Abstract

Simple application of the Einstein model combined with the elastic description of solid state is developed. The frequency of quantum oscillators has been assumed as volume dependent and, furthermore, elastic energy terms of static character have been included to complete the description. Such an extension enables to construct the complete thermodynamics. In particular, the model yields practical equation of state and describes the thermal expansion coefficient as well as the isothermal compressibility of solids. The thermodynamic properties resulting from the Gibbs free-energy analysis have been calculated and illustrated in figures. Some comparison of the theoretical results with experimental data for solid argon has been made.

Keywords: Einstein model, equation of state, thermodynamic properties

## 1 Introduction

The Einstein model of the solid state has been known for many years as the first model able to describe qualitatively the low-temperature behaviour of the specific heat [2, 3, 4]. This model often serves as an approximation for studies of optical phonons, or soft modes in some intermetallic compounds [5] and, as well, is found suitable to describe thermal properties of some modern low-dimensional structures [6]. It is also referred to as a prototype for which the more sophisticated Debye model has been developed [7]. Nonetheless, it has been shown that in many cases the Einstein model provides better results than the Debye model [8]. It is known that in spite of its usefulness and simplicity, the pure Einstein model lacks the ability to describe fully the thermodynamic properties of the solid state. For instance, neither the proper equation of state, nor the thermal expansion or compressibility can be obtained. Thus, apart from the celebrated specific heat behaviour, which is an important consequence of the quantum nature of harmonic oscillators, other thermodynamic properties have not successfully been described within this model.

There are numerous attempts in the literature aimed at improving the pure Einstein model [9, 10, 11, 12, 13].
For instance, a ”variational Einstein model” has been developed in Ref. [9] for describing low temperature solids from the Feynman path integral perspective. The formalism has been applied to a specific system, consisting of solid hydrogen with lithium impurities. However, one of the theoretical results of Ref. [9] seems to be controversial, namely the dependence of the free energy, which is presented there as an increasing function of temperature (see Fig.16 of Ref. [9]), giving the indication that the corresponding entropy is negative.
Cankurtaran and Askerov [10] have introduced the so called Einstein-Debye model for which the calculations of the thermal expansion coefficient were possible by introducing the Grüneisen parameter.
However, neither the Gibbs free-energy, nor the compressibility could be calculated within that approach. Thus, the complete thermodynamic description has not been achieved there.
In Ref. [11] the ”nonlocal Einstein crystal” has been considered, for which the thermodynamic functions have been constructed. The considerations presented in Ref. [11] seem to be of pure theoretical character when the Einstein crystal is replaced by a harmonic sublattice with a size that is coupled to the total volume.
Additionally, the Mie-Grüneisen approach has been developed by Holzapfel et al. for the description of simple metals [12] as well as some simple solids, such as NaCl and MgO [13]. The optimized pseudo-Debye-Einstein model has been used there, and good agreement with the experimental data has been achieved. However, the method presented in Refs. [12, 13] does not seem to be straightforward enough to use since it requires many fitting parameters to be introduced for calculations of the free-energy.

On the other hand, there exists a rich literature on the description of solid state properties within the elastic models [14, 15, 16, 17, 18, 19, 20, 21, 22, 23].
These studies are supported by the experimental measurements [24, 25].
For instance, Birch [20] considers the single-crystal and polycrystalline NaCl at high pressure using the finite strain (Eulerian) theory. In this approach, the equation of state contains coefficients which have to be fitted in each particular temperature. Despite the usefulness of the model for studies of isothermal compression, it does not involve the thermal (vibrational) energy and therefore is not suitable for description of the specific heat. Similarly, in the paper of Vinet et al. [21, 22] the universal form of the equation of state has been found based on the scaling argumentation. The results obtained there have been discussed in the context of the Birch-Murnaghan equation of state [19] and show very good agrement with the experimental data. However, the Gibbs free-energy has not been derived in that works (Ref. [21, 22]), so that the full thermodynamic description of the system has not been constructed. As a matter of fact, the same deficiency can be risen for some other papers devoted to the application of the equation of state based on the elastic models. A comparative review of the experimental equations of state existing in the literature is given in Refs. [26, 27].

Among other modern methods one should mention the first principles computations within the density-functional theory (DFT). This includes, for instance, the equation of state, elastic constants and phonon dispersion relation calculations [28]. On the other hand, by the first-principles molecular-dynamics simulation (FPMD) method the thermal expansivity and the specific heat have been calculated [29]. A good agreement with the experimental data is often achieved. However, the analytic description remains always desirable from the point of view of understanding the physics behind the model.

Such a situation motivated us to combine the Einstein and elastic models and to complete them in a way that enables to overcome the above-mentioned problems. Looking in the literature, such idea can already be found in some papers [30, 31], however, only particular thermodynamic characteristics have been studied there. Thus, our aim is to construct the self-consistent and full thermodynamic description. In order to make it possible, the following aspects should be taken into account: Firstly, the common frequency of quantum oscillators should not be treated as a constant value but should be volume (and, in consequence, temperature) dependent. Obviously, such idea is not new and for the first time emerged in the paper of Grüneisen [32], where a specific frequency/volume dependence has been assumed as a hypothesis. This assumption is supported by the argument that the wavelength of collective excitations (phonons) depends on the crystal size.
Secondly, the vibrational Einstein Hamiltonian necessarily needs to be completed by the static, elastic part [33]. The elastic energy accounts for the mutual interactions of the atoms even if they are not in a thermal vibrating state [19, 23].
It is known that the static energy is responsible for crystal compressibility [19], and owing to its magnitude and nature, cannot be inferred from the model of independent oscillators. It turns out that for the purpose of the combined model the elastic free-energy should also be modified by taking into account the linear term (which is normally neglected in the elastic theory) as well as the static entropy. The details of this modification will be explained in the next section. Owing to analytical simplicity we will neglect the electronic excitation term. It is argued that this term is generally a small correction to the solid equation of state [31].

Taking this into account, we propose improving the model and constructing the Gibbs energy of the system from the beginning. We assume that the balance between the internal and external pressures keeps the system in mechanical equilibrium.
According to our knowledge, the detailed balance between the expanding pressure of quantum oscillators and the compressive elastic pressure has not been discussed in the literature.
Having obtained the expression for the Gibbs energy, a full thermodynamic description of the crystal can be achieved.
The main parameters of the theory, which can be extracted from experiment, are the Einstein frequency in the ground state (or the Einstein temperature), the volume elastic modulus in the ground state - supplemented by the structural space-filling coefficient - and the Grüneisen parameter. We will show that these parameters, in principle, allow us to calculate other properties within a reasonable range of experimental values.

The paper is organized as follows: In the following, theoretical section, the outline of the model is given and the method of derivation the Gibbs energy is presented. In particular, for the combined model the new equation of state is obtained. In the last section, some representative numerical results concerning the thermodynamic properties are illustrated in figures. These calculations are compared with the experimental data for solid argon. A critical discussion of the presented approach is also included.

## 2 Theoretical model

### 2.1 The model Hamiltonian

The Hamiltonian is assumed in the form of:

(1) |

where the elastic, volume dependent part, , can be written as:

(2) |

whereas the oscillatory part, , is given in the form:

(3) |

In (2) and (3) is the number of atoms which are treated as the three-dimensional oscillators, and is the excitation number operator connected with the -th oscillator. The form of the elastic Hamiltonian (2) assumed in this paper is different (much simpler) from that presented in Refs. [30] or [31]. The elastic Hamiltonian, , consists of several terms. The most important is the harmonic term , where -constant is the volume elastic modulus in the ground state. We found that the linear term is also necessary, where -constant is responsible for the internal (compressive) pressure. Although the linear, anharmonic part is assumed to be very small in comparison to the harmonic one (), its role is very important in balancing the internal (expanding) pressure produced by the anharmonic oscillators, i.e., by the part. A certain requirement for -values is imposed as the equilibrium condition for the whole energy, which is discussed later. The relative elastic deformation is defined by the relation:

(4) |

where is the crystal volume in a current thermodynamic state, and is the volume at external pressure (vacuum) and temperature . Equation (4) indicates that a non-zero value of volume deformation occurs when the external pressure and/or the temperature is applied. Thus, a possibility of introducing the isothermal compressibility and thermal expansion coefficient has been opened within the presented approach.

A composition of the Hamiltonian (1) consisting of the classical () and quantum () parts is justified for the purpose of constructing the thermodynamics [23]. Obtaining the free energy of the full system is one of the main goals and that can be conveniently done by summing up the free energies corresponding to those two parts. A similar modus operandi has also been presented in the paper Ref. [33] albeit for a different approach. It should be stressed that by introducing the elastic Hamiltonian, , the static non-vibrational energy has been taken into account, which is crucial for the equation of state. It is worth mentioning that the method can be further generalized by introducing also the electronic part of the free energy [23].

### 2.2 The elastic free-energy

In order to calculate the elastic (Helmholtz) free energy, one has to evaluate the static internal energy and the static entropy. The static internal energy can be immediately found from the classical part of the Hamiltonian:

(5) |

From Eq.(2) it is seen that this energy is explicitly volume-dependent. On the other hand, the static entropy is connected with filling up the volume of the system in the absence of thermal movement. Let us denote by the own volume occupied by atoms in the crystal state. For a given crystallographic lattice we can introduce the space filling coefficient . For instance, for FCC and HCP structures and the model of hard spheres . Thus, inside each primitive cell is the probability of finding an atom, whereas is the probability of opposite event. The static entropy connected with the bimodal distribution (where ) can be calculated from the mean value of , using the general formula for the exact differential of entropy: [34]

(6) |

Hence, in our case:

(7) |

where is a constant of integration. Since is arbitrary, as a matter of convenience we can choose . This arbitrariness is due to the fact that only changes in entropy are experimentally measurable and have thermodynamic significance [34, 35].

The entropy is additive with respect to the number of primitive cells . For the purpose of this paper it is assumed in the first approximation that the space-filling coefficient is volume independent and, in consequence, the corresponding entropy is constant. Thus, the static entropy (for ) results in the residual entropy of the crystal for , which is due to the fact that the whole volume is not perfectly filled. Moreover, the important role of
(being dependent on )
may manifest itself in the structural (1st order) phase transitions [23].

The elastic, Helmholtz free-energy can now be found from the formula:

(8) |

Hence, we finally obtain:

(9) | |||||

For the full thermodynamic description this energy should be completed by the vibrational, Einstein part.

### 2.3 The vibrational free-energy

For the Hamiltonian (3) it is convenient to employ the canonical ensemble. The Helmholtz free-energy can be found from the formula:

(10) |

where the statistical sum can be calculated exactly as:

(11) |

where . Thus we obtain:

(12) |

On the other hand, the vibrational internal energy is given by the statistical mean value of the corresponding Hamiltonian:

(13) |

Hence, the vibrational entropy can be found from the relationship:

(14) | |||||

The specific heat at constant is given by the well-known formula:

(15) | |||||

where .

Further, in order to improve the Einstein model it is assumed that is not a constant but volume dependent according to the Grüneisen assumption [32]:

(16) |

where is the Grüneisen constant. Thus, we can write:

(17) |

where we made use of the relation (4), and in (17) is the Einstein frequency in the ground state (i.e., at , , and ). Then, the variable can be presented in the form of:

(18) |

where

(19) |

is the so-called Einstein temperature. The correction of the Einstein model presented by Eq.(18) implies that all the thermodynamic potentials will depend on the volume via elastic deformation .

### 2.4 Thermodynamic functions of the combined model

The total Helmholtz free-energy is given by the sum of expressions (9) and (12), with the help of (16)-(19):

(20) | |||||

The total Gibbs energy is then given by:

(21) |

The equation of state can be obtained from the variational principle:

(22) |

which, from (21), is equivalent to:

(23) |

where is given by (20). From (23), after differentiation of (20), we obtain:

(24) | |||||

For the high-temperature limit equation (24) takes a form:

(25) |

which is the classical lattice dynamics pressure [23]. On the other hand, in the low-temperature limit () and for small pressure the deformation is small too, and Eq.(24) can be linearized as:

(26) |

Naturally, when then in the ground state is as well. Consequently, we have to demand that the internal pressure of oscillators, which then is equal to , cancels out the pressure arising from the anharmonic, linear term, i.e., . Then the system is in equilibrium without any external force. The condition for this demand leads to the specific choice for the anharmonic constant , namely:

(27) |

where

(28) |

The energy constant presents a convenient unit for introduction of the dimensionless external pressure defined by:

(29) |

and the dimensionless temperature :

(30) |

With the above quantities, and Eqs.(27) - (28), the equation of state (24) takes the simple form:

(31) | |||||

Eq.(31), obtained within the approximations assumed in this paper, presents a dimensionless equation of state from which the elastic deformation can be calculated for arbitrary temperature and external pressure .

The Gibbs free-energy (21) in the dimensionless form can be written as follows:

(32) | |||||

where is given by (31). It is easy to show that from the formula (32) all thermodynamic properties can be derived self-consistently. For instance, the total entropy fulfills the relationship:

(33) |

and, with the use of (31), it can be obtained in the dimensionless form:

(34) | |||||

On the other hand, the volume of the sample satisfies the equation:

(35) |

which again leads to the equation of state (31). Other thermodynamic properties follow from the second derivatives of the Gibbs energy. For instance, the heat capacity at constant is given by:

(36) |

whereas the heat capacity at constant can be calculated from the relationship:

(37) |

On the other hand, the internal energy, , can be written in the dimensionless form:

(38) | |||||

Application of (38) in (37) leads to the dimensionless expression for :

(39) |

It should be noted that for the particular case, when , Eq.(39) is equivalent to formula
(15). However, in general, a non-zero should be determined from the equation of state (31).

Among other response functions which can be calculated either from the Gibbs free energy, or from the equation of state is, for example, the thermal expansion coefficient:

(40) |

and the isothermal compressibility:

(41) |

In particular, for and , i.e., in the ground state, these response functions can be obtained from the approximate equation of state (26). For this purpose, Eq.(26) can be re-written in the form:

(42) |

It can be observed that Eq.(42) for presents Hooke’s law at . Further, we can formally make use of the linear expansion of volume, namely:

(43) |

It can be noticed that Eq.(43) can equivalently be written as:

(44) |

where and are the response functions (40) and (41), respectively, and these functions are taken in the ground state. Now, by comparison of Eqs.(42) and (44) we obtain the result:

(45) |

and

(46) |

It is apparent that the both elastic coefficients, and , are important for determination of isothermal compressibility at . Since , where is the zero-point energy of Einstein oscillators, it is obvious that for this compressibility an elastic energy plays a dominant role. Finally, let us note that the equation of state derived by Vinet et al. [22] for reduces for small to the form: , which agrees with our equation(42) in the limit .

### 2.5 A note on the Grüneisen equation

The Grüneisen parameter is introduced by the definition [3, 4]:

(47) |

which is equivalent to our relation (16). The typical experimental values of the Grüneisen parameter belong to the range [4]. It can be mentioned that some particular definitions of that parameter, depending on the model, have been discussed in Ref. [36]. The experimental determination of is connected with the so-called Grüneisen equation:

(48) |

which has been tested extensively, especially in the low temperature region. It is worth noticing that the derivation of (48) can be based on the exact thermodynamic identity:

(49) |

In the case in hand, in order to calculate , the equation of state should be taken in the form of (24). It can easily be checked that after the calculation of the derivative in (49), eq.(48) is satisfied. It should also be noted that the calculations of all quantities in (48) using the present method require prior solution of the equation of state (31). In particular, for the low-temperature region, where tends to zero according to the 3rd law of thermodynamics, we obtain from eq.(48) that . This limit is in agreement with the previous result (45).

## 3 Numerical results and discussion

### 3.1 Exemplary calculations for the model

In order to perform the numerical calculations, based on the formulas from the preceding section, it is necessary to estimate the energy constants, from which and are the most important. The volume elastic modulus in the ground state, , can be found, for instance, from the sound velocity measurements and its experimental value is of the order () J. In turn, -coefficient is given by Eq.(28), where is the Einstein temperature. Assuming that is typically of the order K we can estimate as J. Thus, the -constant is about 2 orders smaller than the elastic bulk modulus , and the linear anharmonic term in the Hamiltonian (2) can be considered as a small correction to the harmonic one. This confirms the fact that such quantities as the thermal expansion coefficient (Eq. 40) or the isothermal compressibility (Eq. 41) are mainly determined by -constant, not by . The above estimations of and constants, together with the assumed , yield from Eq.(40) 1/K for the temperatures near the Einstein temperature, which is a physically reasonable order of value. In the same token assuming a realistic amount of volume per atom, i.e., , we obtain on the basis of Eq.(41) 1/Pa, which is also a correct order for the isothermal compressibility.

Although the anharmonic elastic energy is a small correction in comparison to the harmonic one, it can be of the same order of magnitude as the vibrational energy resulting from Einstein oscillators. Therefore, this anharmonic term is vital in balancing the expanding pressure resulting from oscillators, in the case when the frequency is volume dependent (17). The expanding pressure arises owing to the decrease of the energy (frequency) of oscillators when they experience collective excitations in increasing volume. The total equilibrium of the system requires balancing of the internal (stretching) forces resulting from expanding quantum oscillators, internal compressive forces of the elastic medium and the external pressure . The resulting deformation is temperature and pressure dependent and can be calculated from the equation of state (31).
The other energy constants, and , play merely a correcting role in modeling the static potential, and are important for high temperatures where the elastic deformation is significant.
The exemplary numerical results, presented in this section in Figs.1-8, are obtained for the following set of parameters: , and .

In Fig.1 we present the dependence of upon for different reduced temperatures 0, 0.5 and 1, as well as for two selected Grüneisen parameters: and . The calculations are based on the equation of state (31) for . It is seen that these dependencies are almost linear in character, and the slopes of the curves correspond to the isothermal compressibility. It can be deduced from Fig.1 that the isothermal compressibility coefficient should be weakly dependent on pressure in a wide range of temperatures. This observation is in accordance with the analogous assumption in the Murnaghan theory of the equation of state [19]. At higher temperatures (for ) the Grüneisen parameter has remarkable influence on the magnitude of elastic deformation, leading to the increase of when increases.

In Fig.2 the dependence of upon for various external pressures ( 0, 0.5 and 1) is presented. The rest of parameters are the same as in Fig.1. In this case the dependencies are not linear and their local slopes are attributed to the thermal expansion coefficient. It is remarkable that by increasing the Grüneisen parameter the relative deformation increases, which is evidently pronounced at high temperatures. Such behaviour is in agreement with the previous figure (Fig.1). The dependence presented in this figure is qualitatively similar to the lattice constant dependence on temperature, as can be seen in the Ref. [37] for the case of AlN. The effect of external pressure on the temperature dependence of relative deformation in Fig.2 resembles, for example, the results obtained in the Ref. [27].

In Fig.3 the dimensionless thermal expansion coefficient (see Eq.40) is plotted vs. temperature for external pressure and Grüneisen parameter . Different curves correspond to various parameters. It is seen that has influence mainly at high temperatures, where it causes the increase of the thermal expansion coefficient. Let us note, on the basis of Eq. (28) and (30), that corresponds to the Einstein temperature. In the low-temperature region the thermal expansion coefficient tends to zero independently on , which is in agreement with the Grüneisen equation (48). A qualitatively similar behaviour for the thermal expansivity has been found in Ref. [13] for the case of MnO, in Ref. [19] for Ti, Al, NaCl and Na, or in Ref. [37] for AlN. The conclusion that the anharmonic term becomes important only in the high-temperature region is in agreement with the free energy calculations from first principles [29]. The similar conclusion can be drawn from the paper Ref. [38] where the influence of intrinsic anharmonicity on the thermodynamic functions has been studied.

The dimensionless isothermal compressibility (see Eq.41) is plotted in Fig.4 vs. for the external pressure and Grüneisen parameter . As in Fig.3, different curves correspond to various parameters. A remarkable influence of -values on the compressibility is seen for the temperatures . On the other hand, for the compressibility tends to some non-zero value, , being in agreement with Eq.(46). A qualitatively similar experimental results have been obtained for Pb and NaCl. [19] The dependence of isothermal compressibility on temperature corresponds to the fact that the elastic constants are temperature dependent (the exemplary studies can be found in Refs. [33, 39]).

In Fig.5 the Gibbs free-energy per 1 atom is presented in -units vs. reduced temperature . The three curves in Fig.5 correspond to different external pressures: 0, 0.5, and 1, and are plotted by the solid, dashed and dashed-dotted lines, respectively. In this case and . We see that application of the external pressure causes the increase of the Gibbs energy. It is demonstrated in Fig.5 that the Gibbs energy is a concave function of temperature with a monotonously decrease vs. . Such behaviour indicates thermodynamically stable solution when the entropy (defined by Eq.33) is positive. Let us note that the initial slope of the Gibbs energy, at corresponds to the residual entropy.

The entropy is illustrated in Figs.6 and 7 as a function of temperature . In both figures . In Fig.6 we plotted entropy per 1 lattice site in the Boltzmann constant units, for external pressure . The solid and dashed curves, which are for and , respectively, correspond to our model, where is a function of temperature and is calculated from the equation of state (31). Due to our choice of integration constant ( in Eq.(7)) we see that the residual (configurational) entropy is present, which does not depend on the Grüneisen parameter . Moreover, one can notice that for the entropy change vs. temperature tend to zero (), which is in accordance with the third law of thermodynamics. On the other hand, the increase of causes some small increase of entropy at high temperatures. For comparison, the dashed-dotted curve presents the entropy for the pure Einstein model, where is put equal to 0 for all temperatures. As stated before, the condition is equivalent to the assumption that the frequency of oscillators is kept constant and does not depend on volume. For the pure Einstein model there is no residual entropy in the ground state, since the configurational ordering of atoms is not taken into account.

It turned out that the entropy is very weakly sensitive to the external pressure, and in Fig.7 we want to study this property in more detail. The relative changes of the entropy are plotted vs. for and 1. The dashed and solid lines are for and , respectively. We note that the sensitivity of on is greater for higher Grüneisen parameter . It is also seen that the influence of high external pressure on the ratio is most remarkable for some intermediate temperatures , whereas for the influence of becomes negligible. Again, the behaviour of near is in agreement with the third law of thermodynamics. The lowering of entropy for , when the external pressure is applied, corresponds to the positive thermal expansion coefficient, in accordance with the Maxwell relation: .

The specific heat per 1 lattice site in the Boltzmann constant units and at constant volume is plotted vs. temperature in Fig.8. The anharmonic parameters are . The main plot has been obtained for and . We note that the value of corresponds to the Einstein temperature. The low-temperature behaviour of the specific heat is typical for the Einstein oscillators, being in agreement with the 3rd law of thermodynamics. On the other hand, the high-temperature part of the curve is classical. In order to see the difference between the present model and the classical Einstein result in more detail, in the inset, we plot this difference for and , as well as for two Grüneisen parameters: (dashed line) and (solid). corresponds to the specific heat of the present model, whereas is the Einstein result, when does not depend on . It is worth noticing that for some small increase of the specific heat occurs for temperatures . Similarly to the entropy, the specific heat is only weakly sensitive to the external pressure . In particular, one can see that the pressure of causes a small decrease of at some restricted range of low temperatures. When the Grüneisen parameter increases all the above changes become enhanced.

### 3.2 A comparison with experimental data for solid argon

For solid argon, which forms FCC structure below the melting temperature of 84 K, the Debye temperature is =85 K. From the approximate formula , which can be derived on the basis of Ref. [3], the Einstein temperature, =68.51 K, can be estimated. Hence, on the basis of eq.(28) we find that the -constant amounts to J. On the other hand, the isothermal compressibility at zero temperature, , can be found from the experiment [40], and for polycrystalline samples it amounts to . The volume per atom can be estimated as . The Grüneisen parameter can be taken as a mean value of the experimental data from various temperature ranges, [40] which yields . Having obtained the data given above, on the basis of eq.(46) the ratio can be estimated, yielding . The other elastic energy constants, which are of higher order, i.e., and , can be treated as the theoretical fitting parameters. We found that the optimal fit to the experimental data (taken from table III, page 561 of Ref. [40]) is obtained for =1250 and =5000, valid for . Taking into account the above set of parameters, we have calculated and for solid argon in the full temperature range K. To complete the comparison, we supplemented the results with two example isotherms. Since the isotherms cover the range of both positive and negative deformation , we decided to allow some asymmetry in our model elastic potential. Thus, for we adopted the anharmonic parameter value of =3000, i.e. the potential is more steep in this range. The results are presented in Figs.9-13 by the solid lines, whereas the experimental data [40] are shown by the open symbols. In addition, the figures have been supplemented with the smoothed experimental data from the Refs. [41, 42, 43].

In Fig.9 the calculation of elastic deformation vs. absolute temperature is presented.
The derivative of elastic deformation over the pressure, i.e., the isothermal compressibility coefficient, , is plotted vs. in Fig.10. It can be noted that in both figures (9 and 10) the experimental dependencies are nonlinear, in agreement with the present theory.
In turn, the thermal expansion coefficient, , at constant pressure () is plotted in Fig.11 for the same temperature range as in Figs. 9 and 10.
Finally, the molar specific heat at constant is presented in Fig.12. For comparison, in Fig.12 the dashed line is plotted for the pure Einstein model when we impose the constraint .

Our calculated isothermal compressibility, , can be related to the adiabatic compressibility, , via the formula:

(50) |

It can be checked that the above formula is equivalent to the Grüneisen equation (48) which is also satisfied in our case. As pointed out in Ref. [40], the values of and can be measured independently; by the ultrasonic method from isentropic sound velocity, and by the piston-displacement method of Bridgman. Using the formula (50) a satisfactory agreement between those two measurements has been obtained in Ref. [40]. Such consistency directly transfers to our calculations, therefore in Fig.10 only curve has been presented.

Regarding the isothermal compressibility (at , ) we adopted here the value for polycrystalline samples after Dobbs and Jones [40]. On this basis the classical (atomic volume) bulk modulus, , can be estimated at K: . On the other hand, for single cubic crystals the bulk modulus can be obtained from the formula [3] , if the elastic constants and are known. These constants for argon single crystals have been deduced from measured isentropic sound velocity [44]. The values extrapolated to K are: and in units of [44]. Hence, the bulk modulus obtained on this basis is . It can be concluded that this value of for single crystal is about larger than for polycrystalline samples given in [40].

In the Fig. 13 two selected isotherms are shown (low- and high-temperature one), presenting elastic deformation as a function of external pressure . It is visible that using a model asymmetric potential, a good fit to experimental data in whole pressure range studied is obtained, both for high- and low-temperature results. This emphasizes the pronounced importance of the exact form of the static lattice potential for pressure studies, since the deformations here exceed 10%. The selected potential provides reliable description for pressures up to 10 kbar. In light of the assumption that our static elastic potential originates from the expansion at the point and , further adjusting the potential would improve the consistency between calculations and experimental data for higher pressures.

It can be concluded, on the basis of Figs. 9-13, that the agreement between calculated lines and the experimental data is quite satisfactory. One should take into account that all these curves have been calculated selfconsistently for the same set of parameters, as described above in this subsection. It can be supposed that the fitting could even be better if we would allow the Grüneisen parameter to vary with temperature and pressure, as it has been suggested from the experimental measurements [40, 45]. As far as the specific heat is concerned, the new result (solid line) differs only insignificantly from the Einstein result (dashed), and the difference is mainly noteworthy at high temperatures. However, it should be noted that the dashed line in Fig.12 is the only result of the pure Einstein model which can be compared with all these experimental data. The full description of temperature dependencies of remaining quantities, such as and has been possible in the improved approach to the Einstein model, when the equation of state (31) is taken into account.

### 3.3 Final remarks

In the present paper a simple combination of the Einstein and elastic models of solid state is presented. First of all, the dependency of frequency of quantum oscillators on the volume has been introduced in a simplified way via Grüneisen assumption. Simultaneously, the elastic properties of the crystal have been taken into account via classical Hamiltonian, containing anharmonic terms.

An idea that the free-energy is a sum of several components has been presented, for instance, in the book of Wallace [23]. However, we have shown that exploitation of this idea led us to the new form of the equation of state (31). Contrary to the Birch-Murnaghan equation of state, which presents only isothermal description [23], in our equation the temperature plays a role equivalent to the rest of variables, i.e., and .
In particular, the thermal expansion coefficient, as well as the isothermal compressibility have simultaneously been derived from this new equation of state. It is well-known that these quantities could not be inferred from the sole Einstein model. On the other hand, our equation of state includes also the pressure resulting from expanding quantum oscillators. Moreover, comparing our method with the papers based on the Wallace approach, one should notice that one of our anharmonic parameters in the static potential, namely , is not independent, but has been related to the Einstein temperature (via Eqs. 29 and 30). As we pointed out, such relationship assures that the equilibrium condition for the total free energy at and is satisfied and the system of quantum oscillators becomes stable. When this point is not discussed, the functional form of the static lattice energy can also be adopted in the form given in Ref. [31] (after Vinet et al.). That form is much more complicated than our polynomial approach, nevertheless it has successfully been applied for very high pressures.

The present method requires several constant parameters, such as: the Einstein temperature, isothermal compressibility (or the volume elastic modulus) at the absolute zero temperature, as well as the Grüneisen parameter, which should be taken from experiment.

The region of applicability of the model (its temperature and pressure range of validity) is mainly connected with the number of terms which are taken into account in the polynomial form of elastic Hamiltonian (2). This form is based on the assumption that the equilibrium central point is for and . Thus, the theory is best applicable for the expansion around this equilibrium. However, as we have seen on the example of solid argon, merely first four terms were sufficient to describe deformation in a wide range; starting from (near the melting temperature and , Fig.9) down to (for K and kbar, Fig.13). Of course, validity of the model for lower values of (for very high pressure, where anharmonicity of static potential plays a role) would require higher order terms in (2). One should also remember that our assumption concerning the space filling coefficient () would not be fulfilled at extremely high pressures. As far as the vibrational anharmonic effects are considered, they all are taken into account by the effective Grüneisen parameter . We think that the approach is useful in the range of temperatures up to the melting point.
In conclusion, the presented formulation allows for a complete thermodynamic description of the system in a relatively wide range of external pressure and temperature.
Our considerations are related to the quasistatic processes. Therefore, the shock-wave experiments cannot be described within this model.

It should be noted that the method is relatively simple, gives analytical form of the equation of state, and therefore it can serve as a first approximation for more advanced approaches.
For instance, in the prototype calculations we have used only single variable for description of the volume elastic deformation. However, it seems possible to generalize the approach for anisotropic deformations and anisotropic external pressures, involving also the Poisson coefficient.
The presented method can be potentially extended basing on the Debye approximation, in which the linear dispersion relation for collective excitations together with the proper density of states are taken into consideration. For instance, in Ref. [31] the vibrational free energy has been taken into account in the high temperature Debye model. In that approach each moment of the density of states function requires a separate Grüneisen parameter, which makes the method much more complicated. On the other hand, for metallic systems the electronic part of the free-energy should also be included [31]. However, such extensions of the method need further studies and should be a subject of separate assignment.

The paper has been partly supported by the grant VEGA 1/0431/10.

## References

- [1]
- [2] A. Einstein, Annalen der Physik 22, 180 (1907).
- [3] C. Kittel, Introduction to Solid State Physics, (John Wiley & Sons, Inc., 1996).
- [4] L. A. Girifalco, Statistical Mechanics of Solids, (Oxford University Press, 2000).
- [5] A. D. Caplin, G. Grüner, J. B. Dunlop, Phys. Rev. Lett. 30, 1138 (1973).
- [6] I. Avramov, M. Michailov, J. Phys.: Condens. Matter 20, 295224 (2008).
- [7] P. Debye, Annalen der Physik 39, 789 (1912).
- [8] V. I. Zubov, Cryst. Res. Technol. 30, 149 (1995).
- [9] D. Li, G. A. Voth, J. Chem. Phys. 96, 5340 (1992).
- [10] M. Cankurtaran, B. M. Askerov, Phys. Stat. Sol. (b) 194, 499 (1996).
- [11] M. T. Lusk, J. Chem. Phys. 121, 11208 (2004); Phys. Rev. B 70, 174103 (2004).
- [12] W. B. Holzapfel, M. Hartwig, W. Sievers, J. Phys. Chem. Ref. Data 30, 515 (2001).
- [13] U. Ponkratz, W. B. Holzapfel, J. Phys.: Condens. Matter 16, S963 (2004).
- [14] F. D. Murnaghan, Amer. J. Math. 59, 235 (1937).
- [15] F. Birch, J. Appl. Phys. 9, 279 (1938).
- [16] J. C. Slater, Introduction to Chemical Physics, (McGraw-Hill Book Co., 1939).
- [17] P. W. Bridgman, The Physics of High Pressure, (G. Bell and Sons Ltd., 1949).
- [18] G. R. Barsch, J. Appl. Phys. 39, 3780 (1968).
- [19] D. C. Wallace, Thermodynamics of Crystals, (John Wiley & Sons, Inc., 1972).
- [20] F. Birch, J. Geophys. Res. 83, 1257 (1978).
- [21] P. Vinet, J. R. Smith, J. Ferrante, J. H. Rose, Phys. Rev. B 35, 1945 (1987).
- [22] P. Vinet, J. H. Rose, J. Ferrante, J. R. Smith, J. Phys.: Condens. Matter 1, 1941 (1989).
- [23] D. C. Wallace, Statistical Physics of Crystals and Liquids, (World Scientific, 2002).
- [24] M. S. Anderson, C. A. Swenson, Phys. Rev. B 10, 5184 (1974).
- [25] M. S. Anderson, C. A. Swenson, Phys. Rev. B 28, 5395 (1983).
- [26] P. B. Roy, S. B. Roy, J. Phys.: Condens. Matter 17, 6193 (2005).
- [27] J. Hama, K. Suito, J. Phys.: Condens. Matter 8, 67 (1996).
- [28] M. Foata-Prestavoine, G. Robert, M.-H. Nadal, S. Bernard, Phys. Rev. B 76, 104104 (2007).
- [29] Z. Wu, Phys. Rev. B 81, 172301 (2010).
- [30] E. D. Chisolm, S. D. Crockett, D. C. Wallace, Phys. Rev. B 68, 104103 (2003).
- [31] C. W. Greeff, M. J. Graf, Phys. Rev. B 69, 054107 (2004).
- [32] E. Grüneisen, Annalen der Physik 39, 257 (1912); Handbuch der Physik, (Verlag Julius Springer, Berlin, 1926).
- [33] Tai-Hyung Kwon, Journal of the Korean Physical Society 26, 349 (1993).
- [34] T. L. Hill, Statistical Mechanics, (McGraw-Hill Book Co., 1956).
- [35] F. Reif, Fundamentals of Statistical and Thermal Physics, (McGraw-Hill Book Co., 1985).
- [36] J. J. Gilvarry, Phys. Rev. 102, 331 (1956).
- [37] H. Kröncke, S. Figge, B. M. Epelbaum, D. Hommel, Acta Phys. Polon. A 114, 1193 (2008).
- [38] A. R. Oganov, P. I. Dorogokupets, J. Phys.: Condens. Matter 16, 1351 (2004).
- [39] S.-C. Kim, T. H. Kwon, Phys. Rev. B 45, 2105 (1992).
- [40] E. K. Dobbs, G. O. Jones, Rep. Prog. Phys. 20, 516 (1957).
- [41] O. G. Peterson, D. N. Batchelder, R. O. Simmons, Phys. Rev. 150, 703 (1966).
- [42] M. S. Anderson, C. A. Swenson, J. Phys. Chem. Solids 36, 145 (1975).
- [43] C. R. Tilford, C. A. Swenson, Phys. Rev. B 5, 719 (1972).
- [44] G. J. Keeler, D. N. Batchelder, J. Phys. C: Solid State Phys. 3, 510 (1970).
- [45] M. Ross, H. K. Mao, P. M. Bell, J. A. Xu, J. Chem. Phys. 85, 1028 (1986).