Virialgleichungen

Im großkanonischen Ensemble lässt sich die Virialentwicklung herleiten. Die großkanonischen Zustandssumme ${\displaystyle \Xi =\Xi (T,V,\mu )}$ ist definiert als

${\displaystyle \Xi =\sum _{N=0}^{\infty }Z_{N}z^{N}=1+Z_{1}z+Z_{2}z^{2}+Z_{3}z^{3}+\dotsb }$

Dabei ist ${\displaystyle Z_{N}:=Z(T,V,N)}$ die kanonische Zustandssumme für ${\displaystyle N}$ Teilchen im Volumen ${\displaystyle V}$ bei Temperatur ${\displaystyle T}$, ${\displaystyle z:=e^{\beta \mu }}$ die Fugazität, ${\displaystyle \beta$ :=(kT)^{-1}} die inverse Temperatur und ${\displaystyle \mu }$ das chemische Potential. Das großkanonische Potential ${\displaystyle \Omega =-pV}$ ist mit dem Logarithmus der Zustandssumme verknüpft: ${\displaystyle \Omega =-kT\,\ln \Xi }$. Somit erhält man eine Zustandsgleichung als Potenzreihe in der Fugazität (${\displaystyle \ln \Xi }$ wird um ${\displaystyle z=0}$ entwickelt):

${\displaystyle {\begin{aligned}{\frac {pV}{kT}}&=-{\frac {\Omega }{kT}}=\ln \Xi =\ln \!\left(1+Z_{1}z+Z_{2}z^{2}+Z_{3}z^{3}+Z_{4}z^{4}+\dotsb \right)\\&=Z_{1}z+\left(Z_{2}-{\frac {1}{2}}Z_{1}^{2}\right)z^{2}+\left(Z_{3}-Z_{1}Z_{2}+{\frac {1}{3}}Z_{1}^{3}\right)z^{3}+\left(Z_{4}-Z_{1}Z_{3}-{\frac {1}{2}}Z_{2}^{2}+Z_{1}^{2}Z_{2}-{\frac {1}{4}}Z_{1}^{4}\right)z^{4}+\dotsb \end{aligned}}}$

Um diese Zustandsgleichung in eine thermische Zustandsgleichung – diese setzt die Zustandsgrößen Druck ${\displaystyle p}$, Volumen ${\displaystyle V}$, Temperatur ${\displaystyle T}$ und Teilchenzahl ${\displaystyle N}$ in Beziehung – umzuwandeln, muss die Fugazität ${\displaystyle z}$ durch die Teilchenzahl ${\displaystyle N}$ ausgedrückt werden.

Durch Ableiten des großkanonischen Potentials nach dem chemischen Potential ${\displaystyle \mu }$ erhält man die negative mittlere Teilchenzahl ${\displaystyle N}$:

${\displaystyle {\begin{aligned}N&=-\left.{\frac {\partial \Omega }{\partial \mu }}\right|_{T,V}=-{\frac {\partial z}{\partial \mu }}\left.{\frac {\partial \Omega }{\partial z}}\right|_{T,V}=-\beta z\left.{\frac {\partial (-\beta ^{-1}\ln \Xi )}{\partial z}}\right|_{T,V}=z{\frac {1}{\Xi }}\left.{\frac {\partial \Xi }{\partial z}}\right|_{T,V}\end{aligned}}}$

Damit folgt

${\displaystyle N\Xi =z\left.{\frac {\partial \Xi }{\partial z}}\right|_{T,V}\quad \Longleftrightarrow \quad N\left(1+Z_{1}z+Z_{2}z^{2}+Z_{3}z^{3}+\dotsb \right)=Z_{1}z+2Z_{2}z^{2}+3Z_{3}z^{3}+\dotsb }$

das mit dem Potenzreihenansatz ${\displaystyle z=\sum _{n=1}^{\infty }\alpha _{n}N^{n}}$ mittels Koeffizientenvergleich für Potenzen in ${\displaystyle N}$ gelöst werden kann. Die ersten Summanden sind

