# Thermodynamic phase transitions and shock singularities

Giuseppe De Nittis , Antonio Moro

## Abstract

We show that, under rather general assumptions on the form of the entropy function, the energy balance equation for a system in thermodynamic equilibrium is equivalent to a set of nonlinear equations of hydrodynamic type. This set of equations is integrable via the method of characteristics and it provides the equation of state for the gas. The shock wave catastrophe set identifies the phase transition. A family of explicitly solvable models of non-hydrodynamic type such as the classical plasma and the ideal Bose gas is also discussed.

## 1. Introduction

A physical system in thermodynamic equilibrium is specified by a certain number of thermodynamic variables and state functions. The occurrence of a phase transition can be interpreted as the fact that there exists a thermodynamic function possessing a certain number of singular (critical) points.

For instance, the energy balance equation for a system in thermodynamic equilibrium can be written as follows (Landau et al. 1980) 1.1 where E is total energy of the system, T is temperature, P is pressure, S is entropy, V is volume and {τj}j=1,…,m is a set of additional parameters that determine the state of the system (Landau et al. 1980; Callen 1985). For a gas mixture of m species, τj=Nj would be the number of particles of each species and the conjugated functions (Λj=μj) are interpreted as the associated chemical potentials. In some other cases, the pairs (τj,Λj) would play the role of a set of generalized fluxes such as mass, heat flux, together with their conjugated generalized potentials (Jou et al. 1993; Muller & Ruggeri 1993).

For the sake of simplicity, let us assume that the state of the system is specified by the thermodynamic variables V , T and P only, and no extra variables are involved. If, from microscopic considerations, one can deduce the explicit form of the Helmholtz free energy F=ETS, then the equation of state is given by the well-known formula (Landau et al. 1980) 1.2 The equation of state can be viewed as a stationary point of the Gibbs potential Φ=F+PV as a function of V . The state at constant pressure and temperature such that the Gibbs potential has a minimum is a state of stable equilibrium. The existence of degenerate critical points of Φ such that the second and possibly higher order derivatives of Φ w.r.t. V vanish detects a phase transition. From the point of view of the singularity theory, one would like to look at the Gibbs potential as a stable unfolding of the map VΦ(V,Tc,Pc), where Tc and Pc are suitable constant parameters at which Φ has a degenerate singularity (Lu 1976). The phase transition corresponds to the catastrophe set of this stable unfolding. The relevance of catastrophe theory in the description of both thermodynamic phase transition and breaking waves near the critical point is well known and there are chapters in a number of textbooks on the subject (Lu 1976; Poston & Stewart 1978; Gaite 1990). The present analysis shows that this intimate connection is more general and extends to regions of the space of variables that are far from the criticality.

A classical example of phase transition is the change of state from gas to liquid. Figure 1 shows an isothermal curve (solid line) of a real gas at a certain temperature T below the critical temperature Tc. The phase transition takes place between the points A and B, where the pressure remains constant as the volume increases. On the left, one observes the liquid state whereas the gas state is observed on the right. The ideal gas model turns out to be too simple to predict the occurrence of phase transitions as the isotherms do not possess singular points (figure 2). The van der Waals model is the simplest model that predicts a phase transition. The dashed line in figure 1 shows the behaviour within the phase transition region according to the van der Waals model. It corresponds to a metastable state and, generically, it is not observed in the experiments. The correct behaviour can be recovered by applying the so-called Maxwell principle (Callen 1985). At the phase transition, the pressure remains constant and its value is determined in such a way that the area below the upper portion of the van der Waals curve is equal to the area above the lower portion.

Figure 1.

Real gas isotherm (solid line) below the critical temperature Tc. The phase transition takes place between the points A and B where the pressure is constant. The dashed line describes a metastable state within the critical region AB as predicted by the van der Waals model. Outside the critical region the van der Waals curve describes with good accuracy several real gases.

Figure 2.

Ideal gas isotherms are given by a family of hyperbolas and there are no singular points at finite pressure and volume.

