Kirkwood-Buff or Fluctuation Solution Theory can be used to provide experimental pair fluctuations, and/or integrals over the pair distribution functions, from experimental thermodynamic data on liquid mixtures. Here, this type of approach is used to provide triplet and quadruplet fluctuations, and the corresponding integrals over the triplet and quadruplet distribution functions, in a purely thermodynamic manner that avoids the use of structure factors. The approach is then applied to binary mixtures of water + methanol and benzene + methanol over the full composition range under ambient conditions. The observed correlations between the different species vary significantly with composition. The magnitude of the fluctuations and integrals appears to increase as the number of the most polar molecule involved in the fluctuation or integral also increases. A simple physical picture of the fluctuations is provided to help rationalize some of these variations.

## I. INTRODUCTION

From a theoretical point of view, liquids and liquid mixtures can be characterized using relative probability distribution functions. Experimentally, the pair distribution function can sometimes be obtained via scattering studies. However, very little experimental information has been made available concerning triplet and higher distribution functions for real liquids and liquid mixtures despite the great interest in these quantities; for example, to understand the deficits of pairwise additive potentials.^{1–9} For pure monatomic liquids, the works of Schofield, Buff and Brout, Gray, Gubbins, and Egelstaff, among others, have shown that some information concerning the triplet distribution functions is available—in terms of integrals over the triplet and combinations of the pair distribution functions—for cases in which the derivatives of the radial distribution functions (RDFs) or structure factors, *S*(*Q*), with respect to pressure or density can be obtained numerically.^{7,10–19} These studies have been crucial for determining the range of applicability of the pair potential and for testing models of the three-body distribution function for the few (and relatively simple) systems that are amenable to this type of study.^{3,5,7,10,17–20}

Unfortunately, the situation is significantly worse for mixtures, especially those composed of polyatomic molecules, and hence even less is known experimentally about the triplet distribution function in these systems, primarily because it has only been possible to obtain RDFs for a very small subset of all possible mixtures.^{6,8,21,22} Nevertheless, a few mixture-based studies have attempted to decipher the importance of three body correlations. For example, Winter and coworkers recently studied derivatives of the structure factor in concentrated lysozyme solutions; however, this involved significant approximations.^{5} Furthermore, osmotic system studies using McMillan-Mayer (MM) theory provide the osmotic virial coefficients,^{23} which do provide information concerning pair and higher solute distributions, but with the restriction of non-volatile solutes at low concentrations.^{23}

We recently illustrated (for pure liquids only) how to bypass the dependence on scattering data to obtain information about triplet and quadruplet correlations from bulk thermodynamic data alone using Fluctuation Solution Theory (FST).^{24} Egelstaff and others previously obtained similar information for some simple fluids (e.g., pure Ne, Ar, Kr, and Rb) at a limited number of state points using scattering studies in which the pressure or density dependence of the long-wavelength limit of the structure factor, *S*(*Q* → 0), was sometimes considered in addition to the *S*(*Q*) derivatives.^{10,17,18,25,26} However, unlike the present study, their main focus was to use scattering data to obtain the higher order correlation data. Therefore, the systems and state points studied were very limited.

Here, we extend the approach presented in our previous study of pure liquids to include liquid mixtures across the full composition range. The fluctuations and integrals are then obtained for two binary mixtures of small molecules. It should be noted that the approach described here is not limited to binary mixtures and can actually be applied to any stable mixture regardless of the type, complexity, and concentration of the components and/or the state point. The approach is essentially the inverse of several expansions that have expressed higher order derivatives of thermodynamic properties in terms of higher order fluctuations and their distribution functions.^{13,27–30} However, previous approaches have been less general;^{13,24–27} as they were typically developed for osmotic systems and/or for infinite dilution only. Unlike the previously mentioned MM theory, which is only applicable to systems with very low solute concentrations,^{23} the FST-based approach presented here is relevant for any concentration. In fact, MM theory is actually a limiting form of FST.^{27}

The work of Bhatia and Ratti, and also Parrinello and Tosi, most closely parallels the current approach.^{31–35} They used a structure factor formulation in the long wavelength (thermodynamic) limit to provide higher order structure factors in terms of purely thermodynamic data. Their approach was then applied to a binary mixture of a liquid Na + K alloy to obtain the triplet fluctuations, albeit with a series of approximations for the thermodynamic data. In principle, this approach could be modified to study any liquid mixture. Unfortunately, to our knowledge, it has not been used to study mixtures of more common interest, e.g., mixtures composed of polyatomic molecules such as methanol + water. In our opinion, this is probably due to the use of a structure factor formulation, which obscures the relationship to the otherwise standard thermodynamic properties of liquid mixtures, and also (incorrectly) suggests that some form of scattering data is required. Here, we follow a similar path. The current approach differs, however, in that fluctuations and integrals are provided, no approximations are used, the theory is extended to include quadruplet correlations in a simple manner, the dependence of the results on the dataset and correlating equation is examined, and the theory is presented in terms of a general matrix formulation employing common thermodynamic quantities that can be easily extended to include any number of components. Furthermore, a simple physical model of the fluctuations is also provided.

