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.

## I. INTRODUCTION

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 constant^{7} and recently in the development of new temperature^{8–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 theory^{34,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 $r12\u2212lr23\u2212mr13\u2212n$, with *l* + *m* + *n* ≤ 9, where *r*_{ij} 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 Gelbart^{37} 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}

## II. THREE-BODY POLARIZABILITY OF HELIUM AND ITS ASYMPTOTICS AT LARGE INTERATOMIC DISTANCES

The components *α*_{μν} of a polarizability tensor ** α** are defined as

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 **r**_{i}, *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

where *α*_{0} is the atomic polarizability, **e** is the unit tensor, *α*_{2}(**r**_{i}, **r**_{j}) is the two-body interaction-induced polarizability tensor of the pair of atoms at **r**_{i} and **r**_{j},

and *α*_{3}(**r**_{1}, **r**_{2}, **r**_{3}) is the three-body interaction-induced polarizability tensor defined essentially in Eq. (2). By ** α**(

**r**

_{i},

**r**

_{j}), we denote the tensor of the total polarizability of two atoms located at

**r**

_{i}and

**r**

_{j}. The tensor

*α*_{3}(

**r**

_{1},

**r**

_{2},

**r**

_{3}), representing the pair-wise non-additive part of the total interaction-induced polarizability tensor

**(**

*α***r**

_{1},

**r**

_{2},

**r**

_{3}) − 3

*α*

_{0}

**e**, can be calculated from Eq. (2) using precomputed two-body polarizability tensors

*α*_{2}(

**r**

_{i},

**r**

_{j}). In this work, it has been calculated directly from the formula

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,

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 *r*_{12}, *r*_{23}, and *r*_{13}. 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) model^{31,33} and takes the form

where *θ*_{3} is the angle between vectors **r**_{1} − **r**_{3} and **r**_{2} − **r**_{3}, *P*_{n}(*x*) is the Legendre polynomial of order *n*, and $S123$ denotes the operator that symmetrizes the expression it precedes with respect to the variables **r**_{1}, **r**_{2}, and **r**_{3}. The classical expression for the third dielectric virial coefficient requires the integration of Eq. (6) over *r*_{13}, *r*_{23}, and *θ*_{3}. It is clear that this integration is ill-defined as it is divergent both at small and at large values of *r*_{13} or *r*_{23}. 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:

where *r*, *R*, and *θ* are the Jacobi coordinates: *r* and *R* are the lengths of the vectors **r** = **r**_{2} − **r**_{1} and **R** = **r**_{3} − (**r**_{1} + **r**_{2})/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 **r**_{1} and **r**_{2}, the symmetrizer $S123$ can be replaced by the operator $2+2P13+2P23$, where the operator $Pij$ permutes the vectors **r**_{i} and **r**_{j}.

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\alpha 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.

## III. COMPUTATIONAL DETAILS AND RESULTS

### A. Polarizability calculations

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 *r*_{12}, *r*_{23}, and *r*_{13}. 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 method^{44} 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 d*X*Z, 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 d*X*Z basis sets and the standard Dunning’s d-aug-cc-pV*X*Z basis sets for helium^{55,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 d*X*Z 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.

r_{12}
. | r_{23}
. | r_{13}
. | X = T
. | X = Q
. | X = 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%) |

