Simulations of Surface Forces in Polyelectrolyte Solutions - The

5116

J. Phys. Chem. B 2008, 112, 5116-5125

Simulations of Surface Forces in Polyelectrolyte Solutions Martin Turesson,*,† Clifford E. Woodward,‡ Torbjo1 rn Åkesson,† and Jan Forsman† Theoretical Chemistry, Chemical Center, Lund UniVersity, P.O. Box 124, S-221 00 Lund, Sweden, and UniVersity College, ADFA, Canberra ACT 2600, Australia ReceiVed: January 22, 2008; In Final Form: February 13, 2008

We have simulated interactions between charged surfaces in the presence of oppositely charged polyelectrolytes by coupling perturbations in the isotension ensemble to a free energy variance minimization scheme. For polymeric systems, this method completely outperforms configurationally biased versions of grand canonical simulations. Proper diffusive equilibrium between bulk and slit has been established for polyelectrolytes with up to 60 monomers per chain. A consequence of imposing diffusive equilibrium conditions, in contrast to previous more restricted models, is the possibility of surface charge inversion; ion-ion correlation and the cooperativity of monomer adsorption drive the formation of a polyion layer close to the surface, that overcompensates the nominal surface charge. This is observed even at modest surface charge densities, and leads to a build up of a long ranged electrostatic barrier. In addition, the onset of charge inversion requires very low bulk polymer densities. Due to screening effects, this leads to a higher and more long-ranged free energy barrier at low, compared to high, bulk densities. Oscillatory forces, reminiscent of those found in simple hard sphere systems, are resolved in the high concentration regime. As a consequence of a second surface charge inversion, the system “stratifies” to form a stable polyelectrolyte layer in the central part of the slit, stabilized by the adsorbed surface layers.

1 Introduction Polyelectrolytes are used in numerous industrial processes, e.q., in wastewater treatment and paper making. The purpose is to modify the colloidal stability of the various suspensions. In this work, we have used simulation methods to study the interactions between charged particles immersed in a polyelectrolyte solution. Our ultimate goal is to predict colloidal stability and relate these predictions to experimental situations. To treat these systems properly, one needs to account for the diffusion of polyelectrolyte molecules between the bulk and the region where their configurations are constrained by the colloidal particles. At true equilibrium, the bulk and constrained molecules have the same chemical potentials. In reality, experimental systems may not be at true equilibrium, due to the slow diffusion of the polyelectrolyte molecules, a property which is sometimes deliberately exploited. Nonequilibrium is also observed in more careful experiments. For example, surface force measurements in the presence of polyelectrolytes are often conducted too quickly, to allow establishment of proper diffusive equilibrium between the bulk solution and the intersurface region.1-4 This can be avoided, to some extent, by the use of relatively short chains.5 In this work, we will simulate surface forces under conditions of true diffusive equilibrium (constant chemical potential). This boundary condition is much more difficult to simulate than that of restricted equilibrium (where the number of polyelectrolyte molecules in the slit is kept constant). In a recent publication,6 we simulated a constant chemical potential polyelectrolyte model, using grand canonical ensemble Monte Carlo (GCMC), augmented with an efficient insertion * Corresponding author. E-mail: [email protected]. † Lund University. ‡ University College, ADFA.

Figure 1. Schematic picture of the model system. Positively charged polymers (curly lines) and negatively charged counterions (omitted for clarity) are allowed to diffuse between a narrow pore and the surrounding bulk solution. The pore is described by two negatively charged planar surfaces bearing a uniform surface charge density (σ) separated with a distance (h). The enhanced region depicts three adjacent monomers in a polymer. Each monomer consists of a point charge, centered in a hard sphere of diameter (d). The polymer has fixed bond lengths (b) and the angle between two adjacent bonds is denoted by (θ). The complete system, including the surface charge, is electroneutral.

algorithm. However, that work was limited to flexible, and rather short, polyions. In this paper, we present a novel version of an established surface force simulation technique utilizing perturbations7 and free energy variance minimizations8-10 in the isotension ensemble (IE). This has enabled us to simulate considerably longer chains than the GCMC would allow. The new method is also well suited for studies of polyions with more complex chain architecture, such as block structures or rings.

10.1021/jp800632e CCC: $40.75 © 2008 American Chemical Society Published on Web 04/03/2008

Surface Force Simulations in Polyelectrolytes

Figure 2. Illustration of the Rosenbluth technique for inserting the fourth monomer, and its counterion, in a trial polyelectrolyte. k number of trial insertions (here 5) for the monomer and its counterion are generated (see panel a). The Rosenbluth factor W is determined, as outlined in the text, and one set of coordinates is chosen from a Boltzmann distribution. In this example, the set of test coordinates corresponding to k ) 2 was chosen, as illustrated in panel b.

