Monte Carlo Simulation Study of Concentration ... - ACS Publications

2922

Langmuir 2002, 18, 2922-2932

Monte Carlo Simulation Study of Concentration Effects and Scattering Functions for Polyelectrolyte Wormlike Micelles Luigi Cannavacciuolo,† Jan Skov Pedersen,‡ and Peter Schurtenberger*,§ Institut fu¨ r Polymere, Eidgeno¨ ssische Technische Hochschule, CH-8092 Zu¨ rich, Switzerland, Department of Chemistry, University of Aarhus, DK-4000 Aarhus, Denmark, and Department of Physics, University of Fribourg, CH-1700 Fribourg, Switzerland Received June 12, 2001. In Final Form: December 21, 2001 We present an extensive Monte Carlo study on systems of many semiflexible chains with excluded volume and electrostatic interactions within the Debye-Hu¨ckel approximation. The model is tuned to mimic charged wormlike micelles in solution under different conditions. Simulations have been performed at different ionic strengths of added salt, charge densities, chain lengths, and volume fractions (η), covering the dilute to concentrated regime. The simple model used is able to reproduce the structural peak of the scattering function S(q), observed in many experiments, and other important features of polyelectrolytes in solution. Universal behavior of S(0) is established after a rescaling of η. Single-chain scattering functions, radii of gyration, and end-to-end distances have been sampled and compared with data from previous simulation studies and theories. The persistence length has also been analyzed and compared with the Odijk-Skolnick-Fixman predictions. The behavior of the quantities studied is in general found to be more complex than scaling theory predictions.

I. Introduction Despite the considerable amount of work that has been done in the last two decades, the understanding of polyelectrolytes (PEL) in solution is still far from being satisfactory, particularly if compared to the knowledge of uncharged polymers.1 The main difficulties arise from two fundamental features of polyelectrolytes: the long-range nature of the Coulomb interactions and the influence not only of the solvent, as in the case of uncharged polymers, but also of counterions and, eventually, added salt ions. This makes both theory and experiments particularly challenging, even for the study of static properties. From a theoretical point of view, the Coulomb potential can only be treated analytically in a mean field manner. Even for weakly charged chains, the electrostatic contribution to the total energy is in general nonperturbative, because due to their long range, interactions can sum up and dominate the thermal energy at some length scale, unless screened out by counterions and, when present, by added salt ions. With respect to the solvent molecules, it is reasonable to imagine the effects of their discreteness to be confined on a much shorter length scale compared to the polyion length. Therefore, apart from a few ab initio works,2-4 it is normal practice to average over the solvent * To whom correspondence should be addressed. † Eidgeno ¨ ssische Technische Hochschule. ‡ University of Aarhus. § University of Fribourg. (1) For a review, see for example: Dautzenberg, H.; Jaeger, W.; Ko¨tz, J.; Philipp, B.; Seidel, Ch.; Stscherbina, D. Polyelectrolytes: Formation, Characterization and Applications; Hanser Publishers: Munich, 1994. Macro-ion Characterization: From Dilute Solution to Complex Fluid; Schmitz, K. S., Ed.; ACS Symposium Series Vol. 548; American Chemical Society: Washington, DC, 1994. Fo¨rster, S.; Schmidt, M. Adv. Polym. Sci. 1995, 120, 50. Barrat, J. L.; Joannny, J. F. Adv. Chem. Phys. 1996, XCIV, 1. (2) Clementi, E. J. Phys. Chem. 1985, 74, 4426. (3) Clementi, E.; Lie, G. C. In Supercomputer Research in Chemistry and Chemical Engineering; Jensen, K. F.; Truhlar, D. G., Eds.; ACS Symposium Series Vol. 353; American Chemical Society: Washington, DC, 1987; p 237. (4) Levesque, D.; Weiss, J. J. In Monte Carlo Method in Condensed Matter Physics; Binder, K., Ed.; Topics in Applied Physics Vol. 71; Springer: Berlin, 1995; p 121.

degrees of freedom and to treat it as a dielectric continuum. This is in general not done for co-ions and counterions, whose discrete nature can produce significant fluctuations of the charge density, which must be taken into account in the framework of a general theory. Scaling theories and renormalization group calculations, so successfully applied for the study of uncharged polymers,5 are based on the assumption that the range of the interactions is much smaller than the length scale determining the macroscopic properties of the system. This is, of course, not the case for polyelectrolytes, for which the Coulomb interactions are long ranged and, even within the Debye-Hu¨ckel approximation, can give rise to a Debye length of the same order of magnitude as the chain length. Due to the presence of additional characteristic length scales, scaling theory for polyelectrolytes is much more difficult and speculative and it is not really clear how a coarse-graining procedure should be carried out. The electrostatic repulsion between charged beads produces a swelling of the chain if compared to the corresponding uncharged case, and therefore a decrease of the overlap concentration should be expected. As a consequence, it usually happens that the corresponding dilute regime falls well beyond the detection sensitivity of the presently available scattering instruments. This makes it almost impossible to determine the structure of polyelectrolytes in the dilute regime where true singlechain properties are probed. This problem is particularly relevant considering that most of the theoretical works deal with approximations valid only in that regime. It turns out to be extremely difficult to perform experiments for checking the validity of the fundamental theoretical assumptions and approximations. To overcome this difficulty, a new model system for polyelectrolytes in solution, based on charged wormlike micelles, has recently been proposed and employed.6-8 The stronger scattering power (5) de Gennes, P. G. Scaling Concepts in Polymer Physics; Cornell University Press: Ithaca, NY, 1979. (6) Magid, L. J. J. Phys. Chem. B 1988, 102, 4064. (7) Jerke, G.; Pedersen, J. S.; Egelhaaf, S. U.; Schurtenberger, P. Langmuir 1998, 14, 6013.

10.1021/la010884f CCC: $22.00 © 2002 American Chemical Society Published on Web 03/09/2002

Simulation of Polyelectrolyte Wormlike Micelles

