Phase Transition in the Boltzmann–Vlasov Equation

In this paper we revisit the problem of explaining phase transition by a study of a form of the Boltzmann equation, where inter-molecular attraction is included by means of a Vlasov term in the evolution equation for the one particle distribution function. We are able to show that for typical gas densities, a uniform state is unstable if the inter-molecular attraction is large enough. Our analysis relies strongly on the assumption, essential to the derivation of the Boltzmann equation, that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu \ll 1,$$\end{document}ν≪1, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu =d/l$$\end{document}ν=d/l is the ratio of the molecular diameter to the mean inter-particle distance; in this case, for fluctuations on the scale of the molecular spacing, the collision term is small, and an explicit approximate solution is possible. We give reasons why we think the resulting approximation is valid, and in conclusion offer some possibilities for extension of the results to finite amplitude.


Introduction
The application of the Boltzmann equation to the fluid mechanical properties of matter has been extraordinarily successful. The Chapman-Enskog expansion [8] of its solution provides a mechanistic explanation of viscosity and thermal conductivity, and equally Boltzmann's Htheorem provides an explicit basis for the concept of entropy and the laws of thermodynamics.
In this context, one of the intriguing properties of materials is that of phase transition, which on the face of it is inconsistent with the H -theorem, which provides a Lyapunov function for the velocity distribution which drives it to a unique (Maxwellian) equilibrium. On the other hand, phase changes correspond to the existence of multiple steady states, as instanced for example by the van der Waals equation, which seems prima facie contradictory to the H -theorem. Actually, it is not as straightforward as this, since the equilibrium Maxwellian distribution depends on number density and temperature, and assumes these quantities are spatially uniform; if one allows spatially varying number density, other possibilities can occur, and indeed this is the subject of the present paper.
Of course, the resolution of this apparent conundrum lies in the fact that the H -theorem assumes a one particle distribution function f (r, v, t) which is independent of the spatial variable r. If this restriction is lifted, the possibility of spatial variation arising through an instability occurs, with the interpretation of such an instability being the onset of condensation. It is therefore of interest to examine the stability of the spatially uniform equilibrium of the Boltzmann equation, with a view to understanding how phase transition occurs.
There is, however, nothing in the Boltzmann equation which seems to provide a mechanism for instability; to provide such a mechanism, we will interpret the intermolecular potential as consisting of two parts: a steep repulsive part, which we model as the usual Boltzmann collisional interaction, and a longer range attractive part, and it is this which provides the instability mechanism: attractive forces between molecules induce a tendency to form clusters, and this tendency is enhanced in conditions of high density or low temperature.
A number of authors have investigated this problem. The original approach uses methods of classical statistical mechanics (e. g., Born and Fuchs [4], Mermin [16] and van Kampen [20]), which allow the derivation of the van der Waals equation of state from a virial expansion of the grand partition function (e. g., Schwabl [19], p. 236 ff). Later authors address the problem directly. For example, Grmela [13] enunciates essentially the same problem which is of concern here. And indeed, his approach is similar: he writes down a Boltzmann-Vlasov equation similar to that which we present below, where the reference to Vlasov reflects the additional attractive term which is used in theories of plasma dynamics [22]. But in common with much of the literature on kinetic theory, the presentation is discursive, and the difficulty of dealing with the linearised collision operator prevents any clear conclusion.
Liboff [15] also presents similar stability results to those we derive below, for a range of different attractive power law potentials, but largely ignores the effects of collisions, only modifying his analysis by including a simplified (algebraic) collision term due to Bhatnagar et al. [3], yielding the so-called BGK approximation, which is analytically tractable. A similar approach has recently been adopted by Benilov and Benilov [2].
De Sobrino [10] takes the application of the Boltzmann equation further. Insofar as phase transition can be understood by derivation of the van der Waals equation of state, he shows that by including the two essential constituents of this equation, intermolecular attraction and molecular crowding, in the prescription of the closure for the two-particle distribution function, the van der Waals equation emerges as the equation of state of the equilibrium solutions. Other approaches and related discussions can be found in the papers of Penrose [18], Wisnivesky [21] and Chen [9].
In this paper we consider the problem of condensation from the point of view of kinetic theory, by analysing the Boltzmann-Vlasov equation for the one particle distribution function, which includes the Boltzmann collision integral and the Vlasov term describing the intermolecular attraction. Although the ideas of scale have been adumbrated before, here we explicitly non-dimensionalise the model, and show that the collision term is small, providing the gas is sufficiently rarified, in the sense that ν = d/l 1: here d is molecular diameter and l is the mean inter-particle spacing; in practice this assumption is justified, except near the critical point. This allows us to derive an explicit condition for instability, much in the manner of Liboff [15]; however, we go further by considering explicitly the corrective rôle of the collision integral, and we show that it is small, although we surmise that this is not the case in the condensed phase equilibrium.
The rest of the paper proceeds as follows. In Sect. 2, we present the Boltzmann equation and its modification by the Vlasov term, and we non-dimensionalise the system, which introduces two important dimensionless parameters: ν, as described above, and β, which measures the strength of the inter-molecular attraction. In Sect. 3, we linearise the model to study spatial instability of the Maxwellian equilibrium, and derive an explicit instability criterion when the density parameter ν is small. We then study the perturbative effect of the collision integral, and show that its effect is uniformly small, thus validating the accuracy of the approximate stability criterion. Finally, in the concluding Sect. 4, we offer some conjectures about the subsequent nonlinear evolution of the system.

