User:Airman72

For molecules that are non attracting hard spheres, $a=0$, the vdW virial expansion becomes simply $Z=(1-b\rho )^{-1}=\sum _{i=0}^{\infty }(b\rho )^{i}$, which illustrates the effect of the excluded volume alone. It was recognized early on that this was in error beginning with the term $(b\rho )^{2}$. Boltzmann calculated its correct value as $(5/8)(b\rho )^{2}$, and used the result to propose an enhanced version of the vdW equation

$$(p+a/v^{2})(v-b/3)=RT[1+2b/(3v)+7b^{2}/(24v^{2})].$$

On expanding $(v-b/3)^{-1}$, this produced the correct coefficients thru $(b/v)^{2}$ and also gave infinite pressure at $b/3$, which is approximately the close packing distance for hard spheres. This was one of the first of many equations of state proposed over the years that attempted to make quantitative improvements to the remarkably accurate explainations of real gase behavior produced by the vdW equation.

Derivation

As Goodstein has noted, "Many derivations, pseudo-derivations, and plausibility arguments have been given for this equation." The derivation presented here uses the canonical ensemble of statistical mechanics applied to a moderately dense gas. This approach makes explicit the assumptions required in order to obtain the van der Waals equation.

The canonical partition function is defined as

$${\cal {Z}}=\sum _{j}\,\exp \left(-{\frac {E_{j}}{kT}}\right)=\sum _{E_{j}}\rho _{E_{j}}\,\exp \left(-{\frac {E_{j}}{kT}}\right)$$

where the first sum is over all the states of the system while the second sum is over the discrete (quantized) energies, and $\rho _{E_{j}}$ is the number of states (degeneracy) with energy $E_{j}$. The function ${\cal {Z}}$ specifies the thermodynamic Helmholtz free energy function by $F(V,T)=-kT\ln[{\cal {Z}}(V,T)]$ (here $V$ is the variable that specifies the $E_{j}$), and through it all the other macroscopic thermodynamic functions,

$$p=-\partial _{V}F|_{T}\qquad S=\partial _{T}F|_{V}\qquad U=F+TS$$

On a macroscopic scale the discrete energy levels are closely spaced so little error is introduced by replacing them with a continuous variable and write

$${\cal {Z}}=\int \,\rho (E)\exp \left(-{\frac {E}{kT}}\right)\,dE.$$

The integration is taken over all energies with little additional error incurred because the integrand is sharply peaked about its average value. Furthermore on this scale the system energy, apart from its internal molecular structure, can be expressed in terms of the momenta and positions of the $N$ molecules that comprise it

$$E=\sum _{i=1}^{N}{\frac {{\bf {p}}_{i}^{2}}{2m}}+\Phi ({\bf {r^{N}}})$$

Here ${\bf {p}}_{i}$ is the vector momentum of the $i$th particle, and $\Phi ({\bf {r^{N}}})$ is the potential energy of the $N$ particles relative to one another where ${\bf {r^{N}}}$ is a shorthand way of designating the $N$ vector locations $({\bf {r}}_{1},{\bf {r}}_{2}\ldots {\bf {r}}_{N})$ of the particles. When the individual particles move in three dimensional space the system state is specified by $6N$ variables so ${\cal {Z}}$ can be represented by a point in this phase space whose elemental volume, $d{\bf {p^{N}}}d{\bf {r^{N}}}=\prod _{i=1}^{N}d{\bf {p}}_{i}d{\bf {r}}_{i}$, has a dimension which is a power of action, [Et]^3N=[mL²/t]^3N. However, in order to make this calculation of ${\cal {Z}}$ correspond to quantum statistical mechanical calculations, the elemental volume must be made dimensionless $d{\bf {p}}^{N}\,d{\bf {r}}^{N}/(N!h^{3N})$, where $h$ is Plank's constant. This results finally in the expression,

$${\cal {Z}}={\frac {1}{N!\,h^{3N}}}\int \int \,\exp \left(-\sum _{i=1}^{N}{\frac {{\bf {p}}_{i}^{2}}{2mkT}}\right)\exp \left[-{\frac {\Phi ({\bf {r^{N}}})}{kT}}\right]\,d{\bf {p^{N}}}d{\bf {r^{N}}}$$

Of the 6$N$ integrals in this expression, $3N$ of them, corresponding to the momenta, can be evaluated. Since

