Improved First-Principles Calculation of the Third Virial Coefficient of Helium

We employ state-of-the-art pair and three-body potentials with path-integral Monte Carlo (PIMC) methods to calculate the third density virial coefficient C(T) for helium. The uncertainties are much smaller than those of the best experimental results, and approximately one-fourth the uncertainty of our previous work. We have extended our results in temperature down to 2.6 K, incorporating the effect of spin statistics that become important below approximately 7 K. Results are given for both the 3He and 4He isotopes. We have also performed PIMC calculations of the third acoustic virial coefficient γa; our calculated values compare well with the limited experimental data available. A correlating equation for C(T) of 4He is presented; differentiation of this equation provides a reliable and simpler way of calculating γa.

garberoglio@fbk.eu michael.moldover@nist.gov allan.harvey@nist.gov We employ state-of-the-art pair and three-body potentials with path-integral Monte Carlo (PIMC) methods to calculate the third density virial coefficient C(T ) for helium. The uncertainties are much smaller than those of the best experimental results, and approximately one-fourth the uncertainty of our previous work. We have extended our results in temperature down to 2.6 K, incorporating the effect of spin statistics that become important below approximately 7 K. Results are given for both the 3 He and 4 He isotopes. We have also performed PIMC calculations of the third acoustic virial coefficient γ a ; our calculated values compare well with the limited experimental data available. A correlating equation for C(T) of 4 He is presented; differentiation of this equation provides a reliable and simpler way of calculating γ a .
Accepted: June 21, 2011 Available online: http://www.nist.gov/jres representation of the pair potential [5] and a three-body potential [6] that were state-of-the-art (or nearly so) at the time the work was performed. The uncertainties obtained in Ref. [4] were smaller than those of the best experimental results by approximately an order of magnitude, with the majority of the uncertainty coming from the three-body potential. Further improvement is desirable-for example, in a proposed pressure standard based on capacitance measurement at 273.16 K, the uncertainty in C is still the largest contributor to the uncertainty budget between approximately 8 MPa and 20 MPa [7]. It is also desirable to extend the results to lower temperatures, where helium plays an important role in temperature metrology.
Recently, a state-of-the-art pair potential for helium has been developed [8,9]. The new potential incorporates not only extremely accurate results for the potential energy in the Born-Oppenheimer (BO) approximation [10], but also accurate calculations for the most important post-BO effects (adiabatic, relativistic, and quantum electrodynamics). The claimed uncertainty of the new pair potential is at least one order of magnitude smaller than that of the potential we used in Ref. [4]. This pair potential has been used to obtain highly accurate values for the second virial coefficient B(T ) and for the low-density limits of the viscosity and thermal conductivity [9]. In addition, a new three-body potential has been developed at the full-configurationinteraction (FCI) level [11], reducing the uncertainty of the three-body potential by approximately a factor of five compared to that used in Ref. [4].
In this work, we take advantage of the availability of these better potentials, and of increased computing capabilities and algorithmic improvements, to recalculate C(T ) with lower uncertainty than could be obtained in Ref. [4] and to extend our calculations to lower temperatures. We also extend our work to the third acoustic virial coefficient, and give some results for the 3 He isotope. In this paper, we will focus on aspects that differ from Ref. [4], such as the calculation of acoustic virials and the low-temperature results. The reader is referred to Ref. [4] for further background, a literature review, and details of the uncertainty analysis. Some additional details of the PIMC calculations, especially at low temperatures where spin statistics become important, are discussed in Ref. [12].