The Boltzmann Equation
The basic equation of statistical mechanics is the Liouville equation, from which we can derive the BBGKY hierarchy: Here f s is the s-particle distribution function, and is a function of the positions r i ∈ V and velocities v i ∈ U of the s particles, as well as time t, P = V × U is the six-dimensional configuration space of position and velocity, and dγ i is the volume element of that space, g is the external acceleration (typically gravity), and a i j is the inter-particle acceleration on particle i due to particle j. If the inter-particle potential is W i j = W (|r i − r j |), then A typical example of such a potential is the Lennard-Jones potential given by where d is the molecular diameter; the potential is portrayed in Fig. 1. Particularly, the one-particle distribution function f 1 = f (r, v, t) satisfies the first equation in the hierarchy: where a = a 12 . The Boltzmann equation now follows from the assumption that the locations of two particles in P are independent of each other, so that we take (2.5) and in this case we obtain where the collision integral takes the form

Intermolecular Forces
Henceforth we ignore the external body force, thus g = 0. We are now motivated by the form of the potential in Fig. 1 to consider W to consist of two parts: an increasing part for r > d corresponding to long range attractive forces, and a vertical line at r = d which corresponds to collisions between particles. In more detail, we consider W = W C (r ) + W LR (r ), with both potentials being zero for r < d, and W LR < 0 is an increasing function, tending to 0 as r → ∞, while the collisional potential W C is 0 for r > d, but ranges from 0 to ∞ at r = d, representing perfectly elastic collision. In this case, the collision integral can be taken to be the sum of two terms, where Q C is the usual Boltzmann form of the collision integral representing perfectly elastic collisions, here U = R 3 is velocity space, v and w are the velocities of two colliding particles, V their relative velocity,k a unit vector along the line between their centres at the point of collision, and d = d 2 V ·k dω, with the solid angle element dω being taken over all solid angle subtended at v such that V ·k > 0 (so that the particles are actually colliding, not separating). A more accurate expression due to Enskog allows for the finite size of the particles, but this is not necessary in the present discussion. Further discussion of the Enskog modification is provided in Sect. 4. Q C is associated with the vertical (repulsive) part of the potential, while Q LR represents the longer range attractive force, where for an inter-particle potential as in (2.2), we have from (2.7) is the acceleration associated with the attractive part of the potential [thus we take a = 0 for ξ < d: equivalently, the integral in (2.11) is taken over ξ ∈ R 3 , ξ > d], and V is spatial volume. As mentioned above, the integral in (2.9) over + is with respect to solid angle in velocity space. Specifically, as shown in Fig. 2, integration is over the solid angle subtended at v over the antipodal sphere through the antipodal points v and w, and more precisely where d is molecular diameter, dω is the element of solid angle indicated in Fig. 2 Note that d has units of velocity. For the Lennard-Jones potential in (2.3), we would take (2.15) and this will be assumed later [after (3.14)]; for now a(ξ ) is quite general. The Boltzmann equation (2.6) with g = 0 thus takes the form and we write ∇ r = ∇. This equation forms the basis of our study.

