We present the first-principles determination of the three-body polarizability and the third dielectric virial coefficient of helium. Coupled-cluster and full configuration interaction methods were used to perform electronic structure calculations. The mean absolute relative uncertainty of the trace of the polarizability tensor, resulting from the incompleteness of the orbital basis set, was found to be 4.7%. Additional uncertainty due to the approximate treatment of triple and the neglect of higher excitations was estimated at 5.7%. An analytic function was developed to describe the short-range behavior of the polarizability and its asymptotics in all fragmentation channels. We calculated the third dielectric virial coefficient and its uncertainty using the classical and semiclassical Feynman–Hibbs approaches. The results of our calculations were compared with experimental data and with recent Path-Integral Monte Carlo (PIMC) calculations [Garberoglio et al., J. Chem. Phys. 155, 234103 (2021)] employing the so-called superposition approximation of the three-body polarizability. For temperatures above 200 K, we observed a significant discrepancy between the classical results obtained using superposition approximation and the ab initio computed polarizability. For temperatures from 10 K up to 200 K, the differences between PIMC and semiclassical calculations are several times smaller than the uncertainties of our results. Except at low temperatures, our results agree very well with the available experimental data but have much smaller uncertainties. The data reported in this work eliminate the main accuracy bottleneck in the optical pressure standard [Gaiser et al., Ann. Phys. 534, 2200336 (2022)] and facilitate further progress in the field of quantum metrology.

Interactions of atoms induce small changes in the polarizability of atomic clusters. In turn, these small changes affect the measured dielectric and optical properties of gases.1–3 In particular, the isotropic interaction-induced polarizability determines the polarized Raman scattering spectrum,4 while the polarizability anisotropy determines the depolarized spectrum and contributes to the density dependence of the Kerr constant.5 The mean polarizability is also related to the dielectric virial coefficients of gases through the Clausius–Mossotti equation.6 The dielectric virial expansion has been extensively used in metrology, formerly in the determination of the Boltzmann constant7 and recently in the development of new temperature8–11 and pressure standards.12,13 At present, the dielectric-constant gas thermometry (DCGT),8–10 the refractive-index gas thermometry (RIGT),14–16 and capacitance measurements enable the determination of thermodynamic temperature with an accuracy of a few parts per million. These experiments, performed mostly on atomic gases, require accurate knowledge of the dielectric virial coefficients that can be obtained from experiments and, more recently, from theory.

In the case of helium, the polarizability of a single atom can be calculated with unprecedented accuracy.17,18 In effect, a vast majority of the uncertainty in the determination of thermophysical properties comes from the error in the second-order and higher order virial coefficients.13 The second dielectric virial coefficient requires the knowledge of the two-body interaction potential and the two-body interaction-induced polarizability. The third dielectric virial coefficient requires the same quantities and, additionally, the knowledge of three-body properties, such as the non-additive three-body interaction potential and the three-body interaction-induced polarizability.

Two- and three-body interaction potentials and two-body interaction-induced polarizabilities for noble gas atoms with their contributions to the virial coefficients were already reported in the literature.19–25 On the contrary, the effects of three-body interaction-induced polarizabilities were studied only experimentally for various systems,26–30 and no complete theoretical description of the three-body polarizability surface is available in the literature. The first attempt to develop a theoretical model describing the three-body polarizabilities was made by Buckingham and Hands.31 They used the dipole-induced-dipole (DID) model but were unable to assess its accuracy as no sufficiently accurate three-body polarizability calculations were available at that time. Their comparison to the self-consistent electron pair calculations (SCEP) of Pérez et al.32 showed a substantial difference between the DID model and the ab initio results. A few years later, Champagne et al.33 used the non-local response theory34,35 and derived analytical behavior of the three-body interaction-induced polarizability for inert gases and simple molecules valid when all distances between interacting atoms or molecules are large. Accounting for the DID interaction, quadrupolar induction, dispersion, and induction–dispersion effects, they obtained all asymptotic terms behaving at large distances as r12lr23mr13n, with l + m + n ≤ 9, where rij are interatomic or intermolecular distances. Unfortunately, they did not compare their results with finite distance calculations of Pérez et al.,32 which were, to the best of our knowledge, the only ab initio calculations of the three-body polarizability of helium available at that time in the literature. Other ab initio calculations of the three-body polarizability were performed by Aleman et al.36 at the SCF level of theory and concerned neon atoms.

For the purpose of calculating the third dielectric virial coefficient, Heller and Gelbart37 suggested the so-called superposition approximation to obtain an estimate of the three-body polarizability using only two-body polarizabilities. However, they used crude two-body polarizability and interaction potential for helium in their third dielectric virial estimation. Recently, Garberoglio et al.38 used this approximation in Path-Integral Monte Carlo (PIMC)39–42 calculations to obtain the third dielectric virial coefficient for helium employing the more accurate two-body polarizability reported by Cencek et al.21 The values of the third dielectric virial coefficient for helium were also given in the paper of Alder et al.43 However, as noted in Ref. 38, these authors arbitrarily adjusted the two-body polarizability in order to improve the agreement with the experimental data.

In this work, we report accurate ab initio calculations of the three-body polarizability for helium performed using the coupled-cluster theory with an approximate account of three-electron excitations (the CC3 method).44 We also present an accurate analytical fit of the calculated three-dimensional isotropic polarizability surface. For metrology applications, we provide a three-dimensional fit of the corresponding local uncertainties of this surface. Finally, the reported polarizability surface is employed in classical and semiclassical calculations to obtain the third dielectric virial coefficient of helium. When applied in fully quantum path-integral calculations,45 our three-body polarizability surface was utilized in the development of a new primary pressure standard based on capacitance measurements, which recently passed the stress test with 2 parts per million (ppm) accuracy level achieved at 7 MPa.46 

The components αμν of a polarizability tensor α are defined as

αμν=2EFμFνF=0,
(1)

where μ, ν ∈ {x, y, z} and E is the energy of a considered quantum state in the presence of a static and uniform electric field F with components Fμ. In this work, we consider a system of three identical closed-shell atoms at fixed positions in space specified by vectors ri, i = 1, 2, 3. The energy and, consequently, the components of the polarizability tensor depend on these atomic positions. For weakly interacting systems, like helium atoms, it is useful to represent this dependence in terms of the so-called many-body expansion

α(r1,r2,r3)=3α0e+α2(r1,r2)+α2(r2,r3)+α2(r1,r3)+α3(r1,r2,r3),
(2)

where α0 is the atomic polarizability, e is the unit tensor, α2(ri, rj) is the two-body interaction-induced polarizability tensor of the pair of atoms at ri and rj,

α2(ri,rj)=α(ri,rj)2α0e,
(3)

and α3(r1, r2, r3) is the three-body interaction-induced polarizability tensor defined essentially in Eq. (2). By α(ri, rj), we denote the tensor of the total polarizability of two atoms located at ri and rj. The tensor α3(r1, r2, r3), representing the pair-wise non-additive part of the total interaction-induced polarizability tensor α(r1, r2, r3) − 3α0e, can be calculated from Eq. (2) using precomputed two-body polarizability tensors α2(ri, rj). In this work, it has been calculated directly from the formula

α3(r1,r2,r3)=α(r1,r2,r3)α(r1,r2)α(r2,r3)α(r1,r3)+3α0e.
(4)

The assumption that we are dealing with closed-shell systems is essential since for three interacting open-shell atoms, the quantum states of subsystems needed to determine the two-body contributions are not well defined (the simplest example is the ground, doublet state of three hydrogen atoms). Formula (2) is particularly useful when the three-body part is much smaller than the sum of two-body polarizabilities. This is always true when the interatomic distances are large, but in the short-range region, the three-body contribution can be large or even dominant.

For three atoms, the three-body polarizability tensor contains four independent components that can be used to define their combinations of physical interest.47 In this work, we restrict ourselves to the mean (isotropic) polarizabilities αn, n = 2, 3,

αn=13Trαn,
(5)

because of their relevance for calculations of the dielectric virial coefficients.

While the components of the polarizability tensor depend on the spatial orientation of the triangle formed by the atoms, the isotropic polarizability α3 is invariant to rotations of the system as a whole and depends only on the interatomic distances. When two or three interatomic distances are large, one can derive asymptotically valid analytic expressions for the dependence of α3 on r12, r23, and r13. While only one fragmentation is possible for two-body systems, a three-body system can separate through two distinct fragmentation channels. The first channel that we shall refer to as the atomic channel corresponds to the situation when all three atoms move away independently to infinity (all three interatomic distances are arbitrarily large). The second channel, referred to as the atom–diatom channel, describes the situation when one atom moves to infinity while the remaining two atoms stay at a limited, small distance from each other (only two interatomic distances are arbitrarily large).

The leading term of the asymptotic expansion of the isotropic three-body polarizability in the first fragmentation channel is predicted by the dipole-induced-dipole (DID) model31,33 and takes the form

αDID=2α03S123P2(cosθ3)r133r233,
(6)

where θ3 is the angle between vectors r1r3 and r2r3, Pn(x) is the Legendre polynomial of order n, and S123 denotes the operator that symmetrizes the expression it precedes with respect to the variables r1, r2, and r3. The classical expression for the third dielectric virial coefficient requires the integration of Eq. (6) over r13, r23, and θ3. It is clear that this integration is ill-defined as it is divergent both at small and at large values of r13 or r23. In Sec. III C and in the  Appendix, we show how this divergence is eliminated by an appropriate regularization. As a result, the inclusion of Eq. (6) in the analytic representation of the three-body polarizability leads to a finite result.