Intermolecular Potentials
For the pair potential, we use the potential first presented by Przybytek et al. [8] and described in more detail by Cencek et al. [9]. A function for the uncertainty of this potential is given in the Supplemental Material for Ref. [8], so that upper-and lower-bound potentials can be obtained by adding or subtracting the uncertainty function from the recommended pair potential. While Przybytek et al. do not attach a rigorous statistical interpretation to their "uncertainty," we believe that it is reasonable to treat it as an expanded uncertainty with coverage factor k = 2, which corresponds to a 95 % confidence limit. For calculations with 3 He, a small adjustment (negligible in the context of this work) was made to scale the adiabatic correction to the pair potential [9] to account for the different mass.
For the nonadditive three-body potential of helium, we use the FCI potential of Cencek et al. [11]. This potential is stated to have a relative uncertainty of 2 %, which again we interpret as an expanded uncertainty at the k = 2 level. For our uncertainty analysis, we construct perturbed three-body potentials FCI-(obtained by multiplying the corresponding potential by 1.02 where it is negative and by 0.98 where it is positive) and FCI+ (multiplying by 0.98 where it is negative and by 1.02 where it is positive).

Third Density Virial Coefficient
It has been shown in Refs. [4] and [12] that the second virial coefficient B(T ) and the third virial coefficient C(T ) for a quantum system obeying Boltzmann statistics can be written as (2) (3) where N A is Avogadro's constant, β = 1/k B T, and the effective two-body and three-body potentials U -2 (r) and V -3 (r 1 , r 2 ) are given by the following path-integral expressions: where U 2 (r) and V 3 (r 1 ,r 2 ,r 3 ) are the two-body and three-body potential energies, respectively, and m is the particle mass. In Eq. (5), the position of particle 3 in the τ = 0 imaginary-time slice has been fixed at the origin of the coordinate system, due to the translational invariance of the integrand. The variables r 1 and r 2 reduce to the position of particles 1 and 2, respectively, in the classical limit, that is, when T is so high that the paths x k (τ) (k = 1, 2, 3) contributing most to the path integrals shrink to a point. A similar procedure was followed for Eq. (4), where we used the translational invariance of the integrand to fix the position of particle 2 in the τ = 0 imaginary-time slice at the origin. The variable r in Eq. (4) denotes the position of particle 1 in the classical limit.
The three-body potential energy is given by (6) where U 3 (r 1 , r 2 , r 3 ) denotes the non-additive part of the three-body potential energy and U 2 (r) is the pair potential. In Eqs. (4) and (5), the path integrals are performed over all closed paths that have the origin as endpoints, that is, paths x(τ) fulfilling the conditions x(0) = x(β ) = 0. The path integrals are normalized in such a way that their discretized form reads [12] ( where x (P+1) = x (1) = 0 de Broglie wavelength of a particle with mass m. The integrand of the middle expression in Eq. (7) can be interpreted as the probability of having a ring polymer with P beads in the positions given by the coordinates x (1) ,..., x (P) [13]. In the following subsection, this quantity will be denoted as F ring . At low temperatures where Boltzmann statistics is no longer a good approximation, the above equations must be extended to incorporate Bose-Einstein or Fermi-Dirac statistics. The details of this extension are given in Ref. [12], where it was shown that the incorporation of spin statistics is necessary for both isotopes of helium at temperatures below approximately 7 K.
The number of beads P was chosen as a function of temperature according to P = 7 + (1200 K)/T for 4 He and P = 7 + (1800 K)/T for 3 He, with the resulting P rounded to the nearest integer. Preliminary calculations showed that these choices for P provide converged results (well within the statistical uncertainty of the calculation) for C(T ) throughout the range of temperatures spanned by the present work. The ring polymers were generated as described in Ref. [4]. As in Ref. [4], we used the VEGAS algorithm [14] for the numerical integration of Eq. (7). We averaged 256 independent calculations, each with 10 6 integration points, in order to obtain the final result and its statistical uncertainty. For low temperatures where spin statistics are significant, the additional non-Boltzmann terms required were calculated as described in Ref. [12] with 128 independent integrations.

Third Acoustic Virial Coefficient
The square of the speed of sound in a gas as a function of pressure on isotherms has the low-density expansion: (8) Here, u is the speed of sound; β a , γ a , and δ a are the temperature-dependent second, third, and fourth acoustic virial coefficients; M is the molar mass; and γ 0 ≡ C p / C v is the ratio of the constant-pressure heat capacity to the constant-volume heat capacity in the ideal-gas state, which is exactly 5/3 for a monatomic gas. [The ideal-gas heat-capacity ratio γ 0 should not be confused with the acoustic virial coefficient γ a .] Insofar as γ a is a second-order correction for non-ideality, it is analogous to the third density virial coefficient C. We choose to discuss RTγ a instead of γ a because RTγ a has both the units and the order of magnitude of the more familiar third density virial coefficient C. Exact thermodynamic relations connect γ a to the density virial coefficients B and C and their first two temperature derivatives [15]. These relations are: (9) where β a is related to the second density virial coefficient B(T ) and its temperature derivatives by (10) and (11) (12) For PIMC calculation of the acoustic virial coefficients, it is necessary to derive path-integral expressions for the temperature derivatives of the density virial coefficients. For this purpose, we use the derivatives of Eqs. (2) and (3) with respect to β = 1/k B T, together with the identities (13) (14) In Eqs. (2) and (3), the second and third density virial coefficients are given as a sum of terms that involve the integral of products of ring-polymer probability distributions with Boltzmann factors of the interaction potential averaged along the path. Making use of the probability distribution can be written as (15) where the last equality defines the quantity a. The Boltzmann factors of the potential along the paths quite generally have the form (16) whose β derivative is given by (17) where the last equality defines the quantity b. .
, the derivatives of the ring-polymer Path-integral expressions for the temperature derivatives of B(T ) and C(T ) can then be derived with the use of Eqs. (15) and (17). For example, the first temperature derivative of B(T ) can be written as (18) where F ring (1) is the probability distribution for the configurations of the first ring polymer, with a corresponding definition for the second, and a 1 and a 2 are the quantity a defined in Eq. (15) for the first and second ring polymer, respectively.
In the case of C(T ), Eq. (18) can be modified to calculate the temperature derivatives of the terms appearing in Eq. (3), taking into account the fact that the potential in the definition of W is actually a threebody potential and that three ring polymers must be considered. As a consequence, there is another distribution probability F ring (3) for the third particle, as well as an integration over these ring-polymer configurations measure). Moreover, a 1 + a 2 must be replaced by a 1 + a 2 + a 3 , and the integral over dr in Eq. (18) becomes an integration over dr 1 dr 2 when calculating dC/dT.
To calculate the second temperature derivatives of the virial coefficients, we use Eq. (14) together with the relations The final result for the second derivative of B(T ) is (21) The second derivative of C(T ) is given by a similar expression, after performing on Eq. (21) the same substitutions described above for the first derivative, together with replacing the term 3(P -1) by 9(P -1)/2.
The temperature derivatives of C(T ) were used to calculate L(T) according to Eq. (11). The same PIMC methodology was used as described in Sec. 3.1; the values were obtained from 256 independent calculations with 10 6 integration steps each. Table 1 shows our calculated C(T ) for 4 He. In addition to all the temperatures given in Ref. [4], Table 1 includes lower temperatures (including some corresponding to fixed points on the ITS-90 temperature scale) and a few additional intermediate temperatures.