$$\exp \left(-\sum _{i=1}^{N}{\frac {{\bf {p}}_{i}^{2}}{2mkT}}\right)=\exp \left(-\sum _{i=1}^{3N}{\frac {p_{i}^{2}}{2mkT}}\right)=\prod _{i=1}^{3N}\exp \left(-{\frac {p_{i}^{2}}{2mkT}}\right)$$

they can all be written as a single definite integral with a well known value

$$\left[(2mkT)^{1/2}\int _{-\infty }^{\infty }\,e^{-x^{2}}\,dx\right]^{3N}=(2\pi mkT)^{3N/2}$$

so that ${\cal {Z}}$ can be written simply as

$${\cal {Z}}={\frac {Q_{N}}{N!\,\Lambda ^{3N}}}\qquad {\mbox{where}}\qquad Q_{N}=\int \,\exp \left[-{\frac {\Phi ({\bf {r^{N}}})}{kT}}\right]\,d{\bf {r^{N}}}$$

Here $\Lambda =[h^{2}/(2\pi mkT)]^{1/2}$ is the thermal de Broglie wavelength and $Q_{N}$ is the configuration integral. The limits of integration of these integrals are specified by the volume occupied by the molecules. Since $F=-kT\ln {\cal {Z}}=kT\{\ln[(N!\,\Lambda ^{3N})]-\ln Q_{N}\}$, and subsequently $p=-\partial _{V}F|_{T}=kT\partial _{V}\ln Q_{N}|_{T}$, it is the configuration integral alone that specifies the equation of state $p=p(V,T)$.

When the molecules do not interact, $\Phi =0$ so $Q_{N}=V^{N}$. Then taking its natural logarithm, differentiating with respect to $V$, and multiplying by $kT$, produces $p=NkT/V$, the ideal gas law.

At this point the system potential energy function, $\Phi ({\bf {r^{N}}})$, must be specified in order to evaluate the $3N$ integrals that make up $Q_{N}$. This is difficult to do in general, but on making the approximation of pairwise additivity, $\Phi$ takes the form

$$\Phi =\sum _{1\leq i<j\leq N}\,\varphi (r_{ij}).$$

When $r_{ij}=|{\bf {r}}_{j}-{\bf {r}}_{i}|$, the distance between the two molecules, this applies to symmetric molecules. More fundamentally this approximation neglects the effect on the force exerted by one molecule on another when another is brought into their vicinity. Although the error created by this and other similar neglects (more than one additional molecule) is unknown, it surely becomes smaller as the molecular density becomes smaller. Using pairwise additivity $Q_{N}$ becomes

$$Q_{N}=\int \,\exp \left[-{\frac {\sum _{1\leq i<j\leq N}\varphi (r_{ij})}{kT}}\right]\,d{\bf {r^{N}}}=\int \,\prod _{1\leq i<j\leq N}\exp \left[-{\frac {\varphi (r_{ij})}{kT}}\right]\,d{\bf {r^{N}}},$$

but the equality is only exactly true in the limit $\rho \rightarrow 0$.

Defining $f_{ij}=\exp[-\varphi (r_{ij})/(kT)]-1$ the integrand can be written as

$${\begin{array}{lll}\prod _{i<j}\exp[-\varphi (r_{ij})(kT)]&=&1+\\&&f_{12}+f_{13}+f_{23}+f_{14}+\cdots +\\&&f_{12}f_{23}+f_{13}f_{23}+f_{12}f_{24}+\cdots +\\&&{\mbox{etc.}}\end{array}}$$

where the second line is a sum of terms that denote an interaction of two molecules, the third line is the interaction among three molecules, and so on. For a dilute gas, $\rho$ small enough, only interactions between two molecules are important , and in this case the partition function simplifies to

$$Q_{N}=\int \,\left(1+\sum _{i<j}f_{ij}\,\right)\,d{\bf {r}}_{ij}d{\bf {r}}^{N-1}$$

Here the differential volume has been separated to emphasize that $f_{ij}$ is a function of ${\bf {r}}_{ij}$ only. Now all $3N$ integrals can be evaluated for the first element of the integrand, and $3(N-1)$ integrals can be evaluated for the second giving

$$Q_{N}=V^{N}+V^{N-1}\sum _{i<j}\int \,f_{ij}\,d{\bf {r}}_{ij}$$

The intermolecular potential is the same for all pairs of molecules, so since the first molecule can be chosen in any one of $N$ ways, while the second can then be chosen in only $(N-1)/2$ ways (since $j>i$), this becomes a single integral in which, because the molecules are spherical, the angular coordinates in $d{\bf {r}}=r^{2}\sin \theta drd\theta d\phi$ have also been integrated