This description is apparently similar to the evolution of a shock wave after the interchange of the independent variable V with the dependent variable P. Figure 3 shows the evolution of the van der Waals curve (up to reflections on the plane) at different temperatures on a (P,V) diagram before and after the critical value Tc. In particular, at the critical temperature Tc, a gradient catastrophe occurs after which the profile of the van der Waals curve becomes multi-valued. This analogy suggests that the volume V as a function of pressure and temperature behaves as the solution to a hyperbolic partial differential equation (PDE) of the form (Whitham 1974) 1.3 The solution to this equation after the critical point has a jump, and consequently exists in the weak sense only. The position of the jump is prescribed by the shock fitting procedure (Whitham 1974) in such a way that the discontinuity cuts off lobes of equal area. This is nothing but the Maxwell principle mentioned above. In the present study, we show that the function V (T,P) satisfies an equation of the form (1.3) under the assumption that the entropy can be decomposed into the sum of a function of V plus a function of T. It turns out that several cases of physical interest belong to this class. Moreover, the above assumption on the functional form of the entropy can also be justified using the Boltzmann principle (Landau et al. 1980) for systems in the semiclassical regime.

Figure 3.

Evolution at different temperatures of a van der Waals curve (up to reflections on the plane) as a nonlinear wave solution to a hyperbolic partial differential equation. Beyond the critical temperature Tc, the solution becomes multi-valued (dashed line) and the shock takes place. The position of the jump can be recovered by the shock fitting procedure (Whitham 1974).

We would like to emphasize that the assumption of general separability of the entropy function covers a broad family of gas models and even solids (Landau et al. 1980). The further restriction to entropies that can be written as the sum of a function of V and T applies to semiclassical weakly interacting gases admitting a general virial expansion with linear coefficients in T. The advantage of this restriction is that such models can simultaneously be treated within the framework of a classical theory of integrable hyperbolic PDEs. The theory of hydrodynamic-type systems of equations of the form (1.3) is well established and provides us with a solid general mathematical framework where the present results are derived and could be generalized. The scalar equation (1.3) is integrable as it possesses infinitely many commuting flows and it is solvable by the classical method of characteristics that provides a local implicit representation of the solution via an algebraic relation of the form This relation turns out to be the equation of state of the thermodynamic system. We mention that the description of general composite thermodynamic systems involves multi-component systems of hydrodynamic type. In this case, the complete integrability requires that some additional and highly non-trivial conditions such as diagonalizability and semi-Hamiltonian property have to be satisfied (Tsarev theorem; Tsarev 1990). The study of such systems within a proper physical setting and the analysis of phase transitions is beyond the scope of the present study and would deserve a separate study.

This approach, based on the study of a set of nonlinear hyperbolic PDEs (conservation laws) for the thermodynamic variables, allows us to construct the equation of state from the macroscopic features of the system such as the knowledge of a certain number of isothermal or isobaric curves, provided the functional form of the entropy fulfils certain general assumptions. This point of view shows some analogies with the one recently proposed by Cooper & Russell (2011). The importance of the theory of PDEs for the study and the interpretation of thermodynamics is well established and goes back to the seminal work of Carathéodory (1909). More recent developments in this direction involve formulations in terms of contact geometry (Arnold 1990). Hence, it is not surprising to come across a possible relevance in thermodynamics of the theory of integrable PDEs.

The paper is organized as follows. In §2, we derive a set of hyperbolic PDEs of hydrodynamic type under the assumption that the entropy is separable as the sum of a function of V and a function of the set of variables (T,τ1,…,τm). The system of equations so obtained is completely integrable and its solution provides the equation of state for the system. The specific form of the entropy function and the arbitrary functions resulting from the integration can be determined on the macroscopic level in terms of a sufficient number of particular isobaric and/or isothermal curves. In §3, we apply the procedure described in §2 to the virial expansion and discuss the ideal and van der Waals gas as special cases. In §4, we discuss the interpretation of the phase transition as a singular sector of the system of hydrodynamic type. A family of explicitly solvable cases for a class of more general separable entropy functions is discussed in §5. Section 6 is devoted to some concluding remarks.

## 2. Equations of state