Third Virial Coefficients
Spin statistics significantly affect C(T ) below about 7 K; this is discussed in detail in Ref. [12], where we describe the method of incorporating these effects and show the size of the various exchange contributions at low temperatures for both 3 He and 4 He.
The low-temperature values in Table 1 differ slightly from those given in Ref. [12]; we discovered a small error in our earlier implementation of the three-body potential and the values in Table 1 supersede those in Ref. [12]. These differences are smaller than the uncertainties of the calculated C(T ), and the conclusions of Ref. [12] are not affected.  in the integration

Third Acoustic Virial Coefficients
In Table 2, we show results for the third acoustic virial coefficient γ a . Because the quantity actually calculated by PIMC is L (see Sec. 3.2 and Eqs. (9) and (11)), and because of the variation of the magnitude of γ a with temperature, we tabulate the quantity RTγ a (and its expanded uncertainty as discussed below). Calculation of RTγ a via Eq. (9) requires values of B and β a , which we obtain from the work of Cencek et al. [9] and which have such small uncertainties that they can be considered exact in the context of these calculations of γ a .   The uncertainty in the path-integral calculation of γ a becomes quite large at low temperatures, which is why PIMC values at lower temperatures are not reported in Table 2. This is due to the statistical uncertainty of the Monte Carlo integration for L; the convergence behavior of this integration is much worse than that for C(T ). The reason for this difference is not completely clear, but it may be due to the fact that the quantity a that is averaged in the calculation of dC/dT [see Eq. (15)] is the so-called thermodynamic estimator of the kinetic energy. This estimator is known to be characterized by a large variance, and therefore long computations are needed to evaluate its average value with small uncertainty [16]. Moreover, the second temperature derivative needed for calculation of the acoustic virial coefficients might add further statistical noise to the calculation.
It is also possible to calculate γ a from a correlation of C(T ), as given in Sec. 4.4, by differentiating the correlation to produce dC/dT and d 2 C/dT 2 as required in Eq. (11). The quantities involving B can again be obtained from the work of Cencek et al. [9]. Values of RTγ a calculated in this manner are also shown in Table 2. They are consistent with the values calculated directly by PIMC.

