A self-consistent thermodynamic model of metallic systems.Application for the description of gold

A self-consistent thermodynamic model of metallic systems. Application for the description of gold

T. Balcerzak    K. Szałowski Department of Solid State Physics, Faculty of Physics and Applied Informatics,
University of Łódź, ulica Pomorska 149/153, 90-236 Łódź, Poland
M. Jaščur Department of Theoretical Physics and Astrophysics, Faculty of Science,
P. J. Šafárik University, Park Angelinum 9, 041 54 Košice, Slovak Republic
July 4, 2019
Abstract

A self-consistent thermodynamic model of metallic system is presented. The expression for the Gibbs energy is derived, which incorporates elastic (static) energy, vibrational energy within the Debye model, and electronic part in Hartee-Fock approximation. The elastic energy is introduced by a volume-dependent anharmonic potential. From the Gibbs energy all thermodynamic quantities, as well as the equation of state, are self-consistently obtained. The model is applied for the description of bulk gold in temperature range K and external pressure up to 30 GPa. The calculated thermodynamic properties are illustrated in figures and show satisfactory agreement with experimental data. The advantages and opportunities for further development of the method are discussed.

Gibbs energy, thermodynamic properties, equation of state, gold
pacs:
64.10.+h; 65.40.+w; 71.10.Ca; 63.10.+a

I Introduction

The properties of noble metals and, in particular, of gold have intensively been studied, both experimentally Neighbours et al. (1958); Daniels et al. (1958); Zimmerman at al. (1962); Isaacs (1965); Martin (1966); Cohen et al. (1966); Skelskey et al. (1970); McLean et al. (1972); White et al. (1972); Lynn et al. (1973); Barin et al. (1973b); Touloukian et al. (1975b); Heinz et al. (1984); Anderson et al. (1989); Collard et al. (1991); Takemura et al. (2008); Yokoo et al. (2009); Kusaba et al. (2009) and theoretically Hiki et al. (1967); Gupta (1968); Rosén et al. (1983); Sundqvist et al. (1985); Shim et al. (2002); Saxena (2004); Baria et al. (2003); Baria (2004); Greeff et al. (2004); Godwal et al. (1989); Moriarty (1995); Boettger (2006); Souvatzis et al. (2006); Boettger at al. (2012); Sun et al. (2013); Garai et al. (2007); Brian et al. (2008); Matsui (2010); Kunc et al. (2010); Jin et al. (2011); Karbasi et al. (2011); Sokolova et al. (2013). The studies of gold are connected with its many technical applications, for instance, as a protective coating material, in electronics, as well as in metrology and medicine. For instance, gold provides an accurate pressure calibration for experiments conducted at high temperatures and pressures Heinz et al. (1984); Shim et al. (2002); Takemura et al. (2008); Yokoo et al. (2009); Jin et al. (2011); Sokolova et al. (2013). In recent years, gold nanoparticles attracted much attention as a class of multifunctional materials for biomedical application Cobley et al. (2011); Bond (2012 b); Guo et al. (2014).

As far as the experimental studies of gold are concerned, they include equation of state as well as the lattice and thermodynamic properties, both in low-temperature regime Zimmerman at al. (1962); Isaacs (1965); Martin (1966); McLean et al. (1972) and in high temperatures and pressures as well Neighbours et al. (1958); Daniels et al. (1958); Skelskey et al. (1970); Takemura et al. (2008); Kusaba et al. (2009); Yokoo et al. (2009); Matsui (2010). At the same time, the thermodynamic and vibrational properties have been studied by the theoretical methods, including specific heat calculations Gupta (1968); Baria (2004), phonon dispersion relation Brian et al. (2008); Lynn et al. (1973), interatomic interaction and density of states Baria et al. (2003), Debye temperatures Rosén et al. (1983); Baria (2004), Grüneisen parameter Hiki et al. (1967), lattice constants Kunc et al. (2010); Matsui (2010) and equation of state Anderson et al. (1989); Shim et al. (2002); Saxena (2004); Godwal et al. (1989); Moriarty (1995); Boettger (2006); Souvatzis et al. (2006); Boettger at al. (2012); Sun et al. (2013). The theoretical methods include studies of anisotropic-continuum model Hiki et al. (1967), classical Mie-Grüneisen approach and Birch-Murnaghan equation of state Anderson et al. (1989); Yokoo et al. (2009), pseudopotential Baria et al. (2003) and embedded-atom Brian et al. (2008) models, quasiharmonic Hiki et al. (1967) and anharmomic Rosén et al. (1983) approach, as well as first principles Godwal et al. (1989); Moriarty (1995); Boettger (2006); Souvatzis et al. (2006) and density functional theory Kunc et al. (2010) calculations.

