3748
Langmuir 1995, 21, 3748-3756
Modeling of Nonionic Micelles Christopher M. Wijmans* and Per Linse Physical Chemistry I, Chemical Center, Lund University, P.O.Box 124, S-221 00 Lund, Sweden Received May 17, 1995. I n Final Form: July 5, 1995@ The equilibrium among differently sized micelles of nonionic linear surfactant chains, as well as their shape and density profiles, have been studied using Monte Carlo (MC) simulations and self-consistentfield (SCF)calculations. The SCF theory predicts a far lower critical micellizationconcentration (cmc),a larger aggregation number at the cmc, and a sharper interface between the lyophobic and lyophilic segments in the micelle, as compared to the MC results. The mean-field approximation used in the SCF theory, and in particular its effect on the chemical potential of the bulk state, is the main factor contributing to the discrepancies in the cmc. The difference in aggregation number and sharpness of the lyophobidyophilic interface are caused by the neglect of micellar shape fluctuations in the SCF theory, Such fluctuations are shown to be more important for surfactants with equal head and tail sizes than for surfactants with a larger head than tail group. The equilibrium micelle size distribution is used to study how the standard chemical potential of surfactants in aggregates of different sizes depends on the aggregation number.
1. Introduction self-consistent field (SCF) theory in order to study the self-assembly of surfactant molecules into micelles, using Nonionic surfactants form micelles in aqueous solution, a lattice formalism to generate all possible conformations and various approaches exist to predict micelle properties of these molecules. Thus, they were able to compute the from the molecular composition of the surfactant. Tanvolume fraction profiles of the head and tail groups in the ford’,2pioneered the “opposingforces”approach, a method micelle which correspond to the minimum in the free whereby the free energy of aggregation is decomposed energy of the system. Bohmer et al.14J5showed that this into an attractive and a repulsive component. The theory closely follows many experimental observations. attractive component results from the system’s urge to In a series of papers Linse16-19demonstrated how the SCF minimize water-hydrocarbon contact, whereas the retheory can be used to model solutions of (impure, polypulsive component accounts for the repulsion between disperse, and commercially available) poly(ethy1ene oxhead groups in a micelle. Using expressions for these forces, the formation of micelles has been ~ t u d i e d , ~ - ~ ide)-poly(propy1ene oxide) copolymers. Useful as it is, the SCF theory still contains four major making it possible to calculate quantities such as the assumptions. First, a lattice formalism is used to describe critical micellization concentration (cmc) and the agthe system. Space is divided into concentric, equidistant gregation number. However, this approach gives no shells, each of which is further divided into lattice cells. information about the molecular structure of a micelle, as Second, a Markov approximation is used to enumerate the calculation is based on some assumption about this the conformations of the molecules. This means that the structure. This approach has been extended to deal with chains are modeled as self-intersecting walks. Each micellar growth as well as phase ~ e p a r a t i o n . ~ conformation is weighted by a potential of mean force Molecular modeling methods have also been developed generated by all other conformations. This potential field which explicitly take many different conformations of the is a mean-field, which constitutes the third important molecules in a micelle into account.6-10 The problem is assumption. The potential of mean force (and, therefore, reduced to that of a single molecule which interacts with also the volume fractions) vary only in the radial direction. a field generated by all other components. A drawback No fluctuations are allowed for that take place within a of these calculations is that they are all based on some layer. Fourth, fluctuations in size and shape ofthe micelle important restriction concerning the conformations of the are neglected. However, it is also possible to apply the molecules. For example, the head groups are pinned to SCF method in a cylindrical or planar (lamellar)geometry. a sphere or given another predetermined distribution. By comparing the critical aggregation concentrations in In the approach pioneered by Leermakers and Scheutthe different geometries one can then, for example, jensl1-l3 no such restrictions are made. They applied a determine under what conditions the formation of spherical micelles is favored above, for example, (lamellar) Abstract published i n Advance A C S Abstracts, September 15, membrane f ~ r m a t i o n . ’ ~But J ~ even if it has thus been 1995. established that a certain surfactant prefers to form (l)Tanford, C. J. Phys. Chem. 1974,24,2469. (2) Tanford, C. The Hydrophobic Effect. Formation of micelles and spherical aggregates, that does not automatically mean biological membranes; 2nd ed.; Wiley: New York, 1980. that more complicated shapes may be totally neglected. (3)Nagarajan, R.; Ruckenstein, E. J. ColloidlnterfaceSei. 1979,71, These drawbacks of the SCF theory can all, in principle, 580. ( 4 )Nagarajan, R. Adu. Colloid Interface Sci. 1986,26,205. be circumvented by performing simulations. To follow @
(5) P u w a d a , S.; Blankschtein, D. J. Chem. Phys. 1990,92,3710. (6) Marcelja, S. Biochim. Biophys. Acta 1974,367,165. (7) Dill, K. A. J . Phys. Chem. 1982,86, 1498. ( 8 ) Gruen, D. W. R. J. Phys. Chem. 1985,89,146. (9) Ben-Shad, A.; Szleifer, I.; Gelbart, W. M. J. Chem. Phys. 1985, 83, 3597. (10) Ben-Shad, A.; Szleifer, I.; Gelbart, W. M. J. Chem. Phys. 1986, 85, 5345. (11)Leermakers, F.A. M.; Scheutjens, J. M. H. M. J. Phys. Chem. 1989,93, 7417. (12) Leermakers, F. A. M.; Scheutjens, J. M. H. M. J . Colloidlnterface Sci. 1990,136, 231.
(13) Leermakers, F. A. M.; van der Schoot, P. P. A. M.; Scheutjens, J. M. H. M.; Lyklema, J. In Surfactants in Solution; K. L. Mittal, Ed.; Plenum: New York, 1990; Vol. 7; pp 25. (14)Bohmer, M. R.; Koopal, L. K. Langmuir 1990,6, 1478. (15) Bohmer, M. R.; Koopal, L. K.; Lyklema, J. J . Phys. Chem. 1991, 89,9569. (16) Linse, P. J . Phys. Chem. 1993,97, 13896. (17) Linse, P. Macromolecules 1993,26,4437. (18) Linse, P. Macromolecules 1994,27,2685. (19) Linse, P. Macromolecules 1994,27,6404.
0743-7463/95/2411-3748$09.00/0 0 1995 American Chemical Society
Modeling of Nonionic Micelles the self-assembly of surfactants using a realistic atomic model is a computationally very demanding task. Molecular dynamics (MD) studies, such as those reported in refs 20-26, therefore tend to consider micelles which are constructed “a priori”. Smit et al.27~28 introduced a simplified model for short surfactant molecules with which they could observe spontaneous formation of micelles. However, it is debatable whether a real thermodynamic equilibrium can be reached within a realistic simulation time.29 As a n alternative several authors have published results of Monte Carlo (MC) simulations of micelle formation using even simpler m ~ d e l s . ~ OIn- ~this ~ context we will especially direct our attention to the model of Mattice and c o - w o r k e r ~ ,who ~ ~ -studied ~~ surfactant selfassembly by treating the surfactant molecules as lattice chains. Among other things, they calculated segmental distribution profiles of the micelles and the dependence of the cmc on the surfactant composition. Their approach has many similarities with that of the SCF model. The aim of this paper is to make a comparison of the results of these two approaches and to understand the origin of the differences obtained. The remaining part of this paper is organized as follows. In the next section we present the three-dimensional lattice model used in the MC simulations, followed by a brief description of the SCF theory and a discussion of the differences that might be expected between the results of these two approaches. In the Results section, data for two different model systems are given. Most attention is paid to a 5 vol % solution ofsurfactant molecules consisting of ten lyophobic and ten lyophilic segments. In order to look a t the influence of the molecular composition of the surfactant, some results are also given for molecules which have only five lyophobic segments. In this case the lyophobicity of the individual segments is increased by a factor 2. The Discussion section consists of two parts. First, we examine the differences between the results of both approaches, given the different assumptions underlying each of them. In the second part, we consider in more detail the micellar size distributions obtained from the three-dimensional lattice model. Conclusions are drawn regarding thermodynamic parameters which govern these distributions. The paper ends with a n Appendix containing more technical aspects of the SCF theory. (20) Jonsson, B.; Eldholm, 0.;Teleman, 0.J. Chem. Phys. 1986,85, 2259. (21) Watanabe, K.; Ferrario, M.; Klein, M. J. Phys. Chem. 1988,92, 819. .~ (22) Wendoloski, J. J.; Kimatian, S.J.; Schutt, C. E.; Salemme, F. R. Science 1989, 243, 636. (23) Shellev. J.:Watanabe. K.; Klein, M. L. Int. J.Quantum Chem.. Quantum SioL Symp. 1990, 17, 103. (24) Shelley, J. C.; Sprik, M.; Klein, M. L. Langmuir 1993, 9, 916. (25) Laaksonen, L.; Rosenholm, J. B. Chem. Phys. Lett. 1993,216, 429. (26) Bocker, J.; Brickmann, J.; Bopp, P. J. J. Phys. Chem. 1994,98, 712. (27) Smit, B.; Hilbers, P. A. J.; Esselink, K.; Rupert, L. A. M.; van Os. N.M. Nature 1990.348. 624. i28) Smit, B.; Hilbers, P. A. J.; Esselink, K.; Rupert, L. A. M.; van Os, N. M. J. Phys. Chem. 1991, 95, 6361. (29) Leermakers, F. A. M.; Lyklema, J. Colloids Surfaces 1992, 67, 239. -.. (30) Haan, S. W.; Pratt, L. R. Chem. Phys. Lett. 1981, 79, 436. (31) Owenson, B.; Pratt, L. R. J . Phys. Chem. 1984,88, 2905. (32) Pratt, L. R.; Owenson, B.; Sun, S. Adv. Colloid Znterface Sci. 1986, 26, 69. (33) Larson, R. G. J . Chem. Phys. 1992, 96, 7904. (34) Rodrigues, K.; Mattice, W. L. J. Chem. Phys. 1991, 95, 5341. (35) Wang, Y.; Mattice, W. L.; Napper, D. H. Langmuir 1993,9,66. (36) Haliloglu, T.; Mattice, W. L. Chem. Eng. Res. 1994, 49, 2851. (37) Adriani, P.; Wang, Y.; Mattice, W. L. J. Chem. Phys. 1994,100, 7718. (38) Zhan, Y.; Mattice, W. L. Macromolecules 1994, 27, 677.
Langmuir, Vol. 11, No. 10, 1995 3749 2. Theoretical Modeling Lattice Model and Monte Carlo Simulation Technique. A three-dimensional simple cubic lattice with dimensions L x L x L is used, and periodic boundary conditions are applied in all three directions. Chain like surfactant molecules with composition A N ~ B are N ~ introduced into the lattice, where N Ais the number of segments in the A block and N B the number of segments in the B block. Each segment occupies one lattice site, and empty lattice sites are considered to hold a solvent molecule S. The MC simulation procedure is as follows. A segment is chosen a t random. If the segment is a n end segment, then with equal probability a flip movement of this segment or a reptation movement of the whole chain is attempted. If the segment is not situated a t the end of the chain, depending on the local chain conformation either a Verdier-Stockmayer “flip” movement39 or a “kink jump”40is attempted. The new configuration is accepted according to the Metropolis algorithm4Iwith a probability minimum of [l, exp(-pAU)l, where AU is the energy difference associated with the change in configuration. The sum of the energies of all nearest-neighbor contacts constitutes the energy of a configuration. We define a contact energy, cpq,withp and q = A, B, or S. We always choose cpp = 0, so that xpq = 6cpq. Both in the MC simulations and the SCF calculations we have used XAB = XAS > 0 and XBS = 0. We define one iteration step as the procedure wherein on average one MC step is attempted for each segment in each chain. The simulation is started by placing all chains parallel to each other on the lattice. First, typically 0.05 x lo6 iterations are performed with all contact energies equal to zero in order to get a homogeneous and isotropic distribution of the surfactants. Then up to 7 x lo6 iterations (system 2) with the interactions turned on are performed to obtain an equilibrated system. Subsequently, the actual simulation is performed for Nt iterations (t for total). Typically, after every 25 iterations the system is sampled and relevant parameters are calculated. If two A segments (on different chains) are nearestneighbors of each other, we consider them to be in the same aggregate. We measure the size of a n aggregate by its aggregation number (Nag.&, and compute numberaverage ( n) and weight-average ( w) aggregation numbers. Isolated surfactant molecules, i.e., “aggregates” with size unity, also called unimers, are included in the calculation of the average aggregation numbers. Large aggregates are referred to as micelles. We determine the shape of a n aggregate by calculating its three principal moments of inertia with respect to its center of mass: I1, Iz, and 13. We then characterize the average shape of such an aggregate by averaging the ratio of the largest and the medium moment of inertia (IllIz) and of the medium and smallest components (Id13). For a perfect sphere I1 = 1 2 = 13,and for a cylinder I I = 1 2 > 13.
In order to compare micelle structures with the SCF model, radial density profiles of the micelles are calculated as follows. The space around the center ofmass of a micelle is divided into concentric spherical shells with a spacing equal to the lattice spacing in the simple cubic lattice. All segments in the first shell are considered to be in layer one, etc. Thus, for the A and B segments, the volume fractions # ( z ) are calculated for z = 1, 2, 3, ... . Self’-ConsistentField Theory. In this approach space is divided into concentric spheres. The spacing between (39) Verdier, P. H.; Stockmayer, W. H. J. Chem. Phys. 1962,36,227. (40) Mehmet, T. G. Macromolecules 1983, 16, 398. (41) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 83, 1087.
3750 Langmuir, Vo2. 11, No. 10, 1995 two spheres equals the diameter of the solvent molecules and the segments which compose the surfactant molecule. Each segment or solvent molecule is located in one of the so-formed layers. All details on a smaller scale than the lattice spacing are automatically neglected. The meanfield approximation is applied within each layer, so that gradients occur in the radial direction only. Within this one-dimensional geometry, all possible conformations of the surfactant chains are generated by treating them as random walks. Each conformation is weighted by a n appropriate Boltzmann factor in order to obtain the equilibrium distribution. Thus, the volume fraction profile is computed of a micelle in equilibrium with a given bulk volume fraction 4b(or chemical potential) ofthe surfactant (eqs A.2-A.5 in the Appendix define the micellar density profile). All details of these calculations have previously been published in the 1iterature.l1-l3 The equations relevant to the results presented in this paper are also given in the Appendix. Once the (radial)density profile of a micelle is known, its excess free energy with respect to the bulk solution, Aexc,can be calculated (using eq A.8). The condition of thermodynamic equilibrium is met if this excess free energy balances the mixing entropy ofthe mice1le.l’ This implies that for a thermodynamically stable micelle Aexc > 0. Furthermore, realistic results may not be expected when A“”‘x 0, as interactions between different micelles are neglected in the SCF model. By applying the theory of small system thermodynamic^,^^ it was shownll that the thermodynamic stability criterion leads to the constraint aAexcl/aOexc > 0, where Oexc is the excess amount of surfactant segments in the micelle with respect to the bulk solution [we define Nagg= Oexc/(NA+ NB)]. Comparison between Lattice Simulationand SCF Theory. The MC results are based on calculations using a simple cubic lattice, whereas the SCF results are based on a spherical variant of the same lattice. So we expect that in both cases the lattice-based nature of the computations will have the same effect on the results. One might expect that the spherical geometry of the lattice used in the SCF model (which is needed to reduce the system to a one-dimensional problem, the mean-field approximation being applied in two dimensions) leads to different results from those found using the threedimensional lattice. However, SCF calculations on polymer chains tethered to a sphere using both the (onedimensional) spherical geometry and a (two-dimensional) cylindrical lattice gave virtually identical volume fraction profiles.43 This reassures us in our assumption that the lattices used in the SCF and MC approaches should both have a very similar effect on the outcome of the calculations. The Markov approximation in the SCF theory is another factor that leads to differences with the simulations. In the MC simulations all chains are modeled as nonoverlapping (self-avoiding) walks. It must be said that more elaborate procedures have been developed and used to evaluate the chain conformation probabilities (the “rotational isomeric state”schemell and the “anisotropicfield” extensionz9of the SCF theory). The random-mixing approximation is a n essential constituent in the SCF theory. This approximation affects both surfactant molecules in the bulk solution and the aggregated state of the surfactant. All spatial density correlations are neglected in the bulk, and in a micelle density correlations within the layers are neglected. As (42)Hall, D. G.; Pethica, B. A. InNonionic Surfactants; M. J. Schick, Ed.; Marcel Dekker: New York, 1976. (43)Wijmans, C. M.; Leermakers, F. A. M.; Fleer, G. J. Langmuir 1994,10, 4514.
Wijmans and Linse we will argue in the Discussion, this approximation makes the prediction of the bulk chemical potential of the surfactant poor, which hence also affects the predicted cmc. Finally, the SCF theory neglects fluctuations in the size and shape of the micelles. For a given surfactant concentration just one spherical micelle, which gives the minimal value for the free energy, is predicted. The theory does not directly provide any information on the width of the micelle size distribution. Neither can distortions from the ideal spherical shape be accounted for. By comparing the SCF calculations with the simulations we hope to get an idea of the effect of this assumption. From the simulations one can directly extract information on the micellar size distribution as well as on micellar shape fluctuations. Up to now we have mainly pointed out the advantages of the MC simulations over the SCF theory. However, the SCF theory has one very important advantage; it has a much smaller computational demand. For low molecular weight systems, as those discussed in this paper, the required central processing unit time is typically 5 orders of magnitude larger for the MC simulations than for the SCF computations. Accurate simulations of polymeric micelles, as studied using the SCF t h e ~ r y , ~ ’ -are ’ ~ ,not ~~ yet possible. Surfactant Model. Two different surfactant models have been used for the comparison between the threedimensional lattice theory and the SCF theory. In the first one (system 11, a surfactant contains 10 lyophobic and 10 lyophilic segments and is characterized by the following parameters: NA= 10, N B= 10, XAB = XAS = 2.7, andXBs = 0. In the second system (system 21, the lyophobic block length is decreased by a factor of 2 and the interaction parameters of the A segments have simultaneously been doubled. Hence,NA = 5, NB= 10, XAB = X A S = 5.4, and XBS = 0. In the simulations, the total number of surfactant molecules are 200 and they are located in a box with L = 44. Our system 1is equal to a model system previously investigated by Wang et al.35 3. Results System 1. In MC simulations of chain molecules, and in particular for systems where self-association can take place, it is necessary to ensure that the system is equilibrated before the sampling is started. Figure 1 shows how the number-average and weight-average aggregation numbers develop during the MC simulation after the initial equilibration with zero contact energies. Going from A to C, the simulation time increases by more than 3 orders of magnitude. From panel A i t is tempting to conclude that the system is equilibrated after 0.05 x lo6iteration steps. This is however not the case; in panel B w is considerably larger. From the even more extended simulation in panel C, we estimate the system to be equilibrated after ca. 0.25 x lo6 iterations. In addition, the normalized weight distributions of the aggregation number of the aggregates (using 0.25 x lo6 iterations for equilibration) are shown in Figure 2 for two simulations differing by 1order of magnitude in length. Although the longer simulation time leads to a curve with higher precision, no systematic shift between the distributions is discerned. From Figure 1 we obtain at equilibrium n =Z 4 and w * 20. This gives a polydispersity ratio wln x 5, which indicates a very broad distribution of the aggregation numbers. The weight (44)Leermakers, F.A.M.; Wijmans, C. M.; Fleer, G. J. Macromolecules 1995,28, 3434.
Modeling of Nonionic Micelles
Langmuir, Vol. 11, No. 10, 1995 3751
10 I
I
30 A
(A)
z
m
t
0 m V
zm V
A
.
- 20
m
z 5
I
5c
0.03 -
a 0.02 -
0.01 -
I
0 '
0 10
5 io4
2.25 104
+
I
, 30
. A
-0
0
v 1.75 I
b
0
C
m 01
5 io5
0
1
--
lo6
10
V
t
1.5
1
o
0 OOo
0
30 A (D
Zm v
+L+
+++
++++++++ t$+
-H
+t+
+.+++?J
++
, " t
'5 +
+++ "++
V
11
+*#a++ ++: +++ 2o
0
20
40
Figure 3. Average ratios of the largest and the medium ( )
0
1 10'
iteration Figure 1. Number-average( ,,) (0)and weight-average (,) (+) aggregation number of the aggregates as a function of simulation length for system 1, starting from a nonequilibrated system (see the text). The total simulation length increases from (A) 0.05 x lo6to (B) 1.25 x lo6t o (C) 12.5 x lo6 iteration steps. Each point represents an average over V ~ OofO the simulation length in each panel. lo6
distribution of the aggregation number (equivalent with the probability of a surfactant to reside in a n aggregate of a given Nagg!is displayed in Figure 2. It shows a wide distribution with maxima a t Nagg= 1and Nagg> > 1,which is typical for micellar solutions. Thus, there are two classes of aggregates: small aggregates including unimers and large aggregates separated by a low probability of finding intermediate-sized aggregates. The aggregation size a t which the minimum probability occurs is 9, and 37% ofthe surfactant molecules are in aggregates smaller than 9. Moreover, the fraction of unimers is 0.17. This corresponds to a volume fraction in the system of 0.008, a value which is the same as that found by Wang et and which can be identified as a measure of the cmc (neglecting the small variation of the concentration of the free surfactant above the cmc). The distribution of micellar aggregation numbers is broad, and the distribution around its maximum a t Nagg= 30 is slightly asymmetric, falling
off more quickly for larger values than for lower values. The size range 10-50 includes 67% of all molecules, and less than 1%of the surfactants are located in aggregates larger than 50. The reason for ? and w being far smaller than 30 is the contribution from the small aggregates (especially unimers). We now consider the average shape of these aggregates. In Figure 3 the average values of the ratios 1142 and I d 3 are given as a function of the micellar aggregation number. Over a large size range is fairly constant a t a value of approximately 1.1. For aggregates in the size range 15-45, that is the range around the maximum micelle fraction, is slightly less than 1.3. This means that these micelles are not perfect spheres, but neither does their shape deviate very drastically from a spherical form (for comparison, a perfect cylinder with a length 1.10 times its diameter has 12/13= 1.30, and an ellipsoid has the same I& ratio when its major axis is 1.26timesitsminor axis). ForN,,, > 50, increases strongly, while remains constant or even decreases slightly. This indicates that these larger micelles have a more rodlike shape. The increased scatter a t large aggregation numbers is due to the low probability of these large aggregates (cf. Figure 2). It should always be remembered that we are only dealing with average shapes. A value of < I I / I 2 > larger than unity can be the result of averaging over a narrow distribution of nonspherical aggregates as well as of averaging over a wider distribution including spherical and nonspherical shapes.
Wijmans and Linse
3752 Langmuir, Vol. 11, No. 10, 1995 1
1
,
,
,
,
,
/
2
3
4
5
6
,
I
7
8
,
,
,
Nz) 0.75
0.5
0.25
0
1
2
3
4
5
6
7
8
z
9
1
0
1
100
k $75
50
25
0
Qb
4 lo5
2 105
0
45
9 1 0 1 1
Z
Figure 4. Volume fractions of A and B segments in micelles (+(z)) as a function of the radial distance ( z ) from the threedimensional lattice model (filled symbols)and the SCF theory (open symbols) for system 1. In the former case the profiles were obtained from micelles with Nagg= 30 by averaging over IO4 micelles, and in the latter case the bulk volume fraction corresponds to the cmc, which leads to Nagg= 45. The curves through the data points are only given as a 'guide to the eye'.
a
'1
1
50 Nagg
Figure 5. (A) Excess surface energy per micelle (Aexc)and (B) the corresponding value of the bulk volume fraction of the surfactant as a functionof the micellar aggregationnumber (Nagg) from the SCF theory for system 1. Figure 4 shows the average radial volume fraction profiles of A and B segments in micelles with Nagg= 30 obtained from the simulation (filledsymbols). As expected, the lyophobicA segments form a dense core around which the B segments are distributed. These B segments are spread out over a large volume, some of them being found fairly near the core center, others protruding far away from the micelle into the solution (up to 10 layers from the micelle center, that is, halfthe contour length of the whole surfactant molecule). The same surfactant system was also studied using the SCF theory. In Figure 5A the excess free energy, A""', is given as a function of the micellar aggregation number. = 45, the excess free energy displays a maximum, For Nagg
Figure 6. Volume fraction of A and B segments in micelles ( ~ z )withNagg ) = 30 (filled symbols)and Nagg = 45 (opensymbols) as a functionof the radial distance (2)from simulation of system 1.
which corresponds to the cmc of the surfactant solution within the framework of the SCF theory. In Figure 5B the surfactant bulk volume fraction q5b is shown for the same micelle sizes. The minimum in this curve, which corresponds to the maximum in Figure 5A, is again the cmc. Thus, from Figure 5 it follows that the critical micelle volume fraction equals and a t this volume fraction micelles are formed with an aggregation number of 45. In Figure 4 the volume fraction profiles a t the cmc predicted by the SCF theory for micelles with Nagg = 45 are also shown (open symbols). The profiles for Nagg= 30 were also for a micelle in a solution close to the cmc.Besides the larger micelle volume (which is due to the larger aggregation number), there are several qualitative differences between the two sets ofprofiles. The SCF predicts a sharper interface of the lyophobic core: the A profile falls off over a smaller distance than in the simulated micelle. Furthermore, there is a (small)dip in the A profile in the center of the micelle, which does not exist in the profile obtained from the three-dimensional lattice. In both theories the B segments are located on the outer side of the lyophobiccore, but the shapes ofthe two B segmental profiles do show a significant difference. The SCF profile is clearly asymmetric. Starting a t the maximum in the profile, the profile falls off far more quickly moving inward (into the core) than outward. The MC simulated profile decays in a similar way in both directions, thus being practically symmetric. A complicating factor when comparing the micellar volume fraction profiles is caused by the different aggregation numbers. An alternative would be to use the same aggregation number for the comparison. Now, the choiceNagg = 30 is not possible, since in the SCF formalism such a n aggregate is too small to be stable (Figure 5). On the other hand, it is possible to use a larger micelle (from the simulation) than the one corresponding t o the maximum in the distribution in Figure 2. In Figure 6 simulated volume fraction profiles of the A and B segments are shown for N,, = 45 as well for Nagg = 30. The occurrence of the larger micelle in the simulation is 1 order of magnitude lower than for the smaller one. The larger micelle, which contains 50%more surfactant molecules, displays a Zower density of lyophobic A segments in the core, and the additional A segments are mainly accommodated in a lyophobic tail sticking out from the center of the aggregate. Moreover, there is a slight increase in the volume fraction of the lyophilic B segments in the core and a considerably enhanced volume fraction of them a t larger distances from the center of the core.
Modeling of Nonionic Micelles
Langmuir, Vol. 11, No. 10, 1995 3753
,
0.08
I 00-
m
-
Zm0.06 -
'", No
.
0.
%
a
0
0..
0.75
0.5
Figure 7. Normalized weight distribution of the aggregation number (Nagg)for system 2. Simulation length Nt = 2.0 x lo6 (open circles) and Nt = 13 x lo6 (filled circles). 3 0
30
40
50 Nagg
Figure 8. Average ratios of the largest and the medium (), and the medium and the smallest ( < I f i s > ) principal for moments of inertia as a function of the aggregate size (Nagg) system 2. Nt as in Figure 7.
Figure 4 showed that there are considerable differences between the volume fraction profiles ofthe micelles at the cmc as predicted by the SCF theory and those of the most likely micelle obtained from the three-dimensional lattice model a t a total volume fraction not too far from the cmc. A comparison between similarly sized micelles (Nagg = 45) from the SCF theory (Figure 4, open symbols) and the simulation (Figure 6, open symbols) shows that the differences between the volume fraction profiles from the two theories seem to be somewhat enhanced compared to the previous comparison a t the cmc. System 2. The results of the second system investigated are given in Figures 7-9, which are the counterparts of Figures 2-4. The larger value of the A-interaction parameters lowers the probability of a MC move wherein a n A-A bond is broken. This makes the motion in the configuration space slower. Therefore, system 2 requires a larger number of iterations in order to become equilibrated (7 x lo6iterations after the interactions are turned on) and to give averages with the same statistical precision. The large scatter in Figure 7 is a n effect of the slower sampling of the configurational space. Disregarding this scatter, one can conclude by comparing the two curves in Figure 7 that one has found the equilibrium distribution. A striking feature of Figure 7 is that there is a clear bimodal distribution of aggregate sizes. There are free unimers (2% of the total surfactant), very low molecular weight aggregates (with aggregation numbers up to four),
t
I
kegments
P
Figure 9. Volume fraction of A and B segments in micelles (+(z)) as a function of the radial distance (z) from the three-
dimensional lattice model (filled symbols) and the SCF theory (open symbols) for system 2. In the former case the profiles were obtained from micelles with Nagg= 45, and in the latter case the bulk volume fraction corresponds to the cmc which = 49.The curves through the data points are only leads t o Nagg given as a 'guide to the eye'. and micelles with aggregation numbers in the range 3050. In contradistinction to system 1 (cf. Figure 2), no aggregates are found in the intermediate regime. Moreover, the volume fraction offree unimers (=cmc)is reduced by 1 order of magnitude and the distribution of micellar sizes is narrower. The cmc from the SCF theory amounts to 2 x 10-8. Figure 8 shows that the average values of IlNz and Id13 are closer to unity ( < 1.1and < 1.15)and depend less on the aggregation number. Hence the shapes of the micelles formed by the A5B10 molecules deviate less from spheres than those formed by the AloBlomolecules. In Figure 9 the average radial volume fraction profiles are given for simulated &B10 micelles withiVag,= 49 (filled symbols). Compared to the simulated profiles of the micelles formed by AloBlo molecules (cf. Figure 41, the core-corona interface of the simulated A5B10 micelle is far sharper. Furthermore, the B profile is clearly asymmetric with respect to its maximum (falling off far more quickly into the core than to the solution side). For the A5B10 surfactant, the SCF theory predicts stable micelles in the size range 49-53, and for comparison, the volume fraction profile forN,,, = 49 is also given in Figure 9 (open symbols). 4. Discussion Comparison of Results from Simulation and the SCF Theory. How can we understand the differences in (i) the cmc, (ii) the aggregation number, and (iii) the micellar structure between the three-dimensional lattice model and the SCF model, and what conclusions can be drawn from these differences? Those are the two questions we want to address in this section. We start by limiting the discussion to system 1. The critical micellar volume fraction predicted by the SCF theory is nearly 3 orders of magnitude lower than and the value obtained from the simulation respectively). Within the framework of the SCF model, the cmc follows directly from the lowest value of the chemical potential of the surfactant molecules for which these molecules can form a (thermodynamically) stable micelle. Concerning the three-dimensional lattice model, Wang et al.35determined the cmc (for our system 1as well as for other systems) by varying the total surfactant volume fraction, and they identified the cmc as the concentration a t the breakpoint of the free unimer
3754 Langmuir, Vol. 11, No. 10, 1995 concentration versus the overall surfactant concentration. The value they reported agrees well with the free unimer concentration in our simulation (9 = 0.008). This should be the case, since the free unimer concentration is only weakly dependent on the total surfactant concentration above the cmc, as the simulations in ref 35 also showed. The breakpoint method is not a very accurate way of determining the cmc, but even if one would allow for a n uncertainty of, say, 20%, this still would not in any way account for the far larger difference as compared to the SCF model. We conclude that the approximations made in the SCF model lead to a far too low cmc, that is, the stability of a micelle is strongly overestimated as compared to free surfactant molecules in solution. The SCF model predicts aggregation numbers between 45 (at the cmc) and 50 (far above the cmc), which is ca. 50%larger than the “optimal” aggregation number of the size distribution obtained from the simulation. In view of the large overestimation of the micellar stability, it is not surprising to find larger aggregation numbers in the SCF model than in the full three-dimensional lattice model. More surprising is the still relatively small increase of the aggregation number. Qualitatively, the segment density profiles of a micelle corresponding to the maximum in the weight distribution of the aggregation number are also rather similar to those obtained from the SCF model. The smaller lyophobiccore is, a t least partly, explained by the lower aggregation number. However, the sharper interface between the micellar core and corona does form a significant feature of a micelle within the SCF model. In contrast, the A profile falls off more gradually and the B profile has a more symmetric shape in the three-dimensional model. In the Introduction several approximations of the SCF model were mentioned which might lead to different results than those of a model with less approximations. First, it does not seem likely that the lattice character of the model can offer a reasonable explanation, as we are comparing the SCF results with results from a lattice model. Second, in the SCF model the surfactant chains are described as (self-intersecting) random walks. This means that these chains have a larger flexibility than a “real” self-avoiding chain. However, this larger flexibility would sooner be expected to counteract than to promote aggregation. The same would be expected for the neglect of correlations between neighboring segment bonds belonging to different molecules. Indeed, Leermakers and L ~ k l e m ashowed ~ ~ that incorporating such correlations in the SCF model leads to a lowering of the cmc and an increase of the aggregation number. It must be noted that the system of ref 29 cannot directly be compared to ours. A different lattice geometry was used, and shortrange interactions between segments belonging to the same chain (“rotational isomeric state’’ scheme) were incorporated in the calculation of chain conformation probabilities. Nevertheless, their findings, which are intuitively very acceptable, indicate that the use of a selfintersecting model for surfactant chains opposes rather than promotes self-assembly. Neither does the constraint imposed by the spherical geometry and the neglect of fluctuations seem to form a reasonable explanation for the differences in the cmc. We expect these approximations both to reduce the tendency of the surfactant to self-assemble and to decrease the aggregation number, but the opposite is observed. Thus, the lower cmc and larger micelles in the SCF model must have a different origin. On the other hand, we do expect the different treatment of the fluctuations to contribute to the differences in the density profiles. The smoother profiles obtained in the simulations arise, since (i) the
Wijmans and Linse average shape deviates from the spherical one (ii) shape fluctuations about the average shape are present and (iii) for a given shape there are radial density fluctuations of the A and B segments. Since the micelles become more nonspherical as their aggregation numbers increase from 30 to 45 (Figure 3), we anticipate that the first effect becomes important for larger aggregation numbers. Indeed, one sees a stronger smoothening of the A profile for the larger micelle (Figure 6). Finally, the mean-field approximation remains as the last important factor determining the differences between the SCF model and the full three-dimensional lattice model. One can easily calculate how the mean-field approximation affects the energy (i.e. the sum of the nearest-neighbor interactions) of a molecule in the bulk solution and in a micelle. In the mean-field approximation, a molecule in the bulk solution a t the critical micelle volume fraction (x < < 1) interacts predominantly ( > ’99%) with solvent. In the micellar state, 24% of the 60 nearest-neighbor contacts of a lyophobic segment are either with a solvent molecule or with a lyophilic B segment. (Due to the neglect of spatial correlations, one also neglects the fact that two nearest-neighbor interactions of a segment are with covalently bonded segments; hence, each A segment can form a maximum of six repulsive nearest-neighbor interactions). In the full threedimensional model, the lyophobic segments of an AloBlo surfactant molecule can in total form 41 nearest-neighbor contacts. A free unimer can lower its energy by forming intramolecular contacts between lyophobic segments, and on average, 5.6 ofthe nearest-neighbor contacts are ofthe A-A type, and the remaining 35.4 (86%)are repulsive contacts with either solvent or B segments. In a micelle with Nagg= 30 there are on average 17.2 repulsive interactions per surfactant molecule (42%). If one only looks a t the energetic contributions to the micellization process, a surfactant molecule has therefore far more to gain from forming micelles in the SCF model than in the three-dimensional lattice model. The transfer of one molecule from the bulk solution to a micelle corresponds to a n enthalpy change of -20.5 and -8.2 kT, respectively. This difference in enthalpy change originates mainly from the different enthalpies in the bulk state and corresponds to a ratio in volume fractions of 2 x lo5, which is even larger than the difference in cmc values that is actually found. As was already indicated, it is reasonable to assume that other factors (the Markov approximation, spherical geometry constraint, and neglect of fluctuations) have a n opposing effect compared to the mean-field nearestneighbor interaction approximation. The conclusions that have been drawn thus far were all based on calculations for one specific system. The calculations for system 2 were performed in order to explore in what manner the chain composition affects the outcome of the modeling. The ASB10 micelles clearly deviate less from the ideal spherical shape than the AloBlo micelles. This trend is expected since a larger head group (B) to tail (A) ratio promotes the formation of a more spherical structure. Furthermore, this larger tendency to form spheres suggests that the SCF model should give a better agreement with the simulation data for the AsBIO molecules than for the AloBlo ones. The micelle size = 49)is indeed closer predicted by the SCF theory (Nagg to the size range found in the simulation than it was for system 1,and the volume fraction profiles are also in closer agreement. Finally, the prediction of the cmc becomes even worse and respectively). However, no firm conclusion can be drawn since the uncertainty in the cmc from the simulation is considerable (the total concentration lies far above the cmc). Nevertheless, one may
Modeling of Nonionic Micelles
Langmuir, Vol. 11, No. 10, 1995 3755
conclude that the large discrepancy in the cmc values between the two models certainly does not disappear in system 2. Micellar Size Distribution and Surfactant Chemical Potential. In the simulation the equilibrium distribution of the aggregation numbers was computed. Figure 2 showed that the normalized weight distribution of the aggregation number for system 1 has a maximum a t Nagg= 30. What factors determine this size distribution? A possible explanation for the reduced probability of micelles with Nagg> 30 could be that in a larger micelle a surfactant molecule is, on average, forced to form less favorable (nearest-neighbor) interactions. One might be inclined to support this hypothesis because of Figure 6, where the core-corona interface is sharper for Nagg= 30 than for Nagg= 45. However, this hypothesis turns out to be false! By counting the number of repulsive interactions during the simulation, one finds that there are 15.6 repulsive interactions per molecule (38%)for Nagg= 45, which is slightly less than the value (17.2)found for the smaller micelle. The fact that micelles of 30 molecules are more favored than those of 45 must therefore have a n entropic cause (there are no other energetic terms). This is an interesting conclusion in the light of a previous discussion in the literature, whether long range (with respect to the attractive forces)repulsive forces are needed for the formation of stable micelles. Pratt et al.31~32 concluded from their simulations that this is the case. However, other studies (SCF c a l c ~ l a t i o n sMD , ~ ~simu~~~ l a t i o n ~ ,and ~ ~ MC , ~ ~~ i m u l a t i o n s do ~ ~ not ) support this conclusion. Our results indicate that the entropy is the most important factor limiting micelle growth. We shall now study the thermodynamics of the micelle distribution in more detail. It is a common practice t o treat a solution of different size micelles as an ideal solution, so for each aggregation number i, the following equation can be written for the chemical potential per surfactant molecules in micelles with Nagg= i:l
1
where pup is the standard chemical potential per surfactant with NJN being the number fraction of aggregates with Nagg= i. The application of this equation can be justified by approximating the aggregates as hard spheres in order to calculate the mixing free energy of the solution. Using the P e r c u s - Y e v i ~ k ~-Carnahan-Stirli~~&~ ~,~~ approximation, a similar equation can be derived for a solution of monodisperse hard spheres with a volume fraction 4 in the limit $r < < l.48 In this case the largest term which has been neglected is linear in 4. With this motivation we apply eq 1to the distributionP,(N,,,) which we calculated for systems 1 and 2. From the micellar size distribution function, the change in the standard chemical potential, pup - p! can be evaluated, and Figure 10 shows pup - p! as a function of the aggregate size. The increase ofthe distribution curves (Figures 2 and 7) as a function of i is clearly accompanied by a decrease of p:. Beyond the maximum in the size distribution, ,uy does level off, but clearly it is not necessary for ,u: to become an increasing function of i in (45) Percus, J. K.; Yevick, G. J. Phys. Reu. 1958, 110, 1. (46) Thiele, E. J. Chem. Phys. 1963, 39, 474. (47) Carnahan, N. F.; Starling, K. E. J. Chem. Phys. 1969,51,635. (48)Israelachvilli, J. N.; Mitchell, D. J.; Ninham, B. W. J. Chem. SOC. Faraday Trans. 2 1976, 72, 1525.
-0.5
-1
I
0
10
20
30
i
40
50
J
60
Figure 10. Change in the standard chemical potential (uy p:) as a function of the aggregation number (i) for system 1and 2 calculated from the micellar size distribution functions using eq 1. O r
-15 I
0
.
'
,
'
10
I
'
'
20
30
i
40
50
I
60
Figure 11. Enthalpic component ofthe change in the chemical potential (hp - h!) as a function of the aggregation number (i) for systems 1and 2 evaluated as the energy difference between a molecule in an aggregate of size n and a free unimer.
order to limit the growth of the aggregates. Different theories of how pe depends on i have been p r o p o ~ e d . ~ , ~ ~ ~ ~ ~ The enthalpic contribution of the change in the standard chemical potential, h: - h!, can easily be determined by counting the number of repulsive interactions of an aggregate. Figure 11shows h: - h! as a function of the aggregate size. Both the standard chemical potential and its enthalpic component are clearly decreasing functions of the aggregation number, over nearly the whole range of the distribution for both systems. However, the enthalpic component shows a far stronger dependence on the aggregation number than the standard chemical potential does. This implies that the entropic contribution of pu," (per surfactant) becomes less favorable with increasing aggregate size. This entropic term originates from the chain conformations, and thus, the chain packing conditions becomes less favorable as the aggregate becomes larger. 5. Conclusions We performed a lattice Monte Carlo simulation to study the equilibrium between different size micelles. Nonionic surfactant molecules were represented as linear chains of segments ("beads") which only interact with each other via nearest-neighbor contacts. This model leads to a polydisperse distribution of aggregate sizes. The great (49)Missel, P. J.; Mazer, N. A,; Benedik, G . B.;Young, C. Y.; Carey, M. C. J. Phys. Chem. 1980,84, 1044.
Wijmans and Linse
3756 Langmuir, Vol. 11, No. 10, 1995 majority of these sizes are significantly smaller than the aggregation number predicted for the same system by a self-consistent field theory. There are also clear differences in the density profile of a micelle as found in the MC simulation compared to the SCF theory. Shape fluctuations of the micelle, which are neglected in the SCF theory, form a reasonable explanation for these differences. Surfactant molecules with a head group twice as large as the tail were seen to form more spherical micelles than surfactants with equal head and tail lengths. As might be expected, the density profiles of the more spherical micelles show a greater resemblance to the profiles calculated within the framework of the SCF theory. A thermodynamic analysis of the size distribution of the aggregates shows that over a wide range of aggregation numbers the relative free energy of a fixed micelle (that is, a micelle without mixing entropy) decreases with increasing aggregation number. This decrease of the micelle free energy is the net result of a (significantly larger) decrease of the enthalpy, which is for a large part (although not completely) compensated by a concomitant decrease of the entropy.
Acknowledgment. We thank Frans Leermakers for a helpful discussion on the numerical details of the SCF calculations. This work was financially supported by the Swedish Research Council for Engineering Science (TFR). Appendix As in the simulation, a surfactant molecule is modeled as a chain of segments, which we numbers = 1,2, NA,..., N A 1, ...,N A Np, =N. The aim is to find the equilibrium segment distribution in a micelle for a given mixture of these surfactant and solvent molecules. Within a selfconsistent field approach, the composition law50r51(known from polymer physics) can be used to calculate the volume fraction of the sth surfactant segment a t r:
+
+
where u,(r) is the potential of mean force with respect to the bulk solution of segment s. The function G plays the role of Green’s function in the diffusion equation. In order to solve eq A.l two simplifylng approximations are introduced. First, the dependence of @, u,and G on r is replaced by a dependence on the scalar r. In other words, it is assumed that a micelle may be described using radial symmetry. In the angular directions a mean-field approximation is now applied. Second, a lattice approximation is used in order to limit the number of surfactant conformations that must be evaluated to find G. The two integrals in eq A . l now become summations:
Here &(r) is the a priori step probability to connect a segment in layer r to one in layer r d . These quantities, A d , are determined by the lattice geometry. Introducing V ( r )= 4/37c$, S ( r )= 4nr2,and the number of sites in layer r (each with a volume equal to the segment diameter squared), one can write
+
1/6 x V ( r - l ) / L ( r ) for d = -1 1/6 x V ( r ) / L ( r ) f o r d = 1 (A.4) = 1 - L 1 ( r )- A+,(?-) for d = 0 lo otherwise The factor C in eq A.2 is proportional to the total amount 0 present in the system, and to the bulk volume fraction @b of surfactant:
For each segmental species p ( p = A, B (surfactant segments), or S (solvent))the segmental potential ofmean force can be written as
up(r)= u’(r) +
d
4s(r)= ( 1 - 4 ~exp(-Pus(r)) ) (A.7) These three segment distribution profiles, in turn, determine the potential profile. Using standard numerical techniques, that potential profile can be found which is consistent with the volume fraction profile to which it gives rise. This must be done under the constraint @A(r) #qJr.) &(r) = 1, for all r. Numerical details can be found in ref 52. For a given segment distribution the excess free energy Aexccan be calculated from the following expression:ll
+
+
r
\
P
d
rl r N
(50) Edwards, S. F. Proc. Phys. SOC.1966,85, 613. (51) De Gennes, P. G. Scaling Concepts i n Polymer Physics; Cornel1 University Press: Ithaca, NY,1979, pp 249.
( X p q x l d ( r4q(r+d)) ) (A.6)
where u’k) is a segment-independent “hard-core”potentia1.52 For a given set of { u p ( r )one } can use eqs A.2-A.5 to calculate the surfactant segment distribution, &(r) and @&I. The solvent distribution is:
4(r;s)= C exp(Bu,(r))CCG(r;str,)G(r;slrN) (A.2) and the discrete functions G obey the following relations:
C q=A,B,S
(A.8) LA9503865 (52) Evers, 0.A.; Scheutjens, J. M. H. M.; Fleer, G. J. Macromolecules 1990,23,5221.
