Exact solution of the van der Waals model in the critical region

The celebrated van der Waals model describes simple fluids in the thermodynamic limit and predicts the existence of a critical point associated to the gas-liquid phase transition. However the behaviour of critical isotherms according to the equation of state, where a gas-liquid phase transition occurs, significantly departs from experimental observations. The correct critical isotherms are heuristically re-established via the Maxwell equal areas rule. A long standing open problem in mean field theory is concerned with the analytic description of van der Waals isotherms for a finite size system that is consistent, in the thermodynamic limit, with the Maxwell prescription. Inspired by the theory of nonlinear conservation laws, we propose a novel mean field approach, based on statistical mechanics, that allows to calculate the van der Waals partition function for a system of large but finite number of particles $N$. Our partition function naturally extends to the whole space of thermodynamic variables, reproduces, in the thermodynamic limit $N\to \infty$, the classical results outside the critical region and automatically encodes Maxwell's prescription. We show that isothermal curves evolve in the space of thermodynamic variables like nonlinear breaking waves and the criticality is explained as the mechanism of formation of a classical hydrodynamic shock.

In late 19th century, due to the outstanding contribution of its founding fathers, L. Boltzmann, J.C. Maxwell and J.W. Gibbs, the statistical thermodynamics has been successfully introduced as the general conceptual framework for understanding equilibrium thermodynamic phenomena by means of statistical mechanics [1,2].The early success of the kinetic theory of gases and the discovery of the mean field approach for the derivation of the classical van der Waals equation [3] nurtured the hope that an equally neat and clear description of critical phenomena would be as effective as it was away from the critical region.Although second order phase transitions, as for example the ergodicity breakdown of real gases at the triple point, have nowadays been completely framed into a rigorous scaffold [6] -also partially guided by techniques from the near field theory [4,5] -first order phase transitions, such as gas-liquid phase transitions, turned out to be more elusive.Indeed, despite the remarkable accuracy of the celebrated van der Waals equation for the description of real gases, the behavior predicted within the critical region, where a real gas turns into a liquid, significantly departs from the experimental observations.Fig. 1 shows the typical isothermal curves of a real gas: above the critical temperature (for T > T c ) the pressure decreases as a strictly monotonic function of the volume; the critical point corresponds to the temperature T = T c where the critical isotherm develops an inflection point; below the critical temperature (at T < T c ) real isotherms are constant within a certain volume interval in spite of the oscillating behavior predicted by the classical van der Waals equation.Interestingly, the behavior of real isothermal curves within the critical region turns out to be intimately connected to the theoretical one via the celebrated Maxwell rule stating that the constant pressure plateau is placed in such a way it cuts lobes of equal areas on the associated van der Waals isotherm.As it is well known, the Maxwell rule corresponds to the condition of thermodynamic equilibrium such that, below the critical temperature, the Gibbs free energy develops two minima of equal value [7].
The remarkable validity, although heuristic, of Maxwell's approach stimulated countless studies aimed at a rigorous statistical mechanical description of first order phase transitions as for instance in the works of Lebowitz and Penrose [8] and van Kampen [9], where large classes of pairwise interaction potentials for particles (continuous and hard-sphere-like respectively) are considered or the work by Griffiths [10] that focusses on the study of analyticity properties of thermodynamic functions.Alternative methods to analyze phase transitions have also been developed based on a macroscopic approach to thermodynamics.For instance, the Landau theory allows to construct suitable asymptotic expansions of the free energy in the order parameters to obtain information of the critical exponents in the vicinity of the critical point (see e.g.[11]); the Widom approach relies on the construction of effective free energy functions based on the analysis of its scaling properties [12].Further recent developments in this direction led to the formulation of the thermodynamic limit as the semiclassical limit of nonlinear conservation laws where phase transitions are associated to shock solutions of a hyperbolic nonlinear PDE, in the class of conservation laws [13][14][15][16].Such nonlinear PDEs can be obtained in mean field theories from the analysis of differential identities of the free energy as showed in [13,14,17] for the Curie-Weiss and the Sherrington-Kirkpatrick models, or from the analysis of thermodynamic Maxwell relations as showed in [15,18] for the van der Waals model.Both the microscopic statistical mechanical approach, via the study of correlation functions asymptotics, and the macroscopic thermodynamic approach, based on the expansion of the free energy in the vicinity of the critical point, reveal the intimate connection with the singularity and catastrophe theory -since the very first pioneering contributions by Arnold -and the Hopf bifurcation theory (see e.g.[19]).
Despite the numerous progresses made in understanding phase transitions in a variety of contexts, from thermodynamics to the classical and quantum field theory [20,21], complex and biological systems [22,23], and the discovery of their intrinsic universality, a global analyti- cal description of phase transitions for the van der Waals gas is still missing.In this Letter, inspired by the theory of nonlinear PDEs, in particular in the class of nonlinear conservation laws, we propose a novel approach that allows to construct the partition function of the van der Waals model for a system of finite number of particles N .
Based on the mean field assumption that configurations of equal volume are equally weighted (Isochoric Weights Thermodynamic Ansatz), we can write the general functional form of the partition function.For finite N , the partition function is, as expected, analytic in the space of thermodynamic variables but it develops a singularity in the thermodynamic limit N → ∞.We use the Laplace method for the asymptotic evaluation of the partition function for large N and require that above the critical point, where the Laplace integral admits one single critical point, the leading asymptotics of the volume expectation value satisfies the classical van der Waals equation, which is assumed to be accurate away from the critical region for the class of systems under consideration.Remarkably, this condition allows to fix uniquely the functional form of the partition function in such a way that the logarithm of the probability density gives the correct Gibbs free energy of the van der Waals model above the critical point.Finally, we prove that in the critical region, defined as the region in the space of thermodynamic variables where the Laplace integral admits multiple critical points, the leading asymptotics for the volume develops a discontinuous behavior, providing the exact analytical description of the first order phase transition.The model.Let us consider one mole of fluid of identical particles of mass m and volume v described by the Hamiltonian of the form where p l is the momentum of the particles, N is Avogadro's number, ψ(r l , r m ) is a two-body interaction potential, and the last term models the interaction with an external field where P > 0 is a real positive mean field coupling constant and the volume v(r 1 , . . ., r N ) is defined as the minimum convex hull associated to the configuration {r 1 , . . ., r N }.The partition function is given by the standard formula for a canonical ensemble Z = d N p i d N r i e −βH N where β = (K B T ) −1 with K B the Boltzmann constant and T the temperature.Integration over the moment variables p l gives that Z = (2πmK B T ) 3N/2 Z where We have introduced the rescaled variables t = 1/(N K B T ) and x = −P/(N K B T ) so that Z is defined for negative x and positive t only.Let us define the density of free energy α N = −N −1 log Z (the quantity −α N is sometimes referred to as mathematical pressure, see e.g.[24]).The expectation value of a given observable O is defined in the usual manner, i.e.
In particular, let us observe that v = −∂α N /∂x.Isochoric Weights Thermodynamic (IWT) ansatz.We now assume that, in the thermodynamic limit, for fixed values of the thermodynamic variables x and t, configurations of equal volume occur with the same probability density, so that there exists a probability measure µ(v) such that the partition function ( 2) is of the form Z = ∞ b dµ(v) where b denotes the hard core volume per particle.This assumption gives a nonlinear generalization of the standard mean field approximation introduced for the statistical mechanical derivation of the van der Waals equation of state (see e.g.[3]).We also note that, from a formal perspective, the IWT ansatz is equivalent to the request that the moments v n for the model (1) are such that the measure µ(v) is the solution to the Stieltjes moments problem [25], that is Expressing the differential as dµ(v) = µ (v)dv, the function µ (v) gives the weight associated to a given volume configuration that, for fixed values of x and t, is the same (isochoric) for all configurations of equal volume.As for the canonical ensemble the logarithm of the probability density in (2) is linear in the variables x and t (this ensures entropy maximization at equilibrium), we have that the probability density µ (v) is such that log µ (v) = N xv + 1  2 tφ(v) + σ(v) for certain functions φ(v) and σ(v).Hence, the partition function takes the form In the following we will prove that the functions φ(v) and σ(v) can be uniquely determined via the request that the expectation value v evaluated, away from the critical region, according to the partition function (3) satisfies, in the thermodynamic limit, the celebrated van der Waals equation The study of the connection between the form of the twobody interaction potential ψ(r l , r n ) and the mean field potential given in terms of the function φ(v) and σ(v) is beyond the scope of the present work.We should also stress however that, in the present context, this information is not essential to fix the functional form of ϕ(v) and σ(v).
Let us now proceed by evaluating the leading order asymptotics, for N → ∞, of the partition function (3) in the region of thermodynamic variables x and t where log µ (v) admits one single critical point.Laplace's formula gives where In particular, the formula (5) implies that α = lim N →∞ α N .Identifying the external field coupling constant P in the Hamiltonian (1) with the physical pressure in Eq. ( 4) and choosing φ(v) = −2a/ v and σ(v) = − log ( v − b), Eq. ( 6) coincides with the van der Waals equation (4).As mentioned above, the matching conditions with the van der Waals equation above the critical region fix the partition function (3).We note that α(v) = G/T , where G is the Gibbs free energy density of the van der Waals model.A direct calculation shows that the partition function so obtained satisfies the Klein-Gordon equation in the light-cone variables Let us also observe that the integral expression ( 7) can be explicitly evaluated at finite N for t = 0 and gives v (x, 0) = b − 1/x − 1/N x that coincides with the equation of state (4) in the limit T → ∞.
Using the definition v = N −1 Z −1 ∂Z/∂x, the Klein-Gordon equation (8) implies that the volume density satisfies the nonlinear viscous conservation law of the type studied in [26] and that is related to the viscous analog of the Camassa-Holm equation.In the thermodynamic limit, above the critical temperature where the gradient of v is bounded, the term of order O(N −1 ) in ( 9) is negligible and the volume density satisfies the Riemann-Hopf type equation ) whose solution develops a gradient catastrophe in finite "time" t.As illustrated in Fig. 2.a, the volume v evolves in the space of thermodynamic parameters just like a nonlinear hyperbolic wave and the gradient catastrophe is associated to the critical point Beyond the critical time t c , the physical solution develops a shock discontinuity, corresponding to a first order phase transition, whose position at fixed t > t c is determined by the equal area rule and its speed U is given the Rankine-Hugoniot condition where v l and v r are the limiting values of v respectively to the left and to the right of the jump.It was observed in [18] that the Rankine-Hugonoit condition is equivalent to the Clausius-Clapeyron equation implying that the shock speed is proportional to the latent heat associated to the first order phase transition and the trajectory of the shock is interpreted as the coexistence curve of the gas-liquid phase as shown in Fig. 2.b.Such connection between phase transitions and scalar shock waves was first observed in the context of magnetic models (see e.g.[16] and reference therein), and in the classical thermodynamic setting in [15,18] where the notion of universality has been also discussed.
The critical region.The subset of the space of thermodynamic variables x and t where the free energy α(v) admits multiple critical (stationary) points defines the critical region associated to the gas-liquid phase transition.In this case, the leading asymptotics at large N of the partition function (7) is given by the formula where v i (x, t) are the three distinct real solutions to the equation of state α (v i ) = 0. Hence, consistently with the classical description of the van der Waals phase transition, below the critical temperature the Gibbs free energy develops three stationary points, two of which are local minima.In the limit N → ∞ the leading contribution  to the partition function is given by the point of local minimum v m such that α(v m ) ≤ α(v i ), for all i = m.Hence, within the critical region, the solution is given by ) is a root of the equation of state α (v m ) = 0 such that the Gibbs free energy has the lowest local minimum.The subset of the (x, t)-plane, such that α takes two equal minima α(v i (x, t)) = α(v j (x, t)), represents the curve of resonance of the exponential contributions in ( 5) and identifies the shock line shown in Fig. 2.b.As already known from the theory of classical shocks for the viscous Burgers equation, such resonance condition is equivalent to the equal areas rule [27].
In Fig. 3 we plot the isothermal curves evaluated using the partition function (7) .As N increases the exact isothermal curves develop an inflection point and rapidly converge to the asymptotic behavior predicted by the Laplace formula.We should emphasize that the partition function (7) provides a global description of isothermal curves in the space of thermodynamic variables and the description of the phase transition is apparently accurate already for relatively small N 10 4 if compared with Avogadro's number.
Concluding remarks.In this letter, we have derived a global mean field partition function, at finite number of particles, for the van der Waals model providing a description of gas-liquid phase transitions as observed in experiments.By construction, the partition function (7) reproduces the van der Waals isotherms away from the critical region and automatically encodes the Maxwell rule.The approach hereby presented, based on the IWT ansatz, can be extended to include the larger class of models associated to partition functions of the form (3) which should give the natural as well as unique extension to the critical region of systems modeled in the thermodynamic limit by equations of state of the form (6).

FIG. 1 .
FIG. 1. Real gas isothermal curves: within the critical region, between the points A and B the behavior predicted by the van der Waals equation (dashed line) departs from the experimentally observed (solid line).The actual critical isotherms are constructed starting from the theoretical ones via the Maxwell equal areas rule.

FIG. 2 .
FIG. 2. a) van der Waals isothermal curves (for the choice of parameters a = 1 and b = 3) above and below the critical temperature Tc = (N KBtc) −1 .b) Shock trajectory (solid line) and critical sector (delimited by the dashes lines) associated to multivalued isotherms.

FIG. 3 .
FIG. 3. Solution for different values of N and comparison with the classical shock at t = 1.15tc for the choice of parameters a = 1 and b = 3.