Recent advances regarding the interplay between *ab initio* calculations and metrology are reviewed, with particular emphasis on gas-based techniques used for temperature and pressure measurements. Since roughly 2010, several thermophysical quantities – in particular, virial and transport coefficients – can be computed from first principles without uncontrolled approximations and with rigorously propagated uncertainties. In the case of helium, computational results have accuracies that exceed the best experimental data by at least one order of magnitude and are suitable to be used in primary metrology. The availability of *ab initio* virial and transport coefficients contributed to the recent SI definition of temperature by facilitating measurements of the Boltzmann constant with unprecedented accuracy. Presently, they enable the development of primary standards of thermodynamic temperature in the range 2.5–552 K and pressure up to 7 MPa using acoustic gas thermometry, dielectric constant gas thermometry, and refractive index gas thermometry. These approaches will be reviewed, highlighting the effect of first-principles data on their accuracy. The recent advances in electronic structure calculations that enabled highly accurate solutions for the many-body interaction potentials and polarizabilities of atoms – particularly helium – will be described, together with the subsequent computational methods, most often based on quantum statistical mechanics and its path-integral formulation, that provide thermophysical properties and their uncertainties. Similar approaches for molecular systems, and their applications, are briefly discussed. Current limitations and expected future lines of research are assessed.

## 1. Introduction

On May 20, 2019, the base SI units of mass (kilogram), electric current (ampere), temperature (kelvin) and amount of substance (mole) were redefined by assigning fixed values to fundamental constants of nature: the Planck constant, the electron charge, the Boltzmann constant, and the Avogadro constant, respectively.^{1–3} By decoupling the base units from specific material artifacts, this new redefinition is expected to lead to improved scientific instruments, reducing the degradation in accuracy when measuring quantities at larger or smaller magnitudes than a predefined unit standard. Additionally, the most accurate experimental technique available at each scale can be used to implement a primary standard, resulting in easier calibrations, increased accuracies of measuring devices, and further technological advancements.

At many conditions, gas-based techniques provide unparalleled performance for primary measurements of temperature and pressure. These involve acoustic, dielectric, or refractivity measurements, because frequency and electromagnetic measurements can be made with very high accuracy. A model, typically expressed as the ideal-gas behavior with a series of corrections in powers of density, is used to relate the measured quantity to the temperature or pressure; in the case of dielectric or refractivity measurements, one set of corrections relates the measured quantity to the gas density and the familiar virial expansion is used to relate the density to the pressure and temperature.

These gas-based methods have been greatly facilitated in recent years by the ability to perform *ab initio* calculations of the thermophysical properties (such as the polarizability and the density, dielectric, and refractivity virial coefficients) of the working gases with no uncontrolled approximations and rigorously defined uncertainties. These calculated properties often have much smaller uncertainties than the best experimental determinations, especially when the gas considered is helium. These techniques have been successfully applied for pressures up to 7 MPa and for thermodynamic temperatures in the range (2.5–552) K (with extension to 1000 K or more progressing^{4}).

These achievements have been facilitated by the increase in supercomputing power and advances in numerical techniques for electronic structure calculations. For example, state-of-the-art calculations for up to three He atoms even include relativistic and quantum electrodynamics effects. In particular, these numerical investigations produce pair and three-body potentials, as well as single-atom, pair, and three-body polarizabilities, with unprecedented accuracy.

Building on these results, the exact quantum statistical mechanics formulation enabled rigorous calculations of the coefficients appearing in the density (virial) expansion of the equation of state, the speed of sound, the dielectric constant, and the refractive index. The path-integral Monte Carlo (PIMC) method has been shown to provide sufficient accuracy for these quantities. As a consequence, it has been possible to devise a fully first-principles chain of calculations with rigorous uncertainty propagation to compute virial coefficients of helium gas.

As a result of these endeavors, since about 2010 thermophysical properties of gaseous helium have been known from theory with an accuracy that in most cases surpasses that of the most precise experimental determinations. Currently, the uncertainties of the *ab initio* second and third virial coefficients of helium are at least one order of magnitude smaller than the experimental ones. The situation is similar for the density dependence of the speed of sound, the dielectric constant, and the refractive index, where it leads to improved accuracy in Acoustic Gas Thermometry (AGT), Dielectric Constant Gas Thermometry (DCGT), and Refractive Index Gas Thermometry (RIGT), respectively.

Section 2 describes these gas-based experimental techniques for temperature and pressure measurement, including their operating principles, temperature and pressure ranges, recent technological improvements, and the sources of uncertainty. We highlight the ways in which theoretical knowledge, in the form of *ab initio* polarizabilities and virial coefficients, has improved these measurements by reducing significant components of the uncertainty.

First-principles calculations of virial coefficients involve two steps: the *ab initio* electronic structure calculation of interatomic potentials and/or polarizabilities, followed by the solution of the exact quantum statistical equations describing virial coefficients.

We therefore present in Sec. 3 a critical review of the state of the art of non-relativistic, relativistic, and quantum electrodynamic electronic structure calculations, with particular emphasis on the determination of uncertainties. Our primary focus will be on helium – which is currently the only substance for which computations can be performed that consistently exceed the accuracy of the best experiments – but other noble gases will be briefly covered due to their importance in metrology. For the sake of completeness, we will recall the hierarchy of physical theories involved in quantum chemical calculations, with particular emphasis on the Full Configuration Interaction (FCI) approach, which is exact within a given orbital basis set and is currently feasible for systems with up to ten electrons. Relativistic and quantum electrodynamic effects (expressed as expansions in powers of the fine-structure constant) have been crucial for achieving the extremely low uncertainty of the latest helium calculations, and are also progressively important in describing larger atoms (notably, neon and argon). Additionally, the evaluation of electronic polarizabilities and magnetic susceptibilities will be discussed. All of these theoretical advances will be exemplified for the case of helium, where we will present the current state of the art regarding interaction potentials and many-body polarizabilities.

Knowledge of interaction potentials and polarizabilities enables calculation of the coefficients appearing in the virial expansion of the equation of state, the speed of sound, the dielectric constant, and the refractive index, which are crucial ingredients in the uncertainty budgets of AGT, DCGT, and RIGT. In the past 15 years, the path-integral approach to quantum statistical mechanics has been successfully applied in calculating virial coefficients without uncontrolled approximations. The main features of this method are reviewed in Sec. 4, with particular attention to the question of uncertainty propagation from the potentials and the polarizabilities. In the case of pair properties, an alternative method based on the solution of the Schrödinger equation is available and provides mutual validation of the path-integral results, as well as enabling the calculation of transport properties. Most of this review is focused on thermodynamic properties, but *ab initio* calculations also provide viscosity and thermal conductivity. We briefly review how this leads to improvements in flow-rate measurements.

Although most efforts have been devoted to noble gases, highly accurate theoretical calculations are also available for molecular systems and have the potential to enable a similar paradigm shift in some metrological applications. We describe in Sec. 5 the present situation in the first-principles calculation of molecular properties, and point out a few areas where computational contributions are expected to have an increasing impact in the near future, namely humidity metrology, measurements of very low pressures, and atmospheric science. We end our review in Sec. 6, where future perspectives and an overview of the status of highly accurate *ab initio* property calculations will be presented.

## 2. Primary Metrology and Thermophysical Properties

### 2.1. Paradigm reversal in temperature metrology

Traditionally, accurate measurements of temperature-dependent thermophysical properties of gases [such as: second density virial coefficient *B*(*T*), viscosity *η*(*T*), thermal conductivity *λ*(*T*)] have been used to determine parameters in evermore-refined models for interatomic and intermolecular potentials. This tradition/paradigm can be traced back to the 18th century when “… Bernoulli had proposed that in Boyle’s law the specific volume *v* be replaced by (*v* − *b*), where *b* was thought to be the volume of the molecules.”^{5} During the past 25 years, the accuracy of the calculated thermophysical properties of the noble gases (particularly helium) has increased dramatically. An example is shown in Fig. 1, which shows how the accuracy of the second virial coefficient *B*(*T*) of ^{4}He improved with time. The data plotted are for temperatures near *T*_{Ne}. (*T*_{Ne} ≡ 24.5561 K is the defined temperature of the triple point of neon on the international temperature scale, ITS-90.^{6}) Since the year 2012, the uncertainty of *B*(*T*_{Ne}), as calculated *ab initio*, has been smaller than the uncertainty of the best measurements of *B*(*T*_{Ne}).

The paradigm reversal (replacing measured thermophysical properties of helium with calculated thermophysical properties) applies to zero-density values of the viscosity *η*(*T*), thermal conductivity *λ*(*T*), and ^{3}He–^{4}He mutual diffusion coefficient as well as to the density and acoustic virial coefficients, relative dielectric permittivity (dielectric constant) *ɛ*_{r}(*p*, *T*), relative magnetic permittivity *μ*_{r}(*p*, *T*), and refractive index $n(p,T)=\epsilon r\mu r$. For many of these properties, the values calculated for helium are standards that are used to calibrate apparatus that measures the same properties for other gases.

The paradigm reversals for *ɛ*_{r}(*p*, *T*) and *n*(*p*, *T*) have been combined with technical advances in the measurement of *ɛ*_{r}(*p*, *T*) and *n*(*p*, *T*) to develop novel pressure standards. One standard operating at optical frequencies and low pressures (100 Pa $\u2264p\u2264$ 100 kPa) is more accurate than manometers based on liquid columns (see Sec. 2.3.1 and Ref. 18). Other standards operating at audio and microwave frequencies and higher pressures (100 kPa $\u2264p\u2264$ 7 MPa) have enabled exacting tests of mechanical pressure generators based on the dimensions of a rotating piston in a cylinder (see Sec. 2.3.2 and Refs. 19 and 20). At still higher pressures (up to 40 MPa), the values of helium’s density calculated from the virial equation of state (VEOS) have been used to calibrate magnetic suspension densimeters.^{21} A more accurate high-pressure scale may result. In Sec. 2.5, we will comment on *ab initio* calculations of transport properties and their contribution to improved flow metrology.

During the past 25 years, the accurate calculations of the thermophysical properties of the noble gases have strongly interacted with gas-based measurements of the thermodynamic temperature *T*. To put this in context, we compare in Fig. 2 the evolution of “consensus” temperature metrology with “thermodynamic” temperature metrology.^{22}

In Fig. 2, the squares represent estimates of the relative uncertainties *u*_{r}(*T*_{scale}) of the consensus temperature scales disseminated by National Metrology Institutes (NMIs). We plot the values of *u*_{r}(*T*_{scale}) near the boiling point of water at intervals of roughly 20 years. Most of the points are at years when the NMIs agreed to disseminate a new consensus scale that was either a better approximation of thermodynamic temperatures and/or an extension of the consensus scale to higher and lower temperatures. The most recent scale is the “International Temperature Scale of 1990” (ITS-90), and temperatures measured using ITS-90 are denoted *T*_{90}.^{6} The data underlying ITS-90 are constant-volume gas thermometry (CVGT) and spectral radiation thermometry linked to CVGT.^{23} The pre-1990 CVGT was based on the ideal-gas equation of state, as corrected by virial coefficients either taken from the experimental literature or measured during the CVGT. Post-1990 thermometry, together with *ab initio* calculations of virial coefficients, revealed that the authors of ITS-90 were unaware that errors in ITS-90 exceeded their expanded (*k* = 2) uncertainty by roughly a factor of two. (See Fig. 3 and the discussion at the end of Sec. 2.2.1).

In Fig. 2, the circles represent the relative uncertainty of determinations of the Boltzmann constant *u*_{r}(*k*_{B}). To determine *k*_{B}, one measures the mean energy *k*_{B}*T* per degree-of-freedom of a system in thermal equilibrium at the thermodynamic temperature *T*. During the interval 1960–2019, the thermodynamic temperature of the triple point of water was defined as *T*_{TPW} ≡ 273.16 K, exactly. Thus, measurements of *k*_{B}*T* that were conducted near *T*_{TPW} had a negligible uncertainty from *T* and *u*_{r}(*k*_{B}*T*) was an excellent proxy for *u*_{r}(*T*), the uncertainty of measurements of *T* under the most favorable conditions.

As displayed in Fig. 2, *u*_{r}(*T*_{scale}) decreased from $\u223c10$ to $\u223c2$ ppm (1 ppm ≡ 1 part in 10^{6}) during the 20th century. Also during the 20th century, the relative uncertainty *u*_{r}(*k*_{B}) decreased from $\u223c20000$ to $\u223c2$ ppm. Thus, *u*_{r}(*k*_{B}) ≫ *u*_{r}(*T*_{scale}) for most of the 20th century, even though *k*_{B} was a “fundamental” constant and, therefore, a worthy challenge for metrology.

Between 1973 and 2017, AGT measurements decreased the uncertainty of *u*_{r}(*k*_{B}) 100-fold from $\u223c40$ to $\u223c0.4$ ppm.^{24,25} By 2017, DCGT achieved the uncertainty *u*_{r}(*k*_{B}) = 1.9 ppm and Johnson noise thermometry achieved *u*_{r}(*k*_{B}) = 2.7 ppm.^{25}

In 1995, Aziz *et al.*^{7} argued that the values of the thermal conductivity *λ*(*T*), viscosity *η*(*T*), and second density virial coefficient *B*(*T*) of helium, as calculated using *ab initio* input, were more accurate than the best available measurements of these quantities. Subsequently, helium-based AGT measurements of *k*_{B} relied on *ab initio* values of *λ*(*T*) to account for the thermoacoustic boundary layer. Just before the Boltzmann constant was defined in 2019, the lowest-uncertainty measurements of *k*_{B} used either the *ab initio* value of thermal conductivity of helium *λ*_{He}(273.16 K) or the value of *λ*_{Ar}(273.16 K) that was deduced from ratio measurements using *λ*_{He}(273.16 K) as a standard.^{26,27}

In 2019, the unit of temperature, the kelvin, was redefined by assigning the fixed numerical value 1.380 649 × 10^{−23} to the Boltzmann constant, *k*_{B}, when *k*_{B} is expressed in the unit J K^{−1}.^{2,3} Thus, the Boltzmann constant can no longer be measured. However, the thermodynamic temperature of the triple point of water now has an uncertainty of a few parts in 10^{7}, although the best current value is still 273.16 K.^{20}

As discussed in the next section, the techniques for measuring thermodynamic temperatures are evolving rapidly. They are becoming more accurate and easier to implement. We anticipate NMIs will disseminate thermodynamic temperatures instead of ITS-90 at temperatures below 25 K. This would not be possible without the accurate *ab initio* values of the thermophysical properties of helium.

### 2.2. Gas thermometry

#### 2.2.1. Acoustic gas thermometry

During the past two decades, AGT has emerged as the most accurate primary thermometry technique over the temperature range 7–552 K, achieving uncertainties as low as 10^{−6}*T*. AGT experiments were instrumental in measuring the Boltzmann constant for the redefinition of the kelvin,^{28} and have revealed small, systematic errors in the ITS-90.^{20,23} The construction of ITS-90 and the definition of *k*_{B} force *T*_{90} and *T* to be essentially equal at the temperature of the triple point of water *T*_{TPW}; however, the derivative d*T*_{90}/d*T* ≈ 1.0001 at *T*_{TPW}. Figure 3 provides evidence that ITS-90 has errors of $\u223c25\xd710\u22126T$ near water’s boiling point and ∼−35 × 10^{−6}*T* near 173 K. This section is necessarily brief; for an in-depth review of AGT, the reader is referred to Ref. 29.

*T*, and the thermodynamic speed of sound,

*w*, in a gas:

^{30}

^{,}

*k*

_{B}is the Boltzmann constant,

*R*=

*k*

_{B}

*N*

_{A}is the molar gas constant,

*N*

_{A}is the Avogadro constant,

*m*is the average molecular mass of the gas,

*γ*

_{0}is the limiting low-pressure value of

*c*

_{p}/

*c*

_{v}where

*c*

_{p}and

*c*

_{v}are the isobaric and isochoric heat capacities, respectively (this ratio is exactly 5/3 for a monatomic gas),

*p*is the gas pressure, and

*β*

_{a}and

*γ*

_{a}are the temperature-dependent acoustic virial coefficients. Helium-4 or argon gas is typically used, as these are considerably less expensive than other noble gases and available in ultra-pure forms, although xenon has also been used.

^{31}

Most modern realizations of primary AGT determine the speed of sound from the resonance frequencies of the acoustic normal modes in a cavity resonator of fixed and stable dimensions. Resonators have been manufactured from copper, aluminum, and stainless steel, with internal volumes between 0.5 and 3 l. Cavity shapes have either been spherical, quasi-spherical (with smooth, deliberate deviations from sphericity), or cylindrical. The use of diamond turning to produce quasi-spherical resonators (QSRs) with extremely accurate forms (∼1 μm) and smooth surfaces (average surface roughness on the order of 3 nm) has significantly improved performance.^{32} In spherical geometries, the best results are obtained from the radially symmetric acoustic modes, since these possess high quality factors and are relatively insensitive to imperfections in the cavity shape. In cylindrical geometries, the longitudinal plane-wave modes are typically favored.

*T*is determined by using the defined value of

*k*

_{B}and by fitting Eq. (1) to measurements of

*w*

^{2}(

*p*,

*T*) to obtain

*T*and

*β*

_{a},

*γ*

_{a}, …. Alternatively, when accurate,

*ab initio*values of

*β*

_{a},

*γ*

_{a}, … are available,

*T*can be determined from Eq. (1) using a measurement of

*w*

^{2}(

*p*,

*T*) at a single pressure. The terms

*γ*

_{0}and

*k*

_{B}are known exactly;

*m*must be determined by an auxiliary experiment; and

*w*is calculated from the radial acoustic mode frequencies,

*f*

_{a}, of the QSR:

*z*

_{a}are the acoustic eigenvalues, Δ

*f*

_{a}is the sum of the acoustic corrections, and

*V*is the cavity volume. If the longitudinal mode frequencies of a cylindrical cavity are used, the term proportional to

*V*

^{1/3}is replaced with a multiple of the cylinder length.

*f*

_{m}of microwave resonances in the cavity, which are related to the volume through the equation

*c*is the speed of light in vacuum,

*n*in the refractive index of the gas in the cavity, Δ

*f*

_{m}is the sum of the electromagnetic corrections, and

*z*

_{m}are the microwave eigenvalues. The microwave modes do not occur in isolation, being at least three-fold degenerate in perfectly spherical cavities. The smooth deformations of the QSR shape lift these degeneracies, enabling accurate measurement of the individual mode frequencies. A key theoretical result is that (to first order) the mean frequency of these mode groups is unaffected by volume-preserving shape deformations.

^{33}

In diamond-turned QSRs, the relative uncertainty in *V* from the microwave method can be less than 1 × 10^{−6}.^{34} This was made possible by improvements in theory,^{35} resonator shape accuracy, and studies of small perturbations due to probes.^{34} Recently, it has been demonstrated that comparable uncertainties can be achieved with low-cost microwave equipment.^{36,37} Accurate microwave dimensional measurements have also been performed in cylindrical acoustic resonators.^{38}

*w*

_{ref}is the measured speed of sound at a known reference temperature

*T*

_{ref}. Most AGT determinations of (

*T*−

*T*

_{90}) use the relative method. The main advantages are that the molecular mass term,

*m*, cancels in the ratio, and that only the relative volume

*V*/

*V*

_{ref}need be measured. Also, many small perturbations to the acoustic and microwave frequencies (e.g., due to shape deformations) either fully or partially cancel in the ratio. As a result, excellent results can be obtained using resonators with modest form accuracies that would be unsuited to absolute AGT. The disadvantages are that relative AGT propagates underlying errors and uncertainty in

*T*

_{ref}, and can require the apparatus to operate over a wide temperature range when no suitable reference points are nearby.

In both absolute and relative primary AGT, maintaining gas purity is of critical importance. Impurities will shift the average molecular mass of the gas, and hence the speed of sound, by an amount that depends on the mass contrast between the bulk gas and impurity. For example, the speed of sound in helium is ∼16 times more sensitive to water vapor than it is in argon. Impurities can either be present in the gas source or arise from outgassing or leaks in the apparatus itself.

Relative AGT requires only that *m* remain unchanged between the measurements at *T* and *T*_{ref}. Temperature dependence in *m* can arise through several mechanisms: impurities such as water, hydrocarbons, or heavy noble gases can be condensed out at low temperatures; higher temperatures ($>500$ K) cause significant outgassing from the walls of steel resonators.^{39} Gas purity is vastly improved by maintaining a flow of gas (typically <50 μmol/s) through the resonator and supply manifold.

Absolute AGT has more stringent requirements on gas purity than relative AGT. To determine an accurate value for *m*, both the isotopic abundance of the gas and any residual impurities must be quantified. Reactive impurities, including water, can be removed from the source gas using gas purifiers, and noble gas impurities can be removed from helium using a cold trap.^{26} The isotopic ratios ^{36}Ar/^{40}Ar and ^{38}Ar/^{40}Ar in argon, and ^{3}He/^{4}He in helium, have been determined by mass spectrometry, and vary significantly from source to source.^{40} Alternatively, isotopically pure ^{40}Ar gas can be used, although this is only available in small quantities and at great expense.^{41}

The low uncertainty of the AGT technique arises from the excellent agreement between acoustic theory and experiment. The simplicity of Eq. (2) hides a number of temperature-, pressure-, and mode-dependent corrections that constitute the term Δ*f*_{a}. The largest of these are the thermoacoustic boundary layer corrections, which arise from an irreversible heat exchange between the oscillating gas and resonator walls.^{41,42} This effect both lowers the frequency of the acoustic resonances and broadens them; a valuable cross-check of experiment and theory can be made by comparing the predicted and measured resonance widths. The radial-mode boundary layer correction in QSRs is approximately proportional to the square root of the gas thermal conductivity – in cylinders, the gas viscosity also features in the correction.^{43} For most temperature ranges, the uncertainty in these parameters can be considered negligibly small for both helium and argon due to improved *ab initio* calculations (see Sec. 2.5).

AGT measurements are typically conducted on isotherms in a pressure range between 25 and 500 kPa, with the optimum pressure range depending on several factors such as the type of gas, temperature, and particular details of the apparatus.^{44} At low pressures, the accuracy in determining *f*_{a} is compromised by weak acoustic signals, interference from neighboring modes due to resonance broadening, and the need to account for details of the interaction of the gas with the resonator’s walls.^{45} At high pressures, higher-order virial terms are required to account for molecular interactions, and the elastic recoil of the resonator walls becomes increasingly significant. The shell recoil effect, which shifts *f*_{a} in proportion to gas density,^{46} is difficult to predict in real resonators^{47,48} because of the complex mechanical properties of the joint(s) formed when the cavity resonator is assembled.

For this and other reasons, it is not common practice to use Eq. (1) to determine *T* from *w*; instead, the measured data are fitted to low-order polynomials that account for the virial coefficients and perturbations that are proportional to pressure. Isotherm measurements have the advantage of data redundancy and reduced uncertainty, but are very slow to execute, with each pressure point taking several hours. Single-state AGT,^{49} which utilizes low-uncertainty *ab initio* calculations of *β*_{a} and *γ*_{a} in helium, offers a much faster means of primary thermometry.

Figure 3 compares AGT measurements from five countries with ITS-90. The AGT data indicate that ITS-90 has an error of $\u223c25\xd710\u22126T$ near water’s boiling point and ∼−35 × 10^{−6}*T* near 173 K. Near *T*_{TPW}, the derivative d*T*_{90}/d*T* ≈ 1 + 1.0 × 10^{−4}. This implies that heat-capacity measurements made using ITS-90 will generate values of the heat capacity that are 0.01% larger than the true heat capacity. However, we are not aware of heat capacity measurement uncertainties as low as 0.01%.

Prior to the AGT publications shown in Fig. 3, Astrov *et al.* corrected an estimate used in their CVGT. They had used measurements of the linear thermal expansion of a metal sample to estimate the thermal expansion of the volume of their CVGT “bulb.” Using additional expansion measurements, Astrov *et al.* corrected their *T* − *T*_{90} results. They now agree, within combined uncertainties, with the AGT data.^{56} (Because AGT uses microwave resonances to measure the cavity’s volume *in situ*, it is not subject to errors from auxiliary measurements of thermal expansion.)

#### 2.2.2. Dielectric constant gas thermometry

^{57,58}and later improved by PTB,

^{59,60}is now a well-established method of primary thermometry. The basic idea of DCGT is to replace the density in the equation of state of a gas by the relative permittivity (dielectric constant)

*ɛ*

_{r}and to measure it by the relative capacitance changes at constant temperature:

*C*

_{c}(

*p*) is the capacitance of the capacitor at pressure

*p*and

*C*

_{c}(0) that at

*p*= 0 Pa, and

*κ*

_{eff}is the effective isothermal compressibility which accounts for the dimensional change of the capacitor due to the gas pressure. In the low-pressure (ideal gas) limit, the working equation can be simply derived by combining the classical ideal-gas law and the Clausius–Mossotti equation:

*A*

_{ɛ}. For a real gas in a general formulation including electric fields, both input equations are power series:

*B*(

*T*),

*C*(

*T*), and

*D*(

*T*) are the second, third, and fourth density virial coefficient, respectively,

*ρ*is the molar density, and

*b*,

*c*,

*d*and

*B*

_{ɛ},

*C*

_{ɛ},

*D*

_{ɛ}are both called the second, third, and fourth dielectric virial coefficient, respectively. The form used in Eq. (8) comes from the tradition of DCGT

^{57,59}of factoring out

*A*

_{ɛ}so that

*b*,

*c*, and

*d*have the same units as

*B*,

*C*, and

*D*. Conversely,

*ab initio*calculations naturally provide the quantities

*B*

_{ɛ},

*C*

_{ɛ}, and

*D*

_{ɛ}.

*ɛ*

_{r}with the relative capacitance change corresponding to Eq. (5). This leads to a power expansion in terms of Ξ = (Δ

*C*

_{c}/

*C*

_{c})/(Δ

*C*

_{c}/

*C*

_{c}+ 3):

The higher-order terms contain combinations of both the dielectric and density virial coefficients and the compressibility. Equation (10) up to the fourth order can be found in Ref. 61.