The leading term of the asymptotic expansion in the second (atom–diatom) fragmentation channel can be obtained from the results of Hunt et al.48 concerning the asymptotics of the interaction-induced polarizability of atom–diatom pairs. It suffices to substitute the value of the polarizability anisotropy β of the diatom with the polarizability anisotropy of the helium dimer and to subtract the isotropic interaction-induced polarizability of the two remaining helium pairs. This procedure leads to the following formula:

αa-d=23α0S123β(r)P2(cosθ)R3,
(7)

where r, R, and θ are the Jacobi coordinates: r and R are the lengths of the vectors r = r2r1 and R = r3 − (r1 + r2)/2, respectively, and θ is the angle between the vectors r and R. By β(r), we denote the polarizability anisotropy of two helium atoms separated by a distance r. Since the function being symmetrized here [and also in Eq. (6)] is already symmetric with respect to the permutation of r1 and r2, the symmetrizer S123 can be replaced by the operator 2+2P13+2P23, where the operator Pij permutes the vectors ri and rj.

It is easy to see that when one of the interatomic distances is small, Eq. (6) cannot represent the exact long-range behavior of the three-body polarizability in the atom–diatom fragmentation channel given by Eq. (7). If the long-range terms included in the analytical representation of the three-body polarizability were given only by Eq. (6) (and the higher-order terms derived in Ref. 33), then the long-range behavior exhibited in Eq. (7) would have to be modeled by fitting with short-range functions. This is a very impractical, if not impossible, procedure. Therefore, the asymptotics of Eq. (7) must be included in the analytic function used to represent the three-body polarizability. However, it must be realized that there are regions in the configuration space where Eqs. (6) and (7) give the same, or very close, results. Indeed, since for a pair of rare gas atoms, the anisotropy β(r) behaves at large r as 6α02/r3,49 Eqs. (6) and (7) give the same values when both r and R become large. Thus, we conclude that the inclusion of both Eqs. (6) and (7) in an analytic representation of the three-body polarizability would lead to double counting in some parts of the configuration space. A solution of this difficulty is presented in Sec. III B.

We performed ab initio calculations of the three-body isotropic polarizability α3 of helium for 748 triangles parameterized by grid points in the three-dimensional region defined by the triangle conditions for the coordinates r12, r23, and r13. The procedure by which the set of grid points was generated is described in detail in Sec. III B. The calculated polarizabilities were obtained using the CC3 method44 as implemented in the Dalton 2018 package.50,51 For any given atomic configuration, individual components of the three-body polarizability tensor were calculated using Eq. (4). All quantities in this formula were obtained using the basis set of the trimer, which is equivalent to applying the so-called counterpoise scheme to remove the basis set superposition error (BSSE).52,53 Note that, within this approach, the polarizability tensor of a single helium atom is no longer diagonal and depends slightly on the position of the atom in the triangle. Therefore, the atomic term 3αμν in Eq. (4) must be replaced by the sum of three, in general different, terms coming from each atom. Similarly, the two-body polarizabilities depend slightly on the shape of the triangle.

All calculations were performed using a family of doubly augmented correlation-consistent Gaussian basis sets, referred to further as dXZ, where X is the cardinal number, developed in Ref. 54 specifically for interactions of helium atoms in their ground state. The most important difference between the dXZ basis sets and the standard Dunning’s d-aug-cc-pVXZ basis sets for helium55,56 is that, in the former, the augmenting functions, optimized for the coefficients in the asymptotic expansion of the pair potential, are more diffuse (have smaller exponents) than in the latter, especially for angular momenta l ≥ 2. Therefore, the dXZ basis sets are more suitable for calculations of molecular properties, such as polarizability, which require a proper description of the electronic wave function at large distances from the nuclei. The highest cardinal number of the employed basis sets was X = 5. Unfortunately, the polarizability calculations using the d5Z basis set were computationally too demanding to be routinely used for all atomic configurations. This is especially true for general geometries where the plane determined by the positions of three non-colinear helium atoms is the sole symmetry element of the system. Therefore, we calculated only 45 points within the d5Z basis set mostly for geometries with a higher symmetry where three atoms either form an isosceles (including an equilateral) triangle or are located on a straight line.

In Table I, we present results obtained for several representative geometries. For a vast majority of atomic configurations, the dQZ values are essentially converged and are very close to the d5Z results. For example, for the equilateral triangle with side lengths equal to 5.73 a.u. (a geometry close to the minimum of the interaction potential of the helium trimer), we observe a difference of only 0.7% between the two basis sets. A similar behavior was found for other atomic configurations.

TABLE I.

Basis set convergence of the isotropic three-body polarizability α3 of helium calculated at the CC3 level of theory using the dXZ family of basis sets. α3[TQ] is the extrapolated polarizability defined in Eq. (8). Interatomic distances, rij, are in a.u., and the polarizability unit is in 10−3 a.u. The estimated absolute percentage uncertainty, 100% × |(α3[TQ] − α3[Q])/α3[TQ]|, is given in parentheses.

r12r23r13X = TX = QX = 5α3[TQ]
5.73 5.73 5.73 −0.040 65 −0.046 49 −0.046 83 −0.050 76 (8.40%) 
5.87 5.87 5.87 −0.059 41 −0.063 95 −0.064 16 −0.067 26 (4.93%) 
6.00 6.00 6.00 −0.067 47 −0.071 08 −0.071 22 −0.073 72 (3.57%) 
6.50 6.50 6.50 −0.061 80 −0.063 39 −0.063 41 −0.064 56 (1.80%) 
3.60 3.87 6.89 −0.228 58 −0.217 22 −0.209 85 −0.208 94 (3.96%) 
2.20 5.30 7.50 1.234 10 1.336 80 1.356 30 1.411 70 (5.31%) 
5.00 5.00 10.00 0.590 28 0.599 34 0.598 77 0.605 95 (1.09%) 
7.95 8.23 10.10 −0.014 51 −0.014 52 −0.014 50 −0.014 53 (0.04%) 
r12r23r13X = TX = QX = 5α3[TQ]
5.73 5.73 5.73 −0.040 65 −0.046 49 −0.046 83 −0.050 76 (8.40%) 
5.87 5.87 5.87 −0.059 41 −0.063 95 −0.064 16 −0.067 26 (4.93%) 
6.00 6.00 6.00 −0.067 47 −0.071 08 −0.071 22 −0.073 72 (3.57%) 
6.50 6.50 6.50 −0.061 80 −0.063 39 −0.063 41 −0.064 56 (1.80%) 
3.60 3.87 6.89 −0.228 58 −0.217 22 −0.209 85 −0.208 94 (3.96%) 
2.20 5.30 7.50 1.234 10 1.336 80 1.356 30 1.411 70 (5.31%) 
5.00 5.00 10.00 0.590 28 0.599 34 0.598 77 0.605 95 (1.09%) 
7.95 8.23 10.10 −0.014 51 −0.014 52 −0.014 50 −0.014 53 (0.04%) 

The results obtained with finite basis sets were extrapolated to the complete basis set (CBS) limit. We assumed that the basis set truncation error of the calculated polarizability vanishes with the increasing value of the basis sets cardinal number X as 1/Xn. For the value of the exponent n, we adopted n = 3, as recommended for the extrapolation of the correlation energies.57,58 This choice can be justified by the fact that the polarizability is defined as the energy derivative, see Eq. (1), and its calculation involves no singular operators that could modify the value of n. Such modifications appear, for example, in the extrapolation of relativistic corrections.54,59,60 With this choice of the exponent n, the CBS limit can be obtained using the following two-point formula:

α3[X1,X]=α3[X]+(X1)3α3[X]α3[X1]X3(X1)3,
(8)

where α3[X − 1] and α3[X] are polarizabilities calculated using basis sets with cardinal numbers X − 1 and X, respectively, and α3[X − 1, X] is the extrapolated value. As the X = 5 results were not available for all considered geometries, we consistently assumed the α3[TQ] extrapolants as our recommended values in the CBS limit.

To each geometry of the trimer, we assigned a theoretical uncertainty of the recommended value of the polarizability. This uncertainty was estimated as the magnitude of the difference between the extrapolated value and the value computed with the largest basis set used in the extrapolation, i.e., σ[TQ] = |α3[TQ] − α3[Q]|. The results obtained using this approach are given in Table I. In some cases, however, this simple procedure may lead to an underestimation of the uncertainty. This occurs when the dTZ value is accidently close to the dQZ one, as can be observed, for example, for the triangle with side lengths (7.95, 8.23, 10.10) a.u. In this case, the estimated uncertainty is only about 0.04% of the extrapolated polarizability, which is half of the difference between the X = 4 and X = 5 results. To avoid underestimation, any uncertainty estimated to be lower than 0.1% of the extrapolated value was scaled up to a conservative value of 1.0%. From now on, we use the notation σ[TQ] to refer to the modified uncertainty. On average, the modified estimated uncertainty is 4.7%.

To obtain a more reliable estimation of the theoretical uncertainty, we assessed the importance of the error of the CC3 model that results from an approximate treatment of triple excitations, as well as the neglect of quadruple and higher excitations. To this end, the full configuration interaction (FCI) calculations of the polarizability were performed for 71 grid points using a program written for the purposes of this project.61 Due to the high computational demands of the FCI method, the calculations were restricted to geometries of high symmetry and were performed using a single basis set of triple-zeta quality. As the dTZ basis set from Ref. 54 was still too large to accomplish the calculations, we prepared its singly augmented version, further referred to as aTZ. The results for six selected geometries are presented in Table II. We found that the error resulting from an approximate character of the CC3 model is almost negligible for small triangles, only about 0.03% of α3. By contrast, for larger triangles, it becomes comparable to the extrapolation errors, σ[TQ]. The average absolute FCI contribution across all 71 grid points is 5.7% of α3[TQ]. Therefore, the value 5.7% was taken as a global, geometry independent uncertainty due to the error of the CC3 model. The total uncertainty σ was obtained by multiplying α3[TQ] by this value of the relative error and adding the result in squares to σ[TQ], namely,