Uncertainty Analysis
The analysis of uncertainty in C(T ) was similar but not identical to that described in Ref. [4]. The contributing factors are the uncertainty in the pair potential, the uncertainty in the three-body potential, and the uncertainty in the convergence of the PIMC calculation.
The standard uncertainty due to PIMC convergence was estimated as the standard deviation of the mean from the 256 independent Monte Carlo runs.
The contributions due to the uncertainties in the potentials were evaluated by calculating C(T ) with perturbed upper-and lower-bound versions of the potentials as described in Sec. 2. In order to avoid noise introduced by the PIMC convergence uncertainty, these calculations were performed with the semiclassical method described in Sec. 3.1 of Ref. [4], which was shown to be fairly accurate down to about 50 K. This approach to estimating uncertainty is adequate even at lower temperatures where the semiclassical values of C(T ) are no longer very accurate, since the needed quantity is not C(T ) itself but rather the difference between C(T ) calculated from the upper perturbed potential and C(T ) calculated from the lower perturbed potential.
Below 20 K, the semiclassical results deviate sufficiently from reality that we no longer trust them for uncertainty analysis. Instead, we observe that the uncertainty due to the potentials increases slowly and smoothly as the temperature is reduced, while the statistical uncertainty of the PIMC calculation increases more quickly. Because of these trends, the contribution of the potential uncertainty, which is our largest uncertainty component above 40 K, is roughly 60 % as large as the PIMC convergence component at 20 K. It is reasonable to assume that this trend continues, so that the uncertainty from the potentials will be less than 60 % of that from the PIMC convergence at lower temperatures. Therefore, we make the conservative estimate that the potential component of the uncertainty is 60 % of that from the PIMC convergence at all temperatures below 20 K.
The last column of Table 1 shows the resulting expanded uncertainties U(C) with coverage factor k = 2. For purposes of illustration, we summarize the uncertainty calculation for the point at 273.16 K. The standard uncertainty of the PIMC integration is 0.0069 cm 6 · mol -2 . The standard uncertainty due to the uncertainty of the two-body potential is 0.0042 cm 6 · mol -2 (1/4 of the difference between C(T )  [11] used in each case). The standard uncertainty due to the three-body potential, computed analogously with perturbed three-body potentials, is 0.0311 cm 6 · mol -2 . These are combined in quadrature to yield a standard uncertainty u(C) = 0.0321 cm 6 · mol -2 , which when multiplied by two yields an expanded uncertainty U(C) of 0.064 cm 6 · mol -2 .
The uncertainty in the acoustic third virial coefficient, shown in Table 2, is computed analogously. In this case, the analysis is less rigorous because the perturbed potentials do not necessarily define upper and lower bounds for the temperature derivatives of C that contribute to γ a . However, at all but the highest temperatures (above 2000 K), the uncertainty in γ a is dominated by the convergence uncertainty in the PIMC calculations, so the expanded uncertainties shown in Table 2 should be good estimates.