Non-dimensionalisation
It is convenient to non-dimensionalise the model. To do this, we define the thermal velocity scale an acceleration scale and a mean inter-molecular distance where n 0 is a typical value of the number density. We then scale the variables 1 as 20) and this leads us to the non-dimensional form of the equation, in which all the variables are dimensionless (in particular d = d S/V ), we have dropped the asterisk from a * , and the dimensionless parameters are defined by This form of scaling the equation is that used by Keller [14], although in his case the interpretation of the parameter ν is rather different.
In the absence of any spatial variation, A = 0 and the left hand side of the equation is simply ∂ f /∂t. As is well known, the existence of Boltzmann's H -function assures us that f approaches an equilibrium on a time scale of O 1 ν 2 , and this is the Maxwellian distribution, which in its dimensionless form is given by where for convenience we will assume that the mean velocity is zero. It is the stability of this state to spatial perturbations which is our concern.

Linear Stability
Before we linearise the equation, it is convenient first to define so that (2.21) 1 takes the form where for ψ(v), The Maxwellian corresponds to and 0 = 0. We now linearise the equation about the steady state by writing and neglecting nonlinear terms; the linearised form of (3.2) is where the linearised form of A is and the linearised collision operator is We now seek normal mode solutions to this equation of the form where k is the wave vector. This leads to the eigenvalue problem for ψ(v) in the form where and we have assumed a normalisation in which To evaluate B, we take Cartesian coordinates in V = R 3 with the z axis in the k direction. By symmetry, the x and y components are zero, and therefore we can write and by changing to spherical polar coordinates, we find that (writing K = νk) This expression was derived by Liboff [15]. The function C(K ) is plotted in Fig. 3; its asymptotic limits for small and large K are [taking a = 6/r 7 , corresponding to (2.15), and with a 0 given in (2.18)] The function C(νk) = C(K ) defined by (3.14). The thin lines give the asymptotic limits from (3.15) Some further detail of the expansion for small K is given in the Appendix. As a consequence, (3.10) is Our aim is to solve this equation to determine σ (k). Eigenfunctions ψ for which Re σ > 0 are unstable.

The Limit → 0
Before we begin, we summarise some of what is known about the linear integral operator L (e. g., Cercignani [7], p. 159 ff). Under the definition of an inner product