σ=σ[TQ]2+0.057×α3[TQ]2.
(9)
TABLE II.

Comparison of the three-body polarizability of helium calculated at the CC3 (α3CC3) and FCI (α3FCI) levels of theory using the aTZ basis set. Interatomic distances, rij, are in a.u. and the polarizability unit is in 10−3 a.u. δ=100%×|(α3FCIα3CC3)/α3CC3|.

r12r23r13α3CC3α3FCIδ
2.50 2.50 2.50 210.809 2 210.768 6 0.02% 
2.75 2.75 2.75 151.949 8 151.908 1 0.03% 
2.60 5.40 8.00 0.876 54 0.922 74 5.27% 
7.56 7.56 5.43 −0.074 77 −0.077 45 3.58% 
6.20 6.20 5.50 −0.070 47 −0.075 38 6.97% 
6.50 6.50 6.50 −0.059 93 −0.062 46 4.21% 
r12r23r13α3CC3α3FCIδ
2.50 2.50 2.50 210.809 2 210.768 6 0.02% 
2.75 2.75 2.75 151.949 8 151.908 1 0.03% 
2.60 5.40 8.00 0.876 54 0.922 74 5.27% 
7.56 7.56 5.43 −0.074 77 −0.077 45 3.58% 
6.20 6.20 5.50 −0.070 47 −0.075 38 6.97% 
6.50 6.50 6.50 −0.059 93 −0.062 46 4.21% 

The final average total uncertainty of our polarizability surface is 9.1% and its median is 5.8%. Using the local, three-dimensional analytical fits of α3[TQ] and σ[TQ], see Sec. III B, the total uncertainty σ obtained from Eq. (9) was employed in Sec. III C to calculate the uncertainty of the thermophysical quantities such as the third dielectric virial coefficient.

To the best of our knowledge, the only theoretical determination of the three-body polarizability of helium has been reported by Pérez et al.32 using the SCEP method. Our results differ substantially from Ref. 32. For example, the ratios of the SCEP and our CC3 results are 1.2, 0.5, and 13 for the equilateral triangles with side lengths equal to 2, 4, and 6 a.u., respectively. For larger isosceles and linear triangles, the differences between the SCEP and CC3 results decrease, especially for linear configurations. The SCEP results of Pérez et al. for larger equilateral triangles do not provide sufficient accuracy for comparison.

In Fig. 1, we compare our ab initio data for equilateral triangles with the theoretical asymptotics of Champagne et al.33 and the superposition approximation.37,38 We omit the results of Pérez et al.32 as their values lie outside of the range of the figure. It is known that the superposition approximation includes the DID term from Eq. (6) at long range. This is confirmed in Fig. 1 where the two models give identical results for large triangles. The DID term and the superposition approximation start to diverge from each other for triangles with side lengths below 6.5 a.u. when exponentially decaying terms in the two-body polarizabilities start contributing significantly to the superposition approximation of the polarizability. It is worth noting that, in contrast to the pure DID term, the superposition approximation predicts a minimum of the three-body polarizability. The differences between the positions of the minima obtained from the ab initio calculations and from the superposition approximation are about 1 a.u. for equilateral triangles. When the interatomic distances are not sufficiently large, the DID model and, consequently, the superposition approximation differ significantly from the results of the ab initio calculations. For example, a reasonable agreement for equilateral triangles is found only when the side lengths are above 9 a.u. (see Fig. 1). Clearly, additional asymptotic terms derived in Ref. 33 are required to correctly describe the polarizability down to 7 a.u. or so. For equilateral triangles with side lengths below 7 a.u., the asymptotic expansion fails completely due to the increasing importance of the exchange and overlap effects. We observed a similar behavior also for isosceles configurations. The superposition approximation fails for isosceles triangles with the base side shorter than 7 a.u. However, it is surprisingly accurate for isosceles triangles with the base longer than 7 a.u. even if the leg sides are shorter than 7 a.u.; see Fig. 2 for representative examples. Based on these observations, we can assume that the transition to the long-range regime takes place for triangles with the shortest side length equal to ∼7 a.u. This is the distance where the pair interaction potential of helium reaches ∼40% of its well depth located at the interatomic separation of 5.6 a.u. We expect similar behavior for other noble gases.

FIG. 1.

Comparison of asymptotic models and extrapolated α3[TQ] ab initio data for the equilateral triangles with side R.

FIG. 1.

Comparison of asymptotic models and extrapolated α3[TQ] ab initio data for the equilateral triangles with side R.

Close modal
FIG. 2.

Comparison of asymptotic models and extrapolated α3[TQ] ab initio data for the isosceles triangles with leg side R for two values of fixed base side length Rbase. (a) Rbase = 6 a.u. (b) Rbase = 8 a.u.

FIG. 2.

Comparison of asymptotic models and extrapolated α3[TQ] ab initio data for the isosceles triangles with leg side R for two values of fixed base side length Rbase. (a) Rbase = 6 a.u. (b) Rbase = 8 a.u.

Close modal

The three-body polarizability surface α3(r12, r23, r13) was fitted with an analytic function αfit(r12, r23, r13) expressed as the sum of three components, one of a short-range and two of a long-range character,

αfit(r12,r23,r13)=αsr(r12,r23,r13)+α3a(r12,r23,r13)+αa-d(r12,r23,r13).
(10)

The short-range function αsr(r12, r23, r13), representing the exponentially decaying contributions to α3(r12, r23, r13), has the form

αsr(r12,r23,r13)=l=14ijkAijk,lS123eαlr12βlr23γlr13r12ir23jr13k,
(11)

where Aijk,l are linear and αl, βl, and γl are nonlinear parameters to be optimized. Possible values of the indices i, j, and k are restricted by the inequalities 0 ≤ ijk ≤ 5 and the additional condition imposed on their sum, namely, i + j + k ≤ 7. Overall, αsr includes 112 linear parameters and 12 nonlinear ones. The form of Eq. (11) differs from the exponentially decaying terms employed in Refs. 20 and 25, where the perimeter of the triangle was used as the main argument of the exponential function. The resulting fitting function incorrectly describes linear geometries since it remains constant when one atom moves between two atoms that have fixed positions. By contrast, the exponential terms used in Eq. (11) are able to describe this movement correctly. Another difference is that we use the interatomic distances and their powers to describe the angular dependence of the potential, rather than the Legendre polynomials of the angles between triangle sides used in Refs. 20 and 25.

The α3a term in Eq. (10) provides the correct asymptotics of α3 in the three-atom fragmentation channel. It is given by the damped asymptotic expansion derived in Ref. 33,

α3a(r12,r23,r13)=2Z1S123P2(cosθ1)r123r133f4(η1r12)f4(η1r13)+10Z2S123P3(cosθ1)r124r134f5(η2r12)f5(η2r13)6Z3(1+3cosθ1cosθ2cosθ3)r123r233r133f4(η3r12)f4(η3r23)f4(η3r13)+4Z4S123P2(cosθ1)r123r136f4(η4r12)f7(η4r13)6πZ5(1+3cosθ1cosθ2cosθ3)r123r233r133f4(η5r12)f4(η5r23)f4(η5r13)+4πZ6S123P2(cosθ1)r123r136f4(η6r12)f7(η6r13),
(12)

where fn are the Tang–Toennies damping functions62 

fn(x)=1ex1+x+x22!++xnn!.
(13)

The Zn coefficients can be calculated with the knowledge of several atomic properties: the static dipole and quadrupole polarizabilities (Z1, …, Z4), the dynamic dipole polarizabilities (Z5 and Z6), and the second dipole hyperpolarizabilities (Z5 and Z6) (see Ref. 33). Their values, calculated using the data published in Refs. 33 and 63, were fixed during the fitting procedure, while the nonlinear parameters η1, …, η6 were fully optimized. The individual terms in Eq. (12) have a clear physical interpretation.33 These arise from the second-order induction interaction (the first two terms), the third-order induction interaction (the next two terms), the dispersion interaction (the fifth term), and the induction–dispersion interaction (the sixth term).

The exponentially decaying short-range part of the fitting function, αsr, is unable to satisfactorily reproduce the calculated data points (within the estimated uncertainty) in the region when two atoms are relatively close to each other and the third atom is at a large distance. To properly describe such configurations, we include an additional function αa−d in Eq. (10). This term is an adaptation of Eq. (7) representing the asymptotics of α3 in the atom–diatom fragmentation channel. To avoid the problem of possible double-counting mentioned in Sec. II, we substituted the exact diatomic polarizability anisotropy β(r) in Eq. (7) with a short-range exponentially decaying polynomial expansion,

αa-d(r12,r23,r13)=S̃123P2(cosθ)(C0+C1r+C2r2)eη7rR3f4(η8R),
(14)