Correlation for Results
We correlated the results for C(T ) in Table 1 as a function of temperature: (22) where T* = T/(100 K) and the parameters a i and b i are given in Table 3. Equation (22) reproduces the values in Table 1 within tolerances smaller than their expanded uncertainties U(C); the fit is much closer than U(C) at high temperatures where the uncertainty is dominated by the Type B contribution from the uncertainty of the potentials. It covers the entire range from 2.6 K to 10 000 K, but should not be extrapolated outside this range. It may be differentiated to obtain dC/dT and d 2 C/dT 2 , which can be used for calculation of acoustic virial coefficients as discussed in Sec. 4

.2.
It is also convenient to have a continuous function for the expanded uncertainty U(C). This will naturally be approximate, as there is significant noise in the uncertainties given in Table 1. U(C) can be represented reasonably well as a function of temperature over the entire range of Table 1 by (23) where τ = log 10 (T/K).

Comparison With Experiment for C(T )
Extensive comparisons with experimental C(T ) data were given in Ref. [4], demonstrating that the uncertainties of calculated C(T ) were much smaller than those obtained from experiment. Since our new values are within the expanded uncertainties of those calculated previously, we do not repeat all the comparisons because the figures would look nearly identical to those in Ref. [4]. Instead, we limit our comparisons to the important range near room temperature and to the lowtemperature range that was not covered in Ref. [4].
In Fig. 1, our results are compared to those from the two most widely used experimental sources for C(T ) [17,18] at temperatures from 250 K to 325 K. The error bars on the experimental points represent expanded uncertainties with coverage factor k = 2, while the expanded uncertainties on our calculated points (see Table 1) are smaller than the size of the symbols. Our calculations are fully consistent with the experimental data, but have smaller uncertainties by factors of approximately 50.
In Fig. 2, we compare our calculated C(T ) with the available experimental data below 40 K [19][20][21][22][23][24]. Error bars drawn on the points from this work represent expanded uncertainties U(C) from Table 1; they are not drawn above 5 K because they would be smaller than the size of the symbols. For clarity, we do not draw error bars for the experimental points; in some sources [19,20] these were not reported and in the others they were usually quite large (on the order of hundreds of cm 6 · mol -2 ), often extending off the scale of Fig. 2. Gaiser et al. [24] described their data obtained by dielectric-constant gas thermometry from 3.7 K to 36 K with a smooth function for C(T ), which we show as a dashed line on Fig. 2. From a figure in Ref. [24], it appears that their expanded (k = 2) uncertainties would be on the order of 20 cm 6 · mol -2 over most of this range, becoming somewhat larger at the lowest temperatures.  Our results are generally consistent with the older experimental sources [19][20][21][22][23] within their scatter and uncertainties. We are for the most part in good agreement with the recent results of Gaiser et al. [24], with moderate disagreement at the low end of their temperature range. To examine this more closely, in Fig. 3 we plot the difference between the function of Gaiser et al. and our results as correlated by Eq. (22); Fig. 3 also shows our calculated PIMC points to demonstrate that Eq. (22) reproduces our results within their uncertainties. Our results agree with Gaiser et al. within mutual expanded uncertainties except between approximately 4 K and 8 K. It is possible that the functional form assumed for C(T ) by Gaiser et al. [24] is not the right shape to represent this system.