$$Q_{N}=V^{N}\left[1+{\frac {N(N-1)}{2V}}(4\pi )\int _{0}^{\infty }\,f(r)r^{2}\,dr\right]$$

With $N\gg 1$, then $N(N-1)\approx N^{2}$, and this is written finally as

$$Q_{N}=V^{N}\left[1-{\frac {N^{2}B(T)}{N_{A}V}}\right]\quad {\mbox{where}}\quad B(T)=-2\pi N_{A}\int _{0}^{\infty }\,f(r)r^{2}\,dr$$

This form for $Q_{N}$ produces

$$p=kT\partial _{V}\left\{N\ln V+\ln \left[1-{\frac {N^{2}B(T)}{N_{A}V}}\right]\right\}$$

Now the infinite series $\ln(1+x)=\sum _{i=1}^{\infty }(-1)^{i+1}x^{i}/i$ converges for$-1<x\leq 1$ so for small enough molar density, $\rho =N/(N_{A}V)$, this becomes more simply

$$p=NkT\partial _{V}\left[\ln V-{\frac {NB(T)}{N_{A}V}}\right]$$

Here only the first term of the logarithmic series, linear in $\rho =N/(N_{A}V)$, has been written. All the remaining terms are $O(\rho ^{2})$ meaning they approach $0$ at the same rate as $\rho ^{2}$. They are not included because terms of this order have already been dropped, namely those that represent the interaction of three or more molecules. Carrying out the differentiation gives the pressure as,

$$p=NkT/V+N^{2}kTB(T)/(N_{A}V^{2})=\rho RT[1+\rho B(T)],$$

This is just two terms of a virial equation of state, and $B(T)$ is called the second virial coefficient. Retaining the $O(\rho ^{2})$ terms that were dropped would have produced the entire virial equation of state, and this shows that the $i$th term of the expansion contains $i$ molecule force interactions. For the dilute gas described here $\rho B(T)\ll 1$, and the higher order terms are negligible.

Recall that

$$B(T)=-2\pi N_{A}\int _{0}^{\infty }\,f(r)r^{2}\,dr$$

where $f(r)=\exp[-\varphi (r)/(kT)]-1$. A characteristic $\varphi (r)$ is shown in dimensionless form in the accompanying plot. It is positive for $r<\sigma$, and negative for $r>\sigma$ with minimum $\varphi (r_{0})=-\epsilon$ at some $r_{0}>\sigma$. Furthermore $\varphi$ increases so rapidly that whenever $r<\sigma$ then $\exp[-\varphi (r)/(kT)]\approx 0$. In addition for $\epsilon /(kT)\ll 1$, the normal case except when $T$ is near $0$, the exponential can be approximated for $r>\sigma$ by two terms of its power series expansion. In these circumstances $f(r)$ can be approximated as