r_{12}
. | r_{23}
. | r_{13}
. | X = T
. | X = Q
. | X = 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/*X*^{n}. 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:

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,

r_{12}
. | r_{23}
. | r_{13}
. | $\alpha 3CC3$ . | $\alpha 3FCI$ . | $\delta $ . |
---|---|---|---|---|---|

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% |

r_{12}
. | r_{23}
. | r_{13}
. | $\alpha 3CC3$ . | $\alpha 3FCI$ . | $\delta $ . |
---|---|---|---|---|---|

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.

### B. Analytic fit of the three-body polarizability surface

The three-body polarizability surface *α*_{3}(*r*_{12}, *r*_{23}, *r*_{13}) was fitted with an analytic function *α*_{fit}(*r*_{12}, *r*_{23}, *r*_{13}) expressed as the sum of three components, one of a short-range and two of a long-range character,

The short-range function *α*_{sr}(*r*_{12}, *r*_{23}, *r*_{13}), representing the exponentially decaying contributions to *α*_{3}(*r*_{12}, *r*_{23}, *r*_{13}), has the form

where *A*_{ijk,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 ≤ *i* ≤ *j* ≤ *k* ≤ 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,

where *f*_{n} are the Tang–Toennies damping functions^{62}

The *Z*_{n} coefficients can be calculated with the knowledge of several atomic properties: the static dipole and quadrupole polarizabilities (*Z*_{1}, …, *Z*_{4}), the dynamic dipole polarizabilities (*Z*_{5} and *Z*_{6}), and the second dipole hyperpolarizabilities (*Z*_{5} and *Z*_{6}) (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,

where $S\u0303123=1+P13+P23$ generates the sum over three possible ways of defining the Jacobi coordinates *r*, *R*, and *θ* for a triangle with sides *r*_{12}, *r*_{23}, and *r*_{13}. 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 *C*_{0} and *C*_{1} different from zero improves the quality of the fit. The parameters *C*_{0}, *C*_{1}, *C*_{2}, *η*_{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(*r*_{12}, *r*_{23}, *r*_{13}) > 1.2 a.u. and max(*r*_{12}, *r*_{23}, *r*_{13}) < 14 a.u., was created through the Sobol sequence. In each step of the iterative procedure, ∼10^{6} 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 ∼10^{6} independent fitting functions. The optimization was performed using the trust-region dogleg procedure^{65} 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 ∼10^{4} 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.

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 *σ*(*r*_{12}, *r*_{23}, *r*_{13}) 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}(*r*_{12}, *r*_{23}, *r*_{13}) ± *σ*(*r*_{12}, *r*_{23}, *r*_{13}). The function *σ*(*r*_{12}, *r*_{23}, *r*_{13}) 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 ≤ *i* ≤ *j* ≤ *k* ≤ 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(*r*_{12}, *r*_{23}, *r*_{13}) > 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(*r*_{12}, *r*_{23}, *r*_{13}) < 2 a.u., which are of little importance in calculations of the physical properties. The final set of points used in the construction of *σ*(*r*_{12}, *r*_{23}, *r*_{13}) comprised almost 65% of the total number of configurations. The average ratio of the values of the fitting function *σ*(*r*_{12}, *r*_{23}, *r*_{13}) to the estimated uncertainties for the whole set of points is about 1.3. The values of all linear and nonlinear parameters of *α*_{3}(*r*_{12}, *r*_{23}, *r*_{13}) and *σ*(*r*_{12}, *r*_{23}, *r*_{13}), along with a Fortran 2003 implementation of the fitting functions, are included in the supplementary material.

### C. Classical calculations of the third dielectric virial coefficient

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}

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,

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

and

In the above formulas, *U*(**r**_{1}, **r**_{2}, **r**_{3}) = *U*_{3}(**r**_{1}, **r**_{2}, **r**_{3}) + *U*_{2}(**r**_{1}, **r**_{2}) + *U*_{2}(**r**_{2}, **r**_{3}) + *U*_{2}(**r**_{1}, **r**_{3}) is the sum of the two-body, *U*_{2}, and three-body, *U*_{3}, interaction potentials, and *β* = 1/*k*_{B}*T*. Note that the integral $C\epsilon 2body(T)$ involves not only the two-body quantities (polarizability and interaction potential) but also the three-body interaction potential *U*_{3}. 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 $V\u22121dr1dr2dr3\u21928\pi r132r232dr13dr23dcos\theta 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 *r*_{13}, *r*_{23}, and *θ*_{3}. However, the resulting integral is ill-defined as the integrals over *r*_{13} and *r*_{23} are divergent at infinity, even if the singularities at *r*_{13} = 0 and *r*_{23} = 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 *r*_{12}. 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\epsilon 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 negligible^{19} 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 Harvey^{74} for calculations of the pressure viral coefficients. The uncertainties due to the polarizabilities were estimated employing the following formula:

where $\alpha 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\epsilon (\alpha nacc)\u2212C\epsilon (\alpha n)|\u2264\sigma (\alpha 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 *U*_{3}. In this case, the equality in Eq. (19) is valid only approximately with an error quadratic in *σ*(*U*_{n}). The total uncertainty was computed as the square root of the sum of squares of the uncertainties of *α*_{2}, *α*_{3}, and *U*_{3}.

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\epsilon 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\epsilon 3body(T)$ cross zero and the two-body terms dominate. However, this trend suggests an increasing importance of $C\epsilon 3body(T)$ for temperatures higher than 2500 K.

From the inspection of the contributions to $C\epsilon 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\epsilon 3body(T)$ dominates for the entire temperature range. The largest contribution to $C\epsilon 3body(T)$ then comes from the uncertainty of the three-body polarizability.

The calculated values of $C\epsilon 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 cm^{9} 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 cm^{9} mol^{−3},^{38} which is close to the value *C*_{ɛ}(300 K) = −0.535 cm^{9} 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.

T
. | $C\epsilon 3body$ . | $4BB\epsilon +C\epsilon 2body$ . | C_{ɛ}
. |
---|---|---|---|

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) |

T
. | $C\epsilon 3body$ . | $4BB\epsilon +C\epsilon 2body$ . | C_{ɛ}
. |
---|---|---|---|

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) |

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 virials^{54,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 *k*_{B}*T*. 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.

## IV. SEMICLASSICAL QUANTUM CORRECTION

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}

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

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\epsilon 2body$) in Eq. (16). The resulting semiclassical values of the third dielectric virial coefficients are presented in Fig. 7.

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.

## V. CONCLUSIONS

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 calculations^{38} 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 Gugan^{90} 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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX: REGULARIZATION OF INTEGRALS CONTRIBUTING TO *C*_{ɛ}(*T*)

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

where the **r**_{1} and **r**_{2} integrations extend over the whole 6D space, *γ* is the angle between the vectors **r**_{1} and **r**_{2}, 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\u2212\beta U(rij)\u22121$, see Eq. (18). The integrals (A1) are ill-defined since they are divergent when the integration is performed first over *r*_{1} or *r*_{2} 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:

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 **r**_{1} − **a** and **r**_{2} − **a**. In fact, the parameter *η* > 0 can be different for *f*_{4} and *f*_{n+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)$,

over the products of one-particle (3D) integrals *I*_{n}(*R*, *a*),

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 **r** − **a** and the *z* axis. Note that due to the cylindrical symmetry of the integration over **r**_{1} and **r**_{2}, the spherical harmonics *Y*_{2m}(*θ*_{a}, *ϕ*_{a}), with *m* = ±1 and *m* = ±2 do not contribute to $In$. Representing *P*_{2}(cos *θ*_{a}) in terms of Cartesian coordinates, we can write

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

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

where

and *g*_{n}(*w*, *a*) is the *n*th derivative of *g*_{0}(*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 Mathematica^{91} 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

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*| ≤ (*a*^{2} + *b*^{2})/2, valid for any real numbers *a* and *b*, we can make the following estimate:

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 |**r** − **a**| ≤ *ρ* of a radius

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

where the integration extends now over the intersection of the exterior of the “small” ball and the interior of the “large” ball. For *a* ≤ *R* − *ρ*, the integral *I*_{n}(*R*, *a*) can be estimated as

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

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

Since *ρ* ∼ *R*^{1/(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 *f*_{4}(*x*) and *f*_{n+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 *f*_{4}(*x*) is replaced by any *f*_{k}(*x*), with *k* ≥ 2 or *f*_{n+1}(*x*) is replaced by any *f*_{k}(*x*), *k* ≥ *n* − 1.