Comparison With Experiment for γ γ a
In this section, we compare our results for RTγ a with measurements spanning the temperature range 3 K to 423 K. In this range, RTγ a has a strong temperature dependence; it increases from approximately -13 000 cm 6 · mol -2 at 3 K to approximately 900 cm 6 · mol -2 near 14 K and then decreases to -9 cm 6 · mol -2 at 423 K. Throughout most of this range, the uncertainty of our PIMC values of RTγ a is on the order of 0.5 % to 3 %. Because of the wide range and precision required, we do not compare our results with measurements on a conventional graph. Instead, we have plotted the quantity 10 -3 (T/K) 1.5 × RTγ a , where the exponent 1.5 was chosen so that the range of the product 10 -3 (T/K) 1.5 × RTγ a in the interval from 3 K to 423 K is much smaller than the range of RTγ a (see Fig. 4).
We examined the speed-of-sound data for 4 He published in archival journals and found three publications from which we could determine accurate values of γ a [25][26][27]. Remarkably, the most recent of these publications is 35 years old. Thus, these studies did not benefit from the dramatic reduction in the uncertainty of speed-of-sound measurements that acoustic thermometry has achieved during the past 20 years [3,[28][29][30].
Gammon [25] measured the speed of sound of 4 He on 14 isotherms spanning the temperature range 98 K to 423 K at intervals of 25 K and spanning the pressure range 10 atm to 150 atm at intervals of 10 atm He calculated in this work with experimental values at near-ambient temperatures. Error bars on experimental points represent expanded uncertainties with coverage factor k = 2; uncertainties for this work (see Table 1) are not shown because the error bars would be smaller than the symbols. 4 He calculated in this work with experimental values at low temperatures. Error bars on calculated points represent expanded uncertainties with coverage factor k = 2 (see Table 1), and are not shown when they would be smaller than the symbol size. Uncertainties for experimental points are not shown for clarity (see text). (1 atm = 0.101 325 MPa). Gammon correlated his data using a classical model He-He pair potential. From his correlation, he identified the isotherm at 348 K and 10 other isolated measurements as outliers. We ignored Gammon's isolated outliers and fit the remaining data on each isotherm, including 348 K, with the function (24) In Eq. (24), we define Z a ≡ Mu 2 / (γ 0 RT lab ) as the square of the speed of sound divided by its ideal-gas value, and we calculated it from Gammon's tabulated values of u 2 . The values of Z a have the narrow range 1 to 1.38; therefore, we weighted every data point on each isotherm equally. We took β a from Cencek et al. [9], and we fitted the parameters T/T lab , γ a / (RT lab ), and, when they are statistically unequal to zero, δ a / (RT lab ) and ε a / (RT lab ). [Here, T is the thermodynamic tempera-ture of each isotherm and T lab is the temperature reported by Gammon after adjustment for the changes in the internationally accepted temperature scale. The parameter T/T lab also accounts for possible changes in M / γ 0 that would occur if there were impurities in the helium.] Seven of Gammon's 14 isotherms (98 K, 148 K, 173 K, 198 K, 248 K, 273 K, and 373 K) were very well behaved; that is, the deviations from a fit with an appropriate number of terms had no obvious pressure dependence. The standard deviation of Z a -β a p / (RT lab ) from the fits, averaged over these isotherms, was 0.000 014. For the remaining isotherms, the deviations from the fits are neither random nor have single, outlying points. Therefore, we are unable to rigorously estimate the uncertainties of γ a / RT. [We confirmed Gammon's identification of the 348 K isotherm as anomalous because a satisfactory fit required the term δ a p 3 / (RT lab ) even though the adjacent isotherms (323 K and 373 K) did not.] For the 13 isotherms (excluding 348 K), we estimated the standard uncertainty u (γ a / RT) by multiplying the result of the fitting routine by (χ 2 /N) 1/2 , where N is the number of degrees of freedom and χ 2 is the sum of the squares of the deviations of the data from the function fitted to it. All of the tabulated (Table 4) uncertainties reflect multiplication of u (RTγ a ) by an additional factor of two to approximate a 95 % confidence limit.

