# First-principles calculation of the third virial coefficient of helium.

Knowledge of the pair and three-body potential-energy surfaces of
helium is now sufficient to allow calculation of the third density
virial coefficient, C(T), with significantly smaller uncertainty than
that work, we employ the best available pair and three-body potentials
for helium and calculate C(T) with path-integral Monte Carlo (PIMC)
calculations supplemented by semiclassical calculations. The values of
C(T) presented extend from 24.5561 K to 10 000 K. In the important
metrological range of temperatures near 273.16 K, our uncertainties are
smaller than the best experimental results by approximately an order of
magnitude, and the reduction in uncertainty at other temperatures is at
least as great. For convenience in calculation of C(T) and its
derivatives, a simple correlating equation is presented.

Key words: calibration; density; helium; metrology; path integral; thermodynamic properties; virial coefficients.

1. Introduction

Accurate knowledge of the thermophysical properties of helium is desirable for many applications in metrology. Examples include microwave resonance measurements for development of a pressure standard (1), measuring the Boltzmann constant (2), and acoustic gas thermometry (3).

In many cases, the nonideal behavior of the gas is a major component of the overall uncertainty. That non-ideality is expressed in the virial expansion

[p/[[rho]RT]] = 1 + B(T)[rho] + C(T)[[rho].sup.2]+ ..., (1)

where p is the pressure, [rho] the molar density, R the molar gas constant, and T the absolute temperature. B(T) is the second virial coefficient, representing the lowest-order deviation from ideal-gas behavior. As the density increases, the contribution from the third virial coefficient C(T) becomes significant. B(T) is a function only of the interactions between pairs of molecules, while C(T) is determined by interactions among three molecules.

Because the helium atom has only two electrons, and because of advances in methodology, algorithms, and computing power for ab initio quantum calculations, it is now possible to calculate properties of individual helium atoms and pairs of atoms with very high accuracy. The use of such calculated values to develop standards for thermophysical properties was first proposed by Aziz et al. (4), after which Hurly and Moldover (5) calculated helium's second virial coefficient B(T), dilute-gas viscosity, and dilute-gas thermal conductivity with uncertainties smaller than those of the best experiments. Hurly and Mehl (6) recently improved on this work, primarily by using a more accurate pair potential. Bich et al. (7) performed similar calculations of pair quantities with another high-accuracy pair potential (8). They also calculated C(T) with quantum effects considered at a first-order perturbation level, but (as explained in Sec. 2.2) the Axilrod-Teller term they used for the three-body potential is not a good representation of three-body effects for helium.

For the third virial coefficient C(T), rigorous first-principles calculations are more difficult. While classical calculation of C(T) for spherically symmetric species is straightforward (9), quantum effects must be considered for a light gas such as helium. Semiclassical perturbation approaches have been derived (10), (11), but they will be insufficient at low temperatures and their uncertainty is difficult to quantify. What is needed is a fully quantum approach as is available for B(T) (9), but this problem has not been solved in closed form for C(T).

Early quantum calculations, aimed to elucidate how C(T) depends on the inter-molecular potentials, were performed at the beginning of the 1960s, using the general theoretical framework developed by Lee and Yang (12). Pais and Uhlenbeck (13) as well as Larsen (14) discussed various approximations to the full quantum problem, and these results were used to estimate the binding energy. Other calculation attempts tried to imitate the exact quantum calculations of B(T) based on scattering phase shifts. Larsen and Mascheroni (15) were able to obtain a rigorous expression for the third virial coefficient under the (unrealistic) assumption of the absence of bound states. Reiner (16) developed a calculation method based on the Faddeev equations for the quantum-mechanical three-body problem, but the difficulties in solving the coupled integral equations could not be overcome, and approximations were required.

The first rigorous quantum calculation of virial coefficients was developed by Fosdick and Jordan (17-19) who derived a path-integral expression for the third virial coefficient in the presence of two-body pairwise additive forces, and showed in detail how to evaluate it numerically in the case of a simple Lennard-Jones model for He.

Fosdick and Jordan's approach was independently rediscovered by Diep and Johnson (20) who devised, by analogy to the classical expression, a path-integral formula for the second virial coefficient of a quantum gas. They used their expression to calculate the second virial coefficient for a new [H.sub.2]-[H.sub.2] potential that they had computed using ab initio calculations, neglecting the rotational degrees of freedom. Their expression was generalized to the case of asymmetric rotors by Schenter (21) who used it to calculate the second virial coefficient of a model of water.

Recently, we developed a rigorous path-integral Monte Carlo procedure to calculate the second virial coefficient of molecular hydrogen, extending the approach pioneered by Fosdick and Jordan (17) to take into account the rotational degrees of freedom of linear molecules. We used this method together with state-of-the-art ab initio calculations of the pair potential to obtain good agreement with experimental data in a wide range of temperatures (22).

In the present work, we further extend the path-integral approach employed by Jordan and Fosdick for spherical particles (18) to calculate the third virial coefficient of [.sup.4]He in the presence of nonadditive three-body interactions. We use this method to calculate the third virial coefficient of [.sup.4]He from 24.5561 K to 10 000 K, using a recent ab initio derived three-body potential (23).

2. Intermolecular Potentials

2.1 Pair Potential

We write the pair potential as [U.sub.2](r), where r is the center-to-center distance between the atoms. We use the pair potential known as [[empty set].sub.07], which was developed by Hurly and Mehl (6) based on the best ab initio calculations available in 2007. For uncertainty analysis, we also use their potentials [[empty set].sub.07.sup.-] and [[empty set].sub.07.sup.+], which represent uncertainty limits (expanded uncertainty with coverage factor k = 2) for the [[empty set].sub.07] potential.

While this work was in progress, a pair potential of higher accuracy (consistent within mutual uncertainties with [[empty set].sub.07]) was published by Jeziorska et al. (24). Because our uncertainties are dominated by the three-body potential (see Sec. 4.2), our results would not have been significantly different had we used that pair potential.

2.2 Three-Body Potential

Calculation of C(T) also requires knowledge of the nonadditive three-body contribution to the potential energy in a triplet of atoms. We write this as [U.sub.3]([r.sub.12], [r.sub.13], [r.sub.23]), where [r.sub.ij] is the distance between atoms i and j. The three-body potential of helium has been studied by Cencek et al. (23), who developed separate ab initio potentials based on symmetry-adapted perturbation theory (SAPT) and coupled-cluster (CC) calculations. The two potentials were estimated to have a maximum uncertainty of 10 %, but very recent work [25] has shown that the SAPT potential is significantly less accurate, while the CC potential is in good agreement (better than 10 %) with higher-level calculations. We therefore employ the CC potential in this work. For purposes of uncertainty analysis, we interpret the estimated 10 % uncertainty as an expanded uncertainty at the k = 2 level. In analogy to the procedure of Hurly and Mehl [6] for the pair potential, we define for our uncertainty analysis three-body potentials CC- (obtained by multiplying the potential by 1.1 where it is negative and by 0.9 where it is positive) and CC+ (multiplying by 0.9 where it is negative and by 1.1 where it is positive).

It is worth noting that, at most temperatures of interest, the main three-body effect for helium is the non-additivity in repulsive forces. The commonly used Axilrod-Teller three-body term [26] accounts only for dispersion interactions. At all but the lowest temperatures considered in this work, the reduced temperature is quite high (compared to the characteristic temperature for the pair dispersion interaction of approximately 11 K), and three-body dispersion effects are less important than three-body repulsion effects. Therefore, an Axilrod-Teller term would seriously misrepresent the three-body potential, even giving a contribution to C(T) of the wrong sign at moderate and high temperatures.

3. Calculation Methods

Let us denote by [Q.sub.N](V,T) the partition function of N particles in a volume V at a temperature T. By defining the quantities [Z.sub.N] as

[[Z.sub.N]/[N!]][equivalent to][[[Q.sub.N](V, T)[V.sup.N]]/[[Q.sub.1][(V, T).sup.N]]], (2)

then the expressions for the second and third virial coefficients, B(T) and C(T), become [9]

B(T) = -[1/[2V]]([Z.sub.2] - [Z.sub.1.sup.2]), (3)

C(T) = 4[B.sup.2](T) - 1/3V([Z.sub.3] - 3[Z.sub.2][Z.sub.1] + 2[Z.sub.1.sup.3]). (4)

3.1 Classical and Semiclassical Calculations

The explicit expression for the partition functions and hence the coefficients [Z.sub.k] appearing in Eqs. (3) and (4) depends on the framework in which the calculations are performed. In classical statistical mechanics (including the correct Boltzmann counting) one has

[Z.sub.class.sup.N] = [integral][d.sup.3][x.sub.1] ... [d.sup.3][x.sub.N]exp( - [beta]U([x.sub.1], ...,[x.sub.N])), (5)

where [beta] = [1/[k.sub.B]]T, [k.sub.B] denotes the Boltzmann constant and U([x.sub.1], ..., [x.sub.N]) is the total potential energy of a configuration with N particles at the positions [x.sub.1], ..., [x.sub.N].

In this case, the third virial coefficient [C.sup.class] (T) is given by the sum of a term depending on the two-body potential [U.sub.2](r), [C.sub.[2-body].sup.class] (T), and a term depending on the three-body potential, [C.sub.[3-body].sup.class], given by

[C.sup.class](T) = [C.sub.[2-body].sup.class](T) + [C.sub.[3-body].sup.class](T), (6)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

We have denoted by [theta] the angle between the vectors [r.sub.12] and [r.sub.13]. The distance between the particles labeled 2 and 3 is therefore given by

[r.sub.23] = [square root of [r.sub.12.sup.2] + [r.sub.13.sup.2] - 2[r.sub.12][r.sub.13]cos[theta]].(9)

The classical formulae are accurate enough for heavy particles at high temperature. If this limit is not attained, quantum diffraction effects (Heisenberg uncertainty) become appreciable, as is the case for hydrogen and helium at and below room temperature. As long as the quantum effects can be considered a small correction to the classical behavior, the expressions given above can be corrected by including the first term in the expansion of the full quantum expression in even powers of h.

The expression for the first quantum correction to the third virial coefficient has been evaluated by Yokota (10). Setting [a.sup.2] = [h.sup.2][beta]/m, where m is the mass of the particles under consideration, the semiclassical expression for the third virial coefficient turns out to be

[C.sup.semi](T) = [C.sup.class](T) + [a.sup.2][C.sub.1](T), (10)

[C.sub.1](T) = 8[B.sub.1](T)[B.sup.class](T) - 2V[B.sub.1](T) + [[[beta].sup.2]/36][integral][d.sup.3][r.sub.12][d.sup.3][r.sub.13][e.sup[beta]U](KUU) (11)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

KUU = [[nabla].sub.12]U * [[nabla].sub.12]U + [[nabla].sub.13]U * [[nabla].sub.13]U + [1/2]([[nabla].sub.12]U * [[nabla].sub.13]U + [[nabla].sub.13]U * [[nabla].sub.12]U). (13)

The classical and semiclassical values of C(T) have been obtained by direct numerical integration of Eqs. (6) and (10). We have used the QAG adaptive algorithm together with the Gauss-Kronrod 15-point rule as implemented in the GNU Scientific Library (27). The interaction has been neglected beyond a cutoff length of L = 4 nm. We have checked that using a larger cutoff value does not appreciably affect the results, and the same cutoff was also used in performing the path-integral Monte Carlo calculations described below.

3.2 Path-Integral Monte Carlo

In the framework of non-relativistic quantum statistical mechanics, the expression for the quantity [Z.sub.N] of Eq. (2) becomes

[Z.sub.N] = [[LAMBDA].sub.m.sup.3N][summation over k][summation over [pi]]<k | exp(-[beta][[^.H].sub.N])[[^.P].sub.[pi]] | k>, (14)