of such micellar solutions, compared with that of standard polyelectrolytes, allows one to explore the very dilute regime and to extrapolate single-coil properties. Moreover, by simply adding small amounts of ionic surfactants at various salt concentrations to solutions containing giant polymer-like nonionic micelles, it is possible to tune the charge density and screening effects to reproduce polyelectrolytes under analogous conditions. With the present work, we try to give more insight into the analogy between PEL and charged wormlike micelles and on how the micelles can be modeled. A powerful tool to manage the specific problems posed by polyelectrolytes turned out to be computer simulations, which have already extensively and successfully been applied to uncharged polymers. In fact, computer simulations provide us, at least in principle, with results under very controlled situations, which can never be reached in experiments. Polyelectrolyte simulations are, however, more difficult and expensive in CPU time compared to studies of neutral polymers. The calculation of the electrostatic energy takes the longest time per Monte Carlo (MC) step, time in principle being proportional to the square of the total number of particles in the system. A complete description of polyelectrolytes in solutions with added salt should include all interactions between polyions, counterions, co-ions, and solvent molecules. Moreover, the study of universal properties and the determination of critical exponents requires the use of large systems and long chains. The fulfilment of all these requirements is beyond the capability of present computers, and therefore only idealized situations can presently be treated. The literature on many-chain polyelectrolyte simulations begins only in the past decade,9-12 and even if progress has been made, it is still modest, and the fundamental issues of polyelectrolytes in solutions are still obscure. Stevens and Kremer have performed molecular dynamics simulations of many-chain systems using a continuum model in the Debye-Hu¨ckel approximation13 and also including the full Coulomb interactions between monomers and counterions.14-18 The work reported in the present paper is an extension of our previous simulations of single chains of charged beads.23 We present an extensive MC study of systems of many chains with excluded volume effects (EV) and with electrostatic interactions treated in the Debye-Hu¨ckel approximation. The parameters defining the model are chosen to reproduce charged wormlike micelles. In particular, the numerical values have been chosen to closely reproduce the situations as in our previous experiments on a mixture of hexaethylene glycol mono-n-hexadecylether (C16E6) and a small amount of the ionic surfactant 1-hexadecane sulfonic acid (C16SO3Na) in very dilute solutions, with added salt (NaCl), which give rise to (8) Sommer, C.; Cannavacciuolo, L.; Egelhaaf, S. U.; Pedersen, J. S.; Schurtenberger, P. Prog. Colloid Polym. Sci. 2000, 115, 347. Sommer, C.; Pedersen, J. S.; Egelhaaf, S. U.; Cannavacciuolo, L.; Kohlbrecher, J.; Schurtenberger, P. Langmuir, in press. (9) Hasalam, A. J.; Jackson, G.; McLeish, T. C. B. J. Chem. Phys. 1999, 111, 416. (10) Ullner, M.; Staikos, G.; Theodorou, D. N. Macromolecules 1998, 31, 7921. (11) Scha¨fer, H.; Seidel, C. Macromolecules 1997, 30, 6658. (12) Christos, G. A.; Carnie, S. L. Chem. Phys. Lett. 1990, 172, 249. (13) Stevens, M. J.; Kremer, K. J. Phys. II 1996, 11, 1607. (14) Stevens, M. J.; Kremer, K. J. Chem. Phys. 1995, 103, 1669. (15) Stevens, M. J.; Kremer, K. In Macro-ion Characterization: From Dilute Solution to Complex Fluid; Schmitz, K. S., Ed.; ACS Symposium Series Vol. 548; American Chemical Society: Washington, DC, 1994. (16) Stevens, M. J.; Kremer, K. Phys. Rev. Lett. 1993, 71, 2228. (17) Stevens, M. J.; Kremer, K. Macromolecules 1993, 26, 4717. (18) Kremer, K.; Du¨nweg, B.; Stevens, M. J. Physica A 1993, 194, 321.

Langmuir, Vol. 18, No. 7, 2002 2923

charged polyelectrolyte-like micelles,8 so-called equilibrium polyelectrolytes. To achieve an exhaustive and satisfactory set of results, according with the large variety of possible conditions for polyelectrolytes in solution, we have simulated systems at various volume fractions η, ionic strengths I, charge densities Zb, and different chain lengths. To better understand the effects of electrostatic interactions, uncharged chains have been simulated as well, at the same volume fractions. Following the same guidelines as in our previous work, we mainly focus on the analysis of the scattering function, because of its importance in the interpretation of experimental data and on size and flexibility properties. A further reason which has motivated us to produce a quite complete set of simulations is the idea to end up with a full parametrization of the scattering function, as already achieved for uncharged systems,19,20 thus providing us with a fundamental model to analyze experimental scattering data. Despite the severe approximations of the continuous model used, we have been able to reproduce wellestablished experimental results like the appearance of structural peaks in the scattering function when interchain interactions are strong enough. Moreover, since most of the theories are based on the Debye-Hu¨ckel approximation, we believe that the present work provides a wide set of data which can be compared with previous theoretical approaches. Of course, other important and still less understood features of polyelectrolytes in solution, like the counterion condensation and the charge distribution around the polyion, deeply involve the exact nature of the particles and need different approaches to be treated, and the study of these effects is beyond the scope of the present work. II. Model and Simulation Technique The model for the individual chains we use in the present work is based on a discrete representation of the wormlike chain (WLC) model by Kratky and Porod.21 It consists of a freely rotating chain of N points separated by the distance a and equal valence angle θ. The chain contour length is thus L ) (N - 1)a, and the intrinsic Kuhn length b is given by

1 + cos θ b)a 1 - cos θ

(1)

