(4)He Thermophysical Properties: New Ab Initio Calculations.

Since 2000, atomic physicists have reduced the uncertainty of the helium-helium "ab initio" potential; for example, from approximately 0.6 % to 0.1 % at 4 bohr, and from 0.8 % to 0.1 % at 5.6 bohr. These results led us to: (1) construct a new inter-atomic potential ϕ 07, (2) recalculate values of the second virial coefficient, the viscosity, and the thermal conductivity of (4)He from 1 K to 10,000 K, and (3), analyze the uncertainties of the thermophysical properties that propagate from the uncertainty of ϕ 07 and from the Born-Oppenheimer approximation of the electron-nucleon quantum mechanical system. We correct minor errors in a previous publication [J. J. Hurly and M. R. Moldover, J. Res. Nat. Inst. Standards Technol. 105, 667 (2000)] and compare our results with selected data published after 2000. The ab initio results tabulated here can serve as standards for the measurement of thermophysical properties.

1. Introduction

In 2000, Hurly and Moldover published a comprehensive report on the application of fundamental physics to the calculation of the thermophysical properties of low-density helium [1]. The present paper is an extension to and update of parts of that paper. We developed a new model potential for the interaction of helium atoms, ϕ07(r), based on the most recent theoretical values of ϕ(r). This potential was used to calculate several important properties of 4He: The density virial coefficient B(T) and its first two temperature derivatives, the zero-density viscosity, and the zero-density thermal conductivity.

Our improved potential and calculations have significantly reduced the uncertainty of the thermophysical properties of helium. For example, at 300 K, the uncertainty of the second virial coefficient is now 1/7 of that reported in Ref. [1] and the uncertainty of the thermal conductivity is 1/3 of that reported in Ref. [1].

The new potential includes the diagonal correction to the Born-Oppenheimer model (DBOC). In addition to the use of this correction, recent discussions of the adiabatic model [2,3] recommend the use of atomic, rather than nuclear, masses in the calculations of atomic interactions. We have examined the sensitivity of the thermophysical properties to this replacement, as well as to the DBOC and to the uncertainties of the theoretical calculations of ϕ.

In the temperature range 3 K to 933 K, helium gas thermometry [4] played a leading role in the formation of the internationally accepted temperature scale, ITS-90. Subsequently, improved gas thermometry has measured T − T90, the differences between the thermodynamic temperature and ITS-90. Thus improved gas thermometry [5] may lead to a future, improved tem perature scale. Each form of gas thermometry (constant volume, dielectric, acoustic) requires the extrapolation of measured gas properties to zero pressure, where gases become “ideal.” In this work, we use fundamental principles to calculate the second density virial coefficient of helium B(T) and the second acoustic virial coefficient of helium β(T) with smaller uncertainties than can be achieved by direct measurements. Our tabulated values for B(T) and β(T) can be used to constrain the extrapolations to zero pressure; thereby leading to more accurate values of the thermodynamic temperature. Acoustic gas thermometry also requires accurate values of the thermal conductivity, which we have tabulated for helium. Recently, May et al. [6] have shown how to combine ab initio values of the viscosity of helium with comparatively simple viscosity-ratio measurements to obtain values of the thermal conductivity of argon that are more accurate than can be achieved by direct measurements. Thus, our tabulation of the viscosity of helium will also facilitate more accurate argon-based acoustic thermometry. Finally, we mention programs to redetermine the Boltzmann constant [7] and to develop an atomic standard of pressure [8] based on accurate measurements of the dielectric constant of helium at the temperature of the triple point of water (TTPW = 273.16 K). Both of these programs will benefit from the reduced uncertainties of B(T).

In the following sections, we first describe the potential and the way it was developed. We summarize the quantum-statistical formulas used for calculating the thermophysical properties of interest, then describe the numerical procedures used for the calculations. We conclude with some comparisons of our theoretical thermophysical properties with recent experimental results.

Standard notation conventions are followed in this paper. All quantum-mechanical formalism is expressed in atomic units except when noted otherwise. Interaction potentials are expressed in hartrees in the formalism, but converted to temperature units (K) for comparison with relevant literature. The CODATA-2002 values of the fundamental constants [9] were used in all calculations.

2. Model Potential ϕ07

The potential model is expressed as the sum of a repulsive term and an attractive term In these equations the cut-off radius r0 = 0.3 bohr is chosen to exclude the unphysical behavior of the potential model at small r; A = 1 hartree (E) defines the units; the a and δ are fit parameters; the C2 are fixed parameters; and the functions f2 account for relativistic retardation. The attractive part of the potential is the sum of multipole attractive terms multiplied by the universal damping functions of Tang and Toennies [10].