where [[^.H].sub.N] is the N-body Hamiltonian operator, |k> denotes a complete set of N-particle states and [[^.P].sub.[pi]] is a permutation operator, with the proper sign to take into account the bosonic or fermionic nature of the particles. [[LAMBDA].sub.m] = [h/[square root of [2[pi]m[k.sub.B]T]] denotes the thermal de Broglie wavelength for a particle of mass m. In the following, we will not be concerned with temperatures so low that exchange effects play a relevant role [28] and hence we will approximate the sum over [[^.P].sub.[pi]] with the identity operator (Boltzmann statistics).

Equation (14) is the starting point to derive the path-integral expressions for the second and third virial coefficients, using Eqs. (3) and (4). In order to avoid cumbersome notation, we will present the derivation in detail in the case of the second virial coefficient. This will allow us to establish some useful notation that will be used to present in the most compact form possible the path-integral formulae for the third virial coefficient.

3.2.1 Second Virial Coefficient

The path-integral formula for the second virial coefficient of Eq. (3), using the quantum-mechanical expression of Eq. (14) in the case of Boltzmann statistics, is readily obtained by first performing a canonical transformation to the center of mass R [equivalent to] [([x.sub.1] + [x.sub.2)]/2] and relative r [equivalent to] [r.sub.1] [equivalent to] [x.sub.2] - [x.sub.1] coordinates. In the equation for [Z.sub.2], the kinetic energy relative to the center-of-mass motion commutes with the kinetic and potential energy of the relative motion and can be integrated out, obtaining a factor [V/[[LAMBDA].sub.M.sup.3]], where M= 2m. As a consequence, B(T) is proportional to the trace of the Hamiltonian describing the relative motion [[^.H].sub.r] = [[^.p].sub.r.sup.2]/2[micro]] + [[^.U].sub.2] (|r|) [equivalent to] [^.T] + [[^.U].sub.2], where [mu] = [m/2] is the reduced mass of the pair and [[^.p].sub.r] is the momentum conjugated to the relative coordinate r.

The resulting expression [Z.sub.2] = [[[LAMBDA].sub.m.sup.6]/[2[[LAMBDA].sub.M.sup.3]]][integral][d.sup.3]r<r | exp(-[beta][[^.H].sub.r]) | r>, (15)

can then be evaluated by performing a Trotter expansion of the kinetic and potential energies of the relative motion,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

keeping a finite value of the Trotter index P and inserting P--1 completeness relations 1 = [integral] [dr.sub.i]|[r.sub.i] > <[r.sub.i]| (i = 2, ...,P) between each of the P factors in Eq. (16). The matrix elements of the kinetic energy operator can be evaluated explicitly, obtaining

<[r.sub.i] | exp(-[[beta][[^.p].sup.2]/[2[micro]P]]) | [r.sub.j]> = [[P.sup.3/2]/[[LAMBDA].sub.[mu].sup.3]]exp(-[K/2][([r.sub.i] - [r.sub.j]).sup.2]), (17)

where K = [2[pi]P/[[LAMBDA].sup.[mu].sub.2]]. The final outcome of this chain of equivalences is to map the calculation of the quantum partition function of Eq. (15) to the calculation of the classical partition function of ring polymers with P beads each (29). The resulting expression can be simplified by introducing the coordinates r = [r.sub.1], [DELTA][r.sub.i] = [r.sub.[i+1]] - [r.sub.i] (i = 1, ..., P--1) and letting

[[bar.U].sub.2](r) = [1/P][P.summation over (i = 1)][U.sub.2]([r.sub.i]), (18)

indicate the average of the two-body potential over the positions occupied by the P beads of a given ring polymer. Performing analogous manipulations in the rather trivial case of [Z.sup.2.sub.1] produces the following expression for the second virial coefficient:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

where [F.sub.ring] is the probability of finding a ring polymer configuration in the ideal gas phase, as shown in Ref. (30), and is given by

[F.sub.ring]([mu];[DELTA][r.sub.1], ...,[DELTA][r.sub.[P-1]]) = [[LAMBDA].sub.[mu].sup.3][([[P.sup.3/2]/[[LAMBDA].sub.[mu].sup.3]]).sup.P] x exp[-[K/2][P.summation over (i = 1)][DELTA][r.sub.i.sup.2]], (20)

where [DELTA] [r.sub.p] = [r.sub.p] - [r.sub.1].

The second virial coefficient can also be written as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

where the effective potential is defined as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (22)

that is, by averaging the factor exp (- [beta] [[bar.U].sub.2]) over ring polymer configurations drawn from an ideal gas distribution and having one bead at the point r. In the classical limit h [right arrow] 0, the ring polymers shrink to a point, and hence Eq. (21) recovers the classical result.

3.2.2 Third Virial Coefficient

The same reasoning that was followed to derive the path-integral expression for the second virial coefficient B(T) in Eq. (19) can be applied to the case of the third virial coefficient. In this case it is useful to evaluate the expectation values over three-body operators after having performed a canonical transformation to the Jacobi coordinates R, r, [rho] and the corresponding momenta P, p, [pi], defined as:

R = [1/3]([x.sub.1] + [x.sub.2] + [x.sub.3]); P = [p.sub.1] + [p.sub.2] + [p.sub.3] (23)

r = [x.sub.2] - [x.sub.1]; p = [1/2]([p.sub.2] - [p.sub.1]) (24)

[rho] = [x.sub.3] - [1/2]([x.sub.1] + [x.sub.2]); [pi] = [1/3](2[p.sub.3] - [p.sub.1] - [p.sub.2]). (25)

As in the case of B(T) the center-of-mass motion can be integrated out, but the Trotter factorization introduces two different ring polymers, corresponding to the coordinates r and [rho]; we note in passing that the masses associated with these degrees of freedom are [M.sub.r] = m/2 and [M.sub.[rho]] = 2m/3, respectively, and that the total potential energy of a three-body configuration,

U([x.sub.1], [x.sub.2], [x.sub.3])[equivalent to][U.sub.3](|[x.sub.2] - [x.sub.1]|,|[x.sub.3] - [x.sub.2]|,|[x.sub.3] - [x.sub.1]|) + [3.summation over (i< j = 1)][U.sub.2](|[x.sub.i] - [x.sub.j]|), (26)

is a function of the coordinates r and [rho] only. As a final result, the expression for C(T) can be written as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

where we have defined, analogously to Eq. (18),

[bar.U](r; [rho]) = [1/P][P.summation over (i = 1)]U([r.sub.i]; [[rho].sub.i]), (28)

as the average of the potential energy of the three bodies over the positions indicated by the beads of the two given ring polymers corresponding to the Jacobi coordinates r = [r.sub.1] and [rho] = [[rho].sub.1]. The three exponentials of [[bar.U].sub.2] appearing in Eq. (27) come from the three terms [Z.sub.2] [Z.sub.1.sup.2] of Eq. (4), when the two-body integral is written as a function of the coordinates

[r.sub.21] = [x.sub.2] - [x.sub.1] = r, (29)

[r.sub.31] = [x.sub.3] - [x.sub.1] = [rho] + [r/2], (30)

[r.sub.32] = [x.sub.3] - [x.sub.2] = [rho] - [r/2], (31)

respectively.

The integrals over the variables [DELTA][r.sub.i] and [DELTA][[rho].sub.i] allow one to define effective two- and three-body potentials as averages over the two kinds of ring polymers, analogously to Eq. (22)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

where x can be either r or [rho] + r/2 or [rho]-r/2. The final expression for C(T) is hence

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (34)

Given the rotational symmetry of the system, the volume of integration can be written as

[d.sup.3][rd.sup.3][rho] = 4[pi][r.sup.2]dr2[pi][[rho].sup.2]d[rho]d cos[THETA], (35)

where r and [rho] are the moduli of the vectors r and [rho], respectively, while [THETA] is the angle between them.

In the actual calculation, it is useful to write the square of the second virial coefficient as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (36)

so that the difference of two integrals in Eq. (34) can be evaluated as the integral of the difference. Note that in Eq. (36) the coordinate [rho] is associated with a particle having mass [M.sub.r] and not [M.sub.p] as in Eqs. (32) and (33).

The calculation of the third virial coefficient in the path-integral formalism follows directly from Eq. (34), which shows that C(T) is given as a three-dimensional integral using effective two- and three-body potentials, given by Eqs. (33) and (32), respectively. Since for each of the atomic positions the effective potentials are obtained as a costly average over configurations of ring polymers, we chose to evaluate Eq. (34) using a Monte Carlo integration procedure, namely the VEGAS algorithm (27), (31), (32). We used N = 5 x [10.sup.5] integration points for the production stage of the algorithm and half as many for the equilibration stage.

For each of the atomic configurations considered in the course of the Monte Carlo integration, we generate n = 200 ring polymers for the r and [rho] coordinates, distributed according to the probability [F.sub.ring]. This can be done very efficiently using an interpolation formula due to Levy (18), (33). For each of the ring polymers, we evaluate the corresponding average potentials [[bar.U].sub.2] and [bar.U] and accumulate their Boltzmann factors to calculate the effective potentials for the given configuration according to Eqs. (32) and (33).

In order to estimate the statistical uncertainty of the values of the third virial coefficient so obtained, we perform 16 of these calculations for each value of T, starting with different seeds for the random number generator.

The final value for C(T) is obtained as the average of the values coming from the 16 independent calculations, and its Type A uncertainty is estimated as the standard error of the mean from the same set of values. Notice that this uncertainty takes into account the statistical error resulting from use of both a finite N and a finite n in the calculations.

Finally, let us discuss the choice of the Trotter index P. Since the path-integral method is exact in the limit of large P, we fixed this number by calculating B(T) with our method for progressively increasing values of the Trotter index P, until our results matched those performed with the phase-shift method (6). We found that the choice P = 7 + 2400 K/T was enough to reach convergence in a wide range of temperatures, and we similarly checked that the same value of P was sufficient in the case of C(T) at 273.16 K. The offset of 7 in the expression for P is due to the approximation inherent in Levy's interpolation formula (18), (33).

4. Results

4.1 Third Virial Coefficients

We calculated C(T) as described in Sec. 3.2 for the potential-energy surface obtained by combining the [[empty set].sub.07] pair potential referenced in Sec. 2.1 with the three-body potential (CC) referenced in Sec. 2.2. Calculations were performed over a wide range of temperatures, as shown in Table 1. In addition to round values, some temperatures were chosen due to their importance in metrology. For example, our lowest temperature of 24.5561 K is the value assigned to the triple point of neon in the International Temperature Scale of 1990 (ITS-90) (34). It should be noted that T in our calculations is the thermodynamic temperature, which may differ on the order of 0.005 % from the corresponding ITS-90 temperature (3). We included temperatures up to 10 000 K, in order to match the range covered by Hurly and Mehl (6) for B(T). Temperatures below those shown in Table 1 were not achievable with available computing resources (the calculations for 24.5561 K took approximately 700 hours of CPU time with 2.2 GHz processors for each combination of two- and three-body potentials.)

4.2 Uncertainty Analysis

The largest contribution to the uncertainty of C(T) comes from imperfect knowledge of the three-body part of the potential energy, but there are smaller contributions that must also be considered. The following three uncertainty contributions are identified:

* Uncertainty in the pair potential. The pair potential [[empty set].sub.07] (see Sec. 2.1) used for our calculations has an unknown systematic error (difference from the true pair potential) that will produce a corresponding uncertainty in C(T). In Ref. (6), it is stated that the pair potentials [[empty set].sub.07.sup.+] and [[empty set].sub.07.sup.+] provide lower and upper bounds, respectively, to the true potential at a level of confidence corresponding to an expanded uncertainty at the k = 2 level. C(T) calculated from these perturbed potentials will represent an uncertainty at the same k = 2 level. Therefore, we take as the standard uncertainty (k= 1) from this source one-fourth of the difference between C(T) calculated with [[empty set].sub.07.sup.+] and with [[empty set].sub.07.sup.-].

* Uncertainty in the three-body potential. As explained in Sec. 2.2, we consider the CC+ and CC- potentials to bound the true three-body potential at a level of confidence corresponding to an expanded uncertainty at the k = 2 level. Therefore, we take as the standard uncertainty (k = 1) from this source one-fourth of the difference between C(T) calculated with CC+ and with CC-. At all temperatures studied, this is the largest of the three uncertainty contributions.