One of the interesting general findings shown below is that FST indicates that the composition and pressure dependence of particle-particle fluctuations both involve the third cumulants of the particle number probability distributions. Specifically, if the third cumulants were zero (meaning the distribution was symmetric), the particle-particle fluctuations would be independent of both composition and pressure. Clearly, this is not the case and is in agreement with Greene and Callen’s acknowledgement that the distribution is not truly Gaussian,^{36} although that is a common assumption or approximation. Thus, in an effort to probe the underlying true distribution, we illustrate how to extract these cumulants from an analysis of the bulk thermodynamic experimental data for two binary mixtures over their full composition ranges, and then we show how these cumulants provide integrals over higher order distribution functions.

## II. THEORY

The most common application of Fluctuation Solution Theory, i.e., traditional Kirkwood-Buff (KB) Theory,^{27} relates second derivatives in the Grand Canonical ensemble (GCE) to second derivatives in closed isothermal isobaric Gibbs ensemble systems.^{29,37–43} Under isothermal conditions, the second derivatives in the GCE only involve the particle number fluctuations between all possible pairs of species. These fluctuations are then related to second derivatives of the Gibbs free energy for an equivalent closed system. The particle number fluctuations can also be related to integrals over pair distributions defined in the GCE. Here, we will extend the usual FST approach to include third and fourth isothermal derivatives in the GCE, and then relate these to third and fourth derivatives of the Gibbs free energy. The fluctuations, and the related integrals over triplet and quadruplet distribution functions, are then determined from available experimental data.

Second, third, and fourth isothermal derivatives with respect to the chemical potentials ({βμ}) in the GCE provide a series of particle number fluctuation densities,^{11,44,45} which we define by

where $\delta N \alpha = N \alpha \u2212 N \alpha $, etc., denote a fluctuation in the number of α particles, *N*_{α} is the instantaneous number of α particles for an open system of volume *V* at a pressure *p*, and β = (*RT*)^{−1} where *R* is the Gas constant and *T* the absolute temperature. The above quantities correspond to the cumulants of the multivariate particle probability distribution for the open system. They are also closely related to structure factors of various order.^{34,35,46} The *B*’s were already defined by Kirkwood and Buff.^{27} For a multivariate symmetric (e.g., Gaussian) distribution, the *C*’s and *D*’s would be zero. For real solutions this is not the case, and the *B*’s, *C*’s, and *D*’s can be used to quantify the pair, triplet and quadruplet correlations between the various particle types.

The fluctuating densities defined in Eq. (1) can be related to a series of corresponding distribution functions. For multicomponent systems, the relationships between the fluctuating quantities and the corresponding distribution functions are provided by^{24}

where *δ*_{αβ} is the Kronecker delta function. The above expressions involve integrals over the *n*-body spatial distribution functions $ g \alpha \beta \u2026 ( n ) ( r 1 , r 2 , \u2026 ) $ that are similar in form to the integrals appearing in the theory of imperfect gases or the McMillan-Mayer theory of (osmotic) solutions,^{23,27}

where the spatial dependences are implied and have been omitted for simplicity. Note that there is, by definition, no angular dependence of the integrals even for molecules. This is not an approximation. The integrals over the spatial distribution functions are valid for any liquid density and are obtained after averaging over the orientations of the molecules explicitly involved in the integral, and averaging over the positions and orientations of all the other molecules in the system.

Hence, if one can obtain the fluctuating quantities in terms of experimental data, then the corresponding integrals over the center of mass based two, three, and four body distribution functions can also be obtained. The fluctuation densities and corresponding integrals provide alternative, but complimentary, descriptions of any correlations that might occur between particles within the system.

The above expressions are restricted to open systems. The next step is to provide a connection to equivalent closed systems, which are of more common interest. The traditional application of FST involves the following second derivatives of the Gibbs free energy,^{27,47,48}

where κ_{T} is the isothermal compressibility, *V*_{m} is the molar volume, and $ V \u0304 \alpha $ is the partial molar volume of species α. In the original KB approach, expressions for the derivatives in Eq. (4) were provided in terms of the pair fluctuations given in Eq. (1), and thereby integrals over the pair distribution functions provided in Eq. (3).^{27} Alternatively, if one knows the composition dependent experimental properties displayed in Eq. (4), one can invert this approach and obtain the fluctuation densities and corresponding integrals. This is the so-called KB inversion approach.^{48} Here, we wish to provide higher fluctuations and correlations in terms of higher derivatives of the expressions in Eq. (4).

The analysis is simplified by defining a corresponding set of dimensionless fluctuating properties such that,