Fig. 2. Comparison of C(T ) for
The values of RTγ a and their expanded uncertainties U(RTγ a ) resulting from fitting Gammon's data are displayed in Fig. 4 and Table 4. Our calculations and Gammon's data agree within combined uncertainties, even though Gammon's values of RTγ a are more negative than our calculated values near the upper end of his temperature range. The uncertainties from fitting acoustic data will be underestimated if they do not account for the bias introduced by truncating the virial expansion. We crudely estimate the effect of truncation by comparing the 2nd and 4th columns in Table 4. For Gammon's isotherms at 98 K, 123 K, 148 K, and 173 K, Table 4 compares RTγ a obtained with and without the term ε a p 4 / (RT lab ) in Eq. (24). The two values of RTγ a agree within combined uncertainties. Table 4 also compares values of RTγ a obtained with and without the term δ a p 3 / (RT lab ) on the isotherms 198 K through 298 K. Except at 298 K, the two values of RTγ a on each isotherm are mutually consistent. Above 298 K, δ a was zero, within its uncertainty; however, the values of RTγ a could be influenced by contributions from δ a . We verified that the uncertainties of the values of a β a from Cencek et al. [9] did not contribute significantly to the uncertainty of γ a .  Plumb and Cataland [27] measured the speed of sound in 4 He on 21 isotherms from 2.3 K to 20 K for the purpose of determining the thermodynamic temperature T. As is often done in acoustic thermometry, Plumb and Cataland deliberately restricted their data to low densities; therefore, they were unable to determine meaningful values of γ a . Because we have β a from Cencek et al. [9], we were able to determine values of γ a on all of the isotherms except those at 5 K, 2.8 K, and 2.3 K. For each isotherm, the values of δ a and ε a in Eq. (24) were set equal to zero, and the data were weighted equally. The results are displayed on Fig. 4 and in Table 5.
Grimsrud and Werntz [26] measured the speed of sound in 4 He on eight isotherms from 2.13 K to 3.816 K. They fitted their data by, in effect, adjusting T/T lab , β a , and γ a . When analyzing their data, we weighted each measurement using the uncertainties that they tabulated. Because we fixed values of β a from Cencek et al. [9], we were able to determine values of γ a with roughly 1/5 the uncertainty achieved by Grimsrud and Werntz. These values are also shown in Table 5 and on Fig. 4. The values of RTγ a from Grimsrud and Werntz are systematically more negative than our calculated values, particularly at the lowest temperatures where, for reasons discussed in the next section, the measurements may be more accurate than our calculations.
In the narrow region of overlap, the data of Grimsrud and Werntz [26] are consistent with the data of Plumb and Cataland [27]. As the temperature decreases, both sets of data tend towards values of RTγ a smaller than those derived from our Eq. (22).
Additional values for γ a derived from acoustic experiments between 2.3 K and 34 K were reported in a conference proceeding by Plumb [31]. Unfortunately, the actual measured data were never reported, so we were not able to apply new high-accuracy values of β a [9] to obtain values of γ a consistent with the best current knowledge, as we did for Refs. [25][26][27]. We therefore do not show the data from Ref. [31] in  He and its expanded (k = 2) uncertainty from the data of Gammon [25], excluding the isotherm at 348 K. Asterisks indicate fits including additional terms, as discussed in the text

Results for 3 He
While the primary focus of this work was on the common isotope 4 He, the same methods can be used for 3 He, which is of interest for cryogenic temperature metrology. Table 6 presents values of C(T) for 3 He, along with their expanded uncertainties. More extensive discussion of the 3 He calculations at low temperatures, along with comparison with the limited experimental data, is given in Ref. [12]. We note that, in addition to the data sources below 10 K examined in Ref. [12], values of C(T ) for 3 He between 14 K and 60 K were reported by Karnus [20]; these seem to be systematically high below about 30 K, similar to the data for 4 He from the same study shown in Fig. 2.
For the same reason discussed for 4 He in Sec. 4.1, the values in Table 6 differ slightly from those reported in Ref. [12], and the new values in Table 6 should be preferred.