Conventional GCMC simulations require random insertions and deletions into (and from) the simulation box. This works reasonably well for fluids consisting of small molecules at intermediate density. However, if the fluid particles are polymers, the insertion probability is substantially lowered. This can be alleviated, to some extent, by using a biased insertion/deletion technique,11-13 where a chain is constructed by adding or removing monomers successively. Surface forces in solutions containing hard sphere polymers with up to 30 monomers per chain have been successfully modeled using this method.10 Furthermore, interfacial tensions of Lennard-Jones polymer fluids with up to 16 monomers per chain was recently calculated by MacDowell and Bryk.14 This was accomplished by a method that utilizes area probability distributions and biased grand canonical sampling. In principle, it should be possible to calculate surface interactions with this technique, via eq 11 in their paper, but the efficiency of this route was not explicitly tested. Greater difficulties arise if the monomers are charged. In this case, conventional GCMC is virtually impossible to perform for chains of even a few monomers. Using the configurationally biased insertion technique, mentioned above, one can successfully simulate polyelectrolyte chains of up to about 15 monomers. In this work, we will present full equilibrium force curves for charged surfaces immersed in a solution containing polyelectrolytes with up to 60 charged monomers per chain. In order to achieve this, we had to abandon the idea of altering the equilibrium concentration of polymers by changing the actual number of chains in the simulation box. Instead, we will look at an approach that varies the solution Volume. We will utilize an optimized rates method due to Bennett,8 to calculate free energy changes as they occur in the isotension ensemble. Bennett’s method aims to minimize the variance in the calculated free energy differences between two simulated systems. This technique has been successfully used in studies of simple liquids9,15,16 as well as neutral polymer solutions.10,17 We will generalize it to treat cases where the monomers and walls are charged. The basic idea of our method is to apply an external pressure parallel to the surfaces. This pressure regulates the particle density in the slit for a given separation. As the surface separation is varied, the parallel pressure is adjusted in such a way that the polyelectrolyte chemical potential is conserved. This is efficiently achieved by Bennett’s method. 2 Model and Simulations We shall adopt the Derjaguin approach, whereby the curved and charged surfaces of two large colloidal particles are modeled

J. Phys. Chem. B, Vol. 112, No. 16, 2008 5117

Figure 3. In the IE, we apply an external parallel pressure. During a simulation, the area is allowed to fluctuate. Only one of the surfaces, with initial area, S0, is shown in the figure. The process S0 f S′ describes increasing the (negatively charged) surface area by an amount ∆S ) +|σ-1|/2. Due to electroneutrality, this change is coupled with a decrease of the number of negatively charged counterions (denoted by (-) in the picture). For the process, S0 f S′′, the surface area decreases by ∆S and the number of counterions increases.

as two parallel flat infinite surfaces. Each surface carries a uniform negative charge density, σ. This is schematically depicted in Figure 1. We expect this to be a reasonable model, as long as the surface separation, h, is smaller than the diameter of the colloidal particles. We define the z direction to be perpendicular to the walls. The impenetrable surfaces are located at z ) 0 and z ) h . A polyion is described by r connected positive point charges separated with a fixed bond length, b. In this work, we have set b ) 7 Å. The system also contains simple monovalent ions of negative charge (not shown in the figure). Both monomers and counterions are centered in hard spheres of diameter d ) 4 Å. The rigidity of the chains is further modified by imposing a nonelectrostatic stiffness potential, EB. This is given by

βEB(i, i + 1) ) cos θ(i, i + 1)

(1)

where β is the inverse thermal energy and θ is the angle between the bonds i and i + 1, see Figure 1. Thus the chain stiffness is regulated by both the monomer charges and the parameter . We will use the restricted primitive model, with a uniform dielectric response of r ) 78.3. The interaction energy, U(R), between any two charges i and j separated with a distance R is thus given by:

U(R) )

{

zizje2/4π0rR R > d ∞ Red

(2)

where zi is the valency of species i and e is the elementary charge. The permittivity of vacuum is denoted by 0. All species interact with the uniformly charged planar surfaces. We denote the number of simulated polyions by N. These are confined in a simulation box with a side length L and a surface area of S ) L2. Periodic boundary conditions are applied in the x-y directions (parallel to the surfaces). A correction for the long ranged electrostatic potential is included in a standard fashion.18 In addition, all charges sense a soft repulsive potential, Vex, from the walls:

Vex ) w(z) + w(h - z) where

βw(z) ) 405 where we have set τ ) 1 Å.

( ) e-z/τ (z/τ)3

(3)

(4)

5118 J. Phys. Chem. B, Vol. 112, No. 16, 2008

Turesson et al.

Figure 4. Surface interactions for systems containing (charged) 14-mers. In the GCMC and IE cases, the system is in equilibrium with a bulk solution, for which Fm ) 10 mM. The “canonical” system is a reference, where the number of polyions are chosen to exactly counteract the surface charge, i.e., in this case, Fm ) 0. The IE simulations were performed with 100 14-mers at all separations (this number varies in the GCMC approach). (a) Osmotic pressure. With the IE approach, this can be obtained either via the virial expression, eq 5, or via differentiation of the grand potential, eq 16. (b) Interaction free energy, as obtained by integration of the GCMC data (eq 8) of the osmotic pressure, together with a typical van der Waals interaction as a comparison.

Additionally, we assume that the region between the surfaces is in equilibrium with an infinitely large bulk solution, i.e., we assume that there are always more chains present in the system than required to neutralize the surfaces. This is not the case in many flocculation studies, in which the polyelectrolyte is gradually added to a suspension of charged particles. In this sense, our model system is more reminiscent of other experimental setups, such as the surface force apparatus, atomic force microscopy, or the thin film pressure balance. Our primary aim in this study is to calculate how the osmotic pressure, Posm varies with surface separation. If we consider forces acting across the surface at z ) h, we can formulate Posm in terms of two separate contributions:

Posm ) Psoft + PMaxwell

(5)

where the first term stems from the nonelectrostatic soft wall repulsion:

Psoft ) -

∫0h Ft(z) (

)