The only degrees of freedom are rotations of the points around the connecting bonds, that is, the choice of the torsion angles. The valence angle is throughout the present simulations chosen so that there are six spheres per Kuhn length, that is, θ ) arccos(5/7). Note that in contrast to the work described in refs 19, 22, and 23, we did not perform an extrapolation of the results to the continuous limit of the model, as this is not possible due to the relatively large size of the studied systems and due to the relatively long range of the interactions (see further below). Excluded volume effects are included by placing N hard spheres of radius F ) 0.1b ) 0.6a at the N points. For the individual chains, configurations with overlap of spheres (19) Pedersen, J. S.; Schurtenberger, P. Macromolecules 1996, 29, 7602. (20) Pedersen, J. S.; Schurtenberger, P. Europhys. Lett. 1999, 45, 666. (21) Kratky, O.; Porod, G. Rev. Trav. Chim. 1949, 68, 1106. (22) Pedersen, J. S.; Laso, M.; Schurtenberger, P. Phys. Rev. E 1996, 54, R5917. (23) Cannavacciuolo, L.; Sommer, C.; Pedersen, J. S.; Schurtenberger, P. Phys. Rev. E 2000, 62, 5409.

2924

Langmuir, Vol. 18, No. 7, 2002

Cannavacciuolo et al.

separated by more than two bonds are disregarded. The value of F has been chosen on the basis of previous experimental observations on the local structure of wormlike micelles24 and on the fact that the uncharged model describes the structure and excluded volume effects of both wormlike micelles and polystyrene in a good solvent.19,20,24 The chains have 60-540 spheres, and the number of chains depends on the length of the chains, so that the number of chains is lower for long chains. The total number of spheres is typically 3000-10000. Note that for simplicity we have only considered strictly monodisperse systems. With only hard-sphere interactions, the model provides a completely athermal representation of semiflexible chains. This was the reason for using it in previous works.19,20,22 However, other representations, for example, with a bending potential for the valence angle, are also possible. The electrostatic interactions are included by associating a charge ZbL/(Nb) to each sphere, where Zb is the number of charges per Kuhn length. We have chosen to use 25, 50, and 75 elementary charges per Kuhn length, which are the values that should closely correspond to the situation encountered in our experiments on polymerlike micelles.7,8 We do not take into account explicitly the presence of co-ions and counterions in this study; the only electrostatic interactions are those between fixed charges on the polyions, screened by clouds of co-ions of the added salt and counterions. This is mathematically expressed by the “crude” approximations of the (Debye-Hu¨ckel) Coulomb interaction energy:

{

(ZbLe/Nb)2 e2F/λD e-r/λD for r > 2F 2 U(r) ) 4π0r (1 + F/λ ) r D ∞ otherwise

(

)

0rkBT 2e2I

Dee ) [〈(RN - R1)2〉]1/2

1/2

(2)

where kB is Boltzmann’s constant, T is the absolute temperature, and I is the ionic strength in mol/L, including the contributions of both added salt and counterions. To relate the distances in the potential and the hard-core radius, we have used the Kuhn length b ) 300 Å which is close to the experimental value.7,8 The electrostatic energy for the individual chains has been calculated for spheres separated by more than two bonds along the contour. Formally, the system is placed in a box with periodic boundary conditions. The side length of the box is chosen so that the desired volume fraction η is obtained: Lbox ) (V/η)1/3, where V ) NcN4πF3/3 is an estimate of the volume of the Nc chains. In the calculations of the electrostatic energy and in the check for overlap, the minimum image convention was used.26 For the electrostatic interactions, (24) Jerke, G.; Pedersen, J. S.; Egelhaaf, S. U.; Schurtenberger, P. Phys. Rev. E 1997, 56, 5772. Stradner, A.; Glatter, O.; Schurtenberger, P. Langmuir 2000, 16, 5345. Pedersen, J. S.; Egelhaaf, S. U.; Schurtenberger, P. J. Phys. Chem. 1995, 99, 1299. (25) Manning, G. S. J. Chem. Phys. 1969, 51, 924. (26) Allen, M. P.; Tildesly, D. J. Computer Simulation of Liquids; Oxford Science Publications: Oxford, 1986.

(3)

with Ri being the positions of the beads, and the radius of gyration

Rg )

In the above expression, e is the elementary charge, 0 is the permittivity of vacuum, r is the dielectric constant of the solvent, and λD is the Debye screening length. The solvent is treated as a dielectric continuum with a dielectric constant r ) 78.5 (H2O at 25 °C), assumed constant at all values of η. The Debye screening length is defined by

λD(I, Zb, η) )

we used a cutoff at one per mille of the thermal energy. The cutoff allows us to use “zippering”27 in connection with the energy calculations and the minimum image convention. The nonoverlapping initial configuration was generated by growing the chains one by one and at each step checking for sphere overlap. At each step, 10 attempts were allowed. If it was not possible to continue the growth after 10 attempts, the chain was deleted and started at a new position. The trial configurations were generated by a reptation method. A polymer was chosen at random, and one end was selected at random. nrep spheres were deleted from this end and then attached to the other end by generating a random configuration. The value of nrep was chosen after some preliminary tests to ensure that the acceptance fraction f of the moves is about 50% as this gives an efficient sampling. The value of nrep therefore strongly depends on the volume fraction. The trial configurations were accepted according to the Metropolis criterion.28 The initial biased configuration was thermalized for neq ) 3Nc(N2/nrep2 + 1) moves before sampling started. Examples of chain configurations after thermalizations are shown in Figure 1. During the simulations, the end-to-end distance

[〈

1

〉]

1/2

N

∑

N i)1

(Ri - RCM)2

(4)

where RCM is the position of the center of mass of the chain, given by

RCM )

1

N

∑ Ri

Ni)1

(5)

were sampled. In these expressions, 〈‚‚‚〉 denotes ensemble average. We have also sampled the full system scattering function S(q) and the single-chain scattering function P(q), where q is the modulus of the scattering vector. For the calculations, we used a method similar to the one described by Frenkel et al.,29 who calculated the scattering in 001, 110, and 111 directions using fast Fourier transformation (FFT). The functions are calculated in the reciprocal lattice points as this agrees with the periodic boundary conditions. However, we have found that the discretization involved in the FFT is significant in the range of qb e 10 which we use. We have therefore used an alternative approach in which the addition rules for sines and cosine are employed as this reduces the number of time-consuming calculations of these functions. Note furthermore that the fact that we are using high-symmetry directions allows the functions in the 110 and 111 directions to be calculated from the sine and cosines calculated for the 001 directions. The minimum q value sampled in the simulations is q ) 2π/Lbox. As the periodic boundary conditions are irrelevant for the single-chain properties, we also calculated the single-chain scattering functions P(q) with the (27) Stellman, S. D.; Froimowitz, M.; Ganz, P. J. J. Comput. Phys. 1971, 7, 178. (28) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, A. H.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (29) Frenkel, D.; Vos, R. J.; de Kruif, C. J.; Vrij, A. J. Chem. Phys. 1986, 84, 4625.