${\displaystyle {\begin{aligned}z=&{\frac {N}{Z_{1}}}+{\frac {N^{2}}{Z_{1}^{2}}}\left({\frac {Z_{1}^{2}-2Z_{2}}{Z_{1}}}\right)+{\frac {N^{3}}{Z_{1}^{3}}}\left({\frac {Z_{1}^{4}-5Z_{1}^{2}Z_{2}+8Z_{2}^{2}-3Z_{1}Z_{3}}{Z_{1}^{2}}}\right)+\\&+{\frac {N^{4}}{Z_{1}^{4}}}\left({\frac {Z_{1}^{6}-9Z_{1}^{4}Z_{2}-11Z_{1}^{3}Z_{3}+32Z_{1}^{2}Z_{2}^{2}+30Z_{1}Z_{2}Z_{3}-40Z_{2}^{3}-4Z_{1}^{2}Z_{4}}{Z_{1}^{3}}}\right)+\dotsb \end{aligned}}}$

Dies ergibt eingesetzt in obige Zustandsgleichung nach Ordnen von Potenzen von ${\displaystyle N}$ die thermische Zustandsgleichung

${\displaystyle {\begin{aligned}{\frac {pV}{kT}}=&N+N^{2}\left({\frac {1}{2}}-{\frac {Z_{2}}{Z_{1}^{2}}}\right)+N^{3}\left({\frac {1}{3}}-{\frac {2Z_{2}}{Z_{1}^{2}}}+{\frac {4Z_{2}^{2}}{Z_{1}^{4}}}-{\frac {2Z_{3}}{Z_{1}^{3}}}\right)+\\&+N^{4}\left({\frac {1}{4}}-{\frac {3Z_{2}}{Z_{1}^{2}}}+{\frac {27Z_{2}^{2}}{2Z_{1}^{4}}}-{\frac {20Z_{2}^{3}}{Z_{1}^{6}}}-{\frac {6Z_{3}}{Z_{1}^{3}}}+{\frac {18Z_{2}Z_{3}}{Z_{1}^{5}}}-{\frac {3Z_{4}}{Z_{1}^{4}}}\right)+\dotsb \end{aligned}}}$

Auf der rechten Seite lässt sich nach dem Äquipartitionstheorem die Ensemble-gemittelte Größe (Summe aller Ortskoordinaten mal Kraft) identifizieren:

${\displaystyle {\frac {pV}{kT}}=N-{\frac {1}{3}}\sum _{i=1}^{3N}\left\langle q_{i}{\frac {\partial (\Phi /kT)}{\partial q_{i}}}\right\rangle =N+{\frac {1}{3kT}}\sum _{i=1}^{3N}\left\langle q_{i}F_{i}\right\rangle =N+{\frac {1}{3kT}}\sum _{i=1}^{N}\left\langle {\boldsymbol {q}}_{i}\cdot {\boldsymbol {F}}_{i}\right\rangle }$

Führt man die Teilchenzahldichte ${\displaystyle \rho =N/V}$ ein, so erhält man die Virialentwicklung, d. h. eine Entwicklung der thermische Zustandsgleichung nach Potenzen der Teilchendichte:

${\displaystyle {\frac {p}{kT}}=\rho \underbrace {1} _{B_{1}}+\rho ^{2}\underbrace {V\left({\frac {1}{2}}-{\frac {Z_{2}}{Z_{1}^{2}}}\right)} _{B_{2}(T)}+\rho ^{3}\underbrace {V^{2}\left({\frac {1}{3}}-{\frac {2Z_{2}}{Z_{1}^{2}}}+{\frac {4Z_{2}^{2}}{Z_{1}^{4}}}-{\frac {2Z_{3}}{Z_{1}^{3}}}\right)} _{B_{3}(T)}+\dotsb }$

wobei die Virialkoeffizienten ${\displaystyle B_{i}(T)}$ identifiziert werden konnten. Im Grenzfall verschwindender Teilchendichte (${\displaystyle \rho \to 0}$) erhält man in führender Ordnung das ideale Gasgesetz

${\displaystyle {\frac {p}{kT}}=\rho \equiv {\frac {N}{V}}}$

Die kanonische Zustandssumme ${\displaystyle Z_{N}}$ für ${\displaystyle N}$ Teilchen ist im klassischen Grenzfall (hier vertauscht der kinetische mit dem potentiellen Energieterm) mit der Hamiltonfunktion ${\displaystyle H=\sum _{i=1}^{N}{\frac {{\boldsymbol {p}}_{i}^{2}}{2m}}+U_{N}}$ geben durch