where ρ = ρ_{1} + ρ_{2} + ⋯ = 1/*V*_{m} is the total number density. We then note that the intensive GCE averages are a function of the intensive thermodynamic variables associated with the GCE, and therefore one can write the following general isothermal differential:

where the sum is over all *m* components in the liquid mixture, and *A* is an intensive property—specifically *ρ*, *ρ*_{α}, *B*_{αβ}, or *C*_{αβγ}. These relationships, coupled with the Gibbs-Duhem (GD) expression at constant temperature $ \u2211 \alpha N \alpha d\beta \mu \alpha =\beta Vdp$, are all that is required to extract the *B*’s, *C*’s, and *D*’s from the relevant thermodynamic data for liquid mixtures. Using Eq. (6), and the fact that ∂*ρ*_{α}/∂*p* = *ρ*_{α}*κ _{T}* and $\u2202 \rho \alpha /\u2202 x \beta =\rho ( \delta \alpha \beta \u2212 \rho \alpha V \u0304 \beta ) / ( 1 \u2212 x \beta ) $, one can write the following matrix expression for a binary solution of 1 and 2:

in terms of two symmetric dimensionless square matrices of size *m* + 1. One matrix contains just GCE fluctuating quantities

where *x*_{α} is the mole fraction of species α, while the other contains the corresponding isothermal experimental data,

where $ \mu \alpha \beta \u2032 =\beta N \u2202 2 G / \u2202 N \beta \u2202 N \alpha p , T = ( 1 \u2212 x \beta ) \mu \alpha \beta = ( 1 \u2212 x \alpha ) \mu \beta \alpha $. This corresponds to the traditional KB theory approach for a binary mixture. In this form, the theory involves just two matrices, albeit with one extra dimension. It should be noted that the determinants of the above matrices for binary solutions are given by $ a =\u2212 \mu 22 / x 1 =\u2212 \mu 11 / x 2 = \mu 21 / x 1 = \mu 12 / x 2 $ and $ b =\u2212 x 1 x 2 \eta /\rho $, respectively, where *η* = *ρ*_{1} + *ρ*_{2} + *ρ*_{1}*ρ*_{2}(*G*_{11} + *G*_{22} − 2*G*_{12}). Hence, one can show that Eq. (7) produces the same results as traditional KB theory.^{27,49} The extension to include additional components is straightforward.

To obtain the higher fluctuations in terms of experimental data for closed systems, we use the expressions in Eq. (6) and write them in matrix form as,

where

The final matrix involves derivatives of the dimensionless pair fluctuations

where the superscript *i* or *p* indicates an isothermal isobaric derivative with respect to *x _{i}*, or an isothermal derivative with respect to pressure at constant composition, respectively. The central row could be removed from the

**c**and

**b**matrices, if desired, as they provide additional expressions for the

_{x}*c*’s that are already contained in rows one and three. Solving Eq. (10) symbolically provides rather complicated expressions. Hence, we have simply solved the equations numerically in the current approach.

Finally, to obtain the quadruplet fluctuations one can write an additional matrix relationship

in terms of a matrix of fluctuations

and a matrix of triplet fluctuations and their derivatives

Clearly, one could continue indefinitely. However, it is unclear if the experimental data can support such manipulations.

In summary, the previous expressions illustrate exactly how to extract quantitative information concerning pair and triplet correlations from readily available experimental data. Traditional KB theory relates second derivatives in the Gibbs and GCEs and can be represented by the following scheme,

where the forward direction indicates the KB inversion procedure,^{48} while the backward direction corresponds to the results obtained in the original KB study.^{27} Here, we have extended this approach to include third derivatives in the Gibbs and GCEs as illustrated by the corresponding scheme

The forward direction corresponds to the analysis performed here and is equivalent to the KB inversion approach for triplet correlations. The reverse direction provides expressions for the corresponding thermodynamic data in terms of the triplet fluctuations and integrals. More details are provided in Appendix A and B, and emphasize the important role that triplet correlations play in determining the thermodynamic properties. Clearly, the latter derivatives can be used to rationalize the composition and pressure dependence of the chemical potentials and partial molar volumes, via simple series expansions, and illustrate how these are related to both the pair and triplet correlations.

## III. METHODS

To apply the above expressions, one requires derivatives of the pair and triplet fluctuations in terms of the experimental isothermal isobaric data. This involves additional composition and pressure derivatives of the original experimental data provided in Eq. (4). The experimental data of interest for traditional KB theory include the composition derivatives of the chemical potentials, usually obtained from derivatives of the excess Gibbs free energy of mixing, the partial molar volumes, usually obtained from derivatives of the excess volume of mixing, and the isothermal compressibility of the solution, all as a function of composition. For a binary solution, there are three unique *B*’s (and therefore *G*_{αβ}’s), four unique *C*’s, and five unique *D*’s. Useful relationships for binary mixtures are provided in Appendix B, and a more detailed description of the general approach for multicomponent systems is provided in Appendix C.

Here, we have chosen to study two binary mixtures. Mixtures of water + methanol were selected as there is substantial data concerning this system. In particular, multiple raw datasets and several correlating equations are available and have been studied in detail. Furthermore, all the required pressure derivatives for the system are also available. Mixtures of benzene + methanol were also examined. This system was chosen as our previous studies have shown interesting mixing behavior where methanol correlations lead to strong micro-heterogeneity at low methanol mole fractions.^{50,51} For this system, some of the pressure derivatives were approximated (see Appendix C). However, based upon the real values of these terms from the water + methanol analysis, the terms approximated for the benzene + methanol system are negligible under the ambient conditions studied here.

Sources of data for the water + methanol system at 298.15 K and 1 bar were as follows. Several data sets and fitting functions for the excess Gibbs free energy of mixing were explored; they are referred to in the accompanying figures by the following abbreviations: RK1, RK2, W1, W2, and W3. The specific datasets and fitting functions are: RK1, 3-parameter Redlich-Kister (RK) equation^{52} fit of the data from Butler *et al.*;^{53} RK2, 3-parameter modified-RK equation^{54} fit of the data from Hu *et al.*;^{54} W1, Wilson equation^{55} parameters provided by Soujanya *et al.*;^{56} W2, Wilson parameters provided on page 61 of Gmehling *et al.*;^{57} and finally W3, Wilson parameters provided on page 62 of Gmehling *et al.*^{57} All analyses were performed using the same density data. The pressure derivatives for the pure liquid volumes, $ \u2202 n V \alpha o /\u2202 p n $ with *n* = 1–3, were obtained from REFPROP version 9.1 (NIST Standard Reference Database 23, NIST Reference Fluid Thermodynamic and Transport Properties).^{58} The excess volume of mixing was fitted to a 4-parameter RK equation using the data of Douheret *et al.*^{59} The pressure derivatives of the excess volume of mixing, $ \u2202 n V m E /\u2202 p n $ with *n* = 1–3, were each fitted to 3-parameter RK equations using the data of Kubota *et al.*^{60}

Sources of data for the benzene + methanol system at 308.15 K and 1 bar were as follows. The excess Gibbs free energy of mixing was obtained from the fit provided by Wilson.^{55} The pressure derivatives for the pure liquid volumes, $ \u2202 n V \alpha o /\u2202 p n $ with *n* = 1–3, were obtained from REFPROP version 9.1.^{58} The excess volume of mixing was taken from the 4-parameter RK coefficients provided by Šerbanović *et al.*^{61} The higher pressure derivatives of the excess volume of mixing were all approximated as zero as explained above.

## IV. RESULTS

### A. Infinite dilution behavior

The main focus of this study is liquid mixtures at finite concentrations. However, it is useful to establish the limiting behavior of the fluctuations and integrals as the concentration of one species disappears in a binary mixture. As the concentration of species 2 tends to zero, the fluctuating quantities involving species 2 also tend to zero. The finite fluctuating quantities that remain only involve species 1. We then have the results described in our earlier study of pure liquids.^{24} Hence, the fluctuations are related to pressure derivatives of the pure liquid density according to

The corresponding *G*’s can then be obtained from Eq. (2). The above expressions also provide some interesting limiting cases. In the incompressible limit (IL), all the pressure derivatives are zero and consequently there are no fluctuations. This corresponds to the results expected for a closed system. The limiting behavior of the integrals is then given by Eq. (2) to provide

Hence, the integrals are only finite when all indices are identical. The integrals obtained here for real liquids do not follow this limit as they are defined for equivalent open systems. A second limiting case, the Gaussian limit (GL), occurs when the particle number fluctuations are assumed to be Gaussian or symmetric in nature. Here, the third (*VC*_{111}) and fourth (*VD*_{1111}) cumulants are then zero, and the higher integrals can be expressed in terms of the pair integrals using Eq. (2).

The results for the pure liquids studied here are presented in Table I as obtained from Eq. (16) and the two limiting cases. The pure liquids all display positive values for the pair fluctuations (required for all systems), with negative values for the triplet fluctuations and positive values for the quadruplet fluctuations. Negative values of *VC*_{111} indicate that a depletion of molecules in a reference volume of the liquid is favored over an addition, while positive values for *VD*_{1111} indicate that the corresponding particle number distribution is more peaked than the normal distribution.^{24} The signs change when examining the pair, triplet, and quadruplet distribution function integrals. However, the deviation of the integrals from the IL is consistent with the signs displayed by the fluctuating quantities. The GL appears to be a good approximation for pure water. However, this will not be true for liquid mixtures (see later). The trends observed between the various liquids follow the variation in the value of *b*_{11} (or ρ_{1}*RT*κ_{T}). Hence, benzene is closer to the IL, even though the compressibility of benzene is much higher than water because the number density of benzene is relatively low.

Property . | Methanol . | Methanol . | Benzene . | Water . | Water/GL . | IL . |
---|---|---|---|---|---|---|

T | 298.15 | 308.15 | 308.15 | 298.15 | 298.15 | |

ρ_{1} | 24.540 34 | 24.246 47 | 11.049 34 | 55.344 56 | 55.344 56 | |

pκ_{T}/10^{−5} | 12.640 | 13.542 | 10.434 | 4.525 | 4.525 | 0 |

b_{11} | 0.076 894 | 0.084 123 | 0.029 538 | 0.062 076 | 0.062 076 | 0 |

c_{111} | −0.051 96 | −0.062 38 | −0.007 426 | −0.014 19 | 0 | 0 |

d_{1111} | 0.106 | 0.138 | 0.004 87 | 0.005 61 | 0 | 0 |

ρ_{1}G_{11} | −0.923 107 | −0.915 877 | −0.970 462 | −0.937 924 | −0.937 924 | −1 |

ρ_{1}^{2}G_{111} | 1.717 4 | 1.685 3 | 1.904 0 | 1.799 6 | 1.813 8 | 2 |

ρ_{1}^{3}G_{1111} | −4.737 | −4.562 | −5.626 | −5.226 | −5.317 | −6 |

Property . | Methanol . | Methanol . | Benzene . | Water . | Water/GL . | IL . |
---|---|---|---|---|---|---|

T | 298.15 | 308.15 | 308.15 | 298.15 | 298.15 | |

ρ_{1} | 24.540 34 | 24.246 47 | 11.049 34 | 55.344 56 | 55.344 56 | |

pκ_{T}/10^{−5} | 12.640 | 13.542 | 10.434 | 4.525 | 4.525 | 0 |

b_{11} | 0.076 894 | 0.084 123 | 0.029 538 | 0.062 076 | 0.062 076 | 0 |

c_{111} | −0.051 96 | −0.062 38 | −0.007 426 | −0.014 19 | 0 | 0 |

d_{1111} | 0.106 | 0.138 | 0.004 87 | 0.005 61 | 0 | 0 |

ρ_{1}G_{11} | −0.923 107 | −0.915 877 | −0.970 462 | −0.937 924 | −0.937 924 | −1 |

ρ_{1}^{2}G_{111} | 1.717 4 | 1.685 3 | 1.904 0 | 1.799 6 | 1.813 8 | 2 |

ρ_{1}^{3}G_{1111} | −4.737 | −4.562 | −5.626 | −5.226 | −5.317 | −6 |

The properties of an infinitely dilute solute are also of considerable interest. The limiting forms of the pair integrals arising from the solution of Eq. (7) are given by

to order *x*_{2} and where *f*_{22} = ∂ln*γ*_{2}/∂*x*_{2}. Using these expressions in Eq. (A1) provides the following limiting expressions for the triplet integrals:

where the derivatives of the pair correlation integrals are given by

and the compressibility derivatives are,

### B. Symmetric ideal solutions

Symmetric ideal (SI) solutions provide a useful reference point for describing real solution behavior.^{47,65,66} Hence, we have also included the corresponding SI results for the two liquid mixtures studied here. SI solutions are characterized by $ G m E = V m E =0$ and give rise to the following general expression for the pair distribution integrals:^{65}

where $ \kappa T = \u2211 \alpha \rho \alpha V \alpha o \kappa T , \alpha o $. Using the appropriate derivatives of Eq. (22) in Eq. (A1) provide general SI expressions for the triplet integrals

The above expression is valid for SI mixtures containing any number of components at any composition. Both the pair and triplet integrals are finite and composition dependent for SI solutions. The values are determined by the volumes of the pure liquids and their pressure dependence.

The SI results for the two binary mixtures, water + methanol and benzene + methanol, are provided in Figures 1 and 2. Even for SI solutions, the distribution of molecules within a fixed volume of the mixture is not Gaussian in nature. Furthermore, the fluctuating quantities display a significant dependence on composition. This appears to be driven by the changes in the number densities, as the corresponding *G*’s display a far simpler dependence on composition. The *b*’s display a single maximum or minimum, the *c*’s display two, and the *d*’s display three. This is to be expected, as they are essentially derivatives of each other. Mathematically, the pair, triplet, and quadruplet fluctuations are second, third, and fourth order polynomials in the composition, respectively. This helps to explain the observed oscillating behavior of these properties. Physically, one can also explain some of these variations using the approach described in Sec. IV F. The *G*’s generally alternate in sign as one goes from *G*_{αα} to *G*_{ααα} to *G*_{αααα} for all compositions.

### C. The effects of different models and experimental datasets

The theory presented in Sec. II is exact. Hence, there are no assumptions or approximations in the equations leading to the triplet and quadruplet fluctuations and integrals. However, the results may be sensitive to the quality of the experimental data used in the analysis. In particular, it is well known that the pair fluctuations and integrals can depend significantly on the quality of the raw excess Gibbs free energy of mixing data and the type of correlating equation used to represent these data.^{42} Hence, it is important to test the degree to which the experimental data provide reproducible values of the cumulants and integrals.

The results of such a test are shown in Figures 3 and 4 for water (1) + methanol (2) mixtures. This system was specifically chosen as there are many determinations of $ G m E $ available. For example, at the time this analysis was performed, there were 126 total water + methanol binary vapor-liquid equilibrium (*T*-*p*-*x*-*y*) determinations (from which $ G m E $ is obtained) in the NIST ThermoData Engine^{67–69} within Aspen Properties^{®}, part of the Aspen Plus Product Family within the aspenONE^{70} Engineering CAD software from Aspen Technology, Inc. Five examples are described in Sec. II and are presented in Figures 3 and 4. The $ G m E $ determinations chosen here do not necessarily pass all the possible thermodynamic consistency tests.^{67,71–78} For the present study, we decided not to eliminate data sets based upon the results of these tests, primarily due to the controversy and lack of clarity concerning how exactly these tests should be performed, and because many of the tests typically require the use of a model (such as Wilson or NRTL) such that they are no longer purely tests of the raw data itself.^{76,78} Furthermore, for most other mixtures, there are going to be significantly fewer raw *T*-*p*-*x*-*y* data sets available from which to obtain $ G m E $. Thus, the seemingly *ad hoc* choices made here, while not fully satisfactory, are actually more representative of the options that might be available for most other systems. Nevertheless, a thorough investigation of the dependence of the water + methanol results on the raw data and correlating equation, incorporating thermodynamic consistency checks and using more data sets and correlating equations, is certainly called for in the future. We anticipate that less disparity between the results would be observed if only thermodynamically consistent data set/correlating equation combinations were used for comparison. Therefore, the results in Figures 3 and 4 may actually tend toward the worst-case scenario.

In Figures 3 and 4, there are clearly quantitative differences between the experimental data sets, especially when using different fitting functions. In particular, use of the RK equation typically resulted in more extreme variations with composition. Due to the well-known recommendation to use the Wilson equation for mixtures containing alcohols,^{79–81} and due to the absence of a physical meaning for the RK parameters, we anticipate that the Wilson analyses are better representations of reality than are the RK analyses. However, the trends in the properties with composition appear to be reasonably clear for all five of the determinations, and consistent results are obtained for the *B*’s and *C*’s for most compositions. The *D*’s clearly show the largest variations between data sets and models. It does appear that the present type of analysis can provide reasonably reliable data for the *B*’s and *C*’s, while care would have to be used when drawing conclusions concerning the *D*’s. A similar story is found for the corresponding *G*’s provided in Figure 4. Clearly, the *G*’s also become less reliable as the relevant species’ concentration approaches infinite dilution.

We have also investigated the effect of using ideal mixture data for both the compressibility and volume of mixing (data not shown). The assumption of an ideal compressibility resulted in essentially no changes to the fluctuations or integrals, while the assumption that $ V m E =0$ resulted in visible differences in the data, but these effects were small compared to the differences observed on changing the dataset and/or correlating equation. It should be noted that these results are probably general for most liquid mixtures under ambient conditions.^{42}

### D. Water and methanol mixtures

Further examination of the results for water + methanol mixtures presented in Figures 3 and 4 leads to several important observations. The results are only slightly different from those predicted for a SI mixture. Indeed, it could be argued that the difference from SI is less than the dispersion observed between datasets and correlating equations. The magnitude of the fluctuating quantities decreases as the number of methanol indices increases. This appears to be due to the higher number density for water, which contributes significantly to the fluctuations, as the corresponding integrals display little change or even a slight increase in magnitude as the number of methanol indices increases.

The fluctuating quantities oscillate in value with composition. A maximum or minimum in the *b*’s occurs at a composition where the *c*’s are close to zero and the *d*’s display a maximum or minimum. Again, this is to be expected from the relationships in Eqs. (1) and (6), and this issue is elaborated upon further in Sec. IV F. The water pair fluctuations (*b*_{11}) are always positive, as are the methanol pair fluctuations (*b*_{22}). The water triplet fluctuations (*c*_{111}) start off negative at low methanol mole fractions and then turn positive as the methanol mole fraction dominates. The opposite trend is observed for the methanol triplet correlations (*c*_{222}). The magnitude of the fluctuations is clearly largest for fluctuations that contain a larger number of water molecules.

The corresponding integrals are displayed in Figure 4. It is clear that fairly consistent results are obtained for the *G*_{11}, *G*_{12}, *G*_{111}, *G*_{112}, *G*_{1111}, and *G*_{1112} integrals at low *x*_{2} compositions. The same is true for the corresponding (1↔2) integrals at low *x*_{1} compositions. As one reaches infinite dilution of one component, the results for integrals containing that component depend substantially on the dataset and correlating equation. Hence, care must be taken when analyzing data in this region. The two limiting behaviors, IL and GL, only appear to be approximately valid for pure liquids, e.g., for ρ^{2}*G*_{111} when *x*_{2} approaches zero, and not for the integrals containing mixed species at finite concentrations (data not shown).

### E. Benzene and methanol mixtures

The results for mixtures of benzene (1) + methanol (2) are displayed in Figures 5 and 6. This mixture clearly deviates substantially from that predicted by SI behavior. In particular, one observes a large positive correlation between methanol molecules, indicated by the large *b*_{22} and *G*_{22} values, for low mole fractions of methanol. Our previous studies have suggested that this indicates a significant deviation from a random distribution of methanol molecules, attributed to the strong desire to form intermolecular hydrogen bonds, as the environment becomes more non-polar.^{50,51} This large positive correlation between methanol molecules extends from the pair to the triplet and quadruplet fluctuations and integrals. The magnitude of the fluctuations and integrals increases as the number of methanol indices increases. Hence, not all the triplet terms will contribute equally to the thermodynamic properties provided by Eqs. (B2)–(B5). In particular, the *c*_{222} and *c*_{221} terms clearly dominate for benzene + methanol mixtures at low methanol mole fractions. The symmetry between the different fluctuations and integrals observed for the SI and water + methanol mixtures is also present for benzene + methanol.

### F. General features of the fluctuating quantities

The expressions provided in Eqs. (1) and (6) also indicate the conditions for which there is a maximum or minimum in the *B*’s or *C*’s as a function of composition, as is commonly observed. Using these expressions, followed by application of the GD expression for binary solutions at constant temperature and pressure provides

where *m*_{2} = *ρ*_{2}/*ρ*_{1}. The chemical potential derivative has to be positive or negative for stable solutions and hence the derivative of *B*_{αβ} can only be zero when

where we have made the substitution, *n _{i}* =

*N*/〈