DCGT works as a primary thermometer if the molar polarizability *A*_{ɛ} and virial coefficients contained in Eq. (10) are known from fundamental principles or independent measurements with sufficient accuracy. The effective compressibility *κ*_{eff} is also required. For classical DCGT, where isotherms are measured and the data are extrapolated to zero pressure via least-squares fitting, only *A*_{ɛ} and *κ*_{eff} are mandatory. This was the way thermodynamic temperature was determined for decades.^{57,59,60} Consequently, in classical DCGT, *ab initio* data on virial coefficients serve as a consistency check or conversely DCGT is used for determination of virial coefficients to check theory.^{61} Since the theoretical calculations of the virial coefficients for helium improved drastically, it is now possible to use higher-order virial coefficients from theory to reduce the number of fitting coefficients or even to use the working equation directly without fitting and to determine temperature at each pressure point via the rearranged working equation. Recently, all three approaches have been tested and compared.^{62} Especially, the point-by-point evaluation is a shift of paradigm and at the moment only possible for helium, where the uncertainty of the *ab initio* calculations, especially of the second density virial coefficient, is small enough. Nevertheless, for other gases not only the virial coefficients but also the molar polarizabilities determined via DCGT have comparable or smaller uncertainties than *ab initio* calculations.^{63} This is a field of potential improvement of theory already started with calculations of *A*_{ɛ} for neon^{64,65} and for argon.^{66}

DCGT was operated in the temperature range from 2.5 K to about 273 K using helium-3, helium-4, neon,^{67} and argon.^{68} All noble gases have the advantage that the molar polarizability is independent of temperature at a level of precision far beyond that of state-of-the-art experiments.^{69}

Besides the use of dielectric measurements in primary thermometry, accurate determinations of polarizability and virial coefficients of noble gases and molecules using gas-filled capacitors have a much longer tradition. These setups, very similar to DCGT, use thermodynamic temperature as one of the input parameters. A complete overview of measurements cannot be given here. Already a very broad overview of existing data, partly at radio frequencies, was summarized by NBS in the 1950s.^{70} In the following decades,^{71} different institutes with changing teams performed measurements until the early 1990s.^{72} In the year 2000, NIST started measurements on gases using capacitors resulting in the most accurate values for the measured molecules.^{73,74} Very recently, PTB established a setup for separate measurement of dielectric and density virial coefficients using a combination of Burnett expansion techniques and DCGT.^{75} The focus of this setup is the determination of properties of energy gases such as hydrogen-methane mixtures in the context of the transition to renewable energy. The setup will also provide lower-uncertainty tests of the *ab initio* calculations of the dielectric and density virial coefficients of the noble gases.

For primary thermometry, most significant recent improvements in DCGT have been achieved by independent determination of *κ*_{eff} using resonant ultrasound spectroscopy around 0 °C and an optimal choice of capacitor materials.^{76} For the Boltzmann experiment with measuring pressures of up to 7 MPa, tungsten carbide was the ideal choice, while at low temperatures beryllium copper was used together with an extrapolation method. Relative uncertainties for *κ*_{eff} in terms of temperature on the level of 1 ppm near 0 °C have been achieved. Equally important are the improvements in pressure measurement. In contrast to AGT, where pressure is a second-order effect, in DCGT *ɛ*_{r} is directly linked to pressure. Therefore, the relative uncertainty in pressure is transferred to a relative uncertainty in temperature. The major steps here are discussed in Sec. 2.3.2 regarding the mechanical pressure standard developed at PTB in the framework of the Boltzmann constant determination.^{77} These systems with relative uncertainties on the level of 1 ppm at pressures up to 7 MPa have been used to calibrate commercially available systems for pressures up to 0.3 MPa with relative uncertainties between 3 and 4 ppm. The dominant uncertainty component in DCGT measurements is the standard deviation of the capacitance measurement. The typical relative uncertainty in terms of temperature connected to this component is on the order of 5 ppm for the low temperature range but was reduced to the 1 ppm level in the case of the Boltzmann experiment at about 0 °C.^{78} Finally, one problem in DCGT using helium is the very small molar polarizability compared to all other gases and molecules. Therefore, special care must be taken concerning impurities and here an especially severe issue is contamination with water.

The polarizability of water at frequencies of capacitance bridges and microwave resonators (see Sec. 5.1.2) is about a factor of 160 larger than that of helium. At cryogenic temperatures, water contamination in the gas phase is naturally reduced by outfreezing, but especially at room temperature the whole measuring setup as well as the gas purifying system must be highly developed. Furthermore, pollution with other noble gases must be treated carefully because they cannot be extracted by getters and filters. Ideally, a mass-spectrometer should be used for the detection of noble gas impurities to allow for an upper estimate of the uncertainty due to gas purity. In summary, with DCGT in the low temperature range from 4 to 25 K uncertainties near 0.2 mK for thermodynamic temperature are achievable. At around 0 °C, the smallest uncertainty for DCGT was achieved during the determination of the Boltzmann constant.^{78} Converted to an uncertainty for thermodynamic temperature, this becomes about 0.5 mK.

In the intermediate range, the uncertainties are larger (between 1 and 2 mK at 200 K^{68}). The main restriction of the present low-temperature setup is the limited pressure range at intermediate temperatures. A measurement of high-pressure isotherms in this range is planned. Together with improved *ab initio* calculations for the second virial coefficients of argon and neon, a single-state version of DCGT might be possible, in analogy with single-state AGT. This could result in a significant reduction in both uncertainty and measurement time.

#### 2.2.3. Refractive index gas thermometry

*ɛ*

_{r}or of the refractive index

*n*in powers of the molar density

*ρ*, that is Eq. (9) in the case of DCGT, and the Lorentz–Lorenz equation

*A*

_{μ}/

*A*

_{ɛ}≈ −1.53 × 10

^{−5}for He,

*B*

_{ɛ}=

*B*

_{R},

*C*

_{ɛ}=

*C*

_{R}, etc.

^{79}Except for the small magnetic-permeability term

*A*

_{μ}(which is well-known from theory for helium

^{80}), low-frequency measurements of

*n*and of

*ɛ*

_{r}are analyzed using the same

*ab initio*constants. RIGT determines the thermodynamic temperature

*T*by combining measurements of the pressure

*p*with the density virial equation of state, Eqs. (7) and (11). The density is eliminated from both equations, either numerically or by iteration, to obtain

The constants *B*, *B*_{R}, *C*, *C*_{R}, etc. that appear in the higher-order terms of Eq. (12) are obtained either from theory or from fitting measurements of *n*^{2}(*p*) on isotherms. [DCGT determines *T* using a version of Eq. (12) in which *ɛ*_{r} replaces *n*^{2}.]

Here, we focus on RIGT conducted at microwave frequencies as developed by Schmidt *et al.*^{81} and as recently reviewed by Rourke *et al.*^{82} These authors determined *n* from measurements of the microwave resonance frequencies *f*_{m} of a gas-filled, metal-walled, quasi-spherical cavity. Typical frequencies ranged from 2.5 to 13 GHz; for this range, the frequency dependence of *n* in the noble gases is negligible. As discussed in Sec. 2.3.1, RIGT has also been realized at optical frequencies in the context of pressure standards.^{83} For helium, the corrections of *A*_{ɛ} and *B*_{R} from zero frequency to optical frequencies have been calculated *ab initio*.^{79,84}

*n*is:

*g*accounts for the penetration of the microwave fields into the cavity’s walls. Usually,

*g*is determined from measurements of the half-widths of the resonances; its contribution to uncertainties is small. The term

*κ*

_{eff}

*p*accounts for the temperature-dependent change of the cavity’s volume in response to the gas pressure

*p*. Often, the uncertainty of

*κ*

_{eff}is the largest contributor to the uncertainty of RIGT. To make this explicit, we manipulate Eqs. (12) and (13) to obtain:

*κ*

_{eff}

*RTn*

^{2}/[3(

*A*

_{ɛ}+

*A*

_{μ})] ≈ 0.007 for a copper-walled cavity immersed in helium near

*T*

_{TPW}. (This estimate assumes that the cavity’s walls are homogeneous and isotropic; therefore,

*κ*

_{eff}=

*κ*

_{T}/3 where

*κ*

_{T}is the isothermal compressibility of copper.) Thus, a relative uncertainty

*u*

_{r}(

*κ*

_{eff}) = 0.01 contributes the relative uncertainty

*u*

_{r}(

*T*

_{TPW}) = 70 × 10

^{−6}to a RIGT determination of

*T*

_{TPW}. In the approximation

*n*

^{2}≈ 1, this uncertainty contribution is a function of

*T*×

*κ*

_{eff}(

*T*), but it is not a function of the pressures measured on an isotherm. Because

*T*×

*κ*

_{eff}(

*T*) decreases with

*T*, RIGT is more attractive at cryogenic temperatures than near or above

*T*

_{TPW}.

Recently, two independent groups explored a two-gas method for measuring *κ*_{eff} of assembled RIGT resonators.^{17,85} Ideally, two-gas measurements would replace measurements of *κ*_{T} of samples of the material comprising the resonator’s wall and also models for the cavity’s deformation under pressure. Both groups relied on new, accurately measured and/or calculated values of the density and refractivity virial coefficients of neon or argon.^{61,79} Using helium and argon, Rourke determined *κ*_{eff} at *T*_{TPW} with the remarkably low uncertainty *u*_{r}(*κ*_{eff}) = 9.6 × 10^{−4}.^{85} Madonna Ripa *et al.* combined helium and neon data to reduce the uncertainty contribution from *κ*_{eff} to their determinations of *T* at the triple points of O_{2} ($\u224854$ K), Ar ($\u224884$ K), and Xe ($\u2248161$ K).^{17} They reported “partial success” and suggested that a revised apparatus using both gases and operating at higher pressures (*p* > 500 kPa) would obtain lower-uncertainty determinations of *T*. They also noted that the two-gas method requires twice as much RIGT data, accurate pressure measurements, and dimensional stability between gas fillings.

Rourke’s review of RIGT^{82} noted five groups implementing RIGT using microwave technology. In contrast, we are aware of only one group (at PTB) implementing DCGT.^{62} The relative popularity of RIGT results from the commercial availability of vector analyzers that can measure microwave frequency ratios with resolutions of 10^{−9}. To our knowledge, using commercially available capacitance bridges, the best attainable capacitance ratio resolution is 70 × 10^{−9}.^{86} To attain higher resolution for DCGT, PTB developed a unique bridge that measures capacitance ratios with a resolution of order 10^{−8} in a 1 s averaging time. To achieve this specification, the PTB bridge must operate at 1 kHz and both the standard (evacuated) capacitor and the unknown (gas-filled) capacitor must have identical construction and be located in the same thermostat.^{87}

Figure 4 illustrates the several strategies being explored for acquiring RIGT data. Absolute RIGT acquires many (*p*, *n*) data on an isotherm and determines *T* via Eq. (14). This method requires state-of-the-art, absolute pressure measurements; therefore, the pressure gradient between the gas-filled cavity and the manometer (normally at ambient temperature) is required.^{88} Uncertainty budgets for absolute RIGT can be found in Refs. 17 and 85.

*T*

_{ref}for which the thermodynamic temperature is already well known, and (2) an unknown isotherm for which

*T*will be determined. As suggested in the lower panel of Fig. 4, one flavor of rRIGT determines

*T*/

*T*

_{ref}by determining the low-pressure limit of the ratio of slopes

^{81}

*T*

_{ref}and

*T*are low temperatures, where the pressure deformation of the cavity

*κ*

_{eff}

*p*is small, this strategy circumvents the problem of accurately determining

*κ*

_{eff}.

Single-pressure RIGT (spRIGT) measures (*p*, *n*, *T*) and (*p*, *n*, *T*_{ref}) and determines *T* from $T/Tref\u2248(nT2\u22121)/(nTref2\u22121)$. This strategy entirely avoids accurate pressure measurements; instead, the pressure in the cavity is required to be identical when *n* is measured at *T* and *T*_{ref} and the pressure (actually, the density of the gas) must be sufficiently low that an approximate pressure is adequate for making the virial corrections. This strategy was used by Gao *et al.* for RIGT between the triple point of neon (*T*_{ref} ≈ 24.5 K) and 5 K.^{89} After establishing *T*_{ref} by acoustic thermometry, they claimed the uncertainties of this implementation of RIGT were smaller than the uncertainties of ITS-90.^{90}

When constant-frequency RIGT (cfRIGT) is implemented, the pressure in the cavity is changed to keep the refractive index constant as the temperature is changed from *T*_{ref} to *T*. In this case, *T*/*T*_{ref} ≈ *p*(*T*, *n*)/*p*(*T*_{ref}, *n*).^{91} This scheme minimizes the frequency-dependent effects of the coaxial cables on the microwave determination of *T*/*T*_{ref}.

To economically search for measurement or modeling errors, one can obtain three redundant values of *T*/*T*_{ref} by measuring microwave frequencies at four judiciously chosen values of (*p*, *n*). Two measurements are made on the isotherm *T*_{ref} at the values (*p*_{1}, *n*_{1}) and (*p*_{2}, *n*_{2}). Two other measurements are on the isotherm *T* at (*p*_{2}, *n*_{3}) and (*p*_{3}, *n*_{1}). spRIGT connects the points (*p*_{2}, *n*_{2}) and (*p*_{2}, *n*_{3}). cfRIGT connects the points (*p*_{1}, *n*_{1}) and (*p*_{3}, *n*_{1}). All four points are used to approximately implement rRIGT via Eq. (15).

Compared with other forms of gas thermometry, relative RIGT has significant advantages at low temperatures. We have already emphasized the availability of microwave network analyzers and the possibility of avoiding state-of-the art pressure measurements. By measuring several microwave resonance frequencies at each state, certain imperfections of the measurements and modeling can be detected. Comparisons of the frequencies of transverse electric (TE) and transverse magnetic (TM) microwave modes might detect the presence of dielectric films such as oxides, oil deposits, or adsorbed water on the cavity’s walls.^{92} Because relative RIGT relies on microwave frequency ratios, the precise shape of the cavity is unimportant. Cavity shapes other than quasispheres may be advantageous in particular applications.

RIGT is simpler and more rugged than relative AGT (rAGT) because RIGT requires neither delicate acoustic transducers nor acoustic ducts. However, RIGT is unlikely to replace rAGT at ambient and higher temperatures because RIGT is more sensitive to the cavity’s dimensions than rAGT by the factor 1/(*ɛ*_{r} − 1), which typically ranges from 200 to 20 000. Furthermore, microwave RIGT is especially sensitive to polar impurities. Adding 1 ppm (mole fraction) of water vapor to dilute argon gas at 293 K will increase the dielectric constant of the gas by 18 ppm and increase the square of the speed of sound by 0.12 ppm. If the water vapor were undetected, these changes would reduce argon’s apparent RIGT temperature by 18 ppm and increase argon’s apparent rAGT temperature by 0.12 ppm. For helium, the corresponding temperatures are reduced by 145 ppm and 4 ppm.

#### 2.2.4. Constant volume gas thermometry

The website of the International Bureau of Weights and Measures includes a document (“*Mise en pratique*…”) that indicates how the SI base unit, the kelvin, may be realized in practice using four different versions of gas thermometry.^{93} Surprisingly, this document omits CVGT, the version of gas thermometry that was the primary basis of ITS-90. In this section, we briefly describe the operation of a particular realization of CVGT and the inconsistent results it generated. This may explain why CVGT was omitted from the *Mise en pratique*. We mention the post-1990 theoretical and experimental developments that suggest an updated realization of CVGT might generate very accurate realizations of the kelvin.

_{NIST90}.” The heart of CV

_{NIST90}was a metal-walled, cylindrical cavity (“gas bulb”;

*V*≈ 407 cm

^{3}) attached to a “dead space” comprised of a capillary leading from the bulb to a constant-volume valve at ambient temperature. The valve separated the gas bulb from a pressure-measurement system. A typical temperature measurement using CV

_{NIST90}began by admitting

*N*

_{r}≈ 0.0023 mol of helium into the gas bulb at a measured reference pressure (

*p*

_{r}≈ 13 kPa) and a measured reference temperature (

*T*

_{r}≈

*T*

_{TPW}).

^{94,95}Then, the valve was closed to seal the helium in the gas bulb and dead space. The bulb was moved into a furnace that was maintained at the unknown temperature

*T*to be determined by CVGT. After the gas bulb equilibrated, the valve was opened to measure the pressure

*p*again. The temperature ratio

*T*/

*T*

_{r}was determined by applying the virial equation at each temperature:

*T*/

*T*

_{r}is determined, in leading order, by the three ratios:

*p*/

*p*

_{r},

*V*

_{T}/

*V*

_{r}, and

*N*

_{r}/

*N*

_{T}. For CV

_{NIST90},

*N*

_{r}/

*N*

_{T}≠ 1 because a tiny quantity of helium flows from the bulb into the capillary when the bulb is moved into the furnace. This quantity was calculated using the measured temperature distribution along the capillary. For CV

_{NIST90},

*V*

_{T}/

*V*

_{r}was calculated using auxiliary measurements of the linear thermal expansion of samples of the platinum–rhodium alloy comprising the gas bulb. These samples had been cut out of the gas bulb after completing all the CVGT measurements.

The simplicity of Eq. (16) hides the many complications of CVGT. We mention three examples. (1) During pressure measurements, helium outside the gas bulb was maintained at the same pressure as the helium inside the gas bulb. (2) Thermo-molecular and hydrostatic pressure gradients in the capillary were taken into account. (3) At high temperatures, creep in the gas bulb’s volume was detected by time-dependent pressure changes; the pressure was extrapolated back in time to its value when the bulb was placed in the furnace.

We denote the second most recent realization of NBS/NIST’s relative CVGT by “CV_{NBS76}.”^{96} Both CV_{NIST90} and CV_{NBS76} shared apparatus and many procedures. However, Ref. 94 lists 11 significant changes. Here, we mention only one. CV_{NIST90}’s two cylindrical gas bulbs had been fabricated entirely from sheets of (80 wt. % Pt + 20 wt. % Rh) alloy. The sides and bottom of CV_{NBS76}’s gas bulb were fabricated from the same alloy; however, the top of the bulb was inadvertently fabricated from (88 wt. % Pt + 12 wt. % Rh) alloy. Perhaps the slight differences in thermal expansions of these alloys led to an anomalous thermal expansion of the volume of CV_{NBS76}’s gas bulb.

Unfortunately, the results from CV_{NIST90} and CV_{NBS76} were inconsistent, within claimed uncertainties, in the range of temperature overlap (505 K $\u2264T\u2264$ 730 K). An approximate expression for the differences is: *T*_{NIST90} − *T*_{NBS76} ≈ 0.090 × (*T*/K − 400) mK. This inconsistency was not explained by the authors of CV_{NIST90} nor by the authors of CV_{NBS76}. Furthermore, the authors did not assert the more recent CV_{NIST90} results were more accurate than the earlier CV_{NBS76} results. The working group that developed ITS-90 had no other data, from NIST or elsewhere, that were suitable for resolving the inconsistency. Therefore, the working group required ITS-90 to be the average of *T*_{NIST90} and *T*_{NBS76} in the overlap range.^{97}

In the range 2.5–308 K, ITS-90 relied, in part, on another realization of CVGT that had a troubled history. Astrov *et al.* deduced the thermal expansion of their copper gas bulb’s volume from measurements of the linear thermal expansion of copper samples taken from the block used to manufacture their bulb.^{98} However, the thermal expansion data were inconsistent with other data for copper. Astrov’s group repeated the thermal expansion measurements using another (better) dilatometer. The more recent expansion data, published in 1995, changed the values of *T* by more than 50 × 10^{−6}*T* in the range 130 K $<T<$ 180 K, where the uncertainties had been estimated as $\u226426\xd710\u22126T$.^{56}

Recently, a working group of the Consultative Committee for Thermometry reviewed primary thermometry below 335 K.^{20} Astrov’s revised CVGT values are close to the current consensus, which is primarily based on AGT and DCGT. The working group retained three other low-temperature realizations of CVGT. Post-1990 AGT measurements of *T* − *T*_{90} near 470 and 552 K indicate that CV_{NIST90} is indeed more accurate than CV_{NBS76}.^{50} Despite the fact that CVGT was the primary basis for the ITS-90, the *Mise en pratique* does not include CVGT. We speculate that no temperature metrology group is pursuing CVGT because: (1) CVGT is complex, (2) Astrov *et al.*’s thermal expansion problem, (3) unexplained problems with NBS/NIST’s CVGT, and (4) rapid advances in other versions of gas thermometry.

We now ask: is CVGT a viable method of primary thermometry today? The gas bulb of a modern CVGT would incorporate feedthroughs to enable measuring microwave resonance frequencies of the bulb’s cavity. The resonance frequencies would determine the bulb’s volumetric thermal expansion, thereby avoiding auxiliary measurements of linear thermal expansion and also avoiding the assumption of isotropic expansion. If the bulb incorporated a valve and a differential-pressure-sensing diaphragm, the dead-space corrections would vanish. (The diaphragm’s motion could be detected using optical interferometry.) Today, the *ab initio* values of *B*(*T*) would reduce the uncertainty component from *B*(*T*) to near zero. A contemporary CVGT could operate at $\u223c5\xd7$ higher helium densities than published experiments without generating significant uncertainties from either the virial coefficients or from pressure-ratio measurements. The higher density, together with simultaneous pressure and microwave measurements, might enable separation of the bulb’s creep from contamination by outgassing. Most outgassing contaminants affect helium’s dielectric constant, refractivity, and speed of sound much more than they affect helium’s pressure, an advantage of CVGT. However, CVGT inherently uses fixed aliquots of gas. Therefore, CVGT cannot benefit from flowing gas techniques that have been used, for example, in high-temperature AGT.^{50} In summary, contemporary CVGT could be competitive with other forms of primary gas thermometry, with a possible exception at the highest temperatures, where flowing gas might be required to maintain gas purity.

### 2.3. Pressure metrology

Traditionally, standards based on the realization of the *mechanical* definition of pressure, the normal force applied per unit area onto the surface of an artifact, include pressure balances and liquid column manometers. The combined overall pressure working range of these instruments extends over seven orders of magnitude, roughly between 10 Pa and 100 MPa. Liquid column manometers achieve their best performance, with relative standard uncertainty as low as 2.5 ppm, near their upper working limit at a few hundred kPa.^{99} With a few notable exceptions, the typical relative standard uncertainty of pressure balances spans between nearly 1 × 10^{−3} at 10 Pa, the lowest end of their utilization range, down to 2–3 ppm in the range between 100 kPa and 3 MPa.^{99,100} One such exception is the remarkable achievement of a relative standard uncertainty as low as 0.9 ppm for the determination of helium pressures up to 7 MPa,^{77} though this achievement required the extensive dimensional characterization, and the cross-float comparison, of the effective areas of six piston–cylinder sets manufactured to extraordinarily tight specifications, with a research effort lasting several years. In spite of this outstanding result, the accurate characterization of pressure balances is challenging, due to the complexity of the dimensional characterization of the cross-sectional area of piston–cylinder assemblies, which includes finite-element modeling of their deformation under pressure.^{101,102} International comparisons periodically provide realistic estimates of the average uncertainty of realization of primary standards among NMIs. In 1999, a comparison of primary mechanical pressure standards in the range 0.62 MPa $<p<$ 6.8 MPa, involving five NMIs leading in pressure metrology exchanging a selected piston–cylinder set, was completed.^{103} The resulting differences Δ*A*_{eff} ≡ 10^{6}(*A*_{eff}/⟨*A*_{eff}⟩ − 1) of the effective area *A*_{eff} of the piston from the reference value ⟨*A*_{eff}⟩ spanned beyond their combined uncertainties with such significant spread to show that the pressure standards realized by different NMIs were mutually inconsistent.

These inconsistencies strengthened the motivation for the development of standards realizing a *thermodynamic* definition of pressure by the experimental determination of a physical property of a gas having a calculable thermodynamic dependence on density, combined with accurate thermometry. This possibility was initially proposed in 1998 by Moldover,^{104} who envisaged, already at that time, the potential of first-principles calculation to accurately predict the thermodynamics and electromagnetic properties of helium and the maturity of experiments determining the dielectric constant using calculable capacitors. The metrological performance of thermodynamic pressure standards has continuously improved over the last two decades to become increasingly competitive in terms of accuracy, providing important alternatives that may test the exactness of the mechanical standards discussed above and eventually replace some of them. Also, due to their reduced complexity and bulkiness, simplified versions of thermodynamics-based standards may be more flexibly adapted to specific technological and scientific applications of pressure metrology. The best-performing recent realizations of gas-based pressure standards include measurements of the dielectric constant using capacitors and of the refractive index at microwave and optical frequencies, respectively using resonant cavities and Fabry–Pérot refractometers. In Secs. 2.3.1 and 2.3.2, we separately discuss the most notable of these developments depending on the pressure range of their application.

#### 2.3.1. Low pressure standards (100 Pa to 100 kPa)

In the low vacuum regime, several experimental methods are available which may provide alternative routes for traceability to the pascal. For the cases involving optical measurements, these methods include: (1) refractometry (interferometry), implemented in various configurations that employ single or multiple cavities or cells with fixed or variable path lengths; (2) line-absorption methods. The achievements and perspectives of all these methods were recently reviewed.^{105}

At present, Fabry–Pérot refractometry with fixed length optical cavities (FLOC) has demonstrated the lowest uncertainty for the realization of pressure standards near atmospheric pressure and down to 100 Pa. In principle, the uncertainty of this method is limited by several optical and mechanical effects, most importantly by the change in the length of the cavity due to compression by the test gas, with the same sensitivity to the imperfect estimate of the compressibility *κ*_{T} that affects RIGT. However, this major uncertainty contribution may be drastically reduced, though not completely eliminated, by measuring the pressure-induced length change of a second reference FLOC monolithically built on the same spacer, which is kept continuously evacuated. In 2015, a dual-cavity FLOC achieved an extremely accurate determination of the refractive index of nitrogen at *λ* = 632.9908 nm, *T* = 302.9190 K and 100.0000 kPa by reference to the pressure realized by a primary standard mercury manometer, and using refractive index measurements in helium to determine the compressibility.^{106} A comparison of the pressures determined by the nitrogen refractometer with the mercury manometer below the primary calibration point at 100 kPa down to 100 Pa showed relative differences within 10 ppm. A direct comparison between laser refractometry with nitrogen and a mercury manometer was realized one year later also at NIST.^{18} The comparison showed relative differences between these instruments within 10 ppm over the range between 100 Pa and 180 kPa. The laser refractometer outperforms the precision and repeatability of the liquid manometer and demonstrates a pressure transfer standard below 1 kPa that is more accurate than its current primary realization. Such remarkably low uncertainty also favorably compares to the best dimensional characterization and modeling of non-rotating piston–cylinder assemblies.^{107}