$$f(r)=\left\{{\begin{array}{lll}-1&&r<\sigma \\-\varphi /(kT)=-\epsilon /(kT){\bar {\varphi }}(r)&&r>\sigma \end{array}}\right.$$

where ${\bar {\varphi }}=\varphi /\epsilon$ has the minimum value of $-1$. Then, on evaluating the integral between $0$ and $\sigma$, the second virial coefficient can be written as,

$$B(T)=2\pi \sigma ^{3}N_{A}/3-2\pi \sigma ^{3}N_{A}\epsilon /(kT)\int _{1}^{\infty }\,|{\bar {\varphi }}(x)|x^{2}\,dx$$

Now $2\pi \sigma ^{3}/3N_{A}=4[4\pi /3(\sigma /2)^{3}]N_{A}=b$ so $B(T)=b-a/(RT)$ where $a=IN_{A}\epsilon b$, and $I$ is a finite numerical factor depending on the dimensionless intermolecular potential function

$$I=3\int _{1}^{\infty }\,|{\bar {\varphi }}(x)|x^{2},dx$$

With these definitions the two term virial equation of state is

$$p=\rho RT[1+\rho b-/\rho a/(RT)]=\rho RT[1+b/v-a/(vRT)]$$

The Taylor expansion of $(1-\rho b)^{-1}$ is given by $(1-\rho b)^{-1}=1+\rho b+O(\rho ^{2}b^{2})$, so when $\rho b=b/v\ll 1$ the terms $O(\rho ^{2}b^{2})$ are ignorable, as has been done throughout this derivation. In that case the expression for $p$ can be written equivalently as,

$$p=\rho RT[(1+\rho b)^{-1}-\rho a/(RT)]=RT/(v-b)-a/v^{2},$$

which is the vdW equation,

According to this derivation the vdW equation is an equivalent of the two term virial equation of statistical mechanics when $\rho b\ll 1$. Consequently it has been shown to be valid in a region where the gas is dilute, $\rho B(T)\ll 1$, or specifically $\rho b=b/v\ll 1$ [the requirement $\rho a/(RT)=I\epsilon /(kT)(b/v)\ll 1$ is true whenever $b/v\ll 1$]. However,as Goodstein has noted, the most interesting behavior of the vdW equation occurs in the vicinity of the critical point where $\rho _{c}b=1/3$, namely in a region where its validity is questionable. Thus he wrote that "Obviously the value of the van der Waals equation rests principally on its empirical behavior rather than its theoretical foundations."

Yet this very remarkable empirical behavior ,which has been described in earlier sections of this article, provides irreplaceable insights; as Boltzmann noted, "...van der Waals has given us such a valuable tool that it would cost us much trouble to obtain by the subtlest deliberations a formula that would really be more useful than the one that van der Waals found by inspiration, as it were."

Mixtures

In 1890 van der Waals published an article that began the study of fluid mixtures. It was subsequently included as Part III of a later published version of his thesis. His essential idea was that in a binary mixture of vdw fluids described by the equations

$$p_{1}={\frac {RT}{v-b_{11}}}-{\frac {a_{11}}{v^{2}}}\quad {\mbox{and}}\quad p_{2}={\frac {RT}{v-b_{22}}}-{\frac {a_{22}}{v^{2}}}$$

the mixture is also a vdW fluid given by

$$p={\frac {RT}{v-b_{x}}}-{\frac {a_{x}}{v^{2}}}$$

where

$$a_{x}=a_{11}x_{1}^{2}+2a_{12}x_{1}x_{2}+a_{22}x_{2}^{2}\quad {\mbox{and}}\quad b_{x}=b_{11}x_{1}^{2}+2b_{12}x_{1}x_{2}+b_{22}x_{2}^{2}$$

Here $x_{1}=N_{1}/N$, and $x_{2}=N_{2}/N$, with $N=N_{1}+N_{2}$ (so that $x_{1}+x_{2}=1$) are the mole fractions of the two fluid substances. Adding the equations for the two fluids shows that $p\neq p_{1}+p_{2}$, although for $v$ sufficiently large $p\approx p_{1}+p_{2}$. The quadratic forms for $a_{x}$ and $b_{x}$ are a consequence of the forces between molecules. This was first shown by Lorentz, and was credited to him by van der Waals. The quantities $a_{11},\,a_{22}$ and $b_{11},\,b_{22}$ in these expressions characterize collisions between two molecules of the same fluid component while $a_{12}=a_{21}$ and $b_{12}=b_{21}$ represent collisions between one molecule of each of the two different component fluids. This idea of van der Waals was later called a one fluid model of mixture behavior.

For molecules that are hard spheres, $b_{12}=(b_{11}+b_{22})/2$, the arithmetic mean of $b_{11}$ and $b_{22}$. Substituting this into the quadratic form, and noting that $x_{1}+x_{2}=1$ then produces

$$b=b_{11}x_{1}+b_{22}x_{2}$$

Van der Waals wrote this relation, but did not make use of it initially. However, it has been used frequently in subsequent studies, and its use is said to produce good agreement with experimental results at high pressure.

Van der Waals also in this article made use of the Helmholtz Potential Minimum Principle, one of a number of alternate extremum princlple of thermodynamics. This one states that in a system in diathermal contact with a heat reservoir $T=T_{R}$, $DF=0$ and $D^{2}F>0$, namely at equilibrium the Helmholtz potential is a minimimum. This leads to a requirement $\partial ^{2}f/\partial v^{2}|_{T}=-\partial p/\partial v|_{T}>0$, which is the same as before. For stable states of a single fluid the curvature of the molar Helmholtz potential, $f=FN_{\mbox{A}}/N$, is positive.

For a single substance the definition of the molar Gibbs free energy can be written in the form $f=g-pv$. Thus when $p$ and $g$ are constant the function $f(v)$ represents a straight line with slope $-p$, and intercept $g$. For a given subcritical, $T_{R}$, with $p=p_{s}(T_{R})$ and a suitable value of $g$ this line is tangent to $f(T_{R},v)$ at the molar volume of each coexisting phase, saturated liquid, $v_{f}(T_{R})$, and saturated vapor, $v_{g}(T_{R})$. Furthermore, each of these points is characterized by the same value of $g$ as well as the same values of $p$ and $T_{R}$. These are the same three specifications for coexistence that were used previously.

Figure 8: The straight line (dotted-solid black) is tangent to the curve $f_{r}(0.875,v_{r})$ (solid-dashed green, dotted gray) at the two points $v_{rf}=0.576$ and $v_{rg}=2.71$. The slope of the straight line, given by $\partial _{v_{r}}f_{r}=-p_{c}v_{c}/(RT_{c})p_{rs}$, is $-0.215$ corresponding to $p_{rs}=0.5730$. All this is consistent with the data of the green curve, $T_{r}=7/8$, of Fig. 1. The intercept on the line is $g$, but its numerical value is arbitrary due to a constant of integration.

As depicted in Fig. 8, the region on the green curve $f(T_{R},v)$ for $v\leq v_{f}$ ($v_{f}$ is designated by the left green circle) is the liquid. As $v$ increases past $v_{f}$ the curvature of $f$ (proportional to $\partial _{v^{2}}^{2}f=-\partial _{v}p$) continually decreases. The point characterized by $\partial _{v^{2}}^{2}f=-\partial _{v}p=0$, is a spinodal point, and between these two points is the metastable superheated liquid. For further increases in $v$ the curvature decreases to a minimum then increases to another spinodal point; between these two spinodal points is the unstable region in which the fluid cannot exist in a homogeneous equilibrium state. With a further increase in $v$ the curvature increases to a maximum at $v_{g}$, where the slope is $p_{s}$; the region between this point and the second spinodal point is the metastable subcooled vapor. Finally, the region $v\geq v_{g}$ is the vapor. In this region the curvature continually decreases until it is zero at infinitely large $v$. The double tangent line is rendered solid between its saturated liquid and vapor values to indicate that states on it are stable, as opposed to the metastable and unstable states, above it (with higher Helmholtz free energy), but black, not green, to indicate that these states are heterogeneous, not homogeneous solutions of the vdW equation.

In the case of a vdW fluid the molar Helmholtz potential can be cast in the form

$$f_{r}=u_{r}-T_{r}s_{r}={\mbox{C}}_{u}+T_{r}\{c-{\mbox{C}}_{s}-\ln[T_{r}^{c}(3v_{r}-1)]\}-9/(8v_{r})$$

where $f_{r}=f/(RT_{c})$. A plot of this function for the isotherm $T_{r}=7/8$ is shown in Fig. 8 along with the line tangent to it at its two coexisting saturation points.

Van der Waals introduced this function because its properties could be easily extended to the binary fluid situation. In a binary mixture of vdW fluids the Helmholtz potential is a function of 2 variables, $f(T_{R},v,x)$, where $x$ is a composition variable, for example $x=x_{2}$ so $x_{1}=1-x$. In this case there are three stability conditions

$${\frac {\partial ^{2}f}{\partial v^{2}}}>0\qquad {\frac {\partial ^{2}f}{\partial x^{2}}}>0\qquad {\frac {\partial ^{2}f}{\partial v^{2}}}{\frac {\partial ^{2}f}{\partial x^{2}}}-\left({\frac {\partial ^{2}f}{\partial x\partial v}}\right)^{2}>0$$

and the Helmholtz potential is a surface (of physical interest in the region $0\leq x\leq 1$). The first two stability conditions show that the curvature in each of the directions $v$ and $x$ are both positive for stable states while the third condition states that stable states correspond to elliptic points on this surface; the limiting condition, $\partial _{x^{2}}^{2}f\partial _{v^{2}}^{2}f-(\partial _{xv}^{2}f)^{2}=0$, specifies the spinodal points.

For a binary mixture the Euler equation, can be written in the form

$$f=-pv+\mu _{1}x_{1}+\mu _{2}x_{2}=-pv+(\mu _{2}-\mu _{1})x+\mu _{1}$$

Here $\mu _{j}=\partial _{x_{j}}f$ are the molar chemical potentials of each substance, $j=1,2$. For $p$, $\mu _{1}$ and $\mu _{2}$, all constant this is the equation of a plane with slopes $-p$ in the $v$ direction, $\mu _{2}-\mu _{1}$ in the $x$ direction, and intercept $\mu _{1}$. As in the case of a single substance, here the plane can be made tangent to the surface, and the locus of the coexisting phase points forms a curve on it. The coexistence conditions are that the two phases have the same $T$, $p$, $\mu _{2}-\mu _{1}$, and $\mu _{1}$; the last two are equivalent to having the same $\mu _{1}$ and $\mu _{2}$ individually which are the Gibbs conditions for material equilibrium in this case. Although this case is similar to the previous one of a single component, here the geometry can be more complex. The surface can develop a wave in the $x$ direction as well as the one in the $v$ direction so there can be two liquid phases, that can be either miscible, or wholly or partially immiscible, as well as the vapor phase. Despite a great deal of work on this problem by van der Waals and his successors, work which produced a great deal of knowledge about fluid mixtures, complete solutions to the problem were only first obtained in 1968 when the availability of modern computers made computation of mathematical problems of this complexity feasible. But the results that were obtained were, in Rowlinson's words,

a spectacular vindication of the essential physical correctness of the ideas behind the van der Waals equation, for almost every kind of critical behavior found in practice can be reproduced by the calculations, and the range of parameters that correlate with the different kinds of behavior are intelligible in terms of the expected effects of size and energy.

In order to obtain these numerical results the values of the constants of the individual component fluids $a_{11},a_{22},b_{11},b_{22}$ must be known. In addition, the effect of collisions between molecules of the different components, given by $a_{12}$ and $b_{12}$, must also be specified. In he absence of experimental data, or computer modelling results to estimate their value, the empirical combining laws,

$$a_{12}=(a_{11}a_{22})^{1/2}\qquad {\mbox{and}}\qquad b_{12}^{1/3}=(b_{11}^{1/3}+b_{22}^{1/3})/2,$$

the geometric and algebraic means respectively can be used. These relations correspond to the empirical combining laws for the intermolecular force constants,

$$\epsilon _{12}=(\epsilon _{11}\epsilon _{22})^{1/2}\qquad {\mbox{and}}\qquad \sigma _{12}=(\sigma _{11}+\sigma _{22})/2,$$

the first of which follows from a simple interpretation of the dispersion forces in terms of polarizabilities of the individual molecules while the second is exact for rigid molecules. Then, generalizing for $n$ fluid components, and using these empirical combinig laws, the expressions for the material constants are:

$$a_{x}=\sum _{i=1}^{n}\sum _{j=1}^{n}(a_{ii}a_{jj})^{1/2}x_{i}x_{j}\quad {\mbox{which can be written as}}\quad a_{x}={\Big (}\sum _{i=1}^{n}a_{ii}^{1/2}x_{i}{\Big )}^{2}$$

$$b_{x}=\sum _{i=1}^{n}b_{ii}x_{i}$$

Using these expressions in the vdW equation is apparently helpful for divers, as well as being important for physicists, physical chemists, and chemical engineers in their study and management of the various phase equilibria and critical behavior observed in fluid mixtures.

Another method of specifying the vdW constants pioneered by W.B. Kay, and known as Kay's rule. specifies the effective critical temperature and pressure of the fluid mixture by

$$T_{cx}=\sum _{i=1}^{n}T_{ci}x_{i}\qquad {\mbox{and}}\qquad p_{cx}=\sum _{i=1}^{n}\,p_{ci}x_{i}$$

In terms of these quantities the vdW mixture constants are then,

$$a_{x}={\frac {27}{64}}{\frac {(R_{x}T_{cx})^{2}}{p_{cx}}}\qquad \qquad b_{x}={\frac {1}{8}}{\frac {R_{x}T_{cx}}{p_{cx}}}$$

where $R_{x}=\sum _{i=1}^{n}\,R_{i}x_{i}$. Kay used these specifications of the mixture critical constants as the basis for calculations of the thermodynamic properties of mixtures.

Kay's idea was adopted by T. W. Leland, who applied it to the molecular parameters, $\epsilon ,\sigma$, which are related to $a,b$ through $T_{c},p_{c}$ by $a\propto \epsilon \sigma ^{3}$ and $b\propto \sigma ^{3}$ (see the introduction to this article or the section on its derivation). Using these together with the quadratic form of $a,b$ for mixtures produces

$$\sigma _{x}^{3}=\sum _{i=i}^{n}\sum _{j=1}^{n}\,\sigma _{ij}^{3}x_{i}x_{j}\qquad {\mbox{and}}\qquad \epsilon _{x}=\left[\sum _{i=1}^{n}\sum _{j=1}^{n}\epsilon _{ij}\sigma _{ij}^{3}x_{i}x_{j}\right]\left[\sum _{i=i}^{n}\sum _{j=1}^{n}\,\sigma _{ij}^{3}x_{i}x_{j}\right]^{-1}$$

which is the van der Waals approximation expressed in terms of the intermolecular constants. This approximation, when compared with computer simulations for mixtures, were in good agreement for a range $1/2<(\sigma _{11}/\sigma _{22})^{3}<2$, namely for molecules of not too different diameters. In fact Rowlinson said of this approximation, "It was, and indeed still is, hard to improve on the original van der Waals recipe when expressed in [this] form...""

Mathematical details

In terms of the right ascension of the Sun, α, and that of a mean Sun moving uniformly along the celestial equator, α_M, the equation of time is defined as the difference, Δt = α_M - α. In this expression Δt is the time difference between apparent solar time (time measured by a sundial) and mean solar time (time measured by a mechanical clock). The left side of this equation is a time difference while the right side terms are angles; however, astronomers regard time and angle as quantities that are related by conversion factors such as; 2π radian = 360° = 1 day = 24 hour. The difference, Δt, is measureable because α can be measured and α_M, by definition, is a linear function of mean solar time.

The equation of time can be calculated based on Newton's theory of celestial motion in which the earth and sun describe elliptical orbits about their common mass center. In doing this it is usual to write α_M = 2πt/t_Y = Λ where

• t is dynamical time, the independent variable in the theory,
• t_Y is the length of time in a tropical year,
• Λ is the ecliptic longitude of a dynamical mean Sun; the angle from the vernal equinox to a fictitious Sun moving uniformly along the ecliptic and coinciding with the apparent sun at apoapsis and periapsis.

Substituting α_M into the equation of time, it becomes

$\Delta t=\Lambda -\alpha =M+\lambda _{p}-\alpha$

The new angles appearing here are:

• M is the mean anomaly; the angle from the periapsis to the dynamical mean Sun,
• λ_p = Λ - M = 4.9412 = 283.11° is the ecliptic longitude of the periapsis written with its value on 1 Jan 2010 at 12 noon.

However, the displayed equation is approximate; it is not accurate over very long times because it ignores the distinction between dynamical time and mean solar time. In addition, an elliptical orbit formulation ignores small perturbations due to the moon and other planets. Another complication is that the orbital parameter values change significantly over long times, for example λ_p increases by about 1.7 degrees per century. Consequently, calculating Δt using the displayed equation with constant orbital parameters produces accurate results only for sufficiently short times (decades). It is possible to write an expression for the equation of time that is valid for centuries, but it is necessarily much more complex.

In order to calculate α, and hence Δt, as a function of M, three additional angles are required; they are

• E the Sun's eccentric anomaly,
• ν the Sun's true anomaly,
• λ = ν + λ_p the Sun's true longitude on the ecliptic.

The celestial sphere and the Sun's elliptical orbit as seen by a geocentric observer looking normal to the ecliptic showing the 6 angles (M, λ_p, α, ν, λ, E) needed for the calculation of the equation of time. For the sake of clarity the drawings are not to scale.

All these angles are shown in the figure on the right, which shows the celestial sphere and the Sun's elliptical orbit seen from the Earth (the same as the Earth's orbit seen from the Sun). In this figure ε = 0.40907 = 23.438° is the obliquity, while e = [1 − (b/a)²]^1/2 = 0.016705 is the eccentricity of the ellipse.