dw(h - z) dz dh

(6)

where Ft(z) is the total particle concentration at z. The other contribution is calculated from electrostatic interactions between the surface charges of the wall and the rest of the system. For a uniform surface charge, this so-called Maxwell term is given by:

PMaxwell ) -

σ2 2r0

(7)

From Posm, we can obtain the interaction free energy per unit area, gs, as

gs(h) )

∫h∞ (Posm(h′) - Pb) dh′

(8)

Here, Pb is the osmotic pressure in the bulk, which can be estimated from Posm at large surface separations. Furthermore, at large separations, gs equals twice the surface tension at a single surface, γs. Hence, we shall report the net free energy per unit area as ∆gs(h) ≡ gs(h) - 2γs. In this model, strong adsorption of polyions to the negatively charged surfaces can lead to an overcompensation (charge inversion) of the bare surface charge. It is therefore instructive to introduce an apparent surface charge density, σapp, defined as:

σapp(z) ) σ +

∫0z ∑ FR(z′)qRdz′, z e h/2

(9)

R

where FR and qR are the density and charge, respectively, of species R. The sum runs over all monomers and counterions in the system. The symmetry of the system allows us to limit the integration to the left half cell. As a consequence, σapp(h/2) ) 0. If σapp > 0, we shall denote its maximum value as σapp(max). In the absence of charge inversion, σapp e 0. The chain configurations are sampled by three different types of Monte Carlo moves: translation of a whole polymer, reptation, and crank shaft moves. In the reptation steps, one end monomer of a randomly chosen chain is detached and then placed, at a random angle, at the other end of the polymer. In the presence of the bond angle potential, we efficiently sample the new angle from a Boltzmann distribution,19 where the weight function is given by eq 1. It should be noted that we also tried two other types of trial moves, namely, the configurationally biased trial moves20,21 and the so-called worm-hole moves.22 We found no evidence that either of these moves provided statistical improvement. Hence, they were discarded. With open boundary conditions, the chemical potential of all species are to remain constant as the separation between the surfaces is varied. As noted earlier, the challenging part is to keep the chemical potential of the polyelectrolyte constant. In order to accomplish this, the simulation needs to include moves that will vary the concentration of polyelectrolyte molecules in the simulation box. In this work, we will introduce a modified version of an optimized rates method7,9,10 applied within the isotension ensemble (IE). Comparisons will be made, for short chains, with results obtained via GCMC simulations.6 The GCMC method utilizes the configurationally biased stepwise insertion/deletion technique, which was originally developed by Rosenbluth,11 and further developed by Siepmann and Frenkel.13 We will describe both methods in more detail below. 2.1. Grand Canonical Ensemble Monte Carlo. Conventional GCMC has successfully been used for open systems containing simple particles and electrolytes. Random insertions and deletions are accepted or rejected in order to establish the correct density in the slit, as determined by the input chemical potential. The problem with using GCMC on chain molecules is that the insertion probability becomes vanishingly small even for fairly low degrees of polymerization. A method to partly circumvent this problem was first reported by Rosenbluth and Rosenbluth.11 A chain is incrementally grown via a number of



Figure 5. (a) Interaction free energies between two charged surfaces immersed in 30-mer (r ) 30) polyelectrolyte solutions, for which Fm ≈ 10 mM. The data represented by circles are obtained using flexible chains ( ) 0), while the squares are data from simulations of semiflexible chains ( ) 6). (b) The variation of the maximum apparent surface charged density (normalized by the nominal surface charge density) with separation for the system described in (a).

k trial positions for each monomer, from which a Boltzmann weighted choice is made. This bias is corrected for at the end of the insertion step. The theory has been thoroughly described and implemented for neutral chain molecules.10 The algorithm was recently generalized for insertion of polyelectrolytes.6 The idea is to insert an electroneutral pair, formed by a monomer and a counterion, until the entire polyelectrolyte consisting of r monomeric units and r monovalent counterions has been created. Chain removals are performed by reversing the process. The method was successfully tested for relatively short polyions (r < 15). A schematic picture of the insertion process is depicted in Figure 2a. Here, k ()5) sets of “virtual” pairs consisting of a monomer and a counterion are generated. The positions of the virtual monomers are randomly generated on the surface of a sphere of radius b centered at the coordinate of the last successfully inserted monomer. For semiflexible polymers, the Boltzmann weights also include the bending energy, EB, eq 1. The k virtual coordinates for the counterion are generated randomly in the simulation box. For each set of (monomer/counterion) pairs, the Boltzmann factor, e-βUi, is calculated, where Ui is the change of the potential energy of the system, due to the virtual insertion of the pair. One set of pairs is then selected with a probability

Pi )

e-βUi wi

(10)

where

A ) µN - P|| Sh

k

wi )

e-βU ∑ j)1

j

(11)

At the end of an insertion step (N f N + 1), when the whole trial polyelectrolyte of r monomers and r counterions has been constructed, the Rosenbluth factor is calculated: r

WN+1 )

wi

∏ i)1 k

(12)

The polyelectrolyte insertion is then being accepted or rejected according to the following criterion:

acc(N f N + 1) )