Most of the methods mentioned above are very specialized and aimed at characterizing only some of thermodynamic properties, for instance, equation of state, but not at obtaining a full and self-consistent thermodynamic description. Only in few papers an attempt has been made to construct the full thermodynamic description based on the expression for the free energy. The thermodynamic potential considered there contains the elastic (static) energy, vibrational energy, as well as the electronic contribution Greeff et al. (2004); Souvatzis et al. (2006). In this context it is worth mentioning that recently the computational method for the Gibbs energy has been presented, which is suitable for solids under the quasi-harmonic approximation Otero-de-la-Roza et al. (2011). In all these papers the electronic energy is taken in the simplest (Sommerfeld) approximation, and the vibrational energy is calculated from the Debye model in quasi-harmonic approximation. On the other hand, the elastic (static) energy calculations involve more elaborate methods, based on the density functional theory (DFT).

We are aware that DFT methods are not always accessible (and convenient) since they require specialized software and sufficient computational resources. Therefore, apart from referring to such methods it is useful to have also the analytical description of the system in question. Such description should be based on the Gibbs energy, so giving correct thermodynamic relationships, and being uncomplicated enough for practical use. In particular, the analytical form of equation of state (EOS) is very useful for understanding mutual relations between physical parameters, as well as contributions to the total pressure from various energy components. Such EOS should be obtained by demanding minimization of the Gibbs energy functional whereas the system is in equilibrium.

Taking the above needs into account, in the present paper we developed a simple approach, based on the derivation of the Gibbs energy, which yields self-consistent thermodynamic description of metallic system using analytical solutions. The Gibbs energy consists of Helmholtz free-energy and grand potential () term. In turn, the Helmholtz free-energy is constructed from the elastic (static) energy, vibrational energy, and electronic part. The elastic energy is described by the method presented in Ref. Balcerzak et al. (2010), which is based on the expansion of anharmonic potential in the power series with respect to the volume deformation. The vibrational energy is taken within the Debye model in the analytical approximation, which has been extended for high temperatures. In turn, the electronic energy, apart from the kinetic (Sommerfeld) term, has been completed by the ground state energy, as well as by the exchange energy in Hartree-Fock approximation. From the Gibbs energy the full thermodynamic description is obtained and presented in a form of analytical expressions. In particular, new EOS has been found for arbitrary temperature and external pressure . However, the price for the simplicity of the method is that some of initial parameters, defined for and , should be taken from experimental data.

The formalism has been presented in the Theoretical model Section in detail and then applied for bulk gold in the Numerical results and discussion Section. The numerical calculations have been performed in the range of temperatures from 0 K up to the melting point, and external pressure up to 30 GPa. The calculated thermodynamic properties, for instance, specific heat, compressibility and thermal expansion, have been presented in figures and compared with experimental data. Some advantages and weak points of the method have been discussed.

Ii Theoretical model

ii.1 General formulation

The Gibbs free energy of a metallic system is assumed in the form of:

 G=Fε+FD+Fel+pV, (1)

where is the elastic energy without lattice vibrations, is the vibrational energy in Debye approximation and is the energy of electronic subsystem. The elastic energy can be presented in a form of a power series:

 Fε=N(Aε+12Bε2+13!Cε3+14!Dε4+15!Eε5+…) (2)

where is a parameter characterizing volume deformation, are the coefficients and is the number of atoms in the lattice. The volume deformation is defined by , where is the volume at absolute zero temperature () and in vacuum conditions (). This energy is a source of static pressure:

 pε=−(∂Fε∂V)T=−NV0(A+Bε+12Cε2+13!Dε3+14!Eε4+…). (3)

The vibrational energy is taken in the Debye approximation Wallace (1972 b):

 FD=N[98kBTD+9kBT(TTD)3∫yD0y2ln(1−e−y)dy] (4)

where

 yD=TD/T (5)

and is the Debye temperature. Eq. (4) can be transformed into a more convenient form:

 FD=N[98kBTD+3kBTln(1−e−yD)−3kBT1y3D∫yD0y3ey−1dy]. (6)

In further considerations we introduce the lattice Grüneisen parameter in the following way Dorogokupets et al. (2007); Dorogokupets (2010):

 γD=γ0D(V/V0)q=γ0D(1+ε)q, (7)

where is a constant parameter. The Debye temperature is connected with the Grüneisen parameter by the relationship Matsui (2010):

 TD=T0De(γ0D−γD)/qg(T)=T0Df(ε)g(T) (8)

where:

 f(ε)=e(γ0D−γD)/q=eγ0D[1−(1+ε)q]/q. (9)

and are the Debye temperature and Grüneisen parameter, respectively, which are taken at and . In Eq. (8) we introduced function which weakly depends on temperature and will be specified latter. This function reflects the fact that the Debye temperature depends not only on the volume, but also on temperature itself. It takes into account the so-called ”intrinsic anharmonicity” which leads to higher order terms in thermodynamic functions and has also been discussed in several papersOganov et al. (2004); Dorogokupets et al. (2004); Holzapfel (2005); Jacobs et al. (2005); Zhang et al. (2012)
It can be easily checked that the above relationship for satisfies the classical Grüneisen assumption Grüneisen (1912, 1926b):

 γD=−VωD(∂ωD∂V)T=−VTD(∂TD∂V)T. (10)