${\displaystyle Z_{N}={\frac {1}{N!h^{3N}}}\underbrace {\int \!d^{3}{\boldsymbol {p}}_{1}\dotso d^{3}{\boldsymbol {p}}_{N}\,\exp \!\left(-\beta \sum _{i=1}^{N}{\frac {{\boldsymbol {p}}_{i}^{2}}{2m}}\right)} _{=(h/\Lambda )^{3N}}\underbrace {\int \!d^{3}{\boldsymbol {r}}_{1}\dotso d^{3}{\boldsymbol {r}}_{N}\,\exp \!\left(-\beta U_{N}\right)} _{=:Q_{N}}={\frac {Q_{N}}{N!\Lambda ^{3N}}}}$

wobei die Impulsintegrationen ausgeführt werden konnten (${\displaystyle \Lambda }$ ist die thermische Wellenlänge) und das Konfigurationsintegral ${\displaystyle Q_{N}}$ eingeführt wurde. Die gesamte potentielle Energie ${\displaystyle U_{N}}$ von ${\displaystyle N}$ Teilchen setzt sich aus einem externen Potential (${\displaystyle u_{1}}$) und internen Potentialen zwischen den Gasmolekülen, resultierend aus Zweiteilchen- und Mehrteilchenwechselwirkungen (${\displaystyle u_{2},u_{3},\dotsc }$), zusammen:

${\displaystyle U_{N}=U_{N}({\boldsymbol {r}}_{1},\dotsc ,{\boldsymbol {r}}_{N})=\sum _{i=1}^{N}u_{1}({\boldsymbol {r}}_{i})+\sum _{\underset {i<j}{i,j=1}}^{N}u_{2}({\boldsymbol {r}}_{i},{\boldsymbol {r}}_{j})+\sum _{\underset {i<j<k}{i,j,k=1}}^{N}u_{3}({\boldsymbol {r}}_{i},{\boldsymbol {r}}_{j},{\boldsymbol {r}}_{k})+\dotsb +u_{N}({\boldsymbol {r}}_{1},\dotsc ,{\boldsymbol {r}}_{N})}$

Mehrteilchenwechselwirkungen können z. B. durch Austauschwechselwirkung, durch Induktion und Dispersion (etwa dem Axilrod-Teller Dreifach-Dipol-Effekt) verursacht werden.

Insbesondere gilt für die Zustandssummen ${\displaystyle Z_{N}}$ mit den niedrigsten Teilchenzahlen ${\displaystyle N}$:

${\displaystyle Z_{0}=1\ ,\quad Z_{1}={\frac {Q_{1}}{\Lambda ^{3}}}\ ,\quad Z_{2}={\frac {Q_{2}}{2\Lambda ^{6}}}\ ,\quad Z_{3}={\frac {Q_{3}}{6\Lambda ^{9}}}}$

wobei die Konfigurationsintegrale wie folgt lauten (die Integrationen über die Ortskoordinaten erstrecken sich auf das zur Verfügung stehende Volumen ${\displaystyle V}$, somit ${\displaystyle \int \!d^{3}{\boldsymbol {r}}=V}$):

${\displaystyle {\begin{aligned}Q_{1}&=\int \!d^{3}{\boldsymbol {r}}_{1}\,e^{-\beta u_{1}({\boldsymbol {r}}_{1})}\\Q_{2}&=\int \!d^{3}{\boldsymbol {r}}_{1}\int \!d^{3}{\boldsymbol {r}}_{2}\,e^{-\beta \left[u_{1}({\boldsymbol {r}}_{1})+u_{1}({\boldsymbol {r}}_{2})+u_{2}({\boldsymbol {r}}_{1},{\boldsymbol {r}}_{2})\right]}\\Q_{3}&=\int \!d^{3}{\boldsymbol {r}}_{1}\int \!d^{3}{\boldsymbol {r}}_{2}\int \!d^{3}{\boldsymbol {r}}_{3}\,e^{-\beta \left[u_{1}({\boldsymbol {r}}_{1})+u_{1}({\boldsymbol {r}}_{2})+u_{1}({\boldsymbol {r}}_{3})+u_{2}({\boldsymbol {r}}_{1},{\boldsymbol {r}}_{2})+u_{2}({\boldsymbol {r}}_{1},{\boldsymbol {r}}_{3})+u_{2}({\boldsymbol {r}}_{2},{\boldsymbol {r}}_{3})+u_{3}({\boldsymbol {r}}_{1},{\boldsymbol {r}}_{2},{\boldsymbol {r}}_{3})\right]}\end{aligned}}}$

Damit schreibt sich die Virialentwicklung als