[

of monovalent counterions, and µ the chemical potential of the polyelectrolyte salt. The deletion step involves the analogous reverse procedure.10 However, for chains with r > 15, this method becomes computationally intractable. We have therefore developed an alternative approach to the problem based on the isotension ensemble. This is discussed in the next section. 2.2. The Isotension Ensemble. The use of Bennett’s method, in the isotension ensemble (IE),7-9 has proven to be a very efficient way to simulate equilibrium surface interactions in polymer solutions.10 In the IE, the area is allowed to fluctuate, subject to an externally applied pressure acting parallel to the surfaces. Our approach is to establish how this parallel pressure (and thus the density) should change with surface separation, by coupling surface area fluctuations and distributions from virtual separation perturbations, to a variance minimization scheme. For a thorough and more general description of the technique, we refer to earlier works.9,10 In order to treat a system with charges, it is useful to work through the thermodynamics and show that one can map charged system variables onto equivalent counterparts in the neutral case. One can then apply our established methods for neutral systems, with appropriate modifications, in order to accommodate these new variables. Consider a simple uncharged system of N polymers, at a temperature T, confined between two walls, each with a surface area S and a separation h. A simple scaling argument shows that the Helmholtz free energy, A(N, S, T, h), can be written as

min 1, WN+1e βµ

V

r

∏ N + 1 i)1 N

V

-

]

+i

(13)

where V is the volume of the simulation box, N- the number

(14)

where µ is the polymer chemical potential. The area averaged pressure acting parallel to the surfaces, P||, is given by

P|| ) -

1 ∂A h ∂S T,N,h

( )

(15)

In the isotension ensemble, the area S fluctuates against a fixed (input) value of P||. Hence, the free energy in the IE is given by the Legendre transform of the Helmholtz free energy, A + P||hS ) µN. This is the appropriate generalization of the Gibbs free energy for this nonuniform system. Thus, for constant N and at fixed chemical potential, the free energy of the IE remains constant as h varies. This fact allows us to use Bennett’s method to estimate free energy differences within the IE in order to simulate open systems without changing the number of polymer molecules. The quantity -P||hS is the grand potential, and thus the equilibrium osmotic pressure in the open system can be expressed as:


Turesson et al. and eq 21 can now be written as

A ) µ˜ pNp - P ˜ ||Sh

Figure 6. Density profiles for monomers, Fm, (opaque symbols) and counterions, Fci, (filled symbols) outside the left surface (at z ) 0). Distributions for two different intrinsic chain stiffnesses, ) 0 and 6, are shown. The inset more clearly reveals the enhanced adsorption of counterions when the polyions are flexible.

Posm )

d(P||(h)h) dh

(16)

The osmotic pressure can also be evaluated from an expression similar to eq 5.10 In addition, using eqs 8 and 16, the interaction free energy (per unit area) between the surfaces can be obtained as

gs(h) ) (Pb - P|| (h))h

(17)

The bulk pressure, Pb, can be estimated from the osmotic pressure at large h. In a system with charged walls and charged species, we first note that the Helmholtz free energy is given by

A ) 2σSψs + N-µ- + Npµp - P||hS

(18)

where N- is the number of counterions, while the chemical potential of a counterion is denoted µ-. The number of polyions is denoted Np and their chemical potential by µp. The surface potential is given by ψs. Electroneutrality implies that

2σS + Npzp + N-z- ) 0

(19)

where zp and z- are the valencies of monomers and counterions, respectively. For fixed σ, the electroneutrality condition means that N- and Np are dependent variables. In our simulations, Np will be fixed. Hence, solving for N- and inserting in eq 18 gives:

A ) 2σSψs -

2σSµ- Npzpµ+ Npµp - P||hS (20) zz-

Equation 20 can be written as

where

µ˜ p ) µp -

zp µ z- -

(22)

In eq 21, we have introduced an “apparent” parallel pressure, P ˜ ||, and polymer chemical potential, µ˜ p. Thus the charged system can be effectively mapped onto an equivalent neutral system,

(23)

Note that µ-/z- - ψs is constant for a system in equilibrium with a bulk solution. Simulations of the charged system in the IE allows fluctuations in S against a fixed (input) value for P ˜ ||. Note that both surfaces fluctuate together so as to maintain the simulation box symmetry. Changing S at fixed σ will alter the total surface charge in the simulation box. This will cause a change in the free energy, due to the interaction with the surface potential. Additionally, electroneutrality implies that the change in the total surface charge must be balanced by a change in the number of counterions (as Np is constant), giving rise to a chemical potential contribution to the free energy. These two contributions to the free energy fluctuations are in additional to the intrinsic pressure volume term; this is reflected in the expression for P ˜ ||. During an actual simulation, we are not free to change the surface area in a continuous manner, due to the integral nature of the counterion valency. Thus, we have only allowed the surface areas to fluctuate in discrete steps, leading to unit changes in the total surface charge. That is, each surface fluctuates in steps of ∆S ) (σ-1/2. A schematic of the area fluctuations are depicted in Figure 3. In an open system, polyions will diffuse to and from the region between the surfaces. The electroneutrality condition means that they will always be accompanied by a chargebalancing number of counterions. Thus, we require the chemical potential combination, µ˜ p, to remain constant in this system. In order to maintain a fixed µ˜ p as the surface separation is changed, we used an IE-adapted version of Bennett’s optimized rates method.9 During a typical simulation, the surface separation is changed, h f [h ( ∆h] and the IE free energy difference is calculated via Bennett’s method. We use the sampled histograms to efficiently obtain the corresponding change in P ˜ || required to make this free energy difference zero.10 The separation is then set to h ( ∆h and the procedure is repeated until the desired range of separations has been spanned. In this work, ∆h varied between 1 and 4 Å. Additional simple electrolyte could in principle be included in the model, in which case they would exchange with a bulk solution via conventional grand canonical steps. In such cases, the electroneutrality condition associated with an area fluctuation can be effected by either adding or deleting simple mobile ions. In this work, we shall restrict ourselves to cases where no simple salt is present, apart from counterions. Effects from salt addition were investigated (for short chains using GCMC) in an earlier work.6 2.3. Comments on the Methods. The IE method completely outperforms the GCMC for long polyelectrolytes. The GCMC approach (even with the Rosenbluth method for insertions/ deletions) becomes impractical for chain lengths exceeding about r ) 15 in the systems we have investigated. With the combination of IE and Bennett’s optimized rates method, however, it is possible to obtain equilibrium surface interactions in the presence of 60-mers. Furthermore, more complex chain architectures, such as ring polymers and block copolymers, can easily be modeled using the IE approach. Another advantage with the IE method is that there is a useful consistency check via the osmotic pressure, which can be obtained using either of eqs 16 or 5. Even minor programming errors generally lead to significant discrepancies between these expressions. A drawback with the IE approach is that it is an iterative process, where the maximum possible step length in separation, ∆h, is limited by



Figure 7. (a) Surface free energies for polyelectrolyte solutions with two different intrinsic chain stiffnesses ( ) 0 and 6) equilibrated with a bulk with Fm ≈ 100 mM (r ) 30). (b) Maximum apparent surface charged density observed as a function of separation for the system in (a).