In the low temperature approximation and we can make use of the integral: . Thus, Eq. (6) can be transformed to the form:

 FD=N[98kBTD−15π4kBT(TTD)3]. (11)

The above energy gives the vibrational pressure for low temperatures:

 pD=−(∂FD∂V)T=3NV0kBTDγD[38+15π4(TTD)4]11+ε (12)

where is given by Eq. (8), and is expressed on the basis of Eq. (10).
On the other hand, for high temperatures, in Eq. (6) we can make use of the following series expansion:

 yey−1=∞∑k=0Bkykk! (13)

where are the Bernoulli numbers: , , , , , , , , , , , , , etc.
This expansion enables us to calculate the integral in the form of a series:

 ∫yD0y3ey−1dy=∞∑k=0Bkk!yk+3Dk+3. (14)

For sufficiently high temperatures, when is small and satisfies the condition: , the series is convergent Dubinov et al. (2008). We found that for practical applications, it is sufficient to include in the series only the terms up to fourth order. Thus, using the high temperature expansion, the vibrational free energy is given by:

 FD=N{98kBTD+3kBTln(1−e−yD) −3kBT[13B0+14B1yD+110B2y2D+1168B4y4D]}. (15)

Such energy gives the following vibrational pressure for high temperatures:

 pD=NV0kBTDγD[32tanh−1(yD/2)−110(1−142y2D)yD]11+ε. (16)

It is worth mentioning that without the second term in the square bracket of Eq. (16), the remaining function reproduces the pressure of oscillators in the Einstein model Balcerzak et al. (2010) with rescaled temperature (, with being the Einstein characteristic temperature).
The electronic free energy in Eq. (1) can be presented in approximate form as:

 Fel=Ne[−Cex√EF+35EF−π241EF(kBT)2], (17)

where is the number of electrons and is the Fermi energy, while the exchange constant is equal to:

 Cex=32πe24πϵ0√2mℏ. (18)

The first term in Eq. (17) corresponds to the exchange energy in the Hartree-Fock approximation Suffczyński (1985 b). The second term is the kinetic energy for , and the last term describes kinetic energy in the low-temperature region, where , denoting the Fermi temperature. The Fermi energy can be presented as a function of the volume:

 EF=E0F1(1+ε)γF (19)

where

 E0F=ℏ22m[3π2NeV0]2/3 (20)

is the Fermi energy at and . In analogy to the Grüneisen assumption for the lattice parameter (see Eq. (10)), -exponent satisfies the equation:

 γF=−VEF(∂EF∂V)T. (21)

In this paper we assume that is a constant parameter (contrary to , which is volume dependent). This parameter can be related to the so-called electronic Grüneisen parameter defined by the formula White et al. (1972):

 γe=−(∂lnDOS(EF)∂lnV)T (22)

where is the density of states at Fermi surface. From Eqs. (21) and (22) one can obtain the relationship:

 γe=1−γF2. (23)

For , which is the free-electron case, we obtain .
From the expression (17) the electronic part of the pressure can be found as:

 pel=−(∂Fel∂V)T=NeV0γF[−12Cex√EF+35EF+π241EF(kBT)2]11+ε, (24)

which is valid for . Under the external pressure , the equation of state follows from the equilibrium condition:

 p=pε+pD+pel, (25)

where is given by Eq. (3), - by Eq. (24), and is given for the low or high temperature regions by Eqs. (12) or (16), respectively.

ii.2 The Grüneisen relationship

For the systems described by -variables we can make use of the exact thermodynamic relationship:

 αpκT=(∂p∂T)V, (26)

where is the thermal volume expansion coefficient:

 αp=1V(∂V∂T)p (27)

and is the isothermal compressibility:

 κT=−1V(∂V∂p)T. (28)

For pressure given in the form of Eq. (25) the temperature partial derivatives can be calculated as follows:

 (∂pε∂T)V=0, (29)
 (∂pD∂T)V=γDCDVV, (30)

where is the phononic heat capacity, and

 (∂pel∂T)V=γFCelVV, (31)

where is the electronic heat capacity. The phononic heat capacity at constant volume for low temperatures can be found from the relationship:

 CDV = −T(∂2FD∂T2)V = 125π4NkB(TTD)3w(T)−9NkBTD[18+π415(TTD)4]Tg(T)∂2g(T)∂T2

where is taken from Eq. (11). The new function is defined by:

 w(T)=[1−Tg(T)∂g(T)∂T]2. (33)

On the other hand, for the high temperature region we have on the basis of Eq. (II.1):

 CDV=3NkB(TDT)2⎡⎣eTD/T(eTD/T−1)2+130−1420(TDT)2⎤⎦w(T) −3NkBTD[1eTD/T−1+12−130TDT+11260(TDT)3]Tg(T)∂2g(T)∂T2. (34)