Figure 1. Example of chain configurations in the simulation box at two different ionic strengths (I). Chains have an overall length of L/b ) 10, with Zb ) 75 charges per Kuhn length. The volume fraction is η ) 0.005. (a) I ) 0.1 M, (b) I ) 0.001 M. Note the difference in the chain elongation and in the average distance between different chains due to the electrostatic repulsion.

same method in an arbitrarily large unit cell to make sure that the function was also sampled in the Guinier region, q2Rg2 , 1. The end-to-end distance and the radius of gyration were sampled at every move for the chain on which the move was attempted. The scattering functions were only sampled for every 100 moves. Statistical uncertainties on the sampled quantities were determined by block analysis.30 III. Monte Carlo Results Simulations have been performed covering the dilute, semidilute, and concentrated regimes, with η in the range of 10-4-0.2,31 at three ionic strengths I ) 0.1, 0.01, and 0.001 M and three charge densities Zb ) 25, 50, and 75. (30) Flyvbjerg, H.; Petersen, H. G. J. Chem. Phys. 1989, 91, 461.

Langmuir, Vol. 18, No. 7, 2002 2925

Note that for the values Zb ) 50 and 75 one obtains a Manning parameter ξ ) λBZb/b25 that is slightly larger than unity and therefore, in the case of classical thin rodlike PEL counterions, condensation occurs. However, this criterion is not directly applicable to wormlike micelles as they possess a finite cross section. As the charge is distributed on a cylindrical surface and not on a onedimensional structure, the threshold for counterion condensation is enhanced and, under the conditions used in our simulations, such a phenomenon is not expected to occur. This is also confirmed by the agreement between experimental data and MC results. The typical chain length used in the simulations was N ) 61 and nb ) 6, where nb is the number of spheres per Kuhn length, given by nb ) b/a. At some setting of the other parameters, values of N ) 7, 19, 181, 541 have been used, whereas nb has been kept fixed. The number of MC steps and of chains used ranges from 5 × 107 steps and 50 chains for the shortest chains to 2 × 108 steps and 10 chains for the longest ones. For the values of the various parameters used, the resulting Debye length and box size fulfill the condition Lbox/λD J 50 . 1 that guarantees the system to be consistently defined in the simulations. A. Size and Shape. As mentioned above, we have sampled the root-mean-square values of the radius of gyration Rg and of the end-to-end distance Dee. In addition, we sampled the root-mean-square end-to-middle distance Dem and the root-mean-square distance between two inner points Dii situated at L/4 and 3L/4 along the chain, as well as the respective distribution functions (only data of Rg and Dee are reported in the present paper). Figure 2 reports the results of the simulations for Rg and Dee at different ionic strengths and charge densities. Both of these quantities are almost constant at low volume fractions and decrease at high volume fractions due to screening effects. The screening of the Coulomb potential has two origins: the presence of salt ions and of the counterions, where the latter increase as the volume fraction is increased. Therefore, the Debye length λD turns out to be dependent on η. The behavior of λD for the various values of the ionic strength and charge density is shown in Figure 3. Returning to Figure 2, we see that as η f 0, the screening decreases, and Rg and Dee slowly approach their maximum values. The chains reach their largest stretching, consistent with the behavior of λD (see Figure 3). In the other limit, where η increases, also the screening effect increases; consequently, Rg and Dee decrease, and the chain conformations are coil-like. Note that for the highest charge density we find the steepest variation of Rg and Dee which again is consistent with the behavior of the Debye length. At a volume fraction of about ∼0.1, electrostatics are totally screened out for all sets of simulations, and chains behave as uncharged. Note that at this concentration the screening length at I ) 0.1 and 0.01 M is less than the radius of the hard spheres (30 Å) for all values of Zb. However, looking at Figure 3, it seems that at I ) 10-3 M complete screening should not occur yet. We expect therefore that an additional screening of the interactions on any given chain is produced by the other chains surrounding it. It has previously been shown by Muthukumar that both the bare EV and the electrostatic interactions between any two segments of a given chain are screened by the presence of other chains at nonzero concentrations.32 The same (31) We are aware that the Debye-Hu¨ckel approximation ultimately breaks down at high concentrations. Nevertheless, we think that it is useful to have simulation data also for this regime that can be compared with full Coulomb simulations in the attempt to better understand the range of validity of the approximations. Of course, our results obtained at high values of η should be considered very carefully.

2926



Figure 3. Debye length vs η at different I and Zb. From above to below: I ) 0.001 M, Zb ) 25, 50, 75; I ) 0.01 M, Zb ) 25, 50, 75; I ) 0.1 M, Zb ) 25, 50, 75.

of screening. The concentration threshold below which “saturation” of Dee occurs can therefore be identified with the absence of counterions and other chains inside the individual polyelectrolyte volumes. At higher concentrations (η ∼ 0.2), an interesting effect occurs. Charged chains are less expanded than uncharged ones. The effect is very small but systematic and is demonstrated in the insets of Figure 2a,b.31 Odijk33 has proposed a scaling theory for polyelectrolyte solutions in the presence of salt. The main results are the scaling of the overlap concentration, given by

c* = Lλc-1RF-3

(6)

and the scaling of the end-to-end distance

Dee ∼ L1/2λB-1/8λD3/8(λcc)-1/8