Let us consider the Gibbs free energy Φ=ETS+V P. Then the balance equation (1.1) takes the form 2.1 where S,V,Λ1,…,Λm are functions of the thermodynamic variables . Equation (2.1) implies 2.2 The system of equations above is compatible if and only if ΦC2(𝒟), i.e. the following Maxwell relations hold 2.3 We observe that, according to the second equation in (2.2), the jump in the volume function V =V (T,P), as in figure 3, can be viewed as a discontinuity of the first derivative of the Gibbs free energy w.r.t. the pressure.

Let us assume that the entropy function S and the potentials Λi depend on the pressure P through the volume only, i.e. 2.4 Equations (2.3), together with the assumptions (2.4), give the following set of equations 2.5 Introducing the notation and equations (2.5) are as follows 2.6 The system of equations (2.6) is compatible, i.e. ∂τiPV =∂PτiV , and it is solvable via the method of characteristics. The general solution is a scalar field V (T,P,τ1,…,τm) constant along the family of characteristic curves on the space of variables (T,P,τ1,…,τm). Characteristics are defined via the system 2.7 where s parametrizes the curves.

Under the extra assumption that the dependence of the entropy function is separable with respect to the volume V and the set of variables (T,τ1,…,τm) as follows the set of equations (2.6) takes the so-called hydrodynamic form 2.8 The characteristic speeds ϕ0,…,ϕm depend in fact on the variable V only. The general solution to the system (2.8) is locally given by the formula 2.9 where f(V) is an arbitrary function of its argument (for a review on the general theory of hydrodynamic type systems, see Dubrovin & Novikov (1989) and references therein). We refer to equation (2.9) as the equation of state. If the characteristic speeds ϕ1,…,ϕm are known, then the function f(V) can be specified either by fixing a particular isothermal curve (or isotherm), T=T0, τi=τ0i or an isobaric curve (or isobar) P=P0, τi=τ0i in the space of thermodynamic variables (V,P,T,τ1,…,τm). We note that assigning the volume on a particular isobaric curve is equivalent to assigning an initial value problem for the system (2.8). Alternatively, if an isothermal curve is known instead of an isobaric one, the equation of state is recovered by solving the same system (2.8) written in evolutionary form with respect to the set of variables T, τ1,…,τm. For instance, choosing the initial datum V 0(P)=V (T0,P,τ01,…,τ0m), on the isothermal curve T=T0, τ1=τ01,…,τm=τ0m, we have, at least locally, where is the inversion of the initial datum V 0(P). The relation above specifies the form of the function f(V) in terms of the initial datum. We also note that the isochoric hypersurfaces V =const. on the m+2 dimensional space 𝒟 are, by definition, the characteristics of the system (2.8). We also note that if, apart from f(V), a certain number of characteristic speeds is not known, these can be obtained by measuring the same number of additional isotherms and/or isobars. For the sake of simplicity, assume that all the functions ϕi, i=0,1,…,m have to be determined. Given m+2 distinct isotherms in the form V i(P)=V (Ti,P,τi1,…,τim), one can locally compute the set of inverse functions . Then, the functions f(V) and ϕi(V), i=0,1,…,m are given by the solution to the following linear system: and

For the particular choice ϕj(V)=cjV j+1 with j=0,…,m, where cjs are arbitrary constants, the system (2.9) coincides with the first m+1 members of the Burgers–Hopf hierarchy (Zakharov 1980) 2.10 In this case, the equation of state is 2.11 The whole hierarchy is simply obtained by taking .

## 3. Virial expansion

In this section, we discuss the class of equations of state in the form of the virial expansion. For the sake of simplicity, we discuss first the ideal gas example and then the van der Waals model as a particular case of the more general virial expansion in the case where no extra variables τj are involved. Applying the procedure outlined above we show how, given either two isobars or equivalently two isotherms, one can recover the form of the characteristic speeds and uniquely determine the equation of state as a particular solution to the equation of hydrodynamic type.

### (a) Ideal gas