Now given a value of 0≤M≤2π, one can calculate α(M) by means of the following procedure:

First, knowing M, calculate E from Kepler's equation

$M=E-e\sin E$

A numerical value can be obtained from an infinite series, graphical, or numerical methods. Alternatively, note that for e = 0, E = M, and for small e, by iteration, E ~ M + e sin M. This can be improved by iterating again, but for the small value of e that characterises the orbit this approximation is sufficient.

Next, knowing E, calculate the true anomaly ν from an elliptical orbit relation

$\nu =2\tan ^{-1}\left[{\sqrt {\frac {1+e}{1-e}}}\tan {\frac {E}{2}}\right]$

The correct branch of the multiple valued function tan⁻¹x to use is the one that makes ν a continuous function of E(M) starting from ν(E=0) = 0. Thus for 0≤ E < π use tan⁻¹x = Tan⁻¹x, and for π < E ≤ 2π use tan⁻¹x = Tan⁻¹x + π. At the specific value E = π for which the argument of tan is infinite, use ν = E. Here Tan⁻¹x is the principal branch, |Tan⁻¹x| < π/2; the function that is returned by calculators and computer applications. Alternatively, note that for e = 0, ν = E and for small e, from a one term Taylor expansion, ν ~ E+e sin E ~ M +2 e sin M.