_{i}*N*〉. Hence, for binary solutions, the maximum or minimum in

_{i}*B*

_{αβ}occurs when the corresponding (reduced) triplet fluctuations are equal for both components. Similar relationships can be obtained for maxima or minima for the

*C*’s in terms of the

*D*’s. Conditions that lead to a maximum or minimum in the

*G*’s can also be developed, but are significantly more complicated and also appear less informative.

A somewhat reasonable physical picture of the two body (and higher) particle number fluctuations can be developed for binary systems. In Figures 1, 3, and 5, the SI and experimental *b*_{11} and *b*_{22} values appear correlated, while the *b*_{12} value appears to be anti-correlated with the other *b*’s. In pure liquids of component 1, the *b*_{11} values describe the fluctuations of *N*_{1} particles in and out of a fixed volume of liquid. The magnitude of these fluctuations is determined by the size of the region and the compressibility of the pure liquid, which is normally small under ambient conditions. When component 1 is very dilute, the value of *b*_{11} tends to zero simply because there are very few *N*_{1} molecules in a given region of the liquid. For intermediate compositions, one can observe particle number fluctuations due to the compressibility of the liquid (Mechanism 1), but also due to exchange of 1 and 2 molecules such that the volume of liquid remains essentially constant for all compositions (Mechanism 2),