In 2017, more accurate measurements in helium and nitrogen were performed between 320 and 420 kPa using a triple-cell heterodyne interferometer referenced to a carefully calibrated piston gauge, showing relative differences within 5 ppm with uncertainties on the order of 10 ppm.^{83} Some pressure distortion errors affecting FLOC might in principle be eliminated by refractive index measurement with a variable length optical cavity (VLOC). The realization of this technique requires extremely challenging dimensional measurements, with displacements on the order of 15 cm that must be determined with picometer uncertainty.^{108} Gas modulation techniques, with the measuring cavity frequently and repeatedly switched between a filled and evacuated condition, have been recently developed,^{109,110} aiming at the reduction of the effects of dimensional instabilities and other short- and long-term fluctuations that affect Fabry–Pérot refractometers. A novel realization of an optical pressure standard, based on a multi-reflection interferometry technique, has also been recently developed, demonstrating the possible realization of the pascal with a relative standard uncertainty of 10 ppm between 10 and 120 kPa.^{111} Optical refractometry for pressure measurement is also being pursued at other NMIs.^{112,113}

With accurate pressure measurement, these optical methods can yield the thermodynamic temperature, becoming another approach to RIGT. This was demonstrated in Ref. 83, where a refractometer was used to measure the Boltzmann constant (albeit with higher uncertainty than AGT and DCGT measurements) prior to the 2019 SI redefinition.

At microwave frequencies, the realization of a low-pressure standard requires a substantial enhancement in frequency resolution. Recently, it was demonstrated by Gambette *et al.* that by coating the internal surface of a copper cavity with a layer of niobium, and working at temperatures below 9 K where niobium becomes superconducting, pressures between 500 Pa and 20 kPa can be realized very precisely.^{114,115} The overall relative standard uncertainty of this method is currently 0.04%, with the largest contribution from non-state-of-the-art thermometry, which is likely to be substantially reduced in future work.

#### 2.3.2. Intermediate pressure standards (0.1–7 MPa)

Differently than initially envisaged, the first realization of a thermodynamic pressure standard was not obtained by capacitance measurements, but using a microwave resonant cavity working in the GHz frequency range, i.e., by a RIGT method. A main motivation for this choice was the development of quasi-spherical microwave resonators, whose internal triaxial ellipsoidal shape slightly deviates from that of a perfect sphere.^{92} This particular geometry resolved the intrinsic degeneracy of microwave modes, allowing enhanced precision in the determination of resonance frequencies.

By 2007, Schmidt *et al.*^{81} demonstrated a pressure standard based on the measurement of the refractive index of helium to achieve overall relative pressure uncertainty *u*_{r}(*p*) within 9 × 10^{−6} between 0.8 and 7 MPa. At the upper limit of the pressure range, the uncertainty was dominated by the uncertainty of the isothermal compressibility *κ*_{T} of maraging steel, which was determined using resonance ultrasound spectroscopy (RUS).^{76} Recently, Gaiser *et al.*^{19} realized Moldover’s original proposal of a capacitance pressure standard using DCGT techniques that they had refined during their measurements of the Boltzmann constant. They achieved the remarkably low uncertainty *u*_{r}(*p*) = 4.4 × 10^{−6} near 7 MPa. Recently, the same experimental data were re-analyzed to take advantage of the increased accuracy of the *ab initio* calculation of the second density virial coefficient *B* of He,^{11} reducing the overall uncertainty of the capacitance pressure standard to *u*_{r}(*p*) = 2.2 × 10^{−6}.^{116}

At pressures below 1 MPa, the uncertainty of the realization of a pressure standard based on DCGT or RIGT with helium is limited by the resolution of relative capacitance or frequency measurements. This limit would be immediately reduced by up to one order of magnitude by using, instead of helium, a more polarizable gas like neon or argon. However, while a significant improvement of the interaction potential, and hence of the *ab initio* calculated *B*, has recently been achieved for neon^{117} and for Ar,^{118} it is not likely that the best available calculations of the molar polarizability *A*_{ɛ} of neon^{65} or argon^{66} can be improved sufficiently to replace experiment in the near future. However, an experimental estimate of *A*_{ɛ} of both neon and argon was obtained by comparative DCGT measurements relative to helium, with relative uncertainty of 2.4 ppm,^{63} and may now be used for the realization of pressure standards with other apparatus. For similar purposes, the ratio of the refractivity of several monatomic and molecular gases, namely Ne, Ar, Xe, N_{2}, CO_{2}, and N_{2}O, to the refractivity of helium was determined at *T* = 293.15 K, *λ* ∼ 633 nm, with standard uncertainty within 16 × 10^{−6}, using interferometry.^{119} At pressures higher than a few MPa, the imperfect determination of the deformation of the cavity under pressure would impact the overall uncertainty of a pressure standard based on RIGT or DCGT. One way to overcome this limit would be to measure the refractivity of two gases at identical values of an unknown pressure using a single apparatus at a known temperature. If the refractivity of both gases were known (either from *ab initio* calculations or reference measurements), the two measurements would determine both the effective compressibility *κ*_{T} of the apparatus and the unknown pressure. The same strategy is also applied to increase the upper pressure range where refractometry methods like FLOC can be applied, though use of helium for the determination of distortion effects requires correcting for diffusion within the glasses used for the construction of these apparatuses.^{120}

### 2.4. High pressures and equation of state

Up to this point, we have considered interactions between temperature and pressure standards and the rigorously calculated, low-density properties of the noble gases including the polarizability and second and third density and dielectric virial coefficients. We now compare *ab initio* calculations with measurements at pressures above 7 MPa and at correspondingly higher densities. The literature includes temperature-dependent values of 6 density virial coefficients of helium,^{121} 7 acoustic virial coefficients of krypton,^{122} and 6 density virial coefficients of argon.^{123} These calculations used the best *ab initio* two-body and nonadditive three-body potentials that were available at the time of publication. Many-body non-additive potentials involving four or more bodies, which are needed for the exact calculation of virial coefficients from the fourth onwards, are not available and are generally neglected, resulting in an uncontrolled approximation. Here, we compare measurements of the density of helium *ρ*_{meas}(*p*, *T*) with values calculated *ab initio*. This comparison avoids fitting *ρ*_{meas}(*p*, *T*) to the VEOS because such fits yield highly correlated values for the separate virial coefficients, each with large uncertainties. Later in this section, we comment on comparisons using speed-of-sound data.

Measurements of gas densities with uncertainties below 0.1% are expensive and rare because they are not required for chemical and mechanical engineering. The uncertainties of most process models are dominated by imperfect models of equipment (heat exchangers, compressors, distillation columns, etc.) and/or imperfect knowledge of the composition of feedstocks and products. An example of a demanding application of gas density and composition measurements is custody transfer of natural gas as it flows through large pipelines near ambient temperature and at high pressures (e.g., 7 MPa). An international comparison among NMIs achieved a *k* = 2 volumetric flow uncertainty of only 0.22%.^{124} In this context, density and composition measurements with uncertainties of order 0.1% are satisfactory for converting volumetric flows into mass flows and heating values.

In Fig. 5, the remarkable data of McLinden and Lösch-Will are used to test the *ab initio* VEOS of helium in the ranges 1 MPa $<p<$ 38 MPa and 223 K $<T<$ 323 K.^{125} These data were acquired using a magnetic suspension densimeter. A weigh scale determined the buoyant forces on two “sinkers” immersed in the helium. The data are precise, well-documented, and traced to SI standards with a claimed, *k* = 2, density uncertainty of 0.015% + 0.001 kg/m^{3} at the temperature extremes and at the highest density. These features attracted previous comparisons with theory.^{21,121,126}

For the present comparison, where recently published theoretical values of the virial coefficients are used, we converted the measured temperatures from the ITS-90 to thermodynamic temperatures using Ref. 20 and we converted the measured mass densities to molar densities using the defined value of the universal gas constant and the molar mass for McLinden and Lösch-Will’s helium sample. At densities below $\u223c4000$ mol/m^{3}, the uncertainties and the values of (*ρ*_{meas}/*ρ*_{calc} − 1) diverge on isotherms as *ρ*^{−1} and/or *p*^{−1}. (See Fig. 5.) These low-density divergences result from time-dependent drifts in the zeros of the densimeter and/or pressure transducer. Because the divergences contain more information about the apparatus than about helium’s VEOS, we do not discuss them.

^{3}, we compared the

*ρ*

_{meas}(

*p*,

*T*) data of McLinden and Lösch-Will

^{125}with the values of $\rho calc*(p,T)$ that are implicitly defined by the truncated VEOS:

*B*

_{calc},

*C*

_{calc}, and

*D*

_{calc}(the latter computed neglecting four-body interactions) were taken from Refs. 11, 127, and 126, respectively. The top panel of Fig. 5 shows that the differences trend downward as the densities increase above about 4000 mol/m

^{3}. This trend, as a function of (

*p*,

*T*), was noted in Ref. 126, together with the suggestion “there may have been a small error in the calibration for the sinkers….” However, the trend (Fig. 5, top) plotted as a function of density suggests that

*ρ*

_{meas}is sensitive to some of the truncated virial coefficients. The truncation suggestion is confirmed by the middle panel of Fig. 5, which includes in

*ρ*

_{calc}(

*p*,

*T*) the two additional terms

*E*

_{calc}(

*T*)

*ρ*

^{4}and

*F*

_{calc}(

*T*)

*ρ*

^{5}calculated semi-classically in Ref. 121. Additional terms [e.g.,

*G*

_{calc}(

*T*)

*ρ*

^{6}from Ref. 121] are less than 1.3 ppm, too small to be visible in Fig. 5.

The claimed *k* = 2 uncertainty of *ρ*_{meas} is 150 ppm;^{125} the span of the upper panels of Fig. 5 is ±150 ppm. The dashed curves (- -) in the middle panel of Fig. 5 represent upper bounds to the uncertainty of *ρ*_{calc}(*T*, *ρ*) at 223 K. For these upper bounds, we used the *k* = 2 uncertainties of the virial coefficients *U*(*B*), *U*(*C*), … provided by their authors. In Eq. (17) we replaced *B* with *B* + *U*(*B*); we replaced *C* with *C* + *U*(*C*), etc. The uncertainties of *ρ*_{calc}(*T*, *ρ*) are smaller at higher temperatures. We conclude that *ρ*_{calc} agrees with *ρ*_{meas} well within combined uncertainties.

At densities above $\u223c4000$ mol/m^{3}, the differences (*ρ*_{meas}/*ρ*_{calc} − 1) are nearly independent of the density; however, the average densities are 34 ppm larger than their expected values *ρ*_{calc}. These offsets are well within the claimed measurement uncertainties (*k* = 2, $\u223c150$ ppm). However, as shown in the lower panel of Fig. 5, the offsets have both a random and a systematic dependence on the temperature. The systematic temperature dependence can be treated as a correction to the calibration of the sinkers’ densities *ρ*_{sinker}(*p*, *T*). Such a correction does not remove the spread (±14 ppm) among the four isotherms at 273 K. Possible causes of this spread are changes between runs of temperature (±3.8 mK) and/or of impurity content (e.g., ±2.3 ppm of N_{2}). In any case, the offsets are smaller than the claimed uncertainties of *ρ*_{meas}(*p*, *T*).

Moldover and McLinden^{21} extended McLinden and Lösch-Will’s data^{125} to 500 K. The extended data are a less-stringent test of the VEOS than Fig. 5 because they span the same pressure range (*p* < 38 MPa) at higher temperatures; therefore, they span a smaller density range. If McLinden’s data could be extended to lower temperatures with comparable uncertainties, they would test helium’s VEOS in greater detail and they might reach a regime where *U*(*ρ*_{meas}) < *U*(*ρ*_{calc}). Schultz and Kofke conducted much more detailed tests of McLinden and Lösch-Will’s data.^{121} We agree with their conclusion that the data are consistent with the VEOS calculated *ab initio*.

It may be possible to significantly reduce the uncertainty of *ρ*_{meas} by improving magnetic suspension densimeters, as suggested by Kayukawa *et al.*^{128} They fabricated sinkers from single crystals of silicon and germanium because these materials have outstanding isotropy, stability, and well-known physical properties. Also, they refined the model and the functioning of their magnetic suspension so that it was independent of the magnetic properties of the fluid under study at the level of 1 ppm. They measured the density of a liquid near ambient temperature and pressure with a claimed *k* = 1 relative uncertainty of 5.4 × 10^{−6}. To date, they have not demonstrated this uncertainty far from ambient temperature and pressure. Even if *ρ*_{meas} achieved such low uncertainties, tests of the VEOS would have to solve problems arising from impure gas samples and imperfect temperature and pressure measurements.

Alternative methods of measuring equations of state have been reviewed by McLinden.^{129} Several methods require filling a container of known volume *V*_{cont}(*p*_{0}, *T*_{0}) with a known quantity of gas and then measuring the pressure as the temperature is changed. These methods resemble the CVGT method discussed in Sec. 2.2.4. Like CVGT, they require accurate values of *V*_{cont}(*p*, *T*); however, unlike CVGT, testing a VEOS requires much higher pressures. Determining *V*_{cont}(*p*, *T*) over wide ranges is complex because: (1) containers comprised of metal alloys have anisotropic elastic and thermal expansions; (2) containers have seals and joints or welds which have complicated mechanical properties; (3) alloys creep and/or anneal under thermal and mechanical stresses. In summary, volumetric methods are unlikely to replace Archimedes-type densimeters because *V*_{cont}(*p*, *T*) is an assembled object subjected to complicated stresses; in contrast, the densimeter’s sinkers are single objects subjected to hydrostatic pressure.

Remarkably, the Burnett method^{130} of measuring the equation of state requires neither determining *V*_{cont}(*p*, *T*) nor measuring quantities of gas. This method uses two pressure vessels with stable volumes *V*_{a} and *V*_{b}. On each isotherm, gas is admitted into *V*_{a} and the pressure is measured. The gas is allowed to expand so that it fills both *V*_{a} and *V*_{b} and the pressure is measured again. *V*_{b} is evacuated and the process is repeated several times. The measured pressures on each isotherm are fitted to the VEOS and an apparatus parameter: the volume ratio at zero pressure (*V*_{a,0} + *V*_{b,0})/*V*_{a,0}. The pressure dependences of *V*_{a} and *V*_{b} must also be known. Usually, they are estimated from elastic constants and models of the pressure vessels; therefore, precise estimates encounter complications of estimating *V*_{cont}(*p*, *T*). Perhaps this explains the large scatter in Burnett determinations of *D*(*T*).^{126} A fairly recent Burnett measurement of the equations of state of nitrogen and hydrogen (353–473 K; 1–100 MPa) claimed *k* = 2 uncertainties of *ρ*_{meas} ranging from 0.07% to 0.24%.^{131}

In addition to *ρ*_{meas}, measurements of the squared speed of sound *w*^{2}(*p*, *T*) in gases have been used to critically test either the VEOS^{122} of Eq. (7) or its acoustic analog, Eq. (1). Accurate values of *w*^{2}(*p*, *T*) in gases are readily available. At the low gas pressures used for acoustic thermometry, the relative expanded uncertainties *U*_{r}(*w*^{2}(*p*, *T*)) measured using quasi-spherical cavity resonators are a few parts in 10^{6} and are dominated by thermometry problems and/or impurities. However, uncertainties grow approximately linearly in pressure because of imperfect models of the recoil of the cavity’s walls in response to the resonating gas. In one study of argon, *U*_{r}(*w*^{2}) ≈ 1.2 × 10^{−4}(*p*/20 MPa) except near the critical point.^{132} At pressures above $\u223c5$ to $\u223c10$ MPa, pulse-echo techniques achieve uncertainties comparable to or smaller than resonance techniques.^{122,133} Remarkably, *w*^{2} from the two techniques agreed within 60–200 ppm within a range of overlap (argon, 250–400 K, $\u223c10$ to $\u223c20$ MPa^{133}).

It is more complex to compare $wmeas2(p,T)$ to a calculated VEOS than to compare *ρ*_{meas} to the same VEOS. To calculate the *n*-th acoustic virial coefficient from the *n*-th density virial coefficient, one also needs the first and second temperature derivatives of the *n*-th virial coefficient as well as all the lower-order density virial coefficients and their temperature derivatives. There are several routes to conduct such a comparison, which are completely equivalent. First, the temperature derivatives of the density virial coefficients can be calculated from *ab initio* potentials using, e.g., the Mayer sampling Monte Carlo (MC) method. Second, the temperature derivatives can be obtained from fits of the theoretically calculated temperature-dependent density virial coefficients. Third, the virial equation of state can be transformed by thermodynamic identities into an acoustic virial equation of state or it can be integrated to formulate a Helmholtz energy equation, from which the speed of sound can be calculated. Speeds of sound calculated by either of the two resulting equations contain contributions from terms with higher acoustic virial coefficients than those used in the density virial equation of state, i.e., it can be expected that the region of convergence of this virial equation of state for the speed of sound extends to higher pressures than that of the acoustic virial equation of state with virial coefficients derived directly from density virial coefficients. These terms describe contributions of configurations of particles which are contained in the low-order density virial coefficients to the higher-order acoustic virial coefficients. Fourth, densities can be calculated from $wmeas2(p,T)$ by the method of thermodynamic integration^{134} and directly compared to the density virial equation of state. As initial conditions for the integration, the density and heat capacity on an isobar must be known. There are subtleties to integrating $wmeas2(p,T)$.^{135} In the first method the uncertainties of the virial coefficients and their temperature derivatives follow from the MC simulation and can be propagated into an uncertainty of the acoustic virial equation of state, while in the other methods the uncertainty of the density virial coefficients or the experimental speeds of sound can be propagated into the acoustic virial equation of state or calculated densities, respectively.

For helium, Gokul *et al.*^{136} calculated the acoustic virial coefficients through the seventh order by the second method outlined above from density virial coefficients. They used the second density virial coefficients reported by Czachorowski *et al.*,^{11} which are based on the pair potential reported in the same work. The higher virial coefficients were taken from the work of Schultz and Kofke.^{121} They are based on the pair potential of Przybytek *et al.*^{137} and the three-body potential of Cencek *et al.*^{138} Uncertainties in the density virial coefficients were propagated into uncertainties in the acoustic virial coefficients by the MC method recommended in Supplement 1 to the “Guide to the Expression of Uncertainty in Measurement.”^{139} Gokul *et al.*^{136} formulated the acoustic virial equation of state as expansion in terms of density or pressure. The uncertainty of speeds of sound calculated with the acoustic VEOS was estimated from the uncertainty of the acoustic virial coefficients.

The density expansion of Gokul *et al.* was compared to the experimental data of Gammon,^{140} Kortbeek *et al.*,^{141} and Plumb and Cataland.^{142} The data of Gammon were measured with a variable-path interferometer operating at 0.5 MHz. They cover the temperature range between 98 and 423 K with pressures up to 15 MPa, and according to the author have an uncertainty of 0.003% of *w*^{2}. For these data, we estimated the expanded (*k* = 2) relative uncertainty *U*_{r}(*w*^{2}) = 0.000 09 by adding uncertainties of 0.003% (for the distance between the crystals), 0.001% (for the precision), and 0.005% (for sample impurities and/or temperature errors, based on the inconsistencies among the 14 isotherms). Gammon’s data agree with the acoustic virial equation of state within 0.01% with a few exceptions. The data of Kortbeek *et al.* were measured with a double-path-length pulse-echo technique, cover the temperature range from 98 to 298 K at pressures between 100 MPa and 1 GPa, and, according to the authors, have an uncertainty of 0.08%. They deviate from the acoustic VEOS between a few tenths of a percent at 100 MPa up to about 4% at 298 K and 1 GPa. These rather large deviations are due to the fact that the acoustic virial equation of state is not converged at such high pressures. The measurements of Plumb and Cataland cover the low temperature range between 2.3 and 20 K at pressures up to 150 kPa. They agree with the acoustic virial equation of state of Gokul *et al.* to within 0.05% except at the lowest measured pressures of about 1.5 kPa, where the deviations reach up to 0.18%. Gokul *et al.* also assessed the pressure range in which the acoustic VEOS is more accurate than the available experimental data for the speed of sound. At low pressures, they observed that speeds of sound calculated with the acoustic VEOS are more accurate than the experimental data of Gammon. Gokul *et al.* further noticed that speeds of sound calculated with the acoustic VEOS are more accurate than the experimental data of Kortbeek *et al.* up to about 300 MPa depending on temperature. At higher pressures, they considered the experimental data of Kortbeek *et al.* to be more accurate than the computed virial equation of state. This conclusion appears to be too optimistic in light of the low uncertainty of 0.08% in the experimental data and the rather large deviations of up to 2% from the virial equation of state below 300 MPa.

Gokul *et al.* also examined the convergence behavior of the acoustic virial equation of state more closely for the expansions in density and pressure. They considered a virial equation of state converged if the value of the speed of sound calculated with it agrees with all higher orders of the expansion within a certain tolerance. They observed that the pressure range in which the expansion in density converges is extended when the tolerance is increased. However, the expansion in pressure hits a pressure limit in the supercritical region, above which increasing the tolerance does not extend the region of convergence farther. Above this pressure limit, the expansion in pressure completely fails. Recently, Wedler and Trusler measured the speed of sound in supercritical helium with a dual-path pulse-echo technique in the temperature range between 273 and 373 K up to 100 MPa with an expanded uncertainty (*k* = 2) of 0.02%–0.04%.^{143} Their data agree with the seventh-order acoustic VEOS in density of Gokul *et al.*^{136} with a few exceptions in the whole range of the measurements within 0.025%, which shows that this form of the VEOS is converged in the region of the measurements.

The first calculation of the third virial coefficient of argon using a first-principles three-body potential was performed by Mas *et al.*^{144} using the empirical pair potential developed by Aziz.^{145} The results agreed almost to within combined uncertainties with the third virial coefficients extracted from experimental data (with theoretical constraints) by Dymond and Alder.^{146} Jäger *et al.* calculated density virial coefficients up to seventh order for argon with their pair and nonadditive three-body potentials.^{123} The calculated virial coefficients were fitted by polynomials in temperature. The seventh-order VEOS was compared with the very accurate (*p*, *ρ*, *T*) data of Gilgen *et al.*,^{147} which were measured with a magnetic suspension densimeter. These data are characterized by a relative uncertainty (*k* = 2) in density of 0.02%. Pressures calculated with the theoretical virial equation of state agree with these data at the highest temperature of the measurements, 340 K, within 0.01%.

In further work, Jäger^{148} used thermodynamic identities to calculate several properties of argon including the speed of sound from the virial equation of state and compared the results with the accurate experimental data of Estrada-Alexanders and Trusler^{132} and Meier and Kabelac.^{133} The data of Estrada-Alexanders and Trusler^{132} were measured with a spherical resonator and cover the temperature range between 110 and 450 K at pressures up to 19 MPa, while the data of Meier and Kabelac were measured with a dual-path-length pulse-echo technique and cover the temperature range between 200 and 420 K with pressures between 9 and 100 MPa. The expanded (*k* = 2) uncertainty of these datasets was estimated to be 0.001%–0.007% and 0.011%–0.036%, respectively. At 300 and 400 K, the calculated speeds of sound agree with both experimental datasets up to 100 MPa within 0.04% and 0.08%, respectively. At the near-critical temperature 146 K and supercritical temperature 250 K, the deviations of the calculated values from the experimental data of Ref. 132 increase with pressure from essentially zero in the ideal-gas limit to about 0.3% at 3.7 MPa and about 0.02% at 12.2 MPa.

In another paper, Jäger *et al.* presented calculations of the second and third density virial coefficient of krypton.^{149} They developed a very accurate pair potential for the krypton dimer, and nonadditive three-body interactions were described by an *ab initio* extended Axilrod–Teller–Muto potential, which was fitted to quantum chemical calculations of the interaction energy of equilateral triangle configurations of three krypton atoms. El Hawary *et al.*^{122} calculated density virial coefficients from the fourth to the eighth using the pair potential and extended Axilrod–Teller–Muto potential of Jäger *et al.* The calculated virial coefficients were fitted to polynomials in temperature, and the virial equation of state was integrated to formulate it as a fundamental equation of state in terms of the Helmholtz energy. Furthermore, El Hawary *et al.* measured the speed of sound in liquid and supercritical krypton between 200 and 420 K at pressures from 6.1 to 100 MPa with an uncertainty (*k* = 2) of 0.005%–0.018%. At 240, 320, and 420 K, the seventh-order and eighth-order virial equations of state agree with each other within 0.02% up to 7, 17, and 38 MPa, respectively. In the region where the virial equation of state is sufficiently converged, the calculated speeds of sound are systematically about 0.08% lower than the experimental data. This small difference is probably due to the uncertainty of the pair potential and the simplified treatment of nonadditive three-body interactions with the extended Axilrod–Teller–Muto model.

At high density in the supercritical region where the virial equation of state does not converge and in the liquid region, thermodynamic properties can be calculated by MC or molecular-dynamics (MD) simulations.^{150} Since the generation of Markov chains in MC simulations avoids some of the numerical errors of algorithms used to integrate the equations of motion in MD simulations, MC simulations are the preferred method for calculating accurate values for thermodynamic properties. In statistical mechanics, there are eight basic ensembles for performing MC or MD simulations,^{151} which are characterized by a thermodynamic potential, three independent variables, and a weight factor, which describes the distribution of systems in the ensemble, in which MC simulations of fluids can be performed. Ströker *et al.*^{152} pointed out that the *NpT* ensemble, in which the number of particles, the pressure, and the temperature are the independent variables, is best suited for the calculation of thermodynamic properties because only ensemble averages involving the enthalpy and volume, but no derivatives of the potential energy with respect to volume, appear in the equations for thermodynamic properties. This means that no derivatives of the potentials are needed in a simulation.

The argon calculations of Mas *et al.*^{144} described above were later extended by performing *NVT*, *NpT*, and Gibbs ensemble MC simulations^{153} along the vapor–liquid coexistence curve. The parameters of the critical point agreed with experiments within 0.8% or better.^{154}