Figure 2. (a) Radius of gyration and (b) end-to-end distance vs η at different ionic strengths and charge densities. Symbols: I ) 0.001 M, Zb ) 75 (]), Zb ) 50 (cross), Zb ) 25 (4); I ) 0.01 M, Zb ) 25 (0); I ) 0.1 M, Zb ) 25 (O). Inset: Zb ) 25, I ) 0.001 M (4), I ) 0.01 M (0), I ) 0.1 (O), uncharged (f). At η ∼ 0.2, charged chains are more contracted than uncharged ones. The broken lines are values for an ideal wormlike chain of the same contour length (L/b ) 10). (c) Shape parameter G vs η. The symbols are the same as in (a) and (b). The two horizontal lines correspond to the limiting values G ) 6.3 (good solvent) and G ) 6 (ideal).

qualitative behavior of Dee is found in the salt-free simulations of Stevens and Kremer,17 where counterions are explicitly taken into account and are the only source

(7)

where λB is the Bjerrum length, RF is the Flory radius, λc is the distance between charges, and c is the concentration. Equations 6 and 7 are subject to the restrictions that an excess of added salt is present and that the c range is confined to the semidilute regime defined by λc-3 . c . c* (for a complete description of the assumption, see ref 33). It is easily seen that our data do not agree with (7). In our data, the exponents of both λD and η vary continuously and interdependently, for example, the exponent of η is found to vary continuously with I. One could think the disagreement to be a finite size (L/b) effect due to the relatively short chains used. However, we are more inclined to believe that our data contain indications for more complex laws governing the behavior of the quantities under consideration. Our perceptions are based mainly on two reasons. First, in their simulations of saltfree polyelectrolyte solutions Stevens and Kremer17 found Dee to follow a scaling behavior with respect to N that is compatible with the Odijk predictions even with shorter chains than those used in this work, and therefore they had negligible finite size effects. Second, the fact that Dee does not scale with a simple power law qualitatively agrees with the conclusions of the renormalization group calculations of Kholodenko and Freed.34 In their work, the (32) See, for example: Muthukumar, M. J. Chem. Phys. 1996, 105, 5183 and references therein. (33) Odijk, T. Macromolecules 1979, 12, 688. (34) Kholodenko, A. L.; Freed, K. F. J. Chem. Phys. 1983, 78, 7412. Bawendi, M. G.; Freed, K. F. J. Chem. Phys. 1986, 84, 449.


Langmuir, Vol. 18, No. 7, 2002 2927

Figure 4. Effects of concentrations and of ionic strength on the scattering function S(q) at fixed Zb ) 25: (a) I ) 0.001 M, from above to below η ) 10-4, 0.001, 0.002 0.003, 0.005, 0.008, 0.01, 0.02, 0.04, 0.1, 0.2; (b) I ) 0.01 M, η ) 10-4, 0.002, 0.003, 0.005, 0.008, 0.01, 0.02, 0.04, 0.1, 0.2; (c) I ) 0.1 M, at the same concentrations as in (b); (d) uncharged at the same concentrations as in (b). The symbols are sampled data, and the full lines are fits.

authors point out that no universality can be established for PEL solutions with added salt. However, a generalized scaling expression for the end-to-end distance is derived in which the scaling exponent is a function of two other scaling variables. The conformation of a chain can be characterized by means of the shape parameter G ) 〈Dee2〉/〈RG2〉. This quantity is very sensitive to the stretching of a chain; it varies from 6 (ideal chain) to 6.3 (flexible chain in good solvent) to 12 (rigid rod). Figure 2c shows the values of G obtained from the simulations. At low concentration, the chains are more stiff, but even for the most interacting setting (at I ) 10-3 M, Z ) 75), they are still very far from being rodlike. At a high volume fraction, chains gradually approach the good solvent and ideal condition, and at even higher concentration they behave as ideal chains due to the screening of excluded volume effects.35 B. Scattering Function and Related Quantities. Figure 4 shows the results of the scattering function at different ionic strengths and volume fractions. Two main effects are evident. At all ionic strengths, as well as for the uncharged system, S(0) decreases with increasing concentration due to interchain interactions. As expected, this effect is much more pronounced when the electrostatic interactions are stronger, that is, at lower ionic strength (35) Flory, P. J. Principles of Polymer Chemistry; Cornell University Press: Ithaca, NY, 1954.

and at higher charge density (see also Figure 5). Physically, this corresponds to an increase of the osmotic compressibility of the system. For strongly interacting chains, a peak in S(q) is observed. This is a well-known feature also observed in experiments on both PEL and micellar solutions. It is a sign of formation of a local order in the system. In real space, a corresponding maximum of the pair correlation function g(r) is observed at distance r ) d, which can be interpreted as an average distance for a shell of neighboring spheres. On increasing concentrations, the peak position moves toward higher scattering vectors, as the chains are on average closer to each other. A similar shift of the peak position is observed at fixed η and with increasing I. In this case, the decrease of d is due to the weakening of the electrostatic repulsion as the screening is increased, and therefore chains are allowed to come closer. The chain length dependence of the scattering function is shown in Figure 6. For L/b = 3, S(q) turns out to depend only slightly on the chain length at relatively high η values. Moreover, in the same regime, the peak position seems to be independent of L/b. An expression for the scattering function of uncharged polymers has been derived by Curro and Schweizer, who developed the polymer reference interaction site model (PRISM),36 based on the liquid-state RISM theory of Chandler and Andersen37 for solutions of diatomic mol(36) Schweizer, K. S.; Curro, J. G. J. Chem. Phys. 1987, 87, 1842.

2928



Figure 6. Effects of the chain length on S(q) at fixed η ) 0.02, I ) 0.001 M, and Zb ) 50. Symbols: L/b ) 1 (O), L/b ) 3 (0), L/b ) 5 (4), L/b ) 10 (]), L/b ) 30 (cross). The symbols are sampled data, and the full lines are fits.

transform of the direct correlation function c(r) of OZ, and P(q) is the single-chain scattering function normalized to P(0) ) 1.39 Due to its shorter range compared to the full correlation function, c(r) is a more convenient function to approximate. In the simulations, we have sampled both S(q) and P(q) at the same q values. One can, therefore, calculate c(q) from eq 8 as

c(q) )

(

1 1 1 β S(q) P(q)

)

