Thermodynamic properties of polydisperse fluid mixtures

We present a systematic study of the thermodynamic properties of a poly-disperse fluid mixture. The size of the particles, σ , is assumed to be distributed according to a continuous distribution function f Σ (σ) , for which we have chosen a Γ-distribution. The interatomic potentials are given by a hard core repulsion plus an adjacent attractive tail in the form of a square-well or a Yukawa potential; for the size-dependence of the attraction strength we have assumed different models. The properties of the mixture are calculated using the optimized random phase approximation (ORPA), a ther-modynamic perturbation theory which is known to give reliable results in the case of simple liquids. To take into account polydispersity we combine the ORPA with the orthogonal decomposition technique where all σ-dependent functions (i.e., the correlation functions and the interatomic potentials) are expanded in terms of orthogonal polynomials p i (σ) associated with the weight function f Σ (σ) .


Introduction
For the past years considerable effort has been devoted to the concepts that make it possible to determine structural and thermodynamic properties of polydisperse fluid mixtures [1].Certainly enough, a very important reason for these activities lies in the fact that suspensions of mesoscopic particles (with diameters ranging from 1 µm to 1 nm), as they are often realized in colloidal suspensions, are in general -as a consequence of their production process -polydisperse.This means that the properties of the particles (size, interatomic potentials, etc.) are continuously distributed according to a distribution function.Such colloidal suspensions play a fundamental role in our every day life: ink, milk, as well as fog and smoke being only a few examples.Further, they are of great biological importance and have vast industrial applications [2].
The first concepts of describing the structural and thermodynamic properties of polydisperse systems date back to the late 1970s [3,4].In subsequent contributions, the concept of polydispersity was cast in a mathematically rigorous form by Salacuse and Stell [5,6] and by Briano and Glandt [7] in the 1980s.Within this framework such systems are viewed as mixtures with formally infinite number of components, each of them characterized by a continuous index X, which replaces the discrete index i in a mixture with a finite number of components.The most convenient choice for X in the systems of spherically symmetric particles is the size of the particles, σ.The amount of particles of size σ in the polydisperse mixture is then determined from a continuous, positive, and normalized size-distribution function f Σ (σ) which replaces the discrete set of concentrations in a mixture with a finite number of components.
The frameworks proposed in recent years in order to describe the properties of polydisperse mixtures using the tools of statistical mechanics were developed by merging this concept for polydispersity with standard liquid state theories that are known to be successful for simple liquids and their mixtures [8]: (i) On the one hand, benefit was taken from the availability of analytic solutions for a certain class of model systems, such as hard spheres (HS), adhesive hard spheres, hard sphere Yukawa, or charged hard sphere systems [4,[9][10][11][12][13][14].Generalized to polydisperse mixtures these models are called truncatable free energy models [1].Their thermodynamic properties can be expressed by a limited number of generalized moments of the distribution function which permits to determine the phase diagram [13,14].(ii) On the other hand, standard liquid state theories, such as integral-equations techniques or perturbation theories, have been successfully generalized to polydisperse mixtures, thus opening access to the systems with arbitrary interactions [15][16][17].Again, the structure and thermodynamic properties can be rather easily calculated.However, the solution of the coexistence equations to determine phase diagrams is -due to the numerical complexity of the problem -still out of reach [18,19].
In this contribution we present a systematic study of the thermodynamic properties of a polydisperse mixture of particles interacting via square-well (SW) and hard sphere Yukawa (HSY) potentials.In both cases, the repulsive part of the interactions are represented by simple hard spheres, a fact which substantially facilitates -without loss of generality -the calculations.The size of the particles is distributed according to a Γ-(or Schulz-) distribution [20].For the size-dependence of the interaction strength (i.e.potential amplitude) we assume different relations.In particular, we want to study the effect of polydispersity on the thermodynamic properties of the system: these are calculated via the optimized random phase approximation (ORPA) [21], a thermodynamic perturbation theory that has proven to give reliable results for one-and two-component systems [22,23].Generalization of this concept to the polydisperse case was recently presented in [17], where the ORPA was merged with a very efficient concept proposed by Lado [16], the so-called orthogonal decomposition technique.The basic idea of the latter concept is to decompose all σ-dependent functions, in particular, the correlation functions and the interatomic potentials, in terms of a set of orthogonal polynomials associated with the distribution functions f Σ (σ).
Since the determination of the phase diagram for the two model systems considered in this contribution is out of reach, we have taken reference -when "locating" states of systems -to the phase diagram of the corresponding one-component system.It can be considered as a particular polydisperse mixture with an infinitely sharp distribution function, i.e. f Σ (σ) = δ(σ − σ), σ being the mean value (or first moment) of σ.
We have calculated the isothermal compressibility, the virial pressure, and the free energy and have studied the effect of polydispersity on these properties for different densities and temperatures of the system.
The paper is organized as follows: in the subsequent section we summarize the theoretical basis of the study (i.e. the systems and the ORPA).In section III we present our investigations on the effect of polydispersity on thermodynamic properties of the systems.The paper is closed with concluding remarks.

Theory
In what follows we introduce the systems we have investigated and briefly summarize those expressions that are required in order to calculate the structure and the thermodynamic properties of a polydisperse mixture within the framework of the ORPA.For a more extensive presentation of the formalism we refer the reader to [17].

The systems
Our approach to the investigation of the properties of a polydisperse system is based on the concept introduced by Salacuse and Stell [5,6], where such a mixture is defined as a "system in which each particle is uniquely associated with a value of some characteristic parameter X, distributed according to a probability distribution function f X (x); X is a continuous random variable".Here, as in many other applications, we have chosen the random variable X to be the diameter of the particles, σ, which is distributed according to a normalized, positive size-distribution function f Σ (σ).In case of non-ambiguity, this function will be denoted in the following as f (σ).The fraction of particles with diameter σ [σ, σ + dσ] is, therefore, given by f (σ)dσ; the system is characterized by a number density ρ and a temperature T The (spherically symmetric) pair potentials acting between particles of size σ 1 and σ 2 and separated by a distance r are denoted by Φ(r; σ 1 , σ 2 ).Here we consider only hard-core potentials, i.e.Φ(r; σ 1 , σ 2 ) can be split into a repulsive HS reference part Φ 0 (r; σ i , σ j ) and a perturbation part Φ 1 (r; σ i , σ j ), In the following we will use indices "0" and "1" for the properties related to the reference and to the perturbation part.The HS potentials Φ 0 (r; σ 1 , σ 2 ) are characterized by additive diameters (Lorentz-rule), d(σ 1 , σ 2 ), i.e., with d(σ, σ) = σ.The structural and thermodynamic properties of the reference system are explicitly known within the Percus-Yevick (PY) approximation [4].
For the attractive tail in Φ(r; σ 1 , σ 2 ) we have considered two model interactions that are frequently used in liquid state physics: the SW and the HSY potential.
The SW potential is given by here, the (σ 1 , σ 2 ) are the parameters for the well-depth and the λ(σ 1 , σ 2 ) characterize the range of the well.This parametrization permits us to take into accountindependently from each other -polydispersity in particle size, interaction strength, and potential range with distribution functions f (σ), f E (σ), and f Λ (σ).
The HSY potential is given by the K(σ 1 , σ 2 ) being the contact values and the z(σ 1 , σ 2 ) being the screening lengths, distributed according to f K (σ) and f z (σ).
For simplicity we assume the particle diameter σ to be distributed according to a Γ-(or Schulz-)distribution [20], where σ is the mean diameter and α is a positive parameter; it is related to the polydispersity index, I, via I = (α+2)/(α+1).The Γ-distribution, being frequently used to describe polydisperse fluids, has the attractive feature that the orthogonal polynomials p i (σ) associated with the Γ-distribution are known explicitly (see below): they are proportional to the associated Laguerre polynomials Further, for the moments of this distribution, m i = ∞ 0 dσf (σ)σ i , closed expressions can be presented (see, for instance, [18]); in particular, the third moment, m 3 , is related to the packing fraction, η, via In figure 1 we have depicted the Γ-distribution function for four different α-values.
The OZ equation is supplemented by the random phase (RPA) closure relation or equivalently, and the core condition (optimization criterion), Using the orthogonal expansion technique proposed by Lado [16], we introduce orthogonal polynomials p i (σ) that are associated with the distribution function f (σ), i.e.
δ ij is the Kronecker-delta.Starting from p 0 (σ) = 1, the polynomials can be constructed -in case they are not known explicitly -via a Gram-Schmidt orthogonalization procedure.All size-dependent functions, x(r; σ) and y(r; σ 1 , σ 2 ), are expanded in terms of these polynomials, i.e., with the coefficient functions Truncating the expansions in ( 13) at some suitably chosen value n, the coefficient functions x i (r) and y ij (r) are most conveniently calculated in practical applications via a Gaussian quadrature algorithm, using the set of n roots, {σ k;n , k = 1, . . ., n}, of the polynomial p n (σ).Detailed numerical studies have shown (for details we refer the reader to [17,18]), that n = 5 is a sufficiently large value to truncate the expansions and simultaneously maintaining a high accuracy for the results.
To be more specific, we can then represent the direct correlation functions of the polydisperse system, c(r; σ 1 , σ 2 ), via similar relations hold for the g(r; σ 1 , σ 2 ) and the h(r; σ 1 , σ 2 ).
It is now convenient to collect the coefficient functions of the correlation functions to matrices, e.g., C ij = c ij (r), and similar relations for the h ij (r) and g ij (r).Further we use the tilde for functions in q-space.We also introduce the coefficient functions for the static structure factor of the HS reference system, S 0;ij (q); they are given by S 0;ij (q) = 1 + ρ h0;ij (q) and are collected to a matrix S 0 (q).
As also outlined in detail in [17], one can show that the solution of the OZ equation along with the ORPA closure relation ( 9), (10), and ( 11) is equivalent to minimizing the functional 3 dq Tr ρ C1 (q)S 0 (q) + ln det 1 − ρ C1 (q)S 0 (q) (17) with respect to variations of the C 1;ij (r) inside the core.In the above equation 1 is the unit matrix and "Tr" denotes the trace of a matrix.For technical and numerical details of the solution of the ORPA we refer the reader to [17].

Thermodynamic properties
Once the ORPA has been solved, several thermodynamic properties can be easily calculated from closed expressions; for the virial pressure P v and the isothermal compressibility χ T they are given by where P c is the compressibility pressure.Within the ORPA and using the coupling constant formalism [16], one can derive for the Helmholtz free energy A the expression where A 0 is the free energy of the polydisperse HS reference system [4] and the other two contributions are found to be N is the total number of particles of all species and V (required below) is the volume of the system.Following the route of Høye and Stell [25] in order to calculate the chemical potential in a one-component system for the ORPA (or the mean spherical approximation), we could derive an explicit expression for the chemical potential of particles of size σ, µ(σ), as follows: we have first generalized their expressions for the µ i of a mixture with a finite number of components M and a set of concentrations {c i }, and have then made the transition to a polydisperse mixture, replacing the c i by f (σ).We finally obtain -for details see [18] -the following expression: (23) or, in terms of the coefficient functions, Index "id" refers to the ideal gas contribution, index "ex;0" to the excess (over ideal gas) contribution of the HS reference system.

Polydisperse square-well mixture
Using the parametrization (3) for the SW interaction we have used the following relation between the well-depth of (σ 1 , σ 2 ) and the particle size: For the well-range we put λ(σ 1 , σ 2 ) ≡ 1.5, i.e. a size independent constant.Relation (25) guarantees that the Berthelot rule [8,12] is satisfied for the interatomic strength.Note that the above equation induces a distribution function for the well-depth f E (σ).For the parameter z in (25) we have chosen z = 0, 1, and 2: for z = 0, i.e., all potentials have -irrespective of the particle diameter -the same attraction, the system is polydisperse only in size, while for z = 1 and 2 our model shows both sizeand amplitude-polydispersity.In figure 2 we show the phase diagram (in terms of the binodal and the spinodal) of the corresponding one-component system [i.e.f (σ) = δ(σ − σ)]; this case has been realized by α = 1000 which leads to a very sharp Γ-distribution (c.f.figure 1).Density is measured in units of ρ * = ρσ 3 .The two dot-dashed, horizontal lines indicate those reduced temperatures t (t = k B T/¯ ) along which the trends in the thermodynamic properties have been investigated.In the following figures of the thermodynamic properties the two vertical dotted lines locate ρ * -values on the spinodal of the onecomponent system for the corresponding isotherms (t = 1.25 and t = 1.12).
As we first look at the characteristic trends in these figures we observe, that, in particular, if the amplitude-polydispersity is strong (z = 2) and the α-values are rather small, the ORPA breaks down close to the spinodal of the one-component case.This can be understood as follows: firstly, as α becomes smaller, the tails of the distribution functions extend to larger σ-values and hence the probability to find large (and therefore strongly attractive) particles in the mixture increases; secondly, the high z-value introduces an additional strong attraction for the large particles even though they might only represent a small fraction of the particles in the system.The combination of these two trends leads to a breakdown of the ORPA since the basic idea of a perturbation of reasonable size due to the square-well is no longer justified.We point out that a similar scenario was encountered for the polydisperse HSY mixture within the MSA [13].
In figures 3 and 4 we show the reduced isothermal compressibility χ * T = ρk B T χ T for the two t-values considered, for z = 0, 1, and 2, and for different degrees of polydispersity; note that α = 1000 represents the one-component case.First we observe the general tendency that the curves are shifted to smaller ρ-values as the polydispersity is increased.Furthermore, we see that size-polydispersity leads to a small shift to lower values in the χ T -curves as the distribution becomes broader.If we take into account both size-and amplitude-polydispersity, the compressibility curves show larger differences with respect to the monodisperse case.The trends outlined above hold for both temperatures considered.
In figures 5 and 6 we show the trends in the virial pressure as density and polydispersity varies.Again we find qualitatively similar trends for the different temperatures.Size-polydispersity alone leaves the pressure in the gas phase unaffected; decreasing α, monotonically increases the pressure, in particular for high densities.As we introduce amplitude-polydispersity, we observe an effect in the pressure of the gas phase while in the fluid phase we also note a decrease in the pressure values for the intermediate densities.In addition, for the case z = 2 the pressure can also become negative.
Further, we have also calculated the pressure directly by differentiating the free energy A with respect to the density, P a .Comparison with the virial pressure gives some insight into the thermodynamic consistency of the ORPA.As space is limited, we have not displayed the curves here and rather discuss the results: the trends observed for P v and P a are qualitatively similar.For small and intermediate densities the internal consistency is very satisfactory, while for larger densities differences of the order of several percent are observed.We point out that this situation is similar to the one observed in simple one-component fluids.
Finally in figure 7 we display the results for the free energy for t = 1.25; for t = 1.12 we observe, again, a qualitatively similar behaviour.Polydispersity in size alone leads to an increase in the free energy for rather high densities only, leaving the small-        and intermediate-density range unaffected.As polydispersity in amplitude comes into play, we observe a decrease as α becomes smaller which leads -in particular for intermediate and large densities -to very pronounced effects for z = 2.

Polydisperse hard sphere Yukawa mixture
We now proceed to the results of the polydisperse HSY mixture.For the sizedependence of the parameters of the interatomic potential (4) we have chosen the following relations which were motivated by the van der Waals study in [24] In contrast to the SW case, disentanglement of size-and amplitude-polydispersity is with this parametrization no longer possible due to the exponent in (4).z has been put to 1.8, the reduced temperature t is given by t = k B T/ K, again ρ * = ρσ 3 .In figure 8 we show the phase diagram (in terms of the binodal and the spinodal) of the corresponding one-component system [i.e.f (σ) = δ(σ − σ)].The two dot-dashed, horizontal lines indicate those reduced temperatures t (t = k B T/¯ ) along which the trends in the thermodynamic properties have been investigated.In the following figures of the thermodynamic properties the two vertical dotted lines locate ρ * -values on the spinodal of the one-component system for the corresponding isotherms (t = 2.16 and t = 1.98).
In figures 9 and 10 we show the results for the isothermal compressibility.With respect to the SW system, we see for a = 0 (amplitude-polydispersity via the exponent) an enormous effect on the compressibility: going from the one-component case  to α = 20 we observe a decrease by a factor of three (in particular for t = 1.98).This strong effect is somewhat reduced as we introduce an additional amplitudepolydispersity via the contact value.We also observe that including polydispersity in the HSY case always makes the compressibility smaller, while we have found opposite tendencies in the SW system.
For the virial pressure P v , shown in figures 11 and 12, we observe, that already a = 0 (i.e.amplitude-polydispersity via the exponent) leads to a strong effect on the pressure even for intermediate densities (while in the SW case only the high-density data were affected).Again, including amplitude-polydispersity via the contact value, brings the curves for different α-values rather close together, while in the SW system we observe the opposite tendency.The region where the pressure becomes negative is now larger than in the SW system.For the internal consistency between P v and P a we have found that it is now slightly worse than in the SW case.
Finally, we investigate the free energy curves in figure 13 for the reduced temperature t = 2.16.Again, we observe that an additional amplitude-polydispersity  brings the curves for different α's rather together with respect to the case that only amplitude-polydispersity via the exponent is present.

Conclusions
We have presented a systematic study of the effect of size-and amplitudepolydispersity on the thermodynamic properties (isothermal compressibility, pressure, and free energy) in a polydisperse fluid mixture.Assuming a Γ-distribution for the size of the particles and different relations between the size and the interaction strength of the particles of the system we have studied mixtures where the particles interact via a SW or a HSY potential.The liquid state theory we have used is a merger of the ORPA and the orthogonal decomposition technique, where the size dependent functions (i.e.correlation functions and interatomic potentials) are decomposed into the orthogonal polynomials associated with the distribution function.Taking the one-component system as a particular polydisperse reference system, we could point out the effect of size-and amplitude-polydispersity on the thermodynamic properties as functions of the density; in these studies we have chosen two particular temperatures.We have observed different tendencies in the SW and in the HSY systems which are obviously due to the longer range of the latter one.Similar to the MSA-based truncatable free energy models we have found that a distribution function, such as the Γ-distribution with an unlimited carrier in σ-space must lead to a breakdown of the ORPA when the distribution has exceeded some value of broadness: then, even though with a very small probability, large particles with a very strong attraction are encountered in the mixture and the premise of the weak perturbative character underlying the concept of the ORPA is violated.

Figure 2 .
Figure 2.Phase diagram of the monodisperse SW system, i.e. f (σ) = δ(σ − σ), calculated within the ORPA: dotted line -spinodal, full line -binodal.The two dot-dashed lines for t = 1.12 and t = 1.25 denote the isotherms along which the trends in the thermodynamic properties in the following figures have been investigated.

Figure 3 .
Figure 3. Reduced dimensionless isothermal compressibility χ * T = ρk B T χ T for the polydisperse SW mixture (as defined in the text) for a reduced temperature t = 1.25; α-values as indicated, z = 0, 1, and 2 (from top to bottom).

Figure 4 .
Figure 4. Reduced dimensionless isothermal compressibility χ *T = ρk B T χ T for the polydisperse SW mixture (as defined in the text) for a reduced temperature t = 1.12; α-values as indicated, z = 0, 1, and 2 (from top to bottom).

Figure 5 .
Figure 5. Reduced virial pressure P v, * = βP v for the polydisperse SW mixture (as defined in the text) for a reduced temperature t = 1.25; α-values as indicated, z = 0, 1, and 2 (from top to bottom).

Figure 6 .
Figure 6.Reduced virial pressure P v, * = βP v for the polydisperse SW mixture (as defined in the text) for a reduced temperature t = 1.12; α-values as indicated, z = 0, 1, and 2 (from top to bottom).

Figure 7 .
Figure 7. Reduced free energy A + = βA/V for the polydisperse SW mixture (as defined in the text) for a reduced temperature t = 1.25; α-values as indicated, z = 0, 1 and 2 (from top to bottom).

Figure 8 .
Figure 8. Phase diagram of the monodisperse HSY system, i.e. f (σ) = δ(σ − σ), calculated within the ORPA: dotted line -spinodal, full line -binodal.The two dot-dashed lines for t = 1.98 and t = 2.16 denote the isotherms along which the trends in the thermodynamic properties in the following figures have been investigated.

Figure 9 .
Figure 9. Reduced dimensionless isothermal compressibility χ * T = ρkBT χ T for the polydisperse HSY mixture (as defined in the text) for a reduced temperature t = 2.16; α-values as indicated, a = 0, and 1 (left and right).

Figure 10 .
Figure 10.Reduced dimensionless isothermal compressibility χ * T = ρkBT χ T for the polydisperse HSY mixture (as defined in the text) for a reduced temperature t = 2.16; α-values as indicated, a = 0, and 1 (left and right).

Figure 11 .
Figure 11.Reduced virial pressure P v, * = βP v for the polydisperse HSY mixture (as defined in the text) for a reduced temperature t = 2.16; α-values as indicated, a = 0 and 1 (from top to bottom).

Figure 12 .
Figure 12.Reduced virial pressure P v, * = βP v for the polydisperse HSY mixture (as defined in the text) for a reduced temperature t = 1.98; α-values as indicated, a = 0 and 1 (from top to bottom).

Figure 13 .
Figure 13.Reduced free energy A + = βA/V for the polydisperse HSY mixture (as defined in the text) for a reduced temperature t = 2.16; α-values as indicated, a = 0 and 1 (left and right).