${\displaystyle {\frac {p}{kT}}=\rho \underbrace {1} _{B_{1}}+\rho ^{2}\underbrace {V\left({\frac {1}{2}}-{\frac {1}{2}}{\frac {Q_{2}}{Q_{1}^{2}}}\right)} _{B_{2}(T)}+\rho ^{3}\underbrace {V^{2}\left({\frac {1}{3}}-{\frac {Q_{2}}{Q_{1}^{2}}}+{\frac {Q_{2}^{2}}{Q_{1}^{4}}}-{\frac {1}{3}}{\frac {Q_{3}}{Q_{1}^{3}}}\right)} _{B_{3}(T)}+\dotsb }$

und die ersten Virialkoeffizienten lauten:

${\displaystyle {\begin{aligned}B_{1}&=1\\B_{2}&={\frac {V}{2}}\left(1-{\frac {2Z_{2}}{Z_{1}^{2}}}\right)={\frac {V}{2}}\left(1-{\frac {Q_{2}}{Q_{1}^{2}}}\right)\\B_{3}&=V^{2}\left({\frac {1}{3}}\left[1-{\frac {6Z_{3}}{Z_{1}^{3}}}\right]-{\frac {2Z_{2}}{Z_{1}^{2}}}\left[1-{\frac {2Z_{2}}{Z_{1}^{2}}}\right]\right)=V^{2}\left({\frac {1}{3}}\left[1-{\frac {Q_{3}}{Q_{1}^{3}}}\right]-{\frac {Q_{2}}{Q_{1}^{2}}}\left[1-{\frac {Q_{2}}{Q_{1}^{2}}}\right]\right)\end{aligned}}}$

Der erste Virialkoeffizient ist identisch 1, der zweite hängt von der Wechselwirkung mit einem äußeren Potential ${\displaystyle u_{1}}$ und Paarwechselwirkungen ${\displaystyle u_{2}}$, der dritte von der Wechselwirkung mit einem äußeren Potential, Paarwechselwirkungen und nicht-additiven 3-Teilchen-Wechselwirkungen ${\displaystyle u_{3}}$ ab, usw.

Der zweite Virialkoeffizient vereinfacht sich für den Spezialfall, dass kein äußeres Potential vorliegt (${\displaystyle u_{1}=0}$) und die Zweiteilchenwechselwirkungen nur vom Abstand ${\displaystyle |{\tilde {\boldsymbol {r}}}|=|{\boldsymbol {r}}_{1}-{\boldsymbol {r}}_{2}|}$ der Teilchen abhängen ${\displaystyle u_{2}({\boldsymbol {r}}_{1},{\boldsymbol {r}}_{2})=u_{2}(|{\tilde {\boldsymbol {r}}}|)}$ mit ${\displaystyle Q_{1}=\int d^{3}{\boldsymbol {r}}=V}$ und ${\displaystyle Q_{2}=V\int \!d^{3}{\tilde {\boldsymbol {r}}}\,e^{-\beta \,u_{2}(|{\tilde {\boldsymbol {r}}}|)}}$ zu

${\displaystyle {\begin{aligned}B_{2}&={\frac {1}{2}}\left[V-{\frac {V^{2}}{V^{2}}}\int \!d^{3}{\tilde {\boldsymbol {r}}}\,e^{-\beta \,u_{2}(|{\tilde {\boldsymbol {r}}}|)}\right]={\frac {1}{2}}\int \!d^{3}{\tilde {\boldsymbol {r}}}\,\underbrace {\left[1-e^{-\beta \,u_{2}(|{\tilde {\boldsymbol {r}}}|)}\right]} _{=:-f(|{\tilde {\boldsymbol {r}}}|)}\\&=-{\frac {1}{2}}\int \!d^{3}{\tilde {\boldsymbol {r}}}\,f(|{\tilde {\boldsymbol {r}}}|)=-{\frac {4\pi }{2}}\int _{0}^{\infty }d{\tilde {r}}\,{\tilde {r}}^{2}\,f({\tilde {r}})\end{aligned}}}$

Dabei wurde im zweiten Schritt die Mayer-Funktion ${\displaystyle f({\boldsymbol {r}}):=e^{-\beta \,u_{2}({\boldsymbol {r}})}-1}$ eingeführt und im dritten aufgrund der Symmetrie Kugelkoordinaten verwendet.