The dipole-dipole and higher multipole parameters C for n ≤ 10 (Table 1) are fixed at the values calculated by Zhang et al. [11]. The coefficients C for helium with the mass of 4He were used in all calculations except those which investigated corrections to the Born-Oppenheimer model. The coefficients C2 for n > 5 were estimated using the three-term recursion formula of Thakkar [12].

Table 1

Attractive interaction coefficients [11] for helium atoms with 4He and infinite mass nucleii

	4He	∞He
C6 (hartree-bohr6)	1.462122853192	1.460977837725
C8 (hartree-bohr8)	14.12578806	14.11785737
C10 (hartree-bohr10)	183.781468	183.691075
C12 (hartree-bohr12)	3267.13274	3265.256092
C14 (hartree-bohr14)	76501.2887	76571.26764
C16 (hartree-bohr16)	2277412.86	2276292.717

Zhang et al. [11] include an extensive tabulation of previous calculations for comparison. The fractional differences between the fixed-nucleon parameters of Zhang et al. and those of Bishop and Pipin [13] are 2.6 × 10−8 for C6, 1.7 × 10−7 for C8, and −9.5 × 10−7 for C10. If these fractional differences are taken as estimates of uncertainties, the total uncertainty in the potential is less than 3 × 10−8 K, and the total fractional uncertainty is less than 5 × 10−8, for r > 10 bohr.

A further, and more significant, source of uncertainty is the extrapolation formula used to estimate C2 for n > 5 from the lower-n values of C2. Thakkar [12] recommends the use of either his Eqs. (29) or (33), with the latter more appropriate for helium (based on the value of ). With the alternative formula, the estimated values of C12, C14, and C16 are 1.3 %, 5.2 %, and 13 % larger. If these differences are used as estimates of the uncertainties of the corresponding potential contributions, the total uncertainty in the potential is less than 4 × 10−5 K, and the total fractional uncertainty is less than 6 × 10−5, for r > 10 bohr.

In principle, the use of the Tang-Toennies damping terms [10] is an additional source of uncertainty. However, these functions differ from unity only for r ~ 10 bohr and below, where the quality of the fit potential can be judged directly by comparison with theoretical potentials. (See Fig. 2.)

Fig. 2

Differences between the theoretical potential values used in fitting ϕ07 and the potential, divided by the uncertainties of the potential values (See Table 3).

The retardation functions f6, f8, and f10 have been calculated by Chen and Chung [14]. Their results for f6 are in excellent agreement with the calculations of Jamieson et al. [15], whose results differ from those of Chen and Chung by a maximum fraction 1.5 × 10−5. The retardation functions satisfy f2(0) = 1; f6 decreases to ½ for r ≈ 500 and approaches 328.47/r for large r; f8 decreases to ½ for r ≈ 660 and approaches 420.62/r for large r; and f10 decreases to ½ for r ≈ 810 and approaches 508.43/r for large r. The functions f2 have the effect, for example, of converting the dipole-dipole interaction from a 1/r6 dependence to a 1/r7 dependence. Retardation has, at most, a marginal effect on all terms except the dipole-dipole term. At r = 660 bohr, the ratio C8f8/r8 to ϕ07 is less than 3 × 10−5; similarly, at r = 810 bohr, the ratio of C10 f10/r10 to ϕ07 is less than 4 × 10−10. Accordingly, the factors f12, f14, and f16, were safely approximated as unity. The code for computing the potential uses cubic spline interpolation of the results of Chen and Chung [14] for f6, f8, and f10.

The parameters δ and a, −2 ≤ j ≤ 2 were determined by fitting the potential model (1)–(3) to selected theoretical values weighted to account for their estimated uncertainties. The retardation functions f2 were set to unity in these fits. Several fits were made. The first was to the selected data set described below, and will be referred to as ϕ07. The second and third fits accounted for the uncertainties of the theoretical values, also described below. The corresponding potentials are designated ϕ07±. An additional fit was made to the potential values without applying the diagonal Born-Oppenheimer correction [16]. The corresponding potential is ϕnboc. The values of δ and the a determined by these fits are listed in Table 2.

Table 2

Variable (fit) potential coefficients

Potential	a−2 (bohr2)	a−1 (bohr)	a0 (−)	a1 (bohr−1)	a2 (bohr−2)	δ (bohr−1)
ϕ07	0.081212	−0.28755	2.14735	−1.97272	−0.051787	1.992657
ϕ07−	0.097486	−0.32441	2.17654	−1.98206	−0.050505	2.006175
ϕ07+	0.065002	−0.25089	2.11837	−1.96343	−0.053050	1.980020
ϕnboc	0.072490	−0.26814	2.13133	−1.96754	−0.052524	1.985551

2.1 Theoretical Values of ϕ