Hence, the magnitude of the fluctuations in mixtures can increase dramatically due to this second mechanism. This helps to explain why the *b*_{11} (and *b*_{22}) fluctuations are essentially anti-correlated with the *b*_{12} values. This correlation/anti-correlation relationship also persists for the higher order fluctuations involving pairs of *c*’s and *d*’s that only differ by a single species index (*c*_{111} and *c*_{112}, for example). Furthermore, if one multiplies Eq. (26) by *δN*_{α}*δN*_{β}/*V* and then ensemble averages, one finds the approximate relationship $ C \alpha \beta 1 V \u0304 1 + C \alpha \beta 2 V \u0304 2 \u22480$, which is true when $RT B \alpha \beta p $ is small. Using the condition given by Eq. (25) implies that not only are the relevant fluctuating quantities equal at the maxima and minima, but they are both close to zero. This type of behavior is observed in Figures 1, 3, and 5.

## V. CONCLUSIONS

We have provided quantitative data concerning the role of triplet and higher correlations for complex liquid mixtures from experiment. The above approach avoids the need for information from scattering data, which is difficult to obtain for complex molecules and/or liquid mixtures. Therefore, data for a wide range of liquid mixtures should be relatively easy to access, albeit at the expense of an absence in spatial resolution. The initial systems studied here suggest the data are sensitive to the correlating equation used for the excess properties, although the trends with composition remain consistent. The correlations vary significantly with composition, even for symmetric ideal solutions, and the degree of variation steadily increases the higher the correlation function. The fluctuations and integrals appear to be dominated by correlations between the most polar molecules in the system. This is to be expected. However, and more importantly, these correlations can now be quantified using the approach presented here. It is envisioned that the type of data provided here will enable more rigorous tests of models and theories for liquids and their mixtures, will help with the development of improved correlating equations, and will provide insights into the pressure, temperature, and composition dependence of the pair correlations in a variety of systems. In particular, extra confidence in theoretical and simulation studies of the triplet correlation functions, $ g \alpha \beta \gamma ( 3 ) $, for example, can be provided when one has agreement with experiment for the corresponding integrals.

## Acknowledgments

The project described was supported by Grant No. R01GM079277 from the National Institute of General Medical Sciences to PES. EAP was supported by the NSF GRF program under grant NSF DGE-0750823. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institute of General Medical Science, the National Institutes of Health, or the National Science Foundation.

### APPENDIX A: GENERAL FST RELATIONSHIPS FOR MULTICOMPONENT SOLUTIONS

In this appendix, we outline a few of the expressions resulting from the theory section in an attempt to clarify what has been accomplished and how this has been achieved. The following general relationships are apparent from Eqs. (10) and (13):

for any number of components. The first expression summarizes the relationships derived by Bhatia and Ratti for the third order structure factors at long wavelengths in binary mixtures.^{34} The expressions for the *b*_{αβ} quantities depend on the number of components, although a general expression can be written in terms of matrix cofactors as indicated by Eq. (7). Nevertheless, the pair fluctuations can be expressed in terms of the experimental data provided in Eq. (4), and the above expressions indicate how additional pressure and composition derivatives lead to the triplet (and quadruplet) fluctuations. The expressions clearly involve many terms, even for binary mixtures. Unfortunately, at present it is unclear which terms may contribute the most to the final results.

The theory also provides FST relationships for several thermodynamic quantities related to third derivatives of the Gibbs free energy. In particular, if one takes derivatives of Eq. (7) one finds

where ∂**a**/∂*X* indicates a matrix of partial derivatives containing the derivative of each of the original elements with respect to *X*, etc., and *X* = *N*_{γ} or *p*. This matrix expression can be solved to provide general expressions for the thermodynamic properties of any multicomponent mixture given by derivatives of the **a** matrix. Derivatives with respect to composition provide the composition-composition-composition (*N*-*N*-*N*) derivatives