17)
L is self-adjoint. Its null-space N is the span of the functions 1, v and 1 2 v 2 ; all other eigenfunctions have negative eigenvalues. Something is known about the action of L on tensors; for example, and such results can be extended to higher order tensors; they are due to a rotational invariance of the operator L. However, the ability to directly solve an equation such as (3.16), despite all this structure, is not clear, and we therefore resort to an approximate method.
, we would hope that as the ambient temperature decreases, the instability that we seek will appear at a sufficiently large value of β. On the other hand, while molecular spacings in liquids are of order d, those in gases are significantly larger (because gas densities are typically much lower at normal temperatures and pressures). 2 This suggests that we consider a gas in which l d, and thus ν 1. This immediately allows us to provide an approximate solution of (3.16). The neglect of the collision integral requires consideration, in case it allows a singular perturbation; in the sequel we give some consideration to this possibility, although we might suspect, since the neglected term is an integral, that the induced perturbation is in fact regular.
If we neglect the term of O(ν 2 ) in (3.16), we simply have (3.19) and the normalisation condition (3.12) implies and it is this which determines σ.
To evaluate (3.20), we take Cartesian axes (v x , v y , ζ ) in the velocity space U , with the ζ axis in the direction of k, so that k · v = kζ. Carrying out the integrals in v x and v y , this leads to Defining σ = √ 2kη, (3.22) this can be manipulated to the form and then This integral is related to a form of the 'plasma dispersion function' (Fettis et al. [12], Abramowitz and Stegun [1], p. 297). Specifically, if we define w(iη) ≡ W (η) = e η 2 erfc η, (3.25) then, providing Re η > 0, the integral in (3.24) is proportional to W (η), and If η is a root of (3.23), then so is −η (by taking the complex conjugate of the equation): it follows that we can take Re η ≥ 0 without loss of generality, since for Re η < 0, we simply consider −η. By taking the imaginary part of the second integral in (3.24), it follows that in fact either η is real or purely imaginary. The latter case corresponds to neutral stability, and is discussed later. For real η, we can then suppose η > 0. In this case, (3.24) shows that a root only exists if βC > 1, and this is shown in Fig. 4. This is the principal result that we obtain: for values βC > 1, the uniform state is unstable, and we associate this with a transition to condensation. With C(νk) as shown in Fig. 3, we see that the criterion for instability is that , (3.27) or in dimensional terms, , (3.28) Fig. 4 The scaled growth rate η = σ √ 2k as a function of βC, given by (3.26) where for illustrative purposes we have used the perfect gas law p = n 0 kT 0 , even though this may not be appropriate for a Boltzmann-Vlasov gas. We can see from Fig. 3 that in the typical case where β ∼ O(1) (and β > 3 8π ), βC decreases to one as k increases from zero to a value of O(1/ν), and thus Fig. 4 shows that η decreases monotonically to zero as k increases. Hence σ is a concave unimodal function which first increases from zero and then decreases back to zero. In particular, the maximum growth rate occurs for k ∼ 1 ν , and is of order 1 ν .
Instability first sets in when β increases through 3 8π , and near this value, long wavelength disturbances are amplified since K 1. We use the expansion in (3.15) for small K , and expand (3.26) for small η, and this leads to the approximation for which the maximum growth rate occurs at k ≈ 1 3ν 10 8 In practice, it is only such long wavelength modes which are of interest.