Figure 8. Apparent surface charge density curves for the 30-mer polyelectrolyte solution with ) 0 and Fm ) 100 mM. (a) At a separation of about 35 Å, there is a crossover from zero to a single surface charge inversion. (b) At h ) 120 Å, a transition from a single to a double charge inversion is observed.

the criterion of configurational overlap between the two systems at adjacent wall separations. Furthermore, one needs to simulate a complete surface force curve up to separations at which bulk conditions can be reliably estimated. A complete surface force curve typically takes about 1 week to obtain on a dual (2.5 GHz PowerPC) processor. On the other hand, GCMC simulations at different surface separations, and in the bulk, can be carried out independent of each other. Such simulations can be (trivially) run in parallel. This notwithstanding, for long chains, GCMC is not a realistic option. An expanded ensemble approach, whereby a chain is allowed to grow or shrink during the simulation,23 might facilitate the establishment of slit-bulk equilibrium, even for the chain lengths we have investigated in this work. However, counting only those configurations which have complete chains is likely to be inefficient. If configurations with partial chains are allowed to contribute, this will in principle generate errors of thermodynamic quantities, such as the osmotic pressure. This error will decrease as the number of chains increase, with the risk of a considerable system size dependence. A possible way to circumvent this problem would be to use an expanded ensemble method to generate the correct density, at a given separation, and then calculate the corresponding osmotic pressure in a separate canonical simulation. However, this would be a cumbersome and, most likely, slow method. As we shall see, accurate osmotic pressure calculations alone require a considerable computational effort. Expression 16, which is directly obtained with the IE method, has superior statistical properties.

3 Results All results reported are given for systems for which the surface charge density is e|σ-1| ) 200 Å2, where e is the elementary charge. The concentration of monomers in the bulk is denoted Fm. Let us first recall a typical interaction curve between charged surfaces in the presence of neutralizing polyions only. Such systems are computationally much easier to run, as there is no exchange with a bulk. The number of polyelectrolyte molecules is fixed (at all separations), due to electroneutrality. Thus, the complete surface force curve can be simply obtained via a sequence of canonical simulations. In part, for this reason, this system has been the focus of several previous simulation studies.1,24,25 An example is provided in Figure 4a, where the canonical approach leads to attraction (within the statistical precision). At very short separations (≈15 Å), there is a steep repulsion. The observed attraction is (for flexible chains) at first hand due to chains forming bridges across the slit from one surface to the other.24,26,27 This is primarily an entropic effect. As the surfaces approach, the chains are configurationally relaxed. If the chain rigidity is increased, the net attraction is gradually dominated by a strong, and more short-ranged, ion-ion correlation attraction, concomitant with strong parallel alignment of polymer chains at the surfaces.25 When the polyelectrolyte salt is added, and diffusive equilibrium is maintained between a bulk and the slit region, the picture changes quite dramatically. Instead of attraction, a repulsive free energy barrier develops, as shown in Figure 4a.


Figure 9. Cartoon of the extrema in the free energy curve corresponding to a vanishing net osmotic pressure, Pnet. At free energy maxima, we observe crossovers from zero to single and single to double charge inversion (denoted by 0, 1, 2). At free energy minima, we find corresponding maxima in σapp(max).

Figure 10. Interaction free energy for a polyelectrolyte system (r ) 60). The chains are flexible, i.e., ) 0. A closer inspection of the curve (see inset) reveals oscillations in the free energy.