In further calculations we will assume that function only weakly depends on temperature and has the simple linear form:

 g(T)=1+rTT0D (35)

where is a constant parameter () over the whole temperature region. Then, in right-hand side of Eqs. (II.2) and (II.2) the second terms containing vanish, whereas in the first terms takes the form of:

 w(T)=(1−rT/T0D1+rT/T0D)2. (36)

We see that function for can enforce some decrease of the specific heat vs. temperature in comparison with the case when . We found that such a possibility can be useful in order to reproduce better the experimental data.

The electronic heat capacity at constant volume is given by the formula:

 CelV=−T(∂2Fel∂T2)V=Neπ221EFk2BT (37)

which is valid for low temperatures in the electronic scale.
Substituting the pressure derivatives (Eqs. (29-31)) into Eq. (26) we obtain the Grüneisen relationship for complex (phononic and electronic) system which can be presented as:

 αpκT=γDCDV+γFCelVV=γeffVCV (38)

where

 CV=CDV+CelV, (39)

and

 γeff=γDCDV+γFCelVCDV+CelV. (40)

Therefore, is the effective Grüneisen parameter for the electronic and phononic complex system.

ii.3 The dimensionless equation of state (EOS)

Equation of state (25) can be presented in a dimensionless form which is convenient for numerical calculations. First, we introduce the reference energy for normalization of various energy coefficients. We define:

 A0D=kBT0D, (41)

where is the Debye temperature at 0K. Then, we can introduce the dimensionless pressure:

 ~p=V0N1A0Dp, (42)

and dimensionless temperature:

 ~T=TT0D. (43)

By the same token, the elastic constants can be normalised as follows:

 ~A=AA0D,~B=BA0D,~C=CA0D,~D=DA0D,~E=EA0D. (44)

Similarly, the dimensionless Fermi and exchange energies can be written as

 ~E0F=E0FA0D (45)

and

 ~E0ex=−Cex√E0FA0D, (46)

respectively.
With the above notation the equation of state (Eq. (25)) in the low temperature limit can be presented in the dimensionless form:

 ~p+~A+~Bε+12~Cε2+13!~Dε3+14!~Eε4= 3γDf(ε)g(T)⎡⎣38+15π4(~Tf(ε)g(T))4⎤⎦11+ε +NeNγF⎡⎣12~E0ex1(1+ε)γF/2+35~E0F1(1+ε)γF+π24~T2~E0F(1+ε)γF⎤⎦11+ε. (47)

From the equilibrium condition for and we must have . This condition allows to determine the -coefficient, namely:

 ~A=98γ0D+NeNγF[12~E0ex+35~E0F]. (48)

Now, substituting Eq. (48) into (II.3) we finally obtain the equation of state in the low temperature region as:

 ~p+98γ0D+~Bε+12~Cε2+13!~Dε3+14!~Eε4+12NeNγF~E0ex[1−1(1+ε)γF/2+1]= 3γDf(ε)g(T)⎡⎣38+15π4(~Tf(ε)g(T))4⎤⎦11+ε +35NeNγF~E0F[1(1+ε)γF+1−1]+NeNπ24γF~T2~E0F(1+ε)γF−1. (49)

When , this equation can be linearized with respect to . Then, by comparison with the approximate formula , which is valid in the limit , the coefficient can be determined:

 ~B=V0N1A0Dκ0−98γ0D(γ0D−q+1)−35γF(γF+1)NeN~E0F−12γF(γF2+1)NeN~E0ex, (50)

where is the isothermal compressibility at and . All parameters necessary to calculate from Eq. (50) can be taken from experimental data.

In the high temperature region the vibrational pressure should be taken from Eq. (16) instead of Eq. (12). This leads to the dimensionless EOS for high temperatures:

 ~p+98γ0D+~Bε+12~Cε2+13!~Dε3+14!~Eε4+12NeNγF~E0ex[1−1(1+ε)γF/2+1]= γDf(ε)g(T)1+ε⎡⎣32tanh−1(f(ε)g(T)2~T)−110f(ε)g(T)~T+1420(f(ε)g(T)~T)3⎤⎦ +35NeNγF~E0F[1(1+ε)γF+1−1]+NeNπ24γF~T2~E0F(1+ε)γF−1 (51)

By comparison of the Eqs. (II.3) and (II.3) one can see that only the phononic part on the r.h.s. of these equation has been modified.

ii.4 Calculation of other thermodynamic quantities

From the equation of state, the isotherms ( for ) or isobars ( for ) can be found directly. Other thermodynamic properties result from differentiation of the equation of state. For instance, the thermal volume expansion coefficient is given by Eq. (27), whereas the isothermal compressibility is given by Eq. (28).
The adiabatic compressibility is measured in the experiment. It is given by the definition:

 κS=−1V(∂V∂p)S (52)