Next knowing ν calculate λ from its definition above

$\lambda =\nu +\lambda _{p}$

The value of λ varies non-linearly with M because the orbit is elliptical, from the approximation for ν, λ ~ M + λ_p + 2 e sin M.

Next, knowing λ calculate α from a relation for the right triangle on the celestial sphere shown above

$\alpha =\tan ^{-1}[\cos \varepsilon \,\tan \lambda ]$

Like ν previously, here the correct branch of tan⁻¹x to use makes α a continuous function of λ(M) starting from α(λ=0)=0. Thus for (2k-1)π/2 < λ < (2k+1)π/2, use tan⁻¹x = Tan⁻¹x + kπ, while for the values λ = (2k+1)π/2 at which the argument of tan is infinite use α = λ. Since λ_p ≤ λ ≤ λ_p+ 2π when M varies from 0 to 2π, the values of k that are needed, with λ_p = 4.9412, are 2, 3, and 4. Although an approximate value for α can be obtained from a one term Taylor expansion like that for ν, it is more efficatious to use the equation sin(α - λ) = - tan²(ε/2) sin(α + λ). Note that for ε = 0, α = λ and for small ε, by iteration, α ~ λ - tan²(ε/2) sin 2λ ~ M + λ_p + 2e sin M - tan²(ε/2) sin(2M+2λ_p).