Ströker *et al.*^{152} carried out semiclassical MC simulations of thermodynamic properties of argon in the *NpT* ensemble at the subcritical isotherm 100 K and the supercritical isotherm 300 K at pressures up to 100 MPa. The interactions between argon atoms were described by the pair potential of Jäger *et al.*^{155} and the nonadditive three-body potential of Jäger *et al.*^{123} Quantum effects were accounted for by the Feynman–Hibbs corrections to the pair potential. Calculated densities agree with the accurate data measured by Gilgen *et al.*^{147} and Klimeck *et al.*^{156} within less than 0.01%, while calculated speeds of sound agree within less than 0.1% with the accurate experimental data of Estrada-Alexanders and Trusler^{132} at low pressure in the supercritical region and Meier and Kabelac^{133} at high pressure in the liquid and supercritical region.

Ströker *et al.*^{157} also performed MC simulations for liquid and supercritical krypton. They employed the accurate pair potential and an extended Axilrod–Teller–Muto potential of Jäger *et al.*^{149} to account for nonadditive three-body interactions. Quantum effects were again accounted for semiclassically. Since the potential models for krypton are not as accurate as those for argon, the deviations of the results for the density and speed of sound from experimental data were larger than for argon, about 0.2% and 0.36%, respectively.

### 2.5. Transport properties and flow metrology

In this section, we describe the impact of the *ab initio* calculations of the zero-density limit of helium’s thermal conductivity *λ*_{He} and viscosity *η*_{He}. First, we mention the impact of *λ*_{He} and *λ*_{Ar} on temperature metrology. Then, we describe how accurate values of *η*_{He} have been used as standards to reduce the uncertainty of viscosity measurements of many gases by a factor of 10. We conclude by briefly considering the impact of accurate viscosity data on metering process gases, for example, during the manufacture of semiconductor chips.

As discussed in Sec. 2.2.1, AGT requires accurate values of *λ* of the working gas at low densities to account for the effect of the thermo-acoustic boundary layer on the measured resonance frequencies. For example, in 2010, Gavioso *et al.* used helium at $\u223c410$ kPa in a single-state, AGT determination of the Boltzmann constant *k*_{B} prior to its definition in 2019.^{158} They reported that a relative standard uncertainty *u*_{r}(*λ*_{He}) = 0.015 generated a relative standard uncertainty of the Boltzmann constant *u*_{r}(*k*_{B}) = (1–3) × 10^{−6}.

Today, an uncertainty of (1–3) × 10^{−6} would be the largest contributor to a state-of-the-art determination of the thermodynamic temperature *T* near 273 K. At low temperatures, the uncertainty of measured values of *λ*_{He} is much larger. Below 20 K, the *λ*_{He} data span a range on the order of ±6%.^{159} This large an uncertainty would lead to *u*_{r}(*T*) > 10^{−5} for acoustic determinations of *T*. Fortunately, the values of *λ*_{He} calculated *ab initio* have extraordinarily small uncertainties, e.g., *u*_{r}(*λ*_{He}) = 9.6 × 10^{−6} at 273 K and *u*_{r}(*λ*_{He}) = 7.3 × 10^{−5} at 10 K.^{10} In essence, the calculated values of *λ*_{He} removed *u*_{r}(*λ*_{He}) from the uncertainty budgets of acoustic thermometers based on helium-filled quasi-spherical cavities.

Cylindrical, argon-filled cavities are being developed for high-temperature acoustic thermometry.^{4,55} These projects require low-uncertainty values of both *λ*_{Ar} and *η*_{Ar}. Low-uncertainty values of *η*_{Ar} were generated from accurate measurements of the ratios *η*_{Ar}/*η*_{He} in the range 200–653 K and the *ab initio* values of *η*_{He}. Then *λ*_{Ar}(*T*) was obtained by combining the ratio-deduced values of *η*_{Ar}(*T*) with values of the Prandtl number Pr_{Ar} calculated from model pair potentials. (Pr = *C*_{p}*η*/*λ*, where *C*_{p} is the constant-pressure heat capacity per mass. For the noble gases, Pr is only weakly sensitive to the potential.)^{160–162} The measured ratios *η*_{Ar}/*η*_{He} were consistent, within a few tenths of a percent, with highly accurate measurements made with an oscillating-disk viscometer^{163} and with calculations of *η*_{Ar} based on *ab initio* Ar–Ar potentials.^{164} Thus, the needs of argon-based acoustic thermometry are now met at all useable temperatures. To put this achievement in context, we note that measuring the thermal conductivity of dilute gases is difficult, even for noble gases near ambient temperature and pressure. Evidence for this appears in Lemmon and Jacobsen’s correlation of the “best” measurements of *λ*_{Ar} and *η*_{Ar} near ambient temperature (270–370 K) and pressure.^{165} The average absolute deviations of selected measurements from their correlation ranged from 0.24% to 1.0%. Lemmon and Jacobsen estimated the uncertainty of the correlated values of *λ*_{Ar} was 2% and the uncertainty of *η*_{Ar} was 0.5%. (With the benefit of *ab initio* calculations and ratio measurements, we now know their correlation overestimated *λ*_{Ar} by 0.54% at 270 K and by 0.45% at 370 K.)

In 2012, Berg and Moldover reviewed measurements of the viscosity of 11 dilute gases near 25 °C.^{166} These measurements were made using 18 different instruments that used five different operating principles and produced 235 independent viscosity ratios during the years 1959–2012. Using the *ab initio* value of *η*_{He} at 25 °C as a reference, the viscosities of the ten other gases (Ne, Ar, Kr, Xe, H_{2}, N_{2}, CH_{4}, C_{2}H_{6}, C_{3}H_{8}, SF_{6}) were determined with low uncertainties *u*_{r}(*η*) ranging from 0.000 27 to 0.000 36. These ratio-derived uncertainties are less than 1/10 the uncertainties claimed for absolute viscosity measurements, such as the measurements of *η*_{Ar} correlated by Lemmon and Jacobsen.^{165} Now, any one of these gases can be used to calibrate a viscometer within these uncertainties. Such ratio-based calibrations have reduced uncertainties of *η* for many other gases^{167} and they have been extended to a wide range of temperatures.^{168} During their study of viscosity ratios, Berg and Moldover observed that the viscosity ratios determined using one instrument (a magnetically suspended, rotating cylinder) were anomalous. Their observation led to an improved theory of the instrument, thereby illustrating the power of combining a reliable standard *η*_{He}(*T*) with precise ratio measurements.^{169}

_{3})

_{3}, WF

_{6}] under conditions differing from the calibration conditions. An accurate transition between gases and conditions can be made using laminar flow meters for which there is a physical model (similar to the model of a capillary tube). Also needed are data for the virial coefficients of the process gas and the viscosity ratio

^{170}

_{6}.

The initial density expansion of the viscosity has the form *η*/*η*_{0} = 1 + *η*_{1}*ρ*, where the low-density limit of the viscosity *η*_{0} depends entirely on pair interactions and the virial-like coefficient *η*_{1} depends on the interactions among two and three molecules. Unfortunately, unlike the density and dielectric virial coefficients and *η*_{0}, no rigorous theory exists for *η*_{1}(*T*). An approximate theory was developed by Rainwater and Friend,^{171,172} who presented quantitative results based on the Lennard-Jones potential. It was later extended with more accurate pair potentials for noble gases.^{173} While the results from the Rainwater–Friend model are in reasonable agreement with the limited experimental data available for the initial density dependence of the viscosity for noble gases,^{173} the error introduced by its approximations is unknown. We note that it is a classical theory, which introduces another source of error for light gases (such as helium) where quantum effects might be important, even at ambient temperatures.

## 3. *Ab Initio* Electronic Structure Calculations

### 3.1. Methodology of electronic structure calculations

In principle, solutions of the equations of relativistic quantum mechanics, possibly including quantum electrodynamics (QED) corrections, can predict all properties of matter to a precision sufficient for thermal metrology applications. In practice, if the goal is to match or exceed the accuracy of experiments, the range of systems reduces to few-particle ones. The first quantum mechanical calculations challenging experimental measurements for molecules appeared only in the 1960s (e.g., Ref. 174), while the first calculations relevant to metrology were published in the mid-1990s.^{7,175,176} Currently, the branches of metrology discussed in this review are becoming increasingly dependent on theoretical input, as discussed in Sec. 2.

Theory improvements leading to results with decreased uncertainties proceed along three main, essentially orthogonal directions: level of physics, truncation of many-electron expansions, and basis set size. There exists an extended hierarchy of approaches in each direction. For the first direction, there exists a set of progressively more accurate physical theories that can be used in calculations relevant for metrology, from Schrödinger’s quantum mechanics for electrons’ motion in the field of nuclei fixed in space to relativistic quantum mechanics and to QED. The second direction is relevant for any many-electron system: one has to choose a truncation of the expansion of the many-electron wave function in terms of virtual-excitation operators at the double, triple, quadruple, etc. level or, equivalently in methods that use explicitly correlated bases (depending explicitly on interelectronic distances), to take into account only correlations of two, three, four, etc. electrons simultaneously. Third, for any given theory and many-electron expansion level, there are several methods of solving quantum equations specific for this level; in particular different types of basis sets are used to expand wave functions, resulting in different magnitudes of uncertainties from such calculations.

The lowest theory level is Schrödinger’s quantum mechanics for electrons moving in the field of nuclei fixed in space, i.e., quantum mechanics in the Born–Oppenheimer (BO) approximation. At the next level, one usually first accounts for the relativistic effects. Post-BO treatment of the Schrödinger equation can be limited to computations of the so-called diagonal adiabatic correction, which is the simplest method of accounting for couplings of electronic and nuclear motions, or it can fully include nonadiabatic effects, i.e., account for the complete couplings of these two types of motion. The highest level of theory applied in calculations relevant to metrology is QED, and it can be implemented at several approximations labeled by powers of the fine-structure constant *α*.

*N*-electron system is represented as a linear combination of Slater determinants constructed from “excitations” of the ground-state HF determinant |Φ

_{0}⟩

*N*-tuple excitations, where $|\Phi ar\u232a$ represents a singly excited Slater determinant formed by replacing spinorbital

*ϕ*

_{a}with

*ϕ*

_{r}. Similarly, $|\Phi abrs\u232a$ represents a doubly excited Slater determinant formed by replacing spinorbital

*ϕ*

_{a}with

*ϕ*

_{r}and spinorbital

*ϕ*

_{b}with

*ϕ*

_{s}, and so on for higher excited determinants. The linear coefficients (CI amplitudes) are computed using the Rayleigh–Ritz variational principle. While the FCI method is conceptually straightforward, the computation time it requires scales with the number of electrons as

*N*!, and therefore it is computationally very costly. One can limit the expansion in Eq. (18) to a subset of excitations (for example, retaining only single and double excitations leads to a method denoted CISD), but truncated expansions are not size extensive. This means that the CISD energy computed for very large separations between two monomers (atoms or molecules) is not equal to the sum of monomers’ energies computed at the CISD level. Only FCI is free of this problem. Thus, truncated CI expansions are not appropriate for calculations of interaction energies.

_{2}, a four-electron system, the expansion can be written as

^{177}

*c*

_{k}are variational parameters, and

*ϕ*

_{k},

*k*> 0 are ECG basis functions. The function

*ϕ*

_{0}is the product of ECG functions for the two helium atoms. The explicit form of

*ϕ*

_{k},

*k*> 0, functions is

*α*

_{ki},

*β*

_{kij}, and

*A*_{ki}= (

*X*

_{ki},

*Y*

_{ki},

*Z*

_{ki}) are nonlinear variational parameters. For a given set of nonlinear parameters, the linear parameters are obtained using the Rayleigh–Ritz variational method. The simplest way to optimize the nonlinear ones is to use the steepest-descent method, recalculating the linear parameters in each step of this method. In actual applications, significantly more advanced optimization methods are used.

*T*

_{i}can be written in terms of pairs of creation $r\u0302\u2020$, $s\u0302\u2020$ and annihilation $a\u0302$, $b\u0302$ operators replacing the occupied spinorbitals

*ϕ*

_{a}and

*ϕ*

_{b}by the virtual ones

*ϕ*

_{r}and

*ϕ*

_{s}. For the two lowest ranks, we have

*t*are different from the amplitudes

*c*. The former amplitudes are obtained by using the expansion (21) in the Schrödinger equation and projecting this equation with subsequent determinants from Eq. (18). Since the resulting set of equations is nonlinear, the solution is obtained in an iterative way. If all the excitation operators are kept in Eq. (22), the method is equivalent to the FCI method, but this expansion is almost always truncated. The simplest CC approach is that of CC doubles (CCD), in which $T\u0302$ is truncated to

*t*of single and double excitations in Eqs. (23) and (24) are computed iteratively while those for triple excitations are evaluated using perturbation theory. A similar approximation, denoted CCSDT(Q), can be made for the CCSDTQ method. In contrast to the truncated CI expansions, the truncated CC expansions are always size extensive. This results from the fact that the exponential ansatz of Eq. (21) can be factored for large separations between subsystems into a product of exponential operators for subsystems. The CC method is applied to interatomic or intermolecular interactions in the supermolecular fashion, i.e., subtracting monomers’ total energies from the total energy of a cluster. Due to size extensivity, the resulting potential-energy surface (PES) dissociates correctly.

^{178–181}The basic assumption of SAPT is the partitioning of the total Hamiltonian

*H*of a cluster into the sum of the Hamiltonians of separated monomers

*V*that collects Coulomb interactions of the electrons and nuclei of a given monomer with those of the other monomers:

*H*

_{0}

_{0}does not satisfy Pauli’s exclusion principle. For large intermonomer separations

*R*, one can ignore this problem and use the Rayleigh–Schrödinger perturbation theory (RSPT), the simplest form of intermolecular perturbation theory. Unfortunately, RSPT leads to unphysical behavior of the interaction energy at short

*R*as it fails to predict the existence of the repulsive walls on potential-energy surfaces. This failure is the result of the lack of correct symmetry of the wave function under exchanges of electrons between interacting monomers. Thus, to describe interactions everywhere in the intermonomer configuration space, one has to perform symmetry adaptation, i.e., antisymmetrization, and this is the origin of the phrase “symmetry-adapted.” There are several ways to do it, but the simplest is to (anti)symmetrize the wave functions of the RSPT method. This leads to the symmetrized Rayleigh–Schrödinger (SRS) approach,

^{182}which is the only SAPT method used in practice. For a dimer, the interaction energy is then expressed as the following series in powers of

*V*:

*V*(orders of perturbation theory) and different terms of the same order can be identified as resulting from different physical interactions: electrostatic (elst), exchange (exch), induction (ind), and dispersion (disp). When SAPT is applied to many-electron systems, monomers can be described at various levels of electronic structure theory: from the HF level to the FCI level. This leads to a hierarchy of SAPT levels of approximations depending on treatment of intramonomer electron correlation. If the monomers are approximated at an order

*n*of many-body perturbation theory (MBPT) with the Møller–Plesset (MP) partition of the Hamiltonians

*H*

_{A}and

*H*

_{B}, denoted as MP

*n*, we can write

*n*→ ∞.

The third direction determining the accuracy of electronic structure calculations involves the size of the basis sets used to expand wave functions. In the CC and CI approaches, the standard technique is to use products of orbital (one-electron) basis sets. Many such basis sets are available; the ones most often used in metrology-related calculations are the correlation consistent (cc) basis sets introduced by Dunning.^{183} These basis sets are denoted by cc-pV*X*Z: cc polarized Valence (i.e., optimized using a frozen-core approximation), *X*-Zeta, where *X* = D, T, Q, 5, … is the so-called cardinal number determining the maximum angular momentum of orbitals. Such basis sets can be augmented by an additional set of diffused functions and are then denoted as aug-cc-pV*X*Z, or two such sets: daug-cc-pV*X*Z. Another option is to use explicitly correlated basis sets in the CC method or to expand the whole many-electron wave function in such a basis set. Explicitly correlated basis sets provide a much faster convergence than products of orbital basis sets, but in most cases require optimizations of a large number of nonlinear parameters.

In order to achieve some target size of uncertainties, one has to choose a proper level in each of the three hierarchies defined earlier. For example, it is possible to perform an FCI calculation for a ten-electron system such as Ne. However, since FCI calculations scale factorially with the number of orbitals, only very small basis sets can be used, resulting in a large uncertainty of the results. Consequently, a better strategy is to use the CCSD(T) method which allows applications of the largest orbital bases available for a system like Ne_{2}. The computed interaction energy will be accurate to about four significant digits relative to the CCSD(T) limit, but will have a fairly large error, of the order of 1%–2%, with respect to the exact interaction energy at the non-relativistic BO level. In contrast, FCI calculations for Ne_{2} employing the smallest sensible basis set of augmented double-zeta size, which would be extremely difficult to perform, would have an error of the order of 40% (such calculations might still be useful in hybrid approaches discussed below).

The orbital basis sets consist of families of bases of varying size. One usually carries out calculations in two or more such basis sets and then performs approximate extrapolations to the complete basis set (CBS) limit. In addition to the standard extrapolations, which assume the *X*^{−3} decay of errors, extrapolations using very accurate ECG results can be performed.^{184–186} CCSD(T)/CBS results may have sufficiently small uncertainties to make calculations of relativistic and diagonal adiabatic corrections necessary, i.e., these corrections may be of the same order of magnitude as the uncertainties of the CCSD(T)/CBS results. To reduce the errors resulting from the truncation of the many-electron expansion, one can follow CCSD(T)/CBS calculations by CCSDT(Q) or FCI ones in smaller basis sets. These effects are then included in an incremental way, i.e., by adding the difference between FCI and CCSD(T) energies computed in the same (small) basis set.

Accurate solutions of quantum equations are followed by estimates of uncertainties, absolutely necessary for metrology purposes. The latter step is often more time-consuming than the former. One should emphasize that theoretical estimates of uncertainties are different from statistical estimates of uncertainties of measurements and in particular one cannot assign a rigorous confidence level to them, although for purposes of metrology one usually assumes that theoretical uncertainties are equivalent to *k* = 2 expanded uncertainties (95% confidence level).

A theoretical estimate of uncertainty consists of several elements. The most rigorous and reliable estimates are those of basis set truncation errors derived from the observed patterns of convergence in basis set. Much more difficult are estimates of uncertainties resulting from truncations of many-electron expansions. Such estimates can sometimes be made by performing higher-level calculations at a single point on a PES, but one most often uses analogy to similar systems for which higher-level calculations have been performed. The same approach can be used to estimate the neglected physical effects, for example, to estimate the uncertainty due to relativistic effects.

Solutions of the electronic Schrödinger equation for a given nuclear configuration of a dimer or a larger cluster, providing accurate quantum mechanical descriptions of such systems, are only the first step in theoretical work of relevance to metrology, as most measured quantities discussed in this review are either bulk properties or response properties of atoms and molecules. In the former case, i.e., to predict properties of gases or liquids relevant for metrology, one needs to know energies of such systems for a large number configurations, i.e., for different geometries of clusters. This issue is approached by using the many-body expansion, where here the bodies are atoms or molecules forming the cluster, starting from two-body (pair) interactions, followed by three-body (pairwise nonadditive) interactions. The approach can be continued to higher-level many-body interactions, but so far this has not been done. The *ab initio* energies are usually fitted to analytic forms only for the two- and three-body interactions.

In addition to energies, metrology applications often require knowledge of accurate values of various properties of atoms and molecules, mainly the static and dynamic polarizabilities and magnetic susceptibilities. These quantities can be computed as analytic energy derivatives with respect to appropriate perturbations. Properties of a single atom or molecule change in condensed phases and the so-called interaction-induced corrections to properties of isolated atoms or molecules are of interest to metrology.

As already mentioned above, although Schrödinger’s quantum mechanics at the BO level provides the bulk of the physical values of interest to metrology, computations of various effects beyond this level are often needed to reduce uncertainties of these properties to the magnitude needed for metrology standards. We will refer to these as post-BO effects. It should be stressed that we really have in mind here the post-nonrelativistic-BO level, since both the relativistic and QED corrections for molecules are usually computed using the BO approximation. One goes beyond this approximation when computing adiabatic and nonadiabatic corrections. Any reasonably detailed description of methodologies used in post-BO calculations would be too voluminous for the present review. Therefore, we refer the reader to the original papers, in particular Refs. 10, 11, 84, 137, 177, and 187–192.

Systems of interest to thermodynamics-based precision metrology are mainly noble-gas atoms and their clusters, and this section will be restricted to such systems, with the majority of text devoted to helium. Apart from being the substance whose behavior is closest to the ideal gas, it is also the only system where theory can currently provide values of physical quantities that are generally more accurate than the measured ones. Nevertheless, neon and argon are also of significant interest since they may be used in secondary standards to improve instrument sensitivity or ease of use. Although for most properties computations for neon and argon have larger uncertainties than the best measurements, such results are still useful as independent checks of experimental work and to guide extrapolation beyond the measured range.

#### 3.1.1. Importance of explicitly correlated basis sets

The current theoretical results for helium owe their very small uncertainties mostly to the use of explicitly correlated basis sets. The calculations involving helium atoms are probably one of the best examples where an important science problem was solved using these bases. To clearly show where this field would be without the use of such basis sets, we discuss in this subsection numerical comparisons of ECG and orbital calculations for He_{2}, performed recently in Ref. 193. The majority of molecular electronic structure calculations are carried out using orbital basis sets. This means that many-electron wave functions are expanded in products of orbitals. The simplest example is the CI method discussed earlier, where the wave function is a linear combination of Slater determinants built of orbitals that are usually obtained by solutions of HF equations. However, expansions in orbital products converge slowly due to the difficulty of reproducing the electron cusps in wave functions.

A way around this difficulty is to use bases that depend explicitly on *r*_{12} = |*r*_{2} − *r*_{1}|, the distance between electrons. Bases of this type are called explicitly correlated. For few-electron systems, such bases are mostly used to directly expand the *N*-electron wave functions of the nonrelativistic BO approximation. The explicitly correlated bases are also often used for many-electron systems within a perturbative or CC approach.^{194,195} For two-electron systems, one uses Hylleraas bases^{196} with polynomial-only dependence on *r*_{12} or, recently more frequently, expansions in purely exponential functions of *r*_{12}, called Slater geminals.^{84,192} Bases combining both types of dependence on *r*_{12} are also sometimes employed.^{197} For more than two electrons, integrals needed for such bases become very expensive and bases involving Gaussian correlation factors $e\u2212\gamma rij2$, i.e., ECG bases, are mostly used. For a review of the ECG approach, see Refs. 194 and 198.

Since expansions in explicitly correlated bases of the type described above approach solutions of the Schrödinger equation, the equivalent orbital calculations should be performed at the FCI level. As already mentioned, FCI calculations scale as *N*! with the number of electrons and therefore are the most expensive of all orbital calculations. Even for He_{2}, FCI calculations cannot be performed using the largest available orbital basis sets. Therefore, the optimal orbital-based strategy is a hybrid one consisting of performing calculations in the largest basis sets at a lower level of theory, for example, at the CCSD(T) level, and adding to these results FCI corrections computed in smaller basis sets.

The BO energies computed in ECG basis sets in Ref. 177 established a new accuracy benchmark for the helium dimer; see the description of these calculations in Sec. 3.3.1. These ECG interaction energies were compared in Ref. 193 to those computed in orbital bases at the hybrid CCSD(T) plus FCI level. The largest available basis sets were applied. For most points, the CCSD(T)+ΔFCI approach gives errors nearly two orders of magnitude larger than the ECG estimated uncertainties. For a couple of points, the CCSD(T)+ΔFCI results are fairly close to the ECG results, but this is mainly due to the former method overestimating the magnitude of the interaction energy at small *R* and underestimating at large *R*. Since these points are near the van der Waals minimum, some previous evaluations of the performance of orbital methods restricted to this region might have been too optimistic. When the whole range of *R* is considered, CCSD(T)+ΔFCI is no match for the ECG approach. One should also realize that any improvements of accuracy of the CCSD(T)+ΔFCI approach would require a huge effort; in particular, one would have to develop quadruple-precision versions of all needed orbital electronic structure codes.

### 3.2. Helium atom polarizability

One of the properties of helium required by precision measurement standards^{199–201} is the helium atom polarizability, both static and dynamic (frequency-dependent). Nonrelativistic calculations of the static polarizability date back to the 1930s and reached an accuracy of 0.1 ppb in 1996 calculations using Hylleraas basis sets.^{202} However, the relativistic correction, which is proportional to *α*^{2}, could be expected to contribute at the 60 ppm level relative to the total polarizability. Unfortunately, the values of these corrections published before 2001 differed significantly from one another. These discrepancies were resolved by accurate calculations of Refs. 187 (using GTGs) and 203 (using Slater geminals) with uncertainties of 20 ppb relative to the total polarizability. This work used the Breit–Pauli operator,^{204} whose expectation values were computed with the ground-state wave function for the nonrelativistic Hamiltonian.

The authors of Ref. 203 also computed the QED correction of order *α*^{3}, which turned out to be significant, amounting to about 20 ppm relative to the total polarizability. However, a part of the *α*^{3} QED correction to the polarizability, resulting from the so-called Bethe logarithm, was only roughly estimated due to the very difficult to compute second electric-field derivative of a second-order-type perturbation theory expression involving the logarithm of the Hamiltonian. The first complete calculation of Bethe-logarithm contribution to helium polarizability was reported in Ref. 188. As such a calculation had never been done before, the algorithms and their numerical implementations had to be developed from scratch. The term containing the electric-field derivative of Bethe’s logarithm turned out to be unexpectedly small, representing only about 0.6% of the total *α*^{3} QED correction. Thus, this correction still makes a contribution of about 20 ppm to the total polarizability.

Further improvement of the accuracy of helium’s static polarizability was achieved in Ref. 192, which concentrated on the second derivative of the Bethe logarithm with respect to the electric field. This quantity can be obtained in a couple of ways, with completely different algorithms. The goal was to achieve agreement between two such approaches and also with Ref. 188. This goal was met, providing a reliable cross-validation for both approaches. The results of Ref. 192, providing currently the most accurate theoretical determination of the polarizability of helium, are shown in Table 1. The value of the *α*^{3} QED contribution computed in Ref. 188 differs from the current one by only 0.03% or 7 ppb relative to the total polarizability. This error is much smaller than the current uncertainty of the *α*^{4} QED contribution, which is estimated to amount to 0.1 ppm; see Table 1.

Contribution . | Reference 192 . |
---|---|

Nonrelativistic | 1.383 809 986 4 |