* Uncertainty in the convergence of the PIMC calculation. This is estimated as the standard deviation of the mean from 16 independent Monte Carlo samples, as described near the end of Sec. 3.2.2.

The first two of these contributions are systematic (Type B) errors, while the third is strictly statistical (Type A). The three contributions are combined in quadrature, and the resulting standard uncertainty is multiplied by a coverage factor k = 2 to produce the expanded uncertainty U(C) in Table 1.

For purposes of illustration, we review the uncertainty calculation (sometimes keeping insignificant digits for clarity) for the point at 273.16 K. The difference between C calculated with the [[empty set].sub.07.sup.+] and [[empty set].sub.07.sup.-] potentials is 0.062 [cm.sup.6] * [mol.sup.-2], so this component of the uncertainty is 0.0155 [cm.sup.6] * [mol.sup.-2]. The difference between calculations with the CC+ and CC- three-body potentials is 0.665 [cm.sup.6] * [mol.sup.-2], so this component of the standard uncertainty is 0.166 [cm.sup.6] * [mol.sup.-2]. The PIMC integration with the CC three-body potential yields C= 112.73 [cm.sup.6] * [mol.sup.-2] with integration uncertainty of 0.03 [cm.sup.6] * [mol.sup.-2]. Combining the three contributions in quadrature yields [u.sub.c] (C) = 0.1696 [cm.sup.6] * [mol.sup.-2], which when multiplied by two yields an expanded uncertainty with coverage factor k = 2 of U(C) = 0.339 [cm.sup.6] * [mol.sup.-2].

4.3 Correlation for Results

For convenience in interpolation, we have correlated the results for C(T) in Table 1 as a function of absolute temperature T with a simple empirical expression:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (37)

where T* = T/(100 K) and the parameters [a.sub.i] and [b.sub.i] are given in Table 2. Equation (37) reproduces all the C(T) values in Table 1 to within better than 0.11 %, which is smaller than the standard uncertainty [u.sub.c] (C) and similar to the portion of the uncertainty due to Type A sources of error. For temperatures above 2000 K, where our calculated values of C(T) are sparse, a finer grid of semiclassical calculations (which effectively coincide with the PIMC calculations at these temperatures; see Sec. 4.5) was used to guide the interpolation.

Equation (37) may be easily differentiated to obtain [dC/dT] and [[d.sup.2]C/d[T.sup.2]], which are useful in the context of acoustic measurements. We know of no rigorous way to estimate the uncertainty of these derivatives, but some idea of their quality may be obtained from the uncertainty in C(T) and the knowledge that the majority of that uncertainty comes from Type B (systematic) sources that would mostly cancel out in the computation of derivatives.

It is important to note that Eq. (37) is strictly for interpolation. It is only valid for the range of temperatures covered in Table 1 (24.5561 K to 10 000 K). The behavior of Eq. (37) outside this range is physically reasonable for short distances, but for example it does not reproduce the maximum in C(T) that is believed to exist near 4 K (35).

4.4 Comparison With Experiment

In this section, we compare our results with selected experimental data. This is not a comprehensive review of experimental work; the purpose of the comparisons is simply to show the scatter of the existing data and to demonstrate the improvement in uncertainty of C(T) obtained in our calculations. Therefore, we have selected for comparison the most recent data sets and those with well-documented uncertainties, adding other data sets in a few cases in order to cover the entire temperature range of the experimental data.

We begin with the near-ambient temperature range, which is of the most interest for metrology and has been the subject of the most experimental study. Figure 1 compares our results with four sources of experimental data (36-39). The error bars for the experimental points in this and subsequent figures represent an expanded uncertainty with coverage factor k = 2, corresponding approximately to a 95 % confidence interval. The appropriate comparable uncertainty corresponding to our calculated results is the last column of Table 1; these uncertainties are not shown on Fig. 1 because they are approximately the size of the symbols themselves.

[FIGURE 1 OMITTED]

These sources of experimental data are fairly consistent; in particular the data of Blancett et al. (36) and of McLinden and Losch-Will (39) have found use in metrology because their expanded uncertainties of approximately 3 [cm.sup.6] * [mol.sup.-2] to 5 [cm.sup.6] * [mol.sup.-2] have been the best available. Our results are fully consistent with these measurements, but have an uncertainty smaller by approximately one order of magnitude (for example, our U(C) at 273.16 K is 0.34 [cm.sup.6] * [mol.sup.-2]). Values of C(T) calculated classically are also shown in Fig. 1. Semiclassical results are not shown; they lie slightly above the full PIMC results but not enough to be clearly visible on the scale of the figure. Figures 2-4 show similar comparisons with selected data at low temperatures (40-43), moderately low temperatures (36), (37), (39), (41), (42), (44), and high temperatures (38), (45-47), respectively. In cases where experimental points are shown without error bars, the original paper did not report uncertainties. In all three figures, the uncertainty of the present calculations (see Table 1) is smaller than the size of the symbols. Classical results are shown on all three figures; only on the low-temperature Fig. 2 are the semiclassical results distinguishable from the PIMC calculations on the scale of the figure. For all these temperature ranges, our results are generally consistent with the available data. Given the scatter and uncertainties of the experimental data, it is reasonable to say that our calculations reduce the uncertainty in these temperature ranges by more than an order of magnitude. There appear to be no experimental data for C(T) available for helium at temperatures higher than those shown in Fig. 4. While estimates for these higher temperatures can be made based on extrapolation, Lennard-Jones potentials, etc., undoubtedly the present results are more reliable than any such estimates.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

After this work was completed, Gaiser and Fellmuth (48) published new experimental data for C(T) of helium from 3.7 K up to 26 K. Only our lowest-temperature point (24.5561 K) overlaps with their data. Interpolating their results to this temperature produces a value of approximately 281 [cm.sup.6] * [mol.sup.-2] with standard uncertainty u(C) of approximately 11 [cm.sup.6] * [mol.sup.-2]. This agrees within its uncertainty with our result given in Table 1.

4.5 Accuracy of Classical and Semiclassical C(T)

Because of the high computational demands of the PIMC calculation, it is natural to consider whether, at sufficiently high temperatures, adequate results may be obtained from the simpler semiclassical calculation or even the much simpler classical calculation. One might consider such a calculation "adequate" if it differed from the full PIMC calculation by less than the expanded uncertainty U(C) in Table 1.

The classical calculation significantly underestimates C(T) even at quite high temperatures. For example, at 1000 K it differs from the PIMC result by more than 0.6 [cm.sup.6] * [mol.sup.-2], more than twice the expanded uncertainty of the PIMC result. Only above approximately 2000 K is the classical calculation adequate in the sense defined above. The degree to which the classical calculation is in error at lower temperatures can be seen in Figs. 1-4. While this error is comparable to the scatter of the experimental data in many cases, it is much larger than the uncertainty in our calculations.

The semiclassical calculation is significantly better, closely reproducing the full PIMC calculations down to about 200 K and producing adequate results at temperatures as low as 120 K. In Fig. 5, we plot the deviation of the semiclassical calculation from Eq. (37) as a function of temperature. Figure 5 also shows the classical calculations and the individual PIMC calculations. Because systematic errors from the potential contribute similarly to both PIMC and semiclassical calculations, the error bars in Fig. 5 represent (at a k = 2 level) only the Type A uncertainty due to the convergence of the PIMC calculations.

[FIGURE 5 OMITTED]

4.6 Contribution of Three-Body Potential to C(T)

We can examine the contribution of the three-body potential to C(T). While this is the most important factor in the overall uncertainty (see Sec. 4.2). it actually contributes only a small fraction of the value of C(T) itself, on the order of 1 % or 2 %

Interestingly, the sign of the three-body contribution changes; it contributes positively to C(T) at low temperatures and negatively at moderate and high temperatures. This reflects the fact that dispersion forces (for which the net three-body contribution to the energy is generally positive) are more important at low temperatures, while exchange repulsion forces (for which the net three-body effect is generally negative) are more important at high temperatures. The "crossover" where the net three-body contribution to C(T) is zero occurs at approximately 170 K.

5. Discussion

We have computed C(T) for helium from state-of-the-art pair and three-body potentials, producing values with uncertainties smaller than those from experiment by at least an order of magnitude. For temperatures at and above the neon triple point, these values (as represented by Eq. (37)) provide a significant improvement for this useful quantity in metrology.

The largest source of uncertainty in our calculations of C(T) is the uncertainty in the three-body potential; therefore, this provides the greatest opportunity for improvement. As mentioned in Sec. 2.2, some very recent work (25), completed after our calculations were finished, has analyzed the three-body potential at the full-configuration-interaction (FCI) level of theory. This FCI potential is significantly more accurate than the CC and SAPT potentials developed earlier in Ref. (23). Given sufficient computing resources (and characterization of the uncertainty of the FCI three-body potential), our calculations here could be repeated with the FCI potential and smaller uncertainties obtained. This is planned for future work.

To give some indication of the improvement expected, in Table 3 we show a few values of [C.sub.FCI] calculated with the new three-body potential of Ref. (25). We are not yet in a position to assign uncertainties to [C.sub.FCI], but our preliminary estimate is that the uncertainties will be reduced by a factor of approximately three compared to our current values. It is clear from Table 3 that the values of [C.sub.FCI] are fully consistent with the values of C(T) computed in this work (so the uncertainty assigned here to the CC three-body potential of Ref. (23) was reasonable), but that our values are systematically lower by a small amount. Table 3 also shows values [C.sub.SAPT] calculated with the SAPT three-body potential of Ref. (23). These values are much further away from the accurate FCI results than those from the CC potential, demonstrating at the level of C(T) the superiority of the CC three-body potential, which the authors of Ref. (25) observed at the level of the potential itself. This provides further justification for our use of only the CC potential from Ref. (23) in our calculation of C(T).

Work is also underway to reduce the uncertainty in the pair potential (49). While improving the pair potential will significantly reduce the uncertainty in the pair quantities such as B(T) calculated by Hurly and Mehl (6), it will not be as helpful for C(T) where the pair potential contributes a relatively small portion of the overall uncertainty. When an improved pair potential (including characterization of its uncertainty) is completed, that would be an opportune time to recalculate C(T) with the FCI three-body potential of Ref. (25) in order to provide a complete and consistent set of state-of-the-art ab initio property values for helium.

It would, of course, be desirable to extend this work to lower temperatures. However, the required computing time per point scales approximately as 1/T, and the computation at 24.5561 K already took over 4000 CPU-hours. Extending the range much below 20 K would require an improvement in the technique or greatly increased computing power. At very low temperatures, around 7 K and below, (28) the use of Boltzmann statistics will begin to be in error. Extension of the PIMC method to incorporate the correct statistics for [.sup.4]He (a boson) is feasible but would introduce additional complexity and computing time.

There would also be some interest in C(T) for the isotope [.sup.3]He. Calculating C(T) for [.sup.3]He would be a straightforward extension of the current work, although the smaller mass (and therefore larger quantum effects) would somewhat increase the computational requirements. Unfortunately, the greatest interest in [.sup.3]He in metrology is for thermometry at very low temperatures that are currently beyond our reach with this method.

For metrology at moderately high pressures, the fourth virial coefficient D(T) might also be of interest. The potential energy should be adequately described by the sum of pair and three-body potentials; the relatively small size of the three-body contribution compared to the pair contribution suggests that nonadditivity of the potential at the four-body level should be tiny (this could be checked by ab initio calculations at selected configurations). PIMC calculation of D(T) would be computationally prohibitive, except perhaps at quite high temperatures where quantum effects are small. However, as discussed in Sec. 4.5, semiclassical calculations are accurate for C(T) above approximately 200 K; it is reasonable to assume that this would also be true for D(T). In the paper deriving the semiclassical perturbation approach to C(T), Yokota (10) indicates that extension of the approach to D(T) is feasible, but to our knowledge it has not been done. Such an extension of the semiclassical approach could provide highly accurate values of D(T) for helium in the important metrological range near room temperature.