where the sums are over all components.

Derivatives with respect to pressure provide a series of relationships. The *N*-*N*-*p* derivatives are given by

which also provide the composition dependence of the partial molar volumes through the relationship $ \u2202 \mu \alpha \beta \u2032 / \u2202 p T , { N} =\beta N \u2202 V \u0304 \alpha / \u2202 N \beta p , T , { N} \u2032 $. The *N*-*p*-*p* derivatives are provided by

which also provide the composition dependence of the compressibility through the relationship $ \u2202 \rho V \u0304 \alpha / \u2202 p T , { N} = \kappa T ( 1 \u2212 \rho V \u0304 \alpha ) \u2212V \u2202 2 \rho / \u2202 p \u2202 N \alpha T $. Finally, the *p*-*p*-*p* derivative is given by

To develop these relationships further requires derivatives of the pair fluctuations obtained from using Eq. (6) to provide

Hence, the pressure and composition dependence of the pair fluctuations are related to the corresponding triplet fluctuations. Using these expressions in Eqs. (A3)–(A6) provides the final results for the thermodynamic derivatives in terms of pair and triplet fluctuations, which we write in the most condensed form as

In the above expressions, we have chosen to use the second derivative quantities provided by the pair fluctuations for simplicity, and because the sign of the quantities is then easier to determine. The above expressions clearly indicate the role of triplet correlations in determining the thermodynamic properties of liquid mixtures. The above relationships can also be used as the starting point for providing additional derivatives in terms of the quadruplet fluctuations, if desired.