In this figure, we show three osmotic pressure curves for the counterion system calculated using three different approaches. One of these results from GCMC simulations and the other two stem from IE simulations using either eq 5 or eq 16 to calculate the osmotic pressure. Note that in the charged case, the latter expression is given by d(P ˜ ||h)/dh. The excellent agreement between the approaches confirms the validity of all the calculated force curves. Figure 4b shows how the interaction free energy varies with separation, as obtained by using eq 8, and the GCMC data in Figure 4a. A van der Waals (vdW) curve with a typical Hamaker constant of 10-20 J is included in order to illustrate the magnitude of the repulsive barrier. The origin of this barrier is electrostatic in nature. In an open electrolyte system with charged surfaces in contact with a (large) bulk solution, it is in principle always possible to accumulate oppositely charged particles in the vicinity of the surfaces to an extent that exceeds the mere neutralization criterion.28,29 This implies a change of sign of σapp(z), as defined in eq 9. When the mobile charges are connected to form a polymer, cooperativity effects and ion-ion correlations enhance the propensity of the system to display overcharging. Even with relatively short chains, the adsorption is generally strong enough to effectively reverse the sign of the surfaces (as measured some distance away).30 When two such overcharged surfaces approach, from large separations, they are repelled by the simple counterions to the polyions via an ordinary double layer mechanism. At

Turesson et al.

Figure 11. h ) 30 Å. The chain bridges efficiently and creates a strong attraction.

short separations, the surfaces are no longer overcharged, and an attractive “well”, similar to that found in the counterionfree case, is observed. The position of the maximum of the free energy barrier is found to almost perfectly match the onset of charge reversal. This separation appears, for the cases studied so far, to be independent of chain length and insensitive to the addition of a simple salt. However, it does vary with surface charge density. Once the repulsive electrostatic contribution becomes small, analysis of the pressure contributions at the mid plane indicate that electrostatic correlation provides an important contribution to the onset of the attraction. All of these effects have been observed in an earlier work.6 Here, we shall expand the scope of that study to include longer chains and also to investigate the effects of nonelectrostatic chain stiffness as well as polyelectrolyte concentration. 3.1. 30-Mers. We will begin by exploring the free energy of interaction displayed by a system containing polyions with r ) 30. Two different bulk monomer concentrations have been investigated, Fm ≈ 10 and 100 mM. 3.1.1. Low Density. First, we consider the case of flexible chains ( ) 0) at a bulk monomer concentration of roughly 10 mM. The number of polyions was fixed at 80. As shown in Figure 5a, there is a long-ranged double layer repulsion, extending to about 300 Å. The maximum of σapp(z) is plotted as a function of separation in Figure 5b. At short distances, we observe no surface charge inversion. As noted above (and in a previous study6), the onset of charge inversion occurs close to the separation where the free energy barrier has a maximum. With flexible chains, the charge inversion seems to level off at a value corresponding to roughly 15% of the absolute value of the nominal surface charge density. The stiffness parameter introduced in eq 1 allows us to regulate the intrinsic rigidity of the polyelectrolyte. The results for ) 6 are also shown in Figure 5. The average chain end to end distance is about 30% larger for the more rigid chains. We see that the maximum repulsion is diminished compared to the flexible case and the position of the maximum is shifted to a slightly shorter separation. This behavior is in agreement with earlier results on neutral systems.25 It is also consistent with the lower degree of charge reversal obtained with stiffer chains, as seen in Figure 5b. Here, we note that the maximum charge reversal is reduced to about 10% of the absolute value of the nominal surface charge density. In order to better understand the behavior of the apparent charge, we plotted in Figure 6 monomer and counterion density profiles for systems containing flexible and semi-flexible chains. The thickness of the adsorbed



Figure 12. h ) 98 Å. The chain lies in a flat configuration at one of the charged surfaces. Left picture, front view; right picture, side view.

monomer density layer is smaller for the stiffer polymers, though the maximum in the profile is larger. This is consistent with the tendency for stiffer polymers to form “flat” configurations at the surface. As indicated by the counterion profile, it is also harder for counterions to screen the monomer interactions in stiff chains, as they must do so in a narrower region, which is electrostatically and entropically unfavorable. In the case of flexible chains, we see a much broader monomer profile, albeit with a smaller maximum. However, more efficient screening by counterions leads to stronger net adsorption and a higher degree of charge inversion. 3.1.2. High Density. Upon increasing the bulk monomer concentration to about 100 mM, a slightly different picture emerges. In Figure 7a, we see that the maximum height of the electrostatic barrier has been reduced by a factor of roughly two compared to the low-density calculations. In addition, the range of the repulsive interaction is reduced as well. From the charge inversion analysis, shown in Figure 7b, we observe that the bulk values of the apparent surface charges are similar to those at lower concentration. Also, as in the low-density case, the onset of surface charge inversion occurs at the position of the maximum of the first free energy barrier. As seen in Figure 7a, there is a small shift between the maxima for flexible and stiff chains. This is manifested as a slight difference between the onsets of their surface charge inversions. This is clearly seen in the inset of Figure 7b, where the crossover to charge inversion occurs at h ) 34 Å for flexible chains and at h ) 29 Å for stiff chains. Complete apparent charge curves for h ) 33 Å and h ) 35 Å ( ) 0) are plotted in Figure 8a. In contrast to Figure 5b, the value of the maximum charge inversion does not monotonically increase toward a bulk value (Figure 7b). This is most evident for the stiff chains, where we detect a distinct maximum in σapp(max). It turns out that the position of this maximum, perfectly match the transition from a repulsive to an attractive net osmotic pressure. This corresponds to the free energy minimum (at h ≈ 80 Å). At larger separations, σapp(max) decreases until the net osmotic pressure once again turns repulsive. Here, the system starts to display a double charge inversion. Thus, the surface charge effectively has the same sign as the bare surface (but smaller in magnitude). In this case, this occurs at about 120 Å for flexible chains (see Figure 8b). Figure 9 is a schematic picture, highlighting the relations between maxima in the free energy and the corresponding transitions of the apparent surface charge. These findings have been confirmed by calculations using a correlation corrected density functional theory (not shown).