Table 3 lists values of the potential ϕ and their uncertainties based on our review of the recent literature [17-26]. The values selected for determination of ϕ07(r) represent a compromise based on availability of calculations at each r, the uncertainty claimed by the authors, and the internal agreement of various calculations for nearby r. The theoretical values were obtained within the Born-Oppenheimer model for fixed nuclear separations. Uncertainties were assigned to each of the selected values. When only a single datum was available, the authors’ uncertainty estimate was used, provided that it was consistent with neighboring values; otherwise the uncertainty was adjusted upward. When several values were available at an r-value, generally the unweighted mean and standard deviation of the more recent calculations was used. The upper-bound potentials of Komasa [19] were used only at small r, where they are in excellent agreement with the quantum-Monte-Carlo calculations of Ceperley and Partridge [17], which have much larger uncertainties.

Table 3

Theoretical potential values used to determine ϕ07. The model potential ϕ07 was fit to the sum of the theoretical potential ϕ and the diagonal Born-Oppenheimer correction ΔϕDBOC with weighting equal to the inverse square of the uncertainty U(ϕ). When a single source is listed, the uncertainty is generally that stated in the source. When multiple sources are cited, the unweighted mean and standard deviation of the set is used. In some cases, indicated by an asterisk, the uncertainty was adjusted upward to account for disagreement with neighboring values.

r (bohr)	ϕ (K)	U(ϕ) (K)	ΔϕDBOC (K)	Source(s)
1	286435	25	158	[19]
1.5	104320	20	36	[19]
2	36144.6	10	11.8	[19]
2.5	11962.0	1.0	4.1	[19,23]*
3	3786.0	20	1.37	[17,19,23,24]
3.5	1111.0	1.0	0.41	[19,20,23]*
4	292.64	0.10	0.10	[20–23,25,26]
4.5	58.400	0.10	0.009	[19,20,23,24]
5	−0.500	0.10	−0.013	[19,20,23,24]
5.1	−4.534	0.025	−0.014	[20,23]
5.6	−10.991	0.011	−0.012	[20–26]
6	−9.671	0.009	−0.011	[20,23]*
6.5	−6.887	0.005	−0.008	[23]
6.6	−6.340	0.020	−0.007	[20,24]
7	−4.619	0.007	−0.005	[26]
7.5	−3.073	0.005	−0.004	[20,23,25]*
8	−2.066	0.002	−0.002	[23]*
9	−0.989	0.001	−0.002	[25]
10	−0.5130	0.0002	−0.001	[23]*
12	−0.166	0.0010	0.000	[25]
15	−0.0423	0.0002	0.000	[25]

The diagonal Born-Oppenheimer correction calculations of Komasa, Cencek and Rychlewski [16] were interpolated using a cubic spline and added to the fixed-nucleon potentials.

Relativistic [27] (+15.4 mK) and radiative [28] (−1.3 mK) corrections to the potential have recently been evaluated only at r = 5.6 bohr. Without additional results at other r we decided, for consistency, to omit these corrections from the determination of ϕ07. The sum of these corrections is small compared with the scatter of the r = 5.6 bohr potentials in Fig. 3, but of the same order as the assigned uncertainty.

Fig. 3

Comparisons of theoretical values of ϕ with the (unretarded) model potential ϕ07. The plotted values, arranged chronologically by date of publication, are from the following sources: * [18], × [19] (upper bound), + [20], □ [21], ▲ [22], Δ [23], ▼ [25], ∇ [26]. The dotted lines represent ϕ07±. These bounds encompass the eight values of ϕ(r) published since 1999 [19–23,25,26], or overlap the authors’ k = 1 uncertainty estimates. The values of ϕ07 less the diagonal Born-Oppenheimer correction are (292.64 ± 0.13) K at 4 bohr and (−10.996 ± 0.015) K at 5.6 bohr.

The model potential defined by Eqs. (1)–(3) was fit to the sum of two quantities, the selected potentials and the corresponding DBOC. The input potentials were weighted by the inverse squares of the uncertainties U(ϕ) in the fit. The coefficients determined in the fit are listed in Table 2. The variance of the fit residuals in the determination of ϕ07 was 0.6.

The upper part of Fig. 1 shows the potential ϕ07 and the selected data used in its determination. The lower part of Fig. 1 and Fig. 3 show fractional differences between many recent theoretical potentials and ϕ07. Figure 2 shows the normalized residuals (ϕ − ϕ07)/U(ϕ).

Fig. 1

Top: The model potential ϕ07 (solid line) and theoretical values of ϕ (open circles) used in its determination. (The vertical scale is proportional to sinh−1(5ϕ/K), which is approximately logarithmic for large ϕ and linear for small |ϕ|.) Bottom: Fractional differences between theoretical values of ϕ and the model potential ϕ07, with error bars as assigned by the authors (when available). The data sources are □ [17], * [18], × [19], + [20], Δ [23], ⋄ [24], ▼ [25], ∇ [26]. The potentials ϕ07± and ϕ00 [1] are shown as solid lines.