Finally, it should be noted that the relationship between the two sets (primed and unprimed) of chemical potential derivatives is provided by

and the above expressions obey the following relationships:

which can be obtained from the GD equation at constant *T*.

### APPENDIX B: FST RELATIONSHIPS FOR BINARY SOLUTIONS

In this appendix, we focus on the most common mixture of interest, namely binary mixtures. For binary mixtures, we have general expressions for the pair fluctuations in terms of experimentally accessible data. Furthermore, the GD expression allows one to focus on a single chemical potential derivative. Hence, Eq. (A1) can be written as

Substitution of the derivatives of the first expression into the second expression, etc., could be performed. However, there are six terms for the pressure derivative and eight terms for the composition derivative of *b*_{αβ} and this route is therefore undesirable. Hence, the approach outlined below in Appendix C is to be preferred. However, in our experience the composition derivatives of the excess Gibbs free energy will dominate the results under ambient conditions.

The expressions for the triplet based thermodynamic properties relevant for binary mixtures can also be simplified by use of the GD expression to give

Other derivatives can be obtained from the relationships provided in Eqs. (A12) and (A13). The relationships for the thermodynamic properties for binary mixtures provided by the pair correlations are given by^{47}

and can be used in Eqs. (B2)–(B5) if desired. The above quantities are typically positive, although the partial molar volume can be negative in some cases, and hence the sign of the above derivatives is directly related to the sign and magnitude of the triplet correlations.

