Path-integral calculation of the fourth virial coefficient of helium isotopes
FFourth virial coefficient of helium isotopes
Path-integral calculation of the fourth virial coefficient of helium isotopes
Giovanni Garberoglio
1, 2, a) and Allan H. Harvey b) European Centre for Theoretical Studies in Nuclear Physics and Related Areas (FBK-ECT*), Trento, I-38123,Italy Trento Institute for Fundamental Physics and Applications (TIFPA-INFN), Trento, I-38123,Italy Applied Chemicals and Materials Division, National Institute of Standards and Technology, Boulder, CO 80305,USA (Dated: 8 January 2021)
We use the path-integral Monte Carlo (PIMC) method and state-of-the-art two-body and three-body potentials to cal-culate the fourth virial coefficients D ( T ) of He and He as functions of temperature from 2.6 K to 2000 K. We deriveexpressions for the contributions of exchange effects due to the bosonic or fermionic nature of the helium isotope;these effects have been omitted from previous calculations. The exchange effects are relatively insignificant for Heat the temperatures considered, but for He they are necessary for quantitative accuracy below about 4 K. Our resultsare consistent with previous theoretical work (and with some of the limited and scattered experimental data) for He;for He there are no experimental values and this work provides the first values of D ( T ) calculated at this level. Theuncertainty of the results depends on the statistical uncertainty of the PIMC calculation, the estimated effect of omittingfour-body and higher terms in the potential energy, and the uncertainty contribution propagated from the uncertainty ofthe potentials. At low temperatures, the uncertainty is dominated by the statistical uncertainty of the PIMC calculations,while at high temperatures the uncertainties related to the three-body potential and to omitted higher-order contributionsbecome dominant. I. INTRODUCTION
State-of-the-art metrology for temperature and pressure in-creasingly relies on properties of helium calculated fromfirst principles. Helium is unique among noble gases inthat its small number of electrons allows the pair interac-tion between atoms to be computed almost exactly; the non-additive three-body interaction can also be calculated accu-rately. Statistical mechanics (sometimes along with otheratomic properties such as the polarizability that can also beaccurately calculated) can then provide properties of heliumgas more accurately than they can be measured. Examplesof this approach include modern acoustic gas thermometry, dielectric-constant gas thermometry, refractive-index gasthermometry, and the recent development of a primary pres-sure standard based on dielectric measurements of helium. Helium also has the advantage of remaining in the gaseousstate to low temperatures, making it the only feasible gas formetrology below about 25 K. Most metrology uses natural he-lium which is predominantly the He isotope, but the rare iso-tope He may be used at very low temperatures due to its evenlower liquefaction temperatures.The most important quantities in these applications are thevirial coefficients that define the equation of state of a gas ofmolar density ρ at temperature T in an expansion around theideal-gas (zero-density) limit, p ρ RT = + B ( T ) ρ + C ( T ) ρ + D ( T ) ρ + · · · , (1)where p is the pressure and R is the molar gas constant. Thesecond virial coefficient B ( T ) depends on the interaction be- a) Electronic mail: [email protected] b) Electronic mail: [email protected] tween two molecules, the third virial coefficient C ( T ) dependson interactions among three molecules, and so on. While formany substances the most accurate values of these coefficientsare those obtained from careful analysis of density data, forsmall molecules first-principles calculations may be able toobtain smaller uncertainties than experiment. Calculation ofvirial coefficients from intermolecular potentials will be de-scribed in the next section; here we note that, due to the lowmass of helium, classical virial coefficient calculations are in-sufficient and quantum effects must be included (including ex-change effects at very low temperatures).For helium, in 2012 Cencek et al. reported values of B ( T ) of unprecedented accuracy calculated from a pair po-tential that incorporated higher-order effects (adiabatic cor-rection to the Born-Oppenheimer approximation, relativisticeffects, quantum electrodynamics). The potential was furtherimproved in 2017 by Przybytek et al. , and still further in 2020by Czachorowski et al. The 2020 work reports B ( T ) for both He and He with uncertainties 5-10 times smaller than theuncertainties obtained by Cencek et al. This accuracy resultsboth from the highly accurate pair potential and from the factthat an exact quantum calculation of B ( T ) is possible via aphase-shift method. For the third virial coefficient C ( T ) , no exact solution isknown, but the fully quantum result can be approached nu-merically with the path-integral Monte Carlo (PIMC) method.The most accurate first-principles values of C ( T ) come fromGarberoglio et al. , who used the same pair potentialas Cencek et al. and the three-body potential reported byCencek et al. They found that it was necessary to accountfor non-Boltzmann statistics (exchange) below approximately5 K for He and 6 K for He.
They also estimated the un-certainty of C ( T ) , which primarily arose from the uncertaintyof the three-body potential and (at low temperatures) from thestatistical uncertainty of the PIMC calculations. These results, a r X i v : . [ phy s i c s . a t m - c l u s ] J a n ourth virial coefficient of helium isotopes 2at least in the Boltzmann case, were confirmed independentlyby the group of Kofke. For many purposes, knowledge of B ( T ) and C ( T ) is suffi-cient. However, at higher pressures, terms containing D ( T ) become significant. Such terms contributed to the uncer-tainty budgets of a recent pressure standard and recent single-pressure refractive-index gas thermometry measurements atcryogenic temperatures. Efforts to calculate D ( T ) for he-lium that include quantum effects have been rare. Garberoglioperformed approximate calculations above 100 K using acentroid-based method. The group of Kofke came the clos-est to a rigorous calculation, employing PIMC (consider-ing only Boltzmann statistics) to calculate D ( T ) for He from2.6 K to 1000 K based on accurate pair and three-body po-tentials. However, no uncertainties were reported apart fromthe statistical uncertainty of the PIMC calculations, and the re-sults of Garberoglio and Harvey for C ( T ) suggest that theneglect of exchange effects will cause these Boltzmann resultsto be in error at the lowest temperatures.In this work, we use state-of-the-art pair and three-body potentials to compute D ( T ) , fully incorporating exchange ef-fects in our PIMC calculations. We also provide the first rig-orous calculations of D ( T ) for He, and provide uncertaintyestimates that include not only statistical uncertainty but alsothe contribution of uncertainties in the potentials.
II. VIRIAL COEFFICIENTS AND EXCHANGE EFFECTS
In Eq. (1), the virial coefficients are known from statisticalmechanics and depend on the N -body partition functions Q N ( V , T ) according to B ( T ) N A = − Z − Z V (2) C ( T ) N A2 = ( Z − Z ) V − Z − Z Z + Z V (3) D ( T ) N A3 = − Z − Z Z − Z + Z Z − Z V + ( Z − Z )( Z − Z Z + Z ) V − ( Z − Z ) V , (4)where N A is the Avogadro constant and the auxiliary functions Z N are defined as Z N N ! = Q N ( V , T ) V N Q ( V , T ) N . (5)The appearance of powers of the Avogadro constant in thedefinition of virial coefficients is due to the fact that Eq. (1)has been written in terms of molar quantities, as the virial co-efficients are usually reported; for the sake of conciseness wewill omit these factors in subsequent formulae for B , C , and D with the understanding that they must be applied as in Eqs.(2)-(4) to produce the virial coefficients in molar units. In general, the partition functions Q N ( V , T ) are given by Q N ( V , T ) = ∑ i (cid:48) (cid:104) i | e − β H N | i (cid:105) (6) = N ! ∑ i , σ (cid:104) i | e − β H N P σ | i (cid:105) , (7)where H N = K N + V N is the Hamiltonian of an N particlesystem, where K N the total kinetic energy and V N the to-tal potential energy of N particles, respectively. The non-additive part of the N body potential will be denoted u N in-stead, so that V = u , V = u + ∑ i < j u ( i , j ) , and V = u + ∑ i < j < k u ( i , j , k ) + ∑ i < j u ( i , j ) . The primed sum in Eq. (6) isover the many-body states | i (cid:105) with the correct symmetry uponparticle exchange in the case of bosonic or fermionic parti-cles. Equation (7) is an equivalent expression for the partitionfunction, where the sum is over a complete set of many-bodystates irrespective of the symmetry upon exchange and σ arethe permutations of N objects. P σ is an operator performingthe permutation in the Hilbert space, which is multiplied bythe sign of the permutation in the case of fermions.In general, not all the degrees of freedom of the state of thesystem that we have denoted by | i (cid:105) appear in the Hamiltonian H N . We will denote by x the set of the degrees of freedomappearing in H N (atomic coordinates, in our case) and with s the other ones (nuclear spins), so that we can write | i (cid:105) = | s (cid:105)| x (cid:105) .In the following we will use Eqs. (2–7) to derive explicitexpressions for the quantum statistical contributions to virialcoefficients that apply at low temperatures. The analogousformula for the second virial coefficient has been known for along time, but its derivation will be useful to fix the notationused in the remainder of the paper. In the case of the thirdvirial coefficient, we will provide derivation of the formulaereported in Refs. 9 and 10. The fourth virial coefficient isnewly developed in this work. A. The second virial coefficient
Let us begin with the calculation of Q . There is only one(trivial) permutation appearing in Eq. (7), and we will denoteit by P · . Performing the sum over the states we obtain Q ( V , T ) = n V Λ , (8)where n is the number of internal states of the atoms we areconsidering and Λ = h / √ π mk B T is the de Broglie thermalwavelength of a particle of mass m , which here can be eitherthe mass of the He atom or the mass of the He atom. Forhelium isotopologues, n = I + I is the nuclear spinstate; I = He (so n =
1) and I = / He ( n = Z = V .In the evaluation of Q ( V , T ) , we must consider two permu-tations. The first one is the identity, which we will denote by P : , whereas the other exchanges the labels of the two particlesand we will denote it by P | . The latter permutation is odd, andhence its contribution is weighted with a + sign in the case ofbosons and a − sign in the case of fermions. We haveourth virial coefficient of helium isotopes 3 Q ( V , T ) = (cid:34)(cid:32) ∑ s , s (cid:104) s , s | P : | s , s (cid:105) (cid:33) (cid:104) x , x | e − β H P : | x , x (cid:105) ± (cid:32) ∑ s , s (cid:104) s , s | P | | s , s (cid:105) (cid:33) (cid:104) x , x | e − β H P | | x , x (cid:105) (cid:35) = (cid:104) n (cid:104) x , x | e − β H P : | x , x (cid:105) ± n (cid:104) x , x | e − β H P | | x , x (cid:105) (cid:105) , (9) ≡ (cid:104) n Q :2 ( V , T ) ± nQ | ( V , T ) (cid:105) , (10)where the last equation defines the Boltzmann and exchangedpartition function, given by Q :2 and Q | , respectively. The fac-tor n in front of Q :2 comes from the fact that P : is the identityand the weight n in front of Q | from the fact that ∑ s , s (cid:104) s , s | P | | s , s (cid:105) = ∑ s , s (cid:104) s , s | s , s (cid:105) = ∑ s , s ( δ s , s ) = n . As detailed in our previous work, the Boltzmann parti-tion function can be evaluated in the path-integral frameworkand its final expression is equivalent to the partition func-tion of a system where each of the two quantum particles isreplaced by a classical ring polymer with P monomers (theequivalence being exact in the P → ∞ limit). The theory de-scribes the form of the probability for a ring-polymer config-uration F ( ∆ r , . . . , ∆ r P − ; m , P , T ) , where ∆ r i = r i + − r i isthe separation between the position of a bead and the subse-quent one. Since the distribution F is translationally invari-ant, it is convenient to assume that r is the null vector. No-tice that the distance between the last bead and the first oneis given by the condition of having a closed polymer, hence r P − r = − ∑ P − i = ∆ r i .Analogously, the exchanged partition function is equiva-lent to the partition function of a single ring polymer of 2 P monomers. Finally, from Eqs. (10), (5), and (2), one obtains,for any value of the nuclear spin I , B ( T ) = B Boltz ( T ) + ( − ) I I + B xc ( T ) (11) B Boltz ( T ) = − V (cid:0) Z :2 − V (cid:1) (12) B xc ( T ) = − Z | V (13)where, denoting by V :2 and V | the potential energies of theequivalent classical systems for the identity and the swap per-mutations, one has Z :2 V = (cid:90) d r (cid:68) e − β V :2 (cid:69) (14) Z | V = Λ / (cid:28) e − β V | (cid:29) . (15)The average in Eq. (14) is on the configurations of two ringpolymers of P beads, whereas the average in Eq. (15) is on justone ring polymer of 2 P beads (and mass m / F described above. The average potentials are defined as V :2 ( r ) = P P ∑ i = u (cid:16)(cid:12)(cid:12)(cid:12) r + r ( ) i − r ( ) i (cid:12)(cid:12)(cid:12)(cid:17) (16) V | = P P ∑ i = u (cid:16)(cid:12)(cid:12)(cid:12) r ( ) i − r ( ) i + P (cid:12)(cid:12)(cid:12)(cid:17) , (17)where r ( ) i and r ( ) i in Eq. (16) denote the coordinates oftwo independent ring polymers of P beads each, and r ( ) i inEq. (17) denotes the coordinates of a ring polymer of 2 P beads. The factor Λ / / in Eq. (15) comes from how theaction of the permutation operator P | on the product of theprobability distributions F of two P -monomer polymers pro-duces the probability distribution of a 2 P -monomer coalescedpolymer. B. The third virial coefficient
The same considerations leading to Eq. (10) can be appliedto the three-particle partition function appearing in the expres-sion of C ( T ) . In this case, the six permutations of three parti-cles can be conveniently divided into three subsets. The firstset includes only the identity permutation, whose representa-tion we will denote as P ∴ . The second set includes the threepermutations that swap two particles (which are odd in char-acter), whose representation will be denoted by P ·| , whereasthe third set includes the two remaining cyclic permutations,which are even, and will be denoted by P (cid:52) . The analogousexpression to Eq. (10) is then Q ( V , T ) = (cid:16) n Q ∴ ± n Q ·| + nQ (cid:52) (cid:17) . (18)Using Eqs. (18) and (10) together with (5), one can writeEq. (3) as C ( T ) = C Boltz + ( − ) I I + C odd ( T ) + C even ( I + ) (19) C Boltz = ( Z :2 − V ) V − Z ∴ − V Z :2 + V V (20) C odd ( T ) = − V (cid:16) Z ·| − V Z | (cid:17) + V (cid:0) Z :2 − V (cid:1) Z | (21) C even ( T ) = V (cid:18) Z | − V Z (cid:52) (cid:19) , (22)where, denoting as V ·| the total three-body energy when twoparticles are coalesced into a single ring polymer and by V (cid:52) ourth virial coefficient of helium isotopes 4the total three-body energy when three particles are coalescedin a single ring polymer, we have Z ∴ V = (cid:90) d r d r (cid:68) exp (cid:16) − β V ∴ ( r , r ) (cid:17)(cid:69) (23) Z ·| V = Λ / (cid:90) d r (cid:28) exp (cid:18) − β V ·| ( r ) (cid:19)(cid:29) (24) Z (cid:52) V = Λ / (cid:68) exp (cid:16) − β V (cid:52) (cid:17)(cid:69) , (25)where, analogously to Eqs. (16) and (17), we have defined V ∴ ( r , r ) = P P ∑ i = V ( r + r ( ) i , r + r ( ) i , r ( ) i + P ) (26) V ·| ( r ) = P P ∑ i = V ( r ( ) i , r + r ( ) i , r + r ( ) i + P ) (27) V (cid:52) = P P ∑ i = V ( r ( ) i , r ( ) i + P , r ( ) i + P ) . (28)In Eqs. (26)–(28), r ( k ) i for k = , , P -bead ring polymers, r ( ) i are the coordinates ofa ring polymer with 2 P beads and mass m /
2, and r ( ) i denotethe coordinates of a ring polymer with 3 P beads and mass m /
3. Note that we have slightly changed the notation fromRefs. 9 and 10, by collecting together all the terms with anodd or even character upon particle exchange.
C. The fourth virial coefficient
The permutation group of 4 particles has an even richerstructure. For our purposes, it is sufficient to recall the pres- ence of the following subsets:• The identity element, whose representation we will de-note as P :: . In this case the sum over the internal statesgives a factor n .• The swapping of a single pair. This subset has odd par-ity, and includes 6 elements. Its representation will bedenoted by P : | . The sum over the internal states resultsin a factor n .• The cyclic permutation on subsets of 3 particles. Thissubset has an even parity and includes 8 elements. Itsrepresentation will be denoted by P ·(cid:52) . The sum on theinternal states produces a factor of n .• The swapping of two distinct pairs. This subset haseven parity and includes 3 elements. Its representationwill be denoted by P || . Also in this case the sum overthe internal states produces a factor n .• The 6 remaining permutations all produce a single ringpolymer. This subset has odd parity and includes, forexample, the cyclic permutations. Its representationwill be denoted by P (cid:3) , and the sum over the internalstates produces a factor of n .We can then write the 4-particle partition function as Q ( V , T ) = (cid:16) n Q ::4 ± n Q : | + n Q ·(cid:52) + n Q || ± nQ (cid:3) (cid:17) , (29)and the expression for D ( T ) turns out to be, after lengthy butstraightforward calculations, D ( T ) = D Boltz ( T ) + D xc ( T ) (30) D Boltz ( T ) = − Z ::4 − V Z ∴ − ( Z :2 ) + V Z :2 − V V + ( Z :2 − V )( Z ∴ − V Z :2 + V ) V − ( Z :2 − V ) V (31) D xc ( T ) = ( − ) I I + D o1 ( T ) + ( I + ) D e1 ( T ) + ( − ) I ( I + ) D o2 ( T ) (32) = D o1 ( T ) + D e1 ( T ) + D o2 ( T ) [ for He ] (33) = − D o1 ( T ) + D e1 ( T ) + − D o2 ( T ) [ for He ] (34) D o1 ( T ) = − V (cid:16) Z : | − V Z ·| − Z :2 Z | + V Z | (cid:17) + V ( Z :2 − V )( Z ·| − V Z | ) + V (cid:16) Z ∴ − V Z :2 + V (cid:17) Z | − V (cid:0) Z :2 − V (cid:1) Z | (35) D e1 ( T ) = − ( Z ·(cid:52) − V Z (cid:52) ) V − ( Z || − Z | ) V + ( Z :2 − V ) Z (cid:52) V + Z | ( Z ·| − V Z | ) V − (cid:0) Z :2 − V (cid:1) Z | V (36) D o2 ( T ) = − V Z (cid:3) + V Z | Z (cid:52) − V Z | , (37)where, denoting again with V σ the total four-body potential energy for the equivalent classical system obtained by apply-ourth virial coefficient of helium isotopes 5ing the permutation σ , we have defined Z ::4 V = (cid:90) d r d r d r (cid:10) exp (cid:0) − β V ::4 (cid:1)(cid:11) (38) Z : | V = Λ / (cid:90) d r d r (cid:28) exp (cid:18) − β V : | (cid:19)(cid:29) (39) Z || V = Λ (cid:90) d r (cid:28) exp (cid:18) − β V || (cid:19)(cid:29) (40) Z ·(cid:52) V = Λ / (cid:90) d r (cid:68) exp (cid:16) − β V ·(cid:52) (cid:17)(cid:69) (41) Z (cid:3) V = Λ (cid:68) exp (cid:16) − β V (cid:3) (cid:17)(cid:69) , (42)with V ::4 ( r , r , r ) = P P ∑ i = V (cid:16) r + r ( ) i , r + r ( ) i , r + r ( ) i , r ( ) i (cid:17) (43) V : | ( r , r ) = P P ∑ i = V (cid:16) r + r ( ) i , r + r ( ) i , r ( ) i , r ( ) i + P (cid:17) (44) V || ( r ) = P P ∑ i = V (cid:16) r ( ) i , r ( ) i + P , r + r ( ) i , r + r ( ) i + P (cid:17) (45) V ·(cid:52) ( r ) = P P ∑ i = V (cid:16) r + r ( ) i , r ( ) i , r ( ) i + P , r ( ) i + P (cid:17) (46) V (cid:3) = P P ∑ i = V (cid:16) r ( ) i , r ( ) i + P , r ( ) i + P , r ( ) i + P (cid:17) . (47)In Eqs. (43)–(47), r ( k ) i with k = , , , P -bead ring polymers, r ( k ) i with k = , P -bead ring polymers and mass m / r ( ) i are the coordinates of a 3 P -bead ring polymer andmass m /
3, and r ( ) i are the coordinates of a 4 P -bead ring poly-mer and mass m / III. NUMERICAL CALCULATIONS
In order to find the optimal number of monomers P as afunction of temperature in the ring-polymer representation ofthe quantum problem, we reanalyzed highly accurate calcula-tions of B ( T ) to find the values P ( T ) above which the calcu-lated B did not change within about one part in 10 . We foundthat convergence in the investigated range of temperature wasassured by choosing P equal to nint ( + / ( T / ) . ) inthe case of He and nint ( + / ( T / ) . ) in the case of He, where nint ( x ) denotes the nearest integer to x . With thischoice, the values of P are close to what we have used in ourprevious work for T >
10 K, but optimize the utilization ofnumerical resources in the low-temperature regime where de-generacy is important.In our calculations of D ( T ) , we performed the multidimen-sional integrations in Eqs. (30)–(42) using the parallel imple-mentation of the VEGAS algorithm. The six-dimensionalintegrations leading to D Boltz were performed using 4 × Monte Carlo calls, the three-dimensional integrations leadingto D o1 using 5 × calls, and the one-dimensional integra-tions leading to D e1 using 5000 calls. The averages appearingin Eqs. (38)–(42) were evaluated using independent ring poly-mers generated by the prescription of Levy. The numberof ring polymers depends on the specific contribution to D ( T ) :we used 16 in the case of D Boltzmann , 128 for D o1 and D e1 ,and 3 × for the evaluation of D o2 . Following Kofke andcoworkers, we found it convenient to evaluate separatelythe two-body contribution to the various components of D ( T ) and the difference leading to the full calculation involving thethree-body potential.We used the very recent two-body potential by Cza-chorowski et al. and the three-body potential by Cencek etal. To the best of our knowledge, no four-body potential isavailable for helium, so in this work we set it to zero. Theadditional uncertainty due to this assumption is discussed inSec. IV C.The evaluation of D ( T ) using the procedure outlined aboveis quite CPU-intensive. At the lowest temperature investigatedhere, where the calculations are more demanding, we neededroughly 20 000 CPU-hours on a modern 2.5 GHz processor,half of which are dedicated to the evaluation of the exchangecontributions. The requirements are roughly inversely propor-tional to the temperature, hence we needed 4000 CPU-hoursat T = 10 K, 650 CPU-hours at 120 K, and so on. IV. ESTIMATION OF UNCERTAINTIESA. PIMC statistical uncertainty
The statistical uncertainty of the PIMC calculations wasevaluated as the standard error of the mean of a set of inde-pendent calculations. The number of independent calculationsthat we employed varied according to the component of D ( T ) that was calculated. We used 8 for the Boltzmann component,24 for D o1 , 64 for D e1 and D o2 . The VEGAS algorithm pro-duces its own estimation of standard uncertainty, which weverified was in very good agreement with the estimate basedon independent calculations. B. Uncertainty propagated from potentials
The usual way to evaluate the contribution to the overalluncertainty of virial coefficients due to the potentials is to per-form the calculation with perturbed two- and three-body po-tentials (where the size of the perturbation is the estimateduncertainty of the potential) and take the difference. Giventhe considerable computing requirements for the evaluationof D ( T ) , we examined a more efficient way to propagate theuncertainty from the potentials to the fourth virial coefficient,starting from the functional derivative of Eq. (4) with respectto the k -body irreducible potentials u k ( r , . . . , r k ) that D ( T ) depends on. We begin by noticing that in the classical limit,ourth virial coefficient of helium isotopes 6the functions Z N become Z N = (cid:90) e − β V N ( r ,..., r N ) N ∏ i = d r i , (48)and that the variation of the n -th virial coefficient, B n , can bewritten as δ B ( k ) n = (cid:90) δ B n δ u k ( r , . . . , r k ) δ u k ( r , . . . , r k ) k ∏ i = d r i , (49)where δ B n / δ u k is the functional derivative of the virial coef-ficient with respect to the k -body irreducible potential. FromEq. (48) one has δ Z δ u ( r , r ) = − β exp [ − β u ( r , r )] (50) δ Z δ u ( r , r ) = − β (cid:90) exp [ − β V ( r , r , r )] d r (51) δ Z δ u ( r , r , r ) = − β exp [ − β V ( r , r , r )] (52) δ Z δ u ( r , r ) = − β (cid:90) exp [ − β V ( r , . . . , r )] d r d r (53) δ Z δ u ( r , r , r ) = − β (cid:90) exp [ − β V ( r , . . . , r )] d r (54) δ Z δ u ( r , r , r , r ) = − β exp [ − β V ( r , . . . , r )] , (55) enabling evaluation of δ D by functional differentiation ofEq. (4). Actually, one can evaluate the variation of D ( T ) withrespect to the two- and three-body potentials separately; usingas an example the variation due to a perturbation of the pairpotential, an expression is obtained of the form δ D ( ) ( T ) = (cid:90) δ u ( r , r ) δ D δ u ( r , r ) d r d r . (56)A naive evaluation of the uncertainty using (the absolutevalue of) Eq. (56) when only δ u is assumed to be at anypoint the absolute value of the estimated uncertainty due tothe pair potential is at best a lower bound on the actual uncer-tainty on D ( T ) . Actually, the other factor in the integrand,that is the functional derivative δ D / δ u , is a function thathas negative and positive regions whose relative importancechanges with temperature. Our calculations show that the es-timated uncertainty using Eq. (56) would cross zero some-where around T =
30 K. Conversely, a more conservative es-timate (possibily resulting in an upper bound) to the propa-gated uncertainty from the potential to the fourth virial co-efficient can be obtained by using the absolute value of theintegrand in Eq. (56); we also decided to consider the high-est absolute value of δ D / δ u when calculated with the mostpositive ( u + = u + δ u ) or the most negative ( u − = u − δ u )pair potential, so that our final expression of the propagateduncertainty due to the k -body potential is δ D ( k ) ( T ) = (cid:90) (cid:12)(cid:12)(cid:12)(cid:12) δ u k δ D δ u k (cid:12)(cid:12)(cid:12)(cid:12) ∏ i = d r i . (57)We considered the uncertainties associated to the potentials asexpanded uncertainties with coverage factor k =
2, so we ob-tained the standard uncertainties by dividing by 2 the valuesobtained using Eq. (57). The standard uncertainties from thepair and three-body potentials were then added in quadratureto obtain the standard uncertainty due to the imperfect knowl-edge of the potential-energy surfaces.The final formulae that we obtain are: δ D ( ) ( T ) = β V (cid:90) δ u ( , ) (cid:12)(cid:12)(cid:12) − β V ( , , , ) − − β V ( , , ) − − β ( V ( , )+ V ( , )) + − β V ( , ) − − β V ( , ) (cid:16) e − β V ( , , ) − − β V ( , ) + (cid:17) − (cid:16) e − β V ( , , ) − e − β V ( , ) (cid:17) (cid:104) e − β V ( , ) − (cid:105) + − β V ( , ) (cid:16) e − β V ( , ) − (cid:17) (cid:12)(cid:12)(cid:12)(cid:12) ∏ i = d r i , (58) δ D ( ) ( T ) = β V (cid:90) δ u ( , , ) (cid:12)(cid:12)(cid:12) − β V ( , , , ) − − β V ( , , ) + − β V ( , , ) (cid:16) e − β V ( , ) − (cid:17)(cid:12)(cid:12)(cid:12) ∏ i = d r i , (59)where in actual numerical integration it can be useful to con-sider that the final result must be invariant under any permu-tation of the four particles’ coordinates so that, for exam- ple, one can use 4e − β V ( , , ) = e − β V ( , , ) + e − β V ( , , ) + e − β V ( , , ) + e − β V ( , , ) . Our final estimate for the propa-ourth virial coefficient of helium isotopes 7gated uncertainty due to the potential is then u V ( D ) = (cid:113)(cid:0) δ D ( ) ( T ) (cid:1) + (cid:0) δ D ( ) ( T ) (cid:1) . (60)In order to take into account the quantum nature of helium,we used the unchanged three-body potential together with thefourth-order Feynman–Hibbs semiclassical pair potential toevaluate Eq. (57). u FH4 ( r ) = u ( r ) + β ¯ h m (cid:18) u (cid:48)(cid:48) ( r ) + u (cid:48) ( r ) r (cid:19) + (cid:18) β ¯ h m (cid:19) (cid:18) u (cid:48)(cid:48)(cid:48)(cid:48) ( r ) + u (cid:48)(cid:48)(cid:48) ( r ) r (cid:19) . (61)In Eq. (61), the number of primes indicates the order of thederivative which, in our calculations, have been evaluated nu-merically using the central-difference formulae with a gridspacing δ r = − Å. We found that this approach providedvery good estimates for the uncertainty of B ( T ) and C ( T ) ,when compared with the traditional way of calculating thiscontribution, down to roughly T = u V ( D ) were still good by run-ning calculations with perturbed potentials and taking onefourth of the difference of the values of D ( T ) obtained in thisway, as in our previous work. Since the statistical uncertaintyis large in this regime, we took the conservative approach ofadding one fourth of the combined uncertainties of D ( T ) withthe perturbed potentials to generate an estimate of u V ( D ) . Wefound this estimate to be compatible with the one obtained us-ing the procedure described in this section in the case of He.In the case of He, the semiclassical estimation of u V ( D ) pro-duces a value of 2 × cm mol − at T = . D ( T ) , hence the approachdescribed in this section is definitely questionable at the low-est temperatures. Analogously to the case of He, we per-formed calculations with the perturbed potentials, finding agood agreement above T = u V ( D ) the values obtained using the difference of thePIMC simulation with perturbed potentials combined with thestatistical uncertainty. C. Uncertainty due to four-body potential
We performed these calculations with high-accuracy pair and three-body potentials. While this truncation of themany-body expansion is rigorous for the third virial coeffi-cient, the fourth virial coefficient is affected by any four-bodynonadditivity that might exist. Unfortunately, there have beenfew attempts to calculate four-body interactions for helium, and to our knowledge no four-body potential has been de-veloped. Our analysis of this uncertainty contribution musttherefore involve some guesswork and approximations.Because helium is not very polarizable, the multibodyforces are weak and the multibody expansion should convergequickly. This is reflected in the behavior of C ( T ) , where al-most all of the quantity is given by the pair potential and thecontribution of the three-body potential is on the order of 1%to 2%. We find similar behavior for D ( T ) , where the three-body effect is much smaller than the two-body contribution atmost temperatures. At high temperatures, however, the rela-tive contribution of three-body effects becomes large (on theorder of 30% at 2000 K); this is probably an artifact of theabsolute value of D becoming small. Because of the apparentrapid convergence of the multibody expansion, a rough esti-mate for the effect of four-body forces might be 10% of thesize of the three-body effect.In order to be more quantitative, we considered the four-body potential developed by Bade in the framework of theDrude model of dispersion. This model takes into accountonly long-range dispersion, so that it results in a pair potentialwith a r − dependence and in the case of three-body forcesit reproduces the Axilrod–Teller form. It has been used toestimate the contribution of four-body forces to D ( T ) in thecases of neon and argon. Bade’s potential depends on thesingle-atom polarizability α , which is known both experi-mentally and theoretically with high precision for He, andan unknown parameter – named ¯ h ω or V in the original pa-pers – that we fixed from the knowledge of the coefficient of r − term in the most recent pair potential for He. With thischoice, Bade’s model produces a value for the Axilrod–Tellerparameter for He that is within 4% of the actual value. Wethen proceeded to evaluate the contribution to D ( T ) from thisfour-body potential, using a semiclassical approach with the4th-order Feynmann–Hibbs expression for the pair potential.As might be expected based on the convergence of the many-body expansion, this contribution is on the order of 10% of thethree-body contribution, except near 5 K where the three-bodycontribution crosses zero.However, Bade’s model is purely for dispersion interac-tions; it does not take into account the repulsive part of thefour-body potential. For C ( T ) , it is known that the Axilrod-Teller three-body dispersion effect is only accurate at low tem-peratures; the repulsive nonadditivity has the opposite signand becomes a larger contribution to C ( T ) above approxi-mately 170 K. The magnitude of this repulsive contributionis similar to the size of the Axilrod-Teller contribution, so sim-ilar behavior might be expected for the contribution of repul-sive nonadditivity to D ( T ) . We therefore estimated an uncer-tainty contribution due to the missing four-body potential asthe maximum of 15% of the absolute value of the three-bodycontribution (interpolated at those temperatures where we didnot calculate it) and the absolute value of the dispersion con-tribution estimated from Bade’s potential. We consider this anexpanded uncertainty with coverage factor k =
2. The corre-sponding standard uncertainty is reported in Table I for Heand in Table III for He in the column labeled u ( D ) .ourth virial coefficient of helium isotopes 8 V. RESULTS AND DISCUSSIONA. Calculated D ( T ) We report the results of our calculations for He in Tables Iand II, where we also compare our results with the values re-ported by the group of Kofke.
In general, we obtain verygood agreement for the Boltzmann contribution to D ( T ) , al-though the statistical uncertainties of our calculations are gen-erally larger. In order to produce Boltzmann values of loweruncertainty based on all data, we combined their results withours at the temperatures where both studies reported results,with statistical weighting according to the uncertainty fromeach source: D = D TW u ( D TW ) + D KG u ( D KG ) u ( D TW ) − + u ( D KG ) − (62) u ( D ) = (cid:115) u ( D TW ) − + u ( D KG ) − , (63)where the subscript TW indicates the values calculated in thiswork, and the subscript KG indicates the value from Kofke’sgroup.At temperatures above about 8 K, the uncertainty budgetis dominated by the uncertainty due to the imperfect knowl-edge of the potentials. At the lower end of this temperaturerange, inspection of the contributions to u V ( D ) shows that thisis mainly due to the uncertainty of the three-body potential.Above about 20 K, the uncertainty is dominated by u ( D ) ,the uncertainty due to possible four-body interactions. At thelowest temperatures, the statistical uncertainty from the PIMCcalculations (primarily the calculation of D Boltz ) becomes asizable contribution to the overall uncertainty budget.The effect of the exchange terms on D ( T ) is minimal, atleast in the case of He. Inspection of Table I shows that thevalues of D xc ( T ) have an uncertainty comparable to their ab-solute value. The breakdown of all the contributions, reportedin Table II, shows that D xc ( T ) is obtained as the sum of termswith opposite signs and comparable magnitude, resulting inconsistent cancellations.The situation for He, reported in Tables III and IV, is qual-itatively similar to that of the heavier isotope, with two notabledifferences. First, the fourth virial coefficient increases withdecreasing temperature in the whole range considered, differ-ent from the case of He where a maximum is observed around T = . He is such thatthe various contributions to D xc ( T ) in Eq. (34) are all posi-tive, resulting in a net positive value of the exchange terms inthe whole range considered. Because the group of Kofke didnot report calculations for He, the values of D Boltz ( T ) inTable III are only those resulting from our calculations. B. Comparison with experimental data
While the virial coefficients in Eq. (1) can be extracted fromexperimental measurements, in practice this becomes quite
Temperature (K)100 200 300 400 500 D ( c m m o l - ) FIG. 1. Comparison of calculated values of D ( T ) for He with exper-imental results and with correlation fitted to first-principlescalculations by Schultz and Kofke. Error bars represent expandeduncertainties with coverage factor k =
2. Expanded uncertainties forthis work are smaller than the size of the symbols; see Table I. difficult for the higher coefficients because they require ex-tracting higher-order effects (a third derivative in the case of D ). Relatively high pressures are also necessary to reach den-sities where D ( T ) is significant.The Burnett method of systematic expansions between twovessels of similar size can reduce the uncertainties in mea-suring virial coefficients. This method was used by Canfieldand coworkers in a series of papers on helium and its mixtureswith other gases, several of which reported virial coefficientsup to the fourth for helium at various temperatures. Noneof these papers reported a full uncertainty analysis; when plot-ting these points we use the parameter uncertainties reportedin the papers, which appear to be only the statistical uncer-tainties from their fits to their data.Highly accurate measurements on helium were performedby McLinden and Lösch-Will with a two-sinker magnetic-suspension densimeter. The same method was used byMoldover and McLinden, who reanalyzed both their dataand the earlier measurements. This analysis was constrainedwith ab initio values of B ( T ) and C ( T ) , allowing them toreduce the uncertainty of D on each isotherm by a factor of 6compared to unconstrained analysis.In Fig. 1, our results from Table I are shown along withthe experimental data for He. The data from the Burnett ex-periments of 50 years ago are in qualitative agreementwith our results, but are too scattered for meaningful quanti-tative comparison. The results of Moldover and McLinden are in excellent agreement with our calculations above 300 K.At lower temperatures, there is a disagreement which is smallbut outside the mutual uncertainties. We note that the pointsof Moldover and McLinden that agree well are all from their“2007 isotherms,” while their points from “2005 isotherms”show a systematic offset. This may suggest an unrecognizedproblem with the 2005 experiments. The 2007 isotherms weremeasured with new sinkers whose volumes were newly cali-ourth virial coefficient of helium isotopes 9 TABLE I. The fourth virial coefficient of He and its contributions: D Boltz , TW ( T ) , Boltzmann contribution calculated in this work; D Boltz , KG ( T ) , Boltzmann contribution from Ref. 14; D Boltz ( T ) , combined Boltzmann values from this work and Ref. 14; D xc ( T ) exchangecontribution; u V ( D ) , propagated uncertainty from the pair and three-body potential; u ( D ) , estimated uncertainty due to the neglect of thefour-body potential; D ( T ) , final values. The numbers in parentheses report standard uncertainty on the last digits. In the last column, thecombined expanded uncertainty of D ( T ) is reported at k = D Boltz , TW ( T ) D Boltz , KG ( T ) D Boltz ( T ) D xc ( T ) u V ( D ) u ( D ) D ( T ) U ( D ) (K) (cm mol − ) (cm mol − ) (cm mol − ) (cm mol − ) (cm mol − ) (cm mol − ) (cm mol − ) (cm mol − )2 . . .
25 -17000(3000) -17000(3000) 1100(600) 3000 300 -16000 80003 . .
75 12400(1500) 12000(1500) 0(300) 1100 200 12000 40004 19000(2000) 19500(1100) 19400(1000) 160(80) 800 180 20000 30004 .
25 20800(900) 20800(900) 50(120) 600 160 21000 20004 . .
75 21000(500) 21000(500) -49(40) 400 120 21000 13005 21100(900) 19200(400) 19500(400) -17(20) 300 110 19500 10005 .
25 17300(400) 17300(400) -1(20) 300 100 17300 10005 . . . .
15 2025.61(12) 2025.61(12) 1.4 8 2030 1580 1781.3(7) 1781.3(7) 1.1 6 1781 1383 .
15 1742.22(7) 1742.22(7) 1.1 6 1742 12103 .
15 1526.67(4) 1526.67(4) 0.9 5 1527 10120 1381.7(5) 1381.7(5) 0.8 4 1382 9123 .
15 1358.59(4) 1358.59(4) 0.8 4 1359 9143 .
15 1224.22(3) 1224.22(3) 0.7 4 1224 8173 .
15 1066.096(16) 1066.096(16) 0.6 3 1066 6200 955.4(4) 955.4(4) 0.6 3 955 6223 .
15 876.991(16) 876.991(16) 0.6 3 877 5250 800.1(3) 800.1(3) 0.5 2 800 5273 .
15 744.051(12) 744.051(12) 0.5 2 744 5273 .
16 743.5(3) 743.5(3) 0.5 2 744 5300 687.0(3) 687.0(3) 0.5 2 687 4323 .
15 645.170(10) 645.170(10) 0.5 2 645 4400 534.0(2) 534.115(7) 534.115(7) 0.5 2 534 4500 433.9(2) 434.029(5) 434.029(5) 0.4 2 434 4600 363.438(5) 363.438(5) 0.4 2 363 4700 310.80(18) 310.920(4) 310.920(4) 0.4 2 311 4800 270.314(4) 270.314(4) 0.4 1.9 270 4900 238.021(4) 238.021(4) 0.4 1.9 238 41000 211.50(14) 211.702(3) 211.702(3) 0.4 1.9 212 41500 130.43(10) 130.43(10) 0.3 1.7 130 42000 89.14(8) 89.14(8) 0.3 1.6 89 3 brated; it is plausible that there may have been a small errorin calibration for the sinkers used in 2005. Figure 1 alsoshows the equation that Schultz and Kofke fit to their resultsat 20 K and above, which were calculated using Boltzmann statistics only. It is not surprising that this equation is in ex-cellent agreement with our results, since exchange effects arenegligible at these temperatures.ourth virial coefficient of helium isotopes 10 TABLE II. The exchange contributions to D ( T ) for He, see Eqs. (30)–(37). Uncertainties here are k = D o1 ( T ) D e1 ( T ) D o2 ( T ) D xc ( T ) (K) (cm mol − ) (cm mol − ) (cm mol − ) (cm mol − )2 . . .
25 -286(585) 1945(149) -596(5) 1062(603)3 . .
75 -491(298) 580(54) -113.9(13) -25(303)4 -198(74) 407(33) -51.9(9) 158(81)4 .
25 -77(116) 151(14) -24.6(4) 50(117)4 . .
75 -91(40) 47(6) -5.64(11) -49(40)5 -47(21) 32(3) -1.39(4) -17(21)5 .
25 -21(22) 21(2) -0.73(2) -1(22)5 . He and its contributions. D xc is defined in Eq. (32) and the other quantities are as in Table I. Theexpanded uncertainty in the last column is at k = D Boltz ( T ) D xc ( T ) u V ( D ) u ( D ) D ( T ) U ( D ) (K) (cm mol − ) (cm mol − ) (cm mol − ) (cm mol − ) (cm mol − ) (cm mol − )2 . .
16 757.6(3) 0.5 2 758 5300 699.3(3) 0.5 2 699 4400 541.1(2) 0.5 2 541 4500 438.61(17) 0.4 2 439 4700 313.22(16) 0.4 1.9 313 41000 212.84(15) 0.4 1.9 213 41500 130.88(5) 0.3 1.8 131 42000 89.41(4) 0.3 1.6 89 3
C. Low-temperature behavior
Since one of the novel features of this work is the incorpo-ration of exchange effects that become important at low tem-perature, in Figs. 2 and 3 we show low-temperature results for He and He, respectively. No low-temperature experimentaldata are available for comparison in either case.In Fig. 2, we see that D ( T ) for He goes through a maxi-mum near 4.5 K, turning sharply negative at the lowest tem- peratures (the three lowest temperatures in Table I are not plot-ted; they would be far below the bottom of the graph). We donot show the effect of exchange, since as discussed above theexchange contributions are relatively small compared to theiruncertainty due to terms of opposite sign. The correlation thatSchultz and Kofke fitted to their Boltzmann results is in ex-cellent agreement with the results obtained here down to itslower temperature limit of 20 K. It extrapolates well down to15 K, but becomes increasingly inaccurate when extrapolatedourth virial coefficient of helium isotopes 11 TABLE IV. The exchange contributions to D ( T ) for He, see Eqs. (30)–(37). Uncertainties here are k = D o1 ( T ) D e1 ( T ) D o2 ( T ) D xc ( T ) (K) (cm mol − ) (cm mol − ) (cm mol − ) (cm mol − )2 . Temperature (K)0 10 20 30 40 D ( c m m o l - ) -20000-1000001000020000 Schultz & Kofke (2019)This workExtrapolated Schultz & Kofke FIG. 2. Calculated values of D ( T ) for He at low temperatures.The correlation fitted to first-principles calculations by Schultz andKofke is extrapolated below 20 K. Error bars represent expandeduncertainties with coverage factor k =
2, and are smaller than thesymbol size at higher temperatures. further.Figure 3 shows the low-temperature results for He. Nomaximum is evident in the temperature range investigated. Wealso show the values of D ( T ) obtained when only Boltzmannstatistics are considered. Incorporation of exchange effects isnecessary for quantitative accuracy below roughly 4 K. Thesignificant exchange effects on D ( T ) for He, compared to asmaller effect for He due to competing terms, is similar tothe behavior found in our previous work for C ( T ) . VI. CONCLUSIONS
The values of D ( T ) presented in Table I (for He) andTable III (for He) represent the first such values calculatedfrom the current state-of-the-art pair and three-body poten-tials. They are also the first to include exchange effects, andthe first to include complete uncertainty estimates. For He,Table III presents the first fully quantum calculations of D ( T ) .Our He calculations at the level of Boltzmann statisticsconfirm previous calculations (with a slightly different two-
Temperature (K)2 3 4 5 6 7 8 D ( c m m o l - ) FIG. 3. Calculated values of D ( T ) for He at low temperatures, alsoshowing results from Boltzmann statistics alone. Error bars representexpanded uncertainties with coverage factor k = body potential) from the group of Kofke. Because thetwo-body potential used in that work differs negligibly fromthe current state of the art for the purpose of calculating D ( T ) ,we could consider those calculations to be additional data thatsupplement ours, making use of their somewhat smaller statis-tical uncertainties to improve our estimates of the Boltzmanncontributions.The exchange contributions to D for He are relativelysmall, due to terms of opposite sign that mostly cancel eachother. Because of this cancellation, our statistical uncertaintyfor the exchange contribution is similar in magnitude to thecontribution itself. In contrast, for the fermion He, the ex-change contributions all contribute in the positive direction,producing a significant effect on D ( T ) at temperatures ofroughly 4 K and below.Several factors contribute to the uncertainty of the results.At low temperatures, the statistical uncertainty of the PIMCcalculations is significant, not only in its own right but also inits contribution to the uncertainty u V ( D ) due to the uncertaintyin the potential, which at the lowest temperatures is estimatedby taking differences between PIMC calculations with shiftedpotentials. This could be somewhat improved with more com-puter resources, and also by employing some of the optimizedourth virial coefficient of helium isotopes 12methods for integrating virial coefficients developed by thegroup of Kofke. Uncertainty due to imperfect knowledge of the potentialcannot be reduced simply with more computer time; this as-pect, which dominates the uncertainty at higher temperatures,requires careful development of potential-energy surfaces.The pair potential is now known with extraordinary accuracy; further improvement may be desirable for other reasons but itwill not reduce the uncertainty of virial coefficients beyond B ( T ) . The three-body potential has an estimated expanded( k =
2) uncertainty of 2% at all configurations, which couldbe reduced somewhat with a concerted effort. Finally, thefour-body potential is unknown, and is the largest uncertaintyabove roughly 10 K. Even an approximate four-body nonaddi-tive potential, for example with a 20% uncertainty, would havea large impact on reducing the uncertainty of D ( T ) . The mostimportant aspect of such a surface for metrology near roomtemperature would be its behavior at short-distance configu-rations dominated by repulsion, since those are the most im-portant configurations at those temperatures. For use at cryo-genic temperatures where the four-body uncertainty is stillsignificant but the dispersion contribution is likely more im-portant, its long-range behavior should approach that derivedby Bade. ACKNOWLEDGMENTS
The authors thank David Kofke and Andrew Schultz of theUniversity at Buffalo for useful discussions on the calculationof the Boltzmann part of D ( T ) , Mark McLinden of NIST forinsight into the experiments reported in Refs. 39 and 40, andDan Friend for suggestions that improved clarity.G.G. acknowledges support from Real-K project 18SIB02,which has received funding from the EMPIR programme co-financed by the Participating States and from the EuropeanUnion’s Horizon 2020 research and innovation programme.The calculations for T ≥
20 K have been performed onthe computing cluster KORE at Fondazione Bruno Kessler.We acknowledge the CINECA award IscraC-FAVOHLA un-der the ISCRA initiative, for the availability of high perfor-mance computing resources and support.
DATA AVAILABILITY
The data that support the findings of this study are availablewithin this article. M. R. Moldover, R. M. Gavioso, J. B. Mehl, L. Pitre, M. de Podesta, andJ. T. Zhang, “Acoustic gas thermometry,” Metrologia , R1–R19 (2014). C. Gaiser, T. Zandt, and B. Fellmuth, “Dielectric-constant gas thermome-try,” Metrologia , S217–S226 (2015). P. M. C. Rourke, C. Gaiser, B. Gao, D. M. Ripa, M. R. Moldover, L. Pitre,and R. J. Underwood, “Refractive-index gas thermometry,” Metrologia ,032001 (2019). C. Gaiser, B. Fellmuth, and W. Sabuga, “Primary gas-pressure standardfrom electrical measurements and thermophysical ab initio calculations.”Nature Phys. , 177–180 (2020). W. Cencek, M. Przybytek, J. Komasa, J. B. Mehl, B. Jeziorski, and K. Sza-lewicz, “Effects of adiabatic, relativistic, and quantum electrodynamics in-teractions on the pair potential and thermophysical properties of helium,” J.Chem. Phys. , 224303 (2012). M. Przybytek, W. Cencek, B. Jeziorski, and K. Szalewicz, “Pair potentialwith submillikelvin uncertainties and nonadiabatic treatment of the halostate of the helium dimer,” Phys. Rev. Lett. , 123401 (2017). P. Czachorowski, M. Przybytek, M. Lesiuk, M. Puchalski, and B. Jeziorski,“Second virial coefficients for He and He from an accurate relativisticinteraction potential,” Phys. Rev. A , 042810 (2020). J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird,
Molecular Theory of Gasesand Liquids (John Wiley & Sons, New York, 1954). G. Garberoglio and A. H. Harvey, “Path-integral calculation of the thirdvirial coefficient of quantum gases at low temperatures,” J. Chem. Phys. , 134106 (2011). G. Garberoglio and A. H. Harvey, “Erratum: Path-integral calculation ofthe third virial coefficient of quantum gases at low temperatures,” J. Chem.Phys. , 199903 (2020). G. Garberoglio, M. R. Moldover, and A. H. Harvey, “Improved first-principles calculation of the third virial coefficient of helium,” J. Res. Natl.Inst. Stand. Technol. , 729–742 (2011). G. Garberoglio, M. R. Moldover, and A. H. Harvey, “Erratum: Improvedfirst-principles calculation of the third virial coefficient of helium,” J. Res.Natl. Inst. Stand. Technol. , 125019 (2020). W. Cencek, K. Patkowski, and K. Szalewicz, “Full-configuration-interaction calculation of three-body nonadditive contribution to helium in-teraction potential,” J. Chem. Phys. , 064105 (2009). K. R. S. Shaul, A. J. Schultz, and D. A. Kofke, “Path-integral Mayer-sampling calculations of the quantum Boltzmann contribution to virial co-efficients of helium-4,” J. Chem. Phys. , 184101 (2012). A. J. Schultz and D. A. Kofke, “Virial coefficients of helium-4 from abinitio -based molecular models,” J. Chem. Eng. Data , 3742–3754 (2019). B. Gao, H. Zhang, D. Han, C. Pan, H. Chen, Y. Song, W. Liu, J. Hu,X. Kong, F. Sparasci, M. Plimmer, E. Luo, and L. Pitre, “Measurement ofthermodynamic temperature between 5 K and 24.5 K with single-pressurerefractive-index gas thermometry,” Metrologia , 065006 (2020). G. Garberoglio, “Quantum effects on virial coefficients: a numerical ap-proach using centroids,” Chem. Phys. Lett. , 19–23 (2012). T. L. Hill,
An Introduction to Statistical Thermodynamics (Dover, NewYork, 1987). T. L. Hill,
Statistical Mechanics (Dover, New York, 1987). M. E. Boyd, S. Y. Larsen, and J. E. Kilpatrick, “Quantum mechanical sec-ond virial coefficient of a Lennard-Jones gas. Helium,” J. Chem. Phys. ,4034–4055 (1969). G. Lepage, “VEGAS: An adaptive multi-dimensional integration program,”Tech. Rep. (Cornell preprint CLNS 80-447, 1980). R. Kreckel, “Parallelization of adaptive MC integrators,” Comput. Phys.Commun. , 258–266 (1997). P. Levy,
Memorial des Sciences Mathematiques (Gauthier Villars, Paris,1954) fascicule 126. H. F. Jordan and L. D. Fosdick, “Three-particle effects in the pair distribu-tion function for He gas,” Phys. Rev. , 128–149 (1968). R. Rodríguez-Cantano, R. Pérez de Tudela, M. Bartolomei, M. I. Hernán-dez, J. Campos-Martínez, T. González-Lezana, P. Villarreal, J. Hernández-Rojas, and J. Bretón, “Examination of the Feynman–Hibbs approach inthe study of NeN-coronene clusters at low temperatures,” J. Phys. Chem. A , 5370–5379 (2016). G. Garberoglio and A. H. Harvey, “First-principles calculation of the thirdvirial coefficient of helium,” J. Res. Natl. Inst. Stand. Technol. , 249–262 (2009). W. Cencek, G. Garberoglio, A. H. Harvey, M. O. McLinden, and K. Sza-lewicz, “Three-body nonadditive potential for argon with estimated un-certainties and third virial coefficient,” J. Phys. Chem. A , 7542–7552(2013). W. L. Bade, “Drude-model calculation of dispersion forces. I. General the-ory,” J. Chem. Phys. , 1280–1284 (1957). W. L. Bade, “Drude-model calculation of dispersion forces. III. The fourth-order contribution,” J. Chem. Phys. , 282–284 (1958). J. Wiebke, E. Pahl, and P. Schwerdtfeger, “Up to fourth virial coeffi-cients from simple and efficient internal-coordinate sampling: Application ourth virial coefficient of helium isotopes 13 to neon,” J. Chem. Phys. , 014508 (2012). J. Wiebke, E. Pahl, and P. Schwerdtfeger, “Sensitivity of the thermal andacoustic virial coefficients of argon to the argon interaction potential,” J.Chem. Phys. , 064702 (2012). C. Gaiser and B. Fellmuth, “Polarizability of helium, neon, and argon: Newperspectives for gas metrology,” Phys. Rev. Lett. , 123203 (2018). M. Puchalski, K. Szalewicz, M. Lesiuk, and B. Jeziorski, “QED calculationof the dipole polarizability of helium atom,” Phys. Rev. A , 022505(2020). A. J. Thakkar, “The generator coordinate method applied to variational per-turbation theory. Multipole polarizabilities, spectral sums, and dispersioncoefficients for helium,” J. Chem. Phys. , 4496–4501 (1981). A. L. Blancett, K. R. Hall, and F. B. Canfield, “Isotherms for the He–Arsystem at 50 ◦ C, 0 ◦ C and -50 ◦ C up to 700 atm,” Physica , 75–91 (1970). K. R. Hall and F. B. Canfield, “A least-squares method for reduction ofBurnett data to compressibility factors and virial coefficients,” Physica ,99–108 (1970). K. R. Hall and F. B. Canfield, “Isotherms for the He–N system at -190 ◦ C,-170 ◦ C and -160 ◦ C up to 700 atm,” Physica , 219–226 (1970). J. A. Provine and F. B. Canfield, “Isotherms for the He–Ar system at -130,-115, and -90 ◦ C up to 700 atm,” Physica , 79–91 (1971). M. O. McLinden and C. Lösch-Will, “Apparatus for wide-ranging, high-accuracy fluid ( p , ρ , T ) measurements based on a compact two-sinker den-simeter,” J. Chem. Thermodyn. , 507–530 (2007). M. R. Moldover and M. O. McLinden, “Using ab initio “data” to accu-rately determine the fourth density virial coefficient of helium,” J. Chem.Thermodyn. , 1193–1203 (2010).41