where the partial derivative is taken at constant entropy , where .
In turn, the heat capacity at constant pressure is defined as:

 Cp=−T(∂2G∂T2)p. (53)

Having calculated , and on the basis of EOS, the quantities defined above ( and ) can be conveniently obtained from the exact thermodynamic relationships:

 κS=κT/(1+TVα2pCVκT) (54)

and

 Cp=CV(1+TVα2pCVκT). (55)

With the help of the generalized Grüneisen relationship (Eq. (38)) the above formulas can also be presented in a more elegant form:

 κS=κT/(1+Tαpγeff) (56)

and

 Cp=CV(1+Tαpγeff) (57)

where the effective Grüneisen parameter is defined by Eq. (40).
It is worth mentioning that using the above equations the exact thermodynamic identity

 CpCV=κTκS (58)

remains fulfilled.
In order to illustrate the method, some exemplary numerical calculations based on the proposed formalism will be presented in the next Section.

Iii Numerical results and discussion

As an application of the theory presented in previous Section we shall describe thermodynamic properties of bulk gold. For such metallic system many experimental results are available which can be compared with calculations performed within the present model. Some of the experimental data at 0 K can serve as the input parameters for present formalism. For instance, the Debye temperature at 0 K amounts to K Zimmerman at al. (1962). On this basis we can calculate the reference energy (see Eq.( 41)), namely J. Another parameter is the atomic volume at 0 K, , which is estimated from Ref. McLean et al. (1972). The isothermal compressibility at 0 K can be calculated as the inverse of the bulk modulus. On the basis of Refs. McLean et al. (1972); Neighbours et al. (1958) we obtained the value . As far as the Grüneisen coefficient is concerned, we assumed the value of , which is an average value from the data reported in Refs. White et al. (1972) and McLean et al. (1972). Then, the parameter is assumed as , which is one of the possible values considered in Ref. Matsui (2010). Regarding electronic properties, we assumed that the parameter in the whole temperature region, which corresponds to the free electron model. The electron density per atom is equal to . Taking into account the above input data, we calculated the normalized Fermi energy at 0K as and the exchange energy . Moreover, the -parameter can be calculated on the basis of Eq. (50), and yields the value .

Other theoretical parameters, which are necessary for further calculations, are connected with the coefficients in the expression for the elastic energy (Eq. (2)). They are treated as the fitting parameters in our theory. We found that the best fit is obtained for the following set of coefficients: , , for , and , , for . These coefficients reflect the asymmetry of the elastic energy with respect to the sign of deformation . Finally, the -coefficient in - function (defined by Eq. (35)) is assumed as: . Having the above set of starting parameters all thermodynamic properties can be calculated for arbitrary and . For the gold crystal we explored the temperature range , where is the melting temperature, and the pressure was from the range up to several tens of GPa.

We start our calculations with the low-temperature region. In particular, for , on the basis of EOS the volume dependence vs. pressure can easily be obtained. Hence, the lattice constant vs. can be calculated. The result is presented in Fig. 1. We see from that figure that the present result fits well DFT (LDA) calculations Kunc et al. (2010). It also agrees with the experimental point for and , namely McLean et al. (1972). It is worth mentioning that the slope of the curve at (, ) can be related to the isothermal compressibility and well reproduces the experimental value McLean et al. (1972).

In Fig. 2 the thermal volume expansivity of Au, is plotted vs. in the low-temperature region. The results are compared with the experimental data which have been averaged from two sources: Ref. White et al. (1972) and Ref. McLean et al. (1972). An excellent agreement of numerical results with the experimental points is seen.

In Fig. 3 the specific heat at constant pressure is presented vs. temperature in the same low-temperature region (). The present result is compared with the experimental data from Ref. Martin (1966). We also present there the results of calculations based on the semi-empirical formula: . This formula has been found in Ref. Isaacs (1965) and reflects the anomalous behaviour of the specific heat (i.e., the negative coefficient at term) in very low temperature region. From Fig. 3 one can conclude that the agreement of the theoretical results with the experimental data is very satisfactory. In this figure, by the diamond symbols the electron contribution to the calculated specific heat is also shown. This contribution is very small in the range of intermediate temperatures, however, it exceeds the phononic specific heat in the extremely low temperature region, namely when K (at K, J/mole K). On the other hand, at K the electron contribution constitutes only 0.56% of the total specific heat. Due to linear increase vs. temperature the electron fraction of the specific heat will increase again for very high temperatures, where the phononic part tends to saturate.

All calculations in the low-temperature region have been done on the basis of EOS in the form of Eq. (II.3). For higher temperatures, when , the appropriate EOS is given by Eq. (II.3). Fig. 4 presents the thermal expansion coefficient vs. in the range of high temperatures, limited by the melting temperature of bulk gold (K for Linde (1994b)). A comparison of calculations with the experimental data taken from Ref. Anderson et al. (1989) and Ref. Yokoo et al. (2009) (after Ref. Touloukian et al. (1975b)) is made. One can see that present result fits well both sets of experimental data, which have been obtained in different temperature regions. Moreover, for , the calculated thermal expansion coefficient shows better agreement with the experimental data than, for instance, in Ref. Brian et al. (2008).

