J. Phys. Chem. B 2001, 105, 9581-9585
9581
Nonideality in Binary Mixtures: Correlations between Excess Volume, Excess Viscosity, and Diffusion Coefficients Arnab Mukherjee and Biman Bagchi* Solid State and Structural Chemistry, Unit Indian Institute of Science, Bangalore-12, India ReceiVed: April 9, 2001; In Final Form: July 3, 2001
Although many experimental studies over the years have shown strong (inverse) correlation between excess viscosity and excess volume of a binary mixture, there does not seem to exist any microscopic theoretical study of such behavior. In this work, we introduce and study two new models of binary mixture with the aim of removing this lacuna. Our models consist of a mixture of two molecular species (1 and 2) with the same diameter and mass but varying solute-solvent Lennard-Jones interactions. In model I, the two different species are strongly attractive, while in model II, the attraction is weaker than that between the pure components. That is, model I encourages structure formation, while the reverse is true for model II. Extensive constant pressure-constant temperature (NPT) molecular dynamics (MD) simulations have been carried out which show that these models can indeed reproduce the strong inverse correlation between the calculated excess volume and excess viscosity. The two models show opposite trends in the composition dependence. Another interesting result is the observation of a re-entrant behavior in the excess volume dependence of total viscosity of the mixturesthe same re-entrance is also observed in the dependence of excess viscosity on excess volume. In addition, we derive the first two terms in the virial expansion for the composition dependence for the infinite frequency shear modulus, G∞. These terms can explain the composition dependence at low composition in terms of partial radial distribution functions and interaction potentials. The calculated self-diffusion coefficients have nonmonotonic behavior and show considerable deviations from Stokes-Einstein relation.
1. Introduction Binary mixtures constitute very important class of solvents which have drawn considerable attention of theoreticians and experimentalists1-9 alike over many years. These mixtures show diverse behavior, including many which can be justly termed anomalous. At a microscopic level, the two components of the mixture consist of molecules which have different molecular sizes and shapes and also different interaction strengths. These properties allow one to achieve the desired tunability in binary mixture solution properties. If we make the simplifying assumption that molecules of the mixture can be approximated as spheres, even then two diameters (σ1 and σ2) and three interaction energies (11, 22, 12) are needed to model the interaction potentials involved in a binary mixture. Therefore, when compared with one component systems, it is not surprising that the binary mixtures show wide variety of behaviors which are difficult to understand within one single, consistent, microscopic picture. In fact, both thermodynamic and transport properties of binary mixtures have remained one of the least understood problems of liquid-phase chemistry. One such oft-discussed (even in the text books of chemistry) anomalous property is nonideality. Nonideality is best described by deviations from the ideal behavior given by Raoult’s law,8 well-known in classical physical chemistry. For a given property, P, of 1 mol of mixture, the latter predicts the following simple dependence on the composition:
Pid ) x1P1 + x2P2
(1)
where xis are the mole fractions and Pis are the values of the property P of the pure (single component) liquids. Nonideality
or deviation from Raoult’s law is decsribed by an excess function, Pex defined by
Pex ) P - (x1P1 + x2P2)
(2)
Most discussed properties are the volume, Vex, specific heat, and the viscosity, ηex. For example, following the above equation, excess volume Vex is defined by
Vex ) V - (x1V1 + x2V2)
(3)
There exist a large body of experimental evidences of the inverse correlation between excess viscosity and excess volume of the liquid mixture. Thus, if the excess volume is positive, then the excess viscosity becomes negative, and vice versa.10,11 It has also been observed that when the constituent species (say, 1 and 2) are like each other, then the excess volume is negative, and excess viscosity is positive, and vice versa. In addition, Vex and ηex show strong composition dependence, with a pronounced extremum, not always at x1 ) x2 ) 0.5. In recent years, several theoretical and computer simulation studies on Lennard-Jones binary mixtures have been carried out.12,13 However, these studies have concentrated mainly on the glass transition in binary mixtures which are known to be good glass formers. In addition, these studies have considered only one particular composition and a unique interaction strength. The absence of any microscopic study of the anomalous composition dependence of viscosity is not surprising because a microscopic calculation of viscosity is quite difficult. In the absence of any microscopic theory, the experimental results have often been fitted to several empirical forms. Prominent among
10.1021/jp011313z CCC: $20.00 © 2001 American Chemical Society Published on Web 09/08/2001
9582 J. Phys. Chem. B, Vol. 105, No. 39, 2001
Mukherjee and Bagchi
them is Eyring’s theory of viscosity extended to treat binary mixtures.14 However, the very basis of Eyring’s theory has been questioned, as this theory is based on creation of holes of the size of the molecules which is energetically unfavorable. Here we have studied two models (referred to as I and II) of binary mixtures in which the solute-solvent interaction strength is varied by keeping all the other parameters unchanged. All the three interactions are described by the modified LennardJones potential. In model I, specific structure formation between solute and solvent molecules is mimiced by strong solutesolvent attractive interaction. This interaction is stronger than that between solvent-solvent and solute-solute interactions. The second model (model II) involves structure-breaking by weak solute-solvent interaction. These two models are perhaps the simplest model to mimic the “structure-making” and “structure-breaking” in binary mixtures. For convenience, we denote the solvent molecules as 1 and the solute molecules as 2. In both the models, 1 and 2 have the same radius and the same mass. We have shown elsewhere that these models can reproduce the observed trend in the composition dependencies of the excess viscosity.15 But no study has been done on excess volume because all the previous simulations were carried out at constant volume. In this work, we carry out extensive constant numberpressure-temperature (NPT) molecular dynamics simulations to evalute the composition dependence of excess volume and viscosity of both the models I and II. Note that only a constant pressure simulation can address the problem of excess volume. We find that model I shows a pronounced negative deviation in excess volume and a concomitant postive deviation in excess viscosity from ideality at the intermediate composition, precisely of the type observed in many experimental situations. Model II, on the other hand, shows a positive deviation in excess volume and negative deviation in viscosity from ideality. In addition, we have formulated a virial expansion of G∞ in the solute composition where we derived the microscopic expression for the first two terms. The virial expansion can explain the signature of deviation very well at low composition changes, which suggets that one can indeed propose a simple explanation of the nonideality of viscosity in terms of interactions among the same and between the two species. The simulation and the theoretical studies suggest that much of the observed behavior can be explained in terms of the “structure-making” or “structure-breaking” tendencies of liquid mixtures arising from the mixing of the two constituent species. Thus, the anomalies have their origin in the interaction energies, most importantly in the value of 12. This is evident from the fact that the Lorentz-Berthelot mixing rule fails to reproduce the anomalous composition dependence of viscosity or diffusion. Organization of the rest of the paper is as follows. In the next section, we present the basic definitions and the main equations. In section 3, we present both the simulation details and the models used in this study. Detailed description of the virial expansion is given in section 4. Section 5 contains the results and discussion. We close the paper with a few conclusions in section 6. 2. Basic Definitions 2.1. Viscosity. Microscopic expression for the time-dependent shear viscosity is formulated in terms of stress autocorrelation function and is given by16,17
η(t) ) (VkBT)-1〈σxz(0)σxz(t)〉
(4)
where σxz is the off-diagonal element of the stress tensor N
σxz )
[(pjxpjz/m) + Fjzxj] ∑ j)1
(5)
Here Fjz is the z-component of the force acting on the jth particle and the corresponding position of the jth particle is xj, and pjz is the z-component of the momentum of jth particle, m being the mass of the particle. Among the total N number of particles present in the system, N1 is the number of solvent particles, and N2 is the number of solute particles, where N1 + N2 ) N. High-frequency shear modulus is given by16,17
G∞ ) (VkBT)-1〈(σxz(0))2〉
(6)
Frequency-dependent viscosity is obtained by Laplace transforming η(t)
η(z) )
∫0∞ dt exp(-zt)η(t)
(7)
Experimentally observed viscosity is given by the zero frequency limit of η(z). 2.2. Self-Diffusion. Self-diffusion is the property of a single tagged particle which is obtained from the velocity autocorrelation function18
Dsi )
∫0∞ dt 〈vi(0)‚vi(t)〉
1 3
(8)
where vi is the velocity of the ith particle. Another definition of self-diffusion is by mean square displacement (MSD) given by Einstein is
Dsi ) lim tf∞
1 〈|r (t) - ri(0)|2〉 6t i
(9)
Diffusion Dis of the ith species is related to friction ζi by Einstein’s relation
Dsi )
kB T mζi
(10)
Stokes-Einstein expression is obtained when the friction is approximated by
ζi ) CηRi
(11)
where η is the viscosity of the liquid, Ri is the radius of atom i, and C is the constant determined by the hydrodynamic boundary condition. 3. Simulation Details We have carried out a series of molecular dynamic simulations at constant pressure (P), temperature (T), and total number of particles (N) of binary mixtures by varying the solute mole fraction from 0 to 1.19-21 Pressure is kept constant by Anderson’s piston method, while in the case of temperature, a damped oscillator method has been adopted which keeps temperature constant at each step.19 The piston mass involved here is 0.0027(m/σ4), which is regarded as optimum.19 Our model binary system consists of a total 500 (solvent (1) + solute (2)) particles. Instead of using simple Lennard-Jones 12 - 6 potential, interaction between any two particles is given by the modified Lennard-Jones potential which sets a cutoff radius rc outside which the potential energy is 0. The particular form of the
Nonideality in Binary Mixtures
J. Phys. Chem. B, Vol. 105, No. 39, 2001 9583
potential is given by22
{[( ) ( ) ] [ ( ) ( ) ] ( ) ( )}
Uij ) 4ij
σ rij
12
-
σ rij
6
+ 6
σ rc
12
-3
σ rc
6
σ rc
12
+4
7
(r/rc)2 σ rc
6
(12)
where the cutoff distance rc in this particular case has been taken as equal to 2.5σ. Use of above potential form takes care of the fact that both potential and force are continuous at the cutoff distance. i and j denote two different particles. We set the diameter (σ) and mass (m) of both the solute and the solvent to unity, for simplicity. The solute-solvent interaction strength lies in the potential well depth 12, where 1 and 2 represent the solvent and solute particles, respectively. Throughout this study, we keep the interaction strengths 11 ) 1.0 (solvent-solvent) and 22 ) 0.5 (solute-solute). We dealt with two very different specific solvent-solute interaction strength values, referred as model I and model II. In model I, 12 ) 2.0, and in model II, 12 ) 0.3. While the former accounts for the situation in which the solute solvent attract each other stronger than they do to their own species, the later describes the opposite scenario. Henceforth, we will be referring the situation in which 12 ) 2.0 and 12 ) 0.3 as attractive and repulsive systems, respectively. For both the models, we have set the pressure equal to 2.0(/σ3). We set the reduced temperature T* () kBT/) as 1.0 °C in model I and 1.24 °C in model II. After many trial runs to verify the existing results on viscosity23 of one component liquids, we selected a time step ∆t* ) 0.002τ for model I and ∆t* ) 0.001τ for model II for the integration of Newtonian equations of motion. The scaled time has been denoted as τ ) σxm/. We have dealt with eight different solute compositions, namely 0.0, 0.1, 0.2, 0.4, 0.6, 0.8, 0.9, and 1.0. For each solute composition, we have equilibrated the system up to 1.5 × 105 time steps. Simulations carried out for another 1 × 106 steps after equilibration during which the stress tensor has been calculated. We have done three runs for each composition point where each one of the last two runs have been started from the particle positions and velocites stored in the last time step by the previous simulation. We have also calculated the radial distribution function in each case to make sure that the clustering and phase separation (especially among the similar species) is avoided. 4. Virial Expansion of Excess Viscosity in Binary Mixtures
2π 15
2
∑
i,j)1
FiFj
d
∫0∞ dr gij(r) dr
1 G∞(x2) ) G∞(x2 ) 0) + G′∞(x2 ) 0)x2 + G′′∞ (x2 ) 0)x22 + 2 ... (14) where G′∞ ) (d/dx)G∞(x2) and G′′∞ ) (d2/dx2)G∞(x2). With the help of eq 13, the virial expression can be written as
1 G∞(x2) ) G∞(x2 ) 0) + [A12 - A11]x2 + [A11 - 2A12]x22 2 (15) where
1 G∞(x2 ) 0) ) FkBT + A11 2
(16)
and
Aij )
4πF2 15
∫0∞ dr gij(r) drd
[
r4
]
dVij(r) dr
(17)
Equation 15 is a new virial expansion presented here for the first time. As expected, the coefficients of this expansion are integrals over the partial radial distribution functions (eqs 16 and 17). We now analyze the above virial expressions to probe into the solute dependence of G∞. As the virial expansion is valid for small solute composition (x2), the nature will depend on Aij’s, which in turn depend on radial distribution functions and interaction potentials. In the case of model I, A12 is greater than A11, and the opposite goes for model II. So, from eq 15, we can see that with addition of small solute (small x2), G∞ will show an increment for model I and decrement for model II over pure solvent. On the other hand, if we can consider any of the ingredients of binary mixture as a solute, the above trend will be observable from both end of the solute axis. Simulations results also support the above prediction. 5. Numerical Results and Discussion
Virial expansions are very useful for understanding the effects of interactions. These expansions are valid at low concentration. Here we have formulated a virial expansion of the highfrequency shear modulus (G∞) against composition starting from pure liquid. The coefficients of the expansion are derived from microscopic definitions presented earlier. These coefficients give correctly the deviation of G∞ of pure solvent under the addition of solute in low concentration and vice versa. With the knowledge at hand about the radial distribution function or the interaction potential, we show that one could predict from the virial expansion the trend of the excess high-frequency shear modulus. Using eq 5 and eq 6, the expression for G∞ for binary mixture can be written as15
G∞ ) (F1 + F2)kBT +
where i, j ) 1 indicates solvent particles and i, j ) 2 denotes solute particles. So F1 is the number density for the solvent particles, and F2 denotes the same for the solute particles. gij(r) is the radial distribution function between particles labeled i and j. Note that Vij includes three different interaction potentials present between the solute and the solvent molecules. The virial expansion of G∞ can be written as
[ ] dVij(r)
r4
dr (13)
In panels a and b of Figure 1, we plot the composition dependence of the excess volume (Vex) and excess viscosity (ηex), respectively, for model I. Note the striking (and the opposite) nonmonotonic dependencies of ηex and Vex on solute composition for the two cases. In Figure 2, we plot ηex against Vex for the same model where the variation is obtained by changing composition. This shows a re-entrant behavior which essentially describes the correlation between the excess volume and the excess viscosity in a dramatic fashion. The composition dependencies of excess volume and excess viscosity are plotted for model II, respectively, in panels a and b of Figure 3. Model II shows completely opposite composition dependence (to that of model I). Note again the nonmonotonic dependencies of excess viscosity and excess volume on composition for both the two cases. The results presented in Figures 1-3 can be partly understood by analyzing the microscopic structure. In panels a and b of Figure 4, we plot all the partial radial distribution functions
9584 J. Phys. Chem. B, Vol. 105, No. 39, 2001
Figure 1. Composition dependence of excess viscosity and excess volume for model I. In panel a, the calculated excess viscosity (ηex) is plotted against solute compostion (x2) for model I. In panel b, the variation of excess volume (Vex) with solute composition (x2) is plotted for the same model. The solid lines are to guide the eye. Error bars are given.
Mukherjee and Bagchi
Figure 3. Composition dependence of excess viscosity and excess volume for model II. In panel a, the excess viscosity (ηex) is plotted against solute composition (x2) for model II. In panel b, excess volume (Vex) is plotted against the solute coposition (x2) for model II. Again solid lines are only for guidelines. Error bars are given.
Figures 5 and 6 show the composition dependence of selfdiffusion coefficients (Dsi ), for models I and II, respectively. These values have also been obtained by using the present (NPT) simulations. Self-diffusion coefficients also show nonmonotonic dependencies on solute composition. 6. Conclusion
Figure 2. Re-entrant nature of excess viscosity (ηex) is plotted against excess volume (Vex) for model I. Solid line is to guide the eye, and the arrows give the direction of the increase in solute composition, showing the re-entrant behavior. Error bars are given.
gij(r) obtained from the simulation for models I and II, respectively. As the solute-solvent interaction strength affects the structure surrounding a solute/solvent to a great extent,24 the above observed features are reflected in the increment of correlation among the unlike species, g12(r), in Figure 4a for model I, while the reverse is seen in Figure 4b for model II. We can refer to such behavior as the “structure forming” for model I and “structure breaking” in model II.
The composition dependence of viscosity, diffusion, and excess volume of binary mixtures is a problem of great importance in solution chemistry. Unfortunately, however, no model-based theoretical study of this problem has been carried out which shows the diverse behavior exhibited by these mixtures. Much of the theoretical studies have been directed to supercooled liquids or to phase separation kinetics. The work presented here is perhaps the first detailed study of the anomalous composition depedence of viscosity, diffusion, and excess volume. The novel aspect of this study is the introduction of two new models which are suffciently simple to allow detailed investigation. To capture the commonly held wisdom that the anomalies in composition dependence arise from the “structure-making” or “structure-breaking” among the two constituent species, the two models are made to represent the two extremes. However, we needed to carry out (NPT) MD simulations in order to allow for the fluctuation in the total volume of the system. In this paper, we have reported results of such constant temperature and pressure (NPT) MD simulations and searched for the correlations between excess volume and excess viscosity. The simulations show that these two are strongly but inversely
Nonideality in Binary Mixtures
J. Phys. Chem. B, Vol. 105, No. 39, 2001 9585
Figure 6. Diffusion coefficients (Dsi ) are plotted against solute compostion (x2) for model II. Filled circles indicate Ds1, while open circles denote Ds2. Error bars are given.
Future work shall consider variation of solute solvent size ratio and also pressure dependence of viscosity and excess volume of nonideal binary mixtures.
Figure 4. Radial distribution function for model I and II. All the partial radial distribution functions g11 (solvent-solvent), g12 (solute-solvent), and g22 (solute-solute) are plotted for 0.4 solute compostion. Panel a shows the results for model I, and panel b shows the same for model II. In both the cases, dashed lines show g11, short dashed lines show g22, and solid lines show g12.
Figure 5. Calculated diffusion coefficients (Dsi ) are plotted against solute compostion (x2) for model I. The filled circles indicate diffusion coefficients for solvent particles (Ds1), and open circles denote diffusion coefficients for solute particles (Ds2). Error bars are given.
correlated for both the attractive and repulsive model liquids, with the trends themselves opposite for the two models. We have also shown the re-entrant nature of excess viscosity with excess volume in the direction of changing solute composition. We have calculated diffusion coefficients (Ds1 and Ds2) for both the models which show nonmonotonic behavior against composition. We have also noticed a breakdown of Stokes-Einstein law in these systems where the product Diη against η gives a scattered distribution.
Acknowledgment. It is a pleasure to thank Dr. Srikanth Sastry and Mr. G. Sinivas for many useful discussions, about both NPT simulations and binary mixtures. We also thank Dr. Sarika Bhattacharyya and Ms. R. Vasanthi for discussions. This work was supported in parts by grants from the Department of Science and Technology (DST) and the Council of Scientific and Industrial Research, India. References and Notes (1) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1951, 19, 774. (2) Ben-Naim, A. J. Chem. Phys. 1977, 67, 4884. (3) Beddard, G. S.; Doust, T.; Hudales, J. Nature 1981, 294, 145. (4) Fleming, G. R. Chemical Application of Ultrafast Spectroscopy; The International Series of Monograph on Chemistry; Oxford University Press: New York, 1986. (5) Cohen, C.; Sutherland, J. W. H.; Deutch, J. M. Phys. Chem. Liq. 1971, 2, 213. (6) Knapp, E. W.; Fischer, S. F. J. Chem. Phys. 1982, 76, 4730. (7) Morillo, M.; Denk, C.; Sanchez-Burgos, F.; Sanchez, A. J. Chem. Phys. 1982, 76, 4730. (8) Berry, R. S.; Rice, S. A.; Ross, J. Physical Chemistry; Oxford University Press: New York, 2000. (9) Rowlinson, J. S.; Swinton, F. L. Liquids and Liquid Mixtures; Butterworth Scientific: Woburn, MA, 1982. (10) Pal, A.; Daas, G. J. Mol. Liq. 2000, 84, 327. (11) Nikam, P. S.; Shirsat, N.; Mehdi Hasan J. Chem. Eng. Ref. Data 1998, 43, 732. (12) Kob, W.; Andersen, H. C. Phys. ReV. E 1995, 51, 4626; Phys. ReV. Lett. 1994, 73, 1376. (13) Sastry, S. Phys. ReV. Lett. 2000, 85, 590. (14) Qunfang, L.; Yu-Chun, H. Fluid Phase Equlib. 1999, 154, 153. (15) Srinivas, G.; Mukherjee, A.; Bagchi, B. J. Chem. Phys. 2001, 114, 6220. Mukherjee, A.; Srinivas, G.; Bagchi, B. Phys. ReV. Lett. 2001, 86, 5926. (16) Zwanzig, R. Annu. ReV. Phys. Chem. 1965, 16, 67. (17) Zwanzig, R.; Mountain, R. D. J. Chem. Phys. 1965, 43, 4464. (18) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (19) Brown, D.; Clarke, J. H. R. Mol. Phys. 1984, 51, 1243. (20) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (21) Haile, J. M.; Graben, H. W. J. Chem. Phys. 1980, 73, 2412. (22) S. Stoddard, D.; Ford, J. Phys. ReV. A 1973, 8, 1504. (23) Heyes, D. M. J. Chem. Phys. 1988, 37, 5677. (24) Lyndel-Bell, R. M.; Rasaiah, J. C. J. Chem. Phys. 107, 1981 (1997).