Let us consider a monoatomic gas of constant number of particles N and equation of state Assume that we are given two particular isobaric curves V 1(T)=V (T,P1) and V 2(T)=V (T,P2) with P1P2 of the form 3.1 where c1 and c2 are two given constants. The equation of state along these isobars can be written as where are obtained by inverting (3.1). Using equations (3.1), we get the following system of linear equations for the functions ϕ0(V) and f(V) 3.2 Subtracting the two equations (3.2), we get 3.3 where is a constant that has to be independent on the particular pair of the chosen isobars. Equation (3.3) implies that the entropy is a logarithmic function of the volume V . We will show, at the end of this section, how this property can be heuristically understood in terms of microscopic considerations. Moreover, one can see that the constant α should be proportional to the total number of particles N, in such a way that α=nR, where n=N/NA is the molar number, NA is the Avogadro number and R is the universal gas constant. Using equation (3.3) in (3.2), one obtains 3.4 The constants α and s specify the physical properties of the gas. The physical constraint implies that s=0 or equivalently This uniquely fixes the solution and provides the equation of state for the ideal gas 3.5 One can easily check that, proceeding as shown in the previous section, the same result can be obtained by fixing two isotherms of the form 3.6

### (b) Real gas

Let us assume that two particular isotherms on the (V,P) plane have the following form 3.7 for an arbitrary positive integer k. In the relation (3.7), it is required that where ν is interpreted as the intrinsic volume occupied by the gas molecules. We observe that the physical conditions are satisfied. The particular choice gives the ideal gas isotherms (3.6). Proceeding as in the previous example, solving the system of linear equations together with the assumption that ϕ0(V) is a function of the volume only, we have that the virial coefficients αi(Tj) are at most linear functions, i.e. and 3.8 3.9 The equation of state reads as 3.10

### (c) van der Waals approximation

Virial expansion of the order k=2 with and provides the van der Waals equation of state This is a particular implicit solution to the equation 3.11 that is equivalent to the Burgers–Hopf equation up to a Galilean change of variables.

### Remark 3.1 (entropy function of a semiclassical gas)

In the examples discussed above, the assumption of separability of the entropy function w.r.t. the variables T and V together with the two particular isobars (or isotherms) of the form (3.7) implies that the entropy has to be logarithmic in the volume V . This result can be heuristically understood, taking into account the definition of the entropy function for a semiclassical gas (Boltzmann principle; Landau et al. 1980) 3.12 where is the statistical weight, N is the number of particles and r is the number of degrees of freedom of a particle. ΔΓ counts the number of distinct microscopic states within the allowed region Δq Δp of the phase space.

If the gas is assumed to be confined in a macroscopic volume V and occupies an intrinsic volume ν, then Δq=(Vν)N. Moreover, as the temperature is a measure of the mean kinetic energy of the system, we can assume in the first instance that where ψ0 is a certain function. Hence, we end up with the following form of the entropy function 3.13 where κ is identified with the Boltzmann constant κ=R/NA and NA is Avogadro's number. Finally, one can reasonably assume that νN as in the van der Waals approximation. In particular, one can set ν=(N/NA)b, where b measures the physical volume occupied by a gas mole. As the ideal gas is made of point-like particles, it is b=0.

## 4. Singular sector and phase transitions

Let us introduce the function 4.1 Hence, the equation locally gives the solution V (T,P,τ1,…,τm) to the system (2.8), provided that the invertibility condition ∂Ω/∂V ≠0 is satisfied. In the following, we will adopt the notation Following Kodama & Konopelchenko (2002), let us introduce the singular sector of co-dimension sm+2 for the system of equations (2.8) 4.2 where and |𝒰j| is the number of solutions considered in 𝒰j. Let us introduce the definition of critical sector that will be useful in the following.

### Definition 4.1 (critical sector)

The critical sector of a system of equations as in (2.8) is the singular sector of highest co-dimension. The critical sector is said to be maximal if its co-dimension is m+2.

The critical sector identifies the occurrence of the phase transition as well as the shock formation of the solution to the system (2.8).

### (a) van der Waals case

For a van der Waals gas, we have ϕ1=ϕ2=⋯=ϕm≡0, and the system of equations (2.8) reduces to equation (3.11). In this case, we have 4.3 where (a polynomial function of V) is defined as follows: The singular sector Z1 is specified by the conditions These conditions give the set 4.4 where V ±(T,P) are the real roots of the quadratic polynomial , i.e.