Regarding compressibility, it should be said that in low-temperature region (K) it is almost constant () and therefore has not been presented. Adiabatic compressibility vs. temperature in the high-temperature region is shown in Fig. 5. The numerical results are compared with the experimental data obtained in Refs. Neighbours et al. (1958) and Yokoo et al. (2009) (after Ref. Collard et al. (1991)). A satisfactory agreement between the present theory and experiment can be noted, although near the melting point some differences are more noticeable. The same remark can be made on Fig. 6, where the specific heat, , is presented vs. temperature. In this case the numerical results are compared with the experimental data taken from Refs. Anderson et al. (1989) and Yokoo et al. (2009) (after Ref. Barin et al. (1973b)). One can see that the agreement between theory and experiment is worse for the highest temperatures, although the difference does not exceed . In particular, for , the calculated specific heat is slightly higher than the experimental one. The same tendency has been observed in Ref. Yokoo et al. (2009). The electron contribution to the specific heat is too small to be presented in Fig. 6 as a separate curve. For instance, at temperatures 100 K, 500 K, 1000 K and 1300 K it amounts to 0.3%, 1.3%, 2.6% and 3.5% of the total specific heat, respectively.

The isotherm curve, describing the volume deformation vs. external pressure, is shown in Fig. 7. The constant temperature amounts to K. The numerical calculations are compared with the experimental data which have been re-calculated from lattice constant measurements Takemura et al. (2008). A non-linear decrease of vs. can be noted, and the negative values of correspond to volume compression. The positive value of for , , is connected with thermal expansion in the temperature range from K (where ) up to K. For higher pressures the calculations become less accurate. The similar discrepancy between theory and experiment for high pressures has been observed in Ref. Karbasi et al. (2011).

In the last figure (Fig. 8) the Debye temperature is plotted as a function of two variables: pressure and temperature. The calculations are based on Eq. 8, whereas is calculated from EOS as a function of and . For and the Debye temperature starts from the value K and decreases when increases. The dependence of on is just opposite; the increasing pressure causes increase of the Debye temperature. Such a behaviour allows to deduce that decrease of the Debye temperature can be obtained by increasing the atomic volume. This conclusion is in agreement with the experimental results of Ref. Kusaba et al. (2009), where the Debye temperature was calculated from Debye-Waller factor by an X-ray diffraction method, and a similar dependence has been reported. However, due to approximate expression for the elastic energy, and specific form of Eq. 8, some of details concerning the Debye temperature behaviour, like an anomaly at low temperatures Gupta (1968); Lynn et al. (1973), have not been reproduced in our calculations.

Iv Summary and final conclusions

In the paper we developed the self-consistent model for thermodynamic description of metallic systems. The idea presented in Ref. Balcerzak et al. (2010), where the Einstein model was combined with the elastic one, has been extended here for the Debye approximation. What is more, the electronic subsystem with its kinetic and exchange energy has been taken into account. The electronic energy has been considered in better approximation than in previous works Greeff et al. (2004); Souvatzis et al. (2006). We have shown that the ground state (Fermi) energy, as well as the exchange energy, both of them being volume dependent, can contribute to the electronic pressure. The regions of low and high temperatures are described by different EOS, however, the whole temperature range has been covered by these equations. Contrary to Birch-Murnaghan equation of state, which presents only an isothermal description Wallace (1972 b), in our EOS all of the variables ( and ) are treated equivalently. The numerical results have been obtained for gold crystal showing satisfactory agreement with the experimental data, as well as with some DFT calculations. The difference between our results and the experimental data does not exceed for the best fit of theoretical parameters.

The greatest difference between the numerical results and experiment turns out to be near the melting point. A possible source of such inaccuracy is the Debye approximation. This approximation has been used for the sake of simplicity, however, it is rather coarse for the real systems at high temperatures. Another reason is that our starting point for the series expansion of elastic energy with respect to is (, ). Thus, our theory is most applicable around the equilibrium point ( and ). For the case of gold it is quite far from the melting point. For instance, in high-pressure equation of state (Birch-Murnaghan) K and bar has been assumed as a reference temperature and pressure point. In our case, expanding the range of pressures up to some extreme values would require higher order terms vs. and new fitting parameters in the elastic potential to be taken into account, as well as other necessary improvements on the presented approach. For instance, the assumption that the elastic coefficients , , , etc., are constants in the whole temperature and pressure region is only an approximation. It has been shown that the elastic coefficients of gold are, to some extent, pressure Daniels et al. (1958) and temperature Neighbours et al. (1958); Garai et al. (2007) dependent. For the above reasons, the application of the method in its present form for the range of extreme pressures would be, in our opinion, rather limited in practice.