where S̃123=1+P13+P23 generates the sum over three possible ways of defining the Jacobi coordinates r, R, and θ for a triangle with sides r12, r23, and r13. When r is small, Eq. (12) gives a negligible contribution due to the strong Tang–Toennies damping and Eq. (14) provides the correct atom–diatom asymptotics. Conversely, when r is large, Eq. (14) gives a negligible contribution and Eq. (12) guarantees the correct asymptotics in the three-atom fragmentation channel. It may be noted that the right-hand-side of Eq. (14) has a singularity at r = 0. This singularity is irrelevant from the physical point of view and we found that keeping C0 and C1 different from zero improves the quality of the fit. The parameters C0, C1, C2, η7, and η8, adjusted in the fitting procedure, are controlling the switching between Eqs. (12) and (14). The total number of fitted parameters in Eqs. (11), (12), and (14) is 135, including 115 linear parameters and 20 nonlinear parameters.

Effective selection of grid points is an important problem in fitting potential energy surfaces.64 In our work, the set of grid points was generated using an iterative approach. An initial set of 250 points, satisfying the conditions min(r12, r23, r13) > 1.2 a.u. and max(r12, r23, r13) < 14 a.u., was created through the Sobol sequence. In each step of the iterative procedure, ∼106 sets of 20 nonlinear parameters were randomly generated in a suitably chosen 20D box. Subsequently, the current set of grid points was used in a local optimization of the nonlinear parameters to obtain ∼106 independent fitting functions. The optimization was performed using the trust-region dogleg procedure65 as implemented in the GSL library.66 During the optimization, the linear parameters (for fixed values of the nonlinear ones) were obtained using the weighted linear least-squares minimization with weights equal to the inverse square of the calculated uncertainties, 1/σ[TQ]2. The fits prepared this way were sorted according to their average relative errors and maximum relative errors, and the top 20 fits were picked for further analysis. The analysis consisted of evaluating the polarizability for ∼104 randomly distributed trimer configurations and finding regions of the configuration space where the largest differences between the preselected 20 fits were observed. A new batch of 50–100 grid points covering the problematic regions were constructed, the corresponding polarizabilities were computed, and the whole procedure was repeated with the enlarged set of grid points. The generation process was finished after six iterations when the differences between trial fits were acceptably small. The final batch of points was added to probe two particular regions: isosceles triangles for which the DID term in Eq. (6) is zero and configurations where the atom–diatom contribution is important, i.e., triangles with one short side. Altogether, the set of 748 grid points was generated. The final fit was prepared using a training set of 672 points (∼90% of the total number of points). The remaining 76 points (∼10% of the total number of points) formed a testing set and were used to assess the accuracy of the fit for points outside of the training set.

The final fit has a mean absolute relative error with respect to uncertainties due to the extrapolation of 0.42σ[TQ] and its median is 0.28σ[TQ]. The average relative error with respect to the values of the polarizability is 1.4% and its median is 0.2%. A majority of the configurations have a relative error smaller than σ[TQ] (see Fig. 3). Most of the points with relative errors larger than σ[TQ] are encountered for triangles with at least one short side (<2 a.u.). In such cases, the polarizability is several orders of magnitude larger than in other regions. For example, the polarizability calculated for the triangle (1.4, 4.6, 6.0) is −91.48 × 10−3 a.u. and the relative error is 4.8σ[TQ]. On the other hand, for a similar triangle (2.0, 5.0, 6.9), the polarizability is −0.73 × 10−3 a.u. Fortunately, these configurations were found to have a minimal influence on the quality of the final results. For example, in the calculation of the third dielectric virial coefficient (see Sec. III C), triangles with one side shorter than 2 a.u. were not important even at the temperature of 2400 K.

FIG. 3.

Absolute differences between calculated and fitted three-body polarizabilities α3 vs the estimated uncertainties σ[TQ]. The straight dotted lines correspond to the relative errors.

FIG. 3.

Absolute differences between calculated and fitted three-body polarizabilities α3 vs the estimated uncertainties σ[TQ]. The straight dotted lines correspond to the relative errors.

Close modal

The second group of triangles with a relative error larger than σ[TQ] are triangles with all sides longer than 16 a.u. For such large interatomic distances, the values of the polarizability are tiny and, hence, are plagued by a large cancellation of significant digits when the supermolecular approach of Eq. (4) is used. However, these regions of the surface are completely determined by the asymptotic expansion. Therefore, the analytical expansion coefficients used in our fit are expected to provide accurate results.

In order to estimate uncertainties of the thermophysical properties of gaseous helium, which can be obtained using our potential, we provide a fit of σ(r12, r23, r13) reflecting the estimated local uncertainties of our ab initio calculations. Except for a small number of outliers, the exact values of the polarizability are expected to be contained within the range α3(r12, r23, r13) ± σ(r12, r23, r13). The function σ(r12, r23, r13) is not intended to precisely fit the uncertainties, but rather to follow the trend with respect to the size of the triangles and provide upper bounds to the calculated uncertainties for the thermophysical properties. The analytical function used to obtain the fit of the uncertainties is analogous to the short-range part of Eq. (11), where now the summation over l goes from 1 to 2, and the allowed values of the indices i, j, and k are restricted by conditions 1 ≤ ijk ≤ 4 and i + j + k ≤ 8. The fit contains 26 linear and 6 nonlinear parameters.

The fitting was performed for the total uncertainties as defined in Eq. (9). Due to a large variation of the estimated uncertainties across the three-dimensional space of the interatomic distances, we selected a smaller subset of the uncertainties for the fitting purposes. First, we discarded configurations with the uncertainties significantly smaller than for the neighboring ones. This is justified by the observation that the resulting fit then bounds these uncertainties from above, which is sufficient for our purposes. Next, we discarded triangles with interatomic distances min(r12, r23, r13) > 16 a.u. as the asymptotic expansion is expected to provide the correct results in this case. Finally, we discarded the smallest triangles with min(r12, r23, r13) < 2 a.u., which are of little importance in calculations of the physical properties. The final set of points used in the construction of σ(r12, r23, r13) comprised almost 65% of the total number of configurations. The average ratio of the values of the fitting function σ(r12, r23, r13) to the estimated uncertainties for the whole set of points is about 1.3. The values of all linear and nonlinear parameters of α3(r12, r23, r13) and σ(r12, r23, r13), along with a Fortran 2003 implementation of the fitting functions, are included in the supplementary material.

As mentioned in Sec. I, the third dielectric coefficient Cɛ plays an important role in the precise determination of various metrological quantities. It is defined by the expansion of the Clausius–Mossotti function for the dielectric constant ɛ in powers of the gas density ρ,67,68

ε1ε+2=ρAε+ρ2Bε+ρ3Cε+,
(15)

where Aɛ is proportional to the atomic polarizability and Bɛ is the second dielectric virial coefficient.

The classical expression for the third dielectric virial coefficient Cɛ(T) was first published in Ref. 68. In this work, we use a factorized expression analogous to the one used in Refs. 40, 69, and 70 for the third pressure virial calculation. Within this framework, the third virial dielectric coefficient is given by the sum of three terms,

Cε(T)=4B(T)Bε(T)+Cε2body(T)+Cε3body(T),
(16)

where B(T) is the second pressure virial coefficient,

Cε2body(T)=2π9VS̃123α2(r1,r2)eβU(r1,r2,r3)eβU2(r1,r2)dr1dr2dr3,
(17)

and

Cε3body(T)=2π9Vα3(r1,r2,r3)eβU(r1,r2,r3)dr1dr2dr3.
(18)

In the above formulas, U(r1, r2, r3) = U3(r1, r2, r3) + U2(r1, r2) + U2(r2, r3) + U2(r1, r3) is the sum of the two-body, U2, and three-body, U3, interaction potentials, and β = 1/kBT. Note that the integral Cε2body(T) involves not only the two-body quantities (polarizability and interaction potential) but also the three-body interaction potential U3. The calculation of the first term in Eq. (16) is trivial since accurate values of B(T) and Bɛ(T) are available in the literature.71 The computation of the second and third terms is somewhat more involved. The translational and rotational invariance can be used to reduce the nine-dimensional integrals in Eqs. (17) and (18) to three-dimensional integrals over two sides of the triangle and the angle between them using the substitution V1dr1dr2dr38πr132r232dr13dr23dcosθ3. Once this exchange of variables is performed, the evaluation of Eq. (18) reduces to the integration of the right-hand-side of Eq. (6) over r13, r23, and θ3. However, the resulting integral is ill-defined as the integrals over r13 and r23 are divergent at infinity, even if the singularities at r13 = 0 and r23 = 0 are removed using suitable damping factors. In the  Appendix, we show that, after an appropriate regularization, this divergence is eliminated and the required integral is actually equal to zero. This does not mean that the DID contribution to Cϵ(T) vanishes. The inclusion of Eq. (6) in the fitting function is necessary to obtain an accurate fit at intermediate distances, which are important in evaluation of Cϵ(T). Moreover, Eq. (6) contributes directly to Cϵ(T) when multiplied by Boltzmann factors dependent on r12. As the classical integration is valid only for higher temperatures, we also calculated the semiclassical results using the Feynman–Hibbs modified two-body potential and polarizabilities.39–41 The results of the semiclassical calculations and the details of the calculations are presented in Sec. IV.

The integrations were performed for all possible triangle configurations with sides up to 8 nm using the adaptive Gauss–Kronrod quadrature of degree 7 and 15.72 In the calculations of the two-body contributions we used the two-body interaction potential reported by Czachorowski et al.19 and the two-body polarizability of Cencek et al.21 For the three-body contribution, we used the three-body polarizability presented in this work and the most recent three-body interaction potential in Ref. 73. Note that all of these potentials have known local uncertainties. We also tested the three-body potential from Ref. 20, but the difference of the computed Cε3body(T) using the alternative three-body interaction potential was several orders of magnitude smaller than the uncertainty originating from the accuracy of the three-body polarizability. The error of the numerical integration was assessed using the Gauss–Kronrod method and was much smaller than the effect of the uncertainties of the polarizability and the potential energy surfaces.