(9)

We found that the function

sin(qRc) (-q2σ2) e qRc

βc(q) ) A

(10)

which is qualitatively similar to the direct correlation function derived in PRISM theory,42 is able to reproduce c(q) with reasonable accuracy. In (10), A, Rc, and σ are fit parameters. It turned out that only two of them are independent. Indeed, we have found that the following relation holds within the accuracy of the data for all the sets of data investigated, independently of η, I, and Zb:

Rc ) 2 log σ + 2

Figure 5. Effects of the charge density on S(q) at fixed η ) 0.02 and different Zb ) 25, 50, 75: (a) I ) 0.001 M, from above to below Zb ) 25, 50, 75; (b) I ) 0.01 M, same Zb as in (a); (c) I ) 0.1 M, same Zb as in (a). The symbols are sampled data, and the full lines are fits.

ecules.38 The PRISM approach consists of a modified Ornstein-Zernike equation (OZ), which has a particularly simple expression in Fourier space,

S(q) )

P(q) 1 + βc(q) P(q)

(8)

where β ) 1/S(0) - 1, c(q) is the three-dimensional Fourier

(11)

Examples of Rc are given in Figure 7. Note that Rc is poorly defined for η f 0. One can thus construct a fit function for S(q) by means of eq 8 with c(q) given by (10) and (11) using the sampled P(q) functions. We have succeeded in fitting this function to the sampled scattering functions. Examples of such fits are reported in Figures 4, 5, and 6. We want to stress that the procedure adopted for fitting S(q) is not based on rigorous arguments, and therefore it should by no means be seen as an attempt to check existing (37) Chandler, D.; Andersen, H. C. J. Chem. Phys. 1972, 57, 1930. (38) Recent works (refs 40 and 41) extend PRISM theory to PEL solutions. Unfortunately, those results cannot be employed for our purpose because of the different length scales involved due to the fact that we are treating a micelle (semiflexible cylinder with finite cross section) in salt solution and not a conventional PEL and because these previous studies consider limiting situations (salt free or rigid rod) only. (39) Note that in the figures we use a normalization with S(q ) 0) ) N. (40) Shew, C.-Y.; Yethiraj, A. J. Chem. Phys. 2000, 113, 8841. (41) Hoffmann, T.; Winkler, R. G.; Reineker, P. J. Chem. Phys. 2001, 114, 10181. (42) Schweizer, K. S.; Curro, J. G. Chem. Phys. 1990, 149, 105.


Figure 7. Values of the fit parameter Rc vs η at I ) 10-3 M, Zb ) 75 (open f), Zb ) 50 (0), Zb ) 25 (O); I ) 0.01 M, Zb ) 75 (]), Zb ) 50 (f), Zb ) 25 (cross); I ) 0.1 M, Zb ) 75 (4); uncharged (1).

Langmuir, Vol. 18, No. 7, 2002 2929

Figure 9. Position of the peak of S(q) vs η from simulation data at I ) 0.001 M, Zb ) 25 (0), 50 (4), and 75 (O), and experimental data on C16E6 doped with 3% ionic surfactant C16SO3Na (9) and with 6% of ionic surfactant (2). The straight line has a slope of 1/2.

This is predicted from most theories43 and also reported in many experimental works. In Figure 9, we also compare the MC data with our previous experimental results on an a mixture of nonionic surfactant and 3% or 6% of ionic surfactant.8 The MC data reproduce the correct scaling exponent and also the position of the crossover at low concentration. Only the nonuniversal prefactor is, probably, slightly different. The agreement shows that the simple model adopted in the simulations is able to reproduce general properties of the investigated system. D. Forward Scattering. An important quantity to consider is the forward scattering, because it is directly related to the osmotic compressibility Π of the system through43

(dΠ dc )

S(0) ) kBT Figure 8. Example of c(q) (symbols) defined by eq 9 with S(q) and P(q) sampled, at η ) 0.01, I ) 0.001 M, and Zb ) 25. The full line is a fit of function 10 to the data.

theories or to provide new closure relations. Indeed, we disregard smaller deviations in reproducing the c(q) but only focus on the high-quality fits obtained for the S(q) data, keeping in mind that only S(q) is an experimentally accessible quantity. The fitted S(q) curve allows one to determine S(0) and the position of the peak, when present. Figure 8 shows an example of c(q) data at I ) 0.001 M, η ) 0.01, and Z ) 25, plotted together with the values given by eq 10, where the fit parameters are those obtained by fitting the corresponding S(q) data with the procedure described above. We have also been able to reproduce the S(q), using for the P(q) in (9) not the sampled function but a fitting function, namely, the numeric parametrization of the wormlike chain model19 with EV. This is a first step in achieving a full parametrization of the S(q, η, I, Z, L). C. Position of the Peak of S(q). A quantitative analysis of the peak position q* is worthwhile, since it can be directly compared with experimental results and theories. Figure 9 contains plots of q* derived from MC data as described in the previous section. q* is displayed versus η at I ) 10-3 M and at Z ) 25, Z ) 50, and Z ) 75. At higher concentrations, q* is found to scale as ∼η1/2.

-1

(12)

S(0) has been determined after fitting the total scattering function as described in section III.B. To find a possible scaling law for S(0), we note that in the case of uncharged systems with different values of b and L, universal behavior is found if S(0) is plotted against the reduced concentration c˜ ∝ c/c*, where c* ∼ Rg-3 is the overlap concentration. In presence of Debye-Hu¨ckel interactions, the Debye length appears as a characteristic length, in addition to the effective hard-sphere radius F which mimics the EV effects. Two effects now have to be taken into account: the swelling of the chains due to increased stiffness and the increase of the excluded volume effect. In rescaling the volume fraction, we have tried to take into account the first effect by replacing Rg with the sampled radius of gyration Rg(I, Zb) in the presence of charges and added salt and at infinite dilution, while the second effect is included by replacing F with F + λD. We have therefore found it reasonable to define a rescaled volume fraction for charged systems as

ηr )

( )( λD + F F

δ

)