To assess the uncertainty of ϕ07 and the propagation of this uncertainty into computed thermophysical properties, the potentials ϕ07+ and ϕ07− were developed. The potential was refitted to theoretical potentials shifted by their uncertainties, that is, to ϕ + ΔϕDBOC ± U(ϕ). Similarly, the effects of the diagonal Born-Oppenheimer correction were assessed by determining the potential ϕnboc through fits to the theoretical ϕ values without adding the correction.

The uncertainty of ϕ07 is difficult to quantify. Figure 2 shows that all but one of the theoretical potentials used in fitting ϕ07 differs from ϕ07 by less than the corresponding uncertainty, consistent with the fit variance of 0.6. Figure 3 shows that all of the theoretical values at r = 4 bohr and r = 5.6 bohr that were used in the fit either fall in the range between ϕ07− and ϕ07+ or have uncertainties overlapping this range. These observations suggest that the uncertainty in ϕ07 should be interpreted as having a large coverage factor [29] k ≈ 2. Table 4 summarizes the properties of the potentials used in this work, and the bound state energies (for angular momentum index ℓ = 0) determined from the potential.

Table 4

The potential minima ϕmin = ϕ(r) for the potentials used in this work, and the corresponding bound-state energies. (The retardation corrections f2 were included in these calculations.)

Potential	ϕmin (K)	rm (bohr)	He mass	Ebound (mK)
ϕ07	−10.999	5.608	Atomic	−1.555
ϕ07−	−10.983	5.608	Atomic	−1.667
ϕ07+	−11.014	5.607	Atomic	−1.438
ϕnboc	−10.985	5.608	Atomic	−1.550
ϕ07	−10.999	5.608	Nuclear	−1.520

3. Atomic Interactions

The thermophysical properties of helium can be evaluated using the formalism of quantum statistical mechanics. In particular, the virial coefficient of the equation of state, the viscosity, and the thermal conductivity can be expressed in terms of the phase shifts associated with the interaction of a pair of helium atoms. The theory and equations used in determining the thermophysical properties are summarized in Sec. 3.1. The following section 3.2 describes the computational techniques used to determine the thermophysical properties.

3.1 Formalism

The interaction of two atoms with a spherically symmetric potential ϕ(r) is described by a quantum mechanical wave function Ψ (r)Y/r, where r is the separation distance and Y is a spherical harmonic. The radial function satisfies where µ is the reduced mass of the He-He system, m is the electron mass, lengths are expressed in units of the Bohr radius a0, and energies are expressed in units of hartree (E).

The solutions to Eq. (4) fall into three ranges. For small r, where the potential is much larger than the angular-momentum term, the solutions must be determined numerically. In the second region of intermediate r, the potential is negligible but the angular momentum term is significant, so Eq. (4) takes the form where that is, κ is just the wave number k = κ/a0 in atomic units. The general solution of Eq. (5) is where j(ξ) and y(ξ) are spherical Bessel and Neumann functions. For large κr the asymptotic form of χ. is which can be recognized as the solution to Eq. (4) in the third region, where both the potential and the angular momentum term are negligible. The thermophysical properties of interest depend on the phase shifts δ(E). The virial coefficient of 4He depends on the sum The convergence of this sum is discussed in the next section. The viscosity and thermal conductivity depend on the quantum cross-sections [30] which are expressed in terms of much more rapidly converging sums.

3.2 Numerical Techniques

Numerical solutions of the radial equation (4) were determined with Numerov’s method using an integration step size This step size was determined empirically to insure that phase shifts obtained with step sizes h0/2 or h0/4 did not differ from those determined with step size h0 within the error criterion |Δδ| < 10−9. Calculations were made for a series of discrete energies in the range 10−11 ≤ E/E ≤ 1. The discrete energies were distributed uniformly on a logarithmic scale.

For each discrete energy, Eq. (4) was integrated upward in r, for ℓ = 0, 2, … ℓ1. A series of nodes of Ψ (r) were found at coordinates r, n = 1, 2, …. The phase shifts at node n, δ, defined by were determined successively. The asymptotic phase shift as r→∞ was obtained when the phase shifts evaluated at a series of nodes agree to within the preset convergence criterion. Convergence was accelerated by using the semi-classical (JWKB) approximation [31,32]. The convergence criterion was that the standard deviation of three successive values of δ was less than 10–9. The maximum angular momentum index ℓ1 was the minimum of either 1000 or the index when |δ| became less than 10–9.

Equation (11) only determines the phase shift within an additive multiple of π. Two conditions were used to get the total phase shifts needed in the sum (9). (1) The limiting values were lim→0δ0(E) = π and lim→0δ (E) = 0 for ℓ > 0; and (2) δ(E) is a continuous function of E [33].