The uncertainties of the third dielectric virial coefficient were estimated through the propagation of errors for both the two-body and three-body polarizabilities and the three-body interaction potential. The errors of the two-body potential are negligible19 and were not included in the propagation. The uncertainties due to the errors in the polarizabilities and the three-body potential were computed using the approach developed by Garberoglio and Harvey74 for calculations of the pressure viral coefficients. The uncertainties due to the polarizabilities were estimated employing the following formula:

|Cε(αnacc)Cε(αn)|=(αnaccαn)δCεδαndXσ(αn)δCεδαndX,
(19)

where αnacc and αn are the accurate and approximate n-body polarizability functions, respectively, δCɛ/δαn is the functional derivative of Cɛ with respect to αn, σ(αn) is our estimate of the k = 2 level uncertainty of αn, and dX is the volume element in the space of coordinates on which αn depends. The equality in Eq. (19) is exact since Cɛ depends linearly on αn. The inequality in Eq. (19) holds at the 95% confidence level since |Cε(αnacc)Cε(αn)|σ(αn) and since σ(αn) is interpreted as the k = 2 level uncertainty of αn. The rightmost side of Eq. (19) was also used to estimate a minor source of uncertainty due to the approximate nature of the three-body potential energy U3. In this case, the equality in Eq. (19) is valid only approximately with an error quadratic in σ(Un). The total uncertainty was computed as the square root of the sum of squares of the uncertainties of α2, α3, and U3.

The results of the classical calculations of Cɛ(T) and of the contributions identified in Eq. (16) are presented in Fig. 4. In contrast to the third pressure virial coefficient,40 the three-body contribution Cε3body(T) to the third dielectric virial coefficient has a comparable magnitude to the combined two-body terms. For temperatures of about 1270 K, the values of Cε3body(T) cross zero and the two-body terms dominate. However, this trend suggests an increasing importance of Cε3body(T) for temperatures higher than 2500 K.

FIG. 4.

Comparison of the classical two-body and three-body contributions to Cɛ.

FIG. 4.

Comparison of the classical two-body and three-body contributions to Cɛ.

Close modal

From the inspection of the contributions to Cε3body(T) originating from the three terms included in the fitting function, see Eq. (10), we conclude that for higher temperatures, the short-range contributions dominate. This is an expected behavior since the helium atoms have enough kinetic energy to penetrate the repulsive wall of the potential. As these short-range contributions are mostly positive, they compensate the (typically negative) long-range contributions. Regarding the uncertainties, the uncertainty of Cε3body(T) dominates for the entire temperature range. The largest contribution to Cε3body(T) then comes from the uncertainty of the three-body polarizability.

The calculated values of Cε3body(T), as well as of the individual contributions defined in Eq. (16), are presented in Table III for a wide range of temperatures. In Figs. 5 and 6, we compare our calculations with theoretical results and experimental values found in the literature. To the best of our knowledge, there are only two relevant theoretical calculations of the third dielectric virial coefficient available in the literature. Both of them use the superposition approximation to obtain the three-body polarizability. Heller and Gelbart provided only one value of Cɛ(300 K) = −0.719 cm9 mol−3.37 In their work, they used relatively simple models for both the two-body potential and the two-body polarizability. More recently, Garberoglio et al.38 used the superposition approximation for the three-body polarizability and the same two-body potential and two-body polarizability as in this work. Their PIMC value of Cɛ(T) at 300 K is −0.556 cm9 mol−3,38 which is close to the value Cɛ(300 K) = −0.535 cm9 mol−3 obtained by us. Overall, our values of Cɛ(T) are negative over the whole range of temperatures, in agreement with the results of Garberoglio et al.38 The classical results are similar to their values up to T = 400 K, but the two curves start to diverge for higher temperatures. According to our results, a local minimum around 1500 K is present, see Fig. 5, while this feature is absent in the results of Garberoglio et al.38 However, note that the rate of decrease of Cɛ(T) with the temperature computed using the superposition approximation slows down around 2000 K. Therefore, it is possible that this may result in a minimum for higher temperatures. One can conclude that the difference between the results of Ref. 38 and our work is due to the incorrect description of the short-range polarizability within the superposition approximation.

TABLE III.

The temperature dependence of the third dielectric virial coefficient Cɛ(T) for helium and its many-body contributions. The temperature is in kelvin and the unit of a dielectric virial coefficient is cm9 mol−3. The uncertainties at k = 2 are shown in parentheses.

TCε3body4BBε+Cε2bodyCɛ
10 −0.158(56) −0.254(14) −0.412(57) 
11 −0.150(49) −0.163(11) −0.313(50) 
12 −0.145(44) −0.114 0(90) −0.259(45) 
13 −0.142(40) −0.087 2(76) −0.229(41) 
14 −0.139(37) −0.072 5(65) −0.212(38) 
15 −0.138(35) −0.064 9(57) −0.202(35) 
20 −0.133(29) −0.066 8(34) −0.199(29) 
25 −0.130(26) −0.084 7(24) −0.215(26) 
30 −0.129(24) −0.103 4(19) −0.232(24) 
35 −0.128(23) −0.120 9(16) −0.249(23) 
45 −0.126(22) −0.151 7(12) −0.277(22) 
50 −0.125(22) −0.165 2(11) −0.290(22) 
75 −0.121(22) −0.220 95(73) −0.342(22) 
100 −0.118(22) −0.263 49(58) −0.381(22) 
125 −0.114(23) −0.298 00(50) −0.412(23) 
150 −0.111(23) −0.327 02(44) −0.438(23) 
175 −0.108(24) −0.352 02(40) −0.460(24) 
200 −0.105(25) −0.373 96(37) −0.479(25) 
273.15 −0.097(26) −0.425 73(32) −0.523(26) 
273.16 −0.097(26) −0.425 73(32) −0.523(26) 
300 −0.094(27) −0.441 34(30) −0.535(27) 
340 −0.090(28) −0.462 11(28) −0.552(28) 
400 −0.084(29) −0.488 86(26) −0.573(29) 
500 −0.074(31) −0.524 83(24) −0.599(31) 
1000 −0.032(38) −0.625 47(18) −0.657(38) 
1200 −0.017(41) −0.647 61(16) −0.664(41) 
1500 0.004(44) −0.671 34(15) −0.667(44) 
1700 0.018(46) −0.682 81(14) −0.665(46) 
2000 0.036(49) −0.695 51(14) −0.659(49) 
2400 0.059(53) −0.706 57(13) −0.647(53) 
TCε3body4BBε+Cε2bodyCɛ
10 −0.158(56) −0.254(14) −0.412(57) 
11 −0.150(49) −0.163(11) −0.313(50) 
12 −0.145(44) −0.114 0(90) −0.259(45) 
13 −0.142(40) −0.087 2(76) −0.229(41) 
14 −0.139(37) −0.072 5(65) −0.212(38) 
15 −0.138(35) −0.064 9(57) −0.202(35) 
20 −0.133(29) −0.066 8(34) −0.199(29) 
25 −0.130(26) −0.084 7(24) −0.215(26) 
30 −0.129(24) −0.103 4(19) −0.232(24) 
35 −0.128(23) −0.120 9(16) −0.249(23) 
45 −0.126(22) −0.151 7(12) −0.277(22) 
50 −0.125(22) −0.165 2(11) −0.290(22) 
75 −0.121(22) −0.220 95(73) −0.342(22) 
100 −0.118(22) −0.263 49(58) −0.381(22) 
125 −0.114(23) −0.298 00(50) −0.412(23) 
150 −0.111(23) −0.327 02(44) −0.438(23) 
175 −0.108(24) −0.352 02(40) −0.460(24) 
200 −0.105(25) −0.373 96(37) −0.479(25) 
273.15 −0.097(26) −0.425 73(32) −0.523(26) 
273.16 −0.097(26) −0.425 73(32) −0.523(26) 
300 −0.094(27) −0.441 34(30) −0.535(27) 
340 −0.090(28) −0.462 11(28) −0.552(28) 
400 −0.084(29) −0.488 86(26) −0.573(29) 
500 −0.074(31) −0.524 83(24) −0.599(31) 
1000 −0.032(38) −0.625 47(18) −0.657(38) 
1200 −0.017(41) −0.647 61(16) −0.664(41) 
1500 0.004(44) −0.671 34(15) −0.667(44) 
1700 0.018(46) −0.682 81(14) −0.665(46) 
2000 0.036(49) −0.695 51(14) −0.659(49) 
2400 0.059(53) −0.706 57(13) −0.647(53) 
FIG. 5.

Comparison of the classical values of Cɛ(T) calculated in this work with the results obtained using the superposition approximation by Garberoglio et al.38 

FIG. 5.

Comparison of the classical values of Cɛ(T) calculated in this work with the results obtained using the superposition approximation by Garberoglio et al.38 

Close modal
FIG. 6.

Comparison of the classical values of Cɛ(T) calculated in this work with the experimental data.75–78 

FIG. 6.

Comparison of the classical values of Cɛ(T) calculated in this work with the experimental data.75–78 

Close modal

