Langmuir 1999, 15, 437-446
437
Thermodynamics of Aggregation of Amphiphiles in Solution from Lattice Monte Carlo Simulations L. A. Rodrı´guez-Guadarrama, Sameer K. Talsania, Kishore K. Mohanty, and Raj Rajagopalan* Department of Chemical Engineering, University of Houston, Houston, Texas 77204-4792 Received June 5, 1998. In Final Form: October 19, 1998 Monte Carlo simulations of a three-dimensional lattice model of an amphiphile-solvent system are presented and examined for their usefulness in predicting the thermodynamics of self-association of amphiphiles in solution. The size distribution of the aggregates obtained from the simulations is usedsin combination with the mass action modelsto obtain the thermodynamic properties of the solution of aggregates and to obtain expressions for the coefficients of an empirical equation for the Gibbs energy available in the literature for aqueous nonionic surfactant solutions. The results are in good agreement with the empirical model, and the essential microstructural features and thermodynamic properties of the micellar solution can be extracted reasonably well using such an empirical fitting of the simulation data. The effects of the lengths of the solvophobic and solvophilic sections of the amphiphile on critical micelle concentration extracted from the simulations are consistent with experimental results reported for several nonionic surfactants in water.
1. Introduction Interest in the aggregation (self-assembly), phase behavior, and properties of solutions of amphiphilic molecules is well-established and is motivated by numerous applications in industry and science.1,2 The equilibrium structures of the amphiphilic molecules in solution are determined by a balance between the internal energy and the entropy of aggregation. The major force that governs the self-assembly derives from the solvophobic effect, which induces the amphiphilic molecules to associate with each other in order to minimize the number of solvophobic contacts.2 Theoretical efforts developed to understand the aggregation of amphiphiles in solution may be classified into two major categories, namely, phenomenological approximations and molecular simulations. In the phenomenological approximations,3-6 the aggregation of the amphiphiles is described by a multiple equilibrium model. The size distribution of the aggregates can be calculated from the difference between the standard chemical potentials of an amphiphilic molecule in an aggregate and of one in the free state. The variation in the standard chemical potential with the aggregate size is evaluated employing a phenomenological approximation, which is usually split into two major contributions, namely, the solvophobic effect described above and the repulsive forces arising from the nature of the headgroup as well as the distance of separation among the headgroups in the aggregate. These theories have been successfully applied in the study of micellization and phase behavior for aqueous ionic and nonionic surfactant solutions. However, * To whom correspondence should be addressed. Present address: Department of Chemical Engineering, University of Florida, Gainesville, FL 32611-6005. E-mail:
[email protected]. (1) Hiemenz, P. C.; Rajagopalan, R. Principles of Colloid and Surface Chemistry, 3rd ed.; Dekker: New York, 1997. (2) Israelachvili, J. N. Intermolecular and Surface Forces, 2nd ed.; Academic Press: New York, 1991. (3) Tanford, C. The Hydrophobic Effect; Wiley: New York, 1980. (4) Nagarajan, R.; Ruckenstein, E. J. Colloid Interface Sci., 1983, 91, 500. (5) Goldstein, R. E. J. Chem. Phys., 1986, 84, 3367. (6) Puvvada, S.; Blankschtein, D. J. Chem. Phys. 1990, 92, 3710.
a major problem with these treatments is that they contain parameters which cannot always be determined experimentally, thus requiring (often crude) approximations. Moreover, these treatments require a priori assumptions concerning the shape of the aggregates. It is convenient to classify the molecular simulation studies into three categories, namely, molecular dynamics, Brownian dynamics, and Monte Carlo simulation methods, on the basis of the molecular and temporal coarse-graining used. In molecular dynamics simulations, detailed atomistic potentials are used and the equations of motion are solved for each amphiphile. Even here, in view of the computational resources required, one is restricted to considering the amphiphilic molecules to be preassembled into an aggregate shape;7 even with such drastic simplifications, only a single aggregate or a few aggregates may be simulated for short periods of time. By simplifying the interactions, it is possible to simulate the spontaneous aggregation of amphiphiles,8 but this also requires considerable computer resources and is often beyond the scope of these techniques. Aggregation of short amphiphiles in solution has also been investigated recently using Brownian dynamics simulations.9 In this method, the amphiphilic molecules are represented by a coarse-grained bead-rod model and the solvent effects enter through a stochastic term in the equation of motion. Unshifted Lennard-Jones interactions are often employed between the beads, and differences in the nature of each bead are accommodated by changing the cutoff distance. Monte Carlo simulations of amphiphilic solutions in lattices have been studied extensively in the last 15 years. Larson10 has studied the structure of aggregates for different types of amphiphiles in solution. However, the majority of Larson’s work has focused on three-component amphiphile-oil-water systems. Self-assembly of diblock11 (7) Shelly, J. C.; Spik, M.; Klein, M. Langmuir, 1993, 9, 916. (8) Rector, D. R.; van Swol, F.; Henderson, J. R. Mol. Phys., 1994, 82, 1009. (9) von Gottberg, F. K.; Smith, K. A.; Hatton, T. A. J. Chem. Phys., 1997, 106, 9850. (10) Larson, R. G. J. Chem. Phys. 1989, 91, 2479.
10.1021/la9806597 CCC: $18.00 © 1999 American Chemical Society Published on Web 12/31/1998
438 Langmuir, Vol. 15, No. 2, 1999
and triblock12,13 copolymers in solution has also been studied by Mattice and associates. They have established scaling laws for the dependence of the cmc (the critical micelle concentration representing the onset of micellization) on the dimensionless interaction energy and the length of the solvophobic section of the diblock and triblock copolymers. Desplat and Care14 have studied a short amphiphilic chain and used the aggregate size distribution obtained from the simulations to study the thermodynamics of aggregation of amphiphiles in solution. A primary objective of this paper is to examine the adequacy of simple lattice Monte Carlo simulations (i.e., ones employing only the most essential interactions, e.g., the solvophobicity of the tail of the surfactant) for the calculation of the thermodynamics of aggregation of amphiphiles. In particular, we examine the correspondence between the results of the simulation and those of a typical phenomenological approximation reported in the literature5 for nonionic surfactants, so that the results of the simulations can be reduced to convenient analytical expressions for the thermodynamic and structural properties of the solution. Another objective is to study the effect of the lengths of the head and tail groups on the critical micelle concentration. The simulation results are compared with an empirical correlation reported in the literature15 for aqueous alkyl ethoxylate (a nonionic surfactant) solutions. The overall goal of the paper is to explore how simple, but sufficient, simulation models can be developed and to examine how the results can be reduced to convenient predictive correlations. The paper is organized as follows. First we discuss, in section 2, the basic methods employed in the study. This includes a brief description of the Monte Carlo simulation technique and the characterization of the aggregates. This is followed in section 3 by an extensive discussion of the transformation of the results (e.g., the cmc and the aggregate size distribution) to thermodynamic properties and of the reduction of the results to predictive correlations. A brief summary and some concluding remarks follow in section 4. 2. Background 2.1. Details of the Monte Carlo Simulations. The model uses a three-dimensional L × L × L cubic lattice with a coordination number of z ) 6; that is, only the nearest neighbor interactions are taken into account. To minimize any possible size effects, we have chosen a minimum box size of L ) 50. Excluded volume and periodic boundary conditions are used. A lattice site can be occupied by one of three types of beads: solvent, denoted by s; head, h; or tail, t. Amphiphilic chains, hitj, are composed of the appropriate numbers (i, j g 1) of h and t beads. Each empty site in the lattice represents a solvent molecule. The interaction energy mm′ between beads m and m′ is assumed to be nonzero and positive if and only if m ) t and m′ ) h or s. In this study, we have selected th/kBT ) ts/kBT ) /kBT, where kB is the Boltzmann constant and T is the absolute temperature. Therefore, the net energy associated with any configuration is a multiple of /kBT. This interaction energy is related to the Flory-Huggins (11) Wang, Y.; Mattice, W. L.; Napper, D. H. Langmuir 1993, 9, 66. (12) Wang, Y.; Mattice, W. L.; Napper, D. H. Macromolecules, 1992, 25, 4073. (13) Nguyen-Misra, M.; Mattice, W. L. Macromolecules 1995, 28, 1444. (14) Desplat, J. C.; Care, C. M. Mol. Phys. 1996, 87, 441. (15) Huibers, P. D. T.; Lobanov, V. S.; Katritzky, A. R.; Shah, D. O.; Karelson, M. Langmuir 1996, 12, 1462.
Rodrı´guez-Guadarrama et al.
parameter χ by χ ≡ z/kBT. The total internal energy for the model is given by
E ) (n + nth) kBT kBT ts
(1)
where nts and nht are the numbers of tail-solvent and headtail pairs, respectively. All simulations are started from a random initial configuration. The only type of move used for the modification of the configuration is the reptation move.16 In this move, the chain to be moved and the chain end to be the lead end for the reptation are chosen randomly. A move is chosen such that the self-avoidance and mutual avoidance conditions are not violated. The probability of acceptance of this move is calculated according to the standard Metropolis algorithm.17 The Metropolis probability is given by
{
(
P ) min 1, exp -
)}
∆E kBT
(2)
where ∆E is the difference in the total internal energy between the trial and old configurations. If a move is not accepted, the old configuration is maintained and a new move is tried. More than 107 moves are done in order to reach a region where the total energy is almost constant. These states are taken to be the equilibrium configurations over which the aggregates are characterized. To avoid metastable states, the simulations for each condition are started at athermal conditions and the interaction energy is increased in small steps (typically in steps of 0.1), after equilibrium is attained at each step, until the desired interaction energy is reached. We have also repeated the simulations several times for each condition in order to verify that the same results are obtained. 2.2. Basic Theory of Surfactant Aggregation. An amphiphilic chain is considered to belong to an aggregate if any of its tail beads occupies a neighboring site of a tail bead from another chain. Nonaggregated amphiphiles are called free or single amphiphiles. The size of an aggregate is described by the aggregation number n, which is the number of amphiphilic chains in the aggregate. The aggregation of amphiphilic molecules in dilute solutions can be represented by the multiple-equilibria model described by kn
An-1 + A1 y\z An
(3)
where A1 refers to the free amphiphilic molecule, An refers to an aggregate with aggregation number n, and kn is the association equilibrium constant for the reaction.4 The amphiphilic solution is composed of N0 solvent molecules, N1 single or free amphiphiles, and Nn aggregates of size n (2 e n e nmax, where nmax is the size of the largest aggregate). If the solution is assumed to be very dilute, then the mutual interactions between different species can be neglected. Aggregates of different sizes can be considered as distinct chemical species, each characterized by its own chemical potential. The standard chemical potentials are represented by µn° for an aggregate of size n, by µ0° for the solvent, and by µ1° for a single amphiphilic molecule. The standard state corresponds to an infinitely (16) Talsania, S. K.; Wang, Y.; Rajagopalan, R.; Mohanty, K. K. J. Colloid Interface Sci. 1997, 190, 92. (17) Metropolis, N. Rosenbluth,A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys., 1953, 21, 1087.
Aggregation of Amphiphiles in Solution
Langmuir, Vol. 15, No. 2, 1999 439
dilute solution. Then the Gibbs free energy of the solution is given by nmax
G)
∑ n)0
nmax
Nnµn° + kBT
∑ n)0
Nn ln
Nn Nt
(4)
where Nt is the total number of species on the lattice, that nmax Nn. The equilibrium state of the system is, Nt ) ∑n)0 corresponds to the minimum in the Gibbs free energy. This condition gives the following expression for the equilibrium size distribution of the aggregate of type n:
(
)
n∆Gn° ) exp ) Kn n kBT X1 Xn
(5)
where X1 and Xn are the mole fractions of singly dispersed amphiphiles and aggregates of size n, respectively, and ∆Gn° ) µn°/n - µ1° is the difference in the standard Gibbs free energy between an amphiphilic molecule in an aggregate of size n and a free amphiphile in the solvent. The overall association equilibrium constant Kn is related to the stepwise association equilibrium constants kn by nmax
Kn )
kn ∏ n)2
(6)
where
ln kn ) ln
Kn n∆Gn° - (n - 1)∆Gn-1° )Kn-1 kBT
(7)
Since the difference n∆Gn° - (n - 1)∆Gn-1° can be written as d(n∆Gn°)/dn, one has
ln kn ) -
(
)
d n∆Gn° dn kBT
(8)
The weight-average aggregation number (i.e., “molecular weight”) Nw is used as a measure of the average size of the aggregates and is defined as nmax
Nw )
n2Xn ∑ n)1 nmax
(9)
∑ nXn
n)1
3. Results and Discussion We have divided the results of the simulations into three major sections. In the first, properties such as the cmc, shape, and size of the aggregates are obtained for a base case amphiphile, h4t4, to see if the results are physically consistent. In addition, these results are used in subsequent sections in the analysis of the thermodynamic properties predicted by the model. In the second section, a method for determining thermodynamic properties such as Gibbs free energy, enthalpy, and entropy of aggregation from the simulations is presented. The values of the Gibbs free energy of aggregation obtained from the simulations are fitted to a phenomenological approximation reported in the literature for aqueous nonionic surfactant solutions5 in order to examine the validity of the empirical equation and the usefulness of combining the simulation results with empirical functions such as the one examined. The parameters obtained from such a fitting are compared
Table 1. Solubility Limits for Various Amphiphiles Predicted Using the Quasichemical Approximation for E/kBT ) 0.7 amphiphile
% by volume
h1t4 h2t4 h3t4 h4t4 h4t3 h4t5 h4t6 h4t7
0.16 0.74 2.50 7.99 100 1.44 0.33 0.08
a A solubility limit of 100% (in the case of h t ) indicates complete 4 3 solubility.
with those from the literature. In the third section, the effect of the lengths of the head and tail of the amphiphile on the cmc is studied. The results are consistent with a correlation reported in the literature for aqueous solutions of alkyl ethoxylate surfactants.15 The results and the comparisons are used to illustrate how simulations can be used as predictive tools for the thermodynamics of aggregation in amphiphile solutions. 3.1. Properties of the Aggregate. For the selfassembly of amphiphiles in solution, it is necessary to have a minimum degree of solvophobicity of the amphiphile.3 In the simulations, the degree of solvophobicity is determined by a combination of two factors: j, the length of the solvophobic tail of the amphiphile hitj, and /kBT, the dimensionless interaction energy. Hence, the selection of these two parameters is crucial in order to obtain aggregates. In the present paper, we have first chosen the amphiphile structure which we would like to study and then selected the appropriate energy parameter, as described below. The structures of the amphiphiles studied were selected so as to be consistent with typical nonionic surfactants (with the tail size equal to or larger than the size of the headgroup, i.e., j g i). To achieve a solvophobicity of the tail consistent with that of a typical hydrocarbon in water, we have used a minimum tail length of four beads and an interaction energy of /kBT ) 0.7. Using the quasichemical approximation,18 we have found that the solubility limit of the chain t4 is 0.005% by volume for /kBT ) 0.7; this is comparable to the solubility limit for n-octane in water at 25 °C. For the results reported in this section, we have used the amphiphile h4t4 with a dimensionless interaction energy of /kBT ) 0.7 (χ ) 4.2). A study of the effects of the lengths of the head and tail groups is presented later in section 3.3. For each amphiphile studied, there exists a minimum interaction energy below which the amphiphiles do not aggregate. For example, for h4t4, this minimum value is /kBT ) 0.56.19 Moreover, there is also a minimum energy above which the amphiphiles phase separate when the concentration exceeds the solubility limit. For h4t4, this value is roughly /kBT ) 0.67. Since the interaction energy chosen for examination here, that is, /kBT ) 0.7, is larger than 0.67, it is necessary to first determine the solubility limit for /kBT ) 0.7 and use amphiphile concentrations below the solubility limit to get true aggregates and not large clusters representing a separate phase. We have determined the solubility limit of h4t4 for /kBT ) 0.7 to be 8% by volume. The solubility limits for all of the amphiphiles used in the simulations have also been estimated using the quasichemical approximation and are reported in Table 1. (18) Guggenheim, E. A. Mixtures; Clarendon Press: Oxford, 1952.
440 Langmuir, Vol. 15, No. 2, 1999
Rodrı´guez-Guadarrama et al.
Figure 1. Calculation of the cmc for the amphiphile h4t4 (for the dimensionless interaction energy /kBT ) 0.7).
3.1.1. Critical Micelle Concentration. Once the solubility limit is established for a given amphiphile, the next step is to find the cmc, that is, the amphiphile concentration at which the formation of aggregates in the solution begins. We define the cmc as the amphiphile concentration at which the number of aggregated amphiphiles equals the nmax nXn. In number of free amphiphiles, that is, X1 ) ∑n)2 Figure 1, which shows the free amphiphile mole fraction X1 against the total h4t4 mole fraction Xamph, the above definition of the cmc corresponds to the intersection of the straight line defined by X1 ) Xamph/2 and the simulation data obtained for X1. At amphiphile concentrations well above the cmc, a decline in free amphiphile concentration is observed. This behavior has been explained by von Gottberg et al.9 and is due to the nonideality of the micellar solution at high concentrations (where interactions among the aggregates cannot be neglected; that is, the concentrations X1 and Xn in eq 5 need to be replaced with the corresponding activities). The cmc of h4t4 occurs at a mole fraction of 0.0025, equal to 2% by volume, for /kBT ) 0.7. The cmc’s for a set of amphiphiles with different architectures are presented in section 3.3 on the basis of the same definition of the cmc. 3.1.2. Shape of Aggregates. An approximate idea of the shape of the aggregates can be obtained from the interpretation of the principal moments of inertia.16 The three principal moments of inertia, I1, I2, and I3, are obtained from the eigenvalues of the matrix of the radii of gyration of the aggregates. The matrix of the radii of gyration is defined by
Rx2R,xR′ )
1
q
∑ (xRp - xRcm)(xR′p - xR′cm)
q p)1
(10)
where R and R′ can take values of 1, 2, or 3, corresponding to the three Cartesian coordinates, xRp represents the coordinate of a site p (in the aggregate) in the R direction, the subscript cm stands for the center of mass, and q is the total number of sites occupied by the aggregate. The characteristic length is defined as lR ) (IR)1/2. A measure of the aggregate asphericity can be obtained from the calculation of l1/l3 and l2/l3 as functions of the aggregate size. For a spherical aggregate l1/l3 ) l2/l3 ) 1, and for a cylindrical one, l1/l3 . 1 and l2/l3 ) 1. For the various amphiphiles and concentration ranges studied in the simulations, we find that the aggregates are roughly spherical given the geometrical constraints of the cubic lattice. The characteristic length ratios, shown
Figure 2. Characteristic lengths of the micelle as functions of the aggregation number for the amphiphile h4t4 at Xamph ) 0.004 and /kBT ) 0.7.
in Figure 2 for h4t4 for a concentration above its cmc (Xamph ) 0.004), are less than 2, indicating that the aggregates are not significantly aspherical. Similar results have been obtained for all the other amphiphiles studied here.20 The shapes of the aggregates have also been gleaned through visualization of the lattice at various configurations generated from the simulation. Figure 3 shows a snapshot of the whole simulation box for an equilibrium configuration and one of a randomly selected aggregate from the configuration, for h4t4 at a concentration above its cmc. This figure illustrates the roughly “spherical” shape of the aggregates. This result is used in section 3.2.2 to establish a simple model for the calculation of the entropy of aggregation. 3.1.3. Aggregate Size. If the average aggregate size is determined by only the reduction of the internal energy per amphiphile due to the formation of aggregates ( i.e., the solvophobic effect), one would expect a continuous increase in the aggregate size as the total amphiphile concentration in the solution is increased. Figure 4, which presents the weight-average aggregation number Nw as a function of the total mole fraction of amphiphile Xamph, shows that the aggregate size actually reaches a constant value. At low amphiphile concentrations, Nw increases as one increases the amphiphile concentration; however, beyond Xamph ) 0.014, the average aggregate size reaches a plateau of about 34 amphiphiles. Although the amphiphiles in the aggregate do not suffer any kind of repulsive interaction among their headgroups (other than the excluded-volume interaction), the aggregates do not grow indefinitely in view of a number of competing interactions, as discussed below. 3.2. Thermodynamics of Aggregation. In the first part of this section, we discuss a method for the evaluation of the Gibbs energy, enthalpy, and entropy of micellization directly from the simulation data. In the second part, the entropy of aggregation is estimated using a theoretical approximation and compared with the results from the simulations. In the third part, the Gibbs energy of aggregation calculated from the simulations is fitted to a phenomenological approximation proposed by Goldstein5 in order to examine Goldstein’s equation and to examine how reliably a combination of simulation results and (19) Rodrı´guez-Guadarrama, L.; Talsania, S. K.; Mohanty, K. K.; Rajagopalan, R. Unpublished work. (20) In calculating the prinicipal moments of intertia for determining the shapes of the micelles, we have included the head beads in eq 10 in the calculations of the radii of gyration. In general, the inclusion of the head beads in eq 10 is likely to influence the shape observed when the headgroup is long. However, similar calculations for h1t4 show that the shape of the core of the micelle is still roughly spherical.
Aggregation of Amphiphiles in Solution
Langmuir, Vol. 15, No. 2, 1999 441
Figure 3. Snapshot of the amphiphilic solution (h4t4) at Xamph ) 0.004 and /kBT ) 0.7. The gray circles are head beads, and the dark ones are tail beads.
442 Langmuir, Vol. 15, No. 2, 1999
Rodrı´guez-Guadarrama et al.
Figure 5. Entropy of aggregation as a function of the aggregation number for h4t4 at Xamph ) 0.004 and /kBT ) 0.7. Figure 4. Weight-average aggregation number as a function of the surfactant concentration for h4t4 for a dimensionless interaction energy of /kBT ) 0.7.
equations such as Goldstein’s can be used to employ simulations as predictive tools. The last part of this section concludes with a derivation for the cmc in terms of the empirical parameters determined from the Goldstein approximation in the previous part. 3.2.1. Estimation of ∆G, ∆H, and ∆S from the Simulations. For a given amphiphile concentration (above the cmc) and dimensionless interaction energy, the size distribution of the aggregates obtained from the simulations can be used directly to calculate the Gibbs free energy of aggregation by rearranging eq 5 as follows:
∆Gn° 1 ) ln X1 - ln Xn kBT n
(11)
(12)
Then the entropy of aggregation ∆Sn°/kB follows from
∆Sn° ∆Hn° ∆Gn° ) kB kBT kBT
4 3 πR ) jna3 3
(13)
Using the size distribution data obtained from the simulations, we have calculated the entropy of aggregation by the method described above. Figure 5 shows ∆Sn°/kB as a function of the aggregate size for h4t4 at a concentration of Xamph ) 0.004 and reveals that as the size of the aggregates increases, the entropy of aggregation decreases. The relation between this reduction in the entropy and the loss of conformational freedom experienced by the amphiphile chain upon transfer from the solution into an aggregate is examined in the next section using a simple theoretical model. 3.2.2. Comparison with Predictions Based on Conformational Statistics. The entropy of aggregation can be calculated from the number of conformations available to an amphiphile inside an aggregate and in the solvent, as has been suggested by Goldstein.5 If one can assume that the aggregates are spherical with a solvophobic core of radius R containing the tail segments of the amphiphiless an assumption roughly valid for n g 10, as confirmed by
(14)
where a is the characteristic length of a lattice site and j is the number of tail beads in the amphiphile. To estimate the entropy of formation of an aggregate of size R, we first assume that the end-to-end distance r of the amphiphilic molecules in the core can be described by a random walk with the Gaussian probability distribution function P(r),
P(r) )
If ∆Gn°/kBT is calculated for a range of interaction energies, the standard enthalpy of aggregation ∆Hn°/kBT can be determined at any */kBT within the range from the GibbsHelmholtz equation
d(∆Gn°/kBT) ∆Hn° |*/kBT ) | kBT d ln(/kBT) */kBT
simulationssthen the volume of the core can be related to the aggregation number as
( ) ( ) 3 2πr02
3/2
exp -
3r2 2r02
(15)
where r0 ) j1/2a is the average extension of an amphiphilic tail if it is assumed to be an ideal chain. The entropy for the configuration with r ) R can then be calculated from the probablity P(R) for r ) R, that is, Sn°(R) ) kB ln P(R). The change in the entropy of n chains due to aggregation is then given by
( )
∆Sn° Sn° - S1° P(R) 3 ) ) ln ) kB kB 2 P(r0)
[( ) ( )] n πj1/2
2
243 128
1/3
(16) One observes from Figure 5 that the magnitudes of the standard entropy of aggregation obtained from the simulations (eqs 11-13) and from eq 16 are of the same order. Note that ∆S1°/kB (i.e., for n ) 1) obtained from this theory is not equal to zero, as the assumption of spherical aggregates (and hence eqs 14 and 16) is not valid for n < 10. The differences observed between the results of the simulation and those of eq 16 can be attributed to a number of approximations that have been used in arriving at eq 16, including the use of short amphiphilic chains, which are not correctly described using the simple ideal chain model, and the assumption that the core of the aggregates is perfectly spherical. However, the general agreement between the simulations and eq 16 indeed suggests that the reduction in the entropy of aggregation is due to the decrease in the number of possible conformations of the amphiphile when it is transferred from the solvent environment to an aggregate. 3.2.3. Phenomenological Model of Goldstein for ∆G. It is possible to relate the dependence of the standard Gibbs free energy of aggregation ∆Gn°/kBT on the aggregation number n to intermolecular interactions.3 A plausible form
Aggregation of Amphiphiles in Solution
Langmuir, Vol. 15, No. 2, 1999 443
results and ours in the dependence of β, ξ, and η on i and j, and it is instructive to discuss the physical significance of the above parameters and their dependence on the lengths of the head and tail groups. The term β represents the main driving force for aggregation and consists of two contributions, namely, (i) the favorability of transferring the tail segments from the solvent to the core of the aggregate and (ii) the free energy change due to the proximity of the head segments to each other in the aggregate. This term is given by
β ) Cβ(i + j) + constant
(18)
where the value of Cβ (for which Goldstein used 0.9) is typically of the order of unity. Our simulations, on the other hand, lead to the relation Figure 6. Gibbs free energy change as a function of the aggregation number for the amphiphile h4t4 at different amphiphile concentrations for /kBT ) 0.7. Table 2. Parameters for the Goldstein Model from Simulations amphiphile
η
β
ξ
h1t4 h2t4 h3t4 h4t4 h4t3 h4t5 h4t6 h4t7
0.056 0.091 0.115 0.126 0.157 0.117 0.196 0.286
10.8 10.6 10.7 10.7 9.35 12.7 14.9 17.4
10.7 10.5 10.6 10.5 9.19 12.6 14.8 17.3
Table 3. Parameters η, β, and ξ of the Goldstein Model for the Amphiphiles hitj As Reported by Goldstein and from the Simulations Goldstein5 Simulations
η
β
ξ
1.0j-1/3 0.0043(i + j)1.67
(i + j) + 1.5 2.07j + 2.55
3.0j2/3 3.7j0.77
β ) Cβj + constant
(19)
showing that β depends exclusively on the tail length. This deviation from Goldstein’s expression is a consequence of the fact that the interaction energy between neighboring head beads is taken to be zero in our work. The term containing the parameter ξ in eq 17 represents the residual solvophilic-solvophobic interactions occurring at the periphery of the aggregate core and is characterized by a surface free energy per unit area, τ, which is proportional to the surface area of the core 4πR2 (for a spherical aggregate). Goldstein’s result gives the relation
ξ ) Cξ j2/3
(20)
where the value of Cξ (taken to be 3.0 by Goldstein) is of the same order as Cβ in eq 19. The simulations lead to
ξ ) Cξ j0.775
(21)
a
The correlation coefficients for the fits obtained from the simulations are 0.89, 0.99, and 0.99, for η, β, and ξ, respectively.
for the function ∆Gn°/kBT proposed by Goldstein5 for aqueous solutions of nonionic surfactants is given by
∆Gn° ) -β + ξn-1/3 + ηn2/3 kBT
(17)
where the coefficients β, ξ, and η are all positive. The physical significance of these parameters is discussed below. The above phenomenological expression is meant for symmetric amphiphiles (i.e., half are solvophobic (tail) segments, and the other half are solvophilic). As noted by Goldstein, the case of asymmetric amphiphiles can be studied using a simple modification to the equation. In what follows, we discuss such an extension and reformulate the Goldstein results for the case of our simulations. Data for ∆Gn°/kBT obtained from our simulations can be fitted to eq 17 to determine β, ξ, and η for a number of amphiphiles hitj (see Figure 6 for the results for a constant interaction energy of /kBT ) 0.7), and one can then obtain β, ξ, and η as functions of i and j. Table 2 presents the parameters for a certain range of i’s and j’s, and Table 3 presents their dependence on i and j thus obtained from the simulations. Also in Table 3 are the expressions presented by Goldstein.5 While ∆Gn°/kBT obtained in the simulations can be fitted well with eq 17, Table 3 reveals some differences between Goldstein’s
The difference in the exponent between eqs 20 and 21 can be attributed to the fact that the aggregates do not have a perfect spherical shape. Finally, the term containing the parameter η in eq 17 represents the entropic contribution to the Gibbs free energy due to the reduction in the configurational entropy of an amphiphile chain when it is transferred from the solvent environment to an aggregate. The coefficient η is described by the relation
η ) Cηj-1/3
(22)
where Cη is a purely geometric factor of order unity. Its exact value (see eq 16) is given by
Cη )
(
)
243 128π2
1/3
) 0.5772 ≈ 0.6
(23)
for which Goldstein used a value of 1.0, as reported in Table 3. The simulations, however, lead to
η ) Cη(i + j)1.67
(24)
We believe that the difference between the Goldstein expression and ours can be traced to the insufficiency of the following assumptions in Goldstein’s estimation of the change in the configurational entropy. (1) The aggregate is spherical. The approximation implicit in this assumption is evident from the simulations (see, for example, Figure 3b).
444 Langmuir, Vol. 15, No. 2, 1999
Rodrı´guez-Guadarrama et al. Table 4. Xcmc’s for Various Amphiphiles Based on Simulations and Eq 30
Figure 7. Aggregate size distribution for the amphiphile h4t4 at different concentrations for /kBT ) 0.7.
(2) The probability distribution function for the endto-end distance of the tail segment is Gaussian, an assumption that is expected to break down for the small tail segments considered in our work (i.e., 3 e j e 7). (3) The average extension of the tail segment in the solvent is described by the ideal chain model (r0 ) j1/2a). This is, perhaps, the most severe of the three assumptions. One can improve on this result by calculating the average size of the tail segment in the solvent from the result based on the excluded-volume effect (i.e., r0 ) j3/5a). Another option is to calculate the exponent ν in the relation r0 ) jνa from simulations of a dilute polymer solution with the polymer-solvent interaction energy (ps ) 0.7) as the only nonzero interaction energy. A detailed calculation of the above parameters is, however, not needed for the purpose of the present paper. Instead, we present merely the expressions obtained from the results of the simulations for a number of i’s and j’s (1 e i e 4 and 3 e j e 7). These expressions will be used in subsequent sections for predicting the cmc. It is also useful to examine the extent of information loss resulting from the reduction of the data from the simulations to eq 17 and what features of the microstructure are preserved in the reduction of the data. Therefore, we have used the values obtained for β, ξ, and η to regenerate the aggregate size distribution. This can be done by numerically solving the total amphiphile concentration nmax
Xamph )
∑ nXn
(25)
n)1
jointly with
Xn ) X1n exp{-n(-β + ξn-1/3 + ηn2/3)}
(26)
(obtained by combining eqs 11 and 17) for X1 and then substituting the result into eq 5 to get the values for nXn. The aggregate size distributions, obtained from the simulations and calculated using the method described above, are shown for two different amphiphile (h4t4) concentrationssabove the cmcsin Figure 7. The differences between the two distributions indicate that the general shapes of the size distribution and the maximum aggregate sizes are preserved although some of the details are lost in the transformation. 3.2.4. cmc: Empirical Equation from Simulations. The cmc’s for various amphiphiles hitj (based on /kBT ) 0.7), as determined from the simulations, are reported in Table
amphiphile
Xcmc simulations
Xcmc (eq 30 + Table 2)
Xcmc (eq 30 + Table 3)
L
h1t4 h2t4 h3t4 h4t4 h4t3 h4t5 h4t6 h4t7
0.0012 0.0018 0.0021 0.0026 0.0060 0.0007 0.0003 0.0001
0.0014 0.0028 0.0038 0.0044 0.0160 0.0010 0.0006 0.0002
0.0015 0.0023 0.0033 0.0047 0.0143 0.0015 0.0005 0.0002
50 50 50 50 50 60 80 100
4. The results, which show the expected increase in the cmc with head size and the decrease with tail size, are in qualitative agreement with experimental observations21-25 reported for aqueous solutions of ionic and nonionic surfactants. A quantitative description of the variation of the cmc’s is proposed in section 3.3. Here, we obtain a relation for Xcmc in terms of the parameters η, β, and ξ of the Goldstein equation (i.e., eq 17). A relation for Xcmc can be developed as follows. First, note that a peak in the size distribution of the aggregates develops for Xamph beyond the cmc, as can be seen clearly in Figure 7. The location of this peak, at n ) n* (say), identifies the most probable aggregate size. As has already been shown by Talsania et al.,16 the first appearance of this peak occurs at the cmc (the size distribution being monotonically decreasing with n prior to the cmc). We denote the magnitude of n* at the cmc by ncmc. Since, as also shown by Talsania et al., the size distribution is fairly flat around ncmc, one can take Xncmc ≈ Xncmc-1. Equation 5 therefore implies that
ln(X1)cmc ≈ ≈
ncmc∆Gn° - (ncmc - 1)∆Gncmc-1° kBT ∆Gn° kBT
(27)
where (X1)cmc is the mole fraction of free amphiphile at the cmc. Note that we have also assumed that ∆Gncmc° ≈ ∆Gncmc-1° in arriving at eq 27. As noted by Nagarajan and Ruckenstein,4 ncmc is usually taken as the aggregation number at which ∆Gn°/kBT reaches its minimum, that is,
( )|
d ∆Gn° dn kBT
)0
(28)
n)ncmc
Thus, one can obtain an estimate for ncmc in terms of η, β, and ξ from eq 17
ncmc )
1ξ 2η
(29)
(21) Ambrosono, L.; Costatino, L.; D’Errico, G.; Vitagliano, V. J. Colloid Interface Sci. 1997, 190, 286. (22) Mukerjee, P.; Mysels, K. J. Critical Micelle Concentration; National Standard Data Series, National Bureau of Standards: Washington, DC, 1971. (23) Ameri, M.; Attwood, D.; Collet, J. H.; Booth, C. J. Chem. Soc., Faraday Trans. 1997, 93, 2545. (24) Jiding, X.; Zhengyu, H. In Surfactants in Solution; Mittal, K. L., Bothorel, P., Eds.; Plenum: New York, 1986; Vol. 5, p 1055. (25) Hofland, H. E. J.; Bouwtra, J. A.; Gooris, G. S.; Spies, F.; Talsma, H.; Junginger, H. E. J. Colloid Interface Sci. 1993, 161, 366.
Aggregation of Amphiphiles in Solution
Langmuir, Vol. 15, No. 2, 1999 445
which, in combination with eqs 27 and 17, leads to
ln
Xcmc 3 ) 1/3(ξ2η)1/3 - β 2 4
Table 5. Coefficients c1, c2, and c3 in the Correlation for log Xcmc Given by Eq 33
(30)
where Xcmc is twice (X1)cmc, the mole fraction of free amphiphile at the cmc, consistent with our previous definition of the cmc. The results obtained for the cmc directly from the simulations and the ones based on eq 30 (with the values of η, β, and ξ determined as in section 3.2.3) are shown in Table 4. Despite the approximations involved in the above derivation and in obtaining η, β, and ξ by nonlinear curve fitting, the comparison in Table 4 shows that eq 30 can serve as a reasonable predictive tool. Finally, notice thatsin view of eqs 7 and 27sthe stepwise association constant kn at n ) ncmc is given by
kncmc )
1 (X1)cmc
(31)
a relation previously identified by Nagarajan and Ruckenstein. An expression for kncmc in terms of β, ξ, and η follows readily from eq 30. Expressions analogous to eqs 30 and 31 for Xcmc and kncmc, respectively, have been presented by Nagarajan and Ruckenstein4 using an empirical relation similar to Goldstein’s suggested by Tanford:3,26
∆Gn° ) -β + ξn-1/3 + η′n1/3 kBT
(32)
where the last term (with the coefficient η′) now represents headgroup repulsion, while the other two have the same physical significance as in Goldstein’s equation, that is, eq 17. We have not examined the use of the Tanford equation here, since the simulation model as employed in this paper does not include any interaction between neighboring head beads. 3.3. cmc: Effect of Head and Tail Lengths. Accurate experimental measurements of the cmc for several linear alkyl ethoxylate surfactants have been reported in the literature.21-25,27 Empirical relationships have also been proposed relating the structural features of the homologous surfactant molecules to the cmc. Here, we shall consider the nonionic surfactants of the alkyl ethoxylate family, which can be represented by CjEi, where i represents the number of ethylene oxide (E ) CH2CH2O) units and j the number of alkyl (C ) CH2) units. Huibers et al.15 have established, on the basis of regression, a linear relationship between the experimental values of the cmc and the number of alkyl and ethylene oxide units. This relation is given by
log Xcmc ) c1 + c2i + c3j
(33)
where c1, c2, and c3 are empirical regression coefficients (and log stands for logarithm to the base 10). We have used the same form of equation to correlate the cmc data obtained from the simulations for hitj. Table 5 reports the values for the coefficients obtained experimentally for alkyl ethoxylates in water at a temperature of 25 °C and for the simulation data. To get a physical interpretation for the coefficients c2 and c3, it is convenient to write the cmc in terms of the (26) Tanford, C. J. Phys. Chem. 1974, 78, 2469. (27) Rosen, M. J. Surfactants and Interfacial Phenomena; Wiley: New York, 1978.
c1
c2
c3
CjEi (experimental15) -0.099 0.044 -0.496 hitj (simulations) -1.14 0.080 -0.454
correlation coeff 0.994 0.993
Gibbs energy change, as derived in section 3.2.4. Thus, from eq 27 one has
ln (X1)cmc ) ln
Xcmc ∆Gn°cmc ∆Gagg° ) ) 2 kBT kBT
(34)
where we have denoted ∆Gncmc°/kBT by ∆Gagg°/kBT, for convenience. For an alkyl ethoxylate surfactant in water, we separate ∆Gagg°/kBT, for the Gibbs energy change for the transfer of a surfactant from water to aggregates, into three main contributions, namely, one for the transfer of the terminal methyl group (CH3), another for intermediate methyl groups in the chain (CH2), and one for the headgroup. Assuming that the contribution from the terminal methyl group is equal to that from a chain methyl group and that the headgroup contribution is not affected by the length of the surfactant tail, one obtains
∆GCH2° Xcmc ∆Gagg° )j + constant ) log 2.303kBT 2.303kBT 2
(35)
The term ∆GCH2°/2.303kBT (which represents the coefficient c3 in eq 33) is a constant that reflects the free energy change involved in transferring the tail group of the surfactant from an aqueous environment to the aggregate. Negative values of c3 favor aggregation and account for the fact that the cmc decreases with an increase in the length of the tail group. A similar argument can be used to explain the relatively small dependence of the cmc on the length of the head of the amphiphile. Thus c2, which is represented by ∆Ghead°/2.303kBT, reflects the free energy change involved in transferring the headgroup from a solvent environment to the aggregate. This free energy change is positive and therefore opposes aggregation. The analysis of the results presented in Table 5 using the above arguments suggests that the free energy changes involved in transferring the tail of the alkyl ethoxylate surfactant from an aqueous environment to the aggregates are very similar to those obtained from the simulations presented here. Thus, the solvophobicity used in the description of the amphiphiles in the present lattice model is very similar to that exhibited by alkyl ethoxylate surfactants in water. The analysis in this section also shows how one can approach the establishment of the correspondence between the model parameters and real systems. 4. Concluding Remarks We have investigated the usefulness of a simple lattice model in describing the thermodynamics of aggregation of short amphiphiles and have tested a phenomenological model proposed by Goldstein for aqueous nonionic surfactant solutions using the results of the simulation. Goldstein’s model takes into account the existence of three principal forces involved in the process of aggregation, namely, the solvophobic effect, the reduction in configurational entropy, and, finally, the creation of a solvophobic interface. Using the simulation results and Goldstein’s model, we have outlined a way to formulate generalized correlations for the principal driving forces for micellization as functions of the lengths of the headgroup and the tail of the amphiphiles, so that one can predict the properties of the aggregates for any amphiphile for a given
446 Langmuir, Vol. 15, No. 2, 1999
dimensionless interaction energy. We establish this point by generating the size distribution for a given amphiphile using the correlation and by comparing the results with those obtained directly from the simulations. We have also used the information obtained from the simulations for the cmc of different amphiphiles to determine the approximate Gibbs free energy changes involved in transferring the head and tail of the amphiphile to the aggregates. These terms were compared with the experimental results reported in the literature for alkyl ethoxylate surfactants, and it is seen that the free energy changes are very similar to these obtained from the simulations. Although our model does not take into account the complicated details of aqueous nonionic surfactant solutionsssuch as hydrogen bonding, hydration of the headgroup, and detailed molecular effectssthe results indicate that it can still be used as a simple model system for studying and testing theories of amphiphilic solutions. It should be noted that it is not a simple matter to map the architecture of a real amphiphilic molecule or a polymer chain or the energies of interactions onto lattice models. Lattice models are most useful to establish some of the characteristic universal properties and behavior of
Rodrı´guez-Guadarrama et al.
macroscopic systems, that is, properties and behavior that depend on global features such as the temperature, solvent quality, and the like.28 The work presented here demonstrates that simple lattice models exhibit many features of self-associating systems in agreement with theory and experiment. In a forthcoming publication, we shall extend the proposed model to study solutions of mixed amphiphiles. Acknowledgment. This work was partially funded by the Texas Higher Education Coordinating Board and the Gulf Coast Hazardous Substance Research Center. We thank Anand Jaganathan, Jorge Jimenez, and Claudia Marin of the University of Florida, Dang Nhan of the University of Houston, and Daniel Schmidt of the Rheinisch-Westfaelische Technische Hochschule Aachen for reading the manuscript critically and for their constructive criticisms of the contents and the presentation. We also extend special thanks to Nikhil Joshi for his assistance with a number of simulations reported here. LA9806597 (28) Vanderzande, C. Lattice Models of Polymers, Cambridge University Press, Cambridge, U.K., 1998.