Finally, Δt can be calculated using the starting value of M and the calculated α(M). The result is usually given as either a set of tabular values, or a graph of Δt as a function of the number of days past periapsis, n, where 0≤n≤ 365.242 (365.242 is the number of days in a tropical year); so that

$M={\frac {2\pi \,n}{365.242}}$

Using the approximation for α(M), Δt can be written as a simple explicit expression, which is designated Δt_a because it is only an approximation.

$\Delta t_{a}=-2e\sin M+\tan ^{2}{\frac {\varepsilon }{2}}\,\sin(2M+2\lambda _{p})=[-7.657\sin M+9.862\sin(2M+3.599)]{\mbox{min}}$

The equation of time as calculated by the exact procedure for Δt described in the text and the asymptotic expression for Δt_a given there.

This equation was first derived by Milne, who wrote it in terms of Λ = M + λ_p. The numerical values written here result from using the orbital parameter values for e, ε, and λ_p given previously in this section. When evaluating the numerical expression for Δt_a as given above, a calculator must be in radian mode to obtain correct values. Note also that the date and time of periapsis (perihelion of the Earth orbit) varies from year to year; a table giving the connection can be found in perihelion.

A comparative plot of the two calculations is shown in the figure on the right. The approximate calculation is seen to be close to the exact one, the absolute error, Err = |(Δt− Δt_a)|, is less than 45 seconds throughout the year; its largest value is 44.8 sec and occurs on day 273. More accurate approximations can be obtained by retaining higher order terms , but they are necessarily more time consuming to evaluate. At some point it is simpler to just evaluate Δt, but Δt_a as written above is easy to evaluate, even with a calculator, and has a nice physical explanation as the sum of two terms, one due to obliquity and the other to eccentricity. This is not true either for Δt considered as a function of M or for higher order approximations of Δt_a.