α^{2} relativistic | −0.000 080 359 9 |

α^{2}/m relativistic recoil | −0.000 000 093 5(1) |

α^{3} QED $\u2212\u2202\u03f52lnk0$ term | 0.000 030 473 8 |

$\u2202\u03f52lnk0$ term | 0.000 000 182 2 |

α^{3}/m QED recoil | 0.000 000 011 12(1) |

α^{4} QED | 0.000 000 56(14) |

Finite nuclear size | 0.000 000 021 7(1) |

Total | 1.383 760 78(14) |

Contribution . | Reference 192 . |
---|---|

Nonrelativistic | 1.383 809 986 4 |

α^{2} relativistic | −0.000 080 359 9 |

α^{2}/m relativistic recoil | −0.000 000 093 5(1) |

α^{3} QED $\u2212\u2202\u03f52lnk0$ term | 0.000 030 473 8 |

$\u2202\u03f52lnk0$ term | 0.000 000 182 2 |

α^{3}/m QED recoil | 0.000 000 011 12(1) |

α^{4} QED | 0.000 000 56(14) |

Finite nuclear size | 0.000 000 021 7(1) |

Total | 1.383 760 78(14) |

Calculations of Refs. 187 and 188 were extended to frequency-dependent polarizabilities.^{84,191} This polarizability was expanded in inverse powers of the wavelength *λ* up to *λ*^{−8}. Different levels of theory were used for each power of *λ*: up to *α*^{4} for the static term, *α*^{2} for inverse powers 2 through 6 (only even powers contribute), and nonrelativistic for 8. The dynamic polarizability at the He–Ne laser wavelength of 632.9908 nm had an uncertainty of 0.1 ppm. This uncertainty results entirely from the uncertainty of the static polarizability. The latter was reduced compared to Ref. 188 mainly because the work of Ref. 205 has shown that the error of the so-called one-loop approximation used to evaluate the *α*^{4} terms is smaller than previously expected, amounting to only about 5% when applied to the excitation energies of helium. Another small change in the static polarizability was due to a slightly improved value of the Bethe-logarithm contribution; see also Ref. 192.

The polarizabilities computed in Refs. 84, 187, 188, and 191 had uncertainties orders of magnitude smaller than the best experimental results. However, recently a new, very accurate measurement of this quantity was published.^{63} The measured value of the molar polarizability, 0.517 254 4(10) cm^{3}/mol, is consistent with the theoretical molar polarizability computed from the atomic one listed in Table 1 and equal to 0.517 254 08(5) cm^{3}/mol. The combined uncertainty is more than three times the difference while the experimental uncertainty is 20 times larger than the theoretical one.

When a helium atom is in a gas or condensed phase, one can expect that its polarizability changes due to interactions with other atoms. More precisely, the polarizability of a helium cluster is not equal to the sum of polarizabilities of helium atoms. This change is often referred to as collision-induced polarizability and for atoms is a function of interatomic distance. For a pair of helium atoms, reliable values of this quantity were computed in Ref. 206, reconciling previously published inconsistent calculations. The results of Ref. 206 were used to compute the second^{79} and third^{207} dielectric virial coefficients of helium. Very recently, the collision-induced three-body polarizability of helium was computed.^{208}

Due to inversion symmetry, a system consisting of one or two helium atoms cannot have a dipole moment in the BO approximation. However, configurations of three or more atoms may have a non-zero dipole moment, which in turn influences the value of the third dielectric virial coefficient.^{207} Presently, the only *ab initio* description of the three-body dipole moment of noble gases is the one developed by Li and Hunt.^{209} However, the results of Ref. 209 apply only at large separations, and do not have associated uncertainties. A dipole-moment surface for the helium trimer with rigorously defined uncertainty is currently being developed.^{210}

### 3.3. Helium dimer potential

#### 3.3.1. BO level

The interest in the helium dimer potential is nearly as old as quantum mechanics. In 1928, Slater^{211} developed the first potential for this system, which gave the interaction energy of −8.8 K at the internuclear distance *R* = 5.6 bohrs (1 bohr $\u224852.91772109$ pm). There is a wide range of helium dimer potentials available in the literature; see Ref. 212 for a comparison of bound-state calculations using a large number of potentials. Figure 6 illustrates the remarkable progress in accuracy of predictions achieved since 1979. Empirical potentials dominated the field until the end of the 1980s; the two most widely used ones, HFDHE2^{213} and HFD-B,^{214} were developed by Aziz *et al.* The first really successful *ab initio* one was the LM-2 potential (published only in a tabular form) developed by Liu and McLean.^{215} Those authors performed CI calculations and, by analyzing the configuration space and basis set convergence, obtained extrapolated interaction energies with estimated uncertainties. Although these estimates were rather crude and do not embrace the current best values for most values of *R*, cf. Fig. 6, they are reasonable.

Aziz and Slaman^{216} used the HFD-B functional form with refitted parameters to “mimic” the behavior of the LM-2 potentials, of the unpublished *ab initio* data computed by Vos *et al.*,^{217} and of the small-*R* Green's-function MC (GFMC) data^{218} to obtain potentials denoted as LM2M1 and LM2M2, differing by assuming, respectively, the smallest and the largest well depth of the LM-2 potential as determined by the estimates of uncertainty. The parameters of these potentials were not fitted directly to *ab initio* data, but chosen by trial and error to reproduce both theoretical data and measured quantities to within their error bars. The LM2M2 potential was considered to be the best helium potential until the mid-1990s, when purely *ab initio* calculations took the lead. Among the latter ones, the TTY potential developed by Tang *et al.*^{219} has a remarkably simple analytical form based on perturbation theory. The HFD-B3-FCI1 potential was obtained by Aziz *et al.*,^{7} who used the HFD-B functional form with its original parameters adjusted so that the new potential runs nearly through the *ab initio* data points. These points were GFMC results of Ref. 218 and the FCI results of van Mourik and van Lenthe.^{220} No uncertainties were assigned to HFD-B3-FCI1, and Fig. 6 shows that it was about as accurate as LM2M2.

The SAPT96 potential^{175,176} opens an era of helium potentials based mostly on calculations with explicitly correlated functions. It was the first fully first-principles He_{2} potential with a systematic estimation of uncertainties. The potential was obtained using a two-level incremental strategy. The leading SAPT corrections (the complete first-order and the bulk of the second-order interaction energies) were computed using GTG basis sets. The GTG-based variant of SAPT was developed in Refs. 221–225. Higher-order SAPT corrections were computed using the general SAPT program based on orbital expansions.^{226–230} Large orbital basis sets including up to *g*-symmetry functions and midbond functions (placed between the nuclei)^{231} were used. The remaining many-electron effects were computed using both SAPT based on FCI-level monomers, with summations to a very high order of perturbation theory (using He_{2}-specific codes), and supermolecular FCI calculations in small orbital basis sets. It is interesting to note that the actual errors of the SAPT96 potential relative to the current best results turned out to be completely dominated by the residual orbital (rather than GTG) contributions. For instance, at *R* = 5.6 bohrs, the orbital part constitutes only −1.81 K out of −11.00 K, but its error was −0.05 K out of the total SAPT96 error of −0.06 K. The factor of 2 underestimation of the uncertainties seen in Fig. 6 for *R* = 5.6 bohrs was entirely due to this issue. SAPT96 is about as accurate as LM2, except for large *R* where it is more accurate, with SAPT96 overestimating and LM2 underestimating the magnitude of interaction energy. With an added retardation correction, SAPT96 was used (under the name SAPT2) by Janzen and Aziz^{232} to calculate properties of helium and found to be the most accurate helium potential at that time.

In 1999, van Mourik and Dunning^{233} calculated CCSD(T) energies in basis sets up to daug-cc-pV6Z, CCSDT − CCSD(T) differences in the daug-cc-pVQZ basis set, and FCI − CCSDT differences in the daug-cc-pVTZ basis set. The CCSD(T) energies were CBS-extrapolated and then refined by adding a correction equal to the *R*-interpolated differences between highly accurate CCSD(T)-R12 results (available at a few distances in Ref. 234) and the obtained CBS limits. The CC-R12 methods are analogous to CC-GTG methods, but the explicit correlation factor enters linearly.^{195} As seen in Fig. 6, the results of Ref. 233 were more accurate than any previously published ones, but no estimates of uncertainties were provided and the computed interaction energies were not fitted.

Supermolecular ECG-based calculations for He_{2} started to appear in the late 1990s,^{235,236} and were initially aimed at providing upper bounds to the interaction energies (by subtracting essentially exact monomer energies), as the authors did not attempt to extrapolate their results to the CBS limits. Another application of explicitly correlated functions to the He–He interaction was a series of papers by Gdanitz,^{237–239} who used the multireference averaged coupled-pair functional method with linear *r*_{12} factors, *r*_{12}-MR-ACPF. The extrapolated results from the last paper of the series, Ref. 239 (denoted “Gdanitz01” in Fig. 6), were among the most accurate results available at that time. However, the reported uncertainties were strongly underestimated at shorter distances (as much as 5 times at 5.6 bohrs and 17 times at 4.0 bohrs).

Another important series of papers was published by Anderson *et al.*,^{240–242} who reported quantum MC energies with progressively reduced statistical uncertainties. Although these results were obtained only for a few internuclear distances, they represented very valuable benchmarks for mainstream electronic structure methods. In fact, until the publication of the CCSAPT07 potential,^{185} the result from Ref. 242, −10.998(5) K (see “Anderson04” in Fig. 6), was the most accurate value available at 5.6 bohrs.

In Refs. 243 and 244, a hybrid supermolecular ECG/orbital method was applied to the helium dimer. The bulk of the correlation effect on the interaction energy, at the CCSD level, was evaluated using GTG functions and the method developed in Refs. 245–256. The nonlinear parameters were optimized at the MP2 level. The effects of noniterative triple excitations [the “(T)” contribution], i.e., the differences between CCSD(T) and CCSD energies, were calculated using large orbital basis sets (up to aug-cc-pV6Z with bond functions and daug-cc-pV6Z) and extrapolated to the CBS limits. Finally, the FCI corrections [differences between FCI and CCSD(T) energies] were obtained in basis sets up to aug-cc-pV5Z with bond functions and daug-cc-pV5Z, and also extrapolated. Results for three distances were reported in Ref. 244 (see “Cencek04” in Fig. 6).

Hurly and Mehl (HM) analyzed the best existing *ab initio* data for the helium dimer and created a new potential^{9} representing a compromise based on uncertainties of existing data and their mutual agreement (for instance, as can be seen in Fig. 6, the result from Ref. 244 was used at *R* = 7.0 bohrs). The diagonal adiabatic corrections from Ref. 257 were added to the final potential, which was then used to calculate the second virial coefficient, viscosity, and thermal conductivity of helium. HM recommended that the values of these thermophysical properties should serve as standards for measurements.

The CCSAPT07 potential^{185} based on the hybrid GTG/orbital method, published in 2007, was a significant improvement over the previous complete potential of this type, i.e., the SAPT96 potential.^{175,176} CCSAPT07 combined three different computational techniques, according to the criterion of the lowest uncertainty available for a given internuclear distance. Variational four-electron ECG calculations were used for *R* ≤ 3.0 bohrs and SAPT+FCI was employed for *R* > 6.5 bohrs. At intermediate distances, the hybrid supermolecular method developed in Refs. 243 and 244 and described above provided the highest accuracy. Compared to Refs. 243 and 244, several computational improvements were introduced,^{184} resulting in significantly reduced uncertainties. The SAPT calculations^{185} of CCSAPT07 followed the SAPT96 recipe, but also with larger basis sets and some computational improvements. The uncertainties of this potential were smaller than some effects that are neglected at the nonrelativistic BO level. Calculations of these effects will be discussed in Sec. 3.3.2.

Another highly accurate potential, by Hellmann, Bich, and Vogel (HBV),^{258} appeared at almost the same time as CCSAPT07. Those authors used very large basis sets (up to daug-cc-pV8Z with added bond functions at the CCSD level, and gradually smaller bases for higher levels of theory up to FCI) followed by CBS extrapolations. After augmenting the HBV potential with adiabatic, approximate relativistic, and retardation corrections, the authors used it to calculate thermophysical properties of helium.^{259} However, the uncertainties of the HBV potential were not estimated, which restricts its usefulness. A direct accuracy comparison between the pure BO component of HBV and CCSAPT07 is now possible because of the much higher accuracy of the present-day benchmark energies,^{177} and we performed such analysis using the values reported in the last column of Table 3 in Ref. 258. Out of 11 distances for which all three energies are available, the largest relative error (with respect to the results of Ref. 177), equal to 0.90%, occurs for the CCSAPT07 energy at 5.0 bohrs, while the error of the HBV energy at this distance is 0.48%. If one excludes this distance, which is close to where the helium potential crosses zero, and calculates the average relative error at the remaining distances, one obtains 0.007% for CCSAPT07 and 0.011% for HBV. Therefore, both potentials exhibit a similar accuracy and represent a significant improvement over all previously published helium dimer potentials.

The current most accurate nonrelativistic BO potential for the helium dimer (labeled as “Przybytek17” in Fig. 6) was published in Ref. 177; see Ref. 193 for details of these calculations. The significant improvement over all previous potentials was achieved by a combination of three factors. First, a pure ECG approach was used, i.e., with all four electrons explicitly correlated and no contributions calculated with orbital methods. Indeed, the residual errors of the older hybrid ECG/orbital potentials, SAPT96^{175,176} and CCSAPT07,^{185} were dominated by insufficient basis set saturation of the relatively small orbital contributions. Second, the use of the monomer contraction method,^{189,260} i.e., the use of the product of helium atoms wave functions as one of the functions in the basis set, dramatically improved the energy convergence with respect to the ECG expansion size. Furthermore, a replacement of the simple product of monomer wave functions by a more compact sum of four-electron functions optimized for two noninteracting helium atoms^{206,261} reduced the computational cost at the nonlinear optimization stage. Third, a near-complete optimization of nonlinear parameters in large basis set expansions was possible due to this reduced cost and due to other improvements of the optimization algorithm.

#### 3.3.2. Physical effects beyond the nonrelativistic BO level

^{185}it became clear that a further reduction of uncertainties required inclusion of post-BO effects. The first calculation of all relevant such effects for the whole potential-energy curve was presented in Ref. 137 and later improved in Refs. 10, 11, and 177. Some post-BO effects for the whole curve were included even earlier in Refs. 258 and 259, but this work omitted non-negligible two-electron terms in the

*α*

^{2}relativistic and

*α*

^{3}QED corrections. The helium dimer potentials of Refs. 11 and 177 include at the post-BO level the diagonal adiabatic correction, relativistic corrections (earlier computed in Ref. 189, but for the minimum separation only), the QED correction, and the retardation effect (a long-range QED correction)

*α*

^{2}Breit–Pauli operator

^{204}for the relativistic correction, and the

*α*

^{3}QED operator

^{262}for the QED correction. In the latter case, one approximation was made in the operator. In the term

*I*is over the nuclei, the value of ln

*k*

_{0}should be computed for each

*R*, but instead a constant value was taken, equal to the value of ln

*k*

_{0}for the helium atom. This is an excellent approximation since ln

*k*

_{0}depends very weakly on

*R*. Calculations for two interacting ground-state hydrogen atoms

^{263}have shown that ln

*k*

_{0}changes by less than 1.15% when

*R*varies from 1.4 bohrs, the distance of the potential minimum, to infinity, where it assumes the atomic value. For H

_{2}, this

*R*-dependence is important since its inclusion changes the dissociation energy by 0.004 cm

^{−1}, while the uncertainty of this quantity is 0.001 cm

^{−1}. This inclusion changes the value of the QED term by 1.8%. The same relative change for He

_{2}would result only in a 0.000 02 K contribution to the interaction energy at the minimum of the potential, negligible compared to the uncertainties coming from other sources.

All post-BO corrections were computed using both four-electron ECG basis sets and orbital basis sets (except for the so-called Araki–Sucher part of the QED operator where only ECG functions were used). The calculations with smaller uncertainties were selected for the final potential. Orbital calculations were performed using a combination of CCSD(T) and FCI approaches or FCI alone. For the adiabatic correction, only FCI was used. The calculations of the average values of the operators listed above with ECG and FCI wave functions are straightforward (although regularization techniques have to be used for singular operators). However, the CCSD(T) wave function needed to compute expectation values is not available (not defined) and instead the CCSD(T) linear response method was used.^{264} The retardation effects of long-range electromagnetic interactions were computed from the Casimir–Polder formula^{265} by subtracting the retardation part of the *α*^{2} relativistic and *α*^{3} QED corrections.^{190}

The calculations of Ref. 177 significantly improved the accuracy of the helium dimer potential, with uncertainties reduced by an order of magnitude compared to those of Refs. 137 and 185. As already discussed, the main improvement was due to the use of larger and better optimized ECG wave functions at the nonrelativistic BO level of theory for all *R* ≤ 9 bohrs. Accuracy of the adiabatic and relativistic corrections was also improved by using larger basis sets than in Refs. 10 and 137. A major theoretical advance was the calculation of the properties of the very weak bound state of He_{2} (the so-called halo state) with full inclusion of nonadiabatic effects.

The accuracy of relativistic and QED contributions was further improved in Ref. 11. The contributions to the interaction energy at the van der Waals minimum are presented in Table 2. Clearly, with the uncertainty of the BO contribution of 0.000 20 K, all the included post-BO contributions are relevant, except for the retardation contribution, but this contribution does become important at very large separations.^{190} One can also see that uncertainties of the adiabatic, relativistic, and QED terms are almost negligible compared to the uncertainty of the BO term. The potential of Ref. 11 was used to compute the second virial coefficient and the second acoustic virial coefficient of helium.

Contribution . | Value . | Uncertainty . |
---|---|---|

V_{BO} | −11.000 71 | 0.000 20 |

V_{ad} | −0.008 904 8 | 0.000 009 7 |

V_{rel} | 0.015 391 1 | 0.000 015 4 |

V_{QED} | −0.001 332 7 | 0.000 001 8 |

V_{ret} | 0.000 012 |

Contribution . | Value . | Uncertainty . |
---|---|---|

V_{BO} | −11.000 71 | 0.000 20 |

V_{ad} | −0.008 904 8 | 0.000 009 7 |

V_{rel} | 0.015 391 1 | 0.000 015 4 |

V_{QED} | −0.001 332 7 | 0.000 001 8 |

V_{ret} | 0.000 012 |

### 3.4. Nonadditive helium potentials

In any fluid, the total interaction energy includes terms beyond pairwise-additive interactions between monomers. These so-called nonadditive contributions begin with three-body nonadditive terms defined as the part of the trimer interaction energy that cannot be recovered by the sum of two-body interactions. The additive and nonadditive interactions form a series called the many-body expansion of interaction energy. Fortunately, for all fluids consisting of monomers interacting via noncovalent forces, this expansion converges very fast and usually it is sufficient to limit calculation to two- and three-body terms. For a review of the many-body expansion, see Ref. 266. For metrology, the three-body potential is needed to calculate the third virial coefficient.

A pairwise-nonadditive potential for helium was developed in Ref. 267 and improved in Ref. 138. In the earlier work, two independent potentials were obtained. One was based on three-body SAPT^{268–272} and the other on the supermolecular CCSD(T) approach. Orbital basis sets up to aug-cc-pV5Z were used. The two potentials were in very good agreement. In Ref. 138, the CCSD(T) potential was improved by calculating the FCI correction in an incremental approach and increasing the number of grid points, with CCSD(T) values taken from Ref. 267 except for the new grid points. Near the minimum of the total potential, the three-body contribution is only −0.0885 K, which should be compared to the total interaction energy of about $\u221233$ K, but the three-body contribution is much larger than the uncertainty resulting from the two-body term, which is 0.0006 K. The uncertainty of the three-body term at the minimum of the total potential was estimated to be 0.002 K.

Recently, the three-body potential for helium was further improved^{273} by adding the relativistic and adiabatic corrections, as well as using a new set of correlation-consistent basis sets specifically developed for helium atoms.^{177} An improved functional form was also used to analytically represent the potential at large distances. In particular, new terms were developed for the case when two atoms remain close while the third is progressively more distant. These refinements resulted in a reduction of the uncertainty by a factor of about 5 overall. In particular, the uncertainty at the minimum was reduced to 0.5 mK, a factor of 4 smaller than that of the previous work.^{138}

### 3.5. Heavier noble-gas atoms

While theory is superior to experiment for the helium atom and helium clusters, this is not the case for neon, and even less so for argon. The simple reason is the number of electrons per atom: 2, 10, and 18, respectively. While for the helium atom and small helium clusters the *N*-electrons explicitly correlated bases can reach ppm or smaller uncertainties, and FCI calculations can be performed in fairly large bases, for neon neither type of calculation can be performed in bases large enough to get meaningful results. To quantify this statement, let us examine the most accurate calculations for the neon dimer,^{117} see Table 3. The calculations at the CCSD(T) level of theory were performed in the largest available basis sets: modified daug-cc-pV8Z with bond functions. The uncertainty of the interaction energy obtained in this way is about 200 ppm, which is only ten times larger than the 20 ppm uncertainty of the He_{2} BO interaction energy. However, uncertainties coming from some excitations of higher rank are significantly larger: the pentuple excitation contribution, Δ(P), increases the uncertainty of the total value of interaction energy to about 1000 ppm. The increase of uncertainties is due to the use of smaller and smaller basis sets as the number of excitations increases: at the CCSDTQ(P) level of theory only the daug-cc-pVDZ basis set could be used. Furthermore, based on the results in Table 3, it is not possible to estimate the uncertainty resulting from neglecting excitations beyond (P). The lower part of Table 3 shows the convergence in the rank of excitation. One can see that while the contribution of the triple excitations is very substantial, a 29% increase in the magnitude of interaction energy relative to the CCSD level, the contribution of quadruple excitations is 57 times smaller than that of triple ones. However, the contribution of pentuple excitations breaks this fairly fast convergence: it is of similar magnitude to that of the quadruple excitations. Note that one cannot blame the noniterative character of the pentuple excitations, as for lower-rank excitations the iterated and noniterated values are fairly similar. One may ask if the value of the pentuple contribution computed in Ref. 117 could be a numerical artifact resulting from the use of a rather small basis set. This issue was investigated in Ref. 117, and the results computed in the aug-cc-pVDZ and aug-cc-pVTZ were 0.0227 and 0.1113 K, respectively. While these results may indicate that even the first digit in the pentuple excitations contribution may be uncertain, they also indicate that the order of magnitude will likely remain the same when going to larger basis sets. This would indicate that for Ne_{2} the coupled-cluster expansion converges very slowly, whereas for other closed-shell systems investigated in the literature CCSDTQ(P) agrees with FCI very well, indicating that effects of higher excitations are negligible. Unfortunately, FCI calculations would be extremely difficult to perform for Ne_{2} even in the aug-cc-pVDZ basis set.

Contribution . | Value . | Uncertainty . |
---|---|---|

CCSD(T) | −41.3301 | 0.0100 |

CCSDT–CCSD(T) | −0.5730 | 0.0115 |

CCSDT(Q)–CCSDT | −0.1602 | 0.0112 |

CCSDTQ–CCSDT(Q) | −0.0043 | 0.0009 |

CCSDTQ(P)–CCSDTQ | 0.1179 | 0.0589 |

CCSD | −32.5355 | |

CCSDT–CCSD | −9.4437 | |

CCSDTQ–CCSDT | −0.1645 | |

CCSDTQ(P)–CCSDTQ | 0.1179 |

Contribution . | Value . | Uncertainty . |
---|---|---|

CCSD(T) | −41.3301 | 0.0100 |

CCSDT–CCSD(T) | −0.5730 | 0.0115 |

CCSDT(Q)–CCSDT | −0.1602 | 0.0112 |

CCSDTQ–CCSDT(Q) | −0.0043 | 0.0009 |

CCSDTQ(P)–CCSDTQ | 0.1179 | 0.0589 |

CCSD | −32.5355 | |

CCSDT–CCSD | −9.4437 | |

CCSDTQ–CCSDT | −0.1645 | |

CCSDTQ(P)–CCSDTQ | 0.1179 |

Similar calculations at the limits of the available technology were reported for Ar_{2} in Ref. 274 (see also earlier calculations^{275} with accurate treatment at very small values of interatomic distances *R*). The value of the interaction energy obtained at the van der Waals minimum is −142.86 K and its uncertainty was estimated at 0.46 K. This uncertainty, representing 3000 ppm (0.3%) of the computed well depth, does not include an estimate of effects beyond CCSDTQ. The results of Ref. 117 for Ne_{2} indicate, however, that the post-CCSDTQ contribution may be not negligible.

The first first-principles three-body potential for argon was developed in Ref. 269 using three-body SAPT. It was then used to compute the third virial coefficient of argon^{144} and to simulate vapor–liquid equilibria.^{153} An improved three-body potential for argon was developed in Ref. 276 using the CCSDT(Q) level of theory and including core correlation and relativistic effects. Uncertainties of the potential were estimated. The authors of Ref. 276 also computed the third virial coefficient, obtaining good overall agreement with experimental data. In particular, in some regions of temperature, theoretical values exhibited smaller uncertainties than experiment and comparisons with theory allowed evaluation of different experiments. When the experimental data were refitted by a new model that included an approximate fourth virial coefficient,^{123} the agreement with theory improved, which can be considered to be a validation of the new model. The work of Ref. 276 shows that despite limitations of accuracy, for some properties of argon theory may provide information relevant for metrology and its accuracy may be competitive with experimental accuracy.

### 3.6. Magnetic susceptibility

Magnetic susceptibilities of noble gas atoms are relevant for RIGT; see Eq. (11). In general, the magnetic susceptibility is several orders of magnitude smaller than its electric counterpart (hence, *A*_{μ} is several orders of magnitude smaller than *A*_{ɛ}). This means that only modest accuracy for the magnetic susceptibility, on the order of 0.1% or even 1%, is sufficient for it to make a negligible contribution to the uncertainty budget of current or planned refractivity-based thermodynamic metrology. Calculations at the BO level are therefore probably sufficient, but it is still desirable to compute additional effects, at least at lowest order, to verify that they are relatively small.