Figure 4 shows the dependence of the phase shifts on ℓ and E. It is clear that for small E, the sum (9) is dominated by the ℓ = 0 term. For larger E many terms contribute to the sum. The Born approximation (see, e.g. Eq. (38.14) of Ref. [34]) for the phase shift is useful in considering the rate of convergence. For small κr the spherical Bessel function can be approximated by the leading term in the Taylor series, (κr)/(2ℓ + 1)!!. The contribution to the integral in Eq. (12) for small r thus decreases rapidly with ℓ. The spherical Bessel function has a maximum for κr near ℓ + 1. For larger ℓ the Born approximation thus becomes dependent mainly on the weaker attractive part of ϕ(r) The contributions from power-law potential terms have a simple form: which has the values and for the dipole-dipole interaction with and without retardation. The infinite sums of these terms are These can be used to get upper and lower limits for the contributions of the C6f6/r6 term to the truncation error of the sum (9). The following test was made to check the summation error. For E ≥ 0.001 hartree, Eq. (12) was used to obtain phase shifts for ℓ1 = 1000 < ℓ ≤ ℓ2 = 3000, and the corresponding contributions to the sum (9) were evaluated numerically. Equations (16) and (17) were then used with ℓ1 → ℓ2 to estimate the truncation error of these numerical sums. The results so obtained were then compared directly with upper and lower limits based on Eqs. (16) and (17). The numerical sums were found to lie very close to the product of C6 and the right-hand side of Eq. (16). The reason is that asymptotic phase shifts for ℓ = 1000 are obtained when r is some multiple of 2π/κ beyond the first zero of the spherical Bessel function j1000(κr), which occurs near . For E > 0.001 hartree this is reached before retardation is significant. For smaller E, the nodes r where Eq. (11) is evaluated occur at larger radii where retardation may be important, but the phase shifts decline sufficiently rapidly with increasing ℓ that convergence is obtained for ℓ ≪ 1000.

Fig. 4

Representative phase shifts as functions of energy. The phase shifts are all positive for small E; δ0 has a zero-E limit of π, otherwise δ(0) = 0. The lines represent, from left to right, ℓ = 0, 2, 4, 10, 20, 40, 100, 200, 400, and 1000.

4. Virial Coefficients

The second virial coefficient of 4He is [33] where and In these equations, , where is the thermal de Broglie wavelength, and −T is the bound state energy in K (Table 4). The temperature derivatives of B(T) can be evaluated directly from Eqs. (18)–(24). Numerical evaluation of the derivatives requires the integrals I2 and I4 in addition to I0.

The thermal contributions Bth(T) require numerical integration over κ. The integrals could formally be written with E as the integration variable. However, the dependence of the sum terms in the integrand for small κ was found to be approximately linear in κ ∝ , so a better spline approximation was obtained by using κ as the independent variable.

Formally, the upper limit of integration is infinite. In practice, the phase shifts become increasingly difficult to calculate at higher energies. Calculations were made only up to E = 1 hartree. The argument of the exponential factor in the integrand, −ακ2/T, has a maximum value at E = 1 hartree equal to −3.16 × 105 K/T, so the integrand is vanishingly small at κmax, even at T = 10000 K (exp(−31.6) ≈ 1.9 × 10−14). The upper limit of integration can thus be safely set at κmax.

Numerical integrations were required for the integrals I0, I2, and I4. For each case, the sum was approximated by cubic splines. The number of knots per decade of energy E was 40 for all except E > 0.1, where 80 knots were required in order to resolve the rapid dependence of the phase shifts. The integrals were calculated as the sum of a series of integrals with κ-limits 0–0.01, 0.01–0.1, 0.1–1, 1–10, and 10−κmax. This procedure insures sufficient sampling of the integrands, whose peak values depend strongly on T.

Figures 5–7 and Table 5 show the virials and the first two temperature derivatives as calculated in this work. Note that the effects of ϕ07± on the results is approximately symmetrical. Half of the difference of each calculated property, as computed with ϕ07+ and ϕ07−, was chosen as a conservative (k ≈ 2) estimate of the uncertainty U(x) of property x. These uncertainty estimates are well-represented by functions of the form with coefficients listed in Table 6. The table also includes an uncertainty calculation for the acoustic virial where γ0 = 5/3 for helium. Equation (25) represents the uncertainties of B, TB′, and T2B″ within 2 %, 3 %, and 2 % (rms), respectively, and with a maximum error less than 10 %. The uncertainty of β is represented within 2 % (rms) with a maximum error of 5 %. As noted previously, the uncertainty of ϕ07 has a large coverage factor k ≈ 2; a similar coverage factor applies to the uncertainties expressed in Eq. (25) and Table 6.

Fig. 5

Virial coefficient of 4He, calculated under various assumptions. The plots of ±B(T) = ±B07(T) were calculated with ϕ07 and atomic masses. The plotted differences are Δ = B − B07, where x designates the way the virials were calculated; x = 00, 07±, and nboc indicates the use of atomic masses and the potentials ϕ00, ϕ07±, and ϕnboc; x = nm indicates calculations with ϕ07 and nuclear, rather than atomic masses.