As shown in figure 4, the singular sector Z1 is a collection of two curves emerging from the critical point The singular sector Z2 is given by the set of equations Note that . More explicitly, one has 4.5 Hence, Z2 is the critical sector of a first-order phase transition for the van der Waals gas. The corresponding critical value of the volume is In this case, the critical sector Z2 is maximal. As shown in figure 5, the critical isotherm has an inflection point owing to the fact that the first non-vanishing derivative is odd or, equivalently, the co-dimension of the critical sector is even.

Figure 4.

Singular sector Z1 for the hydrogen gas according to the van der Waals model: a=2.476×10−2 m6 Pa mol−2, b=2.661× 10−5 m3 mol−1, n=103, R=8.3144 J K mol−1. The intersection of the two lines is the critical sector Z2={Tc=33.159 K, Pc=1.29508×106 Pa}.

Figure 5.

The bold line shows the critical isotherm for the hydrogen gas according to the van der Waals model: a=2.476×10−2 m6 Pa mol−2, b=2.661×10−5 m3 mol−1, n=103, R=8.3144 J K mol−1. The critical point is V c=0.07983 m3, Tc=33.159 K, Pc=1.29508×106 Pa. Below the critical isotherm, there are two stationary points where the isotherm has a local maximum and a local minimum. On the (T,P) plane, stationary points describe the singular sector Z1 in figure 4.

### (b) Deformed virial expansion

Let us consider the function 4.6 The state of the system is specified by the set of variables (T,P,V,τ1,…,τm). Let us define Note that μm+2 as for k>m+2 the system of equations is overdetermined. For V ∈𝒰μ, the system of equations 4.7 can be viewed as a linear system in the set of m+2 variables (T,P,τ1,…,τm).

If μ<m+2, solving the system (4.7) with respect to the set of variables (T,P,τ1,…,τμ−1), then one can parametrize the critical sector Zμ as follows In particular, for μ=m+1, Zμ is a curve in the m+2 dimensional affine space (T,P,τ1,…,τm).

If μ=m+2, then the critical sector is maximal. It is given by a set of isolated points where V k is the kth root of As a particular example, let us consider the following deformed virial expansion 4.8 where γj are arbitrary real constants. The rank of the system (4.7) is maximal and the system is compatible. For the particular choice one obtains a deformed van der Waals model. Let us analyse the one parameter deformation m=1. In this case, μ=2, the conditions Ω=Ω′=Ω′′=0 give 4.9 and The critical sector Z2 is the curve (4.9) in the parameter space (T,P,τ1). Figure 6 shows the critical sector for the hydrogen gas. The plain circle on the curve corresponds to the critical volume V c=3nb for the van der Waals model such that τ1=0.

Figure 6.

Critical sector for the deformed van der Waals model in the case m=1 with a=2.476×10−2 m6 Pa mol−2, b=2.661× 10−5 m−3 mol−1 and n=103, R=8.3144 J K mol−1. Rescaled variables , and are used. The critical point is Tc=33.159 K, Pc=1.29508×106 Pa, τc=0 K Pa J m−3 and corresponds to the critical volume V c=0.07983 m3.

### Remark 4.2 (unfoldings)

Locally, in the vicinity of the critical point, the van der Waals equation is associated with the universal unfolding of co-dimension two of a Riemann–Hugoniot catastrophe (Lu 1976). More precisely, the universal unfolding gives the Gibbs free energy where F=ETS. The deformed virial expansion discussed above can be interpreted as a higher co-dimension unfolding of the van der Waals model.

## 5. Separable entropy functions

Let us consider a system the state of which is completely specified by the set of thermodynamic variables (T,P,V). Assume that we are given a more general separable entropy function of the form 5.1 where σj(V) and ρj(T) are arbitrary functions of their arguments. The system (2.6) reduces to the following equation: 5.2 The characteristic curves P(T) are defined by the equation 5.3 As the solution V to equation (5.2) is constant along the characteristics, integration of (5.3) leads to the following implicit formula of the solution: 5.4 where f(V) is an arbitrary function of its argument. As discussed above, the function f(V) can be fixed by specifying either a particular isobar or a particular isotherm.

### (a) Examples

Few cases of physical interest that can be obtained as a virial expansion with coefficients depending on the temperature fit in the class of separable entropy functions discussed above. We consider here some explicit examples.