Similarly to other models leading to EOS, in our calculations we have used only a single variable for description of the volume elastic deformation. However, the approach can be generalized for anisotropic deformations (also including anisotropic external pressures). It should also be mentioned that our considerations are limited to the quasistatic processes and the shock-wave experiments cannot be described within this model.

In this paper the numerical calculations have been performed for the case of gold only. The analysis for another metal can be done analogously, whereas the numerical calculations of all thermodynamic properties should be performed simultaneously from one set of fitting parameters. Among these parameters there are characteristics of elastic potential (, , , etc.), as well as other parameters (, ) which form an unique set for a given metal. The best fitting of all curves to the experimental data means in practice multiple and time-consuming calculations. However, the numerical calculations and analysis of the properties for other metals is beyond the scope of present paper, which is mainly devoted to the detailed presentation of theoretical model.

In spite of the features mentioned above, in our opinion, the presented model can be useful for thermodynamic description of metallic systems. Its advantage follows from the fact that all thermodynamic properties can be found on the basis of one single expression: the Gibbs energy. Hence, the self-consistency of the theory is preserved and all thermodynamic relationships (like Grüneisen equation) are exactly fulfilled. As mentioned above, the model can be further improved when the volume deformation is treated as anisotropic quantity. On the other hand, the vibrational energy can be taken more accurately than in the Debye approximation, for instance, by better modelling of the phononic dispersion relations and density of states. Also, the electronic energy calculations would benefit from a more realistic model of band structure. For further improvement of the model, electron-phonon interaction might also be taken into account. Of course, such improvements will make the model more accurate, but, at the same time, more complicated for the practical use.

References