Regarding the experimental values, only a few of them can be found in the literature.75–78 The most recent value of Gaiser and Fellmuth comes from the dielectric-constant gas thermometry.75 They obtained the third pressure virial coefficient from the measured isotherms and converted this quantity to the third dielectric virial coefficient with the use of the best available pressure virials54,79,80 and the second dielectric virial coefficient.81 Other experiments performed at room temperature are also consistent with our data. However, for low temperatures, the difference between the experiments and our theoretical calculations is substantial, even if we take large experimental errors into account. This suggests that more accurate measurements are needed to resolve this discrepancy.

It is also important to mention the missing contribution to Cɛ(T) due to the three-body dipole moment. To date, no dipole moment surface has been published and Garberoglio et al.38 used an approximate model containing only the long-range asymptotic description of Hunt and Li.35 They determined that the contribution due to the dipole moment is negligible. Our preliminary study of the three-body dipole moment confirms their findings. When expressed in atomic units, the three-body dipole moments computed by us using the CC3 theory are one to two orders of magnitude smaller than the corresponding three-body polarizabilities. Note that the dipole moments contribute to Cɛ(T) through squares of their length divided by kBT. Therefore, the contribution of the three-body dipoles to Cɛ(T) is expected to be smaller than the present uncertainty of Cɛ(T), except possibly for very low temperatures.

It is well known that the classical calculations of the virial coefficients fail to properly describe the properties of the gas at low temperatures when quantum effects play an important role.2,38,40,82–85 For this reason, we also consider the leading-order quantum correction and provide a semiclassical estimation of the third dielectric virial coefficient. The Wigner–Kirkwood expansion is widely used for this purpose,2,86,87 but the Feynman–Hibbs effective potentials have been recently used as well.82 As noted by Guillot and Guissani,88 the Feynman–Hibbs effective potentials are more accurate and have better convergence properties than the Wigner–Kirkwood expansion. Therefore, we used the Feynman–Hibbs effective potential in this work and modified the two-body interaction potentials with the quadratic Feynman–Hibbs (QFH) correction,39,88,89

U2QFH(rij)=U2(rij)+2β12m2U2(rij)rij2+2rijU2(rij)rij,
(20)

where m is the mass of the helium atom. The modified two-body potential is used instead of the original two-body potential in the calculation of the virial coefficients. The same modification was used by Shaul and collaborators to calculate semiclassical values of the pressure virial coefficients.82 

Similarly, it is possible to introduce the same modification to the two-body polarizability as the polarizability is the second derivative of the potential with respect to the electric field. The QFH modified polarizability reads

α2QFH(rij)=α2(rij)+2β12m2α2(rij)rij2+2rijα2(rij)rij.
(21)

Since there is no three-body Feynman–Hibbs effective potential available, we modified only the two-body potentials and polarizabilities. Using the modified two-body potential and polarizability from Eqs. (20) and (21), we calculated semiclassical results for the two-body terms (BBɛ, Cε2body) in Eq. (16). The resulting semiclassical values of the third dielectric virial coefficients are presented in Fig. 7.

FIG. 7.

Comparison of the classical and semiclassical values of Cɛ(T) calculated in this work with the results from the quantum calculation (PIMC) of Garberoglio et al.38 employing the superposition approximation.

FIG. 7.

Comparison of the classical and semiclassical values of Cɛ(T) calculated in this work with the results from the quantum calculation (PIMC) of Garberoglio et al.38 employing the superposition approximation.

Close modal

Note that this procedure is easily justified for the third pressure virial coefficient where the two-body contribution accounts for over 95% of its value. The ratio of the two- and three-body contributions to the third dielectric virial coefficient is closer to the unity (see Fig. 4), and hence, our semiclassical values may provide incomplete information. Nevertheless, this crude semiclassical approximation to the two-body terms of Cɛ(T) only gives results that agree with the PIMC values of the Garberoglio et al.38 down to 10 K. It is possible that this is a consequence of a fortuitous error cancellation, but it may also indicate that the semiclassical quantum correction to the three-body polarizability surface is small. This requires a further study employing the quadratic Feynman–Hibbs potential for three particles.

We have presented the first complete ab initio calculations of the three-body polarizability of helium and constructed an analytical three-dimensional fit of its isotropic part. The fit is based on ab initio calculations at the CC3 level of theory extrapolated to the complete basis set limit. We have taken into account the exact asymptotic behavior of the three-body polarizability for atomic configurations corresponding to two possible fragmentation channels: three free atoms and atom–diatom. Our uncertainty budget takes into account the uncertainty of the complete basis set extrapolation as well as the approximate nature of the CC3 model. The analytical fit obtained by us has a mean absolute relative error with respect to uncertainties due to the extrapolation of 0.42σ, which corresponds to the mean absolute percentage error of 1.4% with respect to the computed ab initio polarizabilities. Furthermore, we present a three-dimensional fit of local uncertainties of our calculations, which is important from the point of view of applications, e.g., in metrology.

Employing the fits developed in this work and the best available two-body potentials and polarizabilities, we performed classical and semiclassical calculations of the third dielectric virial coefficient Cɛ(T) for gaseous helium. In contrast to the third pressure virial coefficient, the three-body effects give a substantial contribution to the total value of Cɛ(T)—as large as 50% for temperature range from 30 to 75 K.

Our classical results agree with the classical calculations of Garberoglio et al.38 from 10 K up to 400 K but start to differ at higher temperatures. This difference is attributed to the increasing importance of short-range contributions to the three-body polarizability, treated in Ref. 38 using the superposition approximation. For temperatures below 100 K, we observed significant differences in comparison with the results of quantum calculations of Garberoglio et al.38 Nonetheless, our semiclassical calculations employing the ab initio three-body polarizability agree very well with the fully quantum PIMC calculations38 down to 10 K. Above 273 K, our results are consistent with the available experimental data but have much smaller uncertainties. For lower temperatures, we observed differences with the results of White and Gugan90 and of Hout and Bose,77 which have uncertainties more than an order of magnitude larger than the uncertainties of our semiclassical calculations. These differences are significant in view of the excellent agreement between the semiclassical and fully quantum PIMC calculations at the relevant temperature range.

See the supplementary material for the ab initio values of the three-body polarizability of helium and the Fortran code with the implementation of the analytical fit to the polarizability.

We thanked Giovanni Garberoglio and Krzysztof Szalewicz for their useful comments on this manuscript. Professor Robert Moszyński was thanked for his involvement at the early stage of this work. The authors acknowledged support from QuantumPascal Project No. 18SIB04, which had received funding from the EMPIR Programme co-financed by the Participating States and from the European Union’s Horizon 2020 research and innovation program. The authors also thanked for support from the National Science Center, Poland, Project No. 2017/27/B/ST4/02739. This research was supported in part by PL-Grid infrastructure.

The authors have no conflicts to disclose.

J. Lang: Conceptualization (supporting); Data curation (equal); Formal analysis (lead); Investigation (lead); Software (equal); Writing – original draft (equal). M. Przybytek: Conceptualization (supporting); Data curation (equal); Investigation (supporting); Software (equal); Writing – original draft (equal). M. Lesiuk: Conceptualization (supporting); Investigation (supporting); Writing – original draft (equal). B. Jeziorski: Conceptualization (lead); Funding acquisition (lead); Investigation (supporting); Supervision (lead); Writing – original draft (equal).

The data that support the findings of this study are available within the article and its supplementary material.

In evaluating the three-body contribution to Cɛ one encounters ill-defined integrals of the form

f4(ηr1)fn+1(ηr2)r13r2nP2(cosγ)dr1dr2,
(A1)

where the r1 and r2 integrations extend over the whole 6D space, γ is the angle between the vectors r1 and r2, and η is the fitting parameter specifying the strength of the short-range damping affected by the Tang–Toennies functions (13). The integrals with n = 3 and n = 6 are generated by the asymptotics of the three-body polarizability, see Eq. (12), whereas when n ≥ 9, they appear due to the long-range behavior of the two-body Mayer functions eβU(rij)1, see Eq. (18). The integrals (A1) are ill-defined since they are divergent when the integration is performed first over r1 or r2 and vanish if the angular integration over γ is executed first. To eliminate this ambiguity, we note that the integrals (A1) are derived from the well-defined three-atom, 9D integrals in Eq. (18) when the integration over the center of mass is performed, which leads to the factorization of volume V. This procedure is mathematically justified only if the integral (A1) is convergent. Thus, we have to go back to the original three-particle integral over the finite volume V. This leads to the following 9D integral:

In=limV1Vf4(η|r1a|)|r1a|3fn+1(η|r2a|)|r2a|nP2(cosγ)dadr1dr2,
(A2)

where each of the three 3D integrations is over the interior of the sphere of radius R with the center at the origin of the coordinate system, and γ denotes the angle between the vectors r1a and r2a. In fact, the parameter η > 0 can be different for f4 and fn+1, but this complication is irrelevant for further reasoning. Unlike in the conventional statistical mechanics treatment, no center of mass is introduced here and the integration is over the positions of all three particles, the position of the particle 3 being denoted by a.

Employing the spherical harmonics addition theorem and the spherical symmetry of the integrand in Eq. (A2), one can see that In is equal to the R → ∞ limit of the one-dimensional integral In(R),

In(R)=3R30Ra2I3(R,a)In(R,a)da,
(A3)

over the products of one-particle (3D) integrals In(R, a),

In(R,a)=fn+1(η|ra|)|ra|nP2(cosθa)dr,
(A4)

where the integration is over the interior of the sphere of radius R located at the origin of the coordinate system, the vector a is placed on the z axis, and θa is the angle between the vector ra and the z axis. Note that due to the cylindrical symmetry of the integration over r1 and r2, the spherical harmonics Y2m(θa, ϕa), with m = ±1 and m = ±2 do not contribute to In. Representing P2(cos θa) in terms of Cartesian coordinates, we can write