Rg(I, Zb) 3 η Rg,u

(13)

where Rg,u is the radius of gyration of the uncharged chain (43) Fo¨rster, S.; Schmidt, M. Adv. Polym. Sci. 1995, 120, 50.

2930



Figure 10. MC data for S(0) vs rescaled concentration ηr defined by (13), with δ ) 1. Data from all sets of simulations are included. The full line is a renormalization group calculation of Otha and Ono, employed with a generalized exponent. In the semidilute regime, the scaling is S(0) ∼ η-1.8.

at infinite dilution. In Figure 10, we plot the MC data of S(0) versus ηr. The data are found to follow very closely a master curve if δ ) 1 is assumed. At high values of η, deviations are found, analogous to what is reported for similar MC simulations on uncharged systems.20 This is not surprising, considering that eq 13 does not account for screening of excluded volume effects at high concentration. Nevertheless, the simple behavior of (13) provides us with a straightforward interpretation: electrostatic interactions on S(0) can in some respect be accounted for through a simple increase of the excluded volume strength and of the stiffness of the chain. The data have also been compared with the renormalization group calculation of Otha and Ono44 for uncharged polymers,

S(0)-1 ) 1 +

(

) )

2 ln(1 + X) 1 9X - 2 + × 8 X 1 1 exp γ + 1 - 2 ln(1 + X) X X

[( (

)]

(14)

and X is a In the original work of Otha and Ono, γ ) rescaled concentration defined as X ) Rc/c*, where R is a constant of the order of unity and c* is the overlap concentration. If we assume X ) Cηr and treat both γ and C as fit parameters, we can achieve a good agreement with the MC data, as Figure 10 demonstrates. In the semidilute regime, we found 1/ , 4

S(0) ∼ η

-τ

with τ ) 1 + γ = 1.8 and C = 42.1 (15)

E. Single-Chain Scattering Function. The singlechain scattering function P(q) is a particularly useful quantity to describe the chain structure at all length scales. For uncharged polymers, it has been shown that5

P(q) ∝ q-1/ν

for 2π/Dee < q < 2π/b

(16)

where ν equals 0.5 and 0.588 for ideal and self-avoiding chains, respectively. For large q, the rodlike behavior, ν ) 1, shows up. Polyelectrolyte chains have been modeled as an assembly of rigid rods which give rise to a global ideal (or (44) Otha, T.; Ono, Y. Phys. Lett. A 1982, 89, 460.

Figure 11. Single-chain scattering function at I ) 0.001 M, Zb ) 50, and η ) 10-4 (full line), 0.01 (broken line), 0.04 (dotted line), and 0.1 (broken dotted line): (a) log-log plot and (b) Holtzer plot.

self-avoiding) chain.33,45,46 By use of arguments based on Odijk-Skolnick-Fixman (OSF) theory, the following scaling has been proposed:

P(q) ∝

{

q-1 if q > 2π/Lp q-2 if q < 2π/Lp

where Lp is the total persistence length. Due to the screening of excluded volume effects, one should expect to find ν ) 0.5 as one approaches the concentrated region. We have sampled the P(q) values at all sets of simulations performed. Figure 11 shows the effects of concentration on P(q), for I ) 10-3 M and Z ) 50. At qb = 3, the q-1 behavior is found at all concentrations. Only the data at η ) 0.2 seem to be slightly shifted. Stevens and Kremer17 have simulated a salt-free polyelectrolyte with explicit account for counterions. In that work, they emphasize that even if the chains are significantly stretched, they are by no means rodlike. Indeed, the scaling exponent found at high q values is only ν ) 0.90. This important difference could, therefore, mean that at short length scale (45) Hayter, J.; Jannink, G.; Brochard-Wyart, F.; de Gennes, P. G. J. Phys., Lett. 1980, 41, 451. (46) de Gennes, P. G.; Pincus, P.; Velasco, R. J. Phys. 1976, 37, 1461.


Langmuir, Vol. 18, No. 7, 2002 2931

the discrete nature of counterions significantly changes the chain structure, hindering the chains to reach the rod limit. At low qb, the data differ significantly from all theoretical predictions. The scaling exponent, indeed, seems to vary continuously from -1 to about the selfavoiding value of ca. -1.7 in contrast to the picture of only two values. The scaling behavior appears much more complex than theoretical predictions. F. Persistence Length. In a previously published work,23 we found the following scaling for the radius of gyration:

Rg2 ) R(L/btot)2 fBD(L/b, b/btot)

(17)

where btot is the total Kuhn length and

[ ( ) ( )]

R2(x) ) 1 +

x c1

2

+

x c2

3 /3

(18)

is the expansion factor, with c1 ) 3.12, c2 ) 8.67, and ) 0.170 and where

1 b b 2 b3 + 1 + (e-2L/b - 1) fBD(L, b) ) Lb 6 2 4L 2L

[

()

]

(19)

is the Benoit-Doty48 relation for the unperturbed meansquare radius of gyration. Assuming that this scaling also holds in the present case of many-chain systems, one can calculate the total Kuhn length by (numerically) finding the zeros of the function

Φ(L/btot, btot/b, Rg) ) Rg2 - R(L/b)2 fBD(L/btot, btot/b) (20) where Rg is the sampled data of the radius of gyration of a chain of length L/b. The data so calculated are displayed in Figure 12 against η. The Kuhn length can alternatively be determined by fitting the single-chain scattering function using the numerical parametrization for the wormlike chain model with excluded volume.19 We have also tried to compare our results with the predictions of the OSF theory. According to this theory, the Kuhn length is given by

btot ) b + bel

(21)

where

()

λB λD 2 λc

bel ) h(L/λD)

2

(22)

where λc is the distance between two neighboring charges and

h(x) ) 1 -

e-x 8 8 + x+5+ 3x 3 x

(

)

(23)

Since the simulations are in the asymptotic regime L/λD . 1, bel can be approximated by

bel =

()

λB λD 2 λc

2

(24)