Accuracy of Semiclassical Calculations
In Ref. [4], we assessed the accuracy of a first-order semiclassical calculation of C(T ), concluding that the semiclassical results were adequate (in the sense of reproducing the fully quantum C(T ) from PIMC calculations within their expanded uncertainties at the k = 2 level) at temperatures above about 120 K. That conclusion can be reassessed in light of the reduced uncertainties achieved in this work. The semiclassical C(T ) deviates from our new PIMC results by more than the expanded uncertainty of our new results at temperatures below about 280 K. For example, at 273.16 K, the semiclassical calculation yields 112.939 cm 6 · mol -2 , which exceeds the PIMC value by slightly more than the expanded uncertainty given in Table 1.

Discussion
The availability of new, state-of-the-art pair and three-body potentials has allowed us to calculate C(T ) for helium with uncertainties approximately one-fourth that of our previous work [4]. In addition, we have extended the temperature range of our results, which previously had a lower bound of 24.5661 K, to 2.6 K.
We also calculated C(T ) for the 3 He isotope. The incorporation of exchange effects (non-Boltzmann statistics) was necessary to achieve accurate results for both isotopes below about 7 K.
Within the temperature range covered in our previous work [4], our new results given in Table 1 are consistent with our previous results. The present C(T ) are somewhat higher than those calculated previously, typically by an amount near one-half of the expanded (k = 2) uncertainties of the results in Ref. [4]. This change is primarily due to the more accurate three-body potential used here [11], and is consistent with a few preliminary calculations using the potential of Ref. [11] that were reported in Ref. [4].  We extended the PIMC method to calculate the third acoustic virial coefficient γ a ; our results are consistent with values of γ a obtained from experimental acoustic data, and have smaller uncertainties above 18 K. However, the statistical uncertainty in the PIMC calculation of γ a becomes quite large at lower temperatures. We believe more reliable values are obtained by differentiating Eq. (22) to obtain the temperature derivatives of C and then using those values in Eq. (11). Values derived in that way are consistent with experimental results for γ a in the entire range of our C(T ) correlation, with the possible exception of the lowest temperatures (below 4 K), where the experimental data lie slightly below our results. One might expect this approach to be less reliable at the low end of the temperature range of Eq. (22), where the uncertainty in the points to which the C(T ) function was fitted is larger and dC/dT and d 2 C/dT 2 derived from Eq. (22) would have relatively large uncertainties.
At low temperatures, the uncertainty of the present results for C(T ) is dominated by the statistical uncertainty of the PIMC integration. This could, of course, be improved somewhat simply by applying more computing resources. Above about 35 K, the uncertainty from the three-body potential becomes the largest contribution. Therefore, for metrology near room temperature, further improvement in the three-body potential would be desirable; a reduction in by a factor of two in the uncertainty of the three-body potential would produce a reduction by almost that factor in the uncertainty of C(T ) at room temperature.
The present results could also be extended to lower temperatures at the expense of more computer time. This could have some application in primary thermometry at these temperatures.
One could perform similar calculations for the "cross" third virial coefficients that characterize isotopic mixtures; these would be C 334 (T ), representing interactions among two 3 He atoms and one 4 He atom, and the similarly defined C 344 (T ). Because the natural abundance of 3 He is tiny, contributions from these coefficients would be insignificant for experiments with naturally occurring helium. We are not aware of any situation in metrology where these mixture coefficients would be useful, but, if needed, the extension of the methods used here to the mixture coefficients would be straightforward.
Another quantity of interest is the fourth virial coefficient D(T ), whose calculation would be a straightforward extension of our methods. With the reduction in the uncertainty of C(T ) achieved in the present work, D(T ) will become the largest uncertainty in some situations in metrology [7]. In principle, calculating D(T ) requires not only pair and three-body potentials but also the nonadditive four-body potential. Such a potential would be difficult to develop, but the relatively small size of the three-body effects in helium suggests that one might be able to assume the nonadditive four-body effects were negligible. It should be possible to test that assumption by performing a few high-level ab initio calculations for simple assemblies of four helium atoms (such as tetrahedrons or squares). The calculation for D(T ) would require major computing resources because of the increased dimensionality of the integral, but such a calculation may at least be feasible near room temperature where the number of beads in the ring polymers in the PIMC procedure would be relatively small.