Ist zusätzlich ${\displaystyle u_{3}=0}$, so gilt ${\displaystyle Q_{3}=V\int \!d^{3}{\tilde {\boldsymbol {r}}}_{1}d^{3}{\tilde {\boldsymbol {r}}}_{2}\,e^{-\beta \,\left[u_{2}(|{\tilde {\boldsymbol {r}}}_{1}|)+u_{2}(|{\tilde {\boldsymbol {r}}}_{2}|)+u_{2}(|{\tilde {\boldsymbol {r}}}_{1}-{\tilde {\boldsymbol {r}}}_{2}|)\right]}}$ und der dritte Virialkoeffizient wird zu

${\displaystyle {\begin{aligned}B_{3}&=\int \!d^{3}{\tilde {\boldsymbol {r}}}_{1}\int \!d^{3}{\tilde {\boldsymbol {r}}}_{2}\,\left\{{\frac {1}{3}}\left(1-e^{-\beta \,\left[u_{2}(|{\tilde {\boldsymbol {r}}}_{1}|)+u_{2}(|{\tilde {\boldsymbol {r}}}_{2}|)+u_{2}(|{\tilde {\boldsymbol {r}}}_{1}-{\tilde {\boldsymbol {r}}}_{2}|)\right]}\right)-e^{-\beta \,u_{2}(|{\tilde {\boldsymbol {r}}}_{1}|)}\left(1-e^{-\beta \,u_{2}(|{\tilde {\boldsymbol {r}}}_{2}|)}\right)\right\}\\&=-{\frac {1}{3}}\int \!d^{3}{\tilde {\boldsymbol {r}}}_{1}\int \!d^{3}{\tilde {\boldsymbol {r}}}_{2}\ f(|{\tilde {\boldsymbol {r}}}_{1}|)\,f(|{\tilde {\boldsymbol {r}}}_{2}|)\,f(|{\tilde {\boldsymbol {r}}}_{1}-{\tilde {\boldsymbol {r}}}_{2}|)\\&=-{\frac {1}{3V}}\int d^{3}{\boldsymbol {r}}_{1}\int d^{3}{\boldsymbol {r}}_{2}\int d^{3}{\boldsymbol {r}}_{3}\,f_{1,2}f_{2,3}f_{3,1}\end{aligned}}}$

Im letzten Schritt wurde ${\displaystyle f_{i,j}=f(|{\boldsymbol {r}}_{i}-{\boldsymbol {r}}_{j}|)}$ verwendet.

Unter Vernachlässigung nicht-additiver Mehrteilchenwechselwirkungen (${\displaystyle u_{n\geq 3}=0}$) lautet der vierte Virialkoeffizient

${\displaystyle {\begin{aligned}B_{4}&=-{\frac {1}{8}}\int d^{3}{\boldsymbol {r}}_{2}\int d^{3}{\boldsymbol {r}}_{3}\int d^{3}{\boldsymbol {r}}_{4}\,f_{1,2}f_{2,3}f_{3,4}f_{4,1}{\Bigl (}3+3f_{1,3}+3f_{2,4}+f_{1,3}f_{2,4}{\Bigr )}\\&=-{\frac {1}{8V}}\int d^{3}{\boldsymbol {r}}_{1}\int d^{3}{\boldsymbol {r}}_{2}\int d^{3}{\boldsymbol {r}}_{3}\int d^{3}{\boldsymbol {r}}_{4}\,f_{1,2}f_{2,3}f_{3,4}f_{4,1}{\Bigl (}3+6f_{1,3}+f_{1,3}f_{2,4}{\Bigr )}\end{aligned}}}$

und der fünfte

${\displaystyle {\begin{aligned}B_{5}=&-{\frac {1}{30V}}\int d^{3}{\boldsymbol {r}}_{1}\int d^{3}{\boldsymbol {r}}_{2}\int d^{3}{\boldsymbol {r}}_{3}\int d^{3}{\boldsymbol {r}}_{4}\int d^{3}{\boldsymbol {r}}_{5}\,f_{1,2}f_{2,3}f_{3,4}f_{4,5}f_{5,1}\cdot \\&\cdot {\Bigl (}12+60f_{2,4}+10f_{1,4}f_{2,5}/f_{5,1}+60f_{1,3}f_{3,5}+30f_{1,4}f_{2,5}+10f_{1,4}f_{2,5}f_{2,4}/f_{5,1}+\\&\ \ +15f_{1,3}f_{3,5}f_{2,4}+30f_{1,4}f_{2,5}f_{2,4}+10f_{1,3}f_{3,5}f_{1,4}f_{2,5}+f_{1,3}f_{2,4}f_{3,5}f_{1,4}f_{2,5}{\Bigr )}\end{aligned}}}$