### APPENDIX C: ANALYSIS OF EXPERIMENTAL DATA

For a general multicomponent mixture, the molar Gibbs free energy of mixing can be written as $ G m = \u2211 \alpha x \alpha \mu \alpha = G m E + \u2211 \alpha x \alpha \mu \alpha id $, while the molar volume can be written as $ V m = V m E + \u2211 \alpha x \alpha V \alpha o $. The latter expression also provides the required relationship for the isothermal compressibility. Experimental excess Gibbs and volume of mixing data are typically provided in terms of a correlating equation (Wilson, Redlich-Kister, NRTL, etc.). Given such an equation, the corresponding excess partial molar quantities are provided by

Here $\beta \mu \alpha E =ln \gamma \alpha $, where *γ*_{α} is the mole fraction scale activity coefficient. Using the correlating equation, one can then obtain higher derivatives of both the total chemical potentials and the partial molar volumes given that ∂*x*_{α}/∂*x*_{β} = (*δ*_{αβ} − *x*_{α})/(1 − *x*_{β}). The corresponding ideal contributions $\beta \mu \alpha id =\beta \mu \alpha o +ln x \alpha $ to the chemical potential derivatives are given by

for any number of components. In addition, we have also used the following relationships for the cross derivatives:

together with, $\u2202 G m E /\u2202p= V m E $, all of which are valid for any liquid mixture.

As noted above, in the present analysis of the experimental data for binary solutions, we have made one minor approximation for the mixture of benzene + methanol. It has been assumed that

when *n* > 1. This is equivalent to neglecting the pressure dependence of the excess molar volume or excess partial molar volumes. A comparison of the data provided in the Sec. IV indicates that the magnitude of the *c*’s obtained for liquid mixtures are typically much larger than for the pure liquids. Furthermore, neglect of these derivatives for the water + methanol system indicated very minor contributions from these terms. Hence, this approximation should be very reasonable under ambient conditions and will likely be needed for other systems, as the pressure dependence of the solution volume as a function of composition is rarely available in the literature.

It is possible to obtain the derivatives of the fluctuating quantities in terms of the derivatives of the experimental data using the expressions provided in Eq. (6). However, this quickly becomes rather tedious. A simpler, but more opaque, alternative has been adopted here. The first derivatives of the pair fluctuations can be obtained from Eq. (7) by taking derivatives with respect to *X* = *x*_{1}, *x*_{2}, or *p*. This can be continued to provide second derivatives of the pair fluctuations such that

Hence, the first and second derivatives of the *b*’s can be expressed in terms of derivatives of a single matrix, **a**, that contains just experimental data. We note that ∂**b**/∂*X*≠**b _{X}**, but the results from Eq. (C5) can be used to obtain the matrix elements of

**b**. This approach can also be applied to obtain derivatives of the triplet fluctuations, as required by Eq. (15), to give

_{X} where *Y* = *x*_{1}, *x*_{2} or *p*. These manipulations do not necessarily provide clarity to the expressions developed in Sec. II. However, they are relatively easy to implement computationally.

Finally, the fluctuating quantities obtained from the above approach have been verified by comparing each pair fluctuation to the integral over the corresponding triplet fluctuations provided by Eq. (6), and comparing each triplet fluctuation to the integral over the corresponding quadruplet fluctuations also provided by Eq. (6).

## REFERENCES

*CKCSST’98, Qingdao, China*(1998) pp. 30-36.