An Approximate Correction
We now turn to the consideration of the accuracy of this stability criterion. The issue is whether the neglect of the collision integral term ν 2 Lψ in (3.16) is justified. Taking (3.19) as our definition of ψ, we have , (3.31) and it follows that It can be seen that the concession admitted in the last inequality is quite weakening, but we have not found any easy way round this. In a similar way we have (3.34) and thus that (3.35) The integral with respect to can be evaluated (using the geometry of Fig. 2), and is 2 3 where the integral proportional to v·w vanishes due to symmetry considerations. The integrals can be evaluated, and this leads finally to We surmise, however, that in fact this is not the case, and that actually |L(ψ)| ≤ O(v).
Our reasons for this surmise follow from consideration of Lψ when v is large. Because the Maxwellian f 0 decays rapidly as w increases, the integral over U in the definition of L is only significant when w ∼ 1. The integral over + thus splits into two parts: a region near the origin where d S ∼ 1 and thus d ∼ 1/V , where either v or w ∼ 1, ψ ∼ 1, and this part of the integral is O(1/V ). Over the rest of and thus we have (since and is thus of O(v): specifically, where we take the ζ axis in the direction of k, and η is still defined through (3.22). In view of (3.26), this is simply (3.42) and allows a correction to (3.19) as (3.43) and this could be used to provide an improved estimate for σ ; the main point, however, is that the collision term appears to provide a regular correction to our leading order result.

Wave Motion
The question arises, of course, as to what happens if βC < 1. In this case, the only possibility is a pure wave motion in which σ = −iω, and in that case the defining equation for ψ, (3.16), can be written as where again we neglect the collision term on the basis that ν is small , and we again normalise ψ by (3.12), that is, The solution of (3.44) is singular, (3.46) and it is necessary to enquire how one should interpret the resulting integral in (3.45). One way to do this is to consider (3.46) as the limit as ε → 0 of the non-singular solution , and carrying out the integral with respect to v x and v y (and with ζ in the direction of k), (3.45) leads to where c = ω/k is the wave speed. If we take ε > 0, then we can replace the contour in (3.48) by one which is indented by a lower semi-circle of radius δ 1. We then let ε → 0, and subsequently δ → 0, and this leads to taking the limit of (3.48) as which is equivalent to interpreting (3.46) as (3.50) where P indicates that principal value integration should be applied in (3.45), and we have used the result |a|δ(at) = δ(t). However, (3.49) evidently has no solution. What has gone wrong? The problem is that the singularity of (3.44) allows the more general solution to be where A is an arbitrary constant, since in the sense of generalised functions, tδ(t) = 0 (Carrier et al. [5], p. 320). This situation occurs in a number of other examples. Particularly, we may think of the (two-dimensional) stability of inviscid Couette flow [6], which can be described by the eigenvalue problem with boundary conditions φ(0) = φ(1) = 0. In the linearised problem, the basic flow is (u, v) = (y, 0), and the perturbed stream function is given by the normal mode solution There are no regular solutions of (3.52) at all (the discrete spectrum is empty), and solutions are obtained from where c = ω/k is the (real) wave speed. With 0 < c < 1, the solutions are just the Green's functions for the differential operator, and they form the continuous spectrum of the problem. 3 For the instability of more general shear flows (U (y), 0), the inviscid stability equation is the Rayleigh equation [11] (U − c)(φ − k 2 φ) − U φ = 0, (3.55) and the solutions are still singular The normalisation condition now leads, after some algebra, to where we have defined It may be noted that if βC < 1, the possibility of A = 0 is not available. In view of (3.57), the wave speed is √ 2 .

Discussion and Conclusions
It has been previously shown [15] that the Maxwellian equilibrium distribution of a (rarified) gas is spatially unstable at sufficiently low temperatures, providing the collisional operator can be ignored. Efforts to include collisional effects have been limited to the algebraic BGK approximation [2,15], or other such simplifications [13]. In this paper, we have followed a different approach analogous to that of Keller [14], where direct asymptotic solutions of the Boltzmann-Vlasov equation are sought. At leading order this reduces the stability problem to that of the collisionless Vlasov problem, but we have additionally shown that the correction due to the collisional term remains small, so that the resulting stability criterion remains accurate. It remains to consider the nonlinear evolution of the distribution function. This is beyond the scope of the present paper, but the results of de Sobrino [10] and the scaled equation in (2.21) suggest a way forward. First, the closure for the two-particle distribution function in (2.5) is replaced by f 2 (r, v; s, w; t) = C 1 2 (r + s), t f (r, v, t) f (s, w, t), (4.1) where C is a crowding coefficient, designed to represent the reduction of available space in dense conditions; de Sobrino suggests Second, the collision integral in (2.9) is replaced by the more accurate The dimensionless form of this, using (2.20), is where from (4.2), and additionally For small ν, the system reduces to that studied earlier, so that the instability result is the same. It is of at least mathematical interest to ask what state the system then evolves to. Naturally, we would hope that this state would correspond to a condensed liquid, although we recognise that it is generally considered that the Boltzmann (or Boltzmann-Vlasov, or Enskog-Vlasov) equation becomes inapplicable in this case, when the particle spacing is of order d. Despite this, it is of interest to enquire whether the solution of (4.7) evolves to a 'liquid-like' state. The obvious suggestion is that f evolves without the collision term towards a singularity, in which a rescaling becomes necessary when f ∼ 1/ν 3 , r ∼ ν, and thus n ∼ 1/ν 3 , A ∼ 1/ν, t ∼ ν and Q ∼ 1/ν 6 , following which the system is described by (4.5) and (4.7), with ν replaced by one everywhere. Exploration of this possibility is postponed to a future study. but the higher order terms evidently do not become effective until K 6, whereas the first order expression in (3.15) is already accurate for K 6.