An diesem Beispiel sieht man, dass die Berechnung höherer Virialkoeffizienten die Auswertung komplexer Integrale erfordert. Die Integranden lassen sich jedoch systematisch mit Hilfe der Graphentheorie bestimmen. Vernachlässigt man nicht-additive Mehrteilchenwechselwirkungen (${\displaystyle u_{n\geq 3}=0}$), so lässen sich die Virialkoeffizienten graphentheoretisch nach Mayer (Clusterentwicklung) bestimmen:

${\displaystyle {\begin{aligned}B_{i+1}&=-{\frac {i}{i+1}}\beta _{i}=-{\frac {i}{i+1}}{\frac {1}{i!}}\left\{\beta _{i}\right\}_{\text{labelled}}=-{\frac {1}{i+1}}{\frac {1}{(i-1)!}}\left\{\beta _{i}\right\}_{\text{labelled}}\\&=-{\frac {i}{V}}\alpha _{i+1}=-{\frac {i}{V}}{\frac {1}{(i+1)!}}\left\{\alpha _{i+1}\right\}_{\text{labelled}}=-{\frac {1}{i+1}}{\frac {1}{(i-1)!}}{\frac {1}{V}}\left\{\alpha _{i+1}\right\}_{\text{labelled}}\end{aligned}}}$

für ${\displaystyle i\geq 1}$. Dabei sind die ${\displaystyle \beta _{i}}$ irreduzible Clusterintegrale über Graphen mit ${\displaystyle (i+1)}$ Knoten und der Mayer-Funktion als Bindung, wobei über ${\displaystyle i}$ Knoten integriert wird. Die ${\displaystyle \alpha _{i}}$ sind irreduzible Clusterintegrale über Graphen mit ${\displaystyle i}$ Knoten und der Mayer-Funktion als Bindung, wobei über alle ${\displaystyle i}$ Knoten integriert wird. Kann ein zusammenhängender Graph nicht durch Zerschneiden einer Kante in zwei unzusammenhängende Graphen geteilt werden, so heißt er irreduzibel.

Viele realistische Gase zeigen für kleine Molekülabstände eine starke Abstoßung (Pauli-Abstoßung bei Überlappung der Atomhüllen) und für große Abstände eine schwache Anziehung (etwa wie ${\displaystyle r^{-6}}$). Ein einfaches Zweiteilchenwechselwirkungs-Potential kann wie folgt modelliert werden:

${\displaystyle u_{2}(r)={\begin{cases}\infty &{\text{für}}\quad r<\sigma \\v(r)&{\text{für}}\quad r\geq \sigma \end{cases}}}$

Das Potential ist unendlich für Abstände kleiner als der Hartkugelradius ${\displaystyle \sigma =2r_{0}}$ und für Abstände größer als dieser leicht negativ: ${\displaystyle v(r)<0}$ und ${\displaystyle |v(r)|\ll kT}$. Somit lässt sich die Mayer-Funktion nähern als

${\displaystyle f(r)\approx {\begin{cases}-1&{\text{für}}\quad r<\sigma \\-{\frac {v(r)}{kT}}&{\text{für}}\quad r\geq \sigma \end{cases}}}$

und der zweite Virialkoeffizient ist

${\displaystyle B_{2}=-{\frac {4\pi }{2}}\int _{0}^{\infty }dr\,r^{2}\,f(r)=2\pi \underbrace {\int _{0}^{\sigma }dr\,r^{2}} _{\sigma ^{3}/3}-{\frac {2\pi }{kT}}\int _{\sigma }^{\infty }dr\,r^{2}\,\left[-v(r)\right]=:b-{\frac {a}{kT}}}$

Hierbei ist ${\displaystyle b}$ mit dem Eigenvolumen ${\displaystyle V_{0}}$ verknüpft

${\displaystyle b:={\frac {2\pi }{3}}\sigma ^{3}=4{\frac {4\pi }{3}}r_{0}^{3}=4V_{0}>0}$

und ${\displaystyle a}$ ist ein Maß für das mittlere äußere (attraktive) Potential

${\displaystyle a:=2\pi \int _{\sigma }^{\infty }dr\,r^{2}\,\underbrace {\left[-v(r)\right]} _{>0}>0}$