The first comprehensive calculation of the magnetic susceptibility of the helium atom was performed by Bruch and Weinhold.^{277} They added corrections for relativistic effects and nuclear motion to an existing high-accuracy calculation at the BO level. However, their calculation included only some of the relativistic corrections that enter at lowest order. Recently, Puchalski *et al.*^{80} presented a definitive calculation of all effects through order *α*^{4}, along with a more accurately computed value for the nonrelativistic BO limit using Slater geminals. They obtained agreement within mutual uncertainties with the calculations of Bruch and Weinhold for individual terms,^{277} but included some terms that had been omitted in the previous work. When converted from the atomic units used in the paper, the final result for ^{4}He corresponds to *A*_{μ} = −7.921 28(13) × 10^{−6} cm^{3} mol^{−1}. The relative uncertainty of this result, primarily due to neglected QED effects that enter at the *α*^{5} level, was conservatively estimated at 16 ppm. This is far more than sufficient for any conceivable application of refractivity for temperature or pressure measurement, and the agreement with previous work encourages confidence in the result.

As with other properties, the greater number of electrons renders the calculation of magnetic susceptibility much more difficult for neon and especially for argon. The current state-of-the-art calculations for neon^{64} and argon^{66} were performed only at the nonrelativistic BO level, with a rough uncertainty estimate for neglected relativistic effects based on the magnitude of those effects for the electric polarizability. The estimated uncertainty of this calculated quantity was ∼0.2% for neon^{64} and 1% for argon.^{66} The limited experimental information for the magnetic susceptibility is discussed in Sec. 4.5.3.

## 4. From Electronic Structure to Thermophysical Properties

Virial expansions are exact results from quantum statistical mechanics which enable a systematically improvable evaluation of various thermophysical properties as a power series in density starting from the ideal-gas reference system. The coefficients appearing in the *N*-th term of the series can be computed from the knowledge of the interaction of clusters of *N* particles.

*p*as a function of density

*ρ*– one obtains Eq. (7)

^{278,279}together with rigorous expressions for the virial coefficients

*B*(

*T*),

*C*(

*T*),

*D*(

*T*), etc., which turn out to be functions of temperature only and are given by

*Q*

_{N}(

*V*,

*T*) is the partition function of a system of

*N*particles evaluated in the canonical ensemble. These partition functions can be calculated once the interaction potential

*U*

_{N}(

**x**

_{1}, …,

**x**

_{N}) among

*N*particles is known; the potential is generally expressed as

*u*

_{2}is the pair potential,

*u*

_{3}is the non-additive contribution to the three-body potential, and so on. In Eq. (42) we specialized to the case of atomic systems, which is the principal topic of this review; in this case

**x**

_{i}represents the position of the

*i*-th atom. In the case of molecules, which we will discuss in Sec. 5, the various potentials appearing in Eq. (42) depend also on coordinates

**y**

_{i}that describe the intramolecular configuration of molecule

*i*. In particular, a single-body potential

*u*

_{1}(

**y**

_{1}) will also appear in Eq. (42). The potentials and their uncertainties can be computed from first principles using the methods described in Sec. 3. The most general expression for

*Q*

_{N}(

*V*,

*T*) in quantum statistical mechanics is given by

*i*⟩ of the

*N*-body Hamiltonian

*H*

_{N}with the proper symmetry upon particle exchange due to the bosonic or fermionic nature of the particles involved. Equation (44) is an equivalent expression where the sum over the states has no restriction on the symmetry and the operators $Pj$ represent the

*j*-th permutation of particles in the Hilbert space, including the sign of the permutation in the case of fermions. The latter expression will be the most convenient when discussing the path-integral MC approach for the calculation of virial coefficients.

^{280,281}The non-relativistic

*N*-body Hamiltonian is conveniently written as

*π*_{i}and mass

*m*

_{i}for the

*i*-th particle and the second equality defines the

*N*-body kinetic energy

*T*

_{N}.

^{282}

*Q*is

*ɛ*

_{r}is generally given as a generalization of the Clausius–Mossotti equation in one of the two equivalent forms given by Eqs. (8) and (9). Until recently, derivations for the coefficients appearing in these equations would agree on the expression for the second dielectric virial coefficient,

*B*

_{ɛ}, but differ in the case of the higher-order coefficients.

^{283–285}A systematic review of the dielectric expansion showed that the correct expressions are

^{207,285}

*α*

_{1}is the atomic polarizability and the functions

*Z*

_{N}are given by expressions similar to Eq. (41), where the interaction Hamiltonian among the constituent particles of Eq. (45) is extended with two terms in order to include the effect of the interactions of the dipole moment and the electronic polarizability of the system with an external electric field of magnitude

*E*

_{0}. In Eqs. (49)–(51), the derivatives are to be evaluated at

*E*

_{0}= 0. The two additional terms in the Hamiltonian are

**m**

_{n}and

*α*_{n}are the (non-additive) dipole moments and the (non-additive) electronic polarizabilities of a system of

*n*particles. In the case of atoms,

**m**

_{1}and

**m**

_{2}are both zero, but a system of three particles has, in general,

**m**

_{3}≠ 0.

^{286,287}

*n*and is given in Eq. (11). The Lorentz–Lorenz equation (11) is relevant to those experiments where the refractive index is measured by optical methods. In this case, the refractive virial coefficients are a function of the angular frequency

*ω*of the electromagnetic radiation as well as the temperature. Usually, the frequency dependence is approximated as a power-law expansion of the form

*S*(−4).

^{288}

### 4.1. Classical limit

Although the focus of this review is on calculations with no uncontrolled approximation, let us briefly discuss the classical limit of the approach we have outlined. Classical expressions can be computed relatively easily, and provide a useful high-temperature check for the more involved calculations described below.

*N*! in the partition functions.

^{279}In the same limit, the kinetic term in the Hamiltonian (45) commutes with the potential energy

*U*

_{N}as well as with $HNdip$ and $HNpol$. Its contribution can be integrated exactly, resulting in a factor of the form $VN/\Lambda m3N$ where $\Lambda m=h/2\pi mkBT$ is the thermal de Broglie wavelength of the atoms under consideration. Putting all of this together, one obtains

*α*_{1}. Additionally, we have denoted by d

**X**

_{N}the integration element in the space of all the coordinates needed to describe a system of

*N*atoms, e.g., the Cartesian coordinates

**x**

_{1}, …,

**x**

_{N}. Since the system is translationally invariant, the integration produces a factor of

*V*with the understanding that one particle, usually labeled as 1, is fixed at the origin of the coordinate system. Using rotational invariance, one can further write for the integration elements

*r*

_{ij}= |

**r**

_{ij}| = |

**x**

_{i}−

**x**

_{j}| and

*θ*

_{ij}is the angle between the vectors

**r**

_{i1}and

**r**

_{j1}. In Eq. (58), the angle

*ϕ*is the polar angle corresponding to the vector

**r**

_{14}in spherical coordinates.

*α*

_{2}is substituted by the Cauchy moment Δ

*S*(−4).

*γ*

_{a}is more involved and is given in the Appendix.

### 4.2. Quantum calculation of virial coefficients

The classical approach can be expected to be valid when Λ/*σ* ≪ 1, where *σ* is the size of the hard-core repulsive region of atoms (which is around 6 bohrs for the noble gases); this implies that the classical formulas will be asymptotically valid for high temperatures and heavy atoms. However, in the case of helium this approximation is too drastic even at room temperature.

The inclusion of quantum effects in the calculation of virial coefficients (density, acoustic, or dielectric) requires evaluating the *N*-body partition functions *Q*_{N} of Eq. (43) in a quantum framework. A straightforward approach would be to consider in Eq. (43) the eigenstates |*i*⟩ of the *N*-body Hamiltonian, *H*_{N}|*i*⟩ = *E*_{i}|*i*⟩, so that Eq. (43) becomes a simple sum. To the best of our knowledge, this method has been demonstrated to date only in the case of the second dielectric virial coefficient.^{79}

*Q*

_{2}(which enables the calculation of virial coefficients of order 2), a very fruitful approach dating back to the late 1930s

^{278,289}is to rewrite it as the sum of three terms: one depending on the bound-state energies, one depending on the phase shifts of the scattering states, and one depending on the bosonic or fermionic nature of the atoms involved. The expression of

*B*(

*T*) then becomes

*μ*is the reduced mass of the pair of atoms considered, $Enlbound$ is the energy of the

*n*-th bound state with relative angular momentum

*l*, and

*f*(

*I*,

*l*) = 1 + (−1)

^{2I+l}/(2

*I*+ 1) with

*I*the nuclear spin in the case of identical atoms (the case of different atoms can be recovered by letting

*I*→ ∞). The quantity

*δ*

_{l}(

*E*) in Eq. (66) is the

*absolute*scattering phase shift for two particles with relative energy

*E*and angular momentum

*l*. Absolute phase shifts are continuous functions of

*E*that tend, in the limit of

*E*→ 0, to

*π*times the number of bound states at angular momentum

*l*. With the advent of electronic computers, the use of Eqs. (64)–(68) enabled the calculation of accurate numerical values

^{290,291}and it is still the most efficient way to compute the second virial coefficient of atomic species.

^{10,11}One important benefit of this method is that once the energies of all the bound states have been computed and phase shifts are known for a sufficiently high number of total angular momenta and scattering energies, the values of

*B*(

*T*) and its derivatives, and hence

*β*

_{a}(

*T*), can be easily computed at all temperatures; knowledge of the collision-induced pair polarizability also enables the calculation of

*B*

_{ɛ}.

^{79}Additionally, transport properties such as the viscosity and the thermal conductivity – see Sec. 4.6 below – can be computed in a straightforward manner.

Unfortunately, this approach cannot be easily extended to higher-order coefficients. Some attempts in this direction were made in the 1960s,^{292,293} but all of them required the introduction of some uncontrolled approximations and did not take into account the non-additive parts of the many-body potential.

#### 4.2.1. Path integral approach

^{280}was shown by Fosdick and Jordan to provide a systematic way to compute virial coefficients of any order without any uncontrolled approximation.

^{294,295}The path-integral formulation is based on a controlled approximation of the exponential of the

*N*-body Hamiltonian, that is

Equation (70) is the Li–Broughton expansion of the exponential of the sum,^{296} which was independently discovered by Kono *et al.*^{297} based on an initial idea by Takahashi and Imada.^{298} It can be shown that Eq. (70) becomes an exact equality in the case *P* → ∞, although in practice satisfactory convergence is reached for a finite value of the parameter *P*. Actually, Eq. (70) becomes an equality in the *P* → ∞ limit also when *O* is omitted in Eq. (70) (this is the original Trotter–Suzuki approach),^{299,300} although in this case convergence requires higher values of *P*; this approach is called the *primitive approximation*,^{281} and, for the sake of simplicity, will be used throughout this review.

*P*− 1 additional completeness relations between the

*P*factors in Eq. (70). Additionally, one uses as a complete set the (generalized) position eigenstates $|i\u232a=|XN(1)\u232a$, where we have included a superscript (1) for later convenience. In this case, the sum over

*i*in Eq. (44) becomes an integral over the 3

*N*coordinates $XN(1)$ and the

*P*− 1 completeness relations can be written as

*k*= 2, …,

*P*. Notice that in this case the effect of the permutation operators $Pj$ is to exchange atomic coordinates in the rightmost ket. For example, if $P(12)$ denotes the permutation of particles 1 and 2 (assumed to be bosons), one has

*T*≳ 10 K even in the case of helium) and the case of density virials of pure species [so that our Hamiltonian is given by Eq. (45) with

*m*

_{i}=

*m*]. The operators

*U*

_{N}(and, if needed,

*O*) of Eq. (70) are diagonal in the position basis. The matrix elements of the exponential of the kinetic energy operators can be calculated exactly

^{301}and are given by

*Z*

_{N}can be written as

*P*→ ∞ limit) to the original quantum statistical formulation, can be interpreted as the partition function of a

*classical*system.

^{301}For each of the original

*N*particles of coordinates $xi(1)$, one has introduced

*P*− 1 copies of coordinates $xi(k)$, which, as one can see from Eq. (77), are connected via harmonic potentials. The equivalent classical system is then made by

*N*ring polymers of

*P*monomers each. As shown by Eq. (76), these polymers interact with the original potential averaged over all the monomers. It can be shown that the functions

*F*

_{i}of Eq. (77) represent probability distributions.

^{302}Although they are not Gaussian probabilities, because of the ring-polymer condition $xi(P+1)=xi(1)$, they can be sampled exactly using an interpolation formula due to Levy

^{294,295}(also known as “the Brownian bridge”). The harmonic intra-polymer interaction, which ultimately comes from the kinetic energy term

*T*

_{N}of the quantum Hamiltonian (45), has the effect that the average “size” of the ring-polymer corresponding to each particle is of the order of the de Broglie thermal wavelength Λ, thus taking into account quantum diffraction (that is, the Heisenberg uncertainty principle).

In order to compute the functions *Z*_{N} (and, hence, the virial coefficients), it is convenient to separate the *NP* vector coordinates $xi(k)$ as follows: first of all, we notice that the energy of the equivalent classical system is invariant upon an overall rigid rotation or rigid translation. We can use the latter property to extract a factor of *V* and at the same time pin one of the coordinates – conventionally the first monomer of particle 1, that is $x1(1)$ – at the origin of the coordinate system. The rotational invariance can be taken into account by assuming that the first monomer of one particle (particle 2, say) lies along the *x* axis of the coordinate system and that the first monomer of another particle (particle 3) lies in the *xy* plane. This convention brings about a factor of 4*π* when *N* = 2 (corresponding to the integration over the two polar angles describing $x2(1)$) and a factor of 8*π*^{2} (that is the integration over the two polar angles describing $x2(1)$ and the azimuthal angle of $x3(1)$) when *N* ≥ 3. The remaining 3*NP* − 6 coordinates (or 3*NP* − 5 in the case of *N* = 2) can be conveniently divided into

The coordinates of the first bead of all the particles, that is $r12=|x2(1)\u2212x1(1)|$ and, for

*N*≥ 3, $r13=|x3(1)\u2212x1(1)|,cos\theta 23$ and $xi(1)$ (the latter only for*N*≥ 4), where*θ*_{23}is the angle between the position of particles 2 and 3 in the*xy*plane.The relative coordinates $\Delta ri(k)$ (

*k*= 1, …,*P*− 1).

*F*

_{i}depend only on $\Delta ri(k)$, one can rewrite the partition functions

*Z*

_{N}of Eq. (75) in the form

^{79}Explicit expressions for the third acoustic virial coefficient in the path-integral formulation are quite cumbersome, for reasons discussed in the Appendix; they can be found in Ref. 303.

It is important to notice that in the case of *C*(*T*) the terms coming from $Z22$ in Eq. (39) actually involve averages over *four* ring polymers, since these two terms involve two particles each and have to be treated as *independent*, lest spurious correlations be introduced in the calculation of the ⟨⋯⟩ average. In fact, in the last term of Eq. (85) two of these polymers are used to compute $e\u2212\beta U2\u0304(r12)\u22121$ and the other two to compute $e\u2212\beta U2\u0304(r13)\u22121$. Similar considerations also apply when calculating *γ*_{a} and *C*_{ɛ} using path integrals.

Quantum effects are taken into account by averaging over the ring-polymers configurations, and at the same time evaluating the interaction energy as an average over the monomers, as in Eq. (76). We recall that in Eqs. (83) and (85) the radial variables $rij=|xi(1)\u2212xj(1)|$ are the distances between the first monomer of particles *i* and *j*. In the classical limit, the size of the ring polymers shrinks to zero so that one recovers the results of Sec. 4.1.

It is worth noting that one can find several semi-classical approximations of the exact path-integral expressions of Eqs. (83)–(86). In general, they can be obtained by expanding the full quantum-mechanical results in powers of *ℏ*^{2}, where the first term is the classical one. This approach was pioneered by Wigner and Kirkwood^{304,305} and subsequently developed by Feynman and Hibbs,^{280} who put forward the idea of estimating semiclassical values by using the classical expressions with suitably modified (and temperature-dependent) potentials. Although the Feynman–Hibbs approach considered systems with pair potentials only, a systematic derivation of semiclassical expressions in the case of three-body interactions has been developed by Yokota.^{306} Even if semiclassical approaches introduce uncontrolled approximations, they are quite effective in the case of heavier atoms such as argon at high temperatures and provide a useful check for the fully quantum calculations.

#### 4.2.2. Exchange effects

*i*under the action of permutation $P$. This is equivalent to saying that some of the ring polymers would coalesce into larger polymers, depending on the specific permutation that is being considered in the sum of Eq. (44). These larger ring polymers are still described by probability distributions similar to those of the Boltzmann case, that is Eq. (77). As an illustrative example, let us see how the probability distribution for the internal coordinates of particles 1 and 2 is modified in the presence of exchange for bosons of spin 0. Defining $Ri=x1(i)$ and $Ri+P=x2(i)$ for

*i*= 1, …,

*P*as well as Δ

**R**

_{i}=

**R**

_{i+1}−

**R**

_{i}(notice that

**R**

_{2P}=

**R**

_{1}because we are considering the permutation involving only particle 1 and 2), and $\Lambda \mu =2\Lambda $, the kinetic energy terms that would give rise to the probabilities

*F*

_{1}

*F*

_{2}can be written as

*P*monomers describing a particle of mass

*μ*=

*m*/2 at the same temperature [cf. Eq. (77)]. In the case of the second virial coefficient, where this is the only exchange term present, this contribution is just a simple average over the larger polymer, and can then be written as

^{307}

### 4.3. Uncertainty propagation

As is apparent from their definition, the calculation of virial coefficients depends on the knowledge of few-body properties of atoms, namely interaction potentials, polarizabilities, and dipole moments. In a completely *ab initio* calculation of virial coefficients, these quantities – as seen in Sec. 3 – are determined by electronic-structure calculations and are provided with a full uncertainty estimation. In this section, we will show how this uncertainty can be propagated to the uncertainty in virial coefficients, using the third virial coefficient *C*(*T*) as an example.

*C*(

*T*) using perturbed pair and three-body potentials, that is:

*δu*

_{i}for

*i*= 2 or

*i*= 3 in the case of the pair and three-body potential, respectively – are given as expanded (

*k*= 2) uncertainties. Assuming that a (

*k*= 2) perturbation of the potential results in a (

*k*= 2) perturbation of the virial coefficient, one fourth of the absolute value of the difference, that is $\delta C[ui]$ in Eq. (92), is interpreted as a standard (

*k*= 1) uncertainty. The overall standard uncertainty in

*C*(

*T*) due to the uncertainty in the potentials is then obtained as a sum in quadrature

^{127,308}it is unsatisfactory for several reasons. First of all, it considers only rigid shifts of the potentials, while in principle the actual potential can be closer to the upper bound for some configurations and closer to the lower bound for others. Secondly, the uncertainty (92) is obtained as a difference of quantities which are themselves computed with some statistical uncertainty. This requires very long runs to make sure that the difference in Eq. (92) is not influenced by the statistical error in the calculation of $C\xb1[ui]$.

*T*as well as functionals of the potentials.

^{126}A variation

*δu*

_{i}in the potential will then produce a corresponding variation in the value of the virial coefficient, given by

*δC*/

*δu*

_{i}. The absolute value in Eq. (94) comes from the conservative choice of assuming that all the variations will contribute with the same (positive) sign to the final uncertainty. We note in passing that in the case of the second virial coefficient,

*B*(

*T*), Eqs. (94) and (92) produce the same result. As is apparent from Eq. (39), the evaluation of Eq. (94) requires the functional derivative of

*Z*

_{3}and

*Z*

_{2}with respect to the pair and three-body potential. As a first approximation, one can use the classical expression Eq. (55) (possibly augmented by semiclassical corrections).

^{126}More accurate results (especially at low temperatures) are obtained by functional differentiation of the path-integral expressions (79) and (80) so that one has

^{303}

^{,}

^{207,210}

In actual practice, these expressions enable rigorous estimation of the uncertainty propagated from the potentials with a much smaller computational effort than that needed to compute virial coefficients. Additionally, the *a priori* knowledge of a lower bound on the uncertainty and its temperature dependence facilitates the process of finding the optimal set of parameters for the path-integral simulations (cutoff distance, number of beads *P*, number of MC integration points) in order to make the statistical uncertainty of the calculation a minor contributor to the total uncertainty.

### 4.4. Mayer sampling and the virial equation of state

Equations (38)–(40) show that the expressions for the virial coefficients become more involved when the order is increased. Although these expressions can be systematically derived using computer-algebra systems, their subsequent implementation in classical or quantum frameworks becomes more and more time-consuming. Taking also into account the limited availability of *ab initio* many-body potentials (at the time of this writing, these are limited to three bodies and have been developed only for a small set of atoms and molecules), it might seem that a fully *ab initio* calculation of the equation of state using virial expansions could not be feasible. Nevertheless, it is observed that the largest contributions to the value of the virial coefficients come from the many-body potentials of lower orders, as already discussed in Sec. 2.4. As a consequence, even if only pair and three-body potentials are available, a calculation of higher-order virial coefficients can provide useful and reasonably accurate representations of the equation of state.^{309,310}

A very efficient procedure to perform this task is based on the diagrammatic approach by Ursell^{311} and Mayer,^{312,313} who showed how the various terms contributing to the virial coefficients can be related to simpler cluster integrals that can be cataloged using a diagrammatic form. The contributions from the diagrams can be added very efficiently using MC sampling methods.^{314} Although the number of diagrams increases exponentially with the order of the virial coefficient, it has been shown that calculations can be kept within a manageable size up to virial coefficients of order 16,^{315–317} resulting in equations of state with very good accuracy up to the binodal (condensation) density.^{310}

Mayer sampling methods, originally developed for monatomic systems, have been extended to molecules^{318} and therefore can also be used to perform path-integral calculations of density^{121,319} and acoustic virial coefficients.^{136} This approach provides an independent validation of the framework outlined in this review. Virial coefficients calculated using both approaches are found to be compatible within mutual uncertainties.^{303}

### 4.5. Numerical results for virial coefficients

As seen in Sec. 4.2, a fully first-principles calculation of virial coefficients requires the knowledge of many-body potentials and, in the case of dielectric properties, polarizabilities, which can be obtained by *ab initio* electronic structure calculations. Currently, as discussed in Sec. 3, the only system for which these calculations can be made without uncontrolled approximations is helium. Much effort has been devoted to produce high-quality potentials from first principles. At the time of writing, the most accurate pair potential is the one developed by Czachorowski *et al.*,^{11} which includes relativistic and QED effects. This potential was developed using exactly the same approach as the potential of Ref. 177, the only difference being that the relativistic and QED corrections were computed using a larger basis set. As a consequence of including the adiabatic corrections and recoil terms, slightly different pair potentials are available for the ^{4}He–^{4}He, ^{3}He–^{3}He, and ^{4}He–^{3}He interactions.

Recently, a new three-body potential for ^{4}He, including relativistic effects, has been developed,^{273} resulting in a significant increase of accuracy with respect to the previous non-relativistic potential (see Sec. 3.4).^{138} In the case of dielectric properties, the single-atom polarizability has been calculated with outstanding accuracy.^{192} The most accurate pair-induced polarizability currently available is that of Cencek *et al.*^{206} and, recently, fully *ab initio* calculations of the three-body polarizability^{208} and dipole moment^{210} have been performed, enabling a calculation of the third dielectric virial coefficient with well-defined uncertainties completely from first principles.^{210}

In the case of neon, the most recent pair potentials and polarizabilities have been computed by Hellmann and coworkers.^{65,117} Parametrizations of three-body potentials have appeared in the literature,^{320} but no first-principles calculations have been published so far.

Due to its easy accessibility and large measurement effects, argon has been the subject of many theoretical studies. However, the large number of electrons prevents calculations of potentials and polarizabilities with the same accuracy as the lighter noble gases, and some uncontrolled approximations are still necessary. The most accurate pair potential so far has been developed by Lang *et al.*,^{118} while a three-body potential with well-characterized uncertainties was computed and characterized by Cencek and co-workers.^{276} Regarding dielectric properties, the most accurate pair polarizability is the one developed by Vogel *et al.*^{164} In the case of neon and argon, no three-body polarizabilities are available. Calculations have been performed using the *superposition approximation*^{321,322} for the three-body polarizability. Although the results of these calculations compare well with the available experimental data, their uncertainty is to a large extent unknown.^{207}

We report in Table 4 the most up-to-date references regarding *ab initio* calculations of virial coefficients. This table to some extent serves as an update to the table of recommended data presented by Rourke.^{323}

. | Helium . | Neon . | Argon . |
---|---|---|---|

B | Reference 11 | Reference 117 | Reference 118 |

C | References 273 and 303 | Reference 324 | Reference 276 |

D | Reference 126 | Reference 324 | Reference 123 |

β_{a} | Reference 11 | Reference 117 | Reference 118 |

γ_{a} | References 273 and 303 | ⋯ | Reference 325 |

A_{ɛ} | Reference 192 | Reference 65 | Reference 66 |

B_{ɛ} | Reference 79 | Reference 117 | References 79 and 164 |

C_{ɛ} | References 208 and 210 | Reference 207 | Reference 207 |

A_{μ}a | Reference 80 | Reference 64 b | Reference 66 b |

B_{R} | Reference 79 | References 79 and 117 c | Reference 79 |

η | Reference 10 | Reference 117 | Reference 118 |

λ | Reference 10 | Reference 117 | Reference 118 |

. | Helium . | Neon . | Argon . |
---|---|---|---|

B | Reference 11 | Reference 117 | Reference 118 |

C | References 273 and 303 | Reference 324 | Reference 276 |

D | Reference 126 | Reference 324 | Reference 123 |

β_{a} | Reference 11 | Reference 117 | Reference 118 |

γ_{a} | References 273 and 303 | ⋯ | Reference 325 |

A_{ɛ} | Reference 192 | Reference 65 | Reference 66 |

B_{ɛ} | Reference 79 | Reference 117 | References 79 and 164 |

C_{ɛ} | References 208 and 210 | Reference 207 | Reference 207 |

A_{μ}a | Reference 80 | Reference 64 b | Reference 66 b |

B_{R} | Reference 79 | References 79 and 117 c | Reference 79 |

η | Reference 10 | Reference 117 | Reference 118 |

λ | Reference 10 | Reference 117 | Reference 118 |

^{a}

Note that *A*_{R} = *A*_{ɛ} + *A*_{μ}.

^{b}

Improvement in progress; see Ref. 326.

#### 4.5.1. Density virial coefficients