• Neighbours et al. (1958) J. R. Neighbours,  G. A. Alers, Phys. Rev. 111, 707 (1958).
• Daniels et al. (1958) W. B. Daniels, and C. S. Smith, Phys. Rev. 111, 713 (1958).
• Zimmerman at al. (1962) J. E. Zimmerman and  L. T. Crane, Phys. Rev. 126, 513 (1962).
• Isaacs (1965) L. L. Isaacs, J. Chem. Phys. 43, 307 (1965).
• Martin (1966) D.  L. Martin, Phys. Rev. 141, 576 (1966).
• Cohen et al. (1966) L. H. Cohen,  W. Clement, Jr., and G. C. Kennedy, Phys. Rev. 145, 519 (1966).
• Skelskey et al. (1970) D. Skelskey and J. Van den Sype, J. Appl. Phys. 41, 4750 (1970).
• McLean et al. (1972) K. O. McLean, C. A. Swenson and C. R. Rose, J. of Low Temp. Phys. 7, 77 (1972).
• White et al. (1972) G. K. White,  J. G. Collins, J. Low Temp. Phys. 7, 43 (1972).
• Lynn et al. (1973) J. W. Lynn,  H. G. Smith, and R. M. Niclow, Phys. Rev. B 8, 3493 (1973).
• Barin et al. ( 1973b) I. Barin and  O. Knacke, Thermophysical Properties of Inorganic Substances, (Springer, 1973b).
• Touloukian et al. ( 1975b) Y. S. Touloukian, R. K. Kirby, R. E. Taylor, and  P. D. Desai, in Thermophysical Properties of Matter, The TPRC Data Series Vol. 12, (Plenum, 1975b).
• Heinz et al. (1984) D. L. Heinz and R. Jeanloz, J. Appl. Phys. 55, 885 (1984).
• Anderson et al. (1989) O. L. Anderson, D. G. Isaak and  S. Yamamoto, J. Appl. Phys. 65, 1534 (1989).
• Collard et al. (1991) S. M. Collard and R. B. McLellan, Acta Metall. Mater. 39, 3143 (1991).
• Takemura et al. (2008) K. Takemura and  A. Dewaele, Phys. Rev. B 78, 104119 (2008).
• Yokoo et al. (2009) M. Yokoo, N. Kawai, K. G. Nakamura, K. Kondo, Y. Tange, and T. Tsuchiya, Phys. Rev. B 80, 104114 (2009).
• Kusaba et al. (2009) K. Kusaba,  T. Kikegawa, Solid State Comm. 149, 371 (2009).
• Hiki et al. (1967) Y. Hiki,  J. F.  Thomas, Jr., and A. V. Granato, Phys. Rev. 153, 764 (1967).
• Gupta (1968) R. P. Gupta, Phys. Rev. 174, 714 (1968).
• Rosén et al. (1983) J. Rosén and G. Grimvall, Phys. Rev. B 27, 7199 (1983).
• Sundqvist et al. (1985) B. Sundqvist,  J. Neve,  Ö. Rapp, Phys. Rev. B 32, 2200 (1985).
• Shim et al. (2002) S.-M. Shim,  T. S. Duffy,  T. Kenichi, Earth and Planetary Science Letters 203, 729 (2002).
• Saxena (2004) S. K. Saxena, J. Phys. Chem. Solids 65, 1561 (2004).
• Baria et al. (2003) J. K. Baria,  A. R. Jani, Physica B 328, 317 (2003).
• Baria (2004) J. K. Baria, Czechoslovak J. Phys. 54, 575 (2004).
• Greeff et al. (2004) C. W. Greeff and M.  J. Graf, Phys. Rev. B 69, 054107 (2004).
• Godwal et al. (1989) B. K. Godwal and R. Jeanloz, Phys. Rev. B 40, 7501 (1989).
• Moriarty (1995) J. A. Moriarty, High Pressure Res. 13, 343 (1995).
• Boettger (2006) J. C. Boettger, Phys. Rev. B 67, 174107 (2003).
• Souvatzis et al. (2006) P. Souvatzis,  A. Delin,  O. Eriksson, Phys. Rev. B 73, 054110 (2006).
• Boettger at al. (2012) J. Boettger, K.G. Honnell, J.H. Peterson, C. Greeff and S. Crockett, AIP Conf. Proc. 1426, 812 (2012).
• Sun et al. (2013) J.-X. Sun, L.-C. Cai, Q. Wu, and K. Jin, Phys. Scr. 88, 035005 (2013).
• Garai et al. (2007) J. Garai and A. Laugier, J. Appl. Phys. 101, 023514 (2007).
• Brian et al. (2008) Q. Brian,  S. K. Bose,  R. C. Shukla, J. Phys. Chem. Solids 69, 168 (2008).
• Matsui (2010) M. Matsui, J. Phys.: Conf Series 215, 012197 (2010).
• Kunc et al. (2010) K. Kunc and  K. Syassen, Phys. Rev. B 81, 134102 (2010).
• Jin et al. (2011) K. Jin,  Q. Wu,  H. Geng,  X. Li,  L. Cai, and X. Zhou, High Pressure Research 31, 560 (2011).
• Karbasi et al. (2011) A. Karbasi,  S. K. Saxena,  R. Hrubiak, CALPHAD: Computer Coupling of Phase Diagram and Thermochemistry 35, 72 (2011).
• Sokolova et al. (2013) T. S. Sokolova,  P. I. Dorogokupets,  K. D. Litasov, Russian Geology and Geophysics 54, 181 (2013).
• Cobley et al. (2011) C. M. Cobley,  J. Chen,  E. Ch. Cho,  L. W. Wang, and Y. Xia, Chem. Soc. Rev. 40, 44 (2011).
• Bond (2012 b) G. C. Bond, in Gold nanoparticles for physics, chemistry and biology, edited by C.  Louis,  O.  Pluchery (Imperial College Press, 2012b).
• Guo et al. (2014) Y.-J. Guo,  G.-M. Sun,  L. Zhang,  Y.-J. Tang,  J.-J. Luo,  P.-H. Yang, Sensors and Actuators B 191, 741 (2014).
• Otero-de-la-Roza et al. (2011) A. Otero-de-la-Roza,  D. Abbasi-Pérez,  V. Luaña, Computer Phys. Commun. 182, 2232 (2011).
• Balcerzak et al. (2010) T. Balcerzak,  K. Szałowski, and M. Jaščur, J. Phys. Condens. Matter 22, 425401 (2010).
• Wallace (1972 b) D. C. Wallace, Thermodynamic of Crystals, (J. Wiley, 1972b).
• Dorogokupets et al. (2007) P. I. Dorogokupets  and A. Dewaele, High Pressure Res. 27, 431 (2007).
• Dorogokupets (2010) P. I. Dorogokupets, Phys. Chem. Miner. 37, 677 (2010).
• Oganov et al. (2004) A. R. Oganov and P. I. Dorogokupets, J. Phys. Condens. Matter 16, 1351 (2004).
• Dorogokupets et al. (2004) P. I. Dorogokupets and A. R. Oganov, Doklady Earth Sciences 395, 238 (2004).
• Holzapfel (2005) W. B. Holzapfel, High Pressure Res. 25, 187 (2005).
• Jacobs et al. (2005) M. H. G. Jacobs and B. H. W. S.  de Jongh, Phys. Chem. Minerals 32, 614 (2005).
• Zhang et al. (2012) D. Zhang and J.-X. Sun, Chin. Phys. B 21, 080508 (2012).
• Grüneisen (1912) E. Grüneisen, Ann. Phys. 344, 257 (1912).
• Grüneisen (1926b) E. Grüneisen, Handbuch der Physik, (Springer, 1926b).
• Dubinov et al. (2008) A. E. Dubinov, and A. A. Dubinova, Technical Physics Letters 34, 999 (2008).
• Suffczyński (1985 b) M. Suffczyński, Electrons in Solids, (Polish Academy of Sciences, Institute of Physics, Ossolineum, 1985b).
• Linde (1994b) CRC Handbook of Chemistry and Physics 1913-1995, 75th edition, edited by D. R.  Linde (CRC Press, 1994b).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters