Second virial coefficient for the Lennard–Jones potential

Physica A 290 (2001) 92–100

www.elsevier.com/locate/physa

Second virial coecient for the Lennard–Jones potential a Departamento

P. Vargasa; b; ∗ , E. Mu˜noza , L. Rodriguezc

de Fsica, Universidad Tecnica F. Santa Mara, Casilla 110-V, Valparaso, Chile Institut fur Festkorperforchung, D-70569 Stuttgart, Germany c Departamento de Fsica, Universidad de Santiago de Chile, Casilla 307, Santiago-2, Chile b Max-Planck

Received 29 December 1999; received in revised form 1 May 2000

Abstract The second virial coecient for the (12−6) Lennard–Jones potential and its rst and second temperature derivatives are explicitely derived. All coecients can be explicitly written as a linear combination of four modi ed Bessel functions which can be evaluated for every temperature with high degree of numerical accuracy. Detailed comparisons with the classical series expansion for the coecient and its temperature derivatives are made. The use of these analytical, simple and closed formulae can be very valuable when using experimental p–V –T data to t a Lennard– c 2001 Elsevier Science B.V. All rights reserved. Jones potential. PACS: 5.20.−y; 5.70.Ce; 5.20.Jj

1. Introduction Many books on statistical mechanics (see, for instance, Refs. [1– 6]) have a chapter devoted to imperfect gases, in which the virial expansion of the equation of state in terms of density is derived. Explicit evaluation of the second virial coecient is given only in simple cases such as hard spheres and square well potentials [1–5]. However, closed forms for the second virial coecient using dierent potentials can be found in the literature. Buckingham potential formula is modi ed in terms of Whittaker functions [7], Lennard–Jones (2n−n) and Woolley potential in terms of parabolic cylinder functions [8,9], Sutherland and particular Mie potentials in terms of generalized hypergeometric functions [10,11]. Also, extensive numerical calculation for the ∗ Correspondence address. Departamento de Fsica, Universidad T ecnica F. Santa Mara, Casilla 110-V, Valparaso, Chile. Fax: +56-32-797656. E-mail address: pvargas@ s.utfsm.cl (P. Vargas).

c 2001 Elsevier Science B.V. All rights reserved. 0378-4371/01/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 0 ) 0 0 3 6 2 - 9

P. Vargas et al. / Physica A 290 (2001) 92–100

93

Fig. 1. Second virial coecient as a function of reduced temperature, as given by Eq. (3) (full-line). Tabulated data from Ref. [2] are also plotted (open circles). The inset shows the percent error between the values obtained using the analytical formula and the tabulated ones.

Lennard–Jones (12−6) potential [6,12,13] as well as its rst derivatives have been done in the past [6]. Although parabolic cylinder, hypergeometric and Whittaker functions are easily programmable, their numerical evaluation can be more involved than working interpolating values from tabulated functions. However, a simple analytical closed form is very useful, especially at low temperatures where the series expansion of the virial coecient converges very slowly. Also simple formula is desirable, for instance, to study the inversion temperatures for real gases in expansion [14]. In this letter we derive an alternative, simpler, analytical form for the second virial coecient for the Lennard–Jones (12−6) potential, and, for the rst time, derive explicit closed forms for its rst and second temperature derivative from which the zero-pressure Joule–Thompson coecient is evaluated.

2. Method If u, the intermolecular potential between two particles, only depends on the relative separation, r, between them the second virial coecient is commonly given by the following expression [1– 6] (Fig. 1): Z ∞ [e−u(r)=kT − 1]r 2 dr : (1) B(T ) = −2 0

94

P. Vargas et al. / Physica A 290 (2001) 92–100

Fig. 2. First temperature derivative of the second virial coecient as a function of reduced temperature, as given by Eq. (4) (full-line). Tabulated data from Ref. [2] are also plotted (open circles). The inset shows the percent error between the values obtained using the analytical formula and the tabulated ones.