3.2. 60-Mers. In this section, we present results from simulations of a polyelectrolyte solution containing flexible ( ) 0) polyions carrying 60 monomers per chain. The system is in equilibrium with a bulk solution of Fm ≈ 100 mM. The interaction free energy curve is given in Figure 10. From the inset, we see that we are able to resolve a number of free energy maxima/minima. As discussed in the previous section, these are related to corresponding oscillations of the apparent surface charge. The crossovers from zero to a single and a single to a double charge inversion match the positions of the free energy maxima very well. In this case, there is also a transition from the second to a third charge inversion. The exact position is hard to localize due to noise, but at separations above 250 Å, we clearly see it. In a region with a double surface charge inversion, the system “stratifies” into three polyelectrolyte layers. As the surfaces are being compressed, the intermediate layer is eventually squeezed out leading to a “local” depletion attraction. Stratification forces in polyelectrolyte systems have been measured extensively experimentally using atomic force microscopy31 and the thin film pressure balance.32-34 In the latter technique, the external pressure is regulated and the stability of the so-called free-standing films is changed in a stepwise fashion, giving rise to oscillatory forces. One can make comparisons with a simple hard sphere system where oscillatory surface forces, due to excluded volume effects, are seen at short separations. The packing of spheres is, however, in our case replaced with a more long ranged soft electrostatic packing of chains. The period of the observed oscillation is about 100 Å, which corresponds well to the bulk correlation length in a dilute system ()F-1/3 polyion). Figures 11 to 13 show snapshots from a simulation of this system. Imagine that we look at the system in the direction perpendicular to the z axis. One of the chains (arbitrarily chosen) has been labeled by making it thicker and black. We follow the configuration of this labeled chain as the separation changes from h ) 30 to189 Å. At short separations (h ) 30 Å, Figure 11), the chain is seen to coil back and forth between the surfaces, creating a strong bridging attraction. As the separation increases, the propensity to form intersurface bridges is reduced and the chain is found to adsorb in a flat configuration at one surface (h ) 98 Å, Figure 12). In order to highlight this, a side view of the system has been added. Upon increasing the separation further (h ) 125 Å, left picture in Figure 13), we observe that many chains, or significant parts


Turesson et al.

Figure 13. Left picture, h ) 125 Å. One tail of the chain has desorbed from the surface. Right picture, h ) 189 Å. “Stratification” of the system has occurred. The labeled chain now resides in the layer created in the middle of the slit.

Figure 14. A comparison between surface free energy curves obtained with a polyion (r ) 14) and a spherical macroion, with a radius of 9 Å carrying the same charge. The counterion bulk concentration (Fci) was in both cases 10 mM.

of them, have detached from the surfaces. For instance, the configuration of the labeled chain has a tail sticking out from the surface. At the largest distance depicted (h ) 189 Å, right picture in Figure 13), the labeled chain is totally desorbed from the surface and a layer of polymers in the central part of the slit is created. This “stratification” of the system occurs as a consequence of surface charge reversal and causes the free energy to oscillate (see inset of Figure 10). 3.3. Spherical Macroions. If we replace the polymeric 14mer with 14-point charges centered in a sphere of radius 9 Å (roughly the radius of gyration), we observe very similar interaction curves, at least at large separations. A comparison is provided in Figure 14. At close separations the curves differ, as seen in the inset of Figure 14. Snapshots from the simulations reveal strong spatial correlations between the spheres in the spherical macroion system, yielding a strong attraction at short separations. It should be emphasized that our comparison is limited to a single case. Still, the agreement is interesting and might be relevant to explore more thoroughly in a future work. 4 Conclusions In this work, full equilibrium free energy interactions between charged surfaces in polyelectrolyte solutions have been successfully simulated by using a modified version of the isotension ensemble coupled to Bennett’s free energy variance minimization scheme. This technique is superior to grand canonical methods. The major advantage is that insertions and deletions

of polyions and counterions are avoided. This means that we are able to simulate much longer polymers without restrictions on the chain structures. Cooperativity effects may allow polyions to adsorb at oppositely charged surfaces strongly enough to effectively reverse the sign of their charge. The nature of the surface interactions depend in an interesting manner upon the bulk monomer concentration. At low concentrations, a long-ranged repulsive free energy barrier is observed. In a higher concentration regime, multiple surface charge inversions cause the free energy to oscillate. At a double surface charge inversion, the system “stratifies” into three layers: one layer at each surface and one “free-standing” polyion layer in the central portion of the intersurface region. These structure oscillations appear to be simply related to a corresponding behavior of the interaction free energy. In the presence of more rigid chains, the height of the repulsive barrier is reduced, and the repulsive barrier is more short-ranged. While stiffer polyions adsorb in a more flat configuration than corresponding flexible ones, the latter are able to more efficiently pack at the surface, thereby generating a stronger degree of overcharging. A strong similarity was found between surface interactions in the presence of chain-like polyions and spherical macroions, although the latter mediate a stronger correlation-dominated attraction at short range. In all cases studied, we find a repulsive barrier. This seems to contradict the fact that polyelectrolytes sometimes are used as flocculants. However, experimental studies have shown that the concentration range in which aggregation occurs is quite narrow.35-39 In fact, the optimum dosage, where flocculation is the highest, corresponds very well to a perfect neutralized system (i.e., zero electrophoretic mobility).40-42 This is in agreement with our findings, where the largest attraction is found in the canonical system, where there is no possibility of overcharging the surfaces. It should be noted that in practical applications nonequilibrium effects are substantial and the flocculation rate may depend on a lot of different factors. References and Notes (1) Dahlgren, M. A. G.; Waltermo, Å.; Blomberg, E.; Claesson, P. M.; Sjo¨stro¨m, L.; Åkesson, T.; Jo¨nsson, B. J. Phys. Chem. 1993, 97, 11769. (2) Dahlgren, M. A. G.; Hollenberg, H. C. M.; Claesson, P. M. Langmuir 1995, 11, 4480. (3) Biggs, S.; Proud, A. D. Langmuir 1997, 13, 7202.