Fig. 6

First temperature derivative B′(T) for 4He, plotted as TB′. The plotted differences are , where x designates the type of calculation (See caption to Fig. 5.)

Fig. 7

Second temperature derivative B″(T) for 4He, plotted as T2B″. The plotted differences are , where x designates the type of calculation (See caption to Fig. 5.)

Table 5

Thermophysical properties of 4He calculated in this work. Calculated quantities are printed with at least one excess figure as an aid in smooth interpolation; for the uncertainties of B and its derivatives use Eqs. (25) and Table 6; for the uncertainties of η and λ use Eq. (37).

T (K)	B (cm3mol−1)	TB′ (cm3mol−1)	T2B″ (cm3mol−1)	η (µPa·s)	λ (mWm−1K−1)
1.0	−475.05	669.19	−1790.42	0.3292	2.632
1.2	−369.75	495.66	−1294.92	0.3405	2.720
1.4	−301.99	388.71	−986.65	0.3583	2.845
1.6	−254.99	318.11	−783.24	0.3844	3.033
1.8	−220.53	268.92	−642.77	0.4183	3.283
2.0	−194.14	233.07	−542.07	0.4586	3.588
2.25	−168.70	200.19	−451.94	0.5161	4.030
2.5	−148.92	175.86	−387.37	0.5791	4.519
2.75	−133.07	157.15	−339.36	0.6457	5.037
3.0	−120.06	142.28	−302.48	0.7141	5.569
3.5	−99.90	120.06	−249.71	0.8511	6.637
4.0	−84.96	104.14	−213.73	0.9834	7.666
4.5	−73.42	92.10	−187.46	1.1078	8.636
5.0	−64.23	82.63	−167.31	1.2239	9.542
6.0	−50.48	68.63	−138.19	1.4339	11.184
7.0	−40.683	58.72	−117.98	1.6209	12.650
8.0	−33.346	51.328	−103.05	1.7919	13.992
9.0	−27.646	45.582	−91.55	1.9514	15.244
10	−23.090	40.984	−82.39	2.1023	16.427
11	−19.366	37.216	−74.93	2.2463	17.556
12	−16.267	34.069	−68.73	2.3846	18.641
14	−11.407	29.105	−58.988	2.6472	20.699
16	−7.776	25.358	−51.679	2.8947	22.638
18	−4.965	22.423	−45.981	3.1299	24.481
20	−2.729	20.060	−41.408	3.3550	26.244
22	−0.911	18.113	−37.653	3.5715	27.939
23	−0.125	17.262	−36.014	3.6768	28.764
24	0.592	16.479	−34.509	3.7804	29.575
25	1.250	15.757	−33.121	3.8824	30.374
26	1.855	15.088	−31.837	3.9828	31.160
28	2.928	13.888	−29.537	4.1795	32.700
30	3.850	12.841	−27.534	4.3709	34.199
35	5.663	10.729	−23.496	4.8301	37.793
40	6.986	9.123	−20.432	5.2659	41.204
45	7.985	7.860	−18.024	5.6828	44.466
50	8.758	6.838	−16.077	6.0837	47.603
60	9.860	5.286	−13.117	6.8465	53.570
70	10.586	4.1591	−10.967	7.5673	59.208
80	11.0827	3.3033	−9.329	8.2549	64.585
90	11.4314	2.6306	−8.038	8.9149	69.746
100	11.6795	2.0876	−6.994	9.5519	74.726
120	11.9830	1.2650	−5.403	10.7689	84.240
140	12.1311	0.6715	−4.2464	11.9245	93.272
160	12.1903	0.2235	−3.3667	13.0310	101.919
180	12.1956	−0.1262	−2.6743	14.0968	110.248
200	12.1673	−0.4063	−2.1150	15.1284	118.308
225	12.1026	−0.6863	−1.5506	16.3769	128.062
250	12.0183	−0.9096	−1.0953	17.5862	137.510
273.16	11.9301	−1.0791	−0.7458	18.6765	146.027
275	11.9228	−1.0913	−0.7205	18.7620	146.695
298.15	11.8289	−1.2315	−0.4280	19.8245	154.994
300	11.8212	−1.2418	−0.4065	19.9084	155.649
325	11.7167	−1.3680	−0.1398	21.0288	164.400
350	11.6113	−1.4752	0.0893	22.1260	172.970
375	11.5063	−1.5671	0.2883	23.2024	181.375
400	11.4026	−1.6465	0.4626	24.2598	189.633
450	11.2008	−1.7763	0.7531	26.3241	205.753
500	11.0082	−1.8771	0.9850	28.3298	221.413
600	10.6523	−2.0209	1.3306	32.1959	251.596
700	10.3332	−2.1157	1.5740	35.9049	280.550
800	10.0462	−2.1802	1.7529	39.4880	308.517
900	9.7867	−2.2251	1.8887	42.9671	335.670
1000	9.5505	−2.25660	1.9943	46.3584	362.135
1200	9.1354	−2.29376	2.1455	52.924	413.367
1400	8.7804	−2.31003	2.2454	59.256	462.771
1600	8.4715	−2.31430	2.3139	65.402	510.71
1800	8.1990	−2.31131	2.3617	71.395	557.46
2000	7.9559	−2.30378	2.39552	77.259	603.20
2500	7.4448	−2.27457	2.44221	91.470	714.02
3000	7.03314	−2.23916	2.45852	105.185	820.95
3500	6.69073	−2.20239	2.45921	118.527	924.97
4000	6.39902	−2.16616	2.45129	131.578	1026.70
4500	6.14591	−2.13125	2.43847	144.394	1126.60
5000	5.92310	−2.09794	2.42281	157.018	1224.99
6000	5.54615	−2.03623	2.38740	181.810	1418.18
7000	5.23652	−1.98063	2.35026	206.135	1607.72
8000	4.97538	−1.93035	2.31346	230.116	1794.54
9000	4.75070	−1.88461	2.27786	253.835	1979.31
10000	4.55433	−1.84276	2.24378	277.355	2162.52