#### (i) Classical plasma

The case of a completely ionized gas (classical plasma) can be recovered setting σj=ρj=0 for j≥2, f(V)=0 and where c is a certain constant. Equation (5.4) gives the equation of state for a classical plasma (Landau et al. 1980) 5.5

#### (ii) Redlich–Kwong equation

Setting σj=ρj=0 for j≥2, f(V)=0 and where c and b are some constants, equation (5.4) gives the Redlich–Kwong equation of state (Redlich & Kwong 1949) 5.6 Equation (5.6) usually provides a more accurate description of a real gas phase transition than the van der Waals equation.

#### (iii) Peng–Robinson equation

Setting σj=ρj=0 for j≥2, f(V)=0 and and where a, c and Tc are constants, equation (5.4) gives the equation 5.7 which is known as the Peng–Robinson equation of state (Peng & Robinson 1976). This equation arises in the theory of non-polar liquids.

#### (iv) Dieterici equation

The choice with j=0,1,2,… and a a suitable constant, provides the Dieterici equation of state (Dieterici 1899) 5.8 A straightforward computation shows that

#### (v) Ideal Bose gas equation

Setting σj=ρj=0 for j≥1, f(V)=0 and where , α and c are constants depending on the physical properties of the system, Liα+1 is the polylogarithm function, ζ is the Riemann zeta function and Tc is the critical temperature at which a Bose–Einstein condensate begins to form. Equation (5.4) results into the ideal Bose gas (Landau et al. 1980) state equation 5.9

## 6. Concluding remarks

We have shown that a large class of equations of state for real gases can be obtained as solutions of a set of integrable hyperbolic PDEs. This set of equations is equivalent to the balance equation (1.1) under the assumption that the entropy function is separable in the variables V and T. The simplest case leads to a set of equations of hydrodynamic type for which the characteristic speeds are a function of the dependent variable V only. The specific functional form of the unknown characteristic speeds as well as the arbitrary functions resulting from the integration procedure can explicitly be determined using a suitable number of particular isothermal or isobaric curves.

For more general separable entropy functions, the hyperbolic conservation law can still be explicitly integrated. This case contains few examples of equations of state for real gases that can be viewed as a virial expansion with coefficients depending on the temperature.

The equations of state predict the occurrence of a phase transition that is understood as the shock wave dynamics of solutions to a hyperbolic PDE. The critical point where the phase transition starts to develop corresponds to the gradient catastrophe set where the shock starts to form. The singularity and catastrophe theory (Lu 1976; Poston & Stewart 1978; Arnold et al. 1993) turns out to be the natural language for the study of local properties of a system near the critical point. Further analysis of the deformed virial expansion involves the study of stable unfoldings of higher co-dimension.

The approach discussed above can be generalized for the study of general composite systems in thermal contact to be associated with a multi-component system of hydrodynamic PDEs.

A more accurate description of the thermodynamic functions near the phase transition should concern the study of more general entropy functions depending on the first and possibly higher derivatives of V . This will produce diffusive and/or dispersive higher order equations. In fact, it should be emphasized that discontinuous solutions, such as in figure 3, can be obtained as the small diffusion limit of the Burgers equation (Whitham 1974).1 This suggests that in proximity to the phase transition higher order effects should play a role in the description of the real phase transition.

## Acknowledgements

We are pleased to thank Boris Dubrovin, Boris Konopelchenko, Kenneth McLaughlin, Tamara Grava, Andrea Raimondo for useful discussions and comments. We are grateful to Giacomo Dossena whose useful comments helped us to improve the manuscript. We also thank the referees for their stimulating comments. A.M. is supported by the ERC grant from PDE (P.I. Boris Dubrovin) and the grant ‘Giovani Ricercatori della SISSA (P.I. A.M.). G.D.N. is supported by the grant no. ANR-08-BLAN-0261-01. A.M. and G.D.N. are partially supported by the Grant of Istituto Nazionale di Alta Matematica - GNFM, Progetti Giovani 2011 (P.I. A.M.).

## Footnotes

• 1 We thank Kenneth McLaughlin for pointing our attention to the results on the small diffusion limit of the Burgers equation.