Values of [dC/dT] and [d.sup.2]C/d[T.sup.2] are of interest for acoustic measurements, where they are related to the acoustic virial coefficients. These may be obtained by differentiating the interpolating Eq. (37), but the uncertainty in such derived values is difficult to quantify. Direct calculation of [dC/dT] and, possibly, [[d.sup.2]C]/[[d[T.sup.2]] might be feasible through the use of histogram reweighting together with the technique of thermodynamic integration to obtain N-particle functions and, as a consequence, the virial coefficients from Eqs. (3) and (4). Evaluation of the partition functions might also provide an alternative route for the calculation of D(T) and/or the incorporation of quantum statistics (bosonic or fermionic) at low temperatures.

Acknowledgments

We thank Mike Moldover for his encouragement and helpful advice during the course of this project. Karl Johnson for facilitating the collaboration that made this work possible, Greg Schenter and Richard Wheatley for discussions in the formative stages of this work, and Krzysztof Szalewicz for discussion of potentials and for making Refs. (23) and (25) and their computer code available prior to publication. We also thank Kevin Coakley for his assistance with the uncertainty analysis, and Eric Lemmon for assistance with data fitting.

The calculations have been performed on the HPC facility Wiglaf at the Physics Department of the University of Trento.

6. References

(1) J. W. Schmidt, R. M. Gavioso, E. F. May, and M. R. Moldover. Polarizability of helium and gas metrology, Phys. Rev. Lett., 98, 254504 (2007).

(2) B. Fellmuth, Ch. Gaiser, and J. Fischer. Determination of the Boltzmann constant--status and prospects, Meas. Sci. Technol. 17, R145 (2006).

(3) L. Pitre, M. R. Moldover, and W. L. Tew, Acoustic thermometry: new results from 273 K to 77 K and progress towards 4 K, Metrologia 43, 142 (2006).

(4) R. A. Aziz, A. R. Jansen, and M. R. Moldover, Ab Initio calculations for helium: A standard for transport property measurements, Phys. Rev. Lett. 74, 1586 (1995).

(5) J. J. Hurly and M. R. Moldover, Ab Initio values of the thermo-physical properties of helium as standards, J. Res. Natl. Inst. Stand. Technol. 105, 667 (2000).

(6) J. J. Hurly and J. B. Mehl. [.sup.4]He thermophysical properties: new ab initio calculations, J. Res. Natl. Inst. Stand. Technol. 112, 75 (2007).

(7) E. Bich, R. Hellmann, and E. Vogel, Ab initio potential energy curve for the helium atom pair and thermophysical properties of the dilute helium gas. 11. Thermophysical standard values for low-density helium, Mol. Phys. 105, 3035 (2007).

(8) R. Hellmann, E. Bich, and E. Vogel, Ab initio potential energy curve for the helium atom pair and thermophysical properties of the dilute helium gas. I. Helium-helium interatomic potential, Mol. Phys. 105, 3013 (2007).

(9) J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird, Molecular Theory of Gases and Liquids, John Wiley & Sons, New York, (1954).

(10) T. Yokota, Expansion theorem of density matrix, virial expansion and new formula for multiple scattering, J. Phys. Soc. Jap. 15,779(1960).

(11) J. Ram and Y. Singh, On the quantum corrections to the virial coefficients of the equation of state of a fluid, Mol. Phys. 26, 539(1973).

(12) T. D. Lee and C. N. Yang, Many-body problem in quantum statistical mechanics. II. Virial expansion for hard-sphere gas, Phys. Rev. 116,25(1959).

(13) A. Pais and G. E. Uhlenbeck, On the quantum theory of the third virial coefficient, Phys. Rev. 116, 250 (1959).

(14) S. Y. Larsen, Quantum-mechanical calculation of the third virial coefficient of [He.sup.4], Phys. Rev. 130, 1426 (1963).

(15) S. Y. Larsen and P. L. Mascheroni, Quantum-mechanical third virial coefficient and three-body phase shifts, Phys. Rev. A 2, 1018(1970).

(16) A. S. Reiner, Application of Faddeev techiques to the quantum theory of the third virial coefficient, Phys. Rev. 151, 170 (1966).

(17) L. D. Fosdick and H. F. Jordan, Path-integral calculation of the two-particle Slater sum for [He.sup.4], Phys. Rev. 143, 58 (1966).

(18) H. F. Jordan and L. D. Fosdick, Three-particle effects in the pair distribution function for [He.sup.4] gas, Phys. Rev. 171, 128 (1968).

(19) L. D. Fosdick, The Monte Carlo method in quantum statistics, SIAM Review 10, 315 (1968).

(20) P. Diep and J. K. Johnson, An accurate [H.sub.2]-[H.sub.2] interaction potential from first principles, J. Chem. Phys. 112, 4465 (2000). Erratum, ibid., 113, 3480 (2000).

(21) G. K. Schenter, The development of effective classical potentials and the quantum statistical mechanical second virial coefficient of water, J. Chem. Phys. 117, 6573 (2002).

(22) K. Patkowski, W. Cencek, P. Jankowski, K. Szalewicz, J. B. Mehl, G. Garberoglio, and A. H. Harvey, Potential energy surface for interactions between two hydrogen molecules, J. Chem. Phys. 129, 094304 (2008).

(23) W. Cencek, M. Jeziorska, O. Akin-Ojo, and K. Szalewicz. Three-body contribution to the helium interaction potential, J. Phys. Chem. A 111, 11311 (2007).

(24) M. Jeziorska, W. Cencek, K. Patkowski, B. Jeziorski, and K. Szaliewicz, Pair potential for helium from symmetry-adapted perturbation theory calculations and from supermolccular data, J. Chem. Phys. 127, 124303 (2007).

(25) W. Cencek, K. Patkowski. and K. Szalewicz, Full-configuration-interaction calculation of three-body nonadditive contribution to helium interaction potential, J. Chem. Phys. 131, 064105 (2009).

(26) B. M. Axilrod and E. Teller, Interaction of the van der Waals type between three atoms, J. Chem. Phys. 11, 299 (1943).

(27) M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, M. Booth, and F. Rossi, GNU Scientific Library Reference Manual, Network Theory, revised second edition (2006). http://www.gnu.org/software/gsl.

(28) M. E. Boyd, S. Y. Larsen, and J. E. Kilpatrick, Quantum mechanical second virial coefficient of a Lennard-Jones gas. Helium, J. Chem. Phys. 50, 4034 (1969).

(29) D. M. Ceperley, Path integrals in the theory of condensed helium, Rev. Mod. Phys. 67, 279 (1995).

(30) G. Garberoglio, Boltzmann bias grand canonical Monte Carlo, J. Chem. Phys. 128, 134109 (2008).

(31) G. P. Lepage, A new algorithm for adaptive multidimensional integration, J. Comp. Phys. 27, 192 (1978).

(32) G. P. Lepage, VEGAS: An adaptive multi-dimensional integration program, Technical report, Cornell preprint CLNS 80-447, 1980.

(33) P. Levy, Memorial des Sciences Mathematiques, Gauthier Villars, Paris, (1954). Fascicule 126.

(34) H. Preston-Thomas, The international temperature scale of 1990 (ITS-90), Metrologia 27, 3 (1990).

(35) P. P. M. Steur, M. Durieux, and G. T. McConville, Analytic expressions for the virial coefficients B(T) and C(T) of [.sup.4]He between 2.6 K and 300 K, Metrologia 24, 69 (1987).

(36) A. L. Blancett, K. R. Hall, and F. B. Canfield, Isotherms for the He-Ar systems at 50[degrees]C, 0[degrees]C and -50[degrees]C up to 700 atm, Physica 47, 75 (1970).

(37) D. D. Dillard, M. Waxman, and R. L. Robinson, Jr., Volumetric data and virial coefficients for helium, krypton, and helium-krypton mixtures, J. Chem. Eng. Data 23, 269 (1978).

(38) H. Mansoorian, K. R. Hall, J. C. Holste, and P. T. Eubank, The density of gaseous ethane and of fluid methyl chloride, and the vapor pressure of methyl chloride, J. Chem. Thermodyn. 13, 1001 (1981).

(39) M. O. McLinden and C. Losch-Will, Apparatus for wide-ranging, high-accuracy fluid (p, [rho], T) measurements based on a compact two-sinker densimeter, J. Chem. Thermodyn. 39, 507 (2007).

(40) D. White, T. Rubin, P. Camky, and H. L. Johnston, The virial coefficients of helium from 20 to 300[degrees]K, J. Phys. Chem. 64, 1607(1960).

(41) J. A. Sullivan and R. E. Sonntag, Compressibilities for helium at temperatures from 70 to 120[degrees]K and pressures to 690 atm, Cryogenics 7, 13 (1967).

(42) K. R. Hall and F. B. Canfield, Isotherms for the He-[N.sub.2] system at-190[degrees]C, -170[degrees]C, and-160[degrees]C up to 700 atm, Physica 47, 219 (1970).

(43) D. Gugan and G W. Michel, Measurements of the polarizabili-ty and of the second and third virial coefficients of He in the range 4.2 - 27.1 K, Mol. Phys. 39, 783 (1980).

(44) J. A. Provine and F. B. Canfield, Isotherms for the He-Ar system at -130, - 115, and -90[degrees]C up to 700 atm, Physica 52, 79 (1971).

(45) A. Michels and H. Wouters, Isotherms of helium between 0[degrees] and 150[degrees]C up to 200 amagat, Physica 8, 923 (1941).

(46) G S. Kell, G. E. McLaurin, and E. Whalley, Second virial coefficient of helium from 0 to 500[degrees]C by the two-temperature gas-expansion method, J. Chem. Phys. 68, 2199 (1978).

(47) M. R. Patel, L. L. Joffrion, and P. T. Eubank, A simple procedure for estimating virial coefficients from Burnett PVT data, AIChEJ., 34, 1229(1988).

(48) C. Gaiser and B. Fellmuth, Helium virial coefficients--a comparison between new highly accurate theoretical and experimental data, Metrologia 46, 525 (2009).

(49) K. Szalewicz, University of Delaware, 2009. Private communication.

About the authors: Giovanni Garberoglio is a physicist affiliated with CNISM (Consorzio Nazionale Interuniversitario per le Scienze fisiche della Materia) and the Physics Department at the University of Trento, Italy. Allan Harvey is a chemical engineer in the Theory and Modeling of Fluids Group of the Thermophysical Properties Division of the NIST Chemical Science and Technology Laboratory. The National Institute of Standards and Technology is an agency of the U.S. Department of Commerce.

Giovanni Garberoglio

CNISM and Dipartimento di Fisica,

Universita di Trento,

via Sommarive 14,

38100 Povo (TN), Italy

and

Allan H. Harvey

Thermophysical Properties Division,

National Institute of Standards and Technology,

325 Broadway,

Boulder CO 80305, USA

garbcrog@science.unitn.it

allan.harvey@nist.gov

Accepted: June 17, 2009

Available online: http://www.nist.gov/jres

Key words: calibration; density; helium; metrology; path integral; thermodynamic properties; virial coefficients.

1. Introduction

Accurate knowledge of the thermophysical properties of helium is desirable for many applications in metrology. Examples include microwave resonance measurements for development of a pressure standard (1), measuring the Boltzmann constant (2), and acoustic gas thermometry (3).

In many cases, the nonideal behavior of the gas is a major component of the overall uncertainty. That non-ideality is expressed in the virial expansion

[p/[[rho]RT]] = 1 + B(T)[rho] + C(T)[[rho].sup.2]+ ..., (1)

where p is the pressure, [rho] the molar density, R the molar gas constant, and T the absolute temperature. B(T) is the second virial coefficient, representing the lowest-order deviation from ideal-gas behavior. As the density increases, the contribution from the third virial coefficient C(T) becomes significant. B(T) is a function only of the interactions between pairs of molecules, while C(T) is determined by interactions among three molecules.

Because the helium atom has only two electrons, and because of advances in methodology, algorithms, and computing power for ab initio quantum calculations, it is now possible to calculate properties of individual helium atoms and pairs of atoms with very high accuracy. The use of such calculated values to develop standards for thermophysical properties was first proposed by Aziz et al. (4), after which Hurly and Moldover (5) calculated helium's second virial coefficient B(T), dilute-gas viscosity, and dilute-gas thermal conductivity with uncertainties smaller than those of the best experiments. Hurly and Mehl (6) recently improved on this work, primarily by using a more accurate pair potential. Bich et al. (7) performed similar calculations of pair quantities with another high-accuracy pair potential (8). They also calculated C(T) with quantum effects considered at a first-order perturbation level, but (as explained in Sec. 2.2) the Axilrod-Teller term they used for the three-body potential is not a good representation of three-body effects for helium.

For the third virial coefficient C(T), rigorous first-principles calculations are more difficult. While classical calculation of C(T) for spherically symmetric species is straightforward (9), quantum effects must be considered for a light gas such as helium. Semiclassical perturbation approaches have been derived (10), (11), but they will be insufficient at low temperatures and their uncertainty is difficult to quantify. What is needed is a fully quantum approach as is available for B(T) (9), but this problem has not been solved in closed form for C(T).

Early quantum calculations, aimed to elucidate how C(T) depends on the inter-molecular potentials, were performed at the beginning of the 1960s, using the general theoretical framework developed by Lee and Yang (12). Pais and Uhlenbeck (13) as well as Larsen (14) discussed various approximations to the full quantum problem, and these results were used to estimate the binding energy. Other calculation attempts tried to imitate the exact quantum calculations of B(T) based on scattering phase shifts. Larsen and Mascheroni (15) were able to obtain a rigorous expression for the third virial coefficient under the (unrealistic) assumption of the absence of bound states. Reiner (16) developed a calculation method based on the Faddeev equations for the quantum-mechanical three-body problem, but the difficulties in solving the coupled integral equations could not be overcome, and approximations were required.

The first rigorous quantum calculation of virial coefficients was developed by Fosdick and Jordan (17-19) who derived a path-integral expression for the third virial coefficient in the presence of two-body pairwise additive forces, and showed in detail how to evaluate it numerically in the case of a simple Lennard-Jones model for He.

Fosdick and Jordan's approach was independently rediscovered by Diep and Johnson (20) who devised, by analogy to the classical expression, a path-integral formula for the second virial coefficient of a quantum gas. They used their expression to calculate the second virial coefficient for a new [H.sub.2]-[H.sub.2] potential that they had computed using ab initio calculations, neglecting the rotational degrees of freedom. Their expression was generalized to the case of asymmetric rotors by Schenter (21) who used it to calculate the second virial coefficient of a model of water.

Recently, we developed a rigorous path-integral Monte Carlo procedure to calculate the second virial coefficient of molecular hydrogen, extending the approach pioneered by Fosdick and Jordan (17) to take into account the rotational degrees of freedom of linear molecules. We used this method together with state-of-the-art ab initio calculations of the pair potential to obtain good agreement with experimental data in a wide range of temperatures (22).

In the present work, we further extend the path-integral approach employed by Jordan and Fosdick for spherical particles (18) to calculate the third virial coefficient of [.sup.4]He in the presence of nonadditive three-body interactions. We use this method to calculate the third virial coefficient of [.sup.4]He from 24.5561 K to 10 000 K, using a recent ab initio derived three-body potential (23).

2. Intermolecular Potentials

2.1 Pair Potential

We write the pair potential as [U.sub.2](r), where r is the center-to-center distance between the atoms. We use the pair potential known as [[empty set].sub.07], which was developed by Hurly and Mehl (6) based on the best ab initio calculations available in 2007. For uncertainty analysis, we also use their potentials [[empty set].sub.07.sup.-] and [[empty set].sub.07.sup.+], which represent uncertainty limits (expanded uncertainty with coverage factor k = 2) for the [[empty set].sub.07] potential.

While this work was in progress, a pair potential of higher accuracy (consistent within mutual uncertainties with [[empty set].sub.07]) was published by Jeziorska et al. (24). Because our uncertainties are dominated by the three-body potential (see Sec. 4.2), our results would not have been significantly different had we used that pair potential.

2.2 Three-Body Potential

Calculation of C(T) also requires knowledge of the nonadditive three-body contribution to the potential energy in a triplet of atoms. We write this as [U.sub.3]([r.sub.12], [r.sub.13], [r.sub.23]), where [r.sub.ij] is the distance between atoms i and j. The three-body potential of helium has been studied by Cencek et al. (23), who developed separate ab initio potentials based on symmetry-adapted perturbation theory (SAPT) and coupled-cluster (CC) calculations. The two potentials were estimated to have a maximum uncertainty of 10 %, but very recent work [25] has shown that the SAPT potential is significantly less accurate, while the CC potential is in good agreement (better than 10 %) with higher-level calculations. We therefore employ the CC potential in this work. For purposes of uncertainty analysis, we interpret the estimated 10 % uncertainty as an expanded uncertainty at the k = 2 level. In analogy to the procedure of Hurly and Mehl [6] for the pair potential, we define for our uncertainty analysis three-body potentials CC- (obtained by multiplying the potential by 1.1 where it is negative and by 0.9 where it is positive) and CC+ (multiplying by 0.9 where it is negative and by 1.1 where it is positive).

It is worth noting that, at most temperatures of interest, the main three-body effect for helium is the non-additivity in repulsive forces. The commonly used Axilrod-Teller three-body term [26] accounts only for dispersion interactions. At all but the lowest temperatures considered in this work, the reduced temperature is quite high (compared to the characteristic temperature for the pair dispersion interaction of approximately 11 K), and three-body dispersion effects are less important than three-body repulsion effects. Therefore, an Axilrod-Teller term would seriously misrepresent the three-body potential, even giving a contribution to C(T) of the wrong sign at moderate and high temperatures.

3. Calculation Methods

Let us denote by [Q.sub.N](V,T) the partition function of N particles in a volume V at a temperature T. By defining the quantities [Z.sub.N] as

[[Z.sub.N]/[N!]][equivalent to][[[Q.sub.N](V, T)[V.sup.N]]/[[Q.sub.1][(V, T).sup.N]]], (2)

then the expressions for the second and third virial coefficients, B(T) and C(T), become [9]

B(T) = -[1/[2V]]([Z.sub.2] - [Z.sub.1.sup.2]), (3)

C(T) = 4[B.sup.2](T) - 1/3V([Z.sub.3] - 3[Z.sub.2][Z.sub.1] + 2[Z.sub.1.sup.3]). (4)

3.1 Classical and Semiclassical Calculations

The explicit expression for the partition functions and hence the coefficients [Z.sub.k] appearing in Eqs. (3) and (4) depends on the framework in which the calculations are performed. In classical statistical mechanics (including the correct Boltzmann counting) one has

[Z.sub.class.sup.N] = [integral][d.sup.3][x.sub.1] ... [d.sup.3][x.sub.N]exp( - [beta]U([x.sub.1], ...,[x.sub.N])), (5)

where [beta] = [1/[k.sub.B]]T, [k.sub.B] denotes the Boltzmann constant and U([x.sub.1], ..., [x.sub.N]) is the total potential energy of a configuration with N particles at the positions [x.sub.1], ..., [x.sub.N].

In this case, the third virial coefficient [C.sup.class] (T) is given by the sum of a term depending on the two-body potential [U.sub.2](r), [C.sub.[2-body].sup.class] (T), and a term depending on the three-body potential, [C.sub.[3-body].sup.class], given by

[C.sup.class](T) = [C.sub.[2-body].sup.class](T) + [C.sub.[3-body].sup.class](T), (6)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

We have denoted by [theta] the angle between the vectors [r.sub.12] and [r.sub.13]. The distance between the particles labeled 2 and 3 is therefore given by

[r.sub.23] = [square root of [r.sub.12.sup.2] + [r.sub.13.sup.2] - 2[r.sub.12][r.sub.13]cos[theta]].(9)

The classical formulae are accurate enough for heavy particles at high temperature. If this limit is not attained, quantum diffraction effects (Heisenberg uncertainty) become appreciable, as is the case for hydrogen and helium at and below room temperature. As long as the quantum effects can be considered a small correction to the classical behavior, the expressions given above can be corrected by including the first term in the expansion of the full quantum expression in even powers of h.

The expression for the first quantum correction to the third virial coefficient has been evaluated by Yokota (10). Setting [a.sup.2] = [h.sup.2][beta]/m, where m is the mass of the particles under consideration, the semiclassical expression for the third virial coefficient turns out to be

[C.sup.semi](T) = [C.sup.class](T) + [a.sup.2][C.sub.1](T), (10)

[C.sub.1](T) = 8[B.sub.1](T)[B.sup.class](T) - 2V[B.sub.1](T) + [[[beta].sup.2]/36][integral][d.sup.3][r.sub.12][d.sup.3][r.sub.13][e.sup[beta]U](KUU) (11)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

KUU = [[nabla].sub.12]U * [[nabla].sub.12]U + [[nabla].sub.13]U * [[nabla].sub.13]U + [1/2]([[nabla].sub.12]U * [[nabla].sub.13]U + [[nabla].sub.13]U * [[nabla].sub.12]U). (13)

The classical and semiclassical values of C(T) have been obtained by direct numerical integration of Eqs. (6) and (10). We have used the QAG adaptive algorithm together with the Gauss-Kronrod 15-point rule as implemented in the GNU Scientific Library (27). The interaction has been neglected beyond a cutoff length of L = 4 nm. We have checked that using a larger cutoff value does not appreciably affect the results, and the same cutoff was also used in performing the path-integral Monte Carlo calculations described below.

3.2 Path-Integral Monte Carlo

In the framework of non-relativistic quantum statistical mechanics, the expression for the quantity [Z.sub.N] of Eq. (2) becomes

[Z.sub.N] = [[LAMBDA].sub.m.sup.3N][summation over k][summation over [pi]]<k | exp(-[beta][[^.H].sub.N])[[^.P].sub.[pi]] | k>, (14)

where [[^.H].sub.N] is the N-body Hamiltonian operator, |k> denotes a complete set of N-particle states and [[^.P].sub.[pi]] is a permutation operator, with the proper sign to take into account the bosonic or fermionic nature of the particles. [[LAMBDA].sub.m] = [h/[square root of [2[pi]m[k.sub.B]T]] denotes the thermal de Broglie wavelength for a particle of mass m. In the following, we will not be concerned with temperatures so low that exchange effects play a relevant role [28] and hence we will approximate the sum over [[^.P].sub.[pi]] with the identity operator (Boltzmann statistics).

Equation (14) is the starting point to derive the path-integral expressions for the second and third virial coefficients, using Eqs. (3) and (4). In order to avoid cumbersome notation, we will present the derivation in detail in the case of the second virial coefficient. This will allow us to establish some useful notation that will be used to present in the most compact form possible the path-integral formulae for the third virial coefficient.

3.2.1 Second Virial Coefficient

The path-integral formula for the second virial coefficient of Eq. (3), using the quantum-mechanical expression of Eq. (14) in the case of Boltzmann statistics, is readily obtained by first performing a canonical transformation to the center of mass R [equivalent to] [([x.sub.1] + [x.sub.2)]/2] and relative r [equivalent to] [r.sub.1] [equivalent to] [x.sub.2] - [x.sub.1] coordinates. In the equation for [Z.sub.2], the kinetic energy relative to the center-of-mass motion commutes with the kinetic and potential energy of the relative motion and can be integrated out, obtaining a factor [V/[[LAMBDA].sub.M.sup.3]], where M= 2m. As a consequence, B(T) is proportional to the trace of the Hamiltonian describing the relative motion [[^.H].sub.r] = [[^.p].sub.r.sup.2]/2[micro]] + [[^.U].sub.2] (|r|) [equivalent to] [^.T] + [[^.U].sub.2], where [mu] = [m/2] is the reduced mass of the pair and [[^.p].sub.r] is the momentum conjugated to the relative coordinate r.

The resulting expression [Z.sub.2] = [[[LAMBDA].sub.m.sup.6]/[2[[LAMBDA].sub.M.sup.3]]][integral][d.sup.3]r<r | exp(-[beta][[^.H].sub.r]) | r>, (15)

can then be evaluated by performing a Trotter expansion of the kinetic and potential energies of the relative motion,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

keeping a finite value of the Trotter index P and inserting P--1 completeness relations 1 = [integral] [dr.sub.i]|[r.sub.i] > <[r.sub.i]| (i = 2, ...,P) between each of the P factors in Eq. (16). The matrix elements of the kinetic energy operator can be evaluated explicitly, obtaining

<[r.sub.i] | exp(-[[beta][[^.p].sup.2]/[2[micro]P]]) | [r.sub.j]> = [[P.sup.3/2]/[[LAMBDA].sub.[mu].sup.3]]exp(-[K/2][([r.sub.i] - [r.sub.j]).sup.2]), (17)

where K = [2[pi]P/[[LAMBDA].sup.[mu].sub.2]]. The final outcome of this chain of equivalences is to map the calculation of the quantum partition function of Eq. (15) to the calculation of the classical partition function of ring polymers with P beads each (29). The resulting expression can be simplified by introducing the coordinates r = [r.sub.1], [DELTA][r.sub.i] = [r.sub.[i+1]] - [r.sub.i] (i = 1, ..., P--1) and letting

[[bar.U].sub.2](r) = [1/P][P.summation over (i = 1)][U.sub.2]([r.sub.i]), (18)

indicate the average of the two-body potential over the positions occupied by the P beads of a given ring polymer. Performing analogous manipulations in the rather trivial case of [Z.sup.2.sub.1] produces the following expression for the second virial coefficient:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

where [F.sub.ring] is the probability of finding a ring polymer configuration in the ideal gas phase, as shown in Ref. (30), and is given by

[F.sub.ring]([mu];[DELTA][r.sub.1], ...,[DELTA][r.sub.[P-1]]) = [[LAMBDA].sub.[mu].sup.3][([[P.sup.3/2]/[[LAMBDA].sub.[mu].sup.3]]).sup.P] x exp[-[K/2][P.summation over (i = 1)][DELTA][r.sub.i.sup.2]], (20)

where [DELTA] [r.sub.p] = [r.sub.p] - [r.sub.1].

The second virial coefficient can also be written as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

where the effective potential is defined as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (22)

that is, by averaging the factor exp (- [beta] [[bar.U].sub.2]) over ring polymer configurations drawn from an ideal gas distribution and having one bead at the point r. In the classical limit h [right arrow] 0, the ring polymers shrink to a point, and hence Eq. (21) recovers the classical result.

3.2.2 Third Virial Coefficient

The same reasoning that was followed to derive the path-integral expression for the second virial coefficient B(T) in Eq. (19) can be applied to the case of the third virial coefficient. In this case it is useful to evaluate the expectation values over three-body operators after having performed a canonical transformation to the Jacobi coordinates R, r, [rho] and the corresponding momenta P, p, [pi], defined as:

R = [1/3]([x.sub.1] + [x.sub.2] + [x.sub.3]); P = [p.sub.1] + [p.sub.2] + [p.sub.3] (23)

r = [x.sub.2] - [x.sub.1]; p = [1/2]([p.sub.2] - [p.sub.1]) (24)

[rho] = [x.sub.3] - [1/2]([x.sub.1] + [x.sub.2]); [pi] = [1/3](2[p.sub.3] - [p.sub.1] - [p.sub.2]). (25)

As in the case of B(T) the center-of-mass motion can be integrated out, but the Trotter factorization introduces two different ring polymers, corresponding to the coordinates r and [rho]; we note in passing that the masses associated with these degrees of freedom are [M.sub.r] = m/2 and [M.sub.[rho]] = 2m/3, respectively, and that the total potential energy of a three-body configuration,

U([x.sub.1], [x.sub.2], [x.sub.3])[equivalent to][U.sub.3](|[x.sub.2] - [x.sub.1]|,|[x.sub.3] - [x.sub.2]|,|[x.sub.3] - [x.sub.1]|) + [3.summation over (i< j = 1)][U.sub.2](|[x.sub.i] - [x.sub.j]|), (26)

is a function of the coordinates r and [rho] only. As a final result, the expression for C(T) can be written as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

where we have defined, analogously to Eq. (18),

[bar.U](r; [rho]) = [1/P][P.summation over (i = 1)]U([r.sub.i]; [[rho].sub.i]), (28)

as the average of the potential energy of the three bodies over the positions indicated by the beads of the two given ring polymers corresponding to the Jacobi coordinates r = [r.sub.1] and [rho] = [[rho].sub.1]. The three exponentials of [[bar.U].sub.2] appearing in Eq. (27) come from the three terms [Z.sub.2] [Z.sub.1.sup.2] of Eq. (4), when the two-body integral is written as a function of the coordinates

[r.sub.21] = [x.sub.2] - [x.sub.1] = r, (29)

[r.sub.31] = [x.sub.3] - [x.sub.1] = [rho] + [r/2], (30)

[r.sub.32] = [x.sub.3] - [x.sub.2] = [rho] - [r/2], (31)

respectively.

The integrals over the variables [DELTA][r.sub.i] and [DELTA][[rho].sub.i] allow one to define effective two- and three-body potentials as averages over the two kinds of ring polymers, analogously to Eq. (22)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

where x can be either r or [rho] + r/2 or [rho]-r/2. The final expression for C(T) is hence

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (34)

Given the rotational symmetry of the system, the volume of integration can be written as

[d.sup.3][rd.sup.3][rho] = 4[pi][r.sup.2]dr2[pi][[rho].sup.2]d[rho]d cos[THETA], (35)

where r and [rho] are the moduli of the vectors r and [rho], respectively, while [THETA] is the angle between them.

In the actual calculation, it is useful to write the square of the second virial coefficient as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (36)

so that the difference of two integrals in Eq. (34) can be evaluated as the integral of the difference. Note that in Eq. (36) the coordinate [rho] is associated with a particle having mass [M.sub.r] and not [M.sub.p] as in Eqs. (32) and (33).

The calculation of the third virial coefficient in the path-integral formalism follows directly from Eq. (34), which shows that C(T) is given as a three-dimensional integral using effective two- and three-body potentials, given by Eqs. (33) and (32), respectively. Since for each of the atomic positions the effective potentials are obtained as a costly average over configurations of ring polymers, we chose to evaluate Eq. (34) using a Monte Carlo integration procedure, namely the VEGAS algorithm (27), (31), (32). We used N = 5 x [10.sup.5] integration points for the production stage of the algorithm and half as many for the equilibration stage.

For each of the atomic configurations considered in the course of the Monte Carlo integration, we generate n = 200 ring polymers for the r and [rho] coordinates, distributed according to the probability [F.sub.ring]. This can be done very efficiently using an interpolation formula due to Levy (18), (33). For each of the ring polymers, we evaluate the corresponding average potentials [[bar.U].sub.2] and [bar.U] and accumulate their Boltzmann factors to calculate the effective potentials for the given configuration according to Eqs. (32) and (33).

In order to estimate the statistical uncertainty of the values of the third virial coefficient so obtained, we perform 16 of these calculations for each value of T, starting with different seeds for the random number generator.

The final value for C(T) is obtained as the average of the values coming from the 16 independent calculations, and its Type A uncertainty is estimated as the standard error of the mean from the same set of values. Notice that this uncertainty takes into account the statistical error resulting from use of both a finite N and a finite n in the calculations.

Finally, let us discuss the choice of the Trotter index P. Since the path-integral method is exact in the limit of large P, we fixed this number by calculating B(T) with our method for progressively increasing values of the Trotter index P, until our results matched those performed with the phase-shift method (6). We found that the choice P = 7 + 2400 K/T was enough to reach convergence in a wide range of temperatures, and we similarly checked that the same value of P was sufficient in the case of C(T) at 273.16 K. The offset of 7 in the expression for P is due to the approximation inherent in Levy's interpolation formula (18), (33).

4. Results

4.1 Third Virial Coefficients

We calculated C(T) as described in Sec. 3.2 for the potential-energy surface obtained by combining the [[empty set].sub.07] pair potential referenced in Sec. 2.1 with the three-body potential (CC) referenced in Sec. 2.2. Calculations were performed over a wide range of temperatures, as shown in Table 1. In addition to round values, some temperatures were chosen due to their importance in metrology. For example, our lowest temperature of 24.5561 K is the value assigned to the triple point of neon in the International Temperature Scale of 1990 (ITS-90) (34). It should be noted that T in our calculations is the thermodynamic temperature, which may differ on the order of 0.005 % from the corresponding ITS-90 temperature (3). We included temperatures up to 10 000 K, in order to match the range covered by Hurly and Mehl (6) for B(T). Temperatures below those shown in Table 1 were not achievable with available computing resources (the calculations for 24.5561 K took approximately 700 hours of CPU time with 2.2 GHz processors for each combination of two- and three-body potentials.)

Table 1. Third virial coefficients C(T) calculated in this work and our estimates (see Sec. 4.2) of their expanded (k = 2) uncertainties U(C) T(K) C([cm.sup.6] * [mol.sup.-2]) U(C) ([cm.sup.6] * [mol.sup.-2]) 24.5561 273.34 1.26 30.0 248.08 1.54 40.0 221.36 0.96 50.0 205.31 0.78 63.15 190.38 0.83 75.0 180.56 0.64 83.806 174.37 0.63 100.0 164.65 0.55 120.0 154.78 0.47 140.0 146.59 0.45 170.0 136.58 0.39 200.0 128.34 0.33 223.152 122.80 0.33 235.0 120.18 0.37 250.0 117.15 0.29 273.16 112.73 0.34 302.915 107.75 0.31 350.0 100.91 0.30 400.0 94.76 0.31 429.75 91.48 0.27 450.0 89.43 0.26 500.0 84.87 0.26 550.0 80.81 0.25 600.0 77.18 0.24 650.0 73.96 0.26 700.0 71.00 0.23 750.0 68.33 0.24 800.0 65.92 0.23 900.0 61.61 0.24 1000.0 57.87 0.24 1200.0 51.80 0.22 1400.0 46.98 0.23 1600.0 43.08 0.22 1800.0 39.81 0.22 2000.0 37.03 0.21 2500.0 31.59 0.21 5000.0 18.31 0.19 10000.0 9.67 0.16

4.2 Uncertainty Analysis

The largest contribution to the uncertainty of C(T) comes from imperfect knowledge of the three-body part of the potential energy, but there are smaller contributions that must also be considered. The following three uncertainty contributions are identified:

* Uncertainty in the pair potential. The pair potential [[empty set].sub.07] (see Sec. 2.1) used for our calculations has an unknown systematic error (difference from the true pair potential) that will produce a corresponding uncertainty in C(T). In Ref. (6), it is stated that the pair potentials [[empty set].sub.07.sup.+] and [[empty set].sub.07.sup.+] provide lower and upper bounds, respectively, to the true potential at a level of confidence corresponding to an expanded uncertainty at the k = 2 level. C(T) calculated from these perturbed potentials will represent an uncertainty at the same k = 2 level. Therefore, we take as the standard uncertainty (k= 1) from this source one-fourth of the difference between C(T) calculated with [[empty set].sub.07.sup.+] and with [[empty set].sub.07.sup.-].

* Uncertainty in the three-body potential. As explained in Sec. 2.2, we consider the CC+ and CC- potentials to bound the true three-body potential at a level of confidence corresponding to an expanded uncertainty at the k = 2 level. Therefore, we take as the standard uncertainty (k = 1) from this source one-fourth of the difference between C(T) calculated with CC+ and with CC-. At all temperatures studied, this is the largest of the three uncertainty contributions.

* Uncertainty in the convergence of the PIMC calculation. This is estimated as the standard deviation of the mean from 16 independent Monte Carlo samples, as described near the end of Sec. 3.2.2.

The first two of these contributions are systematic (Type B) errors, while the third is strictly statistical (Type A). The three contributions are combined in quadrature, and the resulting standard uncertainty is multiplied by a coverage factor k = 2 to produce the expanded uncertainty U(C) in Table 1.

For purposes of illustration, we review the uncertainty calculation (sometimes keeping insignificant digits for clarity) for the point at 273.16 K. The difference between C calculated with the [[empty set].sub.07.sup.+] and [[empty set].sub.07.sup.-] potentials is 0.062 [cm.sup.6] * [mol.sup.-2], so this component of the uncertainty is 0.0155 [cm.sup.6] * [mol.sup.-2]. The difference between calculations with the CC+ and CC- three-body potentials is 0.665 [cm.sup.6] * [mol.sup.-2], so this component of the standard uncertainty is 0.166 [cm.sup.6] * [mol.sup.-2]. The PIMC integration with the CC three-body potential yields C= 112.73 [cm.sup.6] * [mol.sup.-2] with integration uncertainty of 0.03 [cm.sup.6] * [mol.sup.-2]. Combining the three contributions in quadrature yields [u.sub.c] (C) = 0.1696 [cm.sup.6] * [mol.sup.-2], which when multiplied by two yields an expanded uncertainty with coverage factor k = 2 of U(C) = 0.339 [cm.sup.6] * [mol.sup.-2].

4.3 Correlation for Results

For convenience in interpolation, we have correlated the results for C(T) in Table 1 as a function of absolute temperature T with a simple empirical expression:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (37)

where T* = T/(100 K) and the parameters [a.sub.i] and [b.sub.i] are given in Table 2. Equation (37) reproduces all the C(T) values in Table 1 to within better than 0.11 %, which is smaller than the standard uncertainty [u.sub.c] (C) and similar to the portion of the uncertainty due to Type A sources of error. For temperatures above 2000 K, where our calculated values of C(T) are sparse, a finer grid of semiclassical calculations (which effectively coincide with the PIMC calculations at these temperatures; see Sec. 4.5) was used to guide the interpolation.

Table 2. Coefficients for Eq. (37) for the third virial coefficient of helium i [a.sub.i] [b.sub.i] 1 -13 337.07 -0.77 2 36 155.73 -0.85 3 -50 678.58 -1.00 4 50 673.92 -1.15 5 -23 876.30 -1.25 6 1 226.921 -1.50

Equation (37) may be easily differentiated to obtain [dC/dT] and [[d.sup.2]C/d[T.sup.2]], which are useful in the context of acoustic measurements. We know of no rigorous way to estimate the uncertainty of these derivatives, but some idea of their quality may be obtained from the uncertainty in C(T) and the knowledge that the majority of that uncertainty comes from Type B (systematic) sources that would mostly cancel out in the computation of derivatives.

It is important to note that Eq. (37) is strictly for interpolation. It is only valid for the range of temperatures covered in Table 1 (24.5561 K to 10 000 K). The behavior of Eq. (37) outside this range is physically reasonable for short distances, but for example it does not reproduce the maximum in C(T) that is believed to exist near 4 K (35).

4.4 Comparison With Experiment

In this section, we compare our results with selected experimental data. This is not a comprehensive review of experimental work; the purpose of the comparisons is simply to show the scatter of the existing data and to demonstrate the improvement in uncertainty of C(T) obtained in our calculations. Therefore, we have selected for comparison the most recent data sets and those with well-documented uncertainties, adding other data sets in a few cases in order to cover the entire temperature range of the experimental data.

We begin with the near-ambient temperature range, which is of the most interest for metrology and has been the subject of the most experimental study. Figure 1 compares our results with four sources of experimental data (36-39). The error bars for the experimental points in this and subsequent figures represent an expanded uncertainty with coverage factor k = 2, corresponding approximately to a 95 % confidence interval. The appropriate comparable uncertainty corresponding to our calculated results is the last column of Table 1; these uncertainties are not shown on Fig. 1 because they are approximately the size of the symbols themselves.

[FIGURE 1 OMITTED]

These sources of experimental data are fairly consistent; in particular the data of Blancett et al. (36) and of McLinden and Losch-Will (39) have found use in metrology because their expanded uncertainties of approximately 3 [cm.sup.6] * [mol.sup.-2] to 5 [cm.sup.6] * [mol.sup.-2] have been the best available. Our results are fully consistent with these measurements, but have an uncertainty smaller by approximately one order of magnitude (for example, our U(C) at 273.16 K is 0.34 [cm.sup.6] * [mol.sup.-2]). Values of C(T) calculated classically are also shown in Fig. 1. Semiclassical results are not shown; they lie slightly above the full PIMC results but not enough to be clearly visible on the scale of the figure. Figures 2-4 show similar comparisons with selected data at low temperatures (40-43), moderately low temperatures (36), (37), (39), (41), (42), (44), and high temperatures (38), (45-47), respectively. In cases where experimental points are shown without error bars, the original paper did not report uncertainties. In all three figures, the uncertainty of the present calculations (see Table 1) is smaller than the size of the symbols. Classical results are shown on all three figures; only on the low-temperature Fig. 2 are the semiclassical results distinguishable from the PIMC calculations on the scale of the figure. For all these temperature ranges, our results are generally consistent with the available data. Given the scatter and uncertainties of the experimental data, it is reasonable to say that our calculations reduce the uncertainty in these temperature ranges by more than an order of magnitude. There appear to be no experimental data for C(T) available for helium at temperatures higher than those shown in Fig. 4. While estimates for these higher temperatures can be made based on extrapolation, Lennard-Jones potentials, etc., undoubtedly the present results are more reliable than any such estimates.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

After this work was completed, Gaiser and Fellmuth (48) published new experimental data for C(T) of helium from 3.7 K up to 26 K. Only our lowest-temperature point (24.5561 K) overlaps with their data. Interpolating their results to this temperature produces a value of approximately 281 [cm.sup.6] * [mol.sup.-2] with standard uncertainty u(C) of approximately 11 [cm.sup.6] * [mol.sup.-2]. This agrees within its uncertainty with our result given in Table 1.

4.5 Accuracy of Classical and Semiclassical C(T)

Because of the high computational demands of the PIMC calculation, it is natural to consider whether, at sufficiently high temperatures, adequate results may be obtained from the simpler semiclassical calculation or even the much simpler classical calculation. One might consider such a calculation "adequate" if it differed from the full PIMC calculation by less than the expanded uncertainty U(C) in Table 1.

The classical calculation significantly underestimates C(T) even at quite high temperatures. For example, at 1000 K it differs from the PIMC result by more than 0.6 [cm.sup.6] * [mol.sup.-2], more than twice the expanded uncertainty of the PIMC result. Only above approximately 2000 K is the classical calculation adequate in the sense defined above. The degree to which the classical calculation is in error at lower temperatures can be seen in Figs. 1-4. While this error is comparable to the scatter of the experimental data in many cases, it is much larger than the uncertainty in our calculations.

The semiclassical calculation is significantly better, closely reproducing the full PIMC calculations down to about 200 K and producing adequate results at temperatures as low as 120 K. In Fig. 5, we plot the deviation of the semiclassical calculation from Eq. (37) as a function of temperature. Figure 5 also shows the classical calculations and the individual PIMC calculations. Because systematic errors from the potential contribute similarly to both PIMC and semiclassical calculations, the error bars in Fig. 5 represent (at a k = 2 level) only the Type A uncertainty due to the convergence of the PIMC calculations.

[FIGURE 5 OMITTED]

4.6 Contribution of Three-Body Potential to C(T)

We can examine the contribution of the three-body potential to C(T). While this is the most important factor in the overall uncertainty (see Sec. 4.2). it actually contributes only a small fraction of the value of C(T) itself, on the order of 1 % or 2 %

Interestingly, the sign of the three-body contribution changes; it contributes positively to C(T) at low temperatures and negatively at moderate and high temperatures. This reflects the fact that dispersion forces (for which the net three-body contribution to the energy is generally positive) are more important at low temperatures, while exchange repulsion forces (for which the net three-body effect is generally negative) are more important at high temperatures. The "crossover" where the net three-body contribution to C(T) is zero occurs at approximately 170 K.

5. Discussion

We have computed C(T) for helium from state-of-the-art pair and three-body potentials, producing values with uncertainties smaller than those from experiment by at least an order of magnitude. For temperatures at and above the neon triple point, these values (as represented by Eq. (37)) provide a significant improvement for this useful quantity in metrology.

The largest source of uncertainty in our calculations of C(T) is the uncertainty in the three-body potential; therefore, this provides the greatest opportunity for improvement. As mentioned in Sec. 2.2, some very recent work (25), completed after our calculations were finished, has analyzed the three-body potential at the full-configuration-interaction (FCI) level of theory. This FCI potential is significantly more accurate than the CC and SAPT potentials developed earlier in Ref. (23). Given sufficient computing resources (and characterization of the uncertainty of the FCI three-body potential), our calculations here could be repeated with the FCI potential and smaller uncertainties obtained. This is planned for future work.

To give some indication of the improvement expected, in Table 3 we show a few values of [C.sub.FCI] calculated with the new three-body potential of Ref. (25). We are not yet in a position to assign uncertainties to [C.sub.FCI], but our preliminary estimate is that the uncertainties will be reduced by a factor of approximately three compared to our current values. It is clear from Table 3 that the values of [C.sub.FCI] are fully consistent with the values of C(T) computed in this work (so the uncertainty assigned here to the CC three-body potential of Ref. (23) was reasonable), but that our values are systematically lower by a small amount. Table 3 also shows values [C.sub.SAPT] calculated with the SAPT three-body potential of Ref. (23). These values are much further away from the accurate FCI results than those from the CC potential, demonstrating at the level of C(T) the superiority of the CC three-body potential, which the authors of Ref. (25) observed at the level of the potential itself. This provides further justification for our use of only the CC potential from Ref. (23) in our calculation of C(T).

Table 3. Selected values of third virial coefficients C(T) calculated in this work from the CC three-body potential (23) and their expanded (k = 2) uncertainties U(C), along with values [C.sub.SAPT] calculated from the SAPT three-body potential (23) and [C.sub.FCI] calculated from the new three-body potential of Cencek et al. (25) T (K) C U(C) [C.sub.SAPT] [C.sub.FCI] ([cm.sup.6] * ([cm.sup.6] * ([cm.sup.6] * ([cm.sup.6] * [mol.sup.-2]) [mol.sup.-2]) [mol.sup.-2]) [mol.sup.-2]) 100.0 164.65 0.55 163.98 164.98 200.0 128.34 0.33 127.89 128.60 273.16 112.73 0.34 112.34 112.92 400.0 94.76 0.31 94.41 94.88 1000.0 57.87 0.24 57.54 57.97

Work is also underway to reduce the uncertainty in the pair potential (49). While improving the pair potential will significantly reduce the uncertainty in the pair quantities such as B(T) calculated by Hurly and Mehl (6), it will not be as helpful for C(T) where the pair potential contributes a relatively small portion of the overall uncertainty. When an improved pair potential (including characterization of its uncertainty) is completed, that would be an opportune time to recalculate C(T) with the FCI three-body potential of Ref. (25) in order to provide a complete and consistent set of state-of-the-art ab initio property values for helium.

It would, of course, be desirable to extend this work to lower temperatures. However, the required computing time per point scales approximately as 1/T, and the computation at 24.5561 K already took over 4000 CPU-hours. Extending the range much below 20 K would require an improvement in the technique or greatly increased computing power. At very low temperatures, around 7 K and below, (28) the use of Boltzmann statistics will begin to be in error. Extension of the PIMC method to incorporate the correct statistics for [.sup.4]He (a boson) is feasible but would introduce additional complexity and computing time.

There would also be some interest in C(T) for the isotope [.sup.3]He. Calculating C(T) for [.sup.3]He would be a straightforward extension of the current work, although the smaller mass (and therefore larger quantum effects) would somewhat increase the computational requirements. Unfortunately, the greatest interest in [.sup.3]He in metrology is for thermometry at very low temperatures that are currently beyond our reach with this method.

For metrology at moderately high pressures, the fourth virial coefficient D(T) might also be of interest. The potential energy should be adequately described by the sum of pair and three-body potentials; the relatively small size of the three-body contribution compared to the pair contribution suggests that nonadditivity of the potential at the four-body level should be tiny (this could be checked by ab initio calculations at selected configurations). PIMC calculation of D(T) would be computationally prohibitive, except perhaps at quite high temperatures where quantum effects are small. However, as discussed in Sec. 4.5, semiclassical calculations are accurate for C(T) above approximately 200 K; it is reasonable to assume that this would also be true for D(T). In the paper deriving the semiclassical perturbation approach to C(T), Yokota (10) indicates that extension of the approach to D(T) is feasible, but to our knowledge it has not been done. Such an extension of the semiclassical approach could provide highly accurate values of D(T) for helium in the important metrological range near room temperature.

Values of [dC/dT] and [d.sup.2]C/d[T.sup.2] are of interest for acoustic measurements, where they are related to the acoustic virial coefficients. These may be obtained by differentiating the interpolating Eq. (37), but the uncertainty in such derived values is difficult to quantify. Direct calculation of [dC/dT] and, possibly, [[d.sup.2]C]/[[d[T.sup.2]] might be feasible through the use of histogram reweighting together with the technique of thermodynamic integration to obtain N-particle functions and, as a consequence, the virial coefficients from Eqs. (3) and (4). Evaluation of the partition functions might also provide an alternative route for the calculation of D(T) and/or the incorporation of quantum statistics (bosonic or fermionic) at low temperatures.

Acknowledgments

We thank Mike Moldover for his encouragement and helpful advice during the course of this project. Karl Johnson for facilitating the collaboration that made this work possible, Greg Schenter and Richard Wheatley for discussions in the formative stages of this work, and Krzysztof Szalewicz for discussion of potentials and for making Refs. (23) and (25) and their computer code available prior to publication. We also thank Kevin Coakley for his assistance with the uncertainty analysis, and Eric Lemmon for assistance with data fitting.

The calculations have been performed on the HPC facility Wiglaf at the Physics Department of the University of Trento.

6. References

(1) J. W. Schmidt, R. M. Gavioso, E. F. May, and M. R. Moldover. Polarizability of helium and gas metrology, Phys. Rev. Lett., 98, 254504 (2007).

(2) B. Fellmuth, Ch. Gaiser, and J. Fischer. Determination of the Boltzmann constant--status and prospects, Meas. Sci. Technol. 17, R145 (2006).

(3) L. Pitre, M. R. Moldover, and W. L. Tew, Acoustic thermometry: new results from 273 K to 77 K and progress towards 4 K, Metrologia 43, 142 (2006).

(4) R. A. Aziz, A. R. Jansen, and M. R. Moldover, Ab Initio calculations for helium: A standard for transport property measurements, Phys. Rev. Lett. 74, 1586 (1995).

(5) J. J. Hurly and M. R. Moldover, Ab Initio values of the thermo-physical properties of helium as standards, J. Res. Natl. Inst. Stand. Technol. 105, 667 (2000).

(6) J. J. Hurly and J. B. Mehl. [.sup.4]He thermophysical properties: new ab initio calculations, J. Res. Natl. Inst. Stand. Technol. 112, 75 (2007).

(7) E. Bich, R. Hellmann, and E. Vogel, Ab initio potential energy curve for the helium atom pair and thermophysical properties of the dilute helium gas. 11. Thermophysical standard values for low-density helium, Mol. Phys. 105, 3035 (2007).

(8) R. Hellmann, E. Bich, and E. Vogel, Ab initio potential energy curve for the helium atom pair and thermophysical properties of the dilute helium gas. I. Helium-helium interatomic potential, Mol. Phys. 105, 3013 (2007).

(9) J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird, Molecular Theory of Gases and Liquids, John Wiley & Sons, New York, (1954).

(10) T. Yokota, Expansion theorem of density matrix, virial expansion and new formula for multiple scattering, J. Phys. Soc. Jap. 15,779(1960).

(11) J. Ram and Y. Singh, On the quantum corrections to the virial coefficients of the equation of state of a fluid, Mol. Phys. 26, 539(1973).

(12) T. D. Lee and C. N. Yang, Many-body problem in quantum statistical mechanics. II. Virial expansion for hard-sphere gas, Phys. Rev. 116,25(1959).

(13) A. Pais and G. E. Uhlenbeck, On the quantum theory of the third virial coefficient, Phys. Rev. 116, 250 (1959).

(14) S. Y. Larsen, Quantum-mechanical calculation of the third virial coefficient of [He.sup.4], Phys. Rev. 130, 1426 (1963).

(15) S. Y. Larsen and P. L. Mascheroni, Quantum-mechanical third virial coefficient and three-body phase shifts, Phys. Rev. A 2, 1018(1970).

(16) A. S. Reiner, Application of Faddeev techiques to the quantum theory of the third virial coefficient, Phys. Rev. 151, 170 (1966).

(17) L. D. Fosdick and H. F. Jordan, Path-integral calculation of the two-particle Slater sum for [He.sup.4], Phys. Rev. 143, 58 (1966).

(18) H. F. Jordan and L. D. Fosdick, Three-particle effects in the pair distribution function for [He.sup.4] gas, Phys. Rev. 171, 128 (1968).

(19) L. D. Fosdick, The Monte Carlo method in quantum statistics, SIAM Review 10, 315 (1968).

(20) P. Diep and J. K. Johnson, An accurate [H.sub.2]-[H.sub.2] interaction potential from first principles, J. Chem. Phys. 112, 4465 (2000). Erratum, ibid., 113, 3480 (2000).

(21) G. K. Schenter, The development of effective classical potentials and the quantum statistical mechanical second virial coefficient of water, J. Chem. Phys. 117, 6573 (2002).

(22) K. Patkowski, W. Cencek, P. Jankowski, K. Szalewicz, J. B. Mehl, G. Garberoglio, and A. H. Harvey, Potential energy surface for interactions between two hydrogen molecules, J. Chem. Phys. 129, 094304 (2008).

(23) W. Cencek, M. Jeziorska, O. Akin-Ojo, and K. Szalewicz. Three-body contribution to the helium interaction potential, J. Phys. Chem. A 111, 11311 (2007).

(24) M. Jeziorska, W. Cencek, K. Patkowski, B. Jeziorski, and K. Szaliewicz, Pair potential for helium from symmetry-adapted perturbation theory calculations and from supermolccular data, J. Chem. Phys. 127, 124303 (2007).

(25) W. Cencek, K. Patkowski. and K. Szalewicz, Full-configuration-interaction calculation of three-body nonadditive contribution to helium interaction potential, J. Chem. Phys. 131, 064105 (2009).

(26) B. M. Axilrod and E. Teller, Interaction of the van der Waals type between three atoms, J. Chem. Phys. 11, 299 (1943).

(27) M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, M. Booth, and F. Rossi, GNU Scientific Library Reference Manual, Network Theory, revised second edition (2006). http://www.gnu.org/software/gsl.

(28) M. E. Boyd, S. Y. Larsen, and J. E. Kilpatrick, Quantum mechanical second virial coefficient of a Lennard-Jones gas. Helium, J. Chem. Phys. 50, 4034 (1969).

(29) D. M. Ceperley, Path integrals in the theory of condensed helium, Rev. Mod. Phys. 67, 279 (1995).

(30) G. Garberoglio, Boltzmann bias grand canonical Monte Carlo, J. Chem. Phys. 128, 134109 (2008).

(31) G. P. Lepage, A new algorithm for adaptive multidimensional integration, J. Comp. Phys. 27, 192 (1978).

(32) G. P. Lepage, VEGAS: An adaptive multi-dimensional integration program, Technical report, Cornell preprint CLNS 80-447, 1980.

(33) P. Levy, Memorial des Sciences Mathematiques, Gauthier Villars, Paris, (1954). Fascicule 126.

(34) H. Preston-Thomas, The international temperature scale of 1990 (ITS-90), Metrologia 27, 3 (1990).

(35) P. P. M. Steur, M. Durieux, and G. T. McConville, Analytic expressions for the virial coefficients B(T) and C(T) of [.sup.4]He between 2.6 K and 300 K, Metrologia 24, 69 (1987).

(36) A. L. Blancett, K. R. Hall, and F. B. Canfield, Isotherms for the He-Ar systems at 50[degrees]C, 0[degrees]C and -50[degrees]C up to 700 atm, Physica 47, 75 (1970).

(37) D. D. Dillard, M. Waxman, and R. L. Robinson, Jr., Volumetric data and virial coefficients for helium, krypton, and helium-krypton mixtures, J. Chem. Eng. Data 23, 269 (1978).

(38) H. Mansoorian, K. R. Hall, J. C. Holste, and P. T. Eubank, The density of gaseous ethane and of fluid methyl chloride, and the vapor pressure of methyl chloride, J. Chem. Thermodyn. 13, 1001 (1981).

(39) M. O. McLinden and C. Losch-Will, Apparatus for wide-ranging, high-accuracy fluid (p, [rho], T) measurements based on a compact two-sinker densimeter, J. Chem. Thermodyn. 39, 507 (2007).

(40) D. White, T. Rubin, P. Camky, and H. L. Johnston, The virial coefficients of helium from 20 to 300[degrees]K, J. Phys. Chem. 64, 1607(1960).

(41) J. A. Sullivan and R. E. Sonntag, Compressibilities for helium at temperatures from 70 to 120[degrees]K and pressures to 690 atm, Cryogenics 7, 13 (1967).

(42) K. R. Hall and F. B. Canfield, Isotherms for the He-[N.sub.2] system at-190[degrees]C, -170[degrees]C, and-160[degrees]C up to 700 atm, Physica 47, 219 (1970).

(43) D. Gugan and G W. Michel, Measurements of the polarizabili-ty and of the second and third virial coefficients of He in the range 4.2 - 27.1 K, Mol. Phys. 39, 783 (1980).

(44) J. A. Provine and F. B. Canfield, Isotherms for the He-Ar system at -130, - 115, and -90[degrees]C up to 700 atm, Physica 52, 79 (1971).

(45) A. Michels and H. Wouters, Isotherms of helium between 0[degrees] and 150[degrees]C up to 200 amagat, Physica 8, 923 (1941).

(46) G S. Kell, G. E. McLaurin, and E. Whalley, Second virial coefficient of helium from 0 to 500[degrees]C by the two-temperature gas-expansion method, J. Chem. Phys. 68, 2199 (1978).

(47) M. R. Patel, L. L. Joffrion, and P. T. Eubank, A simple procedure for estimating virial coefficients from Burnett PVT data, AIChEJ., 34, 1229(1988).

(48) C. Gaiser and B. Fellmuth, Helium virial coefficients--a comparison between new highly accurate theoretical and experimental data, Metrologia 46, 525 (2009).

(49) K. Szalewicz, University of Delaware, 2009. Private communication.

About the authors: Giovanni Garberoglio is a physicist affiliated with CNISM (Consorzio Nazionale Interuniversitario per le Scienze fisiche della Materia) and the Physics Department at the University of Trento, Italy. Allan Harvey is a chemical engineer in the Theory and Modeling of Fluids Group of the Thermophysical Properties Division of the NIST Chemical Science and Technology Laboratory. The National Institute of Standards and Technology is an agency of the U.S. Department of Commerce.

Giovanni Garberoglio

CNISM and Dipartimento di Fisica,

Universita di Trento,

via Sommarive 14,

38100 Povo (TN), Italy

and

Allan H. Harvey

Thermophysical Properties Division,

National Institute of Standards and Technology,

325 Broadway,

Boulder CO 80305, USA

garbcrog@science.unitn.it

allan.harvey@nist.gov

Accepted: June 17, 2009

Available online: http://www.nist.gov/jres

Printer friendly Cite/link Email Feedback | |

Author: | Garberoglio, Giovanni; Harvey, Allan H. |
---|---|

Publication: | Journal of Research of the National Institute of Standards and Technology |

Article Type: | Report |

Geographic Code: | 1USA |

Date: | Sep 1, 2009 |

Words: | 8919 |

Previous Article: | A method for assigning priorities to United States Measurement System (USMS) needs: nano-electrotechnologies. |

Next Article: | Perturbations from ducts on the modes of acoustic thermometers. |

Topics: |

## Reader Opinion