Provided that u(r) approaches zero faster than 1=r 3 , necessary condition for existence of the virial, integration by parts of Eq. (1) gives an equivalent expression Z ∞ du(r) −u(r)=kT 3 2 e r dr : (2) B(T ) = − 3kT 0 dr Then, if u(r) is the Lennard–Jones (12−6) potential     12   6 − u(r) = 4 r r and making the following substitution x = r= and T ∗ = kT=, the expression for the second virial coecient of Eq. (2) reads    6 )  12 Z ∞( 3 ∗ 12 6 4 2 1 1 +6 −12 e−4=T {(1=x) −(1=x) } x2 d x : B(T ∗ ) = − 3 T∗ 0 x x By de ning B∗ (T ∗ )=B(T ∗ )=(23 =3), substituting, u=1=x3 , and after some straightforward algebraic manipulations (see Appendix A) we obtain the following dimensionless form for the second virial coecient (Figs. 2 and 3): Z ∞ ∗ ∗ 2 2 8 e1=T (u2 − 1) e−(1=T )(u −1) du : B∗ (T ∗ ) = √ ∗ 2T 0 In Appendix B we demonstrate that the function P(z), de ned as Z ∞ 2 2 (u2 − 1) e−2z(u −1) du ; P(z) ≡ 0

P. Vargas et al. / Physica A 290 (2001) 92–100

95

Fig. 3. Second temperature derivative of the second virial coecient as a function of reduced temperature, as given by Eq. (5) (full-line). Tabulated data from Ref. [2] are also plotted (open circles). The inset shows the percent error between the values obtained using the analytical formula and the tabulated ones.

is a linear combination of four modi ed Bessel functions [15]: P(z) = 18  e−z (I−3=4 (z) − I−1=4 (z) − I1=4 (z) + I3=4 (z)) : Then, it turns out, by replacing z =1=2T ∗ , the second virial coecient for the Lennard– Jones (12−6) potential has the following analytically closed expression: √      2 1=2T ∗ 1 1 ∗ ∗ e I−3=4 + I3=4 B (T ) = 2T ∗ 2T ∗ 2T ∗     1 1 − I−1=4 : (3) −I1=4 2T ∗ 2T ∗ It is now a trivial matter, using the recurrence formulas for the modi ed Bessel function (see 9.6.26 of Abramowitz and Stegun [15]), to derive the following expressions for the temperature derivatives, and zero-pressure Joule–Thompson coecient (Fig. 4). B1∗ (T ∗ ) ≡ T ∗

dB∗ (T ∗ ) ; dT ∗

     1 √ 1 1 1=2T ∗ 2e I3=4 + I−3=4 8T ∗ 2T ∗ 2T ∗     1 1 −3I−1=4 ; −3I1=4 ∗ 2T 2T ∗

B1∗ (T ∗ ) = −

B2∗ (T ∗ ) ≡ T ∗2

∗ ∗ d 2 B∗ (T ∗ ) ∗ dB1 (T ) = T − B1∗ (T ∗ ) ; dT ∗2 dT ∗

(4)

96

P. Vargas et al. / Physica A 290 (2001) 92–100

Fig. 4. B1∗ (T ∗ ) − B∗ (T ∗ ) (reduced zero-pressure Joule–Thompson coecient) as a function of reduced temperature, as given by Eq. (1) and Eq. (2) (full-line). Tabulated data from Ref. [2] are also plotted (open circles). The inset shows the percent error between the values obtained using the analytical formula and the tabulated ones.

B2∗ (T ∗ ) =

      ∗  1 1 √ e1=2T (−4 + 5T ∗ ) I3=4 + I −3=4 2T ∗ 2T ∗ 16 2T ∗2      1 1 ∗ + I1=4 : −(4 + 21T ) I−1=4 2T ∗ 2T ∗

(5)

Zeros for the dierent thermodynamical functions are in T ∗ = 3:417 928 023; 25:152 573 456, 48:289 836 984 and 6:430 798 472 for B∗ ; B1∗ ; B2∗ , and B1∗ − B∗ , respectively.

3. Results Figs. 1– 4 illustrate the tabulated second virial coecients and zero-pressure Joule– Thompson coecient as given by Hirschfelder et al. [6] in their Table I-B (open circles), together with the analytical expressions given in Eqs. (3) – (5) (full line curves). For a given reduced temperature T ∗ , the percent error in the case of the second virial coecient is de ned as 100 × [exact(B∗ ) − tabulated(B∗ )]=exact(B∗ ). Same de nition applies for the other functions and they are shown as insets in each gure. We see that the overall behavior coincides almost perfectly in all temperature ranges. Maximum dierences are found in the derivatives which are nevertheless less than 0:004%. These results clearly demonstrate the high quality of the numerical work performed almost 50 years ago.

P. Vargas et al. / Physica A 290 (2001) 92–100

97

4. Conclusion Second virial coecient for the Lennard–Jones (12−6) potential, as well as its rst and second temperature derivatives have been derived for the rst time from simple and closed formulas. We expect the expression as given in Eqs. (1) – (3) will be of great help in tting a Lennard–Jones (12−6) potential from experimental data. Nowadays they can be easily evaluated in every temperature range with high-numerical accuracy using the modi ed Bessel functions which come precompiled in the majority of numerical library packages or scienti c softwares. Acknowledgements This work was supported by the Direccion de Investigacion Cient ca y Tecnologica (DICYT) from the Universidad de Santiago de Chile, FONDECYT under contrat 1990812 and the Max-Planck Institut fur Festkorperphysik, Stuttgart, Germany. Appendix A According to the text de nition, we have  6 )  12 Z ∞( ∗ 12 6 4 1 1 ∗ ∗ +6 −12 e−(4=T ){(1=x) −(1=x) } x2 d x : B (T ) = − ∗ T 0 x x By making the substitution u = 1=x3 we obtain Z ∞ ∗ 4 2 8 B∗ (T ∗ ) = ∗ (2u2 − 1) e−4=T (u −u ) du ; T 0 √ and now, making u = x= 2 we nally get the desired expression: Z ∞ ∗ ∗ 2 2 8 B∗ (T ∗ ) = √ e1=T (u2 − 1) e−1=T (u −1) du : 2T ∗ 0 Appendix B We recall the integral representation of the modi ed Bessel function I (see formula 9.6.20 of Abramowitz and Stegun [15]). Z Z 1  z cos t sin  ∞ −z cosh t−t e cos t dt − e dt : I =  0  0 Then it follows that Z ∞ e−z cosh t sinh t dt = 0

 1 (I + I− ) − 2 sin  sin 

which we will use for the demonstration.

Z 0



ez cos t cos t dt

(6)

98

P. Vargas et al. / Physica A 290 (2001) 92–100

We have to evaluate analytically P(z) de ned as Z ∞ 2 2 (u2 − 1) e−2z(u −1) du : P(z) ≡ 0

By substituting x = u2 it follows that Z 1 ∞ (u − 1) −2z(u−1)2 √ du : e P(z) = 2 0 u Let us rst separate P(z) as a summation of two terms. The interval splitting [0; 2] and [2; ∞] will be clear later in the demonstration: Z Z 1 ∞ (u − 1) −2z(u−1)2 1 2 (u − 1) −2z(u−1)2 √ √ e e du + du P(z) = 2 0 2 2 u u = P1 (z) + P2 (z) : In the last term, we do the following transformation, which is valid for u ¿ 2: (u − 1)2 =

t 1 1 + cosh t = cosh2 ; 2 2 2

u = 1 + cosh

t (u ¿ 2) ; 2

2(u − 1) du = 12 dt sinh t : Then,

Z

Z ∞ (u − 1) −2z(u−1)2 sinh t −z cosh t 1 √ e du = √ e−z dt ; e cosh t=4 u 8 2 2 0 Z t t 1 −z ∞ sinh cosh e−z cosh t dt = √ e 4 2 2 2 0  Z ∞ t 3t 1 e−z cosh t dt : sinh − sinh = √ e−z 4 4 4 2 0

P2 (z) =

1 2

∞

Now, doing the following transformation in the rst term, P1 (z); valid for 0 ¡ u ¡ 2, (u − 1)2 =

1 1 t + cos t = cos2 ; 2 2 2

u = 1 + cos

t (0 ¡ u ¡ 2) ; 2

2(u − 1)du = − 12 dt sin t ; Then 1 P1 (z) = − 8

Z

e−z = √ 2 2

0

2

Z 0

p 2

sin t 1 + cos t=2

−2z((1=2)+(1=2) cos t)

e

t e−z t sin cos e−z cos t dt = √ 4 2 4 2

e−z dt = √ 8 2

Z 2  0

Z 0

2

sin t −z cos t e dt ; cos t=4

t 3t sin − sin 4 4



e−z cos t dt :

P. Vargas et al. / Physica A 290 (2001) 92–100

99

These expressions lead to the following expression for P(z):  Z 2  3t t e−z e−z cos t dt sin − sin P(z) = √ 4 4 4 2 0   Z ∞ t 3t −z cosh t e dt : sinh − sinh + 4 4 0 From the derived expression for the modi ed Bessel functions of Eq. (6), using as particular cases  = 34 , and 14 we get Z ∞ √ Z  z cos t 1 √ 3t 3t e−z cosh t sinh dt =  2(I3=4 + I−3=4 ) − 2 e cos dt ; 4 2 4 0 0 Z Z ∞  √ 1 √ t t e−z cosh t sinh dt =  2(I1=4 + I−1=4 ) − 2 ez cos t cos dt : 4 2 4 0 0 Therefore by replacing them in the P(z) expression, we have  Z 2  t 1 √ e−z 3t e−z cos t dt +  2(I3=4 + I−3=4 ) P(z) = √ sin − sin 4 4 2 4 2 0  t dt 4 0 0   1 √ e−z 1 √ = √ Q(z) +  2(I3=4 (z) + I−3=4 (z)) −  2(I1=4 (z) + I−1=4 (z)) : 2 2 4 2 To see that Q = 0, let us write    Z 2  √ Z  z cos t t t 3t 3 e−z cos t dt − 2 dt : e sin − sin cos t − cos Q= 4 4 4 4 0 0 √ Z − 2



ez cos t cos

√ 3t 1 √ dt −  2(I1=4 + I−1=4 ) + 2 4 2

Z



ez cos t cos

and substituting t = x +  in the rst integral and then expanding the trigonometrical functions, we have     Z  3 3t 3 3t cos + cos sin sin Q= 4 4 4 4 −    t  t  z cos t cos − cos sin e dt − sin 4 4 4 4   √ Z  z cos t t 3t dt : e cos − cos − 2 4 4 0 Odd functions give zero in the integrals from [−; ], then by transforming the integrals of the even functions in the [ − ; ] interval by twice the same integrals in the [0; ] interval, we get   √ Z  t 3t ez cos t dt cos − cos Q= 2 4 4 0   √ Z  z cos t 3 1 cos t − cos t dt = 0 : e − 2 4 4 0

100

P. Vargas et al. / Physica A 290 (2001) 92–100

Thus, we have proved that P(z) =

 e−z ((I3=4 (z) + I−3=4 (z)) − (I1=4 (z) + I−1=4 (z))) : 8

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

M.C. McQuarrie, Statistical Thermodynamics, Harper & Row, New York, 1976. R.P. Feynman, Statistical Mechanics, a Set of Lectures, Benjamin, New York, 1972. S. Ma, Statistical Mechanics, World Scienti c, Singapore, 1985. C.V. Heer, Statistical Mechanics, Kinetic Theory and Stochastic Processes, Academic Press, New York, 1972. M. Plischke, B. Bergensen, Equilibrium Statistical Physics, 2nd Edition, World Scienti c, Singapore, 1994. J.O. Hirschfelder, C.F. Curtiss, R.B. Bird, Molecular Theory of Gases and Liquids, Wiley, New York, 1954. S.F. Ragab, A.A. Helmy, T.L. Hassanein, M.A. El-Naggar, J. Low-Temp. Phys. 111 (1998) 447. H. Guerin, J. Phys. B 26 (1993) L693. A.J.M. Garrett, J. Phys. A 13 (1980) 379. D. Levi, M. de Llano, J. Chem. Phys. 63 (1975) 4561. E. Ley-Koo, M. de Llano, J. Chem. Phys. 65 (1976) 3802. L.F. Epstein, G.M. Roe, J. Chem. Phys. 19 (1951) 1320. J.A. Barker, P.J. Leonard, A. Pompe, J. Chem. Phys. 44 (1966) 4206. D. Huang, Physica A 256 (1998) 30. M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions, U.S. National Bureau of Standards, 1964; Dover, New York, 1965.