In(R,a)=12fn+1(η|ra|)|ra|n+2(3z2r24az+2a2)dr,
(A5)

where the integration is again over the interior of the sphere |r| ≤ R. Performing this integration in spherical coordinates r, θ, and ϕ, and making the substitution u = cos θ, we find

In(R,a)=π0Rr21+1fn+1[η(r22aru+a2)1/2]r22aru+a2n/2+1×(3r2u2r24aru+2a2)dudr.
(A6)

Let us first consider the case n = 3. It is then possible to express the integration over u and r in closed form. The result reads

I3(R,a)=2πR33a3f4(ηR+ηa)f4(ηRηa)+πη5R8576a3g7(ηR,a)3a2R2+3g5(ηR,a)+3a4R4+2a2R23g3(ηR,a)a2R213g1(ηR,a),
(A7)

where

g0(w,a)=eww(eaweaw)
(A8)

and gn(w, a) is the nth derivative of g0(w, a) with respect to w multiplied by the factor of (−1)n. The evaluation of the final integral over the variable a in Eq. (A3) is now elementary, although rather tedious, and cannot be performed by hand. However, it can by accomplished by the Mathematica91 software using symbolic integration methods. The result is lengthy, but we are interested only in the leading terms in the large R asymptotic expansion, which are

I3(R)=89π2192η1R25π28η21R2+763π2576η31R3+O1R4.
(A9)

It is clear that the R → ∞ limit of I3(R) equals zero. This means that the regularized integral I3 vanishes and the integral (A1) can be assumed to be zero for n = 3.

Let us now consider the general case of arbitrary n ≥ 4. Using the obvious inequality |ab| ≤ (a2 + b2)/2, valid for any real numbers a and b, we can make the following estimate:

|In(R)|12I3(R)+32R30Ra2In2(R,a)da.
(A10)

When n ≥ 4, the second term on the right-hand side of the inequality (A10) can be rigorously estimated without analytic calculation. For this purpose, we define a “small” ball |ra| ≤ ρ of a radius

ρ=4πRn31/(n3)
(A11)

centered around the point a. Note that this “small” ball does not have to be entirely included in the “large” ball, |r| ≤ R. Since |fn+1(x)| < 1 and P2(cos θa) ≤ 1, one can easily show that

|ra|ρ|r|Rfn+1(η|ra|)|ra|nP2(cosθa)dr1R,
(A12)

where the integration extends now over the intersection of the exterior of the “small” ball and the interior of the “large” ball. For aRρ, the integral In(R, a) can be estimated as

In(R,a)=|ra|<ρ|r|Rfn+1(η|ra|)|ra|nP2(cosθa)dr+|ra|ρ|r|Rfn+1(η|ra|)|ra|nP2(cosθa)dr1R,
(A13)

because the first integral vanishes by spherical symmetry. For a > Rρ, we can use the following estimate:

In(R,a)fn+1(η|ra|)|ra|ndr=Cn,
(A14)

where the integration extends over the whole 3D space so that Cn does not depend on R. Splitting the a integration into aRρ and aRρ and using Eqs. (A13) and (A14), the second term on the right-hand-side of the inequality (A10) can be estimated as follows:

32R30Ra2In2(R,a)da12R1ρR3+12Cn21ρR31=OρR.
(A15)

Since ρR1/(2n−6), we see that, for n ≥ 4, the second term in the inequality in Eq. (A10) goes to zero when R → ∞. In view of Eq. (A9), the same applies to In(R). Therefore, the regularized integral of Eq. (A2) is equal to zero and the ill-defined integrals of Eq. (A1) can be neglected in calculations of Cɛ.

Although in the short-range damping of the asymptotic three-body energy, as given in Eq. (12), we used only the f4(x) and fn+1(x) damping functions, the reasoning similar to that of Eqs. (A10)(A15) shows that the regularized integral (A1) vanishes also when the damping function f4(x) is replaced by any fk(x), with k ≥ 2 or fn+1(x) is replaced by any fk(x), kn − 1.

1.
A.
Borysow
and
L.
Frommhold
,
Adv. Chem. Phys.
75
,
439
(
1989
).
2.
R.
Moszynski
,
T. G. A.
Heijmen
, and
A. v. d.
Avoird
,
Chem. Phys. Lett.
247
,
440
(
1995
).
3.
J. P.
McTague
and
G.
Birnbaum
,
Phys. Rev. Lett.
21
,
661
(
1968
).
4.
R.
Moszynski
,
T. G. A.
Heijmen
,
P. E. S.
Wormer
, and
A.
Van der Avoird
,
J. Chem. Phys.
104
,
6997
(
1996
).
5.
W.
Skomorowski
and
R.
Moszynski
,
Mol. Phys.
111
,
1414
(
2013
).
6.
A. D.
Buckingham
and
J. A.
Pople
,
Trans. Faraday Soc.
51
,
1179
(
1955
).
7.
C.
Gaiser
,
B.
Fellmuth
,
N.
Haft
,
A.
Kuhn
,
B.
Thiele-Krivoi
,
T.
Zandt
,
J.
Fischer
,
O.
Jusko
, and
W.
Sabuga
,
Metrologia
54
,
280
(
2017
).
8.
C.
Gaiser
,
T.
Zandt
, and
B.
Fellmuth
,
Metrologia
52
,
S217
(
2015
).
9.
C.
Gaiser
,
B.
Fellmuth
, and
N.
Haft
,
Metrologia
54
,
141
(
2017
).
10.
C.
Gaiser
,
B.
Fellmuth
, and
N.
Haft
,
Metrologia
57
,
055003
(
2020
).
11.
D. M.
Ripa
,
D.
Imbraguglio
,
C.
Gaiser
,
P. P. M.
Steur
,
D.
Giraudi
,
M.
Fogliati
,
M.
Bertinetti
,
G.
Lopardo
,
R.
Dematteis
, and
R. M.
Gavioso
,
Metrologia
58
,
025008
(
2021
).
12.
J. W.
Schmidt
,
R. M.
Gavioso
,
E. F.
May
, and
M. R.
Moldover
,
Phys. Rev. Lett.
98
,
254504
(
2007
).
13.
C.
Gaiser
,
B.
Fellmuth
, and
W.
Sabuga
,
Nat. Phys.
16
,
177
(
2020
).
14.
B.
Gao
,
H.
Zhang
,
D.
Han
,
C.
Pan
,
H.
Chen
,
Y.
Song
,
W.
Liu
,
J.
Hu
,
X.
Kong
,
F.
Sparasci
 et al,
Metrologia
57
,
065006
(
2020
).
15.
P. M. C.
Rourke
,
C.
Gaiser
,
B.
Gao
,
D. M.
Ripa
,
M. R.
Moldover
,
L.
Pitre
, and
R. J.
Underwood
,
Metrologia
56
,
032001
(
2019
).
16.
P. M. C.
Rourke
,
J. Phys. Chem. Ref. Data
50
,
033104
(
2021
).
17.
K.
Piszczatowski
,
M.
Puchalski
,
J.
Komasa
,
B.
Jeziorski
, and
K.
Szalewicz
,
Phys. Rev. Lett.
114
,
173004
(
2015
).
18.
M.
Puchalski
,
K.
Szalewicz
,
M.
Lesiuk
, and
B.
Jeziorski
,
Phys. Rev. A
101
,
022505
(
2020
).
19.
P.
Czachorowski
,
M.
Przybytek
,
M.
Lesiuk
,
M.
Puchalski
, and
B.
Jeziorski
,
Phys. Rev. A
102
,
042810
(
2020
).
20.
W.
Cencek
,
K.
Patkowski
, and
K.
Szalewicz
,
J. Chem. Phys.
131
,
064105
(
2009
).
21.
W.
Cencek
,
J.
Komasa
, and
K.
Szalewicz
,
J. Chem. Phys.
135
,
014301
(
2011
).
22.
E.
Vogel
,
B.
Jäger
,
R.
Hellmann
, and
E.
Bich
,
Mol. Phys.
108
,
3335
(
2010
).
23.
K.
Patkowski
and
K.
Szalewicz
,
J. Chem. Phys.
133
,
094304
(
2010
).
24.
W.
Cencek
,
M.
Jeziorska
,
O.
Akin-Ojo
, and
K.
Szalewicz
,
J. Phys. Chem. A
111
,
11311
(
2007
).
25.
W.
Cencek
,
G.
Garberoglio
,
A. H.
Harvey
,
M. O.
McLinden
, and
K.
Szalewicz
,
J. Phys. Chem. A
117
,
7542
(
2013
).
26.
U.
Bafile
,
L.
Ulivi
,
M.
Zoppi
, and
F.
Barocchi
,
Chem. Phys. Lett.
138
,
559
(
1987
).
27.
U.
Bafile
,
L.
Ulivi
,
M.
Zoppi
,
M.
Moraldi
, and
L.
Frommhold
,
Phys. Rev. A
44
,
4450
(
1991
).
28.
S.
Pestelli
,
U.
Bafile
,
L.
Ulivi
, and
M.
Zoppi
,
Phys. Rev. A
49
,
4602
(
1994
).
29.
F.
Barocchi
,
M.
Celli
, and
M.
Zoppi
,
Phys. Rev. A
38
,
3984
(
1988
).
30.
J.
Van der Elsken
and
R. A.
Huijts
,
J. Chem. Phys.
88
,
3007
(
1988
).
31.
A. D.
Buckingham
and
I. D.
Hands
,
Chem. Phys. Lett.
185
,
544
(
1991
).
32.
J. J.
Pérez
,
J. H. R.
Clarke
, and
A.
Hinchliffe
,
Chem. Phys. Lett.
104
,
583
(
1984
).
33.
M. H.
Champagne
,
X.
Li
, and
K. L. C.
Hunt
,
J. Chem. Phys.
112
,
1893
(
2000
).
34.
X.
Li
and
K. L. C.
Hunt
,
J. Chem. Phys.
105
,
4076
(
1996
).
35.
X.
Li
and
K. L. C.
Hunt
,
J. Chem. Phys.
107
,
4133
(
1997
).
36.
C.
Aleman
,
J. J.
Perez
, and
A.
Hinchliffe
,
Int. J. Mass Spectrom. Ion Processes
122
,
331
(
1992
).
37.
D. F.
Heller
and
W. M.
Gelbart
,
Chem. Phys. Lett.
27
,
359
(
1974
).
38.
G.
Garberoglio
,
A. H.
Harvey
, and
B.
Jeziorski
,
J. Chem. Phys.
155
,
234103
(
2021
).
39.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill
,
New York
,
1965
).
40.
G.
Garberoglio
and
A. H.
Harvey
,
J. Res. Natl. Inst. Stand. Technol.
114
,
249
(
2009
).
41.
G.
Garberoglio
and
A. H.
Harvey
,
J. Chem. Phys.
134
,
134106
(
2011
).
42.
K. R. S.
Shaul
,
A. J.
Schultz
, and
D. A.
Kofke
,
J. Chem. Phys.
137
,
184101
(
2012
).
43.
B. J.
Alder
,
J. C.
Beers
 II