References

• Burington R S 1949 Handbook of Mathematical Tables and Formulas (Sandusky, Ohio: Handbook Publishers)
• Duffett-Smith P 1988 Practical Astronomy with your Calculator Third Edition (Cambridge: Cambridge University Press)
• Heilbron J L 1999 The Sun in the Church, (Cambridge Mass: Harvard University Press|isbn=0-674-85433-0)
• Hinch E J 1991 Perturbation Methods, (Cambridge: Cambridge University Press)
• Hughes D W, et al. 1989, The Equation of Time, Monthly Notices of the Royal Astronomical Society 238 pp 1529-1535
• Meeus, J 1997 Mathematical Astronomy Morsels, (Richmond, Virginia: Willman-Bell)
• Milne R M 1921, Note on the Equation of Time, The Mathematical Gazette 10 (The Mathematical Association) pp 372–375
• Moulton F R 1970 An Introduction to Celestial Mechanics, Second Revised Edition, (New York: Dover)
• Muller M 1995, "Equation of Time - Problem in Astronomy", Acta Phys Pol A 88 Supplement, S-49.
• Roy A E 1978 Orbital Motion, (Adam Hilger|ISBN=0-85274-228-2)
• Whitman A M 2007, "A Simple Expression for the Equation of Time", Journal Of the North American Sundial Society 14 pp 29–33