The most accurate *ab initio* values of the second virial coefficients of helium for both isotopes are those computed by Czachorowski *et al.*^{11} In order to visualize the recent progress in this field, we report in Fig. 7 the evolution of the theoretical uncertainty of *B*(*T*) in the past 20 years. Theoretical and computational improvements enabled a reduction of two orders of magnitude in the relative uncertainty, which is presently on the order of 10^{−4} at low temperatures ($<10$ K) and decreases to less than 10^{−5} at higher temperatures. In general, the current theoretical uncertainties of *B*(*T*) are more than one order of magnitude smaller than the best experimental determinations.

Figure 8 shows the development of the uncertainty in the calculations of *C*(*T*) for helium in the past 12 years, starting from the first calculation with fully characterized uncertainties from 2011,^{127} whose results were independently confirmed a year later using the Mayer sampling approach.^{319} One can clearly see that the subsequent improvement of the pair potential resulted in a reduction of the uncertainty at the lowest temperatures (*T* ≲ 50 K), while the uncertainty at the highest temperatures is dominated by the propagated uncertainty from the three-body potential. Recent improvements resulted in a further reduction of the uncertainty by a factor of $\u223c5$ across the whole temperature range 10–3000 K. The current theoretical uncertainty in *C*(*T*) is a few parts in 10^{4} at high temperature, and increases to a few parts per 10^{3} below 50 K. At temperatures below ∼10 K, the theoretical uncertainty budget is dominated by the propagated uncertainty from the pair potential.

Although no well-characterized four-body potential has yet been published for helium, several groups have performed calculations of the fourth virial coefficient, *D*(*T*). Although initially the effect of the four-body potential was neglected,^{319} more recent work tried to estimate its contribution using known asymptotic values.^{126} These results are in good agreement with the limited experimental information.

In the case of neon, the most recent calculations for *B*(*T*) with a pair potential having well-characterized uncertainties^{117} resulted in a relative uncertainty at *T* = 273.16 K of *u*_{r}(*B*) = 2 × 10^{−3}. As expected, this is larger than the corresponding uncertainty for helium, due to the fact that electronic structure calculations for the heavier atoms are much more computationally demanding. Unfortunately, the three-body potential for neon is only approximately known at the moment. To the best of our knowledge, no first-principles calculation is available in the literature, and only a semi-empirical parametrization is currently known.^{320} As a consequence, no *ab initio* calculation of higher-order coefficients has been performed to date and only approximate values are known.^{324}

The pair potential of argon is well characterized and has been calculated independently by two groups,^{155,274} and hence thermophysical properties at the pair level are well characterized.^{29,164,325} The relative uncertainty of *B*(*T*) at *T* = 273.16 K is *u*_{r} ∼ 0.6%. The pair potential has recently been improved by including relativistic effects, but the uncertainty of the resulting second virial coefficients is still larger than for the best experimental determinations.^{118}

The three-body potential for argon has also been computed independently by two groups^{123,276} and its uncertainty has been rigorously assessed. Therefore, the third virial coefficient of argon is also known with rigorously propagated uncertainties. The relative uncertainty is on the order of *u*_{r} ∼ 1% at *T* = 273.16 K and increases up to *u*_{r} ∼ 6% at *T* = 80 K. Analogously to the other noble gases, the four-body (and higher) non-additive contribution to the potential energy of argon is not known from first principles. Nevertheless, higher-order virial coefficients for argon, up to the seventh, have been computed based on pair and three-body potentials.^{123}

#### 4.5.2. Acoustic virial coefficients

The situation regarding first-principles calculations of acoustic virial coefficients closely follows that of the density virials. In the usual approach using phase shifts, the calculation of *B*(*T*) also provides the temperature derivatives needed to compute *β*_{a}(*T*), and therefore very accurate values for the second acoustic virials for helium,^{11} neon,^{117} and argon^{29,118,164} can be found in the papers where the pair potential and *B*(*T*) calculations are reported.

In the case of the third acoustic virial coefficient, the situation is similar. The most accurate values of *γ*_{a} for helium isotopes are reported in Refs. 273 and 303, which are in very good agreement with the values obtained independently using the Mayer sampling approach.^{136} The current relative uncertainty in *γ*_{a} for helium from *ab initio* calculations is *u*_{r} ∼ 0.02% − 0.2% across the temperature range from 10 to 1000 K.^{303}

As already mentioned, the lack of an accurate three-body potential for neon has prevented a fully first-principles calculation of the third virial coefficient, and hence no *ab initio* values of *γ*_{a} are currently available for neon.

Regarding argon, *ab initio* acoustic virial coefficients up to the fourth, together with a thorough analysis of their associated uncertainties, have been reported by Wiebke *et al.*^{325} The uncertainty of *γ*_{a} at *T* = 273.16 is $\u223c1.4%$.

#### 4.5.3. Dielectric and refractivity virial coefficients

The first dielectric virial coefficient *A*_{ɛ} for helium has been computed in Ref. 192 with an accuracy exceeding the best experimental determination. In the case of neon and argon, the most accurate theoretical results are less accurate than the best experimental determination.^{63} The most accurate computed value for neon can be found in Ref. 65, and a calculation for argon, including the frequency dependence needed for refractivity estimates, has recently appeared.^{66}

Magnetic susceptibilities computed from first principles and the corresponding quantities *A*_{μ} that are used in RIGT are available for helium,^{80} neon,^{64} and argon.^{66} Work in progress will significantly reduce the uncertainties from theory for neon and argon.^{326} As noted by Rourke,^{323} there are some discrepancies between the *ab initio* calculations of the susceptibilities and the experimental values often cited from Barter *et al.*;^{328} the discrepancies are many times larger than the stated uncertainties in the theoretical calculations. This may be due to errors in the 1930s-era argon data used in Barter’s calibration; error in the theoretical value seems unlikely at least for helium, where there is independent verification as discussed in Sec. 3.6. It is noteworthy that the large discrepancy between theory and Barter’s experiments is in the opposite direction for helium than it is for neon and argon, suggesting that Barter might have had an experimental problem specific to helium. A modern experimental determination of *A*_{μ} for helium and argon (perhaps involving measuring the ratio of the two) would be highly desirable. Even a 1% uncertainty for this measurement would be good enough to resolve the existing discrepancies, which are on the order of 7%.

First-principles calculations of *B*_{ɛ}(*T*) for helium have been available for a long time.^{284} Reference values from the latest pair potential and polarizability can be found in Ref. 79. These results have been independently confirmed (except at the lowest temperatures) by semiclassical calculations.^{329} Due to the recent development in three-body polarizabilities^{208} and dipole-moment surfaces,^{210}^{,} *ab initio* values of *C*_{ɛ}(*T*) with well-defined uncertainties are also available for both helium isotopes.^{210} These values agree with the limited experimental data available, but have much smaller uncertainties.

In the case of neon, the most accurate *ab initio* *B*_{ɛ}(*T*) has been computed by Hellmann and co-workers,^{117} who also reported well-characterized uncertainties. The results are in very good agreement with DCGT measurements. The third dielectric virial coefficient of neon is only approximately known from *ab initio* calculations, since the contributions from the three-body polarizability and dipole-moment surfaces can only be estimated with several uncontrolled approximations.^{207}

Regarding argon, the second dielectric virial coefficient has been computed using a fully *ab initio* procedure in Refs. 79, 164, and 329. Analogously to neon, the lack of *ab initio* three-body surfaces for the polarizability and dipole moment has prevented a fully first-principles calculation of *C*_{ɛ}(*T*) for argon. Approximate values were reported in Ref. 207.

Calculations of the second refractivity virial coefficient, *B*_{R}, for helium, neon, and argon were performed by Garberoglio and Harvey^{79} using the best pair potentials and Cauchy moments available at the time, although in many cases a rigorous uncertainty propagation was not possible. In the case of neon, the subsequent improved *B*_{ɛ} from Hellmann *et al.*^{117} can be combined with the frequency-dependent correction from Ref. 79 to provide improved values of *B*_{R}.

### 4.6. Transport properties

When the thermodynamic equilibrium of a gas is perturbed, dynamic processes will tend to restore it. The actual response depends on the specific kind of induced non-homogeneity: density variations will give rise to diffusive processes, relative motions will be damped by internal friction, and temperature gradients will result in heat flowing through the system.

The kinetic theory of gases^{330} provides a theoretical framework to analyze non-equilibrium behavior and transport properties of gases, determining how the flux of matter, momentum, or heat depends on the spatial variation of density, velocity, or temperature. The most accurate description is based on the Boltzmann equation, which describes the evolution of the state of a fluid where simultaneous interactions of three or more particles are neglected; hence, it is valid in the low-density regime only. Despite this limited scope, additional approximations are needed to make the kinetic equations manageable, for example by limiting the strength of the inhomogeneities to the linear or quadratic regime, which are situations that find widespread application.

*η*) and thermal conductivity (

*λ*) describe the linear relation between momentum and temperature inhomogeneities and the resulting internal friction and heat

*π*

_{ij}is the pressure tensor,

*p*the isotropic pressure, $u$ the macroscopic velocity,

**q**the heat flux, and

*T*the temperature. Kinetic theory shows how to compute

*η*and

*λ*from the details of the microscopic interaction between atoms. To this end, it is useful to define

*σ*(

*E*,

*θ*) is the differential cross section for two particles with energy

*E*in the scattering reference frame (

*E*=

*μv*

^{2}/2, where

*μ*=

*m*/2 is the reduced mass and

*v*the modulus of the relative velocity). The quantities defined by Eq. (102) are known as collision integrals. Equation (101) is valid when the cross section is calculated either in the classical or quantum regime; in the latter case one must further consider the fermionic or bosonic nature of the interacting atoms.

^{278}The viscosity and thermal conductivity are given by

*k*of the approximations involved, which in turn involve collision integrals of higher order. In the quantum case, collision integrals cannot be computed using path-integral MC methods, but their value depends on the scattering phase shift (see Sec. 4.2). For example, the expression for

*Q*

^{(2)}(

*E*) becomes

^{331},

*k*= 3 and

*k*= 5, respectively.

As pointed out in Sec. 2.5, the accuracy of *ab initio* calculations of transport properties for helium vastly exceeds that of experiments. We report in Fig. 9 the evolution of the relative uncertainty in the theoretical calculation of *η*_{He} in the past 20 years. The most recent theoretical values, which have an accuracy that is more than enough for several metrological applications, can be found in Ref. 10. It is worth noting that a more accurate pair potential has been published in the meantime,^{11} although no corresponding calculation of transport properties has yet been published.

In the case of neon, the best theoretical estimates of transport coefficients are given in Ref. 117, while for argon they can be found in Ref. 118. For both gases, the best experimental results are obtained from ratio measurements using the *ab initio* value of the viscosity or thermal conductivity of helium.

## 5. Molecular Systems

While the focus of this review is on noble gases, which are the fluids of choice for most *ab initio*-based primary temperature and pressure metrology, first-principles thermophysical properties for molecular species can also be of interest and make significant contributions. Three of the most promising areas are humidity metrology, low-pressure metrology, and atmospheric physics.

There are two main factors that make rigorous *ab initio* calculations of properties much more difficult for molecules than for monatomic species. The first is the increased dimensionality, where interactions depend not only on distance but on the relative orientations of the molecules. This not only complicates the development of potential-energy surfaces between molecules, but also makes the calculation of properties such as virial coefficients a sampling problem in many dimensions. Second, for rigorous calculations the internal degrees of freedom of the molecule must be considered, because properties of interest (such as the mean polarizability) depend on the molecular geometry and a distribution of geometries is sampled for each quantum state of the molecule. In some cases it may be adequate to assume a rigid molecule, but at a minimum an estimate of the uncertainty introduced by this assumption is needed, even though it might be difficult to compute.

In this section, we will describe the calculation of single-molecule quantities and quantities involving two or more molecules, along with their use to calculate properties of interest for metrology. Particular attention will be given to methods for addressing the challenges specific to molecular species. Finally, we will discuss some metrological applications that use properties of molecular species.

### 5.1. Single-molecule calculations

#### 5.1.1. Intramolecular potentials

In order to compute values of a property of a molecule averaged over nuclear motions, it is necessary to have a PES for the molecule. Such surfaces can be developed with *ab initio* calculations, and they can often be refined if accurate spectroscopic measurements are available. Development of the intramolecular potential is relatively straightforward for diatomic molecules such as H_{2}, N_{2}, and CO because the potential is one-dimensional, but the dimensionality and complexity increases quickly with the number of atoms. Surfaces of sufficiently high quality for most purposes have been developed for the triatomic molecules H_{2}O^{335} and CO_{2}.^{336} These intramolecular potential-energy surfaces are also needed in order to sample configurations when considering molecular flexibility for pair calculations as described in Sec. 5.2.2. Except for few-electron diatomic species and two-electron triatomics, pure *ab initio* surfaces are not accurate enough to provide rovibrational spectra competitive with experiments, and the most accurate molecular surfaces are always semiempirical.

#### 5.1.2. Electromagnetic properties

In contrast to noble gases, molecular species have multipole moments in the BO approximation (dipole, quadrupole, etc.). The most significant for metrology is the electric dipole moment. Rigorous *ab initio* calculation of the dipole moment for a molecule such as H_{2}O requires the development of a surface in which the dipole moment vector is given as a function of atomic coordinates, along with the single-molecule PES. The dipole moment for a given rovibrational state can then be computed as the expectation value averaged over the wave function of that state. Because the population of states changes with temperature, the average dipole moment will also change (slightly) with temperature; this has been analyzed for H_{2}O and its isotopologues by Garberoglio *et al.*^{337}

The polarizability is another important quantity, both in the static limit for capacitance-based metrology and at higher frequencies for metrology based on optical refractivity. Unlike a noble gas whose polarizability at a given frequency is a single number, the polarizability of a molecule is a tensor that reflects the variation with direction of the applied field and of the molecular axes. However, the quantity of interest for metrology is the mean polarizability, defined as 1/3 of the trace of the polarizability tensor.

Polarizability reflects the response of the electrons to an electric field. It can be computed *ab initio* in a relatively straightforward way. While for monatomic species (and homonuclear diatomic species) the electronic polarizability is the only contribution, more complicated molecules have an additional contribution in the static limit and at low frequencies; this is usually called the vibrational polarizability. It can be thought of as the electric field distorting the molecule (and therefore its charge distribution) by pushing the negatively and positively charged parts of the molecule in opposite directions.

The molecular dipole moment and polarizability are defined as the first- and second-order response to an externally applied electric field *E*_{0}, respectively. They can be computed by numerical differentiation of the molecular energy computed in the BO approximation as a function of *E*_{0}, or by perturbation theory. Although in principle these two approaches should give the same result, in practice some differences are observed. For atomic systems, the results from perturbation theory are found to be more accurate than numerical differentiation and are generally preferred.^{206} In the case of water, numerical differentiation is considered more accurate for dipole-moment calculations.^{338}

Once intramolecular potential-energy surfaces, polarizability surfaces, and dipole-moment surfaces are available, one can calculate the temperature-dependent electromagnetic response of a molecule, that is the first dielectric virial coefficient *A*_{ɛ} [see Eq. (8)], which is generally given by two contributions:^{337} the first is proportional to the rovibrational and thermal average of the electronic polarizability surface, while the second depends on the squared modulus of the transition matrix element of the dipole-moment surface. Additionally, one can separate the contribution from the dipole-moment transition matrix elements into those transitions where the vibrational state of the molecule changes and those for which the vibrational state of the molecule does not change, but the rotational state does: these two components of the dipole-moment contribution to the molecular polarizability are known as vibrational and rotational polarizabilities, respectively.^{339}

For small molecules (two or three atoms), one can solve directly the many-body Schrödinger equation for nuclear motion^{340} (e.g., using the efficient discrete-variable representation^{341} of the few-body Hamiltonian^{342}) and then perform the appropriate rovibrational and thermal averages to obtain *A*_{ɛ}. It has recently been shown that the path-integral approach outlined in Sec. 4 can be successfully used to compute the first dielectric virial coefficient of water.^{337} It can possibly be generalized to larger molecules, where the direct solution of the many-body Schrödinger equation becomes very demanding in terms of computational power.

In the case of water, computational results using the most accurate intramolecular potential-energy surface,^{335} polarizability surface,^{343} and dipole-moment surface^{338} are within 0.1% of the experimental value for the static dipole moment,^{344} although the theoretical surfaces for water do not yet have rigorously assigned uncertainties.

#### 5.1.3. Spectroscopy

It is now possible, especially for molecules containing only two or three atoms, to compute the positions and intensities of spectroscopic lines *ab initio*. The calculation of line positions requires only the single-molecule PES. The more important quantity for thermodynamic metrology, however, is the intensity of specific lines. This requires both the PES and a surface for the dipole moment as a function of the coordinates. Accurate *ab initio* dipole-moment surfaces have been developed for H_{2}O,^{338} CO_{2},^{345,346} and CO.^{347} The possible use in pressure metrology of intensities calculated from the surfaces for CO and CO_{2} will be discussed in Sec. 5.4.

### 5.2. Calculations for molecular clusters

#### 5.2.1. Interaction potentials

The development of interaction potentials for molecular gases is more difficult than for atomic ones due to the additional degrees of freedom, but much of the description in Sec. 3 is still applicable. A common approximation when developing intermolecular pair potentials is to treat the molecules as rigid rotors, which reduces the dimensionality considerably. For example, the PES of a pair of flexible water molecules has 12 degrees of freedom. By freezing the four OH bond lengths and the two HOH bond angles, only six degrees of freedom, usually taken to be the center-of-mass separation and five angles describing the mutual orientation, remain. To minimize the consequences of freezing the intramolecular degrees of freedom, the zero-point vibrationally averaged structures of the monomers are often used instead of the corresponding equilibrium structures.^{348,349}

However, even a six-dimensional dimer PES requires investigating thousands or even tens of thousands of pair configurations with high-level *ab initio* methods. As discussed in Sec. 3, the most commonly applied level of theory is CCSD(T)^{350} for molecular monomers; this method is usually applied with the frozen-core (FC) approximation. Such a level of theory was only the starting point in the schemes used to develop the most accurate pair potentials for the noble gases beyond helium. For the CCSD(T) method, the computational cost scales with the seventh power of the size of the molecules, and the scaling becomes even steeper for post-CCSD(T) methods.

In recent years, several intermolecular PESs have been developed that go beyond the CCSD(T)/FC level of electronic structure theory. The first step is to include all electrons in the calculations. Examples of all-electron (AE) surfaces are the flexible-monomer water dimer PES of Ref. 351 and the rigid-monomer ammonia dimer PES of Ref. 352. Also, post-CCSD(T)/AE terms were used in the H_{2}–CO flexible-monomer PESs starting in 2012.^{353,354} The T(Q) contributions were shown to have surprisingly large effects on the H_{2}–CO spectra.^{355}

Intermolecular pair potentials can be accurately represented analytically by a number of different base functional forms. Mimicking the anisotropy of the PES is most commonly achieved either by using spherical harmonics expansions or by placing interaction sites at different positions in the molecules, with each site in one molecule interacting with each site in the other molecule through an isotropic function. The site-site form is also often used for the empirical effective pair potentials commonly employed in MD and MC simulations of large molecular systems. The analytic functions used to represent high-dimensional *ab initio* PESs for pairs of small rigid molecules typically have a few tens up to a few hundred fit parameters.

Determination of these parameters, i.e., fitting a PES to a set of grid points in a dimer configurational space and the corresponding interaction energies, was until recently a major task taking often several months of human effort. This bottleneck has recently been removed by computer codes that perform such fitting automatically. In particular, the autoPES program^{351,356} can develop both rigid- and flexible-monomer fits at arbitrary level of electronic structure theory. The automation is complete: a user just inputs specifications of monomers and the program provides on output an analytic PES. This means that the program determines the set of grid points, runs electronic structure calculations for each point, and performs the fit. In addition to developing automation, the autoPES project introduced several improvements in the strategy of generating PESs. In particular, the large-*R* region of a PES is computed *ab initio* from the asymptotic expansion. Such expansion predicts interaction energies well down to *R* about two times larger than the van der Waals minimum distance. This means that no electronic structure calculations are needed in this region and autoPES can develop accurate PESs for dimers of few-atomic monomers using only about 1000 grid points, while most published work used dozens of thousands of points.

Accurate analytic rigid-rotor PESs exist for a large number of both like-species and unlike-species molecule pairs. For metrology, the most noteworthy of these are N_{2}–N_{2},^{357} CO_{2}–CO_{2},^{358,359} H_{2}O–CO_{2},^{360} H_{2}O–N_{2},^{361} and H_{2}O–O_{2}.^{362} Other accurate PESs of this type are: N_{2}–HF,^{363} H_{2}O–H_{2}O,^{351,364} (HF)_{2},^{365} and H_{2}–CO.^{353,355}

Many of these PESs (e.g., those from Refs. 357, 358, and 360–362) are based on nonrelativistic interaction energies corresponding to the frozen-core CCSD(T) level of theory in the CBS limit and are represented analytically by site-site potential functions, with each site-site interaction modeled by a modified Tang–Toennies type potential^{366} with an added Coulomb interaction term. In the case of the N_{2}–N_{2} PES,^{357} corrections to the interaction energies for post-CCSD(T), relativistic, and core-core and core-valence correlation effects were considered. Motivated by the availability of extremely accurate experimental data for the second virial coefficients of N_{2} and CO_{2}, the N_{2}–N_{2}^{357} and CO_{2}–CO_{2}^{358} PESs were additionally fine-tuned such that these data are almost perfectly matched by the values resulting from the PESs. The maximum well depths of the PESs were changed by the fine-tuning by less than 1%. Such fine-tuning does, however, mean that properties such as virial coefficients calculated from these tuned potentials cannot be considered to be truly from first principles for the purpose of metrology.

The second group of PESs listed above was also developed using either CCSD(T), with FC or AE, or SAPT. Post-CCSD(T) terms were considered in some cases, as already mentioned above. A range of different functional forms was used in the fitting; for larger monomers it was most often the site-site form.

While the error introduced by approximating molecules as rigid rotors is believed to be small for the molecules considered here, more rigorous calculations should include the intramolecular degrees of freedom; this has been done for example for the H_{2}–H_{2}, H_{2}–CO, and H_{2}O–H_{2}O potentials.^{351,353,354,367,368} There are several difficulties involved in the generation of fully flexible potentials. The first is the larger number of degrees of freedom. A system of *N* molecules approximated as rigid rotors can be described by *C*_{r} = 6*N* − 6 coordinates, while *C*_{f} = 3*nN* − 6 coordinates are necessary to fully describe a configuration of the same molecules if each of the monomers has *n* atoms. For sampling *c* configurations per degree of freedom, the number of calculations needed to explore the PES grows exponentially as $cCr|f$. In the case of, say, the water trimer (*N* = 3, *n* = 3), even assuming *c* = 3 one goes from 3^{12} ≈ 5 × 10^{5} configurations for rigid models to 3^{21} ≈ 10^{10} configurations for a fully flexible approach. The exponential increase of the number of configurations as a function of the number of degrees of freedom to be considered is sometimes called the dimensionality curse. Not all of these configurations are equally important and there is room for significant pruning and clever sampling strategies: one of the most useful starts from potentials developed for rigid molecules and enables the development of fully flexible versions optimizing the number of additional molecular configurations to be evaluated.^{369,370} More generally, even for a few degrees of freedom, the product of dimensions strategy leading to the *c*^{C} is the worst strategy to follow. Instead, one uses various types of guided MC generation of grid points. In particular, the statistically guided grid generation method of Ref. 371 reduces the number of points needed for a six-dimensional PES to about 300 (assuming the use of *ab initio* asymptotics). Another important issue regards the choice of a suitable form for the analytic potential and the fitting procedure. As in the case of rigid potentials, site-site interaction models (based on exponential functions at short range, inverse powers at long range, and Coulomb potentials) are commonly used for intermolecular flexible potentials. For the intramolecular interactions, Morse functions are often used but polynomial expansions work sufficiently well for molecules in their low-energy rovibrational state.^{351} Nevertheless, the dimensionality curse drastically limits the development of fully flexible potentials and for the time being only pair and three-body potentials involving diatomic and triatomic molecules (notably water^{364,372,373}) have been developed.

#### 5.2.2. Density virial coefficients

The calculation of density virial coefficients for molecular systems can be performed in a way very similar to that for noble gases. The main difference concerns the evaluation of the matrix elements of the free-molecule kinetic energy operator, that is the generalization of Eq. (74) which in turn depends on the specific degrees of freedom considered in the molecular model under consideration.

In the most general case, one considers the translational degrees of freedom of all the atoms in the molecule. Equation (74) remains the same (with the obvious modification of an atom-dependent mass *m*), but one needs an intramolecular potential to keep the molecule bound and, in general, a large number of beads, especially if light atoms (such as hydrogen or one of its isotopes) are to be considered. This approach allows flexibility effects to be fully accounted for and has been applied to investigate the second virial coefficient of hydrogen^{367} and water^{368} isotopologues. As one might expect, flexibility is more important at higher temperatures. On the other hand, this approach requires intramolecular and intermolecular potentials that depend on all the degrees of freedom, which in turn call for very demanding *ab initio* electronic structure calculations.

At sufficiently low temperatures, molecules occupy their vibrational ground state, and rigid-monomer models are expected to be quite (although not perfectly) accurate. In this case, a whole molecule is described as a rigid rotor, that is by three translational and three rotational degrees of freedom (2 in the case of linear molecules). The matrix elements of the kinetic energy operator are, in this case, more complicated than that in Eq. (74), but their expression has been worked out for both linear^{374} and non-linear^{375,376} rotors.

The rigid-rotor approximation of a molecular system is, in principle, an uncontrolled approximation and, consequently, cannot directly provide rigorous data for metrological applications. On the other hand, the associated uncertainties can be partially offset by the fact that potential-energy surfaces can be generated with higher accuracy than in the case of fully flexible models.^{261,377} Validation of the *ab initio* results with experimental data can be used to establish the temperature range in which a rigid model is valid, and provide useful estimates of virial coefficients where experimental data are lacking. Additionally, rigid models can be a stepping stone toward the more accurate fully flexible approaches.

Also, semiclassical approximations of density^{378} or dielectric virial coefficients^{285} for molecular systems are available. They are generally much easier to evaluate than by path-integral calculations, and are quite accurate in many cases.^{337,368,379}

#### 5.2.3. Dielectric and refractivity virial coefficients