Surface Force Simulations in Polyelectrolytes (4) Poptoshev, E.; Rutland, M. W.; Claesson, P. M. Langmuir 1999, 15, 7789. (5) Maurdev, G.; Meagher, L.; Ennis, J.; Gee, L. M. Macromolecules 2001, 34, 4151. (6) Turesson, M.; Åkesson, T.; Forsman, J. Langmuir 2007, 23, 9555. (7) Svensson, B.; Woodward, C. E. J. Chem. Phys. 1994, 100, 4575. (8) Bennett, C. H. J. Comput. Phys. 1976, 22, 245. (9) Forsman, J.; Woodward, C. E. Mol. Sim. 1997, 19, 85. (10) Turesson, M.; Forsman, J.; Åkesson, T. Phys. ReV. E 2007, 76, 021801. (11) Rosenbluth, M. N.; Rosenbluth., A. W. J. Chem. Phys. 1955, 23, 356. (12) Harris, J.; Rice, S. J. Chem. Phys. 1988, 88, 1298. (13) Siepmann, J.; Frenkel, D. Mol. Phys. 1992, 75, 59. (14) MacDowell, L. G.; Bryk, P. Phys. ReV. E. 2007, 75, 061609. (15) Deitrick, G. L.; Scriven, L. E.; Davis, H. T. J. Chem. Phys. 1989, 90, 2370. (16) Errington, J. R.; Kofke, D. A. J. Chem. Phys. 2007, 127, 174709. (17) Czezowski, A.; Woodward, C. E. Comput. Phys. Commmun. 2001, 142, 117. (18) Torrie, G. M.; Vallueau, J. P. J. Phys. Chem. 1982, 86, 3251. (19) Rubinstein, R. Y. Simulation and Monte Carlo methods; Wiley: New York, 1981. (20) Siepmann, J. Mol. Phys. 1990, 70, 1145. (21) de Pablo, J.; Laso, M.; Suter, U. J. Chem. Phys. 1992, 96, 2395. (22) Houdayer, J. J. Chem. Phys. 2002, 116, 1783. (23) Escobedo, F. A.; Suen, J. K. C. J. Chem. Phys. 1995, 103, 2703.

J. Phys. Chem. B, Vol. 112, No. 16, 2008 5125 (24) Åkesson, T.; Woodward, C. E.; Jo¨nsson, B. J. Chem. Phys. 1989, 91, 2461. (25) Turesson, M.; Åkesson, T.; Forsman, J. Langmuir 2006, 22, 5734. (26) van Opheusden, J. H. J. J. Phys. A.: Math. Gen. 1988, 21, 2739. (27) Podgornik, R. Chem. Phys. Lett. 1990, 174, 191. (28) van Megan, W.; Snook, I. J. Chem. Phys. 1980, 73, 4656. (29) Messina, R. Macromolecules 2004, 37, 621. (30) Sjo¨stro¨m, L.; Åkesson, T.; Jo¨nsson, B. Ber. Bunsen-Ges. 1996, 100, 889. (31) Milling, A. J. J. Phys. Chem. 1996, 100, 8986. (32) Klitzing, R.; Espert, A.; Asnacios, A.; Hellweg, T.; Colin, A.; Langevin, D. Colloids Surf. 1999, 149, 131. (33) Klapp, S.; Qu, D.; v. Klitzing, R. J. Phys. Chem. B 2007, 111, 1296. (34) Bergeron, V.; Langevin, D.; Asnacios, A. Langmuir 1996, 12, 1550-1556. (35) Dixon, J. K.; LaMer, V. K.; Li, C.; Messinger, S.; Linford, H. B. J. Colloid Interface Sci. 1967, 23, 465. (36) Vincent, B. AdV. Colloid Interface Sci. 1974, 4, 193. (37) Zhang, J.; Buffle, J. J. Colloid Interface Sci. 1995, 174, 500. (38) Walker, H. W.; Grant, S. B. Colloids Surf., A 1996, 119, 229. (39) Feretti, R.; Zhang, J.; Buffle, J. Colloids Surf., A 1997, 121, 203. (40) Bouyer, F.; Robben, A.; Yu, W.; Borkovec, M. Langmuir 2001, 17, 5225. (41) Kleimann, J.; Gehin-Delval, C.; Auweter, H.; Borkovec, M. Langmuir 2005, 21, 3688. (42) Gregory, J. Colloids Surf. 1988, 31, 231.

Simulations of Surface Forces in Polyelectrolyte Solutions - The

Recommend Documents