,
H. L.
Strauss
, and
J. J.
Weis
,
Proc. Natl. Acad. Sci. U. S. A.
77
,
3098
(
1980
).
44.
K.
Hald
,
F.
Pawłowski
,
P.
Jørgensen
, and
C.
Hättig
,
J. Chem. Phys.
118
,
1292
(
2003
).
45.
G.
Garberoglio
,
A. H.
Harvey
,
J.
Lang
, and
B.
Jeziorski
, “
Path-integral calculation of the third dielectric virial coefficient of helium based on an ab initio three-body dipole and polarizability surfaces
,” (unpublished) (2022).
46.
C.
Gaiser
,
B.
Fellmuth
, and
W.
Sabuga
,
Ann. Phys.
534
,
2200336
(
2022
).
47.
C. G.
Gray
and
K. E.
Gubbins
,
Theory of Molecular Fluids: Volume 1: Fundamentals
(
Oxford University Press,
,
1984
), Vol. 1.
48.
K. L. C.
Hunt
,
Y. Q.
Liang
, and
S.
Sethuraman
,
J. Chem. Phys.
89
,
7126
(
1988
).
49.
P. R.
Certain
and
P. J.
Fortune
,
J. Chem. Phys.
55
,
5818
(
1971
).
50.
K.
Aidas
,
C.
Angeli
,
K. L.
Bak
,
V.
Bakken
,
R.
Bast
,
L.
Boman
,
O.
Christiansen
,
R.
Cimiraglia
,
S.
Coriani
,
P.
Dahle
 et al,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
269
(
2014
).
51.
Dalton, a molecular electronic structure program, release 2018
,”
2018
, see http://daltonprogram.org.
52.
S. F.
Boys
and
F.
Bernardi
,
Mol. Phys.
19
,
553
(
1970
).
53.
B.
Skwara
,
W.
Bartkowiak
, and
D. L.
Da Silva
,
Theor. Chim. Acta
122
,
127
(
2009
).
54.
W.
Cencek
,
M.
Przybytek
,
J.
Komasa
,
J. B.
Mehl
,
B.
Jeziorski
, and
K.
Szalewicz
,
J. Chem. Phys.
136
,
224303
(
2012
).
55.
T. H.
Dunning
and
J.
Chem
,
Phys.
90
,
1007
(
1989
).
56.
D. E.
Woon
and
T. H.
Dunning
,
J. Chem. Phys.
100
,
2975
(
1994
).
57.
A.
Halkier
,
T.
Helgaker
,
P.
Jørgensen
,
W.
Klopper
,
H.
Koch
,
J.
Olsen
, and
A. K.
Wilson
,
Chem. Phys. Lett.
286
,
243
(
1998
).
58.
T.
Helgaker
,
W.
Klopper
, and
D. P.
Tew
,
Mol. Phys.
106
,
2107
(
2008
).
59.
A.
Halkier
,
T.
Helgaker
,
W.
Klopper
, and
J.
Olsen
,
Chem. Phys. Lett.
319
,
287
(
2000
).
60.
M.
Przybytek
,
W.
Cencek
,
B.
Jeziorski
, and
K.
Szalewicz
,
Phys. Rev. Lett.
119
,
123401
(
2017
).
61.
M.
Przybytek
, General FCI program hector,
2014
.
62.
K. T.
Tang
and
J. P.
Toennies
,
J. Chem. Phys.
80
,
3726
(
1984
).
63.
D. M.
Bishop
and
J.
Pipin
,
Chem. Phys. Lett.
236
,
15
(
1995
).
64.
M. P.
Metz
and
K.
Szalewicz
,
J. Chem. Phys.
152
,
134111
(
2020
).
65.
M. J. D.
Powell
,
Numerical Methods for Nonlinear Algebraic Equations
, edited by
P.
Rabinowitz
(
Gordon and Breach Science
,
London
,
1970
), p.
87
144
.
66.
M.
Galassi
,
J.
Davies
,
J.
Theiler
,
B.
Gough
,
G.
Jungman
,
P.
Alken
,
M.
Booth
,
F.
Rossi
, and
R.
Ulerich
,
GNU Scientific Library
(
Network Theory Limited
,
2002
).
67.
T. L.
Hill
,
J. Chem. Phys.
28
,
61
(
1958
).
68.
C. G.
Gray
,
K. E.
Gubbins
, and
C. G.
Joslin
,
Theory of Molecular Fluids: Volume 2: Applications
(
Oxford University Press,
,
2011
), Vol. 2.
69.
P. G.
Kusalik
,
F.
Liden
, and
I. M.
Svishchev
,
J. Chem. Phys.
103
,
10169
(
1995
).
70.
E. M.
Mas
,
V. F.
Lotrich
, and
K.
Szalewicz
,
J. Chem. Phys.
110
,
6694
(
1999
).
71.
G.
Garberoglio
and
A. H.
Harvey
,
J. Res. Natl. Inst. Stand. Technol.
125
,
022
(
2020
).
72.
A.
Kronrod
,
Nodes and Weights of Quadrature Formulas: Sixteen-Place Tables
(
Consultants Bureau
,
New York
,
1965
).
73.
J.
Lang
,
G.
Garberoglio
,
M.
Przybytek
,
M.
Jeziorska
, and
B.
Jeziorski
, “
Three-body potential and third virial coefficients for helium including relativistic and nuclear motion effects
,” (unpublished) (2022).
74.
G.
Garberoglio
and
A. H.
Harvey
,
J. Chem. Phys.
154
,
104107
(
2021
).
75.
C.
Gaiser
and
B.
Fellmuth
,
J. Chem. Phys.
150
,
134303
(
2019
).
76.
M.
Lallemand
and
D.
Vidal
,
J. Chem. Phys.
66
,
4776
(
1977
).
77.
J.
Huot
and
T. K.
Bose
,
J. Chem. Phys.
95
,
2683
(
1991
).
78.
S.
Kirouac
and
T. K.
Bose
,
J. Chem. Phys.
64
,
1580
(
1976
).
79.
G.
Garberoglio
,
M. R.
Moldover
, and
A. H.
Harvey
,
J. Res. Natl. Inst. Stand. Technol.
116
,
729
(
2011
).
80.
E.
Bich
,
R.
Hellmann
, and
E.
Vogel
,
Mol. Phys.
105
,
3035
(
2007
).
81.
A.
Rizzo
,
C.
Hättig
,
B.
Fernández
, and
H.
Koch
,
J. Chem. Phys.
117
,
2609
(
2002
).
82.
K. R. S.
Shaul
,
A. J.
Schultz
,
D. A.
Kofke
, and
M. R.
Moldover
,
Chem. Phys. Lett.
531
,
11
(
2012
).
83.
J.
Ram
and
Y.
Singh
,
Mol. Phys.
26
,
539
(
1973
).
84.
B.
Song
and
Q.-Y.
Luo
,
Metrologia
57
,
025007
(
2020
).
85.
J. O.
Hirschfelder
,
C. F.
Curtiss
, and
R. B.
Bird
,
Molecular Theory of Gases and Liquids
(
Wiley
,
New York
,
1954
).
86.
H. E.
DeWitt
,
J. Math. Phys.
3
,
1003
(
1962
).
87.
T.
Kihara
,
Rev. Mod. Phys.
25
,
831
(
1953
).
88.
B.
Guillot
and
Y.
Guissani
,
J. Chem. Phys.
108
,
10162
(
1998
).
89.
R. P.
Feynman
,
Statistical Mechanics; A Set of Lectures
(
W.A. Benjamin
,
Reading, MA
,
1972
).
90.
M. P.
White
and
D.
Gugan
,
Metrologia
29
,
37
(
1992
).
91.
Wolfram Research, Inc.
, “
Mathematica, version 7.0
,”
Champaign, IL
(
2008
).

Supplementary Material