The calculation of dielectric and refractivity virial coefficients for molecular species is much more difficult than for the monatomic systems discussed in Sec. 4.5.3. In addition to the increased dimensionality, the charge asymmetry creates additional polarization effects in interacting molecules. A complete treatment must therefore include the effect of the molecular interactions not only on the polarizability of the molecules, but also on their charge distribution. Because of this complexity, it seems unlikely that coefficients beyond the second virial will be calculated in the foreseeable future, and quantitatively accurate calculations with realistic uncertainty estimates may be limited to diatomic molecules such as N_{2} or H_{2}.

The only attempt at such calculations we are aware of for realistic (polarizable) molecular models is the work of Stone *et al.*,^{380} who calculated the second dielectric virial coefficient for several small molecules, including CO and H_{2}O. A recent experimental determination of the second dielectric virial coefficient for CO^{381} was in qualitative but not quantitative agreement with the prediction of Stone *et al.*

For rigorous metrology, it would be necessary to characterize the uncertainty of the surfaces describing the mutual polarization and pair polarizability of the molecules. The dimensionality, and therefore the complexity, of these calculations for a diatomic molecule like N_{2} would be similar to that for the three-body polarizability and dipole surfaces for monatomic gases.

#### 5.2.4. Molecular collisions

In some pressure-metrology applications near vacuum conditions, collision rates, which are related to collision integrals, are required. We already introduced collision integrals for atom–atom collisions in Sec. 4.6, but the concept can be generalized to include atom–molecule and molecule–molecule collisions, enabling the calculation of transport properties for dilute molecular gases. While the collision integrals for atom–atom collisions result in a classical treatment from the solution of the linearized Boltzmann equation and in the quantum-mechanical case from the solution of the linearized Uehling–Uhlenbeck equation,^{382} the corresponding classical and quantum-mechanical equations for collisions involving molecules are the linearized Curtiss–Kagan–Maksimov equation^{383–386} and the linearized Waldmann–Snider equation.^{387–390}

The formalism for the calculation of collision integrals involving molecules is much more complex than in the case of atom–atom collisions. Relations for classical collision integrals were derived by Curtiss for rigid linear molecules^{391} and extended to rigid nonlinear molecules by Dickinson *et al.*^{392} The quantum-mechanical calculation of collision integrals involving two molecules has rarely been attempted because of the mathematical complexity and large computational requirements, whereas atom–molecule collisions have been studied quantum-mechanically more often. For collisions between a helium atom and a nitrogen molecule, collision integrals were calculated both classically and quantum-mechanically.^{393,394} The comparison showed that quantum effects are small except at low temperatures. The degree to which the quantum nature of collisions can be neglected for pairs with larger expected quantum effects, such as H_{2}O–H_{2}O, remains an open question, but the agreement with experiment of classically calculated dilute-gas viscosities for H_{2}O^{395} suggests that the classical approximation is adequate for most purposes.

### 5.3. Humidity metrology

Much humidity metrology requires knowledge of humid air’s departure from ideal-gas behavior. Because the densities are low, this can be described by the virial expansion. The second virial coefficient of pure water has been calculated^{368} based on flexible *ab initio* pair potentials computed at a high level of theory.^{364,372,373} It is necessary to take the flexibility of the water molecule into account to obtain quantitative accuracy.^{368}

The most important contribution to the nonideality of humid air comes from the interaction second virial coefficient of water with air. While fairly accurate measurements of this quantity exist near ambient temperatures, it can now be computed with similar or better uncertainty by combining the cross second virial coefficients for water with the main components of dry air.^{396} Good quality pair potentials exist for water with argon,^{397} nitrogen,^{361} and oxygen,^{362} and these have been combined by Hellmann^{362} to produce accurate water–air second virial coefficients between 150 and 2000 K.

For humidity metrology at pressures significantly higher than atmospheric, corrections at the third virial coefficient level become significant. Only very limited data exist for the relevant third virial coefficients (water–water–air and water–air–air),^{398} so *ab initio* calculation of these quantities would be useful. This requires development of three-body potential-energy surfaces for systems such as H_{2}O–N_{2}–N_{2} and H_{2}O–H_{2}O–O_{2}. To our knowledge, no high-accuracy surfaces exist for these three-molecule systems, but their development should be feasible with current technology.

The same framework can be used for humidity metrology in other gases. Hygrometers are typically calibrated with air or nitrogen as the carrier gas, but some error will be introduced if the calibration is used in the measurement of moisture in a different gas. Calibrations can be adjusted if *ab initio* values of the cross second virial coefficient are known for water with the gas of interest. Such values have been developed for several important gases, such as carbon dioxide,^{360} methane,^{399} helium,^{400} and hydrogen.^{401}

Some emerging technologies for humidity metrology can be aided by *ab initio* property calculations. Instruments to measure humidity from the change in dielectric constant with water content of a gas^{402,403} require the first dielectric virial coefficient of water, which depends on its molecular polarizability and dipole moment. These quantities and their temperature dependence have been a subject of recent theoretical study.^{337}

Spectroscopic measurement of humidity has also been proposed;^{404} this requires the intensity of an absorption line for the water molecule. Thus far, work in this area has used measured line intensities due to their smaller uncertainty compared to *ab initio* values. The recent work of Rubin *et al.*^{405} demonstrated mutually consistent sub-percent accuracy for both experimental and theoretical intensities based on a semiempirical PES for an H_{2}O line, offering promise for the future use of calculated intensities to reduce the uncertainty of humidity metrology.

### 5.4. Pressure metrology

Molecular calculations are also promising for pressure metrology at low pressures.^{105} Refractivity-based pressure measurements using noble gases are discussed in Sec. 2.3. Some proposed approaches use ratios of the refractivity of a more refractive gas (such as nitrogen or argon) to that of helium. Use of nitrogen in these systems would be aided by good *ab initio* results for the polarizability of the N_{2} molecule and its second density and refractivity virial coefficients.

For low pressures, on the order of 1 Pa and below, absorption spectroscopy is a promising approach for pressure measurement. The absorption of a gas such as CO or CO_{2} can be used to measure low gas densities (from which the pressure is calculated by the ideal-gas law, perhaps with a second virial correction); this can be a primary pressure standard if the line intensity is calculated from semiempirical potential-energy and dipole-moment *ab initio* surfaces tuned to spectral data. Even if measured intensities are used, theoretical results are valuable to check their accuracy. For CO_{2}, uncertainty of intensity measurements and agreement between theory and experiment below 0.5% has been obtained.^{346,406} The simpler CO molecule is more amenable to accurate theoretical calculations; consistency between experimental and theoretical line intensities on the order of 0.1% has recently been achieved.^{347} In these calculations, the potential-energy curve was purely empirical, but the dipole-moment surface was obtained *ab initio*. An unresolved question in this work so far is the uncertainty of *ab initio* calculated line intensities, which must depend in a complex way on the uncertainties in the intramolecular potential and in the dipole-moment surface. Without reasonable estimates for the uncertainty of calculated intensities, the utility of this spectroscopic method for primary pressure standards is diminished.

For ultrahigh vacuum, gas densities can be measured based on the collision rate between the gas and a collection of trapped ultra-cold atoms. Both lithium and rubidium have been proposed as the trapped species.^{407–413} While in some implementations an apparatus constant is derived from measurements,^{409,410} it has recently been recognized^{414} that the proposed procedure introduces error when light species (such as Li and H_{2}) are involved in the collisions.

It is also possible to determine the relevant proportionality factor for the collision rate from first principles using collision cross sections calculated from *ab initio* pair potentials and quantum collision theory. These calculations have been performed for lithium with H_{2} (the most common gas in metallic vacuum systems) and He;^{415,416} *ab initio* calculations with rubidium are more challenging due to the large number of electrons. A recent paper has reported first-principles collision rate coefficients for both Rb and Li with noble gases, H_{2}, and N_{2}.^{417} It is also possible to measure the ratio of two collision pairs (for example, Rb–H_{2} versus Li–H_{2}) to obtain the coefficient for a system that is more difficult to calculate *ab initio*;^{407,414} in this approach a low uncertainty for the simpler-to-calculate system (that with fewer electrons) is essential.

### 5.5. Atmospheric physics

In atmospheric physics, the interaction of radiation with atmospheric gases, particularly H_{2}O and CO_{2}, has received increasing attention for climate studies; it is also important for Earth-based astronomy where the atmosphere is in the optical path. Scientists in these fields rely on line positions and intensities in the HITRAN database.^{418} Increasingly, *ab initio* calculations are being used to supplement experimental measurements for these quantities, as has recently been summarized for CO_{2}.^{419}

### 5.6. Transport properties

While transport properties of molecular gases are of little relevance in precision metrology, for the sake of completeness we mention briefly the current state of the art for pure molecular gases. Most of the transport property calculations for such gases performed so far are based on classically calculated collision integrals for rigid molecules using the formalism of Curtiss^{383} for linear molecules and of Dickinson *et al.*^{392} for nonlinear molecules (see Sec. 5.2.4).

A representative example of such calculations for gases consisting of small molecules other than H_{2} are the classical shear viscosity and thermal conductivity calculations of Hellmann and Vogel^{395} and Hellmann and Bich,^{420} respectively, for pure H_{2}O. The agreement with the best experimental data is within a few tenths of a percent for the viscosity and a few percent for the thermal conductivity. For both properties, these deviations correspond to the typical uncertainties of the best experimental data. The significant contribution to the thermal conductivity due to the transport of energy “stored” in the vibrational degrees of freedom, which is not directly accounted for by the classical rigid-rotor calculations, was estimated using a scheme that only requires knowledge of the ideal-gas heat capacity in addition to the rigid-rotor collision integrals.^{420} The main assumption in this scheme is that collisions that change the vibrational energy levels of the molecules are so rare that their effects on the collision integrals are negligible.

For pure H_{2}, classical calculations are not accurate enough even at ambient temperature. Fully quantum-mechanical calculations were performed by Mehl *et al.*^{421} using a spherically-averaged modification of a H_{2}–H_{2} PES,^{261} thus reducing the complexity of the collision calculations to that for monatomic gases. Despite this approximation, the calculated shear viscosity and thermal conductivity values for H_{2} agree very well with the best experimental data, particularly in the case of the viscosity where the agreement is within 0.1%.

## 6. Concluding Remarks and Future Perspectives

The outstanding progress achieved during the last three decades by the *ab initio* calculation of the thermophysical properties of pure fluids and mixtures has drastically reduced the uncertainty of the measurement of these properties and of the thermodynamic variables temperature, pressure, and composition.

For example, consider primary thermometry. *Ab initio* calculations directly contributed to the acoustic and dielectric determination of the value of the Boltzmann constant that is used in the new SI definition of the kelvin. The remarkably accurate theoretical calculations of the polarizability and the non-ideality of thermometric gases have also facilitated simplified measurement strategies and techniques.^{29,49,53,54} Consequently, new paths directly disseminating the thermodynamic temperature are now available at temperatures below 25 K, where the realization of ITS-90 is particularly complicated. Various methods of gas thermometry have determined *T* with uncertainties that are comparable to or even lower than the uncertainty of realizations of ITS-90.^{35,62,63} Improved theory has also suggested that primary CVGT could usefully be revisited, as discussed in Sec. 2.2.4.

In the near future, technical achievements will likely further reduce the uncertainty of measurements of the thermodynamic temperature and the thermophysical properties of gases. Efforts include: (1) improving the purity of the thermometric gases at their point of use, (2) implementing two-gas methods to reduce the uncertainties from compressibility of the apparatus, and (3) developing robust microphones (possibly based on optical interferometry) to facilitate cryogenic AGT. In the remainder of this section, we will summarize current limitations and describe some prospects for future contributions.

### 6.1. Current limitations of *ab initio* property calculations

As described in Sec. 3, *ab initio* calculations of properties for individual helium atoms and pairs of atoms have achieved extraordinarily small uncertainties. Even for three-body interactions, the potential energy is now known with small uncertainty, and good surfaces are available for the three-body polarizability and dipole moment. This enables accurate calculations, with no uncontrolled approximations, of the second and third density, acoustic, and dielectric virial coefficients. This high accuracy is due to the small number of electrons involved; electron correlation at the FCI level is still tractable for three helium atoms with a total of six electrons.

For DCGT and RIGT, it would be desirable to have similarly accurate properties for neon and argon, because their higher polarizability (and therefore stronger response) reduces the relative effect of other sources of uncertainty such as imperfect knowledge of the compressibility of the apparatus or the presence of impurities in the gas. Unfortunately, this level of accuracy for neon and argon is unlikely to be obtained in the foreseeable future. The neon atom has ten electrons, as many as five helium atoms, and argon has 18. While recent efforts have (at large computational expense) significantly reduced the uncertainty of single-atom and dimer quantities for neon and argon,^{64–66,117,118} they do not approach the levels of accuracy achieved for helium. For example, the relative uncertainty of the best calculation of the static polarizability of a neon atom^{65} is more than 100 times greater than that of a helium atom.^{192} Similarly, the relative uncertainty of the pair potential minimum energy is about 100 times larger for neon^{117} than for helium.^{11} Therefore, the relative uncertainties of calculated gas-phase thermophysical properties will be much higher for other gases than for helium. In such cases, the most accurate values of properties may be obtained by measuring ratios of properties relative to that of helium. This has already been done for the static polarizability of neon and argon^{63} and for the low-density viscosity of several gases.^{166–168}

Refractivity-based thermal metrology^{82,323} requires *A*_{R}, and preferably also *B*_{R} and *C*_{R}. At microwave frequencies, the static values (*A*_{ɛ}, *B*_{ɛ}, etc.) can be used. At optical frequencies, *A*_{R} has been computed at a state-of-the-art level for helium,^{84} neon,^{64} and argon.^{66}^{,} *B*_{R} has been computed at a state-of-the-art level for helium,^{79} but corresponding calculations for neon and argon rely on values for the Cauchy moment Δ*S*(−4) that could be significantly improved.

Even with state-of-the-art *ab initio* results, it seems likely that ratio measurements using helium, such as those of Egan *et al.* for *A*_{R},^{119} will produce lower uncertainties. To our knowledge, the theory for calculating *C*_{R} at optical frequencies is not available. Therefore, at the moment, it is necessary to take rather uncertain values from experiment or assume (based on the small difference between *B*_{R} and *B*_{ɛ}) that it is equal to *C*_{ɛ}.

As mentioned in Sec. 4.5.3 and also noted by Rourke,^{323} another issue for refractivity methods is the unclear situation surrounding the *A*_{μ} contribution. The best calculations of the magnetic susceptibility for helium,^{80} neon,^{64} and argon^{66} disagree with the old, sparse measurements of these quantities^{328} by amounts much larger than their stated uncertainties. Independent calculations of the magnetic susceptibility for one or more of these species would be helpful in assessing this discrepancy, but what is most needed is a modern measurement of the magnetic susceptibility of a noble gas (probably argon), either as an absolute measurement or as a ratio to a substance with a better-known magnetic susceptibility, such as liquid water.

To reach higher pressures with helium-based apparatus, it would be desirable to have reliable values, with uncertainties, for the fourth virial coefficient *D*(*T*). The most complete first-principles estimate so far^{126} used high-accuracy two-body and three-body potentials, but had a significant uncertainty component due to the unknown four-body potential. Accurate calculations of the nonadditive four-body potential for helium are feasible with modern methods. A four-body PES for helium, even if its relative uncertainty was as large as 10%, would allow reference-quality calculation of *D*(*T*) and enable improved metrology. The fitting of *ab initio* calculations to functional forms with many variables could, in this case, benefit from recent progress in machine-learning-based methods.^{422}

### 6.2. Molecular gases

Nitrogen is an attractive option for gas-based metrology due to its availability in high purity and its longstanding use in traditional apparatus such as piston gauges, but its lack of spherical symmetry and its internal degree of freedom add complication to *ab initio* calculation of its properties. The development of potential-energy surfaces for pair and three-body interactions for rigid molecular models is certainly feasible. This is also possible for flexible models, pending the difficulties discussed in Sec. 5.2.1. Once these surfaces are available, the methods for calculations of density virial coefficients have already been proven.^{261,368,377} (see Sec. 5.2.2). To the best of our knowledge, no fully *ab initio* calculation of dielectric virial coefficients for molecular systems has been performed. This task will require the development of the molecular interaction-induced polarizability function and dipole-moment function. The path-integral approach described in Sec. 4 can certainly be extended to compute these quantities as well as rigorously propagate their uncertainties.

### 6.3. Improved uncertainty estimations

As mentioned in Sec. 4.3, much progress has been made in estimating realistic uncertainties for density and dielectric virial coefficients. The old method of simply displacing the potentials in a “plus” and “minus” direction, while correct for one-dimensional integrations such as *B* and *B*_{ɛ}, is inefficient and can produce inaccurate results for higher coefficients. The functional differentiation approach discussed in Sec. 4.3 provides more rigorous results.

However, it is not entirely clear how to obtain uncertainties for acoustic virial coefficients, because they involve temperature derivatives of *B*(*T*) and *C*(*T*). The rigorous assignment of uncertainty to a derivative of a function computed from uncertain input is an unsolved problem as far as we are aware. Binosi *et al.*^{303} recently applied a statistical method (the Schlessinger Point Method) to the estimation of uncertainties for acoustic virial coefficients; this may provide a way forward.

A similar issue exists for the low-density transport properties. The very low uncertainty of the viscosity of helium shown in Fig. 9 near 40 K, obtained with the traditional method of “plus” and “minus” perturbations to the pair potential, is an artifact of competing effects on the collision integral of perturbations from different parts of the potential. While *B*(*T*), for example, exhibits monotonic behavior with respect to perturbations in the potential, that is not the case for the collision integrals used to compute transport properties, which can cause uncertainties to be artificially underestimated. This was recognized by Hellmann and co-workers, who created potentials perturbed in additional ways to provide a non-rigorous but reasonable estimation method for the uncertainty of low-density transport properties for krypton,^{149} xenon,^{423} and neon.^{117} Further analysis would be welcome to improve the rigor of uncertainty estimates for transport collision integrals.

### 6.4. Transport properties

In addition to the uncertainty issue just mentioned, we see two areas for improvement in the field of transport properties. The first concerns the density dependence beyond the low-density limiting values discussed in this work. As mentioned in Sec. 2.5, for flow metrology it would be desirable to know the viscosity with small and rigorous uncertainties not only at zero density, but at the real densities at which instruments are calibrated. The first correction should be a virial-like term linear in density, but the most successful theory so far^{171–173} relies on some simplifying assumptions. A more rigorous theory would be a significant advance. Even if the initial density dependence were only known for helium, that would enable better metrology for other gases because of the established methods for measuring viscosity ratios.

The second area is the transport properties of molecular species, such as N_{2} or H_{2}O. As mentioned in Sec. 4.6, classical collision integrals can be calculated for these species when they are modeled as rigid rotors. While it is believed that the errors introduced by the assumptions of classical dynamics and rigid molecules are small, it would be desirable to have verification from a more rigorous calculation. One might expect quantum effects to be significant for the dynamics of H_{2}O collisions, since they make a large contribution to *B*(*T*) for H_{2}O.^{367} Since fully quantum calculation of collision integrals is currently intractable for all but the simplest systems, the development of a viable “semiclassical” method for transport properties would be desirable. No such formulation exists to our knowledge.

### 6.5. Simulations of liquid helium

While we have focused on the gaseous systems where *ab initio* properties are already making major contributions to metrology, the thermophysical properties of condensed phases (particularly for helium) are also important in temperature metrology. For example, the vapor pressures of liquid ^{3}He and ^{4}He are part of the definition of ITS-90.^{6} With highly accurate two- and three-body potentials for helium (perhaps eventually supplemented by a four-body potential), high-accuracy simulation of thermodynamic properties of liquid helium may become feasible.

In fact, path-integral simulations of liquid ^{4}He can be performed without uncontrolled approximations,^{281} although, to the best of our knowledge, the most recent *ab initio* potentials have not yet been employed to compute any liquid helium property (e.g., the specific heat – and hence the vapor pressure, via the Clapeyron equation – or the temperature of superfluid transition). Consequently, the accuracy of first-principles many-body potentials for liquid ^{4}He is largely unknown. The use of three-body (or higher, when available) potentials would require considerable computational resources, as has been recently observed in simulations of liquid para-H_{2},^{424} but theoretical developments in efficient simulation methods for degenerate systems^{425} might pave the way for a fully *ab initio* calculation of the thermophysical properties of condensed ^{4}He.

In the case of fermionic systems such as ^{3}He, the path-integral approach suffers in principle from a “sign problem,”^{426} which generally requires some approximations and results in a large statistical uncertainty. However, two research groups have recently claimed to have overcome these limitations,^{427,428} which might result in accurate calculations of thermophysical properties in the liquid phase also for this isotope.

### 6.6. Reproducibility and validation

It is desirable for metrological standards to be based on multiple independent studies, so that they will not be distorted by a single unrecognized error. For example, for the recent redefinition of the SI in which several fundamental physical constants were assigned exact values, it was required that the value assigned to the Boltzmann constant be based on consistent results from at least two independent experiments using different techniques and meeting a low uncertainty threshold.^{28} Similarly, metrological application of the calculated results discussed in this Review would be on a firmer basis if there was independent confirmation of the results.

The danger of an unrecognized error in calculated quantities is not merely hypothetical. For several years, the “best” calculated values of *C* for ^{3}He were in error below about 4.5 K because the effects of nuclear spin on the quantum exchange contribution had been incorporated incorrectly; this was eventually recognized and corrected in Errata.^{126,127} An early quantum calculation of *B*_{ɛ} of argon^{429} disagreed with a later study,^{79} apparently because of inexact handling of resonance states in the earlier work. Ideally, there would be independent confirmation of all the results cited in Table 4 so that any errors could be detected.

One helpful step in this direction would be more complete documentation of calculations, including computer code, so that others can reproduce or check the work. It is common to provide computer code for potential-energy surfaces, but the calculation of virial coefficients has typically been performed with specialized software that is not public.

More important for metrology, however, would be independent verification of the calculated results. Conceptually, this has two parts: validation of the calculated quantities and surfaces described in Sec. 3 (potential-energy, polarizability, and dipole surfaces; atomic and magnetic polarizabilities) and validation of the calculation of virial coefficients from these quantities (described in Sec. 4).

Validation of calculated virial coefficients is probably the easier of the two parts, because it is typically less computationally demanding. This has been done for a few quantities; for example, two groups have performed fully quantum calculations (in one case neglecting exchange effects that become important below 7 K) of *C*^{127,319} and *D*^{126,319} for ^{4}He. Consistency checks can also be made by comparing different calculation methods, including classical and semiclassical approaches that should agree with the quantum calculations at high temperatures. The error in *B*_{ɛ} for argon mentioned above was detected by comparing phase-shift calculations to PIMC and semiclassical calculations, showing the value of multi-method comparisons.

The independent validation of calculated atomic quantities and intermolecular surfaces is more difficult, because these require large amounts of dedicated computer time. There have been a few cases where parallel efforts have produced independent, high-quality results; these include *A*_{ɛ} for neon^{64,65} and the three-body potential of argon.^{123,276} Some validation is also provided when the state of the art advances and new potentials are produced that agree with previous potentials (but have smaller uncertainties); this has been the case with the sequential development of pair potentials for helium (Sec. 3.3). In some cases, however, these are not truly independent verifications because they are developed by the same group and use many of the same methods. While it may be difficult to justify the extensive work required to independently confirm a state-of-the-art calculated surface, there would be value in performing spot checks of a few points. This would require developers of surfaces to make their calculated points available (or at least a subset of them), and also the multiple calculated quantities that typically contribute to each point.

We believe that more attention should be paid to the reproducibility and validation of the calculated results that are increasingly important in precision metrology. Work of this nature may not be very attractive to funding agencies (or graduate students), but it is needed for more confident use of gas-based metrology.

## Acknowledgments

We thank Mark McLinden, Patrick Egan, and Ian Bell of NIST for helpful comments, and Richard Rusby of NPL for valuable discussion regarding the CVGT technique. K.S. acknowledges support from the NSF Grant No. CHE-2154908.

We acknowledge support from *Real-K* Project No. 18SIB02, which has received funding from the EMPIR program co-financed by the Participating States and from the European Union’s Horizon 2020 research and innovation program.

## 7. Author Declarations

### 7.1 Conflict of Interest

The authors have no conflicts to disclose.

## Data Availability

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

### 9. Appendix: Formulae for the Third Acoustic Virial Coefficient, *γ*_{a}

*c*,

*c*

_{T}and

*c*

_{TT}are functions of the temperature

*T*as well as

*r*

_{12},

*r*

_{13}, and

*r*

_{23}through their dependence on

*U*

_{3}. Performing the substitution

*γ*

_{0}= 5/3, we obtain

*γ*

_{a}is more complicated, due to the fact that the ring-polymer distribution function

*F*of Eq. (77) depends on temperature. In particular, defining $U$ so that

*γ*

_{a}. However, this straightforward approach is characterized by large variance in the MC simulations, since Eq. (A9) has a form analogous to the

*thermodynamic*estimator of the kinetic energy.

^{301}It is possible to derive equivalent expressions with smaller variance, using the same ideas that lead to the

*virial*estimator.

^{301,430}The resulting formulas are very cumbersome, and can be found in Ref. 303.

## REFERENCES

*T*−

*T*

_{90}measurements over the temperature range from 430 K to 1358 K under the auspices of the EMPIR InK2 project

*Ab initio*calculations for helium: A standard for transport property measurements

*Ab initio*values of the thermophysical properties of helium as standards

^{4}He thermophysical properties: New

*ab initio*calculations

^{4}He and

^{3}He from an accurate relativistic interaction potential

^{4}He between 13.81 K and 287 K using a constant-volume gas thermometer

*ab initio*‘data’ to accurately determine the fourth density virial coefficient of helium

*h*,

*e*,

*k*, and

*N*

_{A}for the revision of the SI

*k*by acoustic thermometry of helium-4 gas

Note that in the literature one finds multiple and inconsistent definitions of the acoustic virial coefficients, depending on the variable chosen for the expansion of *w*^{2} (the pressure *p* or the molar density *ρ*) and the powers of *RT* included in the definition of the acoustic virials. We used the convention put forward in Ref. 282; in this case *β*_{a} has the same dimensions as the second virial coefficient *B* and *RTγ*_{a} has the same dimensions as the third virial coefficient *C*.