Table 6

Coefficients in Eq. (25) for estimating the uncertainty of B(T) and its temperature derivatives.

Property	c0	c1	c2	c3	c4
B	0.1341	−1.4474	0.0960	−0.00327	–
TB′	0.6612	−1.8415	0.2173	−0.02476	0.00128
T2B″	1.8238	−2.2109	0.3379	−0.04263	0.002166
βa	0.2661	−1.4560	0.1134	−0.00479	–

Figures 5–7 show that that the effects of neglecting the diagonal Born-Oppenheimer correction are no larger than the uncertainties so assigned, and that the effect of using nuclear rather than atomic masses is less than the uncertainties except at the highest temperatures. The differences between values of B(T) calculated with ϕ07 and ϕ00 [1] differ by less than the combined uncertainties (Eq. (25), Table 6 of Ref. [1]).

As noted above, the integrals required for Bth and its temperature derivatives were done numerically. The automatic integration routine was controlled by specifying an error criterion, which was set sufficiently low that errors in B and its derivatives due to the numerical integrations were negligible. This process only insures that the spline-approximated integrand is integrated accurately. It is of more interest to insure that the approximation of the sum term by a spline does not introduce a significant error into the calculation. To estimate this, the number of knots was reduced by eliminating alternate knots, and recalculating B(T) and its derivatives with the cruder spline approximation. The absolute fractional differences of the two evaluations of B(T) was less than 3 × 10−7 except near the zero of B(T). The maximum absolute fractional differences of the two evaluations of TB′(T) was 4 × 10−7, and the maximum absolute fractional differences of the two evaluations of T2B″(T) was 2 × 10−6.

We recommend the use of cubic spline interpolation for estimation of B(T) between temperatures listed in Table 5. Our tests of such interpolations indicate that the fractional interpolation error is generally much less than 10−4 except at the temperature extremes. (Interpolation near the extremes can be improved by using the tabulated higher derivatives to set the end conditions.)

5. Viscosity and Thermal Conductivity

The kinetic coefficients depend on the quantum cross-sections [30] defined by and

Equations (27)–(29) converge rapidly; numerical evaluation was straightforward. The collision integrals needed for computation of kinetic coefficients are expressed in terms of normalized cross sections, defined for even n > 0 by where r (actually an arbitrary length) is the radial position of the potential minimum. The collision integrals are defined as where β ≡ E/(k). Collision integrals with n = 2, s = 2, 4, … 10; n = 4; s = 4, 6, 8; and n = 6, s = 6 are needed for the fifth-order calculation of viscosity and thermal conductivity [35]. The collision integrals were calculated by using cubic spline representations of the collision integrals, and dividing the integrals in Eq. (31) into 11 sub-intervals, with limits 0–10−10, 10−10–10−9, … 0.1–1. This division insured adequate sampling of the integrands, whose peak locations vary rapidly with temperature. (The errors introduced by truncating the infinite integral are neglibible.)

The viscosity is [35,36] where is obtained by solving a set of linear equations for . The components of the symmetric matrix B are listed in Appendix A of Ref. [35]. In particular, since b11 = 4Ω(2,2)*, the viscosity can be expressed as Similarly, the thermal conductivity can be determined from the solution of where the components of A are defined in Appendix B of Ref. [35], and is a column vector with components ζ. The thermal conductivity depends only on ζ1:

To insure that the complicated formulas for the components of B (and the corresponding matrix for the thermal conductivity) were transcribed accurately, the following procedure was followed. The formulas were extracted from an electronic copy of Ref. [35]. These were further edited to conform with Fortran notation. Subsequently, Viehland provided Fortran codes that generated the Fortran code for calculating the matrices directly [37]. Numerical evaluations using the two implementations were identical within machine precision.

Errors in the numerical integrations required for calculating η and λ were estimated by eliminating alternate knots in the spline representations of the collision integrals and repeating the calculations. The two values of η(T), and the two values of λ(T) so determined had an absolute fractional difference of less than 3 × 10−6 at T = 1 K. This difference declined with T and remained below 1.2 × 10−7 for T > 20 K.

The viscosities and thermal conductivities determined in this work are listed in Table 5. Figure 8 and a nearly identical figure for Δλ/λ show the sensitivity of the calculations to the choice of potential and the choice of nuclear instead of atomic masses. The effects of using potential ϕ07± is nearly symmetric. Half of the differences between values of η or λ calculated with these two potentials approaches 0.35 % at low temperature. The differences reverse sign near 42 K. Above this temperature, the half-difference is bound by 0.02 %. A reasonable estimate of the the relative uncertainty U in either η or λ is the minimum of 0.35 % and the equation

Fig. 8

Sensitivity of the viscosity of 4He to various options in the calculations. The fractional difference between η and the value calculated with ϕ07 and atomic masses is plotted as the fraction Δη/η = (η − η07)/η07, where x specifies the type of calculation (See caption to Fig. 5.) A similar plot for the thermal conductivity differs from this plot only in minor details.

Values of the viscosity and thermal conductivity at temperatures between those listed in Table 5 can be obtained by interpolating with cubic splines. Our tests indicate that cubic spline interpolation introduces a fractional error of less than 10−5 except near the temperature extremes.

6. Validation of Computations

The Fortran code used for calculating the phase shifts and for subsequent calculation of the thermo-physical properties was tested by an independent development of new codes by one of us (Mehl) to test the results of Hurly and Moldover [1]. The test demonstrated excellent agreement of the sum (9) and the quantum cross-sections (27)–(29).

The test revealed two errors in the calculation of the thermophysical properties reported in Ref. [1]. The first was an incorrect sign assigned to the bound-state contribution to the published virials, which mainly affected the low temperature results for 4He and for 3He-4He mixtures. The second was due to inconsistent units conversion. The code used by Hurly to calculate the thermal conductivity was based on the equivalent of Eq. (36) in Hirschfelder et al. [36]. Their Eq. (8.2–31) uses a calorie unit in a numerical prefactor. Conversion of this to J using a current definition of the calorie introduced a factor of 1.000545 error in the thermal conductivity results published in Ref. [1]. The published values are high by this factor.

7. Comparisons With Recent Experiments

Hurly and Moldover [1] compared their results with a wide range of experimental results. Here we limit our comparisons to a few accurate experiments published since 2000. Figure 9 compares the recent second virial measurements of McLinden and Lösch-Will [38]. The agreement is excellent.

Fig. 9

Second virial coefficients of 4He measured by McLinden et al. [38], compared with values calculated with ϕ07± (dashed lines). Values calculated with ϕ07 fall between the dashed lines. The error bars indicate the experimental k = 1 uncertainties.

Figure 10 compares the recent measurements of the acoustic virial by Pitre, Moldover and Tew [5]. The measurements fall well within the combined (k = 2) uncertainty of the predicted slope β and the experimental uncertainty except at high temperatures, where the disagreement is on the order of the scatter in the measurements.

Fig. 10

Acoustic virial coefficient of Pitre et al. [5] compared with values calculated with ϕ07; Δβ = β,expt − β,calc. The dashed lines are plots of β,07± − β,07, and indicate the uncertainty of the theoretical calculation. The error bars indicate the experimental (k = 1) uncertainties. Other lines show Δβ corresponding to ϕnm, and ϕnboc. The acoustic virial is clearly sensitive to the differences between the various potentials.

Berg’s high quality measurement of the viscosity [39,40] at 298.15 K (expressed with a k = 2 uncertainty), (19.842 ± 0.014) µPa·s, and the calculated value (19.824 ± 0.004) µPa·s differ by the sum of their k = 2 uncertainties.

8. Concluding Remarks

As shown in Fig. 3, multiple research groups have provided us with very accurate ab initio “data” at 4.0 and 5.6 bohr. In order to fully exploit these data, it would be desirable to have theoretical potentials of similar accuracy at nearby r. The most demanding gas metrology is conducted near 273 K; thus, it would be very desirable to generate ab initio values of the potential at, for example r = 3.89 and 4.13 bohr (corresponding to ϕ = 200 K and 450 K) with uncertainties comparable to those already achieved at 4.0 and 5.6 bohr.