As expected, the data are in reasonably good agreement with OSF theory in the dilute regime, at low ionic strength (47) Odijk, T. J. Polym. Sci., Polym. Phys. 1977, 15, 477. Skolnick, J.; Fixman, M. Macromolecules 1977, 10, 944. (48) Benoit, H.; Doty, P. J. Phys. Chem. 1953, 57, 958.

Figure 12. Total Kuhn length vs η, at Zb ) 75, 50, and 25 (from above to below) calculated by means of the scaling law (17) as described in section III.F. Note the increase of stiffness at low concentration. Data have been divided by the Kuhn length of uncharged systems (bc) at the same concentrations: (a) I ) 0.001 M, (b) I ) 0.01 M, and (c) I ) 0.1 M. The full line represents OSF theory predictions.

and at high charge density, that is, when the chain approach the rodlike behavior, in accordance with general OSF assumptions. Under the above limiting condition, the OSF scaling bel ∼ λD2 is observed (Figure 12a) and only the prefactor is slightly different from predictions.

2932


We do not observe in any of our sets of data the linear scaling as predicted by variational calculations on a single intrinsically flexible chain49-53 or by MC simulations with variable intrinsic stiffness of the chains.52,54 Neither do we observe the so-called sublinear regime evidenced in some simulation study.54 This is probably due to the fact that a fingerprint of the semiflexible nature of our chains is present at all length scales investigated. Note that OSF as well as variational calculations are “single-chain” theories, and therefore bel as given by eq 22 does not account for the presence of other chains. Only the contribution from the counterions is included through λD. For this reason, conclusions from comparison of MC data and theoretical predictions should be taken very carefully. Presently, no theory exists, to our knowledge, that is able to describe the dependence of the Kuhn length on the polyion concentration. IV. Conclusions We have performed Monte Carlo simulations on systems of many semiflexible chains with excluded volume and electrostatic interactions within the Debye-Hu¨ckel approximations. Simulations have been done at various ionic strengths, charge densities, chain lengths, and volume fractions from the dilute to the semidilute and the concentrated regime. The model used has been tailored to mimic charged wormlike micelles, using a similar approach as in our previous work on single chains.23 The radius of gyration and the end-to-end distance were sampled in the simulations. As expected, chains are stretched at low concentrations due to the electrostatic repulsion. As the concentration increases, chains contract because of the screening effects of the counterions, added salt ions, and other chains. The analysis of the shape parameter shows that chain conformations first approach the good solvent limit and, at even higher concentrations, the ideal wormlike behavior. Moreover, data clearly show that screening effects due to the presence of other chains set in at a lower concentration compared to the uncharged case. In fact, we observe that at some higher concentration charged chains are more contracted than uncharged ones under the same conditions. However, the behavior of Rg and Dee is found to be more complex than the existing scaling theories predict. The single-chain and the many-chain scattering functions were also sampled in the simulations. When interchain interactions are strong enough, the many-chain scattering function shows the characteristic structural peak, as observed in many experiments. Moreover, the many-chain scattering functions at high η show only a (49) Ha, B.-Y.; Thirumalai, D. J. Chem. Phys. 1999, 110, 7533. (50) Netz, R. R.; Orland, H. Eur. Phys. J. B 1999, 8, 81. (51) Dobrinin, A. V.; Colby, R. H.; Rubinstein, M. Macromolecules 1995, 28, 1859. (52) Ullner, M.; Jo¨nsson, Bo.; Peterson, C.; Sommelius, O.; So¨derberg, Bo. J. Chem. Phys. 1997, 107, 1279. (53) Barrat, J. L.; Joanny, J.-F. Europhys. Lett. 1993, 24, 333. (54) Micka, U.; Kremer, K. J. Phys.: Condens. Matter 1996, 8, 9463; Phys. Rev. E 1995, 54, 2653; Europhys. Lett. 1997, 38, 279.


slight dependence on the chain length for chains longer than about 3 Kuhn lengths. Therefore, relatively short chains can be used in the simulations. Using a fit function based on the PRISM theory and with an empirical expression for the direct correlation function of OrnsteinZernike, it was possible to fit the S(q) in the full range of the scattering vector and to determine the forward scattering and the position of the peak. After S(0) is plotted versus a rescaled concentration ηr, the data collapse on a single master curve. This remarkable property shows that within the approximation used, the effect of the electrostatic interactions on the osmotic compressibility can be accounted for by an increase of the excluded volume strength and of the stiffness of the chains. In the semidilute regime, S(0) decays with a power law with an exponent of ∼1.8. At higher η, when the chains enter the good solvent or the ideal regime, deviations in the scaling plot show up. The analysis of the position of the peak shows excellent agreement with both experiments and theory. In particular, the MC data reproduce the results of our experiments on charged wormlike micelles very well. In this work, we have employed a so-called continuous model in which the discrete nature of the counterion and of the added salt ions is not taken into account and a mean field potential is used. Despite the severe approximations used, general properties of the simulated systems such as the position of the peak are found to agree with experiments. The dimensions of the chains in our model system are significantly larger than for corresponding traditional polymer polyelectrolyte systems, and the charge density is in general lower. The difference in the dimensions, which is typically a factor of 5 for example for the cross-section size, makes the continuous model a much better approximation than it is for conventional polyelectrolytes. We have in this work taken the first steps in the direction of providing a fitting function which can be used in the analysis of our light scattering and small-angle neutron scattering data on charged wormlike micelles.8 We have now numerical expressions for S(0) which can be used in attempts to determine the growth laws of the micelles from the forward scattering as determined by light scattering. We have furthermore expressions for the direct correlation function c(q) and the single-chain scattering function P(q) which will allow us to aim at fitting the full q range of the measured scattering data. In this procedure, the variation of the Kuhn length with concentration and ionic strength as determined by the simulations should be employed together with the knowledge of S(0) and the growth law which we hope we can determine from the analysis of the light scattering data. Acknowledgment. The financial support of the Swiss National Science Foundation (Grant Nos. 20-53381.98 and 20-46627.96) is gratefully acknowledged. LA010884F

Monte Carlo Simulation Study of Concentration ... - ACS Publications

Recommend Documents