Somit ist in dieser Näherung ${\displaystyle B_{2}=b-a/(kT)}$ für kleine Temperaturen negativ und für große positiv mit dem Grenzwert ${\displaystyle B_{2}=b}$ für ${\displaystyle T\to \infty }$. Die Temperatur ${\displaystyle T=a/(bk)}$, bei der der zweite Virialkoeffizient verschwindet, heißt Boyle-Temperatur.

Die Virialentwicklung lautet also:

${\displaystyle {\frac {p}{kT}}=\rho +\rho ^{2}\left(b-{\frac {a}{kT}}\right)+{\mathcal {O}}(\rho ^{3})}$

Für geringe Dichte (${\displaystyle \rho }$ klein) betrachten wir nur die Entwicklung nur bis zur zweiten Ordnung. Umformen und Ausnutzen von ${\displaystyle (1-\rho b)^{-1}\approx 1+\rho b}$ für ${\displaystyle \rho b\ll 1}$ (das mittlere für jedes Teilchen zur Verfügung stehende Volumen ${\displaystyle \rho ^{-1}}$ ist viel größer als das Eigenvolumen ${\displaystyle b}$) liefert:

${\displaystyle p+a\rho ^{2}=kT\rho \left(1+\rho b\right)\approx kT{\frac {\rho }{1-\rho b}}=kT{\frac {N}{V-Nb}}}$

${\displaystyle \Longrightarrow \quad \left(p+a{\frac {N^{2}}{V^{2}}}\right)\left(V-Nb\right)\approx NkT}$

Letzte Gleichung ist die sog. Van-der-Waals-Gleichung.

Für ein System aus harten Kugeln mit Radius ${\displaystyle R}$ lautet das Wechselwirkungspotential bzw. die Mayer-Funktion (mit ${\displaystyle 2R=\sigma }$)

${\displaystyle u_{2}(r)={\begin{cases}\infty &{\text{für}}\quad r<2R\\0&\mathrm {f{\ddot {u}}r} \quad r\geq 2R\end{cases}}\quad \Rightarrow \quad f(r)={\begin{cases}-1&{\text{für}}\quad r<2R\\0&{\text{für}}\quad r\geq 2R\end{cases}}}$

Obwohl dieses System in der Natur nie realisiert wird, findet es oft Anwendung in der statistischen Mechanik, da die Struktur realer Flüssigkeiten hauptsächlich von abstoßenden Kräften bestimmt wird. Dieses Modell ist somit das einfachste Modell mit flüssigkeitsähnlichen Eigenschaften bei hohen Dichten und wird als Referenzsystem für Störungsrechnungen verwendet. Es zeigt bereits reichhaltige strukturelle und thermodynamische Eigenschaften, wie z. B. einen flüssig-fest Phasenübergang im Bereich 0,494 bis 0,545 der Packungsdichte. Die größtmögliche Packungsdichte ist bei 0,7405 erreicht (siehe dichteste Kugelpackung).

Die ersten vier Virialkoeffizienten können für das Hartkugelmodell analytisch berechnet werden (${\displaystyle V_{0}}$ bezeichne im Folgenden das Kugelvolumen ${\displaystyle V_{0}={\tfrac {4}{3}}\pi R^{3}}$). Der zweite Virialkoeffizient ergibt sich analog zum Van-der-Waals-Gas ohne attraktiven Teil zu ${\displaystyle B_{2}=4V_{0}}$. Der dritte lässt sich wie folgt berechnen:

${\displaystyle {\begin{aligned}B_{3}&=-{\frac {1}{3}}\int \!\mathrm {d} ^{3}{\boldsymbol {r}}\int \!\mathrm {d} ^{3}{\boldsymbol {r}}'\,f(|{\boldsymbol {r}}|)f(|{\boldsymbol {r}}'|)f(|{\boldsymbol {r}}-{\boldsymbol {r}}'|)\quad {\text{mit}}\quad {\boldsymbol {r}}_{2}:={\boldsymbol {r}}'-{\boldsymbol {r}}\,,\ {\boldsymbol {r}}_{1}:={\frac {{\boldsymbol {r}}+{\boldsymbol {r}}'}{2}}\\&=-{\frac {1}{3}}\int \!\mathrm {d} ^{3}{\boldsymbol {r}}_{2}\,f(|{\boldsymbol {r}}_{2}|)\int \!d^{3}{\boldsymbol {r}}_{1}\,f\left(\left|{\boldsymbol {r}}_{1}-{\frac {1}{2}}{\boldsymbol {r}}_{2}\right|\right)f\left(\left|{\boldsymbol {r}}_{1}+{\frac {1}{2}}{\boldsymbol {r}}_{2}\right|\right)\end{aligned}}}$

Das Integral über ${\displaystyle r_{1}}$ lässt sich geometrisch auswerten: Es entspricht dem Überlappvolumen zweier Kugeln mit Radius ${\displaystyle \sigma }$, deren Abstand ${\displaystyle r_{2}}$ beträgt. Dieses Volumen ist das doppelte Volumen einer Kugelkalotte ${\displaystyle \pi h^{2}(3\sigma -h)/3}$ mit Höhe ${\displaystyle h=\sigma -r_{2}/2}$. Die Höhe ist die Hälfte der „Überlappstrecke“ ${\displaystyle 2\sigma -r_{2}}$ zwischen den Kugelmittelpunkten. Anschließend lässt sich über ${\displaystyle r_{2}}$ in Kugelkoordinaten integrieren (das Minuszeichen verschwindet, da ${\displaystyle f=-1}$ innerhalb der Kugel um den Ursprung ist):

${\displaystyle B_{3}=-{\frac {1}{3}}\int _{0}^{\sigma }\!\mathrm {d} r_{2}\,4\pi r_{2}^{2}\,(-1){\frac {2}{3}}\pi \left(\sigma -{\frac {1}{2}}r_{2}\right)^{2}\left[3\sigma -\left(\sigma -{\frac {1}{2}}r_{2}\right)\right]={\frac {5}{18}}\pi ^{2}\sigma ^{6}=10V_{0}^{2}}$

Die analytische Berechnung des vierten Virialkoeffizienten erfordert eine langwierige Rechnung. Das Ergebnis lautet:

${\displaystyle B_{4}={\frac {2707\pi +(438{\sqrt {2}}-4131\arccos {\frac {1}{3}})}{70\pi }}\,V_{0}^{3}\approx 18{,}364768\,V_{0}^{3}}$

Die höheren Virialkoeffizienten wurden numerisch berechnet. Die ersten zehn Terme der Virialentwicklung lauten unter Verwendung der Packungsdichte ${\displaystyle \eta =\rho V_{0}}$:

${\displaystyle {\begin{aligned}{\frac {p}{kT}}=\rho &{\big [}1+4\,\eta +10\,\eta ^{2}+18{,}364768\,\eta ^{3}+28{,}22445\,\eta ^{4}+39{,}81545\,\eta ^{5}+\\&+53{,}3418\,\eta ^{6}+68{,}534\,\eta ^{7}+85{,}805\,\eta ^{8}+105{,}8\,\eta ^{9}+\dotsb {\big ]}\end{aligned}}}$

Die Virialkoeffizienten hängen hier nicht von der Temperatur ab. Das Hartkugelsystem ist wie das ideale Gas ein ‚athermales‘ System. Die kanonischen Zustandssummen und strukturelle Eigenschaften sind unabhängig von der Temperatur, sie hängen nur von der Packungsdichte ab. Die freie Energie kommt rein von der Entropie, nicht von Beiträgen potentieller Energie.

Nähert man die ersten fünf Terme durch die benachbarte ganze Zahl (4, 10, 18, 28, 40), so lassen sich diese durch eine einfache Vorschrift berechnen ${\displaystyle B_{i+1}=(i^{2}+3i)V_{0}^{i}}$ für ${\displaystyle i\geq 1}$. Verwendet man diese zur Berechnung aller Koeffizienten, so erhält man eine genäherte Zustandsgleichung für harte Kugeln:

${\displaystyle {\frac {p}{kT\rho }}\approx 1+\sum _{i=1}^{\infty }(i^{2}+3i)\eta ^{i}={\frac {1+\eta +\eta ^{2}-\eta ^{3}}{(1-\eta )^{3}}}}$

Im letzten Schritt wurden dabei die ersten beiden Ableitungen der geometrischen Reihe verwendet. Diese Zustandsgleichung nach Carnahan und Starling stimmt trotz der Näherung für höhere Virialkoeffizienten sehr gut mit Simulationsergebnissen für die flüssige Phase überein (größte Abweichungen ca. 1 %). Der Phasenübergang und die Zustandsgleichung der festen Phase ist in der Carnahan-Starling Näherung nicht enthalten.